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;
186recompute:
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}
#define SDL_assert(condition)
Definition: SDL_assert.h:169
signed int int32_t
#define SDL_memset
GLint GLint GLint GLint GLint GLint y
Definition: SDL_opengl.h:1574
GLdouble GLdouble GLdouble GLdouble q
Definition: SDL_opengl.h:2087
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1574
const GLfloat * m
GLfloat f
GLbyte nx
GLdouble GLdouble z
GLdouble n
#define SDL_arraysize(array)
Definition: SDL_stdinc.h:115
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
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
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
static const double PIo2[]
Definition: k_rem_pio2.c:135
static const double zero
Definition: k_rem_pio2.c:147
static const double two24
Definition: k_rem_pio2.c:149
static const int init_jk[]
Definition: k_rem_pio2.c:133
static const double one
Definition: k_rem_pio2.c:148
static const double twon24
Definition: k_rem_pio2.c:150
#define scalbn
Definition: math_private.h:46
double floor(double x)
Definition: s_floor.c:33

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().