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