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__
38static const int two_over_pi[] = {
39#else
40static int two_over_pi[] = {
41#endif
420xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
430x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
440x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
450xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
460x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
470x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
480x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
490xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
500x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
510x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
520x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
53};
54
55#ifdef __STDC__
56static const int npio2_hw[] = {
57#else
58static int npio2_hw[] = {
59#endif
600x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
610x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
620x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
630x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
640x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
650x404858EB, 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__
79static const double
80#else
81static double
82#endif
83zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
84half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
85two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
86invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
87pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
88pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
89pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
90pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
91pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
92pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
93
94#ifdef __STDC__
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