计算数论的方法论:两平方定理与复杂度
从Fermat 两平方定理出发,说明存在性证明、格约化、位复杂度和最终精确验证之间的关系。
1. 计算问题的尺度
本课程不是把“理论”与“计算”对立起来,而是研究怎样把一个存在性证明改造成可执行、可验证并具有复杂度控制的算法。输入一个整数 时,输入长度是 ,所以遍历到 的算法对输入长度而言是指数时间。
设 为奇素数。则 有整数解,当且仅当 。解在交换 以及改变符号的意义下唯一。
证明Gaussian 整数证明
在Euclidean 整环 中,范数为 。若 ,则 ,故 是模 的平方,即 。反之,若 ,存在 ,于是 ,但 在 中不再是素元。唯一分解给出 。唯一性来自 的唯一分解,单位只有 。
2. 从证明到算法
证明真正给出的计算接口是:先求一个模 的 的平方根,再把同余关系转化为二维格中的短向量。二维格约化等价于Euclidean 算法,因此整体成本只需要多项式数量级的位运算。
原讲义取 ,计算得到
这个例子强调:先选 再构造 会掩盖算法困难;先给定素数再恢复表示才是实际问题。
3. 线性时间与二次时间
对部分和 ,若每次重新求和,工作量为 ;维护累加器只需 。在数论计算中,同样的差别常出现在重复求 gcd、重复展开级数以及没有缓存的模运算中。
算法部分和的正确实现
s = 0
for n in range(1, N+1):
s += a[n]
plot(n, s)补充复杂度记号与可验证计算
本讲义统一使用位复杂度而非“算术操作次数”。当整数有 位时,一次乘法的成本记为 。任何数值猜测都必须在最后回到精确算术:代入恒等式、检查最小多项式、验证分歧型或核对群作用。
从存在性到可执行算法
课程首先说明其总体目标与精神:所关心的并不只是某个数学对象“存在”,而是如何把存在性证明改写成能够实际执行、能够估计复杂度、并且最终能够用精确算术核验的计算过程。为了在不依赖后续大量背景知识的情形下展示这种思想,原讲义选择了Fermat著名的两平方定理作为第一个例子。
两平方定理断言:素数 能写成两个不同整数平方之和,当且仅当 ;这种表示除交换两个平方项以及改变符号外是唯一的。原讲义取
于是定理保证丢番图方程 有本质唯一的解。真正的计算问题是:怎样把这组 找出来?逐一试验所有 的做法虽然在形式上会在有限时间内终止,但这个“有限”远远不够小;即使计算机速度提高,也可以随时把输入的位数加倍,使这种穷举继续失效。
一个经典证明几乎直接产生高效算法。Cornacchia 在 1908 年提出的思想是:若 ,则 在模 意义下是 的平方根;反过来,给定 ,可以借助二维格约化恢复 。二维格约化本质上就是Euclidean 算法,因此所需时间可控制在 。这里“多项式时间”的变量不是 本身,而是描述 所需的位数,约为 。
原讲义指出,PARI/GP 等软件已经实现了所需的模平方根与 LLL 约化。Victor Miller 在 1992 年给出的一行程序,经改写为较新的 GP 语法后可写成:
代码PARI/GP 实现
fermat(p)=qflll([lift(sqrt(Mod(-1,p))),p;1,0])[1,]对上述大素数调用该程序,几乎瞬间返回
直接平方并相加即可精确验证结果。原讲义特别解释,选择 的前若干位并没有数论上的特殊意义;这样选只是为了避免先挑好 再反向构造 的“作弊”。这个规模的素数足够密集,找到例子并不困难。至于如何知道该整数确实是素数,原讲义说明:对这种位数,整数分解与素性证明早已成为常规工具;有限域或数域上的多项式分解也同样非平凡但已经标准化。尤其值得注意的是,LLL 算法最初的重要应用之一,正是证明 中多项式分解可在多项式时间内完成。
“几乎”高效的原因:模平方根
上面的论证为何只说“几乎”得到高效算法?因为还必须求出模 的 的平方根。计算Legendre 符号很容易,它能判断一个非零元素是否为平方;然而在不作额外假设时,一般并不知道如何确定性地在多项式时间内从“它是平方”构造平方根。若假设与模 Legendre特征相关的广义 Riemann 假设,则可得到确定性多项式时间算法;无条件情形下则有随机多项式时间算法。寻找一个二次非剩余已经足够,而且“找二次非剩余”与“求平方根”在多项式时间归约意义下彼此等价。
原讲义随后给出一个看似迂回但理论上重要的事实:对固定的小整数,模 平方根可借助椭圆曲线算术在确定性多项式时间内求出,尽管这种方法并不实用。Schoof 在 1985 年关于有限域椭圆曲线点数计算的论文中,把这一点作为算法的应用之一。对本例,考虑具有复乘的椭圆曲线
它的复乘环含有 。当 时,点数具有 或 的形式;因此确定性地计算 后,便能恢复两平方表示。
一个基础但关键的复杂度插曲
原讲义提醒:即便极为常规的计算也可能隐藏不必要的低效率。设要画出实数列 的部分和
对应的 个点 。若在每个 处重新计算一次 ,总工作量约为 ;维护一个累加变量只需约 次操作。低于线性时间通常不可能,因为仅仅读取或生成全部 就需要 时间。原讲义说,这个错误虽看似不应出现在研究生课程中,却确实来自一位知名计算数论学者的伪代码;因此任何人都可能在“显而易见”的地方犯复杂度错误。
课程的第一个主任务:显式曲线覆盖
接下来若干周的推动问题,是计算具有给定分歧的曲线覆盖。设 是紧Riemann 曲面之间次数 的映射, 是有限的分支点集。给定 ,可能的 只有有限多个;拓扑上它们对应于 的指数 子群。紧Riemann 曲面又等价于复数域上的光滑射影代数曲线,所以问题可以改述为:给定曲线 、有限点集 与整数 ,构造全部曲线 及其到 的映射。
若 定义在 上,则每个覆盖在某个有限扩张 上定义。许多数论和代数几何中的计算问题,都能编码成从 以及可能的组合数据——例如 ——恢复 的问题。即使 、分支点数与次数都不大,这也已经相当困难。由基本群子群得到覆盖的拓扑构造,本质上依赖Riemann 存在定理,因而具有“超越性”:它能给出例如定义域次数的上界,却不会自动给出代数方程。把这种拓扑存在性转化为显式方程,正是课程最初的一系列目标。
补充:两平方定理的多种证明与算法化差异
两平方定理有多种证明,而它们提供的算法信息并不相同。Wilson 定理与二项式系数证明可以说明 是模 的平方,却不会直接给出较小的整数 。几何数论证明考察格 ,用 Minkowski 定理保证存在非零向量 满足 ;同余条件使 ,因而只能等于 。这条证明已经非常接近算法,因为“找短向量”可由约化实现。
Gaussian 整数证明则揭示唯一性。若 ,多项式 在 中分裂,故理想 的范数为 。在Euclidean 整环中求其生成元 ,得到 。若还有 ,则 与 只差单位与共轭,从而表示唯一。实现时可直接在 做 gcd;Cornacchia 只是把同一运算压缩到普通整数Euclidean 算法。
输入长度与位复杂度
设 有 位。试除所有 需要约 次循环,是指数时间。Euclidean 算法只需 次除法;若 是 位整数乘法成本,则快速 gcd 可在 位操作内完成。Tonelli–Shanks 的主要成本是 次模乘与幂运算。因而整个两平方恢复在随机多项式时间内完成。
复杂度分析还应区分“算出一个答案”与“证明输入条件”。若 未知是否为素数,先运行确定性 AKS 在理论上可行,但实践中通常用 Baillie–PSW 筛查加 ECPP 或 APR-CL 生成可验证素性证书。若输入是合数且仍满足 ,模平方根和格约化可能产生候选,但两平方表示的存在与唯一性不再由素数定理保证。
可复现实验
复现实验时应输出:输入整数、素性证书摘要、所用平方根 、Euclidean 余数序列、首次低于 的余数、最终整数平方检验,以及 的精确值。这样结果不依赖软件黑箱。对巨大整数,可另外给出十六进制或哈希,避免复制过程中的数字错误。
存在性证明的算法设计原则
本例体现四条贯穿全课程的原则。第一,找出证明中真正需要构造的中间对象,这里是模平方根。第二,把代数条件编码成一个具有几何意义的格或簇。第三,选择适合维数和系数规模的约化或提升算法。第四,数值阶段结束后回到精确恒等式。后续的 Belyi 映射、代数数恢复与模曲线公式都遵循同一模板。
进一步:表示数、素数分解与算法边界
两平方恒等式不仅证明可表示数对乘法封闭,还给出组合表示的算法。若已知 的素因子分解,则Fermat 两平方定理断言: 可写为两平方和,当且仅当每个 的指数为偶数。对 的素因子先求一组表示,再用 Brahmagupta–Fibonacci 恒等式合成; 的偶次幂只贡献实整数因子。
表示数的精确个数满足
其中按有序、有符号整数对计数, 分别统计模 4 同余于 1、3 的约数。这个公式来自 中的唯一分解,也可由 theta 级数比较证明。它提醒我们:“找出一个表示”和“枚举全部表示”是不同复杂度问题。
Fermat 大数的例子中,首先观察指数结构和已知因子远比盲目试除重要。若目标只是证明合数,一个非平凡因子即可;若要写成两平方和,则还需掌握所有模 3 素因子的奇偶性。现代整数分解通常分阶段使用试除、Pollard rho、ECM、二次筛和数域筛。算法选择由待找因子大小而非整数总位数单独决定。
复杂度陈述应明确输入模型。整数 的输入长度是 ,所以 的试除对数输入是指数时间;而多项式于 的算法才是位复杂度意义下的多项式时间。算术操作次数还要乘上大整数乘法成本 。忽略位长会把许多表面上“线性”的循环误判为高效算法。
完成计算后可独立验证:检查 的精确整数等式;对每个宣称的素因子做确定性或有证书的素性证明;核对因子乘积;若表示由Gaussian 整数相乘得到,检查共轭选择仅改变符号与次序而不改变范数。
练习与解答两平方与位复杂度
练习 1。证明若 且素数 整除 ,则 。
解。若 ,则 ,与 矛盾;故 ,同理 。因此这类素数在范数中只能出现偶指数。
练习 2。由Gaussian 整数唯一分解推导素数 的表示在符号与交换外唯一。
解。在 中 。任何 给出因子 ,只能与 或 相伴;单位 正对应符号与交换。
练习 3。若试除到 ,以输入位数 表示操作次数。
解。候选约为 ,故对位长是指数复杂度。每次余数运算还需要大整数成本,所以不能称为多项式时间。
深入补充:一般整数的两平方表示数与完整算法
Fermat定理只回答素数何时可写成两个平方;对合数,分解型同时控制存在性、表示数与算法结构。记 为有序带符号解 满足 的个数,并令 分别表示 的正因子中模 同余于 的个数。
对任意正整数 ,有
等价地,若
则当且仅当所有 都为偶数时 可表示为两个平方;在此情形
证明可直接在Gaussian 整数环中完成。素数 分裂为 ,每个指数 可在 与 之间以 种方式分配;素数 在 中仍不可约,若其指数为奇数便不可能成为范数。单位 给出最后的因子 。这个论证不仅给出存在性,还给出枚举全部表示的乘法结构。
有 ,故
除去交换与符号,正整数表示恰有三组:
三组各自产生八个有序带符号解,合计 个;这里不存在坐标为零或两坐标相等所造成的轨道缩短。
先检查每个 的指数是否为偶数。对每个 ,用模平方根加 Cornacchia 算法求一个Gauss素因子 。随后在
中枚举每个 的指数分配,乘出Gaussian 整数 ,再取 。为避免重复,可规定第一非零坐标为正且 ,最后按需要恢复完整轨道。
这个算法揭示了一个重要的复杂度边界:给定素因数分解后,判定与构造都可在输入长度和输出规模的多项式时间内完成;若只给出一般整数 ,最困难的步骤通常不是 Cornacchia,而是分解 。因此“两平方问题是容易的”必须说明输入模型:对素数输入可直接做模平方根,对已分解整数也容易,而对未分解合数,算法成本与整数分解紧密相连。
还可用 theta 级数把公式编码为生成函数。Jacobi 恒等式给出
其中 分别对应 为偶数、、。比较 系数便再次得到 。这说明格点计数、模形式与Gaussian 整数分解在同一公式中汇合。
补充计算:Jacobi 表示数公式、Gaussian 乘法与全整数证书
设 \(r_2(n)=\#\{(x,y)\in\mathbf Z^2:x^2+y^2=n\}\),其中次序与符号均计入。记 \(χ_4\) 为模 4 的非主特征,则 等价地,若 \(n=2^a∏p_i^{α_i}∏q_j^{β_j}\),其中 \(p_i≡1 (mod 4)\),\(q_j≡3 (mod 4)\),则某个 \(β_j\) 为奇数时 \(r_2(n)=0\);否则
证明 由Gaussian 整数的单位与共轭分解计数
在 \(Z[i]\) 中,2 与 1+i 伴随;每个 \(p≡1 (mod 4)\) 分裂为 \(π\bar π\);每个 \(q≡3 (mod 4)\) 仍为Gaussian 素数。若 q 在 n 中的指数为奇数,则任一范数为 n 的Gaussian 整数都要求 q 的指数在元素与共轭中各占一半,矛盾。故可表示性等价于所有 \(β_j\) 偶。
设所有 \(β_j\) 偶。范数为 n 的元素可写成单位 u、固定实因子 \(∏q_j^{β_j/2}\)、\((1+i)^a\),以及对每个 \(p_i\) 在 \(π_i\) 与 \(\bar π_i\) 之间分配 \(α_i\) 个指数的乘积。第 i 个分裂素数有 \(α_i+1\) 种分配,四个单位给出四倍;共轭与交换坐标已包含在单位和指数分配中。因此得到 \(4∏(α_i+1)\)。
另一方面,Dirichlet 卷积恒等式 \(1*χ_4\) 的素数幂值分别为 1、k+1、以及 k 为偶数时 1、奇数时 0,与上式逐素数一致,故 \(r_2=4(1*χ_4)\)。
Python 分解、构造、枚举与计数的统一核验程序
from collections import Counter
from math import isqrt
def factorint(n: int) -> Counter:
"""Trial division with wheel 2,3,6k±1; sufficient for certificates."""
if n < 1:
raise ValueError("n must be positive")
factors = Counter()
for p in (2, 3):
while n % p == 0:
factors[p] += 1
n //= p
f, step = 5, 2
while f * f <= n:
while n % f == 0:
factors[f] += 1
n //= f
f += step
step = 6 - step
if n > 1:
factors[n] += 1
return factors
def divisors_from_factorization(factors: Counter) -> list[int]:
divisors = [1]
for p, e in sorted(factors.items()):
old = divisors[:]
divisors = []
power = 1
for _ in range(e + 1):
divisors.extend(d * power for d in old)
power *= p
return sorted(divisors)
def r2_formula(n: int) -> int:
"""Number of ordered signed representations n=x^2+y^2."""
d1 = d3 = 0
for d in divisors_from_factorization(factorint(n)):
if d % 4 == 1:
d1 += 1
elif d % 4 == 3:
d3 += 1
return 4 * (d1 - d3)
def gaussian_mul(z: tuple[int, int],
w: tuple[int, int]) -> tuple[int, int]:
a, b = z
c, d = w
return a*c - b*d, a*d + b*c
def gaussian_pow(z: tuple[int, int], e: int) -> tuple[int, int]:
out = (1, 0)
while e:
if e & 1:
out = gaussian_mul(out, z)
z = gaussian_mul(z, z)
e >>= 1
return out
def prime_sum_two_squares_brutal(p: int) -> tuple[int, int]:
"""Small-prime base case used only to build exact certificates."""
for x in range(isqrt(p), -1, -1):
y2 = p - x*x
y = isqrt(y2)
if y*y == y2:
return x, y
raise ValueError(f"{p} has no representation")
def one_representation(n: int) -> tuple[int, int] | None:
"""
Construct one representation using Gaussian multiplication.
A prime q≡3 mod 4 must occur to even exponent.
"""
factors = factorint(n)
z = (1, 0)
scalar = 1
for p, e in sorted(factors.items()):
if p == 2:
z = gaussian_mul(z, gaussian_pow((1, 1), e))
elif p % 4 == 1:
base = prime_sum_two_squares_brutal(p)
z = gaussian_mul(z, gaussian_pow(base, e))
else:
if e & 1:
return None
scalar *= p ** (e // 2)
return abs(z[0] * scalar), abs(z[1] * scalar)
def enumerate_representations(n: int) -> set[tuple[int, int]]:
result = set()
bound = isqrt(n)
for x in range(-bound, bound + 1):
y2 = n - x*x
if y2 < 0:
continue
y = isqrt(y2)
if y*y == y2:
result.add((x, y))
result.add((x, -y))
return result
def certificate(n: int) -> dict:
factors = factorint(n)
candidate = one_representation(n)
enum = enumerate_representations(n)
return {
"n": n,
"factorization": dict(factors),
"formula": r2_formula(n),
"enumerated": len(enum),
"candidate": candidate,
"candidate_ok": candidate is None or
candidate[0]**2 + candidate[1]**2 == n,
"count_ok": r2_formula(n) == len(enum),
}
samples = [5, 13, 25, 65, 85, 125, 221, 325, 1105, 1885]
for n in samples:
c = certificate(n)
print(
f"n={n:4d} factors={str(c['factorization']):18s} "
f"r2={c['formula']:2d} rep={str(c['candidate']):10s} "
f"count_ok={c['count_ok']}"
)
print("\nfirst failures of representability below 80:")
bad = [n for n in range(1, 80) if one_representation(n) is None]
print(bad)
print("\nidentity checks:")
for a, b in [(5, 13), (13, 17), (25, 29), (65, 73)]:
z = one_representation(a)
w = one_representation(b)
zw = gaussian_mul(z, w)
print(f"{a}*{b}={a*b}: {z}·{w}={zw}, norm={zw[0]**2+zw[1]**2}")
print("\nJacobi-count and constructive-certificate table, 1 <= n <= 30:")
for n in range(1,31):
fac=factorint(n)
rep=one_representation(n)
r2=r2_formula(n)
enum_count=len(enumerate_representations(n))
fac_text="*".join(f"{p}^{e}" if e>1 else str(p)
for p,e in sorted(fac.items()))
if rep is None:
rep_text="--"
check=(r2==0 and enum_count==0)
else:
rep_text=f"{rep[0]}^2+{rep[1]}^2"
check=(rep[0]**2+rep[1]**2==n and r2==enum_count)
print(f" n={n:3d} factor={fac_text:17s} r2={r2:3d} "
f"enumerated={enum_count:3d} certificate={rep_text:13s} exact={check}")
运行结果 样本分解、表示数与范数恒等式
n= 5 factors={5: 1} r2= 8 rep=(2, 1) count_ok=True
n= 13 factors={13: 1} r2= 8 rep=(3, 2) count_ok=True
n= 25 factors={5: 2} r2=12 rep=(3, 4) count_ok=True
n= 65 factors={5: 1, 13: 1} r2=16 rep=(4, 7) count_ok=True
n= 85 factors={5: 1, 17: 1} r2=16 rep=(7, 6) count_ok=True
n= 125 factors={5: 3} r2=16 rep=(2, 11) count_ok=True
n= 221 factors={13: 1, 17: 1} r2=16 rep=(10, 11) count_ok=True
n= 325 factors={5: 2, 13: 1} r2=24 rep=(1, 18) count_ok=True
n=1105 factors={5: 1, 13: 1, 17: 1} r2=32 rep=(9, 32) count_ok=True
n=1885 factors={5: 1, 13: 1, 29: 1} r2=32 rep=(6, 43) count_ok=True
first failures of representability below 80:
[3, 6, 7, 11, 12, 14, 15, 19, 21, 22, 23, 24, 27, 28, 30, 31, 33, 35, 38, 39, 42, 43, 44, 46, 47, 48, 51, 54, 55, 56, 57, 59, 60, 62, 63, 66, 67, 69, 70, 71, 75, 76, 77, 78, 79]
identity checks:
5*13=65: (2, 1)·(3, 2)=(4, 7), norm=65
13*17=221: (3, 2)·(4, 1)=(10, 11), norm=221
25*29=725: (3, 4)·(5, 2)=(7, 26), norm=725
65*73=4745: (4, 7)·(8, 3)=(11, 68), norm=4745
Jacobi-count and constructive-certificate table, 1 <= n <= 30:
n= 1 factor= r2= 4 enumerated= 4 certificate=1^2+0^2 exact=True
n= 2 factor=2 r2= 4 enumerated= 4 certificate=1^2+1^2 exact=True
n= 3 factor=3 r2= 0 enumerated= 0 certificate=-- exact=True
n= 4 factor=2^2 r2= 4 enumerated= 4 certificate=0^2+2^2 exact=True
n= 5 factor=5 r2= 8 enumerated= 8 certificate=2^2+1^2 exact=True
n= 6 factor=2*3 r2= 0 enumerated= 0 certificate=-- exact=True
n= 7 factor=7 r2= 0 enumerated= 0 certificate=-- exact=True
n= 8 factor=2^3 r2= 4 enumerated= 4 certificate=2^2+2^2 exact=True
n= 9 factor=3^2 r2= 4 enumerated= 4 certificate=3^2+0^2 exact=True
n= 10 factor=2*5 r2= 8 enumerated= 8 certificate=1^2+3^2 exact=True
n= 11 factor=11 r2= 0 enumerated= 0 certificate=-- exact=True
n= 12 factor=2^2*3 r2= 0 enumerated= 0 certificate=-- exact=True
n= 13 factor=13 r2= 8 enumerated= 8 certificate=3^2+2^2 exact=True
n= 14 factor=2*7 r2= 0 enumerated= 0 certificate=-- exact=True
n= 15 factor=3*5 r2= 0 enumerated= 0 certificate=-- exact=True
n= 16 factor=2^4 r2= 4 enumerated= 4 certificate=4^2+0^2 exact=True
n= 17 factor=17 r2= 8 enumerated= 8 certificate=4^2+1^2 exact=True
n= 18 factor=2*3^2 r2= 4 enumerated= 4 certificate=3^2+3^2 exact=True
n= 19 factor=19 r2= 0 enumerated= 0 certificate=-- exact=True
n= 20 factor=2^2*5 r2= 8 enumerated= 8 certificate=2^2+4^2 exact=True
n= 21 factor=3*7 r2= 0 enumerated= 0 certificate=-- exact=True
n= 22 factor=2*11 r2= 0 enumerated= 0 certificate=-- exact=True
n= 23 factor=23 r2= 0 enumerated= 0 certificate=-- exact=True
n= 24 factor=2^3*3 r2= 0 enumerated= 0 certificate=-- exact=True
n= 25 factor=5^2 r2= 12 enumerated= 12 certificate=3^2+4^2 exact=True
n= 26 factor=2*13 r2= 8 enumerated= 8 certificate=1^2+5^2 exact=True
n= 27 factor=3^3 r2= 0 enumerated= 0 certificate=-- exact=True
n= 28 factor=2^2*7 r2= 0 enumerated= 0 certificate=-- exact=True
n= 29 factor=29 r2= 8 enumerated= 8 certificate=5^2+2^2 exact=True
n= 30 factor=2*3*5 r2= 0 enumerated= 0 certificate=-- exact=True
对 Re(s)>1, 素数 p≡1 (mod 4) 的 Euler 因子为 4 的全局常数之外的 (1-p^{-s})^{-2};素数 q≡3 (mod 4) 的因子为 (1-q^{-2s})^{-1};素数 2 的因子为 (1-2^{-s})^{-1}。对 n=1105=5·13·17,有 r_2(n)=4·2^3=32;程序枚举得到 32 个有序带符号解。