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

int32_t attribute_hidden __kernel_rem_pio2 (double *x, double *y, int e0, int nx, const unsigned int prec, const int32_t *ipio2)
 

Variables

static const int init_jk [] = {2,3,4,6}
 
static const double PIo2 []
 
static const double zero = 0.0
 
static const double one = 1.0
 
static const double two24 = 1.67772160000000000000e+07
 
static const double twon24 = 5.96046447753906250000e-08
 

Function Documentation

◆ __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 at line 152 of file k_rem_pio2.c.

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 }

References floor(), i, init_jk, j, k, one, PIo2, scalbn, SDL_arraysize, SDL_assert, SDL_memset, two24, twon24, and zero.

Referenced by __ieee754_rem_pio2().

Variable Documentation

◆ init_jk

const int init_jk[] = {2,3,4,6}
static

Definition at line 133 of file k_rem_pio2.c.

Referenced by __kernel_rem_pio2().

◆ one

const double one = 1.0
static

Definition at line 148 of file k_rem_pio2.c.

Referenced by __kernel_rem_pio2().

◆ PIo2

const 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 135 of file k_rem_pio2.c.

Referenced by __kernel_rem_pio2().

◆ two24

const double two24 = 1.67772160000000000000e+07
static

Definition at line 149 of file k_rem_pio2.c.

Referenced by __kernel_rem_pio2().

◆ twon24

const double twon24 = 5.96046447753906250000e-08
static

Definition at line 150 of file k_rem_pio2.c.

Referenced by __kernel_rem_pio2().

◆ zero

const double zero = 0.0
static

Definition at line 147 of file k_rem_pio2.c.

Referenced by __kernel_rem_pio2().

SDL_memset
#define SDL_memset
Definition: SDL_dynapi_overrides.h:386
one
static const double one
Definition: k_rem_pio2.c:148
zero
static const double zero
Definition: k_rem_pio2.c:147
q
GLdouble GLdouble GLdouble GLdouble q
Definition: SDL_opengl.h:2087
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
x
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1574
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
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:1574
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:29
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