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

Well, OK, I phrased my first post badly.  Let me try again:

I want to make this situation better, as you may have noticed.

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

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

Maybe we should just implement floats in Python.

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

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

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

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?
> 
> I believe Python should raise exceptions in these cases by default,
> 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.

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

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

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

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

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

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

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

[...]

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

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

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

> 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?  How come they
haven't tried to reduce the mess before?

I guess "academia" is an answer to both those questions...

[...]

> >> 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.
> 
> IME the rounding modes in 754 are almost never changed from the
> default to-nearest-even, and I can live with that.  The rounding modes
> in the Decimal module will be changed a lot, though (due to the
> insanity of politically mandated rounding rules in currency
> calculations).

Decimal is a solved problem, here (and also an inspiration for
design).

> >>> 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.
> 
> So far, I think FreeBSD remains the only known Python platform that
> enables some FP traps by default.  So it's OK if you _just_ don't care
> much about embedding on FreeBSD <wink>.

Even less to worry about!

Cheers,
mwh

-- 
  I never realized it before, but having looked that over I'm certain
  I'd rather have my eyes burned out by zombies with flaming dung
  sticks than work on a conscientious Unicode regex engine.
                                             -- Tim Peters, 3 Dec 1998