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.15 2023/08/08 06:31:17 mrg 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
__kernel_rem_pio2(double * x,double * y,int e0,int nx,int prec)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 /* if nx < 2, fq[0] may be accessed uninitialised */
306 if (nx < 2) {
307 fq[0] = 0;
308 }
309
310 /* initialize jk*/
311 jk = init_jk[prec];
312 jp = jk;
313
314 /* determine jx,jv,q0, note that 3>q0 */
315 jx = nx-1;
316 jv = (e0-3)/24; if(jv<0) jv=0;
317 q0 = e0-24*(jv+1);
318
319 /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
320 j = jv-jx; m = jx+jk;
321 for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
322
323 /* compute q[0],q[1],...q[jk] */
324 for (i=0;i<=jk;i++) {
325 for(j=0,fw=0.0;j<=jx;j++)
326 fw += x[j]*f[jx+i-j];
327 q[i] = fw;
328 }
329
330 jz = jk;
331 recompute:
332 /* distill q[] into iq[] reversingly */
333 for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
334 fw = (double)((int32_t)(twon24* z));
335 iq[i] = (int32_t)(z-two24*fw);
336 z = q[j-1]+fw;
337 }
338
339 /* compute n */
340 z = scalbn(z,q0); /* actual value of z */
341 z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */
342 n = (int32_t) z;
343 z -= (double)n;
344 ih = 0;
345 if(q0>0) { /* need iq[jz-1] to determine n */
346 i = (iq[jz-1]>>(24-q0)); n += i;
347 iq[jz-1] -= i<<(24-q0);
348 ih = iq[jz-1]>>(23-q0);
349 }
350 else if(q0==0) ih = iq[jz-1]>>23;
351 else if(z>=0.5) ih=2;
352
353 if(ih>0) { /* q > 0.5 */
354 n += 1; carry = 0;
355 for(i=0;i<jz ;i++) { /* compute 1-q */
356 j = iq[i];
357 if(carry==0) {
358 if(j!=0) {
359 carry = 1; iq[i] = 0x1000000- j;
360 }
361 } else iq[i] = 0xffffff - j;
362 }
363 if(q0>0) { /* rare case: chance is 1 in 12 */
364 switch(q0) {
365 case 1:
366 iq[jz-1] &= 0x7fffff; break;
367 case 2:
368 iq[jz-1] &= 0x3fffff; break;
369 }
370 }
371 if(ih==2) {
372 z = one - z;
373 if(carry!=0) z -= scalbn(one,q0);
374 }
375 }
376
377 /* check if recomputation is needed */
378 if(z==zero) {
379 j = 0;
380 for (i=jz-1;i>=jk;i--) j |= iq[i];
381 if(j==0) { /* need recomputation */
382 for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */
383
384 for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */
385 f[jx+i] = (double) ipio2[jv+i];
386 for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
387 q[i] = fw;
388 }
389 jz += k;
390 goto recompute;
391 }
392 }
393
394 /* chop off zero terms */
395 if(z==0.0) {
396 jz -= 1; q0 -= 24;
397 while(iq[jz]==0) { jz--; q0-=24;}
398 } else { /* break z into 24-bit if necessary */
399 z = scalbn(z,-q0);
400 if(z>=two24) {
401 fw = (double)((int32_t)(twon24*z));
402 iq[jz] = (int32_t)(z-two24*fw);
403 jz += 1; q0 += 24;
404 iq[jz] = (int32_t) fw;
405 } else iq[jz] = (int32_t) z ;
406 }
407
408 /* convert integer "bit" chunk to floating-point value */
409 fw = scalbn(one,q0);
410 for(i=jz;i>=0;i--) {
411 q[i] = fw*(double)iq[i]; fw*=twon24;
412 }
413
414 /* compute PIo2[0,...,jp]*q[jz,...,0] */
415 for(i=jz;i>=0;i--) {
416 for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
417 fq[jz-i] = fw;
418 }
419
420 /* compress fq[] into y[] */
421 switch(prec) {
422 case 0:
423 fw = 0.0;
424 for (i=jz;i>=0;i--) fw += fq[i];
425 y[0] = (ih==0)? fw: -fw;
426 break;
427 case 1:
428 case 2:
429 fw = 0.0;
430 for (i=jz;i>=0;i--) fw += fq[i];
431 STRICT_ASSIGN(double,fw,fw);
432 y[0] = (ih==0)? fw: -fw;
433 fw = fq[0]-fw;
434 for (i=1;i<=jz;i++) fw += fq[i];
435 y[1] = (ih==0)? fw: -fw;
436 break;
437 case 3: /* painful */
438 for (i=jz;i>0;i--) {
439 fw = fq[i-1]+fq[i];
440 fq[i] += fq[i-1]-fw;
441 fq[i-1] = fw;
442 }
443 for (i=jz;i>1;i--) {
444 fw = fq[i-1]+fq[i];
445 fq[i] += fq[i-1]-fw;
446 fq[i-1] = fw;
447 }
448 for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
449 if(ih==0) {
450 y[0] = fq[0]; y[1] = fq[1]; y[2] = fw;
451 } else {
452 y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
453 }
454 }
455 return n&7;
456 }
457