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

[Tim Peters]

[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.

[Tim]
>>>> 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.

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

[Tim]
>> Sure.  I was explaining why I keep saying what you say you don't want
>> me to say:  until 2.5 actually comes out, what purpose would it serve
>> to stop warning people that 754 special-value behavior is a x-platform
>> crapshoot?  Much of it (albeit less so) will remain a crapshoot after
>> 2.5 comes out too.

[Michael] 
> Well, OK, I phrased my first post badly.  Let me try again:
>
> I want to make this situation better, as you may have noticed.

Yup, I did notice that!  I've even pointed it out.  In a positive light <wink>.

>>>>> 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).

>> The "conditions" in IBM's decimal standard map, many-to-one, on to a
>> smaller collection of "signals" in that standard.  It has 8 signals:
>> the 5 I named above from 754, plus "clamped", "rounded", and
>> "subnormal".  Distinctions are excruciatingly subtle; e.g., "rounded"
>> and "inexact" would be the same thing in 754, but, as you suggest, in
>> the decimal standard a result can be exact yet also rounded (if it
>> "rounds away" one or more trailing zeroes), due to the unnormalized
>> model.
 
> Right, yes, that last one confused me for a while.
> 
> Why doesn't 754 have subnormal?

Who cares <0.1 wink>.

>  Actually, I think I'm confused about when Underflow is signalled -- is it
> when a denormalized result is about to be returned or when a genuine
> zero is about to be returned?

Underflow in 754 is involved -- indeed, the definition is different
depending on whether the underflow trap is or is not enabled(!).  On
top of that, it's not entirely defined -- some parts are left to the
implementer's discrection.  See

     http://www2.hursley.ibm.com/decimal/854M8208.pdf

for what, from a brief glance, looks very much like the underflow text
in the final 754 standard.

Note that subnormals and underflow are much less a real concern with
754 doubles than with 754 floats, because the latter have such a small
dynamic range.

>>>>> 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.

>> And, I'll add, "should not" on rounded, clamped and subnormal too.

> Sure, but we already have a conforming implementation of 854 with
> settable traps and flags and rounding modes and all that jazz.

No, we don't, but I assume you're talking about the decimal module. 
If so, the decimal module enables traps on overflow, invalid
operation, and divide-by-0 by default.  A conforming implementation
would have to disable them by default.

Apart from that difference in defaults, the decimal module does intend
to conform fully to the proposed decimal FP standard.

> Maybe we should just implement floats in Python.

Certainly the easiest way to get 754 semantics across boxes!  Been
there, done that, BTW -- it's mondo slow.

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

>>> Later :)

>> That's a problem, though.  754 subsets are barely an improvement over
>> what Python does today:

> Well, my contention is that the consistent application of one
> particular 754 subset would be an improvement.  Maybe I'm wrong!

I wrote more about that later (as in exploiting that infix * doesn't
raise OverlowError, but pow(float, float) probably does).

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

>> Has it been marked as deprecated yet (entered into the PEP for
>> deprecated modules, raises deprecation warnings, etc)?  I don't know.
>> IMO it should become deprecated, but I don't have time to push that.

> A bit of googling suggests that more people pass --with-fpectl to
> configure than I expected, but I doubt more than 1% of those actually
> use the features thus provided (of course, this is a guess).

I expect 1% is way high.  Before we stopped building fpectl by
default, Guido asked and heard back that there were no known users
even at LLNL anymore (the organization that contributed the code).

>>>>> 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

Remember this below.

>>>> 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?

>> The 754 standard mandates "non-stop" mode (all traps disabled) at
>> startup, but who in this chain is _claiming_ to implement the 754
>> standard?  Your platform C may or may not, and your OS may or may not.
 
> Well, fair point, but it's a convenient stick to beat unconforming
> platform vendors with.

Yup.

> I only really care about platforms that make at least a half-hearted
> effort towards 754 compliance.  I'm not looking to improve behaviour
> on a PDP-10 running TOPS (or whatever).

>>>> You're seeing native HW fp behavior then.

>>> But anyway, shouldn't we try to raise exceptions in these cases?

Note that the only cases you could have been talking about here were
the plain * and / examples above.

>> I believe Python should raise exceptions in these cases by default,

And I was talking about the plain * and / examples too.

>> 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).

> Um.  Python raises OverflowError more often that it used to, it seems
> to me, and already _dares_ to raise it from (e.g.) math.exp with no
> way of turning this off I can see.

I was talking about * and / (and, by extension, also + and -).  As I
originally said, libm is a different story (for implementation
reasons).  I don't believe Python could get away with raising
exceptions on overflow or invalid operation for + - * or / unless
there were also a way to get back current behavior (which is non-stop
+ - * / on almost all (and possibly all) current Python platforms,
with the exception of divide-by-0).

>> Many serious numeric programmers would be livid, and justifiably so,
>> if they couldn't get non-stop mode back.

> Maybe I'm missing something (and writing this reply in multiple stages
> makes this more likely) but they don't have it now, do they?

They do for + - * /.

>> The most likely x-platfrom accident so far is that they've been
>> getting non-stop mode in Python since its beginning.

Again referring to + - * /.
 
> I know I'm repeating myself, but didn't this changeset:
> 
> http://fisheye.cenqua.com/changelog/python?cs=MAIN:tim_one:20010905223655
> 
> change that rather?  (It affects 2.2 and later, I think).

That was about libm error reporting, not + - * /.  It also wasn't
trying to _increase_ libm error coverage, it was trying to get back
the error coverage we _had_ before some C compilers started exploiting
the much looser libm error reporting allowed by C99.  C89 required
more from libm, and as people moved to newer C compilers Python
stopped raising exceptions where it had raised them before, because
their libm stopped setting errno in nasty cases.  The patch above
tried to get former behavior back.

Ideally, users would have a choice between exceptions or non-stop mode
for libm functions too.

>>> 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.

>> That's possible, anyway.  Doing it in software is slow, but if you
>> look at the horrid contortions fpectlmodule.c endures to try to
>> exploit HW facilities, it's not clear that the latter is any faster.

> I think it might well be slow, too, at least some of the time; the
> conclusion of this rather long thread:
> 
> http://lists.freebsd.org/pipermail/freebsd-stable/2005-May/014575.html
> 
> was "turn --with-fpectl off".

LOL!   That's a classic -- thanks for pointing it out.  There's indeed
little you can do that's slower than slopping setjmp/longjmp around
every piddly 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?

>> Didn't say the user did care.

> I know you didn't, and I also knew the reason for the different
> behaviours, but that doesn't make me feel that much better.

Maybe we could stick to questions you do want answers to, then?

>> Why doesn't Python already supply a fully 754-conforming arithmetic
>> on 754 boxes?  It's got almost everything to do with implementation
>> headaches, and very little to do with what users care about.
>> Because all the C facilities are a x-platform mess, the difference
>> between calling and not calling libm can be the difference between
>> using the platform libm or Python needing to write its own libm.
>> For example, there's no guarantee that math.sqrt(-1) will raise
>> ValueError in Python, because Python currently relies on the
>> platform libm sqrt to detect _and report_ errors.  The C standards
>> don't require much of anything there.

> Can't we use the stuff defined in Appendix F and header <fenv.h> of
> C99 to help here?  I know this stuff is somewhat optional, but it's
> available AFAICT on the platforms I actually use (doesn't mean it
> works, of course).

It's entirely optional part of C99.  Python doesn't require C99.  The
most important example of a compiler that doesn't support any of that
stuff is Microsoft's, although they have their own MS-specific ways to
spell most of it.

> I'm thinking something like this:
>
> fexcept_t flags;
>
> feclearexcept(FE_ALL_EXCEPT);
>
> /* stuff, e.g. r = exp(PyFloat_AS_DOUBLE(x)) */
>
> fegetexceptflag(&flags, FE_ALL_EXCEPT)
> 
> /* inspect flags to see if any of the flags we're currently trapping
>   are set */

Assuming the platform libm sets 754 flags appropriately, that's a fine
way to proceed on platforms that also support that specific spelling.

...

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

> > There are at least two other causes:  some cases of divide-by-0 (like
> > 1/0 returns +Inf), and non-exceptional production of an infinite
> > result from infinite operands (like sqrt(+Inf) should return +Inf, and
> > there's nothing exceptional about that).

> Yeah, but I think those can be dealt with (if we really wanted to).

They certainly could be.

> OTOH, I think Inexact would be pain itself to detect without hw
> support, and Underflow would be tedious.

In the absence of HW help, there's really no hope for those short of
emulating fp arithmetic in SW.

>>>>> 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.

>> OTOH, I bet I could find people willing to pay for the ability to
>> "turn off" divide-by-0 exceptions (or the people talking my ear off at
>> PyCon this year were just bluffing <wink>).

> Noted.

>> BTW, since there's so little the HW can help us with here in reality
>> (since there's no portable way to get at it),

> In what way does C99's fenv.h fail?  Is it just insufficiently
> available, or is there some conceptual lack?

Just that it's not universally supported.  Look at fpectlmodule.c for
a sample of the wildly different ways it _is_ spelled across some
platforms.  A maze of #ifdefs could work too, provided we defined a
PyWhatever_XYZ API to hide platform spelling details.

>> any predictability is going to be supplied by software, and software
>> for that will work better over the long run if it's designed from
>> the start realizing that it needs to consult some form of context to
>> decide what is and is not "an exception".
 
> Also noted.

> [...]
>>> 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.

>> It's a delicate house of cards right now for sure.  In some ways,
>> experts have learned to live with-- and exploit --it:  when they
>> _want_ an overflow exception they'll use pow(), and when they don't
>> they'll just multiply.

> Another question: where can I find these experts?

I'm one.  I'd expect to find others on the numerical Python lists, not
on c.l.py.  The folks at Enthought would be good to talk with.

> How come they haven't tried to reduce the mess before?

It's tedious, exacting, time-consuming, and an endless source of
x-platform problems.  I'd much rather work on the Python Challenge
riddles now <0.5 wink>, speaking of which ...