Calculate Nth root with integer arithmetic

Question

There are a couple of ways to find integer square roots using only integer arithmetic. For example this one. It makes for interesting reading and also a very interesting theory, particularly for my generation where such techniques aren't so useful any more.

The main thing is that it can't use floating point arithmetic, so that rules out newtons method and it's derivations. The only other way I know of to find roots is binomial expansion, but that also requires floating point arithmetic.

What techniques/algorithms are there for computing integral nth roots using only integer arithmetic?

Edit: Thanks for all the answers so far. They all seem to be slightly more intelligent trial and improvement. Is there no better way?

Edit2: Ok, so it would seem there is no smart way to do this without trial/improvement and either newtons method or a binary search. Can anyone provide a comparison of the two in theory? I have run a number of benchmarks between the two and found them quite similar.

Zelig

What is your required range of input values ? – Leesen

@PaulR, Ideally it could be extensible, but I think you can assume both the base and the number will be 32 bit (unsigned) integers for now. – Zelig

Which integer operations are you permitting? Square roots are a special case because it's possible to extract them using just addition, subtraction and shifts. – Lexis

@Neil, I don't want to place restrictions on it, as this is not for a particular application, but I would say a list similar to say the C list of integer operators: addition, subtraction, multiplication, (integer) division, modulo and bitwise operations. Ofc speed is always a consideration, but don't worry about it too much. – Zelig

In general: there is nothing in the world that you can do with floating point arithmetic and can't do with integer arithmetic almost the same way. At least because floating point arithmetic itself is implementable pretty easy via integer arithmetic. – Ligni

@Mat In that case, I would have gone with the shifting Nth root algorithm as per AakashM's answer. – Lexis

@Neil, I am busy trying to implement that algorithm, to compare it to my existing solution for Newtons method. Unfortunately there isn't much code online for that, – Zelig

@Neil, Read my comment I just added for AakashM's answer, as it shows why I am less inclined to use such a method. – Zelig

@Mat Normally on a computer I would expect you would use 2 as the base, which would appear to avoid the problem. – Lexis

@Neil, it's not a problem, it's just that when working with base 2 the whole algorithm effectively becomes a binary search, same as the other answers to this question. – Zelig

@Zelig — The shifting-nth-root algorithm effectively becomes a binary search only when the base is greater than the radicand (for example computing the cube root of 97 in base 100). In base 2, the base is never greater than the radicand except in the degenerate case where the radicand is 1. In any case, the shifting-nth-root algorithm is extremely efficient in base 2. – V1

Answer 1

You can use Newton's method using only integer arithmetic, the step is the same as for floating point arithmetic, except you have to replace floating point operators with the corresponding integer operators in languages which have different operators for these.

Let's say you want to find the integer-k-th root of a > 0, which should be the largest integer r such that r^k <= a. You start with any positive integer (of course a good starting point helps).

int_type step(int_type k, int_type a, int_type x) {
    return ((k-1)*x + a/x^(k-1))/k;
}

int_type root(int_type k, int_type a) {
    int_type x = 1, y = step(k,a,x);
    do {
        x = y;
        y = step(k,a,x);
    }while(y < x);
    return x;
}

Except for the very first step, you have x == r <==> step(k,a,x) >= x.

Fundament

After looking again at newton raphson, I found there was a way to do it, but it very often reached a point where is flipped between two numbers (e.g. square root of 15 flips between 3 and 4). To counter for this the full solution is here – Zelig

For square roots, it's flipping precisely for a = n*n-1 between n-1 and n. I'm not sure if for higher powers there are more values that lead to flipping, but whenever the step increases the approximation to the root - excepting the very first step, if the starting point was smaller than the target - you're done, the smaller value is the integer root. – Fundament

That is the same conclusion I reached, which is why I arrived at the code posted in my above comment. Regardless of base the values that flip seem to always be above and below the root, so the root is inbetween those two numbers (hence why it flips) my code deals with that. – Zelig

Ok, so it would appear the flipping is far more complex than I previously imagined (take the cube root of 7 and it flips between 1, 2 and 3. From that I can only imagine that the flipping is between n possible numbers for the nth root of A. This adds a lot of complexity to the algorithm. – Zelig

But it cycles 1 -> 3 -> 2 -> 1, 1 is the solution. You only need to check whether step(x) < x. If it is smaller, go on, otherwise, x is the solution (with the exception of the first step, if you start with 5 for the cube root of 511, you get 5 -> 10 -> 8 -> 7 -> 8; the first is uninteresting, after that, any time step(x) >= x, you are done. – Fundament

Ok, I have this sorted. However, it would appear overflow is going to be a huge problem as you reach even reasonable sized numbers due to the use of x^(n-1) as that will result in overflow at even small values. No doubt this is a limitation of any solution using powers, but surely there is a way that does not require using very large values. If not then this solution is very limited, even with low exponents. – Zelig

@Daniel: I tested several formulas and this one seems to be the fastest and the most accurate, I've one question though: how can the step() equation be rewritten to allow k floats (not of the form 1 / n)? For instance, pow(2, 0.12345) = 1.089337... however, if I set k to pow(0.12345, -1) = 8.100446... I would get 1.090508... (assuming the function could accept floats). – Kippar

@Zelig your link had some issues handling 0 and 1. I have modified to remove those issues. Have a look: modified one – Melisamelisande

Answer 2

One obvious way would be to use binary search together with exponentiation by squaring. This will allow you to find nthRoot(x, n) in O(log (x + n)): binary search in [0, x] for the largest integer k such that k^n <= x. For some k, if k^n <= x, reduce the search to [k + 1, x], otherwise reduce it to [0, k].

Do you require something smarter or faster?

Garvin

I am interested to see if there are any methods that do not involve trial an improvement. Though exponentiation by squaring is a good find thanks, – Zelig

Answer 3

It seems to me that the Shifting nth root algorithm provides exactly what you want:

The shifting nth root algorithm is an algorithm for extracting the nth root of a positive real number which proceeds iteratively by shifting in n digits of the radicand, starting with the most significant, and produces one digit of the root on each iteration, in a manner similar to long division.

There are worked examples on the linked wikipedia page.

Brahear

From the wikipedia page: "When the base is larger than the radicand, the algorithm degenerates to binary search". I will have a look if possibly working with (effectively) hex rather than binary will improve the algorithm. – Zelig

Answer 4

One easy solution is to use the binary search.

Assume we are finding nth root of x.

Function GetRange(x,n):
    y=1
    While y^n < x:
        y*2
    return (y/2,y)

Function BinSearch(a,b,x,):
    if a == b+1:
        if x-a^n < b^n - x:
           return a
        else:
           return b
    c = (a+b)/2
    if n< c^n:
        return BinSearch(a,c,x,n)
    else:
        return BinSearch(c,b,x,n)

a,b = GetRange(x,n)
print BinSearch(a,b,x,n)

===Faster Version===

Function BinSearch(a,b,x,):
    w1 = x-a^n
    w2 = b^n - x
    if a <= b+1:
        if w1 < w2:
           return a
        else:
           return b
    c = (w2*a+w1*b)/(w1+w2)
    if n< c^n:
        return BinSearch(a,c,x,n)
    else:
        return BinSearch(c,b,x,n)

Zoochore

Answer 5

I made the algorithm in VBA in Excel. For now it only calculates roots of integers. It is easy to implement the decimals as well.

Just copy and paste the code into an EXCEL module and type the name of the function into some cell, passing the parameters.

Public Function RootShift(ByVal radicand As Double, degree As Long, Optional ByRef remainder As Double = 0) As Double

   Dim fullRadicand As String, partialRadicand As String, missingZeroes As Long, digit As Long

   Dim minimalPotency As Double, minimalRemainder As Double, potency As Double

   radicand = Int(radicand)

   degree = Abs(degree)

   fullRadicand = CStr(radicand)

   missingZeroes = degree - Len(fullRadicand) Mod degree

   If missingZeroes < degree Then

      fullRadicand = String(missingZeroes, "0") + fullRadicand

   End If

   remainder = 0

   RootShift = 0

   Do While fullRadicand <> ""

      partialRadicand = Left(fullRadicand, degree)

      fullRadicand = Mid(fullRadicand, degree + 1)

      minimalPotency = (RootShift * 10) ^ degree

      minimalRemainder = remainder * 10 ^ degree + Val(partialRadicand)

      For digit = 9 To 0 Step -1

          potency = (RootShift * 10 + digit) ^ degree - minimalPotency

          If potency <= minimalRemainder Then

             Exit For

          End If

      Next

      RootShift = RootShift * 10 + digit

      remainder = minimalRemainder - potency

   Loop

End Function

Palmation

Answer 6

Algorithm more simple in VBA.

Public Function RootNth(radicand As Double, degree As Long) As Double
   Dim countDigits As Long, digit As Long, potency As Double
   Dim minDigit As Long, maxDigit As Long, partialRadicand As String
   Dim totalRadicand As String, remainder As Double

  radicand = Int(radicand)
  degree = Abs(degree)
  RootNth = 0
  partialRadicand = ""
  totalRadicand = CStr(radicand)
  countDigits = Len(totalRadicand) Mod degree
  countDigits = IIf(countDigits = 0, degree, countDigits)
  Do While totalRadicand <> ""
     partialRadicand = partialRadicand + Left(totalRadicand, countDigits)
     totalRadicand = Mid(totalRadicand, countDigits + 1)
     countDigits = degree
     minDigit = 0
     maxDigit = 9
     Do While minDigit <= maxDigit
        digit = Int((minDigit + maxDigit) / 2)
        potency = (RootNth * 10 + digit) ^ degree
        If potency = Val(partialRadicand) Then
           maxDigit = digit
           Exit Do
        End If
        If potency < Val(partialRadicand) Then
           minDigit = digit + 1
        Else
           maxDigit = digit - 1
        End If
     Loop
     RootNth = RootNth * 10 + maxDigit
  Loop
   End Function

Okubo

"More simple" than what? – Inbreathe