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

[Michael Hudson]

I doubt anyone else is reading this by now, so I've trimmed quotes
fairly ruthlessly :)

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

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

!

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

Uh, yes.  Apologies for the lack of precision.

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

Right, that's what I meant.

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

No doubt.

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

Interesting.

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

Ah, OK, this distinction passed me by.

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

Hmm, is <fenv.h> optional?  I'm not finding those words.  I know
Appendix F is.

> Python doesn't require C99.

Sure.  But it would be possible to, say, detect C99 floating point
facilities at ./configure time and use them if available.

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

OK, *that's* a serious issue.

If you had to guess, do you think it likely that MS would ship fenv.h
in the next interation of VC++?

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

It even seems to work, on darwin/ppc (so, with GCC) at least.

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

The more I think about it, the less wise I think detecting stuff this
was is sane.

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

C'mon, fpectlmodule.c is _old_.  Maybe I'm stupidly optimistic, but
perhaps in the last near-decade things have got a little better here.

> A maze of #ifdefs could work too, provided we defined a
> PyWhatever_XYZ API to hide platform spelling details.

Hopefully it wouldn't be that bad a maze; frankly GCC & MSVC++ covers
more than all the cases I care about.

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

OK.

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

Heh.

Cheers,
mwh

-- 
  For every complex problem, there is a solution that is simple,
  neat, and wrong.                                    -- H. L. Mencken