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