Limit of digit-by-digit calculation of square roots

Question

I am trying to get as good an estimate of pi as I can using the Chudnovsky algorithm in Python. This algorithm implies getting the square root of 640 320.

After doing some research, I found a quite effective way to compute square roots; the method is called "digit-by-digit calculation" (see here). So, after trying to implement it, I found out that the first 13 decimals are correct, and then I get strange results (the next one is a 0 instead of a 4, and then the next 'digit' is 128, then -1024...)

I tried checking my function, but it looks fine to me (besides, I would probably not find the correct first 13 decimals otherwise). Thus, my question is : are there some limits in this digit-by-digit calculation method?

If, by any chance, you would like to see my code, here it is:

def sqrt(input_number,accuracy):
"""input_number is a list that represents a number we want to get the square root of.
For example, 12.56 would be [[1,2], [5,6], '+']"""
    
    if input_number[2]!="+":
        raise ValueError("Cannot find the real square root of a negative number: '"+pl(input_number)+"'")
    
"""Below creates the right number of elements for the required accuracy of the
square root"""
    if len(input_number[0])%2==1:
        input_number[0].insert(0,0)
        
    if len(input_number[1])<2*accuracy:
        for i in range(2*accuracy-len(input_number[1])):
            input_number[1].append(0)
            
    if len(input_number[1])%2==1:
        input_number[1].append(0)
    
# Below makes the pairs of digits required in the algorithm
    pairs=[[10*input_number[0][2*i]+input_number[0][2*i+1] for i in range(int(len(input_number[0])/2))],[10*input_number[1][2*i]+input_number[1][2*i+1] for i in range(int(len(input_number[1])/2))]]

    
"""Performs the algorithm, where c,x,y and p have the same definition
as on the Wikipedia link above. r is the remainder. pairs[0] is the pairs
of digits before the decimal dot, and pairs[1] represents the pairs of
digits after the dot. square_root is the computed square root of input_number."""
    p=0
    r=0
    square_root=[[],[],"+"]
    for i in range(len(pairs[0])):
        c=100*r+pairs[0][i]
        x=int((-20*p+(400*p**2+4*c)**.5)/2)
        y=20*p*x+x**2
        r=c-y
        p=10*p+x
        square_root[0].append(x)
    
    for i in range(len(pairs[1])):
        print(p,r,c)
        c=100*r+pairs[1][i]
        x=int((-20*p+(400*p**2+4*c)**.5)/2)
        y=20*p*x+x**2
        r=c-y
        p=10*p+x
        square_root[1].append(x)
    
    
    return square_root

Why not comment your code, or at least use descriptive variable names? — Douglas Zare, Apr 19 '15 at 01:55
Yes, you're totally right. I just added more comments and changed some variable names. — AnthonyD973, Apr 19 '15 at 02:10

0xrgb · Accepted Answer · 2015-04-19T07:40:45.063

The problem is this code.

x = int((-20 * p + (400 * p ** 2 + 4 * c) ** .5) / 2)

This code performs Floating-point subtraction. It causes loss of significance, because of two close numbers are subtracted.

>>> p = 10**15
>>> c = 10**15
>>> x = (-20 * p + (400 * p ** 2 + 4 * c) ** .5) / 2
>>> x
0.0

so, you should use integer sqrt instead of **.5. and change loop like this.

for i in range(len(pairs[0])):
    c = 100 * r + pairs[0][i]
    #x = int((-20 * p + (400 * p ** 2 + 4 * c) ** .5) / 2)
    x = (-20 * p + isqrt(400 * p ** 2 + 4 * c)) // 2
    y = 20 * p * x + x ** 2
    r = c - y
    p = 10 * p + x
    square_root[0].append(x)

for i in range(len(pairs[1])):
    #print(p,r,c)
    c = 100 * r + pairs[1][i]
    #x = int((-20 * p + (400 * p ** 2 + 4 * c) ** .5) / 2)
    x = (-20 * p + isqrt(400 * p ** 2 + 4 * c)) // 2
    y = 20 * p * x + x ** 2
    r = c - y
    p = 10 * p + x
    square_root[1].append(x)

and define isqrt -- the integer sqrt function

# from http://stackoverflow.com/questions/15390807/integer-square-root-in-python
def isqrt(n):
    x = n
    y = (x + 1) // 2
    while y < x:
        x = y
        y = (x + n // x) // 2
    return x

Then, you could get an accurated value of sqrt(2).

>>> print( sqrt([[2], [0], '+'], 25) )
[[1], [4, 1, 4, 2, 1, 3, 5, 6, 2, 3, 7, 3, 0, 9, 5, 0, 4, 8, 8, 0, 1, 6, 8, 8, 7], '+']

Thank you very much, that solved the problem. I'll know about this in the future. So, a**.5 =sqrt(a) is quite inaccurate a calculation for large numbers? — AnthonyD973, Apr 19 '15 at 11:34
Yes it is. (And subtraction) If you need more information about this, you should read: http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html — 0xrgb, Apr 19 '15 at 11:58

Limit of digit-by-digit calculation of square roots

1 Answers1