二次剩余

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

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

当 (n) 不为 (p) 的倍数时，如果存在 (x)，使得 (x^2 equiv n (mod 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,yle frac{p-1}{2},x^2 equiv y^2)，那么 ((x+y)(x-y) equiv 0)。又因为 (x+y,x-y) 均不为 (p) 的倍数，且 (x+y ot= 0,x-y ot= 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)（(nequiv g^t)）。若 (n^{frac{p-1}{2}} equiv 1)，那么 (g^{t imes frac{p-1}{2}} equiv 1)。根据原根定义，有 (t imes frac{p-1}{2} | p-1)，即 (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(mod p))。
• (i^p equiv -i(mod p))。证明：(i^pequiv (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(mod p))。证明考虑二项式定理，只有 (i=0) 和 (i=p) 的时候组合数为1.

然后开始正式推式子：

[egin{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
", ans); return ;}
	if (ans > ans_)	swap(ans, ans_);
	printf("%lld %lld
", ans, ans_);
}