diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/math/__invtrigl.c | 42 | ||||
| -rw-r--r-- | src/math/__invtrigl.h | 63 | ||||
| -rw-r--r-- | src/math/acos.c | 1 | ||||
| -rw-r--r-- | src/math/acosl.c | 9 | ||||
| -rw-r--r-- | src/math/asinl.c | 4 | ||||
| -rw-r--r-- | src/math/atan2.c | 1 | ||||
| -rw-r--r-- | src/math/atan2l.c | 32 | ||||
| -rw-r--r-- | src/math/atanl.c | 45 | 
8 files changed, 87 insertions, 110 deletions
diff --git a/src/math/__invtrigl.c b/src/math/__invtrigl.c index a821842c..a51330e2 100644 --- a/src/math/__invtrigl.c +++ b/src/math/__invtrigl.c @@ -28,9 +28,8 @@  #include "__invtrigl.h"  #if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 -/* - * asinl() and acosl() - */ + +/* coefficients used in asinl() and acosl() */  const long double  pS0 =  1.66666666666666666631e-01L,  pS1 = -4.16313987993683104320e-01L, @@ -45,38 +44,9 @@ qS3 = -1.68285799854822427013e+00L,  qS4 =  3.90699412641738801874e-01L,  qS5 = -3.14365703596053263322e-02L; -/* - * atanl() - */ -const long double atanhi[] = { -	 4.63647609000806116202e-01L, -	 7.85398163397448309628e-01L, -	 9.82793723247329067960e-01L, -	 1.57079632679489661926e+00L, -}; - -const long double atanlo[] = { -	 1.18469937025062860669e-20L, -	-1.25413940316708300586e-20L, -	 2.55232234165405176172e-20L, -	-2.50827880633416601173e-20L, -}; - -const long double aT[] = { -	 3.33333333333333333017e-01L, -	-1.99999999999999632011e-01L, -	 1.42857142857046531280e-01L, -	-1.11111111100562372733e-01L, -	 9.09090902935647302252e-02L, -	-7.69230552476207730353e-02L, -	 6.66661718042406260546e-02L, -	-5.88158892835030888692e-02L, -	 5.25499891539726639379e-02L, -	-4.70119845393155721494e-02L, -	 4.03539201366454414072e-02L, -	-2.91303858419364158725e-02L, -	 1.24822046299269234080e-02L, -}; - +const long double pi_hi = 3.1415926535897932384626433832795L;  const long double pi_lo = -5.01655761266833202345e-20L; +const long double pio2_hi = 1.57079632679489661926L; +const long double pio2_lo = -2.50827880633416601173e-20L; +  #endif diff --git a/src/math/__invtrigl.h b/src/math/__invtrigl.h index c3ad3c49..22748b68 100644 --- a/src/math/__invtrigl.h +++ b/src/math/__invtrigl.h @@ -32,15 +32,6 @@  #define BIAS            (LDBL_MAX_EXP - 1)  #define MANH_SIZE       LDBL_MANH_SIZE -/* Approximation thresholds. */ -#define ASIN_LINEAR     (BIAS - 32)     /* 2**-32 */ -#define ACOS_CONST      (BIAS - 65)     /* 2**-65 */ -#define ATAN_CONST      (BIAS + 65)     /* 2**65 */ -#define ATAN_LINEAR     (BIAS - 32)     /* 2**-32 */ - -/* 0.95 */ -#define THRESH  ((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|LDBL_NBIT) -  /* Constants shared by the long double inverse trig functions. */  #define pS0     __pS0  #define pS1     __pS1 @@ -54,56 +45,24 @@  #define qS3     __qS3  #define qS4     __qS4  #define qS5     __qS5 -#define atanhi  __atanhi -#define atanlo  __atanlo -#define aT      __aT +#define pi_hi   __pi_hi  #define pi_lo   __pi_lo +#define pio2_hi __pio2_hi +#define pio2_lo __pio2_lo -#define pio2_hi atanhi[3] -#define pio2_lo atanlo[3] -#define pio4_hi atanhi[1] - -#ifdef STRUCT_DECLS -typedef struct longdouble { -	uint64_t mant; -	uint16_t expsign; -} LONGDOUBLE; -#else -typedef long double LONGDOUBLE; -#endif - -extern const LONGDOUBLE pS0, pS1, pS2, pS3, pS4, pS5, pS6; -extern const LONGDOUBLE qS1, qS2, qS3, qS4, qS5; -extern const LONGDOUBLE atanhi[], atanlo[], aT[]; -extern const LONGDOUBLE pi_lo; - -#ifndef STRUCT_DECLS -static inline long double -P(long double x) -{ -	return (x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 + x * \ -		(pS4 + x * (pS5 + x * pS6))))))); -} - -static inline long double -Q(long double x) -{ -	return (1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * qS5))))); -} +extern const long double pS0, pS1, pS2, pS3, pS4, pS5, pS6; +extern const long double qS1, qS2, qS3, qS4, qS5; +extern const long double pi_hi, pi_lo, pio2_hi, pio2_lo; -static inline long double -T_even(long double x) +static long double P(long double x)  { -	return (aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + x * \ -		(aT[8] + x * (aT[10] + x * aT[12])))))); +	return x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 + +		x * (pS4 + x * (pS5 + x * pS6))))));  } -static inline long double -T_odd(long double x) +static long double Q(long double x)  { -	return (aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + x * \ -		(aT[9] + x * aT[11]))))); +	return 1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * qS5))));  } -#endif  #endif diff --git a/src/math/acos.c b/src/math/acos.c index 54d266ee..0eb15bed 100644 --- a/src/math/acos.c +++ b/src/math/acos.c @@ -38,6 +38,7 @@  static const double  pi      = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */  pio2_hi = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */ +// FIXME  static const volatile double  pio2_lo = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */  static const double diff --git a/src/math/acosl.c b/src/math/acosl.c index cc565336..83857d49 100644 --- a/src/math/acosl.c +++ b/src/math/acosl.c @@ -23,9 +23,7 @@ long double acosl(long double x)  }  #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384  #include "__invtrigl.h" - -static const long double -pi = 3.14159265358979323846264338327950280e+00L; +#define ACOS_CONST      (BIAS - 65)     /* 2**-65 */  long double acosl(long double x)  { @@ -41,7 +39,8 @@ long double acosl(long double x)  			if (expsign > 0)  				return 0.0;  /* acos(1) = 0 */  			else -				return pi + 2.0 * pio2_lo;  /* acos(-1)= pi */ +				// FIXME +				return pi_hi + 2.0 * pio2_lo;  /* acos(-1)= pi */  		}  		return (x - x) / (x - x);  /* acos(|x|>1) is NaN */  	} @@ -60,7 +59,7 @@ long double acosl(long double x)  		s = sqrtl(z);  		r = p / q;  		w = r * s - pio2_lo; -		return pi - 2.0 * (s + w); +		return pi_hi - 2.0 * (s + w);  	} else {                   /* x > 0.5 */  		z = (1.0 - x) * 0.5;  		s = sqrtl(z); diff --git a/src/math/asinl.c b/src/math/asinl.c index ddd807e2..7572767c 100644 --- a/src/math/asinl.c +++ b/src/math/asinl.c @@ -24,6 +24,10 @@ long double asinl(long double x)  #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384  #include "__invtrigl.h"  static const long double huge = 1.000e+300; +static const long double pio4_hi = 7.85398163397448309628e-01L; +#define ASIN_LINEAR     (BIAS - 32)     /* 2**-32 */ +/* 0.95 */ +#define THRESH  ((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|LDBL_NBIT)  long double asinl(long double x)  { diff --git a/src/math/atan2.c b/src/math/atan2.c index 143c3834..825a38fe 100644 --- a/src/math/atan2.c +++ b/src/math/atan2.c @@ -39,6 +39,7 @@  #include "libm.h" +// FIXME  static const volatile double  tiny  = 1.0e-300;  static const double diff --git a/src/math/atan2l.c b/src/math/atan2l.c index 0fc901c8..45cbfcc1 100644 --- a/src/math/atan2l.c +++ b/src/math/atan2l.c @@ -24,10 +24,8 @@ long double atan2l(long double y, long double x)  }  #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384  #include "__invtrigl.h" -static const volatile long double -tiny = 1.0e-300; -static const long double -pi = 3.14159265358979323846264338327950280e+00L; +// FIXME: +static const volatile long double tiny = 1.0e-300;  long double atan2l(long double y, long double x)  { @@ -55,9 +53,9 @@ long double atan2l(long double y, long double x)  	if (expty==0 && ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)==0) {  		switch(m) {  		case 0: -		case 1: return y;        /* atan(+-0,+anything)=+-0 */ -		case 2: return  pi+tiny; /* atan(+0,-anything) = pi */ -		case 3: return -pi-tiny; /* atan(-0,-anything) =-pi */ +		case 1: return y;           /* atan(+-0,+anything)=+-0 */ +		case 2: return  pi_hi+tiny; /* atan(+0,-anything) = pi */ +		case 3: return -pi_hi-tiny; /* atan(-0,-anything) =-pi */  		}  	}  	/* when x = 0 */ @@ -69,15 +67,15 @@ long double atan2l(long double y, long double x)  			switch(m) {  			case 0: return  pio2_hi*0.5+tiny; /* atan(+INF,+INF) */  			case 1: return -pio2_hi*0.5-tiny; /* atan(-INF,+INF) */ -			case 2: return  1.5*pio2_hi+tiny; /*atan(+INF,-INF)*/ -			case 3: return -1.5*pio2_hi-tiny; /*atan(-INF,-INF)*/ +			case 2: return  1.5*pio2_hi+tiny; /* atan(+INF,-INF) */ +			case 3: return -1.5*pio2_hi-tiny; /* atan(-INF,-INF) */  			}  		} else {  			switch(m) { -			case 0: return  0.0;     /* atan(+...,+INF) */ -			case 1: return -0.0;     /* atan(-...,+INF) */ -			case 2: return  pi+tiny; /* atan(+...,-INF) */ -			case 3: return -pi-tiny; /* atan(-...,-INF) */ +			case 0: return  0.0;        /* atan(+...,+INF) */ +			case 1: return -0.0;        /* atan(-...,+INF) */ +			case 2: return  pi_hi+tiny; /* atan(+...,-INF) */ +			case 3: return -pi_hi-tiny; /* atan(-...,-INF) */  			}  		}  	} @@ -95,11 +93,11 @@ long double atan2l(long double y, long double x)  	else                     /* safe to do y/x */  		z = atanl(fabsl(y/x));  	switch (m) { -	case 0: return z;              /* atan(+,+) */ -	case 1: return -z;             /* atan(-,+) */ -	case 2: return pi - (z-pi_lo); /* atan(+,-) */ +	case 0: return z;               /* atan(+,+) */ +	case 1: return -z;              /* atan(-,+) */ +	case 2: return pi_hi-(z-pi_lo); /* atan(+,-) */  	default: /* case 3 */ -		return (z-pi_lo) - pi; /* atan(-,-) */ +		return (z-pi_lo)-pi_hi; /* atan(-,-) */  	}  }  #endif diff --git a/src/math/atanl.c b/src/math/atanl.c index 36072c17..33ecff0d 100644 --- a/src/math/atanl.c +++ b/src/math/atanl.c @@ -23,8 +23,53 @@ long double atanl(long double x)  }  #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384  #include "__invtrigl.h" + +#define ATAN_CONST      (BIAS + 65)     /* 2**65 */ +#define ATAN_LINEAR     (BIAS - 32)     /* 2**-32 */  static const long double huge = 1.0e300; +static const long double atanhi[] = { +	 4.63647609000806116202e-01L, +	 7.85398163397448309628e-01L, +	 9.82793723247329067960e-01L, +	 1.57079632679489661926e+00L, +}; + +static const long double atanlo[] = { +	 1.18469937025062860669e-20L, +	-1.25413940316708300586e-20L, +	 2.55232234165405176172e-20L, +	-2.50827880633416601173e-20L, +}; + +static const long double aT[] = { +	 3.33333333333333333017e-01L, +	-1.99999999999999632011e-01L, +	 1.42857142857046531280e-01L, +	-1.11111111100562372733e-01L, +	 9.09090902935647302252e-02L, +	-7.69230552476207730353e-02L, +	 6.66661718042406260546e-02L, +	-5.88158892835030888692e-02L, +	 5.25499891539726639379e-02L, +	-4.70119845393155721494e-02L, +	 4.03539201366454414072e-02L, +	-2.91303858419364158725e-02L, +	 1.24822046299269234080e-02L, +}; + +static long double T_even(long double x) +{ +	return aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + +		x * (aT[8] + x * (aT[10] + x * aT[12]))))); +} + +static long double T_odd(long double x) +{ +	return aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + +		x * (aT[9] + x * aT[11])))); +} +  long double atanl(long double x)  {  	union IEEEl2bits u;  | 
