二次剩余之 Cipolla 算法学习笔记

前言

曾经，在数个月以前，我学过一波二次剩余……

只不过当时是被 SaInT7 所胁迫，而且时间仓促，精神不集中，最终只解得皮毛。

现在我重新学了一遍，略作总结。

二次剩余问题

已知同余方程 $x^2\equiv n(\mathrm{m}\mathrm{o}\mathrm{d}p)$，给定 $n,p$，判断该方程是否有解，如果有解则输出所有非负整数解。

注意到 $n=0$ 的情况是平凡的，只有 $x=0$ 一解。如果该方程有解且 $n\ne 0$，则称 $n$ 为模 $p$ 意义下的二次剩余。

在此文中，只讲解 $p$ 为奇素数的一种解法：Cipolla 算法，其他情况的解法详见其他文章。

一些定理

以下均默认 $n\ne 0$，同余方程为模 $p$ 意义，$p$ 为奇素数。

定理 1

如果 $n$ 是二次剩余，则原方程恰有两不同解。

证：

令一解为 $x_1$，则 $p-x_1$ 显然也为一解。注意到 $p$ 是奇素数，因此 $x_1\ne p-x_1$。因此原方程至少有两解。

我们从原方程的解集中任选两不同解，记为 $u,v$。

则有

$$u^2-v^2\equiv 0$$

平方差得

$$(u+v)(u-v)\equiv 0$$

因为 $u\ne v$，所以 $u+v\equiv 0$。所以对于每个 $u$，$v$ 是唯一的，然而我们又是任选的 $u,v$，所以解集大小只能恰好为 2。

定理 2

对于一个同余方程组 $x^2\equiv n$，有解的 $n$ 恰有 $\frac{p-1}{2}$ 个。

证：

根据定理 1，无序模意义相反数对集合与二次剩余集合形成一一映射，而无序模 $p$ 意义相反数对恰有 $\frac{p-1}{2}$ 个，因此得证。

定理 3

$$(a+b{)}^p\equiv a^p+b^p$$

证：

二项式定理展开左侧，可得到系数分别为 $C_p^0,C_p^1,\cdots ,C_p^p$。

展开组合数，发现除 $C_p^0$ 项与 $C_p^p$ 项外，每项的分子上均有乘子 $p$，而又因为 $p$ 是奇素数无法约去，因此这些项都是 $p$ 的倍数，即与 0 同余，不考虑。

所以只需考虑第一项与最后一项，则证毕。

定理 4（欧拉判定准则）

$n^{\frac{p-1}{2}}\equiv 1$ 与 $n$ 是二次剩余等价，$n^{\frac{p-1}{2}}\equiv p-1$ 与 $n$ 是非二次剩余等价。

证：

若 $n$ 是二次剩余:

设 $x$ 为一解，费马小定理有 $x^{p-1}\equiv 1$，因此 $(x^2{)}^{\frac{p-1}{2}}\equiv n^{\frac{p-1}{2}}\equiv 1$。

若 $n^{\frac{p-1}{2}}\equiv 1$：

设 $n\equiv g^k$，其中 $g$ 为模 $p$ 意义下的原根。则 $g^{k\times \frac{p-1}{2}}\equiv 1$。因此 $p-1|k\times \frac{p-1}{2}$，得 $k$ 为偶数。那么 $g^{\frac{k}{2}}$ 显然为一解，因此 $n$ 是二次剩余。

根据费马小定理，$n^{\frac{p-1}{2}}\equiv 1$ 或 $n^{\frac{p-1}{2}}\equiv p-1$。对于一个 $n$ 只存在其 是/不是 二次剩余两种情况，第一种与 $n^{\frac{p-1}{2}}\equiv 1$ 等价，则第二种必与 $n^{\frac{p-1}{2}}\equiv p-1$ 等价。

Cipolla 算法与说明

Cipolla 算法

对于一组 $n,p$，先用欧拉判定准则判断无解，并特判 $n=0$。

找到一个数 $a$ 使得 $a^2-n$ 不是二次剩余（用欧拉判定准则判断）。根据定理 2，可以随机选取 $a$，期望 2 次即能找到非二次剩余。

然后我们扩域，令单位 $i$ 满足 $i^2\equiv a^2-n$。

$(a+i{)}^{\frac{p+1}{2}}$ 为一个解，其模意义相反数为另一个解。

说明

一、取 $(a+i{)}^{\frac{p+1}{2}}$ 为解的正确性。

我们只需证明 $(a+i{)}^{p+1}\equiv n$。

运用定理 3，等价于 $(a+i)(a^p+i^p)\equiv n$。

然后发现 $i^p\equiv (i^2{)}^{\frac{p-1}{2}}\times i\equiv (a^2-n{)}^{\frac{p-1}{2}}\times i\equiv p-i$。

再运用费马小定理，得知 $a^p+i^p\equiv a-i$。

因此 $(a+i)(a^p+i^p)\equiv (a+i)(a-i)\equiv a^2-a^2+n\equiv n$。

二、为什么 $(a+i{)}^{\frac{p+1}{2}}$ 的结果可以输出，即结果不含 $i$？

把结果用 $A+Bi$ 表示，则 $(A+Bi{)}^2\equiv n$，即 $A^2+2ABi+B^2{a}^2-B^{2}n\equiv n$。

式子右侧虚部系数为 0，因此 $AB\equiv 0$。若 $B\not\equiv 0$，则 $A\equiv 0$，即 $(Bi{)}^2\equiv n$。 化成 $i^2\equiv n{B}^{-2}$。

因为 $B^{-2}$ 是二次剩余（$B^{-1}$ 肯定为一解），而 $n$ 又是二次剩余，所以 $n{B}^{-2}$ 也是二次剩余，即 $i^2$ 是二次剩余，矛盾！

因此 $B\equiv 0$。

三、怎么实现？

大可以参照复数类的实现方法，实现一个“扩域类”。

不过据说也有不需要的写法？

四、时空复杂度？

时间复杂度期望 $\Theta (\log p)$，空间复杂度 $O(1)$。

代码实现

inline ll mul(ll x,ll y)
{
	ll ret=0;
	while(y)
	{
		if(y&1) ret=(ret+x)%p;
		x=(x+x)%p;
		y>>=1;
	}
	return ret;
}
inline ll qpow(ll x,ll y)
{
	ll ret=1;
	while(y)
	{
		if(y&1) ret=mul(ret,x);
		x=mul(x,x);
		y>>=1;
	}
	return ret;
}
struct Extend
{
	ll r,v;
	Extend(ll r=0,ll v=0): r(r),v(v) {}
	Extend operator * (const Extend &x) const
	{return Extend((mul(r,x.r)+mul(mul(v,x.v),v2))%p,(mul(r,x.v)+mul(v,x.r))%p);} 
};
inline Extend qpow(ext x,ll y)
{
	Extend ret(1,0);
	while(y)
	{
		if(y&1) ret=ret*x;
		x=x*x;
		y>>=1;
	}
	return ret;
}
bool check(ll x){return qpow(x,(p-1)>>1)==1;}
ll sqroot(ll x)
{
	ll tmp=mul(mul(rand(),rand()),mul(rand(),rand()));
	while(check((mul(tmp,tmp)-x+p)%p)) tmp=mul(mul(rand(),rand()),mul(rand(),rand()));
	v2=(mul(tmp,tmp)-x+p)%p;
	return qpow(Extend(tmp,1),(p+1)>>1).r;
}