[Jython-checkins] jython: Re-work cmath.cos and cmath.cosh for accuracy and corner cases.
jeff.allen
jython-checkins at python.org
Tue Jan 20 23:07:15 CET 2015
https://hg.python.org/jython/rev/1bcc508b734a
changeset: 7547:1bcc508b734a
user: Jeff Allen <ja.py at farowl.co.uk>
date: Sat Jan 10 23:25:02 2015 +0000
summary:
Re-work cmath.cos and cmath.cosh for accuracy and corner cases.
This includes changes to cmath_testcases.txt, and its generator,
so that reference results that are ostensibly real have -0.0 as their
imaginary part where cmath should. (Only cos and cosh affected so far.)
files:
Lib/test/cmath_testcases.txt | 30 ++--
Misc/make_cmath_testcases.py | 14 +-
src/org/python/modules/cmath.java | 119 +++++++++++++++--
3 files changed, 131 insertions(+), 32 deletions(-)
diff --git a/Lib/test/cmath_testcases.txt b/Lib/test/cmath_testcases.txt
--- a/Lib/test/cmath_testcases.txt
+++ b/Lib/test/cmath_testcases.txt
@@ -1750,15 +1750,15 @@
cosh0053 cosh 0.0003 0.0 -> 1.0000000450000003375 0.0
cosh0054 cosh 0.2 0.0 -> 1.0200667556190758485 0.0
cosh0055 cosh 1.0 0.0 -> 1.5430806348152437785 0.0
-cosh0056 cosh -1e-18 0.0 -> 1.0 0.0
-cosh0057 cosh -0.0003 0.0 -> 1.0000000450000003375 0.0
-cosh0058 cosh -1.0 0.0 -> 1.5430806348152437785 0.0
+cosh0056 cosh -1e-18 0.0 -> 1.0 -0.0
+cosh0057 cosh -0.0003 0.0 -> 1.0000000450000003375 -0.0
+cosh0058 cosh -1.0 0.0 -> 1.5430806348152437785 -0.0
cosh0059 cosh 1.3169578969248168 0.0 -> 2.0000000000000001504 0.0
-cosh0060 cosh -1.3169578969248168 0.0 -> 2.0000000000000001504 0.0
+cosh0060 cosh -1.3169578969248168 0.0 -> 2.0000000000000001504 -0.0
cosh0061 cosh 17.328679513998633 0.0 -> 16777216.000000021938 0.0
cosh0062 cosh 18.714973875118524 0.0 -> 67108864.000000043662 0.0
cosh0063 cosh 709.7827 0.0 -> 8.9883497833190073272e+307 0.0
-cosh0064 cosh -709.7827 0.0 -> 8.9883497833190073272e+307 0.0
+cosh0064 cosh -709.7827 0.0 -> 8.9883497833190073272e+307 -0.0
-- special values
cosh1000 cosh 0.0 0.0 -> 1.0 0.0
@@ -2071,20 +2071,20 @@
cos0023 cos 0.45124351152540226 1.6992693993812158 -> 2.543477948972237 -1.1528193694875477
-- Additional real values (Jython)
-cos0050 cos 1e-150 0.0 -> 1.0 0.0
-cos0051 cos 1e-18 0.0 -> 1.0 0.0
-cos0052 cos 1e-09 0.0 -> 0.9999999999999999995 0.0
-cos0053 cos 0.0003 0.0 -> 0.9999999550000003375 0.0
-cos0054 cos 0.2 0.0 -> 0.98006657784124162892 0.0
-cos0055 cos 1.0 0.0 -> 0.5403023058681397174 0.0
+cos0050 cos 1e-150 0.0 -> 1.0 -0.0
+cos0051 cos 1e-18 0.0 -> 1.0 -0.0
+cos0052 cos 1e-09 0.0 -> 0.9999999999999999995 -0.0
+cos0053 cos 0.0003 0.0 -> 0.9999999550000003375 -0.0
+cos0054 cos 0.2 0.0 -> 0.98006657784124162892 -0.0
+cos0055 cos 1.0 0.0 -> 0.5403023058681397174 -0.0
cos0056 cos -1e-18 0.0 -> 1.0 0.0
cos0057 cos -0.0003 0.0 -> 0.9999999550000003375 0.0
cos0058 cos -1.0 0.0 -> 0.5403023058681397174 0.0
-cos0059 cos 1.0471975511965976 0.0 -> 0.50000000000000009945 0.0
-cos0060 cos 2.5707963267948966 0.0 -> -0.84147098480789647357 0.0
+cos0059 cos 1.0471975511965976 0.0 -> 0.50000000000000009945 -0.0
+cos0060 cos 2.5707963267948966 0.0 -> -0.84147098480789647357 -0.0
cos0061 cos -2.5707963267948966 0.0 -> -0.84147098480789647357 0.0
-cos0062 cos 18 0.0 -> 0.66031670824408014482 0.0
-cos0063 cos 18.0 0.0 -> 0.66031670824408014482 0.0
+cos0062 cos 18 0.0 -> 0.66031670824408014482 -0.0
+cos0063 cos 18.0 0.0 -> 0.66031670824408014482 -0.0
-- special values
cos1000 cos -0.0 0.0 -> 1.0 0.0
diff --git a/Misc/make_cmath_testcases.py b/Misc/make_cmath_testcases.py
--- a/Misc/make_cmath_testcases.py
+++ b/Misc/make_cmath_testcases.py
@@ -103,14 +103,24 @@
}
def generate_cases() :
- fmt = "{}{:04d} {} {!r} 0.0 -> {} 0.0"
+ fmt = "{}{:04d} {} {!r} 0.0 -> {} {!r}"
for fn in sorted(cases_to_generate.keys()):
print "-- Additional real values (Jython)"
count, xlist = cases_to_generate[fn]
for x in xlist:
+ # Compute the function (in the reference library)
func = getattr(mpmath, fn)
y = func(x)
- print fmt.format(fn, count, fn, x, mpmath.nstr(y, 20) )
+ # For the benefit of cmath tests, get the sign of imaginary zero right
+ zero = 0.0
+ if math.copysign(1., x) > 0.:
+ if fn=='cos' :
+ zero = -0.0
+ else :
+ if fn=='cosh' :
+ zero = -0.0
+ # Output one test case at sufficient precision
+ print fmt.format(fn, count, fn, x, mpmath.nstr(y, 20), zero )
count += 1
def test_main():
diff --git a/src/org/python/modules/cmath.java b/src/org/python/modules/cmath.java
--- a/src/org/python/modules/cmath.java
+++ b/src/org/python/modules/cmath.java
@@ -142,16 +142,105 @@
return (PyComplex)half.__mul__(log(one.__add__(x).__div__(one.__sub__(x))));
}
- public static PyComplex cos(PyObject in) {
- PyComplex x = complexFromPyObject(in);
- return new PyComplex(Math.cos(x.real) * math.cosh(x.imag), -Math.sin(x.real)
- * math.sinh(x.imag));
+ /**
+ * Return the cosine of z.
+ *
+ * @param z
+ * @return cos <i>z</i>
+ */
+ public static PyComplex cos(PyObject z) {
+ return cosOrCosh(complexFromPyObject(z), false);
}
- public static PyComplex cosh(PyObject in) {
- PyComplex x = complexFromPyObject(in);
- return new PyComplex(Math.cos(x.imag) * math.cosh(x.real), Math.sin(x.imag)
- * math.sinh(x.real));
+ /**
+ * Return the hyperbolic cosine of z.
+ *
+ * @param z
+ * @return cosh <i>z</i>
+ */
+ public static PyComplex cosh(PyObject z) {
+ return cosOrCosh(complexFromPyObject(z), true);
+ }
+
+ /**
+ * Helper to compute either cos <i>z</i> or cosh <i>z</i>.
+ *
+ * @param z
+ * @param h <code>true</code> to compute cosh <i>z</i>, <code>false</code> to compute cos
+ * <i>z</i>.
+ * @return
+ */
+ private static PyComplex cosOrCosh(PyComplex z, boolean h) {
+ double x, y, u, v;
+ PyComplex w;
+
+ if (h) {
+ // We compute w = cosh(z). Let w = u + iv and z = x + iy.
+ x = z.real;
+ y = z.imag;
+ // Then the function body computes cosh(x+iy), according to:
+ // u = cosh(x) cos(y),
+ // v = sinh(x) sin(y),
+ // And we return w = u + iv.
+ } else {
+ // We compute w = sin(z). Unusually, let z = y - ix, so x + iy = iz.
+ y = z.real;
+ x = -z.imag;
+ // Then the function body computes cosh(x+iy) = cosh(iz) = cos(z) as before.
+ }
+
+ if (y == 0.) {
+ // Real argument for cosh (or imaginary for cos): use real library.
+ u = math.cosh(x); // This will raise a range error on overflow.
+ // v is zero but follows the sign of x*y (in which y could be -0.0).
+ v = Math.copySign(1., x) * y;
+
+ } else if (x == 0.) {
+ // Imaginary argument for cosh (or real for cos): imaginary result at this point.
+ u = Math.cos(y);
+ // v is zero but follows the sign of x*y (in which x could be -0.0).
+ v = x * Math.copySign(1., y);
+
+ } else {
+
+ // The trig calls will not throw, although if y is infinite, they return nan.
+ double cosy = Math.cos(y), siny = Math.sin(y), absx = Math.abs(x);
+
+ if (absx == Double.POSITIVE_INFINITY) {
+ if (!Double.isNaN(cosy)) {
+ // w = (inf,inf), but "rotated" by the direction cosines.
+ u = absx * cosy;
+ v = x * siny;
+ } else {
+ // Provisionally w = (inf,nan), which will raise domain error if y!=nan.
+ u = absx;
+ v = Double.NaN;
+ }
+
+ } else if (absx > ATLEAST_27LN2) {
+ // Use 0.5*e**x approximation. This is also the region where we risk overflow.
+ double r = Math.exp(absx - 2.);
+ // r approximates 2cosh(x)/e**2: multiply in this order to avoid inf:
+ u = r * cosy * HALF_E2;
+ // r approximates 2sinh(|x|)/e**2: put back the proper sign of x in passing.
+ v = Math.copySign(r, x) * siny * HALF_E2;
+ if (Double.isInfinite(u) || Double.isInfinite(v)) {
+ // A finite x gave rise to an infinite u or v.
+ throw math.mathRangeError();
+ }
+
+ } else {
+ // Normal case, without risk of overflow.
+ u = Math.cosh(x) * cosy;
+ v = Math.sinh(x) * siny;
+ }
+ }
+
+ // Compose the result w = u + iv.
+ w = new PyComplex(u, v);
+
+ // If that generated a nan, and there wasn't one in the argument, raise a domain error.
+ return exceptNaN(w, z);
}
/**
@@ -391,6 +480,7 @@
/**
* Return the sine of z.
+ *
* @param z
* @return sin <i>z</i>
*/
@@ -400,6 +490,7 @@
/**
* Return the hyperbolic sine of z.
+ *
* @param z
* @return sinh <i>z</i>
*/
@@ -420,7 +511,7 @@
PyComplex w;
if (h) {
- // We compute w = sinh(z). Let w = u + iv and z = x + iy. We compute w = sinh(z) from:
+ // We compute w = sinh(z). Let w = u + iv and z = x + iy.
x = z.real;
y = z.imag;
// Then the function body computes sinh(x+iy), according to:
@@ -428,13 +519,11 @@
// v = cosh(x) sin(y),
// And we return w = u + iv.
} else {
- // We compute w = sin(z). Unusually, let z = y - ix.
+ // We compute w = sin(z). Unusually, let z = y - ix, so x + iy = iz.
y = z.real;
x = -z.imag;
- // Then the function body computes sinh(x+iy) = sinh(iz) = i sin(z), according to:
- // u = sinh(x) cos(y),
- // v = cosh(x) sin(y).
- // But we finally return w = v - iu = sin(z)
+ // Then the function body computes sinh(x+iy) = sinh(iz) = i sin(z) as before,
+ // but we finally return w = v - iu = sin(z).
}
if (y == 0.) {
@@ -458,7 +547,7 @@
if (!Double.isNaN(cosy)) {
// w = (inf,inf), but "rotated" by the direction cosines.
u = x * cosy;
- v = Math.abs(x) * siny;
+ v = absx * siny;
} else {
// Provisionally w = (inf,nan), which will raise domain error if y!=nan.
u = x;
--
Repository URL: https://hg.python.org/jython
More information about the Jython-checkins
mailing list