啃了这么久 oi wiki，我决定写一个正常人稍微能理解的东西。

二次剩余

定义：对于两个整数 $a,p$ 满足 $\gcd(a,p)=1$，若存在整数 $x$ 满足 $x^2\equiv a\ (\text{mod}\ p)$，其中 $0<x<p$，那么我们称 $a$ 为模 $p$ 下的二次剩余（后文在模数显然为 $p$ 的情况下也可以简写为‘$a$ 为二次剩余’），否则为二次非剩余。

于是，题目就变成了判断 $y$ 是否为二次剩余并求出一个 $x$。

如何判断二次剩余？

Euler 判别法：

当模数为奇素数时，我们有如下定理：

对于上述的 $a,p$，$a$ 为二次剩余当且仅当 $a^{\frac{p-1}{2}}\equiv 1(\text{mod}\ p)$。

证明如下：

显然，$p$ 为质数且 $\gcd(a,p)=1$，根据费马小定理，有 $a^{p-1}\equiv 1(\text{mod}\ p)$。同时 $x\le p$，所以 $\gcd(x,p)=1$，因此 $x^{p-1}\equiv 1(\text{mod}\ p)$。

我们进一步考虑第一个式子，易得 $a^{p-1}-1\equiv 0(\text{mod}\ p)$，因为 $p$ 为奇数，所以运用平方差公式得到：

$$(a^{\frac{p-1}{2}}+1)(a^{\frac{p-1}{2}}-1)\equiv 0(\text{mod}\ p)$$

因此，$a$ 必定满足：

$$a^{\frac{p-1}{2}}\equiv \pm 1(\text{mod}\ p)$$

这里出现了 $\pm 1$，这让我们不禁想起费马小定理中的那个 $1$。

为了利用这么一个性质，我们看这么一个式子：

$$x^{p-1}-a^{\frac{p-1}{2}}$$

同理，我们用 $p$ 为奇素数的性质把前面改一下：

$$(x^2)^{\frac{p-1}{2}}-a^{\frac{p-1}{2}}$$

我们知道，对于任意的 $x^n-y^n$，在 $n\ge 2$ 时总能分解成 $x-y$ 与一个多项式的乘积，运用一下，得到：

$$(x^2-a)P(x)$$

其中 $P(x)$ 是一个整数系多项式。

看回对 $x$ 使用费马小定理的那个式子，两边同时乘一个 $x$ 可得：

$$x^p-x\equiv 0(\text{mod}\ p)$$

把左边添加几项，再搬到右边，变成：

$$x^p-x\equiv x(x^{p-1}-a^{\frac{p-1}{2}})+x(a^{\frac{p-1}{2}}-1)(\text{mod}\ p)$$

带入上面的式子，得到：

$$x^p-x\equiv x(x^2-a)P(x)+x(a^{\frac{p-1}{2}}-1)(\text{mod}\ p)$$

这就是证明要用的式子。

首先证明充分性：令 $a^{\frac{p-1}{2}}\equiv1(\text{mod}\ p)$，则 $x(x^2-a)P(x)\equiv x^p-x(\text{mod}\ p)$。

我们将右边记为 $F(x)$，显然对于任意的 $0<x<p$，$F(x)$ 在模 $p$ 意义下均为 $0$，所以所有的形如 $x-a$ 的式子都是 $F(x)$ 的因式。

带进去得到：

$$x(x^2-a)P(x)\equiv (x-0)(x-1)(x-2)\dots(x-(p-1))(\text{mod}\ p)$$

此时就能得到 $(x^2-a)|(x-0)(x-1)(x-2)\dots(x-(p-1))$，也就是：

$$\frac{(x-0)(x-1)(x-2)\dots(x-(p-1))}{(x^2-a)}\in \mathbb{Z}$$

但下面的式子中有一个二次项，如果这个式子没办法继续分解，整个分数就不可能是整数，所以下面的式子一定是可以因式分解的。

我们假定 $(x^2-a)\equiv(x-r_1)(x-r_2)(\text{mod}\ p)$，此时直接取 $x=r_1$，显然右边变成了 $0$，所以：

$$r^2_1-a\equiv 0(\text{mod}\ p)$$

综上，$r^2_1\equiv a$，我们找到了一个解。

接下来证明必要性。

我们设 $r^2=a$，显然 $a^{\frac{p-1}{2}}=r^{p-1}$，又 $r^{p-1}\equiv 1(\text{mod}\ p)$，所以证完了。

因此，在题目中，快速幂判断一下就可以了，用不上 Legendre 符号。

求解二次剩余

此处讲解 Cipolla 算法，该算法的核心就是找到一个二次非剩余，然后虚构一个它的解（就像虚数单位 $i$）。进而求出答案。

首先，我们通过随机的方法找到一个 $r$ 使得 $r^2-a$ 为二次非剩余（对于 $a=0$ 不存在，所以需要特判一下），期望两次就能找到，判断就可以使用快速幂就可以。

为什么期望两次就够？

我们考虑什么时候 $r^2-a$ 是二次剩余，此时存在一个 $x$ 使得 $$r^2-a\equiv x^2(\text{mod}\ p)$$移项之后得到：$$r^2-x^2\equiv a(\text{mod}\ p)$$ 平方差公式就有：$$(r+x)(r-x)\equiv a(\text{mod}\ p)$$ 在模 $p$ 意义下，$(r+x)$ 只有 $[1,p-1]$ 共 $p-1$ 种取值。因为 $p$ 是质数，所以对于每一种取值，总有对应的 $(r-x)$，也就有对应的 $r$ 和 $x$。所以原方程共有 $p-1$ 组解。

分类讨论一下就能发现，$x\equiv 0(\text{mod}\ p)$ 时 $r$ 和 $-r$ 都为原方程的解，而其他的 $p-1-2=p-3$ 种情况中，每一个 $r+x$ 都有对应的 $r+(-x+p)$，所以一个 $r$ 就对应两个解，因此合法的 $r$ 的取值共有 $2+\frac{p-3}{2}=\frac{p+1}{2}$ 个。

那么每一次得到二次剩余的概率就是：$$\frac{\frac{p+1}{2}}{p}=\frac{p+1}{2p}$$ 二次非剩余的就是：$$1-\frac{p+1}{2p}=\frac{p-1}{2p}$$ 期望次数就是：$$\frac{2p}{p-1}\approx2$$

得到 $r$ 之后，我们定义一个 $w$ 使得 $w^2\equiv r^2-a(\text{mod}\ p)$，这样，我们就创造出了模 $p$ 意义下的”复数“，这里面的所有数都可以被写成 $x+yw$ 的形式，其中 $x,y\in\mathbb{Z}_p$。

接下来我们定义 $\beta=r+w$，在模 $p$ 意义下考虑这么一个式子： 于是：

$$(\beta^{\frac{p+1}{2}})^2=a$$

令 $\gamma=\beta^\frac{p+1}{2}=g_0+g_1w(g_0,g_1\in\mathbb{Z}_p)$，则：

$$g_0^2+2g_0g_1w+g_1^2\times(r^2-a)\equiv a(\text{mod}\ p)$$

但 $a$ 是二次剩余，不可能含有 $w$，所以 $2g_0g_1\equiv0(\text{mod}\ p)$，所以 $g_0\equiv0(\text{mod}\ p)$ 或 $g_1\equiv0(\text{mod}\ p)$。

但 $g_0\equiv 0$ 时，上式变为：

$$g_1^2\times(r^2-a)\equiv a(\text{mod}\ p)$$

因为 $r^2-a$ 为二次非剩余，所以 $g_1^2\times(r^2-a)$ 不可能是二次剩余，这与 $a$ 是二次剩余相矛盾。

所以只有 $g_1=0$，那么 $g_0$，也就是 $(r+w)^{\frac{p+1}{2}}$ 的”实部“就是我们要求的答案。

代码：

#include<bits/stdc++.h>
using namespace std;
long long d,p;//d即为r^2-a，我们定义的复数单位w满足w^2=d,p为模数 
struct C{//定义模意义下的复数 
	long long a,b;//a为实部，b为虚部
	C operator *(const C &ano)const{
		return {(a*ano.a+d*(b*ano.b%p)%p)%p,(a*ano.b+b*ano.a)%p};//(a+bw)(x+yw)=ax+ayw+bxw+byw^2=(ax+byd)+(ay+bx)w
	}
};
mt19937 rd(99573);
C power(C a,long long b){//复数的快速幂，算结果用 
	C ans;
	ans.a=1;
	ans.b=0;
	while(b){
		if(b&1){
			ans=ans*a;
		}
		a=a*a;
		b>>=1;
	}
	return ans;
}
long long power(long long a,long long b){//快速幂，快速判断二次剩余 
	long long ans=1;
	while(b){
		if(b&1){
			ans=ans*a%p;
		}
		a=a*a%p;
		b>>=1;
	}
	return ans;
}
long long cip(long long a){
	if(a==0)return 0;
	if(p==2)return 1;
	if(power(a,(p-1)/2)==p-1)return -1;
	
	//构建w^2
	long long r;
	for(r=rd()%p;;r=rd()%p)if(power((r*r-a+p)%p,(p-1)/2)==p-1)break;
	d=(r*r-a+p)%p; 
	
	//快速幂求答案 
	return power({r,1},(p+1)/2).a;
}
int main(){
	int t;
	cin>>t;
	while(t--){
		long long y;
		cin>>y>>p;
		cout<<cip(y)<<"\n";
	}
	return 0;
}

Jacobi 符号

在动手求平方根之前，必须先确认方程 $x^2 \equiv Y \pmod P$ 到底有没有解。

根据欧拉判别法，$Y$ 有平方根当且仅当 $Y^{\frac{P-1}{2}} \equiv 1 \pmod P$。虽然这没错，但直接算这个幂次太慢了。我们需要一个像“辗转相除法”一样快的判定工具，这就是 Jacobi 符号 $\left(\frac{a}{n}\right)$。

它不需要分解质因数，只靠以下三条规则就能在 $O(\log n)$ 时间内算出结果：

1. 取模性质：$\left(\frac{a}{n}\right) = \left(\frac{a \bmod n}{n}\right)$

2. 提取因子2：$\left(\frac{2}{n}\right) = (-1)^{\frac{n^2-1}{8}}$

3. 二次互反律：若 $a, n$ 都是正奇数且互质，则$\left(\frac{a}{n}\right)\left(\frac{n}{a}\right) = (-1)^{\frac{a-1}{2}\cdot\frac{n-1}{2}}$

当 $n$ 是合数时，$\left(\frac{a}{n}\right)=1$ 不代表 $a$ 一定是二次剩余；但 $\left(\frac{a}{n}\right)=-1$ 一定代表 $a$ 不是二次剩余。在 Cipolla 算法中，主判断时 $P$ 是素数，完全等价于 Legendre 符号；找辅助参数时只用它做排除，逻辑也是安全的。

讲解 Cipolla 算法

直接求 $\sqrt{Y} \pmod P$ 很难，Cipolla 的天才想法是：既然在整数里找不到 $\sqrt{d}$，我们就人为造一个出来！ 这和虚数单位 $i=\sqrt{-1}$ 的思路一模一样。

随机选一个整数 $b$，计算 $d = b^2 - Y \pmod P$。反复尝试直到 $\left(\frac{d}{P}\right) = -1$（即 $d$ 不是二次剩余）。（期望次数为 $2$ ,其实就是一半的数字都符合条件） 此时 $\sqrt{d}$ 在模 $P$ 整数中不存在，我们把它当作新的“虚数单位” $\omega$。 所有数都写成 $x + y\omega$ 的形式，其中 $x,y \in \mathbb{Z}_P$。

记住 $\omega^2 = d$，展开后合并同类项即可： $$(x_1 + y_1\omega)(x_2 + y_2\omega) = (x_1x_2 + y_1y_2 d) + (x_1y_2 + x_2y_1)\omega$$

在普通整数模 $P$ 中，我们有 $x^P \equiv x$。在这个新数系中，有一个极其重要的性质：

$$(x + y\omega)^P = x - y\omega$$

证明：

1. 二项式展开：$(x+y\omega)^P = x^P + \binom{P}{1}x^{P-1}(y\omega) + \dots + (y\omega)^P$。因为 $P$ 是素数，中间所有组合数 $\binom{P}{k}$ 都能被 $P$ 整除，模 $P$ 后全部消失！只剩首尾两项：$x^P + (y\omega)^P$。
2. 费马小定理+欧拉判别：$x^P \equiv x$，$y^P \equiv y$；而 $(\omega)^P = \omega \cdot (\omega^2)^{\frac{P-1}{2}} = \omega \cdot d^{\frac{P-1}{2}}$。因为我们特意选了 $d$ 为非二次剩余，所以 $d^{\frac{P-1}{2}} \equiv -1$，故 $\omega^P = -\omega$。
3. 合并：$(x+y\omega)^P = x + y(-\omega) = x - y\omega$。证毕！

性质说明在新数系里，$P$ 次幂就等于“共轭”（把 $\omega$ 变成 $-\omega$），就像复数中 $\bar{z}$ 的作用一样。

考虑元素 $\beta = b + \omega$，利用上面的引理计算 $\beta^{P+1}$：

$$\beta^{P+1} = \beta^P \cdot \beta = (b - \omega)(b + \omega) = b^2 - \omega^2$$

代入 $\omega^2 = d = b^2 - Y$，得：

$$\beta^{P+1} = b^2 - (b^2 - Y) = Y$$

既然 $\beta^{P+1} = Y$，那么两边开方（指数除以2）：

$$\left( \beta^{\frac{P+1}{2}} \right)^2 = Y$$

令 $\gamma = \beta^{\frac{P+1}{2}} = g_0 + g_1\omega$，则 $\gamma^2 = Y$。展开左边：

$$g_0^2 + g_1^2 d + 2g_0g_1\omega = Y$$

注意右边 $Y$ 是纯整数（没有 $\omega$ 项），所以左边的“虚部”必须为 0：

$$2g_0g_1 \equiv 0 \pmod P$$

因为 $P$ 是奇素数，$2 \not\equiv 0$，所以要么 $g_0=0$，要么 $g_1=0$。

• 若 $g_0=0$，则 $g_1^2 d = Y$，意味着 $\left(\frac{Y}{P}\right) = \left(\frac{d}{P}\right) = -1$，但这与我们一开始判断 $Y$ 是二次剩余矛盾！
• 因此必然 $g_1=0$，即 $\gamma = g_0$ 是一个纯整数，且 $g_0^2 \equiv Y \pmod P$。

结论：$\beta^{\frac{P+1}{2}}$ 的实部 $g_0$ 就是我们要找的 $\sqrt{Y}$！

时间复杂度

步骤	复杂度
Jacobi 符号	$O(\log P)$
找非二次剩余	期望 $O(1)$
扩域快速幂	$O(\log P)$
总计	$O(\log P)$

代码

呼~，终于到代码了。

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

// 扩域元素: x + y*sqrt(d)
struct node{int x,y;};

// 扩展欧几里得求逆元 (迭代版，避免递归开销)
inline int modInv(int a,int m){
	int b=m,x=1,y=0,t;
	while(true){
		t=a/b;
		a-=t*b;
		if(!a){
			if(b==-1)y=-y;
			return(y<0)?y+m:y;
		}
		x-=t*y;
		t=b/a;
		b-=t*a;
		if(!b){
			if(a==-1)x=-x;
			return(x<0)?x+m:x;
		}
		y-=t*x;
	}
}

// Jacobi 符号判定 (O(log n))
// 返回 1: 可能是二次剩余, -1: 非二次剩余, 0: 整除
inline int jacobi(int a,int m){
	int s=1;
	if(a<0)a=a%m+m;
	while(m>1){
		a%=m;
		if(!a)return 0;
		// 提取因子 2: (2/m) = (-1)^((m^2-1)/8)
		int r=__builtin_ctz(a);
		if((r&1)&&((m+2)&4))s=-s;
		a>>=r;
		// 二次互反律调整符号
		if(a&m&2)s=-s;
		swap(a,m);
	}
	return s;
}

// Cipolla 算法求模平方根
inline int modSqrt(int a,int p){
	if(p==2)return a&1; // 处理素数 2 的边界情况
	
	int j=jacobi(a,p);
	if(!j)return 0;     // a 是 p 的倍数
	if(j==-1)return -1; // a 不是二次剩余，无解

	int b,d;
	// 随机寻找 b 使得 d = b^2 - a 为非二次剩余 (期望尝试 2 次)
	while(true){
		b=rand()%p;
		d=(b*b-a)%p;
		if(d<0)d+=p;
		if(jacobi(d,p)==-1)break;
	}

	// 在扩域 Z_p[sqrt(d)] 中计算 (b + sqrt(d))^((p+1)/2)
	// f 为底数 (b + 1*sqrt(d)), g 为结果累加器 (初始为 1 + 0*sqrt(d))
	node f={b,1},g={1,0};
	int tmp;
	
	// 快速幂
	for(int e=(p+1>>1);e;e>>=1){
		if(e&1){
			// g = g * f
			tmp=(g.x*f.x%p+d*(g.y*f.y%p)%p)%p; // 实部: x1x2 + y1y2*d
			g.y=(g.x*f.y%p+g.y*f.x%p)%p;       // 虚部: x1y2 + x2y1
			g.x=tmp;
		}
		// f = f * f
		tmp=(f.x*f.x%p+d*(f.y*f.y%p)%p)%p;
		f.y=(2*f.x*f.y)%p;
		f.x=tmp;
	}
	
	// 根据推导，结果的虚部 g.y 必为 0，实部 g.x 即为 sqrt(a)
	return min(g.x,p-g.x);
}

signed main(){
	srand(time(0));
	ios::sync_with_stdio(false);
	cin.tie(0),cout.tie(0);
	int T;
	cin>>T;
	while(T--){
		int Y,P;
		cin>>Y>>P;
		cout<<modSqrt(Y,P)<<"\n";
	}
}