math.nroot [was Re: A brief question.]
tim.peters at gmail.com
Sun Jul 3 19:09:34 CEST 2005
>>> I'd rather like to see a well implemented math.nthroot. 64**(1/3.0)
>>> gives 3.9999999999999996, and this error could be avoided.
>> >>> math.exp(math.log(64)/3.0)
> 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
>>> 64**(1./3.) # also too small
>>> exp(log(64)/3) # happens to be on the nose
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
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
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
> 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
> 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
>>> nan < 1.0
>>> nan > 1.0
>>> nan == nan
>>> nan < nan
>>> nan > nan
>>> nan != nan
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.
More information about the Python-list