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