added expm1l logp1l and improved acoshl atanhl

Author: ZERICO2005

src/libc/acoshl.c

@@ -1,9 +1,6 @@
 #include <math.h>
 #include "__float64_constants.h"
 
-/* enable when log1pl is implemented */
-#if 0
-
 /**
  * @remarks Minimum relative precision of:
  * 2^-42.37 at +1.000000358e+00 with ideal log1pl
@@ -15,14 +12,3 @@ long double acoshl(long double x) {
     }
     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/atanhl.c

@@ -1,23 +1,13 @@
 #include <math.h>
 
-/* enable when log1pl is implemented */
-#if 0
-
 /**
- * @remarks Relative precision of:
- * 2^-51 at +1.277396381e-01 with ideal log1pl
+ * @remarks Minimum ulp:
+ * ulp of -2    at +0x1.cf33c905a6324p-4 with ideal log1pl
+ * ulp of +4098 at +0x1.ffff1c95280b8p-15 with current log1pl and ideal logl
+ * ulp of +4135 at -0x1.f319bfec310c4p-14 with current log1pl and logl
  */
 long double atanhl(long double arg) {
     long double x = fabsl(arg);
     x = 0.5L * log1pl((2.0L * x) / (1.0L - x));
     return copysignl(x, arg);
 }
-
-#else
-
-long double atanhl(long double x)
-{
-    return 0.5L * logl((1.0L + x) / (1.0L - x));
-}
-
-#endif

src/libc/expm1l.c

@@ -0,0 +1,14 @@
+#include <math.h>
+#include "__float64_constants.h"
+
+/**
+ * @remarks Approximate minimum ulp:
+ * ulp of -4098 at +0x1.3683c27b9c856p-13 with ideal expl
+ */
+long double expm1l(long double x) {
+    if (fabsl(x) < 0x1.0p-13L) {
+        long double x_sqr = x * x;
+        return (x_sqr * x) * F64_1_DIV_3 + 0.5L * (x_sqr) + x;
+    }
+    return expl(x) - 1.0L;
+}

src/libc/log1pl.c

@@ -0,0 +1,14 @@
+#include <math.h>
+#include "__float64_constants.h"
+
+/**
+ * @remarks Approximate minimum ulp:
+ * ulp of +4097 at +0x1.0003239899b0ap-13 with ideal expl
+ */
+long double log1pl(long double x) {
+    if (fabsl(x) < 0x1.2p-13L) {
+        long double x_sqr = x * x;
+        return (x_sqr * x) * F64_1_DIV_3 - 0.5L * (x_sqr) + x;
+    }
+    return logl(1.0L + x);
+}

test/floating_point/float64_math/src/main.cpp

@@ -46,16 +46,24 @@ static size_t run_test(void) {
     if (test_result(std::atanh(  -0.9L), -1.472219489583220230004513715944L)) { return __LINE__; }
     if (test_result(std::atan (  12.0L),  1.487655094906455389320653376989L)) { return __LINE__; }
     if (test_result(std::atan (  -1.1L), -0.832981266674431705417693561837L)) { return __LINE__; }
+
     // if (test_result(std::sinh (   5.0L),  74.20321057778875897700947199607L)) { return __LINE__; }
     // if (test_result(std::sinh (  -9.0L), -4051.541902082789960515223595898L)) { return __LINE__; }
     // if (test_result(std::cosh (   8.0L),  1490.479161252178088627715460421L)) { return __LINE__; }
     // if (test_result(std::cosh (  -2.0L),  3.762195691083631459562213477774L)) { return __LINE__; }
     // if (test_result(std::tanh (   0.8L),  0.664036770267848963684844656400L)) { return __LINE__; }
     // if (test_result(std::tanh (  -3.2L), -0.996682397839651156180968063061L)) { return __LINE__; }
+
     if (test_result(std::exp  (   6.3L),  544.5719101259290330593886677332L)) { return __LINE__; }
     if (test_result(std::exp  (  -4.2L),  0.014995576820477706211984360229L)) { return __LINE__; }
     if (test_result(std::log  (1.0e-6L), -13.81551055796427410410794872811L)) { return __LINE__; }
     if (test_result(std::log  (3.0e+8L),  19.51929303262047516353917687440L)) { return __LINE__; }
+
+    if (test_result(std::expm1 (    0.2L),  0.221402758160169833921071994640L)) { return __LINE__; }
+    if (test_result(std::expm1 (-1.0e-8L), -9.999999950000000166666666250e-9L)) { return __LINE__; }
+    if (test_result(std::log1p (   -0.2L), -0.223143551314209755766295090310L)) { return __LINE__; }
+    if (test_result(std::log1p ( 1.0e-8L),  9.999999950000000333333330833e-9L)) { return __LINE__; }
+
     if (test_result(std::erfc (  -1.3L),  1.934007944940652436603893327504L)) { return __LINE__; }
     if (test_result(std::erfc (   3.1L),  0.000011648657367199596033706164L)) { return __LINE__; }
     if (test_result(std::erf  (  -2.2L), -0.998137153702018108556548243971L)) { return __LINE__; }