-
-
Notifications
You must be signed in to change notification settings - Fork 31.6k
cmath is numerically unsound #45722
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Comments
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 Tim Peters remarks in [1] that basically all of cmath is unsound. I just wanted to make sure that this issue remains on the radar. |
Can you propose a patch? |
On Samstag 03 November 2007, Martin v. Löwis wrote:
Other than point at how boost.math does things, I don't have the time to work Andreas |
I have to review a few complex math topics for some of my current course |
I took a look at this a while back, and got as far as writing a pure As Tim Peters pointed out, the real problem here is a lack of decent |
It would be ok if a test is only run on a system with IEEE floats, and |
Here is (quite a large) patch, cmath.patch, that fixes a variety of
The patch includes documentation updates, but there are still no tests. |
Here's an updated patch for cmath, complete with tests and documentation changes. Any feedback would be much appreciated. |
A note to any potential reviewers for this patch (cmath3.patch): I'm especially
And it would be really great if someone with access to a Windows box could just |
#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
|
Okay: would it be okay for me to move the #ifdef MS_WINDOWS sequence There are replacements for log1p and asinh in the patch, so it shouldn't |
It may even be a good idea to move asinh, copysign and log1p to a new This is beyond a simple fix for the cmath module and should be discussed |
Guido, how do you like the idea of Include/pymath.h and Python/pymath.c |
Are you crazy? Adding new features to 2.5? No way! |
See bpo-1640 and svn+ssh://pythondev@svn.python.org/python/branches/trunk-math |
Substantial fixes for the cmath module went into the trunk and the py3k |
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: >>> import cmath
>>> cmath.asinh(-2j)
(1.3169578969248166-1.5707963267948966j) cheers, |
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? |
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. |
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: cheers, |
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, |
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. :) |
ah! |
Note: these values reflect the state of the issue at the time it was migrated and might not reflect the current state.
Show more details
GitHub fields:
bugs.python.org fields:
The text was updated successfully, but these errors were encountered: