SDL 2.0
math_private.h File Reference
#include "SDL_endian.h"
+ Include dependency graph for math_private.h:
+ This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Data Structures

union  ieee_double_shape_type
 
union  ieee_float_shape_type
 

Macros

#define _IEEE_LIBM
 
#define attribute_hidden
 
#define libm_hidden_proto(x)
 
#define libm_hidden_def(x)
 
#define strong_alias(x, y)
 
#define atan   SDL_uclibc_atan
 
#define __ieee754_atan2   SDL_uclibc_atan2
 
#define copysign   SDL_uclibc_copysign
 
#define cos   SDL_uclibc_cos
 
#define __ieee754_exp   SDL_uclibc_exp
 
#define fabs   SDL_uclibc_fabs
 
#define floor   SDL_uclibc_floor
 
#define __ieee754_fmod   SDL_uclibc_fmod
 
#define __ieee754_log   SDL_uclibc_log
 
#define __ieee754_log10   SDL_uclibc_log10
 
#define __ieee754_pow   SDL_uclibc_pow
 
#define scalbln   SDL_uclibc_scalbln
 
#define scalbn   SDL_uclibc_scalbn
 
#define sin   SDL_uclibc_sin
 
#define __ieee754_sqrt   SDL_uclibc_sqrt
 
#define tan   SDL_uclibc_tan
 
#define EXTRACT_WORDS(ix0, ix1, d)
 
#define GET_HIGH_WORD(i, d)
 
#define GET_LOW_WORD(i, d)
 
#define INSERT_WORDS(d, ix0, ix1)
 
#define SET_HIGH_WORD(d, v)
 
#define SET_LOW_WORD(d, v)
 
#define GET_FLOAT_WORD(i, d)
 
#define SET_FLOAT_WORD(d, i)
 

Typedefs

typedef unsigned int u_int32_t
 

Functions

double __ieee754_sqrt (double) attribute_hidden
 
double __ieee754_acos (double) attribute_hidden
 
double __ieee754_acosh (double) attribute_hidden
 
double __ieee754_log (double) attribute_hidden
 
double __ieee754_atanh (double) attribute_hidden
 
double __ieee754_asin (double) attribute_hidden
 
double __ieee754_atan2 (double, double) attribute_hidden
 
double __ieee754_exp (double) attribute_hidden
 
double __ieee754_cosh (double) attribute_hidden
 
double __ieee754_fmod (double, double) attribute_hidden
 
double __ieee754_pow (double, double) attribute_hidden
 
double __ieee754_lgamma_r (double, int *) attribute_hidden
 
double __ieee754_gamma_r (double, int *) attribute_hidden
 
double __ieee754_lgamma (double) attribute_hidden
 
double __ieee754_gamma (double) attribute_hidden
 
double __ieee754_log10 (double) attribute_hidden
 
double __ieee754_sinh (double) attribute_hidden
 
double __ieee754_hypot (double, double) attribute_hidden
 
double __ieee754_j0 (double) attribute_hidden
 
double __ieee754_j1 (double) attribute_hidden
 
double __ieee754_y0 (double) attribute_hidden
 
double __ieee754_y1 (double) attribute_hidden
 
double __ieee754_jn (int, double) attribute_hidden
 
double __ieee754_yn (int, double) attribute_hidden
 
double __ieee754_remainder (double, double) attribute_hidden
 
int32_t __ieee754_rem_pio2 (double, double *) attribute_hidden
 
double __ieee754_scalb (double, double) attribute_hidden
 
double __kernel_sin (double, double, int) attribute_hidden
 
double __kernel_cos (double, double) attribute_hidden
 
double __kernel_tan (double, double, int) attribute_hidden
 
int32_t __kernel_rem_pio2 (double *, double *, int, int, const unsigned int, const int32_t *) attribute_hidden
 

Macro Definition Documentation

◆ __ieee754_atan2

#define __ieee754_atan2   SDL_uclibc_atan2

Definition at line 35 of file math_private.h.

◆ __ieee754_exp

#define __ieee754_exp   SDL_uclibc_exp

Definition at line 38 of file math_private.h.

◆ __ieee754_fmod

#define __ieee754_fmod   SDL_uclibc_fmod

Definition at line 41 of file math_private.h.

◆ __ieee754_log

#define __ieee754_log   SDL_uclibc_log

Definition at line 42 of file math_private.h.

◆ __ieee754_log10

#define __ieee754_log10   SDL_uclibc_log10

Definition at line 43 of file math_private.h.

◆ __ieee754_pow

#define __ieee754_pow   SDL_uclibc_pow

Definition at line 44 of file math_private.h.

◆ __ieee754_sqrt

#define __ieee754_sqrt   SDL_uclibc_sqrt

Definition at line 48 of file math_private.h.

◆ _IEEE_LIBM

#define _IEEE_LIBM

Definition at line 24 of file math_private.h.

◆ atan

#define atan   SDL_uclibc_atan

Definition at line 34 of file math_private.h.

◆ attribute_hidden

#define attribute_hidden

Definition at line 25 of file math_private.h.

◆ copysign

#define copysign   SDL_uclibc_copysign

Definition at line 36 of file math_private.h.

◆ cos

#define cos   SDL_uclibc_cos

Definition at line 37 of file math_private.h.

◆ EXTRACT_WORDS

#define EXTRACT_WORDS (   ix0,
  ix1,
  d 
)
Value:
do { \
ieee_double_shape_type ew_u; \
ew_u.value = (d); \
(ix0) = ew_u.parts.msw; \
(ix1) = ew_u.parts.lsw; \
} while (0)
SDL_PRINTF_FORMAT_STRING const char int SDL_PRINTF_FORMAT_STRING const char int SDL_PRINTF_FORMAT_STRING const char int SDL_PRINTF_FORMAT_STRING const char const char SDL_SCANF_FORMAT_STRING const char return SDL_ThreadFunction const char void return Uint32 return Uint32 SDL_AssertionHandler void SDL_SpinLock SDL_atomic_t int int return SDL_atomic_t return void void void return void return int return SDL_AudioSpec SDL_AudioSpec return int int return return int SDL_RWops int SDL_AudioSpec Uint8 ** d

Definition at line 99 of file math_private.h.

◆ fabs

#define fabs   SDL_uclibc_fabs

Definition at line 39 of file math_private.h.

◆ floor

#define floor   SDL_uclibc_floor

Definition at line 40 of file math_private.h.

◆ GET_FLOAT_WORD

#define GET_FLOAT_WORD (   i,
  d 
)
Value:
do { \
ieee_float_shape_type gf_u; \
gf_u.value = (d); \
(i) = gf_u.word; \
} while (0)
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

Definition at line 166 of file math_private.h.

◆ GET_HIGH_WORD

#define GET_HIGH_WORD (   i,
  d 
)
Value:
do { \
ieee_double_shape_type gh_u; \
gh_u.value = (d); \
(i) = gh_u.parts.msw; \
} while (0)

Definition at line 109 of file math_private.h.

◆ GET_LOW_WORD

#define GET_LOW_WORD (   i,
  d 
)
Value:
do { \
ieee_double_shape_type gl_u; \
gl_u.value = (d); \
(i) = gl_u.parts.lsw; \
} while (0)

Definition at line 118 of file math_private.h.

◆ INSERT_WORDS

#define INSERT_WORDS (   d,
  ix0,
  ix1 
)
Value:
do { \
ieee_double_shape_type iw_u; \
iw_u.parts.msw = (ix0); \
iw_u.parts.lsw = (ix1); \
(d) = iw_u.value; \
} while (0)

Definition at line 127 of file math_private.h.

◆ libm_hidden_def

#define libm_hidden_def (   x)

Definition at line 27 of file math_private.h.

◆ libm_hidden_proto

#define libm_hidden_proto (   x)

Definition at line 26 of file math_private.h.

◆ scalbln

#define scalbln   SDL_uclibc_scalbln

Definition at line 45 of file math_private.h.

◆ scalbn

#define scalbn   SDL_uclibc_scalbn

Definition at line 46 of file math_private.h.

◆ SET_FLOAT_WORD

#define SET_FLOAT_WORD (   d,
  i 
)
Value:
do { \
ieee_float_shape_type sf_u; \
sf_u.word = (i); \
(d) = sf_u.value; \
} while (0)

Definition at line 175 of file math_private.h.

◆ SET_HIGH_WORD

#define SET_HIGH_WORD (   d,
  v 
)
Value:
do { \
ieee_double_shape_type sh_u; \
sh_u.value = (d); \
sh_u.parts.msw = (v); \
(d) = sh_u.value; \
} while (0)
const GLdouble * v
Definition: SDL_opengl.h:2064

Definition at line 137 of file math_private.h.

◆ SET_LOW_WORD

#define SET_LOW_WORD (   d,
  v 
)
Value:
do { \
ieee_double_shape_type sl_u; \
sl_u.value = (d); \
sl_u.parts.lsw = (v); \
(d) = sl_u.value; \
} while (0)

Definition at line 147 of file math_private.h.

◆ sin

#define sin   SDL_uclibc_sin

Definition at line 47 of file math_private.h.

◆ strong_alias

#define strong_alias (   x,
  y 
)

Definition at line 28 of file math_private.h.

◆ tan

#define tan   SDL_uclibc_tan

Definition at line 49 of file math_private.h.

Typedef Documentation

◆ u_int32_t

typedef unsigned int u_int32_t

Definition at line 31 of file math_private.h.

Function Documentation

◆ __ieee754_acos()

double __ieee754_acos ( double  )

◆ __ieee754_acosh()

double __ieee754_acosh ( double  )

◆ __ieee754_asin()

double __ieee754_asin ( double  )

◆ __ieee754_atan2()

double __ieee754_atan2 ( double  y,
double  x 
)

Definition at line 50 of file e_atan2.c.

51{
52 double z;
53 int32_t k,m,hx,hy,ix,iy;
54 u_int32_t lx,ly;
55
56 EXTRACT_WORDS(hx,lx,x);
57 ix = hx&0x7fffffff;
58 EXTRACT_WORDS(hy,ly,y);
59 iy = hy&0x7fffffff;
60 if(((ix|((lx|-(int32_t)lx)>>31))>0x7ff00000)||
61 ((iy|((ly|-(int32_t)ly)>>31))>0x7ff00000)) /* x or y is NaN */
62 return x+y;
63 if(((hx-0x3ff00000)|lx)==0) return atan(y); /* x=1.0 */
64 m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */
65
66 /* when y = 0 */
67 if((iy|ly)==0) {
68 switch(m) {
69 case 0:
70 case 1: return y; /* atan(+-0,+anything)=+-0 */
71 case 2: return pi+tiny;/* atan(+0,-anything) = pi */
72 case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
73 }
74 }
75 /* when x = 0 */
76 if((ix|lx)==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
77
78 /* when x is INF */
79 if(ix==0x7ff00000) {
80 if(iy==0x7ff00000) {
81 switch(m) {
82 case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */
83 case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
84 case 2: return 3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
85 case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
86 }
87 } else {
88 switch(m) {
89 case 0: return zero ; /* atan(+...,+INF) */
90 case 1: return -zero ; /* atan(-...,+INF) */
91 case 2: return pi+tiny ; /* atan(+...,-INF) */
92 case 3: return -pi-tiny ; /* atan(-...,-INF) */
93 }
94 }
95 }
96 /* when y is INF */
97 if(iy==0x7ff00000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
98
99 /* compute y/x */
100 k = (iy-ix)>>20;
101 if(k > 60) z=pi_o_2+0.5*pi_lo; /* |y/x| > 2**60 */
102 else if(hx<0&&k<-60) z=0.0; /* |y|/x < -2**60 */
103 else z=atan(fabs(y/x)); /* safe to do y/x */
104 switch (m) {
105 case 0: return z ; /* atan(+,+) */
106 case 1: {
107 u_int32_t zh;
108 GET_HIGH_WORD(zh,z);
109 SET_HIGH_WORD(z,zh ^ 0x80000000);
110 }
111 return z ; /* atan(-,+) */
112 case 2: return pi-(z-pi_lo);/* atan(+,-) */
113 default: /* case 3 */
114 return (z-pi_lo)-pi;/* atan(-,-) */
115 }
116}
signed int int32_t
GLint GLint GLint GLint GLint GLint y
Definition: SDL_opengl.h:1574
GLint GLint GLint GLint GLint x
Definition: SDL_opengl.h:1574
const GLfloat * m
GLdouble GLdouble z
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 pi_o_4
Definition: e_atan2.c:45
static const double pi
Definition: e_atan2.c:47
static const double tiny
Definition: e_atan2.c:43
static const double zero
Definition: e_atan2.c:44
static const double pi_o_2
Definition: e_atan2.c:46
static const double pi_lo
Definition: e_atan2.c:48
unsigned int u_int32_t
Definition: math_private.h:31
#define SET_HIGH_WORD(d, v)
Definition: math_private.h:137
#define EXTRACT_WORDS(ix0, ix1, d)
Definition: math_private.h:99
#define GET_HIGH_WORD(i, d)
Definition: math_private.h:109
double atan(double x)
Definition: s_atan.c:71
double fabs(double x)
Definition: s_fabs.c:22

References atan(), EXTRACT_WORDS, fabs(), GET_HIGH_WORD, k, pi, pi_lo, pi_o_2, pi_o_4, SET_HIGH_WORD, tiny, and zero.

◆ __ieee754_atanh()

double __ieee754_atanh ( double  )

◆ __ieee754_cosh()

double __ieee754_cosh ( double  )

◆ __ieee754_exp()

double __ieee754_exp ( double  x)

Definition at line 100 of file e_exp.c.

101{
102 double y;
103 double hi = 0.0;
104 double lo = 0.0;
105 double c;
106 double t;
107 int32_t k=0;
108 int32_t xsb;
109 u_int32_t hx;
110
111 GET_HIGH_WORD(hx,x);
112 xsb = (hx>>31)&1; /* sign bit of x */
113 hx &= 0x7fffffff; /* high word of |x| */
114
115 /* filter out non-finite argument */
116 if(hx >= 0x40862E42) { /* if |x|>=709.78... */
117 if(hx>=0x7ff00000) {
118 u_int32_t lx;
119 GET_LOW_WORD(lx,x);
120 if(((hx&0xfffff)|lx)!=0)
121 return x+x; /* NaN */
122 else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */
123 }
124 #if 1
125 if(x > o_threshold) return huge*huge; /* overflow */
126 #else /* !!! FIXME: check this: "huge * huge" is a compiler warning, maybe they wanted +Inf? */
127 if(x > o_threshold) return INFINITY; /* overflow */
128 #endif
129
130 if(x < u_threshold) return twom1000*twom1000; /* underflow */
131 }
132
133 /* argument reduction */
134 if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
135 if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
136 hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
137 } else {
138 k = (int32_t) (invln2*x+halF[xsb]);
139 t = k;
140 hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */
141 lo = t*ln2LO[0];
142 }
143 x = hi - lo;
144 }
145 else if(hx < 0x3e300000) { /* when |x|<2**-28 */
146 if(huge+x>one) return one+x;/* trigger inexact */
147 }
148 else k = 0;
149
150 /* x is now in primary range */
151 t = x*x;
152 c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
153 if(k==0) return one-((x*c)/(c-2.0)-x);
154 else y = one-((lo-(x*c)/(2.0-c))-hi);
155 if(k >= -1021) {
156 u_int32_t hy;
157 GET_HIGH_WORD(hy,y);
158 SET_HIGH_WORD(y,hy+(k<<20)); /* add k to y's exponent */
159 return y;
160 } else {
161 u_int32_t hy;
162 GET_HIGH_WORD(hy,y);
163 SET_HIGH_WORD(y,hy+((k+1000)<<20)); /* add k to y's exponent */
164 return y*twom1000;
165 }
166}
GLdouble GLdouble t
Definition: SDL_opengl.h:2071
const GLubyte * c
static const double halF[2]
Definition: e_exp.c:84
static const double invln2
Definition: e_exp.c:93
static const double P1
Definition: e_exp.c:94
static const double twom1000
Definition: e_exp.c:86
static const double ln2LO[2]
Definition: e_exp.c:91
static const double ln2HI[2]
Definition: e_exp.c:89
static const double u_threshold
Definition: e_exp.c:88
static const double P3
Definition: e_exp.c:96
static const double o_threshold
Definition: e_exp.c:87
static const double huge
Definition: e_exp.c:85
static const double P2
Definition: e_exp.c:95
static const double one
Definition: e_exp.c:83
static const double P4
Definition: e_exp.c:97
static const double P5
Definition: e_exp.c:98
#define GET_LOW_WORD(i, d)
Definition: math_private.h:118

References GET_HIGH_WORD, GET_LOW_WORD, halF, huge, invln2, k, ln2HI, ln2LO, o_threshold, one, P1, P2, P3, P4, P5, SET_HIGH_WORD, twom1000, and u_threshold.

◆ __ieee754_fmod()

double __ieee754_fmod ( double  x,
double  y 
)

Definition at line 23 of file e_fmod.c.

24{
25 int32_t n,hx,hy,hz,ix,iy,sx,i;
26 u_int32_t lx,ly,lz;
27
28 EXTRACT_WORDS(hx,lx,x);
29 EXTRACT_WORDS(hy,ly,y);
30 sx = hx&0x80000000; /* sign of x */
31 hx ^=sx; /* |x| */
32 hy &= 0x7fffffff; /* |y| */
33
34 /* purge off exception values */
35 if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */
36 ((hy|((ly|-(int32_t)ly)>>31))>0x7ff00000)) /* or y is NaN */
37 return (x*y)/(x*y);
38 if(hx<=hy) {
39 if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
40 if(lx==ly)
41 return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/
42 }
43
44 /* determine ix = ilogb(x) */
45 if(hx<0x00100000) { /* subnormal x */
46 if(hx==0) {
47 for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
48 } else {
49 for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
50 }
51 } else ix = (hx>>20)-1023;
52
53 /* determine iy = ilogb(y) */
54 if(hy<0x00100000) { /* subnormal y */
55 if(hy==0) {
56 for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
57 } else {
58 for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
59 }
60 } else iy = (hy>>20)-1023;
61
62 /* set up {hx,lx}, {hy,ly} and align y to x */
63 if(ix >= -1022)
64 hx = 0x00100000|(0x000fffff&hx);
65 else { /* subnormal x, shift x to normal */
66 n = -1022-ix;
67 if(n<=31) {
68 hx = (hx<<n)|(lx>>(32-n));
69 lx <<= n;
70 } else {
71 hx = lx<<(n-32);
72 lx = 0;
73 }
74 }
75 if(iy >= -1022)
76 hy = 0x00100000|(0x000fffff&hy);
77 else { /* subnormal y, shift y to normal */
78 n = -1022-iy;
79 if(n<=31) {
80 hy = (hy<<n)|(ly>>(32-n));
81 ly <<= n;
82 } else {
83 hy = ly<<(n-32);
84 ly = 0;
85 }
86 }
87
88 /* fix point fmod */
89 n = ix - iy;
90 while(n--) {
91 hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
92 if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
93 else {
94 if((hz|lz)==0) /* return sign(x)*0 */
95 return Zero[(u_int32_t)sx>>31];
96 hx = hz+hz+(lz>>31); lx = lz+lz;
97 }
98 }
99 hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
100 if(hz>=0) {hx=hz;lx=lz;}
101
102 /* convert back to floating value and restore the sign */
103 if((hx|lx)==0) /* return sign(x)*0 */
104 return Zero[(u_int32_t)sx>>31];
105 while(hx<0x00100000) { /* normalize x */
106 hx = hx+hx+(lx>>31); lx = lx+lx;
107 iy -= 1;
108 }
109 if(iy>= -1022) { /* normalize output */
110 hx = ((hx-0x00100000)|((iy+1023)<<20));
111 INSERT_WORDS(x,hx|sx,lx);
112 } else { /* subnormal output */
113 n = -1022 - iy;
114 if(n<=20) {
115 lx = (lx>>n)|((u_int32_t)hx<<(32-n));
116 hx >>= n;
117 } else if (n<=31) {
118 lx = (hx<<(32-n))|(lx>>n); hx = sx;
119 } else {
120 lx = hx>>(n-32); hx = sx;
121 }
122 INSERT_WORDS(x,hx|sx,lx);
123 x *= one; /* create necessary signal */
124 }
125 return x; /* exact output */
126}
GLdouble n
static const double Zero[]
Definition: e_fmod.c:21
static const double one
Definition: e_fmod.c:21
#define INSERT_WORDS(d, ix0, ix1)
Definition: math_private.h:127

References EXTRACT_WORDS, i, INSERT_WORDS, one, and Zero.

◆ __ieee754_gamma()

double __ieee754_gamma ( double  )

◆ __ieee754_gamma_r()

double __ieee754_gamma_r ( double  ,
int *   
)

◆ __ieee754_hypot()

double __ieee754_hypot ( double  ,
double   
)

◆ __ieee754_j0()

double __ieee754_j0 ( double  )

◆ __ieee754_j1()

double __ieee754_j1 ( double  )

◆ __ieee754_jn()

double __ieee754_jn ( int  ,
double   
)

◆ __ieee754_lgamma()

double __ieee754_lgamma ( double  )

◆ __ieee754_lgamma_r()

double __ieee754_lgamma_r ( double  ,
int *   
)

◆ __ieee754_log()

double __ieee754_log ( double  x)

Definition at line 85 of file e_log.c.

86{
87 double hfsq,f,s,z,R,w,t1,t2,dk;
88 int32_t k,hx,i,j;
89 u_int32_t lx;
90
91 EXTRACT_WORDS(hx,lx,x);
92
93 k=0;
94 if (hx < 0x00100000) { /* x < 2**-1022 */
95 if (((hx&0x7fffffff)|lx)==0)
96 return -two54/zero; /* log(+-0)=-inf */
97 if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
98 k -= 54; x *= two54; /* subnormal number, scale up x */
99 GET_HIGH_WORD(hx,x);
100 }
101 if (hx >= 0x7ff00000) return x+x;
102 k += (hx>>20)-1023;
103 hx &= 0x000fffff;
104 i = (hx+0x95f64)&0x100000;
105 SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */
106 k += (i>>20);
107 f = x-1.0;
108 if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */
109 if(f==zero) {if(k==0) return zero; else {dk=(double)k;
110 return dk*ln2_hi+dk*ln2_lo;}
111 }
112 R = f*f*(0.5-0.33333333333333333*f);
113 if(k==0) return f-R; else {dk=(double)k;
114 return dk*ln2_hi-((R-dk*ln2_lo)-f);}
115 }
116 s = f/(2.0+f);
117 dk = (double)k;
118 z = s*s;
119 i = hx-0x6147a;
120 w = z*z;
121 j = 0x6b851-hx;
122 t1= w*(Lg2+w*(Lg4+w*Lg6));
123 t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
124 i |= j;
125 R = t2+t1;
126 if(i>0) {
127 hfsq=0.5*f*f;
128 if(k==0) return f-(hfsq-s*(hfsq+R)); else
129 return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
130 } else {
131 if(k==0) return f-s*(f-R); else
132 return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
133 }
134}
GLdouble s
Definition: SDL_opengl.h:2063
GLfloat f
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
GLubyte GLubyte GLubyte GLubyte w
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
static const double Lg7
Definition: e_log.c:81
static const double Lg6
Definition: e_log.c:80
static const double two54
Definition: e_log.c:74
static const double zero
Definition: e_log.c:83
static const double ln2_hi
Definition: e_log.c:72
static const double Lg4
Definition: e_log.c:78
static const double Lg5
Definition: e_log.c:79
static const double ln2_lo
Definition: e_log.c:73
static const double Lg3
Definition: e_log.c:77
static const double Lg2
Definition: e_log.c:76
static const double Lg1
Definition: e_log.c:75

References EXTRACT_WORDS, GET_HIGH_WORD, i, j, k, Lg1, Lg2, Lg3, Lg4, Lg5, Lg6, Lg7, ln2_hi, ln2_lo, SET_HIGH_WORD, two54, and zero.

◆ __ieee754_log10()

double __ieee754_log10 ( double  x)

Definition at line 61 of file e_log10.c.

62{
63 double y,z;
64 int32_t i,k,hx;
65 u_int32_t lx;
66
67 EXTRACT_WORDS(hx,lx,x);
68
69 k=0;
70 if (hx < 0x00100000) { /* x < 2**-1022 */
71 if (((hx&0x7fffffff)|lx)==0)
72 return -two54/zero; /* log(+-0)=-inf */
73 if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
74 k -= 54; x *= two54; /* subnormal number, scale up x */
75 GET_HIGH_WORD(hx,x);
76 }
77 if (hx >= 0x7ff00000) return x+x;
78 k += (hx>>20)-1023;
79 i = ((u_int32_t)k&0x80000000)>>31;
80 hx = (hx&0x000fffff)|((0x3ff-i)<<20);
81 y = (double)(k+i);
82 SET_HIGH_WORD(x,hx);
84 return z+y*log10_2hi;
85}
static const double two54
Definition: e_log10.c:54
static const double log10_2lo
Definition: e_log10.c:57
static const double zero
Definition: e_log10.c:59
static const double ivln10
Definition: e_log10.c:55
static const double log10_2hi
Definition: e_log10.c:56
#define __ieee754_log
Definition: math_private.h:42

References __ieee754_log, EXTRACT_WORDS, GET_HIGH_WORD, i, ivln10, k, log10_2hi, log10_2lo, SET_HIGH_WORD, two54, and zero.

◆ __ieee754_pow()

double __ieee754_pow ( double  x,
double  y 
)

Definition at line 103 of file e_pow.c.

104{
105 double z,ax,z_h,z_l,p_h,p_l;
106 double y1,t1,t2,r,s,t,u,v,w;
107 int32_t i,j,k,yisint,n;
108 int32_t hx,hy,ix,iy;
109 u_int32_t lx,ly;
110
111 EXTRACT_WORDS(hx,lx,x);
112 /* x==1: 1**y = 1 (even if y is NaN) */
113 if (hx==0x3ff00000 && lx==0) {
114 return x;
115 }
116 ix = hx&0x7fffffff;
117
118 EXTRACT_WORDS(hy,ly,y);
119 iy = hy&0x7fffffff;
120
121 /* y==zero: x**0 = 1 */
122 if((iy|ly)==0) return one;
123
124 /* +-NaN return x+y */
125 if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
126 iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
127 return x+y;
128
129 /* determine if y is an odd int when x < 0
130 * yisint = 0 ... y is not an integer
131 * yisint = 1 ... y is an odd int
132 * yisint = 2 ... y is an even int
133 */
134 yisint = 0;
135 if(hx<0) {
136 if(iy>=0x43400000) yisint = 2; /* even integer y */
137 else if(iy>=0x3ff00000) {
138 k = (iy>>20)-0x3ff; /* exponent */
139 if(k>20) {
140 j = ly>>(52-k);
141 if((j<<(52-k))==ly) yisint = 2-(j&1);
142 } else if(ly==0) {
143 j = iy>>(20-k);
144 if((j<<(20-k))==iy) yisint = 2-(j&1);
145 }
146 }
147 }
148
149 /* special value of y */
150 if(ly==0) {
151 if (iy==0x7ff00000) { /* y is +-inf */
152 if (((ix-0x3ff00000)|lx)==0)
153 return one; /* +-1**+-inf is 1 (yes, weird rule) */
154 if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
155 return (hy>=0) ? y : zero;
156 /* (|x|<1)**-,+inf = inf,0 */
157 return (hy<0) ? -y : zero;
158 }
159 if(iy==0x3ff00000) { /* y is +-1 */
160 if(hy<0) return one/x; else return x;
161 }
162 if(hy==0x40000000) return x*x; /* y is 2 */
163 if(hy==0x3fe00000) { /* y is 0.5 */
164 if(hx>=0) /* x >= +0 */
165 return __ieee754_sqrt(x);
166 }
167 }
168
169 ax = fabs(x);
170 /* special value of x */
171 if(lx==0) {
172 if(ix==0x7ff00000||ix==0||ix==0x3ff00000){
173 z = ax; /*x is +-0,+-inf,+-1*/
174 if(hy<0) z = one/z; /* z = (1/|x|) */
175 if(hx<0) {
176 if(((ix-0x3ff00000)|yisint)==0) {
177 z = (z-z)/(z-z); /* (-1)**non-int is NaN */
178 } else if(yisint==1)
179 z = -z; /* (x<0)**odd = -(|x|**odd) */
180 }
181 return z;
182 }
183 }
184
185 /* (x<0)**(non-int) is NaN */
186 if(((((u_int32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x);
187
188 /* |y| is huge */
189 if(iy>0x41e00000) { /* if |y| > 2**31 */
190 if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */
191 if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
192 if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
193 }
194 /* over/underflow if x is not close to one */
195 if(ix<0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
196 if(ix>0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
197 /* now |1-x| is tiny <= 2**-20, suffice to compute
198 log(x) by x-x^2/2+x^3/3-x^4/4 */
199 t = x-1; /* t has 20 trailing zeros */
200 w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
201 u = ivln2_h*t; /* ivln2_h has 21 sig. bits */
202 v = t*ivln2_l-w*ivln2;
203 t1 = u+v;
204 SET_LOW_WORD(t1,0);
205 t2 = v-(t1-u);
206 } else {
207 double s2,s_h,s_l,t_h,t_l;
208 n = 0;
209 /* take care subnormal number */
210 if(ix<0x00100000)
211 {ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); }
212 n += ((ix)>>20)-0x3ff;
213 j = ix&0x000fffff;
214 /* determine interval */
215 ix = j|0x3ff00000; /* normalize ix */
216 if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */
217 else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */
218 else {k=0;n+=1;ix -= 0x00100000;}
219 SET_HIGH_WORD(ax,ix);
220
221 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
222 u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
223 v = one/(ax+bp[k]);
224 s = u*v;
225 s_h = s;
226 SET_LOW_WORD(s_h,0);
227 /* t_h=ax+bp[k] High */
228 t_h = zero;
229 SET_HIGH_WORD(t_h,((ix>>1)|0x20000000)+0x00080000+(k<<18));
230 t_l = ax - (t_h-bp[k]);
231 s_l = v*((u-s_h*t_h)-s_h*t_l);
232 /* compute log(ax) */
233 s2 = s*s;
234 r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
235 r += s_l*(s_h+s);
236 s2 = s_h*s_h;
237 t_h = 3.0+s2+r;
238 SET_LOW_WORD(t_h,0);
239 t_l = r-((t_h-3.0)-s2);
240 /* u+v = s*(1+...) */
241 u = s_h*t_h;
242 v = s_l*t_h+t_l*s;
243 /* 2/(3log2)*(s+...) */
244 p_h = u+v;
245 SET_LOW_WORD(p_h,0);
246 p_l = v-(p_h-u);
247 z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */
248 z_l = cp_l*p_h+p_l*cp+dp_l[k];
249 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
250 t = (double)n;
251 t1 = (((z_h+z_l)+dp_h[k])+t);
252 SET_LOW_WORD(t1,0);
253 t2 = z_l-(((t1-t)-dp_h[k])-z_h);
254 }
255
256 s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
257 if(((((u_int32_t)hx>>31)-1)|(yisint-1))==0)
258 s = -one;/* (-ve)**(odd int) */
259
260 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
261 y1 = y;
262 SET_LOW_WORD(y1,0);
263 p_l = (y-y1)*t1+y*t2;
264 p_h = y1*t1;
265 z = p_l+p_h;
267 if (j>=0x40900000) { /* z >= 1024 */
268 if(((j-0x40900000)|i)!=0) /* if z > 1024 */
269 return s*huge*huge; /* overflow */
270 else {
271 if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */
272 }
273 } else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */
274 if(((j-0xc090cc00)|i)!=0) /* z < -1075 */
275 return s*tiny*tiny; /* underflow */
276 else {
277 if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */
278 }
279 }
280 /*
281 * compute 2**(p_h+p_l)
282 */
283 i = j&0x7fffffff;
284 k = (i>>20)-0x3ff;
285 n = 0;
286 if(i>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
287 n = j+(0x00100000>>(k+1));
288 k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */
289 t = zero;
290 SET_HIGH_WORD(t,n&~(0x000fffff>>k));
291 n = ((n&0x000fffff)|0x00100000)>>(20-k);
292 if(j<0) n = -n;
293 p_h -= t;
294 }
295 t = p_l+p_h;
296 SET_LOW_WORD(t,0);
297 u = t*lg2_h;
298 v = (p_l-(t-p_h))*lg2+t*lg2_l;
299 z = u+v;
300 w = v-(z-u);
301 t = z*z;
302 t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
303 r = (z*t1)/(t1-two)-(w+z*w);
304 z = one-(r-z);
306 j += (n<<20);
307 if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */
308 else SET_HIGH_WORD(z,j);
309 return s*z;
310}
GLdouble GLdouble GLdouble r
Definition: SDL_opengl.h:2079
GLfixed y1
static const double lg2_h
Definition: e_pow.c:93
static const double cp_h
Definition: e_pow.c:97
static const double L2
Definition: e_pow.c:82
static const double L3
Definition: e_pow.c:83
static const double ovt
Definition: e_pow.c:95
static const double cp
Definition: e_pow.c:96
static const double ivln2
Definition: e_pow.c:99
static const double two
Definition: e_pow.c:76
static const double dp_h[]
Definition: e_pow.c:72
static const double tiny
Definition: e_pow.c:79
static const double zero
Definition: e_pow.c:74
static const double P1
Definition: e_pow.c:87
static const double L4
Definition: e_pow.c:84
static const double ivln2_h
Definition: e_pow.c:100
static const double two53
Definition: e_pow.c:77
static const double L1
Definition: e_pow.c:81
static const double bp[]
Definition: e_pow.c:71
static const double L6
Definition: e_pow.c:86
static const double L5
Definition: e_pow.c:85
static const double P3
Definition: e_pow.c:89
static const double cp_l
Definition: e_pow.c:98
static const double lg2_l
Definition: e_pow.c:94
static const double ivln2_l
Definition: e_pow.c:101
static const double huge
Definition: e_pow.c:78
static const double P2
Definition: e_pow.c:88
static const double dp_l[]
Definition: e_pow.c:73
static const double one
Definition: e_pow.c:75
static const double P4
Definition: e_pow.c:90
static const double lg2
Definition: e_pow.c:92
static const double P5
Definition: e_pow.c:91
#define scalbn
Definition: math_private.h:46
#define __ieee754_sqrt
Definition: math_private.h:48
#define SET_LOW_WORD(d, v)
Definition: math_private.h:147

References __ieee754_sqrt, bp, cp, cp_h, cp_l, dp_h, dp_l, EXTRACT_WORDS, fabs(), GET_HIGH_WORD, huge, i, ivln2, ivln2_h, ivln2_l, j, k, L1, L2, L3, L4, L5, L6, lg2, lg2_h, lg2_l, one, ovt, P1, P2, P3, P4, P5, scalbn, SET_HIGH_WORD, SET_LOW_WORD, tiny, two, two53, and zero.

◆ __ieee754_rem_pio2()

int32_t __ieee754_rem_pio2 ( double  x,
double *  y 
)

Definition at line 69 of file e_rem_pio2.c.

70{
71 double z=0.0,w,t,r,fn;
72 double tx[3];
73 int32_t e0,i,j,nx,n,ix,hx;
74 u_int32_t low;
75
76 GET_HIGH_WORD(hx,x); /* high word of x */
77 ix = hx&0x7fffffff;
78 if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */
79 {y[0] = x; y[1] = 0; return 0;}
80 if(ix<0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */
81 if(hx>0) {
82 z = x - pio2_1;
83 if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */
84 y[0] = z - pio2_1t;
85 y[1] = (z-y[0])-pio2_1t;
86 } else { /* near pi/2, use 33+33+53 bit pi */
87 z -= pio2_2;
88 y[0] = z - pio2_2t;
89 y[1] = (z-y[0])-pio2_2t;
90 }
91 return 1;
92 } else { /* negative x */
93 z = x + pio2_1;
94 if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */
95 y[0] = z + pio2_1t;
96 y[1] = (z-y[0])+pio2_1t;
97 } else { /* near pi/2, use 33+33+53 bit pi */
98 z += pio2_2;
99 y[0] = z + pio2_2t;
100 y[1] = (z-y[0])+pio2_2t;
101 }
102 return -1;
103 }
104 }
105 if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
106 t = fabs(x);
107 n = (int32_t) (t*invpio2+half);
108 fn = (double)n;
109 r = t-fn*pio2_1;
110 w = fn*pio2_1t; /* 1st round good to 85 bit */
111 if(n<32&&ix!=npio2_hw[n-1]) {
112 y[0] = r-w; /* quick check no cancellation */
113 } else {
114 u_int32_t high;
115 j = ix>>20;
116 y[0] = r-w;
117 GET_HIGH_WORD(high,y[0]);
118 i = j-((high>>20)&0x7ff);
119 if(i>16) { /* 2nd iteration needed, good to 118 */
120 t = r;
121 w = fn*pio2_2;
122 r = t-w;
123 w = fn*pio2_2t-((t-r)-w);
124 y[0] = r-w;
125 GET_HIGH_WORD(high,y[0]);
126 i = j-((high>>20)&0x7ff);
127 if(i>49) { /* 3rd iteration need, 151 bits acc */
128 t = r; /* will cover all possible cases */
129 w = fn*pio2_3;
130 r = t-w;
131 w = fn*pio2_3t-((t-r)-w);
132 y[0] = r-w;
133 }
134 }
135 }
136 y[1] = (r-y[0])-w;
137 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
138 else return n;
139 }
140 /*
141 * all other (large) arguments
142 */
143 if(ix>=0x7ff00000) { /* x is inf or NaN */
144 y[0]=y[1]=x-x; return 0;
145 }
146 /* set z = scalbn(|x|,ilogb(x)-23) */
147 GET_LOW_WORD(low,x);
148 SET_LOW_WORD(z,low);
149 e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */
150 SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20)));
151 for(i=0;i<2;i++) {
152 tx[i] = (double)((int32_t)(z));
153 z = (z-tx[i])*two24;
154 }
155 tx[2] = z;
156 nx = 3;
157 while((nx > 0) && tx[nx-1]==zero) nx--; /* skip zero term */
159 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
160 return n;
161}
GLbyte nx
static const int32_t two_over_pi[]
Definition: e_rem_pio2.c:24
static const double invpio2
Definition: e_rem_pio2.c:61
static const double pio2_1t
Definition: e_rem_pio2.c:63
static const double zero
Definition: e_rem_pio2.c:58
static const double two24
Definition: e_rem_pio2.c:60
static const double pio2_2t
Definition: e_rem_pio2.c:65
static const double half
Definition: e_rem_pio2.c:59
static const double pio2_3t
Definition: e_rem_pio2.c:67
static const int32_t npio2_hw[]
Definition: e_rem_pio2.c:38
static const double pio2_1
Definition: e_rem_pio2.c:62
static const double pio2_3
Definition: e_rem_pio2.c:66
static const double pio2_2
Definition: e_rem_pio2.c:64
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

References __kernel_rem_pio2(), fabs(), GET_HIGH_WORD, GET_LOW_WORD, half, i, invpio2, j, npio2_hw, pio2_1, pio2_1t, pio2_2, pio2_2t, pio2_3, pio2_3t, SET_HIGH_WORD, SET_LOW_WORD, two24, two_over_pi, and zero.

Referenced by cos(), sin(), and tan().

◆ __ieee754_remainder()

double __ieee754_remainder ( double  ,
double   
)

◆ __ieee754_scalb()

double __ieee754_scalb ( double  ,
double   
)

◆ __ieee754_sinh()

double __ieee754_sinh ( double  )

◆ __ieee754_sqrt()

double __ieee754_sqrt ( double  x)

Definition at line 87 of file e_sqrt.c.

88{
89 double z;
90 int32_t sign = (int)0x80000000;
91 int32_t ix0,s0,q,m,t,i;
92 u_int32_t r,t1,s1,ix1,q1;
93
94 EXTRACT_WORDS(ix0,ix1,x);
95
96 /* take care of Inf and NaN */
97 if((ix0&0x7ff00000)==0x7ff00000) {
98 return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
99 sqrt(-inf)=sNaN */
100 }
101 /* take care of zero */
102 if(ix0<=0) {
103 if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
104 else if(ix0<0)
105 return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
106 }
107 /* normalize x */
108 m = (ix0>>20);
109 if(m==0) { /* subnormal x */
110 while(ix0==0) {
111 m -= 21;
112 ix0 |= (ix1>>11); ix1 <<= 21;
113 }
114 for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
115 m -= i-1;
116 ix0 |= (ix1>>(32-i));
117 ix1 <<= i;
118 }
119 m -= 1023; /* unbias exponent */
120 ix0 = (ix0&0x000fffff)|0x00100000;
121 if(m&1){ /* odd m, double x to make it even */
122 ix0 += ix0 + ((ix1&sign)>>31);
123 ix1 += ix1;
124 }
125 m >>= 1; /* m = [m/2] */
126
127 /* generate sqrt(x) bit by bit */
128 ix0 += ix0 + ((ix1&sign)>>31);
129 ix1 += ix1;
130 q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
131 r = 0x00200000; /* r = moving bit from right to left */
132
133 while(r!=0) {
134 t = s0+r;
135 if(t<=ix0) {
136 s0 = t+r;
137 ix0 -= t;
138 q += r;
139 }
140 ix0 += ix0 + ((ix1&sign)>>31);
141 ix1 += ix1;
142 r>>=1;
143 }
144
145 r = sign;
146 while(r!=0) {
147 t1 = s1+r;
148 t = s0;
149 if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
150 s1 = t1+r;
151 if(((t1&sign)==sign)&&(s1&sign)==0) s0 += 1;
152 ix0 -= t;
153 if (ix1 < t1) ix0 -= 1;
154 ix1 -= t1;
155 q1 += r;
156 }
157 ix0 += ix0 + ((ix1&sign)>>31);
158 ix1 += ix1;
159 r>>=1;
160 }
161
162 /* use floating add to find out rounding direction */
163 if((ix0|ix1)!=0) {
164 z = one-tiny; /* trigger inexact flag */
165 if (z>=one) {
166 z = one+tiny;
167 if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;}
168 else if (z>one) {
169 if (q1==(u_int32_t)0xfffffffe) q+=1;
170 q1+=2;
171 } else
172 q1 += (q1&1);
173 }
174 }
175 ix0 = (q>>1)+0x3fe00000;
176 ix1 = q1>>1;
177 if ((q&1)==1) ix1 |= sign;
178 ix0 += (m <<20);
179 INSERT_WORDS(z,ix0,ix1);
180 return z;
181}
GLdouble GLdouble GLdouble GLdouble q
Definition: SDL_opengl.h:2087
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat s1
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat s0
static const double tiny
Definition: e_sqrt.c:85
static const double one
Definition: e_sqrt.c:85

References EXTRACT_WORDS, i, INSERT_WORDS, one, and tiny.

◆ __ieee754_y0()

double __ieee754_y0 ( double  )

◆ __ieee754_y1()

double __ieee754_y1 ( double  )

◆ __ieee754_yn()

double __ieee754_yn ( int  ,
double   
)

◆ __kernel_cos()

double __kernel_cos ( double  x,
double  y 
)

Definition at line 59 of file k_cos.c.

60{
61 double a,hz,z,r,qx;
62 int32_t ix;
63 GET_HIGH_WORD(ix,x);
64 ix &= 0x7fffffff; /* ix = |x|'s high word*/
65 if(ix<0x3e400000) { /* if x < 2**27 */
66 if(((int)x)==0) return one; /* generate inexact */
67 }
68 z = x*x;
69 r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))));
70 if(ix < 0x3FD33333) /* if |x| < 0.3 */
71 return one - (0.5*z - (z*r - x*y));
72 else {
73 if(ix > 0x3fe90000) { /* x > 0.78125 */
74 qx = 0.28125;
75 } else {
76 INSERT_WORDS(qx,ix-0x00200000,0); /* x/4 */
77 }
78 hz = 0.5*z-qx;
79 a = one-qx;
80 return a - (hz - (z*r-x*y));
81 }
82}
GLboolean GLboolean GLboolean GLboolean a
static const double C5
Definition: k_cos.c:56
static const double C2
Definition: k_cos.c:53
static const double C1
Definition: k_cos.c:52
static const double C6
Definition: k_cos.c:57
static const double C3
Definition: k_cos.c:54
static const double one
Definition: k_cos.c:51
static const double C4
Definition: k_cos.c:55

References C1, C2, C3, C4, C5, C6, GET_HIGH_WORD, INSERT_WORDS, and one.

Referenced by cos(), and sin().

◆ __kernel_rem_pio2()

int32_t __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
#define SDL_memset
#define SDL_arraysize(array)
Definition: SDL_stdinc.h:115
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
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().

◆ __kernel_sin()

double __kernel_sin ( double  x,
double  y,
int  iy 
)

Definition at line 52 of file k_sin.c.

53{
54 double z,r,v;
55 int32_t ix;
56 GET_HIGH_WORD(ix,x);
57 ix &= 0x7fffffff; /* high word of x */
58 if(ix<0x3e400000) /* |x| < 2**-27 */
59 {if((int)x==0) return x;} /* generate inexact */
60 z = x*x;
61 v = z*x;
62 r = S2+z*(S3+z*(S4+z*(S5+z*S6)));
63 if(iy==0) return x+v*(S1+z*r);
64 else return x-((z*(half*y-v*r)-y)-v*S1);
65}
static const double S3
Definition: k_sin.c:47
static const double S6
Definition: k_sin.c:50
static const double S4
Definition: k_sin.c:48
static const double half
Definition: k_sin.c:44
static const double S5
Definition: k_sin.c:49
static const double S1
Definition: k_sin.c:45
static const double S2
Definition: k_sin.c:46

References GET_HIGH_WORD, half, S1, S2, S3, S4, S5, and S6.

Referenced by cos(), and sin().

◆ __kernel_tan()

double __kernel_tan ( double  x,
double  y,
int  iy 
)

Definition at line 69 of file k_tan.c.

70{
71 double z,r,v,w,s;
72 int32_t ix,hx;
73 GET_HIGH_WORD(hx,x);
74 ix = hx&0x7fffffff; /* high word of |x| */
75 if(ix<0x3e300000) /* x < 2**-28 */
76 {if((int)x==0) { /* generate inexact */
77 u_int32_t low;
78 GET_LOW_WORD(low,x);
79 if(((ix|low)|(iy+1))==0) return one/fabs(x);
80 else return (iy==1)? x: -one/x;
81 }
82 }
83 if(ix>=0x3FE59428) { /* |x|>=0.6744 */
84 if(hx<0) {x = -x; y = -y;}
85 z = pio4-x;
86 w = pio4lo-y;
87 x = z+w; y = 0.0;
88 }
89 z = x*x;
90 w = z*z;
91 /* Break x^5*(T[1]+x^2*T[2]+...) into
92 * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
93 * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
94 */
95 r = T[1]+w*(T[3]+w*(T[5]+w*(T[7]+w*(T[9]+w*T[11]))));
96 v = z*(T[2]+w*(T[4]+w*(T[6]+w*(T[8]+w*(T[10]+w*T[12])))));
97 s = z*x;
98 r = y + z*(s*(r+v)+y);
99 r += T[0]*s;
100 w = x+r;
101 if(ix>=0x3FE59428) {
102 v = (double)iy;
103 return (double)(1-((hx>>30)&2))*(v-2.0*(x-(w*w/(w+v)-r)));
104 }
105 if(iy==1) return w;
106 else { /* if allow error up to 2 ulp,
107 simply return -1.0/(x+r) here */
108 /* compute -1.0/(x+r) accurately */
109 double a,t;
110 z = w;
111 SET_LOW_WORD(z,0);
112 v = r-(z - x); /* z+v = r+x */
113 t = a = -1.0/w; /* a = -1.0/w */
114 SET_LOW_WORD(t,0);
115 s = 1.0+t*z;
116 return t+a*(s+t*v);
117 }
118}
static const double T[]
Definition: k_tan.c:53
static const double pio4lo
Definition: k_tan.c:52
static const double pio4
Definition: k_tan.c:51
static const double one
Definition: k_tan.c:50

References fabs(), GET_HIGH_WORD, GET_LOW_WORD, one, pio4, pio4lo, SET_LOW_WORD, and T.

Referenced by tan().