Possible bug (was Re: numpy, overflow, inf, ieee, and rich comparison)

[Edward Jason Riedy]

Wow.  I almost missed an ieee754 discussion... ;)  Note that I'd
like to see ieee support in Python, but it's going to take a _lot_
of work.  Until then, consistent behavior makes sense.  It works
for other languages, for some definition of `works.'

And Tim Peters writes:
 - 
 - He is an entertaining writer, although his rhetoric is often more 
 - extreme than either of ours <wink>.

He's just as entertaining (and comprehensive) in person.  ;)  
There was a recent discussion about `merging' Rexx's decimal 
FP and ieee854 with the hopes of hardware implementation [1].  
It was educational.  Rexx's is still mostly a subset, although 
it does have an exceptional situation that may be useful (can't 
recall it at the moment, and my notes are in my office).  I
meant to poke at a software 854 implementation for Python 
(the I/O's easier than with 754), but other things have 
caught up with me.

 - Have you read the 754 std?  He's in favor of the 754 std.

Only because so few people are interested in ieee854 hardware
support.  (And there are a few bits that he wouldn't mind seeing
modified, as far as I can tell.  Requirements like an input NaN 
always being a signalling NaN aren't so useful after all...)

 - I agree with him there wholeheartedly:  without the 754-mandated flags,
 - *and* without raising exceptions, the subset of 754 remaining sucks.

You don't mean what you think you mean.  ;)  More later...

 - [matching experiences]:  a vast majority still want divide-by-0, 
 - invalid operation, and overflow to raise an exception the vast 
 - majority of the time.

That doesn't make them right.

 - No.  I want exactly the choices 754 lays out:  user-settable sticky flags
 - and user-maskable exceptions.

Note that ieee754 traps are not the same as the standard view of
exceptions in programming languages; they're much lower level.
Exceptions destroy things as they pass upwards, and that doesn't
help debugging an odd FP problem.  

The traps in ieee754 and 854 are of the style seen in hardware 
(and, notably, many Lisp-based languages).  You jump to a handler, 
passing it the address of the failed instruction and some stored
state (more or less a continuation).  The handler can then do 
whatever it wants, including substituting another value and 
continuing after the insn that trapped.

 - Kinda, except that each of the 5 754-defined traps must be individually
 - maskable.  754 spells out what's needed.

Not fully.  Maskable on a per-process or per-thread basis?  And
do threads and processes inherit their parent's defaults or the
ieee754 defaults?

This gets more interesting when you add in things like the 
`microthreads' people seem to like (threads as programming
language construct rather than system construct).

CPython's currently at a point where people can make fairly 
arbitrary choices on these issues and simply enforce them.  
That's because it's slow.  As it gets faster and faster, these
issues will keep coming up.  Amusingly enough, this is similar
to how Rexx slipped decimal fp into common use.  It was added
when Rexx was slow, so no one cared about speed, but they all
liked the functionality.  (Whether or not it's sufficiently
used beyond fixed-point is another question.)

[On `presubstitution']
 - The basic idea is that the programmer gets to specify, via library 
 - calls (whatever), which particular fp value each 754 exceptional 
 - case returns.

It's a straight-forward application of ieee{7,8}54 traps.  That's
essentially how software i87 emulation works, too...  

 - Then set it back to +- Inf (or anything else they needed) later.

This is the other problem with traps and the like.  When is later?
If a thread of execution is weaving between different libraries,
as is common in OO systems, when do you set up the trap handlers?
(Consider an ODE solver with plug-in functions.)  And can the
user override all lower-level choices to insert `debugging' traps?
How about some of them, possibly on a per-module basis, to keep 
from being drowned (or surprised by finding non-error traps in
an unrelated module)?

The current next implementation of decimal arithmetic in Java [2]
requires an FP environment object to be passed into every operation
to deal with some of these issues.  Obviously, that's a touch
painful.  AFAIK, no one has a good solution yet.

There's a related issue that would be nice to solve at the same
time:  Deciding in which precision to compute a quantity.  Python's
decision that all floating point numbers are doubles isn't always
the best.  I'm tired, so I'm failing to cook up an example, but
here's the general idea...  

Consider what happens when you string together different libraries 
written by different people and you need to shuttle FP numbers from 
one to the next.  How do you know you have a reasonably correct 
answer after churning through all these calculations?  Assume that 
the libraries do the best, most accurate job they can individually.
Annoyingly enough, you can't say much about your final answer.

Take an extreme reduction:  The result of every single ieee754
operation is ``the correct result correctly rounded,'' but not
all algorithms that use these operations produce an accurate 
answer to the same precision as the individual operations.

The idea that's been presented (by someone I mention far too 
often in these messages) is to pin down some intermediate 
results.  Make sure you know them, and you know how well you 
know them.  One way that works for strings of basic operations 
is the old K&R C way: floats come in, floats go out, but all 
intermediate results are doubles.  This doesn't generalize 
perfectly, as we don't have an infinite recursion of precisions 
to use, but it seems to give a good engineering guideline:  
Compute intermediate results to high precision, but use the
end result to a lower precision to have more confidence that
the answer means something.  That statement's very fuzzy,
rather like my thinking after getting up to go to a farmer's
market way too early this morning, but I think you can get 
what I mean.  The devil's in the details, as always.

How is this related to trap handlers and exceptional behavior?
Both rely on some FP context being passed across expressions,
and they need to FP context to depend on program context.  If 
an elegant solution can be found, it could work for both.  But
that's a big if.  I think it might be sensible in a dynamic
language, but I think most things might be sensible in dynamic
languages.

 - If you can find two NumPy users in favor of disabling them all by 
 - default, I'll be mildly astonished.

Do he and I count?  ;)

 - Scared the hell out of language stds committees too, which goes a lot 
 - farther than the 754 folks would like to admit toward explaining why 
 - language committees stayed away from 754 bindings in droves.  

Dr. Kahan has said as much, but he has stuck to his guns.  It is 
not a trivial point, it's a different way of thinking.  If all 
people had stuck to the overflow-will-trap way of thinking, they 
would still be looking for an optimal symmetric eigenpair algorithm.
Likewise, lazy language people would likely still be thinking too 
hard about exceptions in lazy languages (Jones, et al. seem to 
credit ieee fp's model for their inspiration).

Not having to stick silly tests around virtually every FP operation
takes a huge mental weight off.  You don't have to worry that some
odd data will cause some inconsequential calculation to go nuts
and stop the whole program.

 - C99 is 15 years(!) after the fact.

And ``common wisdom'' is that changes in thinking take about 20
years to trickle out to the world at large, so it's a bit ahead of
schedule.

 - Even $3 business calculators have a sqrt key these days, and you do
 - not want to be on the other end of email when trying to explain to 
 - a businessperson why their crappy std deviation program returned a 
 - complex number <0.1 wink>.

Because its core algorithms weren't designed by Dr. Kahan (HP71B's 
arithmetic and equation solver)?  ;)

Jason, who has negative time right now and is going more negative
by babbling...

[1] http://www2.hursley.ibm.com/decimal/
[2] http://www2.hursley.ibm.com/decimalj/