Talk:Nth root algorithm

I've been playing around with finding (integer) nth roots for large n. Unfortunately, the following implementation of Newton's method (in Python) is ridiculously slow:

def nthroot(y, n):
    x, xp = 1, -1
    while abs(x - xp) > 1:
        xp, x = x, x - x/n + y/(n * x**(n-1))
    while x**n > y:
        x -= 1
    return x

For example, nthroot(12345678901234567890123456789012345**100,100) takes several minutes to finish.

Using binary search, the solution can be found in seconds.

Am I just doing something wrong, or is Newton's method inherently inefficient for integers and such large n? - Fredrik | talk 23:22, 24 May 2005 (UTC)

Brief answer (let me know if you need more information): Are you sure that the standard datatype of Python can handle 125-digit integers? I think you need a special library for this. Besides, it may be necessary to start the iteration with a more sensible initial guess than x = 1 (generally, Newton's method may not converge unless the initial guess is close to the answer, but extracting N-th roots may be a special case). -- Jitse Niesen 19:39, 6 Jun 2005 (UTC)

No, no special library is required. I have tried various different initial guesses, but they offer no improvement. Fredrik | talk 20:55, 6 Jun 2005 (UTC)

I thought a bit more about it now. Basically, you are right: Newton's method is inefficient for large n. There are two reasons for it. Firstly, the computation of x**(n-1) in exact arithmetics is expensive if x is a large integer, so every iteration takes a long time. Secondly, the method needs a lot of iterations. For instance, nthroot(1234567890123456789012345**10, 10) requires 4725 iterations (note that I replaced 100 in your example with 10), while a binary search requires about 100 iterations. However, the method improves if you give a good initial guess. After replacing the first line of the nthroot function with x, xp = int(y**(1.0/n)), -1, the above example gets the correct answer in only 3 iterations.

The current article is rather bad: it only presents the algorithm, which is trivial if you know Newton's method. There is no analysis, no discussion on how to get the initial guess and no references. I am not convinced that the algorithm is used that often; I thought the common approach is to use exponential and logarithms. Bignum arithmetics is a separate issue. I hope you can find the time to improve the article. -- Jitse Niesen 23:37, 7 Jun 2005 (UTC)

The catch with that solution is that y**(1.0/n) overflows for large y. Fortunately this can be avoided by using B**int(log(y, B)/n) with some integer B instead. However, this still leaves Newton (as implemented above) significantly slower than bsearch for very large n (a test with n=300 is nearly instantaneous with bsearch, but takes several seconds with Newton). I'll do some more detailed benchmarking when I get time, hopefully soon. For reference, I'll provide my bsearch implementation here: - Fredrik | talk 00:24, 8 Jun 2005 (UTC)

# Returns a tuple of the (floor, ceil) values of the exact solution
def nthroot_bsearch(x, n):
    guess = 1 # Initial guess must be less than the result
    step = 1
    while 1:
        w = (guess+step)**n
        if w == x:
            return (guess+step,) * 2
        elif w < x:
            step <<= 1
        elif step == 1:
            return guess, guess+1
        else:
            guess += step >> 1
            step = 1