xref: /freebsd/lib/msun/src/k_rem_pio2.c (revision 3fb8f1272b50cb87cb624b321f7b81e76627c437)
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