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