[Python-Dev] calculator module

[Tim Peters]

[Raymond Hettinger]
> ...
> Also, I'll clarify the part about numerical accuracy.  The other
> respondent took "pretty good" as meaning sloppy.  What was meant was
> using the best algorithms and approximations that can be implemented
> cleanly and portably.  This is a high standard, but much lower than
> would be expected of a system library (down to the last ULP or pretty
> close to it).

It's also very much lower than the HP calculators you (& I) remember so
fondly.  A large part of why they were so pleasant to use is that HP went to
*extraordinary* lengths to make them numerically robust.  Indeed, they hired
Kahan to design and implement robust algorithms for the HP Solver, the HP
adaptive numerical Romberg integration, and the HP financial calculators'
IRR solvers (note that an IRR solver is essentially finding the roots of an
N-degree polynomial, so there are N roots, and it takes extreme care to
isolate the one that makes physical sense; & if cash flows are both positive
and negative, there can even be more than one root that makes physical
sense).  They were all solid advances in the numeric state of the art.
Kahan wrote some articles about them which are worth tracking down.  Alas, I
think they appeared in the HP Systems Journal, and offhand I wasn't able to
find them online just now.

> A practical and important consequence of "pretty good" is that we can
> ignore bug reports that say "last two digits of gamma(3.1) don't agree
> with Matlab."

Another professional thing HP did is document worst-case error bounds on all
their function implementations.  It's appalling that most C vendors don't
(and Python inherits this sad story).  It's fine by me if, for example, a
numeric function is guaranteed good to only one digit, provided that the
documentation says so.

BTW, it should be much easier to write high-quality implementations using
Decimal than using native binary fp.  The library writer building on Decimal
can easily request a temporary increase in working precision, then round the
final extended-precision result back.  With limitations, even the simplest
numerically naive algorithms can often be coaxed into delivering
high-quality results that way.  95% of the battle when using native binary
fp is finding clever ways to rephrase computation so that native-precision
operations don't go insane.

Simple example:  hypot(x, y) = sqrt(x**2 + y**2).  The naive implementation
in native fp is useless, because it suffers spurious overflow and underflow
across a huge section of the input domain.  Use an extended-range,
extended-precision arithmetic internally, though, and the naive
implementation is both accurate and robust across the entire domain.