Предел вычисления цифр квадратными корнями по цифрам

Question

Я пытаюсь получить как хорошую оценку pi, так как я могу использовать алгоритм Чудновского в Python. Этот алгоритм подразумевает получение квадратного корня из 640 320.

После некоторых исследований я нашел довольно эффективный способ вычисления квадратных корней; метод называется «вычислением по цифрам» (см. here). Поэтому, пытаясь реализовать его, я обнаружил, что первые 13 десятичных знаков верны, и затем я получаю странные результаты (следующий - это 0 вместо 4, а затем следующая цифра - 128, затем -1024 ...)

Я попытался проверить свою функцию, но это выглядит хорошо для меня (кроме того, я бы, вероятно, не нашел правильные первые 13 десятичных знаков в противном случае). Таким образом, мой вопрос: существуют ли какие-то ограничения в этом методе расчета по цифре?

Если, случайно, вы хотели бы видеть мой код, вот он:

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

Answer 1

Проблема заключается в код.

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

Этот код выполняет с плавающей запятой вычитание. Это приводит к потере значимости, поскольку два близких числа вычитаются.

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

так, вы должны использовать целочисленный SQRT вместо **.5. и измените петлю следующим образом.

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) 

и определить isqrt - функция целое SQRT

# 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 

Тогда вы могли бы получить accurated значение 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], '+']