SDL  2.0
k_rem_pio2.c
Go to the documentation of this file.
1 /* @(#)k_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: k_rem_pio2.c,v 1.7 1995/05/10 20:46:25 jtc Exp $";
16 #endif
17 
18 /*
19  * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
20  * double x[],y[]; int e0,nx,prec; int ipio2[];
21  *
22  * __kernel_rem_pio2 return the last three digits of N with
23  * y = x - N*pi/2
24  * so that |y| < pi/2.
25  *
26  * The method is to compute the integer (mod 8) and fraction parts of
27  * (2/pi)*x without doing the full multiplication. In general we
28  * skip the part of the product that are known to be a huge integer (
29  * more accurately, = 0 mod 8 ). Thus the number of operations are
30  * independent of the exponent of the input.
31  *
32  * (2/pi) is represented by an array of 24-bit integers in ipio2[].
33  *
34  * Input parameters:
35  * x[] The input value (must be positive) is broken into nx
36  * pieces of 24-bit integers in double precision format.
37  * x[i] will be the i-th 24 bit of x. The scaled exponent
38  * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
39  * match x's up to 24 bits.
40  *
41  * Example of breaking a double positive z into x[0]+x[1]+x[2]:
42  * e0 = ilogb(z)-23
43  * z = scalbn(z,-e0)
44  * for i = 0,1,2
45  * x[i] = floor(z)
46  * z = (z-x[i])*2**24
47  *
48  *
49  * y[] ouput result in an array of double precision numbers.
50  * The dimension of y[] is:
51  * 24-bit precision 1
52  * 53-bit precision 2
53  * 64-bit precision 2
54  * 113-bit precision 3
55  * The actual value is the sum of them. Thus for 113-bit
56  * precison, one may have to do something like:
57  *
58  * long double t,w,r_head, r_tail;
59  * t = (long double)y[2] + (long double)y[1];
60  * w = (long double)y[0];
61  * r_head = t+w;
62  * r_tail = w - (r_head - t);
63  *
64  * e0 The exponent of x[0]
65  *
66  * nx dimension of x[]
67  *
68  * prec an integer indicating the precision:
69  * 0 24 bits (single)
70  * 1 53 bits (double)
71  * 2 64 bits (extended)
72  * 3 113 bits (quad)
73  *
74  * ipio2[]
75  * integer array, contains the (24*i)-th to (24*i+23)-th
76  * bit of 2/pi after binary point. The corresponding
77  * floating value is
78  *
79  * ipio2[i] * 2^(-24(i+1)).
80  *
81  * External function:
82  * double scalbn(), floor();
83  *
84  *
85  * Here is the description of some local variables:
86  *
87  * jk jk+1 is the initial number of terms of ipio2[] needed
88  * in the computation. The recommended value is 2,3,4,
89  * 6 for single, double, extended,and quad.
90  *
91  * jz local integer variable indicating the number of
92  * terms of ipio2[] used.
93  *
94  * jx nx - 1
95  *
96  * jv index for pointing to the suitable ipio2[] for the
97  * computation. In general, we want
98  * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
99  * is an integer. Thus
100  * e0-3-24*jv >= 0 or (e0-3)/24 >= jv
101  * Hence jv = max(0,(e0-3)/24).
102  *
103  * jp jp+1 is the number of terms in PIo2[] needed, jp = jk.
104  *
105  * q[] double array with integral value, representing the
106  * 24-bits chunk of the product of x and 2/pi.
107  *
108  * q0 the corresponding exponent of q[0]. Note that the
109  * exponent for q[i] would be q0-24*i.
110  *
111  * PIo2[] double precision array, obtained by cutting pi/2
112  * into 24 bits chunks.
113  *
114  * f[] ipio2[] in floating point
115  *
116  * iq[] integer array by breaking up q[] in 24-bits chunk.
117  *
118  * fq[] final product of x*(2/pi) in fq[0],..,fq[jk]
119  *
120  * ih integer. If >0 it indicates q[] is >= 0.5, hence
121  * it also indicates the *sign* of the result.
122  *
123  */
124 
125 
126 /*
127  * Constants:
128  * The hexadecimal values are the intended ones for the following
129  * constants. The decimal values may be used, provided that the
130  * compiler will convert from decimal to binary accurately enough
131  * to produce the hexadecimal values shown.
132  */
133 
134 #include "math_libm.h"
135 #include "math_private.h"
136 
137 #include "SDL_assert.h"
138 
141 #ifdef __STDC__
142  static const int init_jk[] = { 2, 3, 4, 6 }; /* initial value for jk */
143 #else
144  static int init_jk[] = { 2, 3, 4, 6 };
145 #endif
146 
147 #ifdef __STDC__
148 static const double PIo2[] = {
149 #else
150 static double PIo2[] = {
151 #endif
152  1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
153  7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
154  5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
155  3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
156  1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
157  1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
158  2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
159  2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
160 };
161 
162 #ifdef __STDC__
163 static const double
164 #else
165 static double
166 #endif
167  zero = 0.0, one = 1.0, two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
168  twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
169 
170 #ifdef __STDC__
172 __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec,
173  const int32_t * ipio2)
174 #else
176 __kernel_rem_pio2(x, y, e0, nx, prec, ipio2)
177  double x[], y[];
178  int e0, nx, prec;
179  int32_t ipio2[];
180 #endif
181 {
182  int32_t jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih;
183  double z, fw, f[20], fq[20], q[20];
184 
185  /* initialize jk */
186  SDL_assert((prec >= 0) && (prec < SDL_arraysize(init_jk)));
187  jk = init_jk[prec];
188  SDL_assert((jk >= 2) && (jk <= 6));
189  jp = jk;
190 
191  /* determine jx,jv,q0, note that 3>q0 */
192  SDL_assert(nx > 0);
193  jx = nx - 1;
194  jv = (e0 - 3) / 24;
195  if (jv < 0)
196  jv = 0;
197  q0 = e0 - 24 * (jv + 1);
198 
199  /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
200  j = jv - jx;
201  m = jx + jk;
202  for (i = 0; i <= m; i++, j++)
203  f[i] = (j < 0) ? zero : (double) ipio2[j];
204 
205  /* compute q[0],q[1],...q[jk] */
206  for (i = 0; i <= jk; i++) {
207  for (j = 0, fw = 0.0; j <= jx; j++) {
208  const int32_t idx = jx + i - j;
209  SDL_assert(idx >= 0);
210  SDL_assert(idx < 20);
211  SDL_assert(idx <= m);
212  fw += x[j] * f[idx];
213  }
214  q[i] = fw;
215  }
216 
217  jz = jk;
218  recompute:
219  /* distill q[] into iq[] reversingly */
220  for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) {
221  fw = (double) ((int32_t) (twon24 * z));
222  iq[i] = (int32_t) (z - two24 * fw);
223  z = q[j - 1] + fw;
224  }
225 
226  /* compute n */
227  z = scalbn(z, q0); /* actual value of z */
228  z -= 8.0 * floor(z * 0.125); /* trim off integer >= 8 */
229  n = (int32_t) z;
230  z -= (double) n;
231  ih = 0;
232  if (q0 > 0) { /* need iq[jz-1] to determine n */
233  i = (iq[jz - 1] >> (24 - q0));
234  n += i;
235  iq[jz - 1] -= i << (24 - q0);
236  ih = iq[jz - 1] >> (23 - q0);
237  } else if (q0 == 0)
238  ih = iq[jz - 1] >> 23;
239  else if (z >= 0.5)
240  ih = 2;
241 
242  if (ih > 0) { /* q > 0.5 */
243  n += 1;
244  carry = 0;
245  for (i = 0; i < jz; i++) { /* compute 1-q */
246  j = iq[i];
247  if (carry == 0) {
248  if (j != 0) {
249  carry = 1;
250  iq[i] = 0x1000000 - j;
251  }
252  } else
253  iq[i] = 0xffffff - j;
254  }
255  if (q0 > 0) { /* rare case: chance is 1 in 12 */
256  switch (q0) {
257  case 1:
258  iq[jz - 1] &= 0x7fffff;
259  break;
260  case 2:
261  iq[jz - 1] &= 0x3fffff;
262  break;
263  }
264  }
265  if (ih == 2) {
266  z = one - z;
267  if (carry != 0)
268  z -= scalbn(one, q0);
269  }
270  }
271 
272  /* check if recomputation is needed */
273  if (z == zero) {
274  j = 0;
275  for (i = jz - 1; i >= jk; i--)
276  j |= iq[i];
277  if (j == 0) { /* need recomputation */
278  for (k = 1; iq[jk - k] == 0; k++); /* k = no. of terms needed */
279 
280  for (i = jz + 1; i <= jz + k; i++) { /* add q[jz+1] to q[jz+k] */
281  f[jx + i] = (double) ipio2[jv + i];
282  for (j = 0, fw = 0.0; j <= jx; j++)
283  fw += x[j] * f[jx + i - j];
284  q[i] = fw;
285  }
286  jz += k;
287  goto recompute;
288  }
289  }
290 
291  /* chop off zero terms */
292  if (z == 0.0) {
293  jz -= 1;
294  q0 -= 24;
295  while (iq[jz] == 0) {
296  jz--;
297  q0 -= 24;
298  }
299  } else { /* break z into 24-bit if necessary */
300  z = scalbn(z, -q0);
301  if (z >= two24) {
302  fw = (double) ((int32_t) (twon24 * z));
303  iq[jz] = (int32_t) (z - two24 * fw);
304  jz += 1;
305  q0 += 24;
306  iq[jz] = (int32_t) fw;
307  } else
308  iq[jz] = (int32_t) z;
309  }
310 
311  /* convert integer "bit" chunk to floating-point value */
312  fw = scalbn(one, q0);
313  for (i = jz; i >= 0; i--) {
314  q[i] = fw * (double) iq[i];
315  fw *= twon24;
316  }
317 
318  /* compute PIo2[0,...,jp]*q[jz,...,0] */
319  for (i = jz; i >= 0; i--) {
320  for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++)
321  fw += PIo2[k] * q[i + k];
322  fq[jz - i] = fw;
323  }
324 
325  /* compress fq[] into y[] */
326  switch (prec) {
327  case 0:
328  fw = 0.0;
329  for (i = jz; i >= 0; i--)
330  fw += fq[i];
331  y[0] = (ih == 0) ? fw : -fw;
332  break;
333  case 1:
334  case 2:
335  fw = 0.0;
336  for (i = jz; i >= 0; i--)
337  fw += fq[i];
338  y[0] = (ih == 0) ? fw : -fw;
339  fw = fq[0] - fw;
340  for (i = 1; i <= jz; i++)
341  fw += fq[i];
342  y[1] = (ih == 0) ? fw : -fw;
343  break;
344  case 3: /* painful */
345  for (i = jz; i > 0; i--) {
346  fw = fq[i - 1] + fq[i];
347  fq[i] += fq[i - 1] - fw;
348  fq[i - 1] = fw;
349  }
350  for (i = jz; i > 1; i--) {
351  fw = fq[i - 1] + fq[i];
352  fq[i] += fq[i - 1] - fw;
353  fq[i - 1] = fw;
354  }
355  for (fw = 0.0, i = jz; i >= 2; i--)
356  fw += fq[i];
357  if (ih == 0) {
358  y[0] = fq[0];
359  y[1] = fq[1];
360  y[2] = fw;
361  } else {
362  y[0] = -fq[0];
363  y[1] = -fq[1];
364  y[2] = -fw;
365  }
366  }
367  return n & 7;
368 }
GLdouble GLdouble z
static double PIo2[]
Definition: k_rem_pio2.c:150
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1574
signed int int32_t
GLdouble GLdouble GLdouble GLdouble q
Definition: SDL_opengl.h:2087
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 int return Display Window Cursor return Display Window return Display Drawable GC int int unsigned int unsigned int return Display Drawable GC int int _Xconst char int return Display Drawable GC int int unsigned int unsigned int return Display return Display Cursor return Display GC return XModifierKeymap return char Display Window int return Display return Display Atom return Display Window XWindowAttributes return Display Window return Display XEvent Bool(*) XPointer return Display Window Bool unsigned int int int Window Cursor Time return Display Window int return KeySym return Display _Xconst char Bool return Display _Xconst char return XKeyEvent char int KeySym XComposeStatus return Display int int int XVisualInfo return Display Window int int return _Xconst char return Display XEvent return Display Drawable GC XImage int int int int unsigned int unsigned int return Display Window Window Window int int int int unsigned int return Display Window Window int int return Display Window unsigned int unsigned int return Display Window Bool long XEvent return Display GC unsigned long return Display Window int Time return Display Window Window return Display Window unsigned long return Display Window XSizeHints Display Colormap XColor int return char int XTextProperty return XFontStruct _Xconst char int int int int XCharStruct return Display Window return Display Time return Display Colormap return Display Window Window int int unsigned int unsigned int int int return Display Window int return XExtensionInfo Display char XExtensionHooks int XPointer return XExtensionInfo XExtensionInfo Display return Display return Display unsigned long Display GC Display char long Display xReply int Bool return Display Bool return Display int SDL_X11_XESetEventToWireRetType return Display Window Window Window Window unsigned int return Display XShmSegmentInfo return Display Drawable GC XImage int int int int unsigned int unsigned int Boo k)
Definition: SDL_x11sym.h:213
int attribute_hidden __kernel_rem_pio2(x, y, int e0, int nx, int prec, ipio2)
Definition: k_rem_pio2.c:176
const GLfloat * m
GLfloat f
#define attribute_hidden
Definition: math_private.h:24
static double one
Definition: k_rem_pio2.c:167
#define scalbn
Definition: math_private.h:40
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
static double two24
Definition: k_rem_pio2.c:167
GLbyte nx
static double twon24
Definition: k_rem_pio2.c:168
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
#define SDL_assert(condition)
Definition: SDL_assert.h:169
static double zero
Definition: k_rem_pio2.c:167
GLdouble n
#define SDL_arraysize(array)
Definition: SDL_stdinc.h:93
#define floor
Definition: math_private.h:37
libm_hidden_proto(scalbn)
Definition: k_rem_pio2.c:139