1 2 /* @(#)e_rem_pio2.c 1.4 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 15 /* __ieee754_rem_pio2(x,y) 16 * 17 * return the remainder of x rem pi/2 in y[0]+y[1] 18 * use __kernel_rem_pio2() 19 */ 20 21 #include "fdlibm.h" 22 23 /* 24 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 25 */ 26 #ifdef __STDC__ 27 static const int two_over_pi[] = { 28 #else 29 static int two_over_pi[] = { 30 #endif 31 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 32 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 33 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 34 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 35 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 36 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 37 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 38 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 39 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 40 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 41 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 42 }; 43 44 #ifdef __STDC__ 45 static const int npio2_hw[] = { 46 #else 47 static int npio2_hw[] = { 48 #endif 49 0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, 50 0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C, 51 0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A, 52 0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C, 53 0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB, 54 0x404858EB, 0x404921FB, 55 }; 56 57 /* 58 * invpio2: 53 bits of 2/pi 59 * pio2_1: first 33 bit of pi/2 60 * pio2_1t: pi/2 - pio2_1 61 * pio2_2: second 33 bit of pi/2 62 * pio2_2t: pi/2 - (pio2_1+pio2_2) 63 * pio2_3: third 33 bit of pi/2 64 * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) 65 */ 66 67 #ifdef __STDC__ 68 static const double 69 #else 70 static double 71 #endif 72 zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ 73 half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ 74 two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ 75 invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ 76 pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ 77 pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ 78 pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */ 79 pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */ 80 pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ 81 pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ 82 83 #ifdef __STDC__ 84 int __ieee754_rem_pio2(double x, double *y) 85 #else 86 int __ieee754_rem_pio2(x,y) 87 double x,y[]; 88 #endif 89 { 90 double z,w,t,r,fn; 91 double tx[3]; 92 int e0,i,j,nx,n,ix,hx; 93 94 hx = __HI(x); /* high word of x */ 95 ix = hx&0x7fffffff; 96 if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ 97 {y[0] = x; y[1] = 0; return 0;} 98 if(ix<0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */ 99 if(hx>0) { 100 z = x - pio2_1; 101 if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */ 102 y[0] = z - pio2_1t; 103 y[1] = (z-y[0])-pio2_1t; 104 } else { /* near pi/2, use 33+33+53 bit pi */ 105 z -= pio2_2; 106 y[0] = z - pio2_2t; 107 y[1] = (z-y[0])-pio2_2t; 108 } 109 return 1; 110 } else { /* negative x */ 111 z = x + pio2_1; 112 if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */ 113 y[0] = z + pio2_1t; 114 y[1] = (z-y[0])+pio2_1t; 115 } else { /* near pi/2, use 33+33+53 bit pi */ 116 z += pio2_2; 117 y[0] = z + pio2_2t; 118 y[1] = (z-y[0])+pio2_2t; 119 } 120 return -1; 121 } 122 } 123 if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */ 124 t = ieee_fabs(x); 125 n = (int) (t*invpio2+half); 126 fn = (double)n; 127 r = t-fn*pio2_1; 128 w = fn*pio2_1t; /* 1st round good to 85 bit */ 129 if(n<32&&ix!=npio2_hw[n-1]) { 130 y[0] = r-w; /* quick check no cancellation */ 131 } else { 132 j = ix>>20; 133 y[0] = r-w; 134 i = j-(((__HI(y[0]))>>20)&0x7ff); 135 if(i>16) { /* 2nd iteration needed, good to 118 */ 136 t = r; 137 w = fn*pio2_2; 138 r = t-w; 139 w = fn*pio2_2t-((t-r)-w); 140 y[0] = r-w; 141 i = j-(((__HI(y[0]))>>20)&0x7ff); 142 if(i>49) { /* 3rd iteration need, 151 bits acc */ 143 t = r; /* will cover all possible cases */ 144 w = fn*pio2_3; 145 r = t-w; 146 w = fn*pio2_3t-((t-r)-w); 147 y[0] = r-w; 148 } 149 } 150 } 151 y[1] = (r-y[0])-w; 152 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 153 else return n; 154 } 155 /* 156 * all other (large) arguments 157 */ 158 if(ix>=0x7ff00000) { /* x is inf or NaN */ 159 y[0]=y[1]=x-x; return 0; 160 } 161 /* set z = ieee_scalbn(|x|,ilogb(x)-23) */ 162 __LO(z) = __LO(x); 163 e0 = (ix>>20)-1046; /* e0 = ieee_ilogb(z)-23; */ 164 __HI(z) = ix - (e0<<20); 165 for(i=0;i<2;i++) { 166 tx[i] = (double)((int)(z)); 167 z = (z-tx[i])*two24; 168 } 169 tx[2] = z; 170 nx = 3; 171 while(tx[nx-1]==zero) nx--; /* skip zero term */ 172 n = __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi); 173 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 174 return n; 175 } 176
