feat(solutions/0057): 添加欧拉计划第57题连分数平方根解法

包含正算递推实现、约分优化及详细数学推导文档

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
+    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()

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$ 的特殊情况。