e_rem_pio2.c revision 2116
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 142116Sjkhstatic char rcsid[] = "$Id: e_rem_pio2.c,v 1.5 1994/08/18 23:05:56 jtc 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;} 1022116Sjkh if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */ 1032116Sjkh t = fabs(x); 1042116Sjkh n = (int32_t) (t*invpio2+half); 1052116Sjkh fn = (double)n; 1062116Sjkh r = t-fn*pio2_1; 1072116Sjkh w = fn*pio2_1t; /* 1st round good to 85 bit */ 1082116Sjkh if(n<32&&ix!=npio2_hw[n-1]) { 1092116Sjkh y[0] = r-w; /* quick check no cancellation */ 1102116Sjkh } else { 1112116Sjkh u_int32_t high; 1122116Sjkh j = ix>>20; 1132116Sjkh y[0] = r-w; 1142116Sjkh GET_HIGH_WORD(high,y[0]); 1152116Sjkh i = j-((high>>20)&0x7ff); 1162116Sjkh if(i>16) { /* 2nd iteration needed, good to 118 */ 1172116Sjkh t = r; 1182116Sjkh w = fn*pio2_2; 1192116Sjkh r = t-w; 1202116Sjkh w = fn*pio2_2t-((t-r)-w); 1212116Sjkh y[0] = r-w; 1222116Sjkh GET_HIGH_WORD(high,y[0]); 1232116Sjkh i = j-((high>>20)&0x7ff); 1242116Sjkh if(i>49) { /* 3rd iteration need, 151 bits acc */ 1252116Sjkh t = r; /* will cover all possible cases */ 1262116Sjkh w = fn*pio2_3; 1272116Sjkh r = t-w; 1282116Sjkh w = fn*pio2_3t-((t-r)-w); 1292116Sjkh y[0] = r-w; 1302116Sjkh } 1312116Sjkh } 1322116Sjkh } 1332116Sjkh y[1] = (r-y[0])-w; 1342116Sjkh if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 1352116Sjkh else return n; 1362116Sjkh } 1372116Sjkh /* 1382116Sjkh * all other (large) arguments 1392116Sjkh */ 1402116Sjkh if(ix>=0x7ff00000) { /* x is inf or NaN */ 1412116Sjkh y[0]=y[1]=x-x; return 0; 1422116Sjkh } 1432116Sjkh /* set z = scalbn(|x|,ilogb(x)-23) */ 1442116Sjkh GET_LOW_WORD(low,x); 1452116Sjkh SET_LOW_WORD(z,low); 1462116Sjkh e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ 1472116Sjkh SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20))); 1482116Sjkh for(i=0;i<2;i++) { 1492116Sjkh tx[i] = (double)((int32_t)(z)); 1502116Sjkh z = (z-tx[i])*two24; 1512116Sjkh } 1522116Sjkh tx[2] = z; 1532116Sjkh nx = 3; 1542116Sjkh while(tx[nx-1]==zero) nx--; /* skip zero term */ 1552116Sjkh n = __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi); 1562116Sjkh if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 1572116Sjkh return n; 1582116Sjkh} 159