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