• 一、连分数简介
  • 1、实数的连分数
  • 2、二次无理数的连分数
• 二、渐近分数列
• 三、丢番图方程简介
  • 1、佩尔方程的解
  • 2、Diophantine Equation

一、连分数简介

在数学中，连分数 即如下表达式：

其中 $a_i\in \mathbb{Z},i\in \mathbb{N}$ 。

连分数具有如下性质：

• 有理数的连分数表示只有有限的项，并且当没有尾随的 1 时，其表示方法唯一；
• 无理数的连分数表示是唯一的；
• 勒让德定理：二次无理数 $y$ 的连分数表示会产生循环，即 $y=[a_0;\overline{a_1,a_2,\cdots,a_2,a_1,2a_0}],a_0=[\,y\,]$；
• 数 $x$ 的截断连分数表示是在特定意义上最佳可能的有理数逼近。

1、实数的连分数

表 1：求 3.245 的连分数

1、Lisp 的求法

(defun frac (n)
	(write (floor n))(format t " ")
	(setq a (- n (floor n))) 
	(if (/= a 0) (frac (/ 1 a)))
)(frac (/ 3245 1000))

3 4 12 4

2、数值灾难

Lisp 不会自动求出 3245/1000 的结果，因此计算过程是分数的计算而非浮点运算，可以得到精确结果。但如果用浮点数实现的话，可能导致 数值灾难。

由于浮点数被用以在计算机中对实数进行近似表示，因此浮点运算与数学运算有所差异。为了得到精确的结果，需要用 欧几里得 GCD 算法 的变体作为替代。

3、欧几里得 GCD 算法（辗转相除法）

def GCD(a,b,ans=[]): # 实数 a/b 的连分数
	if b:	return GCD(b, a%b, ans+[int(a/b)])
	else:	return ans
print(GCD(3245,1000))

[3, 4, 12, 4]

2、二次无理数的连分数

上述算法对无理数的连分数有些无能为力。下面着重讨论无理数 $\sqrt{n}$ 的连分数展开技巧¹：

表 2：求 $\sqrt{7}$ 的连分数($[\,\sqrt{7}\,]=2$)

将上述算法一般化：设 $k=[\,\sqrt{n}\,]$，并令 $n_c=n-c^2$。

其中，若求得 $\frac{i_{l-1}+i_l}{q}$，则取满足 $i_{l+1}\leq k$ 的最大值，使 $i_l+i_{l+1}$ 能被 $n_l/q$ 整除。

如果考虑将 $\frac{i+j}{q}=a \Rightarrow i-\frac{a}{q}-j$ 的形式，有：

其中，$\frac{n_{i_1}}{n_k}=q_1,\,\frac{n_{i_2}}{q_1}=q_2,\cdots$

二、渐近分数列

如下表所列，渐近分数常用于无理数的逼近。

表 3：连分数与渐近分数列

数学上可以证明：由连分数得到的渐近分数，其值是在分子或分母小于下一个渐近分数的分数中，最接近精确值的近似值。

对渐进分数列 $\{\frac{p_i}{q_i},\,i\in \mathbb{N}^+\}$，有 $\frac{p_1}{q_1}=\frac{a_0}{1},\frac{p_1}{q_1}=\frac{a_1a_0+1}{a_1}$，其递推式为：

$$\begin{bmatrix}p_{i+1}\\q_{i+1}\end{bmatrix}=\begin{bmatrix}p_i&p_{i-1}\\q_i&q_{i-1}\end{bmatrix}\begin{bmatrix}a_{i}\\1\end{bmatrix}$$

显然，渐近分数列的奇数项小于原值，偶数项大于原值。

三、丢番图方程简介

丢番图方程 又名整系数多项式方程，其变量与系数均为整数。若一个丢番图方程具有以下的形式：

则称此二元二次不定方程为佩尔方程。

1、佩尔方程的解

• 若 $n$ 是完全平方数，佩尔方程只有平凡解 $(\pm 1, 0)$。
• 对于其余情况，拉格朗日证明了佩尔方程总有非平凡解，而这些解可由 $\sqrt{n}$ 的连分数求出。

设 $\{\frac{p_i}{q_i}\}$ 是 $\sqrt{n}$ 的的渐近分数列，则存在 $i$ 使得 $(p_i,q_i)$ 为佩尔方程的解。取 $\min i$ 对应的 $(p_i,q_i)$ 作为佩尔方程的基本解（最小解），记作$(x_1,y_1)$，则全部的非平凡解解 $\{x,y\}$ 可表示为 $x_i + y_i\sqrt n = (x_1 + y_1\sqrt n)^i$，或由递推关系式求得：

$$\begin{bmatrix}x_{i+1}\\y_{i+1}\end{bmatrix}=\begin{bmatrix}x_1&ny_1\\y_1&x_1\end{bmatrix}\begin{bmatrix}x_{i}\\y_{i}\end{bmatrix}$$

规律：记 $\sqrt{n}$ 的连分数循环节长度为 $l$。

1. 当 $l$ 为偶数，$(p_l,q_l)$ 即是基本解；
2. 当 $l$ 为奇数，$(p_{2l},q_{2l})$ 即是基本解。

Example：求方程 $x^2 - 7y^2= 1$ 的解

由 $\sqrt{7}$ 的渐近分数列 $\{\frac{p_i}{q_i}\}=\frac{2}{1},\frac{3}{1},\frac{5}{2},\frac{8}{3},\cdots$，求得 $x^2 - 7y^2= 1$ 的基本解 $(8,3)$，其所有非平凡解

2、Diophantine Equation²

题意：求最小的非完全平方数 $n\in(1,1000]$ 使佩尔方程 $x^2 - ny^2= 1$ 的基本解 $x_1$ 最大。

from math import *
def iqj(i0,q0,j0,d):
	i = j0
	q = (d - j0*j0)//q0
	max_num = j0 + floor(sqrt(d))
	j = max([x for x in range(max_num,0,-1) if x % q == 0]) - j0
	a = (i + j)//q
	return i,q,j,a
def pell (d,h,k):
	if h*h - d*k*k == 1:	return True
	else:		        return False 
def frac(d):
	i0 = 0
	q0 = 1
	j0 = floor(sqrt(d))
	a0 = (i0 + j0)//q0
	h0 = a0
	k0 = 1
	if pell(d,h0,k0):	return h0
	else:
		i0,q0,j0,a1=iqj(i0,q0,j0,d)
		h1 = a1*a0 + 1
		k1 = a1
		while(not pell(d,h1,k1)):
			i0,q0,j0,a=iqj(i0,q0,j0,d)
			h = a*h1 + h0
			k = a*k1 + k0
			h0,k0 = h1,k1
			h1,k1 = h,k
		return h1
mh = 0
ans = 0
for d in [i+1 for i in range(1000) if sqrt(i+1)!=int(sqrt(i+1))]:
	h = frac(d)
	if h > mh:
		mh = h
		ans = d
print(ans)

661

如果您喜欢这篇文章，欢迎转载、演绎或用于商业目的，但请务必保留作者署名以及本文链接！

Copyright © Pandaman