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