连续分数和佩尔方程-数值问题

Question

数学背景

Continued fractions是一种表示数字（有理或无理）的方法，用basic recursion formula进行计算。给定数字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是最终表示。我们还可以通过

定义一系列的convergents

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。

Pell's equation是形式为x^2-D*y^2=1的问题，其中所有数字都是整数，在我们的情况下D不是理想的平方。给定D的解决方案，可将x is given by continued fractions最小化。基本上，对于上述方程式，可以保证此（基本）解为x=h[n]，对于找到的最低y=k[n]为n，它可以解决sqrt（D）]的连续分数展开式>。

---

问题

我无法使D=61的这种简单算法起作用。我首先注意到它不能解决100个系数的Pell方程，因此我将它与Wolfram Alpha的convergents和continued fraction representation进行了比较，并注意到第20个元素失败-与3相比，表示是4，得出不同的会聚点-Wolfram上的h[20]=335159612与我的425680601相比。

我在以下两个系统上测试了下面的代码，两种语言（虽然公平地说，我猜是Python是C），并且得到相同的结果-循环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
一个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;
}
数学背景连续分数是一种表示数字（有理还是非有理数）的方法，使用基本的递归公式进行计算。给定数字r，我们定义r [0] = r并具有：对于...中的n，]]

mpmath，可以使用python的多精度库。请注意所有重要数字均为mp格式。

在下面的代码中，x, y and i是标准的多精度整数。 r和f是多精度实数。请注意，初始计数设置为高于20。

from mpmath import mp, mpf

mp.dps = 50  # precision in number of decimal digits

def continued_fraction(D, count=22, thresh=mpf(1E-12), verbose=False):
    cf = []
    h = (0, 1)
    k = (1, 0)
    r = start = mp.sqrt(D)
    initial_count = count
    x = 0 # some dummy starting values, they will be overwritten early in the while loop
    y = 1
    while abs(x/y - start) > thresh and count > 0:
        i = int(mp.floor(r))
        cf.append(i)
        x, y = i*h[-1] + h[-2], i*k[-1] + k[-2]
        if verbose or 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
        f = r - i
        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=22, verbose=True, thresh=mpf(1e-100))

输出类似于Wolfram的：

...
335159612²-61x42912791² = 3
1431159437²-61x183241189² = -12
1766319049²-61x226153980² = 1
[7, 1, 4, 3, 1, 2, 2, 1, 3, 4, 1, 14, 1, 4, 3, 1, 2, 2, 1, 3, 4, 1]

Answer 1

mpmath，可以使用python的多精度库。请注意所有重要数字均为mp格式。