Added the remaining long double <math.h> routines

Author: ZERICO2005

src/libc/acosl.c

@@ -0,0 +1,21 @@
+
+
+#include <errno.h>
+#include <math.h>
+#include <stdbool.h>
+#include "__float64_constants.h"
+
+/**
+ * @remarks Minimum ulp:
+ * ulp of -13337 at +0x1.ffffffe84b398p-1 assuming ideal atanl
+ */
+long double acosl(long double x) {
+	if (x < 0.5L) {
+		return F64_PI2 - asinl(x);
+	}
+	if (x <= 1.0L) {
+		return atanl((sqrtl(1.0L - x) * sqrtl(1.0L + x)) / x);
+	}
+	errno = EDOM;
+	return __builtin_nanl("");
+}

src/libc/asinl.c

@@ -0,0 +1,38 @@
+/* Copyright (c) 2000-2008 Zilog, Inc. */
+
+/**
+ * asin(arg) return the arcsin of arg.
+ * arctan is called after appropriate range reduction.
+ */
+
+#include <errno.h>
+#include <math.h>
+#include <stdbool.h>
+#include "__float64_constants.h"
+
+/**
+ * @remarks Minimum ulp:
+ * ulp of -805 at -0x1.ffffffc97db1bp-1 assuming ideal atanl
+ */
+long double asinl(long double arg) {
+    bool arg_sign;
+    long double temp;
+    arg_sign = signbit(arg);
+    arg = fabsl(arg);
+
+    if(arg > 1.0L) {
+        errno = EDOM;
+        return 0.0L;
+    }
+
+    temp = sqrtl(1.0L - arg * arg);
+    if(arg > 0.7L) {
+        temp = F64_PI2 - atanl(temp / arg);
+    } else {
+        temp = atanl(arg / temp);
+    }
+    if (arg_sign) {
+        temp = -temp;
+    }
+    return temp;
+}

src/libc/cbrtf.c

@@ -1,8 +1,18 @@
 #include <math.h>
 
+/**
+ * @remarks Minimum ulp:
+ * ulp of -1  at +0x1.000006p-4   with ideal expf logf (1/16 < |x| < 16)
+ * ulp of -5  at +0x1.cb149cp-4   with current expf logf (1/16 < |x| < 16)
+ * ulp of -70 at +0x1.e9b19ep-100 with current expf logf (1.0e-30 < |x| < 1.0e+30)
+ * ulp of +81 at +0x1.e40000p-142 with current expf logf (entire range)
+ */
 float cbrtf(float x)
 {
-    return copysignf(powf(fabsf(x), 1.f/3.f), x);
+    if (x == 0.0f) {
+        return x;
+    }
+    return copysignf(expf(logf(fabsf(x)) / 3.0f), x);
 }
 
 double cbrt(double) __attribute__((alias("cbrtf")));

src/libc/cbrtl.c

@@ -0,0 +1,13 @@
+#include <math.h>
+
+/**
+ * @remarks Minimum ulp:
+ * ulp of -257 at +0x1.ffe3f201d9f88p-808 assuming ideal logl expl
+ */
+long double cbrtl(long double x)
+{
+    if (x == 0.0L) {
+        return x;
+    }
+    return copysignl(expl(logl(fabsl(x)) / 3.0L), x);
+}

src/libc/coshf.c

@@ -1,45 +1,16 @@
-/************************************************************************/
-/*                                                                      */
-/*                      Copyright (C)1987-2008 by                       */
-/*                             Zilog, Inc.                              */
-/*                                                                      */
-/*                         San Jose, California                         */
-/*                                                                      */
-/************************************************************************/
-/*
-    sinh(arg) returns the hyperbolic sine of its floating-
-    point argument.
-
-    The exponential function is called for arguments
-    greater in magnitude than 0.5.
-
-    A series is used for arguments smaller in magnitude than 0.5.
-    The coefficients are #2029 from Hart & Cheney. (20.36D)
-
-    cosh(arg) is computed from the exponential function for
-    all arguments.
-*/
-
 #include <math.h>
 
 /**
  * @remarks Minimum ulp:
  * ulp of -1  at +0x1.39a3fep-12 with ideal expf (|x| < 21.0f)
  * ulp of -18 at +0x1.0a049cp+4  with current expf (|x| < 21.0f)
  */
-float _coshf_c(float arg) {
-    float val;
-    arg = fabsf(arg);
-
-    val = expf(arg);
-
-    if (arg > 21.0f) {
-        return val / 2.0f;
+float _coshf_c(float x) {
+    x = fabsf(x);
+    if (x < 88.5f) {
+        return (expf(x) + expf(-x)) / 2.0f;
     }
-
-    val += expf(-arg);
-    val /= 2.0f;
-    return val;
+    return expf(x - (float)M_LN2);
 }
 
 double _cosh_c(double) __attribute__((alias("_coshf_c")));

src/libc/coshl.c

@@ -0,0 +1,15 @@
+#include <math.h>
+#include "__float64_constants.h"
+
+/**
+ * @remarks Minimum ulp:
+ * ulp of -1 at -0x1.f9134e7749c00p+0 with ideal expl (|x| < 709.0)
+ * ulp of +496 at +0x1.62e4249341dc8p+9 with ideal expl (|x| >= 709.0)
+ */
+long double coshl(long double x) {
+    x = fabsl(x);
+    if (x < 709.0L) {
+        return (expl(x) + expl(-x)) / 2.0L;
+    }
+    return expl(x - F64_LN2);
+}

src/libc/exp2f.c

@@ -2,7 +2,8 @@
 
 float exp2f(float x)
 {
-    return expf(x * (float)M_LOG2E);
+    // this can be made much more accurately using frexp and ldexp
+    return expf(x * (float)M_LN2);
 }
 
 double exp2(double) __attribute__((alias("exp2f")));

src/libc/exp2l.c

@@ -0,0 +1,8 @@
+#include <math.h>
+#include "__float64_constants.h"
+
+long double exp2l(long double x)
+{
+    // this can be made much more accurately using frexp and ldexp
+    return expl(x * F64_LN2);
+}

src/libc/include/__math_def.h

@@ -172,9 +172,9 @@ double      fmod(double, double);
 float       fmodf(float, float);
 long double fmodl(long double, long double);
 
-double      frexp(double, int *);
-float       frexpf(float, int *);
-long double frexpl(long double, int *);
+double      frexp(double, int *) __attribute__((nonnull(2)));
+float       frexpf(float, int *) __attribute__((nonnull(2)));
+long double frexpl(long double, int *) __attribute__((nonnull(2)));
 
 double      hypot(double, double);
 float       hypotf(float, float);
@@ -268,9 +268,9 @@ double      remainder(double, double);
 float       remainderf(float, float);
 long double remainderl(long double, long double);
 
-double      remquo(double, double, int *);
-float       remquof(float, float, int *);
-long double remquol(long double, long double, int *);
+double      remquo(double, double, int *) __attribute__((nonnull(3)));
+float       remquof(float, float, int *) __attribute__((nonnull(3)));
+long double remquol(long double, long double, int *) __attribute__((nonnull(3)));
 
 double      rint(double);
 float       rintf(float);

src/libc/include/math.h

@@ -33,7 +33,6 @@
     float: _isnanf \
 )(x))
 
-
 #define isinf(x) ((int)_Generic(__math_promote(x), \
     long double: _isinfl, \
     default: _isinff, \