Issue1640

Created on **2007-12-17 15:48** by **christian.heimes**, last changed **2008-05-01 21:18** by **christian.heimes**. This issue is now **closed**.

Files | ||||
---|---|---|---|---|

File name | Uploaded | Description | Edit | |

trunk_float_math_combined.patch | christian.heimes, 2007-12-17 20:17 | |||

trunk_float_math_combined2.patch | christian.heimes, 2007-12-18 10:06 | |||

trunk_pymath_hyberbolic_complex2.patch | christian.heimes, 2008-01-21 18:04 | |||

invhyp.c | mark.dickinson, 2008-01-21 20:42 |

Messages (56) | |||
---|---|---|---|

msg58694 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2007-12-17 15:48 | |

The patch adds several small enhancements to the math module and pyport.h. * Py_MATH_PI and Py_MATH_E in long double precision * Py_IS_NAN and Py_IS_INFINITY use isnan() and isinf() functions were available (checked by configure) * isnan and isinf for the math module * Bessel (1st and 2nd kind), erf, erfc, lgamma function for math module. They are defined in almost (or all?) math libraries. Together with http://bugs.python.org/issue1635 it implements most of PEP 754. |
|||

msg58698 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2007-12-17 17:40 | |

Cool! If there's a move to add functions to the math module, there are some others that are part of C99 (but not C89), would be good to have, and that I'd consider more fundamental than the Bessel, error, gamma functions; for example, the inverse hyperbolic trig functions (acosh, asinh, atanh), log1p, expm1, copysign. |
|||

msg58699 - (view) | Author: Adam Olsen (Rhamphoryncus) | Date: 2007-12-17 17:53 | |

Minor typo. Should be IEEE: "Return the sign of an int, long or float. On platforms with full IEE 754\n\" |
|||

msg58702 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2007-12-17 19:27 | |

Mark Dickinson wrote: > Mark Dickinson added the comment: > > Cool! If there's a move to add functions to the math module, there are > some others that are part of C99 (but not C89), would be good to have, and > that I'd consider more fundamental than the Bessel, error, gamma > functions; for example, the inverse hyperbolic trig functions (acosh, > asinh, atanh), log1p, expm1, copysign. I've added the inverse hyperbolic, log1p and expm1. copysign is too low level but I've added sign(x) -> -1/0/+1. It uses copysign() where available so you can emulate copysign(x, y) with sign(x) * y. Do you want some more functions? http://www.dinkumware.com/manuals/?manual=compleat&page=math.html Feel free to provide a patch :) Christian |
|||

msg58706 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2007-12-17 20:17 | |

Combined patch for #1635 and #1640. UPDATES: * Fixed Py_NAN macro * Fixed typo * Added inverse hyperbolic trigonometric functions * Added expm1 and log1p Missing: * doc updates * some real unit tests for the a*h and other new functions |
|||

msg58735 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2007-12-18 10:06 | |

The patch fixes some compatibility problems on Windows and implements several methods like acosh in C when the system's math lib doesn't provide the functions. |
|||

msg58749 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2007-12-18 14:36 | |

Unfortunately, implementing asinh, acosh, atanh correctly is not as simple as just using the formulas (this is one of the reasons that it's useful to be able to get at the library definitions). For example, the formula asinh(x) = log(x + sqrt(1.+(x*x))) has significant error (i.e. much more than just a few ulps) when x is small (so that the argument of the log is close to 1), and when x is large and negative (so that the addition involves adding two numbers of almost equal magnitude but opposite sign). It also overflows unnecessarily at around 1e154 (on an IEEE754 system), as a result of the intermediate calculation x*x overflowing, and it fails to give the correct signs on zero arguments. (asinh(0.) = 0., asinh(-0.) = -0.) So either some serious work is required here, or code should be borrowed in an appropriately legal fashion from some other library, or those few people who don't have asinh, acosh, etc. already in libm will have to live without. Mark |
|||

msg58758 - (view) | Author: Guido van Rossum (gvanrossum) * | Date: 2007-12-18 19:36 | |

i suggest abandoning any attempts at implementing math ourselves. I also suggest not combining this with #1635 but reviewing and (eventually) applying the latter first. |
|||

msg58770 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2007-12-18 21:05 | |

Guido van Rossum wrote: > Guido van Rossum added the comment: > > i suggest abandoning any attempts at implementing math ourselves. I > also suggest not combining this with #1635 but reviewing and > (eventually) applying the latter first. How do you like a compromise? I rip out the new functions for the math module except isnan(), isinf() and sign(). The new patch would contain the inf and nan roundtrip plus the helper methods. Mark is right. My naive functions are a numerical nightmare for small and large numbers. I copied them from my math handbook 'cause I couldn't find solid implementations in my numerical handbook. It contains all sort of algorithms for approximations, ODE, PDE, FFT but no sample of the inverse hyperbolic functions. We could reuse the functions from the BSD libm. They are used in numerous C runtime libraries (BSD, Mac OS X, uclibc). http://packages.e.kth.se/common/src/os/NetBSD/1.3.2/lib/libm/src/ Christian |
|||

msg58786 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2007-12-18 23:28 | |

The trunk_math_sign_inf_nan patch contains just three new method isnan(), isinf() and sign(). It also fixes a minor numerical issue with one function that did small / (small / large) instead of small * (large / small). |
|||

msg59115 - (view) | Author: Guido van Rossum (gvanrossum) * | Date: 2008-01-03 00:11 | |

One nit: you added a blank line to the end of test_math.py. This will cause the checkin to fail. :-) One bigger issue: the sign() function doesn't seem to work properly for nans. E.g. on Linux I get: >>> inf = 1e1000 >>> nan = inf/inf >>> mnan = -nan >>> math.sign(nan) -1 >>> math.sign(mnan) 1 IOW a positive nan is considered negative and vice versa. (This is probably due to the way nans defy testing, always returning false.) I'm also curious why math.sign(0.0) returns 1 -- this is going to cause a lot of confusion. This is also missing doc patches. |
|||

msg59137 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-01-03 16:57 | |

Guido van Rossum wrote: > Guido van Rossum added the comment: > > One nit: you added a blank line to the end of test_math.py. > This will cause the checkin to fail. :-) *grr* stupid white space check hook > One bigger issue: the sign() function doesn't seem to work properly for > nans. E.g. on Linux I get: > >>>> inf = 1e1000 >>>> nan = inf/inf >>>> mnan = -nan >>>> math.sign(nan) > -1 >>>> math.sign(mnan) > 1 > > IOW a positive nan is considered negative and vice versa. (This is > probably due to the way nans defy testing, always returning false.) If I recall the definition correctly NaNs don't have a sign. The content of the sign bit is not defined for NaNs. I could fix the sign but it's just eye candy and a waste of CPU cycles. IMO it would be more appropriate to return 0 for NaNs instead of +1 or -1. > I'm also curious why math.sign(0.0) returns 1 -- this is going to > cause a lot of confusion. math.sign(0.0) == 1 and math.sign(-0.0) == -1 is the main purpose of the sign() function. Most ordinary users are still going to use x > 0.0 or x < 0.0 instead of math.sign(). math.sign() is for the power users who need to determinate whether an operation hits the 0 from the left or right side of the number line. Because 0.0 == -0.0 it's not possible the distinguish a positive from a negative zero with comparison operations. if x == 0.0 and str(x).startswith("-") was the only existing way to detect a negative zero. Christian |
|||

msg59138 - (view) | Author: Guido van Rossum (gvanrossum) * | Date: 2008-01-03 17:28 | |

> > One nit: you added a blank line to the end of test_math.py. > > This will cause the checkin to fail. :-) > > *grr* stupid white space check hook No, you edited a line that didn't need editing. :-) > > One bigger issue: the sign() function doesn't seem to work properly for > > nans. E.g. on Linux I get: > > > >>>> inf = 1e1000 > >>>> nan = inf/inf > >>>> mnan = -nan > >>>> math.sign(nan) > > -1 > >>>> math.sign(mnan) > > 1 > > > > IOW a positive nan is considered negative and vice versa. (This is > > probably due to the way nans defy testing, always returning false.) > > If I recall the definition correctly NaNs don't have a sign. The content > of the sign bit is not defined for NaNs. I could fix the sign but it's > just eye candy and a waste of CPU cycles. IMO it would be more > appropriate to return 0 for NaNs instead of +1 or -1. Perhaps you recall wrong, as negating the nan returns one with the opposite sign? This seems to indicate that there are positive and negative nans. > > I'm also curious why math.sign(0.0) returns 1 -- this is going to > > cause a lot of confusion. > > math.sign(0.0) == 1 and math.sign(-0.0) == -1 is the main purpose of the > sign() function. Most ordinary users are still going to use x > 0.0 or x > < 0.0 instead of math.sign(). math.sign() is for the power users who > need to determinate whether an operation hits the 0 from the left or > right side of the number line. > > Because 0.0 == -0.0 it's not possible the distinguish a positive from a > negative zero with comparison operations. if x == 0.0 and > str(x).startswith("-") was the only existing way to detect a negative zero. Hm, OK, but then passing a zero of some other type (e.g. int) should also return +1 as the sign. I also think the function's name should be changed, because I (and I assume many others) have grown up with a sign() function that essentially returns cmp(x, 0.0). Perhaps it would be better to have a function math.isneg() that returns True for -0.0 and anything smaller and False for +0.0 and anything larger. It could also return the proper sign of a nan. I suggest that you check in the isinf() and isnan() functions and their tests once you have docs for them, and hold off on the sign() function until we've got agreement on it. |
|||

msg59152 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-01-03 20:21 | |

Guido van Rossum wrote: > Hm, OK, but then passing a zero of some other type (e.g. int) should > also return +1 as the sign. I also think the function's name should be > changed, because I (and I assume many others) have grown up with a > sign() function that essentially returns cmp(x, 0.0). > > Perhaps it would be better to have a function math.isneg() that > returns True for -0.0 and anything smaller and False for +0.0 and > anything larger. It could also return the proper sign of a nan. I'm fine with a isneg() function but I wouldn't "fix" it for NaN. It has probably some kind of obscure meaning. The best explanation I was able to find, is http://www.cisl.ucar.edu/docs/trap.error/errortypes.html "Note: Since NaN is "not a number," there really isn't a "+" or "-" sign associated with it." Christian |
|||

msg59153 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-03 20:29 | |

Why not implement copysign? It's a standard, familiar function with well- documented meaning all over the web, and it can easily be used to create whatever sign test is necessary, letting the user decide what results (s)he wants for +/-0, +/-nan, etc. |
|||

msg59154 - (view) | Author: Guido van Rossum (gvanrossum) * | Date: 2008-01-03 20:35 | |

> Why not implement copysign? It's a standard, familiar function with well- > documented meaning all over the web, and it can easily be used to create > whatever sign test is necessary, letting the user decide what results > (s)he wants for +/-0, +/-nan, etc. Good idea. Since you seem to like providing patches, can you create one for math.copysign()? |
|||

msg59155 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-01-03 20:42 | |

Guido van Rossum wrote: > Good idea. Since you seem to like providing patches, can you create > one for math.copysign()? Don't forget it's copysign() on Unix but _copysign() on Windows. #if defined(MS_WINDOWS) || defined(HAVE_COPYSIGN) #ifdef MS_WINDOWS FUNC2(copysign, _copysign, #else FUNC2(copysign, copysign, #endif "doc"); #endif should work on all systems. Christian |
|||

msg59156 - (view) | Author: Guido van Rossum (gvanrossum) * | Date: 2008-01-03 20:50 | |

> Guido van Rossum wrote: > > Good idea. Since you seem to like providing patches, can you create > > one for math.copysign()? > > Don't forget it's copysign() on Unix but _copysign() on Windows. > > #if defined(MS_WINDOWS) || defined(HAVE_COPYSIGN) > #ifdef MS_WINDOWS > FUNC2(copysign, _copysign, > #else > FUNC2(copysign, copysign, > #endif > "doc"); > #endif > > should work on all systems. Well, the Python API in the math module should always be called copysign(). :-) And what to do if neither is present? Are there any systems without it? |
|||

msg59165 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-01-03 22:03 | |

Guido van Rossum wrote: > Well, the Python API in the math module should always be called copysign(). :-) > > And what to do if neither is present? Are there any systems without it? takes care of it. It's a macro to define a function which accepts two floats and returns a float: FUNC2(funcname, func, docstring). On Windows _copysign is always defined and on other systems we can use HAVE_COPYSIGN. I added it to configure.in a while ago. Christian |
|||

msg59166 - (view) | Author: Guido van Rossum (gvanrossum) * | Date: 2008-01-03 22:04 | |

OK, just check it in then, but do add docs! On Jan 3, 2008 2:03 PM, Christian Heimes <report@bugs.python.org> wrote: > > Christian Heimes added the comment: > > Guido van Rossum wrote: > > Well, the Python API in the math module should always be called copysign(). :-) > > > > And what to do if neither is present? Are there any systems without it? > > takes care of it. It's a macro to define a function which accepts two > floats and returns a float: FUNC2(funcname, func, docstring). > > On Windows _copysign is always defined and on other systems we can use > HAVE_COPYSIGN. I added it to configure.in a while ago. > > Christian > > > __________________________________ > Tracker <report@bugs.python.org> > <http://bugs.python.org/issue1640> > __________________________________ > |
|||

msg59211 - (view) | Author: Tim Peters (tim.peters) * | Date: 2008-01-04 03:51 | |

The functionality of what was called (and I agree confusingly so) "sign()" here is supplied by "signbit()" in C99. That standard isn't free, but the relevant part is incorporated in the free Open Group standards: http://www.opengroup.org/onlinepubs/000095399/functions/signbit.html The 754 standard says NaNs have sign bits, but assigns no meaning to them. In particular, the value of the sign bit of a NaN resulting from a string->double routine applied to the string "nan" isn't defined by 754, or (AFAICT) by C99 either. |
|||

msg60257 - (view) | Author: George Castillo (gmcastil) | Date: 2008-01-20 02:12 | |

Is there still interest in implementing the inverse hyperbolic trig functions for real numbers? I would be willing to explore this if there is. |
|||

msg60259 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-20 03:07 | |

George: I'm certainly still interested in having asinh, acosh and atanh in math---I'm not sure about anyone else, but I consider these three functions to be basic ingredients in any math library. They even appear in most calculus texts. And the complex versions are already in cmath. The right thing to do would be to wrap the libm functions when they exist, but provide fallbacks for the platforms where they don't. I think these functions are present on OS X and in the GNU libm, but they might be missing on Windows. There's a fallback version of asinh already provided in the patch to issue #1381, which I believe to be reasonably sound (but since I wrote it, a second opinion would be good). For acosh and atanh, you'd have to find the right way to deal with domain errors (e.g. acosh(0.5), atanh(2.0)) and singularities (atanh(1.0), atanh(-1.0)). I'd suggest trying to follow the details in Appendix F (section F.9.2) of the C99 standard where possible. |
|||

msg60264 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-01-20 08:57 | |

I'm +1 in adding fallbacks for important functions like copysign, asinh, acosh and atanh. expm1 and log1p may be worth adding, too. Windows doesn't have any of the functions except of _copysign(). But why write our own version if we can reuse existing implementations? I found a set of very well written and documented implementations in the uclibc sources of libm. The sources are under a BSDish license: * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== In #1381 I suggested two new files Python/pymath.c and Include/pymath.h. We could stick all the replacement implementations in the files and maybe move some of the code from Include/pyport.h and Python/hypot.c to the new files. |
|||

msg61290 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-01-20 13:56 | |

The patch implements replacements for copysign, log1p (taken from Mark's patch #1381) asinh, acosh, atanh (taken from uclibc's libm). I modified the a*h functions for our needs. They set an errno and use our macros. The patch also adds the three hyberbolic arc functions to the math module. Math related code is now in Include/pymath.h and Python/pymath.c. copysign is now defined on Windows, too. I also found a bug in configure.in (my fault). The C99 function is called finite(), not isfinite(). Christian |
|||

msg61306 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-20 15:17 | |

Excellent! I'll take a look. |
|||

msg61308 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-20 16:11 | |

One question: is there a policy on what should happen for singularities and domain errors? If not, I think it would be useful to have one. Following the policy may not be achievable on all platforms, but it should be easy to do on major platforms, and at least we'll know what we're aiming for in general. Maybe it already exists, and I missed it :) For domain errors (e.g. sqrt(-1), log(-1), acosh(0)), the obvious two options are: - raise an exception (e.g. OverflowError), or - return a NaN For singularities (e.g. log(0), atanh(1)), the options are basically the same: - raise an exception (preferably something different from OverflowError), or - return the IEEE-754 recommended value (usually +/-Inf) I suspect there are use-cases for both types of behaviour here. Of course, the *right* thing to do, in some sense, would be to have a thread-local floating-point environment that makes it possible for the user to choose whether he/she wants an exception or a special value, much like the way Decimal behaves at the moment. But that would be a big change, almost certainly requiring a PEP and a lot of work. A few months ago I definitely would have said that an exception should be raised in both cases, as already happens (mostly); but since then NaNs and Infinities have acquired greater legitimacy within Python. Tim, if you're listening: any words of wisdom? Should I ask this on python-dev? |
|||

msg61309 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-20 16:33 | |

whoops. OverflowError should be something else in the previous post; of course, OverflowError is inappropriate for domain errors (but entirely appropriate for something like exp(1000)). Currently, on Linux I get: - overflow (exp(1000)) -> OverflowError - domain error (sqrt(-1)) -> ValueError - singularity (log(0)) -> OverflowError On OS X I get: - overflow -> OverflowError - domain error -> NaN - singularity -> OverflowError |
|||

msg61310 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-01-20 16:53 | |

Mark Dickinson wrote: > Currently, on Linux I get: > - overflow (exp(1000)) -> OverflowError > - domain error (sqrt(-1)) -> ValueError > - singularity (log(0)) -> OverflowError Windows raises the same exceptions as Linux. |
|||

msg61311 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-20 17:09 | |

Okay: for now, I guess we just follow the pattern that already exists on Linux and Windows. I think the OS X sqrt(-1) behaviour is a bug. |
|||

msg61366 - (view) | Author: George Castillo (gmcastil) | Date: 2008-01-20 22:24 | |

Just a quick addition here regarding the singularities to these functions. The atanh(x) is only defined for |x| < 1, so atanh(1) or atanh(-1) isn't singular there so much as simply isn't defined. So, even though the function approaches infinite as x -> 1, it wouldn't really be correct to return a value at |x| = 1. I think raising an exception at those points would be more correct. |
|||

msg61367 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-20 22:54 | |

No: IEEE-754r and the C99 standard both say clearly that atanh(1) should be infinity and atanh(-1) should be -infinity, and furthermore that the 'divide-by-zero' exception should be raised rather than the 'invalid' exception. It's a singularity, just like log(0). (This makes even more sense viewed from the perspective of complex arithmetic, where atanh is defined at all points in the complex plane except -1 and 1, where it has log-type singularities.) The general idea is that it's meaningful to set atanh(1) = infinity because that's what the limit of atanh(x) is as x approaches 1 from below; similarly for atanh(-1) and log(0). |
|||

msg61370 - (view) | Author: George Castillo (gmcastil) | Date: 2008-01-21 00:30 | |

I misunderstood the rationale for the function returning infinite at those points - I didn't realize that C99 was the governing force behind the implementation of these functions, rather than mathematical rigor. Thanks for pointing it out. In that case, I agree with you that, in order to conform, these functions would need to return the values required by those documents. |
|||

msg61371 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-21 01:02 | |

George: I think my last post was a bit rude. I apologize if it came across that way. Mathematical rigor and IEEE-754 recommendations aren't necessarily in conflict here, though. For example, the natural log function from (0, infinity) to (-infinity, infinity) extends naturally and uniquely to a continuous function on the closed subset [0, infinity] of the extended real line---i.e., the real line together with the two extra points -infinity and infinity. With the appropriate topology, the extended real line is a perfectly well-defined and well-behaved mathematical object, though of course it's no longer a field. Since IEEE-754 floats include infinities, it's reasonable, and sometimes useful, to regard the set of IEEE-floats as a computational model of the extended real line instead of the reals. At any rate, I agree with you that log(0) and atanh(1) should raise Python exceptions, at least for now. But these calculations are qualitatively different from log(-1) and atanh(2), and it wouldn't be at all unreasonable if they raised a different exception--- e.g. ZeroDivisionError instead of ValueError. |
|||

msg61372 - (view) | Author: Tim Peters (tim.peters) * | Date: 2008-01-21 01:05 | |

Mark, these are the "spirit of 754" rules: 1. If the mathematical result is a real number, but of magnitude too large to approximate by a machine float, overflow is signaled and the result is an infinity (with the appropriate sign). 2. If the mathematical result is a real number, but of magnitude too small to approximate by a machine float, underflow is signaled and the result is a zero (with the appropriate sign). 3. At a singularity (a value x such that the limit of f(y) as y approaches x exists and is an infinity), "divide by zero" is signaled and the result is an infinity (with the appropriate sign). This is complicated a little by that the left-side and right-side limits may not be the same; e.g., 1/x approaches +inf or -inf as x approaches 0 from the positive or negative directions. In that specific case, the sign of the zero determines the result of 1/0. 4. At a point where a function has no defined result in the extended reals (i.e., the reals plus an infinity or two), invalid operation is signaled and a NaN is returned. And these are what Python has historically /tried/ to do (but not always successfully, as platform libm behavior varies a lot): For #1, raise OverflowError. For #2, return a zero (with the appropriate sign if that happens by accident ;-)). For #3 and #4, raise ValueError. It may have made sense to raise Python's ZeroDivisionError in #3, but historically that's only been raised for division by zero and mod by zero. You're right that you can't make everyone happy, so may as well stick with historical consistency, and wait for a PEP to introduce a more flexible mechanism. As W. Kahan wrote long ago, the reason they're called "exceptions" is that no matter which fixed policy you adopt, someone will take strident exception to it ;-) |
|||

msg61376 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-21 03:34 | |

Thanks, Tim! Dare I suggest extending these rules to encompass things like sqrt(NaN), log(inf), etc., as follows: - return a special value in Python where IEEE-754r/C99 specifies a special value, but doesn't raise any of the three divide-by-zero, invalid, or overflow exceptions, and - raise OverflowError or ValueError as appropriate where IEEE-754r specifies raising one of these exceptions. So e.g. cos(infinity) should give a ValueError, while log(infinity) and exp(infinity) should not raise any Python exception, but should return an infinity instead. And most single variable operations should return an input NaN unaltered, without raising an exception. |
|||

msg61385 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-01-21 12:24 | |

Mark Dickinson wrote: > So e.g. cos(infinity) should give a ValueError, while log(infinity) and exp(infinity) > should not raise any Python exception, but should return an infinity instead. And most > single variable operations should return an input NaN unaltered, without raising an > exception. The matter should be discussed in a proper PEP and targeted for Python 3.0. Python 3.0 is the right place for subtle changes which may break code. For Python 2.6 we must not change the exception or outcome of a function and new functions should be as consistent with existing ones as possible. I still don't like the idea of math.atanh(1) == inf. Why? See for yourself: 18.714973875118524 >>> math.atanh(.99999999999999999) Traceback (most recent call last): File "<stdin>", line 1, in <module> ValueError: math domain error (Linux) Christian |
|||

msg61386 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-01-21 12:31 | |

The mail interface screwed up the message: >>> math.atanh(.9999999999999999) 18.714973875118524 >>> math.atanh(.99999999999999999) Traceback (most recent call last): File "<stdin>", line 1, in <module> ValueError: math domain error |
|||

msg61394 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-01-21 13:52 | |

I've added your complex patch and its tests to my patch. The tests are showing some flaws in the atanh (or log1p) function on Windows: AssertionError: atanh0022:atanh(-1.000000) returned -18.36840028483855, expected -18.714973875118524 On Linux the tests are passing just fine. |
|||

msg61398 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-01-21 14:35 | |

It's most probably the fault of log1p(): AssertionError: atanh0022:atanh(-0.99999999999999989) returned -18.36840028483855, expected -18.714973875118524 |
|||

msg61399 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-21 14:53 | |

Christian: I'm definitely not proposing atanh(1) = inf: it should raise ValueError. I'm proposing that we follow Tim's rules for now; this means no change for finite inputs. The new thing here is that since you've made inf and nan more accessible and consistent across platforms, I think we should make sure that the math functions do the right thing for an *input* of +/-inf or nan. I'm almost sure that the current behavior of e.g. exp(float("inf")) is more-or-less accidental rather than designed. I think I'm missing the point of your math.atanh(.999...) example. .99999999999999999 *is* already exactly equal to 1.0, so you're just proving that math.atanh(1.0) currently gives a ValueError. (Which, again, I think is the right thing to do.) >>> x = .99999999999999999 >>> x == 1.0 True The atanh0022 result is definitely a bug: it looks like either asinh or log1p is buggy. I'll try to figure it out. |
|||

msg61401 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-21 15:11 | |

Christian: would it be possible for you to tell me what the argument of the log1p call is on Windows in that last branch of c_atanh, for the testcase atanh0022. On OS X I get: Input to log1p is 3.2451855365842669075e+32 It's hard to imagine that there's anything wrong with log1p here, since all it does for a large input like this is compute log(1+x). |
|||

msg61402 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-21 15:25 | |

Okay: here's an attempted guess at what's happening on Windows: Near the end of c_atanh, there's a line that reads: r.real = log1p(4.*z.real/((1-z.real)*(1-z.real) + ay*ay))/4.; Somehow, when z.real is 0.99999999999999989 and ay is 0, the argument to log1p is ending up one-quarter of the size that it should be. I suspect that the 1-z.real calculation is producing, for reasons that are beyond me, the float 2**-52 instead of the correct value of 2**- 53. Christian: if you have any time to play with this, could you try replacing this line with something like: double temp = 1-z.real printf("z.real is %.19e\n", z.real); r.real = log1p(4.*z.real/(temp*temp + ay*ay))/4.; and see if either the problem goes away or if you can confirm that temp is coming out to be 2.2204460492503131e-16 rather than 1.1102230246251565e-16. |
|||

msg61403 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-21 15:31 | |

Sorry: those lines should have been: double temp = 1-z.real; printf("temp is %.19e\n", temp); r.real = log1p(4.*z.real/(temp*temp + ay*ay))/4.; |
|||

msg61414 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-01-21 16:56 | |

Mark Dickinson wrote: > Mark Dickinson added the comment: > > Christian: would it be possible for you to tell me what the argument of the > log1p call is on Windows in that last branch of c_atanh, for the testcase > atanh0022. > > On OS X I get: > > Input to log1p is 3.2451855365842669075e+32 > > It's hard to imagine that there's anything wrong with log1p here, since all > it does for a large input like this is compute log(1+x). You got me wrong (and I didn't explain it properly). All complex functions pass the test. math.atanh() fails. I think my implementation of Python/pymath.c:atanh() doesn't return the right value for arguments almost 1.0. Christian |
|||

msg61419 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-21 17:14 | |

Sorry: I should have read more carefully. So math.atanh works on Linux but is producing some strange results on Windows. It's still rather puzzling though. I still suspect that it's the argument to log1p that's coming out wrong rather than the result. |
|||

msg61423 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-01-21 17:24 | |

Mark Dickinson wrote: > Sorry: I should have read more carefully. So math.atanh works on Linux > but is producing some strange results on Windows. > > It's still rather puzzling though. I still suspect that it's the > argument to log1p that's coming out wrong rather than the result. It uses t = 0.5 * log1p((x + x) / (1.0 - x)) for t > 0.5. I presume the culprit is in "2 / x" where x is almost 0. Do you have an idea how we can increase the accuracy for values nearly 1.? Christian |
|||

msg61431 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-01-21 18:04 | |

I disabled the tests for atanh0022 and atanh0023. Other changes: Added math.log1p + tests Added SUN license to Doc/licenses.rst Added docs to Doc/library/math.rst |
|||

msg61441 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-21 19:29 | |

The problem with atanh is that it should be using absx throughout, rather than x. So in "if (absx < 0.5)" branch and the following branch, replace all occurrences of x with absx, and it should work. One other comment: asinh shouldn't directly set errno for a NaN. It should do the same as acosh and atanh: return x+x. This makes asinh(float("nan")) return a nan, which makes it consistent with acosh and atanh, consistent with the way that Linux and OS X behave, and consistent with the other single-argument functions in math. I think asinh should also return x+x for x an infinity. This again should make it consistent with the way that the libm asinh works on OS X and Linux. |
|||

msg61442 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-21 19:39 | |

Also, for the C-level routines, atanh(1.0) and atanh(-1.0) should definitely return infinity and -infinity (and probably set errno as well.) Note that this is not an argument about what Python should do: Python will still raise a ValueError for atanh(1.0) and atanh(-1.0). But the atanh is supposed to be a drop-in replacement for the libm atanh, on those platforms where it's missing. And the C99 standard is clear about return values, even though it's not useful when it comes to deciding whether to set errno or not. |
|||

msg61443 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-21 20:13 | |

One more comment: the Py_IS_NAN test in acosh is likely never reached, since the comparison x >= two_pow_p28 will probably evaluate to 0 when x is a NaN. |
|||

msg61444 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-21 20:42 | |

Also: -1.0 shouldn't be returned at this level to indicate an error; these are pure C functions we're writing---the Python wrappers, and Python error conventions, apply further up the chain somewhere. Just so I'm doing something a little more constructive than yelling criticisms from the sideline, I've attached a file invhyp.c that's how I think the C functions should look. Some notes: - it doesn't touch errno, but lets the platform decide how to handle errors (i.e. produce a special value/set errno/signal a floating-point exception/some combination of these). This will make the asinh, acosh, atanh functions behave in the same way that the regular libm functions behave on any platform. So e.g. if a particular platform is used to setting errno for domain errors like sqrt(-1), it'll do so for atanh/asinh/acosh. And another platform that signals a floating-point exception for sqrt(-1) will do the same for atanh(3). - I've left in the "huge+x > 1.0" test in asinh; it's kinda pointless, but also fairly harmless. I think its only point is to make sure that the inexact flag gets set where appropriate, and Python doesn't much care about the inexact flag. - I'm reasonably sure that the test "x == 1." in acosh is safe. |
|||

msg61474 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-01-22 01:41 | |

Mark Dickinson wrote: > - it doesn't touch errno, but lets the platform decide how to handle errors (i.e. produce a > special value/set errno/signal a floating-point exception/some combination of these). This will > make the asinh, acosh, atanh functions behave in the same way that the regular libm functions > behave on any platform. So e.g. if a particular platform is used to setting errno for domain > errors like sqrt(-1), it'll do so for atanh/asinh/acosh. And another platform that signals a > floating-point exception for sqrt(-1) will do the same for atanh(3). I tried your patch. It fixes the problem with atanh0022 and 0023 test but it breaks in other places. math.acosh(0) returns NaN and does NOT raise an exception. math.atanh(1) raises OverflowError instead of ValueError. The libm of uclibc *does* set errno for signaling NaN. The relevant code is in w_acos.c and k_standard.c. I don't know how emit a signal for a NaN but setting errno sounds reasonable for me and it gives the desired output. w_acosh.c: double acosh(double x) /* wrapper acosh */ { #ifdef _IEEE_LIBM return __ieee754_acosh(x); #else double z; z = __ieee754_acosh(x); if(_LIB_VERSION == _IEEE_ || isnan(x)) return z; if(x<1.0) { return __kernel_standard(x,x,29); /* acosh(x<1) */ } else return z; #endif } kernel_standard: case 129: /* acosh(x<1) */ exc.type = DOMAIN; exc.name = type < 100 ? "acosh" : "acoshf"; exc.retval = zero/zero; if (_LIB_VERSION == _POSIX_) errno = EDOM; else if (!matherr(&exc)) { if (_LIB_VERSION == _SVID_) { (void) WRITE2("acosh: DOMAIN error\n", 20); } errno = EDOM; } break; |
|||

msg61479 - (view) | Author: Mark Dickinson (mark.dickinson) * | Date: 2008-01-22 02:10 | |

Hmmm. For atanh(1): raising OverflowError is actually consistent with what currently happens. The only singularity that's already present in math is log(0), and it seems that that raises OverflowError on OS X, Linux and Windows... I wonder whether this is what Tim meant to say? For acosh(0): you're right, and I'm wrong---this should definitely return a NaN and set errno. I guess that dividing 0 by 0 doesn't set errno on Windows. Okay: let's set it directly there. I do still think that asinh(nan), atanh(nan) and acosh(nan) should return nan and not raise any exceptions, just for the sake of consistency with Linux/OS X and with the other libm functions. I guess I don't really care about asinh(+/-inf), etc: an infinite return value will be caught by the stuff in math_1 anyway. |
|||

msg61498 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-01-22 11:10 | |

Mark Dickinson wrote: > For atanh(1): raising OverflowError is actually consistent with what currently happens. > The only singularity that's already present in math is log(0), and it seems that that > raises OverflowError on OS X, Linux and Windows... I wonder whether this is what Tim > meant to say? But on Linux atanh(1) sets errno = EDOM and raises a ValueError. Either the outcome is not consistent or the libm developers have a different philosophy than Time and you. *g* > For acosh(0): you're right, and I'm wrong---this should definitely return a NaN and set > errno. I guess that dividing 0 by 0 doesn't set errno on Windows. Okay: let's set it > directly there. I'm not an expert on IEEE754 and NaN but I guess "Signaling NaN" also means it sets errno. Personally I'm not interested in reimplementing or fixing Windows' mathematical library (although it's tempting *g*) but to get the results right for our purpose. You are right, we should try to mimic the C99 standard whenever it is possible. But it's more important to get the functions of our math module right. > I do still think that asinh(nan), atanh(nan) and acosh(nan) should return nan and not > raise any exceptions, just for the sake of consistency with Linux/OS X and with the other > libm functions. Yes, it makes perfectly sense and glibc's libm already does the right thing. However Windows screws up again so I've to fix the function call manually. Damn you Windows! I inserted the four lines in mathmodule.c:math_1() and something similar in math_2() to fix Windows. #ifndef __GNUC__ /* Windows et al */ if (Py_IS_NAN(x)) return PyFloat_FromDouble(x+x); #endif > I guess I don't really care about asinh(+/-inf), etc: an infinite return value will be > caught by the stuff in math_1 anyway. Are you fine with OverflowError? I'm going to create a branch for your work. It's getting tedious sending large patches back and forth. Christian |
|||

msg66047 - (view) | Author: Christian Heimes (christian.heimes) * | Date: 2008-05-01 21:18 | |

The patch was part of the merge of Mark's and my trunk-math branch. It was merged into the trunk and 3.0 a while ago. |

History | |||
---|---|---|---|

Date | User | Action | Args |

2008-05-01 21:18:02 | christian.heimes | set | status: open -> closed resolution: accepted messages: + msg66047 |

2008-01-22 11:10:55 | christian.heimes | set | messages: + msg61498 |

2008-01-22 02:10:29 | mark.dickinson | set | messages: + msg61479 |

2008-01-22 01:41:59 | christian.heimes | set | messages: + msg61474 |

2008-01-21 20:42:14 | mark.dickinson | set | files:
+ invhyp.c messages: + msg61444 |

2008-01-21 20:13:59 | mark.dickinson | set | messages: + msg61443 |

2008-01-21 19:39:53 | mark.dickinson | set | messages: + msg61442 |

2008-01-21 19:29:33 | mark.dickinson | set | messages: + msg61441 |

2008-01-21 18:10:23 | gvanrossum | set | nosy: - gvanrossum |

2008-01-21 18:05:05 | christian.heimes | set | files: - trunk_pymath_hyberbolic.patch |

2008-01-21 18:04:59 | christian.heimes | set | files: - trunk_mathmodule.patch |

2008-01-21 18:04:53 | christian.heimes | set | files: - trunk_math_sign_inf_nan.patch |

2008-01-21 18:04:45 | christian.heimes | set | files: - trunk_pymath_hyberbolic_complex.patch |

2008-01-21 18:04:33 | christian.heimes | set | files:
+ trunk_pymath_hyberbolic_complex2.patch messages: + msg61431 |

2008-01-21 17:24:19 | christian.heimes | set | messages: + msg61423 |

2008-01-21 17:14:49 | mark.dickinson | set | messages: + msg61419 |

2008-01-21 16:56:35 | christian.heimes | set | messages: + msg61414 |

2008-01-21 15:31:42 | mark.dickinson | set | messages: + msg61403 |

2008-01-21 15:25:01 | mark.dickinson | set | messages: + msg61402 |

2008-01-21 15:11:26 | mark.dickinson | set | messages: + msg61401 |

2008-01-21 14:53:14 | mark.dickinson | set | messages: + msg61399 |

2008-01-21 14:35:34 | christian.heimes | set | messages: + msg61398 |

2008-01-21 13:52:49 | christian.heimes | set | files:
+ trunk_pymath_hyberbolic_complex.patch messages: + msg61394 |

2008-01-21 12:31:09 | christian.heimes | set | messages: + msg61386 |

2008-01-21 12:24:50 | christian.heimes | set | messages: + msg61385 |

2008-01-21 03:34:24 | mark.dickinson | set | messages: + msg61376 |

2008-01-21 01:05:35 | tim.peters | set | messages: + msg61372 |

2008-01-21 01:02:13 | mark.dickinson | set | messages: + msg61371 |

2008-01-21 00:30:44 | gmcastil | set | messages: + msg61370 |

2008-01-20 22:54:34 | mark.dickinson | set | messages: + msg61367 |

2008-01-20 22:24:41 | gmcastil | set | messages: + msg61366 |

2008-01-20 17:09:36 | mark.dickinson | set | messages: + msg61311 |

2008-01-20 16:53:30 | christian.heimes | set | messages: + msg61310 |

2008-01-20 16:33:26 | mark.dickinson | set | messages: + msg61309 |

2008-01-20 16:11:36 | mark.dickinson | set | messages: + msg61308 |

2008-01-20 15:17:47 | mark.dickinson | set | messages: + msg61306 |

2008-01-20 13:56:08 | christian.heimes | set | files:
+ trunk_pymath_hyberbolic.patch messages: + msg61290 |

2008-01-20 08:57:07 | christian.heimes | set | messages: + msg60264 |

2008-01-20 03:07:15 | mark.dickinson | set | messages: + msg60259 |

2008-01-20 02:12:07 | gmcastil | set | nosy:
+ gmcastil messages: + msg60257 |

2008-01-04 03:51:25 | tim.peters | set | nosy:
+ tim.peters messages: + msg59211 |

2008-01-03 22:04:28 | gvanrossum | set | messages: + msg59166 |

2008-01-03 22:03:31 | christian.heimes | set | messages: + msg59165 |

2008-01-03 20:50:25 | gvanrossum | set | messages: + msg59156 |

2008-01-03 20:42:19 | christian.heimes | set | messages: + msg59155 |

2008-01-03 20:35:13 | gvanrossum | set | messages: + msg59154 |

2008-01-03 20:29:15 | mark.dickinson | set | messages: + msg59153 |

2008-01-03 20:21:15 | christian.heimes | set | messages: + msg59152 |

2008-01-03 17:28:15 | gvanrossum | set | messages: + msg59138 |

2008-01-03 16:57:46 | christian.heimes | set | messages: + msg59137 |

2008-01-03 00:11:00 | gvanrossum | set | messages: + msg59115 |

2007-12-18 23:28:28 | christian.heimes | set | files:
+ trunk_math_sign_inf_nan.patch messages: + msg58786 |

2007-12-18 21:05:41 | christian.heimes | set | messages: + msg58770 |

2007-12-18 19:36:45 | gvanrossum | set | nosy:
+ gvanrossum messages: + msg58758 |

2007-12-18 14:36:33 | mark.dickinson | set | messages: + msg58749 |

2007-12-18 10:06:24 | christian.heimes | set | files:
+ trunk_float_math_combined2.patch messages: + msg58735 |

2007-12-17 20:17:25 | christian.heimes | set | files:
+ trunk_float_math_combined.patch messages: + msg58706 |

2007-12-17 19:54:27 | christian.heimes | link | issue1635 superseder |

2007-12-17 19:27:44 | christian.heimes | set | messages: + msg58702 |

2007-12-17 17:53:35 | Rhamphoryncus | set | nosy:
+ Rhamphoryncus messages: + msg58699 |

2007-12-17 17:40:43 | mark.dickinson | set | nosy:
+ mark.dickinson messages: + msg58698 |

2007-12-17 16:04:29 | christian.heimes | set | dependencies: + Float patch for inf and nan on Windows (and other platforms) |

2007-12-17 15:48:44 | christian.heimes | create |