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 cancellation 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