src/fdlibm/e_rem_pio2.c

/* -*- Mode: C; tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
 *
 * ***** BEGIN LICENSE BLOCK *****
 * Version: MPL 1.1/GPL 2.0/LGPL 2.1
 *
 * The contents of this file are subject to the Mozilla Public License Version
 * 1.1 (the "License"); you may not use this file except in compliance with
 * the License. You may obtain a copy of the License at
 * http://www.mozilla.org/MPL/
 *
 * Software distributed under the License is distributed on an "AS IS" basis,
 * WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
 * for the specific language governing rights and limitations under the
 * License.
 *
 * The Original Code is Mozilla Communicator client code, released
 * March 31, 1998.
 *
 * The Initial Developer of the Original Code is
 * Sun Microsystems, Inc.
 * Portions created by the Initial Developer are Copyright (C) 1998
 * the Initial Developer. All Rights Reserved.
 *
 * Contributor(s):
 *
 * Alternatively, the contents of this file may be used under the terms of
 * either of the GNU General Public License Version 2 or later (the "GPL"),
 * or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
 * in which case the provisions of the GPL or the LGPL are applicable instead
 * of those above. If you wish to allow use of your version of this file only
 * under the terms of either the GPL or the LGPL, and not to allow others to
 * use your version of this file under the terms of the MPL, indicate your
 * decision by deleting the provisions above and replace them with the notice
 * and other provisions required by the GPL or the LGPL. If you do not delete
 * the provisions above, a recipient may use your version of this file under
 * the terms of any one of the MPL, the GPL or the LGPL.
 *
 * ***** END LICENSE BLOCK ***** */

/* @(#)e_rem_pio2.c 1.4 95/01/18 */
/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * ====================================================
 *
 */

/* __ieee754_rem_pio2(x,y)
 * 
 * return the remainder of x rem pi/2 in y[0]+y[1] 
 * use __kernel_rem_pio2()
 */

#include "fdlibm.h"

/*
 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi 
 */
#ifdef __STDC__
static const int two_over_pi[] = {
#else
static int two_over_pi[] = {
#endif
0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 
0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 
0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 
0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 
0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 
0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 
0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 
0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 
0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 
0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 
0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, 
};

#ifdef __STDC__
static const int npio2_hw[] = {
#else
static int npio2_hw[] = {
#endif
0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
0x404858EB, 0x404921FB,
};

/*
 * invpio2:  53 bits of 2/pi
 * pio2_1:   first  33 bit of pi/2
 * pio2_1t:  pi/2 - pio2_1
 * pio2_2:   second 33 bit of pi/2
 * pio2_2t:  pi/2 - (pio2_1+pio2_2)
 * pio2_3:   third  33 bit of pi/2
 * pio2_3t:  pi/2 - (pio2_1+pio2_2+pio2_3)
 */

#ifdef __STDC__
static const double 
#else
static double 
#endif
zero =  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
two24 =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
invpio2 =  6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
pio2_1  =  1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
pio2_1t =  6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
pio2_2  =  6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
pio2_2t =  2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
pio2_3  =  2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
pio2_3t =  8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */

#ifdef __STDC__
        int __ieee754_rem_pio2(double x, double *y)
#else
        int __ieee754_rem_pio2(x,y)
        double x,y[];
#endif
{
        fd_twoints u, ux, uz;
        double z = 0;
        double w,t,r,fn;
        double tx[3];
        int e0,i,j,nx,n,ix,hx;

        u.d = x;
        hx = __HI(u);           /* high word of x */
        ix = hx&0x7fffffff;
        if(ix<=0x3fe921fb)   /* |x| ~<= pi/4 , no need for reduction */
            {y[0] = x; y[1] = 0; return 0;}
        if(ix<0x4002d97c) {  /* |x| < 3pi/4, special case with n=+-1 */
            if(hx>0) { 
                z = x - pio2_1;
                if(ix!=0x3ff921fb) {    /* 33+53 bit pi is good enough */
                    y[0] = z - pio2_1t;
                    y[1] = (z-y[0])-pio2_1t;
                } else {                /* near pi/2, use 33+33+53 bit pi */
                    z -= pio2_2;
                    y[0] = z - pio2_2t;
                    y[1] = (z-y[0])-pio2_2t;
                }
                return 1;
            } else {    /* negative x */
                z = x + pio2_1;
                if(ix!=0x3ff921fb) {    /* 33+53 bit pi is good enough */
                    y[0] = z + pio2_1t;
                    y[1] = (z-y[0])+pio2_1t;
                } else {                /* near pi/2, use 33+33+53 bit pi */
                    z += pio2_2;
                    y[0] = z + pio2_2t;
                    y[1] = (z-y[0])+pio2_2t;
                }
                return -1;
            }
        }
        if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
            t  = fd_fabs(x);
            n  = (int) (t*invpio2+half);
            fn = (double)n;
            r  = t-fn*pio2_1;
            w  = fn*pio2_1t;    /* 1st round good to 85 bit */
            if(n<32&&ix!=npio2_hw[n-1]) {       
                y[0] = r-w;     /* quick check no cancellation */
            } else {
                j  = ix>>20;
                y[0] = r-w;
                u.d = y[0];
                i = j-(((__HI(u))>>20)&0x7ff);
                if(i>16) {  /* 2nd iteration needed, good to 118 */
                    t  = r;
                    w  = fn*pio2_2;     
                    r  = t-w;
                    w  = fn*pio2_2t-((t-r)-w);  
                    y[0] = r-w;
                    u.d = y[0];
                    i = j-(((__HI(u))>>20)&0x7ff);
                    if(i>49)  { /* 3rd iteration need, 151 bits acc */
                        t  = r; /* will cover all possible cases */
                        w  = fn*pio2_3; 
                        r  = t-w;
                        w  = fn*pio2_3t-((t-r)-w);      
                        y[0] = r-w;
                    }
                }
            }
            y[1] = (r-y[0])-w;
            if(hx<0)    {y[0] = -y[0]; y[1] = -y[1]; return -n;}
            else         return n;
        }
    /* 
     * all other (large) arguments
     */
        if(ix>=0x7ff00000) {            /* x is inf or NaN */
            y[0]=y[1]=x-x; return 0;
        }
    /* set z = scalbn(|x|,ilogb(x)-23) */
        ux.d = x; uz.d = z;
        __LO(uz) = __LO(ux);
        z = uz.d;
        e0      = (ix>>20)-1046;        /* e0 = ilogb(z)-23; */
        uz.d = z;
        __HI(uz) = ix - (e0<<20);
        z = uz.d;
        for(i=0;i<2;i++) {
                tx[i] = (double)((int)(z));
                z     = (z-tx[i])*two24;
        }
        tx[2] = z;
        nx = 3;
        while(tx[nx-1]==zero) nx--;     /* skip zero term */
        n  =  __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi);
        if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
        return n;
}
1	siliconforks	2	/* -- Mode: C; tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 --
2			*
3			* *** BEGIN LICENSE BLOCK ***
4			* Version: MPL 1.1/GPL 2.0/LGPL 2.1
5			*
6			* The contents of this file are subject to the Mozilla Public License Version
7			* 1.1 (the "License"); you may not use this file except in compliance with
8			* the License. You may obtain a copy of the License at
9			* http://www.mozilla.org/MPL/
10			*
11			* Software distributed under the License is distributed on an "AS IS" basis,
12			* WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
13			* for the specific language governing rights and limitations under the
14			* License.
15			*
16			* The Original Code is Mozilla Communicator client code, released
17			* March 31, 1998.
18			*
19			* The Initial Developer of the Original Code is
20			* Sun Microsystems, Inc.
21			* Portions created by the Initial Developer are Copyright (C) 1998
22			* the Initial Developer. All Rights Reserved.
23			*
24			* Contributor(s):
25			*
26			* Alternatively, the contents of this file may be used under the terms of
27			* either of the GNU General Public License Version 2 or later (the "GPL"),
28			* or the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
29			* in which case the provisions of the GPL or the LGPL are applicable instead
30			* of those above. If you wish to allow use of your version of this file only
31			* under the terms of either the GPL or the LGPL, and not to allow others to
32			* use your version of this file under the terms of the MPL, indicate your
33			* decision by deleting the provisions above and replace them with the notice
34			* and other provisions required by the GPL or the LGPL. If you do not delete
35			* the provisions above, a recipient may use your version of this file under
36			* the terms of any one of the MPL, the GPL or the LGPL.
37			*
38			* *** END LICENSE BLOCK *** */
39
40			/* @(#)e_rem_pio2.c 1.4 95/01/18 */
41			/*
42			* ====================================================
43			* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
44			*
45			* Developed at SunSoft, a Sun Microsystems, Inc. business.
46			* Permission to use, copy, modify, and distribute this
47			* software is freely granted, provided that this notice
48			* is preserved.
49			* ====================================================
50			*
51			*/
52
53			/* __ieee754_rem_pio2(x,y)
54			*
55			* return the remainder of x rem pi/2 in y[0]+y[1]
56			* use __kernel_rem_pio2()
57			*/
58
59			#include "fdlibm.h"
60
61			/*
62			* Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
63			*/
64			#ifdef __STDC__
65			static const int two_over_pi[] = {
66			#else
67			static int two_over_pi[] = {
68			#endif
69			0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
70			0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
71			0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
72			0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
73			0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
74			0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
75			0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
76			0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
77			0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
78			0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
79			0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
80			};
81
82			#ifdef __STDC__
83			static const int npio2_hw[] = {
84			#else
85			static int npio2_hw[] = {
86			#endif
87			0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
88			0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
89			0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
90			0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
91			0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
92			0x404858EB, 0x404921FB,
93			};
94
95			/*
96			* invpio2: 53 bits of 2/pi
97			* pio2_1: first 33 bit of pi/2
98			* pio2_1t: pi/2 - pio2_1
99			* pio2_2: second 33 bit of pi/2
100			* pio2_2t: pi/2 - (pio2_1+pio2_2)
101			* pio2_3: third 33 bit of pi/2
102			* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
103			*/
104
105			#ifdef __STDC__
106			static const double
107			#else
108			static double
109			#endif
110			zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
111			half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
112			two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
113			invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
114			pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
115			pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
116			pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
117			pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
118			pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
119			pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
120
121			#ifdef __STDC__
122			int __ieee754_rem_pio2(double x, double *y)
123			#else
124			int __ieee754_rem_pio2(x,y)
125			double x,y[];
126			#endif
127			{
128			fd_twoints u, ux, uz;
129			double z = 0;
130			double w,t,r,fn;
131			double tx[3];
132			int e0,i,j,nx,n,ix,hx;
133
134			u.d = x;
135			hx = __HI(u); /* high word of x */
136			ix = hx&0x7fffffff;
137			if(ix<=0x3fe921fb) /* \|x\| ~<= pi/4 , no need for reduction */
138			{y[0] = x; y[1] = 0; return 0;}
139			if(ix<0x4002d97c) { /* \|x\| < 3pi/4, special case with n=+-1 */
140			if(hx>0) {
141			z = x - pio2_1;
142			if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */
143			y[0] = z - pio2_1t;
144			y[1] = (z-y[0])-pio2_1t;
145			} else { /* near pi/2, use 33+33+53 bit pi */
146			z -= pio2_2;
147			y[0] = z - pio2_2t;
148			y[1] = (z-y[0])-pio2_2t;
149			}
150			return 1;
151			} else { /* negative x */
152			z = x + pio2_1;
153			if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */
154			y[0] = z + pio2_1t;
155			y[1] = (z-y[0])+pio2_1t;
156			} else { /* near pi/2, use 33+33+53 bit pi */
157			z += pio2_2;
158			y[0] = z + pio2_2t;
159			y[1] = (z-y[0])+pio2_2t;
160			}
161			return -1;
162			}
163			}
164			if(ix<=0x413921fb) { /* \|x\| ~<= 2^19(pi/2), medium size /
165			t = fd_fabs(x);
166			n = (int) (t*invpio2+half);
167			fn = (double)n;
168			r = t-fn*pio2_1;
169			w = fnpio2_1t; / 1st round good to 85 bit */
170			if(n<32&&ix!=npio2_hw[n-1]) {
171			y[0] = r-w; /* quick check no cancellation */
172			} else {
173			j = ix>>20;
174			y[0] = r-w;
175			u.d = y[0];
176			i = j-(((__HI(u))>>20)&0x7ff);
177			if(i>16) { /* 2nd iteration needed, good to 118 */
178			t = r;
179			w = fn*pio2_2;
180			r = t-w;
181			w = fn*pio2_2t-((t-r)-w);
182			y[0] = r-w;
183			u.d = y[0];
184			i = j-(((__HI(u))>>20)&0x7ff);
185			if(i>49) { /* 3rd iteration need, 151 bits acc */
186			t = r; /* will cover all possible cases */
187			w = fn*pio2_3;
188			r = t-w;
189			w = fn*pio2_3t-((t-r)-w);
190			y[0] = r-w;
191			}
192			}
193			}
194			y[1] = (r-y[0])-w;
195			if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
196			else return n;
197			}
198			/*
199			* all other (large) arguments
200			*/
201			if(ix>=0x7ff00000) { /* x is inf or NaN */
202			y[0]=y[1]=x-x; return 0;
203			}
204			/* set z = scalbn(\|x\|,ilogb(x)-23) */
205			ux.d = x; uz.d = z;
206			__LO(uz) = __LO(ux);
207			z = uz.d;
208			e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */
209			uz.d = z;
210			__HI(uz) = ix - (e0<<20);
211			z = uz.d;
212			for(i=0;i<2;i++) {
213			tx[i] = (double)((int)(z));
214			z = (z-tx[i])*two24;
215			}
216			tx[2] = z;
217			nx = 3;
218			while(tx[nx-1]==zero) nx--; /* skip zero term */
219			n = __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi);
220			if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
221			return n;
222			}