从两平方和分解到拉格朗日四平方和定理

on 2025-10-02

起因

最近刷到一道ctf题目，例题简化如下

import Crypto.Util.*

p = getPrime(512)
q = getPrime(512)
n0 = p*p + q*q
print(n0)

核心考点在于将分解n0分解为p，q。

思路1

很直观的来看，我们可以进行以下推导

$$ q * q = n_0 - p * p $$

$$ q = \sqrt{n_0 - p * p} $$

写成Python代码就是

from gmpy2 import isqrt, is_square

for p in range(2, n0):
    qq = n0 - p * p
    if is_square(qq):
        q = isqrt(qq)
        print(f"p = {p}, q = {q}")
        break

但是这里的复杂度是O(n)，n0是一个1024bit的数，也就是我们要计算$2^{1024}$次！

然后我们可以想到，那p和q都是素数。我们的搜索下界是$2^{511}$，上界是$2^{513}$。然后有素数定理（Prime Number Theorem）可以计算素数个数：

$$ \pi(x) \approx Li(x) = \int_{2}^{x} \frac{dt}{ln(t)} $$

简化估算可得，在$2^{512}$附近大约每355个数就有一个素数。即使只枚举素数，复杂度仍然是$O(2^{512} / 355) \approx O(2^{504})$，依然不可行。

思路2：两平方和分解

既然暴力枚举不可行，我们需要利用数论中的定理。

我们的问题是 $n_0 = p^2 + q^2$ 的一个特殊情况。虽然拉格朗日四平方和定理告诉我们任何正整数都可以表示为至多四个平方数之和，但我们这里只需要两个平方数。

对于我们的问题：$n_0 = p^2 + q^2$ 显然可以表示为两平方和（这是构造方式），关键在于如何高效地找到 $p$ 和 $q$。

分解算法

我们可以使用以下方法来分解 $n_0$：

高斯整数分解法：利用高斯整数环的性质
Cornacchia 算法：专门用于求解 $x^2 + y^2 = n$ 的算法

下面介绍实用的 Cornacchia 算法。

Cornacchia 算法

Cornacchia 算法用于求解方程 $x^2 + dy^2 = n$，其中 $d$ 是正整数。对于我们的问题，$d = 1$。

算法步骤

检查 $n$ 是否为素数且 $n \equiv 1 \pmod{4}$（或 $n = 2$）
找到 $x_0$ 使得 $x_0^2 \equiv -1 \pmod{n}$
使用欧几里得算法：
- 令 $a = x_0, b = n$（即初始时 $b = n$）
- 重复执行 $a, b = b \bmod a, a$ 直到 $a < \sqrt{n}$
检查 $(n - a^2)$ 是否为完全平方数

为什么这个算法有效？

这是一个基于高斯整数理论的优美结果。

核心原理

如果 $n = a^2 + b^2$，且 $r^2 \equiv -1 \pmod{n}$，则在高斯整数环 $\mathbb{Z}[i]$ 中：

$$ n = (a + bi)(a - bi) $$

其中 $r^2 + 1 \equiv 0 \pmod{n}$ 意味着：

$$ (r + i)(r - i) \equiv 0 \pmod{n} $$

关键洞察：$\gcd(n, r + i)$ 在高斯整数中是 $n$ 的一个非平凡因子，恰好形式为 $(a + bi)$ 或其共轭。

为什么辗转相除有效

在整数上对 $(n, r)$ 做欧几里得算法，实际上是在模拟高斯整数中的 $\gcd(n, r+i)$：

每步 $n = qr + s$ 对应高斯整数中 $n = q(r+i) + (s-qi)$
余数的模 $|s-qi| = \sqrt{s^2 + q^2}$ 严格递减
当首次出现 $s \leq \sqrt{n}$ 时，必有 $s^2 + q^2 \approx n$
此时 $s$ 就是 $a$ 或 $b$

严格证明要点

引理：若 $\gcd(n, r+i) = a + bi$，则 $|a+bi|^2 = a^2 + b^2$ 整除 $n$。

在高斯整数环 $\mathbb{Z}[i]$ 中，由于 $r^2 \equiv -1 \pmod{n}$，我们有：

$$ n \mid (r^2 + 1) = (r+i)(r-i) $$

而 $n = (a+bi)(a-bi)$，所以 $\gcd(n, r+i)$ 必然包含 $n$ 的某个因子。欧几里得算法不断缩小余数的模，当首次出现余数 $s$ 满足 $s < \sqrt{n}$ 时，如果 $n - s^2$ 是完全平方数，则找到了分解。

这个过程的正确性来自：辗转相除在高斯整数中找到了 $\gcd(n, r+i)$ 的整数部分投影，而这个投影恰好对应 $n$ 在实部和虚部上的分解。

效率分析：

这个方法的效率来自：

不需要试除
只需 $O(\log n)$ 次辗转相除即可提取平方根
基于高斯整数的优雅性质

复杂度为 $O(\log^2 n)$（包括模运算），远优于暴力枚举的 $O(\sqrt{n})$。

Python 实现：

from gmpy2 import isqrt, is_square, powmod, invert
import random

def tonelli_shanks(n, p):
    """求解 x^2 ≡ n (mod p)"""
    if powmod(n, (p - 1) // 2, p) != 1:
        return None

    # 找 Q 和 S 使得 p - 1 = Q * 2^S
    Q, S = p - 1, 0
    while Q % 2 == 0:
        Q //= 2
        S += 1

    # 找一个非二次剩余 z
    z = 2
    while powmod(z, (p - 1) // 2, p) != p - 1:
        z += 1

    M = S
    c = powmod(z, Q, p)
    t = powmod(n, Q, p)
    R = powmod(n, (Q + 1) // 2, p)

    while True:
        if t == 0:
            return 0
        if t == 1:
            return R

        # 找最小的 i 使得 t^(2^i) = 1
        i = 1
        temp = (t * t) % p
        while temp != 1:
            temp = (temp * temp) % p
            i += 1

        b = powmod(c, 1 << (M - i - 1), p)
        M = i
        c = (b * b) % p
        t = (t * c) % p
        R = (R * b) % p

def cornacchia(n):
    """使用 Cornacchia 算法求解 x^2 + y^2 = n"""
    if n % 4 == 3:
        return None

    # 对于合数，先尝试找 x0 使得 x0^2 ≡ -1 (mod n)
    # 这里使用启发式方法
    x0 = tonelli_shanks(n - 1, n) if n > 2 else None

    if x0 is None:
        # 如果 n 不是素数，尝试随机搜索
        for _ in range(1000):
            x0 = random.randint(1, n - 1)
            if (x0 * x0) % n == n - 1:
                break
        else:
            return None

    # 确保 x0 > n/2
    if x0 < n // 2:
        x0 = n - x0

    # 欧几里得算法
    a, b = n, x0
    sqrt_n = isqrt(n)

    while a > sqrt_n:
        a, b = b, a % b

    # 检查 n - a^2 是否为完全平方数
    c2 = n - a * a
    if is_square(c2):
        c = isqrt(c2)
        return (int(a), int(c))

    return None

# 使用示例
def solve_two_squares(n0):
    """分解 n0 = p^2 + q^2"""
    result = cornacchia(n0)
    if result:
        p, q = result
        print(f"找到分解: {n0} = {p}^2 + {q}^2")
        print(f"p = {p}")
        print(f"q = {q}")
        print(f"验证: {p}^2 + {q}^2 = {p*p + q*q}")
        return p, q
    else:
        print("未找到分解")
        return None

if __name__ == "__main__":
    # 小数测试
    n = 5  # 5 = 1^2 + 2^2
    solve_two_squares(n)

更简单的方法：使用 SageMath

如果你有 SageMath 环境，可以直接使用内置函数：

# SageMath
result = two_squares(n0)
print(result)  # (p, q)

拉格朗日四平方和定理

拉格朗日四平方和定理（Lagrange's four-square theorem）：

任何正整数都可以表示为至多四个平方数之和。

$$ n = a^2 + b^2 + c^2 + d^2 $$

这个定理更加一般，对所有正整数都成立。分解算法也类似，可以使用：

Rabin-Shallit 算法：用于四平方和分解
递归分解：利用恒等式进行分解

在这里不给出相应的Python实现。在SageMath中，可以通过相关函数实现。