[Python-Dev] trunk-math

[Steve Holden]

Mark Dickinson wrote:
> Dear all,
> 
> I'd like to draw your attention to some of the work that's been going on 
> in the trunk-math branch.  Christian Heimes and I have been working on 
> various aspects of Python mathematics, and we're hoping to get at least 
> some of this work into Python 2.6/3.0.  Most of the changes are 
> completed or nearly complete, and 2.6/3.0 isn't very far away, so it 
> seems like a good time to try to get some feedback from python-dev.
> 
> Here's an overview of the changes (overview originally written by 
> Christian, edited and augmented by me.  I hope Christian will step in 
> and correct me if I misrepresent him or his work here.)  Many of the 
> changes were motivated by Christian's work (already in the trunk) in 
> making infinities and nans more accessible and portable for Python 
> users. (See issue #1640.)
> 
> * Structural reorganization: there are new files Include/pymath.h and 
> Objects/pymath.c with math-related definitions and replacement functions 
> for platforms without copysign, log1p, hypot and inverse hyperbolic 
> functions.
> 
> * New math functions: inverse hyperbolic functions (acosh, asinh, atanh).
> 
> * New float methods: is_finite, is_inf, is_integer and is_nan.
> 
> * New cmath functions: phase, polar and rect, isinf and isnan.
> 
> * New complex method: is_finite.
> 
> * Work on math and cmath functions to make them handle special values 
> (infinities and nans) and floating-point exceptions according to the C99 
> standard.   The general philosophy follows the ideas put forth by Tim 
> Peters and co. many moons ago. and repeated in the issue #1640 thread 
> more recently:  where the C99 standard (or IEEE 754) specifies raising 
> 'invalid' or 'divide-by-zero' Python should raise a ValueError.  Where 
> the C99 standard specifies raising 'overflow' Python should raise 
> OverflowError.  'underflow' and 'inexact' flags are ignored.  From a 
> user's perspective, this means that infinities and nans are never 
> produced by math or cmath operations on finite values (but see below). 
>  sqrt(-1) will always raise ValueError, instead of returning a NaN.  See 
> issue #711019 and the resulting warning at the bottom of the math module 
> documentation.   Although complex_abs doesn't live in cmathmodule.c, it 
> was also fixed up this way.
> 
> * The cmath module has been rescued:  it's no longer numerically unsound 
> (see issue #1381).  For the majority of functions (sin, cos, tan, sinh, 
> cosh, tanh, asin, acos, atan, asinh, acosh, atanh, exp, sqrt) the real 
> and imaginary parts of the results are always within a few ulps of the 
> true values.  (In extensive but non-exhaustive testing I haven't found 
> errors of more than 5 ulps in either the real or imaginary parts for 
> these functions.)  For log and log10 the errors can be larger when the 
> argument has absolute value close to 1; this seems pretty much 
> unavoidable without using multiple-precision arithmetic.  pow and 
> two-argument log are less accurate; again, this is essentially 
> unavoidable without adding hundreds of extra lines of code.
> 
> * Many more tests. In addition to a handful of extra test_* methods in 
> test_math and test_cmath, there are over 1700 testcases (in a 
> badly-named file Lib/test/cmath.ctest) for the cmath and math functions, 
> with a testcase format inspired in no small part by the decimal .decTest 
> file format.  Most of the testcase values were produced independently of 
> Python using MPFR and interval arithmetic (C code available on request); 
>  some were created by hand.
> 
> * There's a per-thread state for division operator. In IEEE 754 mode 
> 1./0. and 1.%0. return INF and 0./0. NAN. The contextlib has a new 
> context "ieee754" and the math lib set_ieee754/get_ieee754 (XXX better 
> place for the functions?)
> 
> Some notes:
> 
> * We've used the C99 standard (especially Annex F and Annex G) as a 
> reference for deciding what the math and cmath functions should do, 
> wherever possible. It seems to make sense to choose to follow some 
> standard, essentially so that all the hard decisions have been thought 
> through thoroughly by a group of experts.  Two obvious choices are the 
> C99 standard and IEEE 754(r); for almost all math issues the two say 
> essentially the same thing.  C99 has the advantage that it includes 
> specifications for complex math, while IEEE 754(r) does not.  (Actually, 
> I've been using draft version N1124 of the C99 standard, not the 
> standard itself, since I'm too cheap to pay up for a genuine version. 
> Google 'N1124' for a copy.)
> 
> * I'm offering to act as long-term maintainer for the cmath module, if 
> that's useful.
> 
> * The most interesting and exciting feature, by far (in my opinion) is 
> the last one.  By way of introduction, here's a snippet from Tim Peters, 
> in a comp.lang.python posting 
> (http://mail.python.org/pipermail/python-list/2005-July/330745.html), 
> answering a question from Michael Hudson about 1e300*1e300 and inf/inf.
> 
> "I believe Python should raise exceptions in these cases by default, 
> because, as above, they correspond to the overflow and invalid-operation 
> signals respectively, and Python should raise exceptions on the 
> overflow, invalid-operation, and divide-by-0 signals by default. But I 
> also believe Python _dare not_ do so unless it also supplies sane 
> machinery for disabling traps on specific signals (along the lines of 
> the relevant standards here). Many serious numeric programmers would be 
> livid, and justifiably so, if they couldn't get non-stop mode back. The 
> most likely x-platfrom accident so far is that they've been getting 
> non-stop mode in Python since its beginning."
> 
> Christian has found a simple, elegant and practical way to make it 
> possible for Python to raise these exceptions by default, while also 
> allowing serious numeric users access to non-stop mode---i.e., a mode 
> that generates inf from 1/0 instead of raising a Python exception.  (I 
> had a much more elaborate plan in mind, involving a thread-local context 
> similar to Decimal's, with control over individual traps and flags. 
>  Christian's solution is far more practical.)  The idea is that any 
> arithmetic operating under a "with ieee754:" acts like arithmetic on an 
> IEEE 754 platform with no traps enabled:  invalid operations like 
> sqrt(-1) produce a nan, division by zero produces an infinity, etc.  No 
> Python exceptions related to floating-point are raised.
> 
> See the thread started by Neal Becker 
> at http://mail.python.org/pipermail/python-list/2008-February/477064.html, 
> entitled "Turn off ZeroDivisionError?" for a recent discussion of these 
> issues.
> 
> I fear that the per-thread state change seems like something where a PEP 
> might be necessary; it's also not clear right now that Christian and I 
> have exactly the same goals here (discussion is ongoing).  But I hope 
> that the rest of the changes are uncontroversial enough to merit 
> consideration for possible inclusion in 2.6/3.0.
> 
> Thoughts?
> 
Only one: if this stunning plethora of improvements is likely to benefit 
by you and Christian having access to published standards I will happily 
champion a PSF-sponsored purchase, and would anticipate little argument.

Plus, thank you both very much for putting the effort in to satisfy 
numerical calculation requirements in an elegant and practical way. The 
immense amount of hard work involved should not go unrecognized.

regards
  Steve
-- 
Steve Holden        +1 571 484 6266   +1 800 494 3119
Holden Web LLC              http://www.holdenweb.com/