1 2 /* 3 * ==================================================== 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5 * 6 * Developed at SunSoft, a Sun Microsystems, Inc. business. 7 * Permission to use, copy, modify, and distribute this 8 * software is freely granted, provided that this notice 9 * is preserved. 10 * ==================================================== 11 */ 12 13 /* 14 * __kernel_rem_pio2(x,y,e0,nx,prec) 15 * double x[],y[]; int e0,nx,prec; 16 * 17 * __kernel_rem_pio2 return the last three digits of N with 18 * y = x - N*pi/2 19 * so that |y| < pi/2. 20 * 21 * The method is to compute the integer (mod 8) and fraction parts of 22 * (2/pi)*x without doing the full multiplication. In general we 23 * skip the part of the product that are known to be a huge integer ( 24 * more accurately, = 0 mod 8 ). Thus the number of operations are 25 * independent of the exponent of the input. 26 * 27 * (2/pi) is represented by an array of 24-bit integers in ipio2[]. 28 * 29 * Input parameters: 30 * x[] The input value (must be positive) is broken into nx 31 * pieces of 24-bit integers in double precision format. 32 * x[i] will be the i-th 24 bit of x. The scaled exponent 33 * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 34 * match x's up to 24 bits. 35 * 36 * Example of breaking a double positive z into x[0]+x[1]+x[2]: 37 * e0 = ilogb(z)-23 38 * z = scalbn(z,-e0) 39 * for i = 0,1,2 40 * x[i] = floor(z) 41 * z = (z-x[i])*2**24 42 * 43 * 44 * y[] output result in an array of double precision numbers. 45 * The dimension of y[] is: 46 * 24-bit precision 1 47 * 53-bit precision 2 48 * 64-bit precision 2 49 * 113-bit precision 3 50 * The actual value is the sum of them. Thus for 113-bit 51 * precision, one may have to do something like: 52 * 53 * long double t,w,r_head, r_tail; 54 * t = (long double)y[2] + (long double)y[1]; 55 * w = (long double)y[0]; 56 * r_head = t+w; 57 * r_tail = w - (r_head - t); 58 * 59 * e0 The exponent of x[0]. Must be <= 16360 or you need to 60 * expand the ipio2 table. 61 * 62 * nx dimension of x[] 63 * 64 * prec an integer indicating the precision: 65 * 0 24 bits (single) 66 * 1 53 bits (double) 67 * 2 64 bits (extended) 68 * 3 113 bits (quad) 69 * 70 * External function: 71 * double scalbn(), floor(); 72 * 73 * 74 * Here is the description of some local variables: 75 * 76 * jk jk+1 is the initial number of terms of ipio2[] needed 77 * in the computation. The minimum and recommended value 78 * for jk is 3,4,4,6 for single, double, extended, and quad. 79 * jk+1 must be 2 larger than you might expect so that our 80 * recomputation test works. (Up to 24 bits in the integer 81 * part (the 24 bits of it that we compute) and 23 bits in 82 * the fraction part may be lost to cancellation before we 83 * recompute.) 84 * 85 * jz local integer variable indicating the number of 86 * terms of ipio2[] used. 87 * 88 * jx nx - 1 89 * 90 * jv index for pointing to the suitable ipio2[] for the 91 * computation. In general, we want 92 * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8 93 * is an integer. Thus 94 * e0-3-24*jv >= 0 or (e0-3)/24 >= jv 95 * Hence jv = max(0,(e0-3)/24). 96 * 97 * jp jp+1 is the number of terms in PIo2[] needed, jp = jk. 98 * 99 * q[] double array with integral value, representing the 100 * 24-bits chunk of the product of x and 2/pi. 101 * 102 * q0 the corresponding exponent of q[0]. Note that the 103 * exponent for q[i] would be q0-24*i. 104 * 105 * PIo2[] double precision array, obtained by cutting pi/2 106 * into 24 bits chunks. 107 * 108 * f[] ipio2[] in floating point 109 * 110 * iq[] integer array by breaking up q[] in 24-bits chunk. 111 * 112 * fq[] final product of x*(2/pi) in fq[0],..,fq[jk] 113 * 114 * ih integer. If >0 it indicates q[] is >= 0.5, hence 115 * it also indicates the *sign* of the result. 116 * 117 */ 118 119 120 /* 121 * Constants: 122 * The hexadecimal values are the intended ones for the following 123 * constants. The decimal values may be used, provided that the 124 * compiler will convert from decimal to binary accurately enough 125 * to produce the hexadecimal values shown. 126 */ 127 128 #include <float.h> 129 130 #include "math.h" 131 #include "math_private.h" 132 133 static const int init_jk[] = {3,4,4,6}; /* initial value for jk */ 134 135 /* 136 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 137 * 138 * integer array, contains the (24*i)-th to (24*i+23)-th 139 * bit of 2/pi after binary point. The corresponding 140 * floating value is 141 * 142 * ipio2[i] * 2^(-24(i+1)). 143 * 144 * NB: This table must have at least (e0-3)/24 + jk terms. 145 * For quad precision (e0 <= 16360, jk = 6), this is 686. 146 */ 147 static const int32_t ipio2[] = { 148 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 149 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 150 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 151 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 152 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 153 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 154 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 155 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 156 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 157 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 158 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 159 160 #if LDBL_MAX_EXP > 1024 161 #if LDBL_MAX_EXP > 16384 162 #error "ipio2 table needs to be expanded" 163 #endif 164 0x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6, 165 0xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2, 166 0xDE4F98, 0x327DBB, 0xC33D26, 0xEF6B1E, 0x5EF89F, 0x3A1F35, 167 0xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30, 168 0x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C, 169 0x467D86, 0x2D71E3, 0x9AC69B, 0x006233, 0x7CD2B4, 0x97A7B4, 170 0xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770, 171 0xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7, 172 0xCB2324, 0x778AD6, 0x23545A, 0xB91F00, 0x1B0AF1, 0xDFCE19, 173 0xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522, 174 0x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16, 175 0xDE3B58, 0x929BDE, 0x2822D2, 0xE88628, 0x4D58E2, 0x32CAC6, 176 0x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E, 177 0x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48, 178 0xD36710, 0xD8DDAA, 0x425FAE, 0xCE616A, 0xA4280A, 0xB499D3, 179 0xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF, 180 0xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55, 181 0x36D9CA, 0xD2A828, 0x8D61C2, 0x77C912, 0x142604, 0x9B4612, 182 0xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929, 183 0x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC, 184 0xC3E7B3, 0x28F8C7, 0x940593, 0x3E71C1, 0xB3092E, 0xF3450B, 185 0x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C, 186 0x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4, 187 0x9794E8, 0x84E6E2, 0x973199, 0x6BED88, 0x365F5F, 0x0EFDBB, 188 0xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC, 189 0x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C, 190 0x90AA47, 0x02E774, 0x24D6BD, 0xA67DF7, 0x72486E, 0xEF169F, 191 0xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5, 192 0x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437, 193 0x10D86D, 0x324832, 0x754C5B, 0xD4714E, 0x6E5445, 0xC1090B, 194 0x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA, 195 0x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD, 196 0x6AE290, 0x89D988, 0x50722C, 0xBEA404, 0x940777, 0x7030F3, 197 0x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3, 198 0xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717, 199 0x3BDF08, 0x2B3715, 0xA0805C, 0x93805A, 0x921110, 0xD8E80F, 200 0xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61, 201 0xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB, 202 0xAA140A, 0x2F2689, 0x768364, 0x333B09, 0x1A940E, 0xAA3A51, 203 0xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0, 204 0x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C, 205 0x5BC3D8, 0xC492F5, 0x4BADC6, 0xA5CA4E, 0xCD37A7, 0x36A9E6, 206 0x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC, 207 0xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED, 208 0x306529, 0xBF5657, 0x3AFF47, 0xB9F96A, 0xF3BE75, 0xDF9328, 209 0x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D, 210 0x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0, 211 0xA8654F, 0xA5C1D2, 0x0F3F0B, 0xCD785B, 0x76F923, 0x048B7B, 212 0x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4, 213 0xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3, 214 0xDA4886, 0xA05DF7, 0xF480C6, 0x2FF0AC, 0x9AECDD, 0xBC5C3F, 215 0x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD, 216 0x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B, 217 0x2A1216, 0x2DB7DC, 0xFDE5FA, 0xFEDB89, 0xFDBE89, 0x6C76E4, 218 0xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761, 219 0x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31, 220 0x48D784, 0x16DF30, 0x432DC7, 0x356125, 0xCE70C9, 0xB8CB30, 221 0xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262, 222 0x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E, 223 0xC4F133, 0x5F6E13, 0xE4305D, 0xA92E85, 0xC3B21D, 0x3632A1, 224 0xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C, 225 0x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4, 226 0xCBDA11, 0xD0BE7D, 0xC1DB9B, 0xBD17AB, 0x81A2CA, 0x5C6A08, 227 0x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196, 228 0xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9, 229 0x4F6A68, 0xA82A4A, 0x5AC44F, 0xBCF82D, 0x985AD7, 0x95C7F4, 230 0x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC, 231 0x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C, 232 0xD0C0B2, 0x485551, 0x0EFB1E, 0xC37295, 0x3B06A3, 0x3540C0, 233 0x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C, 234 0xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0, 235 0x3C3ABA, 0x461846, 0x5F7555, 0xF5BDD2, 0xC6926E, 0x5D2EAC, 236 0xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22, 237 0x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893, 238 0x745D7C, 0xB2AD6B, 0x9D6ECD, 0x7B723E, 0x6A11C6, 0xA9CFF7, 239 0xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5, 240 0x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F, 241 0xBEFDFD, 0xEF4556, 0x367ED9, 0x13D9EC, 0xB9BA8B, 0xFC97C4, 242 0x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF, 243 0x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B, 244 0x9C2A3E, 0xCC5F11, 0x4A0BFD, 0xFBF4E1, 0x6D3B8E, 0x2C86E2, 245 0x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138, 246 0xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E, 247 0xCC2254, 0xDC552A, 0xD6C6C0, 0x96190B, 0xB8701A, 0x649569, 248 0x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34, 249 0xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9, 250 0x9B5861, 0xBC57E1, 0xC68351, 0x103ED8, 0x4871DD, 0xDD1C2D, 251 0xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F, 252 0xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855, 253 0x382682, 0x9BE7CA, 0xA40D51, 0xB13399, 0x0ED7A9, 0x480569, 254 0xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B, 255 0x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE, 256 0x5FD45E, 0xA4677B, 0x7AACBA, 0xA2F655, 0x23882B, 0x55BA41, 257 0x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49, 258 0xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F, 259 0xAE5ADB, 0x86C547, 0x624385, 0x3B8621, 0x94792C, 0x876110, 260 0x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8, 261 0x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365, 262 0xB1933D, 0x0B7CBD, 0xDC51A4, 0x63DD27, 0xDDE169, 0x19949A, 263 0x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270, 264 0x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5, 265 0x4D7E6F, 0x5119A5, 0xABF9B5, 0xD6DF82, 0x61DD96, 0x023616, 266 0x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B, 267 0x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, 0x8071E0, 268 #endif 269 270 }; 271 272 static const double PIo2[] = { 273 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */ 274 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */ 275 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */ 276 3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */ 277 1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */ 278 1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */ 279 2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */ 280 2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */ 281 }; 282 283 static const double 284 zero = 0.0, 285 one = 1.0, 286 two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ 287 twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */ 288 289 int 290 __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec) 291 { 292 int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih; 293 double z,fw,f[20],fq[20],q[20]; 294 295 /* initialize jk*/ 296 jk = init_jk[prec]; 297 jp = jk; 298 299 /* determine jx,jv,q0, note that 3>q0 */ 300 jx = nx-1; 301 jv = (e0-3)/24; if(jv<0) jv=0; 302 q0 = e0-24*(jv+1); 303 304 /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ 305 j = jv-jx; m = jx+jk; 306 for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j]; 307 308 /* compute q[0],q[1],...q[jk] */ 309 for (i=0;i<=jk;i++) { 310 for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw; 311 } 312 313 jz = jk; 314 recompute: 315 /* distill q[] into iq[] reversingly */ 316 for(i=0,j=jz,z=q[jz];j>0;i++,j--) { 317 fw = (double)((int32_t)(twon24* z)); 318 iq[i] = (int32_t)(z-two24*fw); 319 z = q[j-1]+fw; 320 } 321 322 /* compute n */ 323 z = scalbn(z,q0); /* actual value of z */ 324 z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */ 325 n = (int32_t) z; 326 z -= (double)n; 327 ih = 0; 328 if(q0>0) { /* need iq[jz-1] to determine n */ 329 i = (iq[jz-1]>>(24-q0)); n += i; 330 iq[jz-1] -= i<<(24-q0); 331 ih = iq[jz-1]>>(23-q0); 332 } 333 else if(q0==0) ih = iq[jz-1]>>23; 334 else if(z>=0.5) ih=2; 335 336 if(ih>0) { /* q > 0.5 */ 337 n += 1; carry = 0; 338 for(i=0;i<jz ;i++) { /* compute 1-q */ 339 j = iq[i]; 340 if(carry==0) { 341 if(j!=0) { 342 carry = 1; iq[i] = 0x1000000- j; 343 } 344 } else iq[i] = 0xffffff - j; 345 } 346 if(q0>0) { /* rare case: chance is 1 in 12 */ 347 switch(q0) { 348 case 1: 349 iq[jz-1] &= 0x7fffff; break; 350 case 2: 351 iq[jz-1] &= 0x3fffff; break; 352 } 353 } 354 if(ih==2) { 355 z = one - z; 356 if(carry!=0) z -= scalbn(one,q0); 357 } 358 } 359 360 /* check if recomputation is needed */ 361 if(z==zero) { 362 j = 0; 363 for (i=jz-1;i>=jk;i--) j |= iq[i]; 364 if(j==0) { /* need recomputation */ 365 for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */ 366 367 for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */ 368 f[jx+i] = (double) ipio2[jv+i]; 369 for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; 370 q[i] = fw; 371 } 372 jz += k; 373 goto recompute; 374 } 375 } 376 377 /* chop off zero terms */ 378 if(z==0.0) { 379 jz -= 1; q0 -= 24; 380 while(iq[jz]==0) { jz--; q0-=24;} 381 } else { /* break z into 24-bit if necessary */ 382 z = scalbn(z,-q0); 383 if(z>=two24) { 384 fw = (double)((int32_t)(twon24*z)); 385 iq[jz] = (int32_t)(z-two24*fw); 386 jz += 1; q0 += 24; 387 iq[jz] = (int32_t) fw; 388 } else iq[jz] = (int32_t) z ; 389 } 390 391 /* convert integer "bit" chunk to floating-point value */ 392 fw = scalbn(one,q0); 393 for(i=jz;i>=0;i--) { 394 q[i] = fw*(double)iq[i]; fw*=twon24; 395 } 396 397 /* compute PIo2[0,...,jp]*q[jz,...,0] */ 398 for(i=jz;i>=0;i--) { 399 for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k]; 400 fq[jz-i] = fw; 401 } 402 403 /* compress fq[] into y[] */ 404 switch(prec) { 405 case 0: 406 fw = 0.0; 407 for (i=jz;i>=0;i--) fw += fq[i]; 408 y[0] = (ih==0)? fw: -fw; 409 break; 410 case 1: 411 case 2: 412 fw = 0.0; 413 for (i=jz;i>=0;i--) fw += fq[i]; 414 STRICT_ASSIGN(double,fw,fw); 415 y[0] = (ih==0)? fw: -fw; 416 fw = fq[0]-fw; 417 for (i=1;i<=jz;i++) fw += fq[i]; 418 y[1] = (ih==0)? fw: -fw; 419 break; 420 case 3: /* painful */ 421 for (i=jz;i>0;i--) { 422 fw = fq[i-1]+fq[i]; 423 fq[i] += fq[i-1]-fw; 424 fq[i-1] = fw; 425 } 426 for (i=jz;i>1;i--) { 427 fw = fq[i-1]+fq[i]; 428 fq[i] += fq[i-1]-fw; 429 fq[i-1] = fw; 430 } 431 for (fw=0.0,i=jz;i>=2;i--) fw += fq[i]; 432 if(ih==0) { 433 y[0] = fq[0]; y[1] = fq[1]; y[2] = fw; 434 } else { 435 y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw; 436 } 437 } 438 return n&7; 439 } 440