SDL  2.0
e_rem_pio2.c
Go to the documentation of this file.
1 /* @(#)e_rem_pio2.c 5.1 93/09/24 */
2 /*
3  * ====================================================
4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5  *
6  * Developed at SunPro, a Sun Microsystems, Inc. business.
7  * Permission to use, copy, modify, and distribute this
8  * software is freely granted, provided that this notice
9  * is preserved.
10  * ====================================================
11  */
12 
13 #if defined(LIBM_SCCS) && !defined(lint)
14 static const char rcsid[] =
15  "$NetBSD: e_rem_pio2.c,v 1.8 1995/05/10 20:46:02 jtc Exp $";
16 #endif
17 
18 /* __ieee754_rem_pio2(x,y)
19  *
20  * return the remainder of x rem pi/2 in y[0]+y[1]
21  * use __kernel_rem_pio2()
22  */
23 
24 #include "math_libm.h"
25 #include "math_private.h"
26 
27 #include "SDL_assert.h"
28 
30 
31 /*
32  * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
33  */
34 #ifdef __STDC__
35  static const int32_t two_over_pi[] = {
36 #else
37  static int32_t two_over_pi[] = {
38 #endif
39  0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
40  0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
41  0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
42  0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
43  0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
44  0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
45  0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
46  0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
47  0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
48  0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
49  0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
50  };
51 
52 #ifdef __STDC__
53 static const int32_t npio2_hw[] = {
54 #else
55 static int32_t npio2_hw[] = {
56 #endif
57  0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
58  0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
59  0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
60  0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
61  0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
62  0x404858EB, 0x404921FB,
63 };
64 
65 /*
66  * invpio2: 53 bits of 2/pi
67  * pio2_1: first 33 bit of pi/2
68  * pio2_1t: pi/2 - pio2_1
69  * pio2_2: second 33 bit of pi/2
70  * pio2_2t: pi/2 - (pio2_1+pio2_2)
71  * pio2_3: third 33 bit of pi/2
72  * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
73  */
74 
75 #ifdef __STDC__
76 static const double
77 #else
78 static double
79 #endif
80  zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
81  half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
82  two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
83  invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
84  pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
85  pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
86  pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
87  pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
88  pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
89  pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
90 
91 #ifdef __STDC__
93 __ieee754_rem_pio2(double x, double *y)
94 #else
97  double x, y[];
98 #endif
99 {
100  double z = 0.0, w, t, r, fn;
101  double tx[3];
102  int32_t e0, i, j, nx, n, ix, hx;
103  u_int32_t low;
104 
105  GET_HIGH_WORD(hx, x); /* high word of x */
106  ix = hx & 0x7fffffff;
107  if (ix <= 0x3fe921fb) { /* |x| ~<= pi/4 , no need for reduction */
108  y[0] = x;
109  y[1] = 0;
110  return 0;
111  }
112  if (ix < 0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */
113  if (hx > 0) {
114  z = x - pio2_1;
115  if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */
116  y[0] = z - pio2_1t;
117  y[1] = (z - y[0]) - pio2_1t;
118  } else { /* near pi/2, use 33+33+53 bit pi */
119  z -= pio2_2;
120  y[0] = z - pio2_2t;
121  y[1] = (z - y[0]) - pio2_2t;
122  }
123  return 1;
124  } else { /* negative x */
125  z = x + pio2_1;
126  if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */
127  y[0] = z + pio2_1t;
128  y[1] = (z - y[0]) + pio2_1t;
129  } else { /* near pi/2, use 33+33+53 bit pi */
130  z += pio2_2;
131  y[0] = z + pio2_2t;
132  y[1] = (z - y[0]) + pio2_2t;
133  }
134  return -1;
135  }
136  }
137  if (ix <= 0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
138  t = fabs(x);
139  n = (int32_t) (t * invpio2 + half);
140  fn = (double) n;
141  r = t - fn * pio2_1;
142  w = fn * pio2_1t; /* 1st round good to 85 bit */
143  if (n < 32 && ix != npio2_hw[n - 1]) {
144  y[0] = r - w; /* quick check no cancellation */
145  } else {
146  u_int32_t high;
147  j = ix >> 20;
148  y[0] = r - w;
149  GET_HIGH_WORD(high, y[0]);
150  i = j - ((high >> 20) & 0x7ff);
151  if (i > 16) { /* 2nd iteration needed, good to 118 */
152  t = r;
153  w = fn * pio2_2;
154  r = t - w;
155  w = fn * pio2_2t - ((t - r) - w);
156  y[0] = r - w;
157  GET_HIGH_WORD(high, y[0]);
158  i = j - ((high >> 20) & 0x7ff);
159  if (i > 49) { /* 3rd iteration need, 151 bits acc */
160  t = r; /* will cover all possible cases */
161  w = fn * pio2_3;
162  r = t - w;
163  w = fn * pio2_3t - ((t - r) - w);
164  y[0] = r - w;
165  }
166  }
167  }
168  y[1] = (r - y[0]) - w;
169  if (hx < 0) {
170  y[0] = -y[0];
171  y[1] = -y[1];
172  return -n;
173  } else
174  return n;
175  }
176  /*
177  * all other (large) arguments
178  */
179  if (ix >= 0x7ff00000) { /* x is inf or NaN */
180  y[0] = y[1] = x - x;
181  return 0;
182  }
183  /* set z = scalbn(|x|,ilogb(x)-23) */
184  GET_LOW_WORD(low, x);
185  SET_LOW_WORD(z, low);
186  e0 = (ix >> 20) - 1046; /* e0 = ilogb(z)-23; */
187  SET_HIGH_WORD(z, ix - ((int32_t) (e0 << 20)));
188  for (i = 0; i < 2; i++) {
189  tx[i] = (double) ((int32_t) (z));
190  z = (z - tx[i]) * two24;
191  }
192  tx[2] = z;
193  nx = 3;
194 
195  /* If this assertion ever fires, here's the static analysis that warned about it:
196  http://buildbot.libsdl.org/sdl-static-analysis/sdl-macosx-static-analysis/sdl-macosx-static-analysis-1101/report-8c6ccb.html#EndPath */
197  SDL_assert((tx[0] != zero) || (tx[1] != zero) || (tx[2] != zero));
198 
199  while (nx && tx[nx - 1] == zero)
200  nx--; /* skip zero term */
201  n = __kernel_rem_pio2(tx, y, e0, nx, 2, two_over_pi);
202  if (hx < 0) {
203  y[0] = -y[0];
204  y[1] = -y[1];
205  return -n;
206  }
207  return n;
208 }
int32_t attribute_hidden __ieee754_rem_pio2(double x, y)
Definition: e_rem_pio2.c:96
#define GET_HIGH_WORD(i, d)
Definition: math_private.h:103
GLdouble GLdouble GLdouble r
Definition: SDL_opengl.h:2079
GLdouble GLdouble z
static double invpio2
Definition: e_rem_pio2.c:83
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1574
signed int int32_t
int attribute_hidden __kernel_rem_pio2(x, y, int e0, int nx, int prec, ipio2)
Definition: k_rem_pio2.c:176
static double half
Definition: e_rem_pio2.c:81
static double two24
Definition: e_rem_pio2.c:82
static double pio2_1t
Definition: e_rem_pio2.c:85
#define attribute_hidden
Definition: math_private.h:24
#define SET_HIGH_WORD(d, v)
Definition: math_private.h:131
static double pio2_2
Definition: e_rem_pio2.c:86
static double pio2_3t
Definition: e_rem_pio2.c:89
return Display return Display Bool Bool int int int return Display XEvent Bool(*) XPointer return Display return Display Drawable _Xconst char unsigned int unsigned int return Display Pixmap Pixmap XColor XColor unsigned int unsigned int return Display _Xconst char char int char return Display Visual unsigned int int int char unsigned int unsigned int in i)
Definition: SDL_x11sym.h:50
unsigned int u_int32_t
Definition: math_private.h:29
static double pio2_1
Definition: e_rem_pio2.c:84
#define SET_LOW_WORD(d, v)
Definition: math_private.h:141
GLbyte nx
GLubyte GLubyte GLubyte GLubyte w
return Display return Display Bool Bool int int int return Display XEvent Bool(*) XPointer return Display return Display Drawable _Xconst char unsigned int unsigned int return Display Pixmap Pixmap XColor XColor unsigned int unsigned int return Display _Xconst char char int char return Display Visual unsigned int int int char unsigned int unsigned int int in j)
Definition: SDL_x11sym.h:50
GLint GLint GLint GLint GLint GLint y
Definition: SDL_opengl.h:1574
static double pio2_3
Definition: e_rem_pio2.c:88
#define fabs
Definition: math_private.h:36
#define SDL_assert(condition)
Definition: SDL_assert.h:169
static double pio2_2t
Definition: e_rem_pio2.c:87
#define GET_LOW_WORD(i, d)
Definition: math_private.h:112
GLdouble n
static int32_t npio2_hw[]
Definition: e_rem_pio2.c:55
libm_hidden_proto(fabs)
Definition: e_rem_pio2.c:29
static double zero
Definition: e_rem_pio2.c:80
GLdouble GLdouble t
Definition: SDL_opengl.h:2071