理论部分

求解类似x^2 \equiv n (\mod p)这样的方程称为二次剩余问题.如果存在x满足以上方程,则称n是在模p意义下的二次剩余.

符号定义

定义勒让德符号如下:

(\frac{n}{p})= \begin{cases} 1& \text {n是模p意义下的二次剩余}\\ -1& \text{n不是模p意义下的二次剩余}\\ 0& \text{p|n} \end{cases}

同时,有如下公式存在:

(\frac{n}{p}) \equiv n^{\frac{p-1}{2}} \ ( \mod p)

根据费马小定理可以轻松证明上述公式.有了这个公式后,就可以在logp的时间内快速判断一个数是不是二次剩余.同时还有一个结论是,模p的剩余系中恰好有一半是二次剩余,一半不是.这点为之后算法的时间复杂度构成了保障.

算法原理

Cipolla求解二次同余方程的流程如下:

1. 判断n是不是二次剩余,如果不是,返回无解
2. 随机出一个数a,使得a^2-n是二次非剩余.根据上面的性质可知,期望随机两次.
3. 令\omega^{2}\equiv a^2-n (\mod p),则y\equiv (a+\omega)^{\frac{p+1}{2}}即为原方程的解

上述流程看起来非常的神奇:一个非二次剩余,强行求解出了对应方程的根,并据此得到了原方程的解.时间复杂度是O(logp)

在证明上面那个算法的正确性之前,首先要对数域进行扩充.类似于实数域无法对负数进行开根,从而补充虚根i得到复数域的情况,把非二次剩余a^2-n在模p意义下的平方根记为\omega,这样,在新的数域下,其中的每个元素都可以表示成a+b\omega的形式,从而完成了对数域的扩充.在模p意义下的数域记为F_{p}的话,新的数域就记为F_{p^2}.

下面这张图证明了这个新数域是对四则运算封闭的:

接下来需要证明构造得到的y一定是原方程的解.首先要证明如下的几个引理:

1. 引理一:\omega^{p}\equiv-\omega(\mod p)

证明:

\omega^{p}\equiv (\omega^{2})^{\frac{p-1}{2}}\omega\equiv-\omega

2. 引理二:(a+b)^{n}\equiv a^n+b^n(\mod n) \ n \in P

证明:

因为n是质数,故组合数(\frac{n}{i})只有在i!=n \ and \ i!= 0的时候才不是n的倍数,故只会保留下首尾两项

接下来根据上面两个引理证明y就是原方程的解(当然,p-y显然也是).

x^{2}\equiv y^2\equiv (a+\omega)^{p+1}

\equiv (a+\omega)^{p}(a+\omega) \equiv(a^p+\omega^p)(a+\omega)

\equiv(a-\omega)(a+\omega)\equiv(a^2-\omega^2)

\equiv(a^2-(a^2-n))\equiv n (\mod p)

剩下一点小问题是计算出来的y是新的数域上的数,也就是说存在"虚部"不为0的情况,此时无法对应到原数域的一个根.根据拉格朗日定理,在模p意义下的多项式方程f(x)\equiv 0(\mod p)至多有deg(f(x))个根.新的数域是对原数域的扩充,故原数域上的根仍是新数域上的根.当原数域上方程x^2-n\equiv 0(\mod p)一定有解时,必有两个根.而之前在新数域下我们也得到了两个根,显然他们同时就是原数域的根.即"虚部"为0.

代码模板

inline LL Mi(LL x,LL y,LL MOD){
    LL p=x,t=1,Res=1;
    for (;t<=y;(t&y)?Res=(Res*p)%MOD:0,p=(p*p)%MOD,t<<=1);
    return Res;
}

inline F Mul(F a,F b,LL MOD){
    F c=F(0,0,a.w);
    c.a=((a.a*b.a)%MOD+(((a.b*b.b)%MOD)*a.w)%MOD)%MOD;
    c.b=(b.a*a.b%MOD+a.a*b.b%MOD)%MOD;
    return c;
}

inline F Pow(F x,LL y,LL MOD){
    LL t=1;
    F p=x,Res=F(1,0,x.w);
    for (;t<=y;Res=(t&y)?Mul(Res,p,MOD):Res,p=Mul(p,p,MOD),t<<=1);
    return Res;
}

inline LL Solve(int x,int p){
    LL a=0,w=0;
    if (!x) return 0;
    if (Mi(x,(p-1)>>1,p)==p-1)  return -1;
    while (1){
        a=rand();   w=(a*a%p-x+p)%p;
        if (Mi(w,(p-1)>>1,p)==p-1){
            F s=Pow(F(a,1,w),(p+1)>>1,p);
            return s.a;
        }
    }
}

相关例题

[模板题]二次剩余

洛谷链接

代码:

#include<cstring>
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<vector>
#include<cmath>
#include<ctime>
#include<cstdlib>
#define l(x) (x<<1)
#define r(x) ((x<<1)+1)
#define LL long long
#define INF 0x3f3f3f3f
using namespace std;

int T;
LL n,p,Ans,cnt;

struct F{
    LL a,b,w;
    F(){};
    F(const LL &x,const LL &y,const LL &c){a=x; b=y; w=c;}
};

inline int Abs(int x){return (x<0)?-x:x;}
inline void Swap(LL &a,LL &b){a^=b^=a^=b;}
inline int Min(int a,int b){return (a<b)?a:b;}
inline int Max(int a,int b){return (a>b)?a:b;}

inline int read(){
    int p=0;    char    c=getchar();
    while (c<48||c>57)  c=getchar();
    while (c>=48&&c<=57)    p=(p<<1)+(p<<3)+c-48,c=getchar();
    return p;
}

inline LL Mi(LL x,LL y,LL MOD){
    LL p=x,t=1,Res=1;
    for (;t<=y;(t&y)?Res=(Res*p)%MOD:0,p=(p*p)%MOD,t<<=1);
    return Res;
}

inline F Mul(F a,F b,LL MOD){
    F c=F(0,0,a.w);
    c.a=((a.a*b.a)%MOD+(((a.b*b.b)%MOD)*a.w)%MOD)%MOD;
    c.b=(b.a*a.b%MOD+a.a*b.b%MOD)%MOD;
    return c;
}

inline F Pow(F x,LL y,LL MOD){
    LL t=1;
    F p=x,Res=F(1,0,x.w);
    for (;t<=y;Res=(t&y)?Mul(Res,p,MOD):Res,p=Mul(p,p,MOD),t<<=1);
    return Res;
}

inline LL Solve(int x,int p){
    LL a=0,w=0;
    if (!x) return 0;
    if (Mi(x,(p-1)>>1,p)==p-1)  return -1;
    while (1){
        a=rand();   w=(a*a%p-x+p)%p;
        if (Mi(w,(p-1)>>1,p)==p-1){
            F s=Pow(F(a,1,w),(p+1)>>1,p);
            return s.a;
        }
    }
}

int main(){
    #ifdef __Marvolo
    freopen("zht.in","r",stdin);
    freopen("zht.out","w",stdout);
    #endif
    srand(time(NULL));
    scanf("%d",&T);
    while (T--){
        scanf("%lld%lld",&n,&p);
        n%=p;
        Ans=Solve(n,p);
        if (Ans==-1)    printf("Hola!\n");
        else {
            cnt=(p-Ans)%p;
            if (cnt==Ans)   printf("%lld\n",Ans);
            else {
                if (cnt>Ans)    Swap(Ans,cnt);
                printf("%lld %lld\n",cnt,Ans);
            }
        }
    }
    return 0;
}

算法推广

注:以下内容都不包括解的存在性的证明.以后有机会可能会补上.

p为质数的次幂的形式

求解方程x^2\equiv a (\mod p^k)

先求解r^2 \equiv a (\mod p)

故有r^2-a\equiv 0(\mod p),则(r^2-a)^k \equiv 0 (\mod p^k)

令r^2-a\equiv t^2-u^2a (\mod p)

然后有(r-\sqrt{a})^k\equiv t-u\sqrt{a},(r+\sqrt{a})^k\equiv t+u\sqrt{a}

这时采取类似上面的扩充数域的思想进行计算.

最终,得到t^2\equiv u^2a(\mod p)

所以,原方程的解为(tu^{-1})^2\equiv a (mod p)

这里求逆元时采用Exgcd的方式.

p为2的次幂的形式

$k\leq2$

特判

$k=3$

同余方程有解当且仅当a\equiv1(\mod 8),此时,解只可能为\pm1,\pm5.

将解记为x\equiv\pm(x_3+t_3 *2^2),x_3=1,5,t_3\in \Z

$k>3$

有递推公式:x_k=\pm(x_{k-1}+\frac{a_k-x_{k-1}^2}{2}+t_{k-1}*2^{k-1}),t\in \Z

其中,a_k=a\%2^k

根据公式从3开始递推计算即可

p为合数的形式

把p进行分解,得到一系列同余方程.按照p为质数次幂的情况求解,之后再根据中国剩余定理进行合并.

参考文章

https://www.cnblogs.com/bztMinamoto/p/10664973.html

https://blog.csdn.net/sizeof_you/article/details/84841962

https://blog.csdn.net/zxyoi_dreamer/article/details/85195819