math.nroot [was Re: A brief question.]
mwh at python.net
Mon Jul 11 14:09:50 CEST 2005
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.
> 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
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:
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:
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:
/* stuff, e.g. r = exp(PyFloat_AS_DOUBLE(x)) */
/* 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>).
> 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".
> > 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
Decimal is a solved problem, here (and also an inspiration for
> >>> 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!
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
More information about the Python-list