arccosh for complex numbers, goofy choice of branch
OK, Which branch do we want to use. As it currently is in numpy and scipy.special arccosh(1.5) = 0.96242365011920694 arccosh(1.5+0j) = -0.96242365011920705 + 0.0j This is consistent with gsl, but inconsistent with Mathematica, NAG, Maple, and probably all sensible implementations which use the generally accepted principal value. I've left this inconsistency raising an error in the ufunc tests until we make a decision. It might be nice to know what FORTRAN and MatLab do with this. Chuck
Matlab is consistent, I'm afraid:
acosh(1.5) ans = 0.9624 acosh(1.5 + 0j) ans = 0.9624
L. On Mon, Mar 17, 2008 at 9:40 AM, Charles R Harris <charlesr.harris@gmail.com> wrote:
OK,
Which branch do we want to use. As it currently is in numpy and scipy.special
arccosh(1.5) = 0.96242365011920694 arccosh(1.5+0j) = -0.96242365011920705 + 0.0j
This is consistent with gsl, but inconsistent with Mathematica, NAG, Maple, and probably all sensible implementations which use the generally accepted principal value. I've left this inconsistency raising an error in the ufunc tests until we make a decision. It might be nice to know what FORTRAN and MatLab do with this.
Chuck
_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
-- Lorenzo Bolla lbolla@gmail.com http://lorenzobolla.emurse.com/
Mon, 17 Mar 2008 08:07:38 -0600, Charles R Harris wrote: [clip]
OK, that does it. I'm going to change it's behavior.
The problem with bad arccosh branch cuts is still present:
import numpy as np numpy.__version__ '1.2.0.dev5436.e45a7627a39d' np.arccosh(-1e-9 + 0.1j) (-0.099834078899207618-1.5707963277899337j) np.arccosh(1e-9 + 0.1j) (0.099834078899207576+1.5707963257998594j) np.arccosh(-1e-9 - 0.1j) (-0.099834078899207618+1.5707963277899337j) np.arccosh(1e-9 - 0.1j) (0.099834078899207576-1.5707963257998594j)
Ticket #854. http://scipy.org/scipy/numpy/ticket/854 I'll write up some tests for all the functions with branch cuts to verify that the cuts and their continuity are correct. (Where "correct" bears some resemblance to "ISO C standard", I think...) -- Pauli Virtanen
On Thu, Jul 17, 2008 at 3:56 PM, Pauli Virtanen <pav@iki.fi> wrote:
Mon, 17 Mar 2008 08:07:38 -0600, Charles R Harris wrote: [clip]
OK, that does it. I'm going to change it's behavior.
The problem with bad arccosh branch cuts is still present:
import numpy as np numpy.__version__ '1.2.0.dev5436.e45a7627a39d' np.arccosh(-1e-9 + 0.1j) (-0.099834078899207618-1.5707963277899337j) np.arccosh(1e-9 + 0.1j) (0.099834078899207576+1.5707963257998594j) np.arccosh(-1e-9 - 0.1j) (-0.099834078899207618+1.5707963277899337j) np.arccosh(1e-9 - 0.1j) (0.099834078899207576-1.5707963257998594j)
Ticket #854. http://scipy.org/scipy/numpy/ticket/854
I'll write up some tests for all the functions with branch cuts to verify that the cuts and their continuity are correct. (Where "correct" bears some resemblance to "ISO C standard", I think...)
Hmm, The problem here is arccosh = log(x + sqrt(x**2 - 1)) when the given numbers are plugged into x**2 - 1, one lies above the negative real axis, the other below and the branch cut [-inf,0] of sqrt introduces the discontinuity. Maybe sqrt(x - 1)*sqrt(x+1) will fix that. I do think the branch cut should be part of the documentation of all the complex functions. I wonder what arccos does here? Ah, here is a reference. <http://www.cs.umbc.edu/help/theory/identities.txt>Note arccosh z = ln(z + sqrt(z-1) sqrt(z+1) ) not sqrt(z**2-1) So I guess that is the fix. Chuck
Hi all, Re: Ticket 854. I wrote tests for the branch cuts for all complex arc* functions in umathmodule. It turns out that all except arccosh were OK. The formula for arcsinh was written in a non-standard form with an unnecessary nc_neg, but this didn't affect the results. I also wrote tests for checking values of the functions at infs and nans. A patch for all of this is attached, with all currently non-passing tests marked as skipped. I'd like to commit this if there are no objections. Another thing I noticed is that the present implementations of the complex functions are naive, so they over- and underflow earlier than necessary:
np.arcsinh(1e8) 19.113827924512311 np.arcsinh(1e8 + 0j) (inf-0j) np.arcsinh(-1e8 + 0j) (-19.113827924512311-0j)
This particular thing in arcsinh occurs because of loss of precision in intermediate stages. (In the version in my patch this loss of precision is still present.) It would be nice to polish these up. BTW, are there obstacles to using the C99 complex functions when they are available? This would avoid quite a bit of drudgework... As an alternative, we can probably steal better implementations from Python's recently polished cmathmodule.c *** Then comes a descent into pedantry: Numpy complex functions are not C99 compliant in handling of the signed zero, inf, and nan. I don't know whether we should comply with C99, but it would be nicer to have these handled in a controlled way rather than as a byproduct of the implementation chosen. 1) The branch cuts for sqrt and arc* don't respect the negative zero:
a = 1 + 0j np.sqrt(-a) 1j np.sqrt(-1 + 0j) 1j
The branch cut of the logarithm however does:
np.log(-a) -3.1415926535897931j np.log(-1 + 0j) 3.1415926535897931j
All complex functions in the C99 standard respect the negative zero in their branch cuts. Do we want to follow? I don't know how to check what Octave and Matlab do regarding this, since I haven't figured out how to place a negative zero in complex numbers in these languages. But at least in practice these languages appear not to respect the sign of zero.
a = 1 + 0j log(-a) ans = 0.000000000000000 + 3.141592653589793i log(-1) ans = 0.000000000000000 + 3.141592653589793i
2) The numpy functions in general don't return C99 compliant results at inf or nan. I wrote up some tests for checking these. Do we want to fix these? -- Pauli Virtanen diff -r 357b6ce2a4bc numpy/core/src/umathmodule.c.src --- a/numpy/core/src/umathmodule.c.src Thu Jul 17 16:17:45 2008 +0300 +++ b/numpy/core/src/umathmodule.c.src Sat Jul 19 15:16:37 2008 +0300 @@ -825,15 +825,16 @@ c@typ@ t; nc_sum@c@(x, &nc_1@c@, &t); + nc_sqrt@c@(&t, &t); nc_diff@c@(x, &nc_1@c@, r); + nc_sqrt@c@(r, r); nc_prod@c@(&t, r, r); - nc_sqrt@c@(r, r); nc_sum@c@(x, r, r); nc_log@c@(r, r); return; /* return nc_log(nc_sum(x, - nc_sqrt(nc_prod(nc_sum(x,nc_1), nc_diff(x,nc_1))))); + nc_prod(nc_sqrt(nc_sum(x,nc_1)), nc_sqrt(nc_diff(x,nc_1))))); */ } @@ -863,12 +864,11 @@ nc_prod@c@(x, x, r); nc_sum@c@(&nc_1@c@, r, r); nc_sqrt@c@(r, r); - nc_diff@c@(r, x, r); + nc_sum@c@(r, x, r); nc_log@c@(r, r); - nc_neg@c@(r, r); return; /* - return nc_neg(nc_log(nc_diff(nc_sqrt(nc_sum(nc_1,nc_prod(x,x))),x))); + return nc_log(nc_sum(nc_sqrt(nc_sum(nc_1,nc_prod(x,x))),x)); */ } diff -r 357b6ce2a4bc numpy/core/tests/test_umath.py --- a/numpy/core/tests/test_umath.py Thu Jul 17 16:17:45 2008 +0300 +++ b/numpy/core/tests/test_umath.py Sat Jul 19 15:16:37 2008 +0300 @@ -1,6 +1,8 @@ from numpy.testing import * import numpy.core.umath as ncu import numpy as np +import nose +from numpy import inf, nan, pi class TestDivision(TestCase): def test_division_int(self): @@ -34,7 +36,6 @@ ncu.sqrt(3+4j)]) assert_almost_equal(x**14, [-76443+16124j, 23161315+58317492j, 5583548873 + 2465133864j]) - class TestLog1p(TestCase): def test_log1p(self): @@ -179,10 +180,10 @@ assert_equal(np.choose(c, (a, 1)), np.array([1,1])) -class TestComplexFunctions(TestCase): +class TestComplexFunctions(object): funcs = [np.arcsin , np.arccos , np.arctan, np.arcsinh, np.arccosh, np.arctanh, np.sin , np.cos , np.tan , np.exp, - np.log , np.sqrt , np.log10] + np.log , np.sqrt , np.log10, np.log1p] def test_it(self): for f in self.funcs: @@ -204,6 +205,205 @@ assert_almost_equal(fcf, fcd, decimal=6, err_msg='fch-fcd %s'%f) assert_almost_equal(fcl, fcd, decimal=15, err_msg='fch-fcl %s'%f) + def test_branch_cuts(self): + # check branch cuts and continuity on them + yield _check_branch_cut, np.log, -0.5, 1j, 1, -1, True + yield _check_branch_cut, np.log10, -0.5, 1j, 1, -1, True + yield _check_branch_cut, np.log1p, -1.5, 1j, 1, -1, True + yield _check_branch_cut, np.sqrt, -0.5, 1j, 1, -1 + + yield _check_branch_cut, np.arcsin, [ -2, 2], [1j, -1j], 1, -1 + yield _check_branch_cut, np.arccos, [ -2, 2], [1j, -1j], 1, -1 + yield _check_branch_cut, np.arctan, [-2j, 2j], [1, -1 ], -1, 1 + + yield _check_branch_cut, np.arcsinh, [-2j, 2j], [-1, 1], -1, 1 + yield _check_branch_cut, np.arccosh, [ -1, 0.5], [1j, 1j], 1, -1 + yield _check_branch_cut, np.arctanh, [ -2, 2], [1j, -1j], 1, -1 + + # check against bogus branch cuts: assert continuity between quadrants + yield _check_branch_cut, np.arcsin, [-2j, 2j], [ 1, 1], 1, 1 + yield _check_branch_cut, np.arccos, [-2j, 2j], [ 1, 1], 1, 1 + yield _check_branch_cut, np.arctan, [ -2, 2], [1j, 1j], 1, 1 + + yield _check_branch_cut, np.arcsinh, [ -2, 2, 0], [1j, 1j, 1 ], 1, 1 + yield _check_branch_cut, np.arccosh, [-2j, 2j, 2], [1, 1, 1j], 1, 1 + yield _check_branch_cut, np.arctanh, [-2j, 2j, 0], [1, 1, 1j], 1, 1 + + def test_branch_cuts_failing(self): + # XXX: signed zeros are not OK for sqrt or for the arc* functions + yield _check_branch_cut, np.sqrt, -0.5, 1j, 1, -1, True + yield _check_branch_cut, np.arcsin, [ -2, 2], [1j, -1j], 1, -1, True + yield _check_branch_cut, np.arccos, [ -2, 2], [1j, -1j], 1, -1, True + yield _check_branch_cut, np.arctan, [-2j, 2j], [1, -1 ], -1, 1, True + yield _check_branch_cut, np.arcsinh, [-2j, 2j], [-1, 1], -1, 1, True + yield _check_branch_cut, np.arccosh, [ -1, 0.5], [1j, 1j], 1, -1, True + yield _check_branch_cut, np.arctanh, [ -2, 2], [1j, -1j], 1, -1, True + test_branch_cuts_failing = dec.skipknownfailure(test_branch_cuts_failing) + + def test_against_cmath(self): + import cmath, sys + + # cmath.asinh is broken in some versions of Python, see + # http://bugs.python.org/issue1381 + broken_cmath_asinh = False + if sys.version_info < (2,5,3): + broken_cmath_asinh = True + + points = [-2, 2j, 2, -2j, -1-1j, -1+1j, +1-1j, +1+1j] + name_map = {'arcsin': 'asin', 'arccos': 'acos', 'arctan': 'atan', + 'arcsinh': 'asinh', 'arccosh': 'acosh', 'arctanh': 'atanh'} + atol = 4*np.finfo(np.complex).eps + for func in self.funcs: + fname = func.__name__.split('.')[-1] + cname = name_map.get(fname, fname) + try: cfunc = getattr(cmath, cname) + except AttributeError: continue + for p in points: + a = complex(func(np.complex_(p))) + b = cfunc(p) + + if cname == 'asinh' and broken_cmath_asinh: + continue + + assert abs(a - b) < atol, "%s %s: %s; cmath: %s"%(fname,p,a,b) + +class TestC99(object): + """Check special functions at special points against the C99 standard""" + # NB: inherits from object instead of TestCase since using test generators + + # + # Non-conforming results are with XXX added to the exception field. + # + + def test_clog(self): + for p, v, e in [ + ((-0., 0.), (-inf, pi), 'divide'), + ((+0., 0.), (-inf, 0.), 'divide'), + ((1., inf), (inf, pi/2), ''), + ((1., nan), (nan, nan), ''), + ((-inf, 1.), (inf, pi), ''), + ((inf, 1.), (inf, 0.), ''), + ((-inf, inf), (inf, 3*pi/4), ''), + ((inf, inf), (inf, pi/4), ''), + ((inf, nan), (inf, nan), ''), + ((-inf, nan), (inf, nan), ''), + ((nan, 0.), (nan, nan), ''), + ((nan, 1.), (nan, nan), ''), + ((nan, inf), (inf, nan), ''), + ((+nan, nan), (nan, nan), ''), + ]: + yield self._check, np.log, p, v, e + + def test_csqrt(self): + for p, v, e in [ + ((-0., 0.), (0.,0.), 'XXX'), # now (-0., 0.) + ((0., 0.), (0.,0.), ''), + ((1., inf), (inf,inf), 'XXX invalid'), # now (inf, nan) + ((nan, inf), (inf,inf), 'XXX'), # now (nan, nan) + ((-inf, 1.), (0.,inf), ''), + ((inf, 1.), (inf,0.), ''), + ((-inf,nan), (nan, -inf), ''), # could also be +inf + ((inf, nan), (inf, nan), ''), + ((nan, 1.), (nan, nan), ''), + ((nan, nan), (nan, nan), ''), + ]: + yield self._check, np.sqrt, p, v, e + + def test_cacos(self): + for p, v, e in [ + ((0., 0.), (pi/2, -0.), 'XXX'), # now (-0., 0.) + ((-0., 0.), (pi/2, -0.), ''), + ((0., nan), (pi/2, nan), 'XXX'), # now (nan, nan) + ((-0., nan), (pi/2, nan), 'XXX'), # now (nan, nan) + ((1., inf), (pi/2, -inf), 'XXX'), # now (nan, -inf) + ((1., nan), (nan, nan), ''), + ((-inf, 1.), (pi, -inf), 'XXX'), # now (nan, -inf) + ((inf, 1.), (0., -inf), 'XXX'), # now (nan, -inf) + ((-inf, inf), (3*pi/4, -inf), 'XXX'), # now (nan, nan) + ((inf, inf), (pi/4, -inf), 'XXX'), # now (nan, nan) + ((inf, nan), (nan, +-inf), 'XXX'), # now (nan, nan) + ((-inf, nan), (nan, +-inf), 'XXX'), # now: (nan, nan) + ((nan, 1.), (nan, nan), ''), + ((nan, inf), (nan, -inf), 'XXX'), # now: (nan, nan) + ((nan, nan), (nan, nan), ''), + ]: + yield self._check, np.arccos, p, v, e + + def test_cacosh(self): + for p, v, e in [ + ((0., 0), (0, pi/2), ''), + ((-0., 0), (0, pi/2), ''), + ((1., inf), (inf, pi/2), 'XXX'), # now: (nan, nan) + ((1., nan), (nan, nan), ''), + ((-inf, 1.), (inf, pi), 'XXX'), # now: (inf, nan) + ((inf, 1.), (inf, 0.), 'XXX'), # now: (inf, nan) + ((-inf, inf), (inf, 3*pi/4), 'XXX'), # now: (nan, nan) + ((inf, inf), (inf, pi/4), 'XXX'), # now: (nan, nan) + ((inf, nan), (inf, nan), 'XXX'), # now: (nan, nan) + ((-inf, nan), (inf, nan), 'XXX'), # now: (nan, nan) + ((nan, 1.), (nan, nan), ''), + ((nan, inf), (inf, nan), 'XXX'), # now: (nan, nan) + ((nan, nan), (nan, nan), '') + ]: + yield self._check, np.arccosh, p, v, e + + def test_casinh(self): + for p, v, e in [ + ((0., 0), (0, 0), ''), + ((1., inf), (inf, pi/2), 'XXX'), # now: (inf, nan) + ((1., nan), (nan, nan), ''), + ((inf, 1.), (inf, 0.), 'XXX'), # now: (inf, nan) + ((inf, inf), (inf, pi/4), 'XXX'), # now: (nan, nan) + ((inf, nan), (nan, nan), 'XXX'), # now: (nan, nan) + ((nan, 0.), (nan, 0.), 'XXX'), # now: (nan, nan) + ((nan, 1.), (nan, nan), ''), + ((nan, inf), (+-inf, nan), 'XXX'), # now: (nan, nan) + ((nan, nan), (nan, nan), ''), + ]: + yield self._check, np.arcsinh, p, v, e + + def test_catanh(self): + for p, v, e in [ + ((0., 0), (0, 0), ''), + ((0., nan), (0., nan), 'XXX'), # now: (nan, nan) + ((1., 0.), (inf, 0.), 'XXX divide'), # now: (nan, nan) + ((1., inf), (inf, 0.), 'XXX'), # now: (nan, nan) + ((1., nan), (nan, nan), ''), + ((inf, 1.), (0., pi/2), 'XXX'), # now: (nan, nan) + ((inf, inf), (0, pi/2), 'XXX'), # now: (nan, nan) + ((inf, nan), (0, nan), 'XXX'), # now: (nan, nan) + ((nan, 1.), (nan, nan), ''), + ((nan, inf), (+0, pi/2), 'XXX'), # now: (nan, nan) + ((nan, nan), (nan, nan), ''), + ]: + yield self._check, np.arctanh, p, v, e + + def _check(self, func, point, value, exc=''): + if 'XXX' in exc: + raise nose.SkipTest + if isinstance(point, tuple): point = complex(*point) + if isinstance(value, tuple): value = complex(*value) + v = dict(divide='ignore', invalid='ignore', + over='ignore', under='ignore') + old_err = np.seterr(**v) + try: + # check sign of zero, nan, etc. + got = complex(func(point)) + got = "(%s, %s)" % (repr(got.real), repr(got.imag)) + expected = "(%s, %s)" % (repr(value.real), repr(value.imag)) + assert got == expected, (got, expected) + + # check exceptions + if exc in ('divide', 'invalid', 'over', 'under'): + v[exc] = 'raise' + np.seterr(**v) + assert_raises(FloatingPointError, func, point) + else: + for k in v.keys(): v[k] = 'raise' + np.seterr(**v) + func(point) + finally: + np.seterr(**old_err) class TestAttributes(TestCase): def test_attributes(self): @@ -216,6 +416,59 @@ assert_equal(add.nout, 1) assert_equal(add.identity, 0) +def _check_branch_cut(f, x0, dx, re_sign=1, im_sign=-1, sig_zero_ok=False, + dtype=np.complex): + """ + Check for a branch cut in a function. + + Assert that `x0` lies on a branch cut of function `f` and `f` is + continuous from the direction `dx`. + + Parameters + ---------- + f : func + Function to check + x0 : array-like + Point on branch cut + dx : array-like + Direction to check continuity in + re_sign, im_sign : {1, -1} + Change of sign of the real or imaginary part expected + sig_zero_ok : bool + Whether to check if the branch cut respects signed zero (if applicable) + dtype : dtype + Dtype to check (should be complex) + + """ + x0 = np.atleast_1d(x0).astype(dtype) + dx = np.atleast_1d(dx).astype(dtype) + + scale = np.finfo(dtype).eps * 1e3 + atol = 1e-4 + + y0 = f(x0) + yp = f(x0 + dx*scale*np.absolute(x0)/np.absolute(dx)) + ym = f(x0 - dx*scale*np.absolute(x0)/np.absolute(dx)) + + assert np.all(np.absolute(y0.real - yp.real) < atol), (y0, yp) + assert np.all(np.absolute(y0.imag - yp.imag) < atol), (y0, yp) + assert np.all(np.absolute(y0.real - ym.real*re_sign) < atol), (y0, ym) + assert np.all(np.absolute(y0.imag - ym.imag*im_sign) < atol), (y0, ym) + + if sig_zero_ok: + # check that signed zeros also work as a displacement + jr = (x0.real == 0) & (dx.real != 0) + ji = (x0.imag == 0) & (dx.imag != 0) + + x = -x0 + x.real[jr] = 0.*dx.real + x.imag[ji] = 0.*dx.imag + x = -x + ym = f(x) + ym = ym[jr | ji] + y0 = y0[jr | ji] + assert np.all(np.absolute(y0.real - ym.real*re_sign) < atol), (y0, ym) + assert np.all(np.absolute(y0.imag - ym.imag*im_sign) < atol), (y0, ym) if __name__ == "__main__": run_module_suite()
day pot-luck invite was a SHAM! The real party is on Saturday, and is not a pot-luck.Remember -- the Sunday pot-luck invite was a SHAM! The real party is on Saturday, and is not a pot-luck.Remember -- the Sunday pot-luck invite was a SHAM! The real party is on Saturday, and is not a pot-luck.Remember -- the Sunday pot-luck invite was a SHAM! The real party is on Saturday, and is not a pot-luck.
Accidental (virus?) post. Humblest apologies for the noise. Please ignore. Gary On Sat, 19 Jul 2008, Gary Strangman wrote:
day pot-luck invite was a SHAM! The real party is on Saturday, and is not a pot-luck.Remember -- the Sunday pot-luck invite was a SHAM! The real party is on Saturday, and is not a pot-luck.Remember -- the Sunday pot-luck invite was a SHAM! The real party is on Saturday, and is not a pot-luck.Remember -- the Sunday pot-luck invite was a SHAM! The real party is on Saturday, and is not a pot-luck. _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
On Sat, Jul 19, 2008 at 7:13 AM, Pauli Virtanen <pav@iki.fi> wrote:
Hi all,
Re: Ticket 854.
I wrote tests for the branch cuts for all complex arc* functions in umathmodule. It turns out that all except arccosh were OK. The formula for arcsinh was written in a non-standard form with an unnecessary nc_neg, but this didn't affect the results. I also wrote tests for checking values of the functions at infs and nans.
A patch for all of this is attached, with all currently non-passing tests marked as skipped. I'd like to commit this if there are no objections.
Another thing I noticed is that the present implementations of the complex functions are naive, so they over- and underflow earlier than necessary:
np.arcsinh(1e8) 19.113827924512311 np.arcsinh(1e8 + 0j) (inf-0j) np.arcsinh(-1e8 + 0j) (-19.113827924512311-0j)
This particular thing in arcsinh occurs because of loss of precision in intermediate stages. (In the version in my patch this loss of precision is still present.)
It would be nice to polish these up. BTW, are there obstacles to using the C99 complex functions when they are available? This would avoid quite a bit of drudgework... As an alternative, we can probably steal better implementations from Python's recently polished cmathmodule.c
The main problem is determining when they are available and if they cover all the needed precisions. Since we will need standalone implementations on some platforms anyway, I am inclined towards stealing from cmathmodule.c if it offers improved code for some of the functions.
***
Then comes a descent into pedantry: Numpy complex functions are not C99 compliant in handling of the signed zero, inf, and nan. I don't know whether we should comply with C99, but it would be nicer to have these handled in a controlled way rather than as a byproduct of the implementation chosen.
Agreed.
1)
The branch cuts for sqrt and arc* don't respect the negative zero:
a = 1 + 0j np.sqrt(-a) 1j np.sqrt(-1 + 0j) 1j
The branch cut of the logarithm however does:
np.log(-a) -3.1415926535897931j np.log(-1 + 0j) 3.1415926535897931j
All complex functions in the C99 standard respect the negative zero in their branch cuts. Do we want to follow?
Hmm. I think so, to the extent we can. This might lead to some unexpected results, but that is what happens for arguments near the branch cuts. Do it and document it.
I don't know how to check what Octave and Matlab do regarding this, since I haven't figured out how to place a negative zero in complex numbers in these languages. But at least in practice these languages appear not to respect the sign of zero.
a = 1 + 0j log(-a) ans = 0.000000000000000 + 3.141592653589793i log(-1) ans = 0.000000000000000 + 3.141592653589793i
2)
The numpy functions in general don't return C99 compliant results at inf or nan. I wrote up some tests for checking these.
Do we want to fix these?
I'd say yes. That way we can refer to the C99 standard to document numpy behavior. Chuck
Pauli Virtanen wrote:
Hi all,
Re: Ticket 854.
I wrote tests for the branch cuts for all complex arc* functions in umathmodule. It turns out that all except arccosh were OK. The formula for arcsinh was written in a non-standard form with an unnecessary nc_neg, but this didn't affect the results. I also wrote tests for checking values of the functions at infs and nans.
Thanks for looking into these. These functions were contributed by Konrad Hinsen (IIRC) many years ago and I don't think they've really been reviewed since then. I'm all for using C99 when it is available and improving these functions with help from cmathmodule. IIRC, the cmathmodule was contributed by Konrad originally also. So +1 on C99 standardization. -Travis
On Sat, 19 Jul 2008 13:29:45 -0500 "Travis E. Oliphant" <oliphant@enthought.com> wrote:
Pauli Virtanen wrote:
Hi all,
Re: Ticket 854.
I wrote tests for the branch cuts for all complex arc* functions in umathmodule. It turns out that all except arccosh were OK. The formula for arcsinh was written in a non-standard form with an unnecessary nc_neg, but this didn't affect the results. I also wrote tests for checking values of the functions at infs and nans.
Thanks for looking into these. These functions were contributed by Konrad Hinsen (IIRC) many years ago and I don't think they've really been reviewed since then.
I'm all for using C99 when it is available and improving these functions with help from cmathmodule. IIRC, the cmathmodule was contributed by Konrad originally also.
So +1 on C99 standardization.
-Travis
====================================================================== ERROR: test_umath.TestC99.test_catanh(<ufunc 'arctanh'>, (nan, nan), (nan, nan), '') ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/lib/python2.4/site-packages/nose-0.10.3-py2.4.egg/nose/case.py", line 182, in runTest self.test(*self.arg) File "/usr/lib/python2.4/site-packages/numpy/core/tests/test_umath.py", line 405, in _check func(point) FloatingPointError: invalid value encountered in arctanh ====================================================================== ERROR: test_umath.TestC99.test_clog(<ufunc 'log'>, (1.0, nan), (nan, nan), '') ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/lib/python2.4/site-packages/nose-0.10.3-py2.4.egg/nose/case.py", line 182, in runTest self.test(*self.arg) File "/usr/lib/python2.4/site-packages/numpy/core/tests/test_umath.py", line 405, in _check func(point) FloatingPointError: invalid value encountered in log ====================================================================== ERROR: test_umath.TestC99.test_clog(<ufunc 'log'>, (nan, 0.0), (nan, nan), '') ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/lib/python2.4/site-packages/nose-0.10.3-py2.4.egg/nose/case.py", line 182, in runTest self.test(*self.arg) File "/usr/lib/python2.4/site-packages/numpy/core/tests/test_umath.py", line 405, in _check func(point) FloatingPointError: invalid value encountered in log ====================================================================== ERROR: test_umath.TestC99.test_clog(<ufunc 'log'>, (nan, 1.0), (nan, nan), '') ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/lib/python2.4/site-packages/nose-0.10.3-py2.4.egg/nose/case.py", line 182, in runTest self.test(*self.arg) File "/usr/lib/python2.4/site-packages/numpy/core/tests/test_umath.py", line 405, in _check func(point) FloatingPointError: invalid value encountered in log ====================================================================== ERROR: test_umath.TestC99.test_clog(<ufunc 'log'>, (nan, nan), (nan, nan), '') ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/lib/python2.4/site-packages/nose-0.10.3-py2.4.egg/nose/case.py", line 182, in runTest self.test(*self.arg) File "/usr/lib/python2.4/site-packages/numpy/core/tests/test_umath.py", line 405, in _check func(point) FloatingPointError: invalid value encountered in log ---------------------------------------------------------------------- Ran 2070 tests in 35.424s FAILED (SKIP=41, errors=5) <nose.result.TextTestResult run=2070 errors=5 failures=0> Nils
On Sun, Jul 20, 2008 at 2:18 AM, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
On Sat, 19 Jul 2008 13:29:45 -0500 "Travis E. Oliphant" <oliphant@enthought.com> wrote:
Pauli Virtanen wrote:
Hi all,
Re: Ticket 854.
I wrote tests for the branch cuts for all complex arc* functions in umathmodule. It turns out that all except arccosh were OK. The formula for arcsinh was written in a non-standard form with an unnecessary nc_neg, but this didn't affect the results. I also wrote tests for checking values of the functions at infs and nans.
Thanks for looking into these. These functions were contributed by Konrad Hinsen (IIRC) many years ago and I don't think they've really been reviewed since then.
I'm all for using C99 when it is available and improving these functions with help from cmathmodule. IIRC, the cmathmodule was contributed by Konrad originally also.
So +1 on C99 standardization.
-Travis
====================================================================== ERROR: test_umath.TestC99.test_catanh(<ufunc 'arctanh'>, (nan, nan), (nan, nan), '') ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/lib/python2.4/site-packages/nose-0.10.3-py2.4.egg/nose/case.py", line 182, in runTest self.test(*self.arg) File "/usr/lib/python2.4/site-packages/numpy/core/tests/test_umath.py", line 405, in _check func(point) FloatingPointError: invalid value encountered in arctanh
<snip> What architecture and OS? Chuck
On Sun, Jul 20, 2008 at 03:34, Charles R Harris <charlesr.harris@gmail.com> wrote:
What architecture and OS?
I get the following on OS X 10.5.3 Intel Core 2 Duo: ====================================================================== FAIL: test_umath.TestC99.test_clog(<ufunc 'log'>, (-0.0, -0.0), (-inf, -0.0), 'divide') ---------------------------------------------------------------------- Traceback (most recent call last): File "/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/nose-0.10.3-py2.5.egg/nose/case.py", line 182, in runTest self.test(*self.arg) File "/Users/rkern/svn/numpy/numpy/core/tests/test_umath.py", line 394, in _check assert got == expected, (got, expected) AssertionError: ('(-inf, 3.1415926535897931)', '(-inf, 0.0)') ====================================================================== FAIL: test_umath.TestC99.test_csqrt(<ufunc 'sqrt'>, (-inf, 1.0), (-0.0, inf), '') ---------------------------------------------------------------------- Traceback (most recent call last): File "/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/nose-0.10.3-py2.5.egg/nose/case.py", line 182, in runTest self.test(*self.arg) File "/Users/rkern/svn/numpy/numpy/core/tests/test_umath.py", line 394, in _check assert got == expected, (got, expected) AssertionError: ('(0.0, inf)', '(-0.0, inf)') ---------------------------------------------------------------------- Ran 1887 tests in 7.716s -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
On Sun, 20 Jul 2008 02:34:48 -0600 "Charles R Harris" <charlesr.harris@gmail.com> wrote:
On Sun, Jul 20, 2008 at 2:18 AM, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
On Sat, 19 Jul 2008 13:29:45 -0500 "Travis E. Oliphant" <oliphant@enthought.com> wrote:
Pauli Virtanen wrote:
Hi all,
Re: Ticket 854.
I wrote tests for the branch cuts for all complex arc* functions in umathmodule. It turns out that all except arccosh were OK. The formula for arcsinh was written in a non-standard form with an unnecessary nc_neg, but this didn't affect the results. I also wrote tests for checking values of the functions at infs and nans.
Thanks for looking into these. These functions were contributed by Konrad Hinsen (IIRC) many years ago and I don't think they've really been reviewed since then.
I'm all for using C99 when it is available and improving these functions with help from cmathmodule. IIRC, the cmathmodule was contributed by Konrad originally also.
So +1 on C99 standardization.
-Travis
====================================================================== ERROR: test_umath.TestC99.test_catanh(<ufunc 'arctanh'>, (nan, nan), (nan, nan), '') ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/lib/python2.4/site-packages/nose-0.10.3-py2.4.egg/nose/case.py", line 182, in runTest self.test(*self.arg) File "/usr/lib/python2.4/site-packages/numpy/core/tests/test_umath.py", line 405, in _check func(point) FloatingPointError: invalid value encountered in arctanh
<snip>
What architecture and OS?
Linux linux 2.6.11.4-21.17-default #1 Fri Apr 6 08:42:34 UTC 2007 i686 athlon i386 GNU/Linux SuSe Linux 9.3 gcc --version gcc (GCC) 3.3.5 20050117 (prerelease) (SUSE Linux) Copyright (C) 2003 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. /usr/bin/python Python 2.4 (#1, Oct 13 2006, 17:13:31) [GCC 3.3.5 20050117 (prerelease) (SUSE Linux)] on linux2 Type "help", "copyright", "credits" or "license" for more information. Nils
Hi, Sorry, Sun, 20 Jul 2008 10:18:47 +0200, Nils Wagner wrote:
ERROR: test_umath.TestC99.test_catanh(<ufunc 'arctanh'>, (nan, nan), (nan, nan), '') FloatingPointError: invalid value encountered in arctanh
Skipped now.
ERROR: test_umath.TestC99.test_clog(<ufunc 'log'>, (1.0, nan), (nan, nan), '') FloatingPointError: invalid value encountered in log
Bug in tests, fixed.
ERROR: test_umath.TestC99.test_clog(<ufunc 'log'>, (nan, 0.0), (nan, nan), '') FloatingPointError: invalid value encountered in log
Bug in tests, fixed.
ERROR: test_umath.TestC99.test_clog(<ufunc 'log'>, (nan, 1.0), (nan, nan), '') FloatingPointError: invalid value encountered in log
Bug in tests, fixed.
ERROR: test_umath.TestC99.test_clog(<ufunc 'log'>, (nan, nan), (nan, nan), '') FloatingPointError: invalid value encountered in log
Skipped. Sun, 20 Jul 2008 03:42:45 -0500, Robert Kern wrote:
FAIL: test_umath.TestC99.test_clog(<ufunc 'log'>, (-0.0, -0.0), (-inf, -0.0), 'divide')
Interesting, there's no test like this in there, only one with positive zeros. Where does this one come from?
FAIL: test_umath.TestC99.test_csqrt(<ufunc 'sqrt'>, (-inf, 1.0), (-0.0, inf), '')
Fails on your platform, skipped. -- Pauli Virtanen
On Sun, 20 Jul 2008 12:10:23 +0000 (UTC) Pauli Virtanen <pav@iki.fi> wrote:
Hi,
Sorry,
Sun, 20 Jul 2008 10:18:47 +0200, Nils Wagner wrote:
ERROR: test_umath.TestC99.test_catanh(<ufunc 'arctanh'>, (nan, nan), (nan, nan), '') FloatingPointError: invalid value encountered in arctanh
Skipped now.
ERROR: test_umath.TestC99.test_clog(<ufunc 'log'>, (1.0, nan), (nan, nan), '') FloatingPointError: invalid value encountered in log
Bug in tests, fixed.
ERROR: test_umath.TestC99.test_clog(<ufunc 'log'>, (nan, 0.0), (nan, nan), '') FloatingPointError: invalid value encountered in log
Bug in tests, fixed.
ERROR: test_umath.TestC99.test_clog(<ufunc 'log'>, (nan, 1.0), (nan, nan), '') FloatingPointError: invalid value encountered in log
Bug in tests, fixed.
ERROR: test_umath.TestC99.test_clog(<ufunc 'log'>, (nan, nan), (nan, nan), '') FloatingPointError: invalid value encountered in log
Skipped.
Sun, 20 Jul 2008 03:42:45 -0500, Robert Kern wrote:
FAIL: test_umath.TestC99.test_clog(<ufunc 'log'>, (-0.0, -0.0), (-inf, -0.0), 'divide')
Interesting, there's no test like this in there, only one with positive zeros. Where does this one come from?
FAIL: test_umath.TestC99.test_csqrt(<ufunc 'sqrt'>, (-inf, 1.0), (-0.0, inf), '')
Fails on your platform, skipped.
-- Pauli Virtanen
_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Ran 2069 tests in 30.300s OK (SKIP=44) <nose.result.TextTestResult run=2069 errors=0 failures=0> Thank you. Cheers, Nils
participants (7)
-
Charles R Harris -
Gary Strangman -
lorenzo bolla -
Nils Wagner -
Pauli Virtanen -
Robert Kern -
Travis E. Oliphant