[Python-Dev] Rewrite of cmath module?

[Tim Peters]

[Mark Dickinson]
> The current cmath module has some significant numerical problems,
> particularly in the inverse trig and inverse hyperbolic trig
> functions.

IIRC, cmath was essentially whipped up over a weekend by a scientist,
transcribing formulas from a math reference.  Such an approach is
bound to suffer "significant numerical problems", but has the
advantage of taking little time to implement ;-)  A few of the worst
have been improved (e.g., see c_quot(), which contains the original
code in a comment), but no systematic effort has been made.

> ...

> I'm wondering whether anyone would be interested in a rewrite of the
> cmath module.  I have a drop-in replacement written in pure Python,
> based largely on the formulas given by Kahan in his `Much Ado about
> Nothing's Sign Bit' article, which I believe eliminates the vast
> majority of the current numerical problems.

I believe that :-)

> Besides the numerical fixes, the only significant difference from the
> current behaviour is that the branch cuts for asinh have been moved.
> (The docs suggest that the current branch cuts for asinh should be
> regarded as buggy).  I'm offering to translate this into C and add
> appropriate documentation and test cases, but I have a couple of
> questions before I start.

Make the Python implementation available first and solicit feedback
from the Numeric community?

> (1) Is there actually any interest in fixing this module?  I guess
> that the cmath module can't be used that much, otherwise these
> problems would have surfaced before.

That's the rub:  most people don't use cmath, while most who do have
no idea what these functions should return.  If there's a gross
problem, most wouldn't even notice.  The few places in cmath that have
been improved were due to users who did notice.  Of course all "should
be" fixed.

> (2) (Disregard if the answer to question 1 is `No'.)  What should be
> done about branch cuts?  The value of a function on a branch cut is
> going to depend on signs of zeros, so it'll be pretty much impossible
> to make any guarantees about the behaviour.  Even so, there are two
> approaches that seem feasible:

...

Follow the C99 standard whenever possible.  C99 added a complex type
to C, and added complex-valued math functions too.  Alas, C99 hasn't
been widely adopted yet, but Python should be compatible with the
platform C implementation to the extent possible.  Of course the
easiest way to do that is to /use/ the platform C implementation when
possible.

...

> (3) Is it worth trying to be consistent in exceptional cases?  The
> current math module documentation notes that math.log(0) might produce
> -inf on one platform and raise ValueError on another, so being consistent
> about whether a non-representable result gives a NaN, an infinity or a
> Python exception seems as though it would be tricky.  I should also note that
> I've made no effort to do the right thing when the argument to any of the
> functions contains a NaN or an infinity in either the real or imaginary part.

Welcome to the club ;-)  To the extent that your code relies on
real-valued libm functions, you inherit uncertainties from the
platform C's implementation of those too.  There's a long debate about
whether that's a feature (Python math functions act the same way as
the platform C's, and likely the platform Fortran's too) or a bug
(Python acts differently on different platforms).  "Feature" is the
traditional answer to that.

...
> [list of specific bad behaviors]
...

> One more thing: since this is my first post to python-dev I should
> probably introduce myself.  I'm a mathematician, working mostly in
> number theory.  I learnt programming and numerical analysis the hard
> way, coding service-life prediction algorithms for rocket motors in
> Fortran 77 during several long summers.  I've been using Python for
> more than ten years, mainly for teaching but also occasionally for
> research purposes.

Of course a strong background in math is very helpful for this.  If
you also spent several long summers cleaning sewers with your bare
hands, your qualifications for working on libm functions are beyond
reproach ;-)