同余方程-5天从入门到入土

· · 算法·理论

upd: 应评论区，我自己查了一些资料之后发现原来写的东西确实有些不对的地方，造成的误解请各位海涵，本人蒟蒻，如果对于某个描述还是有疑问的话欢迎继续指出，谢谢各位 ^_^

特别鸣谢 \color{red}\mathsf{NaVi\_Awson} 学长 以及 \color{orange}\mathsf{\text{落汐}} 大佬/qq

目录

（为什么没了呢？因为我不会了TAT）

\color{brown}\texttt{同步题单戳这里}

概念

1.带余除法

实际上就是小学学的那种“3 除以 2 商 1 余 1”的除法，书面一点的描述就是 a=q\times b+r 其中 a,b,q,r 都是整数，a 叫做被除数，b 叫做除数 (b\ne 0)，q 是商，r 是余数（0\le r <b），C++中，整型之间的除法是整数除法，设 a=q\times b+r ，那么 a/b=q

本文中的除法都是整数除法

2.同余方程

同余方程就是长成这样的柿子： ax\equiv c\pmod b
这个式子什么意思？ 就是字面意思，同余即指余数相同（ax\bmod b=c\bmod b）

算法

Day 1.\mathsf{Exgcd} (扩展欧几里得算法)(求解二元一次不定方程/单个线性同余方程)

不定方程 ax+by=c (或者 ax\equiv c\pmod{b} 设 ax=c+(-y)\times b 移项可得不定方程 ax+by=c )

我们先看一下简单一点的一个方程 ax+by=\gcd(a,b)

证明 \boxed{\gcd(a,b)=\gcd(b,a\bmod b)}

设 g=\gcd(a,b)
则 a=g\times k_1 , b=g\times k_2，
则a\bmod b=(g\times k_1)\bmod(g\times k_2)=(k_1\bmod k_2)\times g， b=g\times k_2

为了使 \gcd(a,b)=\gcd(b,a\bmod b)

假设 $k_1\bmod k_2$ 与 $k_2$ 有公因子为 $r

则k_1\bmod k_2=n_1\times r ,k_2=n_2\times r
所以 k_1=x\times k_2+n_1\times r=(x\times n_2+n_1)\times r

此时 \gcd(k_1,k_2)=r 不符合 \gcd(a,b)=g (否则\gcd(a,b)=r\times g)

所以 k_1\bmod k_2 与 k_2 没有公因子
即 \gcd(a,b)=\gcd(b, a\bmod b)
证毕

证明 \boxed{\texttt{当} a,b\texttt{不全为0时}ax+by=\gcd(a,b)\texttt{一定有解}}

已证:\gcd(a,b)=\gcd(b,a\bmod b)
所以我们可以得到推论：

ax+by=\gcd(a,b)$ 可以转化为 $bx+(a\bmod b)y=\gcd(b,a\bmod b)

所以当 b=0 时(欧几里得算法最后一步)
方程变成 ax=\gcd(a,b)=a
所以一定有一组解 x=1,y=0
证毕

利用这个方法我们可以递归到边界(b=0),然后将 x 设为 1，y 设为 0

然后在递归返回的时候更新 x 和 y 的值

由于方程 ax+by=\gcd(a,b) 转化为了方程 bx+(a\bmod b)y=\gcd(b,a\bmod b)
而 a\bmod b=a-b\times (a/b)除号当然是整除啦

方程左边 bx+(a-b\times (a/b))y=ay+(x-(a/b)\times y)b

所以转化后新的 \boxed{x=y},新的 \boxed{y=x-(a/b)\times y}

但是代码实现可以更简单一些,利用传引用在访问下一层递归的时候就交换 x 和 y
Code

int g,x,y;
void Exgcd(int a,int b,int &x,int &y)
{
    if(!b)
    {
        x=1;y=0;g=a;
        return ;
    }
    Exgcd(b,a%b,y,x);//交换
    y-=(a/b)*x;
    return ;
}

那么解决了方程 ax+by=\gcd(a,b) 之后 ax+by=c 也很好解决了

我们解出 ax+by=\gcd(a,b) 后，如果 c 不是 \gcd(a,b) 的整数倍，那么这个同余方程是没有解的。

如果有解，那么解为 x\times (c/\gcd(a,b))

通解为 x_i=x\times (c/\gcd(a,b))+(b/\gcd(a,b))\times k (k \in \mathbb Z)

利用 Exgcd 可以解决较为基本的不定方程和同余方程问题

时间复杂度 O(\log n) (证明可以自己去了解)

Day 2.\mathsf{CRT}(中国剩余定理/孙子定理)(求解模数两两互质的线性同余方程组)

(线性的意思就是一次的方程) 同余方程组长这样：

\begin{cases} x\equiv r_1\pmod {m_1}\\x\equiv r_2\pmod {m_2} \\ \cdots \\x\equiv r_i\pmod {m_i}\\ \cdots \\x\equiv r_n\pmod {m_n}\end{cases}

其中 m_1,m_2,\cdots,m_n 两两互质 中国剩余定理：对于如上的线性同余方程组，若其中 m_1,m_2,\cdots,m_n 两两互质 设 M=\prod\limits_{i=1}^{n}m_i，M_i=M/m_i 设 K_i 为 M_iK_i\equiv 1\pmod {m_i}的解 则对于任意的 r_1,r_2,\cdots,r_n 方程组有整数解 \boxed{x=\sum\limits_{i=1}^{n}r_iM_iK_i}

证明中国剩余定理
因为 M_i=M/m_i 是除了 m_i 以外的所有的 m 的倍数(模数两两互质)

所以 \forall k \ne i, r_iM_iK_i\equiv 0 \pmod {m_k} (都说了是倍数那么右边肯定是 0 了）

又因为 K_i 为 M_iK_i\equiv 1\pmod {m_i}的解，左右同时乘以 r_i 得到 r_iM_iK_i \equiv r_i\pmod {m_i}

所以当 x=\sum\limits_{i=1}^{n}r_iM_iK_i 时,对于任意的 x\equiv r_i\pmod {m_i},因为\forall k \ne i， r_iM_iK_i\equiv 0 \pmod {m_k} 而且 r_iM_iK_i \equiv r_i\pmod {m_i}，所以方程成立。

即方程组的解为 x=\sum\limits_{i=1}^{n}r_iM_iK_i
证毕

当题目保证给出的模数两两互质的时候，\mathsf{CRT\text{ }it\text{ }is!}
步骤：

  1. 求出M=\prod\limits_{i=1}^{n}m_i
  2. 利用 Exgcd 求出 M_iK_i\equiv 1\pmod {m_i}的解 K_i
  3. 逐步累加 r_iM_iK_i (即利用 x=\sum\limits_{i=1}^{n}r_iM_iK_i)得到解(模 M 意义下)

Code： 

typedef long long LL;
LL CRT()
{
    LL pro=1;
    for(LL i=1;i<=N;i++)
        pro*=m[i];
    LL x,y,Mi;LL ans=0;
    for(LL i=1;i<=N;i++)
    {
        Mi=pro/m[i];
        Exgcd(Mi,m[i],x,y);//x即为求解的K_i
        ans=(ans+Mi*r[i]*x)%pro;
    }
    return (ans+pro)%pro;
}

时间复杂度 O(n\log n)

Day 3.\mathsf{ExCRT}(扩展中国剩余定理/孙子定理)(求解一般的线性同余方程组)

\mathsf{ExCRT}$ 实际上 $\mathsf{superset}$ 了 $\mathsf{CRT.}\space\space\mathcal{By\space connect}

不是所有的方程组模数都两两互质(正如不是所有牛奶...算了) 那么现在所有的 m 之间都没有什么关系了
这样的话方程 M_iK_i\equiv 1\pmod {m_i} 不一定有解，无法通过 x=\sum\limits_{i=1}^{n}r_iM_iK_i 得到解

好的那么现在你面对的是？

\begin{cases} x\equiv r_1\pmod {m_1}\\x\equiv r_2\pmod {m_2} \\ \cdots \\x\equiv r_i\pmod {m_i}\\ \cdots \\x\equiv r_n\pmod {m_n}\end{cases}

Q:现在你会什么？ A:解一个方程

假设我们解出来了第一个方程，将其解设为 ans ,我们知道第一个方程的通解 x_v=ans+k\times m_1,k\in \mathbb{Z}

我们现在的目的就是要求一个同时满足第一第二两个方程的 x ，那么我们把第一个方程的通解代入第二个方程

ans+k\times m_1\equiv r_2\pmod {m_2}$ 移项得$k\times m_1\equiv r_2-ans\pmod {m_2}

即求解不定方程 m_1\times k+m_2\times y=r_2-ans

解出来 k 之后得到第一第二两个方程的公共解 ans^\prime=ans+k\times \text{lcm}(m_1,m_2)

那么前两个同余方程等价为：

x\equiv ans^\prime\pmod {\text{lcm}(m_1,m_2)}

为什么等价？ 因为向方程中代入 x_v=ans^\prime+k\times \text{lcm}(m_1,m_2 ),k\in \mathbb{Z}

你会发现这两个方程和上面的等价方程同时成立，所以这里可以等价

那么你就会发现一件事情，你成功地把两个方程整成了一个方程

那么下一步，继续从方程组里面拿一个方程出来，这时候你手上又有了两个方程

\begin{cases} x\equiv ans\pmod{\text{lcm}}\\x\equiv r_i\pmod{m_i} \end{cases}

由于第一个方程的通解是 x_v=ans+k\times \text{lcm},k\in \mathbb{Z}

所以我们就是要求一个满足条件的 x_v 使 x_v\equiv r_i\pmod{m_i} 成立

所以就是求解 ans+k\times \text{lcm}\equiv r_i\pmod{m_i}

变形得 k\times \text{lcm}\equiv r_i-ans\pmod{m_i}

即求解不定方程 \text{lcm}\times k+m_i\times y=r_i-ans

如果方程无解，则同余方程组无解

方程有解，那么解 ans^\prime=ans+k\times \text{lcm}

循环求解即可。

形式化描述过程就是：

设 ans 为前 i-1 个同余方程的公共解，\text{lcm} 是前 i-1 个方程模数的最小公倍数

那么对于第 i 个方程我们求解不定方程 \text{lcm}\times k+m_i\times y=r_i-ans

不定方程无解则同余方程组无解

得到 k 后我们知道了前 i 个同余方程的公共解 ans^\prime =ans+\text{lcm}\times k

前 i 个方程模数的最小公倍数 \text{lcm}^\prime =\text{lcm}\times(m_i/\gcd(\text{lcm},m_i))

此时可求解第 i+1 个方程，循环即可。

\mathsf{Code} 
inline LL ExCRT()
{
    LL ans=0,lcm=1;
    LL k1,k2,c;
    for(int i=1;i<=N;i++)
    {
        Exgcd(lcm,m[i],k1,k2);
        c=((r[i]-ans)%m[i]+m[i])%m[i]; // c=r2-r1(顺便调整为正的)
        if(c%g) return -1; //无解
        k1=((k1*c/g)%(m[i]/g)+(m[i]/g))%(m[i]/g);
        ans=ans+k1*lcm;
        lcm=lcm*m[i]/g; //lcm=lcm(lcm,m_i)
        ans=(ans%lcm+lcm)%lcm;      
    }
    return (ans%lcm+lcm)%lcm;
}

时间复杂度 O(n\log n) （常数比 \mathsf{CRT} 大一点）

Day 4.\mathsf{BSGS}(\mathsf{Baby Step Giant Step})(解未知数在指数位置且模数为质数的同余方程)

有时候事情并不会都是那么简单，经过之前几个小节内容的学习你应该已经可以熟练地解决线性同余方程的相关问题 但是啊

A^x\equiv B\pmod{C}

其中 A,B,C 已知 且 A,C 互质，需要求解最小的 x 这个该怎么办呢？

\mathsf{BSGS\text{ }it\text{ }is!}

首先我们想一下对于这个问题暴力怎么搞

当然直接枚举指数啊qwq

但是这样枚举的尽头是什么呢？

首先我们知道 A^0\equiv 1\pmod{C} ，然后又由费马小定理(不知道这个的同学可以先去了解一下逆元再了解费马小定理)

A^{C-1}\equiv 1\pmod{C}

所以我们发现实际上枚举指数的时候枚举 C 个就会开始重复，那么暴力枚举 \left[0, C-1\right) 即可。 

int a,m,r;
scanf("%d%d%d",&a,&m,&r);
int base=1;
for(int i=0;i<m-1;i++)
{
    if(base==r) {printf("%d",i);return 0;}
    base=(base*a)%m;
}
printf("no solution");
return 0;

但是假如现在 C=19260817 但是只给你 200ms 这个暴力算法是会超时的（而且题目要是想卡暴力质数可以出的很大）

那么我们要对暴力算法进行改进，核心肯定是缩小枚举的范围

这就是 \mathsf{BSGS} 的雏形了

设 x=i\times \sqrt{C}-j

对于同余方程 A^{i\times \sqrt{C}-j}\equiv B\pmod{C}

左右两边同时除以 B 得 \dfrac{A^{i\times \sqrt{C}}}{A^j\times B}\equiv 1\pmod{C}

预处理分母，然后当模 C 意义下有 A^{i\times \sqrt{C}}=A^j\times B 时\dfrac{A^{i \times \sqrt{C}}}{A^j\times B}=1, 方程有解 i\times \sqrt{C}-j

分母怎么预处理呢？

因为我们的解形式是 i\times \sqrt{C}-j ，所以 j\in \left[0, \sqrt{C}\right)

预处理出 (A^j\times B)\bmod C,j\in\left[0, \sqrt{C}\right) ,将其幂指数存在哈希表里面

暴力枚举 i,i\in\left[ 1,\sqrt{C}\right] ，对于每一个 A^{i\times \sqrt{C}} 在哈希表里面查询是否有对应的 A^j\times B

若存在那么最小解就是 \boxed{i\times \sqrt{C}-\small\texttt{哈希表中存储的幂指数}}

那么现在是不是能理解 \mathsf{BabyStepGiantStep} 这个名字了呢？

我们在预处理哈希表的时候就是在走 \mathsf{BabyStep} (一次加 1)

枚举 i\times\sqrt{C} 中的 i 时就是在走 \mathsf{GiantStep} (一次加 \sqrt{C} )

\mathsf{Code} 
LL BSGS()
{
    Hash.clear();
    //很多选手会使用map但是手写Hash的效率比map高了好几倍
    LL D=B%C,sqrtm=ceil(sqrt(C));
    for(int i=0;i<sqrtm;i++)
    {
        //预处理幂值，建立幂值到幂指数的映射
        Hash.insert(D,i);
        D=D*A%C;
    }
    LL t=ti=fast_pow(A,sqrtm),ans;
    //快速幂求出A^sqrtm
    for(int i=1;i<=sqrtm;i++)
    {
        if((ans=Hash.find(ti))!=-1)
         //如果找到了
            return i*sqrtm-ans;
        ti=ti*t%C;
    }
    return -1;
}

时间复杂度 O(\sqrt{C}) (快速幂的 \log n 在 \sqrt{C} 下忽略不计)

Day5.\mathsf{ExBSGS}(\mathsf{ExtendedBabyStepGiantStep})(解未知数在指数位置的一般模数同余方程)

就像我们在 \mathsf{ExCRT} 的时候说过的一样，不是所有的模数都那么美好让 A,C 互质

不互质的时候分式 \dfrac{A^{i \times \sqrt{C}}}{A^j\times B} 会因为 A^j 没得逆元导致无意义

我们只会互质的解法，怎么办？

没有条件创造条件！

对于同余方程 A^x\equiv B\pmod{C},A,C 不互质

那么我们让他互质！

设 g=\gcd(A,C) ,我们需要对方程执行消去因子

设 A=a\times g,C=c\times g 我们必须要把方程中的幂底数和模数化为互质的，如果途中发现 B\bmod g\ne 0 那么这个方程无解

证明 \boxed {\texttt{若}A=a\times g,C=c\times g,g\ne 1.\texttt{而}B\bmod g\ne0,\texttt{方程无解}}

设 D=A^{x-1} 即有同余方程 A\times D\equiv B\pmod{C}

化为不定方程 ax+by=c 形式：A\times D+C\times y=B

在学习 \mathsf{Exgcd} 的时候我们就已经知道了 上述方程有解的充分必要条件就是 B 可以被 \gcd(A,C) 也就是 g 整除，而现在 B\bmod g\ne 0 所以方程无解

证毕

为了使模数互质我们执行消去因子达到 \gcd(A,C)=1 为止

依据 \mathsf{BSGS} 的经验我们再来

对于方程 A^x\equiv B\pmod{C}

设 A=a\times g,C=c\times g,g\ne 1 ,全除以 g 得到

最后可以得到这样一个柿子 : $D\times A^{x-cnt}\equiv B^\prime \pmod{C^\prime}

其中的 D 是每消去一次因子左边的值除掉一个 g 后得到的商的乘积 （即 D=\prod\limits_{i=1}^{cnt}(A/g_i), g_i为第 i 次消去因子时的 \gcd(A,C) )

设 x=i\times\sqrt{C^\prime}-j+cnt ,两边同时除以 B^\prime 得

\dfrac{D\times A^{i\times \sqrt{C^\prime}}}{A^j\times B^\prime}\equiv 1\pmod{C^\prime}

现在 A,C^\prime 互质了，套上 \mathsf{BSGS} 求解出对应的 i\times\sqrt{C^\prime}-j 然后加上 cnt 即可

但是这样得到的 x>cnt 我们需要对小于 cnt 的解特判，消去因子的时候做即可

\mathsf{Code} 
inline LL BSGS()
{
    LL sqrtm=ceil(sqrt(C));
    H.clear();//手写Hash很香的qwq
    LL t=B%C;
    for(int i=0;i<sqrtm;i++)
    {
        H.insert(t,i);
        t=t*A%C;
    }
    t=fast_pow(A,sqrtm,C);
    LL ti=D*t%C,ans;//初始化成D*A^(i*sqrtm)
    for(int i=1;i<=sqrtm;i++)
    {
        if((ans=H.find(ti))!=-1&&i*sqrtm-ans!=0)
            return i*sqrtm-ans;
        ti=ti*t%C;
    }
    return -1;
}
inline LL ExBSGS()
{
    LL cnt=0;D=1;
    if(B==1) return 0;
    while((g=gcd(A,C))!=1)//消去因子
    {
        if(B%g) return -1;
        B/=g;C/=g;
        D=(D*A/g)%C;
        cnt++;
        if(D==B) return cnt;
        //每消去一次g，D就会乘以A/g 在模C意义下出现D==B时消去了cnt次
        //所以解为cnt
    }
    LL ans=BSGS();
    if(ans==-1) return -1;
    else return ans+cnt;
}

时间复杂度 O(\sqrt{C}) (消去因子部分复杂度是 \log 级别,忽略)

\huge\mathsf{Congratulations!Learning\text{ }Completed}

如果还有任何对算法不理解的问题珂以私信窝，力所能及的范围内窝一定尽力解答 ^_^