【7】同余学习笔记

· · 算法·理论

前言

中国剩余定理是大玄学，我终于懂得了龟速乘有什么用。没什么逻辑

实际上大纲没有专门列出同余本身的整体评级，但是大部分知识点都是 7 级，最后就自己评为 7 级。

UPD on 2025.2.18 ：新增了更深入理解裴蜀定理的内容和一道例题。

长文警告：本文一共 1175 行，请合理安排阅读时间。

欧几里得算法

欧几里得算法

最大公约数：最大公因数，也称最大公约数、最大公因子，指两个或多个整数共有约数中最大的一个，记作 \gcd(a,b) 。

互质数：公因数只有 1 的两个或多个非零自然数，叫做互质数。

求 \gcd(a,b)，使用欧几里得算法。

int gcd(int a,int b)
{
    if(b==0)return a;
    else return gcd(b,a%b);
}

实际上就是辗转相除法，还有辗转相减法，不过没有辗转相除法快。

$1$：对于任意整数 $a\mid b$ ，如果 $d\mid a$ 并且 $d\mid b$ ，则 $d\mid \gcd(a,b) 2$：$\gcd(an,bn)=n\times\gcd(a,b) 3$：对于正整数 $a,b,n$ ，若 $n\mid ab$ 并且 $\gcd(a,n)=1$ ，则 $\gcd(b,n)=1

裴蜀定理

给定一个关于未知数 x,y 的一次不定方程且 x,y 均为整数，求 x,y 的值。

ax+by=c

这个式子也被称之为裴蜀等式。

方程有整数解，当且仅当 c 是 \gcd(a,b) 的倍数，裴蜀等式有解时必然有无穷解，且此时和可以取到 \gcd(a,b) 的任意整数倍。在辗转相除法下任意情况，裴蜀等式均成立。这称为裴蜀定理。下面给出证明。

原问题等价于 ax+by 能得到的最小正整数为 \gcd(a,b)。

首先有一点是显然的，最终的答案一定为 \gcd(a,b) 的倍数。因为 ax+by=ak_1\gcd(a,b)+bk_2\gcd(a,b)=(ak_1+bk_2)\gcd(a,b)，其中 k_1,k_2 均为整数。显然 ak_1+bk_2 为整数。

设 ax+by 能得到的最小正整数为 c，此时有 ax_0+by_0=c，考虑带余除法，假设 a 不能整除 c，设 a=pc+r，其中 p,r 均为整数且 0\le r\lt c。

于是 r=a-pc=a-p(ax_0+by_0)=a(1-px_0)-bpy_0。令 x'=1-px_0,y'=py_0，可以整理为 ax'+by' 的形式，又变为了裴蜀等式。

因为 r\ge 0，所以 ax'+by'\ge 0，而根据条件 ax+by 能得到的最小正整数为 c，所以 ax'+by'\ge c，即 r\ge c。而根据假设 r\lt c，矛盾，所以 a 能整除 c。所以 c 必然为 a 的因数。

同理假设 b 不能整除 c 也能反证出相同的结论，c 也必然为 b 的因数。因此，c 同时为 a,b 的因数，所以 c 是 \gcd(a,b) 的因数。

综上所述，c 是 \gcd(a,b) 的倍数，也是 \gcd(a,b) 的因数，所以 c=\gcd(a,b)。结论得证。

裴蜀定理的一个常见推论是，正整数 k 的任意整数倍在模 p 意义下可以取到的值为 \gcd(k,p) 的整数倍。可以使用以下式子表示，其中 n,m 均为整数。

nk\equiv m\gcd(k,p)\ ({\rm mod}\ p)

扩展欧几里得算法

我们讨论裴蜀等式较为特殊的情况，如下所示：

ax+by=\gcd(a,b)

在欧几里得算法的最后，一定会有 a=\gcd(a,b),b=0 的情况。在这种情况下，一定有 x=1,y=0。

我们逆推扩展欧几里得算法的过程，则有 \gcd(a,b)=\gcd(b,a\bmod b)。则有如下式子：

bx_1+(a\bmod b)y_1=bx_1+(a-\lfloor\frac{a}{b}\rfloor\times b)y_1=ay_1+b(x_1-\lfloor\frac{a}{b}\rfloor\times y_1)=\gcd(a,b) ax+by=\gcd(a,b)

对比二式系数得到 x=y_1,y=x_1-\lfloor\frac{a}{b}\rfloor\times y_1，这样，我们就得到了 x,y 的递推式。求完 a,b 后倒着递推回去，即可得到 x,y 的值。

对于更一般的 ax+by=c 的情况，我们只需要把等式两边同时乘以 \frac{c}{\gcd(a,b)} 即可。也就是说，如果扩展欧几里得算法求出来 x',y'，那么最终的 x=x\times\frac{c}{\gcd(a,b)},y=y\times\frac{c}{\gcd(a,b)}。

扩展欧几里得算法既可以求出 \gcd(a,b) ，又可以求出裴蜀等式的 x,y 。

int exgcd(int a,int b,int &x,int &y)
{
    if(b==0)
       {
        x=1;
        y=0;
        return a;
       }
    int r=exgcd(b,a%b,x,y);
    int d=x;
    x=y;
    y=d-a/b*y;
    return r;
}

中国剩余定理

中国剩余定理

求解如下同余方程组：(最小非负整数解)

\begin{cases} x \equiv b_1\ ({\rm mod}\ a_1) \\ x\equiv b_2\ ({\rm mod}\ a_2) \\ ... \\ x \equiv b_n\ ({\rm mod}\ a_n)\end{cases}

其中 a_i,b_i 为给定的整数，且保证 b_i 两两互质。

令 M=\prod_{i=1}^{n}a_i，则对于每个式子，求出 m_i=\frac{M}{a_i} 在模 a_i 意义下的逆元记作 t_i，原方程组的解为 \sum_{i=1}^{n}b_im_it_i。

扩展中国剩余定理

求解如下同余方程组：(最小非负整数解)

\begin{cases} x \equiv b_1\ ({\rm mod}\ a_1) \\ x\equiv b_2\ ({\rm mod}\ a_2) \\ ... \\ x \equiv b_n\ ({\rm mod}\ a_n)\end{cases}

其中 a_i,b_i 为给定的整数，不保证 b_i 两两互质。

考虑两两合并同余方程组，具体的，对于如下两个式子，我们希望将其合并为一个同余式。

\begin{cases} x \equiv b_1\ ({\rm mod}\ a_1) \\ x\equiv b_2\ ({\rm mod}\ a_2)\end{cases}

由小学学过的知识，我们知道，合并之后的模数 a 等于 \text{lcm}(a_1,a_2)。

我们首先将两个式子用同余式的定义展开得到：

\begin{cases} x = b_1 + a_1k_1 \\ x = b_2 - a_2k_2\end{cases}

在这一步中，下面的式子展开为 x = b_2+a_2k_2 也可以，但是为了方便下一步计算，所以选择拆成减号的形式。

由于两个式子右边都等于 x，代入得到：

b_1 + a_1k_1=b_2 - a_2k_2

移项得到如下式子：

a_1k_1+a_2k_2=b_2-b_1

这一是个经典的裴蜀等式，直接运用 exgcd 求解即可。最后把求出的值带回去合并时最开始方程组两个式子中的第一个式子，得到 x=b_1+a_1k_1。

经过多次合并，最后整个方程组被合并为了一个式子，如下：

x \equiv b_n\ ({\rm mod}\ a_n)

显然，最小非负整数解一定为 b_n。

代码中 mul 为慢速乘，作用是防止溢出，相当于 x=x*(c/q)%(a[i]/q)。

代码中的 nm 相当于上述推导过程中的 a_1，ans 相当于上述推导过程中的 b_1。需要注意的是，当一次合并过程中裴蜀等式无解时，整个方程组也不会有解，返回 -1。

long long excrt()
{
    long long nm=a[1],ans=b[1];
    for(int i=2;i<=n;i++)
        {
            long long c=((b[i]-ans)%a[i]+a[i])%a[i];
            long long q=exgcd(nm,a[i],x,y);
            if(c%q!=0)return -1;
            x=mul(x,c/q,a[i]/q);
            ans+=x*nm,nm*=(a[i]/q),ans=(ans%nm+nm)%nm; 
        }
    return (ans%nm+nm)%nm;
}

应用

带系数的同余方程组(见例 5)

square-free-number 降低模数(见例 10)

欧拉定理

欧拉定理及其推论

欧拉定理：若正整数 a,n 互质，则 a^{\varphi(n)} \equiv 1\ ({\rm mod}\ n)，其中 \varphi(n) 为欧拉函数。

欧拉定理的推论：若正整数 a,n 互质，则对于任意正整数 b，有 a^{b} \equiv a^{b\mod \varphi(n)}\ ({\rm mod}\ n)。

证明：

设 b=q\times \varphi(n)+r，其中 0\le r\lt\varphi(n)，即 r=b\mod\varphi(n)。

a^b\equiv a^{q\times\varphi(n)+r}\equiv(a^{\varphi(n)})^{q}\times a^r\equiv 1^q\times a^r\equiv a^r\equiv a^{b\mod \varphi(n)}\ ({\rm mod}\ n)

扩展欧拉定理

对于任意正整数 a,n,b，有如下式子：

a^b\equiv\begin{cases} a^b (b\lt\varphi(n))\\ a^{b\mod \varphi(n)+\varphi(n)}(b\ge\varphi(n))\end{cases}\ ({\rm mod}\ n)

BSGS

BSGS(大步小步算法)

给定正整数 a,b,p，且保证 a,p 互质，求下列高次同余方程的最小非负整数解：

a^x\equiv b\pmod p

由欧拉定理，得到 a^{x}\equiv a^{x\mod\varphi(n)}\pmod p，因此，该方程的最小非负整数解一定在 [0,\varphi(p)] 之间，否则可以通过欧拉定理得到一个更小的解。又因为 \varphi(p)\le p-1，所有最小非负整数解一定在 [0,p-1] 之间。

考虑根号分治。令 t=\lfloor\sqrt{p}\rfloor，x=i\times t-j，其中 0\le j\lt t。则原方程可化为：

a^{i\times t-j}\equiv b\pmod p

两边同时乘以 a^j 得到如下式子：

a^{i\times t}\equiv a^jb\pmod p

考虑分别枚举左右两边。先枚举右边，0\le j\lt t，通过递推求出每一个 a^jb\mod p 的值，插入一个 Hash 表。接下来，枚举左边。先通过快速幂求出 a^t\mod p 的值，然后每次将上一次的余数乘以这个值，即可求出每一个 a^{i\times t} 的值。由于 x 在 [0,p-1] 之间，所以 i 也只需要枚举 [1,t]。枚举出值之后，在 Hash 表中进行匹配，如果成功，直接返回。如果一次都没有成功，返回无解。

总的时间复杂度为 O(\sqrt{q})，较为优秀。

代码中 h.count(x) 表示数值 x 在 Hash 表中的出现次数，主要用于卡常。

unordered_map<long long,long long>h;
long long bsgs(long long a,long long b,long long mod)
{
    h.clear();
    long long t=sqrt(mod)+1,now=1,inc=power(a,t,mod);
    if(a==0)
       {
        if(b==0)return 1;
        else return -1;
       }
    for(int i=0;i<=t-1;i++)h[now*b%mod]=i+1,now=now*a%mod;
    now=1;
    for(int i=1;i<=t;i++)
        {
        now=now*inc%mod;
        if(h.count(now)>0)return i*t-h[now]+1;
        }
    return -1;
}

exBSGS(扩展大步小步算法)

给定正整数 a,b,p，且不保证 a,p 互质，求下列高次同余方程的最小非负整数解：

a^x\equiv b\pmod p

考虑将原问题转化为 a,p 互质的情况，先将原方程转化为如下形式：

a^{x-1}a+kp=b

设 \gcd(a,p)=d，则根据裴蜀定理，必须有 d\mid b，否则直接返回无解。然后，考虑把 d 除掉，得到如下式子：

a^{x-1}\times\frac{a}{d}+k\times\frac{p}{d}=\frac{b}{d}

令 a^{x-1}\to a^x,\frac{p}{d}\to p,\frac{b}{d}\to b，重复这个过程，每一次单独提出来一个 a，最后使得 \gcd(a,p)=1。

设最后一共进行了这个过程 cnt 次，所有过程 d 的乘积为 d'，则最后可以得出以下式子：

a^{x-cnt}\times\frac{a^{cnt}}{d'}+k\times\frac{p}{d'}=\frac{b}{d'}

根据上述转化，最终 \gcd(a,p)=1，即 \gcd(a,\frac{p}{d'})=1。 转化为同余式，即可得到：

a^{x-cnt}\times\frac{a^{cnt}}{d'}\equiv\frac{b}{d'}\pmod {\frac{p}{d'}}

对这个式子直接求一遍 bsgs，注意左边需要额外乘以一个系数 \frac{a^{cnt}}{d'}。最后求出来的结果还需要增加 cnt，因为我们需要把剪掉的部分加回来。

long long exbsgs(long long a,long long b)
{
    long long d=gcd(a,p),cnt=0,sum=1,ans=0;
    a%=p,b%=p;
    if(b==1||p==1)return 0;
    if(a==0)
       {
        if(b==0)return 1;
        else return -1;
       }
    while(d!=1)
       {
        if(b%d!=0)return -1;
        cnt++,p/=d,sum=sum*(a/d)%p,b/=d;
        d=gcd(a,p);
        if(sum==b)return cnt;
       }
    exgcd(sum%p,p,x,y);
    x=(x%p+p)%p;
    ans=bsgs(a,x*b%p,p);
    if(ans==-1)return -1;
    else return ans+cnt;
}

例题

例题 1 ：

P4549 【模板】裴蜀定理

我们不难把裴蜀定理推广到多个整数的形式，证明过程几乎是完全一样的，最后的和一定是这些整数的 \gcd。注意负数需要取相反数后再求 \gcd。

#include <bits/stdc++.h>
long long n,a,ans;
long long gcd(long long a,long long b)
{
    if(b==0)return a;
    else return gcd(b,a%b);
}

int main()
{
    scanf("%lld",&n);
    scanf("%lld%lld",&a,&ans);
    if(a<0)a=-a;
    if(ans<0)ans=-ans;
    ans=gcd(a,ans);
    for(int i=2;i<n;i++)
        {
            scanf("%lld",&a);
            if(a<0)a=-a;
            ans=gcd(a,ans);
        }
    printf("%lld",ans);
    return 0;
}

例题 2 ：

P1082 [NOIP2012 提高组] 同余方程

ax\equiv 1\pmod b

由同余的性质得到：

ax+by\equiv 1+0\pmod b

合并，去掉模数限制得到：

ax+by=1

然后就转化成了裴蜀等式，扩展欧几里得算法解决。实际上就是求出了 a 在模 b 意义下的逆元。

#include <bits/stdc++.h>
using namespace std;
int n,a,b,x,y; 
int exgcd(int a,int b,int &x,int &y)
{
    if(b==0)
       {
        x=1;
        y=0;
        return a;
       }
    int r=exgcd(b,a%b,x,y);
    int d=x;
    x=y;
    y=d-a/b*y;
    return r;
}

int main()
{
    scanf("%d%d",&a,&b);
    int r=exgcd(a,b,x,y);
    printf("%d\n",(x+b)%b);
    return 0;
}

例题 3 ：

P2421 [NOI2002] 荒岛野人

由于题目保证了答案不会超过 10^6，考虑直接枚举答案 m，判断是否合法。

首先，如果连 \max\{c_i\} 个山洞都没有，肯定不合法，直接跳过这一部分。

然后，考虑枚举每两个野人在有生之年会不会到同一个山洞。抽象成数学，需要解一个同余方程：（其中只有 x 为未知数，表示经过的年数）

c_i+p_ix\equiv c_j+p_jx\pmod m

移项合并得：

(p_i-p_j)x\equiv c_i-c_j\pmod m

采用类似例题 1 的方式得到如下式子：

(p_i-p_j)x+my=c_i-c_j

直接扩展欧几里得算法解决，求出 x。若原方程有解并且 x 小于两个野人的寿命，那么这两个野人在有生之年就会发生冲突，这个答案不合法，枚举下一个。

#include <bits/stdc++.h>
using namespace std;
long long x,y,s,c1[16],p1[16],l1[16],max1=0;
long long gcd(long long a,long long b)
{
    if(b==0)return a;
    else return gcd(b,a%b);
}

long long exgcd(long long a,long long b,long long &x,long long &y)
{
    if(b==0)
       {
        x=1;
        y=0;
        return a;
       }
    long long r=exgcd(b,a%b,x,y);
    long long d=x;
    x=y;
    y=d-a/b*y;
    return r;
}

int check(long long a,long long b,long long m,long long n,long long l)
{
    long long p=n-m,q=a-b;
    if(q%gcd(p,l)!=0)return 0;
    long long r=exgcd(p,l,x,y);
    return ((x*q/r)%(l/r)+(l/r))%(l/r);
}

bool judge(long long a)
{
    for(int i=0;i<s;i++)
        for(int j=i+1;j<s;j++)
            {
                long long b=check(c1[i],c1[j],p1[i],p1[j],a);
                if(b==0)continue;
                if(b<=l1[i]&&b<=l1[j])return 0;
            }
    return 1;
}

int main()
{
    scanf("%lld",&s);
    for(int i=0;i<s;i++)
        {
        scanf("%lld%lld%lld",&c1[i],&p1[i],&l1[i]);
        max1=max(max1,c1[i]);
        }
    for(int i=max1;i<1000000;i++)
        {
            if(judge(i))
               {
                printf("%lld",i);
                return 0;
               }
        } 
    return 0;
}

例题 4 ：

P4777 【模板】扩展中国剩余定理（EXCRT）

扩展中国剩余定理板子题，不多赘述。

#include <bits/stdc++.h>
using namespace std;
long long n,x,y,a[200000],b[200000];
long long mul(long long a,long long p,long long m)
{
    long long x=a,ans=0;
    while(p)
       {
        if(p%2==1)ans=(ans+x)%m;
        p/=2;
        x=(x+x)%m;
       }
    return ans;
}

long long exgcd(long long a,long long b,long long &x,long long &y)
{
    if(b==0)
       {
        x=1;
        y=0;
        return a;
       }
    long long r=exgcd(b,a%b,x,y);
    long long d=x;
    x=y;
    y=d-a/b*y;
    return r;
}

long long excrt()
{
    long long nm=a[1],ans=b[1];
    for(int i=2;i<=n;i++)
        {
            long long c=((b[i]-ans)%a[i]+a[i])%a[i];
            long long q=exgcd(nm,a[i],x,y);
            if(c%q!=0)return -1;
            x=mul(x,c/q,a[i]/q);
            ans+=x*nm,nm*=(a[i]/q),ans=(ans%nm+nm)%nm; 
        }
    return (ans%nm+nm)%nm;
}

int main()
{
    scanf("%lld",&n);
    for(int i=1;i<=n;i++)scanf("%lld%lld",&a[i],&b[i]);
    printf("%lld",excrt());
    return 0;
}

例题 5 ：

P4774 [NOI2018] 屠龙勇士

由于必须按照顺序打每一条龙，每一条龙掉落到剑是一定的，所以如果能打到某一条龙，那么所可以选择的剑的集合是一定的。需要一个可以动态维护前驱，插入，删除元素的数据结构，考虑平衡树。

依题意，巨龙最后会不断以固定的速度回复血量，生命值刚好为 0 时才会死亡。考虑通过一个同余式表示整个过程。

对于第 i 条龙，若其初始血量为 b_i，每次回复 a_i 点血量，打其时剑的攻击力为 c_i，c_i 通过平衡树维护出，则能打死的情况下有如下同余式：

c_ix\equiv b_i\pmod {a_i}

能打死所有的龙，需要所有同余式都有一个解。由题意，这个解是相同的。所以联立所有式子，得到如下方程组：

\begin{cases} c_1x \equiv b_1\ ({\rm mod}\ a_1) \\ c_2x\equiv b_2\ ({\rm mod}\ a_2) \\ ... \\ c_nx \equiv b_n\ ({\rm mod}\ a_n)\end{cases}

接下来，就需要求解带系数的同余方程组。不可以直接将 c_i 除过去，因为有可能逆元不存在。依旧考虑合并方程，方便起见，我们让已经合并的式子系数为 1，则需要合并如下式子：

\begin{cases} x \equiv b_1\ ({\rm mod}\ a_1) \\ c_2x\equiv b_2\ ({\rm mod}\ a_2)\end{cases}

由小学学过的知识，我们知道，合并之后的模数 a 等于 \text{lcm}(a_1,a_2)。

我们首先将两个式子用同余式的定义展开得到：

\begin{cases} x = b_1 + a_1k_1 \\ c_2x = b_2 - a_2k_2\end{cases}

将上面式子乘以 c_2 得到：

\begin{cases} c_2x = b_1c_2 + a_1c_2k_1 \\ c_2x = b_2 - a_2k_2\end{cases}

由于两个式子右边都等于 c_2x，代入得到：

b_1c_2 + a_1c_2k_1=b_2 - a_2k_2

移项得到如下式子：

(a_1c_2)k_1+(a_2)k_2=(b_2-b_1c_2)

这一是个经典的裴蜀等式，直接运用 exgcd 求解即可。最后把求出的值带回去合并时最开始方程组两个式子中的第一个式子，得到 x=b_1+a_1k_1。

经过多次合并，最后整个方程组被合并为了一个式子，如下：

x \equiv b_n\ ({\rm mod}\ a_n)

显然，最小非负整数解一定为 b_n。

注意多使用慢速乘，防止溢出。同时注意，必须要把巨龙的血量砍为负数。如果最后求出的 x 不足以将某只巨龙血量砍为负数，特判将 x 改为可以砍成负数的最小值。

#include <bits/stdc++.h>
using namespace std;
long long t,n,m,x,y,s[200040],a[200040],b[200040],c[200040],z[200040],ch[200040][2],val[200040],key[200040],siz[200040],tol[200040],cnt=0,root=0;
long long create(long long v)
{
    val[++cnt]=v;
    key[cnt]=rand();
    siz[cnt]=1;
    tol[cnt]=1;
    return cnt;
}

void pushup(long long now)
{
    siz[now]=siz[ch[now][0]]+siz[ch[now][1]]+tol[now];
}

void build()
{
    root=create((long long)-1999999999999),ch[root][1]=create((long long)1999999999999);
    pushup(root);
}

void rotate(long long &now,long long to)
{
    long long tmp=ch[now][to^1];
    ch[now][to^1]=ch[tmp][to];
    ch[tmp][to]=now;
    now=tmp;
    pushup(ch[now][to]);pushup(now);
}

void insert(long long &now,long long v)
{
    if(now==0)
       {
       now=create(v);
       return;
       }
    if(v==val[now])tol[now]++;
    else 
       {
        long long to=0;
        if(v<val[now])to=0;
        else to=1;
        insert(ch[now][to],v);
        if(key[now]<key[ch[now][to]])rotate(now,to^1);
       }
    pushup(now);
}

void del(long long &now,long long v)
{
    if(now==0)return ;
    if(v==val[now])
       {
       if(tol[now]>1)
          {
          tol[now]--;
          pushup(now);
          return ;
          }
       else if(ch[now][0]||ch[now][1])
            {
                if(!ch[now][1]||key[ch[now][0]]>key[ch[now][1]])rotate(now,1),del(ch[now][1],v);
                else rotate(now,0),del(ch[now][0],v);
                pushup(now);
            }
       else now=0;
       return ;
       }
    else 
       {
        long long to=0;
        if(v<val[now])to=0;
        else to=1;
        del(ch[now][to],v);
        pushup(now);
       }
}

long long pre(long long x)
{
    long long now=root,ans=0;
    while(now)
       {
        if(x<val[now])now=ch[now][0];
        else ans=val[now],now=ch[now][1];
       }
    return ans;
}

long long nxt(long long x)
{
    long long now=root,ans=0;
    while(now)
       {
        if(x>=val[now])now=ch[now][1];
        else ans=val[now],now=ch[now][0];
       }
    return ans;
}

long long exgcd(long long a,long long b,long long &x,long long &y)
{
    if(b==0)
       {
        x=1,y=0;
        return a;
       }
    long long r=exgcd(b,a%b,x,y),d=x;
    x=y,y=d-a/b*y;
    return r;
}

long long mul(long long a,long long p,long long m)
{
    long long x=a,ans=0;
    while(p)
       {
        if(p&1)ans=(ans%m+x%m)%m;
        p>>=1;
        x=(x%m+x%m)%m;
       }
    return ans;
}

long long excrt()
{
    long long nm=1,ans=0;
    for(long long i=1;i<=n;i++)
        {
            long long d=((b[i]-mul(ans,c[i],a[i]))%a[i]+a[i])%a[i],r=exgcd(mul(nm,c[i],a[i]),a[i],x,y);
            if(d%r!=0)return -1;
            x=(x%a[i]+a[i])%a[i];
            ans=ans+mul(nm,mul(x,(d/r),a[i]/r),a[i]/r*nm),nm=a[i]/r*nm,ans=(ans%nm+nm)%nm;
        }
    return ans;
}

int main()
{
    scanf("%lld",&t);
    while(t--)
       {
       scanf("%lld%lld",&n,&m);
       memset(ch,0,sizeof(ch));
       root=0,cnt=0;build();
       for(long long i=1;i<=n;i++)scanf("%lld",&s[i]),b[i]=s[i];
       for(long long i=1;i<=n;i++)scanf("%lld",&a[i]);
       for(long long i=1;i<=n;i++)scanf("%lld",&z[i]);
       for(long long i=1;i<=m;i++)
           {
           scanf("%lld",&x);
           insert(root,x);
           }
       for(long long i=1;i<=n;i++)
            {
            c[i]=pre(b[i]);
            if(c[i]==(long long)-1999999999999)c[i]=nxt(c[i]);
            del(root,c[i]);
            insert(root,z[i]);
            }
        long long ans=excrt();
        if(ans==0)
           for(long long i=1;i<=n;i++)
               ans=max(ans,(s[i]-1)/c[i]+1);
        printf("%lld\n",ans);    
       }
    return 0;
}

例题 6 ：

P5091 【模板】扩展欧拉定理

根据扩展欧拉定理，本质上就是需要求出次数除以模数的欧拉函数值的余数。我们可以逐位读入，利用余数的可乘性和可加性，每次将当前余数乘以 10 在加上读入的数字，再取余。

求出余数后，直接使用快速幂解决。

注意当 b\lt\varphi(m) 时，不需要再将次数增加 \varphi(m)。

#include <bits/stdc++.h>
using namespace std;
long long a,m,n,now,flag,mod;
char b;
long long power(long long a,long long p)
{
    long long ans=1,x=a;
    while(p)
       {
        if(p&1)ans=ans*x%mod;
        p>>=1;
        x=x*x%mod;
       }
    return ans;
}

int main()
{
    scanf("%lld%lld",&a,&m);
    mod=n=m;
    for(long long i=2;i*i<=m;i++)
        if(m%i==0)
           {
            n=n/i*(i-1);
            while(m%i==0)m/=i;
           }
    if(m!=1)n=n/m*(m-1);
    while(b<'0'||b>'9')b=getchar();
    while(b>='0'&&b<='9')
       {
        long long t=(now*10+b-'0')%n;
        if(t<now*10+b-'0')flag=1;
        now=t;
        b=getchar();
       }
    if(flag)printf("%lld",power(a,now+n));
    else printf("%lld",power(a,now));
    return 0; 
}

例题 7 ：

P4139 上帝与集合的正确用法

实际上，题目就是要求以下式子的值：

2^{2^{2^{\dots}}}\mod p

由于幂的层数可以看作无限大，所以 2^{2^{2^{\dots}}} 一定是大于 \varphi(p)的。根据扩展欧拉定理，将式子化为如下形式：

2^{(2^{2^{\dots}}\mod\varphi(p)+\varphi(p))}\mod p

很显然，中间那个 2^{2^{\dots}}\mod\varphi(p) 可以采取类似的操作求出，所以我们直接递归求解，最后加上 \varphi(p) 进行快速幂，就是这一层的结果。

当模数为 1 时，停止递归。因为 \varphi(1)=1，再递归下去，模数始终为 1，得出的值刚好为题目所要求的定值。

#include <bits/stdc++.h>
using namespace std;
long long t,n,pri[1000010],phi[10000010],cnt,ans;
bool a[10000010];
long long power(long long a,long long p,long long mod)
{
    long long x=a,ans=1;
    while(p)
       {
        if(p&1)ans=ans*x%mod;
        p>>=1;
        x=x*x%mod;
       }
    return ans;
}

void linear(long long n)
{
    a[1]=1;phi[1]=1;
    for(long long i=2;i<=n;i++)
        {
            if(!a[i])pri[++cnt]=i,phi[i]=i-1;
            for(long long j=1;j<=cnt&&i*pri[j]<=n;j++)
                {
                a[i*pri[j]]=1;
                phi[i*pri[j]]=phi[i]*(pri[j]-1);
                if(i%pri[j]==0)
                   {
                   phi[i*pri[j]]=phi[i]*pri[j];
                   break;
                   }
                }
        }
}

long long oula(long long now)
{
    if(now==1)return 0;
    else return power(2,oula(phi[now])%now+phi[now],now);
}

int main()
{
    scanf("%lld",&t);
    linear(10000008);
    while(t--)
        {
        scanf("%lld",&n);
        printf("%lld\n",oula(n));
        }
    return 0;
}

例题 8 ：

P3846 [TJOI2007] 可爱的质数/【模板】BSGS

bsgs 模板题，不多赘述。

#include <bits/stdc++.h>
using namespace std;
long long p,b,n,ans;
unordered_map<long long,long long>h;
long long power(long long a,long long p,long long mod)
{
    long long ans=1,x=a;
    while(p)
       {
        if(p%2==1)ans=ans*x%mod;
        p/=2;
        x=x*x%mod;
       }
    return ans; 
}

long long bsgs(long long a,long long b,long long mod)
{
    h.clear();
    long long t=sqrt(mod)+1,now=1,inc=power(a,t,mod);
    if(a==0)
       {
        if(b==0)return 1;
        else return -1;
       }
    for(int i=0;i<=t-1;i++)h[now*b%mod]=i+1,now=now*a%mod;
    now=1;
    for(int i=1;i<=t;i++)
        {
        now=now*inc%mod;
        if(h.count(now)>0)return i*t-h[now]+1;
        }
    return -1;
}

int main()
{
    scanf("%lld%lld%lld",&p,&b,&n);
    ans=bsgs(b,n,p);
    if(ans==-1)printf("no solution");
    else printf("%lld\n",ans);
    return 0;
}

例题 9 ：

P4195 【模板】扩展 BSGS/exBSGS

exbsgs 模板题，不多赘述。

注意此题数据较强，需要大量特判和卡常技巧，写起来不是很容易。

#include <bits/stdc++.h>
using namespace std;
long long a,p,b,x,y;
map<long long,long long>h;
inline long long read()
{
    long long x=0,f=1;char ch=getchar();
    while (ch<'0'||ch>'9'){if (ch=='-') f=-1;ch=getchar();}
    while (ch>='0'&&ch<='9'){x=x*10+ch-48;ch=getchar();}
    return x*f;
}

long long gcd(long long a,long long b)
{
    if(b==0)return a;
    else return gcd(b,a%b);
}

long long exgcd(long long a,long long b,long long &x,long long &y)
{
    if(b==0)
       {
        x=1,y=0;
        return a;
       }
    long long r=exgcd(b,a%b,x,y),d=x;
    x=y,y=d-a/b*y;
    return r;
}

long long power(long long a,long long p,long long mod)
{
    long long ans=1,x=a;
    while(p)
       {
        if(p&1)ans=ans*x%mod;
        p>>=1;
        x=x*x%mod;
       }
    return ans; 
}

long long bsgs(long long a,long long b,long long mod)
{
    h.clear();
    long long t=sqrt(mod)+1,now=1,inc=power(a,t,mod);
    if(a==0)
       {
        if(b==0)return 1;
        else return -1;
       }
    for(int i=0;i<=t-1;i++)h[now*b%mod]=i+1,now=now*a%mod;
    now=1;
    for(int i=1;i<=t;i++)
        {
        now=now*inc%mod;
        if(h.count(now)>0)return i*t-h[now]+1;
        }
    return -1;
}

long long exbsgs(long long a,long long b)
{
    long long d=gcd(a,p),cnt=0,sum=1,ans=0;
    a%=p,b%=p;
    if(b==1||p==1)return 0;
    if(a==0)
       {
        if(b==0)return 1;
        else return -1;
       }
    while(d!=1)
       {
        if(b%d!=0)return -1;
        cnt++,p/=d,sum=sum*(a/d)%p,b/=d;
        d=gcd(a,p);
        if(sum==b)return cnt;
       }
    exgcd(sum%p,p,x,y);
    x=(x%p+p)%p;
    ans=bsgs(a,x*b%p,p);
    if(ans==-1)return -1;
    else return ans+cnt;
}

int main()
{
    a=read(),p=read(),b=read();
    while(1)
        {
            if(a==0&&b==0&&p==0)break;
            long long ans=exbsgs(a,b);
            if(ans==-1)printf("No Solution\n");
            else printf("%lld\n",ans);
            a=read(),p=read(),b=read();
        }
    return 0;
}

例题 10 ：

P2480 [SDOI2010] 古代猪文

基础数论大杂烩。根据题意，题目就是给定整数 q,n，要求以下式子的值：

q^{\sum_{d\mid n}C_{n}^{d}}\mod 999911659

因为 999911659 是大质数，根据欧拉定理的推论，可以将原式化为如下形式：

q^{\sum_{d\mid n}C_{n}^{d}\mod \varphi(999911659)}\mod 999911659 q^{\sum_{d\mid n}C_{n}^{d}\mod 999911658}\mod 999911659

关键就转化为了求 \sum_{d\mid n}C_{n}^{d}\mod 999911658 的值。很显然，对于求组合数，999911658 作为模数实在是太大了，不方便使用 Lucas 定理，考虑降低模数。

把 999911658 质因数分解，得到 999911658=2\times3\times4679\times35617，它的所有质因子指数均为 1，我们把这种数叫做 square-free-number。

我们可以考虑分别求出 a_1=\sum_{d\mid n}C_{n}^{d}\mod 2，a_2=\sum_{d\mid n}C_{n}^{d}\mod 3，a_3=\sum_{d\mid n}C_{n}^{d}\mod 4679，a_4=\sum_{d\mid n}C_{n}^{d}\mod 35617。由于每一个模数都不大，直接枚举 n 的约数，可以通过 Lucas 定理求出每一个组合数的值，并求和。

然后，我们构造出一个线性同余方程组：

\begin{cases} x \equiv a_1\ ({\rm mod}\ 2) \\ x\equiv a_2\ ({\rm mod}\ 3) \\ x\equiv a_3\ ({\rm mod}\ 4679) \\ x \equiv a_4\ ({\rm mod}\ 35617)\end{cases}

直接使用中国剩余定理求出 x，即为 \sum_{d\mid n}C_{n}^{d}\mod 999911658 的值。这个可以根据 excrt 的合并过程证明，此处略。

最后，用快速幂求出 q^x，就是最后的答案了。

注意，当 q=999911659 时，原式的值必然为 0，需要特判。

#include <bits/stdc++.h>
using namespace std;
long long n,g,x,y,jc[50000],inv[50000],a[5]={0,2,3,4679,35617},b[5],m[5],t[5];
const long long mod=999911659; 
long long mul(long long a,long long p,long long m)
{
    long long x=a,ans=0;
    while(p)
       {
        if(p%2==1)ans=(ans+x)%m;
        p/=2;
        x=(x+x)%m;
       }
    return ans;
}

long long power(long long a,long long p,long long m)
{
    long long ans=1,x=a;
    while(p)
       {
        if(p&1)ans=ans*x%m;
        p>>=1;
        x=x*x%m;
       }
    return ans;
}

long long exgcd(long long a,long long b,long long &x,long long &y)
{
    if(b==0)
       {
        x=1;
        y=0;
        return a;
       }
    long long r=exgcd(b,a%b,x,y);
    long long d=x;
    x=y;
    y=d-a/b*y;
    return r;
}

void cal(long long m)
{
    jc[0]=1,inv[0]=1; 
    for(int i=1;i<=m-1;i++)jc[i]=jc[i-1]*i%m;
    inv[m-1]=power(jc[m-1],m-2,m);
    for(int i=m-2;i>=1;i--)inv[i]=inv[i+1]*(i+1)%m;
}

long long get_c(long long n,long long k,long long m)
{
    if(k>n)return 0;
    return jc[n]*inv[n-k]%m*inv[k]%m;
}

long long lucas(long long n,long long k,long long m)
{
    if(k==0)return 1;
    return get_c(n%m,k%m,m)*lucas(n/m,k/m,m)%m;
}

long long solve(long long m)
{
    long long r=0;
    cal(m);
    for(long long i=1;i*i<=n;i++)
        if(n%i==0)
            {
            if(i!=n/i)r=(r+lucas(n,i,m)+lucas(n,n/i,m))%m;
            else r=(r+lucas(n,i,m))%m;
            }
    return r;
}

long long excrt()
{
    long long nm=a[1],ans=b[1];
    for(int i=2;i<=4;i++)
        {
            long long c=((b[i]-ans)%a[i]+a[i])%a[i];
            long long q=exgcd(nm,a[i],x,y);
            if(c%q!=0)return -1;
            x=mul(x,c/q,a[i]/q);
            ans+=x*nm,nm*=(a[i]/q),ans=(ans%nm+nm)%nm; 
        }
    return (ans%nm+nm)%nm;
}

int main()
{
    scanf("%lld%lld",&n,&g);
    if(g==mod)
       {
       printf("0");
       return 0;
       }
    for(int i=1;i<=4;i++)b[i]=solve(a[i]);
    printf("%lld",power(g,excrt(),mod));
    return 0; 
}

后记

这篇成功取代 【7】Tarjan学习笔记，共 1175 行，成为最长的学习笔记。

upd on 2024.2.14 被 【6】线段树学习笔记 取代。

数论这一块，套路很多，只能说熟能生巧吧。

但是最近几年提高组和国赛好像没怎么考数论。