GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: lib/libm/src/k_rem_pio2.c Lines: 70 92 76.1 %
Date: 2017-11-13 Branches: 53 78 67.9 %

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

7840
	for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
302
303
    /* compute q[0],q[1],...q[jk] */
304
3360
	for (i=0;i<=jk;i++) {
305
12600
	    for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
306
	}
307
308
280
	jz = jk;
309
recompute:
310
    /* distill q[] into iq[] reversingly */
311
3472
	for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
312
1400
	    fw    =  (double)((int32_t)(twon24* z));
313
1400
	    iq[i] =  (int32_t)(z-two24*fw);
314
1400
	    z     =  q[j-1]+fw;
315
	}
316
317
    /* compute n */
318
336
	z  = scalbn(z,q0);		/* actual value of z */
319
336
	z -= 8.0*floor(z*0.125);		/* trim off integer >= 8 */
320
336
	n  = (int32_t) z;
321
336
	z -= (double)n;
322
	ih = 0;
323
336
	if(q0>0) {	/* need iq[jz-1] to determine n */
324
	    i  = (iq[jz-1]>>(24-q0)); n += i;
325
	    iq[jz-1] -= i<<(24-q0);
326
	    ih = iq[jz-1]>>(23-q0);
327
	}
328
336
	else if(q0==0) ih = iq[jz-1]>>23;
329
560
	else if(z>=0.5) ih=2;
330
331
336
	if(ih>0) {	/* q > 0.5 */
332
224
	    n += 1; carry = 0;
333
2352
	    for(i=0;i<jz ;i++) {	/* compute 1-q */
334
952
		j = iq[i];
335
952
		if(carry==0) {
336
224
		    if(j!=0) {
337
224
			carry = 1; iq[i] = 0x1000000- j;
338
224
		    }
339
728
		} else  iq[i] = 0xffffff - j;
340
952
	    }
341
224
	    if(q0>0) {		/* rare case: chance is 1 in 12 */
342
	        switch(q0) {
343
	        case 1:
344
	    	   iq[jz-1] &= 0x7fffff; break;
345
	    	case 2:
346
	    	   iq[jz-1] &= 0x3fffff; break;
347
	        }
348
	    }
349
224
	    if(ih==2) {
350
224
		z = one - z;
351
448
		if(carry!=0) z -= scalbn(one,q0);
352
	    }
353
	}
354
355
    /* check if recomputation is needed */
356
336
	if(z==zero) {
357
	    j = 0;
358
336
	    for (i=jz-1;i>=jk;i--) j |= iq[i];
359
112
	    if(j==0) { /* need recomputation */
360
112
		for(k=1;iq[jk-k]==0;k++);   /* k = no. of terms needed */
361
362
224
		for(i=jz+1;i<=jz+k;i++) {   /* add q[jz+1] to q[jz+k] */
363
56
		    f[jx+i] = (double) ipio2[jv+i];
364
448
		    for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
365
56
		    q[i] = fw;
366
		}
367
		jz += k;
368
56
		goto recompute;
369
	    }
370
	}
371
372
    /* chop off zero terms */
373
280
	if(z==0.0) {
374
56
	    jz -= 1; q0 -= 24;
375
112
	    while(iq[jz]==0) { jz--; q0-=24;}
376
	} else { /* break z into 24-bit if necessary */
377
224
	    z = scalbn(z,-q0);
378
224
	    if(z>=two24) {
379
224
		fw = (double)((int32_t)(twon24*z));
380
224
		iq[jz] = (int32_t)(z-two24*fw);
381
224
		jz += 1; q0 += 24;
382
224
		iq[jz] = (int32_t) fw;
383
224
	    } else iq[jz] = (int32_t) z ;
384
	}
385
386
    /* convert integer "bit" chunk to floating-point value */
387
280
	fw = scalbn(one,q0);
388
3808
	for(i=jz;i>=0;i--) {
389
1624
	    q[i] = fw*(double)iq[i]; fw*=twon24;
390
	}
391
392
    /* compute PIo2[0,...,jp]*q[jz,...,0] */
393
3808
	for(i=jz;i>=0;i--) {
394

20328
	    for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
395
1624
	    fq[jz-i] = fw;
396
	}
397
398
    /* compress fq[] into y[] */
399

560
	switch(prec) {
400
	    case 0:
401
		fw = 0.0;
402
		for (i=jz;i>=0;i--) fw += fq[i];
403
		y[0] = (ih==0)? fw: -fw;
404
		break;
405
	    case 1:
406
	    case 2:
407
		fw = 0.0;
408
3808
		for (i=jz;i>=0;i--) fw += fq[i];
409
		STRICT_ASSIGN(double,fw,fw);
410
280
		y[0] = (ih==0)? fw: -fw;
411
280
		fw = fq[0]-fw;
412
3248
		for (i=1;i<=jz;i++) fw += fq[i];
413
280
		y[1] = (ih==0)? fw: -fw;
414
280
		break;
415
	    case 3:	/* painful */
416
		for (i=jz;i>0;i--) {
417
		    fw      = fq[i-1]+fq[i];
418
		    fq[i]  += fq[i-1]-fw;
419
		    fq[i-1] = fw;
420
		}
421
		for (i=jz;i>1;i--) {
422
		    fw      = fq[i-1]+fq[i];
423
		    fq[i]  += fq[i-1]-fw;
424
		    fq[i-1] = fw;
425
		}
426
		for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
427
		if(ih==0) {
428
		    y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
429
		} else {
430
		    y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
431
		}
432
	}
433
840
	return n&7;
434
280
}