4

I was working with numbers of 200 digits in python. When finding the square root of a number using math.sqrt(n) I am getting a wrong answer.

In[1]: n=9999999999999999999999999999999999999999999999999999999999999999999999
       999999999999999999999999998292000000000000000000000000000000000000000000
       0000000000000000000000000000000000000000000000000000726067

In[2]: x=int(math.sqrt(n))

In[3]: x
Out[1]: 10000000000000000159028911097599180468360808563945281389781327
        557747838772170381060813469985856815104L

In[4]: x*x
Out[2]: 1000000000000000031805782219519836346574107361670094060730052612580
        0264077231077619856175974095677538298443892851483731336069235827852
        3336313169161345893842466001164011496325176947445331439002442530816L

In[5]: math.sqrt(n)
Out[3]: 1e+100

The value of x is coming larger than expected since x*x (201 digits) is larger than n (200 digits). What is happening here? Is there some concept I am getting wrong here? How else can I find the root of very large numbers?

Juan John Mathews
  • 728
  • 1
  • 10
  • 24
  • 1
    change `x=int(math.sqrt(n))` to `x=math.sqrt(n)` and try it – Sailesh Kotha Jan 26 '15 at 13:20
  • Its gives the same answer in exponential form as 1e+100. – Juan John Mathews Jan 26 '15 at 13:32
  • 3
    `math.sqrt()` returns a floating-point value that is subject to IEEE 754 arithmetic, not a bignum. – Frédéric Hamidi Jan 26 '15 at 13:38
  • 1
    http://stackoverflow.com/questions/15390807/integer-square-root-in-python has useful answers. – georg Jan 26 '15 at 13:41
  • 2
    To handle such large numbers you should look at the standard `decimal` module, or the 3rd-party `mpmath` module. Alternatively, if you just want an integer square root, you can write one using the [Babylonian method](http://en.wikipedia.org/wiki/Methods_of_computing_square_roots), aka Hero's method, which is a special case of Newton's method. – PM 2Ring Jan 26 '15 at 13:43

3 Answers3

9

Using the decimal module:

import decimal
D = decimal.Decimal
n = D(99999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999982920000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000726067)

with decimal.localcontext() as ctx:
    ctx.prec = 300
    x = n.sqrt()
    print(x)
    print(x*x)
    print(n-x*x)

yields

9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999145.99999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999983754999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999998612677

99999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999982920000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000726067.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

0E-100
unutbu
  • 842,883
  • 184
  • 1,785
  • 1,677
6

math.sqrt returns an IEEE-754 64-bit result, which is roughly 17 digits. There are other libraries that will work with high-precision values. In addition to the decimal and mpmath libraries mentioned above, I maintain the gmpy2 library (https://code.google.com/p/gmpy/).

>>> import gmpy2
>>> n=gmpy2.mpz(99999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999982920000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000726067)
>>> gmpy2.get_context().precision=2048
>>> x=gmpy2.sqrt(n)
>>> x*x
mpfr('99999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999982920000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000726067.0',2048)
>>>

The gmpy2 library can also return integer square roots (isqrt) or quickly check if an integer is an exact square (is_square).

Frédéric Hamidi
  • 258,201
  • 41
  • 486
  • 479
casevh
  • 11,093
  • 1
  • 24
  • 35
1

Here's an integer square root program using Hero's method that I wrote a while back. For the initial approximation it uses a number of half the bit length of the input value, so it starts converging pretty quickly. However, I haven't timed it to see if it's faster in Python than just using a simpler initial approximation. :)

#! /usr/bin/env python

''' Long integer square roots. Newton's method.

    Written by PM 2Ring. Adapted from C to Python 2008.10.19
'''

import sys

def root(m):
    # Get initial approximation
    n, a, k = m, 1, 0
    while n > a:
        n >>= 1
        a <<= 1
        k += 1
        #print k, ':', n, a

    # Go back one step & average
    a = n + (a>>2)
    #print a

    # Apply Newton's method
    while k:
        a = (a + m // a) >> 1
        k >>= 1
        #print k, ':', a
    return a

def main():
    m = len(sys.argv) > 1 and int(sys.argv[1]) or 2*10L**100
    print "The Square Root of", m
    print root(m)

if __name__ == '__main__':
    main()
PM 2Ring
  • 54,345
  • 6
  • 82
  • 182