二次剩余

给出 \(n,p\)，求解方程

\[x^2\equiv n~(\operatorname{mod}~p) \]

其中保证 \(p\) 是奇素数。

---

二次剩余数量

我们假设 \(n\) 在模意义下有多个不同的根，其中两个是 \(x_1,x_2\)，那么 \(x_1^2\equiv x_2^2\)。

移项加拆括号得：\((x_1-x_2)(x_1+x_2)=0\)。因为 \(x_1\neq x_2\)，所以 \(x_1+x_2=0\) 且此方程只有这两个根。

又可以知道，任意一对相反数都对应一个二次剩余，而且这些二次剩余是两两不同的。
也就说二次剩余的数量恰为 \(\frac{p-1}{2}\)，其他的非 \(0\) 数都是非二次剩余，数量也是 \(\frac{p-1}{2}\)。

---

欧拉准则

由费马小定理我们可以知道：\(n^{p-1}\equiv 1~(\operatorname{mod}~p)\)。

由于 \(p\) 是奇素数，所以：\(\left(n^{\frac{p-1}{2}}\right)^2\equiv 1~(\operatorname{mod}~p)\)。

因为 \(1\) 开根后等于 \(1,-1\)：\(n^{\frac{p-1}{2}}\equiv \pm 1~(\operatorname{mod}~p)\)。

我们假设 \(n^{\frac{p-1}{2}}\equiv 1~(\operatorname{mod}~p)\)，令 \(n=g^k\)，\(g\) 为模 \(p\) 的原根。

则有 \(g^{k\frac{p-1}{2}}\equiv 1\)，又 \(g^{\varphi(p)}\equiv 1\)，所以 \(p-1|k\frac{p-1}{2}\)，\(k\) 为偶数。

那么我们可以得出，当 \(n=g^{\frac{k}{2}}\) 时，\(n\) 为二次剩余。

所以 \(n^{\frac{p-1}{2}} \equiv 1\) 与 \(n\) 是二次剩余是等价的，\(n^{\frac{p-1}{2}} \equiv -1\) 与 \(n\) 是非二次剩余是等价的。

---

Cipolla法求解方程

我们先求出一个 \(a\) 使得 \(a^2-n\) 是非二次剩余。因为它的数量是一半一半，所以随机检验期望 \(2\) 次就能找出这个 \(a\)。

定义 \(i^2=a^2-n\)，因为 \(a^2-n\) 是非二次剩余，所以我们求不出来 \(i\) 。

但我们可以类比虚数，定义这一个 \(i\) ，然后所有数变成 \(A+Bi\)。

引理1：\(i^p=-i\)。

证明：\(i^p=i(i^2)^{\frac{p-1}{2}}=i(a^2-n)^{\frac{p-1}{2}}=-i\)

引理2：\((a+b)^p=a^p+b^p\)。

证明：二项式定理展开，消去模 \(p\) 为 \(0\) 的组合数，\((a+b)^p=\sum\limits_{i=0}^p\binom{p}{i}a^ib^{p-i}=a^p+b^p\)

推论：\((a+i)^{p+1}=n\)。

证明：\((a+i)^{p+1}=a^{p+1}+i^{p+1}=a^2-i^2=a^2-(a^2-n)=n\)

因为 \(p+1\) 是偶数，所以 \(x_1=(a+i)^{\frac{p+1}{2}}\) 就是我们求出来的一个解。

至于虚部，可以证明这里虚部为 \(0\)，证明自行bdfs。

因为两个解是相反数，那么另一个解即为 \(x_2=p-x_1\)。

至此，我们可以在期望复杂度 \(O(\log p)\) 内求解二次剩余。

代码

我这里写了两个快速幂，其实可以只用写扩域的快速幂。

#include<bits/stdc++.h>
#define pc(x) putchar(x)
#define int long long
using namespace std;
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while(ch<'0'||ch>'9'){f=ch=='-'?-1:f;ch=getchar();}
    while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    return x*f;
}
void write(int x)
{
    if(x<0){x=-x;putchar('-');}
    if(x>9)write(x/10);
    putchar(x%10+48);
}
int T,n,p,i2,x0,x1;
struct CP
{
    int x,y;
    CP(int xx=0,int yy=0){x=xx;y=yy;}
    friend CP operator *(CP x,CP y)
    {return CP((x.x*y.x+i2*x.y%p*y.y)%p,(x.x*y.y+x.y*y.x)%p);}
};
CP qpow(CP x,int y)
{
    CP res(1,0);
    while(y)
    {
        if(y&1)res=res*x;
        x=x*x;y>>=1;
    }return res;
}
int qpow(int x,int y)
{
    int res=1;
    while(y)
    {
        if(y&1)res=res*x%p;
        x=x*x%p;y>>=1;
    }return res;
}
bool check(int x){return qpow(x,(p-1)>>1)==1;}
signed main()
{
    srand(time(NULL));T=read();
    while(T-->0)
    {
        n=read(),p=read();
        if(!n){puts("0");continue;}
        if(qpow(n,(p-1)/2)==p-1){puts("Hola!");continue;}
        int a=rand()%p;
        while(!a||check((a*a-n+p)%p))a=rand()%p;
        i2=(a*a%p-n+p)%p;x0=qpow(CP(a,1),(p+1)>>1).x;
        x1=p-x0;if(x0>x1)swap(x0,x1);
        write(x0),pc(' ');if(x1!=x0)write(x1);
        pc('\n');
    }return 0;
}