二次剩余（Cipolla's Algorithm）

当 \(a=0\) 或者 \(p=2\) 的情况，特判，只考虑奇素数的情况。

若 \(x^2\equiv a\pmod p\) 有解，那么 \(a\) 就是 \(p\) 的二次剩余，否则为非二次剩余。

定义勒让德符号 \(\left(\frac{a}{p}\right)=\begin{cases}1,a是p的二次剩余\\-1,a是p的非二次剩余\end{cases}\)

欧拉判别法：\(\left(\frac{a}{p}\right)=a^{\frac{p-1}{2}}(\gcd(a,p)=1)\)

如果 \(x^2\equiv a\pmod p\) 有解，那么 \(x^{p-1}\equiv 1\pmod p\)，那么 \(a^{\frac{p-1}{2}}\equiv 1\pmod p\)

如果无解，考虑将 \(1\sim p-1\) 的数两两配对，使得\(\forall x,y\in[1,p),xy=a\)，如果有解，一定会有 \(x\) 与 \(p-x\) 无法配对。

既然配成了 \(\frac{p-1}{2}\) 对，它们的乘积就是 \((p-1)!\equiv -1\pmod p\) 所以 \(a^{\frac{p-1}{2}}\equiv -1\pmod p\)

然后有两个结论用来证明做法的正确性。

1. 若 \(x^2\equiv a\pmod p\) 有解，那么 \(a\) 的个数为 \(\frac{p-1}{2}\)，相当于 \(a\) 与 \(p-a\) 一个有解，一个无解。
2. \((a+b)^p\equiv a^p+b^p\pmod p\)，直接二项式定理展开 \(\binom{p}{i}\) 当且仅当 \(i=0\) 或者 \(i=p\) 时 \(p\) 才能被消去，否则都会被 \(p\) 取模为0。

第一条结论说明随机出一个非二次剩余的期望次数为两次。

考虑随机出一个数 \(A\)，使得 \(A^2-a\) 是一个非二次剩余，也就是 \((A^2-a)^{\frac{p-1}{2}}\equiv -1\pmod p\)

令 \(\omega=\sqrt{A^2-a}\)，那么它应该是一个虚数，所有的复数都能用 \(x+y\omega\) 表示。

最后的答案 \(x\equiv (A+\omega)^{\frac{p+1}{2}}\pmod p\)

证明就是 \(x^2\equiv (A+\omega)^{p+1}\pmod p\)

后面这个东西因为第二条结论可以变成 \(x^2\equiv (A^p+\omega ^p)(A+\omega)\pmod p\)

因为 \((A^2-a)^{\frac{p-1}{2}}\equiv -1\pmod p\) 也就是 \({\omega}^{p-1}\equiv -1\pmod p\)

最后式子化简成 \(x^2\equiv (A-\omega)(A+\omega)\pmod p\)

因为 \(A^2-{\omega}^2=a\) 所以 \(x^2\equiv a\pmod p\)

那么总结一下做法：

特判 \(a=0\) 或者 \(p=2\) 的情况，然后如果 \(a^{\frac{p-1}{2}}\equiv -1\pmod p\) 那么无解。

然后随机出一个 \(A\)，使得 \(A^2-a\) 是 \(p\) 的非二次剩余，那么答案就是 \((A+\omega)^{\frac{p+1}{2}}\)

时间复杂度 \(O(2\log p)\)

代码：

#include 
#include 
#include 
using namespace std;
typedef long long lll; int n,mod,now;
int mo(int x,int y){return x+y>=mod?x+y-mod:x+y;}
struct CP{
	 lll x,y;
	 inline CP operator *(const CP &t)const{
	     return (CP){mo(x*t.x%mod,y*t.y%mod*now%mod),mo(x*t.y%mod,y*t.x%mod)};
	 }
};
int iut(){
	int ans=0; char c=getchar();
	while (!isdigit(c)) c=getchar();
	while (isdigit(c)) ans=ans*10+c-48,c=getchar();
	return ans;
}
void print(int ans){
    if (ans>9) print(ans/10);
    putchar(ans%10+48);
}
CP KSM(CP x,int y){
	CP ans=(CP){1,0};
	for (;y;y>>=1,x=x*x)
	    if (y&1) ans=ans*x;
	return ans;
}
int ksm(int x,int y){
	int ans=1;
	for (;y;y>>=1,x=1ll*x*x%mod)
	    if (y&1) ans=1ll*ans*x%mod;
	return ans;
}
int Cipolla(int n,int mod){
	if (!n) return 0;
	if (mod==2) return 1;
	if (ksm(n,(mod-1)>>1)==mod-1) return -1;
	for (int A=1;;){
		A=rand()%(mod-1)+1;
		now=mo(1ll*A*A%mod,mod-n);
		if (ksm(now,(mod-1)>>1)==mod-1){
			CP Ans=KSM((CP){A,1},(mod+1)>>1);
			return Ans.x;
		}
	}
}
int main(){
	for (int Test=iut();Test;--Test,putchar(10)){
		n=iut(),mod=iut();
		int Ans=Cipolla(n,mod);
		if (!n||mod==2) print(Ans);
		else if (Ans==-1) printf("Hola!");
		else{
			int _Ans=mod-Ans;
			if (Ans>_Ans) swap(Ans,_Ans);
			print(Ans),putchar(32),print(_Ans);
		} 
	}
	return 0;
}

三次剩余（Peralta Method Extension）

\(a=0\) 和 \(p\leq 3\) 的情况特殊判断。

那么质数就可以在模三意义下分为一和二两种。

众所周知，一个三次单位根为 \(\omega=\frac{\sqrt{3}i-1}{2}\)，

那么如果它是一个实数且方程有解，那么就有三个实数解，\(sqrt{p-3}\) 用二次剩余求出来。

也就是说如果找到一个实数解 \(x\)，那么 \(x\omega\) 和 \(x{\omega}^2\) 也是合法的实数解

如果 \(p\bmod 3=2\) 那么有且仅有一解 \(a^{\frac{2p-1}{3}}\)，另外两个解都不是实数。

考虑 \(p\bmod 3=1\) 的情况，那么有解当且仅当 \(a^{\frac{p-1}{3}}\equiv 1\pmod p\)

引理：如果 \(k|p-1\)，那么 \(a\) 是 \(p\) 的 \(k\) 次剩余当且仅当 \(a^{\frac{p-1}{k}}\equiv 1\pmod p\)

证明：充分性显然，考虑必要性的证明，设 \(p\) 的一个原根为 \(g\) 且 \(g^t=a\)，

那么在指数上，\(t\frac{p-1}{k}\equiv 0\pmod{p-1}\) 那也就是要保证 \(k|p-1\) 必要性也可以证明出来。

随机出三个数 \(A,B,C\) 构造多项式 \(f(x)=Ax^2+Bx+C,x^3=a\) 那么 \(f(x)^{p-1}\equiv 1\pmod p\)

那么如果 \(f(x)^{\frac{p-1}{3}}\) 只有一次项 \(B'\) 非0，那么 \((B'x)^3=1\) 即 \(B'=x^{-1}\)

期望次数为9次，所以时间复杂度为 \(O(9\log p)\)

代码：

#include 
#include 
#include 
using namespace std;
typedef long long lll; int A,mod,now,NOW;
int mo(int x,int y){return x+y>=mod?x+y-mod:x+y;}
struct CP{
	 lll x,y;
	 inline CP operator *(const CP &t)const{
	     return (CP){mo(x*t.x%mod,y*t.y%mod*now%mod),mo(x*t.y%mod,y*t.x%mod)};
	 }
};
int iut(){
	int ans=0; char c=getchar();
	while (!isdigit(c)) c=getchar();
	while (isdigit(c)) ans=ans*10+c-48,c=getchar();
	return ans;
}
void print(int ans){
    if (ans>9) print(ans/10);
    putchar(ans%10+48);
}
CP KSM(CP x,int y){
	CP ans=(CP){1,0};
	for (;y;y>>=1,x=x*x)
	    if (y&1) ans=ans*x;
	return ans;
}
int ksm(int x,int y){
	int ans=1;
	for (;y;y>>=1,x=1ll*x*x%mod)
	    if (y&1) ans=1ll*ans*x%mod;
	return ans;
}
int Cipolla(int n,int mod){
	if (!n) return 0;
	if (mod==2) return 1;
	if (ksm(n,(mod-1)>>1)==mod-1) return -1;
	for (int A=1;;){
		A=rand()%(mod-1)+1;
		now=mo(1ll*A*A%mod,mod-n);
		if (ksm(now,(mod-1)>>1)==mod-1){
			CP Ans=KSM((CP){A,1},(mod+1)>>1);
			return Ans.x;
		}
	}
}
struct Three{
	lll f[3];
	inline Three operator *(const Three &t)const{
		Three F=(Three){0,0,0};
	    for (int i=0;i<3;++i)
	    for (int j=0;j<3;++j)
		if (i+j<3) F.f[i+j]+=f[i]*t.f[j]%mod;
		    else F.f[i+j-3]+=f[i]*t.f[j]%mod*NOW%mod;
		for (int i=0;i<3;++i) F.f[i]%=mod;
		return F;
	}
};
Three ksM(Three x,int y){
	Three ans=(Three){1,0,0};
	for (;y;y>>=1,x=x*x)
	    if (y&1) ans=ans*x;
	return ans;
}
Three Peralta(int A,int mod){
	if (!A) return (Three){0,-1,-1};
	if (mod<4) return (Three){A,-1,-1};
	if (mod%3==2) return (Three){ksm(A,(2*mod-1)/3),-1,-1};
	if (ksm(A,(mod-1)/3)!=1) return (Three){-1,-1,-1};
	int omega=1ll*(Cipolla(mod-3,mod)-1)*ksm(2,mod-2)%mod;
	Three X; NOW=A;
	while (1){
		X=(Three){rand()%mod,rand()%mod,rand()%mod};
		Three Ans=ksM(X,(mod-1)/3);
		if (!Ans.f[0]&&!Ans.f[2]&&Ans.f[1]){
			Ans.f[1]=ksm(Ans.f[1],mod-2);
			Ans.f[0]=Ans.f[1]*omega%mod;
			Ans.f[2]=Ans.f[0]*omega%mod;
			return Ans;
		}
	}
}
int main(){
	for (int Test=iut();Test;--Test,putchar(10)){
		mod=iut(),A=iut();
		Three Ans=Peralta(A,mod);
		if (Ans.f[0]==-1) printf("No Solution");
		else if (Ans.f[1]==-1) print(Ans.f[0]);
		else{
			sort(Ans.f,Ans.f+3);
			for (int i=0;i<3;++i)
			    print(Ans.f[i]),putchar(32);
		}
	}
	return 0;
}

后记

二次剩余来自 beginend的二次剩余Cipolla算法学习小记
三次剩余来自 skywalkert的寻找模质数意义下的二次剩余与三次剩余