Генерация цифр квадратного корня из 2
Я хочу генерировать цифры квадратного корня из двух- трех миллионов цифр.
Я знаю о Ньютоне-Рафсоне, но я не очень понимаю, как реализовать его в C или C++ из-за отсутствия поддержки biginteger. Может ли кто-нибудь указать мне правильное направление?
Кроме того, если кто-нибудь знает, как сделать это в Python (я новичок), я был бы также признателен.
9 ответов
Вы можете попробовать использовать отображение:
a/b -> (a+2b)/(a+b) начиная с a= 1, b= 1, Это сходится к sqrt(2) (фактически дает его непрерывные дробные представления).
Теперь ключевой момент: это можно представить как матричное умножение (аналогично Фибоначчи)
Если a_n и b_n - это n-е числа в шагах, то
[1 2] [a_n b_n]T = [a_ (n + 1) b_ (n + 1)]T
[1 1]
что сейчас дает нам
[1 2]n [a_1 b_1]T = [a_ (n + 1) b_ (n + 1)]T
[1 1]
Таким образом, если матрица 2x2 равна A, нам нужно вычислить An, что можно сделать путем многократного возведения в квадрат и использовать только целочисленную арифметику (поэтому вам не нужно беспокоиться о проблемах точности).
Также обратите внимание, что a/b, который вы получите, всегда будет в сокращенной форме (как gcd(a,b) = gcd(a+2b, a+b)), поэтому, если вы думаете об использовании класса дроби для представления промежуточного результатов нет!
Поскольку n-й знаменатель подобен (1+sqrt(2))^n, чтобы получить 3 миллиона цифр, вам, вероятно, потребуется вычислить до 3671656-го члена.
Обратите внимание, что даже если вы ищете ~3,6-миллионный термин, повторное возведение в квадрат позволит вам вычислить n-й член в умножениях и сложениях O (Log n).
Кроме того, это можно легко сделать параллельным, в отличие от итеративных, таких как Ньютон-Рафсон и т. Д.
РЕДАКТИРОВАТЬ: мне нравится эта версия лучше, чем предыдущая. Это общее решение, которое принимает как целые, так и десятичные дроби; при n = 2 и точности = 100000 это занимает около двух минут. Спасибо Полу Макгуайру за его предложения и другие предложения, добро пожаловать!
def sqrt_list(n, precision):
ndigits = [] # break n into list of digits
n_int = int(n)
n_fraction = n - n_int
while n_int: # generate list of digits of integral part
ndigits.append(n_int % 10)
n_int /= 10
if len(ndigits) % 2: ndigits.append(0) # ndigits will be processed in groups of 2
decimal_point_index = len(ndigits) / 2 # remember decimal point position
while n_fraction: # insert digits from fractional part
n_fraction *= 10
ndigits.insert(0, int(n_fraction))
n_fraction -= int(n_fraction)
if len(ndigits) % 2: ndigits.insert(0, 0) # ndigits will be processed in groups of 2
rootlist = []
root = carry = 0 # the algorithm
while root == 0 or (len(rootlist) < precision and (ndigits or carry != 0)):
carry = carry * 100
if ndigits: carry += ndigits.pop() * 10 + ndigits.pop()
x = 9
while (20 * root + x) * x > carry:
x -= 1
carry -= (20 * root + x) * x
root = root * 10 + x
rootlist.append(x)
return rootlist, decimal_point_index
Что касается произвольных больших чисел, вы можете взглянуть на Многофункциональную арифметическую библиотеку GNU (для C/C++).
Вот краткая версия для вычисления квадратного корня от целого числа до цифр точности. Он работает путем нахождения целочисленного квадратного корня из a после умножения на 10, возведенного в 2 x цифры.
def sqroot(a, digits):
a = a * (10**(2*digits))
x_prev = 0
x_next = 1 * (10**digits)
while x_prev != x_next:
x_prev = x_next
x_next = (x_prev + (a // x_prev)) >> 1
return x_next
Всего несколько предостережений.
Вам нужно будет преобразовать результат в строку и добавить десятичную точку в правильном месте (если вы хотите, чтобы десятичная точка была напечатана).
Преобразование очень большого целого числа в строку не очень быстро.
Деление очень больших целых чисел тоже не очень быстро (в Python).
В зависимости от производительности вашей системы вычисление квадратного корня из 2-3 миллионов знаков после запятой может занять час или более.
Я не доказал, что цикл всегда завершится. Он может колебаться между двумя значениями, отличающимися последней цифрой. Или не может.
Для работы? Используйте библиотеку!
Ради забавы? Повезло тебе:)
Напишите программу, которая будет имитировать то, что вы будете делать с карандашом и бумагой. Начните с 1 цифры, затем 2 цифры, затем 3, ..., ...
Не беспокойся о Ньютоне или о ком-либо еще. Просто сделай это по-своему.
Самый хороший способ, вероятно, использует продолжение расширения дроби [1; 2, 2, ...] квадратный корень из двух.
def root_two_cf_expansion():
yield 1
while True:
yield 2
def z(a,b,c,d, contfrac):
for x in contfrac:
while a > 0 and b > 0 and c > 0 and d > 0:
t = a // c
t2 = b // d
if not t == t2:
break
yield t
a = (10 * (a - c*t))
b = (10 * (b - d*t))
# continue with same fraction, don't pull new x
a, b = x*a+b, a
c, d = x*c+d, c
for digit in rdigits(a, c):
yield digit
def rdigits(p, q):
while p > 0:
if p > q:
d = p // q
p = p - q * d
else:
d = (10 * p) // q
p = 10 * p - q * d
yield d
def decimal(contfrac):
return z(1,0,0,1,contfrac)
decimal((root_two_cf_expansion()) возвращает итератор всех десятичных цифр. t1 а также t2 В алгоритме указаны минимальные и максимальные значения следующей цифры. Когда они равны, мы выводим эту цифру.
Обратите внимание, что это не обрабатывает некоторые исключительные случаи, такие как отрицательные числа в непрерывной дроби.
(Этот код является адаптацией кода на Haskell для обработки непрерывных дробных фрагментов.)
Ну, вот код, который я написал. Он сгенерировал миллион цифр после десятичного числа для квадратного корня из 2 примерно за 60800 секунд для меня, но мой ноутбук спал, когда он запускал программу, это должно быть быстрее. Вы можете попытаться сгенерировать 3 миллиона цифр, но это может занять пару дней.
def sqrt(number,digits_after_decimal=20):
import time
start=time.time()
original_number=number
number=str(number)
list=[]
for a in range(len(number)):
if number[a]=='.':
decimal_point_locaiton=a
break
if a==len(number)-1:
number+='.'
decimal_point_locaiton=a+1
if decimal_point_locaiton/2!=round(decimal_point_locaiton/2):
number='0'+number
decimal_point_locaiton+=1
if len(number)/2!=round(len(number)/2):
number+='0'
number=number[:decimal_point_locaiton]+number[decimal_point_locaiton+1:]
decimal_point_ans=int((decimal_point_locaiton-2)/2)+1
for a in range(0,len(number),2):
if number[a]!='0':
list.append(eval(number[a:a+2]))
else:
try:
list.append(eval(number[a+1]))
except IndexError:
pass
p=0
c=list[0]
x=0
ans=''
for a in range(len(list)):
while c>=(20*p+x)*(x):
x+=1
y=(20*p+x-1)*(x-1)
p=p*10+x-1
ans+=str(x-1)
c-=y
try:
c=c*100+list[a+1]
except IndexError:
c=c*100
while c!=0:
x=0
while c>=(20*p+x)*(x):
x+=1
y=(20*p+x-1)*(x-1)
p=p*10+x-1
ans+=str(x-1)
c-=y
c=c*100
if len(ans)-decimal_point_ans>=digits_after_decimal:
break
ans=ans[:decimal_point_ans]+'.'+ans[decimal_point_ans:]
total=time.time()-start
return ans,total
Вот более эффективная целочисленная функция квадратного корня (в Python 3.x), которая должна завершаться во всех случаях. Он начинается с числа, намного ближе к квадратному корню, поэтому требуется меньше шагов. Обратите внимание, что для int.bit_length требуется Python 3.1+. Ошибка проверки опущена для краткости.
def isqrt(n):
x = (n >> n.bit_length() // 2) + 1
result = (x + n // x) // 2
while abs(result - x) > 1:
x = result
result = (x + n // x) // 2
while result * result > n:
result -= 1
return result
Python уже поддерживает большие целые числа из коробки, и если это единственное, что сдерживает вас в C/C++, вы всегда можете написать быстрый контейнерный класс самостоятельно.
Единственная проблема, которую вы упомянули, это отсутствие больших целых чисел. Если вы не хотите использовать для этого библиотеку, ищите ли вы помощь в написании такого класса?