Python sqrt limit for very large numbers?

Question

I'm working with very large numbers (1,000,000 digits) and I need to calculate their square root. I seem to be hitting on a limit in my code.

y = 10**309
x = y**0.5
print(x)

And I'm getting this error:

x = y**0.5
OverflowError: int too large to convert to float

The code works till 10**308. But beyond that it seems broken. I've checked this in command line as well. Same error. Can someone please help me?

If this is a Python limit, is there an alternate method I could use?

Answer 1

Simplifiy your problem, using a bit of math.

Note that sqrt(a*b) = sqrt(a) * sqrt(b) (for real, positive numbers at least).

So, any number larger than, say, 10^100, divide by 10^100. That's a, and the result of the division is b, so that your original number = a * b. Then use the square root of 10^100 (= 10^50), multiply that by the square root of b, and you have your answer.

With your example:

import math
x = 10**309
a = 1e100
b = 1e209   # Note: you can't calculate this within Python; just use plain math here
y = 1e50 * math.sqrt(1e209)

Example for a not-so-round number:

x = 3.1415 * 1e309
a = 1e100
b = 3.1415e209   # Again, just subtract the exponent: 309 - 100
y = 1e50 * math.sqrt(3.1415e209)

Or for an integer that's not a power of 10, fully written out:

x = 707070
x = 70.707 * 1e4  # note: even number in exponent
x = 70.707e4
a = 1e2  # sqrt(1e2) = 1e1 = 10
b = 70.707e2
y = 10 * sqrt(70.707e2)

A few notes:

• Python handles ridiculously large integer numbers without problems. For floating point numbers, it uses standard (C) conventions, and limits itself to 64 bit precision. You almost always get floating point numbers when taking a square root of something.

• 1e309 means 10**309, and 3.1415e209 means 3.1415 * 10**209. This is a standard programming convention.

Answer 2

You should use the gmpy2 module. It provides very fast multiple-precision arithmetic.

On my system, operations on million digit numbers are very fast.

In [8]: a=gmpy2.mpz('3'*1000000)

In [9]: %timeit gmpy2.isqrt(a)
10 loops, best of 3: 22.8 ms per loop

In [10]: %timeit (a+1)*(a-1)
10 loops, best of 3: 20.9 ms per loop

Working with 100,000,000 digit numbers only takes a few seconds.

In [20]: a.num_digits(10)
Out[20]: 99995229

In [21]: %timeit gmpy2.isqrt(a)
1 loops, best of 3: 5.05 s per loop

In [22]: %timeit (a+1)*(a-1)
1 loops, best of 3: 3.49 s per loop

Disclaimer: I'm the current maintainer of gmpy2.

Answer 3

Based on what I think are similar questions, you can look into using the Decimal class.

Here is an example using what you have

>>> x = 10**309
>>> y =x**.5
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
OverflowError: long int too large to convert to float
>>> import decimal
>>> d = decimal.Decimal(x)
>>> d.sqrt()
Decimal('3.162277660168379331998893544E+154')
>>> float(d.sqrt())
3.1622776601683792e+154

May be helpful, it doesn't give you your error.