Skip to content

算法Algorithm

Cornacchia 算法:两平方和构造

Cornacchia algorithm · Cornacchia sum-of-two-squares algorithm · Cornacchia 两平方和算法

给定−1的模平方根,以截断欧几里得余数链构造本原两平方和,用余数和系数双不变量证明平方根阈值后的结果恰好等于模数。

知道 325 能写成两平方和,并不意味着要从 02,12,22 一直试到接近325。一个更短的输入是同余证书 572+1=10⋅325。Cornacchia算法把325与57做两次带余除法,就得到 325=172+62。

本页完整证明 x2+y2=m 的版本。一般Cornacchia算法还处理 x2+dy2=m,但 d>1 的完整性需要进一步论证;文末会指出本页证明中不能直接照搬的一步。

形式陈述 ​

输入与输出 ​

输入整数 m>1,以及满足模同余

t2≡−1(modm)

的整数 t。先取 tmodm 的最小非负代表,再将它与 m−t 中较小者记为 t,于是 1≤t≤m/2。零不可能是模 m>1 的 −1 的平方根。

输出正整数 x,y,满足

x2+y2=m,gcd(x,y)=1.

模数不必为素数。算法只承诺为所给模平方根构造一份表示,不承诺一次输出所有表示。

余数链在何处停下 ​

设 r0=m,r1=t,执行整数欧几里得带余除法

rj−1=qjrj+rj+1,0≤rj+1<rj,

在第一个满足 rk2<m 的正余数处停止。此时令 x=rk,检查 m−x2 是否为整数平方;它必为平方,正平方根就是 y。

实现时可以同步维护系数,连最后一次开方也省去:令

s0=0,s1=1,sj+1=qjsj+sj−1.

在相同停止位置直接输出 (rk,sk)。保留整数平方检查仍有价值,因为它能发现错误输入、抄错的商或程序故障。

text
检查 m > 1 且 (t*t + 1) mod m = 0
将 t 约简到 1 ≤ t ≤ m/2
(a, b, u, v) ← (m, t, 0, 1)
while b*b ≥ m:
    q ← a div b
    (a, b, u, v) ← (b, a-q*b, v, u+q*v)
检查 b*b + v*v = m 且 gcd(b,v) = 1
返回 (b,v)

所有比较都是整数比较。代码中的元组更新使用旧值同时计算,不能先覆盖 a,b 再用新值更新另一个分量。

直觉

两个一起保存的不变量 ​

余数向下变小,系数向上积累。对每个 j≥1,有

rj−1sj+rjsj−1=m,rj≡(−1)j−1tsj(modm).

第一式在 j=1 时就是 m⋅1+t⋅0=m。若它在当前步成立,下一步的左边为

rj(qjsj+sj−1)+(rj−1−qjrj)sj=rjsj−1+rj−1sj=m.

所以它始终成立。第二式从 r1=t 开始,把相邻两个交替符号的同余代入 rj+1=rj−1−qjrj,就得到新的交替符号与系数 sj+1。

由输入 t2≡−1,第二个不变量立即给出

rj2+sj2≡0(modm).

也就是说,每一对余数和系数的平方和都是 m 的倍数。还需要一个大小界,才能把“某个倍数”缩成“恰好一次”。

第一次越过平方根阈值足够小 ​

输入同余保证 gcd(m,t)=1:任何共同素因子都会同时整除 t2 与 t2+1。因此普通Euclidean链的最后一个正余数是一,在到零之前一定已经满足 rk2<m,算法会以正余数停止。

由于选择的是第一个越过阈值的位置,有 rk−12≥m。第一不变量中的两项非负,于是

0<sk≤mrk−1≤m,0<rk<m.

所以

0<rk2+sk2<2m.

它又是 m 的整数倍,两条信息合起来只允许 rk2+sk2=m。证明中出现平方根只是为了说明界;执行时只需要比较整数平方。

为什么结果还一定互素 ​

从同余不变量可写 rk=hm+εtsk,其中 h∈Z,ε∈{1,−1}。记 Q=(t2+1)/m∈Z,将这个式子代入已经证明的 rk2+sk2=m,除以 m 得

1=h2m+2εhtsk+Qsk2.

若某素数同时整除 rk,sk,它也整除 m=rk2+sk2,因而整除右边全部三项,与左边为一矛盾。这一步证明了本原性,不是观察几个例子后额外猜出的性质。

例子与边界

325与57的完整轨迹 ​

先验证 572+1=3250=10⋅325。余数和系数为

j rj sj 生成信息
0 325 0 初值
1 57 1 初值
2 40 5 325=5⋅57+40
3 17 6 57=1⋅40+17

402>325 而 172<325,所以到第三项停止。最后的不变量可直接检查:

40⋅6+17⋅5=325,17≡57⋅6(mod325).

结果为 172+62=325。若输入另一个平方根 t=18,一开始就有 182<325,直接输出 (18,1)。

一份模根不等于全部表示 ​

模325的 −1 平方根恰为 18,57,268,307。符号配对后只有 18,57 两个输入,给出两类本原表示。可是 325 还有 (10,15) 这一类非本原表示,它不会从这个输入契约中产生。

原因是由本原表示反求模根时,需要求 y−1(modm),再令 t≡xy−1。本原表示保证 gcd(y,m)=1;对 (10,15),15 在模325下不可逆。这不是算法漏解,而是所求对象确实不同。要找非本原表示,可先提取共同因子的平方,再对较小模数处理。

不能无证明地把系数一换成d ​

令 d=5,m=7,t=3,则 t2≡−5(mod7)。余数步骤给出 7=2⋅3+1,系数为二,但

12+5⋅22=21=3⋅7,

不等于七。对于一般 d,同余仍能保证 r2+ds2 是 m 的倍数,可原来的“小于 2m”不再成立。事实上七也不能写成 x2+5y2:y=0 时要求七为平方,|y|=1 时要求二为平方,更大的 |y| 已超出七。

所以一般版本必须保留最终整除与平方检查,并另证相应的完整性。本页的两平方和证明不能替代那项工作。

推论与应用

模平方根怎样得到 ​

对 p≡1(mod4) 的奇素数,可先找一个二次非剩余 a。Euler判据给出 a(p−1)/2≡−1,所以

t=a(p−1)/4(modp)

就是所需平方根。例如 p=29,a=2,27≡12,122≡−1;算法随后由 29=2⋅12+5 输出 (5,2)。

若独立均匀地从非零模 p 元素取样,恰有一半是二次非剩余,因此找到它的期望试验次数为二;每次仍需要模幂计算。对一般合数,不能直接使用素数的Euler判据;可以将模根作为已核验输入,或在已知素因子分解时另行求根并组合。两种情况的成本承诺必须分开。

可复算成本与终点 ​

给定模根之后,算法使用不超过完整Euclidean算法的除法步数,即 O(log⁡m) 步。各步操作的是 O(log⁡m) 位整数,大整数除法和乘法不算固定成本;使用通常的Euclidean位复杂度分析可得这部分的二次位复杂度上界。这个预算不包含求模根或分解 m。

配套整数范数证书单元要求同时保留输入同余、每次商余数、系数递推和最终平方和。任何一条抄错,都能被局部等式或最终检查定位,而不必重新相信整个搜索过程。

参考资料
  • Andrew V. Sutherland,18.783 Elliptic Curves, Problem Set 2,2025,Problem 1,pp.1–2:Cornacchia的平方根停止阈值、一般模数的根分支及最终平方检查。
  • Francesc Fité、Andrew V. Sutherland,Sato–Tate groups of y²=x⁸+c and y²=x⁷−cx,2016,§3的Cornacchia算法段,PDF第6页(印刷页108):算法输入与给定模平方根后的位复杂度。本文针对 d=1 另给余数—系数双不变量及本原性证明。
  • Julius M. Basilla,On the solution of x²+dy²=m,Proceedings of the Japan Academy, Series A,80(5),2004,pp.40–41,doi:10.3792/pjaa.80.40:一般Cornacchia正确性的原始研究参考。本轮未直接取得论文正文,不据此宣称已核对一般 d 的证明细节。
关系图谱5 个相邻概念 · 2 类关系

拖动节点调整位置。

显示关系

显示:依赖

  1. 前置三跳
  2. 前置二跳
  3. 前置一跳
  4. 当前条目
  5. 后续一跳
  6. 后续二跳
  7. 后续三跳
文字版关系按与当前条目的最短距离分组
类型化关系