[Jython-checkins] jython: Re-work cmath.log, and cmath.log10 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/0c240eade776
changeset:   7545:0c240eade776
user:        Jeff Allen <ja.py at farowl.co.uk>
date:        Fri Jan 09 23:20:20 2015 +0000
summary:
  Re-work cmath.log, and cmath.log10 for accuracy and corner cases.

files:
  src/org/python/modules/cmath.java |  158 +++++++++++++----
  src/org/python/modules/math.java  |    2 +-
  2 files changed, 120 insertions(+), 40 deletions(-)


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
@@ -23,6 +23,9 @@
     /** ln({@link Double#MAX_VALUE}) or a little less */
     private static final double NEARLY_LN_DBL_MAX = 709.4361393;
 
+    /** log<sub>10</sub>e (Ref: Abramowitz & Stegun [1972], p3). */
+    private static final double LOG10E = 0.43429448190325182765;
+
     private static PyComplex c_prodi(PyComplex x) {
         return (PyComplex)x.__mul__(i);
     }
@@ -204,18 +207,6 @@
         return exceptNaN(new PyComplex(u, v), zz);
     }
 
-    public static PyComplex log(PyObject in) {
-        PyComplex x = complexFromPyObject(in);
-        if (isNaN(x)) {
-            if (Double.isInfinite(x.real) || Double.isInfinite(x.imag)) {
-                return new PyComplex(Double.POSITIVE_INFINITY, Double.NaN);
-            } else {
-                return PyComplex.NaN;
-            }
-        }
-        return new PyComplex(Math.log(abs(x)), Math.atan2(x.imag, x.real));
-    }
-
     public static double phase(PyObject in) {
         PyComplex x = complexFromPyObject(in);
         return Math.atan2(x.imag, x.real);
@@ -269,36 +260,125 @@
         return Double.isNaN(x.real) || Double.isNaN(x.imag);
     }
 
-    public static PyComplex log10(PyObject in) {
-        PyComplex x = complexFromPyObject(in);
-        if (isNaN(x)) {
-            if (Double.isInfinite(x.real) || Double.isInfinite(x.imag)) {
-                return new PyComplex(Double.POSITIVE_INFINITY, Double.NaN);
+    /**
+     * Returns the natural logarithm of <i>w</i>.
+     *
+     * @param w
+     * @return ln <i>w</i>
+     */
+    public static PyComplex log(PyObject w) {
+        PyComplex ww = complexFromPyObject(w);
+        double u = ww.real, v = ww.imag;
+        // The real part of the result is the log of the magnitude.
+        double lnr = logHypot(u, v);
+        // The imaginary part of the result is the arg. This may result in a nan.
+        double theta = Math.atan2(v, u);
+        PyComplex z = new PyComplex(lnr, theta);
+        return exceptNaN(z, ww);
+    }
+
+    /**
+     * Returns the common logarithm of <i>w</i> (base 10 logarithm).
+     *
+     * @param w
+     * @return log<sub>10</sub><i>w</i>
+     */
+    public static PyComplex log10(PyObject w) {
+        PyComplex ww = complexFromPyObject(w);
+        double u = ww.real, v = ww.imag;
+        // The expression is the same as for base e, scaled in magnitude.
+        double logr = logHypot(u, v) * LOG10E;
+        double theta = Math.atan2(v, u) * LOG10E;
+        PyComplex z = new PyComplex(logr, theta);
+        return exceptNaN(z, ww);
+    }
+
+    /**
+     * Returns the logarithm of <i>w</i> to the given base. If the base is not specified, returns
+     * the natural logarithm of <i>w</i>. There is one branch cut, from 0 along the negative real
+     * axis to -∞, continuous from above.
+     *
+     * @param w
+     * @param b
+     * @return log<sub>b</sub><i>w</i>
+     */
+    public static PyComplex log(PyObject w, PyObject b) {
+        PyComplex ww = complexFromPyObject(w), bb = complexFromPyObject(b), z;
+        double u = ww.real, v = ww.imag, br = bb.real, bi = bb.imag, x, y;
+        // Natural log of w is (x,y)
+        x = logHypot(u, v);
+        y = Math.atan2(v, u);
+
+        if (bi != 0. || br <= 0.) {
+            // Complex or negative real base requires complex log: general case.
+            PyComplex lnb = log(bb);
+            z = (PyComplex)(new PyComplex(x, y)).__div__(lnb);
+
+        } else {
+            // Real positive base: frequent case. (b = inf or nan ends up here too.)
+            double lnb = Math.log(br);
+            z = new PyComplex(x / lnb, y / lnb);
+        }
+
+        return exceptNaN(z, ww);
+    }
+
+    /**
+     * Helper function for the log of a complex number, dealing with the log magnitude, and without
+     * intermediate overflow or underflow. It returns ln <i>r</i>, where <i>r<sup>2</sup> =
+     * u<sup>2</sup>+v<sup>2</sup></i>. To do this it computes
+     * &frac12;ln(u<sup>2</sup>+v<sup>2</sup>). Special cases are handled as follows:
+     * <ul>
+     * <li>if u or v is NaN, it returns NaN</li>
+     * <li>if u or v is infinite, it returns positive infinity</li>
+     * <li>if u and v are both zero, it raises a ValueError</li>
+     * </ul>
+     * We have this function instead of <code>Math.log(Math.hypot(u,v))</code> because a valid
+     * result is still possible even when <code>hypot(u,v)</code> overflows, and because there's no
+     * point in taking a square root when a log is to follow.
+     *
+     * @param u
+     * @param v
+     * @return &frac12;ln(u<sup>2</sup>+v<sup>2</sup>)
+     */
+    private static double logHypot(double u, double v) {
+
+        if (Double.isInfinite(u) || Double.isInfinite(v)) {
+            return Double.POSITIVE_INFINITY;
+
+        } else {
+            // Cannot overflow, but if u=v=0 will return -inf.
+            int scale = 0, ue = Math.getExponent(u), ve = Math.getExponent(v);
+            double lnr;
+
+            if (ue < -511 && ve < -511) {
+                // Both u and v are too small to square, or zero. (Just one would be ok.)
+                scale = 600;
+            } else if (ue > 510 || ve > 510) {
+                // One of these is too big to square and double (or is nan or inf).
+                scale = -600;
+            }
+
+            if (scale == 0) {
+                // Normal case: there is no risk of overflow or log of zero.
+                lnr = 0.5 * Math.log(u * u + v * v);
             } else {
-                return PyComplex.NaN;
+                // We must work with scaled values, us = u * 2**n etc..
+                double us = Math.scalb(u, scale);
+                double vs = Math.scalb(v, scale);
+                // rs**2 = r**2 * 2**2n
+                double rs2 = us * us + vs * vs;
+                // So ln(r) = ln(u**2+v**2)/2 = ln(us**2+vs**2)/2 - n ln(2)
+                lnr = 0.5 * Math.log(rs2) - scale * math.LN2;
+            }
+
+            // (u,v) = 0 leads to ln(r) = -inf, but that's a domain error
+            if (lnr == Double.NEGATIVE_INFINITY) {
+                throw math.mathDomainError();
+            } else {
+                return lnr;
             }
         }
-        double l = abs(x);
-        return new PyComplex(math.log10(new PyFloat(l)), Math.atan2(x.imag, x.real)
-                / Math.log(10.0));
-    }
-
-    public static PyComplex log(PyObject in, PyObject base) {
-        return log(complexFromPyObject(in), complexFromPyObject(base));
-    }
-
-    public static PyComplex log(PyComplex x, PyComplex base) {
-        if (isNaN(x)) {
-            if (Double.isInfinite(x.real) || Double.isInfinite(x.imag)) {
-                return new PyComplex(Double.POSITIVE_INFINITY, Double.NaN);
-            } else {
-                return PyComplex.NaN;
-            }
-        }
-        double l = abs(x);
-        PyComplex log_base = log(base);
-        return (PyComplex)new PyComplex(math.log(new PyFloat(l)), Math.atan2(x.imag, x.real))
-                .__div__(log_base);
     }
 
     public static PyComplex sin(PyObject in) {
diff --git a/src/org/python/modules/math.java b/src/org/python/modules/math.java
--- a/src/org/python/modules/math.java
+++ b/src/org/python/modules/math.java
@@ -24,7 +24,7 @@
     private static final double MINUS_ONE = -1.0;
     private static final double TWO = 2.0;
     private static final double EIGHT = 8.0;
-    private static final double LN2 = 0.693147180559945309417232121458; // Ref OEIS A002162
+    static final double LN2 = 0.693147180559945309417232121458; // Ref OEIS A002162
 
     private static final double INF = Double.POSITIVE_INFINITY;
     private static final double NINF = Double.NEGATIVE_INFINITY;

-- 
Repository URL: https://hg.python.org/jython


More information about the Jython-checkins mailing list