math.nroot [was Re: A brief question.]
Steven D'Aprano
steve at REMOVETHIScyber.com.au
Sun Jul 3 05:37:26 CEST 2005
On Sun, 03 Jul 2005 02:22:23 +0200, Fredrik Johansson wrote:
> On 7/3/05, Tom Anderson <twic at urchin.earth.li> wrote:
>> That's one way. I'd do:
>>
>> root = value ** 0.5
>>
>> Does that mean we can expect Guido to drop math.sqrt in py3k? :)
>
> I'd rather like to see a well implemented math.nthroot. 64**(1/3.0)
> gives 3.9999999999999996, and this error could be avoided.
py> 64**(1/3.0)
3.9999999999999996
Use the fact that if z = x**y, log(z) = y*log(x).
py> import math
py> math.exp(1/3.0*math.log(64))
3.9999999999999991
That's even less accurate.
py> math.exp(math.log(64)/3.0)
4.0
Success!!!
Another (random-ish) example:
py> 17**5
1419857
py> 1419857**(1.0/5)
17.000000000000004
py> math.exp(math.log(1419857)/5.0)
17.0
def nroot(x, y):
"""Returns the yth root of x.
More accurate than using x**(1/y).
Handles most, but not all, edge cases correctly. (Does not check for
y as either negative zero or negative infinity.)
"""
if x < 0.0:
raise ValueError("Complex valued result.")
if y == 0.0: # x**inf
if x > 1.0:
try:
return float("inf")
except:
try:
return 1.0/0.0
except:
raise ValueError("Overflow")
elif x == 1.0: # 1**N = 1 for any positive N
return 1.0
else: # 0**N = 0 for any non-zero N
return 0.0
# check the limiting case where y is positive infinity
if 1.0/y == 0.0:
# x**0 = 1 for all real x (including the limiting case x = 0)
return 1.0
if x == 0.0: # 0**N = 0 for any non-zero N
return 0.0
else:
return math.exp(math.log(x)/y)
Note how much simpler this would be if we could guarantee proper
infinities and NaNs in the code. We could turn a 23-line block to a
one-liner.
--
Steven.
More information about the Python-list
mailing list