下文只在模质数意义下讨论。
即给定
n
>
0
,
p
,
p
n>0,p,p
n>0,p,p是质数,求所有
x
∈
[
0
,
p
)
,
x
2
=
n
(
m
o
d
p
)
x\in[0,p),x^2=n\pmod p
x∈[0,p),x2=n(modp)。
若存在解则
n
n
n称为
p
p
p的二次剩余。
首先将
p
=
2
p=2
p=2的特殊情况判掉。下文
p
p
p是奇质数。
首先显然模奇质数
p
p
p意义下会有恰好
p
−
1
2
\frac{p-1}2
2p−1个二次剩余。
首先如何判定,有一个定理:
(
−
1
)
[
n
不
是
p
的
二
次
剩
余
]
=
n
p
−
1
2
(
m
o
d
p
)
(-1)^{[n不是p的二次剩余]}=n^{\frac{p-1}2}\pmod p
(−1)[n不是p的二次剩余]=n2p−1(modp)
并且若有解
x
x
x,则仅有
x
x
x和
p
−
x
p-x
p−x这两组解。这个定理易证,略去证明。
(上述两个定理都可以通过百度“二次剩余”搜到证明(应该能))
那么考虑一个算法:在
[
1
,
p
)
[1,p)
[1,p)随机一个整数
t
t
t,使得:
t
2
−
n
t^2-n
t2−n不是关于
p
p
p的二次剩余。由于有恰好一半不是,因此随机两步就会随机出一组。那么令
w
=
t
2
−
n
w=\sqrt{t^2-n}
w=t2−n(这里可能不存在在模
p
p
p意义下后半部分的平方根,这里只是定义一个有
w
2
=
t
2
−
n
(
m
o
d
p
)
w^2=t^2-n\pmod p
w2=t2−n(modp)性质的
w
w
w)。
可以显然证明:
w
P
=
−
w
w^P=-w
wP=−w。
那么考虑令
x
=
(
t
+
w
)
p
+
1
2
x=(t+w)^{\frac{p+1}2}
x=(t+w)2p+1
x
2
=
(
t
+
w
)
p
+
1
=
(
t
+
w
)
(
t
+
w
)
p
=
(
t
+
w
)
∑
i
=
0
p
(
p
i
)
t
i
w
p
−
i
=
(
t
+
w
)
(
t
p
+
w
p
)
=
(
t
+
w
)
(
t
−
w
)
=
t
2
−
w
2
=
t
2
−
t
2
+
n
=
n
x^2=(t+w)^{p+1}\\=(t+w)(t+w)^p\\=(t+w)\sum_{i=0}^p\binom pit^iw^{p-i}\\=(t+w)(t^p+w^p)\\=(t+w)(t-w)=t^2-w^2\\=t^2-t^2+n=n
x2=(t+w)p+1=(t+w)(t+w)p=(t+w)∑i=0p(ip)tiwp−i=(t+w)(tp+wp)=(t+w)(t−w)=t2−w2=t2−t2+n=n
因此
x
x
x就是
x
2
=
n
(
m
o
d
p
)
x^2=n\pmod p
x2=n(modp)的一组解,其解集为
{
x
,
p
−
x
}
\{x,p-x\}
{x,p−x}
时间复杂度
O
(
lg
p
)
O(\lg p)
O(lgp)。实现的时候可以将
(
t
+
w
)
k
(t+w)^k
(t+w)k视为
v
+
c
w
v+cw
v+cw,利用
w
2
=
t
2
−
n
w^2=t^2-n
w2=t2−n即可做一次乘法。代码(应该是正确的):
#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define gc getchar()
using namespace std;typedef long long lint;
inline int inn() { int x,ch;while((ch=gc)<'0'||ch>'9');x=ch^'0';while((ch=gc)>='0'&&ch<='9') x=(x<<1)+(x<<3)+(ch^'0');return x; }
inline int fast_pow(int x,int k,int p,int ans=1) { for(;k;k>>=1,x=(lint)x*x%p) (k&1)?ans=(lint)ans*x%p:0;return ans; }
struct RAND{ unsigned int x;RAND() { x=1000000007; }
inline int operator() (int a,int b) { return x=(x<<11)^(x>>3)^x^998244353,int(x%(b-a+1)+a); }
}rnd; struct E{ int t,c;E(int _t=0,int _c=0) { t=_t,c=_c; } };
inline E tms(const E &a,const E &b,int w2,int p) { return E(((lint)a.t*b.t+(lint)a.c*b.c%p*w2)%p,((lint)a.c*b.t+(lint)a.t*b.c)%p); }
inline int calc(E x,int w2,int k,int p) { E ans(1,0);for(;k;k>>=1,x=tms(x,x,w2,p)) if(k&1) ans=tms(ans,x,w2,p);return ans.t; }
inline int solve(int n,int p)
{
if(fast_pow(n,(p-1)/2,p)==p-1) return -1;int t,w;
do{
t=rnd(1,p-1),w=((lint)t*t-n)%p;if(w<0) w+=p;
}while(fast_pow(w,(p-1)/2,p)==1);
return calc(E(t,1),w,(p+1)/2,p);
}
int main()
{
for(int T=inn();T;T--)
{
int n=inn(),p=inn();n%=p;if(!n) { printf("0\n");continue; }
if(p==2) { printf("1\n");continue; }int ans1=solve(n,p),ans2=p-ans1;
if(ans1==-1) { printf("No answer!\n");continue; }
if(ans1>ans2) swap(ans1,ans2);printf("%d %d\n",ans1,ans2);
}
return 0;
}
`