二次剩余

给定 \(n,p\)，求 \(x^2 \equiv n (\bmod p)\) 的解 \(x\)。\(5 \times 10^5\) 组询问。保证 \(p\) 为奇素数。（即模意义下开根号）

做法简述：

1. 通过随机 \(a\) 找到 \((a^2-n)^{\frac{p-1}{2}}=-1 \pmod p\) 的 \(a\)
2. 根据 \(i^2 = a^2-n\) 建立复数域
3. \((a+i)^{\frac{p+1}{2}}\) 即为 \(\sqrt n\)。可以保证其没有虚部。

全文参考oi-wiki 以及 这篇博客

当 \(n\) 不为 \(p\) 的倍数时，如果存在 \(x\)，使得 \(x^2 \equiv n (\bmod p)\) ，那么我们称 \(n\) 为模 \(p\) 的二次剩余；如果不存在，则称 \(n\) 为模 \(p\) 的非二次剩余。

暂不考虑 \(n = 0\) 的情况，且仅考虑 \(p\) 为奇素数的情况。

二次剩余的个数

模 \(p\) 的二次剩余共有 \(\frac{p - 1}{2}\)，非二次剩余也有 \(\frac{p-1}{2}\) 个。

证明：\(0...\frac{p-1}{2}\) 的平方对应互不相同的所有二次剩余，且 \(i,p-i\) 二次剩余相同。考虑证明 \(0...\frac{p-1}{2}\) 的平方互不相同。若存在 \(x,y\le \frac{p-1}{2},x^2 \equiv y^2\)，那么 \((x+y)(x-y) \equiv 0\)。又因为 \(x+y,x-y\) 均不为 \(p\) 的倍数，且 \(x+y \not= 0,x-y\not= 0\)，所以不成立，原命题成立。

欧拉判别准则

若 \(n\) 为 \(p\) 的二次剩余，那么 \(n^{\frac{p-1}{2}} \equiv 1\)；若 \(n\) 为 \(p\) 的非二次剩余，那么 \(n^{\frac{p-1}{2}} \equiv -1\)。

证明：

根据费马小定理，\((n^{\frac{p-1}{2}}-1)(n^{\frac{p-1}{2}}+1) \equiv 0\)。因此 \(n^{\frac{p-1}{2}} \equiv 1 或 -1\) 。显然若 \(n\) 为二次剩余，即存在 \(x^2 \equiv n\) ，那么 \(x^{p-1} \equiv 1\)。

找到 \(p\) 的原根 \(g\)，以及 \(n\) 的对数 \(t\)（\(n\equiv g^t\)）。若 \(n^{\frac{p-1}{2}} \equiv 1\)，那么 \(g^{t \times \frac{p-1}{2}} \equiv 1\)。根据原根定义，有 \(p-1 | t \times \frac{p-1}{2}\)，即 \(t\) 必定为偶数。那么 \(g^{\frac{t}{2}}\) 即为一个解。

Cipolla 算法

首先随机出一个 \(a\) 使得 \(a^2-n\) 为 \(p\) 的非二次剩余。由于可行解很多，期望随机次数为 \(2\)。

定义“复数域”，类似真正的复数域，只不过满足 \(i^2=a^2-n\)。每个数可以表示成 \(A+Bi\) 的形式。

那么有结论:

\[\large \sqrt n = (a+i)^{\frac{p+1}{2}} \]

（注意是 \(p+1\)）

证明：

首先摆出三个比较显然的结论：

• \(x^p \equiv x(\bmod p)\)。
• \(i^p \equiv -i(\bmod p)\)。证明：\(i^p\equiv (i^2)^\frac{p-1}{2}i \equiv (a^2-n)^\frac{p-1}{2}i \equiv -i\)
• \((x+y)^p \equiv x^p+y^p(\bmod p)\)。证明考虑二项式定理，只有 \(i=0\) 和 \(i=p\) 的时候组合数为1.

然后开始正式推式子：

\[\begin{aligned} (a+i)^\frac{p+1}{2} &\equiv ((a+i)^{p+1})^{1/2}\\ &\equiv ((a+i)^{p}(a+i))^{1/2}\\ &\equiv ((a^p+i^p)(a+i))^{1/2}\\ &\equiv ((a-i)(a+i))^{1/2}\\ &\equiv (a^2-a^2+n)^{1/2}\\ &\equiv \sqrt n \end{aligned} \]

复杂度：单次 \(O(\log p)\)。

inline bool check(ll a) {
	return quickpow(a, (P - 1) >> 1) == P - 1;
}
struct Complex {
	ll x, y;
	Complex(ll xx = 0, ll yy = 0) { x = xx, y = yy; }
	inline Complex operator *(const Complex &a) const {
		return Complex((x * a.x + y * a.y % P * w) % P, (x * a.y + y * a.x) % P);
	}
};
inline Complex Pow(Complex x, ll k)
inline void work() {
	read(n), read(P);
	n = (n % P + P) % P;
	if (!n) {
		puts("0");
		return ;
	}
	if (quickpow(n, (P - 1) >> 1) != 1) { puts("Hola!"); return ; }
	a = Rand(P), w = (a * a - n + P) % P;
	while (!check(w))	a = Rand(P), w = (a * a - n + P) % P;
	Complex t(a, 1);
	ll ans = Pow(t, (P + 1) >> 1).x;
	ll ans_ = P - ans;
	if (ans == ans_) { printf("%lld\n", ans); return ;}
	if (ans > ans_)	swap(ans, ans_);
	printf("%lld %lld\n", ans, ans_);
}