From d156541f31fa4095ac6a9ace2f4568285f1e0a6e Mon Sep 17 00:00:00 2001 From: Sidney Zhang Date: Wed, 16 Sep 2026 11:30:19 +0800 Subject: [PATCH] =?UTF-8?q?feat(solutions/0057):=20=E6=B7=BB=E5=8A=A0?= =?UTF-8?q?=E6=AC=A7=E6=8B=89=E8=AE=A1=E5=88=92=E7=AC=AC57=E9=A2=98?= =?UTF-8?q?=E8=BF=9E=E5=88=86=E6=95=B0=E5=B9=B3=E6=96=B9=E6=A0=B9=E8=A7=A3?= =?UTF-8?q?=E6=B3=95?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 包含正算递推实现、约分优化及详细数学推导文档 --- solutions/0057.SquareRoot/euler_57.py | 52 +++- solutions/0057.SquareRoot/readme.md | 337 ++++++++++++++++++++++++++ 2 files changed, 388 insertions(+), 1 deletion(-) create mode 100644 solutions/0057.SquareRoot/readme.md diff --git a/solutions/0057.SquareRoot/euler_57.py b/solutions/0057.SquareRoot/euler_57.py index 52ff584..2c818e9 100644 --- a/solutions/0057.SquareRoot/euler_57.py +++ b/solutions/0057.SquareRoot/euler_57.py @@ -1,10 +1,26 @@ ''' +It is possible to show that the square root of two can be expressed as an infinite continued fraction. +sqrt{2}=1+frac{1}{2+frac{1}{2+frac{1}{2+...}}} +By expanding this for the first four iterations, we get: +begin{aligned} + & 1+frac{1}{2}=frac{3}{2}=1.5 \\ + & 1+frac{1}{2+frac{1}{2}}=frac{7}{5}=1.4 \\ + & 1+frac{1}{2+frac{1}{2+frac{1}{2}}}=frac{17}{12}=1.41666... \\ + & 1+frac{1}{2+frac{1}{2+frac{1}{2+frac{1}{2}}}}=frac{41}{29}=1.41379 ldots +end{aligned} + +The next three expansions are 99/70, 239/169, and 577/408, but the eighth expansion, 1393/985, +is the first example where the number of digits in the numerator exceeds the number of digits in the denominator. + +In the first one-thousand expansions, how many fractions contain a numerator with more digits than the denominator? ''' import time from functools import wraps from typing import Any, Callable, TypeVar +from math import gcd + F = TypeVar("F", bound=Callable[..., Any]) @@ -37,4 +53,38 @@ def benchmark(repeat: int = 1) -> Callable[[F], F]: wrapper.total_time = 0.0 # type: ignore[attr-defined] return wrapper # type: ignore[return-value] - return decorator \ No newline at end of file + return decorator + + +def compute_square_root(n: int = 1000) -> list[tuple[int,int]] : + """ + Compute the square root of 2 up to n expansions. + """ + cf = [1] + [2 for _ in range(0, n)] + p2, p1 = 0, 1 # p_{-2}, p_{-1} + q2, q1 = 1, 0 # q_{-2}, q_{-1} + results = [] + for a in cf: + p = a * p1 + p2 + q = a * q1 + q2 + p2, p1 = p1, p + q2, q1 = q1, q + if gc := gcd(p, q) != 1: + p //= gc + q //= gc + results.append((p, q)) + return results + +@benchmark(20) +def main() -> None: + results = compute_square_root(1000) + count = 0 + for p, q in results: + if len(str(p)) > len(str(q)): + count += 1 + print(f"{count} fractions contain a numerator with more digits than the denominator.") + print(f"sqrt{2} = {results[-1][0]/results[-1][1]:.15f}") + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/solutions/0057.SquareRoot/readme.md b/solutions/0057.SquareRoot/readme.md new file mode 100644 index 0000000..4493d90 --- /dev/null +++ b/solutions/0057.SquareRoot/readme.md @@ -0,0 +1,337 @@ +方法1:倒算 + +- N = a_n ; D = 1 +- k:(n-1 → 0) => N,D = a_k*N+D ,N +- 最后: gcd(N,D) 约分 + +方法2:正算 + +正算如何成立呢?这个我还没有很好的想法。 + +``` +p2, p1 = 0, 1 # p_{-2}, p_{-1} +q2, q1 = 1, 0 # q_{-2}, q_{-1} + +for a in cf: + p = a * p1 + p2 + q = a * q1 + q2 + p2, p1 = p1, p + q2, q1 = q1, q + +# 结果就是 p1 / q1 +``` + + +----- + +广义连分数 + +形式为: +$b_0+\frac{a_1}{b_1+\frac{a_2}{b_2+\cdots+\frac{a_n}{b_n}}}$ + +用递推: +$A_{-1}=1,\quad A_0=b_0$ +$B_{-1}=0,\quad B_0=1$ +对 $k=1,\dots,n$ : +$A_k=b_kA_{k-1}+a_kA_{k-2}$ +$B_k=b_kB_{k-1}+a_kB_{k-2}$ + +最后: +$\text{值}=\frac{A_n}{B_n}$ + +简单连分数只是 $b_0=a_0,\ b_k=a_k,\ a_k=1$ 的特殊情况。 + + +----- + +## 1. 从最直观的“从后往前”看 + +简单连分数: + +$$ +[a_0;a_1,\dots,a_n] +=a_0+\frac1{a_1+\frac1{a_2+\cdots+\frac1{a_n}}} +$$ + +如果从最里面往外算,设当前尾部的值为: + +$$ +T_k=[a_k;a_{k+1},\dots,a_n]=\frac{N_k}{D_k} +$$ + +那么往外加一层 $a_{k-1}$,就变成: + +$$ +T_{k-1}=a_{k-1}+\frac1{T_k} +=a_{k-1}+\frac1{N_k/D_k} +=a_{k-1}+\frac{D_k}{N_k} +$$ + +通分: + +$$ +T_{k-1}=\frac{a_{k-1}N_k+D_k}{N_k} +$$ + +所以: + +$$ +N_{k-1}=a_{k-1}N_k+D_k,\qquad D_{k-1}=N_k +$$ + +这就是从后往前的递推: + +$$ +(N,D)\leftarrow (aN+D,\;N) +$$ + +它的道理非常直接: +**连分数每往外加一层,就是做一次 $x\mapsto a+1/x$。** +如果 $x=N/D$,那么: + +$$ +a+\frac1x=a+\frac DN=\frac{aN+D}{N} +$$ + +所以分子变成 $aN+D$,分母变成原来的分子 $N$。 + +--- + +## 2. 从前往后的 $p,q$ 递推为什么也对? + +从前往后定义收敛项: + +$$ +\frac{p_k}{q_k}=[a_0;a_1,\dots,a_k] +$$ + +递推公式是: + +$$ +p_k=a_kp_{k-1}+p_{k-2} +$$ + +$$ +q_k=a_kq_{k-1}+q_{k-2} +$$ + +初始条件: + +$$ +p_{-2}=0,\quad p_{-1}=1 +$$ + +$$ +q_{-2}=1,\quad q_{-1}=0 +$$ + +这个递推看起来像是突然出现的,但它其实来自矩阵乘法。 + +--- + +## 3. 关键:线性分式变换与矩阵表示 + +考虑变换: + +$$ +f_a(x)=a+\frac1x +$$ + +它可以写成: + +$$ +f_a(x)=\frac{ax+1}{x} +$$ + +这种形式叫线性分式变换。它可以对应一个 $2\times 2$ 矩阵: + +$$ +f_a \longleftrightarrow M(a)= +\begin{pmatrix} +a & 1\\ +1 & 0 +\end{pmatrix} +$$ + +为什么?因为线性分式变换: + +$$ +x\mapsto \frac{Ax+B}{Cx+D} +$$ + +对应矩阵: + +$$ +\begin{pmatrix} +A & B\\ +C & D +\end{pmatrix} +$$ + +而 $a+1/x=(ax+1)/(1x+0)$,所以矩阵就是: + +$$ +\begin{pmatrix} +a & 1\\ +1 & 0 +\end{pmatrix} +$$ + +线性分式变换有个重要性质: +**复合两个线性分式变换,对应矩阵相乘。** + +也就是说: + +$$ +f_{a_0}(f_{a_1}(x)) +$$ + +对应矩阵: + +$$ +M(a_0)M(a_1) +$$ + +于是整个连分数: + +$$ +[a_0;a_1,\dots,a_n] +$$ + +对应矩阵乘积: + +$$ +M(a_0)M(a_1)\cdots M(a_n) +$$ + +--- + +## 4. 矩阵乘积展开就是递推 + +设: + +$$ +M(a_0)M(a_1)\cdots M(a_k) += +\begin{pmatrix} +p_k & p_{k-1}\\ +q_k & q_{k-1} +\end{pmatrix} +$$ + +当 $k=-1$ 时,空乘积是单位矩阵: + +$$ +\begin{pmatrix} +1 & 0\\ +0 & 1 +\end{pmatrix} +$$ + +可以对应: + +$$ +p_{-1}=1,\quad p_{-2}=0,\quad q_{-1}=0,\quad q_{-2}=1 +$$ + +现在右乘一个 $M(a_k)$: + +$$ +\begin{pmatrix} +p_{k-1} & p_{k-2}\\ +q_{k-1} & q_{k-2} +\end{pmatrix} +\begin{pmatrix} +a_k & 1\\ +1 & 0 +\end{pmatrix} += +\begin{pmatrix} +a_kp_{k-1}+p_{k-2} & p_{k-1}\\ +a_kq_{k-1}+q_{k-2} & q_{k-1} +\end{pmatrix} +$$ + +所以新的左上角就是: + +$$ +p_k=a_kp_{k-1}+p_{k-2} +$$ + +新的左下角就是: + +$$ +q_k=a_kq_{k-1}+q_{k-2} +$$ + +这正是从前往后的递推公式。 + +最后,矩阵乘积作用在无穷大上: + +$$ +\begin{pmatrix} +p_n & p_{n-1}\\ +q_n & q_{n-1} +\end{pmatrix} +(\infty) += +\frac{p_n}{q_n} +$$ + +所以: + +$$ +[a_0;a_1,\dots,a_n]=\frac{p_n}{q_n} +$$ + +--- + +## 5. 总结一句话 + +递推之所以成立,是因为: + +$$ +a+\frac1x=\frac{ax+1}{x} +$$ + +每一步都是同一个线性分式变换。 +线性分式变换复合对应矩阵相乘,而矩阵相乘展开后,自然得到: + +$$ +p_k=a_kp_{k-1}+p_{k-2} +$$ + +$$ +q_k=a_kq_{k-1}+q_{k-2} +$$ + +从后往前算,是直接反复做 $N,D\leftarrow aN+D,N$; +从前往后算,是用 $p,q$ 递推维护当前收敛项的分子分母。 +两者本质上是同一件事,只是计算方向不同。 + +广义连分数也一样,如果每层是: + +$$ +b_k+\frac{a_k}{x} +$$ + +它对应矩阵: + +$$ +\begin{pmatrix} +b_k & a_k\\ +1 & 0 +\end{pmatrix} +$$ + +递推就变成: + +$$ +A_k=b_kA_{k-1}+a_kA_{k-2} +$$ + +$$ +B_k=b_kB_{k-1}+a_kB_{k-2} +$$ + +简单连分数只是 $b_k=a_k,\ a_k=1$ 的特殊情况。 \ No newline at end of file