*SDL project

5 files changed, 219 insertions, 6 deletions

diff --git a/Source/3rdParty/SDL2/src/libm/e_exp.c b/Source/3rdParty/SDL2/src/libm/e_exp.c
new file mode 100644
index 0000000..d8cd4a4
--- /dev/null
+++ b/Source/3rdParty/SDL2/src/libm/e_exp.c
@@ -0,0 +1,187 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/* __ieee754_exp(x)
+ * Returns the exponential of x.
+ *
+ * Method
+ *   1. Argument reduction:
+ *      Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
+ *	Given x, find r and integer k such that
+ *
+ *               x = k*ln2 + r,  |r| <= 0.5*ln2.
+ *
+ *      Here r will be represented as r = hi-lo for better
+ *	accuracy.
+ *
+ *   2. Approximation of exp(r) by a special rational function on
+ *	the interval [0,0.34658]:
+ *	Write
+ *	    R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
+ *      We use a special Reme algorithm on [0,0.34658] to generate
+ * 	a polynomial of degree 5 to approximate R. The maximum error
+ *	of this polynomial approximation is bounded by 2**-59. In
+ *	other words,
+ *	    R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
+ *  	(where z=r*r, and the values of P1 to P5 are listed below)
+ *	and
+ *	    |                  5          |     -59
+ *	    | 2.0+P1*z+...+P5*z   -  R(z) | <= 2
+ *	    |                             |
+ *	The computation of exp(r) thus becomes
+ *                             2*r
+ *		exp(r) = 1 + -------
+ *		              R - r
+ *                                 r*R1(r)
+ *		       = 1 + r + ----------- (for better accuracy)
+ *		                  2 - R1(r)
+ *	where
+ *			         2       4             10
+ *		R1(r) = r - (P1*r  + P2*r  + ... + P5*r   ).
+ *
+ *   3. Scale back to obtain exp(x):
+ *	From step 1, we have
+ *	   exp(x) = 2^k * exp(r)
+ *
+ * Special cases:
+ *	exp(INF) is INF, exp(NaN) is NaN;
+ *	exp(-INF) is 0, and
+ *	for finite argument, only exp(0)=1 is exact.
+ *
+ * Accuracy:
+ *	according to an error analysis, the error is always less than
+ *	1 ulp (unit in the last place).
+ *
+ * Misc. info.
+ *	For IEEE double
+ *	    if x >  7.09782712893383973096e+02 then exp(x) overflow
+ *	    if x < -7.45133219101941108420e+02 then exp(x) underflow
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+
+#include "math_libm.h"
+#include "math_private.h"
+
+static const double
+one	= 1.0,
+halF[2]	= {0.5,-0.5,},
+huge	= 1.0e+300,
+twom1000= 9.33263618503218878990e-302,     /* 2**-1000=0x01700000,0*/
+o_threshold=  7.09782712893383973096e+02,  /* 0x40862E42, 0xFEFA39EF */
+u_threshold= -7.45133219101941108420e+02,  /* 0xc0874910, 0xD52D3051 */
+ln2HI[2]   ={ 6.93147180369123816490e-01,  /* 0x3fe62e42, 0xfee00000 */
+	     -6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */
+ln2LO[2]   ={ 1.90821492927058770002e-10,  /* 0x3dea39ef, 0x35793c76 */
+	     -1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */
+invln2 =  1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
+P1   =  1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
+P2   = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
+P3   =  6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
+P4   = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
+P5   =  4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
+
+double __ieee754_exp(double x)	/* default IEEE double exp */
+{
+	double y;
+	double hi = 0.0;
+	double lo = 0.0;
+	double c;
+	double t;
+	int32_t k=0;
+	int32_t xsb;
+	u_int32_t hx;
+
+	GET_HIGH_WORD(hx,x);
+	xsb = (hx>>31)&1;		/* sign bit of x */
+	hx &= 0x7fffffff;		/* high word of |x| */
+
+    /* filter out non-finite argument */
+	if(hx >= 0x40862E42) {			/* if |x|>=709.78... */
+            if(hx>=0x7ff00000) {
+	        u_int32_t lx;
+		GET_LOW_WORD(lx,x);
+		if(((hx&0xfffff)|lx)!=0)
+		     return x+x; 		/* NaN */
+		else return (xsb==0)? x:0.0;	/* exp(+-inf)={inf,0} */
+	    }
+		#if 1
+		if(x > o_threshold) return huge*huge; /* overflow */
+		#else  /* !!! FIXME: check this: "huge * huge" is a compiler warning, maybe they wanted +Inf? */
+		if(x > o_threshold) return INFINITY; /* overflow */
+		#endif
+
+	    if(x < u_threshold) return twom1000*twom1000; /* underflow */
+	}
+
+    /* argument reduction */
+	if(hx > 0x3fd62e42) {		/* if  |x| > 0.5 ln2 */
+	    if(hx < 0x3FF0A2B2) {	/* and |x| < 1.5 ln2 */
+		hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
+	    } else {
+		k  = (int32_t) (invln2*x+halF[xsb]);
+		t  = k;
+		hi = x - t*ln2HI[0];	/* t*ln2HI is exact here */
+		lo = t*ln2LO[0];
+	    }
+	    x  = hi - lo;
+	}
+	else if(hx < 0x3e300000)  {	/* when |x|<2**-28 */
+	    if(huge+x>one) return one+x;/* trigger inexact */
+	}
+	else k = 0;
+
+    /* x is now in primary range */
+	t  = x*x;
+	c  = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
+	if(k==0) 	return one-((x*c)/(c-2.0)-x);
+	else 		y = one-((lo-(x*c)/(2.0-c))-hi);
+	if(k >= -1021) {
+	    u_int32_t hy;
+	    GET_HIGH_WORD(hy,y);
+	    SET_HIGH_WORD(y,hy+(k<<20));	/* add k to y's exponent */
+	    return y;
+	} else {
+	    u_int32_t hy;
+	    GET_HIGH_WORD(hy,y);
+	    SET_HIGH_WORD(y,hy+((k+1000)<<20));	/* add k to y's exponent */
+	    return y*twom1000;
+	}
+}
+
+/*
+ * wrapper exp(x)
+ */
+#ifndef _IEEE_LIBM
+double exp(double x)
+{
+	static const double o_threshold =  7.09782712893383973096e+02; /* 0x40862E42, 0xFEFA39EF */
+	static const double u_threshold = -7.45133219101941108420e+02; /* 0xc0874910, 0xD52D3051 */
+
+	double z = __ieee754_exp(x);
+	if (_LIB_VERSION == _IEEE_)
+		return z;
+	if (isfinite(x)) {
+		if (x > o_threshold)
+			return __kernel_standard(x, x, 6); /* exp overflow */
+		if (x < u_threshold)
+			return __kernel_standard(x, x, 7); /* exp underflow */
+	}
+	return z;
+}
+#else
+strong_alias(__ieee754_exp, exp)
+#endif
+libm_hidden_def(exp)
diff --git a/Source/3rdParty/SDL2/src/libm/e_rem_pio2.c b/Source/3rdParty/SDL2/src/libm/e_rem_pio2.c
index df7c2b8..5e055d6 100644
--- a/Source/3rdParty/SDL2/src/libm/e_rem_pio2.c
+++ b/Source/3rdParty/SDL2/src/libm/e_rem_pio2.c
@@ -154,7 +154,7 @@ int32_t attribute_hidden __ieee754_rem_pio2(double x, double *y)
 	}
 	tx[2] = z;
 	nx = 3;
-	while(tx[nx-1]==zero) nx--;	/* skip zero term */
+	while((nx > 0) && tx[nx-1]==zero) nx--;	/* skip zero term */
 	n  =  __kernel_rem_pio2(tx,y,e0,nx,2,two_over_pi);
 	if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;}
 	return n;
diff --git a/Source/3rdParty/SDL2/src/libm/k_rem_pio2.c b/Source/3rdParty/SDL2/src/libm/k_rem_pio2.c
index 7b04275..393db54 100644
--- a/Source/3rdParty/SDL2/src/libm/k_rem_pio2.c
+++ b/Source/3rdParty/SDL2/src/libm/k_rem_pio2.c
@@ -128,6 +128,8 @@
 #include "math_libm.h"
 #include "math_private.h"
 
+#include "SDL_assert.h"
+
 static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
 
 static const double PIo2[] = {
@@ -147,13 +149,19 @@ one    = 1.0,
 two24   =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
 twon24  =  5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
 
-int attribute_hidden __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int32_t *ipio2)
+int32_t attribute_hidden __kernel_rem_pio2(double *x, double *y, int e0, int nx, const unsigned int prec, const int32_t *ipio2)
 {
 	int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
 	double z,fw,f[20],fq[20],q[20];
 
+	if (nx < 1) {
+		return 0;
+	}
+
     /* initialize jk*/
+	SDL_assert(prec < SDL_arraysize(init_jk));
 	jk = init_jk[prec];
+	SDL_assert(jk > 0);
 	jp = jk;
 
     /* determine jx,jv,q0, note that 3>q0 */
@@ -164,6 +172,9 @@ int attribute_hidden __kernel_rem_pio2(double *x, double *y, int e0, int nx, int
     /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
 	j = jv-jx; m = jx+jk;
 	for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
+	if ((m+1) < SDL_arraysize(f)) {
+        SDL_memset(&f[m+1], 0, sizeof (f) - ((m+1) * sizeof (f[0])));
+    }
 
     /* compute q[0],q[1],...q[jk] */
 	for (i=0;i<=jk;i++) {
@@ -179,6 +190,9 @@ recompute:
 	    iq[i] =  (int32_t)(z-two24*fw);
 	    z     =  q[j-1]+fw;
 	}
+	if (jz < SDL_arraysize(iq)) {
+        SDL_memset(&iq[jz], 0, sizeof (q) - (jz * sizeof (iq[0])));
+    }
 
     /* compute n */
 	z  = scalbn(z,q0);		/* actual value of z */
@@ -238,7 +252,8 @@ recompute:
     /* chop off zero terms */
 	if(z==0.0) {
 	    jz -= 1; q0 -= 24;
-	    while(iq[jz]==0) { jz--; q0-=24;}
+		SDL_assert(jz >= 0);
+	    while(iq[jz]==0) { jz--; SDL_assert(jz >= 0); q0-=24;}
 	} else { /* break z into 24-bit if necessary */
 	    z = scalbn(z,-q0);
 	    if(z>=two24) {
@@ -260,6 +275,9 @@ recompute:
 	    for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
 	    fq[jz-i] = fw;
 	}
+	if ((jz+1) < SDL_arraysize(f)) {
+        SDL_memset(&fq[jz+1], 0, sizeof (fq) - ((jz+1) * sizeof (fq[0])));
+    }
 
     /* compress fq[] into y[] */
 	switch(prec) {
diff --git a/Source/3rdParty/SDL2/src/libm/math_libm.h b/Source/3rdParty/SDL2/src/libm/math_libm.h
index eb7bdd5..3c751c5 100644
--- a/Source/3rdParty/SDL2/src/libm/math_libm.h
+++ b/Source/3rdParty/SDL2/src/libm/math_libm.h
@@ -18,6 +18,10 @@
      misrepresented as being the original software.
   3. This notice may not be removed or altered from any source distribution.
 */
+
+#ifndef math_libm_h_
+#define math_libm_h_
+
 #include "../SDL_internal.h"
 
 /* Math routines from uClibc: http://www.uclibc.org */
@@ -26,6 +30,7 @@ double SDL_uclibc_atan(double x);
 double SDL_uclibc_atan2(double y, double x);    
 double SDL_uclibc_copysign(double x, double y);       
 double SDL_uclibc_cos(double x);         
+double SDL_uclibc_exp(double x);
 double SDL_uclibc_fabs(double x);        
 double SDL_uclibc_floor(double x);
 double SDL_uclibc_fmod(double x, double y);
@@ -37,4 +42,6 @@ double SDL_uclibc_sin(double x);
 double SDL_uclibc_sqrt(double x);
 double SDL_uclibc_tan(double x);
 
+#endif /* math_libm_h_ */
+
 /* vi: set ts=4 sw=4 expandtab: */
diff --git a/Source/3rdParty/SDL2/src/libm/math_private.h b/Source/3rdParty/SDL2/src/libm/math_private.h
index 1c0c8a4..d0ef66a 100644
--- a/Source/3rdParty/SDL2/src/libm/math_private.h
+++ b/Source/3rdParty/SDL2/src/libm/math_private.h
@@ -35,6 +35,7 @@ typedef unsigned int u_int32_t;
 #define __ieee754_atan2 SDL_uclibc_atan2
 #define copysign        SDL_uclibc_copysign
 #define cos             SDL_uclibc_cos
+#define __ieee754_exp   SDL_uclibc_exp
 #define fabs            SDL_uclibc_fabs
 #define floor           SDL_uclibc_floor
 #define __ieee754_fmod  SDL_uclibc_fmod
@@ -206,7 +207,7 @@ __ieee754_sqrt(double)
      extern double __ieee754_jn(int, double) attribute_hidden;
      extern double __ieee754_yn(int, double) attribute_hidden;
      extern double __ieee754_remainder(double, double) attribute_hidden;
-     extern int __ieee754_rem_pio2(double, double *) attribute_hidden;
+     extern int32_t __ieee754_rem_pio2(double, double *) attribute_hidden;
 #if defined(_SCALB_INT)
      extern double __ieee754_scalb(double, int) attribute_hidden;
 #else
@@ -220,7 +221,7 @@ __ieee754_sqrt(double)
      extern double __kernel_sin(double, double, int) attribute_hidden;
      extern double __kernel_cos(double, double) attribute_hidden;
      extern double __kernel_tan(double, double, int) attribute_hidden;
-     extern int __kernel_rem_pio2(double *, double *, int, int, int,
-                                  const int *) attribute_hidden;
+     extern int32_t __kernel_rem_pio2(double *, double *, int, int, const unsigned int,
+                                  const int32_t *) attribute_hidden;
 
 #endif /* _MATH_PRIVATE_H_ */