From 397dff81ee6efb4cca5593426a02db682e8b7e84 Mon Sep 17 00:00:00 2001 From: Akos Kiss Date: Sun, 13 Mar 2016 22:30:47 +0100 Subject: [PATCH] Specialize fdlibm to jerry Keep IEEE code paths only: * removed SVID, XOPEN, POSIX code paths from everywhere. * deleted s_lib_version.c, as version is only useful if multiple standards are supported. * deleted k_standard.c, as it handles non-IEEE exception cases only. * renamed the e_{acos,asin,atan2,exp,fmod,log,pow,sqrt}.c sources as s_.*, dropped the __ieee754_ prefix from the names of the appropriate functions therein, and deleted the w_{acos,asin,atan2,exp,fmod,log,pow,sqrt}.c wrapper code. Keep C99 declaration variants only: * removed old C-style function declaration variants. * removed data declaration variants where const qualifier was not used. Clean unused sources/functions: * removed s_{rint,significand,tanh}.c and the appropriate functions defined therein. JerryScript-DCO-1.0-Signed-off-by: Akos Kiss akiss@inf.u-szeged.hu --- third-party/fdlibm/e_rem_pio2.c | 59 +- third-party/fdlibm/include/fdlibm.h | 174 +---- third-party/fdlibm/k_cos.c | 29 +- third-party/fdlibm/k_rem_pio2.c | 73 +- third-party/fdlibm/k_sin.c | 27 +- third-party/fdlibm/k_standard.c | 733 -------------------- third-party/fdlibm/k_tan.c | 12 +- third-party/fdlibm/{e_acos.c => s_acos.c} | 25 +- third-party/fdlibm/{e_asin.c => s_asin.c} | 25 +- third-party/fdlibm/s_atan.c | 35 +- third-party/fdlibm/{e_atan2.c => s_atan2.c} | 23 +- third-party/fdlibm/s_ceil.c | 15 +- third-party/fdlibm/s_copysign.c | 7 +- third-party/fdlibm/s_cos.c | 13 +- third-party/fdlibm/{e_exp.c => s_exp.c} | 41 +- third-party/fdlibm/s_fabs.c | 7 +- third-party/fdlibm/s_finite.c | 9 +- third-party/fdlibm/s_floor.c | 15 +- third-party/fdlibm/{e_fmod.c => s_fmod.c} | 17 +- third-party/fdlibm/s_isnan.c | 9 +- third-party/fdlibm/s_lib_version.c | 35 - third-party/fdlibm/{e_log.c => s_log.c} | 47 +- third-party/fdlibm/{e_pow.c => s_pow.c} | 51 +- third-party/fdlibm/s_rint.c | 84 --- third-party/fdlibm/s_scalbn.c | 23 +- third-party/fdlibm/s_significand.c | 30 - third-party/fdlibm/s_sin.c | 13 +- third-party/fdlibm/{e_sqrt.c => s_sqrt.c} | 102 ++- third-party/fdlibm/s_tan.c | 13 +- third-party/fdlibm/s_tanh.c | 82 --- third-party/fdlibm/w_acos.c | 39 -- third-party/fdlibm/w_asin.c | 41 -- third-party/fdlibm/w_atan2.c | 40 -- third-party/fdlibm/w_exp.c | 48 -- third-party/fdlibm/w_fmod.c | 39 -- third-party/fdlibm/w_log.c | 39 -- third-party/fdlibm/w_pow.c | 59 -- third-party/fdlibm/w_sqrt.c | 38 - 38 files changed, 266 insertions(+), 1905 deletions(-) delete mode 100644 third-party/fdlibm/k_standard.c rename third-party/fdlibm/{e_acos.c => s_acos.c} (87%) rename third-party/fdlibm/{e_asin.c => s_asin.c} (89%) rename third-party/fdlibm/{e_atan2.c => s_atan2.c} (90%) rename third-party/fdlibm/{e_exp.c => s_exp.c} (85%) rename third-party/fdlibm/{e_fmod.c => s_fmod.c} (92%) delete mode 100644 third-party/fdlibm/s_lib_version.c rename third-party/fdlibm/{e_log.c => s_log.c} (82%) rename third-party/fdlibm/{e_pow.c => s_pow.c} (92%) delete mode 100644 third-party/fdlibm/s_rint.c delete mode 100644 third-party/fdlibm/s_significand.c rename third-party/fdlibm/{e_sqrt.c => s_sqrt.c} (88%) delete mode 100644 third-party/fdlibm/s_tanh.c delete mode 100644 third-party/fdlibm/w_acos.c delete mode 100644 third-party/fdlibm/w_asin.c delete mode 100644 third-party/fdlibm/w_atan2.c delete mode 100644 third-party/fdlibm/w_exp.c delete mode 100644 third-party/fdlibm/w_fmod.c delete mode 100644 third-party/fdlibm/w_log.c delete mode 100644 third-party/fdlibm/w_pow.c delete mode 100644 third-party/fdlibm/w_sqrt.c diff --git a/third-party/fdlibm/e_rem_pio2.c b/third-party/fdlibm/e_rem_pio2.c index 7242bb232..46e2a9e96 100644 --- a/third-party/fdlibm/e_rem_pio2.c +++ b/third-party/fdlibm/e_rem_pio2.c @@ -13,39 +13,31 @@ */ /* __ieee754_rem_pio2(x,y) - * - * return the remainder of x rem pi/2 in y[0]+y[1] + * + * return the remainder of x rem pi/2 in y[0]+y[1] * use __kernel_rem_pio2() */ #include "fdlibm.h" /* - * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi + * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi */ -#ifdef __STDC__ static const int two_over_pi[] = { -#else -static int two_over_pi[] = { -#endif -0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, -0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, -0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, -0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, -0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, -0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, -0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, -0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, -0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, -0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, -0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, +0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, +0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, +0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, +0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, +0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, +0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, +0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, +0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, +0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, +0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, +0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, }; -#ifdef __STDC__ static const int npio2_hw[] = { -#else -static int npio2_hw[] = { -#endif 0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, 0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C, 0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A, @@ -64,11 +56,7 @@ static int npio2_hw[] = { * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) */ -#ifdef __STDC__ -static const double -#else -static double -#endif +static const double zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ @@ -80,12 +68,7 @@ pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */ pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ -#ifdef __STDC__ - int __ieee754_rem_pio2(double x, double *y) -#else - int __ieee754_rem_pio2(x,y) - double x,y[]; -#endif +int __ieee754_rem_pio2(double x, double *y) { double z,w,t,r,fn; double tx[3]; @@ -126,7 +109,7 @@ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ fn = (double)n; r = t-fn*pio2_1; w = fn*pio2_1t; /* 1st round good to 85 bit */ - if(n<32&&ix!=npio2_hw[n-1]) { + if(n<32&&ix!=npio2_hw[n-1]) { y[0] = r-w; /* quick check no cancellation */ } else { j = ix>>20; @@ -134,16 +117,16 @@ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ i = j-(((__HI(y[0]))>>20)&0x7ff); if(i>16) { /* 2nd iteration needed, good to 118 */ t = r; - w = fn*pio2_2; + w = fn*pio2_2; r = t-w; - w = fn*pio2_2t-((t-r)-w); + w = fn*pio2_2t-((t-r)-w); y[0] = r-w; i = j-(((__HI(y[0]))>>20)&0x7ff); if(i>49) { /* 3rd iteration need, 151 bits acc */ t = r; /* will cover all possible cases */ w = fn*pio2_3; r = t-w; - w = fn*pio2_3t-((t-r)-w); + w = fn*pio2_3t-((t-r)-w); y[0] = r-w; } } @@ -152,7 +135,7 @@ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} else return n; } - /* + /* * all other (large) arguments */ if(ix>=0x7ff00000) { /* x is inf or NaN */ diff --git a/third-party/fdlibm/include/fdlibm.h b/third-party/fdlibm/include/fdlibm.h index fcee3cfeb..ecf5e7770 100644 --- a/third-party/fdlibm/include/fdlibm.h +++ b/third-party/fdlibm/include/fdlibm.h @@ -24,56 +24,17 @@ #ifdef __LITTLE_ENDIAN #define __HI(x) *(1+(int*)&x) #define __LO(x) *(int*)&x -#define __HIp(x) *(1+(int*)x) -#define __LOp(x) *(int*)x #else #define __HI(x) *(int*)&x #define __LO(x) *(1+(int*)&x) -#define __HIp(x) *(int*)x -#define __LOp(x) *(1+(int*)x) -#endif - -#ifdef __STDC__ -#define __P(p) p -#else -#define __P(p) () #endif /* * ANSI/POSIX */ -extern int signgam; - #define MAXFLOAT ((float)3.40282346638528860e+38) -enum fdversion {fdlibm_ieee = -1, fdlibm_svid, fdlibm_xopen, fdlibm_posix}; - -#define _LIB_VERSION_TYPE enum fdversion -#define _LIB_VERSION _fdlib_version - -/* if global variable _LIB_VERSION is not desirable, one may - * change the following to be a constant by: - * #define _LIB_VERSION_TYPE const enum version - * In that case, after one initializes the value _LIB_VERSION (see - * s_lib_version.c) during compile time, it cannot be modified - * in the middle of a program - */ -extern _LIB_VERSION_TYPE _LIB_VERSION; - -#define _IEEE_ fdlibm_ieee -#define _SVID_ fdlibm_svid -#define _XOPEN_ fdlibm_xopen -#define _POSIX_ fdlibm_posix - -struct exception { - int type; - char *name; - double arg1; - double arg2; - double retval; -}; - #define HUGE MAXFLOAT /* @@ -93,126 +54,39 @@ struct exception { /* * ANSI/POSIX */ -extern double acos __P((double)); -extern double asin __P((double)); -extern double atan __P((double)); -extern double atan2 __P((double, double)); -extern double cos __P((double)); -extern double sin __P((double)); -extern double tan __P((double)); +extern double acos (double); +extern double asin (double); +extern double atan (double); +extern double atan2 (double, double); +extern double cos (double); +extern double sin (double); +extern double tan (double); -extern double cosh __P((double)); -extern double sinh __P((double)); -extern double tanh __P((double)); +extern double exp (double); +extern double log (double); -extern double exp __P((double)); -extern double frexp __P((double, int *)); -extern double ldexp __P((double, int)); -extern double log __P((double)); -extern double log10 __P((double)); -extern double modf __P((double, double *)); +extern double pow (double, double); +extern double sqrt (double); -extern double pow __P((double, double)); -extern double sqrt __P((double)); +extern double ceil (double); +extern double fabs (double); +extern double floor (double); +extern double fmod (double, double); -extern double ceil __P((double)); -extern double fabs __P((double)); -extern double floor __P((double)); -extern double fmod __P((double, double)); - -extern double erf __P((double)); -extern double erfc __P((double)); -extern double gamma __P((double)); -extern double hypot __P((double, double)); -extern int isnan __P((double)); -extern int finite __P((double)); -extern double j0 __P((double)); -extern double j1 __P((double)); -extern double jn __P((int, double)); -extern double lgamma __P((double)); -extern double y0 __P((double)); -extern double y1 __P((double)); -extern double yn __P((int, double)); - -extern double acosh __P((double)); -extern double asinh __P((double)); -extern double atanh __P((double)); -extern double cbrt __P((double)); -extern double logb __P((double)); -extern double nextafter __P((double, double)); -extern double remainder __P((double, double)); -#ifdef _SCALB_INT -extern double scalb __P((double, int)); -#else -extern double scalb __P((double, double)); -#endif - -extern int matherr __P((struct exception *)); - -/* - * IEEE Test Vector - */ -extern double significand __P((double)); +extern int isnan (double); +extern int finite (double); /* * Functions callable from C, intended to support IEEE arithmetic. */ -extern double copysign __P((double, double)); -extern int ilogb __P((double)); -extern double rint __P((double)); -extern double scalbn __P((double, int)); - -/* - * BSD math library entry points - */ -extern double expm1 __P((double)); -extern double log1p __P((double)); - -/* - * Reentrant version of gamma & lgamma; passes signgam back by reference - * as the second argument; user must allocate space for signgam. - */ -#ifdef _REENTRANT -extern double gamma_r __P((double, int *)); -extern double lgamma_r __P((double, int *)); -#endif /* _REENTRANT */ +extern double copysign (double, double); +extern double scalbn (double, int); /* ieee style elementary functions */ -extern double __ieee754_sqrt __P((double)); -extern double __ieee754_acos __P((double)); -extern double __ieee754_acosh __P((double)); -extern double __ieee754_log __P((double)); -extern double __ieee754_atanh __P((double)); -extern double __ieee754_asin __P((double)); -extern double __ieee754_atan2 __P((double,double)); -extern double __ieee754_exp __P((double)); -extern double __ieee754_cosh __P((double)); -extern double __ieee754_fmod __P((double,double)); -extern double __ieee754_pow __P((double,double)); -extern double __ieee754_lgamma_r __P((double,int *)); -extern double __ieee754_gamma_r __P((double,int *)); -extern double __ieee754_lgamma __P((double)); -extern double __ieee754_gamma __P((double)); -extern double __ieee754_log10 __P((double)); -extern double __ieee754_sinh __P((double)); -extern double __ieee754_hypot __P((double,double)); -extern double __ieee754_j0 __P((double)); -extern double __ieee754_j1 __P((double)); -extern double __ieee754_y0 __P((double)); -extern double __ieee754_y1 __P((double)); -extern double __ieee754_jn __P((int,double)); -extern double __ieee754_yn __P((int,double)); -extern double __ieee754_remainder __P((double,double)); -extern int __ieee754_rem_pio2 __P((double,double*)); -#ifdef _SCALB_INT -extern double __ieee754_scalb __P((double,int)); -#else -extern double __ieee754_scalb __P((double,double)); -#endif +extern int __ieee754_rem_pio2 (double,double*); /* fdlibm kernel function */ -extern double __kernel_standard __P((double,double,int)); -extern double __kernel_sin __P((double,double,int)); -extern double __kernel_cos __P((double,double)); -extern double __kernel_tan __P((double,double,int)); -extern int __kernel_rem_pio2 __P((double*,double*,int,int,int,const int*)); +extern double __kernel_sin (double,double,int); +extern double __kernel_cos (double,double); +extern double __kernel_tan (double,double,int); +extern int __kernel_rem_pio2 (double*,double*,int,int,int,const int*); diff --git a/third-party/fdlibm/k_cos.c b/third-party/fdlibm/k_cos.c index 7fb855d25..ebf1e4a59 100644 --- a/third-party/fdlibm/k_cos.c +++ b/third-party/fdlibm/k_cos.c @@ -15,7 +15,7 @@ * __kernel_cos( x, y ) * kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164 * Input x is assumed to be bounded by ~pi/4 in magnitude. - * Input y is the tail of x. + * Input y is the tail of x. * * Algorithm * 1. Since cos(-x) = cos(x), we need only to consider positive x. @@ -26,14 +26,14 @@ * cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x * where the remez error is * - * | 2 4 6 8 10 12 14 | -58 - * |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2 - * | | - * - * 4 6 8 10 12 14 + * | 2 4 6 8 10 12 14 | -58 + * |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2 + * | | + * + * 4 6 8 10 12 14 * 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then * cos(x) = 1 - x*x/2 + r - * since cos(x+y) ~ cos(x) - sin(x)*y + * since cos(x+y) ~ cos(x) - sin(x)*y * ~ cos(x) - x*y, * a correction term is necessary in cos(x) and hence * cos(x+y) = 1 - (x*x/2 - (r - x*y)) @@ -48,11 +48,7 @@ #include "fdlibm.h" -#ifdef __STDC__ -static const double -#else -static double -#endif +static const double one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */ C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */ @@ -61,12 +57,7 @@ C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */ C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */ C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */ -#ifdef __STDC__ - double __kernel_cos(double x, double y) -#else - double __kernel_cos(x, y) - double x,y; -#endif +double __kernel_cos(double x, double y) { double a,hz,z,r,qx; int ix; @@ -76,7 +67,7 @@ C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */ } z = x*x; r = z*(C1+z*(C2+z*(C3+z*(C4+z*(C5+z*C6))))); - if(ix < 0x3FD33333) /* if |x| < 0.3 */ + if(ix < 0x3FD33333) /* if |x| < 0.3 */ return one - (0.5*z - (z*r - x*y)); else { if(ix > 0x3fe90000) { /* x > 0.78125 */ diff --git a/third-party/fdlibm/k_rem_pio2.c b/third-party/fdlibm/k_rem_pio2.c index ec473ac0d..664fd9335 100644 --- a/third-party/fdlibm/k_rem_pio2.c +++ b/third-party/fdlibm/k_rem_pio2.c @@ -14,12 +14,12 @@ /* * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2) * double x[],y[]; int e0,nx,prec; int ipio2[]; - * - * __kernel_rem_pio2 return the last three digits of N with + * + * __kernel_rem_pio2 return the last three digits of N with * y = x - N*pi/2 * so that |y| < pi/2. * - * The method is to compute the integer (mod 8) and fraction parts of + * The method is to compute the integer (mod 8) and fraction parts of * (2/pi)*x without doing the full multiplication. In general we * skip the part of the product that are known to be a huge integer ( * more accurately, = 0 mod 8 ). Thus the number of operations are @@ -28,10 +28,10 @@ * (2/pi) is represented by an array of 24-bit integers in ipio2[]. * * Input parameters: - * x[] The input value (must be positive) is broken into nx + * x[] The input value (must be positive) is broken into nx * pieces of 24-bit integers in double precision format. - * x[i] will be the i-th 24 bit of x. The scaled exponent - * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 + * x[i] will be the i-th 24 bit of x. The scaled exponent + * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 * match x's up to 24 bits. * * Example of breaking a double positive z into x[0]+x[1]+x[2]: @@ -61,15 +61,15 @@ * * nx dimension of x[] * - * prec an integer indicating the precision: + * prec an integer indicating the precision: * 0 24 bits (single) * 1 53 bits (double) * 2 64 bits (extended) * 3 113 bits (quad) * * ipio2[] - * integer array, contains the (24*i)-th to (24*i+23)-th - * bit of 2/pi after binary point. The corresponding + * integer array, contains the (24*i)-th to (24*i+23)-th + * bit of 2/pi after binary point. The corresponding * floating value is * * ipio2[i] * 2^(-24(i+1)). @@ -80,12 +80,12 @@ * * Here is the description of some local variables: * - * jk jk+1 is the initial number of terms of ipio2[] needed + * jk jk+1 is the initial number of terms of ipio2[] needed * in the computation. The recommended value is 2,3,4, * 6 for single, double, extended,and quad. * - * jz local integer variable indicating the number of - * terms of ipio2[] used. + * jz local integer variable indicating the number of + * terms of ipio2[] used. * * jx nx - 1 * @@ -98,16 +98,16 @@ * * jp jp+1 is the number of terms in PIo2[] needed, jp = jk. * - * q[] double array with integral value, representing the + * q[] double array with integral value, representing the * 24-bits chunk of the product of x and 2/pi. * * q0 the corresponding exponent of q[0]. Note that the * exponent for q[i] would be q0-24*i. * * PIo2[] double precision array, obtained by cutting pi/2 - * into 24 bits chunks. + * into 24 bits chunks. * - * f[] ipio2[] in floating point + * f[] ipio2[] in floating point * * iq[] integer array by breaking up q[] in 24-bits chunk. * @@ -121,25 +121,17 @@ /* * 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 + * 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 "fdlibm.h" -#ifdef __STDC__ static const int init_jk[] = {2,3,4,6}; /* initial value for jk */ -#else -static int init_jk[] = {2,3,4,6}; -#endif -#ifdef __STDC__ static const double PIo2[] = { -#else -static double PIo2[] = { -#endif 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */ 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */ 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */ @@ -150,22 +142,13 @@ static double PIo2[] = { 2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */ }; -#ifdef __STDC__ -static const double -#else -static double -#endif +static const double zero = 0.0, one = 1.0, two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */ -#ifdef __STDC__ - int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int *ipio2) -#else - int __kernel_rem_pio2(x,y,e0,nx,prec,ipio2) - double x[], y[]; int e0,nx,prec; int ipio2[]; -#endif +int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int *ipio2) { int jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih; double z,fw,f[20],fq[20],q[20]; @@ -258,7 +241,7 @@ recompute: while(iq[jz]==0) { jz--; q0-=24;} } else { /* break z into 24-bit if necessary */ z = scalbn(z,-q0); - if(z>=two24) { + if(z>=two24) { fw = (double)((int)(twon24*z)); iq[jz] = (int)(z-two24*fw); jz += 1; q0 += 24; @@ -283,29 +266,29 @@ recompute: case 0: fw = 0.0; for (i=jz;i>=0;i--) fw += fq[i]; - y[0] = (ih==0)? fw: -fw; + y[0] = (ih==0)? fw: -fw; break; case 1: case 2: fw = 0.0; - for (i=jz;i>=0;i--) fw += fq[i]; - y[0] = (ih==0)? fw: -fw; + for (i=jz;i>=0;i--) fw += fq[i]; + y[0] = (ih==0)? fw: -fw; fw = fq[0]-fw; for (i=1;i<=jz;i++) fw += fq[i]; - y[1] = (ih==0)? fw: -fw; + y[1] = (ih==0)? fw: -fw; break; case 3: /* painful */ for (i=jz;i>0;i--) { - fw = fq[i-1]+fq[i]; + fw = fq[i-1]+fq[i]; fq[i] += fq[i-1]-fw; fq[i-1] = fw; } for (i=jz;i>1;i--) { - fw = fq[i-1]+fq[i]; + fw = fq[i-1]+fq[i]; fq[i] += fq[i-1]-fw; fq[i-1] = fw; } - for (fw=0.0,i=jz;i>=2;i--) fw += fq[i]; + for (fw=0.0,i=jz;i>=2;i--) fw += fq[i]; if(ih==0) { y[0] = fq[0]; y[1] = fq[1]; y[2] = fw; } else { diff --git a/third-party/fdlibm/k_sin.c b/third-party/fdlibm/k_sin.c index dfcad764e..959b29942 100644 --- a/third-party/fdlibm/k_sin.c +++ b/third-party/fdlibm/k_sin.c @@ -15,10 +15,10 @@ * kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854 * Input x is assumed to be bounded by ~pi/4 in magnitude. * Input y is the tail of x. - * Input iy indicates whether y is 0. (if iy=0, y assume to be 0). + * Input iy indicates whether y is 0. (if iy=0, y assume to be 0). * * Algorithm - * 1. Since sin(-x) = -sin(x), we need only to consider positive x. + * 1. Since sin(-x) = -sin(x), we need only to consider positive x. * 2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0. * 3. sin(x) is approximated by a polynomial of degree 13 on * [0,pi/4] @@ -26,13 +26,13 @@ * sin(x) ~ x + S1*x + ... + S6*x * where * - * |sin(x) 2 4 6 8 10 12 | -58 - * |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2 - * | x | - * + * |sin(x) 2 4 6 8 10 12 | -58 + * |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2 + * | x | + * * 4. sin(x+y) = sin(x) + sin'(x')*y * ~ sin(x) + (1-x*x/2)*y - * For better accuracy, let + * For better accuracy, let * 3 2 2 2 2 * r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6)))) * then 3 2 @@ -41,11 +41,7 @@ #include "fdlibm.h" -#ifdef __STDC__ -static const double -#else -static double -#endif +static const double half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */ S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */ @@ -54,12 +50,7 @@ S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */ S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */ S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */ -#ifdef __STDC__ - double __kernel_sin(double x, double y, int iy) -#else - double __kernel_sin(x, y, iy) - double x,y; int iy; /* iy=0 if y is zero */ -#endif +double __kernel_sin(double x, double y, int iy) { double z,r,v; int ix; diff --git a/third-party/fdlibm/k_standard.c b/third-party/fdlibm/k_standard.c deleted file mode 100644 index d6e4a82ec..000000000 --- a/third-party/fdlibm/k_standard.c +++ /dev/null @@ -1,733 +0,0 @@ - -/* @(#)k_standard.c 1.3 95/01/18 */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - * - */ - -#include "fdlibm.h" -#include - -#ifndef _USE_WRITE -#include /* fputs(), stderr */ -// #define WRITE2(u,v) fputs(u, stderr) -#else /* !defined(_USE_WRITE) */ -#include /* write */ -// #define WRITE2(u,v) write(2, u, v) -#undef fflush -#endif /* !defined(_USE_WRITE) */ - -static double zero = 0.0; /* used as const */ - -/* - * Standard conformance (non-IEEE) on exception cases. - * Mapping: - * 1 -- acos(|x|>1) - * 2 -- asin(|x|>1) - * 3 -- atan2(+-0,+-0) - * 4 -- hypot overflow - * 5 -- cosh overflow - * 6 -- exp overflow - * 7 -- exp underflow - * 8 -- y0(0) - * 9 -- y0(-ve) - * 10-- y1(0) - * 11-- y1(-ve) - * 12-- yn(0) - * 13-- yn(-ve) - * 14-- lgamma(finite) overflow - * 15-- lgamma(-integer) - * 16-- log(0) - * 17-- log(x<0) - * 18-- log10(0) - * 19-- log10(x<0) - * 20-- pow(0.0,0.0) - * 21-- pow(x,y) overflow - * 22-- pow(x,y) underflow - * 23-- pow(0,negative) - * 24-- pow(neg,non-integral) - * 25-- sinh(finite) overflow - * 26-- sqrt(negative) - * 27-- fmod(x,0) - * 28-- remainder(x,0) - * 29-- acosh(x<1) - * 30-- atanh(|x|>1) - * 31-- atanh(|x|=1) - * 32-- scalb overflow - * 33-- scalb underflow - * 34-- j0(|x|>X_TLOSS) - * 35-- y0(x>X_TLOSS) - * 36-- j1(|x|>X_TLOSS) - * 37-- y1(x>X_TLOSS) - * 38-- jn(|x|>X_TLOSS, n) - * 39-- yn(x>X_TLOSS, n) - * 40-- gamma(finite) overflow - * 41-- gamma(-integer) - * 42-- pow(NaN,0.0) - */ - - -#ifdef __STDC__ - double __kernel_standard(double x, double y, int type) -#else - double __kernel_standard(x,y,type) - double x,y; int type; -#endif -{ - struct exception exc; -#ifndef HUGE_VAL /* this is the only routine that uses HUGE_VAL */ -#define HUGE_VAL inf - double inf = 0.0; - - __HI(inf) = 0x7ff00000; /* set inf to infinite */ -#endif - -#ifdef _USE_WRITE - (void) fflush(stdout); -#endif - exc.arg1 = x; - exc.arg2 = y; - switch(type) { - case 1: - /* acos(|x|>1) */ - exc.type = DOMAIN; - exc.name = "acos"; - exc.retval = zero; - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if(_LIB_VERSION == _SVID_) { - // // (void) WRITE2("acos: DOMAIN error\n", 19); - // } - // errno = EDOM; - // } - break; - case 2: - /* asin(|x|>1) */ - exc.type = DOMAIN; - exc.name = "asin"; - exc.retval = zero; - // if(_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if(_LIB_VERSION == _SVID_) { - // // (void) WRITE2("asin: DOMAIN error\n", 19); - // } - // errno = EDOM; - // } - break; - case 3: - /* atan2(+-0,+-0) */ - exc.arg1 = y; - exc.arg2 = x; - exc.type = DOMAIN; - exc.name = "atan2"; - exc.retval = zero; - // if(_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if(_LIB_VERSION == _SVID_) { - // // (void) WRITE2("atan2: DOMAIN error\n", 20); - // } - // errno = EDOM; - // } - break; - case 4: - /* hypot(finite,finite) overflow */ - exc.type = OVERFLOW; - exc.name = "hypot"; - if (_LIB_VERSION == _SVID_) - exc.retval = HUGE; - else - exc.retval = HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // errno = ERANGE; - // } - break; - case 5: - /* cosh(finite) overflow */ - exc.type = OVERFLOW; - exc.name = "cosh"; - if (_LIB_VERSION == _SVID_) - exc.retval = HUGE; - else - exc.retval = HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // errno = ERANGE; - // } - break; - case 6: - /* exp(finite) overflow */ - exc.type = OVERFLOW; - exc.name = "exp"; - if (_LIB_VERSION == _SVID_) - exc.retval = HUGE; - else - exc.retval = HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // errno = ERANGE; - // } - break; - case 7: - /* exp(finite) underflow */ - exc.type = UNDERFLOW; - exc.name = "exp"; - exc.retval = zero; - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // errno = ERANGE; - // } - break; - case 8: - /* y0(0) = -inf */ - exc.type = DOMAIN; /* should be SING for IEEE */ - exc.name = "y0"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("y0: DOMAIN error\n", 17); - // } - // errno = EDOM; - // } - break; - case 9: - /* y0(x<0) = NaN */ - exc.type = DOMAIN; - exc.name = "y0"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("y0: DOMAIN error\n", 17); - // } - // errno = EDOM; - // } - break; - case 10: - /* y1(0) = -inf */ - exc.type = DOMAIN; /* should be SING for IEEE */ - exc.name = "y1"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("y1: DOMAIN error\n", 17); - // } - // errno = EDOM; - // } - break; - case 11: - /* y1(x<0) = NaN */ - exc.type = DOMAIN; - exc.name = "y1"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("y1: DOMAIN error\n", 17); - // } - // errno = EDOM; - // } - break; - case 12: - /* yn(n,0) = -inf */ - exc.type = DOMAIN; /* should be SING for IEEE */ - exc.name = "yn"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("yn: DOMAIN error\n", 17); - // } - // errno = EDOM; - // } - break; - case 13: - /* yn(x<0) = NaN */ - exc.type = DOMAIN; - exc.name = "yn"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("yn: DOMAIN error\n", 17); - // } - // errno = EDOM; - // } - break; - case 14: - /* lgamma(finite) overflow */ - exc.type = OVERFLOW; - exc.name = "lgamma"; - if (_LIB_VERSION == _SVID_) - exc.retval = HUGE; - else - exc.retval = HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // errno = ERANGE; - // } - break; - case 15: - /* lgamma(-integer) or lgamma(0) */ - exc.type = SING; - exc.name = "lgamma"; - if (_LIB_VERSION == _SVID_) - exc.retval = HUGE; - else - exc.retval = HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("lgamma: SING error\n", 19); - // } - // errno = EDOM; - // } - break; - case 16: - /* log(0) */ - exc.type = SING; - exc.name = "log"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("log: SING error\n", 16); - // } - // errno = EDOM; - // } - break; - case 17: - /* log(x<0) */ - exc.type = DOMAIN; - exc.name = "log"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("log: DOMAIN error\n", 18); - // } - // errno = EDOM; - // } - break; - case 18: - /* log10(0) */ - exc.type = SING; - exc.name = "log10"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("log10: SING error\n", 18); - // } - // errno = EDOM; - // } - break; - case 19: - /* log10(x<0) */ - exc.type = DOMAIN; - exc.name = "log10"; - if (_LIB_VERSION == _SVID_) - exc.retval = -HUGE; - else - exc.retval = -HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("log10: DOMAIN error\n", 20); - // } - // errno = EDOM; - // } - break; - case 20: - /* pow(0.0,0.0) */ - /* error only if _LIB_VERSION == _SVID_ */ - exc.type = DOMAIN; - exc.name = "pow"; - exc.retval = zero; - // if (_LIB_VERSION != _SVID_) exc.retval = 1.0; - // else if (!matherr(&exc)) { - // // (void) WRITE2("pow(0,0): DOMAIN error\n", 23); - // errno = EDOM; - // } - break; - case 21: - /* pow(x,y) overflow */ - exc.type = OVERFLOW; - exc.name = "pow"; - if (_LIB_VERSION == _SVID_) { - exc.retval = HUGE; - y *= 0.5; - if(xzero) ? HUGE : -HUGE); - else - exc.retval = ( (x>zero) ? HUGE_VAL : -HUGE_VAL); - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // errno = ERANGE; - // } - break; - case 26: - /* sqrt(x<0) */ - exc.type = DOMAIN; - exc.name = "sqrt"; - if (_LIB_VERSION == _SVID_) - exc.retval = zero; - else - exc.retval = zero/zero; - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("sqrt: DOMAIN error\n", 19); - // } - // errno = EDOM; - // } - break; - case 27: - /* fmod(x,0) */ - exc.type = DOMAIN; - exc.name = "fmod"; - if (_LIB_VERSION == _SVID_) - exc.retval = x; - else - exc.retval = zero/zero; - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("fmod: DOMAIN error\n", 20); - // } - // errno = EDOM; - // } - break; - case 28: - /* remainder(x,0) */ - exc.type = DOMAIN; - exc.name = "remainder"; - exc.retval = zero/zero; - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("remainder: DOMAIN error\n", 24); - // } - // errno = EDOM; - // } - break; - case 29: - /* acosh(x<1) */ - exc.type = DOMAIN; - exc.name = "acosh"; - exc.retval = zero/zero; - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("acosh: DOMAIN error\n", 20); - // } - // errno = EDOM; - // } - break; - case 30: - /* atanh(|x|>1) */ - exc.type = DOMAIN; - exc.name = "atanh"; - exc.retval = zero/zero; - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("atanh: DOMAIN error\n", 20); - // } - // errno = EDOM; - // } - break; - case 31: - /* atanh(|x|=1) */ - exc.type = SING; - exc.name = "atanh"; - exc.retval = x/zero; /* sign(x)*inf */ - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("atanh: SING error\n", 18); - // } - // errno = EDOM; - // } - break; - case 32: - /* scalb overflow; SVID also returns +-HUGE_VAL */ - exc.type = OVERFLOW; - exc.name = "scalb"; - exc.retval = x > zero ? HUGE_VAL : -HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // errno = ERANGE; - // } - break; - case 33: - /* scalb underflow */ - exc.type = UNDERFLOW; - exc.name = "scalb"; - exc.retval = copysign(zero,x); - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // errno = ERANGE; - // } - break; - case 34: - /* j0(|x|>X_TLOSS) */ - exc.type = TLOSS; - exc.name = "j0"; - exc.retval = zero; - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2(exc.name, 2); - // // (void) WRITE2(": TLOSS error\n", 14); - // } - // errno = ERANGE; - // } - break; - case 35: - /* y0(x>X_TLOSS) */ - exc.type = TLOSS; - exc.name = "y0"; - exc.retval = zero; - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2(exc.name, 2); - // // (void) WRITE2(": TLOSS error\n", 14); - // } - // errno = ERANGE; - // } - break; - case 36: - /* j1(|x|>X_TLOSS) */ - exc.type = TLOSS; - exc.name = "j1"; - exc.retval = zero; - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2(exc.name, 2); - // // (void) WRITE2(": TLOSS error\n", 14); - // } - // errno = ERANGE; - // } - break; - case 37: - /* y1(x>X_TLOSS) */ - exc.type = TLOSS; - exc.name = "y1"; - exc.retval = zero; - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2(exc.name, 2); - // // (void) WRITE2(": TLOSS error\n", 14); - // } - // errno = ERANGE; - // } - break; - case 38: - /* jn(|x|>X_TLOSS) */ - exc.type = TLOSS; - exc.name = "jn"; - exc.retval = zero; - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2(exc.name, 2); - // // (void) WRITE2(": TLOSS error\n", 14); - // } - // errno = ERANGE; - // } - break; - case 39: - /* yn(x>X_TLOSS) */ - exc.type = TLOSS; - exc.name = "yn"; - exc.retval = zero; - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2(exc.name, 2); - // // (void) WRITE2(": TLOSS error\n", 14); - // } - // errno = ERANGE; - // } - break; - case 40: - /* gamma(finite) overflow */ - exc.type = OVERFLOW; - exc.name = "gamma"; - if (_LIB_VERSION == _SVID_) - exc.retval = HUGE; - else - exc.retval = HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = ERANGE; - // else if (!matherr(&exc)) { - // errno = ERANGE; - // } - break; - case 41: - /* gamma(-integer) or gamma(0) */ - exc.type = SING; - exc.name = "gamma"; - if (_LIB_VERSION == _SVID_) - exc.retval = HUGE; - else - exc.retval = HUGE_VAL; - // if (_LIB_VERSION == _POSIX_) - // errno = EDOM; - // else if (!matherr(&exc)) { - // if (_LIB_VERSION == _SVID_) { - // // (void) WRITE2("gamma: SING error\n", 18); - // } - // errno = EDOM; - // } - break; - case 42: - /* pow(NaN,0.0) */ - /* error only if _LIB_VERSION == _SVID_ & _XOPEN_ */ - exc.type = DOMAIN; - exc.name = "pow"; - exc.retval = x; - // if (_LIB_VERSION == _IEEE_ || - // _LIB_VERSION == _POSIX_) exc.retval = 1.0; - // else if (!matherr(&exc)) { - // errno = EDOM; - // } - // break; - } - return exc.retval; -} diff --git a/third-party/fdlibm/k_tan.c b/third-party/fdlibm/k_tan.c index 4c2c4f1ed..c49762fbd 100644 --- a/third-party/fdlibm/k_tan.c +++ b/third-party/fdlibm/k_tan.c @@ -10,7 +10,6 @@ * ==================================================== */ -/* INDENT OFF */ /* __kernel_tan( x, y, k ) * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854 * Input x is assumed to be bounded by ~pi/4 in magnitude. @@ -26,9 +25,9 @@ * tan(x) ~ x + T1*x + ... + T13*x * where * - * |tan(x) 2 4 26 | -59.2 - * |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2 - * | x | + * |tan(x) 2 4 26 | -59.2 + * |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2 + * | x | * * Note: tan(x+y) = tan(x) + tan'(x)*y * ~ tan(x) + (1+x*x)*y @@ -68,10 +67,9 @@ static const double xxx[] = { #define pio4 xxx[14] #define pio4lo xxx[15] #define T xxx -/* INDENT ON */ -double -__kernel_tan(double x, double y, int iy) { +double __kernel_tan(double x, double y, int iy) +{ double z, r, v, w, s; int ix, hx; diff --git a/third-party/fdlibm/e_acos.c b/third-party/fdlibm/s_acos.c similarity index 87% rename from third-party/fdlibm/e_acos.c rename to third-party/fdlibm/s_acos.c index d7c9ed225..4bf4166c6 100644 --- a/third-party/fdlibm/e_acos.c +++ b/third-party/fdlibm/s_acos.c @@ -11,15 +11,15 @@ * ==================================================== */ -/* __ieee754_acos(x) - * Method : +/* acos(x) + * Method : * acos(x) = pi/2 - asin(x) * acos(-x) = pi/2 + asin(x) * For |x|<=0.5 * acos(x) = pi/2 - (x + x*x^2*R(x^2)) (see asin.c) * For x>0.5 - * acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2))) - * = 2asin(sqrt((1-x)/2)) + * acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2))) + * = 2asin(sqrt((1-x)/2)) * = 2s + 2s*z*R(z) ...z=(1-x)/2, s=sqrt(z) * = 2f + (2c + 2s*z*R(z)) * where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term @@ -37,13 +37,9 @@ #include "fdlibm.h" -#ifdef __STDC__ -static const double -#else -static double -#endif -one= 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ -pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ +static const double +one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ +pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */ pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */ @@ -57,12 +53,7 @@ qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */ qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ -#ifdef __STDC__ - double __ieee754_acos(double x) -#else - double __ieee754_acos(x) - double x; -#endif +double acos(double x) { double z,p,q,r,w,s,c,df; int hx,ix; diff --git a/third-party/fdlibm/e_asin.c b/third-party/fdlibm/s_asin.c similarity index 89% rename from third-party/fdlibm/e_asin.c rename to third-party/fdlibm/s_asin.c index 1d503fc22..9a7581460 100644 --- a/third-party/fdlibm/e_asin.c +++ b/third-party/fdlibm/s_asin.c @@ -11,13 +11,13 @@ * ==================================================== */ -/* __ieee754_asin(x) - * Method : +/* asin(x) + * Method : * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ... * we approximate asin(x) on [0,0.5] by * asin(x) = x + x*x^2*R(x^2) * where - * R(x^2) is a rational approximation of (asin(x)-x)/x^3 + * R(x^2) is a rational approximation of (asin(x)-x)/x^3 * and its remez error is bounded by * |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75) * @@ -44,11 +44,7 @@ #include "fdlibm.h" -#ifdef __STDC__ static const double -#else -static double -#endif one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ huge = 1.000e+300, pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */ @@ -66,12 +62,7 @@ qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */ qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */ qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ -#ifdef __STDC__ - double __ieee754_asin(double x) -#else - double __ieee754_asin(x) - double x; -#endif +double asin(double x) { double t = 0,w,p,q,c,r,s; int hx,ix; @@ -80,8 +71,8 @@ qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ if(ix>= 0x3ff00000) { /* |x|>= 1 */ if(((ix-0x3ff00000)|__LO(x))==0) /* asin(1)=+-pi/2 with inexact */ - return x*pio2_hi+x*pio2_lo; - return (x-x)/(x-x); /* asin(|x|>1) is NaN */ + return x*pio2_hi+x*pio2_lo; + return (x-x)/(x-x); /* asin(|x|>1) is NaN */ } else if (ix<0x3fe00000) { /* |x|<0.5 */ if(ix<0x3e400000) { /* if |x| < 2**-27 */ if(huge+x>one) return x;/* return x with inexact if x!=0*/ @@ -109,6 +100,6 @@ qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ p = 2.0*s*r-(pio2_lo-2.0*c); q = pio4_hi-2.0*w; t = pio4_hi-(p-q); - } - if(hx>0) return t; else return -t; + } + if(hx>0) return t; else return -t; } diff --git a/third-party/fdlibm/s_atan.c b/third-party/fdlibm/s_atan.c index 0093eafdf..fecb7fd1e 100644 --- a/third-party/fdlibm/s_atan.c +++ b/third-party/fdlibm/s_atan.c @@ -26,41 +26,29 @@ * [39/16,INF] atan(x) = atan(INF) + atan( -1/t ) * * 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 + * 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 "fdlibm.h" -#ifdef __STDC__ static const double atanhi[] = { -#else -static double atanhi[] = { -#endif 4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */ 7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */ 9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */ 1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */ }; -#ifdef __STDC__ static const double atanlo[] = { -#else -static double atanlo[] = { -#endif 2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */ 3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */ 1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */ 6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */ }; -#ifdef __STDC__ static const double aT[] = { -#else -static double aT[] = { -#endif 3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */ -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */ 1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */ @@ -74,20 +62,11 @@ static double aT[] = { 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */ }; -#ifdef __STDC__ - static const double -#else - static double -#endif +static const double one = 1.0, huge = 1.0e300; -#ifdef __STDC__ - double atan(double x) -#else - double atan(x) - double x; -#endif +double atan(double x) { double w,s1,s2,z; int ix,hx,id; @@ -109,9 +88,9 @@ huge = 1.0e300; x = fabs(x); if (ix < 0x3ff30000) { /* |x| < 1.1875 */ if (ix < 0x3fe60000) { /* 7/16 <=|x|<11/16 */ - id = 0; x = (2.0*x-one)/(2.0+x); + id = 0; x = (2.0*x-one)/(2.0+x); } else { /* 11/16<=|x|< 19/16 */ - id = 1; x = (x-one)/(x+one); + id = 1; x = (x-one)/(x+one); } } else { if (ix < 0x40038000) { /* |x| < 2.4375 */ diff --git a/third-party/fdlibm/e_atan2.c b/third-party/fdlibm/s_atan2.c similarity index 90% rename from third-party/fdlibm/e_atan2.c rename to third-party/fdlibm/s_atan2.c index 4e731baa3..f9445481f 100644 --- a/third-party/fdlibm/e_atan2.c +++ b/third-party/fdlibm/s_atan2.c @@ -12,10 +12,10 @@ * */ -/* __ieee754_atan2(y,x) +/* atan2(y,x) * Method : * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). - * 2. Reduce x to positive by (if x and y are unexceptional): + * 2. Reduce x to positive by (if x and y are unexceptional): * ARG (x+iy) = arctan(y/x) ... if x > 0, * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, * @@ -33,19 +33,15 @@ * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2; * * 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 + * 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 "fdlibm.h" -#ifdef __STDC__ -static const double -#else -static double -#endif +static const double tiny = 1.0e-300, zero = 0.0, pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */ @@ -53,12 +49,7 @@ pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */ pi = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */ pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */ -#ifdef __STDC__ - double __ieee754_atan2(double y, double x) -#else - double __ieee754_atan2(y,x) - double y,x; -#endif +double atan2(double y, double x) { double z; int k,m,hx,hy,ix,iy; diff --git a/third-party/fdlibm/s_ceil.c b/third-party/fdlibm/s_ceil.c index af74592ed..9500f2ab4 100644 --- a/third-party/fdlibm/s_ceil.c +++ b/third-party/fdlibm/s_ceil.c @@ -22,18 +22,9 @@ #include "fdlibm.h" -#ifdef __STDC__ static const double huge = 1.0e300; -#else -static double huge = 1.0e300; -#endif -#ifdef __STDC__ - double ceil(double x) -#else - double ceil(x) - double x; -#endif +double ceil(double x) { int i0,i1,j0; unsigned i,j; @@ -43,7 +34,7 @@ static double huge = 1.0e300; if(j0<20) { if(j0<0) { /* raise inexact if x != 0 */ if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0<0) {i0=0x80000000;i1=0;} + if(i0<0) {i0=0x80000000;i1=0;} else if((i0|i1)!=0) { i0=0x3ff00000;i1=0;} } } else { @@ -62,7 +53,7 @@ static double huge = 1.0e300; if((i1&i)==0) return x; /* x is integral */ if(huge+x>0.0) { /* raise inexact flag */ if(i0>0) { - if(j0==20) i0+=1; + if(j0==20) i0+=1; else { j = i1 + (1<<(52-j0)); if(j 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 + * 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 "fdlibm.h" -#ifdef __STDC__ static const double -#else -static double -#endif one = 1.0, halF[2] = {0.5,-0.5,}, huge = 1.0e+300, @@ -98,12 +94,7 @@ P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ -#ifdef __STDC__ - double __ieee754_exp(double x) /* default IEEE double exp */ -#else - double __ieee754_exp(x) /* default IEEE double exp */ - double x; -#endif +double exp(double x) /* default IEEE double exp */ { double y,hi,lo,c,t; int k = 0,xsb; @@ -125,7 +116,7 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ } /* argument reduction */ - if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ + 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 { @@ -135,7 +126,7 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ lo = t*ln2LO[0]; } x = hi - lo; - } + } else if(hx < 0x3e300000) { /* when |x|<2**-28 */ if(huge+x>one) return one+x;/* trigger inexact */ } @@ -144,7 +135,7 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ /* 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); + if(k==0) return one-((x*c)/(c-2.0)-x); else y = one-((lo-(x*c)/(2.0-c))-hi); if(k >= -1021) { __HI(y) += (k<<20); /* add k to y's exponent */ diff --git a/third-party/fdlibm/s_fabs.c b/third-party/fdlibm/s_fabs.c index 0c4dd6436..e7377d821 100644 --- a/third-party/fdlibm/s_fabs.c +++ b/third-party/fdlibm/s_fabs.c @@ -17,12 +17,7 @@ #include "fdlibm.h" -#ifdef __STDC__ - double fabs(double x) -#else - double fabs(x) - double x; -#endif +double fabs(double x) { __HI(x) &= 0x7fffffff; return x; diff --git a/third-party/fdlibm/s_finite.c b/third-party/fdlibm/s_finite.c index 42e1728d7..9f8613707 100644 --- a/third-party/fdlibm/s_finite.c +++ b/third-party/fdlibm/s_finite.c @@ -18,14 +18,9 @@ #include "fdlibm.h" -#ifdef __STDC__ - int finite(double x) -#else - int finite(x) - double x; -#endif +int finite(double x) { - int hx; + int hx; hx = __HI(x); return (unsigned)((hx&0x7fffffff)-0x7ff00000)>>31; } diff --git a/third-party/fdlibm/s_floor.c b/third-party/fdlibm/s_floor.c index b37c41ba4..2d1fb84ce 100644 --- a/third-party/fdlibm/s_floor.c +++ b/third-party/fdlibm/s_floor.c @@ -22,18 +22,9 @@ #include "fdlibm.h" -#ifdef __STDC__ static const double huge = 1.0e300; -#else -static double huge = 1.0e300; -#endif -#ifdef __STDC__ - double floor(double x) -#else - double floor(x) - double x; -#endif +double floor(double x) { int i0,i1,j0; unsigned i,j; @@ -43,7 +34,7 @@ static double huge = 1.0e300; if(j0<20) { if(j0<0) { /* raise inexact if x != 0 */ if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ - if(i0>=0) {i0=i1=0;} + if(i0>=0) {i0=i1=0;} else if(((i0&0x7fffffff)|i1)!=0) { i0=0xbff00000;i1=0;} } @@ -63,7 +54,7 @@ static double huge = 1.0e300; if((i1&i)==0) return x; /* x is integral */ if(huge+x>0.0) { /* raise inexact flag */ if(i0<0) { - if(j0==20) i0+=1; + if(j0==20) i0+=1; else { j = i1+(1<<(52-j0)); if(j>31; + hx |= (unsigned)(lx|(-lx))>>31; hx = 0x7ff00000 - hx; return ((unsigned)(hx))>>31; } diff --git a/third-party/fdlibm/s_lib_version.c b/third-party/fdlibm/s_lib_version.c deleted file mode 100644 index 40f065e50..000000000 --- a/third-party/fdlibm/s_lib_version.c +++ /dev/null @@ -1,35 +0,0 @@ - -/* @(#)s_lib_version.c 1.3 95/01/18 */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -/* - * MACRO for standards - */ - -#include "fdlibm.h" - -/* - * define and initialize _LIB_VERSION - */ -#ifdef _POSIX_MODE -_LIB_VERSION_TYPE _LIB_VERSION = _POSIX_; -#else -#ifdef _XOPEN_MODE -_LIB_VERSION_TYPE _LIB_VERSION = _XOPEN_; -#else -#ifdef _SVID3_MODE -_LIB_VERSION_TYPE _LIB_VERSION = _SVID_; -#else /* default _IEEE_MODE */ -_LIB_VERSION_TYPE _LIB_VERSION = _IEEE_; -#endif -#endif -#endif diff --git a/third-party/fdlibm/e_log.c b/third-party/fdlibm/s_log.c similarity index 82% rename from third-party/fdlibm/e_log.c rename to third-party/fdlibm/s_log.c index 3798bc802..a11c1a941 100644 --- a/third-party/fdlibm/e_log.c +++ b/third-party/fdlibm/s_log.c @@ -11,28 +11,28 @@ * ==================================================== */ -/* __ieee754_log(x) +/* log(x) * Return the logrithm of x * - * Method : - * 1. Argument Reduction: find k and f such that - * x = 2^k * (1+f), + * Method : + * 1. Argument Reduction: find k and f such that + * x = 2^k * (1+f), * where sqrt(2)/2 < 1+f < sqrt(2) . * * 2. Approximation of log(1+f). * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) * = 2s + 2/3 s**3 + 2/5 s**5 + ....., * = 2s + s*R - * We use a special Reme algorithm on [0,0.1716] to generate - * a polynomial of degree 14 to approximate R The maximum error + * We use a special Reme algorithm on [0,0.1716] to generate + * a polynomial of degree 14 to approximate R The maximum error * of this polynomial approximation is bounded by 2**-58.45. In * other words, * 2 4 6 8 10 12 14 * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s - * (the values of Lg1 to Lg7 are listed in the program) + * (the values of Lg1 to Lg7 are listed in the program) * and * | 2 14 | -58.45 - * | Lg1*s +...+Lg7*s - R(z) | <= 2 + * | Lg1*s +...+Lg7*s - R(z) | <= 2 * | | * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. * In order to guarantee error in log below 1ulp, we compute log @@ -40,14 +40,14 @@ * log(1+f) = f - s*(f - R) (if f is not too large) * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) * - * 3. Finally, log(x) = k*ln2 + log(1+f). + * 3. Finally, log(x) = k*ln2 + log(1+f). * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) - * Here ln2 is split into two floating point number: + * Here ln2 is split into two floating point number: * ln2_hi + ln2_lo, * where n*ln2_hi is always exact for |n| < 2000. * * Special cases: - * log(x) is NaN with signal if x < 0 (including -INF) ; + * log(x) is NaN with signal if x < 0 (including -INF) ; * log(+INF) is +INF; log(0) is -INF with signal; * log(NaN) is that NaN with no signal. * @@ -56,19 +56,15 @@ * 1 ulp (unit in the last place). * * 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 + * 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 "fdlibm.h" -#ifdef __STDC__ static const double -#else -static double -#endif ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ @@ -80,14 +76,9 @@ Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ -static double zero = 0.0; +static const double zero = 0.0; -#ifdef __STDC__ - double __ieee754_log(double x) -#else - double __ieee754_log(x) - double x; -#endif +double log(double x) { double hfsq,f,s,z,R,w,t1,t2,dk; int k,hx,i,j; @@ -118,14 +109,14 @@ static double zero = 0.0; if(k==0) return f-R; else {dk=(double)k; return dk*ln2_hi-((R-dk*ln2_lo)-f);} } - s = f/(2.0+f); + s = f/(2.0+f); dk = (double)k; z = s*s; i = hx-0x6147a; w = z*z; j = 0x6b851-hx; - t1= w*(Lg2+w*(Lg4+w*Lg6)); - t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); + t1= w*(Lg2+w*(Lg4+w*Lg6)); + t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); i |= j; R = t2+t1; if(i>0) { diff --git a/third-party/fdlibm/e_pow.c b/third-party/fdlibm/s_pow.c similarity index 92% rename from third-party/fdlibm/e_pow.c rename to third-party/fdlibm/s_pow.c index 89a077cca..f3981792d 100644 --- a/third-party/fdlibm/e_pow.c +++ b/third-party/fdlibm/s_pow.c @@ -10,14 +10,14 @@ * ==================================================== */ -/* __ieee754_pow(x,y) return x**y +/* pow(x,y) return x**y * * n * Method: Let x = 2 * (1+f) * 1. Compute and return log2(x) in two pieces: * log2(x) = w1 + w2, * where w1 has 53-24 = 29 bit trailing zeros. - * 2. Perform y*log2(x) = n+y' by simulating muti-precision + * 2. Perform y*log2(x) = n+y' by simulating muti-precision * arithmetic, where |y'|<=0.5. * 3. Return x**y = 2**n*exp(y'*log2) * @@ -45,23 +45,19 @@ * Accuracy: * pow(x,y) returns x**y nearly rounded. In particular * pow(integer,integer) - * always returns the correct integer provided it is + * always returns the correct integer provided it is * representable. * * 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 + * 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 "fdlibm.h" -#ifdef __STDC__ -static const double -#else -static double -#endif +static const double bp[] = {1.0, 1.5,}, dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ @@ -94,12 +90,7 @@ ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */ ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ -#ifdef __STDC__ - double __ieee754_pow(double x, double y) -#else - double __ieee754_pow(x,y) - double x, y; -#endif +double pow(double x, double y) { double z,ax,z_h,z_l,p_h,p_l; double y1,t1,t2,r,s,t,u,v,w; @@ -113,12 +104,12 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ ix = hx&0x7fffffff; iy = hy&0x7fffffff; /* y==zero: x**0 = 1 */ - if((iy|ly)==0) return one; + if((iy|ly)==0) return one; /* +-NaN return x+y */ if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) || - iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) - return x+y; + iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) + return x+y; /* determine if y is an odd int when x < 0 * yisint = 0 ... y is not an integer @@ -126,7 +117,7 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ * yisint = 2 ... y is an even int */ yisint = 0; - if(hx<0) { + if(hx<0) { if(iy>=0x43400000) yisint = 2; /* even integer y */ else if(iy>=0x3ff00000) { k = (iy>>20)-0x3ff; /* exponent */ @@ -137,11 +128,11 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ j = iy>>(20-k); if((j<<(20-k))==iy) yisint = 2-(j&1); } - } - } + } + } /* special value of y */ - if(ly==0) { + if(ly==0) { if (iy==0x7ff00000) { /* y is +-inf */ if(((ix-0x3ff00000)|lx)==0) return y - y; /* inf**+-1 is NaN */ @@ -156,7 +147,7 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ if(hy==0x40000000) return x*x; /* y is 2 */ if(hy==0x3fe00000) { /* y is 0.5 */ if(hx>=0) /* x >= +0 */ - return sqrt(x); + return sqrt(x); } } @@ -169,13 +160,13 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ if(hx<0) { if(((ix-0x3ff00000)|yisint)==0) { z = (z-z)/(z-z); /* (-1)**non-int is NaN */ - } else if(yisint==1) + } else if(yisint==1) z = -z; /* (x<0)**odd = -(|x|**odd) */ } return z; } } - + n = (hx>>31)+1; /* (x<0)**(non-int) is NaN */ @@ -193,7 +184,7 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ /* over/underflow if x is not close to one */ if(ix<0x3fefffff) return (hy<0)? s*huge*huge:s*tiny*tiny; if(ix>0x3ff00000) return (hy>0)? s*huge*huge:s*tiny*tiny; - /* now |1-x| is tiny <= 2**-20, suffice to compute + /* now |1-x| is tiny <= 2**-20, suffice to compute log(x) by x-x^2/2+x^3/3-x^4/4 */ t = ax-one; /* t has 20 trailing zeros */ w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25)); @@ -225,7 +216,7 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ __LO(s_h) = 0; /* t_h=ax+bp[k] High */ t_h = zero; - __HI(t_h)=((ix>>1)|0x20000000)+0x00080000+(k<<18); + __HI(t_h)=((ix>>1)|0x20000000)+0x00080000+(k<<18); t_l = ax - (t_h-bp[k]); s_l = v*((u-s_h*t_h)-s_h*t_l); /* compute log(ax) */ @@ -287,7 +278,7 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ n = ((n&0x000fffff)|0x00100000)>>(20-k); if(j<0) n = -n; p_h -= t; - } + } t = p_l+p_h; __LO(t) = 0; u = t*lg2_h; diff --git a/third-party/fdlibm/s_rint.c b/third-party/fdlibm/s_rint.c deleted file mode 100644 index 3095e0d38..000000000 --- a/third-party/fdlibm/s_rint.c +++ /dev/null @@ -1,84 +0,0 @@ - -/* @(#)s_rint.c 1.3 95/01/18 */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -/* - * rint(x) - * Return x rounded to integral value according to the prevailing - * rounding mode. - * Method: - * Using floating addition. - * Exception: - * Inexact flag raised if x not equal to rint(x). - */ - -#include "fdlibm.h" - -#ifdef __STDC__ -static const double -#else -static double -#endif -TWO52[2]={ - 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ - -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */ -}; - -#ifdef __STDC__ - double rint(double x) -#else - double rint(x) - double x; -#endif -{ - int i0,j0,sx; - unsigned i,i1; - double w,t; - i0 = __HI(x); - sx = (i0>>31)&1; - i1 = __LO(x); - j0 = ((i0>>20)&0x7ff)-0x3ff; - if(j0<20) { - if(j0<0) { - if(((i0&0x7fffffff)|i1)==0) return x; - i1 |= (i0&0x0fffff); - i0 &= 0xfffe0000; - i0 |= ((i1|-i1)>>12)&0x80000; - __HI(x)=i0; - w = TWO52[sx]+x; - t = w-TWO52[sx]; - i0 = __HI(t); - __HI(t) = (i0&0x7fffffff)|(sx<<31); - return t; - } else { - i = (0x000fffff)>>j0; - if(((i0&i)|i1)==0) return x; /* x is integral */ - i>>=1; - if(((i0&i)|i1)!=0) { - if(j0==19) i1 = 0x40000000; else - i0 = (i0&(~i))|((0x20000)>>j0); - } - } - } else if (j0>51) { - if(j0==0x400) return x+x; /* inf or NaN */ - else return x; /* x is integral */ - } else { - i = ((unsigned)(0xffffffff))>>(j0-20); - if((i1&i)==0) return x; /* x is integral */ - i>>=1; - if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20)); - } - __HI(x) = i0; - __LO(x) = i1; - w = TWO52[sx]+x; - return w-TWO52[sx]; -} diff --git a/third-party/fdlibm/s_scalbn.c b/third-party/fdlibm/s_scalbn.c index 329be8b89..95d51b0ef 100644 --- a/third-party/fdlibm/s_scalbn.c +++ b/third-party/fdlibm/s_scalbn.c @@ -11,31 +11,22 @@ * ==================================================== */ -/* +/* * scalbn (double x, int n) - * scalbn(x,n) returns x* 2**n computed by exponent - * manipulation rather than by actually performing an + * scalbn(x,n) returns x* 2**n computed by exponent + * manipulation rather than by actually performing an * exponentiation or a multiplication. */ #include "fdlibm.h" -#ifdef __STDC__ static const double -#else -static double -#endif -two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ -twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ +two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ +twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ huge = 1.0e+300, tiny = 1.0e-300; -#ifdef __STDC__ - double scalbn (double x, int n) -#else - double scalbn (x,n) - double x; int n; -#endif +double scalbn (double x, int n) { int k,hx,lx; hx = __HI(x); @@ -45,7 +36,7 @@ tiny = 1.0e-300; if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */ x *= two54; hx = __HI(x); - k = ((hx&0x7ff00000)>>20) - 54; + k = ((hx&0x7ff00000)>>20) - 54; if (n< -50000) return tiny*x; /*underflow*/ } if (k==0x7ff) return x+x; /* NaN or Inf */ diff --git a/third-party/fdlibm/s_significand.c b/third-party/fdlibm/s_significand.c deleted file mode 100644 index 1a2163671..000000000 --- a/third-party/fdlibm/s_significand.c +++ /dev/null @@ -1,30 +0,0 @@ - -/* @(#)s_significand.c 1.3 95/01/18 */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -/* - * significand(x) computes just - * scalb(x, (double) -ilogb(x)), - * for exercising the fraction-part(F) IEEE 754-1985 test vector. - */ - -#include "fdlibm.h" - -#ifdef __STDC__ - double significand(double x) -#else - double significand(x) - double x; -#endif -{ - return __ieee754_scalb(x,(double) -ilogb(x)); -} diff --git a/third-party/fdlibm/s_sin.c b/third-party/fdlibm/s_sin.c index 43394e577..83fa2e558 100644 --- a/third-party/fdlibm/s_sin.c +++ b/third-party/fdlibm/s_sin.c @@ -20,8 +20,8 @@ * __ieee754_rem_pio2 ... argument reduction routine * * Method. - * Let S,C and T denote the sin, cos and tan respectively on - * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 + * Let S,C and T denote the sin, cos and tan respectively on + * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 * in [-pi/4 , +pi/4], and let n = k mod 4. * We have * @@ -39,17 +39,12 @@ * trig(NaN) is that NaN; * * Accuracy: - * TRIG(x) returns trig(x) nearly rounded + * TRIG(x) returns trig(x) nearly rounded */ #include "fdlibm.h" -#ifdef __STDC__ - double sin(double x) -#else - double sin(x) - double x; -#endif +double sin(double x) { double y[2],z=0.0; int n, ix; diff --git a/third-party/fdlibm/e_sqrt.c b/third-party/fdlibm/s_sqrt.c similarity index 88% rename from third-party/fdlibm/e_sqrt.c rename to third-party/fdlibm/s_sqrt.c index 7bd8de803..6d096fb9c 100644 --- a/third-party/fdlibm/e_sqrt.c +++ b/third-party/fdlibm/s_sqrt.c @@ -11,15 +11,15 @@ * ==================================================== */ -/* __ieee754_sqrt(x) +/* sqrt(x) * Return correctly rounded sqrt. * ------------------------------------------ * | Use the hardware sqrt if you have one | * ------------------------------------------ - * Method: - * Bit by bit method using integer arithmetic. (Slow, but portable) + * Method: + * Bit by bit method using integer arithmetic. (Slow, but portable) * 1. Normalization - * Scale x to y in [1,4) with even powers of 2: + * Scale x to y in [1,4) with even powers of 2: * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then * sqrt(x) = 2^k * sqrt(y) * 2. Bit by bit computation @@ -28,9 +28,9 @@ * i+1 2 * s = 2*q , and y = 2 * ( y - q ). (1) * i i i i - * - * To compute q from q , one checks whether - * i+1 i + * + * To compute q from q , one checks whether + * i+1 i * * -(i+1) 2 * (q + 2 ) <= y. (2) @@ -45,7 +45,7 @@ * s + 2 <= y (3) * i i * - * The advantage of (3) is that s and y can be computed by + * The advantage of (3) is that s and y can be computed by * i i * the following recurrence formula: * if (3) is false @@ -57,10 +57,10 @@ * -i -(i+1) * s = s + 2 , y = y - s - 2 (5) * i+1 i i+1 i i - * - * One may easily use induction to prove (4) and (5). + * + * One may easily use induction to prove (4) and (5). * Note. Since the left hand side of (3) contain only i+2 bits, - * it does not necessary to do a full (53-bit) comparison + * it does not necessary to do a full (53-bit) comparison * in (3). * 3. Final rounding * After generating the 53 bits result, we compute one more bit. @@ -70,7 +70,7 @@ * The rounding mode can be detected by checking whether * huge + tiny is equal to huge, and whether huge - tiny is * equal to huge for some floating point number "huge" and "tiny". - * + * * Special cases: * sqrt(+-0) = +-0 ... exact * sqrt(inf) = inf @@ -83,21 +83,14 @@ #include "fdlibm.h" -#ifdef __STDC__ -static const double one = 1.0, tiny=1.0e-300; -#else -static double one = 1.0, tiny=1.0e-300; -#endif +static const double +one = 1.0, +tiny = 1.0e-300; -#ifdef __STDC__ - double __ieee754_sqrt(double x) -#else - double __ieee754_sqrt(x) - double x; -#endif +double sqrt(double x) { double z; - int sign = (int)0x80000000; + int sign = (int)0x80000000; unsigned r,t1,s1,ix1,q1; int ix0,s0,q,m,t,i; @@ -108,7 +101,7 @@ static double one = 1.0, tiny=1.0e-300; if((ix0&0x7ff00000)==0x7ff00000) { return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf sqrt(-inf)=sNaN */ - } + } /* take care of zero */ if(ix0<=0) { if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */ @@ -142,12 +135,12 @@ static double one = 1.0, tiny=1.0e-300; r = 0x00200000; /* r = moving bit from right to left */ while(r!=0) { - t = s0+r; - if(t<=ix0) { - s0 = t+r; - ix0 -= t; - q += r; - } + t = s0+r; + if(t<=ix0) { + s0 = t+r; + ix0 -= t; + q += r; + } ix0 += ix0 + ((ix1&sign)>>31); ix1 += ix1; r>>=1; @@ -155,9 +148,9 @@ static double one = 1.0, tiny=1.0e-300; r = sign; while(r!=0) { - t1 = s1+r; + t1 = s1+r; t = s0; - if((tone) { if (q1==(unsigned)0xfffffffe) q+=1; - q1+=2; + q1+=2; } else q1 += (q1&1); } @@ -195,18 +188,18 @@ static double one = 1.0, tiny=1.0e-300; /* Other methods (use floating-point arithmetic) ------------- -(This is a copy of a drafted paper by Prof W. Kahan +(This is a copy of a drafted paper by Prof W. Kahan and K.C. Ng, written in May, 1986) - Two algorithms are given here to implement sqrt(x) + Two algorithms are given here to implement sqrt(x) (IEEE double precision arithmetic) in software. Both supply sqrt(x) correctly rounded. The first algorithm (in Section A) uses newton iterations and involves four divisions. The second one uses reciproot iterations to avoid division, but requires more multiplications. Both algorithms need the ability - to chop results of arithmetic operations instead of round them, + to chop results of arithmetic operations instead of round them, and the INEXACT flag to indicate when an arithmetic operation - is executed exactly with no roundoff error, all part of the + is executed exactly with no roundoff error, all part of the standard (IEEE 754-1985). The ability to perform shift, add, subtract and logical AND operations upon 32-bit words is needed too, though not part of the standard. @@ -216,7 +209,7 @@ A. sqrt(x) by Newton Iteration (1) Initial approximation Let x0 and x1 be the leading and the trailing 32-bit words of - a floating point number x (in IEEE double format) respectively + a floating point number x (in IEEE double format) respectively 1 11 52 ...widths ------------------------------------------------------ @@ -224,7 +217,7 @@ A. sqrt(x) by Newton Iteration ------------------------------------------------------ msb lsb msb lsb ...order - + ------------------------ ------------------------ x0: |s| e | f1 | x1: | f2 | ------------------------ ------------------------ @@ -249,7 +242,7 @@ A. sqrt(x) by Newton Iteration (2) Iterative refinement - Apply Heron's rule three times to y, we have y approximates + Apply Heron's rule three times to y, we have y approximates sqrt(x) to within 1 ulp (Unit in the Last Place): y := (y+x/y)/2 ... almost 17 sig. bits @@ -274,12 +267,12 @@ A. sqrt(x) by Newton Iteration it requires more multiplications and additions. Also x must be scaled in advance to avoid spurious overflow in evaluating the expression 3y*y+x. Hence it is not recommended uless division - is slow. If division is very slow, then one should use the + is slow. If division is very slow, then one should use the reciproot algorithm given in section B. (3) Final adjustment - By twiddling y's last bit it is possible to force y to be + By twiddling y's last bit it is possible to force y to be correctly rounded according to the prevailing rounding mode as follows. Let r and i be copies of the rounding mode and inexact flag before entering the square root program. Also we @@ -310,7 +303,7 @@ A. sqrt(x) by Newton Iteration I := i; ... restore inexact flag R := r; ... restore rounded mode return sqrt(x):=y. - + (4) Special cases Square root of +inf, +-0, or NaN is itself; @@ -329,7 +322,7 @@ B. sqrt(x) by Reciproot Iteration k := 0x5fe80000 - (x0>>1); y0:= k - T2[63&(k>>14)]. ... y ~ 1/sqrt(x) to 7.8 bits - Here k is a 32-bit integer and T2[] is an integer array + Here k is a 32-bit integer and T2[] is an integer array containing correction terms. Now magically the floating value of y (y's leading 32-bit word is y0, the value of its trailing word y1 is set to zero) approximates 1/sqrt(x) @@ -350,9 +343,9 @@ B. sqrt(x) by Reciproot Iteration Apply Reciproot iteration three times to y and multiply the result by x to get an approximation z that matches sqrt(x) - to about 1 ulp. To be exact, we will have + to about 1 ulp. To be exact, we will have -1ulp < sqrt(x)-z<1.0625ulp. - + ... set rounding mode to Round-to-nearest y := y*(1.5-0.5*x*y*y) ... almost 15 sig. bits to 1/sqrt(x) y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x) @@ -361,7 +354,7 @@ B. sqrt(x) by Reciproot Iteration z := z + 0.5*z*(1-z*y) ... about 1 ulp to sqrt(x) Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that - (a) the term z*y in the final iteration is always less than 1; + (a) the term z*y in the final iteration is always less than 1; (b) the error in the final result is biased upward so that -1 ulp < sqrt(x) - z < 1.0625 ulp instead of |sqrt(x)-z|<1.03125ulp. @@ -408,27 +401,27 @@ B. sqrt(x) by Reciproot Iteration I := 1; ... Raise Inexact flag: z is not exact else { j := 1 - [(x0>>20)&1] ... j = logb(x) mod 2 - k := z1 >> 26; ... get z's 25-th and 26-th + k := z1 >> 26; ... get z's 25-th and 26-th fraction bits I := i or (k&j) or ((k&(j+j+1))!=(x1&3)); } R:= r ... restore rounded mode return sqrt(x):=z. - If multiplication is cheaper then the foregoing red tape, the + If multiplication is cheaper then the foregoing red tape, the Inexact flag can be evaluated by I := i; I := (z*z!=x) or I. - Note that z*z can overwrite I; this value must be sensed if it is + Note that z*z can overwrite I; this value must be sensed if it is True. Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be zero. -------------------- - z1: | f2 | + z1: | f2 | -------------------- bit 31 bit 0 @@ -445,7 +438,6 @@ B. sqrt(x) by Reciproot Iteration 11 01 even ------------------------------------------------- - (4) Special cases (see (4) of Section A). - + (4) Special cases (see (4) of Section A). + */ - diff --git a/third-party/fdlibm/s_tan.c b/third-party/fdlibm/s_tan.c index 1f5564bce..55fd9f3f3 100644 --- a/third-party/fdlibm/s_tan.c +++ b/third-party/fdlibm/s_tan.c @@ -19,8 +19,8 @@ * __ieee754_rem_pio2 ... argument reduction routine * * Method. - * Let S,C and T denote the sin, cos and tan respectively on - * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 + * Let S,C and T denote the sin, cos and tan respectively on + * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 * in [-pi/4 , +pi/4], and let n = k mod 4. * We have * @@ -38,17 +38,12 @@ * trig(NaN) is that NaN; * * Accuracy: - * TRIG(x) returns trig(x) nearly rounded + * TRIG(x) returns trig(x) nearly rounded */ #include "fdlibm.h" -#ifdef __STDC__ - double tan(double x) -#else - double tan(x) - double x; -#endif +double tan(double x) { double y[2],z=0.0; int n, ix; diff --git a/third-party/fdlibm/s_tanh.c b/third-party/fdlibm/s_tanh.c deleted file mode 100644 index 7d77c2eac..000000000 --- a/third-party/fdlibm/s_tanh.c +++ /dev/null @@ -1,82 +0,0 @@ - -/* @(#)s_tanh.c 1.3 95/01/18 */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -/* Tanh(x) - * Return the Hyperbolic Tangent of x - * - * Method : - * x -x - * e - e - * 0. tanh(x) is defined to be ----------- - * x -x - * e + e - * 1. reduce x to non-negative by tanh(-x) = -tanh(x). - * 2. 0 <= x <= 2**-55 : tanh(x) := x*(one+x) - * -t - * 2**-55 < x <= 1 : tanh(x) := -----; t = expm1(-2x) - * t + 2 - * 2 - * 1 <= x <= 22.0 : tanh(x) := 1- ----- ; t=expm1(2x) - * t + 2 - * 22.0 < x <= INF : tanh(x) := 1. - * - * Special cases: - * tanh(NaN) is NaN; - * only tanh(0)=0 is exact for finite argument. - */ - -#include "fdlibm.h" - -#ifdef __STDC__ -static const double one=1.0, two=2.0, tiny = 1.0e-300; -#else -static double one=1.0, two=2.0, tiny = 1.0e-300; -#endif - -#ifdef __STDC__ - double tanh(double x) -#else - double tanh(x) - double x; -#endif -{ - double t,z; - int jx,ix; - - /* High word of |x|. */ - jx = __HI(x); - ix = jx&0x7fffffff; - - /* x is INF or NaN */ - if(ix>=0x7ff00000) { - if (jx>=0) return one/x+one; /* tanh(+-inf)=+-1 */ - else return one/x-one; /* tanh(NaN) = NaN */ - } - - /* |x| < 22 */ - if (ix < 0x40360000) { /* |x|<22 */ - if (ix<0x3c800000) /* |x|<2**-55 */ - return x*(one+x); /* tanh(small) = small */ - if (ix>=0x3ff00000) { /* |x|>=1 */ - t = expm1(two*fabs(x)); - z = one - two/(t+two); - } else { - t = expm1(-two*fabs(x)); - z= -t/(t+two); - } - /* |x| > 22, return +-1 */ - } else { - z = one - tiny; /* raised inexact flag */ - } - return (jx>=0)? z: -z; -} diff --git a/third-party/fdlibm/w_acos.c b/third-party/fdlibm/w_acos.c deleted file mode 100644 index e463eaf9c..000000000 --- a/third-party/fdlibm/w_acos.c +++ /dev/null @@ -1,39 +0,0 @@ - -/* @(#)w_acos.c 1.3 95/01/18 */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -/* - * wrap_acos(x) - */ - -#include "fdlibm.h" - - -#ifdef __STDC__ - double acos(double x) /* wrapper acos */ -#else - double acos(x) /* wrapper acos */ - double x; -#endif -{ -#ifdef _IEEE_LIBM - return __ieee754_acos(x); -#else - double z; - z = __ieee754_acos(x); - if(_LIB_VERSION == _IEEE_ || isnan(x)) return z; - if(fabs(x)>1.0) { - return __kernel_standard(x,x,1); /* acos(|x|>1) */ - } else - return z; -#endif -} diff --git a/third-party/fdlibm/w_asin.c b/third-party/fdlibm/w_asin.c deleted file mode 100644 index e8182857c..000000000 --- a/third-party/fdlibm/w_asin.c +++ /dev/null @@ -1,41 +0,0 @@ - -/* @(#)w_asin.c 1.3 95/01/18 */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - * - */ - -/* - * wrapper asin(x) - */ - - -#include "fdlibm.h" - - -#ifdef __STDC__ - double asin(double x) /* wrapper asin */ -#else - double asin(x) /* wrapper asin */ - double x; -#endif -{ -#ifdef _IEEE_LIBM - return __ieee754_asin(x); -#else - double z; - z = __ieee754_asin(x); - if(_LIB_VERSION == _IEEE_ || isnan(x)) return z; - if(fabs(x)>1.0) { - return __kernel_standard(x,x,2); /* asin(|x|>1) */ - } else - return z; -#endif -} diff --git a/third-party/fdlibm/w_atan2.c b/third-party/fdlibm/w_atan2.c deleted file mode 100644 index 80ad39b35..000000000 --- a/third-party/fdlibm/w_atan2.c +++ /dev/null @@ -1,40 +0,0 @@ - -/* @(#)w_atan2.c 1.3 95/01/18 */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - * - */ - -/* - * wrapper atan2(y,x) - */ - -#include "fdlibm.h" - - -#ifdef __STDC__ - double atan2(double y, double x) /* wrapper atan2 */ -#else - double atan2(y,x) /* wrapper atan2 */ - double y,x; -#endif -{ -#ifdef _IEEE_LIBM - return __ieee754_atan2(y,x); -#else - double z; - z = __ieee754_atan2(y,x); - if(_LIB_VERSION == _IEEE_||isnan(x)||isnan(y)) return z; - if(x==0.0&&y==0.0) { - return __kernel_standard(y,x,3); /* atan2(+-0,+-0) */ - } else - return z; -#endif -} diff --git a/third-party/fdlibm/w_exp.c b/third-party/fdlibm/w_exp.c deleted file mode 100644 index 7819ca133..000000000 --- a/third-party/fdlibm/w_exp.c +++ /dev/null @@ -1,48 +0,0 @@ - -/* @(#)w_exp.c 1.4 04/04/22 */ -/* - * ==================================================== - * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. - * - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -/* - * wrapper exp(x) - */ - -#include "fdlibm.h" - -#ifdef __STDC__ -static const double -#else -static double -#endif -o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ -u_threshold= -7.45133219101941108420e+02; /* 0xc0874910, 0xD52D3051 */ - -#ifdef __STDC__ - double exp(double x) /* wrapper exp */ -#else - double exp(x) /* wrapper exp */ - double x; -#endif -{ -#ifdef _IEEE_LIBM - return __ieee754_exp(x); -#else - double z; - z = __ieee754_exp(x); - if(_LIB_VERSION == _IEEE_) return z; - if(finite(x)) { - if(x>o_threshold) - return __kernel_standard(x,x,6); /* exp overflow */ - else if(x 0.0) return z; - if(x==0.0) - return __kernel_standard(x,x,16); /* log(0) */ - else - return __kernel_standard(x,x,17); /* log(x<0) */ -#endif -} diff --git a/third-party/fdlibm/w_pow.c b/third-party/fdlibm/w_pow.c deleted file mode 100644 index ea548e752..000000000 --- a/third-party/fdlibm/w_pow.c +++ /dev/null @@ -1,59 +0,0 @@ - -/* @(#)w_pow.c 1.3 95/01/18 */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -/* - * wrapper pow(x,y) return x**y - */ - -#include "fdlibm.h" - - -#ifdef __STDC__ - double pow(double x, double y) /* wrapper pow */ -#else - double pow(x,y) /* wrapper pow */ - double x,y; -#endif -{ -#ifdef _IEEE_LIBM - return __ieee754_pow(x,y); -#else - double z; - z=__ieee754_pow(x,y); - if(_LIB_VERSION == _IEEE_|| isnan(y)) return z; - if(isnan(x)) { - if(y==0.0) - return __kernel_standard(x,y,42); /* pow(NaN,0.0) */ - else - return z; - } - if(x==0.0){ - if(y==0.0) - return __kernel_standard(x,y,20); /* pow(0.0,0.0) */ - if(finite(y)&&y<0.0) - return __kernel_standard(x,y,23); /* pow(0.0,negative) */ - return z; - } - if(!finite(z)) { - if(finite(x)&&finite(y)) { - if(isnan(z)) - return __kernel_standard(x,y,24); /* pow neg**non-int */ - else - return __kernel_standard(x,y,21); /* pow overflow */ - } - } - if(z==0.0&&finite(x)&&finite(y)) - return __kernel_standard(x,y,22); /* pow underflow */ - return z; -#endif -} diff --git a/third-party/fdlibm/w_sqrt.c b/third-party/fdlibm/w_sqrt.c deleted file mode 100644 index 4dd589e25..000000000 --- a/third-party/fdlibm/w_sqrt.c +++ /dev/null @@ -1,38 +0,0 @@ - -/* @(#)w_sqrt.c 1.3 95/01/18 */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - -/* - * wrapper sqrt(x) - */ - -#include "fdlibm.h" - -#ifdef __STDC__ - double sqrt(double x) /* wrapper sqrt */ -#else - double sqrt(x) /* wrapper sqrt */ - double x; -#endif -{ -#ifdef _IEEE_LIBM - return __ieee754_sqrt(x); -#else - double z; - z = __ieee754_sqrt(x); - if(_LIB_VERSION == _IEEE_ || isnan(x)) return z; - if(x<0.0) { - return __kernel_standard(x,x,26); /* sqrt(negative) */ - } else - return z; -#endif -}