1 /* 2 * Copyright (c) 1998, 2001, Oracle and/or its affiliates. All rights reserved. 3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. 4 * 5 * This code is free software; you can redistribute it and/or modify it 6 * under the terms of the GNU General Public License version 2 only, as 7 * published by the Free Software Foundation. Oracle designates this 8 * particular file as subject to the "Classpath" exception as provided 9 * by Oracle in the LICENSE file that accompanied this code. 10 * 11 * This code is distributed in the hope that it will be useful, but WITHOUT 12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 14 * version 2 for more details (a copy is included in the LICENSE file that 15 * accompanied this code). 16 * 17 * You should have received a copy of the GNU General Public License version 18 * 2 along with this work; if not, write to the Free Software Foundation, 19 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 20 * 21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA 22 * or visit www.oracle.com if you need additional information or have any 23 * questions. 24 */ 25 26 /* __ieee754_rem_pio2(x,y) 27 * 28 * return the remainder of x rem pi/2 in y[0]+y[1] 29 * use __kernel_rem_pio2() 30 */ 31 32 #include "fdlibm.h" 33 34 /* 35 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 36 */ 37 #ifdef __STDC__ 38 static const int two_over_pi[] = { 39 #else 40 static int two_over_pi[] = { 41 #endif 42 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 43 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 44 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 45 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 46 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 47 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 48 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 49 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 50 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 51 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 52 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 53 }; 54 55 #ifdef __STDC__ 56 static const int npio2_hw[] = { 57 #else 58 static int npio2_hw[] = { 59 #endif 60 0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, 61 0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C, 62 0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A, 63 0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C, 64 0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB, 65 0x404858EB, 0x404921FB, 66 }; 67 68 /* 69 * invpio2: 53 bits of 2/pi 70 * pio2_1: first 33 bit of pi/2 71 * pio2_1t: pi/2 - pio2_1 72 * pio2_2: second 33 bit of pi/2 73 * pio2_2t: pi/2 - (pio2_1+pio2_2) 74 * pio2_3: third 33 bit of pi/2 75 * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) 76 */ 77 78 #ifdef __STDC__ 79 static const double 80 #else 81 static double 82 #endif 83 zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ 84 half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ 85 two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ 86 invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */ 87 pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */ 88 pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */ 89 pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */ 90 pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */ 91 pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ 92 pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ 93 94 #ifdef __STDC__ __ieee754_rem_pio2(double x,double * y)95 int __ieee754_rem_pio2(double x, double *y) 96 #else 97 int __ieee754_rem_pio2(x,y) 98 double x,y[]; 99 #endif 100 { 101 double z,w,t,r,fn; 102 double tx[3]; 103 int e0,i,j,nx,n,ix,hx; 104 105 hx = __HI(x); /* high word of x */ 106 ix = hx&0x7fffffff; 107 if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ 108 {y[0] = x; y[1] = 0; return 0;} 109 if(ix<0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */ 110 if(hx>0) { 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 } else { /* negative x */ 122 z = x + pio2_1; 123 if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */ 124 y[0] = z + pio2_1t; 125 y[1] = (z-y[0])+pio2_1t; 126 } else { /* near pi/2, use 33+33+53 bit pi */ 127 z += pio2_2; 128 y[0] = z + pio2_2t; 129 y[1] = (z-y[0])+pio2_2t; 130 } 131 return -1; 132 } 133 } 134 if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */ 135 t = fabs(x); 136 n = (int) (t*invpio2+half); 137 fn = (double)n; 138 r = t-fn*pio2_1; 139 w = fn*pio2_1t; /* 1st round good to 85 bit */ 140 if(n<32&&ix!=npio2_hw[n-1]) { 141 y[0] = r-w; /* quick check no cancellation */ 142 } else { 143 j = ix>>20; 144 y[0] = r-w; 145 i = j-(((__HI(y[0]))>>20)&0x7ff); 146 if(i>16) { /* 2nd iteration needed, good to 118 */ 147 t = r; 148 w = fn*pio2_2; 149 r = t-w; 150 w = fn*pio2_2t-((t-r)-w); 151 y[0] = r-w; 152 i = j-(((__HI(y[0]))>>20)&0x7ff); 153 if(i>49) { /* 3rd iteration need, 151 bits acc */ 154 t = r; /* will cover all possible cases */ 155 w = fn*pio2_3; 156 r = t-w; 157 w = fn*pio2_3t-((t-r)-w); 158 y[0] = r-w; 159 } 160 } 161 } 162 y[1] = (r-y[0])-w; 163 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 164 else return n; 165 } 166 /* 167 * all other (large) arguments 168 */ 169 if(ix>=0x7ff00000) { /* x is inf or NaN */ 170 y[0]=y[1]=x-x; return 0; 171 } 172 /* set z = scalbn(|x|,ilogb(x)-23) */ 173 __LO(z) = __LO(x); 174 e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ 175 __HI(z) = ix - (e0<<20); 176 for(i=0;i<2;i++) { 177 tx[i] = (double)((int)(z)); 178 z = (z-tx[i])*two24; 179 } 180 tx[2] = z; 181 nx = 3; 182 while(tx[nx-1]==zero) nx--; /* skip zero term */ 183 n = __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi); 184 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} 185 return n; 186 } 187