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

```