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.

Colonel Panic
  • 132,665
  • 89
  • 401
  • 465
  • 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. – Alexandre C. Apr 12 '13 at 19:04
  • @Alexandre Huh, `125**(1/3)` -> `4.999999999999999` – wjandrea Jan 20 '22 at 20:48

11 Answers11

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.

NPE
  • 486,780
  • 108
  • 951
  • 1,012
  • Maybe `==` should be replaced with `<=` since testing for float equality is often treacherous. – unutbu Apr 12 '13 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 '13 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 '13 at 19:01
  • @Andrey: I've added an explanation to the answer. – NPE Apr 12 '13 at 19:02
  • So better use `math.ceil`, it is more elegant that `if`. – Andrey Apr 12 '13 at 19:02
  • because you are getting something like 4.99999 and int( 4.9999 ) = 4 – Joran Beasley Apr 12 '13 at 19:03
  • it wont alway need to ciel sometimes it will be floored ... also ciel returns a float afaik ... – Joran Beasley Apr 12 '13 at 19:03
  • 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 '13 at 19:18
  • 2
    Be careful with big numbers: nth_root((10**20+2)**2,2)=10**20 with the method in this answer. – Peter de Rivaz Apr 12 '13 at 19:19
  • I think the simplest robust solution would be to push everything to `Decimal` and use the precision of `val` to set the needed precision for `1/n`, finally getting the right rounded int from that. – DSM Apr 12 '13 at 19:24
  • @DSM: I am pretty sure that for any given precision setting it would be possible to construct an example that would still fail. – NPE Apr 12 '13 at 19:26
  • @NPE: for any *fixed* precision setting. I'd have to think harder about it than I want to do on a Friday, but for a particular `val`, ISTM you can choose a precision sufficient to get the right answer. – DSM Apr 12 '13 at 19:38
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
Szabolcs Dombi
  • 5,493
  • 3
  • 39
  • 71
user448810
  • 17,381
  • 4
  • 34
  • 59
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
Colonel Panic
  • 132,665
  • 89
  • 401
  • 465
  • 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!) – Peter de Rivaz Apr 12 '13 at 19:29
  • @Peter, so the first guess really can be wrong by more than 1? – Colonel Panic Apr 12 '13 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. – Daniel Fischer Apr 12 '13 at 19:30
  • 1
    int(((10**40)**2)**(1.0/2))=10**40+303786028427003666890752, so it can be off by a long way! – Peter de Rivaz Apr 12 '13 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.

Michele d'Amico
  • 22,111
  • 8
  • 69
  • 76
Natig Aliyev
  • 379
  • 6
  • 18
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
>>> 
Doug Currie
  • 40,708
  • 1
  • 95
  • 119
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
Escualo
  • 40,844
  • 23
  • 87
  • 135
1
def nth_root(n, k):
    x = n**(1./k)
    y = int(x)
    return y + 1 if y != x else y
ecline6
  • 896
  • 8
  • 15
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)

Stochastically
  • 7,616
  • 5
  • 30
  • 58
  • 2
    The problem is with the 1/3, not the 5: 5 is exactly representable by a floating point number. – Alexandre C. Apr 12 '13 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. – Stochastically Apr 12 '13 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 '13 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. – Stochastically Apr 12 '13 at 19:17
  • @Stochastically: For integer inputs it *is* possible to avoid rounding errors. See my answer for one such method. – NPE Apr 12 '13 at 19:19
  • @Stochastically: A well implemented `pow` returns exactly 10 for `pow(100, .5)`. The implementation of a floating-point function is not required to use only normal floating point to produce its answer; internally it may use extended precision or other techniques to ensure a good answer. – Eric Postpischil Apr 12 '13 at 19:25
  • @NPE, just voted for your answer, it's a good answer for integers :-). – Stochastically Apr 13 '13 at 09:55
  • @EricPostpischil I doubt doubt that it's always possible to find better solutions, but however much precision you have, rounding error doesn't go away because you can't have infinite precision. – Stochastically Apr 13 '13 at 10:00
  • @Stochastically: The fact that rounding error is necessary when a mathematical result is not exactly representable is irrelevant to the fact that a well implemented `pow` returns exactly 10 for `pow(100, .5)`. Your assertion that, because logarithms are not exactly representable, results are necessarily inaccurate is false. A properly designed routine could use the natural logarithm of 100, which is not exactly representable, and still get an exactly correct answer even though intermediate arithmetic is inexact. – Eric Postpischil Apr 13 '13 at 12:52
  • @EricPostpischil it sounds like magic, the ability to consistently get an exact result at the end of a calculation, even though an intermediate result is (very slightly) inaccurate. `pow(100,0.5)` is a very simple example of course. For the magic to any use, it would have to be true in all situation where an exact result is possible (because any small subset of situations could be hard coded). So I'm still sceptical, unless you can point me to such an implementation? – Stochastically Apr 13 '13 at 13:19
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
Alexandre C.
  • 55,948
  • 11
  • 128
  • 197
  • 2
    `rtn(4 ** (1./3))` is 2, but `2**3 > 4`. – unutbu Apr 12 '13 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 ! – Alexandre C. Apr 15 '15 at 07:18
0

Do this before everything:

from __future__ import division

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

khan
  • 7,005
  • 15
  • 48
  • 70
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.

ravi tanwar
  • 598
  • 5
  • 16