[Numpy-svn] r6520 - in trunk/numpy/core: src tests
numpy-svn at scipy.org
numpy-svn at scipy.org
Sun Mar 1 12:07:29 EST 2009
Author: ptvirtan
Date: 2009-03-01 11:07:06 -0600 (Sun, 01 Mar 2009)
New Revision: 6520
Modified:
trunk/numpy/core/src/umath_funcs.inc.src
trunk/numpy/core/tests/test_umath.py
Log:
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.
Modified: trunk/numpy/core/src/umath_funcs.inc.src
===================================================================
--- trunk/numpy/core/src/umath_funcs.inc.src 2009-02-27 22:27:14 UTC (rev 6519)
+++ trunk/numpy/core/src/umath_funcs.inc.src 2009-03-01 17:07:06 UTC (rev 6520)
@@ -96,6 +96,17 @@
#c=f,,l#
*/
+/* Perform the operation result := 1 + coef * x * result,
+ * with real coefficient `coef`.
+ */
+#define SERIES_HORNER_TERM at c@(result, x, coef) \
+ do { \
+ nc_prod at c@((result), (x), (result)); \
+ (result)->real *= (coef); \
+ (result)->imag *= (coef); \
+ nc_sum at c@((result), &nc_1 at c@, (result)); \
+ } while(0)
+
/* constants */
static c at typ@ nc_1 at c@ = {1., 0.};
static c at typ@ nc_half at c@ = {0.5, 0.};
@@ -319,15 +330,33 @@
* 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 at typ@ a, *pa=&a;
- nc_prod at c@(x, x, r);
- nc_diff at c@(&nc_1 at c@, r, r);
- nc_sqrt at c@(r, r);
- nc_prodi at c@(x, pa);
- nc_sum at c@(pa, r, r);
- nc_log at c@(r, r);
- nc_prodi at c@(r, r);
- nc_neg at c@(r, r);
+ if (fabs(x->real) > 1e-3 || fabs(x->imag) > 1e-3) {
+ c at typ@ a, *pa=&a;
+ nc_prod at c@(x, x, r);
+ nc_diff at c@(&nc_1 at c@, r, r);
+ nc_sqrt at c@(r, r);
+ nc_prodi at c@(x, pa);
+ nc_sum at c@(pa, r, r);
+ nc_log at c@(r, r);
+ nc_prodi at c@(r, r);
+ nc_neg at 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 at typ@ x2;
+ nc_prod at c@(x, x, &x2);
+
+ *r = nc_1 at c@;
+ SERIES_HORNER_TERM at c@(r, &x2, 25./42);
+ SERIES_HORNER_TERM at c@(r, &x2, 9./20);
+ SERIES_HORNER_TERM at c@(r, &x2, 1./6);
+ nc_prod at c@(r, x, r);
+ }
return;
}
@@ -338,11 +367,29 @@
/*
* return nc_log(nc_sum(nc_sqrt(nc_sum(nc_1,nc_prod(x,x))),x));
*/
- nc_prod at c@(x, x, r);
- nc_sum at c@(&nc_1 at c@, r, r);
- nc_sqrt at c@(r, r);
- nc_sum at c@(r, x, r);
- nc_log at c@(r, r);
+ if (fabs(x->real) > 1e-3 || fabs(x->imag) > 1e-3) {
+ nc_prod at c@(x, x, r);
+ nc_sum at c@(&nc_1 at c@, r, r);
+ nc_sqrt at c@(r, r);
+ nc_sum at c@(r, x, r);
+ nc_log at 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 at typ@ x2;
+ nc_prod at c@(x, x, &x2);
+
+ *r = nc_1 at c@;
+ SERIES_HORNER_TERM at c@(r, &x2, -25./42);
+ SERIES_HORNER_TERM at c@(r, &x2, -9./20);
+ SERIES_HORNER_TERM at c@(r, &x2, -1./6);
+ nc_prod at c@(r, x, r);
+ }
return;
}
@@ -352,12 +399,30 @@
/*
* return nc_prod(nc_i2,nc_log(nc_quot(nc_sum(nc_i,x),nc_diff(nc_i,x))));
*/
- c at typ@ a, *pa=&a;
- nc_diff at c@(&nc_i at c@, x, pa);
- nc_sum at c@(&nc_i at c@, x, r);
- nc_quot at c@(r, pa, r);
- nc_log at c@(r,r);
- nc_prod at c@(&nc_i2 at c@, r, r);
+ if (fabs(x->real) > 1e-3 || fabs(x->imag) > 1e-3) {
+ c at typ@ a, *pa=&a;
+ nc_diff at c@(&nc_i at c@, x, pa);
+ nc_sum at c@(&nc_i at c@, x, r);
+ nc_quot at c@(r, pa, r);
+ nc_log at c@(r,r);
+ nc_prod at c@(&nc_i2 at 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 at typ@ x2;
+ nc_prod at c@(x, x, &x2);
+
+ *r = nc_1 at c@;
+ SERIES_HORNER_TERM at c@(r, &x2, -5./7);
+ SERIES_HORNER_TERM at c@(r, &x2, -3./5);
+ SERIES_HORNER_TERM at c@(r, &x2, -1./3);
+ nc_prod at c@(r, x, r);
+ }
return;
}
@@ -367,12 +432,30 @@
/*
* return nc_prod(nc_half,nc_log(nc_quot(nc_sum(nc_1,x),nc_diff(nc_1,x))));
*/
- c at typ@ a, *pa=&a;
- nc_diff at c@(&nc_1 at c@, x, r);
- nc_sum at c@(&nc_1 at c@, x, pa);
- nc_quot at c@(pa, r, r);
- nc_log at c@(r, r);
- nc_prod at c@(&nc_half at c@, r, r);
+ if (fabs(x->real) > 1e-3 || fabs(x->imag) > 1e-3) {
+ c at typ@ a, *pa=&a;
+ nc_diff at c@(&nc_1 at c@, x, r);
+ nc_sum at c@(&nc_1 at c@, x, pa);
+ nc_quot at c@(pa, r, r);
+ nc_log at c@(r, r);
+ nc_prod at c@(&nc_half at 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 at typ@ x2;
+ nc_prod at c@(x, x, &x2);
+
+ *r = nc_1 at c@;
+ SERIES_HORNER_TERM at c@(r, &x2, 5./7);
+ SERIES_HORNER_TERM at c@(r, &x2, 3./5);
+ SERIES_HORNER_TERM at c@(r, &x2, 1./3);
+ nc_prod at c@(r, x, r);
+ }
return;
}
@@ -463,5 +546,7 @@
return;
}
+#undef SERIES_HORNER_TERM at c@
+
/**end repeat**/
Modified: trunk/numpy/core/tests/test_umath.py
===================================================================
--- trunk/numpy/core/tests/test_umath.py 2009-02-27 22:27:14 UTC (rev 6519)
+++ trunk/numpy/core/tests/test_umath.py 2009-03-01 17:07:06 UTC (rev 6520)
@@ -422,6 +422,81 @@
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
More information about the Numpy-svn
mailing list