summaryrefslogtreecommitdiff
path: root/contrib/SDL-3.2.8/src/libm
diff options
context:
space:
mode:
author3gg <3gg@shellblade.net>2025-12-27 12:03:39 -0800
committer3gg <3gg@shellblade.net>2025-12-27 12:03:39 -0800
commit5a079a2d114f96d4847d1ee305d5b7c16eeec50e (patch)
tree8926ab44f168acf787d8e19608857b3af0f82758 /contrib/SDL-3.2.8/src/libm
Initial commit
Diffstat (limited to 'contrib/SDL-3.2.8/src/libm')
-rw-r--r--contrib/SDL-3.2.8/src/libm/e_atan2.c135
-rw-r--r--contrib/SDL-3.2.8/src/libm/e_exp.c201
-rw-r--r--contrib/SDL-3.2.8/src/libm/e_fmod.c145
-rw-r--r--contrib/SDL-3.2.8/src/libm/e_log.c153
-rw-r--r--contrib/SDL-3.2.8/src/libm/e_log10.c107
-rw-r--r--contrib/SDL-3.2.8/src/libm/e_pow.c348
-rw-r--r--contrib/SDL-3.2.8/src/libm/e_rem_pio2.c162
-rw-r--r--contrib/SDL-3.2.8/src/libm/e_sqrt.c458
-rw-r--r--contrib/SDL-3.2.8/src/libm/k_cos.c83
-rw-r--r--contrib/SDL-3.2.8/src/libm/k_rem_pio2.c315
-rw-r--r--contrib/SDL-3.2.8/src/libm/k_sin.c66
-rw-r--r--contrib/SDL-3.2.8/src/libm/k_tan.c119
-rw-r--r--contrib/SDL-3.2.8/src/libm/math_libm.h50
-rw-r--r--contrib/SDL-3.2.8/src/libm/math_private.h234
-rw-r--r--contrib/SDL-3.2.8/src/libm/s_atan.c119
-rw-r--r--contrib/SDL-3.2.8/src/libm/s_copysign.c30
-rw-r--r--contrib/SDL-3.2.8/src/libm/s_cos.c74
-rw-r--r--contrib/SDL-3.2.8/src/libm/s_fabs.c30
-rw-r--r--contrib/SDL-3.2.8/src/libm/s_floor.c76
-rw-r--r--contrib/SDL-3.2.8/src/libm/s_isinf.c24
-rw-r--r--contrib/SDL-3.2.8/src/libm/s_isinff.c24
-rw-r--r--contrib/SDL-3.2.8/src/libm/s_isnan.c31
-rw-r--r--contrib/SDL-3.2.8/src/libm/s_isnanf.c33
-rw-r--r--contrib/SDL-3.2.8/src/libm/s_modf.c68
-rw-r--r--contrib/SDL-3.2.8/src/libm/s_scalbn.c74
-rw-r--r--contrib/SDL-3.2.8/src/libm/s_sin.c74
-rw-r--r--contrib/SDL-3.2.8/src/libm/s_tan.c68
27 files changed, 3301 insertions, 0 deletions
diff --git a/contrib/SDL-3.2.8/src/libm/e_atan2.c b/contrib/SDL-3.2.8/src/libm/e_atan2.c
new file mode 100644
index 0000000..0f5b2c0
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/e_atan2.c
@@ -0,0 +1,135 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/* __ieee754_atan2(y,x)
14 * Method :
15 * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
16 * 2. Reduce x to positive by (if x and y are unexceptional):
17 * ARG (x+iy) = arctan(y/x) ... if x > 0,
18 * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
19 *
20 * Special cases:
21 *
22 * ATAN2((anything), NaN ) is NaN;
23 * ATAN2(NAN , (anything) ) is NaN;
24 * ATAN2(+-0, +(anything but NaN)) is +-0 ;
25 * ATAN2(+-0, -(anything but NaN)) is +-pi ;
26 * ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
27 * ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
28 * ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
29 * ATAN2(+-INF,+INF ) is +-pi/4 ;
30 * ATAN2(+-INF,-INF ) is +-3pi/4;
31 * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
32 *
33 * Constants:
34 * The hexadecimal values are the intended ones for the following
35 * constants. The decimal values may be used, provided that the
36 * compiler will convert from decimal to binary accurately enough
37 * to produce the hexadecimal values shown.
38 */
39
40#include "math_libm.h"
41#include "math_private.h"
42
43static const double
44tiny = 1.0e-300,
45zero = 0.0,
46pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
47pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
48pi = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */
49pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
50
51double attribute_hidden __ieee754_atan2(double y, double x)
52{
53 double z;
54 int32_t k,m,hx,hy,ix,iy;
55 u_int32_t lx,ly;
56
57 EXTRACT_WORDS(hx,lx,x);
58 ix = hx&0x7fffffff;
59 EXTRACT_WORDS(hy,ly,y);
60 iy = hy&0x7fffffff;
61 if(((ix|((lx|-(int32_t)lx)>>31))>0x7ff00000)||
62 ((iy|((ly|-(int32_t)ly)>>31))>0x7ff00000)) /* x or y is NaN */
63 return x+y;
64 if(((hx-0x3ff00000)|lx)==0) return atan(y); /* x=1.0 */
65 m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */
66
67 /* when y = 0 */
68 if((iy|ly)==0) {
69 switch(m) {
70 case 0:
71 case 1: return y; /* atan(+-0,+anything)=+-0 */
72 case 2: return pi+tiny;/* atan(+0,-anything) = pi */
73 case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */
74 }
75 }
76 /* when x = 0 */
77 if((ix|lx)==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
78
79 /* when x is INF */
80 if(ix==0x7ff00000) {
81 if(iy==0x7ff00000) {
82 switch(m) {
83 case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */
84 case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */
85 case 2: return 3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/
86 case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/
87 }
88 } else {
89 switch(m) {
90 case 0: return zero ; /* atan(+...,+INF) */
91 case 1: return -zero ; /* atan(-...,+INF) */
92 case 2: return pi+tiny ; /* atan(+...,-INF) */
93 case 3: return -pi-tiny ; /* atan(-...,-INF) */
94 }
95 }
96 }
97 /* when y is INF */
98 if(iy==0x7ff00000) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny;
99
100 /* compute y/x */
101 k = (iy-ix)>>20;
102 if(k > 60) z=pi_o_2+0.5*pi_lo; /* |y/x| > 2**60 */
103 else if(hx<0&&k<-60) z=0.0; /* |y|/x < -2**60 */
104 else z=atan(fabs(y/x)); /* safe to do y/x */
105 switch (m) {
106 case 0: return z ; /* atan(+,+) */
107 case 1: {
108 u_int32_t zh;
109 GET_HIGH_WORD(zh,z);
110 SET_HIGH_WORD(z,zh ^ 0x80000000);
111 }
112 return z ; /* atan(-,+) */
113 case 2: return pi-(z-pi_lo);/* atan(+,-) */
114 default: /* case 3 */
115 return (z-pi_lo)-pi;/* atan(-,-) */
116 }
117}
118
119/*
120 * wrapper atan2(y,x)
121 */
122#ifndef _IEEE_LIBM
123double atan2(double y, double x)
124{
125 double z = __ieee754_atan2(y, x);
126 if (_LIB_VERSION == _IEEE_ || isnan(x) || isnan(y))
127 return z;
128 if (x == 0.0 && y == 0.0)
129 return __kernel_standard(y,x,3); /* atan2(+-0,+-0) */
130 return z;
131}
132#else
133strong_alias(__ieee754_atan2, atan2)
134#endif
135libm_hidden_def(atan2)
diff --git a/contrib/SDL-3.2.8/src/libm/e_exp.c b/contrib/SDL-3.2.8/src/libm/e_exp.c
new file mode 100644
index 0000000..f39bb5c
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/e_exp.c
@@ -0,0 +1,201 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/* __ieee754_exp(x)
14 * Returns the exponential of x.
15 *
16 * Method
17 * 1. Argument reduction:
18 * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
19 * Given x, find r and integer k such that
20 *
21 * x = k*ln2 + r, |r| <= 0.5*ln2.
22 *
23 * Here r will be represented as r = hi-lo for better
24 * accuracy.
25 *
26 * 2. Approximation of exp(r) by a special rational function on
27 * the interval [0,0.34658]:
28 * Write
29 * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
30 * We use a special Reme algorithm on [0,0.34658] to generate
31 * a polynomial of degree 5 to approximate R. The maximum error
32 * of this polynomial approximation is bounded by 2**-59. In
33 * other words,
34 * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
35 * (where z=r*r, and the values of P1 to P5 are listed below)
36 * and
37 * | 5 | -59
38 * | 2.0+P1*z+...+P5*z - R(z) | <= 2
39 * | |
40 * The computation of exp(r) thus becomes
41 * 2*r
42 * exp(r) = 1 + -------
43 * R - r
44 * r*R1(r)
45 * = 1 + r + ----------- (for better accuracy)
46 * 2 - R1(r)
47 * where
48 * 2 4 10
49 * R1(r) = r - (P1*r + P2*r + ... + P5*r ).
50 *
51 * 3. Scale back to obtain exp(x):
52 * From step 1, we have
53 * exp(x) = 2^k * exp(r)
54 *
55 * Special cases:
56 * exp(INF) is INF, exp(NaN) is NaN;
57 * exp(-INF) is 0, and
58 * for finite argument, only exp(0)=1 is exact.
59 *
60 * Accuracy:
61 * according to an error analysis, the error is always less than
62 * 1 ulp (unit in the last place).
63 *
64 * Misc. info.
65 * For IEEE double
66 * if x > 7.09782712893383973096e+02 then exp(x) overflow
67 * if x < -7.45133219101941108420e+02 then exp(x) underflow
68 *
69 * Constants:
70 * The hexadecimal values are the intended ones for the following
71 * constants. The decimal values may be used, provided that the
72 * compiler will convert from decimal to binary accurately enough
73 * to produce the hexadecimal values shown.
74 */
75
76#include "math_libm.h"
77#include "math_private.h"
78
79#ifdef __WATCOMC__ /* Watcom defines huge=__huge */
80#undef huge
81#endif
82
83static const double
84one = 1.0,
85halF[2] = {0.5,-0.5,},
86huge = 1.0e+300,
87twom1000= 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0*/
88o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
89u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */
90ln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
91 -6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */
92ln2LO[2] ={ 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
93 -1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */
94invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
95P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
96P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
97P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
98P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
99P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
100
101union {
102 Uint64 u64;
103 double d;
104} inf_union = {
105 SDL_UINT64_C(0x7ff0000000000000) /* Binary representation of a 64-bit infinite double (sign=0, exponent=2047, mantissa=0) */
106};
107
108double __ieee754_exp(double x) /* default IEEE double exp */
109{
110 double y;
111 double hi = 0.0;
112 double lo = 0.0;
113 double c;
114 double t;
115 int32_t k=0;
116 int32_t xsb;
117 u_int32_t hx;
118
119 GET_HIGH_WORD(hx,x);
120 xsb = (hx>>31)&1; /* sign bit of x */
121 hx &= 0x7fffffff; /* high word of |x| */
122
123 /* filter out non-finite argument */
124 if(hx >= 0x40862E42) { /* if |x|>=709.78... */
125 if(hx>=0x7ff00000) {
126 u_int32_t lx;
127 GET_LOW_WORD(lx,x);
128 if(((hx&0xfffff)|lx)!=0)
129 return x+x; /* NaN */
130 else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */
131 }
132 #if 1
133 if(x > o_threshold) return inf_union.d; /* overflow */
134 #elif 1
135 if(x > o_threshold) return huge*huge; /* overflow */
136 #else /* !!! FIXME: check this: "huge * huge" is a compiler warning, maybe they wanted +Inf? */
137 if(x > o_threshold) return INFINITY; /* overflow */
138 #endif
139
140 if(x < u_threshold) return twom1000*twom1000; /* underflow */
141 }
142
143 /* argument reduction */
144 if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
145 if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
146 hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
147 } else {
148 k = (int32_t) (invln2*x+halF[xsb]);
149 t = k;
150 hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */
151 lo = t*ln2LO[0];
152 }
153 x = hi - lo;
154 }
155 else if(hx < 0x3e300000) { /* when |x|<2**-28 */
156 if(huge+x>one) return one+x;/* trigger inexact */
157 }
158 else k = 0;
159
160 /* x is now in primary range */
161 t = x*x;
162 c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
163 if(k==0) return one-((x*c)/(c-2.0)-x);
164 else y = one-((lo-(x*c)/(2.0-c))-hi);
165 if(k >= -1021) {
166 u_int32_t hy;
167 GET_HIGH_WORD(hy,y);
168 SET_HIGH_WORD(y,hy+(k<<20)); /* add k to y's exponent */
169 return y;
170 } else {
171 u_int32_t hy;
172 GET_HIGH_WORD(hy,y);
173 SET_HIGH_WORD(y,hy+((k+1000)<<20)); /* add k to y's exponent */
174 return y*twom1000;
175 }
176}
177
178/*
179 * wrapper exp(x)
180 */
181#ifndef _IEEE_LIBM
182double exp(double x)
183{
184 static const double o_threshold = 7.09782712893383973096e+02; /* 0x40862E42, 0xFEFA39EF */
185 static const double u_threshold = -7.45133219101941108420e+02; /* 0xc0874910, 0xD52D3051 */
186
187 double z = __ieee754_exp(x);
188 if (_LIB_VERSION == _IEEE_)
189 return z;
190 if (isfinite(x)) {
191 if (x > o_threshold)
192 return __kernel_standard(x, x, 6); /* exp overflow */
193 if (x < u_threshold)
194 return __kernel_standard(x, x, 7); /* exp underflow */
195 }
196 return z;
197}
198#else
199strong_alias(__ieee754_exp, exp)
200#endif
201libm_hidden_def(exp)
diff --git a/contrib/SDL-3.2.8/src/libm/e_fmod.c b/contrib/SDL-3.2.8/src/libm/e_fmod.c
new file mode 100644
index 0000000..32c0249
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/e_fmod.c
@@ -0,0 +1,145 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/*
14 * __ieee754_fmod(x,y)
15 * Return x mod y in exact arithmetic
16 * Method: shift and subtract
17 */
18
19#include "math_libm.h"
20#include "math_private.h"
21
22static const double one = 1.0, Zero[] = {0.0, -0.0,};
23
24double attribute_hidden __ieee754_fmod(double x, double y)
25{
26 int32_t n,hx,hy,hz,ix,iy,sx,i;
27 u_int32_t lx,ly,lz;
28
29 EXTRACT_WORDS(hx,lx,x);
30 EXTRACT_WORDS(hy,ly,y);
31 sx = hx&0x80000000; /* sign of x */
32 hx ^=sx; /* |x| */
33 hy &= 0x7fffffff; /* |y| */
34
35 /* purge off exception values */
36 if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */
37 ((hy|((ly|-(int32_t)ly)>>31))>0x7ff00000)) /* or y is NaN */
38 return (x*y)/(x*y);
39 if(hx<=hy) {
40 if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
41 if(lx==ly)
42 return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/
43 }
44
45 /* determine ix = ilogb(x) */
46 if(hx<0x00100000) { /* subnormal x */
47 if(hx==0) {
48 for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
49 } else {
50 for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
51 }
52 } else ix = (hx>>20)-1023;
53
54 /* determine iy = ilogb(y) */
55 if(hy<0x00100000) { /* subnormal y */
56 if(hy==0) {
57 for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
58 } else {
59 for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
60 }
61 } else iy = (hy>>20)-1023;
62
63 /* set up {hx,lx}, {hy,ly} and align y to x */
64 if(ix >= -1022)
65 hx = 0x00100000|(0x000fffff&hx);
66 else { /* subnormal x, shift x to normal */
67 n = -1022-ix;
68 if(n<=31) {
69 hx = (hx<<n)|(lx>>(32-n));
70 lx <<= n;
71 } else {
72 hx = lx<<(n-32);
73 lx = 0;
74 }
75 }
76 if(iy >= -1022)
77 hy = 0x00100000|(0x000fffff&hy);
78 else { /* subnormal y, shift y to normal */
79 n = -1022-iy;
80 if(n<=31) {
81 hy = (hy<<n)|(ly>>(32-n));
82 ly <<= n;
83 } else {
84 hy = ly<<(n-32);
85 ly = 0;
86 }
87 }
88
89 /* fix point fmod */
90 n = ix - iy;
91 while(n--) {
92 hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
93 if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
94 else {
95 if((hz|lz)==0) /* return sign(x)*0 */
96 return Zero[(u_int32_t)sx>>31];
97 hx = hz+hz+(lz>>31); lx = lz+lz;
98 }
99 }
100 hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
101 if(hz>=0) {hx=hz;lx=lz;}
102
103 /* convert back to floating value and restore the sign */
104 if((hx|lx)==0) /* return sign(x)*0 */
105 return Zero[(u_int32_t)sx>>31];
106 while(hx<0x00100000) { /* normalize x */
107 hx = hx+hx+(lx>>31); lx = lx+lx;
108 iy -= 1;
109 }
110 if(iy>= -1022) { /* normalize output */
111 hx = ((hx-0x00100000)|((iy+1023)<<20));
112 INSERT_WORDS(x,hx|sx,lx);
113 } else { /* subnormal output */
114 n = -1022 - iy;
115 if(n<=20) {
116 lx = (lx>>n)|((u_int32_t)hx<<(32-n));
117 hx >>= n;
118 } else if (n<=31) {
119 lx = (hx<<(32-n))|(lx>>n); hx = sx;
120 } else {
121 lx = hx>>(n-32); hx = sx;
122 }
123 INSERT_WORDS(x,hx|sx,lx);
124 x *= one; /* create necessary signal */
125 }
126 return x; /* exact output */
127}
128
129/*
130 * wrapper fmod(x,y)
131 */
132#ifndef _IEEE_LIBM
133double fmod(double x, double y)
134{
135 double z = __ieee754_fmod(x, y);
136 if (_LIB_VERSION == _IEEE_ || isnan(y) || isnan(x))
137 return z;
138 if (y == 0.0)
139 return __kernel_standard(x, y, 27); /* fmod(x,0) */
140 return z;
141}
142#else
143strong_alias(__ieee754_fmod, fmod)
144#endif
145libm_hidden_def(fmod)
diff --git a/contrib/SDL-3.2.8/src/libm/e_log.c b/contrib/SDL-3.2.8/src/libm/e_log.c
new file mode 100644
index 0000000..f935fa7
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/e_log.c
@@ -0,0 +1,153 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13#if defined(_MSC_VER) /* Handle Microsoft VC++ compiler specifics. */
14/* C4723: potential divide by zero. */
15#pragma warning ( disable : 4723 )
16#endif
17
18/* __ieee754_log(x)
19 * Return the logrithm of x
20 *
21 * Method :
22 * 1. Argument Reduction: find k and f such that
23 * x = 2^k * (1+f),
24 * where sqrt(2)/2 < 1+f < sqrt(2) .
25 *
26 * 2. Approximation of log(1+f).
27 * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
28 * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
29 * = 2s + s*R
30 * We use a special Reme algorithm on [0,0.1716] to generate
31 * a polynomial of degree 14 to approximate R The maximum error
32 * of this polynomial approximation is bounded by 2**-58.45. In
33 * other words,
34 * 2 4 6 8 10 12 14
35 * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
36 * (the values of Lg1 to Lg7 are listed in the program)
37 * and
38 * | 2 14 | -58.45
39 * | Lg1*s +...+Lg7*s - R(z) | <= 2
40 * | |
41 * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
42 * In order to guarantee error in log below 1ulp, we compute log
43 * by
44 * log(1+f) = f - s*(f - R) (if f is not too large)
45 * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
46 *
47 * 3. Finally, log(x) = k*ln2 + log(1+f).
48 * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
49 * Here ln2 is split into two floating point number:
50 * ln2_hi + ln2_lo,
51 * where n*ln2_hi is always exact for |n| < 2000.
52 *
53 * Special cases:
54 * log(x) is NaN with signal if x < 0 (including -INF) ;
55 * log(+INF) is +INF; log(0) is -INF with signal;
56 * log(NaN) is that NaN with no signal.
57 *
58 * Accuracy:
59 * according to an error analysis, the error is always less than
60 * 1 ulp (unit in the last place).
61 *
62 * Constants:
63 * The hexadecimal values are the intended ones for the following
64 * constants. The decimal values may be used, provided that the
65 * compiler will convert from decimal to binary accurately enough
66 * to produce the hexadecimal values shown.
67 */
68
69#include "math_libm.h"
70#include "math_private.h"
71
72static const double
73ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
74ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
75two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
76Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
77Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
78Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
79Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
80Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
81Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
82Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
83
84static const double zero = 0.0;
85
86double attribute_hidden __ieee754_log(double x)
87{
88 double hfsq,f,s,z,R,w,t1,t2,dk;
89 int32_t k,hx,i,j;
90 u_int32_t lx;
91
92 EXTRACT_WORDS(hx,lx,x);
93
94 k=0;
95 if (hx < 0x00100000) { /* x < 2**-1022 */
96 if (((hx&0x7fffffff)|lx)==0)
97 return -two54/zero; /* log(+-0)=-inf */
98 if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
99 k -= 54; x *= two54; /* subnormal number, scale up x */
100 GET_HIGH_WORD(hx,x);
101 }
102 if (hx >= 0x7ff00000) return x+x;
103 k += (hx>>20)-1023;
104 hx &= 0x000fffff;
105 i = (hx+0x95f64)&0x100000;
106 SET_HIGH_WORD(x,hx|(i^0x3ff00000)); /* normalize x or x/2 */
107 k += (i>>20);
108 f = x-1.0;
109 if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */
110 if(f==zero) {if(k==0) return zero; else {dk=(double)k;
111 return dk*ln2_hi+dk*ln2_lo;}
112 }
113 R = f*f*(0.5-0.33333333333333333*f);
114 if(k==0) return f-R; else {dk=(double)k;
115 return dk*ln2_hi-((R-dk*ln2_lo)-f);}
116 }
117 s = f/(2.0+f);
118 dk = (double)k;
119 z = s*s;
120 i = hx-0x6147a;
121 w = z*z;
122 j = 0x6b851-hx;
123 t1= w*(Lg2+w*(Lg4+w*Lg6));
124 t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
125 i |= j;
126 R = t2+t1;
127 if(i>0) {
128 hfsq=0.5*f*f;
129 if(k==0) return f-(hfsq-s*(hfsq+R)); else
130 return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
131 } else {
132 if(k==0) return f-s*(f-R); else
133 return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
134 }
135}
136
137/*
138 * wrapper log(x)
139 */
140#ifndef _IEEE_LIBM
141double log(double x)
142{
143 double z = __ieee754_log(x);
144 if (_LIB_VERSION == _IEEE_ || isnan(x) || x > 0.0)
145 return z;
146 if (x == 0.0)
147 return __kernel_standard(x, x, 16); /* log(0) */
148 return __kernel_standard(x, x, 17); /* log(x<0) */
149}
150#else
151strong_alias(__ieee754_log, log)
152#endif
153libm_hidden_def(log)
diff --git a/contrib/SDL-3.2.8/src/libm/e_log10.c b/contrib/SDL-3.2.8/src/libm/e_log10.c
new file mode 100644
index 0000000..b6e736b
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/e_log10.c
@@ -0,0 +1,107 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13#if defined(_MSC_VER) /* Handle Microsoft VC++ compiler specifics. */
14/* C4723: potential divide by zero. */
15#pragma warning ( disable : 4723 )
16#endif
17
18/* __ieee754_log10(x)
19 * Return the base 10 logarithm of x
20 *
21 * Method :
22 * Let log10_2hi = leading 40 bits of log10(2) and
23 * log10_2lo = log10(2) - log10_2hi,
24 * ivln10 = 1/log(10) rounded.
25 * Then
26 * n = ilogb(x),
27 * if(n<0) n = n+1;
28 * x = scalbn(x,-n);
29 * log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))
30 *
31 * Note 1:
32 * To guarantee log10(10**n)=n, where 10**n is normal, the rounding
33 * mode must set to Round-to-Nearest.
34 * Note 2:
35 * [1/log(10)] rounded to 53 bits has error .198 ulps;
36 * log10 is monotonic at all binary break points.
37 *
38 * Special cases:
39 * log10(x) is NaN with signal if x < 0;
40 * log10(+INF) is +INF with no signal; log10(0) is -INF with signal;
41 * log10(NaN) is that NaN with no signal;
42 * log10(10**N) = N for N=0,1,...,22.
43 *
44 * Constants:
45 * The hexadecimal values are the intended ones for the following constants.
46 * The decimal values may be used, provided that the compiler will convert
47 * from decimal to binary accurately enough to produce the hexadecimal values
48 * shown.
49 */
50
51#include "math_libm.h"
52#include "math_private.h"
53
54static const double
55two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
56ivln10 = 4.34294481903251816668e-01, /* 0x3FDBCB7B, 0x1526E50E */
57log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
58log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
59
60static const double zero = 0.0;
61
62double attribute_hidden __ieee754_log10(double x)
63{
64 double y,z;
65 int32_t i,k,hx;
66 u_int32_t lx;
67
68 EXTRACT_WORDS(hx,lx,x);
69
70 k=0;
71 if (hx < 0x00100000) { /* x < 2**-1022 */
72 if (((hx&0x7fffffff)|lx)==0)
73 return -two54/zero; /* log(+-0)=-inf */
74 if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
75 k -= 54; x *= two54; /* subnormal number, scale up x */
76 GET_HIGH_WORD(hx,x);
77 }
78 if (hx >= 0x7ff00000) return x+x;
79 k += (hx>>20)-1023;
80 i = ((u_int32_t)k&0x80000000)>>31;
81 hx = (hx&0x000fffff)|((0x3ff-i)<<20);
82 y = (double)(k+i);
83 SET_HIGH_WORD(x,hx);
84 z = y*log10_2lo + ivln10*__ieee754_log(x);
85 return z+y*log10_2hi;
86}
87
88/*
89 * wrapper log10(X)
90 */
91#ifndef _IEEE_LIBM
92double log10(double x)
93{
94 double z = __ieee754_log10(x);
95 if (_LIB_VERSION == _IEEE_ || isnan(x))
96 return z;
97 if (x <= 0.0) {
98 if(x == 0.0)
99 return __kernel_standard(x, x, 18); /* log10(0) */
100 return __kernel_standard(x, x, 19); /* log10(x<0) */
101 }
102 return z;
103}
104#else
105strong_alias(__ieee754_log10, log10)
106#endif
107libm_hidden_def(log10)
diff --git a/contrib/SDL-3.2.8/src/libm/e_pow.c b/contrib/SDL-3.2.8/src/libm/e_pow.c
new file mode 100644
index 0000000..d1a141e
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/e_pow.c
@@ -0,0 +1,348 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/* __ieee754_pow(x,y) return x**y
14 *
15 * n
16 * Method: Let x = 2 * (1+f)
17 * 1. Compute and return log2(x) in two pieces:
18 * log2(x) = w1 + w2,
19 * where w1 has 53-24 = 29 bit trailing zeros.
20 * 2. Perform y*log2(x) = n+y' by simulating muti-precision
21 * arithmetic, where |y'|<=0.5.
22 * 3. Return x**y = 2**n*exp(y'*log2)
23 *
24 * Special cases:
25 * 1. +-1 ** anything is 1.0
26 * 2. +-1 ** +-INF is 1.0
27 * 3. (anything) ** 0 is 1
28 * 4. (anything) ** 1 is itself
29 * 5. (anything) ** NAN is NAN
30 * 6. NAN ** (anything except 0) is NAN
31 * 7. +-(|x| > 1) ** +INF is +INF
32 * 8. +-(|x| > 1) ** -INF is +0
33 * 9. +-(|x| < 1) ** +INF is +0
34 * 10 +-(|x| < 1) ** -INF is +INF
35 * 11. +0 ** (+anything except 0, NAN) is +0
36 * 12. -0 ** (+anything except 0, NAN, odd integer) is +0
37 * 13. +0 ** (-anything except 0, NAN) is +INF
38 * 14. -0 ** (-anything except 0, NAN, odd integer) is +INF
39 * 15. -0 ** (odd integer) = -( +0 ** (odd integer) )
40 * 16. +INF ** (+anything except 0,NAN) is +INF
41 * 17. +INF ** (-anything except 0,NAN) is +0
42 * 18. -INF ** (anything) = -0 ** (-anything)
43 * 19. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
44 * 20. (-anything except 0 and inf) ** (non-integer) is NAN
45 *
46 * Accuracy:
47 * pow(x,y) returns x**y nearly rounded. In particular
48 * pow(integer,integer)
49 * always returns the correct integer provided it is
50 * representable.
51 *
52 * Constants :
53 * The hexadecimal values are the intended ones for the following
54 * constants. The decimal values may be used, provided that the
55 * compiler will convert from decimal to binary accurately enough
56 * to produce the hexadecimal values shown.
57 */
58
59#include "math_libm.h"
60#include "math_private.h"
61
62#if defined(_MSC_VER) /* Handle Microsoft VC++ compiler specifics. */
63/* C4756: overflow in constant arithmetic */
64#pragma warning ( disable : 4756 )
65#endif
66
67#ifdef __WATCOMC__ /* Watcom defines huge=__huge */
68#undef huge
69#endif
70
71static const double
72bp[] = {1.0, 1.5,},
73dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
74dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
75zero = 0.0,
76one = 1.0,
77two = 2.0,
78two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */
79huge = 1.0e300,
80tiny = 1.0e-300,
81 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
82L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
83L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
84L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
85L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
86L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
87L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
88P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
89P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
90P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
91P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
92P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
93lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
94lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
95lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
96ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */
97cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
98cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
99cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/
100ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
101ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
102ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
103
104double attribute_hidden __ieee754_pow(double x, double y)
105{
106 double z,ax,z_h,z_l,p_h,p_l;
107 double y1,t1,t2,r,s,t,u,v,w;
108 int32_t i,j,k,yisint,n;
109 int32_t hx,hy,ix,iy;
110 u_int32_t lx,ly;
111
112 EXTRACT_WORDS(hx,lx,x);
113 /* x==1: 1**y = 1 (even if y is NaN) */
114 if (hx==0x3ff00000 && lx==0) {
115 return x;
116 }
117 ix = hx&0x7fffffff;
118
119 EXTRACT_WORDS(hy,ly,y);
120 iy = hy&0x7fffffff;
121
122 /* y==zero: x**0 = 1 */
123 if((iy|ly)==0) return one;
124
125 /* +-NaN return x+y */
126 if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
127 iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
128 return x+y;
129
130 /* determine if y is an odd int when x < 0
131 * yisint = 0 ... y is not an integer
132 * yisint = 1 ... y is an odd int
133 * yisint = 2 ... y is an even int
134 */
135 yisint = 0;
136 if(hx<0) {
137 if(iy>=0x43400000) yisint = 2; /* even integer y */
138 else if(iy>=0x3ff00000) {
139 k = (iy>>20)-0x3ff; /* exponent */
140 if(k>20) {
141 j = ly>>(52-k);
142 if(((u_int32_t)j<<(52-k))==ly) yisint = 2-(j&1);
143 } else if(ly==0) {
144 j = iy>>(20-k);
145 if((j<<(20-k))==iy) yisint = 2-(j&1);
146 }
147 }
148 }
149
150 /* special value of y */
151 if(ly==0) {
152 if (iy==0x7ff00000) { /* y is +-inf */
153 if (((ix-0x3ff00000)|lx)==0)
154 return one; /* +-1**+-inf is 1 (yes, weird rule) */
155 if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */
156 return (hy>=0) ? y : zero;
157 /* (|x|<1)**-,+inf = inf,0 */
158 return (hy<0) ? -y : zero;
159 }
160 if(iy==0x3ff00000) { /* y is +-1 */
161 if(hy<0) return one/x; else return x;
162 }
163 if(hy==0x40000000) return x*x; /* y is 2 */
164 if(hy==0x3fe00000) { /* y is 0.5 */
165 if(hx>=0) /* x >= +0 */
166 return __ieee754_sqrt(x);
167 }
168 }
169
170 ax = fabs(x);
171 /* special value of x */
172 if(lx==0) {
173 if(ix==0x7ff00000||ix==0||ix==0x3ff00000){
174 z = ax; /*x is +-0,+-inf,+-1*/
175 if(hy<0) z = one/z; /* z = (1/|x|) */
176 if(hx<0) {
177 if(((ix-0x3ff00000)|yisint)==0) {
178 z = (z-z)/(z-z); /* (-1)**non-int is NaN */
179 } else if(yisint==1)
180 z = -z; /* (x<0)**odd = -(|x|**odd) */
181 }
182 return z;
183 }
184 }
185
186 /* (x<0)**(non-int) is NaN */
187 if(((((u_int32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x);
188
189 /* |y| is huge */
190 if(iy>0x41e00000) { /* if |y| > 2**31 */
191 if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */
192 if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
193 if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
194 }
195 /* over/underflow if x is not close to one */
196 if(ix<0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
197 if(ix>0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
198 /* now |1-x| is tiny <= 2**-20, suffice to compute
199 log(x) by x-x^2/2+x^3/3-x^4/4 */
200 t = x-1; /* t has 20 trailing zeros */
201 w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
202 u = ivln2_h*t; /* ivln2_h has 21 sig. bits */
203 v = t*ivln2_l-w*ivln2;
204 t1 = u+v;
205 SET_LOW_WORD(t1,0);
206 t2 = v-(t1-u);
207 } else {
208 double s2,s_h,s_l,t_h,t_l;
209 n = 0;
210 /* take care subnormal number */
211 if(ix<0x00100000)
212 {ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); }
213 n += ((ix)>>20)-0x3ff;
214 j = ix&0x000fffff;
215 /* determine interval */
216 ix = j|0x3ff00000; /* normalize ix */
217 if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */
218 else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */
219 else {k=0;n+=1;ix -= 0x00100000;}
220 SET_HIGH_WORD(ax,ix);
221
222 /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
223 u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
224 v = one/(ax+bp[k]);
225 s = u*v;
226 s_h = s;
227 SET_LOW_WORD(s_h,0);
228 /* t_h=ax+bp[k] High */
229 t_h = zero;
230 SET_HIGH_WORD(t_h,((ix>>1)|0x20000000)+0x00080000+(k<<18));
231 t_l = ax - (t_h-bp[k]);
232 s_l = v*((u-s_h*t_h)-s_h*t_l);
233 /* compute log(ax) */
234 s2 = s*s;
235 r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
236 r += s_l*(s_h+s);
237 s2 = s_h*s_h;
238 t_h = 3.0+s2+r;
239 SET_LOW_WORD(t_h,0);
240 t_l = r-((t_h-3.0)-s2);
241 /* u+v = s*(1+...) */
242 u = s_h*t_h;
243 v = s_l*t_h+t_l*s;
244 /* 2/(3log2)*(s+...) */
245 p_h = u+v;
246 SET_LOW_WORD(p_h,0);
247 p_l = v-(p_h-u);
248 z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */
249 z_l = cp_l*p_h+p_l*cp+dp_l[k];
250 /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
251 t = (double)n;
252 t1 = (((z_h+z_l)+dp_h[k])+t);
253 SET_LOW_WORD(t1,0);
254 t2 = z_l-(((t1-t)-dp_h[k])-z_h);
255 }
256
257 s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
258 if(((((u_int32_t)hx>>31)-1)|(yisint-1))==0)
259 s = -one;/* (-ve)**(odd int) */
260
261 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
262 y1 = y;
263 SET_LOW_WORD(y1,0);
264 p_l = (y-y1)*t1+y*t2;
265 p_h = y1*t1;
266 z = p_l+p_h;
267 EXTRACT_WORDS(j,i,z);
268 if (j>=0x40900000) { /* z >= 1024 */
269 if(((j-0x40900000)|i)!=0) /* if z > 1024 */
270 return s*huge*huge; /* overflow */
271 else {
272 if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */
273 }
274 } else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */
275 if(((j-0xc090cc00)|i)!=0) /* z < -1075 */
276 return s*tiny*tiny; /* underflow */
277 else {
278 if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */
279 }
280 }
281 /*
282 * compute 2**(p_h+p_l)
283 */
284 i = j&0x7fffffff;
285 k = (i>>20)-0x3ff;
286 n = 0;
287 if(i>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
288 n = j+(0x00100000>>(k+1));
289 k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */
290 t = zero;
291 SET_HIGH_WORD(t,n&~(0x000fffff>>k));
292 n = ((n&0x000fffff)|0x00100000)>>(20-k);
293 if(j<0) n = -n;
294 p_h -= t;
295 }
296 t = p_l+p_h;
297 SET_LOW_WORD(t,0);
298 u = t*lg2_h;
299 v = (p_l-(t-p_h))*lg2+t*lg2_l;
300 z = u+v;
301 w = v-(z-u);
302 t = z*z;
303 t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
304 r = (z*t1)/(t1-two)-(w+z*w);
305 z = one-(r-z);
306 GET_HIGH_WORD(j,z);
307 j += (n<<20);
308 if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */
309 else SET_HIGH_WORD(z,j);
310 return s*z;
311}
312
313/*
314 * wrapper pow(x,y) return x**y
315 */
316#ifndef _IEEE_LIBM
317double pow(double x, double y)
318{
319 double z = __ieee754_pow(x, y);
320 if (_LIB_VERSION == _IEEE_|| isnan(y))
321 return z;
322 if (isnan(x)) {
323 if (y == 0.0)
324 return __kernel_standard(x, y, 42); /* pow(NaN,0.0) */
325 return z;
326 }
327 if (x == 0.0) {
328 if (y == 0.0)
329 return __kernel_standard(x, y, 20); /* pow(0.0,0.0) */
330 if (isfinite(y) && y < 0.0)
331 return __kernel_standard(x,y,23); /* pow(0.0,negative) */
332 return z;
333 }
334 if (!isfinite(z)) {
335 if (isfinite(x) && isfinite(y)) {
336 if (isnan(z))
337 return __kernel_standard(x, y, 24); /* pow neg**non-int */
338 return __kernel_standard(x, y, 21); /* pow overflow */
339 }
340 }
341 if (z == 0.0 && isfinite(x) && isfinite(y))
342 return __kernel_standard(x, y, 22); /* pow underflow */
343 return z;
344}
345#else
346strong_alias(__ieee754_pow, pow)
347#endif
348libm_hidden_def(pow)
diff --git a/contrib/SDL-3.2.8/src/libm/e_rem_pio2.c b/contrib/SDL-3.2.8/src/libm/e_rem_pio2.c
new file mode 100644
index 0000000..851560f
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/e_rem_pio2.c
@@ -0,0 +1,162 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/* __ieee754_rem_pio2(x,y)
14 *
15 * return the remainder of x rem pi/2 in y[0]+y[1]
16 * use __kernel_rem_pio2()
17 */
18
19#include "math_libm.h"
20#include "math_private.h"
21
22/*
23 * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
24 */
25static const int32_t two_over_pi[] = {
260xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62,
270x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A,
280x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
290xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41,
300x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8,
310x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
320x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5,
330xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08,
340x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
350x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880,
360x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
37};
38
39static const int32_t npio2_hw[] = {
400x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C,
410x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C,
420x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A,
430x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
440x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB,
450x404858EB, 0x404921FB,
46};
47
48/*
49 * invpio2: 53 bits of 2/pi
50 * pio2_1: first 33 bit of pi/2
51 * pio2_1t: pi/2 - pio2_1
52 * pio2_2: second 33 bit of pi/2
53 * pio2_2t: pi/2 - (pio2_1+pio2_2)
54 * pio2_3: third 33 bit of pi/2
55 * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
56 */
57
58static const double
59zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
60half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
61two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
62invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
63pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
64pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
65pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
66pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
67pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
68pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
69
70int32_t attribute_hidden __ieee754_rem_pio2(double x, double *y)
71{
72 double z=0.0,w,t,r,fn;
73 double tx[3];
74 int32_t e0,i,j,nx,n,ix,hx;
75 u_int32_t low;
76
77 GET_HIGH_WORD(hx,x); /* high word of x */
78 ix = hx&0x7fffffff;
79 if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */
80 {y[0] = x; y[1] = 0; return 0;}
81 if(ix<0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */
82 if(hx>0) {
83 z = x - pio2_1;
84 if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */
85 y[0] = z - pio2_1t;
86 y[1] = (z-y[0])-pio2_1t;
87 } else { /* near pi/2, use 33+33+53 bit pi */
88 z -= pio2_2;
89 y[0] = z - pio2_2t;
90 y[1] = (z-y[0])-pio2_2t;
91 }
92 return 1;
93 } else { /* negative x */
94 z = x + pio2_1;
95 if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */
96 y[0] = z + pio2_1t;
97 y[1] = (z-y[0])+pio2_1t;
98 } else { /* near pi/2, use 33+33+53 bit pi */
99 z += pio2_2;
100 y[0] = z + pio2_2t;
101 y[1] = (z-y[0])+pio2_2t;
102 }
103 return -1;
104 }
105 }
106 if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
107 t = fabs(x);
108 n = (int32_t) (t*invpio2+half);
109 fn = (double)n;
110 r = t-fn*pio2_1;
111 w = fn*pio2_1t; /* 1st round good to 85 bit */
112 if(n<32&&ix!=npio2_hw[n-1]) {
113 y[0] = r-w; /* quick check no cancellation */
114 } else {
115 u_int32_t high;
116 j = ix>>20;
117 y[0] = r-w;
118 GET_HIGH_WORD(high,y[0]);
119 i = j-((high>>20)&0x7ff);
120 if(i>16) { /* 2nd iteration needed, good to 118 */
121 t = r;
122 w = fn*pio2_2;
123 r = t-w;
124 w = fn*pio2_2t-((t-r)-w);
125 y[0] = r-w;
126 GET_HIGH_WORD(high,y[0]);
127 i = j-((high>>20)&0x7ff);
128 if(i>49) { /* 3rd iteration need, 151 bits acc */
129 t = r; /* will cover all possible cases */
130 w = fn*pio2_3;
131 r = t-w;
132 w = fn*pio2_3t-((t-r)-w);
133 y[0] = r-w;
134 }
135 }
136 }
137 y[1] = (r-y[0])-w;
138 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
139 else return n;
140 }
141 /*
142 * all other (large) arguments
143 */
144 if(ix>=0x7ff00000) { /* x is inf or NaN */
145 y[0]=y[1]=x-x; return 0;
146 }
147 /* set z = scalbn(|x|,ilogb(x)-23) */
148 GET_LOW_WORD(low,x);
149 SET_LOW_WORD(z,low);
150 e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */
151 SET_HIGH_WORD(z, ix - ((int32_t)(e0<<20)));
152 for(i=0;i<2;i++) {
153 tx[i] = (double)((int32_t)(z));
154 z = (z-tx[i])*two24;
155 }
156 tx[2] = z;
157 nx = 3;
158 while((nx > 0) && tx[nx-1]==zero) nx--; /* skip zero term */
159 n = __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi);
160 if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
161 return n;
162}
diff --git a/contrib/SDL-3.2.8/src/libm/e_sqrt.c b/contrib/SDL-3.2.8/src/libm/e_sqrt.c
new file mode 100644
index 0000000..8ac58c6
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/e_sqrt.c
@@ -0,0 +1,458 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/* __ieee754_sqrt(x)
14 * Return correctly rounded sqrt.
15 * ------------------------------------------
16 * | Use the hardware sqrt if you have one |
17 * ------------------------------------------
18 * Method:
19 * Bit by bit method using integer arithmetic. (Slow, but portable)
20 * 1. Normalization
21 * Scale x to y in [1,4) with even powers of 2:
22 * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
23 * sqrt(x) = 2^k * sqrt(y)
24 * 2. Bit by bit computation
25 * Let q = sqrt(y) truncated to i bit after binary point (q = 1),
26 * i 0
27 * i+1 2
28 * s = 2*q , and y = 2 * ( y - q ). (1)
29 * i i i i
30 *
31 * To compute q from q , one checks whether
32 * i+1 i
33 *
34 * -(i+1) 2
35 * (q + 2 ) <= y. (2)
36 * i
37 * -(i+1)
38 * If (2) is false, then q = q ; otherwise q = q + 2 .
39 * i+1 i i+1 i
40 *
41 * With some algebric manipulation, it is not difficult to see
42 * that (2) is equivalent to
43 * -(i+1)
44 * s + 2 <= y (3)
45 * i i
46 *
47 * The advantage of (3) is that s and y can be computed by
48 * i i
49 * the following recurrence formula:
50 * if (3) is false
51 *
52 * s = s , y = y ; (4)
53 * i+1 i i+1 i
54 *
55 * otherwise,
56 * -i -(i+1)
57 * s = s + 2 , y = y - s - 2 (5)
58 * i+1 i i+1 i i
59 *
60 * One may easily use induction to prove (4) and (5).
61 * Note. Since the left hand side of (3) contain only i+2 bits,
62 * it does not necessary to do a full (53-bit) comparison
63 * in (3).
64 * 3. Final rounding
65 * After generating the 53 bits result, we compute one more bit.
66 * Together with the remainder, we can decide whether the
67 * result is exact, bigger than 1/2ulp, or less than 1/2ulp
68 * (it will never equal to 1/2ulp).
69 * The rounding mode can be detected by checking whether
70 * huge + tiny is equal to huge, and whether huge - tiny is
71 * equal to huge for some floating point number "huge" and "tiny".
72 *
73 * Special cases:
74 * sqrt(+-0) = +-0 ... exact
75 * sqrt(inf) = inf
76 * sqrt(-ve) = NaN ... with invalid signal
77 * sqrt(NaN) = NaN ... with invalid signal for signaling NaN
78 *
79 * Other methods : see the appended file at the end of the program below.
80 *---------------
81 */
82
83#include "math_libm.h"
84#include "math_private.h"
85
86static const double one = 1.0, tiny = 1.0e-300;
87
88double attribute_hidden __ieee754_sqrt(double x)
89{
90 double z;
91 int32_t sign = (int)0x80000000;
92 int32_t ix0,s0,q,m,t,i;
93 u_int32_t r,t1,s1,ix1,q1;
94
95 EXTRACT_WORDS(ix0,ix1,x);
96
97 /* take care of Inf and NaN */
98 if((ix0&0x7ff00000)==0x7ff00000) {
99 return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
100 sqrt(-inf)=sNaN */
101 }
102 /* take care of zero */
103 if(ix0<=0) {
104 if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */
105 else if(ix0<0)
106 return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
107 }
108 /* normalize x */
109 m = (ix0>>20);
110 if(m==0) { /* subnormal x */
111 while(ix0==0) {
112 m -= 21;
113 ix0 |= (ix1>>11); ix1 <<= 21;
114 }
115 for(i=0;(ix0&0x00100000)==0;i++) ix0<<=1;
116 m -= i-1;
117 ix0 |= (ix1>>(32-i));
118 ix1 <<= i;
119 }
120 m -= 1023; /* unbias exponent */
121 ix0 = (ix0&0x000fffff)|0x00100000;
122 if(m&1){ /* odd m, double x to make it even */
123 ix0 += ix0 + ((ix1&sign)>>31);
124 ix1 += ix1;
125 }
126 m >>= 1; /* m = [m/2] */
127
128 /* generate sqrt(x) bit by bit */
129 ix0 += ix0 + ((ix1&sign)>>31);
130 ix1 += ix1;
131 q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
132 r = 0x00200000; /* r = moving bit from right to left */
133
134 while(r!=0) {
135 t = s0+r;
136 if(t<=ix0) {
137 s0 = t+r;
138 ix0 -= t;
139 q += r;
140 }
141 ix0 += ix0 + ((ix1&sign)>>31);
142 ix1 += ix1;
143 r>>=1;
144 }
145
146 r = sign;
147 while(r!=0) {
148 t1 = s1+r;
149 t = s0;
150 if((t<ix0)||((t==ix0)&&(t1<=ix1))) {
151 s1 = t1+r;
152 if(((t1&sign)==(u_int32_t)sign)&&(s1&sign)==0) s0 += 1;
153 ix0 -= t;
154 if (ix1 < t1) ix0 -= 1;
155 ix1 -= t1;
156 q1 += r;
157 }
158 ix0 += ix0 + ((ix1&sign)>>31);
159 ix1 += ix1;
160 r>>=1;
161 }
162
163 /* use floating add to find out rounding direction */
164 if((ix0|ix1)!=0) {
165 z = one-tiny; /* trigger inexact flag */
166 if (z>=one) {
167 z = one+tiny;
168 if (q1==(u_int32_t)0xffffffff) { q1=0; q += 1;}
169 else if (z>one) {
170 if (q1==(u_int32_t)0xfffffffe) q+=1;
171 q1+=2;
172 } else
173 q1 += (q1&1);
174 }
175 }
176 ix0 = (q>>1)+0x3fe00000;
177 ix1 = q1>>1;
178 if ((q&1)==1) ix1 |= sign;
179 ix0 += (m <<20);
180 INSERT_WORDS(z,ix0,ix1);
181 return z;
182}
183
184/*
185 * wrapper sqrt(x)
186 */
187#ifndef _IEEE_LIBM
188double sqrt(double x)
189{
190 double z = __ieee754_sqrt(x);
191 if (_LIB_VERSION == _IEEE_ || isnan(x))
192 return z;
193 if (x < 0.0)
194 return __kernel_standard(x, x, 26); /* sqrt(negative) */
195 return z;
196}
197#else
198strong_alias(__ieee754_sqrt, sqrt)
199#endif
200libm_hidden_def(sqrt)
201
202
203/*
204Other methods (use floating-point arithmetic)
205-------------
206(This is a copy of a drafted paper by Prof W. Kahan
207and K.C. Ng, written in May, 1986)
208
209 Two algorithms are given here to implement sqrt(x)
210 (IEEE double precision arithmetic) in software.
211 Both supply sqrt(x) correctly rounded. The first algorithm (in
212 Section A) uses newton iterations and involves four divisions.
213 The second one uses reciproot iterations to avoid division, but
214 requires more multiplications. Both algorithms need the ability
215 to chop results of arithmetic operations instead of round them,
216 and the INEXACT flag to indicate when an arithmetic operation
217 is executed exactly with no roundoff error, all part of the
218 standard (IEEE 754-1985). The ability to perform shift, add,
219 subtract and logical AND operations upon 32-bit words is needed
220 too, though not part of the standard.
221
222A. sqrt(x) by Newton Iteration
223
224 (1) Initial approximation
225
226 Let x0 and x1 be the leading and the trailing 32-bit words of
227 a floating point number x (in IEEE double format) respectively
228
229 1 11 52 ...widths
230 ------------------------------------------------------
231 x: |s| e | f |
232 ------------------------------------------------------
233 msb lsb msb lsb ...order
234
235
236 ------------------------ ------------------------
237 x0: |s| e | f1 | x1: | f2 |
238 ------------------------ ------------------------
239
240 By performing shifts and subtracts on x0 and x1 (both regarded
241 as integers), we obtain an 8-bit approximation of sqrt(x) as
242 follows.
243
244 k := (x0>>1) + 0x1ff80000;
245 y0 := k - T1[31&(k>>15)]. ... y ~ sqrt(x) to 8 bits
246 Here k is a 32-bit integer and T1[] is an integer array containing
247 correction terms. Now magically the floating value of y (y's
248 leading 32-bit word is y0, the value of its trailing word is 0)
249 approximates sqrt(x) to almost 8-bit.
250
251 Value of T1:
252 static int T1[32]= {
253 0, 1024, 3062, 5746, 9193, 13348, 18162, 23592,
254 29598, 36145, 43202, 50740, 58733, 67158, 75992, 85215,
255 83599, 71378, 60428, 50647, 41945, 34246, 27478, 21581,
256 16499, 12183, 8588, 5674, 3403, 1742, 661, 130,};
257
258 (2) Iterative refinement
259
260 Apply Heron's rule three times to y, we have y approximates
261 sqrt(x) to within 1 ulp (Unit in the Last Place):
262
263 y := (y+x/y)/2 ... almost 17 sig. bits
264 y := (y+x/y)/2 ... almost 35 sig. bits
265 y := y-(y-x/y)/2 ... within 1 ulp
266
267
268 Remark 1.
269 Another way to improve y to within 1 ulp is:
270
271 y := (y+x/y) ... almost 17 sig. bits to 2*sqrt(x)
272 y := y - 0x00100006 ... almost 18 sig. bits to sqrt(x)
273
274 2
275 (x-y )*y
276 y := y + 2* ---------- ...within 1 ulp
277 2
278 3y + x
279
280
281 This formula has one division fewer than the one above; however,
282 it requires more multiplications and additions. Also x must be
283 scaled in advance to avoid spurious overflow in evaluating the
284 expression 3y*y+x. Hence it is not recommended uless division
285 is slow. If division is very slow, then one should use the
286 reciproot algorithm given in section B.
287
288 (3) Final adjustment
289
290 By twiddling y's last bit it is possible to force y to be
291 correctly rounded according to the prevailing rounding mode
292 as follows. Let r and i be copies of the rounding mode and
293 inexact flag before entering the square root program. Also we
294 use the expression y+-ulp for the next representable floating
295 numbers (up and down) of y. Note that y+-ulp = either fixed
296 point y+-1, or multiply y by nextafter(1,+-inf) in chopped
297 mode.
298
299 I := FALSE; ... reset INEXACT flag I
300 R := RZ; ... set rounding mode to round-toward-zero
301 z := x/y; ... chopped quotient, possibly inexact
302 If(not I) then { ... if the quotient is exact
303 if(z=y) {
304 I := i; ... restore inexact flag
305 R := r; ... restore rounded mode
306 return sqrt(x):=y.
307 } else {
308 z := z - ulp; ... special rounding
309 }
310 }
311 i := TRUE; ... sqrt(x) is inexact
312 If (r=RN) then z=z+ulp ... rounded-to-nearest
313 If (r=RP) then { ... round-toward-+inf
314 y = y+ulp; z=z+ulp;
315 }
316 y := y+z; ... chopped sum
317 y0:=y0-0x00100000; ... y := y/2 is correctly rounded.
318 I := i; ... restore inexact flag
319 R := r; ... restore rounded mode
320 return sqrt(x):=y.
321
322 (4) Special cases
323
324 Square root of +inf, +-0, or NaN is itself;
325 Square root of a negative number is NaN with invalid signal.
326
327
328B. sqrt(x) by Reciproot Iteration
329
330 (1) Initial approximation
331
332 Let x0 and x1 be the leading and the trailing 32-bit words of
333 a floating point number x (in IEEE double format) respectively
334 (see section A). By performing shifs and subtracts on x0 and y0,
335 we obtain a 7.8-bit approximation of 1/sqrt(x) as follows.
336
337 k := 0x5fe80000 - (x0>>1);
338 y0:= k - T2[63&(k>>14)]. ... y ~ 1/sqrt(x) to 7.8 bits
339
340 Here k is a 32-bit integer and T2[] is an integer array
341 containing correction terms. Now magically the floating
342 value of y (y's leading 32-bit word is y0, the value of
343 its trailing word y1 is set to zero) approximates 1/sqrt(x)
344 to almost 7.8-bit.
345
346 Value of T2:
347 static int T2[64]= {
348 0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,
349 0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
350 0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
351 0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
352 0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
353 0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
354 0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
355 0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd,};
356
357 (2) Iterative refinement
358
359 Apply Reciproot iteration three times to y and multiply the
360 result by x to get an approximation z that matches sqrt(x)
361 to about 1 ulp. To be exact, we will have
362 -1ulp < sqrt(x)-z<1.0625ulp.
363
364 ... set rounding mode to Round-to-nearest
365 y := y*(1.5-0.5*x*y*y) ... almost 15 sig. bits to 1/sqrt(x)
366 y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
367 ... special arrangement for better accuracy
368 z := x*y ... 29 bits to sqrt(x), with z*y<1
369 z := z + 0.5*z*(1-z*y) ... about 1 ulp to sqrt(x)
370
371 Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that
372 (a) the term z*y in the final iteration is always less than 1;
373 (b) the error in the final result is biased upward so that
374 -1 ulp < sqrt(x) - z < 1.0625 ulp
375 instead of |sqrt(x)-z|<1.03125ulp.
376
377 (3) Final adjustment
378
379 By twiddling y's last bit it is possible to force y to be
380 correctly rounded according to the prevailing rounding mode
381 as follows. Let r and i be copies of the rounding mode and
382 inexact flag before entering the square root program. Also we
383 use the expression y+-ulp for the next representable floating
384 numbers (up and down) of y. Note that y+-ulp = either fixed
385 point y+-1, or multiply y by nextafter(1,+-inf) in chopped
386 mode.
387
388 R := RZ; ... set rounding mode to round-toward-zero
389 switch(r) {
390 case RN: ... round-to-nearest
391 if(x<= z*(z-ulp)...chopped) z = z - ulp; else
392 if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp;
393 break;
394 case RZ:case RM: ... round-to-zero or round-to--inf
395 R:=RP; ... reset rounding mod to round-to-+inf
396 if(x<z*z ... rounded up) z = z - ulp; else
397 if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp;
398 break;
399 case RP: ... round-to-+inf
400 if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else
401 if(x>z*z ...chopped) z = z+ulp;
402 break;
403 }
404
405 Remark 3. The above comparisons can be done in fixed point. For
406 example, to compare x and w=z*z chopped, it suffices to compare
407 x1 and w1 (the trailing parts of x and w), regarding them as
408 two's complement integers.
409
410 ...Is z an exact square root?
411 To determine whether z is an exact square root of x, let z1 be the
412 trailing part of z, and also let x0 and x1 be the leading and
413 trailing parts of x.
414
415 If ((z1&0x03ffffff)!=0) ... not exact if trailing 26 bits of z!=0
416 I := 1; ... Raise Inexact flag: z is not exact
417 else {
418 j := 1 - [(x0>>20)&1] ... j = logb(x) mod 2
419 k := z1 >> 26; ... get z's 25-th and 26-th
420 fraction bits
421 I := i or (k&j) or ((k&(j+j+1))!=(x1&3));
422 }
423 R:= r ... restore rounded mode
424 return sqrt(x):=z.
425
426 If multiplication is cheaper then the foregoing red tape, the
427 Inexact flag can be evaluated by
428
429 I := i;
430 I := (z*z!=x) or I.
431
432 Note that z*z can overwrite I; this value must be sensed if it is
433 True.
434
435 Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be
436 zero.
437
438 --------------------
439 z1: | f2 |
440 --------------------
441 bit 31 bit 0
442
443 Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd
444 or even of logb(x) have the following relations:
445
446 -------------------------------------------------
447 bit 27,26 of z1 bit 1,0 of x1 logb(x)
448 -------------------------------------------------
449 00 00 odd and even
450 01 01 even
451 10 10 odd
452 10 00 even
453 11 01 even
454 -------------------------------------------------
455
456 (4) Special cases (see (4) of Section A).
457
458 */
diff --git a/contrib/SDL-3.2.8/src/libm/k_cos.c b/contrib/SDL-3.2.8/src/libm/k_cos.c
new file mode 100644
index 0000000..399a40c
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/k_cos.c
@@ -0,0 +1,83 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/*
14 * __kernel_cos( x, y )
15 * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
16 * Input x is assumed to be bounded by ~pi/4 in magnitude.
17 * Input y is the tail of x.
18 *
19 * Algorithm
20 * 1. Since cos(-x) = cos(x), we need only to consider positive x.
21 * 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
22 * 3. cos(x) is approximated by a polynomial of degree 14 on
23 * [0,pi/4]
24 * 4 14
25 * cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
26 * where the remez error is
27 *
28 * | 2 4 6 8 10 12 14 | -58
29 * |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2
30 * | |
31 *
32 * 4 6 8 10 12 14
33 * 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then
34 * cos(x) = 1 - x*x/2 + r
35 * since cos(x+y) ~ cos(x) - sin(x)*y
36 * ~ cos(x) - x*y,
37 * a correction term is necessary in cos(x) and hence
38 * cos(x+y) = 1 - (x*x/2 - (r - x*y))
39 * For better accuracy when x > 0.3, let qx = |x|/4 with
40 * the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125.
41 * Then
42 * cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)).
43 * Note that 1-qx and (x*x/2-qx) is EXACT here, and the
44 * magnitude of the latter is at least a quarter of x*x/2,
45 * thus, reducing the rounding error in the subtraction.
46 */
47
48#include "math_libm.h"
49#include "math_private.h"
50
51static const double
52one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
53C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
54C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
55C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
56C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
57C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
58C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
59
60double attribute_hidden __kernel_cos(double x, double y)
61{
62 double a,hz,z,r,qx;
63 int32_t ix;
64 GET_HIGH_WORD(ix,x);
65 ix &= 0x7fffffff; /* ix = |x|'s high word*/
66 if(ix<0x3e400000) { /* if x < 2**27 */
67 if(((int)x)==0) return one; /* generate inexact */
68 }
69 z = x*x;
70 r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6)))));
71 if(ix < 0x3FD33333) /* if |x| < 0.3 */
72 return one - (0.5*z - (z*r - x*y));
73 else {
74 if(ix > 0x3fe90000) { /* x > 0.78125 */
75 qx = 0.28125;
76 } else {
77 INSERT_WORDS(qx,ix-0x00200000,0); /* x/4 */
78 }
79 hz = 0.5*z-qx;
80 a = one-qx;
81 return a - (hz - (z*r-x*y));
82 }
83}
diff --git a/contrib/SDL-3.2.8/src/libm/k_rem_pio2.c b/contrib/SDL-3.2.8/src/libm/k_rem_pio2.c
new file mode 100644
index 0000000..3dd5b2b
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/k_rem_pio2.c
@@ -0,0 +1,315 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/*
14 * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
15 * double x[],y[]; int e0,nx,prec; int ipio2[];
16 *
17 * __kernel_rem_pio2 return the last three digits of N with
18 * y = x - N*pi/2
19 * so that |y| < pi/2.
20 *
21 * The method is to compute the integer (mod 8) and fraction parts of
22 * (2/pi)*x without doing the full multiplication. In general we
23 * skip the part of the product that are known to be a huge integer (
24 * more accurately, = 0 mod 8 ). Thus the number of operations are
25 * independent of the exponent of the input.
26 *
27 * (2/pi) is represented by an array of 24-bit integers in ipio2[].
28 *
29 * Input parameters:
30 * x[] The input value (must be positive) is broken into nx
31 * pieces of 24-bit integers in double precision format.
32 * x[i] will be the i-th 24 bit of x. The scaled exponent
33 * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
34 * match x's up to 24 bits.
35 *
36 * Example of breaking a double positive z into x[0]+x[1]+x[2]:
37 * e0 = ilogb(z)-23
38 * z = scalbn(z,-e0)
39 * for i = 0,1,2
40 * x[i] = floor(z)
41 * z = (z-x[i])*2**24
42 *
43 *
44 * y[] ouput result in an array of double precision numbers.
45 * The dimension of y[] is:
46 * 24-bit precision 1
47 * 53-bit precision 2
48 * 64-bit precision 2
49 * 113-bit precision 3
50 * The actual value is the sum of them. Thus for 113-bit
51 * precison, one may have to do something like:
52 *
53 * long double t,w,r_head, r_tail;
54 * t = (long double)y[2] + (long double)y[1];
55 * w = (long double)y[0];
56 * r_head = t+w;
57 * r_tail = w - (r_head - t);
58 *
59 * e0 The exponent of x[0]
60 *
61 * nx dimension of x[]
62 *
63 * prec an integer indicating the precision:
64 * 0 24 bits (single)
65 * 1 53 bits (double)
66 * 2 64 bits (extended)
67 * 3 113 bits (quad)
68 *
69 * ipio2[]
70 * integer array, contains the (24*i)-th to (24*i+23)-th
71 * bit of 2/pi after binary point. The corresponding
72 * floating value is
73 *
74 * ipio2[i] * 2^(-24(i+1)).
75 *
76 * External function:
77 * double scalbn(), floor();
78 *
79 *
80 * Here is the description of some local variables:
81 *
82 * jk jk+1 is the initial number of terms of ipio2[] needed
83 * in the computation. The recommended value is 2,3,4,
84 * 6 for single, double, extended,and quad.
85 *
86 * jz local integer variable indicating the number of
87 * terms of ipio2[] used.
88 *
89 * jx nx - 1
90 *
91 * jv index for pointing to the suitable ipio2[] for the
92 * computation. In general, we want
93 * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
94 * is an integer. Thus
95 * e0-3-24*jv >= 0 or (e0-3)/24 >= jv
96 * Hence jv = max(0,(e0-3)/24).
97 *
98 * jp jp+1 is the number of terms in PIo2[] needed, jp = jk.
99 *
100 * q[] double array with integral value, representing the
101 * 24-bits chunk of the product of x and 2/pi.
102 *
103 * q0 the corresponding exponent of q[0]. Note that the
104 * exponent for q[i] would be q0-24*i.
105 *
106 * PIo2[] double precision array, obtained by cutting pi/2
107 * into 24 bits chunks.
108 *
109 * f[] ipio2[] in floating point
110 *
111 * iq[] integer array by breaking up q[] in 24-bits chunk.
112 *
113 * fq[] final product of x*(2/pi) in fq[0],..,fq[jk]
114 *
115 * ih integer. If >0 it indicates q[] is >= 0.5, hence
116 * it also indicates the *sign* of the result.
117 *
118 */
119
120
121/*
122 * Constants:
123 * The hexadecimal values are the intended ones for the following
124 * constants. The decimal values may be used, provided that the
125 * compiler will convert from decimal to binary accurately enough
126 * to produce the hexadecimal values shown.
127 */
128
129#include "math_libm.h"
130#include "math_private.h"
131
132
133static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
134
135static 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
146static const double
147zero = 0.0,
148one = 1.0,
149two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
150twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
151
152int32_t attribute_hidden __kernel_rem_pio2(const 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*/
162 SDL_assert(prec < SDL_arraysize(init_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 (iq) - (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 SDL_zero(fq);
275 for(i=jz;i>=0;i--) {
276 for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
277 fq[jz-i] = fw;
278 }
279
280 /* compress fq[] into y[] */
281 switch(prec) {
282 case 0:
283 fw = 0.0;
284 for (i=jz;i>=0;i--) fw += fq[i];
285 y[0] = (ih==0)? fw: -fw;
286 break;
287 case 1:
288 case 2:
289 fw = 0.0;
290 for (i=jz;i>=0;i--) fw += fq[i];
291 y[0] = (ih==0)? fw: -fw;
292 fw = fq[0]-fw;
293 for (i=1;i<=jz;i++) fw += fq[i];
294 y[1] = (ih==0)? fw: -fw;
295 break;
296 case 3: /* painful */
297 for (i=jz;i>0;i--) {
298 fw = fq[i-1]+fq[i];
299 fq[i] += fq[i-1]-fw;
300 fq[i-1] = fw;
301 }
302 for (i=jz;i>1;i--) {
303 fw = fq[i-1]+fq[i];
304 fq[i] += fq[i-1]-fw;
305 fq[i-1] = fw;
306 }
307 for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
308 if(ih==0) {
309 y[0] = fq[0]; y[1] = fq[1]; y[2] = fw;
310 } else {
311 y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
312 }
313 }
314 return n&7;
315}
diff --git a/contrib/SDL-3.2.8/src/libm/k_sin.c b/contrib/SDL-3.2.8/src/libm/k_sin.c
new file mode 100644
index 0000000..b7596ae
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/k_sin.c
@@ -0,0 +1,66 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/* __kernel_sin( x, y, iy)
14 * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
15 * Input x is assumed to be bounded by ~pi/4 in magnitude.
16 * Input y is the tail of x.
17 * Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
18 *
19 * Algorithm
20 * 1. Since sin(-x) = -sin(x), we need only to consider positive x.
21 * 2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0.
22 * 3. sin(x) is approximated by a polynomial of degree 13 on
23 * [0,pi/4]
24 * 3 13
25 * sin(x) ~ x + S1*x + ... + S6*x
26 * where
27 *
28 * |sin(x) 2 4 6 8 10 12 | -58
29 * |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2
30 * | x |
31 *
32 * 4. sin(x+y) = sin(x) + sin'(x')*y
33 * ~ sin(x) + (1-x*x/2)*y
34 * For better accuracy, let
35 * 3 2 2 2 2
36 * r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
37 * then 3 2
38 * sin(x) = x + (S1*x + (x *(r-y/2)+y))
39 */
40
41#include "math_libm.h"
42#include "math_private.h"
43
44static const double
45half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
46S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
47S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
48S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
49S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
50S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
51S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
52
53double attribute_hidden __kernel_sin(double x, double y, int iy)
54{
55 double z,r,v;
56 int32_t ix;
57 GET_HIGH_WORD(ix,x);
58 ix &= 0x7fffffff; /* high word of x */
59 if(ix<0x3e400000) /* |x| < 2**-27 */
60 {if((int)x==0) return x;} /* generate inexact */
61 z = x*x;
62 v = z*x;
63 r = S2+z*(S3+z*(S4+z*(S5+z*S6)));
64 if(iy==0) return x+v*(S1+z*r);
65 else return x-((z*(half*y-v*r)-y)-v*S1);
66}
diff --git a/contrib/SDL-3.2.8/src/libm/k_tan.c b/contrib/SDL-3.2.8/src/libm/k_tan.c
new file mode 100644
index 0000000..1405e03
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/k_tan.c
@@ -0,0 +1,119 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/* __kernel_tan( x, y, k )
14 * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
15 * Input x is assumed to be bounded by ~pi/4 in magnitude.
16 * Input y is the tail of x.
17 * Input k indicates whether tan (if k=1) or
18 * -1/tan (if k= -1) is returned.
19 *
20 * Algorithm
21 * 1. Since tan(-x) = -tan(x), we need only to consider positive x.
22 * 2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0.
23 * 3. tan(x) is approximated by a odd polynomial of degree 27 on
24 * [0,0.67434]
25 * 3 27
26 * tan(x) ~ x + T1*x + ... + T13*x
27 * where
28 *
29 * |tan(x) 2 4 26 | -59.2
30 * |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2
31 * | x |
32 *
33 * Note: tan(x+y) = tan(x) + tan'(x)*y
34 * ~ tan(x) + (1+x*x)*y
35 * Therefore, for better accuracy in computing tan(x+y), let
36 * 3 2 2 2 2
37 * r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
38 * then
39 * 3 2
40 * tan(x+y) = x + (T1*x + (x *(r+y)+y))
41 *
42 * 4. For x in [0.67434,pi/4], let y = pi/4 - x, then
43 * tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
44 * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
45 */
46
47#include "math_libm.h"
48#include "math_private.h"
49
50static const double
51one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
52pio4 = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */
53pio4lo= 3.06161699786838301793e-17, /* 0x3C81A626, 0x33145C07 */
54T[] = {
55 3.33333333333334091986e-01, /* 0x3FD55555, 0x55555563 */
56 1.33333333333201242699e-01, /* 0x3FC11111, 0x1110FE7A */
57 5.39682539762260521377e-02, /* 0x3FABA1BA, 0x1BB341FE */
58 2.18694882948595424599e-02, /* 0x3F9664F4, 0x8406D637 */
59 8.86323982359930005737e-03, /* 0x3F8226E3, 0xE96E8493 */
60 3.59207910759131235356e-03, /* 0x3F6D6D22, 0xC9560328 */
61 1.45620945432529025516e-03, /* 0x3F57DBC8, 0xFEE08315 */
62 5.88041240820264096874e-04, /* 0x3F4344D8, 0xF2F26501 */
63 2.46463134818469906812e-04, /* 0x3F3026F7, 0x1A8D1068 */
64 7.81794442939557092300e-05, /* 0x3F147E88, 0xA03792A6 */
65 7.14072491382608190305e-05, /* 0x3F12B80F, 0x32F0A7E9 */
66 -1.85586374855275456654e-05, /* 0xBEF375CB, 0xDB605373 */
67 2.59073051863633712884e-05, /* 0x3EFB2A70, 0x74BF7AD4 */
68};
69
70double attribute_hidden __kernel_tan(double x, double y, int iy)
71{
72 double z,r,v,w,s;
73 int32_t ix,hx;
74 GET_HIGH_WORD(hx,x);
75 ix = hx&0x7fffffff; /* high word of |x| */
76 if(ix<0x3e300000) /* x < 2**-28 */
77 {if((int)x==0) { /* generate inexact */
78 u_int32_t low;
79 GET_LOW_WORD(low,x);
80 if(((ix|low)|(iy+1))==0) return one/fabs(x);
81 else return (iy==1)? x: -one/x;
82 }
83 }
84 if(ix>=0x3FE59428) { /* |x|>=0.6744 */
85 if(hx<0) {x = -x; y = -y;}
86 z = pio4-x;
87 w = pio4lo-y;
88 x = z+w; y = 0.0;
89 }
90 z = x*x;
91 w = z*z;
92 /* Break x^5*(T[1]+x^2*T[2]+...) into
93 * x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
94 * x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
95 */
96 r = T[1]+w*(T[3]+w*(T[5]+w*(T[7]+w*(T[9]+w*T[11]))));
97 v = z*(T[2]+w*(T[4]+w*(T[6]+w*(T[8]+w*(T[10]+w*T[12])))));
98 s = z*x;
99 r = y + z*(s*(r+v)+y);
100 r += T[0]*s;
101 w = x+r;
102 if(ix>=0x3FE59428) {
103 v = (double)iy;
104 return (double)(1-((hx>>30)&2))*(v-2.0*(x-(w*w/(w+v)-r)));
105 }
106 if(iy==1) return w;
107 else { /* if allow error up to 2 ulp,
108 simply return -1.0/(x+r) here */
109 /* compute -1.0/(x+r) accurately */
110 double a,t;
111 z = w;
112 SET_LOW_WORD(z,0);
113 v = r-(z - x); /* z+v = r+x */
114 t = a = -1.0/w; /* a = -1.0/w */
115 SET_LOW_WORD(t,0);
116 s = 1.0+t*z;
117 return t+a*(s+t*v);
118 }
119}
diff --git a/contrib/SDL-3.2.8/src/libm/math_libm.h b/contrib/SDL-3.2.8/src/libm/math_libm.h
new file mode 100644
index 0000000..b7b1614
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/math_libm.h
@@ -0,0 +1,50 @@
1/*
2 Simple DirectMedia Layer
3 Copyright (C) 1997-2025 Sam Lantinga <slouken@libsdl.org>
4
5 This software is provided 'as-is', without any express or implied
6 warranty. In no event will the authors be held liable for any damages
7 arising from the use of this software.
8
9 Permission is granted to anyone to use this software for any purpose,
10 including commercial applications, and to alter it and redistribute it
11 freely, subject to the following restrictions:
12
13 1. The origin of this software must not be misrepresented; you must not
14 claim that you wrote the original software. If you use this software
15 in a product, an acknowledgment in the product documentation would be
16 appreciated but is not required.
17 2. Altered source versions must be plainly marked as such, and must not be
18 misrepresented as being the original software.
19 3. This notice may not be removed or altered from any source distribution.
20*/
21
22#ifndef math_libm_h_
23#define math_libm_h_
24
25#include "SDL_internal.h"
26
27/* Math routines from uClibc: http://www.uclibc.org */
28
29extern double SDL_uclibc_atan(double x);
30extern double SDL_uclibc_atan2(double y, double x);
31extern double SDL_uclibc_copysign(double x, double y);
32extern double SDL_uclibc_cos(double x);
33extern double SDL_uclibc_exp(double x);
34extern double SDL_uclibc_fabs(double x);
35extern double SDL_uclibc_floor(double x);
36extern double SDL_uclibc_fmod(double x, double y);
37extern int SDL_uclibc_isinf(double x);
38extern int SDL_uclibc_isinff(float x);
39extern int SDL_uclibc_isnan(double x);
40extern int SDL_uclibc_isnanf(float x);
41extern double SDL_uclibc_log(double x);
42extern double SDL_uclibc_log10(double x);
43extern double SDL_uclibc_modf(double x, double *y);
44extern double SDL_uclibc_pow(double x, double y);
45extern double SDL_uclibc_scalbn(double x, int n);
46extern double SDL_uclibc_sin(double x);
47extern double SDL_uclibc_sqrt(double x);
48extern double SDL_uclibc_tan(double x);
49
50#endif /* math_libm_h_ */
diff --git a/contrib/SDL-3.2.8/src/libm/math_private.h b/contrib/SDL-3.2.8/src/libm/math_private.h
new file mode 100644
index 0000000..b7d60e4
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/math_private.h
@@ -0,0 +1,234 @@
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 * from: @(#)fdlibm.h 5.1 93/09/24
14 * $Id: math_private.h,v 1.3 2004/02/09 07:10:38 andersen Exp $
15 */
16
17#ifndef _MATH_PRIVATE_H_
18#define _MATH_PRIVATE_H_
19
20/* #include <endian.h> */
21/* #include <sys/types.h> */
22
23#define _IEEE_LIBM
24#define attribute_hidden
25#define libm_hidden_proto(x)
26#define libm_hidden_def(x)
27#define strong_alias(x, y)
28#define weak_alias(x, y)
29
30#if !defined(SDL_PLATFORM_HAIKU) && !defined(SDL_PLATFORM_PSP) && !defined(SDL_PLATFORM_3DS) && !defined(SDL_PLATFORM_PS2) /* already defined in a system header. */
31typedef unsigned int u_int32_t;
32#endif
33
34#define atan SDL_uclibc_atan
35#define __ieee754_atan2 SDL_uclibc_atan2
36#define copysign SDL_uclibc_copysign
37#define cos SDL_uclibc_cos
38#define __ieee754_exp SDL_uclibc_exp
39#define fabs SDL_uclibc_fabs
40#define floor SDL_uclibc_floor
41#define __ieee754_fmod SDL_uclibc_fmod
42#undef __isinf
43#define __isinf SDL_uclibc_isinf
44#undef __isinff
45#define __isinff SDL_uclibc_isinff
46#undef __isnan
47#define __isnan SDL_uclibc_isnan
48#undef __isnanf
49#define __isnanf SDL_uclibc_isnanf
50#define __ieee754_log SDL_uclibc_log
51#define __ieee754_log10 SDL_uclibc_log10
52#define modf SDL_uclibc_modf
53#define __ieee754_pow SDL_uclibc_pow
54#define scalbln SDL_uclibc_scalbln
55#define scalbn SDL_uclibc_scalbn
56#define sin SDL_uclibc_sin
57#define __ieee754_sqrt SDL_uclibc_sqrt
58#define tan SDL_uclibc_tan
59
60/* The original fdlibm code used statements like:
61 n0 = ((*(int*)&one)>>29)^1; * index of high word *
62 ix0 = *(n0+(int*)&x); * high word of x *
63 ix1 = *((1-n0)+(int*)&x); * low word of x *
64 to dig two 32 bit words out of the 64 bit IEEE floating point
65 value. That is non-ANSI, and, moreover, the gcc instruction
66 scheduler gets it wrong. We instead use the following macros.
67 Unlike the original code, we determine the endianness at compile
68 time, not at run time; I don't see much benefit to selecting
69 endianness at run time. */
70
71/* A union which permits us to convert between a double and two 32 bit
72 ints. */
73
74/*
75 * Math on arm is special:
76 * For FPA, float words are always big-endian.
77 * For VFP, floats words follow the memory system mode.
78 * For Maverick, float words are always little-endian.
79 */
80
81#if (SDL_FLOATWORDORDER == SDL_BIG_ENDIAN)
82
83typedef union
84{
85 double value;
86 struct
87 {
88 u_int32_t msw;
89 u_int32_t lsw;
90 } parts;
91} ieee_double_shape_type;
92
93#else
94
95typedef union
96{
97 double value;
98 struct
99 {
100 u_int32_t lsw;
101 u_int32_t msw;
102 } parts;
103} ieee_double_shape_type;
104
105#endif
106
107/* Get two 32 bit ints from a double. */
108
109#define EXTRACT_WORDS(ix0,ix1,d) \
110do { \
111 ieee_double_shape_type ew_u; \
112 ew_u.value = (d); \
113 (ix0) = ew_u.parts.msw; \
114 (ix1) = ew_u.parts.lsw; \
115} while (0)
116
117/* Get the more significant 32 bit int from a double. */
118
119#define GET_HIGH_WORD(i,d) \
120do { \
121 ieee_double_shape_type gh_u; \
122 gh_u.value = (d); \
123 (i) = gh_u.parts.msw; \
124} while (0)
125
126/* Get the less significant 32 bit int from a double. */
127
128#define GET_LOW_WORD(i,d) \
129do { \
130 ieee_double_shape_type gl_u; \
131 gl_u.value = (d); \
132 (i) = gl_u.parts.lsw; \
133} while (0)
134
135/* Set a double from two 32 bit ints. */
136
137#define INSERT_WORDS(d,ix0,ix1) \
138do { \
139 ieee_double_shape_type iw_u; \
140 iw_u.parts.msw = (ix0); \
141 iw_u.parts.lsw = (ix1); \
142 (d) = iw_u.value; \
143} while (0)
144
145/* Set the more significant 32 bits of a double from an int. */
146
147#define SET_HIGH_WORD(d,v) \
148do { \
149 ieee_double_shape_type sh_u; \
150 sh_u.value = (d); \
151 sh_u.parts.msw = (v); \
152 (d) = sh_u.value; \
153} while (0)
154
155/* Set the less significant 32 bits of a double from an int. */
156
157#define SET_LOW_WORD(d,v) \
158do { \
159 ieee_double_shape_type sl_u; \
160 sl_u.value = (d); \
161 sl_u.parts.lsw = (v); \
162 (d) = sl_u.value; \
163} while (0)
164
165/* A union which permits us to convert between a float and a 32 bit
166 int. */
167
168typedef union
169{
170 float value;
171 u_int32_t word;
172} ieee_float_shape_type;
173
174/* Get a 32 bit int from a float. */
175
176#define GET_FLOAT_WORD(i,d) \
177do { \
178 ieee_float_shape_type gf_u; \
179 gf_u.value = (d); \
180 (i) = gf_u.word; \
181} while (0)
182
183/* Set a float from a 32 bit int. */
184
185#define SET_FLOAT_WORD(d,i) \
186do { \
187 ieee_float_shape_type sf_u; \
188 sf_u.word = (i); \
189 (d) = sf_u.value; \
190} while (0)
191
192/* ieee style elementary functions */
193extern double __ieee754_sqrt(double) attribute_hidden;
194extern double __ieee754_acos(double) attribute_hidden;
195extern double __ieee754_acosh(double) attribute_hidden;
196extern double __ieee754_log(double) attribute_hidden;
197extern double __ieee754_atanh(double) attribute_hidden;
198extern double __ieee754_asin(double) attribute_hidden;
199extern double __ieee754_atan2(double, double) attribute_hidden;
200extern double __ieee754_exp(double) attribute_hidden;
201extern double __ieee754_cosh(double) attribute_hidden;
202extern double __ieee754_fmod(double, double) attribute_hidden;
203extern double __ieee754_pow(double, double) attribute_hidden;
204extern double __ieee754_lgamma_r(double, int *) attribute_hidden;
205extern double __ieee754_gamma_r(double, int *) attribute_hidden;
206extern double __ieee754_lgamma(double) attribute_hidden;
207extern double __ieee754_gamma(double) attribute_hidden;
208extern double __ieee754_log10(double) attribute_hidden;
209extern double __ieee754_sinh(double) attribute_hidden;
210extern double __ieee754_hypot(double, double) attribute_hidden;
211extern double __ieee754_j0(double) attribute_hidden;
212extern double __ieee754_j1(double) attribute_hidden;
213extern double __ieee754_y0(double) attribute_hidden;
214extern double __ieee754_y1(double) attribute_hidden;
215extern double __ieee754_jn(int, double) attribute_hidden;
216extern double __ieee754_yn(int, double) attribute_hidden;
217extern double __ieee754_remainder(double, double) attribute_hidden;
218extern int32_t __ieee754_rem_pio2(double, double *) attribute_hidden;
219#if defined(_SCALB_INT)
220extern double __ieee754_scalb(double, int) attribute_hidden;
221#else
222extern double __ieee754_scalb(double, double) attribute_hidden;
223#endif
224
225/* fdlibm kernel function */
226#ifndef _IEEE_LIBM
227extern double __kernel_standard(double, double, int) attribute_hidden;
228#endif
229extern double __kernel_sin(double, double, int) attribute_hidden;
230extern double __kernel_cos(double, double) attribute_hidden;
231extern double __kernel_tan(double, double, int) attribute_hidden;
232extern int32_t __kernel_rem_pio2(const double *, double *, int, int, const unsigned int, const int32_t *) attribute_hidden;
233
234#endif /* _MATH_PRIVATE_H_ */
diff --git a/contrib/SDL-3.2.8/src/libm/s_atan.c b/contrib/SDL-3.2.8/src/libm/s_atan.c
new file mode 100644
index 0000000..ce429d2
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/s_atan.c
@@ -0,0 +1,119 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/* atan(x)
14 * Method
15 * 1. Reduce x to positive by atan(x) = -atan(-x).
16 * 2. According to the integer k=4t+0.25 chopped, t=x, the argument
17 * is further reduced to one of the following intervals and the
18 * arctangent of t is evaluated by the corresponding formula:
19 *
20 * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
21 * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
22 * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
23 * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
24 * [39/16,INF] atan(x) = atan(INF) + atan( -1/t )
25 *
26 * Constants:
27 * The hexadecimal values are the intended ones for the following
28 * constants. The decimal values may be used, provided that the
29 * compiler will convert from decimal to binary accurately enough
30 * to produce the hexadecimal values shown.
31 */
32
33#include "math_libm.h"
34#include "math_private.h"
35
36static const double atanhi[] = {
37 4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */
38 7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */
39 9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */
40 1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */
41};
42
43static const double atanlo[] = {
44 2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */
45 3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */
46 1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */
47 6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */
48};
49
50static const double aT[] = {
51 3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */
52 -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */
53 1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */
54 -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */
55 9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */
56 -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */
57 6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */
58 -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */
59 4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */
60 -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */
61 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
62};
63
64#ifdef __WATCOMC__ /* Watcom defines huge=__huge */
65#undef huge
66#endif
67
68static const double
69one = 1.0,
70huge = 1.0e300;
71
72double atan(double x)
73{
74 double w,s1,s2,z;
75 int32_t ix,hx,id;
76
77 GET_HIGH_WORD(hx,x);
78 ix = hx&0x7fffffff;
79 if(ix>=0x44100000) { /* if |x| >= 2^66 */
80 u_int32_t low;
81 GET_LOW_WORD(low,x);
82 if(ix>0x7ff00000||
83 (ix==0x7ff00000&&(low!=0)))
84 return x+x; /* NaN */
85 if(hx>0) return atanhi[3]+atanlo[3];
86 else return -atanhi[3]-atanlo[3];
87 } if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
88 if (ix < 0x3e200000) { /* |x| < 2^-29 */
89 if(huge+x>one) return x; /* raise inexact */
90 }
91 id = -1;
92 } else {
93 x = fabs(x);
94 if (ix < 0x3ff30000) { /* |x| < 1.1875 */
95 if (ix < 0x3fe60000) { /* 7/16 <=|x|<11/16 */
96 id = 0; x = (2.0*x-one)/(2.0+x);
97 } else { /* 11/16<=|x|< 19/16 */
98 id = 1; x = (x-one)/(x+one);
99 }
100 } else {
101 if (ix < 0x40038000) { /* |x| < 2.4375 */
102 id = 2; x = (x-1.5)/(one+1.5*x);
103 } else { /* 2.4375 <= |x| < 2^66 */
104 id = 3; x = -1.0/x;
105 }
106 }}
107 /* end of argument reduction */
108 z = x*x;
109 w = z*z;
110 /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
111 s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10])))));
112 s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9]))));
113 if (id<0) return x - x*(s1+s2);
114 else {
115 z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
116 return (hx<0)? -z:z;
117 }
118}
119libm_hidden_def(atan)
diff --git a/contrib/SDL-3.2.8/src/libm/s_copysign.c b/contrib/SDL-3.2.8/src/libm/s_copysign.c
new file mode 100644
index 0000000..630f477
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/s_copysign.c
@@ -0,0 +1,30 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/*
14 * copysign(double x, double y)
15 * copysign(x,y) returns a value with the magnitude of x and
16 * with the sign bit of y.
17 */
18
19#include "math_libm.h"
20#include "math_private.h"
21
22double copysign(double x, double y)
23{
24 u_int32_t hx,hy;
25 GET_HIGH_WORD(hx,x);
26 GET_HIGH_WORD(hy,y);
27 SET_HIGH_WORD(x,(hx&0x7fffffff)|(hy&0x80000000));
28 return x;
29}
30libm_hidden_def(copysign)
diff --git a/contrib/SDL-3.2.8/src/libm/s_cos.c b/contrib/SDL-3.2.8/src/libm/s_cos.c
new file mode 100644
index 0000000..ef85e71
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/s_cos.c
@@ -0,0 +1,74 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/* cos(x)
14 * Return cosine function of x.
15 *
16 * kernel function:
17 * __kernel_sin ... sine function on [-pi/4,pi/4]
18 * __kernel_cos ... cosine function on [-pi/4,pi/4]
19 * __ieee754_rem_pio2 ... argument reduction routine
20 *
21 * Method.
22 * Let S,C and T denote the sin, cos and tan respectively on
23 * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
24 * in [-pi/4 , +pi/4], and let n = k mod 4.
25 * We have
26 *
27 * n sin(x) cos(x) tan(x)
28 * ----------------------------------------------------------
29 * 0 S C T
30 * 1 C -S -1/T
31 * 2 -S -C T
32 * 3 -C S -1/T
33 * ----------------------------------------------------------
34 *
35 * Special cases:
36 * Let trig be any of sin, cos, or tan.
37 * trig(+-INF) is NaN, with signals;
38 * trig(NaN) is that NaN;
39 *
40 * Accuracy:
41 * TRIG(x) returns trig(x) nearly rounded
42 */
43
44#include "math_libm.h"
45#include "math_private.h"
46
47double cos(double x)
48{
49 double y[2],z=0.0;
50 int32_t n, ix;
51
52 /* High word of x. */
53 GET_HIGH_WORD(ix,x);
54
55 /* |x| ~< pi/4 */
56 ix &= 0x7fffffff;
57 if(ix <= 0x3fe921fb) return __kernel_cos(x,z);
58
59 /* cos(Inf or NaN) is NaN */
60 else if (ix>=0x7ff00000) return x-x;
61
62 /* argument reduction needed */
63 else {
64 n = __ieee754_rem_pio2(x,y);
65 switch(n&3) {
66 case 0: return __kernel_cos(y[0],y[1]);
67 case 1: return -__kernel_sin(y[0],y[1],1);
68 case 2: return -__kernel_cos(y[0],y[1]);
69 default:
70 return __kernel_sin(y[0],y[1],1);
71 }
72 }
73}
74libm_hidden_def(cos)
diff --git a/contrib/SDL-3.2.8/src/libm/s_fabs.c b/contrib/SDL-3.2.8/src/libm/s_fabs.c
new file mode 100644
index 0000000..782b6f7
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/s_fabs.c
@@ -0,0 +1,30 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/*
14 * fabs(x) returns the absolute value of x.
15 */
16
17/*#include <features.h>*/
18/* Prevent math.h from defining a colliding inline */
19#undef __USE_EXTERN_INLINES
20#include "math_libm.h"
21#include "math_private.h"
22
23double fabs(double x)
24{
25 u_int32_t high;
26 GET_HIGH_WORD(high,x);
27 SET_HIGH_WORD(x,high&0x7fffffff);
28 return x;
29}
30libm_hidden_def(fabs)
diff --git a/contrib/SDL-3.2.8/src/libm/s_floor.c b/contrib/SDL-3.2.8/src/libm/s_floor.c
new file mode 100644
index 0000000..4809af1
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/s_floor.c
@@ -0,0 +1,76 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/*
14 * floor(x)
15 * Return x rounded toward -inf to integral value
16 * Method:
17 * Bit twiddling.
18 * Exception:
19 * Inexact flag raised if x not equal to floor(x).
20 */
21
22/*#include <features.h>*/
23/* Prevent math.h from defining a colliding inline */
24#undef __USE_EXTERN_INLINES
25#include "math_libm.h"
26#include "math_private.h"
27
28#ifdef __WATCOMC__ /* Watcom defines huge=__huge */
29#undef huge
30#endif
31
32static const double huge = 1.0e300;
33
34double floor(double x)
35{
36 int32_t i0,i1,j0;
37 u_int32_t i,j;
38 EXTRACT_WORDS(i0,i1,x);
39 j0 = ((i0>>20)&0x7ff)-0x3ff;
40 if(j0<20) {
41 if(j0<0) { /* raise inexact if x != 0 */
42 if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
43 if(i0>=0) {i0=i1=0;}
44 else if(((i0&0x7fffffff)|i1)!=0)
45 { i0=0xbff00000;i1=0;}
46 }
47 } else {
48 i = (0x000fffff)>>j0;
49 if(((i0&i)|i1)==0) return x; /* x is integral */
50 if(huge+x>0.0) { /* raise inexact flag */
51 if(i0<0) i0 += (0x00100000)>>j0;
52 i0 &= (~i); i1=0;
53 }
54 }
55 } else if (j0>51) {
56 if(j0==0x400) return x+x; /* inf or NaN */
57 else return x; /* x is integral */
58 } else {
59 i = ((u_int32_t)(0xffffffff))>>(j0-20);
60 if((i1&i)==0) return x; /* x is integral */
61 if(huge+x>0.0) { /* raise inexact flag */
62 if(i0<0) {
63 if(j0==20) i0+=1;
64 else {
65 j = i1+(1<<(52-j0));
66 if(j<(u_int32_t)i1) i0 +=1 ; /* got a carry */
67 i1=j;
68 }
69 }
70 i1 &= (~i);
71 }
72 }
73 INSERT_WORDS(x,i0,i1);
74 return x;
75}
76libm_hidden_def(floor)
diff --git a/contrib/SDL-3.2.8/src/libm/s_isinf.c b/contrib/SDL-3.2.8/src/libm/s_isinf.c
new file mode 100644
index 0000000..9486b05
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/s_isinf.c
@@ -0,0 +1,24 @@
1#include "SDL_internal.h"
2/*
3 * Written by J.T. Conklin <jtc@netbsd.org>.
4 * Changed to return -1 for -Inf by Ulrich Drepper <drepper@cygnus.com>.
5 * Public domain.
6 */
7
8/*
9 * isinf(x) returns 1 is x is inf, -1 if x is -inf, else 0;
10 * no branching!
11 */
12
13#include "math.h"
14#include "math_private.h"
15
16int __isinf(double x)
17{
18 int32_t hx,lx;
19 EXTRACT_WORDS(hx,lx,x);
20 lx |= (hx & 0x7fffffff) ^ 0x7ff00000;
21 lx |= -lx;
22 return ~(lx >> 31) & (hx >> 30);
23}
24libm_hidden_def(__isinf)
diff --git a/contrib/SDL-3.2.8/src/libm/s_isinff.c b/contrib/SDL-3.2.8/src/libm/s_isinff.c
new file mode 100644
index 0000000..184c9aa
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/s_isinff.c
@@ -0,0 +1,24 @@
1#include "SDL_internal.h"
2/*
3 * Written by J.T. Conklin <jtc@netbsd.org>.
4 * Public domain.
5 */
6
7/*
8 * isinff(x) returns 1 if x is inf, -1 if x is -inf, else 0;
9 * no branching!
10 */
11
12#include "math.h"
13#include "math_private.h"
14
15int __isinff (float x)
16{
17 int32_t ix,t;
18 GET_FLOAT_WORD(ix,x);
19 t = ix & 0x7fffffff;
20 t ^= 0x7f800000;
21 t |= -t;
22 return ~(t >> 31) & (ix >> 30);
23}
24libm_hidden_def(__isinff)
diff --git a/contrib/SDL-3.2.8/src/libm/s_isnan.c b/contrib/SDL-3.2.8/src/libm/s_isnan.c
new file mode 100644
index 0000000..4831086
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/s_isnan.c
@@ -0,0 +1,31 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/*
14 * isnan(x) returns 1 is x is nan, else 0;
15 * no branching!
16 */
17
18#include "math.h"
19#include "math_private.h"
20
21int __isnan(double x)
22{
23 int32_t hx,lx;
24 EXTRACT_WORDS(hx,lx,x);
25 hx &= 0x7fffffff;
26 hx |= (u_int32_t)(lx|(-lx))>>31;
27 hx = 0x7ff00000 - hx;
28 return (int)(((u_int32_t)hx)>>31);
29}
30weak_alias(__isnan, isnan)
31libm_hidden_def(__isnan)
diff --git a/contrib/SDL-3.2.8/src/libm/s_isnanf.c b/contrib/SDL-3.2.8/src/libm/s_isnanf.c
new file mode 100644
index 0000000..1cb308c
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/s_isnanf.c
@@ -0,0 +1,33 @@
1#include "SDL_internal.h"
2/* s_isnanf.c -- float version of s_isnan.c.
3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
4 */
5
6/*
7 * ====================================================
8 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
9 *
10 * Developed at SunPro, a Sun Microsystems, Inc. business.
11 * Permission to use, copy, modify, and distribute this
12 * software is freely granted, provided that this notice
13 * is preserved.
14 * ====================================================
15 */
16
17/*
18 * isnanf(x) returns 1 is x is nan, else 0;
19 * no branching!
20 */
21
22#include "math.h"
23#include "math_private.h"
24
25int __isnanf(float x)
26{
27 int32_t ix;
28 GET_FLOAT_WORD(ix,x);
29 ix &= 0x7fffffff;
30 ix = 0x7f800000 - ix;
31 return (int)(((u_int32_t)(ix))>>31);
32}
33libm_hidden_def(__isnanf)
diff --git a/contrib/SDL-3.2.8/src/libm/s_modf.c b/contrib/SDL-3.2.8/src/libm/s_modf.c
new file mode 100644
index 0000000..55f83ba
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/s_modf.c
@@ -0,0 +1,68 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/*
14 * modf(double x, double *iptr)
15 * return fraction part of x, and return x's integral part in *iptr.
16 * Method:
17 * Bit twiddling.
18 *
19 * Exception:
20 * No exception.
21 */
22
23#include "math_libm.h"
24#include "math_private.h"
25
26static const double one = 1.0;
27
28double modf(double x, double *iptr)
29{
30 int32_t i0,i1,_j0;
31 u_int32_t i;
32 EXTRACT_WORDS(i0,i1,x);
33 _j0 = ((i0>>20)&0x7ff)-0x3ff; /* exponent of x */
34 if(_j0<20) { /* integer part in high x */
35 if(_j0<0) { /* |x|<1 */
36 INSERT_WORDS(*iptr,i0&0x80000000,0); /* *iptr = +-0 */
37 return x;
38 } else {
39 i = (0x000fffff)>>_j0;
40 if(((i0&i)|i1)==0) { /* x is integral */
41 *iptr = x;
42 INSERT_WORDS(x,i0&0x80000000,0); /* return +-0 */
43 return x;
44 } else {
45 INSERT_WORDS(*iptr,i0&(~i),0);
46 return x - *iptr;
47 }
48 }
49 } else if (_j0>51) { /* no fraction part */
50 *iptr = x*one;
51 /* We must handle NaNs separately. */
52 if (_j0 == 0x400 && ((i0 & 0xfffff) | i1))
53 return x*one;
54 INSERT_WORDS(x,i0&0x80000000,0); /* return +-0 */
55 return x;
56 } else { /* fraction part in low x */
57 i = ((u_int32_t)(0xffffffff))>>(_j0-20);
58 if((i1&i)==0) { /* x is integral */
59 *iptr = x;
60 INSERT_WORDS(x,i0&0x80000000,0); /* return +-0 */
61 return x;
62 } else {
63 INSERT_WORDS(*iptr,i0,i1&(~i));
64 return x - *iptr;
65 }
66 }
67}
68libm_hidden_def(modf)
diff --git a/contrib/SDL-3.2.8/src/libm/s_scalbn.c b/contrib/SDL-3.2.8/src/libm/s_scalbn.c
new file mode 100644
index 0000000..b3a0604
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/s_scalbn.c
@@ -0,0 +1,74 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/*
14 * scalbln(double x, long n)
15 * scalbln(x,n) returns x * 2**n computed by exponent
16 * manipulation rather than by actually performing an
17 * exponentiation or a multiplication.
18 */
19
20#include "math_libm.h"
21#include "math_private.h"
22#include <limits.h>
23
24#ifdef __WATCOMC__ /* Watcom defines huge=__huge */
25#undef huge
26#endif
27
28static const double
29two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
30twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
31huge = 1.0e+300,
32tiny = 1.0e-300;
33
34double scalbln(double x, long n)
35{
36 int32_t k, hx, lx;
37
38 EXTRACT_WORDS(hx, lx, x);
39 k = (hx & 0x7ff00000) >> 20; /* extract exponent */
40 if (k == 0) { /* 0 or subnormal x */
41 if ((lx | (hx & 0x7fffffff)) == 0)
42 return x; /* +-0 */
43 x *= two54;
44 GET_HIGH_WORD(hx, x);
45 k = ((hx & 0x7ff00000) >> 20) - 54;
46 }
47 if (k == 0x7ff)
48 return x + x; /* NaN or Inf */
49 k = (int32_t)(k + n);
50 if (k > 0x7fe)
51 return huge * copysign(huge, x); /* overflow */
52 if (n < -50000)
53 return tiny * copysign(tiny, x); /* underflow */
54 if (k > 0) { /* normal result */
55 SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
56 return x;
57 }
58 if (k <= -54) {
59 if (n > 50000) /* in case integer overflow in n+k */
60 return huge * copysign(huge, x); /* overflow */
61 return tiny * copysign(tiny, x); /* underflow */
62 }
63 k += 54; /* subnormal result */
64 SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
65 return x * twom54;
66}
67libm_hidden_def(scalbln)
68
69
70double scalbn(double x, int n)
71{
72 return scalbln(x, n);
73}
74libm_hidden_def(scalbn)
diff --git a/contrib/SDL-3.2.8/src/libm/s_sin.c b/contrib/SDL-3.2.8/src/libm/s_sin.c
new file mode 100644
index 0000000..511bc21
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/s_sin.c
@@ -0,0 +1,74 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/* sin(x)
14 * Return sine function of x.
15 *
16 * kernel function:
17 * __kernel_sin ... sine function on [-pi/4,pi/4]
18 * __kernel_cos ... cose function on [-pi/4,pi/4]
19 * __ieee754_rem_pio2 ... argument reduction routine
20 *
21 * Method.
22 * Let S,C and T denote the sin, cos and tan respectively on
23 * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
24 * in [-pi/4 , +pi/4], and let n = k mod 4.
25 * We have
26 *
27 * n sin(x) cos(x) tan(x)
28 * ----------------------------------------------------------
29 * 0 S C T
30 * 1 C -S -1/T
31 * 2 -S -C T
32 * 3 -C S -1/T
33 * ----------------------------------------------------------
34 *
35 * Special cases:
36 * Let trig be any of sin, cos, or tan.
37 * trig(+-INF) is NaN, with signals;
38 * trig(NaN) is that NaN;
39 *
40 * Accuracy:
41 * TRIG(x) returns trig(x) nearly rounded
42 */
43
44#include "math_libm.h"
45#include "math_private.h"
46
47double sin(double x)
48{
49 double y[2],z=0.0;
50 int32_t n, ix;
51
52 /* High word of x. */
53 GET_HIGH_WORD(ix,x);
54
55 /* |x| ~< pi/4 */
56 ix &= 0x7fffffff;
57 if(ix <= 0x3fe921fb) return __kernel_sin(x,z,0);
58
59 /* sin(Inf or NaN) is NaN */
60 else if (ix>=0x7ff00000) return x-x;
61
62 /* argument reduction needed */
63 else {
64 n = __ieee754_rem_pio2(x,y);
65 switch(n&3) {
66 case 0: return __kernel_sin(y[0],y[1],1);
67 case 1: return __kernel_cos(y[0],y[1]);
68 case 2: return -__kernel_sin(y[0],y[1],1);
69 default:
70 return -__kernel_cos(y[0],y[1]);
71 }
72 }
73}
74libm_hidden_def(sin)
diff --git a/contrib/SDL-3.2.8/src/libm/s_tan.c b/contrib/SDL-3.2.8/src/libm/s_tan.c
new file mode 100644
index 0000000..dc10a1e
--- /dev/null
+++ b/contrib/SDL-3.2.8/src/libm/s_tan.c
@@ -0,0 +1,68 @@
1#include "SDL_internal.h"
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13/* tan(x)
14 * Return tangent function of x.
15 *
16 * kernel function:
17 * __kernel_tan ... tangent function on [-pi/4,pi/4]
18 * __ieee754_rem_pio2 ... argument reduction routine
19 *
20 * Method.
21 * Let S,C and T denote the sin, cos and tan respectively on
22 * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
23 * in [-pi/4 , +pi/4], and let n = k mod 4.
24 * We have
25 *
26 * n sin(x) cos(x) tan(x)
27 * ----------------------------------------------------------
28 * 0 S C T
29 * 1 C -S -1/T
30 * 2 -S -C T
31 * 3 -C S -1/T
32 * ----------------------------------------------------------
33 *
34 * Special cases:
35 * Let trig be any of sin, cos, or tan.
36 * trig(+-INF) is NaN, with signals;
37 * trig(NaN) is that NaN;
38 *
39 * Accuracy:
40 * TRIG(x) returns trig(x) nearly rounded
41 */
42
43#include "math_libm.h"
44#include "math_private.h"
45
46double tan(double x)
47{
48 double y[2],z=0.0;
49 int32_t n, ix;
50
51 /* High word of x. */
52 GET_HIGH_WORD(ix,x);
53
54 /* |x| ~< pi/4 */
55 ix &= 0x7fffffff;
56 if(ix <= 0x3fe921fb) return __kernel_tan(x,z,1);
57
58 /* tan(Inf or NaN) is NaN */
59 else if (ix>=0x7ff00000) return x-x; /* NaN */
60
61 /* argument reduction needed */
62 else {
63 n = __ieee754_rem_pio2(x,y);
64 return __kernel_tan(y[0],y[1],1-((n&1)<<1)); /* 1 -- n even
65 -1 -- n odd */
66 }
67}
68libm_hidden_def(tan)
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-29 23:35:11 +0000