k_rem_pio2.c revision 176356
1141296Sdas 2141296Sdas/* @(#)k_rem_pio2.c 1.3 95/01/18 */ 32116Sjkh/* 42116Sjkh * ==================================================== 52116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 62116Sjkh * 7141296Sdas * Developed at SunSoft, a Sun Microsystems, Inc. business. 82116Sjkh * Permission to use, copy, modify, and distribute this 9141296Sdas * software is freely granted, provided that this notice 102116Sjkh * is preserved. 112116Sjkh * ==================================================== 122116Sjkh */ 132116Sjkh 14175499Sbde#include <sys/cdefs.h> 15175499Sbde__FBSDID("$FreeBSD: head/lib/msun/src/k_rem_pio2.c 176356 2008-02-17 07:31:59Z das $"); 162116Sjkh 172116Sjkh/* 18176356Sdas * __kernel_rem_pio2(x,y,e0,nx,prec) 19176356Sdas * double x[],y[]; int e0,nx,prec; 20141296Sdas * 21141296Sdas * __kernel_rem_pio2 return the last three digits of N with 222116Sjkh * y = x - N*pi/2 232116Sjkh * so that |y| < pi/2. 242116Sjkh * 25141296Sdas * The method is to compute the integer (mod 8) and fraction parts of 262116Sjkh * (2/pi)*x without doing the full multiplication. In general we 272116Sjkh * skip the part of the product that are known to be a huge integer ( 282116Sjkh * more accurately, = 0 mod 8 ). Thus the number of operations are 292116Sjkh * independent of the exponent of the input. 302116Sjkh * 312116Sjkh * (2/pi) is represented by an array of 24-bit integers in ipio2[]. 322116Sjkh * 332116Sjkh * Input parameters: 34141296Sdas * x[] The input value (must be positive) is broken into nx 352116Sjkh * pieces of 24-bit integers in double precision format. 36141296Sdas * x[i] will be the i-th 24 bit of x. The scaled exponent 37141296Sdas * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 382116Sjkh * match x's up to 24 bits. 392116Sjkh * 402116Sjkh * Example of breaking a double positive z into x[0]+x[1]+x[2]: 412116Sjkh * e0 = ilogb(z)-23 422116Sjkh * z = scalbn(z,-e0) 432116Sjkh * for i = 0,1,2 442116Sjkh * x[i] = floor(z) 452116Sjkh * z = (z-x[i])*2**24 462116Sjkh * 472116Sjkh * 482116Sjkh * y[] ouput result in an array of double precision numbers. 492116Sjkh * The dimension of y[] is: 502116Sjkh * 24-bit precision 1 512116Sjkh * 53-bit precision 2 522116Sjkh * 64-bit precision 2 532116Sjkh * 113-bit precision 3 542116Sjkh * The actual value is the sum of them. Thus for 113-bit 552116Sjkh * precison, one may have to do something like: 562116Sjkh * 572116Sjkh * long double t,w,r_head, r_tail; 582116Sjkh * t = (long double)y[2] + (long double)y[1]; 592116Sjkh * w = (long double)y[0]; 602116Sjkh * r_head = t+w; 612116Sjkh * r_tail = w - (r_head - t); 622116Sjkh * 63176356Sdas * e0 The exponent of x[0]. Must be <= 16360 or you need to 64176356Sdas * expand the ipio2 table. 652116Sjkh * 662116Sjkh * nx dimension of x[] 672116Sjkh * 682116Sjkh * prec an integer indicating the precision: 692116Sjkh * 0 24 bits (single) 702116Sjkh * 1 53 bits (double) 712116Sjkh * 2 64 bits (extended) 722116Sjkh * 3 113 bits (quad) 732116Sjkh * 742116Sjkh * External function: 752116Sjkh * double scalbn(), floor(); 762116Sjkh * 772116Sjkh * 782116Sjkh * Here is the description of some local variables: 792116Sjkh * 802116Sjkh * jk jk+1 is the initial number of terms of ipio2[] needed 812116Sjkh * in the computation. The recommended value is 2,3,4, 822116Sjkh * 6 for single, double, extended,and quad. 832116Sjkh * 84141296Sdas * jz local integer variable indicating the number of 85141296Sdas * terms of ipio2[] used. 862116Sjkh * 872116Sjkh * jx nx - 1 882116Sjkh * 892116Sjkh * jv index for pointing to the suitable ipio2[] for the 902116Sjkh * computation. In general, we want 912116Sjkh * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8 922116Sjkh * is an integer. Thus 932116Sjkh * e0-3-24*jv >= 0 or (e0-3)/24 >= jv 942116Sjkh * Hence jv = max(0,(e0-3)/24). 952116Sjkh * 962116Sjkh * jp jp+1 is the number of terms in PIo2[] needed, jp = jk. 972116Sjkh * 982116Sjkh * q[] double array with integral value, representing the 992116Sjkh * 24-bits chunk of the product of x and 2/pi. 1002116Sjkh * 1012116Sjkh * q0 the corresponding exponent of q[0]. Note that the 1022116Sjkh * exponent for q[i] would be q0-24*i. 1032116Sjkh * 1042116Sjkh * PIo2[] double precision array, obtained by cutting pi/2 105141296Sdas * into 24 bits chunks. 1062116Sjkh * 107141296Sdas * f[] ipio2[] in floating point 1082116Sjkh * 1092116Sjkh * iq[] integer array by breaking up q[] in 24-bits chunk. 1102116Sjkh * 1112116Sjkh * fq[] final product of x*(2/pi) in fq[0],..,fq[jk] 1122116Sjkh * 1132116Sjkh * ih integer. If >0 it indicates q[] is >= 0.5, hence 1142116Sjkh * it also indicates the *sign* of the result. 1152116Sjkh * 1162116Sjkh */ 1172116Sjkh 1182116Sjkh 1192116Sjkh/* 1202116Sjkh * Constants: 121141296Sdas * The hexadecimal values are the intended ones for the following 122141296Sdas * constants. The decimal values may be used, provided that the 123141296Sdas * compiler will convert from decimal to binary accurately enough 1242116Sjkh * to produce the hexadecimal values shown. 1252116Sjkh */ 1262116Sjkh 127175499Sbde#include <float.h> 128175499Sbde 1292116Sjkh#include "math.h" 1302116Sjkh#include "math_private.h" 1312116Sjkh 1322116Sjkhstatic const int init_jk[] = {2,3,4,6}; /* initial value for jk */ 1332116Sjkh 134176356Sdas/* 135176356Sdas * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 136176356Sdas * 137176356Sdas * integer array, contains the (24*i)-th to (24*i+23)-th 138176356Sdas * bit of 2/pi after binary point. The corresponding 139176356Sdas * floating value is 140176356Sdas * 141176356Sdas * ipio2[i] * 2^(-24(i+1)). 142176356Sdas * 143176356Sdas * NB: This table must have at least (e0-3)/24 + jk terms. 144176356Sdas * For quad precision (e0 <= 16360, jk = 6), this is 686. 145176356Sdas */ 146176356Sdasstatic const int32_t ipio2[] = { 147176356Sdas0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 148176356Sdas0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 149176356Sdas0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 150176356Sdas0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 151176356Sdas0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 152176356Sdas0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 153176356Sdas0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 154176356Sdas0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 155176356Sdas0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 156176356Sdas0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 157176356Sdas0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 158176356Sdas 159176356Sdas#if LDBL_MAX_EXP > 1024 160176356Sdas#if LDBL_MAX_EXP > 16384 161176356Sdas#error "ipio2 table needs to be expanded" 162176356Sdas#endif 163176356Sdas0x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6, 164176356Sdas0xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2, 165176356Sdas0xDE4F98, 0x327DBB, 0xC33D26, 0xEF6B1E, 0x5EF89F, 0x3A1F35, 166176356Sdas0xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30, 167176356Sdas0x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C, 168176356Sdas0x467D86, 0x2D71E3, 0x9AC69B, 0x006233, 0x7CD2B4, 0x97A7B4, 169176356Sdas0xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770, 170176356Sdas0xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7, 171176356Sdas0xCB2324, 0x778AD6, 0x23545A, 0xB91F00, 0x1B0AF1, 0xDFCE19, 172176356Sdas0xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522, 173176356Sdas0x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16, 174176356Sdas0xDE3B58, 0x929BDE, 0x2822D2, 0xE88628, 0x4D58E2, 0x32CAC6, 175176356Sdas0x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E, 176176356Sdas0x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48, 177176356Sdas0xD36710, 0xD8DDAA, 0x425FAE, 0xCE616A, 0xA4280A, 0xB499D3, 178176356Sdas0xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF, 179176356Sdas0xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55, 180176356Sdas0x36D9CA, 0xD2A828, 0x8D61C2, 0x77C912, 0x142604, 0x9B4612, 181176356Sdas0xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929, 182176356Sdas0x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC, 183176356Sdas0xC3E7B3, 0x28F8C7, 0x940593, 0x3E71C1, 0xB3092E, 0xF3450B, 184176356Sdas0x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C, 185176356Sdas0x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4, 186176356Sdas0x9794E8, 0x84E6E2, 0x973199, 0x6BED88, 0x365F5F, 0x0EFDBB, 187176356Sdas0xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC, 188176356Sdas0x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C, 189176356Sdas0x90AA47, 0x02E774, 0x24D6BD, 0xA67DF7, 0x72486E, 0xEF169F, 190176356Sdas0xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5, 191176356Sdas0x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437, 192176356Sdas0x10D86D, 0x324832, 0x754C5B, 0xD4714E, 0x6E5445, 0xC1090B, 193176356Sdas0x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA, 194176356Sdas0x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD, 195176356Sdas0x6AE290, 0x89D988, 0x50722C, 0xBEA404, 0x940777, 0x7030F3, 196176356Sdas0x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3, 197176356Sdas0xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717, 198176356Sdas0x3BDF08, 0x2B3715, 0xA0805C, 0x93805A, 0x921110, 0xD8E80F, 199176356Sdas0xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61, 200176356Sdas0xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB, 201176356Sdas0xAA140A, 0x2F2689, 0x768364, 0x333B09, 0x1A940E, 0xAA3A51, 202176356Sdas0xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0, 203176356Sdas0x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C, 204176356Sdas0x5BC3D8, 0xC492F5, 0x4BADC6, 0xA5CA4E, 0xCD37A7, 0x36A9E6, 205176356Sdas0x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC, 206176356Sdas0xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED, 207176356Sdas0x306529, 0xBF5657, 0x3AFF47, 0xB9F96A, 0xF3BE75, 0xDF9328, 208176356Sdas0x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D, 209176356Sdas0x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0, 210176356Sdas0xA8654F, 0xA5C1D2, 0x0F3F0B, 0xCD785B, 0x76F923, 0x048B7B, 211176356Sdas0x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4, 212176356Sdas0xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3, 213176356Sdas0xDA4886, 0xA05DF7, 0xF480C6, 0x2FF0AC, 0x9AECDD, 0xBC5C3F, 214176356Sdas0x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD, 215176356Sdas0x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B, 216176356Sdas0x2A1216, 0x2DB7DC, 0xFDE5FA, 0xFEDB89, 0xFDBE89, 0x6C76E4, 217176356Sdas0xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761, 218176356Sdas0x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31, 219176356Sdas0x48D784, 0x16DF30, 0x432DC7, 0x356125, 0xCE70C9, 0xB8CB30, 220176356Sdas0xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262, 221176356Sdas0x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E, 222176356Sdas0xC4F133, 0x5F6E13, 0xE4305D, 0xA92E85, 0xC3B21D, 0x3632A1, 223176356Sdas0xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C, 224176356Sdas0x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4, 225176356Sdas0xCBDA11, 0xD0BE7D, 0xC1DB9B, 0xBD17AB, 0x81A2CA, 0x5C6A08, 226176356Sdas0x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196, 227176356Sdas0xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9, 228176356Sdas0x4F6A68, 0xA82A4A, 0x5AC44F, 0xBCF82D, 0x985AD7, 0x95C7F4, 229176356Sdas0x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC, 230176356Sdas0x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C, 231176356Sdas0xD0C0B2, 0x485551, 0x0EFB1E, 0xC37295, 0x3B06A3, 0x3540C0, 232176356Sdas0x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C, 233176356Sdas0xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0, 234176356Sdas0x3C3ABA, 0x461846, 0x5F7555, 0xF5BDD2, 0xC6926E, 0x5D2EAC, 235176356Sdas0xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22, 236176356Sdas0x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893, 237176356Sdas0x745D7C, 0xB2AD6B, 0x9D6ECD, 0x7B723E, 0x6A11C6, 0xA9CFF7, 238176356Sdas0xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5, 239176356Sdas0x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F, 240176356Sdas0xBEFDFD, 0xEF4556, 0x367ED9, 0x13D9EC, 0xB9BA8B, 0xFC97C4, 241176356Sdas0x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF, 242176356Sdas0x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B, 243176356Sdas0x9C2A3E, 0xCC5F11, 0x4A0BFD, 0xFBF4E1, 0x6D3B8E, 0x2C86E2, 244176356Sdas0x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138, 245176356Sdas0xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E, 246176356Sdas0xCC2254, 0xDC552A, 0xD6C6C0, 0x96190B, 0xB8701A, 0x649569, 247176356Sdas0x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34, 248176356Sdas0xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9, 249176356Sdas0x9B5861, 0xBC57E1, 0xC68351, 0x103ED8, 0x4871DD, 0xDD1C2D, 250176356Sdas0xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F, 251176356Sdas0xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855, 252176356Sdas0x382682, 0x9BE7CA, 0xA40D51, 0xB13399, 0x0ED7A9, 0x480569, 253176356Sdas0xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B, 254176356Sdas0x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE, 255176356Sdas0x5FD45E, 0xA4677B, 0x7AACBA, 0xA2F655, 0x23882B, 0x55BA41, 256176356Sdas0x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49, 257176356Sdas0xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F, 258176356Sdas0xAE5ADB, 0x86C547, 0x624385, 0x3B8621, 0x94792C, 0x876110, 259176356Sdas0x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8, 260176356Sdas0x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365, 261176356Sdas0xB1933D, 0x0B7CBD, 0xDC51A4, 0x63DD27, 0xDDE169, 0x19949A, 262176356Sdas0x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270, 263176356Sdas0x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5, 264176356Sdas0x4D7E6F, 0x5119A5, 0xABF9B5, 0xD6DF82, 0x61DD96, 0x023616, 265176356Sdas0x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B, 266176356Sdas0x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, 0x8071E0, 267176356Sdas#endif 268176356Sdas 269176356Sdas}; 270176356Sdas 2712116Sjkhstatic const double PIo2[] = { 2722116Sjkh 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */ 2732116Sjkh 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */ 2742116Sjkh 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */ 2752116Sjkh 3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */ 2762116Sjkh 1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */ 2772116Sjkh 1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */ 2782116Sjkh 2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */ 2792116Sjkh 2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */ 2802116Sjkh}; 2812116Sjkh 282141296Sdasstatic const double 2832116Sjkhzero = 0.0, 2842116Sjkhone = 1.0, 2852116Sjkhtwo24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ 2862116Sjkhtwon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */ 2872116Sjkh 288176356Sdasint 289176356Sdas__kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec) 2902116Sjkh{ 2912116Sjkh int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih; 2922116Sjkh double z,fw,f[20],fq[20],q[20]; 2932116Sjkh 2942116Sjkh /* initialize jk*/ 2952116Sjkh jk = init_jk[prec]; 2962116Sjkh jp = jk; 2972116Sjkh 2982116Sjkh /* determine jx,jv,q0, note that 3>q0 */ 2992116Sjkh jx = nx-1; 3002116Sjkh jv = (e0-3)/24; if(jv<0) jv=0; 3012116Sjkh q0 = e0-24*(jv+1); 3022116Sjkh 3032116Sjkh /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ 3042116Sjkh j = jv-jx; m = jx+jk; 3052116Sjkh for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j]; 3062116Sjkh 3072116Sjkh /* compute q[0],q[1],...q[jk] */ 3082116Sjkh for (i=0;i<=jk;i++) { 3092116Sjkh for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw; 3102116Sjkh } 3112116Sjkh 3122116Sjkh jz = jk; 3132116Sjkhrecompute: 3142116Sjkh /* distill q[] into iq[] reversingly */ 3152116Sjkh for(i=0,j=jz,z=q[jz];j>0;i++,j--) { 3162116Sjkh fw = (double)((int32_t)(twon24* z)); 3172116Sjkh iq[i] = (int32_t)(z-two24*fw); 3182116Sjkh z = q[j-1]+fw; 3192116Sjkh } 3202116Sjkh 3212116Sjkh /* compute n */ 3222116Sjkh z = scalbn(z,q0); /* actual value of z */ 3232116Sjkh z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */ 3242116Sjkh n = (int32_t) z; 3252116Sjkh z -= (double)n; 3262116Sjkh ih = 0; 3272116Sjkh if(q0>0) { /* need iq[jz-1] to determine n */ 3282116Sjkh i = (iq[jz-1]>>(24-q0)); n += i; 3292116Sjkh iq[jz-1] -= i<<(24-q0); 3302116Sjkh ih = iq[jz-1]>>(23-q0); 331141296Sdas } 3322116Sjkh else if(q0==0) ih = iq[jz-1]>>23; 3332116Sjkh else if(z>=0.5) ih=2; 3342116Sjkh 3352116Sjkh if(ih>0) { /* q > 0.5 */ 3362116Sjkh n += 1; carry = 0; 3372116Sjkh for(i=0;i<jz ;i++) { /* compute 1-q */ 3382116Sjkh j = iq[i]; 3392116Sjkh if(carry==0) { 3402116Sjkh if(j!=0) { 3412116Sjkh carry = 1; iq[i] = 0x1000000- j; 3422116Sjkh } 3432116Sjkh } else iq[i] = 0xffffff - j; 3442116Sjkh } 3452116Sjkh if(q0>0) { /* rare case: chance is 1 in 12 */ 3462116Sjkh switch(q0) { 3472116Sjkh case 1: 3482116Sjkh iq[jz-1] &= 0x7fffff; break; 3492116Sjkh case 2: 3502116Sjkh iq[jz-1] &= 0x3fffff; break; 3512116Sjkh } 3522116Sjkh } 3532116Sjkh if(ih==2) { 3542116Sjkh z = one - z; 3552116Sjkh if(carry!=0) z -= scalbn(one,q0); 3562116Sjkh } 3572116Sjkh } 3582116Sjkh 3592116Sjkh /* check if recomputation is needed */ 3602116Sjkh if(z==zero) { 3612116Sjkh j = 0; 3622116Sjkh for (i=jz-1;i>=jk;i--) j |= iq[i]; 3632116Sjkh if(j==0) { /* need recomputation */ 3642116Sjkh for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */ 3652116Sjkh 3662116Sjkh for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */ 3672116Sjkh f[jx+i] = (double) ipio2[jv+i]; 3682116Sjkh for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; 3692116Sjkh q[i] = fw; 3702116Sjkh } 3712116Sjkh jz += k; 3722116Sjkh goto recompute; 3732116Sjkh } 3742116Sjkh } 3752116Sjkh 3762116Sjkh /* chop off zero terms */ 3772116Sjkh if(z==0.0) { 3782116Sjkh jz -= 1; q0 -= 24; 3792116Sjkh while(iq[jz]==0) { jz--; q0-=24;} 3802116Sjkh } else { /* break z into 24-bit if necessary */ 3812116Sjkh z = scalbn(z,-q0); 382141296Sdas if(z>=two24) { 3832116Sjkh fw = (double)((int32_t)(twon24*z)); 3842116Sjkh iq[jz] = (int32_t)(z-two24*fw); 3852116Sjkh jz += 1; q0 += 24; 3862116Sjkh iq[jz] = (int32_t) fw; 3872116Sjkh } else iq[jz] = (int32_t) z ; 3882116Sjkh } 3892116Sjkh 3902116Sjkh /* convert integer "bit" chunk to floating-point value */ 3912116Sjkh fw = scalbn(one,q0); 3922116Sjkh for(i=jz;i>=0;i--) { 3932116Sjkh q[i] = fw*(double)iq[i]; fw*=twon24; 3942116Sjkh } 3952116Sjkh 3962116Sjkh /* compute PIo2[0,...,jp]*q[jz,...,0] */ 3972116Sjkh for(i=jz;i>=0;i--) { 3982116Sjkh for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k]; 3992116Sjkh fq[jz-i] = fw; 4002116Sjkh } 4012116Sjkh 4022116Sjkh /* compress fq[] into y[] */ 4032116Sjkh switch(prec) { 4042116Sjkh case 0: 4052116Sjkh fw = 0.0; 4062116Sjkh for (i=jz;i>=0;i--) fw += fq[i]; 407141296Sdas y[0] = (ih==0)? fw: -fw; 4082116Sjkh break; 4092116Sjkh case 1: 4102116Sjkh case 2: 4112116Sjkh fw = 0.0; 412141296Sdas for (i=jz;i>=0;i--) fw += fq[i]; 413175507Sbde STRICT_ASSIGN(double,fw,fw); 414141296Sdas y[0] = (ih==0)? fw: -fw; 4152116Sjkh fw = fq[0]-fw; 4162116Sjkh for (i=1;i<=jz;i++) fw += fq[i]; 417141296Sdas y[1] = (ih==0)? fw: -fw; 4182116Sjkh break; 4192116Sjkh case 3: /* painful */ 4202116Sjkh for (i=jz;i>0;i--) { 421141296Sdas fw = fq[i-1]+fq[i]; 4222116Sjkh fq[i] += fq[i-1]-fw; 4232116Sjkh fq[i-1] = fw; 4242116Sjkh } 4252116Sjkh for (i=jz;i>1;i--) { 426141296Sdas fw = fq[i-1]+fq[i]; 4272116Sjkh fq[i] += fq[i-1]-fw; 4282116Sjkh fq[i-1] = fw; 4292116Sjkh } 430141296Sdas for (fw=0.0,i=jz;i>=2;i--) fw += fq[i]; 4312116Sjkh if(ih==0) { 4322116Sjkh y[0] = fq[0]; y[1] = fq[1]; y[2] = fw; 4332116Sjkh } else { 4342116Sjkh y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw; 4352116Sjkh } 4362116Sjkh } 4372116Sjkh return n&7; 4382116Sjkh} 439