e_rem_pio2l.h revision 7659
12116Sjkh/* @(#)e_rem_pio2.c 5.1 93/09/24 */ 22116Sjkh/* 32116Sjkh * ==================================================== 42116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 52116Sjkh * 62116Sjkh * Developed at SunPro, a Sun Microsystems, Inc. business. 72116Sjkh * Permission to use, copy, modify, and distribute this 82116Sjkh * software is freely granted, provided that this notice 92116Sjkh * is preserved. 102116Sjkh * ==================================================== 112116Sjkh */ 122116Sjkh 132116Sjkh#ifndef lint 147659Sbdestatic char rcsid[] = "$Id: e_rem_pio2.c,v 1.1.1.1 1994/08/19 09:39:44 jkh Exp $"; 152116Sjkh#endif 162116Sjkh 172116Sjkh/* __ieee754_rem_pio2(x,y) 182116Sjkh * 192116Sjkh * return the remainder of x rem pi/2 in y[0]+y[1] 202116Sjkh * use __kernel_rem_pio2() 212116Sjkh */ 222116Sjkh 232116Sjkh#include "math.h" 242116Sjkh#include "math_private.h" 252116Sjkh 262116Sjkh/* 272116Sjkh * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 282116Sjkh */ 292116Sjkh#ifdef __STDC__ 302116Sjkhstatic const int32_t two_over_pi[] = { 312116Sjkh#else 322116Sjkhstatic int32_t two_over_pi[] = { 332116Sjkh#endif 342116Sjkh0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 352116Sjkh0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 362116Sjkh0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 372116Sjkh0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 382116Sjkh0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 392116Sjkh0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 402116Sjkh0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 412116Sjkh0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 422116Sjkh0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 432116Sjkh0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 442116Sjkh0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 452116Sjkh}; 462116Sjkh 472116Sjkh#ifdef __STDC__ 482116Sjkhstatic const int32_t npio2_hw[] = { 492116Sjkh#else 502116Sjkhstatic int32_t npio2_hw[] = { 512116Sjkh#endif 522116Sjkh0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, 532116Sjkh0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C, 542116Sjkh0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A, 552116Sjkh0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C, 562116Sjkh0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB, 572116Sjkh0x404858EB, 0x404921FB, 582116Sjkh}; 592116Sjkh 602116Sjkh/* 612116Sjkh * invpio2: 53 bits of 2/pi 622116Sjkh * pio2_1: first 33 bit of pi/2 632116Sjkh * pio2_1t: pi/2 - pio2_1 642116Sjkh * pio2_2: second 33 bit of pi/2 652116Sjkh * pio2_2t: pi/2 - (pio2_1+pio2_2) 662116Sjkh * pio2_3: third 33 bit of pi/2 672116Sjkh * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) 682116Sjkh */ 692116Sjkh 702116Sjkh#ifdef __STDC__ 712116Sjkhstatic const double 722116Sjkh#else 732116Sjkhstatic double 742116Sjkh#endif 752116Sjkhzero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ 762116Sjkhhalf = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ 772116Sjkhtwo24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ 782116Sjkhinvpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ 792116Sjkhpio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ 802116Sjkhpio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ 812116Sjkhpio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */ 822116Sjkhpio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */ 832116Sjkhpio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ 842116Sjkhpio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ 852116Sjkh 862116Sjkh#ifdef __STDC__ 872116Sjkh int32_t __ieee754_rem_pio2(double x, double *y) 882116Sjkh#else 892116Sjkh int32_t __ieee754_rem_pio2(x,y) 902116Sjkh double x,y[]; 912116Sjkh#endif 922116Sjkh{ 932116Sjkh double z,w,t,r,fn; 942116Sjkh double tx[3]; 952116Sjkh int32_t e0,i,j,nx,n,ix,hx; 962116Sjkh u_int32_t low; 972116Sjkh 982116Sjkh GET_HIGH_WORD(hx,x); /* high word of x */ 992116Sjkh ix = hx&0x7fffffff; 1002116Sjkh if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ 1012116Sjkh {y[0] = x; y[1] = 0; return 0;} 1027659Sbde if(ix<0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */ 1037659Sbde if(hx>0) { 1047659Sbde z = x - pio2_1; 1057659Sbde if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */ 1067659Sbde y[0] = z - pio2_1t; 1077659Sbde y[1] = (z-y[0])-pio2_1t; 1087659Sbde } else { /* near pi/2, use 33+33+53 bit pi */ 1097659Sbde z -= pio2_2; 1107659Sbde y[0] = z - pio2_2t; 1117659Sbde y[1] = (z-y[0])-pio2_2t; 1127659Sbde } 1137659Sbde return 1; 1147659Sbde } else { /* negative x */ 1157659Sbde z = x + pio2_1; 1167659Sbde if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */ 1177659Sbde y[0] = z + pio2_1t; 1187659Sbde y[1] = (z-y[0])+pio2_1t; 1197659Sbde } else { /* near pi/2, use 33+33+53 bit pi */ 1207659Sbde z += pio2_2; 1217659Sbde y[0] = z + pio2_2t; 1227659Sbde y[1] = (z-y[0])+pio2_2t; 1237659Sbde } 1247659Sbde return -1; 1257659Sbde } 1267659Sbde } 1272116Sjkh if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */ 1282116Sjkh t = fabs(x); 1292116Sjkh n = (int32_t) (t*invpio2+half); 1302116Sjkh fn = (double)n; 1312116Sjkh r = t-fn*pio2_1; 1322116Sjkh w = fn*pio2_1t; /* 1st round good to 85 bit */ 1332116Sjkh if(n<32&&ix!=npio2_hw[n-1]) { 1342116Sjkh y[0] = r-w; /* quick check no cancellation */ 1352116Sjkh } else { 1362116Sjkh u_int32_t high; 1372116Sjkh j = ix>>20; 1382116Sjkh y[0] = r-w; 1392116Sjkh GET_HIGH_WORD(high,y[0]); 1402116Sjkh i = j-((high>>20)&0x7ff); 1412116Sjkh if(i>16) { /* 2nd iteration needed, good to 118 */ 1422116Sjkh t = r; 1432116Sjkh w = fn*pio2_2; 1442116Sjkh r = t-w; 1452116Sjkh w = fn*pio2_2t-((t-r)-w); 1462116Sjkh y[0] = r-w; 1472116Sjkh GET_HIGH_WORD(high,y[0]); 1482116Sjkh i = j-((high>>20)&0x7ff); 1492116Sjkh if(i>49) { /* 3rd iteration need, 151 bits acc */ 1502116Sjkh t = r; /* will cover all possible cases */ 1512116Sjkh w = fn*pio2_3; 1522116Sjkh r = t-w; 1532116Sjkh w = fn*pio2_3t-((t-r)-w); 1542116Sjkh y[0] = r-w; 1552116Sjkh } 1562116Sjkh } 1572116Sjkh } 1582116Sjkh y[1] = (r-y[0])-w; 1592116Sjkh if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 1602116Sjkh else return n; 1612116Sjkh } 1622116Sjkh /* 1632116Sjkh * all other (large) arguments 1642116Sjkh */ 1652116Sjkh if(ix>=0x7ff00000) { /* x is inf or NaN */ 1662116Sjkh y[0]=y[1]=x-x; return 0; 1672116Sjkh } 1682116Sjkh /* set z = scalbn(|x|,ilogb(x)-23) */ 1692116Sjkh GET_LOW_WORD(low,x); 1702116Sjkh SET_LOW_WORD(z,low); 1712116Sjkh e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ 1722116Sjkh SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20))); 1732116Sjkh for(i=0;i<2;i++) { 1742116Sjkh tx[i] = (double)((int32_t)(z)); 1752116Sjkh z = (z-tx[i])*two24; 1762116Sjkh } 1772116Sjkh tx[2] = z; 1782116Sjkh nx = 3; 1792116Sjkh while(tx[nx-1]==zero) nx--; /* skip zero term */ 1802116Sjkh n = __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi); 1812116Sjkh if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 1822116Sjkh return n; 1832116Sjkh} 184