[Python-Dev] trunk-math

[Mark Dickinson]

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?

Mark
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.python.org/pipermail/python-dev/attachments/20080215/1cc88a88/attachment-0001.htm