Математическое обоснование
Непрерывные дроби представляют собой способ представления чисел (рациональных или нет) с базисом c формулой рекурсии чтобы рассчитать это. Для числа r мы определяем r[0]=r и имеем:
for n in range(0..N):
a[n] = floor(r[n])
if r[n] == [an]: break
r[n+1] = 1 / (r[n]-a[n])
, где a - окончательное представление. Мы также можем определить ряд сходящихся по
h[-2,-1] = [0, 1]
k[-2, -1] = [1, 0]
h[n] = a[n]*h[n-1]+h[n-2]
k[n] = a[n]*k[n-1]+k[n-2]
, где h[n]/k[n] сходятся к r .
уравнение Пелла - это проблема вида x^2-D*y^2=1, где все числа являются целыми числами, а D не является идеальным квадратом в нашем случае. Решение для данного D, которое минимизирует x , дается непрерывными дробями . В принципе, для приведенного выше уравнения гарантируется, что это (фундаментальное) решение равно x=h[n] и y=k[n] для самого низкого найденного n, что решает уравнение в продолженном расширении дроби sqrt (D) .
Проблема
Мне не удается заставить этот простой алгоритм работать для D=61. Сначала я заметил, что оно не решает уравнение Пелла для 100 коэффициентов, поэтому я сравнил его с сходимостями Вольфрама Альфы и продолжение дробного представления и заметил, что 20-й элемент не работает - представление 3 по сравнению с 4, который я получаю, получая различные конвергенты - h[20]=335159612 на Вольфраме по сравнению с 425680601 для меня.
Я тестировал приведенный ниже код, два языка (хотя, если честно, Python - это C под капотом, я думаю), на двух системах и получаю одинаковый результат - разность на l oop 20. Замечу, что конвергенты все еще точны и сходятся! Почему я получаю другие результаты по сравнению с Wolfram Alpha, и можно ли это исправить?
Для тестирования вот программа Python для решения уравнения Пелла для D=61, печать первых 20 конвергентов и представление непрерывной дроби cf (и некоторый дополнительный ненужный пух):
from math import floor, sqrt # Can use mpmath here as well.
def continued_fraction(D, count=100, thresh=1E-12, verbose=False):
cf = []
h = (0, 1)
k = (1, 0)
r = start = sqrt(D)
initial_count = count
x = (1+thresh+start)*start
y = start
while abs(x/y - start) > thresh and count:
i = int(floor(r))
cf.append(i)
f = r - i
x, y = i*h[-1] + h[-2], i*k[-1] + k[-2]
if verbose is True or verbose == initial_count-count:
print(f'{x}\u00B2-{D}x{y}\u00B2 = {x**2-D*y**2}')
if x**2 - D*y**2 == 1:
print(f'{x}\u00B2-{D}x{y}\u00B2 = {x**2-D*y**2}')
print(cf)
return
count -= 1
r = 1/f
h = (h[1], x)
k = (k[1], y)
print(cf)
raise OverflowError(f"Converged on {x} {y} with count {count} and diff {abs(start-x/y)}!")
continued_fraction(61, count=20, verbose=True, thresh=-1) # We don't want to stop on account of thresh in this example
A c программа делает то же самое:
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
int main() {
long D = 61;
double start = sqrt(D);
long h[] = {0, 1};
long k[] = {1, 0};
int count = 20;
float thresh = 1E-12;
double r = start;
long x = (1+thresh+start)*start;
long y = start;
while(abs(x/(double)y-start) > -1 && count) {
long i = floor(r);
double f = r - i;
x = i * h[1] + h[0];
y = i * k[1] + k[0];
printf("%ld\u00B2-%ldx%ld\u00B2 = %lf\n", x, D, y, x*x-D*y*y);
r = 1/f;
--count;
h[0] = h[1];
h[1] = x;
k[0] = k[1];
k[1] = y;
}
return 0;
}