SDL  2.0
k_rem_pio2.c File Reference
#include "math_libm.h"
#include "math_private.h"
#include "SDL_assert.h"
+ Include dependency graph for k_rem_pio2.c:

Go to the source code of this file.

Functions

 libm_hidden_proto (scalbn)
 
int attribute_hidden __kernel_rem_pio2 (x, y, int e0, int nx, int prec, ipio2)
 

Variables

static double PIo2 []
 
static double zero = 0.0
 
static double one = 1.0
 
static double two24 = 1.67772160000000000000e+07
 
static double twon24 = 5.96046447753906250000e-08
 

Function Documentation

§ __kernel_rem_pio2()

int attribute_hidden __kernel_rem_pio2 ( x  ,
y  ,
int  e0,
int  nx,
int  prec,
ipio2   
)

Definition at line 176 of file k_rem_pio2.c.

References floor, i, j, k, one, scalbn, SDL_arraysize, SDL_assert, two24, twon24, and zero.

Referenced by __ieee754_rem_pio2().

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  fw += x[j] * f[jx + i - j];
209  q[i] = fw;
210  }
211 
212  jz = jk;
213  recompute:
214  /* distill q[] into iq[] reversingly */
215  for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) {
216  fw = (double) ((int32_t) (twon24 * z));
217  iq[i] = (int32_t) (z - two24 * fw);
218  z = q[j - 1] + fw;
219  }
220 
221  /* compute n */
222  z = scalbn(z, q0); /* actual value of z */
223  z -= 8.0 * floor(z * 0.125); /* trim off integer >= 8 */
224  n = (int32_t) z;
225  z -= (double) n;
226  ih = 0;
227  if (q0 > 0) { /* need iq[jz-1] to determine n */
228  i = (iq[jz - 1] >> (24 - q0));
229  n += i;
230  iq[jz - 1] -= i << (24 - q0);
231  ih = iq[jz - 1] >> (23 - q0);
232  } else if (q0 == 0)
233  ih = iq[jz - 1] >> 23;
234  else if (z >= 0.5)
235  ih = 2;
236 
237  if (ih > 0) { /* q > 0.5 */
238  n += 1;
239  carry = 0;
240  for (i = 0; i < jz; i++) { /* compute 1-q */
241  j = iq[i];
242  if (carry == 0) {
243  if (j != 0) {
244  carry = 1;
245  iq[i] = 0x1000000 - j;
246  }
247  } else
248  iq[i] = 0xffffff - j;
249  }
250  if (q0 > 0) { /* rare case: chance is 1 in 12 */
251  switch (q0) {
252  case 1:
253  iq[jz - 1] &= 0x7fffff;
254  break;
255  case 2:
256  iq[jz - 1] &= 0x3fffff;
257  break;
258  }
259  }
260  if (ih == 2) {
261  z = one - z;
262  if (carry != 0)
263  z -= scalbn(one, q0);
264  }
265  }
266 
267  /* check if recomputation is needed */
268  if (z == zero) {
269  j = 0;
270  for (i = jz - 1; i >= jk; i--)
271  j |= iq[i];
272  if (j == 0) { /* need recomputation */
273  for (k = 1; iq[jk - k] == 0; k++); /* k = no. of terms needed */
274 
275  for (i = jz + 1; i <= jz + k; i++) { /* add q[jz+1] to q[jz+k] */
276  f[jx + i] = (double) ipio2[jv + i];
277  for (j = 0, fw = 0.0; j <= jx; j++)
278  fw += x[j] * f[jx + i - j];
279  q[i] = fw;
280  }
281  jz += k;
282  goto recompute;
283  }
284  }
285 
286  /* chop off zero terms */
287  if (z == 0.0) {
288  jz -= 1;
289  q0 -= 24;
290  while (iq[jz] == 0) {
291  jz--;
292  q0 -= 24;
293  }
294  } else { /* break z into 24-bit if necessary */
295  z = scalbn(z, -q0);
296  if (z >= two24) {
297  fw = (double) ((int32_t) (twon24 * z));
298  iq[jz] = (int32_t) (z - two24 * fw);
299  jz += 1;
300  q0 += 24;
301  iq[jz] = (int32_t) fw;
302  } else
303  iq[jz] = (int32_t) z;
304  }
305 
306  /* convert integer "bit" chunk to floating-point value */
307  fw = scalbn(one, q0);
308  for (i = jz; i >= 0; i--) {
309  q[i] = fw * (double) iq[i];
310  fw *= twon24;
311  }
312 
313  /* compute PIo2[0,...,jp]*q[jz,...,0] */
314  for (i = jz; i >= 0; i--) {
315  for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++)
316  fw += PIo2[k] * q[i + k];
317  fq[jz - i] = fw;
318  }
319 
320  /* compress fq[] into y[] */
321  switch (prec) {
322  case 0:
323  fw = 0.0;
324  for (i = jz; i >= 0; i--)
325  fw += fq[i];
326  y[0] = (ih == 0) ? fw : -fw;
327  break;
328  case 1:
329  case 2:
330  fw = 0.0;
331  for (i = jz; i >= 0; i--)
332  fw += fq[i];
333  y[0] = (ih == 0) ? fw : -fw;
334  fw = fq[0] - fw;
335  for (i = 1; i <= jz; i++)
336  fw += fq[i];
337  y[1] = (ih == 0) ? fw : -fw;
338  break;
339  case 3: /* painful */
340  for (i = jz; i > 0; i--) {
341  fw = fq[i - 1] + fq[i];
342  fq[i] += fq[i - 1] - fw;
343  fq[i - 1] = fw;
344  }
345  for (i = jz; i > 1; i--) {
346  fw = fq[i - 1] + fq[i];
347  fq[i] += fq[i - 1] - fw;
348  fq[i - 1] = fw;
349  }
350  for (fw = 0.0, i = jz; i >= 2; i--)
351  fw += fq[i];
352  if (ih == 0) {
353  y[0] = fq[0];
354  y[1] = fq[1];
355  y[2] = fw;
356  } else {
357  y[0] = -fq[0];
358  y[1] = -fq[1];
359  y[2] = -fw;
360  }
361  }
362  return n & 7;
363 }
static double PIo2[]
Definition: k_rem_pio2.c:150
GLdouble n
signed int int32_t
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1567
GLdouble GLdouble GLdouble GLdouble q
Definition: SDL_opengl.h:2080
const GLfloat * m
GLfloat f
static double one
Definition: k_rem_pio2.c:167
#define scalbn
Definition: math_private.h:40
static double two24
Definition: k_rem_pio2.c:167
GLint GLint GLint GLint GLint GLint y
Definition: SDL_opengl.h:1567
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
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 in i)
Definition: SDL_x11sym.h:50
#define SDL_assert(condition)
Definition: SDL_assert.h:167
GLdouble GLdouble z
static double zero
Definition: k_rem_pio2.c:167
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:211
#define SDL_arraysize(array)
Definition: SDL_stdinc.h:90
#define floor
Definition: math_private.h:37

§ libm_hidden_proto()

libm_hidden_proto ( scalbn  )

Definition at line 139 of file k_rem_pio2.c.

References PIo2.

142  { 2, 3, 4, 6 }; /* initial value for jk */
143 #else
144  static int init_jk[] = { 2, 3, 4, 6 };

Variable Documentation

§ one

double one = 1.0
static

Definition at line 167 of file k_rem_pio2.c.

Referenced by __kernel_rem_pio2().

§ PIo2

double PIo2[]
static
Initial value:
= {
1.57079625129699707031e+00,
7.54978941586159635335e-08,
5.39030252995776476554e-15,
3.28200341580791294123e-22,
1.27065575308067607349e-29,
1.22933308981111328932e-36,
2.73370053816464559624e-44,
2.16741683877804819444e-51,
}

Definition at line 150 of file k_rem_pio2.c.

Referenced by libm_hidden_proto().

§ two24

double two24 = 1.67772160000000000000e+07
static

Definition at line 167 of file k_rem_pio2.c.

Referenced by __kernel_rem_pio2().

§ twon24

double twon24 = 5.96046447753906250000e-08
static

Definition at line 168 of file k_rem_pio2.c.

Referenced by __kernel_rem_pio2().

§ zero

double zero = 0.0
static

Definition at line 167 of file k_rem_pio2.c.

Referenced by __kernel_rem_pio2().