uc-sdk
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
e_rem_pio2.c
Go to the documentation of this file.
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__
27 static const int two_over_pi[] = {
28 #else
29 static int two_over_pi[] = {
30 #endif
31 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
32 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
33 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
34 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
35 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
36 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
37 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
38 0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
39 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
40 0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
41 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
42 };
43 
44 #ifdef __STDC__
45 static const int npio2_hw[] = {
46 #else
47 static int npio2_hw[] = {
48 #endif
49 0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
50 0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
51 0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
52 0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
53 0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
54 0x404858EB, 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__
68 static const double
69 #else
70 static double
71 #endif
72 zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
73 half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
74 two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
75 invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
76 pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
77 pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
78 pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
79 pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
80 pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
81 pio2_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 = 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 = scalbn(|x|,ilogb(x)-23) */
162  __LO(z) = __LO(x);
163  e0 = (ix>>20)-1046; /* e0 = 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 }