Improved accuracy of libm

src/libc/__float32_constants.h

@@ -0,0 +1,32 @@
+#ifndef __FLOAT32_CONSTANTS_H
+#define __FLOAT32_CONSTANTS_H
+
+#define F32_1_DIV_3 0.33333333333333333333333333333333f /* 1/3 */
+#define F32_2_DIV_3 0.66666666666666666666666666666667f /* 2/3 */
+#define F32_1_DIV_6 0.16666666666666666666666666666667f /* 1/6 */
+
+#define F32_SQRT2         1.414213562373095048801688724209698f /* sqrt(2) */
+#define F32_SQRT3         1.732050807568877293527446341505872f /* sqrt(3) */
+#define F32_INV_SQRT2     0.70710678118654752440084436210485f  /* 1 / sqrt(2) */
+#define F32_INV_SQRT3     0.577350269189625764509148780501957f /* 1 / sqrt(3) */
+#define F32_SQRT2_PLUS_1  2.4142135623730950488016887242097f   /* sqrt(2) + 1 */
+#define F32_SQRT2_MINUS_1 0.4142135623730950488016887242097f   /* sqrt(2) - 1 */
+
+#define F32_E      2.718281828459045235360287471352662f /* e */
+#define F32_LN2    0.693147180559945309417232121458176f /* ln(2) */
+#define F32_LN10   2.302585092994045684017991454684364f /* ln(10) */
+#define F32_LOG2E  1.442695040888963407359924681001892f /* log2(e) */
+#define F32_LOG10E 0.434294481903251827651128918916605f /* log10(e) */
+
+#define F32_PI           3.141592653589793238462643383279502f /* pi */
+#define F32_PI2          1.5707963267948966192313216916398f   /* 1/2 * pi */
+#define F32_PI4          0.78539816339744830961566084581988f  /* 1/4 * pi */
+#define F32_3PI4         2.3561944901923449288469825374596f   /* 3/4 * pi */
+#define F32_INV_PI       0.318309886183790671537767526745028f /* 1 / pi */
+#define F32_INV_SQRTPI   0.564189583547756286948079451560772f /* 1 / sqrt(pi) */
+#define F32_2_DIV_SQRTPI 1.1283791670955125738961589031215f   /* 2 / sqrt(pi) */
+
+#define F32_EGAMMA 0.577215664901532860606512090082402f /* Euler-Mascheroni constant */
+#define F32_PHI    1.618033988749894848204586834365638f /* Golden Ratio */
+
+#endif /* __FLOAT32_CONSTANTS_H */

src/libc/__float64_constants.h

@@ -0,0 +1,32 @@
+#ifndef __FLOAT64_CONSTANTS_H
+#define __FLOAT64_CONSTANTS_H
+
+#define F64_1_DIV_3 0.33333333333333333333333333333333L /* 1/3 */
+#define F64_2_DIV_3 0.66666666666666666666666666666667L /* 2/3 */
+#define F64_1_DIV_6 0.16666666666666666666666666666667L /* 1/6 */
+
+#define F64_SQRT2         1.414213562373095048801688724209698L /* sqrt(2) */
+#define F64_SQRT3         1.732050807568877293527446341505872L /* sqrt(3) */
+#define F64_INV_SQRT2     0.70710678118654752440084436210485L  /* 1 / sqrt(2) */
+#define F64_INV_SQRT3     0.577350269189625764509148780501957L /* 1 / sqrt(3) */
+#define F64_SQRT2_PLUS_1  2.4142135623730950488016887242097L   /* sqrt(2) + 1 */
+#define F64_SQRT2_MINUS_1 0.4142135623730950488016887242097L   /* sqrt(2) - 1 */
+
+#define F64_E      2.718281828459045235360287471352662L /* e */
+#define F64_LN2    0.693147180559945309417232121458176L /* ln(2) */
+#define F64_LN10   2.302585092994045684017991454684364L /* ln(10) */
+#define F64_LOG2E  1.442695040888963407359924681001892L /* log2(e) */
+#define F64_LOG10E 0.434294481903251827651128918916605L /* log10(e) */
+
+#define F64_PI           3.141592653589793238462643383279502L /* pi */
+#define F64_PI2          1.5707963267948966192313216916398L   /* 1/2 * pi */
+#define F64_PI4          0.78539816339744830961566084581988L  /* 1/4 * pi */
+#define F64_3PI4         2.3561944901923449288469825374596L   /* 3/4 * pi */
+#define F64_INV_PI       0.318309886183790671537767526745028L /* 1 / pi */
+#define F64_INV_SQRTPI   0.564189583547756286948079451560772L /* 1 / sqrt(pi) */
+#define F64_2_DIV_SQRTPI 1.1283791670955125738961589031215L   /* 2 / sqrt(pi) */
+
+#define F64_EGAMMA 0.577215664901532860606512090082402L /* Euler-Mascheroni constant */
+#define F64_PHI    1.618033988749894848204586834365638L /* Golden Ratio */
+
+#endif /* __FLOAT64_CONSTANTS_H */

src/libc/acoshf.c

@@ -1,8 +1,16 @@
 #include <math.h>
 
-float acoshf(float x)
-{
-    return logf(x + sqrtf(x * x - 1));
+/**
+ * @remarks Minimum ulp:
+ * ulp of -2  at +0x1.0007e8p+0 with ideal log1pf
+ * ulp of +10 at +0x1.01c964p+0 with current log1pf
+ */
+float acoshf(float x) {
+    if (x < 0x1.0p+63f) {
+        float t = x - 1.0f;
+        return log1pf(t + sqrtf(2.0f * t + t * t));
+    }
+    return logf(x) + (float)M_LN2;
 }
 
 double acosh(double) __attribute__((alias("acoshf")));

src/libc/acoshl.c

@@ -1,6 +1,28 @@
 #include <math.h>
+#include "__float64_constants.h"
 
-long double acoshl(long double x)
-{
-    return logl(x + sqrtl(x * x - 1));
+/* enable when log1pl is implemented */
+#if 0
+
+/**
+ * @remarks Minimum relative precision of:
+ * 2^-42.37 at +1.000000358e+00 with ideal log1pl
+ */
+long double acoshl(long double x) {
+    if (x < 0x1.0p+511L) {
+        long double t = x - 1.0L;
+        return log1pl(t + sqrtl(2.0L * t + t * t));
+    }
+    return logl(x) + F64_LN2;
 }
+
+#else
+
+long double acoshl(long double x) {
+    if (x < 0x1.0p+511L) {
+        return logl(x + sqrtl(x * x - 1.0L));
+    }
+    return logl(x) + F64_LN2;
+}
+
+#endif

src/libc/asinf.c

@@ -7,31 +7,33 @@
 
 #include <errno.h>
 #include <math.h>
+#include <stdbool.h>
 
-#define pio2  1.57079632679490
-
+/**
+ * @remarks Minimum ulp:
+ * ulp of +8 at +0x1.ffe956p-1
+ */
 float _asinf_c(float arg) {
-    float sign, temp;
+    bool arg_sign;
+    float temp;
+    arg_sign = signbit(arg);
+    arg = fabsf(arg);
 
-    sign = 1.;
-    if(arg < 0) {
-        arg = -arg;
-        sign = -1.;
-    }
-
-    if(arg > 1.) {
+    if(arg > 1.0f) {
         errno = EDOM;
-        return(0.);
+        return 0.0f;
     }
 
-    temp = sqrtf(1. - arg*arg);
-    if(arg > 0.7) {
-        temp = pio2 - atanf(temp/arg);
+    temp = sqrtf(1.0f - arg*arg);
+    if(arg > 0.7f) {
+        temp = (float)M_PI_2 - atanf(temp/arg);
     } else {
         temp = atanf(arg/temp);
     }
-
-    return(sign*temp);
+    if (arg_sign) {
+        temp = -temp;
+    }
+    return temp;
 }
 
 double _asin_c(double) __attribute__((alias("_asinf_c")));

src/libc/asinhf.c

@@ -1,8 +1,20 @@
 #include <math.h>
+#include "__float32_constants.h"
 
-float asinhf(float x)
-{
-    return logf(x + sqrtf(x * x + 1));
+/**
+ * @remarks Minimum ulp:
+ * ulp of -20 at +0x1.ff2afcp-5
+ */
+float asinhf(float arg) {
+    float x = fabsf(arg);
+    if (x < 0.0703125f) {
+        x = x - (x * x * x) * F32_1_DIV_6;
+    } else if (x < 0x1.0p+63f) {
+        x = logf(x + sqrtf(x * x + 1.0f));
+    } else {
+        x = logf(x) + F32_LN2;
+    }
+    return copysignf(x, arg);
 }
 
 double asinh(double) __attribute__((alias("asinhf")));

src/libc/asinhl.c

@@ -1,6 +1,18 @@
 #include <math.h>
+#include "__float64_constants.h"
 
-long double asinhl(long double x)
-{
-    return logl(x + sqrtl(x * x + 1));
+/**
+ * @remarks Relative precision of:
+ * 2^-42.41 at +9.768006857e-04 with ideal logl
+ */
+long double asinhl(long double arg) {
+    long double x = fabsl(arg);
+    if (x < 0x1.0p-10L) {
+        x = x - (x * x * x) * F64_1_DIV_6;
+    } else if (x < 0x1.0p+511L) {
+        x = logl(x + sqrtl(x * x + 1.0L));
+    } else {
+        x = logl(x) + F64_LN2;
+    }
+    return copysignl(x, arg);
 }

src/libc/atan2f.c

@@ -1,27 +1,26 @@
 /* Copyright (c) 2000-2008 Zilog, Inc. */
 #include <math.h>
+#include "__float32_constants.h"
 
-#define pio2 1.57079632679489e0
-
-float _atan2f_c(float arg1,float arg2) {
+float _atan2f_c(float arg1, float arg2) {
     float satan(float);
 
-    if((arg1+arg2)==arg1) {
-        if(arg1 >= 0.) {
-            return(pio2);
+    if ((arg1+arg2)==arg1) {
+        if (arg1 >= 0.0f) {
+            return (F32_PI2);
         } else {
-            return(-pio2);
+            return (-F32_PI2);
         }
-    } else if(arg2 < 0.) {
-        if(arg1 >= 0.) {
-            return(pio2+pio2 - satan(-arg1/arg2));
+    } else if (arg2 < 0.0f) {
+        if(arg1 >= 0.0f) {
+            return (F32_PI - satan(-arg1/arg2));
         } else {
-            return(-pio2-pio2 + satan(arg1/arg2));
+            return (-F32_PI + satan(arg1/arg2));
         }
-    } else if(arg1 > 0) {
-        return(satan(arg1/arg2));
+    } else if (arg1 > 0.0f) {
+        return (satan(arg1/arg2));
     } else {
-        return(-satan(-arg1/arg2));
+        return (-satan(-arg1/arg2));
     }
 }
 

src/libc/atan2l.c

@@ -0,0 +1,25 @@
+/* Copyright (c) 2000-2008 Zilog, Inc. */
+#include <math.h>
+#include "__float64_constants.h"
+
+long double atan2l(long double arg1, long double arg2) {
+    long double f64_satan(long double);
+
+    if ((arg1+arg2)==arg1) {
+        if (arg1 >= 0.0L) {
+            return (F64_PI2);
+        } else {
+            return (-F64_PI2);
+        }
+    } else if (arg2 < 0.0L) {
+        if(arg1 >= 0.0L) {
+            return (F64_PI - f64_satan(-arg1/arg2));
+        } else {
+            return (-F64_PI + f64_satan(arg1/arg2));
+        }
+    } else if (arg1 > 0.0L) {
+        return (f64_satan(arg1/arg2));
+    } else {
+        return (-f64_satan(-arg1/arg2));
+    }
+}

src/libc/atanf.c

@@ -14,34 +14,36 @@
  * coefficients are #5077 from Hart & Cheney. (19.56D)
  */
 #include <math.h>
+#include "__float32_constants.h"
 
-#define sq2p1 2.41421356237309e0
-#define sq2m1 0.414213562373095e0
-#define pio2  1.57079632679489e0
-#define pio4  0.785398163397448e0
-#define p4    0.161536412982230e2
-#define p3    0.268425481955040e3
-#define p2    0.115302935154049e4
-#define p1    0.178040631643320e4
-#define p0    0.896785974036639e3
-#define q4    0.589569705084446e2
-#define q3    0.536265374031215e3
-#define q2    0.166678381488163e4
-#define q1    0.207933497444541e4
-#define q0    0.896785974036639e3
+/**
+ * @note These coefficients are for float64
+ */
+#define p4 0.161536412982230e2f
+#define p3 0.268425481955040e3f
+#define p2 0.115302935154049e4f
+#define p1 0.178040631643320e4f
+#define p0 0.896785974036639e3f
+#define q4 0.589569705084446e2f
+#define q3 0.536265374031215e3f
+#define q2 0.166678381488163e4f
+#define q1 0.207933497444541e4f
+#define q0 0.896785974036639e3f
 
 /**
  * atan makes its argument positive and
  * calls the inner routine satan.
+ *
+ * @remarks Minimum ulp:
+ * ulp of +4 at +0x1.a85846p-2
  */
-
 float _atanf_c(float arg) {
     float satan(float);
 
-    if(arg>0) {
-        return(satan(arg));
+    if (signbit(arg)) {
+        return (-satan(-arg));
     } else {
-        return(-satan(-arg));
+        return (satan(arg));
     }
 }
 
@@ -52,7 +54,6 @@ double _atan_c(double) __attribute__((alias("_atanf_c")));
  * is in and calls atan.
  */
 
-
 /**
  * xatan evaluates a series valid in the
  * range [-0.414...,+0.414...].
@@ -65,7 +66,7 @@ static float xatan(float arg) {
     argsq = arg*arg;
     value = ((((p4*argsq + p3)*argsq + p2)*argsq + p1)*argsq + p0);
     value = value/(((((argsq + q4)*argsq + q3)*argsq + q2)*argsq + q1)*argsq + q0);
-    return(value*arg);
+    return (value*arg);
 }
 
 /**
@@ -74,11 +75,15 @@ static float xatan(float arg) {
  */
 
 float satan(float arg) {
-    if(arg < sq2m1) {
-        return(xatan(arg));
-    } else if(arg > sq2p1) {
-        return(pio2 - xatan(1.0/arg));
+    if (arg < F32_SQRT2_MINUS_1) {
+        return (xatan(arg));
+    } else if (arg > F32_SQRT2_PLUS_1) {
+        if (arg > 0x1.0p+25f) {
+            /* rounds to pi/2 */
+            return F32_PI2;
+        }
+        return (F32_PI2 - xatan(1.0f/arg));
     } else {
-        return(pio4 + xatan((arg-1.0)/(arg+1.0)));
+        return (F32_PI4 + xatan((arg-1.0f)/(arg+1.0f)));
     }
 }

src/libc/atanhf.c

@@ -1,8 +1,15 @@
 #include <math.h>
 
-float atanhf(float x)
-{
-    return .5f * logf((1 + x) / (1 - x));
+/**
+ * @remarks Minimum ulp:
+ * ulp of +2  at +0x1.fe2dd0p-9 with ideal log1pf
+ * ulp of -9  at +0x1.e1e8d0p-5 with current log1pf and ideal logf
+ * ulp of +11 at +0x1.f30abcp-5 with current log1pf and logf
+ */
+float atanhf(float arg) {
+    float x = fabsf(arg);
+    x = 0.5f * log1pf((2.0f * x) / (1.0f - x));
+    return copysignf(x, arg);
 }
 
 double atanh(double) __attribute__((alias("atanhf")));