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