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