Fixed #1008: loss of precision issues in arcsin, arcsinh, arctan, arctanh

The complex-valued arc* functions (that -> 0 for z -> 0) have loss of
precision issues for small arguments. This patch addresses this
by switching to a series expansion (in Horner form) in this regime.

Tests included.

2 files changed, 186 insertions, 26 deletions

diff --git a/numpy/core/src/umath_funcs.inc.src b/numpy/core/src/umath_funcs.inc.src
index ef87bd10a..161e12771 100644
--- a/numpy/core/src/umath_funcs.inc.src
+++ b/numpy/core/src/umath_funcs.inc.src
@@ -96,6 +96,17 @@ npy_Object@Kind@(PyObject *i1, PyObject *i2)
    #c=f,,l#
 */
 
+/* Perform the operation  result := 1 + coef * x * result,
+ * with real coefficient `coef`.
+ */
+#define SERIES_HORNER_TERM@c@(result, x, coef)                  \
+    do {                                                        \
+        nc_prod@c@((result), (x), (result));                    \
+        (result)->real *= (coef);                               \
+        (result)->imag *= (coef);                               \
+        nc_sum@c@((result), &nc_1@c@, (result));                \
+    } while(0)
+
 /* constants */
 static c@typ@ nc_1@c@ = {1., 0.};
 static c@typ@ nc_half@c@ = {0.5, 0.};
@@ -319,15 +330,33 @@ nc_asin@c@(c@typ@ *x, c@typ@ *r)
      * return nc_neg(nc_prodi(nc_log(nc_sum(nc_prod(nc_i,x),
      * nc_sqrt(nc_diff(nc_1,nc_prod(x,x)))))));
      */
-    c@typ@ a, *pa=&a;
-    nc_prod@c@(x, x, r);
-    nc_diff@c@(&nc_1@c@, r, r);
-    nc_sqrt@c@(r, r);
-    nc_prodi@c@(x, pa);
-    nc_sum@c@(pa, r, r);
-    nc_log@c@(r, r);
-    nc_prodi@c@(r, r);
-    nc_neg@c@(r, r);
+    if (fabs(x->real) > 1e-3 || fabs(x->imag) > 1e-3) {
+        c@typ@ a, *pa=&a;
+        nc_prod@c@(x, x, r);
+        nc_diff@c@(&nc_1@c@, r, r);
+        nc_sqrt@c@(r, r);
+        nc_prodi@c@(x, pa);
+        nc_sum@c@(pa, r, r);
+        nc_log@c@(r, r);
+        nc_prodi@c@(r, r);
+        nc_neg@c@(r, r);
+    }
+    else {
+        /*
+         * Small arguments: series expansion, to avoid loss of precision
+         * asin(x) = x [1 + (1/6) x^2 [1 + (9/20) x^2 [1 - ...]]]
+         *
+         * |x| < 1e-3 => |rel. error| < 1e-20
+         */
+        c@typ@ x2;
+        nc_prod@c@(x, x, &x2);
+        
+        *r = nc_1@c@;
+        SERIES_HORNER_TERM@c@(r, &x2, 25./42);
+        SERIES_HORNER_TERM@c@(r, &x2, 9./20);
+        SERIES_HORNER_TERM@c@(r, &x2, 1./6);
+        nc_prod@c@(r, x, r);
+    }
     return;
 }
 
@@ -338,11 +367,29 @@ nc_asinh@c@(c@typ@ *x, c@typ@ *r)
     /*
      * return nc_log(nc_sum(nc_sqrt(nc_sum(nc_1,nc_prod(x,x))),x));
      */
-    nc_prod@c@(x, x, r);
-    nc_sum@c@(&nc_1@c@, r, r);
-    nc_sqrt@c@(r, r);
-    nc_sum@c@(r, x, r);
-    nc_log@c@(r, r);
+    if (fabs(x->real) > 1e-3 || fabs(x->imag) > 1e-3) {
+        nc_prod@c@(x, x, r);
+        nc_sum@c@(&nc_1@c@, r, r);
+        nc_sqrt@c@(r, r);
+        nc_sum@c@(r, x, r);
+        nc_log@c@(r, r);
+    }
+    else {
+        /*
+         * Small arguments: series expansion, to avoid loss of precision
+         * asinh(x) = x [1 - (1/6) x^2 [1 - (9/20) x^2 [1 - ...]]]
+         *
+         * |x| < 1e-3 => |rel. error| < 1e-20
+         */
+        c@typ@ x2;
+        nc_prod@c@(x, x, &x2);
+        
+        *r = nc_1@c@;
+        SERIES_HORNER_TERM@c@(r, &x2, -25./42);
+        SERIES_HORNER_TERM@c@(r, &x2, -9./20);
+        SERIES_HORNER_TERM@c@(r, &x2, -1./6);
+        nc_prod@c@(r, x, r);
+    }
     return;
 }
 
@@ -352,12 +399,30 @@ nc_atan@c@(c@typ@ *x, c@typ@ *r)
     /*
      * return nc_prod(nc_i2,nc_log(nc_quot(nc_sum(nc_i,x),nc_diff(nc_i,x))));
      */
-    c@typ@ a, *pa=&a;
-    nc_diff@c@(&nc_i@c@, x, pa);
-    nc_sum@c@(&nc_i@c@, x, r);
-    nc_quot@c@(r, pa, r);
-    nc_log@c@(r,r);
-    nc_prod@c@(&nc_i2@c@, r, r);
+    if (fabs(x->real) > 1e-3 || fabs(x->imag) > 1e-3) {
+        c@typ@ a, *pa=&a;
+        nc_diff@c@(&nc_i@c@, x, pa);
+        nc_sum@c@(&nc_i@c@, x, r);
+        nc_quot@c@(r, pa, r);
+        nc_log@c@(r,r);
+        nc_prod@c@(&nc_i2@c@, r, r);
+    }
+    else {
+        /*
+         * Small arguments: series expansion, to avoid loss of precision
+         * atan(x) = x [1 - (1/3) x^2 [1 - (3/5) x^2 [1 - ...]]]
+         *
+         * |x| < 1e-3 => |rel. error| < 1e-20
+         */
+        c@typ@ x2;
+        nc_prod@c@(x, x, &x2);
+        
+        *r = nc_1@c@;
+        SERIES_HORNER_TERM@c@(r, &x2, -5./7);
+        SERIES_HORNER_TERM@c@(r, &x2, -3./5);
+        SERIES_HORNER_TERM@c@(r, &x2, -1./3);
+        nc_prod@c@(r, x, r);
+    }
     return;
 }
 
@@ -367,12 +432,30 @@ nc_atanh@c@(c@typ@ *x, c@typ@ *r)
     /*
      * return nc_prod(nc_half,nc_log(nc_quot(nc_sum(nc_1,x),nc_diff(nc_1,x))));
      */
-    c@typ@ a, *pa=&a;
-    nc_diff@c@(&nc_1@c@, x, r);
-    nc_sum@c@(&nc_1@c@, x, pa);
-    nc_quot@c@(pa, r, r);
-    nc_log@c@(r, r);
-    nc_prod@c@(&nc_half@c@, r, r);
+    if (fabs(x->real) > 1e-3 || fabs(x->imag) > 1e-3) {
+        c@typ@ a, *pa=&a;
+        nc_diff@c@(&nc_1@c@, x, r);
+        nc_sum@c@(&nc_1@c@, x, pa);
+        nc_quot@c@(pa, r, r);
+        nc_log@c@(r, r);
+        nc_prod@c@(&nc_half@c@, r, r);
+    }
+    else {
+        /*
+         * Small arguments: series expansion, to avoid loss of precision
+         * atan(x) = x [1 + (1/3) x^2 [1 + (3/5) x^2 [1 - ...]]]
+         *
+         * |x| < 1e-3 => |rel. error| < 1e-20
+         */
+        c@typ@ x2;
+        nc_prod@c@(x, x, &x2);
+        
+        *r = nc_1@c@;
+        SERIES_HORNER_TERM@c@(r, &x2, 5./7);
+        SERIES_HORNER_TERM@c@(r, &x2, 3./5);
+        SERIES_HORNER_TERM@c@(r, &x2, 1./3);
+        nc_prod@c@(r, x, r);
+    }
     return;
 }
 
@@ -463,5 +546,7 @@ nc_tanh@c@(c@typ@ *x, c@typ@ *r)
     return;
 }
 
+#undef SERIES_HORNER_TERM@c@
+
 /**end repeat**/
 
diff --git a/numpy/core/tests/test_umath.py b/numpy/core/tests/test_umath.py
index 1b4ba97ca..8a9be6374 100644
--- a/numpy/core/tests/test_umath.py
+++ b/numpy/core/tests/test_umath.py
@@ -422,6 +422,81 @@ class TestComplexFunctions(object):
 
                 assert abs(a - b) < atol, "%s %s: %s; cmath: %s"%(fname,p,a,b)
 
+    def check_loss_of_precision(self, dtype):
+        """Check loss of precision in complex arc* functions"""
+
+        # Check against known-good functions
+
+        info = np.finfo(dtype)
+        real_dtype = dtype(0.).real.dtype
+
+        eps = info.eps
+
+        def check(x, rtol):
+            d = np.absolute(np.arcsinh(x)/np.arcsinh(x+0j).real - 1)
+            assert np.all(d < rtol), (x[np.argmax(d)], d.max())
+
+            d = np.absolute(np.arcsinh(x)/np.arcsin(1j*x).imag - 1)
+            assert np.all(d < rtol), (x[np.argmax(d)], d.max())
+
+            d = np.absolute(np.arctanh(x)/np.arctanh(x+0j).real - 1)
+            assert np.all(d < rtol), (x[np.argmax(d)], d.max())
+
+            d = np.absolute(np.arctanh(x)/np.arctan(1j*x).imag - 1)
+            assert np.all(d < rtol), (x[np.argmax(d)], d.max())
+        
+        # The switchover was chosen as 1e-3; hence there can be up to
+        # ~eps/1e-3 of relative cancellation error before it
+
+        x_series = np.logspace(np.log10(info.tiny/eps).real, -3, 200,
+                               endpoint=False)
+        x_basic = np.logspace(dtype(-3.).real, -1e-8, 10)
+
+        check(x_series, 2*eps)
+        check(x_basic, 2*eps/1e-3)
+
+        # Check a few points
+
+        z = np.array([1e-5*(1+1j)], dtype=dtype)
+        p = 9.999999999333333333e-6 + 1.000000000066666666e-5j
+        d = np.absolute(1-np.arctanh(z)/p)
+        assert np.all(d < 1e-15)
+
+        p = 1.0000000000333333333e-5 + 9.999999999666666667e-6j
+        d = np.absolute(1-np.arcsinh(z)/p)
+        assert np.all(d < 1e-15)
+
+        p = 9.999999999333333333e-6j + 1.000000000066666666e-5
+        d = np.absolute(1-np.arctan(z)/p)
+        assert np.all(d < 1e-15)
+
+        p = 1.0000000000333333333e-5j + 9.999999999666666667e-6
+        d = np.absolute(1-np.arcsin(z)/p)
+        assert np.all(d < 1e-15)
+
+        # Check continuity across switchover points
+
+        def check(func, z0, d=1):
+            z0 = np.asarray(z0, dtype=dtype)
+            zp = z0 + abs(z0) * d * eps * 2
+            zm = z0 - abs(z0) * d * eps * 2
+            assert np.all(zp != zm), (zp, zm)
+
+            # NB: the cancellation error at the switchover is at least eps
+            good = (abs(func(zp) - func(zm)) < 2*eps)
+            assert np.all(good), (func, z0[~good])
+
+        for func in (np.arcsinh,np.arcsinh,np.arcsin,np.arctanh,np.arctan):
+            pts = [rp+1j*ip for rp in (-1e-3,0,1e-3) for ip in(-1e-3,0,1e-3)
+                   if rp != 0 or ip != 0]
+            check(func, pts, 1)
+            check(func, pts, 1j)
+            check(func, pts, 1+1j)
+
+    def test_loss_of_precision(self):
+        for dtype in [np.complex64, np.complex_, np.longcomplex]:
+            yield self.check_loss_of_precision, dtype
+
 class TestAttributes(TestCase):
     def test_attributes(self):
         add = ncu.add