math.nroot [was Re: A brief question.]

[Michael Hudson]

Tim Peters <tim.peters at gmail.com> writes:

> [Tim Peters]
> >> All Python behavior in the presence of infinities, NaNs, and signed
> >> zeroes is a platform-dependent accident, mostly inherited from that
> >> all C89 behavior in the presence of infinities, NaNs, and signed
> >> zeroes is a platform-dependent crapshoot.
>  
> [Michael Hudson]
> > As you may have noticed by now, I'd kind of like to stop you saying
> > this :) -- at least on platforms where doubles are good old-fashioned
> > 754 8-byte values.
> 
> Nope, I hadn't noticed!  I'll stop saying it when it stops being true,
> though <wink>.  Note that since there's not even an alpha out for 2.5
> yet, none of the good stuff you did in CVS counts for users yet.

Well, obviously.  OTOH, there's nothing I CAN do that will be useful
for users until 2.5 actually comes out.

> > But first, I'm going to whinge a bit, and lay out some stuff that Tim
> > at least already knows (and maybe get some stuff wrong, we'll see).
> >
> > Floating point standards lay out a number of "conditions": Overflow
> > (number too large in magnitude to represent), Underflow (non-zero
> > number to small in magnitude to represent), Subnormal (non-zero number
> > to small in magnitude to represent in a normalized way), ...
> 
> The 754 standard has five of them:  underflow, overflow, invalid
> operation, inexact, and "divide by 0" (which should be understood more
> generally as a singularity; e.g., divide-by-0 is also appropriate for
> log(0)).

OK, the decimal standard has more, which confused me for a bit
(presumably it has more because it doesn't normalize after each
operation).

> > For each condition, it should (at some level) is possible to trap each
> > condition, or continue in some standard-mandated way (e.g. return 0
> > for Underflow).
> 
> 754 requires that, yes.
> 
> > While ignoring the issue of allowing the user to control this, I do
> > wish sometimes that Python would make up it's mind about what it does
> > for each condition.
> 
> Guido and I agreed long ago that Python "should", by default, raise an
> exception on overflow, invalid operation, and divide by 0, and "should
> not", by default, raise an exception on underflow or inexact.

OK.

> Such defaults favor non-expert use.  Experts may or may not be happy
> with them, so Python "should" also allow changing the set.

Later :)

(In the mean time can we just kill fpectl, please?)

> > There are a bunch of conditions which we shouldn't and don't trap by
> > default -- Underflow for example.  For the conditions that probably should
> > result in an exception, there are inconsistencies galore:
> 
> > >>> inf = 1e300 * 1e300 # <- Overflow, no exception
> > >>> nan = inf/inf # <- InvalidOperation, no exception
> 
> Meaning you're running on a 754 platform whose C runtime arranged to
> disable the overflow and invalid operation traps.

Isn't that the standard-mandated start up environment?

> You're seeing native HW fp behavior then.

But anyway, shouldn't we try to raise exceptions in these cases?  I
don't think it's a particularly good idea to try to utilize the fp
hardware's ability to do this at this stage, btw, but to add some kind
of check after each operation.

> > >>> pow(1e100, 100) <- Overflow, exception
> > Traceback (most recent call last):
> >  File "<stdin>", line 1, in ?
> > OverflowError: (34, 'Numerical result out of range')
> > >>> math.sqrt(-1) # <- InvalidOperation, exception
> > Traceback (most recent call last):
> >  File "<stdin>", line 1, in ?
> > ValueError: math domain error
> 
> Unlike the first two examples, these call libm functions.

And the user cares about this why?

> Then it's a x-platform crapshoot whether and when the libm functions
> set errno to ERANGE or EDOM, and somewhat of a mystery whether it's
> better to reproduce what the native libm considers to be "an error",
> or try to give the same results across platforms.  Python makes a
> weak attempt at the latter.

Well, you can at least be pretty sure that an infinite result is the
result of an overflow condition, I guess.

> > At least we're fairly consistent on DivisionByZero...
> 
> When it's a division by 0, yes.  It's cheap and easy to test for that.
>  However, many expert uses strongly favor getting back an infinity
> then instead, so it's not good that Python doesn't support a choice
> about x/0.

Indeed.  But I'd rather work on non-settable predictability first.

> > If we're going to trap Overflow consistently, we really need a way of
> > getting the special values reliably -- which is what pep 754 is about,
> > and its implementation may actually work more reliably in 2.5 since my
> > recent work...
> 
> I don't know what you have in mind.  

Well, the reason I headbutted into this stuff again recently was
writing (grotty) string to float parsing code for PyPy.  If you write
(where 'e' is the exponent parsed from, say '1e309'):

while e > 0:
    result *= 10
    e -= 1

you get an infinite result in the large e case.  If instead you write
the seemingly much more sensible:

result *= pow(10, e)

you don't, you get an overflow error instead.  This make me sad.

Whether my code should have returned an infinite value or raised an
exception for "1e1000" is an open question, but I'd like the answer to
be less accidental.

> For example, checking the result of x*y to see whether it's an
> infinity is not a reliable way to detect overflow, and it fails in
> more than one way (e.g., one of the inputs may have been an infinity
> (in which case OverflowError is inappropriate), 

Well, yeah.  What other ways can it fail?

> and overflow doesn't always result in an infinity either (depends on
> the rounding mode in effect)).

If we're allowing the user to alter rounding mode, I imagine we'll
have a much better idea about this sort of thing.  For now, I think we
can assume a default rounding that doesn't do that, and someone who
has a platform like where it has his or her own problem.

> > On the issue of platforms that start up processes with traps enabled,
> > I think the correct solution is to find the incantation to turn them
> > off again and use that in Py_Initialize(), though that might upset
> > embedders.
> 
> Hard to know.  Python currently has a hack to disable traps on
> FreeBSD, in python.c's main().

Oh, right.  Well, I don't much care about embedding, anyway.

Cheers,
mwh

-- 
  MARVIN:  Oh dear, I think you'll find reality's on the blink again.
                   -- The Hitch-Hikers Guide to the Galaxy, Episode 12