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
143static 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 */
157static const int32_t ipio2[] = {
1580xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
1590x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
1600x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
1610xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
1620x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
1630x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
1640x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
1650xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
1660x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
1670x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
1680x4D7327, 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
1740x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6,
1750xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2,
1760xDE4F98, 0x327DBB, 0xC33D26, 0xEF6B1E, 0x5EF89F, 0x3A1F35,
1770xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30,
1780x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C,
1790x467D86, 0x2D71E3, 0x9AC69B, 0x006233, 0x7CD2B4, 0x97A7B4,
1800xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770,
1810xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7,
1820xCB2324, 0x778AD6, 0x23545A, 0xB91F00, 0x1B0AF1, 0xDFCE19,
1830xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522,
1840x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16,
1850xDE3B58, 0x929BDE, 0x2822D2, 0xE88628, 0x4D58E2, 0x32CAC6,
1860x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E,
1870x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48,
1880xD36710, 0xD8DDAA, 0x425FAE, 0xCE616A, 0xA4280A, 0xB499D3,
1890xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF,
1900xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55,
1910x36D9CA, 0xD2A828, 0x8D61C2, 0x77C912, 0x142604, 0x9B4612,
1920xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929,
1930x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC,
1940xC3E7B3, 0x28F8C7, 0x940593, 0x3E71C1, 0xB3092E, 0xF3450B,
1950x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C,
1960x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4,
1970x9794E8, 0x84E6E2, 0x973199, 0x6BED88, 0x365F5F, 0x0EFDBB,
1980xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC,
1990x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C,
2000x90AA47, 0x02E774, 0x24D6BD, 0xA67DF7, 0x72486E, 0xEF169F,
2010xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5,
2020x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437,
2030x10D86D, 0x324832, 0x754C5B, 0xD4714E, 0x6E5445, 0xC1090B,
2040x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA,
2050x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD,
2060x6AE290, 0x89D988, 0x50722C, 0xBEA404, 0x940777, 0x7030F3,
2070x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3,
2080xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717,
2090x3BDF08, 0x2B3715, 0xA0805C, 0x93805A, 0x921110, 0xD8E80F,
2100xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61,
2110xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB,
2120xAA140A, 0x2F2689, 0x768364, 0x333B09, 0x1A940E, 0xAA3A51,
2130xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0,
2140x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C,
2150x5BC3D8, 0xC492F5, 0x4BADC6, 0xA5CA4E, 0xCD37A7, 0x36A9E6,
2160x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC,
2170xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED,
2180x306529, 0xBF5657, 0x3AFF47, 0xB9F96A, 0xF3BE75, 0xDF9328,
2190x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D,
2200x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0,
2210xA8654F, 0xA5C1D2, 0x0F3F0B, 0xCD785B, 0x76F923, 0x048B7B,
2220x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4,
2230xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3,
2240xDA4886, 0xA05DF7, 0xF480C6, 0x2FF0AC, 0x9AECDD, 0xBC5C3F,
2250x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD,
2260x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B,
2270x2A1216, 0x2DB7DC, 0xFDE5FA, 0xFEDB89, 0xFDBE89, 0x6C76E4,
2280xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761,
2290x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31,
2300x48D784, 0x16DF30, 0x432DC7, 0x356125, 0xCE70C9, 0xB8CB30,
2310xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262,
2320x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E,
2330xC4F133, 0x5F6E13, 0xE4305D, 0xA92E85, 0xC3B21D, 0x3632A1,
2340xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C,
2350x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4,
2360xCBDA11, 0xD0BE7D, 0xC1DB9B, 0xBD17AB, 0x81A2CA, 0x5C6A08,
2370x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196,
2380xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9,
2390x4F6A68, 0xA82A4A, 0x5AC44F, 0xBCF82D, 0x985AD7, 0x95C7F4,
2400x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC,
2410x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C,
2420xD0C0B2, 0x485551, 0x0EFB1E, 0xC37295, 0x3B06A3, 0x3540C0,
2430x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C,
2440xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0,
2450x3C3ABA, 0x461846, 0x5F7555, 0xF5BDD2, 0xC6926E, 0x5D2EAC,
2460xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22,
2470x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893,
2480x745D7C, 0xB2AD6B, 0x9D6ECD, 0x7B723E, 0x6A11C6, 0xA9CFF7,
2490xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5,
2500x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F,
2510xBEFDFD, 0xEF4556, 0x367ED9, 0x13D9EC, 0xB9BA8B, 0xFC97C4,
2520x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF,
2530x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B,
2540x9C2A3E, 0xCC5F11, 0x4A0BFD, 0xFBF4E1, 0x6D3B8E, 0x2C86E2,
2550x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138,
2560xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E,
2570xCC2254, 0xDC552A, 0xD6C6C0, 0x96190B, 0xB8701A, 0x649569,
2580x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34,
2590xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9,
2600x9B5861, 0xBC57E1, 0xC68351, 0x103ED8, 0x4871DD, 0xDD1C2D,
2610xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F,
2620xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855,
2630x382682, 0x9BE7CA, 0xA40D51, 0xB13399, 0x0ED7A9, 0x480569,
2640xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B,
2650x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE,
2660x5FD45E, 0xA4677B, 0x7AACBA, 0xA2F655, 0x23882B, 0x55BA41,
2670x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49,
2680xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F,
2690xAE5ADB, 0x86C547, 0x624385, 0x3B8621, 0x94792C, 0x876110,
2700x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8,
2710x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365,
2720xB1933D, 0x0B7CBD, 0xDC51A4, 0x63DD27, 0xDDE169, 0x19949A,
2730x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270,
2740x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5,
2750x4D7E6F, 0x5119A5, 0xABF9B5, 0xD6DF82, 0x61DD96, 0x023616,
2760x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B,
2770x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, 0x8071E0,
278#endif
279
280};
281
282static 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
293static const double
294zero   = 0.0,
295one    = 1.0,
296two24   =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
297twon24  =  5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
298
299int
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;
331recompute:
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