Fixes to asinh(), atanh(), and sinh(): The 'basic' expansions all exhibited bad cancellation errors near zero (<= 1E-10). This can be easily seen e.g. with x = 1E-30, the results are all zero though they should be extremely close to x. The cutoff values (1E-5, 1E-9) are chosen so that the next term in the Taylor series is negligible (for double). Functions could maybe be optimized to use only first term (x) and a smaller cutoff, just bigger than where the cancellation occurs.

3 files changed, 49 insertions, 5 deletions

diff --git a/nuttx/libc/math/lib_asinh.c b/nuttx/libc/math/lib_asinh.c
index 7b573fcdd..5a30ac1aa 100644
--- a/nuttx/libc/math/lib_asinh.c
+++ b/nuttx/libc/math/lib_asinh.c
@@ -49,6 +49,20 @@
 #ifdef CONFIG_HAVE_DOUBLE
 double asinh(double x)
 {
-  return log(x + sqrt(x * x + 1));
+  double y;
+  double z = x * x;
+
+  if (z < 1E-9)
+    {
+      /* x - 1/6 * x^3 + 3/40 * x^5 - 5/112 * x^7 + ... */
+
+      y = x * (1 - z / 6.0);
+    }
+  else
+    {
+      y = log(x + sqrt(z + 1));
+    }
+
+  return y;
 }
 #endif
diff --git a/nuttx/libc/math/lib_atanh.c b/nuttx/libc/math/lib_atanh.c
index 9932a5e9b..396e80852 100644
--- a/nuttx/libc/math/lib_atanh.c
+++ b/nuttx/libc/math/lib_atanh.c
@@ -49,6 +49,21 @@
 #ifdef CONFIG_HAVE_DOUBLE
 double atanh(double x)
 {
-  return 0.5 * log((1 + x) / (1 - x));
+  double y;
+
+  if (fabs(x) < 1E-5)
+    {
+      double z = x * x;
+
+      /* x + 1/3 * x^3 + 1/5 * x^5 + 1/7 * x^7 + ... */
+
+      y = x * (1 + z / 3.0);
+    }
+  else
+    {
+      y = log((1 + x) / (1 - x)) / 2.0;
+    }
+
+  return y;
 }
 #endif
diff --git a/nuttx/libc/math/lib_sinh.c b/nuttx/libc/math/lib_sinh.c
index da79e7cfd..714a52816 100644
--- a/nuttx/libc/math/lib_sinh.c
+++ b/nuttx/libc/math/lib_sinh.c
@@ -3,7 +3,7 @@
  *
  * This file is a part of NuttX:
  *
- *   Copyright (C) 2012 Gregory Nutt. All rights reserved.
+ *   Copyright (C) 2012, 2015 Gregory Nutt. All rights reserved.
  *   Ported by: Darcy Gong
  *
  * It derives from the Rhombs OS math library by Nick Johnson which has
@@ -41,7 +41,22 @@
 #ifdef CONFIG_HAVE_DOUBLE
 double sinh(double x)
 {
-  x = exp(x);
-  return ((x - (1.0 / x)) / 2.0);
+  double y;
+  double z;
+
+  if (fabs(x) < 1E-5)
+    {
+      /* x + 1/3! * x^3 + 1/5! * x^5 + 1/7! * x^7 + ... */
+
+      z = x * x;
+      y = x * (1.0 + z / 6.0);
+    }
+  else
+    {
+      z = exp(x);
+      y = (z - 1.0 / z) / 2.0;
+    }
+
+  return y;
 }
 #endif