17

I want to find the greatest integer less than or equal to the kth root of n. I tried

int(n**(1/k))

But for n=125, k=3 this gives the wrong answer! I happen to know that 5 cubed is 125.

>>> int(125**(1/3))
4

What's a better algorithm?

Background: In 2011, this slip-up cost me beating Google Code Jam problem Expensive Dinner.

2
  • 3
    Rounding is causing what you observe. Here, this is because 1/3 is rounded down. I don't know about IEEE754 compliance of Python, but on other machines you could have observed 5. Apr 12, 2013 at 19:04
  • @Alexandre Huh, 125**(1/3) -> 4.999999999999999
    – wjandrea
    Jan 20, 2022 at 20:48

11 Answers 11

15

One solution first brackets the answer between lo and hi by repeatedly multiplying hi by 2 until n is between lo and hi, then uses binary search to compute the exact answer:

def iroot(k, n):
    hi = 1
    while pow(hi, k) < n:
        hi *= 2
    lo = hi // 2
    while hi - lo > 1:
        mid = (lo + hi) // 2
        midToK = pow(mid, k)
        if midToK < n:
            lo = mid
        elif n < midToK:
            hi = mid
        else:
            return mid
    if pow(hi, k) == n:
        return hi
    else:
        return lo

A different solution uses Newton's method, which works perfectly well on integers:

def iroot(k, n):
    u, s = n, n+1
    while u < s:
        s = u
        t = (k-1) * s + n // pow(s, k-1)
        u = t // k
    return s
8
  • Thanks. Test them with Peter's examples, both give the correct answers. Apr 13, 2013 at 19:52
  • Would this be guaranteed to find the exact integer root? Sep 11, 2014 at 0:18
  • @Ben: Yes. Both methods do.
    – user448810
    Sep 11, 2014 at 12:46
  • 1
    Why is that? Isn't it possible that it gets stuck on a different integer? Sep 11, 2014 at 13:31
  • @Ben: Do you have an example that fails?
    – user448810
    Sep 11, 2014 at 13:52
15

How about:

def nth_root(val, n):
    ret = int(val**(1./n))
    return ret + 1 if (ret + 1) ** n == val else ret

print nth_root(124, 3)
print nth_root(125, 3)
print nth_root(126, 3)
print nth_root(1, 100)

Here, both val and n are expected to be integer and positive. This makes the return expression rely exclusively on integer arithmetic, eliminating any possibility of rounding errors.

Note that accuracy is only guaranteed when val**(1./n) is fairly small. Once the result of that expression deviates from the true answer by more than 1, the method will no longer give the correct answer (it'll give the same approximate answer as your original version).

Still I am wondering why int(125**(1/3)) is 4

In [1]: '%.20f' % 125**(1./3)
Out[1]: '4.99999999999999911182'

int() truncates that to 4.

12
  • Maybe == should be replaced with <= since testing for float equality is often treacherous.
    – unutbu
    Apr 12, 2013 at 18:59
  • @unutbu: In that expression everything's integer ("Here, both val and n are expected to be integer and positive").
    – NPE
    Apr 12, 2013 at 18:59
  • 3
    @Andrey: Try typing 125**(1./3) into the interpreter. You should get something like 4.999999999999999, not 5. int floors it to 4.
    – unutbu
    Apr 12, 2013 at 19:01
  • 1
    @EricPostpischil: Because that would fail on 4 ** (1./3) and many other inputs, whereas the method in my answer won't.
    – NPE
    Apr 12, 2013 at 19:18
  • 2
    Be careful with big numbers: nth_root((1020+2)**2,2)=1020 with the method in this answer. Apr 12, 2013 at 19:19
4

My cautious solution after being so badly burned:

def nth_root(N,k):
    """Return greatest integer x such that x**k <= N"""
    x = int(N**(1/k))      
    while (x+1)**k <= N:
        x += 1
    while x**k > N:
        x -= 1
    return x
4
  • Sometimes I've found it useful to use this method in conjunction with a variable step size to accelerate convergence (only needed when you are using really large numbers!) Apr 12, 2013 at 19:29
  • @Peter, so the first guess really can be wrong by more than 1? Apr 12, 2013 at 19:30
  • You get of course burned by the N**(1./k) if N is too large to be converted to float. For small enough N, N**(1./k) is close, so the adjustment won't take long. If you're interested in treating large numbers, where the conversion to float loses too much precision, so the adjustment could take long, a few Newton-Raphson steps would get the precise result. Apr 12, 2013 at 19:30
  • 1
    int(((1040)**2)**(1.0/2))=1040+303786028427003666890752, so it can be off by a long way! Apr 12, 2013 at 19:31
3

Why not to try this :

125 ** (1 / float(3)) 

or

pow(125, 1 / float(3))

It returns 5.0, so you can use int(), to convert to int.

1

Here it is in Lua using Newton-Raphson method

> function nthroot (x, n) local r = 1; for i = 1, 16 do r = (((n - 1) * r) + x / (r ^ (n -   1))) / n end return r end
> return nthroot(125,3)
5
> 

Python version

>>> def nthroot (x, n):
...     r = 1
...     for i in range(16):
...             r = (((n - 1) * r) + x / (r ** (n - 1))) / n
...     return r
... 
>>> nthroot(125,3)
5
>>> 
3
  • nthroot(10**4,4)=32 with this method, I would have expected 10? Apr 12, 2013 at 19:39
  • You need to increase the number of iterations for larger numbers; for 10**4 range(32) is plenty. Apr 12, 2013 at 20:00
  • This method seems exceedingly slow to me. It takes 21 times around the loop to compute nthroot(10**4,4) and it takes 423 times around the loop the compute nthroot(10**15,15). Oct 21, 2018 at 4:47
1

I wonder if starting off with a method based on logarithms can help pin down the sources of rounding error. For example:

import math
def power_floor(n, k):
    return int(math.exp(1.0 / k * math.log(n)))

def nth_root(val, n):
    ret = int(val**(1./n))
    return ret + 1 if (ret + 1) ** n == val else ret

cases = [
    (124, 3),
    (125, 3),
    (126, 3),
    (1, 100),
    ]


for n, k in cases:
    print "{0:d} vs {1:d}".format(nth_root(n, k), power_floor(n, k))

prints out

4 vs 4
5 vs 5
5 vs 5
1 vs 1
1
def nth_root(n, k):
    x = n**(1./k)
    y = int(x)
    return y + 1 if y != x else y
0

You can round to nearest integer instead of rounding down / to zero (I don't know what Python specifies) :

def rtn (x):
    return int (x + 0.5)

>>> rtn (125 ** (1/3))
5
2
  • 2
    rtn(4 ** (1./3)) is 2, but 2**3 > 4.
    – unutbu
    Apr 12, 2013 at 19:10
  • @mbomb007: It seems that it works if x is non negative (this needs to be checked against the spec -- if there is round-to-zero in Python, then this needs to be adjusted for negative numbers). But non integer powers are well defined only for non negative numbers ! Apr 15, 2015 at 7:18
0

int(125**(1/3)) should clearly be 5, i.e. the right answer, so this must be standard computer rounding error, i.e internally the result is 4.9999999999 which gets rounded down to 4. This problem will exist with whatever algorithm you use. One simple ad-hoc solution is to add a tiny number e.g. int((125**(1/3)) + 0.00000001)

10
  • 2
    The problem is with the 1/3, not the 5: 5 is exactly representable by a floating point number. Apr 12, 2013 at 19:05
  • Yes, 5 is perfectly represented as a floating point number, but you'd still have this problem with e.g. (100**(1/2)) even though all of 100, (1/2) and the answer 10 are perfectly representable. The problem is that to take powers, computers use logs, and the logs of these numbers aren't perfectly representable. Apr 12, 2013 at 19:11
  • 2
    Adding 0.00000001 is just a hack. For example, it would make the algorithm fail on taking the 20th root of 95367431259155: the result would be 5 instead of 4.
    – NPE
    Apr 12, 2013 at 19:16
  • @NPE, yes true. "ad-hoc" = "hack" in my language LOL. The fact is, it's impossible to beat rounding error all the time. Instead of 0.00000001, a number much closer to the accuracy of the machine you're using would probably be optimal. Apr 12, 2013 at 19:17
  • @Stochastically: For integer inputs it is possible to avoid rounding errors. See my answer for one such method.
    – NPE
    Apr 12, 2013 at 19:19
0

Do this before everything:

from __future__ import division

and then run any of the above specified techniques to have your results.

0
0
def nthrootofm(a,n):
    a= pow(a,(1/n))
    return 'rounded:{},'.format(round(a))
a=125
n=3
q=nthrootofm(a,n)
print(q)

just used a format string , maybe this helps.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.