35

I need a way to compute the nth root of a long integer in Python.

I tried pow(m, 1.0/n), but it doesn't work:

OverflowError: long int too large to convert to float

Any ideas?

By long integer I mean REALLY long integers like:

11968003966030964356885611480383408833172346450467339251 196093144141045683463085291115677488411620264826942334897996389 485046262847265769280883237649461122479734279424416861834396522 819159219215308460065265520143082728303864638821979329804885526 557893649662037092457130509980883789368448042961108430809620626 059287437887495827369474189818588006905358793385574832590121472 680866521970802708379837148646191567765584039175249171110593159 305029014037881475265618958103073425958633163441030267478942720 703134493880117805010891574606323700178176718412858948243785754 898788359757528163558061136758276299059029113119763557411729353 915848889261125855717014320045292143759177464380434854573300054 940683350937992500211758727939459249163046465047204851616590276 724564411037216844005877918224201569391107769029955591465502737 961776799311859881060956465198859727495735498887960494256488224 613682478900505821893815926193600121890632

Improve this question
4
  • As David is implying, pow(n, 1/3) will give you the cubic (i.e. 3rd) root of n. Commented Dec 10, 2008 at 13:59
  • 5
    No it won't, since 1/3 == 0 in python < 3. Commented Dec 10, 2008 at 14:05
  • (But it won't be what the OP wanted, either). Commented Dec 10, 2008 at 14:06
  • Py3 doesn't have integer limitations... they can grow forever until memory runs out. I tested on my installation. That's a solution. Commented Nov 9, 2015 at 4:08

12 Answers 12

Reset to default
30

If it's a REALLY big number. You could use a binary search.

def find_invpow(x,n): """Finds the integer component of the n'th root of x, an integer such that y ** n <= x < (y + 1) ** n. """ high = 1 while high ** n <= x: high *= 2 low = high/2 while low < high: mid = (low + high) // 2 if low < mid and mid**n < x: low = mid elif high > mid and mid**n > x: high = mid else: return mid return mid + 1

For example:

>>> x = 237734537465873465 >>> n = 5 >>> y = find_invpow(x,n) >>> y 2986 >>> y**n <= x <= (y+1)**n True >>> >>> x = 119680039660309643568856114803834088331723464504673392511960931441> >>> n = 45 >>> y = find_invpow(x,n) >>> y 227661383982863143360L >>> y**n <= x < (y+1)**n True >>> find_invpow(y**n,n) == y True >>>
Improve this answer
Sign up to request clarification or add additional context in comments.

4 Comments

It will still fail for numbers > ~10**1000. Changing mid = (low + high) // 2 to mid = int((low + high) // 2) + 1 would fix that.
Downvoted because it has bugs. I tried find_invpow(64, 3) and got 3, even though the result should have been 4.
The line low = high/2 should be low = high // 2 if you actually want an integer out at the end.
Dichotomic search seems a pretty slow method. Prefer Newton's iterations, which double the number of exact digit every time. To find a very accurate starting value, use floating-point.
17

Gmpy is a C-coded Python extension module that wraps the GMP library to provide to Python code fast multiprecision arithmetic (integer, rational, and float), random number generation, advanced number-theoretical functions, and more.

Includes a root function:

x.root(n): returns a 2-element tuple (y,m), such that y is the (possibly truncated) n-th root of x; m, an ordinary Python int, is 1 if the root is exact (x==y**n), else 0. n must be an ordinary Python int, >=0.

For example, 20th root:

>>> import gmpy >>> i0=11968003966030964356885611480383408833172346450467339251 >>> m0=gmpy.mpz(i0) >>> m0 mpz(11968003966030964356885611480383408833172346450467339251L) >>> m0.root(20) (mpz(567), 0)
Improve this answer

3 Comments

With gmpy2 this results in: 'mpz' object has no attribute 'root'.
gmpy2 uses a new mpfr type based on the MPFR library. gmpy2.root(x, n) -> mpfr Return n-th root of x. The result always an 'mpfr'.
@Zelphir gmpy2 has gmpy2.iroot to compute integer roots.
9

If you are looking for something standard, fast to write with high precision. I would use decimal and adjust the precision (getcontext().prec) to at least the length of x.

Code (Python 3.0)

from decimal import * x = '11968003966030964356885611480383408833172346450467339251\ 196093144141045683463085291115677488411620264826942334897996389\ 485046262847265769280883237649461122479734279424416861834396522\ 819159219215308460065265520143082728303864638821979329804885526\ 557893649662037092457130509980883789368448042961108430809620626\ 059287437887495827369474189818588006905358793385574832590121472\ 680866521970802708379837148646191567765584039175249171110593159\ 305029014037881475265618958103073425958633163441030267478942720\ 703134493880117805010891574606323700178176718412858948243785754\ 898788359757528163558061136758276299059029113119763557411729353\ 915848889261125855717014320045292143759177464380434854573300054\ 940683350937992500211758727939459249163046465047204851616590276\ 724564411037216844005877918224201569391107769029955591465502737\ 961776799311859881060956465198859727495735498887960494256488224\ 613682478900505821893815926193600121890632' minprec = 27 if len(x) > minprec: getcontext().prec = len(x) else: getcontext().prec = minprec x = Decimal(x) power = Decimal(1)/Decimal(3) answer = x**power ranswer = answer.quantize(Decimal('1.'), rounding=ROUND_UP) diff = x - ranswer**Decimal(3) if diff == Decimal(0): print("x is the cubic number of", ranswer) else: print("x has a cubic root of ", answer)

Answer

x is the cubic number of 22873918786185635329056863961725521583023133411 451452349318109627653540670761962215971994403670045614485973722724603798 107719978813658857014190047742680490088532895666963698551709978502745901 704433723567548799463129652706705873694274209728785041817619032774248488 2965377218610139128882473918261696612098418

Improve this answer

Comments

8

You can make it run slightly faster by avoiding the while loops in favor of setting low to 10 ** (len(str(x)) / n) and high to low * 10. Probably better is to replace the len(str(x)) with the bitwise length and using a bit shift. Based on my tests, I estimate a 5% speedup from the first and a 25% speedup from the second. If the ints are big enough, this might matter (and the speedups may vary). Don't trust my code without testing it carefully. I did some basic testing but may have missed an edge case. Also, these speedups vary with the number chosen.

If the actual data you're using is much bigger than what you posted here, this change may be worthwhile.

from timeit import Timer def find_invpow(x,n): """Finds the integer component of the n'th root of x, an integer such that y ** n <= x < (y + 1) ** n. """ high = 1 while high ** n < x: high *= 2 low = high/2 while low < high: mid = (low + high) // 2 if low < mid and mid**n < x: low = mid elif high > mid and mid**n > x: high = mid else: return mid return mid + 1 def find_invpowAlt(x,n): """Finds the integer component of the n'th root of x, an integer such that y ** n <= x < (y + 1) ** n. """ low = 10 ** (len(str(x)) / n) high = low * 10 while low < high: mid = (low + high) // 2 if low < mid and mid**n < x: low = mid elif high > mid and mid**n > x: high = mid else: return mid return mid + 1 x = 237734537465873465 n = 5 tests = 10000 print "Norm", Timer('find_invpow(x,n)', 'from __main__ import find_invpow, x,n').timeit(number=tests) print "Alt", Timer('find_invpowAlt(x,n)', 'from __main__ import find_invpowAlt, x,n').timeit(number=tests)

Norm 0.626754999161

Alt 0.566340923309

Improve this answer

3 Comments

Your second function, find_invpowAlt, gives a wildly wrong answer for x=118997879821732370764604711647724283139870175351576755860556891902958645241483485254092600557474860904935286687480039428945219115513349647465379580432922136155355040992635166676363150438436216219094913514982415747153956476970303302126880391024128871557664284712411567099374094385902892603751471822837746770111, n=3
The line low = high/2 should be low = high // 2 if you actually want an integer out at the end.
I wrote it using low = high/2 intentionally, as I directly copied that answer from Markus's Code as a means of benchmarking it when compared to my proposed improvement. I'll note that since then, Markus has updated his code to correct a bug and tzs has replied to my code that it also has a bug. However, I don't actually remember the details of the problem nor of my solution, so I no longer feel competent to address these shortcomings.
4

I may suggest four methods for solving your task. First is based on Binary Search. Second is based on Newton's Method. Third is based on Shifting n-th Root Algorithm. Fourth is called by me Chord-Tangent method described by me in picture here.

Binary Search was already implemented in many answers above. I just introduce here my own vision of it and its implementation.

As alternative I also implement Optimized Binary Search method (marked Opt). This method just starts from range [hi / 2, hi) where hi is equal to 2^(num_bit_length / k) if we're computing k-th root.

Newton's Method is new here, as I see it wasn't implemented in other answers. It is usually considered to be faster than Binary Search, although my own timings in code below don't show any speedup. Hence this method here is just for reference/interest.

Shifting Method is 30-50% faster than optimized binary search method, and should be even faster if implemented in C++, because C++ has fast 64 bit arithemtics which is partially used in this method.

Chord-Tangent Method:

enter image description here

Chord-Tangent Method is invented by me on piece of paper (see image above), it is inspired and is an improvement of Newton method. Basically I draw a Chord and a Tangent Line and find intersection with horizontal line y = n, these two intersections form lower and upper bound approximations of location of root solution (x0, n) where n = x0 ^ k. This method appeared to be fastest of all, while all other methods do more than 2000 iterations, this method does just 8 iterations, for the case of 8192-bit numbers. So this method is 200-300x times faster than previous (by speed) Shifting Method.

As an example I generate really huge random integer of 8192 bits in size. And measure timings of finding cubic root with both methods.

In test() function you can see that I passed k = 3 as root's power (cubic root), you can pass any power instead of 3.

Try it online!

def binary_search(begin, end, f, *, niter = [0]): while begin < end: niter[0] += 1 mid = (begin + end) >> 1 if f(mid): begin = mid + 1 else: end = mid return begin def binary_search_kth_root(n, k, *, verbose = False): # https://en.wikipedia.org/wiki/Binary_search_algorithm niter = [0] res = binary_search(0, n + 1, lambda root: root ** k < n, niter = niter) if verbose: print('Binary Search iterations:', niter[0]) return res def binary_search_opt_kth_root(n, k, *, verbose = False): # https://en.wikipedia.org/wiki/Binary_search_algorithm niter = [0] hi = 1 << (n.bit_length() // k - 1) while hi ** k <= n: niter[0] += 1 hi <<= 1 res = binary_search(hi >> 1, hi, lambda root: root ** k < n, niter = niter) if verbose: print('Binary Search Opt iterations:', niter[0]) return res def newton_kth_root(n, k, *, verbose = False): # https://en.wikipedia.org/wiki/Newton%27s_method f = lambda x: x ** k - n df = lambda x: k * x ** (k - 1) x, px, niter = n, 2 * n, [0] while abs(px - x) > 1: niter[0] += 1 px = x x -= f(x) // df(x) if verbose: print('Newton Method iterations:', niter[0]) mini, minv = None, None for i in range(-2, 3): v = abs(f(x + i)) if minv is None or v < minv: mini, minv = i, v return x + mini def shifting_kth_root(n, k, *, verbose = False): # https://en.wikipedia.org/wiki/Shifting_nth_root_algorithm B_bits = 64 r, y = 0, 0 B = 1 << B_bits Bk_bits = B_bits * k Bk_mask = (1 << Bk_bits) - 1 niter = [0] for i in range((n.bit_length() + Bk_bits - 1) // Bk_bits - 1, -1, -1): alpha = (n >> (i * Bk_bits)) & Bk_mask B_y = y << B_bits Bk_yk = (y ** k) << Bk_bits Bk_r_alpha = (r << Bk_bits) + alpha Bk_yk_Bk_r_alpha = Bk_yk + Bk_r_alpha beta = binary_search(1, B, lambda beta: (B_y + beta) ** k <= Bk_yk_Bk_r_alpha, niter = niter) - 1 y, r = B_y + beta, Bk_r_alpha - ((B_y + beta) ** k - Bk_yk) if verbose: print('Shifting Method iterations:', niter[0]) return y def chord_tangent_kth_root(n, k, *, verbose = False): niter = [0] hi = 1 << (n.bit_length() // k - 1) while hi ** k <= n: niter[0] += 1 hi <<= 1 f = lambda x: x ** k df = lambda x: k * x ** (k - 1) # https://i.sstatic.net/et9O0.jpg x_begin, x_end = hi >> 1, hi y_begin, y_end = f(x_begin), f(x_end) for icycle in range(1 << 30): if x_end - x_begin <= 1: break niter[0] += 1 if 0: # Do Binary Search step if needed x_mid = (x_begin + x_end) >> 1 y_mid = f(x_mid) if y_mid > n: x_end, y_end = x_mid, y_mid else: x_begin, y_begin = x_mid, y_mid # (y_end - y_begin) / (x_end - x_begin) = (n - y_begin) / (x_n - x_begin) -> x_n = x_begin + (n - y_begin) * (x_end - x_begin) // (y_end - y_begin) y_n = f(x_n) tangent_x = x_n + (n - y_n) // df(x_n) + 1 chord_x = x_n + (n - y_n) * (x_end - x_n) // (y_end - y_n) assert chord_x <= tangent_x, (chord_x, tangent_x) x_begin, x_end = chord_x, tangent_x y_begin, y_end = f(x_begin), f(x_end) assert y_begin <= n, (chord_x, y_begin, n, n - y_begin) assert y_end > n, (icycle, tangent_x - binary_search_kth_root(n, k), y_end, n, y_end - n) if verbose: print('Chord Tangent Method iterations:', niter[0]) return x_begin def test(): import random, timeit nruns = 3 bits = 8192 n = random.randrange(1 << (bits - 1), 1 << bits) a = binary_search_kth_root(n, 3, verbose = True) b = binary_search_opt_kth_root(n, 3, verbose = True) c = newton_kth_root(n, 3, verbose = True) d = shifting_kth_root(n, 3, verbose = True) e = chord_tangent_kth_root(n, 3, verbose = True) assert abs(a - b) <= 0 and abs(a - c) <= 1 and abs(a - d) <= 1 and abs(a - e) <= 1, (a - b, a - c, a - d, a - e) print() print('Binary Search timing:', round(timeit.timeit(lambda: binary_search_kth_root(n, 3), number = nruns) / nruns, 3), 'sec') print('Binary Search Opt timing:', round(timeit.timeit(lambda: binary_search_opt_kth_root(n, 3), number = nruns) / nruns, 3), 'sec') print('Newton Method timing:', round(timeit.timeit(lambda: newton_kth_root(n, 3), number = nruns) / nruns, 3), 'sec') print('Shifting Method timing:', round(timeit.timeit(lambda: shifting_kth_root(n, 3), number = nruns) / nruns, 3), 'sec') print('Chord Tangent Method timing:', round(timeit.timeit(lambda: chord_tangent_kth_root(n, 3), number = nruns) / nruns, 3), 'sec') if __name__ == '__main__': test()

Output:

Binary Search iterations: 8192 Binary Search Opt iterations: 2732 Newton Method iterations: 9348 Shifting Method iterations: 2752 Chord Tangent Method iterations: 8 Binary Search timing: 0.506 sec Binary Search Opt timing: 0.05 sec Newton Method timing: 2.09 sec Shifting Method timing: 0.03 sec Chord Tangent Method timing: 0.001 sec
Improve this answer

1 Comment

Why are you assigning niter a default value?
3

Oh, for numbers that big, you would use the decimal module.

ns: your number as a string

ns = "11968003966030964356885611480383408833172346450467339251196093144141045683463085291115677488411620264826942334897996389485046262847265769280883237649461122479734279424416861834396522819159219215308460065265520143082728303864638821979329804885526557893649662037092457130509980883789368448042961108430809620626059287437887495827369474189818588006905358793385574832590121472680866521970802708379837148646191567765584039175249171110593159305029014037881475265618958103073425958633163441030267478942720703134493880117805010891574606323700178176718412858948243785754898788359757528163558061136758276299059029113119763557411729353915848889261125855717014320045292143759177464380434854573300054940683350937992500211758727939459249163046465047204851616590276724564411037216844005877918224201569391107769029955591465502737961776799311859881060956465198859727495735498887960494256488224613682478900505821893815926193600121890632" from decimal import Decimal d = Decimal(ns) one_third = Decimal("0.3333333333333333") print d ** one_third

and the answer is: 2.287391878618402702753613056E+305

TZ pointed out that this isn't accurate... and he's right. Here's my test.

from decimal import Decimal def nth_root(num_decimal, n_integer): exponent = Decimal("1.0") / Decimal(n_integer) return num_decimal ** exponent def test(): ns = "11968003966030964356885611480383408833172346450467339251196093144141045683463085291115677488411620264826942334897996389485046262847265769280883237649461122479734279424416861834396522819159219215308460065265520143082728303864638821979329804885526557893649662037092457130509980883789368448042961108430809620626059287437887495827369474189818588006905358793385574832590121472680866521970802708379837148646191567765584039175249171110593159305029014037881475265618958103073425958633163441030267478942720703134493880117805010891574606323700178176718412858948243785754898788359757528163558061136758276299059029113119763557411729353915848889261125855717014320045292143759177464380434854573300054940683350937992500211758727939459249163046465047204851616590276724564411037216844005877918224201569391107769029955591465502737961776799311859881060956465198859727495735498887960494256488224613682478900505821893815926193600121890632" nd = Decimal(ns) cube_root = nth_root(nd, 3) print (cube_root ** Decimal("3.0")) - nd if __name__ == "__main__": test()

It's off by about 10**891

Improve this answer

3 Comments

Um. This could work, but it's not accurate. Using your terms, and if answer equals done_third, then how much (answer3 - d) should be?
Decimal is about as accurate as you need it to be... my 0.333 string is just for brevity.
tz... you're correct. It' way off... oh well. The newton's method above does rock though!
2

Possibly for your curiosity:

http://en.wikipedia.org/wiki/Hensel_Lifting

This could be the technique that Maple would use to actually find the nth root of large numbers.

Pose the fact that x^n - 11968003.... = 0 mod p, and go from there...

Improve this answer

Comments

1

I came up with my own answer, which takes @Mahmoud Kassem's idea, simplifies the code, and makes it more reusable:

def cube_root(x): return decimal.Decimal(x) ** (decimal.Decimal(1) / decimal.Decimal(3))

I tested it in Python 3.5.1 and Python 2.7.8, and it seemed to work fine.

The result will have as many digits as specified by the decimal context the function is run in, which by default is 28 decimal places. According to the documentation for the power function in the decimal module, "The result is well-defined but only “almost always correctly-rounded”.". If you need a more accurate result, it can be done as follows:

with decimal.localcontext() as context: context.prec = 50 print(cube_root(42))
Improve this answer

Comments

1

You could use SymPy for this:

In [29]: n = 1196800396603096435688561148038340883317234645046733925119609314414104568346308529111567 ...: 74884116202648269423348979963894850462628472657692808832376494611224797342794244168618343965 ...: 22819159219215308460065265520143082728303864638821979329804885526557893649662037092457130509 ...: 98088378936844804296110843080962062605928743788749582736947418981858800690535879338557483259 ...: 01214726808665219708027083798371486461915677655840391752491711105931593050290140378814752656 ...: 18958103073425958633163441030267478942720703134493880117805010891574606323700178176718412858 ...: 94824378575489878835975752816355806113675827629905902911311976355741172935391584888926112585 ...: 57170143200452921437591774643804348545733000549406833509379925002117587279394592491630464650 ...: 47204851616590276724564411037216844005877918224201569391107769029955591465502737961776799311 ...: 85988106095646519885972749573549888796049425648822461368247890050582189381592619360012189063 ...: 2 In [30]: import sympy In [31]: r = sympy.Integer(n) ** sympy.Rational(1, 3) In [32]: print(r) 228739187861856353290568639617255215830231334114514523493181096276535406707619622159719944036700456144859737227246037981077199788136588570141900477426804900885328956669636985517099785027459017044337235675487994631296527067058736942742097287850418176190327742484882965377218610139128882473918261696612098418

This is just using SymPy's general symbolic expression simplification. The above will give an exact result if possible or otherwise an expression e.g.:

In [33]: r2 = sympy.Integer(10) ** sympy.Rational(1, 3) In [34]: r2 Out[34]: 3 ____ ╲╱ 10 In [35]: r2.evalf(50) Out[35]: 2.1544346900318837217592935665193504952593449421921

SymPy also has a function specifically for computing the nth root of an integer as an integer:

In [38]: sympy.integer_nthroot(n, 3) Out[38]: (228739187861856353290568639617255215830231334114514523493181096276535406707619622159719944036700456144859737227246037981077199788136588570141900477426804900885328956669636985517099785027459017044337235675487994631296527067058736942742097287850418176190327742484882965377218610139128882473918261696612098418, True)

The True here indicates that the result is an exact cube root. When an exact root is not possible this function will return the floor of the cube root:

In [39]: sympy.integer_nthroot(10, 3) Out[39]: (2, False)

In other words 2 is the greatest integer such that 2**3 <= 10 but the False indicates that 2**3 != 10.

Improve this answer

Comments

0

In older versions of Python, 1/3 is equal to 0. In Python 3.0, 1/3 is equal to 0.33333333333 (and 1//3 is equal to 0).

So, either change your code to use 1/3.0 or switch to Python 3.0 .

Improve this answer

2 Comments

As far as I can tell, there was no indication that it used Python 2.
@SolomonUckoI This question and answer were both posted only 1 week after python 3 was released.
-1

Try converting the exponent to a floating number, as the default behaviour of / in Python is integer division

n**(1/float(3))

Improve this answer

3 Comments

n**(1.0/3) will do the job too
OverflowError: long int too large to convert to float, it just would not work cause the huge number cannot be converted
As far as I can tell, there was no indication that it used Python 2.
-3

Well, if you're not particularly worried about precision, you could convert it to a sting, chop off some digits, use the exponent function, and then multiply the result by the root of how much you chopped off.

E.g. 32123 is about equal to 32 * 1000, the cubic root is about equak to cubic root of 32 * cubic root of 1000. The latter can be calculated by dividing the number of 0s by 3.

This avoids the need for the use of extension modules.

Improve this answer

4 Comments

I'm worried about precision because i know this number, for example is the third power of another integer. (And of course i need to know this integer) :)
Also, converting to string might become problematic with huge numbers.
@Attila: Not in the definition of huge he gave in the post. Also, not sure exactly what format he is accepting input in, but the basic idea of truncation will work in most formats.
He has working code, just the problem is precision

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.