math - 连分数和 Pell 方程 - 数值问题
问题描述
数学背景
连分数是一种表示数字(有理或无理)的方法,使用基本的递归公式来计算它。给定一个数字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]
并且对于求解sqrt( )的连分数展开中的方程y=k[n]
的最低找到。n
D
问题
我无法让这个简单的算法为D=61
. 我首先注意到它没有解决 Pell 方程的 100 个系数,所以我将它与 Wolfram Alpha 的收敛和连分数表示进行比较,并注意到第 20 个元素失败 - 表示3
与4
我得到的比较,产生不同的收敛 -h[20]=335159612
在 Wolfram 上425680601
与我。
我在两个系统上测试了下面的代码,两种语言(虽然公平地说,Python 是 C 引擎盖下我猜),并得到相同的结果 - 循环 20 上的差异。我会注意到收敛仍然是准确的并且收敛!与 Wolfram Alpha 相比,为什么我得到不同的结果,是否可以修复它?
为了测试,这里有一个 Python 程序来求解 Pell 方程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;
}
解决方案
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]
推荐阅读
- javascript - 对对象数组进行排序,但第一位是固定的
- matlab - 在 3D 矩阵中查找具有特定值的簇
- c# - 如何在选择游戏对象时将材质分配给游戏对象,并在使用 C# 取消选择游戏对象时将其更改回原始材质?
- facebook - 如何让 Facebook 的访问令牌显示我的页面提要?
- google-cloud-platform - 指定 Cloud Vision 请求的输出文件名
- c# - 刚体移动时,刚体速度有时为零
- python - Flask Restplus API 的 Swagger UI 中出现错误 401
- python - 如何使用文件管道使用 Python/Scrapy 下载 (PDF) 文件?
- css - CSS:你能在一个类中引用一个类吗?
- google-apps-script - 是否可以为 Google 表单生成自动填充脚本?