[Jython-checkins] jython (merge default -> default): Merge cmath changes to trunk
jeff.allen
jython-checkins at python.org
Tue Jan 20 23:07:17 CET 2015
https://hg.python.org/jython/rev/9c81bac53d7a
changeset: 7553:9c81bac53d7a
parent: 7542:88209128cde8
parent: 7552:e89bf3a7c8dd
user: Jeff Allen <ja.py at farowl.co.uk>
date: Tue Jan 20 22:06:52 2015 +0000
summary:
Merge cmath changes to trunk
files:
Lib/test/cmath_testcases.txt | 30 +-
Lib/test/test_cmath.py | 144 +-
Lib/test/test_float_jy.py | 7 +
Misc/make_cmath_testcases.py | 14 +-
NEWS | 3 +-
src/org/python/core/PyFloat.java | 5 +-
src/org/python/modules/cmath.java | 1079 ++++++++++++++--
src/org/python/modules/math.java | 8 +-
8 files changed, 1042 insertions(+), 248 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/Lib/test/test_cmath.py b/Lib/test/test_cmath.py
--- a/Lib/test/test_cmath.py
+++ b/Lib/test/test_cmath.py
@@ -1,4 +1,4 @@
-from test.test_support import run_unittest, is_jython
+from test.test_support import run_unittest, verbose
from test.test_math import parse_testfile, test_file
import unittest
import cmath, math
@@ -52,7 +52,6 @@
'cos', 'cosh', 'exp', 'log', 'log10', 'sin', 'sinh',
'sqrt', 'tan', 'tanh']]
# test first and second arguments independently for 2-argument log
-
test_functions.append(lambda x : cmath.log(x, 1729. + 0j))
test_functions.append(lambda x : cmath.log(14.-27j, x))
@@ -244,24 +243,30 @@
unit_interval = test_values + [-x for x in test_values] + \
[0., 1., -1.]
+ # test_values for acosh, atanh
+ ge_one = [1.] + [1./x for x in test_values]
+ unit_open = test_values + [-x for x in test_values] + [0.]
+
# test_values for log, log10, sqrt
- positive = test_values + [1.] + [1./x for x in test_values]
+ positive = test_values + ge_one
nonnegative = [0.] + positive
# test_values for functions defined on the whole real line
real_line = [0.] + positive + [-x for x in positive]
test_functions = {
- # FIXME uncomment tests for Jython
- #'acos' : unit_interval,
- #'asin' : unit_interval,
- #'atan' : real_line,
- #'cos' : real_line,
- #'cosh' : real_line,
+ 'acos' : unit_interval,
+ 'asin' : unit_interval,
+ 'atan' : real_line,
+ 'acosh' : ge_one, # Jython added
+ 'asinh' : real_line, # Jython added
+ 'atanh' : unit_open, # Jython added
+ 'cos' : real_line,
+ 'cosh' : real_line,
'exp' : real_line,
'log' : positive,
'log10' : positive,
- #'sin' : real_line,
+ 'sin' : real_line,
'sinh' : real_line,
'sqrt' : nonnegative,
'tan' : real_line,
@@ -273,7 +278,7 @@
for v in values:
z = complex_fn(v)
self.rAssertAlmostEqual(float_fn(v), z.real)
- self.rAssertAlmostEqual(0., z.imag)
+ self.assertEqual(0., z.imag)
# test two-argument version of log with various bases
for base in [0.5, 2., 10.]:
@@ -282,10 +287,9 @@
self.rAssertAlmostEqual(math.log(v, base), z.real)
self.assertEqual(0., z.imag)
- @unittest.skipIf(is_jython, "FIXME: not working in Jython")
def test_specific_values(self):
if not float.__getformat__("double").startswith("IEEE"):
- return
+ self.skipTest('needs IEEE double')
def rect_complex(z):
"""Wrapped version of rect that accepts a complex number instead of
@@ -297,78 +301,88 @@
two floats."""
return complex(*polar(z))
+ raised_fmt = '\n' \
+ '{}: {}(complex({!r}, {!r}))\n' \
+ 'Raised: {!r}\n' \
+ 'Expected: complex({!r}, {!r})'
+ not_raised_fmt = '\n' \
+ '{} not raised in test {}: {}(complex({!r}, {!r}))\n' \
+ 'Received: complex({!r}, {!r})'
+ value_fmt = '\n' \
+ '{}: {}(complex({!r}, {!r}))\n' \
+ 'Expected: complex({!r}, {!r})\n' \
+ 'Received: complex({!r}, {!r})\n' \
+ 'Received value insufficiently close to expected value.'
+ failures = 0
+
for id, fn, ar, ai, er, ei, flags in parse_testfile(test_file):
arg = complex(ar, ai)
- expected = complex(er, ei)
if fn == 'rect':
function = rect_complex
elif fn == 'polar':
function = polar_complex
else:
function = getattr(cmath, fn)
- if 'divide-by-zero' in flags or 'invalid' in flags:
- try:
+
+ try:
+ # Catch and count failing test cases locally
+ if 'divide-by-zero' in flags or 'invalid' in flags:
try:
actual = function(arg)
except ValueError:
- continue
+ pass
else:
- self.fail('ValueError not raised in test '
- '{}: {}(complex({!r}, {!r}))'.format(id, fn, ar, ai))
- except AssertionError, ex:
- print "Got", function, ex
- except BaseException, ex:
- print "Got", function, ex
+ failures += 1
+ self.fail(not_raised_fmt.format('ValueError',
+ id, fn, ar, ai, actual.real, actual.imag))
- try:
- if 'overflow' in flags:
+ elif 'overflow' in flags:
try:
actual = function(arg)
except OverflowError:
- continue
- except BaseException, ex:
- print "\nGot", function, ex
+ pass
else:
- self.fail('OverflowError not raised in test '
- '{}: {}(complex({!r}, {!r}))'.format(id, fn, ar, ai))
- except AssertionError, ex:
- print "\nGot", function, ex
+ failures += 1
+ self.fail(not_raised_fmt.format('OverflowError',
+ id, fn, ar, ai, actual.real, actual.imag))
- try:
- actual = function(arg)
- except BaseException, ex:
- print "\nGot", function, ex
+ else :
+ actual = function(arg)
- if 'ignore-real-sign' in flags:
- actual = complex(abs(actual.real), actual.imag)
- expected = complex(abs(expected.real), expected.imag)
- if 'ignore-imag-sign' in flags:
- actual = complex(actual.real, abs(actual.imag))
- expected = complex(expected.real, abs(expected.imag))
+ # Make sign of expected conform to actual, where ignored.
+ exr, exi = er, ei
+ if 'ignore-real-sign' in flags:
+ exr = math.copysign(er, actual.real)
+ if 'ignore-imag-sign' in flags:
+ exi = math.copysign(ei, actual.imag)
- # for the real part of the log function, we allow an
- # absolute error of up to 2e-15.
- if fn in ('log', 'log10'):
- real_abs_err = 2e-15
- else:
- real_abs_err = 5e-323
+ # for the real part of the log function, we allow an
+ # absolute error of up to 2e-15.
+ if fn in ('log', 'log10'):
+ real_abs_err = 2e-15
+ else:
+ real_abs_err = 5e-323
- error_message = (
- '{}: {}(complex({!r}, {!r}))\n'
- 'Expected: complex({!r}, {!r})\n'
- 'Received: complex({!r}, {!r})\n'
- 'Received value insufficiently close to expected value.'
- ).format(id, fn, ar, ai,
- expected.real, expected.imag,
- actual.real, actual.imag)
- try:
- self.rAssertAlmostEqual(expected.real, actual.real,
- abs_err=real_abs_err,
- msg=error_message)
- self.rAssertAlmostEqual(expected.imag, actual.imag,
- msg=error_message)
- except AssertionError, ex:
- print "\nGot", ex, error_message
+ error_message = value_fmt.format(id, fn, ar, ai, er, ei,
+ actual.real, actual.imag)
+ self.rAssertAlmostEqual(exr, actual.real,
+ abs_err=real_abs_err,
+ msg=error_message)
+ self.rAssertAlmostEqual(exi, actual.imag,
+ msg=error_message)
+
+ except AssertionError as ex:
+ failures += 1
+ if verbose :
+ print(ex)
+ except BaseException as ex:
+ failures += 1
+ if verbose :
+ print(raised_fmt.format(id, fn, ar, ai, ex, er, ei))
+
+ if failures :
+ self.fail('{} discrepancies'.format(failures))
+
def assertCISEqual(self, a, b):
eps = 1E-7
@@ -446,6 +460,8 @@
self.assertTrue(math.isnan(abs(complex(2.3, NAN))))
self.assertEqual(abs(complex(INF, NAN)), INF)
self.assertTrue(math.isnan(abs(complex(NAN, NAN))))
+
+ # result overflows
if float.__getformat__("double").startswith("IEEE"):
self.assertRaises(OverflowError, abs, complex(1.4e308, 1.4e308))
diff --git a/Lib/test/test_float_jy.py b/Lib/test/test_float_jy.py
--- a/Lib/test/test_float_jy.py
+++ b/Lib/test/test_float_jy.py
@@ -85,6 +85,13 @@
self.assert_(type(float('inf')), float)
self.assertRaises(OverflowError, long, float('Infinity'))
+ def test_minus_zero(self):
+ # Some operations confused by -0.0
+ mz = float('-0.0')
+ self.assertEquals(mz, 0.)
+ self.assertEquals(repr(mz)[0], '-')
+ self.assertEquals(repr(abs(mz))[0], '0')
+
def test_float_none(self):
self.assertRaises(TypeError, float, None)
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/NEWS b/NEWS
--- a/NEWS
+++ b/NEWS
@@ -39,7 +39,8 @@
- [ 2252 ] Args in sys.argv are now unicode if characters > 127
- [ 2092 ] Upgrade to JLine2 and add tab completion support
- [ 2236 ] Interactive parser does not accept try ... except E as e: syntax
- - [ 2237 ] More robust math/cmath support (ongoing)
+ - [ 2237 ] Fully conformant math and cmath support
+ - [ 2244 ] More robust testing of math and cmath modules
New features
- Full support of Python buffer protocol, along with Java ByteBuffer support
diff --git a/src/org/python/core/PyFloat.java b/src/org/python/core/PyFloat.java
--- a/src/org/python/core/PyFloat.java
+++ b/src/org/python/core/PyFloat.java
@@ -839,10 +839,7 @@
@ExposedMethod(doc = BuiltinDocs.float___abs___doc)
final PyObject float___abs__() {
- if (getValue() < 0) {
- return float___neg__();
- }
- return float___float__();
+ return new PyFloat(Math.abs(getValue()));
}
@Override
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
@@ -2,6 +2,7 @@
import org.python.core.Py;
import org.python.core.PyComplex;
+import org.python.core.PyException;
import org.python.core.PyFloat;
import org.python.core.PyInstance;
import org.python.core.PyObject;
@@ -12,48 +13,21 @@
public static final PyFloat pi = new PyFloat(Math.PI);
public static final PyFloat e = new PyFloat(Math.E);
- private static final PyComplex one = new PyComplex(1.0, 0.0);
- private static final PyComplex half = new PyComplex(0.5, 0.0);
- private static final PyComplex i = new PyComplex(0.0, 1.0);
- private static final PyComplex half_i = new PyComplex(0.0, 0.5);
-
/** 2<sup>-½</sup> (Ref: Abramowitz & Stegun [1972], p2). */
private static final double ROOT_HALF = 0.70710678118654752440;
+ /** ln({@link Double#MAX_VALUE}) or a little less */
+ private static final double NEARLY_LN_DBL_MAX = 709.4361393;
+ /**
+ * For x larger than this, <i>e<sup>-x</sup></i> is negligible compared with
+ * <i>e<sup>x</sup></i>, or equivalently 1 is negligible compared with <i>e<sup>2x</sup></i>, in
+ * IEEE-754 floating point. Beyond this, sinh <i>x</i> and cosh <i>x</i> are adequately
+ * approximated by 0.5<i>e<sup>x</sup></i>. The smallest theoretical value is 27 ln(2).
+ */
+ private static final double ATLEAST_27LN2 = 18.72;
+ private static final double HALF_E2 = 0.5 * Math.E * Math.E;
- private static PyComplex c_prodi(PyComplex x) {
- return (PyComplex)x.__mul__(i);
- }
-
- private static boolean isNaN(PyComplex x) {
- return Double.isNaN(x.real) || Double.isNaN(x.imag);
- }
-
- private static double abs(PyComplex x) {
- boolean isNaN = isNaN(x);
- boolean isInfinite = !isNaN && (Double.isInfinite(x.real) || Double.isInfinite(x.imag));
- if (isNaN) {
- return Double.NaN;
- }
- if (isInfinite) {
- return Double.POSITIVE_INFINITY;
- }
- double real_abs = Math.abs(x.real);
- double imag_abs = Math.abs(x.imag);
-
- if (real_abs < imag_abs) {
- if (x.imag == 0.0) {
- return real_abs;
- }
- double q = x.real / x.imag;
- return imag_abs * Math.sqrt(1 + q * q);
- } else {
- if (x.real == 0.0) {
- return imag_abs;
- }
- double q = x.imag / x.real;
- return real_abs * Math.sqrt(1 + q * q);
- }
- }
+ /** log<sub>10</sub>e (Ref: Abramowitz & Stegun [1972], p3). */
+ private static final double LOG10E = 0.43429448190325182765;
private static PyComplex complexFromPyObject(PyObject obj) {
// If op is already of type PyComplex_Type, return its value
@@ -88,74 +62,555 @@
return new PyComplex(obj.asDouble(), 0);
}
- public static PyObject acos(PyObject in) {
- PyComplex x = complexFromPyObject(in);
- return c_prodi(log(x.__add__(i.__mul__(sqrt(one.__sub__(x.__mul__(x))))))).__neg__();
+ /**
+ * Return the arc cosine of w. There are two branch cuts. One extends right from 1 along the
+ * real axis to ∞, continuous from below. The other extends left from -1 along the real
+ * axis to -∞, continuous from above.
+ *
+ * @param w
+ * @return cos<sup>-1</sup><i>w</i>
+ */
+ public static PyComplex acos(PyObject w) {
+ return _acos(complexFromPyObject(w));
}
- public static PyComplex acosh(PyObject in) {
- PyComplex x = complexFromPyObject(in);
- PyComplex a = sqrt(x.__sub__(one));
- PyComplex b = sqrt(x.__add__(one));
- PyComplex c = sqrt(half);
- PyComplex r = log(c.__mul__(b.__add__(a)));
- return ((PyComplex)r.__add__(r));
+ /**
+ * Helper to compute cos<sup>-1</sup><i>w</i>. The method used is as in CPython:
+ * <p>
+ * <i>a = (1-w)<sup>½</sup> = √2</i> sin <i>z/2</i> <br>
+ * <i>b = (1+w)<sup>½</sup> = √2</i> cos <i>z/2</i>
+ * <p>
+ * Then, with <i>z = x+iy</i>, <i>a = a<sub>1</sub>+ia<sub>2</sub></i>, and <i>b =
+ * b<sub>1</sub>+ib<sub>2</sub></i>,
+ * <p>
+ * a<sub>1</sub> / b<sub>1</sub> = tan <i>x/2</i> <br>
+ * a<sub>2</sub>b<sub>1</sub> - a<sub>1</sub>b<sub>2</sub> = sinh <i>y</i>
+ * <p>
+ * and we use {@link Math#atan2(double, double)} and {@link math#asinh(double)} to obtain
+ * <i>x</i> and <i>y</i>.
+ * <p>
+ * For <i>w</i> sufficiently large that <i>w<sup>2</sup></i>≫1, cos<sup>-1</sup><i>w</i>
+ * ≈ -i ln(<i>2w</i>).
+ *
+ * @param w
+ * @return cos<sup>-1</sup><i>w</i>
+ */
+ private static PyComplex _acos(PyComplex w) {
+
+ // Let z = x + iy and w = u + iv.
+ double x, y, u = w.real, v = w.imag;
+
+ if (Math.abs(u) > 0x1p27 || Math.abs(v) > 0x1p27) {
+ /*
+ * w is large: approximate 2cos(z) by exp(i(x+iy)) or exp(-i(x+iy)), whichever
+ * dominates. Hence, z = x+iy = i ln(2(u+iv)) or -i ln(2(u+iv))
+ */
+ x = Math.atan2(Math.abs(v), u);
+ y = Math.copySign(logHypot(u, v) + math.LN2, -v);
+
+ } else if (Double.isNaN(v)) {
+ // Special cases
+ x = (u == 0.) ? Math.PI / 2. : v;
+ y = v;
+
+ } else {
+ // Normal case, without risk of overflow.
+ PyComplex a = sqrt(new PyComplex(1. - u, -v)); // a = sqrt(1-w) = sqrt(2) sin(z/2)
+ PyComplex b = sqrt(new PyComplex(1 + u, v)); // b = sqrt(1+w) = sqrt(2) cos(z/2)
+ // Arguments here are sin(x/2)cosh(y/2), cos(x/2)cosh(y/2) giving tan(x/2)
+ x = 2. * Math.atan2(a.real, b.real);
+ // 2 (cos(x/2)**2+sin(x/2)**2) sinh(y/2)cosh(y/2) = sinh y
+ y = math.asinh(a.imag * b.real - a.real * b.imag);
+ }
+
+ // If that generated a nan, and there wasn't one in the argument, raise a domain error.
+ return exceptNaN(new PyComplex(x, y), w);
}
- public static PyComplex asin(PyObject in) {
- PyComplex x = complexFromPyObject(in);
- PyComplex squared = (PyComplex)x.__mul__(x);
- PyComplex sq1_minus_xsq = sqrt(one.__sub__(squared));
- return (PyComplex)c_prodi(log(sq1_minus_xsq.__add__(c_prodi(x)))).__neg__();
+ /**
+ * Return the hyperbolic arc cosine of w. There is one branch cut, extending left from 1 along
+ * the real axis to -∞, continuous from above.
+ *
+ * @param w
+ * @return cosh<sup>-1</sup><i>w</i>
+ */
+ public static PyComplex acosh(PyObject w) {
+ return _acosh(complexFromPyObject(w));
}
- public static PyComplex asinh(PyObject in) {
- PyComplex x = complexFromPyObject(in);
- PyComplex a = sqrt(x.__add__(i));
- PyComplex b = sqrt(x.__sub__(i));
- PyComplex z = sqrt(half);
- PyComplex r = log(z.__mul__(a.__add__(b)));
- return ((PyComplex)r.__add__(r));
+ /**
+ * Helper to compute z = cosh<sup>-1</sup><i>w</i>. The method used is as in CPython:
+ * <p>
+ * <i>a = (w-1)<sup>½</sup> = √2</i> sinh <i>z/2</i> <br>
+ * <i>b = (w+1)<sup>½</sup> = √2</i> cosh <i>z/2</i>
+ * <p>
+ * Then, with <i>z = x+iy</i>, <i>a = a<sub>1</sub>+ia<sub>2</sub></i>, and <i>b =
+ * b<sub>1</sub>+ib<sub>2</sub></i>,
+ * <p>
+ * a<sub>1</sub>b<sub>1</sub> + a<sub>2</sub>b<sub>2</sub> = sinh <i>x</i> <br>
+ * a<sub>2</sub> / b<sub>1</sub> = tan <i>y/2</i>
+ * <p>
+ * and we use {@link math#asinh(double)} and {@link Math#atan2(double, double)} to obtain
+ * <i>x</i> and <i>y</i>.
+ * <p>
+ * For <i>w</i> sufficiently large that <i>w<sup>2</sup></i>≫1,
+ * cosh<sup>-1</sup><i>w</i> ≈ ln(<i>2w</i>). We do not use this method also to compute
+ * cos<sup>-1</sup><i>w</i>, because the branch cuts do not correspond.
+ *
+ * @param w
+ * @return cosh<sup>-1</sup><i>w</i>
+ */
+ private static PyComplex _acosh(PyComplex w) {
+
+ // Let z = x + iy and w = u + iv.
+ double x, y, u = w.real, v = w.imag;
+
+ if (Math.abs(u) > 0x1p27 || Math.abs(v) > 0x1p27) {
+ /*
+ * w is large: approximate 2cosh(z) by exp(x+iy) or exp(-x-iy), whichever dominates.
+ * Hence, z = x+iy = ln(2(u+iv)) or -ln(2(u+iv))
+ */
+ x = logHypot(u, v) + math.LN2;
+ y = Math.atan2(v, u);
+
+ } else if (v == 0. && !Double.isNaN(u)) {
+ /*
+ * We're on the real axis (and maybe the branch cut). u = cosh x cos y. In all cases,
+ * the sign of y follows v.
+ */
+ if (u >= 1.) {
+ // As real library, cos y = 1, u = cosh x.
+ x = math.acosh(u);
+ y = v;
+ } else if (u < -1.) {
+ // Left part of cut: cos y = -1, u = -cosh x
+ x = math.acosh(-u);
+ y = Math.copySign(Math.PI, v);
+ } else {
+ // -1 <= u <= 1: cosh x = 1, u = cos y.
+ x = 0.;
+ y = Math.copySign(Math.acos(u), v);
+ }
+
+ } else {
+ // Normal case, without risk of overflow.
+ PyComplex a = sqrt(new PyComplex(u - 1., v)); // a = sqrt(w-1) = sqrt(2) sinh(z/2)
+ PyComplex b = sqrt(new PyComplex(u + 1., v)); // b = sqrt(w+1) = sqrt(2) cosh(z/2)
+ // 2 sinh(x/2)cosh(x/2) (cos(y/2)**2+sin(y/2)**2) = sinh x
+ x = math.asinh(a.real * b.real + a.imag * b.imag);
+ // Arguments here are cosh(x/2)sin(y/2) and cosh(x/2)cos(y/2) giving tan y/2
+ y = 2. * Math.atan2(a.imag, b.real);
+ }
+
+ // If that generated a nan, and there wasn't one in the argument, raise a domain error.
+ return exceptNaN(new PyComplex(x, y), w);
}
- public static PyComplex atan(PyObject in) {
- PyComplex x = complexFromPyObject(in);
- return (PyComplex)half_i.__mul__(log(i.__add__(x).__div__(i.__sub__(x))));
+ /**
+ * Return the arc sine of w. There are two branch cuts. One extends right from 1 along the real
+ * axis to ∞, continuous from below. The other extends left from -1 along the real axis to
+ * -∞, continuous from above.
+ *
+ * @param w
+ * @return sin<sup>-1</sup><i>w</i>
+ */
+ public static PyComplex asin(PyObject w) {
+ return asinOrAsinh(complexFromPyObject(w), false);
}
- public static PyComplex atanh(PyObject in) {
- PyComplex x = complexFromPyObject(in);
- return (PyComplex)half.__mul__(log(one.__add__(x).__div__(one.__sub__(x))));
+ /**
+ * Return the hyperbolic arc sine of w. There are two branch cuts. One extends from 1j along the
+ * imaginary axis to ∞j, continuous from the right. The other extends from -1j along the
+ * imaginary axis to -∞j, continuous from the left.
+ *
+ * @param w
+ * @return sinh<sup>-1</sup><i>w</i>
+ */
+ public static PyComplex asinh(PyObject w) {
+ return asinOrAsinh(complexFromPyObject(w), true);
}
- 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));
+ /**
+ * Helper to compute either sin<sup>-1</sup><i>w</i> or sinh<sup>-1</sup><i>w</i>. The method
+ * used is as in CPython:
+ * <p>
+ * <i>a = (1-iw)<sup>½</sup> = √2 </i>sin(<i>π/4-iz/2</i>) <br>
+ * <i>b = (1+iw)<sup>½</sup> = √2 </i>cos(<i>π/4-iz/2</i>)
+ * <p>
+ * Then, with <i>w = u+iv</i>, <i>z = x+iy</i>, <i>a = a<sub>1</sub>+ia<sub>2</sub></i>, and
+ * <i>b = b<sub>1</sub>+ib<sub>2</sub></i>,
+ * <p>
+ * a<sub>1</sub>b<sub>2</sub> - a<sub>2</sub>b<sub>2</sub> = sinh <i>x</i> <br>
+ * v / (a<sub>2</sub>b<sub>1</sub> - a<sub>1</sub>b<sub>2</sub>) = tan <i>y</i>
+ * <p>
+ * and we use {@link math#asinh(double)} and {@link Math#atan2(double, double)} to obtain
+ * <i>x</i> and <i>y</i>.
+ * <p>
+ * For <i>w</i> sufficiently large that <i>w<sup>2</sup></i>≫1,
+ * sinh<sup>-1</sup><i>w</i> ≈ ln(<i>2w</i>). When computing sin<sup>-1</sup><i>w</i>, we
+ * evaluate <i>-i</i> sinh<sup>-1</sup><i>iw</i> instead.
+ *
+ * @param w
+ * @param h <code>true</code> to compute sinh<sup>-1</sup><i>w</i>, <code>false</code> to
+ * compute sin<sup>-1</sup><i>w</i>.
+ * @return sinh<sup>-1</sup><i>w</i> or sin<sup>-1</sup><i>w</i>
+ */
+ private static PyComplex asinOrAsinh(PyComplex w, boolean h) {
+ double u, v, x, y;
+ PyComplex z;
+
+ if (h) {
+ // We compute z = asinh(w). Let z = x + iy and w = u + iv.
+ u = w.real;
+ v = w.imag;
+ // Then the function body computes x + iy = asinh(w).
+ } else {
+ // We compute w = asin(z). Unusually, let w = u - iv, so u + iv = iw.
+ v = w.real;
+ u = -w.imag;
+ // Then as before, the function body computes asinh(u+iv) = asinh(iw) = i asin(w),
+ // but we finally return z = y - ix = -i asinh(iw) = asin(w).
+ }
+
+ if (Double.isNaN(u)) {
+ // Special case for nan in real part. Default clause deals naturally with v=nan.
+ if (v == 0.) {
+ x = u;
+ y = v;
+ } else if (Double.isInfinite(v)) {
+ x = Double.POSITIVE_INFINITY;
+ y = u;
+ } else { // Any other value of v -> nan+nanj
+ x = y = u;
+ }
+
+ } else if (Math.abs(u) > 0x1p27 || Math.abs(v) > 0x1p27) {
+ /*
+ * w is large: approximate 2sinh(z) by exp(x+iy) or -exp(-x-iy), whichever dominates.
+ * Hence, z = x+iy = ln(2(u+iv)) or -ln(-2(u+iv))
+ */
+ x = logHypot(u, v) + math.LN2;
+ if (Math.copySign(1., u) > 0.) {
+ y = Math.atan2(v, u);
+ } else {
+ // Adjust for sign, choosing the angle so that -pi/2 < y < pi/2
+ x = -x;
+ y = Math.atan2(v, -u);
+ }
+
+ } else {
+ // Normal case, without risk of overflow.
+ PyComplex a = sqrt(new PyComplex(1. + v, -u)); // a = sqrt(1-iw)
+ PyComplex b = sqrt(new PyComplex(1. - v, u)); // b = sqrt(1+iw)
+ // Combine the parts so as that terms in y cancel, leaving us with sinh x:
+ x = math.asinh(a.real * b.imag - a.imag * b.real);
+ // The arguments are v = cosh x sin y, and cosh x cos y
+ y = Math.atan2(v, a.real * b.real - a.imag * b.imag);
+ }
+
+ // Compose the result w according to whether we're computing asin(w) or asinh(w).
+ if (h) {
+ z = new PyComplex(x, y); // z = x + iy = asinh(u+iv).
+ } else {
+ z = new PyComplex(y, -x); // z = y - ix = -i asinh(v-iu) = asin(w)
+ }
+
+ // If that generated a nan, and there wasn't one in the argument, raise a domain error.
+ return exceptNaN(z, w);
}
- 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 arc tangent of w. There are two branch cuts. One extends from 1j along the
+ * imaginary axis to ∞j, continuous from the right. The other extends from -1j along the
+ * imaginary axis to -∞j, continuous from the left.
+ *
+ * @param w
+ * @return tan<sup>-1</sup><i>w</i>
+ */
+ public static PyComplex atan(PyObject w) {
+ return atanOrAtanh(complexFromPyObject(w), false);
}
- public static PyComplex exp(PyObject in) {
- PyComplex x = complexFromPyObject(in);
- double l = Math.exp(x.real);
- return new PyComplex(l * Math.cos(x.imag), l * Math.sin(x.imag));
+ /**
+ * Return the hyperbolic arc tangent of w. There are two branch cuts. One extends from 1 along
+ * the real axis to ∞, continuous from below. The other extends from -1 along the real
+ * axis to -∞, continuous from above.
+ *
+ * @param w
+ * @return tanh<sup>-1</sup><i>w</i>
+ */
+ public static PyComplex atanh(PyObject w) {
+ return atanOrAtanh(complexFromPyObject(w), true);
}
- 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);
+ /**
+ * Helper to compute either tan<sup>-1</sup><i>w</i> or tanh<sup>-1</sup><i>w</i>. The method
+ * used is close to that used in CPython. For <i>z</i> = tanh<sup>-1</sup><i>w</i>:
+ * <p>
+ * <i>z = </i>½ln(<i>1 + 2w/(1-w)</i>)
+ * <p>
+ * Then, letting <i>z = x+iy</i>, and <i>w = u+iv</i>,
+ * <p>
+ * <i>x = </i>¼ln(<i>1 + 4u/((1-u)<sup>2</sup>+v<sup>2</sup>)</i>) <i> =
+ * -</i>¼ln(<i>1 - 4u/((1+u)<sup>2</sup>+v<sup>2</sup>)</i>)<br>
+ * <i>y = </i>½tan<sup>-1</sup>(<i>2v / ((1+u)(1-u)-v<sup>2</sup>)</i>)<br>
+ * <p>
+ * We use {@link math#log1p(double)} and {@link Math#atan2(double, double)} to obtain <i>x</i>
+ * and <i>y</i>. The second expression for <code>x</code> is used when <i>u<0</i>. For
+ * <i>w</i> sufficiently large that <i>w<sup>2</sup></i>≫1, tanh<sup>-1</sup><i>w</i>
+ * ≈ 1/w ± <i>iπ/2</i>). For small <i>w</i>, tanh<sup>-1</sup><i>w</i> ≈
+ * <i>w</i>. When computing tan<sup>-1</sup><i>w</i>, we evaluate <i>-i</i>
+ * tanh<sup>-1</sup><i>iw</i> instead.
+ *
+ * @param w
+ * @param h <code>true</code> to compute tanh<sup>-1</sup><i>w</i>, <code>false</code> to
+ * compute tan<sup>-1</sup><i>w</i>.
+ * @return tanh<sup>-1</sup><i>w</i> or tan<sup>-1</sup><i>w</i>
+ */
+ private static PyComplex atanOrAtanh(PyComplex w, boolean h) {
+ double u, v, x, y;
+ PyComplex z;
+
+ if (h) {
+ // We compute z = atanh(w). Let z = x + iy and w = u + iv.
+ u = w.real;
+ v = w.imag;
+ // Then the function body computes x + iy = atanh(w).
+ } else {
+ // We compute w = atan(z). Unusually, let w = u - iv, so u + iv = iw.
+ v = w.real;
+ u = -w.imag;
+ // Then as before, the function body computes atanh(u+iv) = atanh(iw) = i atan(w),
+ // but we finally return z = y - ix = -i atanh(iw) = atan(w).
+ }
+
+ double absu = Math.abs(u), absv = Math.abs(v);
+
+ if (absu >= 0x1p511 || absv >= 0x1p511) {
+ // w is large: approximate atanh(w) by 1/w + i pi/2. 1/w = conjg(w)/|w|**2.
+ if (Double.isInfinite(absu) || Double.isInfinite(absv)) {
+ x = Math.copySign(0., u);
} else {
- return PyComplex.NaN;
+ // w is also too big to square, carry a 2**-N scaling factor.
+ int N = 520;
+ double uu = Math.scalb(u, -N), vv = Math.scalb(v, -N);
+ double mod2w = uu * uu + vv * vv;
+ x = Math.scalb(uu / mod2w, -N);
+ }
+ // We don't need the imaginary part of 1/z. Just pi/2 with the sign of v. (If not nan.)
+ if (Double.isNaN(v)) {
+ y = v;
+ } else {
+ y = Math.copySign(Math.PI / 2., v);
+ }
+
+ } else if (absu < 0x1p-53) {
+ // u is small enough that u**2 may be neglected relative to 1.
+ if (absv > 0x1p-27) {
+ // v is not small, but is not near overflow either.
+ double v2 = v * v;
+ double d = 1. + v2;
+ x = Math.copySign(Math.log1p(4. * absu / d), u) * 0.25;
+ y = Math.atan2(2. * v, 1. - v2) * 0.5;
+ } else {
+ // v is also small enough that v**2 may be neglected (or is nan). So z = w.
+ x = u;
+ y = v;
+ }
+
+ } else if (absu == 1. && absv < 0x1p-27) {
+ // w is close to +1 or -1: needs a different expression, good as v->0
+ x = Math.copySign(Math.log(absv) - math.LN2, u) * 0.5;
+ if (v == 0.) {
+ y = Double.NaN;
+ } else {
+ y = Math.copySign(Math.atan2(2., absv), v) * 0.5;
+ }
+
+ } else {
+ /*
+ * Normal case, without risk of overflow. The basic expression is z =
+ * 0.5*ln((1+w)/(1-w)), which for positive u we rearrange as 0.5*ln(1+2w/(1-w)) and for
+ * negative u as -0.5*ln(1-2w/(1+w)). By use of absu, we reduce the difference between
+ * the expressions fo u>=0 and u<0 to a sign transfer.
+ */
+ double lmu = (1. - absu), lpu = (1. + absu), v2 = v * v;
+ double d = lmu * lmu + v2;
+ x = Math.copySign(Math.log1p(4. * absu / d), u) * 0.25;
+ y = Math.atan2(2. * v, lmu * lpu - v2) * 0.5;
+ }
+
+ // Compose the result w according to whether we're computing atan(w) or atanh(w).
+ if (h) {
+ z = new PyComplex(x, y); // z = x + iy = atanh(u+iv).
+ } else {
+ z = new PyComplex(y, -x); // z = y - ix = -i atanh(v-iu) = atan(w)
+ }
+
+ // If that generated a nan, and there wasn't one in the argument, raise a domain error.
+ return exceptNaN(z, w);
+ }
+
+ /**
+ * Return the cosine of z.
+ *
+ * @param z
+ * @return cos <i>z</i>
+ */
+ public static PyComplex cos(PyObject z) {
+ return cosOrCosh(complexFromPyObject(z), false);
+ }
+
+ /**
+ * 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;
}
}
- return new PyComplex(Math.log(abs(x)), Math.atan2(x.imag, x.real));
+
+ // 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);
+ }
+
+ /**
+ * Return the exponential value e<sup>z</sup>.
+ *
+ * @param z
+ * @return e<sup>z</sup>
+ */
+ public static PyComplex exp(PyObject z) {
+ PyComplex zz = complexFromPyObject(z);
+ double x = zz.real, y = zz.imag, r, u, v;
+ /*
+ * This has a lot of corner-cases, and some of them make little sense sense, but it matches
+ * CPython and passes the regression tests.
+ */
+ if (y == 0.) {
+ // Real value: use a real solution. (This may raise a range error.)
+ u = math.exp(x);
+ // v follows sign of y.
+ v = 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);
+
+ if (x == Double.NEGATIVE_INFINITY) {
+ // w = (0,0) but "signed" by the direction cosines (even in they are nan).
+ u = Math.copySign(0., cosy);
+ v = Math.copySign(0., siny);
+
+ } else if (x == Double.POSITIVE_INFINITY) {
+ if (!Double.isNaN(cosy)) {
+ // w = (inf,inf), but "signed" by the direction cosines.
+ u = Math.copySign(x, cosy);
+ v = Math.copySign(x, siny);
+ } else {
+ // Provisionally w = (inf,nan), which will raise domain error if y!=nan.
+ u = x;
+ v = Double.NaN;
+ }
+
+ } else if (x > NEARLY_LN_DBL_MAX) {
+ // r = e**x would overflow but maybe not r*cos(y) and r*sin(y).
+ r = Math.exp(x - 1); // = r / e
+ u = r * cosy * Math.E;
+ v = r * siny * Math.E;
+ 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.
+ // Compute r = exp(x), and return w = u + iv = r (cos(y) + i*sin(y))
+ r = Math.exp(x);
+ u = r * cosy;
+ v = r * siny;
+ }
+ }
+ // If that generated a nan, and there wasn't one in the argument, raise domain error.
+ return exceptNaN(new PyComplex(u, v), zz);
}
public static double phase(PyObject in) {
@@ -170,25 +625,42 @@
return new PyTuple(new PyFloat(r), new PyFloat(phi));
}
+ /**
+ * Return the complex number x with polar coordinates r and phi. Equivalent to
+ * <code>r * (math.cos(phi) + math.sin(phi)*1j)</code>.
+ *
+ * @param r radius
+ * @param phi angle
+ * @return
+ */
public static PyComplex rect(double r, double phi) {
- // Handle various edge cases
+ double x, y;
+
if (Double.isInfinite(r) && (Double.isInfinite(phi) || Double.isNaN(phi))) {
- return new PyComplex(Double.POSITIVE_INFINITY, Double.NaN);
+ x = Double.POSITIVE_INFINITY;
+ y = Double.NaN;
+
+ } else if (phi == 0.0) {
+ // cos(phi)=1, sin(phi)=phi: finesse oddball r in computing y, but not x.
+ x = r;
+ if (Double.isNaN(r)) {
+ y = phi;
+ } else if (Double.isInfinite(r)) {
+ y = phi * Math.copySign(1., r);
+ } else {
+ y = phi * r;
+ }
+
+ } else if (r == 0.0 && (Double.isInfinite(phi) || Double.isNaN(phi))) {
+ // Ignore any problems (inf, nan) with phi
+ x = y = 0.;
+
+ } else {
+ // Text-book case, using the trig functions.
+ x = r * Math.cos(phi);
+ y = r * Math.sin(phi);
}
- if (phi == 0.0) { // NB this test will succeed if phi is 0.0 or -0.0
- if (Double.isNaN(r)) {
- return new PyComplex(Double.NaN, 0.0);
- } else if (r == Double.POSITIVE_INFINITY) {
- return new PyComplex(r, phi);
- } else if (r == Double.NEGATIVE_INFINITY) {
- return new PyComplex(r, -phi);
- }
- }
- if (r == 0.0 && (Double.isInfinite(phi) || Double.isNaN(phi))) {
- return new PyComplex(0.0, 0.0);
- }
-
- return new PyComplex(r * Math.cos(phi), r * Math.sin(phi));
+ return exceptNaN(new PyComplex(x, y), r, phi);
}
/**
@@ -211,48 +683,231 @@
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
+ * ½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 ½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));
+ /**
+ * Return the sine of z.
+ *
+ * @param z
+ * @return sin <i>z</i>
+ */
+ public static PyComplex sin(PyObject z) {
+ return sinOrSinh(complexFromPyObject(z), false);
}
- 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);
+ /**
+ * Return the hyperbolic sine of z.
+ *
+ * @param z
+ * @return sinh <i>z</i>
+ */
+ public static PyComplex sinh(PyObject z) {
+ return sinOrSinh(complexFromPyObject(z), true);
+ }
+
+ /**
+ * Helper to compute either sin <i>z</i> or sinh <i>z</i>.
+ *
+ * @param z
+ * @param h <code>true</code> to compute sinh <i>z</i>, <code>false</code> to compute sin
+ * <i>z</i>.
+ * @return
+ */
+ private static PyComplex sinOrSinh(PyComplex z, boolean h) {
+ double x, y, u, v;
+ PyComplex w;
+
+ if (h) {
+ // 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:
+ // u = sinh(x) cos(y),
+ // v = cosh(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 as before, the function body computes sinh(x+iy) = sinh(iz) = i sin(z),
+ // but we finally return w = v - iu = sin(z).
+ }
+
+ if (y == 0.) {
+ // Real argument for sinh (or imaginary for sin): use real library.
+ u = math.sinh(x); // This will raise a range error on overflow.
+ // v follows the sign of y (which could be -0.0).
+ v = y;
+
+ } else if (x == 0.) {
+ // Imaginary argument for sinh (or real for sin): imaginary result at this point.
+ v = Math.sin(y);
+ // u follows sign of x (which could be -0.0).
+ u = x;
+
+ } 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 = x * cosy;
+ v = absx * siny;
+ } else {
+ // Provisionally w = (inf,nan), which will raise domain error if y!=nan.
+ u = x;
+ 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:
+ v = r * siny * HALF_E2;
+ // r approximates 2sinh(|x|)/e**2: put back the proper sign of x in passing.
+ u = Math.copySign(r, x) * cosy * HALF_E2;
+ if (Double.isInfinite(u) || Double.isInfinite(v)) {
+ // A finite x gave rise to an infinite u or v.
+ throw math.mathRangeError();
+ }
+
} else {
- return PyComplex.NaN;
+ // Normal case, without risk of overflow.
+ u = Math.sinh(x) * cosy;
+ v = Math.cosh(x) * siny;
}
}
- 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) {
- PyComplex x = complexFromPyObject(in);
- return new PyComplex(Math.sin(x.real) * math.cosh(x.imag), Math.cos(x.real)
- * math.sinh(x.imag));
- }
+ // Compose the result w according to whether we're computing sin(z) or sinh(z).
+ if (h) {
+ w = new PyComplex(u, v); // w = u + iv = sinh(x+iy).
+ } else {
+ w = new PyComplex(v, -u); // w = v - iu = sin(y-ix) = sin(z)
+ }
- public static PyComplex sinh(PyObject in) {
- PyComplex x = complexFromPyObject(in);
- return new PyComplex(Math.cos(x.imag) * math.sinh(x.real), Math.sin(x.imag)
- * math.cosh(x.real));
+ // If that generated a nan, and there wasn't one in the argument, raise a domain error.
+ return exceptNaN(w, z);
}
/**
@@ -365,35 +1020,143 @@
}
}
- public static PyComplex tan(PyObject in) {
- PyComplex x = complexFromPyObject(in);
-
- double sr = Math.sin(x.real);
- double cr = Math.cos(x.real);
- double shi = math.sinh(x.imag);
- double chi = math.cosh(x.imag);
- double rs = sr * chi;
- double is = cr * shi;
- double rc = cr * chi;
- double ic = -sr * shi;
- double d = rc * rc + ic * ic;
-
- return new PyComplex(((rs * rc) + (is * ic)) / d, ((is * rc) - (rs * ic)) / d);
+ /**
+ * Return the tangent of z.
+ *
+ * @param z
+ * @return tan <i>z</i>
+ */
+ public static PyComplex tan(PyObject z) {
+ return tanOrTanh(complexFromPyObject(z), false);
}
- public static PyComplex tanh(PyObject in) {
- PyComplex x = complexFromPyObject(in);
+ /**
+ * Return the hyperbolic tangent of z.
+ *
+ * @param z
+ * @return tanh <i>z</i>
+ */
+ public static PyComplex tanh(PyObject z) {
+ return tanOrTanh(complexFromPyObject(z), true);
+ }
- double si = Math.sin(x.imag);
- double ci = Math.cos(x.imag);
- double shr = math.sinh(x.real);
- double chr = math.cosh(x.real);
- double rs = ci * shr;
- double is = si * chr;
- double rc = ci * chr;
- double ic = si * shr;
- double d = rc * rc + ic * ic;
+ /**
+ * Helper to compute either tan <i>z</i> or tanh <i>z</i>. The expression used is:
+ * <p>
+ * tanh(<i>x+iy</i>) = (sinh <i>x</i> cosh <i>x + i</i> sin <i>y</i> cos <i>y</i>) /
+ * (sinh<sup>2</sup><i>x +</i> cos<sup>2</sup><i>y</i>)
+ * <p>
+ * A simplification is made for x sufficiently large that <i>e<sup>2|x|</sup></i>≫1 that
+ * deals satisfactorily with large or infinite <i>x</i>. When computing tan, we evaluate
+ * <i>i</i> tan <i>iz</i> instead, and the approximation applies to
+ * <i>e<sup>2|y|</sup></i>≫1.
+ *
+ * @param z
+ * @param h <code>true</code> to compute tanh <i>z</i>, <code>false</code> to compute tan
+ * <i>z</i>.
+ * @return tan or tanh <i>z</i>
+ */
+ private static PyComplex tanOrTanh(PyComplex z, boolean h) {
+ double x, y, u, v, s;
+ PyComplex w;
- return new PyComplex(((rs * rc) + (is * ic)) / d, ((is * rc) - (rs * ic)) / d);
+ if (h) {
+ // We compute w = tanh(z). Let w = u + iv and z = x + iy.
+ x = z.real;
+ y = z.imag;
+ // Then the function body computes tanh(x+iy), according to:
+ // s = sinh**2 x + cos**2 y
+ // u = sinh x cosh x / s,
+ // v = sin y cos y / s,
+ // And we return w = u + iv.
+ } else {
+ // We compute w = tan(z). Unusually, let z = y - ix, so x + iy = iz.
+ y = z.real;
+ x = -z.imag;
+ // Then the function body computes tanh(x+iy) = tanh(iz) = i tan(z) as before,
+ // but we finally return w = v - iu = tan(z).
+ }
+
+ if (y == 0.) {
+ // Real argument for tanh (or imaginary for tan).
+ u = Math.tanh(x);
+ // v is zero but follows the sign of y (which could be -0.0).
+ v = y;
+
+ } else if (x == 0. && !Double.isNaN(y)) {
+ // Imaginary argument for tanh (or real for tan): imaginary result at this point.
+ v = Math.tan(y); // May raise domain error
+ // u is zero but follows sign of x (which could be -0.0).
+ u = x;
+
+ } 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 > ATLEAST_27LN2) {
+ // e**2x is much greater than 1: exponential approximation to sinh and cosh.
+ s = 0.25 * Math.exp(2 * absx);
+ u = Math.copySign(1., x);
+ if (s == Double.POSITIVE_INFINITY) {
+ // Either x is inf or 2x is large enough to overflow exp(). v=0, but signed:
+ v = Math.copySign(0., siny * cosy);
+ } else {
+ v = siny * cosy / s;
+ }
+
+ } else {
+ // Normal case: possible overflow in s near (x,y) = (0,pi/2) is harmless.
+ double sinhx = Math.sinh(x), coshx = Math.cosh(x);
+ s = sinhx * sinhx + cosy * cosy;
+ u = sinhx * coshx / s;
+ v = siny * cosy / s;
+ }
+ }
+
+ // Compose the result w according to whether we're computing tan(z) or tanh(z).
+ if (h) {
+ w = new PyComplex(u, v); // w = u + iv = tanh(x+iy).
+ } else {
+ w = new PyComplex(v, -u); // w = v - iu = tan(y-ix) = tan(z)
+ }
+
+ // If that generated a nan, and there wasn't one in the argument, raise a domain error.
+ return exceptNaN(w, z);
}
+
+ /**
+ * Turn a <code>NaN</code> result into a thrown <code>ValueError</code>, a math domain error, if
+ * the original argument was not itself <code>NaN</code>. A <code>PyComplex</code> is a
+ * <code>NaN</code> if either component is a <code>NaN</code>.
+ *
+ * @param result to return (if we return)
+ * @param arg to include in check
+ * @return result if <code>arg</code> was <code>NaN</code> or <code>result</code> was not
+ * <code>NaN</code>
+ * @throws PyException (ValueError) if <code>result</code> was <code>NaN</code> and
+ * <code>arg</code> was not <code>NaN</code>
+ */
+ private static PyComplex exceptNaN(PyComplex result, PyComplex arg) throws PyException {
+ if ((Double.isNaN(result.real) || Double.isNaN(result.imag))
+ && !(Double.isNaN(arg.real) || Double.isNaN(arg.imag))) {
+ throw math.mathDomainError();
+ } else {
+ return result;
+ }
+ }
+
+ /**
+ * Raise <code>ValueError</code> if <code>result</code> is a <code>NaN</code>, but neither
+ * <code>a</code> nor <code>b</code> is <code>NaN</code>. Same as
+ * {@link #exceptNaN(PyComplex, PyComplex)}.
+ */
+ private static PyComplex exceptNaN(PyComplex result, double a, double b) throws PyException {
+ if ((Double.isNaN(result.real) || Double.isNaN(result.imag))
+ && !(Double.isNaN(a) || Double.isNaN(b))) {
+ throw math.mathDomainError();
+ } else {
+ return result;
+ }
+ }
+
}
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;
@@ -515,7 +515,7 @@
*
* @return ValueError("math domain error")
*/
- private static PyException mathDomainError() {
+ static PyException mathDomainError() {
return Py.ValueError("math domain error");
}
@@ -524,7 +524,7 @@
*
* @return OverflowError("math range error")
*/
- private static PyException mathRangeError() {
+ static PyException mathRangeError() {
return Py.OverflowError("math range error");
}
@@ -575,7 +575,7 @@
*/
private static double exceptInf(double result, double arg) {
if (Double.isInfinite(result) && !Double.isInfinite(arg)) {
- throw Py.OverflowError("math range error");
+ throw mathRangeError();
} else {
return result;
}
--
Repository URL: https://hg.python.org/jython
More information about the Jython-checkins
mailing list