summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--numpy/core/src/umath_funcs.inc.src137
-rw-r--r--numpy/core/tests/test_umath.py75
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
View Lorry Controller Status