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

[Tim Peters]

[Fredrik Johansson]
>>> I'd rather like to see a well implemented math.nthroot. 64**(1/3.0)
>>> gives 3.9999999999999996, and this error could be avoided.

[Steven D'Aprano]
>> >>> math.exp(math.log(64)/3.0)
>> 4.0
>>
>> Success!!!

[Tom Anderson]
> Eeeeeeenteresting. I have no idea why this works. Given that math.log is
> always going to be approximate for numbers which aren't rational powers of
> e (which, since e is transcendental, is all rational numbers, and
> therefore all python floats, isn't it?), i'd expect to get the same
> roundoff errors here as with exponentiation. Is it just that the errors
> are sufficiently smaller that it looks exact?

Writing exp(log(x)*y) rather than x**y is in _general_ a terrible
idea, but in the example it happens to avoid the most important
rounding error entirely:  1./3. is less than one-third, so 64**(1./3.)
is less than 64 to the one-third.  Dividing by 3 instead of
multiplying by 1./3. is where the advantage comes from here:

>>> 1./3.  # less than a third
0.33333333333333331
>>> 64**(1./3.)  # also too small
3.9999999999999996
>>> exp(log(64)/3)  # happens to be on the nose
4.0

If we feed the same roundoff error into the exp+log method in
computing 1./3., we get a worse result than pow got:

>>> exp(log(64) * (1./3.))  # worse than pow's
3.9999999999999991

None of this generalizes usefully -- these are example-driven
curiousities.  For example, let's try 2000 exact cubes, and count how
often "the right" answer is delivered:

c1 = c2 = 0
for i in range(1, 2001):
    p = i**3
    r1 = p ** (1./3.)
    r2 = exp(log(p)/3)
    c1 += r1 == i
    c2 += r2 == i
print c1, c2

On my box that prints

    3 284

so "a wrong answer" is overwhelmingly more common either way.  Fredrik
is right that if you want a library routine that can guarantee to
compute exact n'th roots whenever possible, it needs to be written for
that purpose.

...

> YES! This is something that winds me up no end; as far as i can tell,
> there is no clean programmatic way to make an inf or a NaN;

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.

> in code i write which cares about such things, i have to start:
>
> inf = 1e300 ** 1e300
> nan = inf - inf

That would be much more portable (== would do what you intended by
accident on many more platforms) if you used multiplication instead of
exponentiation in the first line.

...

> And then god forbid i should actually want to test if a number is NaN,
> since, bizarrely, (x == nan) is true for every x; instead, i have to
> write:
> 
> def isnan(x):
>        return (x == 0.0) and (x == 1.0)

The result of that is a platform-dependent accident too.  Python 2.4
(but not eariler than that) works hard to deliver _exactly_ the same
accident as the platform C compiler delivers, and at least NaN
comparisons work "as intended" (by IEEE 754) in 2.4 under gcc and MSVC
7.1 (because those C implementations treat NaN comparisons as intended
by IEEE 754; note that MSVC 6.0 did not):

>>> inf = 1e300 * 1e300
>>> nan == nan
>>> nan = inf - inf
>>> nan == 1.0
False
>>> nan < 1.0
False
>>> nan > 1.0
False
>>> nan == nan
False
>>> nan < nan
False
>>> nan > nan
False
>>> nan != nan
True

So at the Python level you can do "x != x" to see whether x is a NaN
in 2.4+(assuming that works in the C with which Python was compiled;
it does under gcc and MSVC 7.1).

> The IEEE spec actually says that (x == nan) should be *false* for every x,
> including nan. I'm not sure if this is more or less stupid than what
> python does!

Python did nothing "on purpose" here before Python 2.4.

> And while i'm ranting, how come these expressions aren't the same:
>
> 1e300 * 1e300
> 1e300 ** 2

Because all Python behavior in the presence of infinities, NaNs and
signed zeroes is a platform-dependent accident.

> And finally, does Guido know something about arithmetic that i don't,

Probably yes, but that's not really what you meant to ask <wink>.

> or is this expression:
>
> -1.0 ** 0.5
>
> Evaluated wrongly?

Read the manual for the precedence rules.  -x**y groups as -(x**y). 
-1.0 is the correct answer.  If you intended (-x)**y, then you need to
insert parentheses to force that order.