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

[Steven D'Aprano]

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.