二次剩余 | Quadratic Residue

求解方程 $x^2\equiv a(\bmod p)$。

如果该方程有正整数解，称 $a$ 是 $p$ 的二次剩余。以下设 $p$ 是及素数。

解的数量

有正整数解的 $a$ 有 $\dfrac{p-1}{2}$ 个。这里的 $a\in[1,p-1]$。

勒让德符号 | Legendre Symbol

$\left(\dfrac{n}{p}\right)=\begin{cases}0, &n \mid p\\-1, &n\nmid p, \nexists x\in \mathbb Z, x^2\equiv n(\bmod p)\\1, &n\nmid p, \exists x\in \mathbb Z, x^2\equiv n(\bmod p)\end{cases}$

欧拉准则

$\left(\dfrac{n}{p}\right)=1 \iff n^{\frac{p-1}{2}}\equiv 1(\bmod p)$

欧拉准则的证明

根据费马定理知 $n^{p-1}\equiv 1(\bmod p)$。 那么显然 $n^\frac{p-1}{2}\equiv \pm1(\bmod p)$。因为把 $1$ 移过去，因式分解可知。

假设存在 $x^2\equiv n(\bmod p)$。那么 $n^\frac{p-1}{2}\equiv x^{p-1}\equiv 1(\bmod p)$，根据定义 $\left(\dfrac{n}{p}\right)=1$。

证毕。

Cipolla 算法

找到一个 $a$ 使得 $\left(\dfrac{a^2-n}{p}\right)=-1$。

定义 $i^2=a^2-n(\bmod p)$，则 $(a+i)^{\frac{p+1}{2}}$ 即为答案。

证明

$x\equiv (a+i)^\frac{p+1}2\equiv ((a+i)^{p+1})^\frac12\equiv ((a+i)^p(a+i))^\frac12(\bmod p)$。

注意到 $(a+i)^p\equiv a^p+i^p$，原因是用二项式定理展开，除了 $a^p$ 和 $i^p$ 两项，其余都存在因子 $p$。

然后呢根据费马定理的推论，$a^p\equiv a(\bmod p)$，并且我们求出了 $i^p\equiv (a^2-p)^\frac{p-1}{2}\cdot i=-1\cdot i=(\bmod p)$，所以 $x\equiv [(a-i)(a+i)]^\frac{1}{2}\equiv (a^2-i^2)^\frac{1}{2}\equiv(a^2-(a^2-n))^\frac{1}{2}\equiv n^\frac{1}{2}(\bmod p)$。

证毕。 ​

实现

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
ll i2, n, p;
struct Complex{
	ll r, i;
	Complex(ll r = 0, ll i = 0): r(r), i(i) {}
};
Complex operator * (Complex x, Complex y){
	Complex ans = Complex(0, 0);
	ans.r = (x.r * y.r % p + x.i * y.i % p * i2 % p + p) % p, ans.i = (x.r * y.i % p + x.i * y.r % p + p) % p;
	return ans;
}
ll qpow(Complex a, ll b, ll p){
	Complex ans = Complex(1, 0); 
	while(b){if(b & 1) ans = ans * a; a = a * a, b >>= 1; }
	return ans.r % p;
}
ll qpow(ll a, ll b, ll p){
	a %= p; ll ans = 1;
	while(b){if(b & 1) ans = ans * a % p; a = a * a % p, b >>= 1; }
	return ans;
}
inline ll rnd(){
	return 1ll * rand() * rand() % p; 
}
ll Cipolla(ll n, ll p){
	assert(p % 2);
	if(n != 0 && qpow(n, (p - 1) >> 1, p) == (p - 1)) return -1;
	ll x = 0; do x = rnd(); while(qpow(i2 = (x * x - n + p) % p, (p - 1) >> 1, p) != (p - 1));
    return (qpow((Complex){x, 1}, (p + 1) >> 1, p) + p) % p;
}
int T;
int main(){
	srand(54255 * 355 + 292);
	cin >> T;
	while(T--){
		cin >> n >> p; if(n == 0){cout << 0; continue;}
	    ll ans = Cipolla(n, p), ans2 = p - ans;
	    if(ans == -1) {cout << -1 << endl; continue; }
	    else{
	    	cout << 1 << ' ';
	    	if(ans > ans2) cout << ans2 << ' ' << ans << endl;
	    	else cout << ans << ' ' << ans2 << endl;
		}
	}
}