From d01075218445d051984c36aae7df65e633b7b04e Mon Sep 17 00:00:00 2001 From: Sidney Zhang Date: Mon, 27 Apr 2026 17:09:27 +0800 Subject: [PATCH] =?UTF-8?q?fix:=20=E7=BB=9F=E4=B8=80=E6=96=87=E4=BB=B6?= =?UTF-8?q?=E7=BC=96=E7=A0=81=E4=B8=BA=20utf-8=20=E5=B9=B6=E5=B0=86=20benc?= =?UTF-8?q?hmark=20=E6=97=A5=E5=BF=97=E6=94=B9=E4=B8=BA=E8=8B=B1=E6=96=87?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- main.py | 18 +-- solutions/0066.DiophantineEq/README.md | 169 +++++++++++++++++++++++ solutions/0066.DiophantineEq/euler_66.py | 92 ++++++++++++ 3 files changed, 266 insertions(+), 13 deletions(-) create mode 100644 solutions/0066.DiophantineEq/README.md create mode 100644 solutions/0066.DiophantineEq/euler_66.py diff --git a/main.py b/main.py index 0ac29e0..d0de7c5 100644 --- a/main.py +++ b/main.py @@ -17,7 +17,7 @@ def add_newproblem(num: int, name: str | None = None) -> None: Path(f"solutions/{title}").mkdir(parents=True, exist_ok=True) file_name = f"solutions/{title}/euler_{num}.py" Path(file_name).touch() - with open(file_name, "w") as f: + with open(file_name, "w", encoding="utf-8") as f: f.write("""''' ''' @@ -29,14 +29,8 @@ F = TypeVar("F", bound=Callable[..., Any]) def benchmark(repeat: int = 1) -> Callable[[F], F]: - ''' - 重复运行目标函数并计算平均耗时。 - - 平均耗时会被存储在 wrapper.avg_time 和 wrapper.total_time 中， - 同时会打印到控制台。函数的返回值不受影响。 - ''' if repeat < 1: - raise ValueError("repeat 必须 >= 1") + raise ValueError("repeat must >= 1") def decorator(func: F) -> F: @wraps(func) @@ -54,18 +48,16 @@ def benchmark(repeat: int = 1) -> Callable[[F], F]: wrapper.total_time = total # type: ignore[attr-defined] print( - f"[Benchmark] {func.__name__} | 重复 {repeat} 次 | " - f"平均: {wrapper.avg_time:.6f}s | 总计: {wrapper.total_time:.6f}s" # type: ignore[attr-defined] + f"[Benchmark] {func.__name__} | repeated {repeat} times | " + f"average: {wrapper.avg_time:.6f}s | total: {wrapper.total_time:.6f}s" # type: ignore[attr-defined] ) return result - # 初始化属性，避免调用前访问报错 wrapper.avg_time = 0.0 # type: ignore[attr-defined] wrapper.total_time = 0.0 # type: ignore[attr-defined] return wrapper # type: ignore[return-value] - return decorator - """) + return decorator""") @app.command("list", help="List all problems in the projecteuler repository.") diff --git a/solutions/0066.DiophantineEq/README.md b/solutions/0066.DiophantineEq/README.md new file mode 100644 index 0000000..51564f9 --- /dev/null +++ b/solutions/0066.DiophantineEq/README.md @@ -0,0 +1,169 @@ +# 佩尔方程（Pell's Equation）—— 基本解的完整数学理论 + +## 一、方程定义与历史背景 + +**佩尔方程**是指一类特殊的二元二次不定方程： + +$$x^2 - Dy^2 = 1$$ + +其中 $D$ 是一个**非平方正整数**（若 $D$ 为完全平方数，则方程退化为 $(x - y\sqrt{D})(x + y\sqrt{D}) = 1$，仅有平凡整数解）。 + +### 历史脉络 +- **公元前400年左右**：古希腊数学家就已研究 $x^2 - 2y^2 = 1$，发现解 $(3,2), (17,12), (99,70), \ldots$ +- **Brahmagupta（628年）**：印度数学家系统研究了"Varga-prakriti"方程，提出"bhāvanā"合成法则，本质上发现了 $(x_1 + y_1\sqrt{D})^n$ 的通解结构 +- **Bhaskara II（1150年）**：在《Līlāvatī》中给出求解 $x^2 - 61y^2 = 1$ 的算法（循环法 cakravāla），得到基本解 $(1766319049, 226153980)$ +- **Fermat（1657年）**：向欧洲数学界挑战求解此方程，Euler 误将解法归功于 Pell，由此得名 + +--- + +## 二、基本理论：解的存在性与结构 + +### 定理 1（Lagrange, 1768）：解的存在性 +对任意非平方正整数 $D$，方程 $x^2 - Dy^2 = 1$ **总有无穷多组正整数解**。 + +### 定义：基本解（最小解） +所有正整数解中使 $x + y\sqrt{D}$ 最小的那组解 $(x_1, y_1)$ 称为**基本解**（或**最小解**、**本源解**）。它是生成全部解的"种子"。 + +### 定理 2：通解结构 +若 $(x_1, y_1)$ 是基本解，则**所有**正整数解 $(x_n, y_n)$ 满足： + +$$x_n + y_n\sqrt{D} = (x_1 + y_1\sqrt{D})^n, \quad n = 1, 2, 3, \ldots$$ + +等价地，解可通过递推得到： +$$\begin{cases} x_{n+1} = x_1 x_n + D y_1 y_n \\ y_{n+1} = y_1 x_n + x_1 y_n \end{cases}$$ + +或写成矩阵形式： +$$\begin{pmatrix} x_{n+1} \\ y_{n+1} \end{pmatrix} = \begin{pmatrix} x_1 & Dy_1 \\ y_1 & x_1 \end{pmatrix} \begin{pmatrix} x_n \\ y_n \end{pmatrix}$$ + +--- + +## 三、核心算法：连分数法 + +求解基本解的标准方法是**连分数展开** $\sqrt{D}$。 + +### 3.1 $\sqrt{D}$ 的连分数展开 + +对非平方正整数 $D$，$\sqrt{D}$ 有**周期性连分数展开**： + +$$\sqrt{D} = [a_0; \overline{a_1, a_2, \ldots, a_L}]$$ + +其中 $a_0 = \lfloor\sqrt{D}\rfloor$，$(a_1, a_2, \ldots, a_L)$ 是纯循环部分，$L$ 为周期长度。 + +### 3.2 计算递推关系 + +展开算法的核心递推（计算周期中各项）： +$$\begin{aligned} +m_0 &= 0, \quad d_0 = 1, \quad a_0 = \lfloor\sqrt{D}\rfloor \\ +m_{k+1} &= d_k \cdot a_k - m_k \\ +d_{k+1} &= \frac{D - m_{k+1}^2}{d_k} \\ +a_{k+1} &= \left\lfloor \frac{a_0 + m_{k+1}}{d_{k+1}} \right\rfloor +\end{aligned}$$ + +当 $(m_k, d_k, a_k)$ 首次重复时，即得到一个完整周期。 + +### 3.3 渐近分数（Convergents） + +记第 $n$ 个渐近分数为 $\frac{p_n}{q_n}$，递推公式： +$$\begin{cases} +p_{-1} = 1, & p_0 = a_0, \quad p_n = a_n p_{n-1} + p_{n-2} \\[4pt] +q_{-1} = 0, & q_0 = 1, \quad q_n = a_n q_{n-1} + q_{n-2} +\end{cases}$$ + +### 3.4 关键定理（Lagrange） + +> **周期 $L$ 为偶数**：基本解 $(x_1, y_1) = (p_{L-1}, q_{L-1})$，即第 $(L-1)$ 个渐近分数满足 $p_{L-1}^2 - D q_{L-1}^2 = 1$ +> +> **周期 $L$ 为奇数**：基本解 $(x_1, y_1) = (p_{2L-1}, q_{2L-1})$，需要走到第 $(2L-1)$ 个渐近分数 + +--- + +## 四、经典算例 + +### 例 1：$D = 2$ + +$$\sqrt{2} = [1; \overline{2}], \quad L = 1 \text{（奇数）}$$ + +渐近分数序列： +$$\frac{1}{1}, \frac{3}{2}, \frac{7}{5}, \frac{17}{12}, \ldots$$ + +检验：$3^2 - 2 \cdot 2^2 = 9 - 8 = 1$ ✓ + +**基本解：$(x_1, y_1) = (3, 2)$** + +通解：$(3 + 2\sqrt{2})^n = x_n + y_n\sqrt{2}$，前几项： +- $n=2$: $(17, 12)$，$17^2 - 2 \cdot 12^2 = 289 - 288 = 1$ +- $n=3$: $(99, 70)$，$99^2 - 2 \cdot 70^2 = 9801 - 9800 = 1$ +- $n=4$: $(577, 408)$ + +### 例 2：$D = 13$ + +$$\sqrt{13} = [3; \overline{1, 1, 1, 1, 6}], \quad L = 5 \text{（奇数）}$$ + +需要计算到第 $2L - 1 = 9$ 个渐近分数： + +$$\frac{p_9}{q_9} = \frac{649}{180}$$ + +验证：$649^2 - 13 \cdot 180^2 = 421201 - 421200 = 1$ ✓ + +**基本解：$(649, 180)$** + +### 例 3：$D = 61$（历史上最著名的难题） + +$$\sqrt{61} = [7; \overline{1, 4, 3, 1, 2, 2, 1, 3, 4, 1, 14}], \quad L = 11 \text{（奇数）}$$ + +周期长达 11，需走到第 $2 \times 11 - 1 = 21$ 个渐近分数： + +**基本解：$(1766319049, 226153980)$** + +验证：$1766319049^2 - 61 \times 226153980^2 = 1$ ✓ + +这正是 Bhaskara II 在 12 世纪手工计算出的惊人结果。 + +--- + +## 五、负佩尔方程：$x^2 - Dy^2 = -1$ + +### 定理 3：可解性判据 +方程 $x^2 - Dy^2 = -1$ 有整数解**当且仅当** $\sqrt{D}$ 的连分数周期 $L$ 为**奇数**。 + +### 定理 4：基本解位置 +当 $L$ 为奇数时，$x^2 - Dy^2 = -1$ 的基本解由第 $(L-1)$ 个渐近分数给出；而 $x^2 - Dy^2 = 1$ 的基本解则由第 $(2L-1)$ 个渐近分数给出。 + +### 例：$D = 5$ + +$$\sqrt{5} = [2; \overline{4}], \quad L = 1 \text{（奇数）}$$ + +- $x^2 - 5y^2 = -1$ 的基本解：$(p_0, q_0) = (2, 1)$，$2^2 - 5 \cdot 1^2 = 4 - 5 = -1$ ✓ +- $x^2 - 5y^2 = 1$ 的基本解：$(p_1, q_1) = (9, 4)$，$9^2 - 5 \cdot 4^2 = 81 - 80 = 1$ ✓ + +对比 $D = 3$（$L = 2$ 偶数）：$x^2 - 3y^2 = -1$ **无解**。 + +--- + +## 六、数值验证汇总 + +| $D$ | 周期 $L$ | 周期奇偶 | $x^2 - Dy^2 = 1$ 基本解 | $x^2 - Dy^2 = -1$ | +|:---:|:---:|:---:|:---|:---| +| 2 | 1 | 奇 | $(3, 2)$ | $(1, 1)$ ✓ | +| 3 | 2 | 偶 | $(2, 1)$ | 无解 | +| 5 | 1 | 奇 | $(9, 4)$ | $(2, 1)$ ✓ | +| 6 | 2 | 偶 | $(5, 2)$ | 无解 | +| 7 | 4 | 偶 | $(8, 3)$ | 无解 | +| 10 | 1 | 奇 | $(19, 6)$ | $(3, 1)$ ✓ | +| 13 | 5 | 奇 | $(649, 180)$ | $(18, 5)$ ✓ | +| 29 | 5 | 奇 | $(9801, 1820)$ | $(70, 13)$ ✓ | +| 61 | 11 | 奇 | $(1766319049, 226153980)$ | $(29718, 3805)$ ✓ | + +--- + +## 七、算法复杂度与注意事项 + +1. **基本解的大小**：$(x_1, y_1)$ 可以**极其巨大**。例如 $D = 991$ 时，$x_1 \approx 3.7 \times 10^{47}$，$y_1 \approx 1.2 \times 10^{46}$。基本解的位数关于 $D$ 可以是指数级增长。 + +2. **计算精度**：实现时必须使用任意精度整数运算（如 Python 的 `int` 类型），标准 64 位整数远不够。 + +3. **与代数数论的联系**：佩尔方程的解群 $\{(x, y) : x^2 - Dy^2 = 1\}$ 同构于实二次域 $\mathbb{Q}(\sqrt{D})$ 中整数环的单位群，基本解对应**基本单位元**（Fundamental Unit）。 + +4. **广义佩尔方程**：$x^2 - Dy^2 = N$（$N \neq \pm 1$）的求解可通过先求 $x^2 - Dy^2 = 1$ 的基本解，再结合特解得到所有解。 + +5. **递推时的细节**：需要注意的一点是，n为偶数时，基本解为 $p_{2n-1}, q_{2n-1}$ ，计算时需要重复a序列，需要重复的是 $a_1$ 开始的a序列，那么 $2n-1$ 项就只需要重复到 $2 a_0$ 的前一项。 diff --git a/solutions/0066.DiophantineEq/euler_66.py b/solutions/0066.DiophantineEq/euler_66.py new file mode 100644 index 0000000..6e1fee3 --- /dev/null +++ b/solutions/0066.DiophantineEq/euler_66.py @@ -0,0 +1,92 @@ +""" +Consider quadratic Diophantine equations of the form: + x^2 - D y^2 = 1 +For example, when D = 13 , the minimal solution in x is $649^2 - 13 * 180^2 = 1$ . +It can be assumed that there are no solutions in positive integers when D is square. +By finding minimal solutions in x for D={2,3,5,6,7} , we obtain the following: + 3^2 - 2 * 2^2 = 1 + 2^2 - 3 * 2^2 = 1 + 9^2 - 5 * 4^2 = 1 + 5^2 - 6 * 2^2 = 1 + 8^2 - 7 * 3^2 = 1 +Hence, by considering minimal solutions in x for D <= 7 , +the largest x is obtained when D = 5 . +Find the value of D <= 1000 in minimal solutions of x for +which the largest value of x is obtained. +""" + +import time +from functools import wraps +from math import isqrt +from typing import Any, Callable, TypeVar + +F = TypeVar("F", bound=Callable[..., Any]) + + +def benchmark(repeat: int = 1) -> Callable[[F], F]: + if repeat < 1: + raise ValueError("repeat >= 1") + + def decorator(func: F) -> F: + @wraps(func) + def wrapper(*args: Any, **kwargs: Any) -> Any: + total = 0.0 + result = None + + for _ in range(repeat): + start = time.perf_counter() + result = func(*args, **kwargs) + end = time.perf_counter() + total += end - start + + wrapper.avg_time = total / repeat # type: ignore[attr-defined] + wrapper.total_time = total # type: ignore[attr-defined] + + print( + f"[Benchmark] {func.__name__} | repeated {repeat} times | " + f"average: {wrapper.avg_time:.6f}s | total: {wrapper.total_time:.6f}s" # type: ignore[attr-defined] + ) + return result + + wrapper.avg_time = 0.0 # type: ignore[attr-defined] + wrapper.total_time = 0.0 # type: ignore[attr-defined] + return wrapper # type: ignore[return-value] + + return decorator + + +def diophantine_eq(D: int) -> int | None: + a = [int(isqrt(D))] + if a[0] * a[0] == D: + return None + p = [0, 1, int(isqrt(D))] + m = 0 + d = 1 + while True: + m = d * a[-1] - m + d = (D - m * m) // d + a.append((a[0] + m) // d) + p.append(p[-1] * a[-1] + p[-2]) + if a[-1] == 2 * a[0]: + if len(a) % 2 == 1: + return p[-2] + else: + for ai in a[1:-1]: + p.append(p[-1] * ai + p[-2]) + return p[-1] + + +@benchmark(repeat=10) +def find_minimal_solution(limit: int) -> int: + min_x: dict[int, int] = {} + for D in range(2, limit + 1): + x = diophantine_eq(D) + if x is not None: + min_x[D] = x + return max(min_x, key=min_x.__getitem__) + + +if __name__ == "__main__": + max_D = 1000 + dd = find_minimal_solution(max_D) + print(f"D <= {max_D} |-> D of max x = {dd}")