[issue1381] cmath is numerically unsound
Sebastien Binet bi...@cern.ch added the comment: 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. -- nosy: +bins versions: +Python 3.1 ___ Python tracker rep...@bugs.python.org http://bugs.python.org/issue1381 ___ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Mark Dickinson dicki...@gmail.com added the comment: 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.9998j) sinh(-1.3169578969248164-1.5707963267948966j) (-1.0605752387249064e-16-1.9993j) Perhaps numpy is using different branch cuts? -- ___ Python tracker rep...@bugs.python.org http://bugs.python.org/issue1381 ___ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Mark Dickinson dicki...@gmail.com added the comment: 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. -- ___ Python tracker rep...@bugs.python.org http://bugs.python.org/issue1381 ___ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Sebastien Binet bi...@cern.ch added the comment:
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.
--
___
Python tracker rep...@bugs.python.org
http://bugs.python.org/issue1381
___
___
Python-bugs-list mailing list
Unsubscribe:
http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Mark Dickinson dicki...@gmail.com added the comment: 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). -- ___ Python tracker rep...@bugs.python.org http://bugs.python.org/issue1381 ___ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Mark Dickinson dicki...@gmail.com added the comment: 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. :) -- ___ Python tracker rep...@bugs.python.org http://bugs.python.org/issue1381 ___ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Sebastien Binet bi...@cern.ch added the comment: Note that in your program, you're feeding complex(-0.0, -2.0) to asinh, not complex(0.0, -2.0). ah! (ducking) -- ___ Python tracker rep...@bugs.python.org http://bugs.python.org/issue1381 ___ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Changes by Jesús Cea Avión [EMAIL PROTECTED]: -- nosy: +jcea ___ Python tracker [EMAIL PROTECTED] http://bugs.python.org/issue1381 ___ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Mark Dickinson [EMAIL PROTECTED] added the comment: 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). -- resolution: - fixed status: open - closed __ Tracker [EMAIL PROTECTED] http://bugs.python.org/issue1381 __ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Christian Heimes added the comment: See #1640 and svn+ssh://[EMAIL PROTECTED]/python/branches/trunk-math -- components: +Extension Modules -Library (Lib) keywords: +patch priority: - normal versions: +Python 3.0 -Python 2.5 __ Tracker [EMAIL PROTECTED] http://bugs.python.org/issue1381 __ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Christian Heimes added the comment: 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? -- nosy: +gvanrossum versions: +Python 2.6 __ Tracker [EMAIL PROTECTED] http://bugs.python.org/issue1381 __ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Guido van Rossum added the comment: Are you crazy? Adding new features to 2.5? No way! __ Tracker [EMAIL PROTECTED] http://bugs.python.org/issue1381 __ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Mark Dickinson added the comment: 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. __ Tracker [EMAIL PROTECTED] http://bugs.python.org/issue1381 __ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Mark Dickinson added the comment: 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. __ Tracker [EMAIL PROTECTED] http://bugs.python.org/issue1381 __ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Christian Heimes added the comment: * 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. -- nosy: +tiran __ Tracker [EMAIL PROTECTED] http://bugs.python.org/issue1381 __ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Mark Dickinson added the comment: 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. Added file: http://bugs.python.org/file9065/cmath3.patch __ Tracker [EMAIL PROTECTED] http://bugs.python.org/issue1381 __ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Mark Dickinson added the comment: 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?) Added file: http://bugs.python.org/file8807/cmath.patch __ Tracker [EMAIL PROTECTED] http://bugs.python.org/issue1381 __Index: configure === --- configure (revision 59184) +++ configure (working copy) @@ -1,5 +1,5 @@ #! /bin/sh -# From configure.in Revision: 58653 . +# From configure.in Revision: 58784 . # Guess values for system-dependent variables and create Makefiles. # Generated by GNU Autoconf 2.61 for python 2.6. # @@ -15226,11 +15226,14 @@ -for ac_func in alarm bind_textdomain_codeset chflags chown clock confstr \ - ctermid execv fork fpathconf ftime ftruncate \ + + + +for ac_func in alarm asinh bind_textdomain_codeset chflags chown clock \ + confstr copysign ctermid execv fork fpathconf ftime ftruncate \ gai_strerror getgroups getlogin getloadavg getpeername getpgid getpid \ getpriority getpwent getspnam getspent getsid getwd \ - kill killpg lchflags lchown lstat mkfifo mknod mktime \ + kill killpg lchflags lchown log1p lstat mkfifo mknod mktime \ mremap nice pathconf pause plock poll pthread_init \ putenv readlink realpath \ select setegid seteuid setgid \ Index: configure.in === --- configure.in(revision 59184) +++ configure.in(working copy) @@ -2303,11 +2303,11 @@ AC_MSG_RESULT(MACHDEP_OBJS) # checks for library functions -AC_CHECK_FUNCS(alarm bind_textdomain_codeset chflags chown clock confstr \ - ctermid execv fork fpathconf ftime ftruncate \ +AC_CHECK_FUNCS(alarm asinh bind_textdomain_codeset chflags chown clock \ + confstr copysign ctermid execv fork fpathconf ftime ftruncate \ gai_strerror getgroups getlogin getloadavg getpeername getpgid getpid \ getpriority getpwent getspnam getspent getsid getwd \ - kill killpg lchflags lchown lstat mkfifo mknod mktime \ + kill killpg lchflags lchown log1p lstat mkfifo mknod mktime \ mremap nice pathconf pause plock poll pthread_init \ putenv readlink realpath \ select setegid seteuid setgid \ Index: Doc/library/cmath.rst === --- Doc/library/cmath.rst (revision 59184) +++ Doc/library/cmath.rst (working copy) @@ -14,6 +14,15 @@ floating-point number, respectively, and the function is then applied to the result of the conversion. +.. note:: + + On platforms with hardware and system-level support for signed + zeros, functions involving branch cuts are continuous on *both* + sides of the branch cut: the sign of the zero distinguishes one + side of the branch cut from the other. On platforms that do not + support signed zeros the continuity is as specified below. + + The functions are: @@ -37,32 +46,37 @@ .. function:: asinh(x) - Return the hyperbolic arc sine of *x*. There are two branch cuts, extending - left from ``±1j`` to ``±âj``, both continuous from above. These branch cuts - should be considered a bug to be corrected in a future release. The correct - branch cuts should extend along the imaginary axis, one from ``1j`` up to - ``âj`` and continuous from the right, and one from ``-1j`` down to ``-âj`` - and continuous from the left. + Return the hyperbolic arc sine of *x*. 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. + ..
[issue1381] cmath is numerically unsound
Mark Dickinson added the comment: 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. -- nosy: +marketdickinson Added file: http://bugs.python.org/file8685/cmath_py.py __ Tracker [EMAIL PROTECTED] http://bugs.python.org/issue1381 __ This module provides exponential, logarithmic, trigonometric, inverse trigonometric, hyperbolic and inverse hyperbolic functions for complex numbers. It is intended as a drop-in replacement for the cmath module in the standard library. For the most part, it uses formulas and methods described by W. Kahan in his `Branch cuts for elementary functions' paper. Design goals - make all functions numerically sound; both the real part and the imaginary part should be within a few ulps of the true result, where this is reasonable. (But see the note on accuracy below.) - avoid unnecessary overflows in intermediate expressions, when the final result is representable. - fix buggy branch cuts in asinh - do the 'right thing' with respect to signed zeros on platforms that follow IEEE754 and C99/IEC 60559. In particular, all branch cuts should be continuous from both sides in this case. - don't do anything unreasonable on platforms that don't support signed zeros, or have signed zero support that doesn't comply fully with the above standards. With no signed zeros, continuity at branch cuts should match the documentation. Behaviour here is untested. Non-design goals - do the right thing with NaNs and infinities. It's hard to do this portably. I believe that many of the routines below do actually do the right thing in the presence of NaNs, but this is mostly accidental. - (related to the above): give sensible and consistent exceptions. Again, it seems difficult to do this across platforms. - produce results that are provably accurate to within 1ulp. Note on accuracy In an ideal world, the complex-valued function f(z) would return the closest representable complex number to the true mathematical value of f(z): that is, both the real and imaginary part of the result would be accurate to within = 0.5ulp. Achieving this level of accuracy is very hard---most C math libraries don't manage it. (But see the crlibm and MPFR libraries.) A slightly more realistic goal is 1ulp. In practice, one might hope that the returned real and imaginary parts are always within a few ulps (say 10 or 20) of those of the closest representable value. But even this is an unrealistic goal in some situtations. For example let z be the complex number: 0.79993338661852249060757458209991455078125 + 0.600088817841970012523233890533447265625j which is exactly representable, assuming IEEE doubles are being used for the real and imaginary parts. The nearest representable complex number to the natural logarithm of z is: 6.16297582203915472977912941627176741932192527428924222476780414581298828125e-33 + 2.4980915447965088560522417537868022918701171875j It would take a lot of effort to get the real part anywhere near correct here. Other problems occur in computing trigonometric functions for complex numbers z with large real part, or hyperbolic trig functions for z with large imaginary part. Notes on signed zeros - There are many subtle difficulties here: for example, the expression 1 + z should add 1 to the real part of z and leave the imaginary part untouched. But in Python, if the imaginary part of z is -0. then the imaginary part of 1+z will be +0: 1 gets coerced to the complex number 1 + 0j, and then the imaginary parts get added to give -0j + 0j = 0j. But z - 1 always does the right thing: subtracting +0. won't change the sign of zero. Similarly, z + (-1.) is fine. So we can either work with the underlying reals directly, or rewrite the erroneous 1+z as -(-z-1), which works. Similarly, 1-z should be rewritten as -(z-1). An alternative fix is to use the complex number 1 - 0j (note negative sign) in place of 1 in the expressions 1+z and 1-z. Similarly, the expression i*z is special-cased so that i*(x+i*y) = -y + i*x; see the function mul_by_j. The code below should `do the right thing' regardless of whether signed zeros are present. In particular: - on a platform (hardware + C math library) that supports signed zeros, so that for example: atan2(-0., positive) gives -0. atan2(0., positive) gives 0. atan2(-0.,
[issue1381] cmath is numerically unsound
Martin v. Löwis added the comment: 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. __ Tracker [EMAIL PROTECTED] http://bugs.python.org/issue1381 __ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
New submission from Andreas Kloeckner: 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. -- components: Library (Lib) messages: 57088 nosy: inducer severity: normal status: open title: cmath is numerically unsound type: behavior versions: Python 2.5 __ Tracker [EMAIL PROTECTED] http://bugs.python.org/issue1381 __ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Martin v. Löwis added the comment: Can you propose a patch? -- nosy: +loewis __ Tracker [EMAIL PROTECTED] http://bugs.python.org/issue1381 __ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Andreas Kloeckner added the comment: 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 __ Tracker [EMAIL PROTECTED] http://bugs.python.org/issue1381 __ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
[issue1381] cmath is numerically unsound
Alan McIntyre added the comment: 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. -- nosy: +alanmcintyre __ Tracker [EMAIL PROTECTED] http://bugs.python.org/issue1381 __ ___ Python-bugs-list mailing list Unsubscribe: http://mail.python.org/mailman/options/python-bugs-list/archive%40mail-archive.com
