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