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