classification
Title: cmath is numerically unsound
Type: behavior Stage:
Components: Extension Modules Versions: Python 3.0, Python 3.1, Python 2.6
process
Status: closed Resolution: fixed
Dependencies: Superseder:
Assigned To: Nosy List: alanmcintyre, bins, christian.heimes, gvanrossum, inducer, jcea, loewis, mark.dickinson
Priority: normal Keywords: patch

Created on 2007-11-03 18:58 by inducer, last changed 2010-05-11 15:01 by bins. This issue is now closed.

Files
File name Uploaded Description Edit
cmath_py.py mark.dickinson, 2007-11-04 15:16
cmath.patch mark.dickinson, 2007-11-26 02:18
cmath3.patch mark.dickinson, 2008-01-05 00:25 Updated patch for cmath module, with docs and tests
Messages (23)
msg57088 - (view) Author: Andreas Kloeckner (inducer) Date: 2007-11-03 18:58
This here basically says it all:

>>> import cmath;[cmath.asinh(i*1e-17).real for i in range(0,20)]
[4.4408920985006257e-16, 4.4408920985006257e-16, 4.4408920985006257e-16,
4.4408920985006257e-16, 4.4408920985006257e-16, 4.4408920985006257e-16,
4.4408920985006257e-16, 4.4408920985006257e-16, 4.4408920985006257e-16,
4.4408920985006257e-16, 4.4408920985006257e-16, 4.4408920985006257e-16,
4.4408920985006257e-16, 4.4408920985006257e-16, 4.4408920985006257e-16,
4.4408920985006257e-16, 4.4408920985006257e-16, 4.4408920985006257e-16,
4.4408920985006257e-16, 4.4408920985006257e-16]

The boost.math toolkit at [2] is an implementation that does better in
the above (real-only) aspect.
[2] http://freespace.virgin.net/boost.regex/toolkit/html/index.html

Tim Peters remarks in [1] that basically all of cmath is unsound.
http://mail.python.org/pipermail/python-bugs-list/2001-February/004126.html

I just wanted to make sure that this issue remains on the radar.
msg57089 - (view) Author: Martin v. Löwis (loewis) * (Python committer) Date: 2007-11-03 19:17
Can you propose a patch?
msg57090 - (view) Author: Andreas Kloeckner (inducer) Date: 2007-11-03 19:42
On Samstag 03 November 2007, Martin v. Löwis wrote:
> Martin v. Löwis added the comment:
>
> Can you propose a patch?

Other than point at how boost.math does things, I don't have the time to work 
on this right now, sorry.

Andreas
msg57092 - (view) Author: Alan McIntyre (alanmcintyre) * (Python committer) Date: 2007-11-04 05:46
I have to review a few complex math topics for some of my current course
work, so I wouldn't mind taking a look into this.  I can't promise I'll
have the time required to make all of cmath correct (assuming it's as
unsound as claimed), but I'll do what I can.
msg57098 - (view) Author: Mark Dickinson (mark.dickinson) * (Python committer) Date: 2007-11-04 15:16
I took a look at this a while back, and got as far as writing a pure 
Python drop-in replacement for cmath, based on Kahan's "branch cuts for 
elementary functions" paper.  This fixes a variety of problems in cmath, 
including the buggy branch cuts for asinh.  File attached, in case it's 
of any use.

As Tim Peters pointed out, the real problem here is a lack of decent 
unit tests for these functions.  I have tests for the file above, but 
they assume IEEE754 floats, which is probably not an acceptable 
assumption in general.
msg57104 - (view) Author: Martin v. Löwis (loewis) * (Python committer) Date: 2007-11-04 18:38
It would be ok if a test is only run on a system with IEEE floats, and
skipped elsewhere. For all practical purposes, Python assumes that all
systems have IEEE float.
msg57835 - (view) Author: Mark Dickinson (mark.dickinson) * (Python committer) Date: 2007-11-26 02:18
Here is (quite a large) patch, cmath.patch, that fixes a variety of 
problems in the current cmath module.  A summary of the changes:

* sqrt, log, acos, acosh, asin, asinh, atan, atanh no longer produce 
overflow errors for very large inputs

* exp, cos, sin, tan, cosh, sinh, tanh produce valid answers in some cases
where they incorrectly overflowed before.

* sqrt and log are more accurate for tiny numbers

* numeric problems in some functions (especially asinh and asin) should
have been fixed

* the branch cuts for asinh have been moved to the standard positions

* the direction of continuity for the branch cuts of tan, tanh have been 
fixed

* on systems with full hardware and system support for signed zeros (most 
modern systems), functions with a branch cut are continuous on both sides 
of the branch cut.  (As recommended by the C99 standard, amongst others.)

The patch includes documentation updates, but there are still no tests.  
(My current tests make heavy use of the MPFR library, and assume IEEE-754 
format doubles, so need to be seriously reworked.)  The tests are on my to-
do list, but I'm unlikely to find time for them before the new year.  In 
the meantime, I'd very much appreciate any feedback on this patch, if 
anyone has the time/energy/inclination to look at it.  (Andreas:  are you 
in a position to give this a test-run?)
msg59268 - (view) Author: Mark Dickinson (mark.dickinson) * (Python committer) Date: 2008-01-05 00:25
Here's an updated patch for cmath, complete with tests and documentation changes.
There's an extra file, Lib/test/cmath.ctest, containing test values for the various functions;  
these test values were obtained using MPFR and interval arithmetic, and (modulo bugs) should 
all be correctly rounded.

Any feedback would be much appreciated.
msg59557 - (view) Author: Mark Dickinson (mark.dickinson) * (Python committer) Date: 2008-01-08 20:55
A note to any potential reviewers for this patch (cmath3.patch):  I'm especially 
looking for non-mathematical feedback here, so please don't be put off by the 
mathematics.  Some questions:

- there's a new file cmath.ctest containing testcases;  should this be put 
directly into Lib/test (where it is right now), or is there a better place.  Is 
the name of this file reasonable?

- is the new stuff in pyport.h (checks for _copysign and copysign, and automatic 
renaming of _copysign to copysign) okay?  Should it be in cmathmodule.c instead?

- is the code in test_cmath that looks for the testcase file okay?

And it would be really great if someone with access to a Windows box could just 
verify that the tests pass with this patch applied;  I've only tested it on OS X 
and Linux.
msg59568 - (view) Author: Christian Heimes (christian.heimes) * (Python committer) Date: 2008-01-08 22:42
* I've added a configure test for copysign a while ago. I'm using this
constructor for Windows compatibility in mathmodule.c:

#ifdef MS_WINDOWS
#  define copysign _copysign
#  define HAVE_COPYSIGN 1
#endif
#ifdef HAVE_COPYSIGN
#endif

I know no platform besides Windows with a _copysign method. You don't
have to add tests for _copysign for Windows. You may want to add the
code to PC/pyconfig.h and define HAVE_COPYSIGN there.

* test_file = testdir + os.sep + 'cmath.ctest
  make that:
  test_file = os.path.join(testdir, 'cmath.ctest')

* Windows doesn't have log1p, the hyperbolic area functions (asinh) and
some other functions, too. IIRC Windows does only implement the Bessel
functions of 1st and 2nd kind but non of the other C99 math functions.
msg59570 - (view) Author: Mark Dickinson (mark.dickinson) * (Python committer) Date: 2008-01-08 22:54
Okay:  would it be okay for me to move the #ifdef MS_WINDOWS sequence 
somewhere where it can be used by both mathmodule.c and cmathmodule.c, 
then?

There are replacements for log1p and asinh in the patch, so it shouldn't 
matter than they're missing on Windows---it might mean that the accuracy 
of some of the functions is slightly reduced on Windows, though.
msg59572 - (view) Author: Christian Heimes (christian.heimes) * (Python committer) Date: 2008-01-08 23:16
It may even be a good idea to move asinh, copysign and log1p to a new
file Python/pymath.c and add compensatory implementations of acosh,
atanh and expm1, too. The math related code in pyport.h could then be
moved to Include/pymath.h. In the past some people have asked me to
implement the reverse hyperbolic functions for the math module.

This is beyond a simple fix for the cmath module and should be discussed
on the python-dev list first.
msg59637 - (view) Author: Christian Heimes (christian.heimes) * (Python committer) Date: 2008-01-10 00:20
Guido, how do you like the idea of Include/pymath.h and Python/pymath.c
containing supplementary mathematical functions and mathematical
constants? Mark's patch for cmath is rather large, can it still be
applied to 2.5?
msg59638 - (view) Author: Guido van Rossum (gvanrossum) * (Python committer) Date: 2008-01-10 00:29
Are you crazy?  Adding new features to 2.5?  No way!
msg61499 - (view) Author: Christian Heimes (christian.heimes) * (Python committer) Date: 2008-01-22 11:16
See #1640 and svn+ssh://pythondev@svn.python.org/python/branches/trunk-math
msg65703 - (view) Author: Mark Dickinson (mark.dickinson) * (Python committer) Date: 2008-04-23 17:08
Substantial fixes for the cmath module went into the trunk and the py3k 
branches as part of the merge of the trunk-math branch.  These fixes 
address the asinh problems in this issue, amongst other things.  See 
r62380 (trunk) and r62384 (py3k).
msg105510 - (view) Author: Sebastien Binet (bins) Date: 2010-05-11 13:29
hi there,

it seems there is still a problem, at least with asinh(-2j)

see:

$ python
Python 2.6.5 (r265:79063, Apr  1 2010, 05:28:39) 
[GCC 4.4.3 20100316 (prerelease)] on linux2
Type "help", "copyright", "credits" or "license" for more information.

>>> import numpy as np
>>> np.__version__
'1.4.0'
>>> import cmath
>>> cmath.asinh(-2j)
 (1.3169578969248166-1.5707963267948966j)

>>> np.arcsinh(-2j)
(-1.3169578969248164-1.5707963267948966j)

same behaviour for python3.1:
$ python3
Python 3.1.2 (r312:79147, Apr  1 2010, 09:12:21) 
[GCC 4.4.3 20100316 (prerelease)] on linux2
Type "help", "copyright", "credits" or "license" for more information.

>>> import cmath
>>> cmath.asinh(-2j)
(1.3169578969248166-1.5707963267948966j)

cheers,
sebastien.
msg105512 - (view) Author: Mark Dickinson (mark.dickinson) * (Python committer) Date: 2010-05-11 13:39
Python's result looks fine to me, as does numpy's: they're both giving a valid inverse hyperbolic sine:

>>> from cmath import sinh
>>> sinh(1.3169578969248166-1.5707963267948966j)  
(1.0605752387249067e-16-1.9999999999999998j)
>>> sinh(-1.3169578969248164-1.5707963267948966j)
(-1.0605752387249064e-16-1.9999999999999993j)

Perhaps numpy is using different branch cuts?
msg105515 - (view) Author: Mark Dickinson (mark.dickinson) * (Python committer) Date: 2010-05-11 14:00
A bit more explanation:  Python takes account of the sign of zero when deciding which side of the branch cut a value lies, which is the proper thing to do when you have signed zeros available (according to the likes of Kahan, anyway);  I suspect that numpy isn't doing that, but is treating all values that lie directly on a branch in the same way.

In this case there's a branch cut from -1j down to -1j*inf.  Values just to the right of that branch cut (i.e., with positive real part) should have a result with positive real part;  values just to the left of it should have negative real part:

Some results (using complex() to create the values, since other ways of creating complex numbers are prone to changing the sign of a zero):

>>> asinh(complex(0.0, -2.0))
(1.3169578969248166-1.5707963267948966j)
>>> asinh(complex(1e-10, -2.0))
(1.3169578969248166-1.5707963267371616j)
>>> asinh(complex(-0.0, -2.0))
(-1.3169578969248166-1.5707963267948966j)
>>> asinh(complex(-1e-10, -2.0))
(-1.3169578969248166-1.5707963267371616j)

So the cmath module is working as intended here.  numpy may or may not be working as intended:  I don't know how much they care about branch cut continuity.
msg105517 - (view) Author: Sebastien Binet (bins) Date: 2010-05-11 14:05
hi Mark,

that may very well be so, but I'd naively standardize on C/Fortran behaviour (but that's probably my physicist bias)

on my platform, the following piece of C-code:

$ cat test_cmath.c
#include <complex.h>
#include <stdio.h>

int main(int argc, char** argv)
{
  complex c = casinh(-2*I);
  printf("asinh(-2j) = %g + %gi\n", creal(c), cimag(c));
  return 0;
}
/* EOF */

gives:
$ ./a.out 
asinh(-2j) = -1.31696 + -1.5708i

cheers,
sebastien.
msg105519 - (view) Author: Mark Dickinson (mark.dickinson) * (Python committer) Date: 2010-05-11 14:08
> that may very well be so, but I'd naively standardize on C/Fortran > behaviour (but that's probably my physicist bias)

Yep, that's exactly what Python does. :)  (Also follows the LISP standard).

Note that in your program, you're feeding complex(-0.0, -2.0) to asinh,
not complex(0.0, 2.0).
msg105521 - (view) Author: Mark Dickinson (mark.dickinson) * (Python committer) Date: 2010-05-11 14:11
> Note that in your program, you're feeding complex(-0.0, -2.0) to asinh,
> not complex(0.0, 2.0).

Bah;  that should be complex(0.0, -2.0) in the second line, of course.

Anyway, try passing conj(2*I) to asinh in your C program and see what happens.  :)
msg105522 - (view) Author: Sebastien Binet (bins) Date: 2010-05-11 15:01
> Note that in your program, you're feeding complex(-0.0, -2.0) to asinh,
> not complex(0.0, -2.0).

ah!
(ducking)
History
Date User Action Args
2010-05-11 15:01:04binssetmessages: + msg105522
2010-05-11 14:11:48mark.dickinsonsetmessages: + msg105521
2010-05-11 14:08:23mark.dickinsonsetmessages: + msg105519
2010-05-11 14:05:16binssetmessages: + msg105517
2010-05-11 14:00:15mark.dickinsonsetmessages: + msg105515
2010-05-11 13:39:44mark.dickinsonsetmessages: + msg105512
2010-05-11 13:29:38binssetnosy: + bins

messages: + msg105510
versions: + Python 3.1
2008-10-13 11:57:09jceasetnosy: + jcea
2008-04-23 17:08:48mark.dickinsonsetstatus: open -> closed
resolution: fixed
messages: + msg65703
2008-01-22 11:17:00christian.heimessetpriority: normal
keywords: + patch
messages: + msg61499
components: + Extension Modules, - Library (Lib)
versions: + Python 3.0, - Python 2.5
2008-01-10 00:29:05gvanrossumsetmessages: + msg59638
2008-01-10 00:20:00christian.heimessetnosy: + gvanrossum
messages: + msg59637
versions: + Python 2.6
2008-01-08 23:16:54christian.heimessetmessages: + msg59572
2008-01-08 22:54:57mark.dickinsonsetmessages: + msg59570
2008-01-08 22:42:17christian.heimessetnosy: + christian.heimes
messages: + msg59568
2008-01-08 20:55:48mark.dickinsonsetmessages: + msg59557
2008-01-05 00:25:28mark.dickinsonsetfiles: + cmath3.patch
messages: + msg59268
2007-11-26 02:18:52mark.dickinsonsetfiles: + cmath.patch
messages: + msg57835
2007-11-04 18:38:37loewissetmessages: + msg57104
2007-11-04 15:16:22mark.dickinsonsetfiles: + cmath_py.py
nosy: + mark.dickinson
messages: + msg57098
2007-11-04 05:46:53alanmcintyresetnosy: + alanmcintyre
messages: + msg57092
2007-11-03 19:42:43inducersetmessages: + msg57090
2007-11-03 19:17:37loewissetnosy: + loewis
messages: + msg57089
2007-11-03 18:58:58inducercreate