二次剩余, 俗称模意义开根。
实际上是求 (x^2equiv n pmod p) 的解。
这有两种解法，一个是利用原根来求，另外一种就是 (Cipolla) 算法（只能适用于奇素数的情况）。
我们下面主要介绍第二种。
1.勒让德符号
对于一个正整数 (n), 他的勒让德符号 (left(dfrac{n} {p} ight))定义为：
根据这个我们可以初步判断一个数是否存在二次剩余。
2.解的数量
对于二次剩余 (x^2equiv npmod p) 有多少个解？
假设有多组解，对于任意不相同的两个解有 (x_0^2 = x_1^2)
移项可得： ((x_0-x_1)(x_0+x_1) equiv 0 pmod p)
因为 (p) 是奇素数，且 (x_0 eq x_1), 所以 (x_0-x_1 eq 0)
那么 (x_0+x_1equiv 0), 也就是说 (x_0) 和 (x_1) 在模 (p) 的意义下互为相反数。
换而言之就是这个方程只有两个解，且这两个解在模 (p) 的意义下互为相反数。
另外，每一对相反数都可以对应一个二次剩余，而且这些二次剩余是两两不同的。
也就是说 二次剩余的数量恰为 (p over 2) ，非二次剩余的数量也为 ({pover 2})
3.欧拉判别
(left(dfrac{n}{p} ight) = n^{{p-1}over 2} mod p)
4.Cipolla
求解二次剩余方程 (x^2equiv n pmod p)
我们先随机出一个 (a) 使得 (a^2-n) 不存在二次剩余， 由于非二次剩余的数量接近于 (pover 2)。
所以期望随机两次就可以随机到一个合适的 (a) 。
然后设 (w = sqrt{a^2-n}) , 但 (a^2-n) 不存在二次剩余，怎么把 (w) 求出来呢？
我们实际上不需要把这个 (w) 求出来。 类似于复数，定义一个 (w) 使得 (w^2 = a^2-n) 。
这样每个数都可以表示成： (a+bw) 的形式， (a) 和 (b) 都是模 (p) 意义下的数。
有一个结论： ((a+w) ^{p+1} = n)
引理1： ((A+B)^p = A^p + B^p)
证明：由二项式定理展开可得： ((A+B)^p = displaystylesum_{i=1}^{p} C_{p}^{i} x^iy^{p-i})
在模 (p) 意义下，组合数 (C_{p}^{i}) 只有 (C_{p}^{0},C_{p}^{p}) 为 (1) ,其他的由于分子上的阶乘被消掉了值为 (0).
所以 ((A+B)^p = A^p + B^p)
引理2： (w^{p} = -w)
证明： 因为 (w^{p} = w(w^2)^{p-1over 2} = w(a^2-n)^{p-1over 2} = -w)
我们来通过两个引理证一下上面那个结论
((a+w)^{p+1} = (a+w) (a^{p} + w^{p}) = (a+w)(a-w) = a^2 - w^2)
又因为： (w^2 = a^2 - n) , 所以： (n = a^2-w^2)
得证： ((a+w)^{p+1} = n)
可以得到 (x = (a+w)^{p+1over 2})， 又因为 (p) 为奇素数，所以 (p+1over2) 肯定为整数，那么 (x) 也一定可以算出来。
对于 ((a+bw)(c+dw) = (ac + dbw^2) + (ad+bc)w) ，重载一下乘法，在进行快速幂即可。
还有一个问题，(x) 虚部一定为 (0) 吗
我们来简单证明一下：
把 (x) 写成复数形式可得： (x^2equiv n Rightarrow (a+bw)^2 equiv n)
展开可得： (a^2 + 2abw + bw^2 equiv n)
移项可得： (a^2+bw^2-n equiv -2abw Rightarrow a^2 + b(a^2-n) - n equiv -2abw)
因为左边虚部为 (0) ,所以右边的虚部肯定也为 (0), 即 (ab = 0) 。
假设说 (b eq 0), 那么只能 (a = 0) 。
所以 ((bw)^2 equiv n Rightarrow w^2 = {nover b^2}),
由于 (nover b^2) 在模 (p) 的意义下是确定的，所以 (w^2) 是二次剩余，这与前面所说的 (w^2) 不是二次剩余相矛盾。
所以 (b = 0), 即 (x) 的虚部肯定为 (0).
最后方程的两个解分别是 (x) 和 (p-x)
code
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<bits/stdc++.h>
using namespace std;
#define int long long
int T,n,p,a,w;
struct Complex
{
int x,y;
Complex(){}
Complex(int a,int b){x = a, y = b;}
};
inline int read()
{
int s = 0,w = 1; char ch = getchar();
while(ch < '0' || ch > '9'){if(ch == '-') w = -1; ch = getchar();}
while(ch >= '0' && ch <= '9'){s = s * 10 + ch - '0'; ch = getchar();}
return s * w;
}
Complex operator * (Complex a,Complex b)
{
Complex c;
c.x = (a.x * b.x % p + a.y * b.y % p * w % p) % p;
c.y = (a.x * b.y % p + a.y * b.x % p) % p;
return c;
}
Complex Pow(Complex a,int b)
{
Complex c = Complex(1,0);
for(; b; b >>= 1)
{
if(b & 1) c = c * a;
a = a * a;
}
return c;
}
int ksm(int a,int b)
{
int res = 1;
for(; b; b >>= 1)
{
if(b & 1) res = res * a % p;
a = a * a % p;
}
return res;
}
bool pd(int n)
{
return ksm(n,(p-1)/2) == 1;
}
void slove(int n,int p)
{
n %= p;
if(!pd(n))
{
printf("Hola!
");
return;
}
srand(time(0)); rand();
while(1)
{
a = rand() % p;
w = (a * a % p - n + p) % p;
if(ksm(w,(p-1)/2) == p-1) break;
}
Complex Ans = Pow(Complex(a,1),(p+1)/2);
int ans1 = Ans.x, ans2 = p - Ans.x;
if(ans1 > ans2) swap(ans1,ans2);
printf("%lld ",ans1);
if(ans1 == ans2) printf("
");
else printf("%lld
",ans2);
}
signed main()
{
T = read();
while(T--)
{
n = read(); p = read();
if(n == 0) printf("%d
",0);
else slove(n,p);
}
return 0;
}