数论总结

数论是 OI 在数学知识板块最重要最基础的一块，它包含许多基础算法，下面分别讲解。

快速幂

快速幂取模

求 $a ^{n}\bmod p$ 的值。

Luogu1226 快速幂

暴力算法时间复杂度为 $O(n)$，当 $n>10 ^{8}$ 时就不能承受。

有没有更快的算法呢？我们考虑把 $n$ 二进制拆分，即 $(n _{t}n _{t-1}\cdots n _{1}n _{0}) _{2}$，那么有：

$$n=n _{t}2 ^{t}+n _{t-1} 2 ^{t-1}+\cdots+n _{1}2 ^{1}+n _{0}2 ^{0}$$

带回 $a ^{n}$，可得：

$$\begin{aligned} a ^{n} & =\left(a ^{n _{t}2 ^{t}+\cdots+n _{0}2 ^{0}}\right)\\ & =a ^{n _{t}2 ^{t}}\times \cdots a ^{n _{0}2 ^{0}} \end{aligned}$$

由于我们可以 $O(1)$ 从 $2 ^{i}$ 项推到 $2 ^{i+1}$ 项，因此这个算法的复杂度为 $O(\log n)$。

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

等比数列求和

已知等比数列首项为 $a _{1}$，公比为 $q$，求前 $n$ 项和 $S _{n}\bmod p$ 的值。

Luogu2044 随机数生成器

这题有很多种方法，首先我们发现等比数列是个递推式，我们可以用矩阵加速递推来做。

其次假如 $p$ 是质数，我们可以直接借助等比数列求和公式：

$$S _{n}=\dfrac{a _{1}(1-q ^{n})}{1-q}$$

把 $1-q$ 的逆元求出来即可。

但当 $p$ 不是质数的时候呢？我们这时候无法保证一定能求出 $1-q$ 关于 $p$ 的逆元。

我们把式子列出来：

$$S _{n}=(a _{1}q+a _{1} q ^{2}+\cdots+a _{1}q ^{n-1})+a _{1}$$

我们令 $b _{n}=a _{1}q+a _{1}q ^{2}+\cdots +a _{1} q ^{n}$，那么 $S _{n}=a _{1}+b _{n-1}$。

我们对 $b _{n}$ 分类讨论，如果 $n$ 是偶数，那么我们提取公因数：

$$b _{n}=\left(a _{1}q+\cdots+a _{1}q ^{\frac{n}{2}}\right)\left(1+q ^{\frac{n}{2}}\right)=b _{\frac{n}{2}}\left(1+q ^{\frac{n}{2}}\right)$$

如果 $n$ 是奇数，我们把 $a _{n}$ 提出来，那么 $n-1$ 是偶数，可以继续处理：

$$b _{n}=b _{n-1}+a _{n}$$

这样同样可以在 $O(\log ^{2} n)$ 的时间里求出等比数列的前 $n$ 项和。

long long solve(long long a1,long long n,long long q)
{
    if(n==1)
        return a1;
    if(n&1)
        return (solve(a1,n-1,q)%p+a1*power(q,n-1,p))%p;
    else
        return (solve(a1,n/2,q)%p*(1+power(q,n/2,p))%p)%p;
}
long long sumn(long long a1,long long n,long long q)
{
    return a1+solve(a1,n-1,q);
}

发现在 solve 函数中，快速幂的指数和 solve 的 $n$ 相同，于是我们把 solve 的返回值改成 pair，即同时记录 $q ^{n}$，即可将时间复杂度优化至 $O(\log n)$。

pair solve(long long a1,long long n,long long q)
{
    if(n==1)
        return make_pair(a1,q);
    if(n&1)
    {
        pair<int,int> p=solve(a1,n-1,q);
        return make_pair(p.first%p+a1*p.second%p,p.second*q%p);
    }
    else
    {
        pair<int,int> p=solve(a1,n/2,q);
        return make_pair(p.first%p*(1+p.second)%p,p.second*p.second%p);
    }
}
long long sumn(long long a1,long long n,long long q)
{
    return a1+solve(a1,n-1,q).first;
}

矩阵快速幂

求矩阵快速幂。

Luogu3390 矩阵快速幂

关于矩阵知识请见线性代数总结。跟快速幂相似，把过程中的乘法换成矩阵乘法即可。

matrix multi(matrix a,matrix b)
{
    matrix c;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            c.a[i][j]=0;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            for(int k=1;k<=n;k++)
                c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]%mod)%mod;
    return c;
}
matrix power(matrix a,long long b)
{
    matrix ans=a;
    b--;
    for(;b;b>>=1)
    {
        if(b&1)
            ans=multi(ans,a);
        a=multi(a,a);
    }
    return ans;
}

---

质数

质数的定义很简单，即大于 $1$ 的自然数中，只被 $1$ 和它本身整除的数叫做质数，其它大于 $1$ 的自然数称作合数。我们用 $\pi(x)$ 代表小于等于 $x$ 的素数个数，随着 $x$ 增大，$\pi(x)\sim \dfrac{x}{\ln (x)}$。

质数判定

判断一个数是否是质数。

暴力做法

时间复杂度为 $O(\sqrt {n})$，即遍历 $1\sim \sqrt{n}$ 里的所有数，判断是否是 $n$ 的约数即可。

bool isprime(int n)
{
    if(n<2)
        return false;
    for(int i=2;i*i<=n;i++)
        if(n%i==0)
            return false;
    return true;
}

若已提前将 $1\sim \sqrt{n}$ 里的素数预先处理出来，那么只用判断这些素数即可。根据素数计数函数可得素数个数约为 $\dfrac{\sqrt {n}}{\log \sqrt{n}}$ 个，故时间复杂度为 $O\left(\dfrac{\sqrt{n}}{\log n}\right)$。

bool isprime(int n)
{
    if(n<2)
        return false;
    for(int i=1;i<=cnt;i++)
        if(n%prime[i]==0)
            return false;
    return true;
}

Miller Rabin 算法

HDU2138 How many prime numbers

时间复杂度为 $O(k\log ^{3} n)$，利用 FFT 可以降到 $O(k\log ^{2} n)$。

Miller Rabin 并不是一个确定性算法，它有一定的错误率，但能极大概率测试其是否为素数。

要学习 Miller Rabin 算法，首先要了解两个基础定理。

---

费马小定理

当 $p$ 为质数时，且 $\gcd(a,p)=1$ 时，有 $a ^{p-1}\equiv 1\pmod p$。但反过来不一定成立，即如果 $a,p$ 互质，且 $a ^{p-1}\equiv 1\pmod p$，不能推出 $p$ 是质数，例如 $\text{Carmichael}$ 数。

它的另一个形式是对于任意整数 $a$，有 $a ^{p}\equiv a\pmod p$。

二次探测定理

当 $p$ 为质数时，且 $0<x<p$，则方程 $x ^{2}\equiv 1\pmod p$ 的解为 $x=1$ 或 $x=p-1$。

---

理解这两个定理后，就可以开始判断质数了。首先对于偶数和 $0,1,2$ 我们可以直接判断是否为质数，否则设要测试的数为 $n$，我们取较小的质数 $a$，设 $s,t$ 满足 $2 ^{s} t=n-1$。然后算出 $a ^{t}$，不断平方 $s$ 次进行二次探测，最后根据费马小定理判断 $a ^{n-1}\bmod n$ 是否等于 $1$ 即可。多次取不同的质数 $a$ 可以使正确性更高。经过测试得，如果 $n<4759123141$，那么只需要测试三个底数 $2,7,61$ 即可；如果用前 $7$ 个质数测试，那么所有不超过 $341550071728320$ 的数都使正确的；如果用 $2,3,7,61,24251$ 作为底数，那么 $10 ^{16}$ 内唯一的强伪质数为 $46856248255981$；如果选取 $30$ 以内的所有质数，那么 int 范围内的数不会出错。

最后需要注意，由于要求 long long 类型的平方，可能会溢出，因此需要快速乘来进行乘法。

long long p[9]={0,2,3,7,61,24251};
inline long long multi(long long a,long long b,long long p)
{
    long long ans=0;
    for(;b;b>>=1)
    {
        if(b&1)
            ans=(ans+a)%p;
        a=(a+a)%p;
    }
    return ans%p;
}
inline long long power(long long a,long long b,long long p)
{
    long long ans=1%p;
    for(;b;b>>=1)
    {
        if(b&1)
            ans=multi(ans,a,p);
        a=multi(a,a,p);
    }
    return ans%p;
}
inline bool Miller_Rabin(long long x)
{
    if(x==2)
        return true;
    if(x<2||(!(x&1))||x==46856248255981)//特判
        return false;
    long long s=0,t=x-1;
    while(!(t&1))
    {
        t>>=1;
        s++;
    }//拆分s,t
    for(int i=1;i<=5&&p[i]<x;i++)
    {
        long long r=power(p[i],t,x);//判断每个质数
        for(int j=1;j<=s;j++)//平方s次
        {
            long long q=multi(r,r,x);
            if(q==1&&r!=1&&r!=x-1)//二次探测定理
                return false;
            r=q;
        }
        if(r!=1)//费马小定理
            return false;
    }
    return true;
}

筛法

筛法的作用其实是将 $1\sim n$ 之内所有的素数筛出来，当然也可以来判断质数(虽然有点大材小用)。

质数筛主要有两种：埃氏筛和线性筛(欧拉筛)

Eratosthenes 筛法

简称埃氏筛。时间复杂度为 $O(n\log \log n)$。就是对于每一个质数就把 $n$ 之内它的倍数标记。

for(int i=2;i<=n;i++)
{
    if(!vis[i])
    {
        prime[++cnt]=i;
        for(int j=i*2;j<=n;j+=i)
            vis[j]=1;
    }
}

Euler 筛法

Luogu3383 线性筛素数

中文名为欧拉筛，也称为线性筛。时间复杂度为 $O(n)$。可以发现在埃氏筛中有些数会被重复标记，例如 $12$ 既会被 $2$ 标记，也会被 $3$ 标记。而线性筛中则保证每个合数只会被它最小的质因数筛去。

for(int i=2;i<=n;i++)
{
    if(!vis[i])
        prime[++cnt]=i;
    for(int j=1;j<=cnt;j++)
    {
        if(i*prime[j]>n)
            break;
        vis[i*prime[j]]=1;
        if(i%prime[j]==0)//这一步保证每个合数只会被标记一次
            break;
    }
}

分解质因数

将一个数 $n$ 分解质因数。

暴力做法

时间复杂度为 $O(\sqrt{n})$，即对于 $2\sim \sqrt{n}$ 的每一个数判断是否能整除，能整除就整除。

void divide(int n)
{
    int q=sqrt(n);
    for(int i=2;i<=q;i++)
        if(n%i==0)
        {
            prime[++cnt]=i;
            while(n%i==0)
            {
                n/=i;
                c[cnt]++;
            }
        }
    if(n>1)
    {
        prime[++cnt]=n;
        c[cnt]=1;
    }
}

同理，若已经将 $1\sim \sqrt{n}$ 预处理好的话，那么时间复杂度为 $\dfrac{\sqrt{n}}{\log n}$。

void divide(int n)
{
    tot=0;
    for(int i=1;i<=cnt;i++)
        if(n%prime[i]==0)
        {
            p[++tot]=prime[i];
      		while(n%prime[i]==0)
            {
                n/=prime[i];
                c[tot]++;
            }
        }
    if(n>1)
    {
        p[++tot]=n;
        c[tot]=1;
    }
}

Pollard Rho 算法

Luogu4718 Pollard-Rho算法

Pollard Rho 是一个随机算法，能快速找到一个数 $n$ 的一个因数，时间复杂度约为 $O(\sqrt[4]{n})$。

我们尝试随机在 $[1,\sqrt{n}]$ 猜一个 $n$ 的约数，概率约为 $\dfrac{1}{\sqrt{n}}$，我们需要改进算法来提高我们的准确率。

---

生日悖论

假如说一个班上有 $k$ 个人，要找到一个人的生日是某月某日，概率是非常低的。但是要找两个生日相同的人，当 $k=23$ 时概率就超过 $50\%$ 了，而当 $k=60$ 时概率为 $0.9999$，极为接近 $1$。但是这似乎与我们的常识并不符合，因此叫做生日悖论。

---

生日悖论给了我们启示：当我们随机选择组合时，满足答案的组合会比满足答案的单体要多得多，因此可以用来提高正确率。 例如在 $\left[1,1000\right]$ 中取得 $42$ 的概率为 $\dfrac{1}{1000}$，但时随机取两个数 $i,j$ 使得 $|i-j|=42$ 的概率却大的多，依次类推，我们选择 $k$ 个数 $x _{1},x _{2},\cdots,x _{k}$ 满足 $|x _{i}-x _{j}|=42$，当 $k=30$ 时成功几率已经超过一半，$k=100$ 时几乎确定能成功。

我们回到原问题，我们同理选择 $k$ 个数，并询问是否存在 $x _{i}-x _{j}$ 能整除 $n$。当 $k=\sqrt{n}$ 时可能性会超过一半。但是对于一个大于 $10 ^{10}$ 的数，我们依然要做大量运算。我们换个思路，询问是否存在 $\gcd(x _{i}-x _{j},n)>1$，这样答案的可能性会大大增加，但我们依然需要选取大约 $\sqrt[4]{n}$ 个数，很明显我们不能随机选择，我们要有规律的生成。

Pollard Rho 算法只用到了两个数 $x,a$，并用到了一个神奇的函数 $f$ 来生成选取的数：

$$f(x)=(x ^{2}+a)\bmod n$$

我们令 $x _{1}=2$，$a$ 由 rand 随机生成开始，令 $x _{2}=f(x _{1}),x _{3}=f(x _{2})\cdots$，即生成规则为 $x _{n+1}=f(x _{n})$。其中 $x _{1},x _{2}\cdots x _{k}$ 即是我们要选取的数。在一定范围内，这个数列是随机的，我们可以取相邻两项做差求最大公约数。

但是这个数列也会进入死循环，因为它的生成数列的轨迹很香希腊字母 $\rho$，因此算法名字叫做 Pollard Rho。在模 $n$ 意义下，函数在迭代过程中总会代入之前的值得到同样结果。

该怎么判断这个环呢？这里有几个不同的方法。

首先就是暴力，拿个数组记录之前所有产生的数判断即可。

另一种方法是 Floyd 发明的判环算法，也被称为“龟兔赛跑”算法。假设我们在一个圆形轨道上行走，我们如何直到我们已经走完了一圈呢？我们设两个人，暂且叫它们 mich 和 clockwhite，我们让 clockwhite 按照 mich 的两倍速度从同一起点往前走，这是一个典型的追及问题，当 clockwhite 第一次追上 mich 时，我们就知道 clockwhite 已经走了至少一圈了。在本问题中我们就设 $s,t$，$s$ 每迭代 $1$ 次 $t$ 就迭代 $2$ 次。当第一次 $x _{s}=x _{t}$ 时，$t-s$ 便是环的长度。

inline long long f(long long x,long long a,long long n)
{
    return (__int128)(x*x+a)%n;
}
inline long long floyd(long long n)
{
    long long a=rand()%(n-1)+1;
    int s=f(2,a,n);
    int t=f(s,a,n);
    while(s!=t)
    {
        long long d=gcd(abs(t-s),n);
        if(d>1)
            return d;//找到约数
        s=f(s,a,n);
        t=f(f(t,a,n),a,n);
    }
    return -1;//没找到,重新调整参数a
}

第三种方法是基于路径倍增，同时还有一些常数优化。由于求 $\gcd$ 的复杂度为 $O(\log n)$，我们频繁调用这个函数会让算法变慢。显然，由于 $\gcd(a,n)>1$，那么对于任意一个正整数 $b$，满足 $\gcd(ab,n)>1$，即 $\gcd(ab\bmod n,n)>1$。因此我们可以把所有选出的测试的数在模 $n$ 意义下乘起来再最后做一次 $\gcd$ 即可。但是要选取多少个呢？我们可以用倍增的思想来做，这样每次累计的样本个数是 $2 ^{k}$ 个，$k$ 依次递增。如此我们可以在一定程度避免环的问题，同时也可以避免取 $\gcd$ 的次数过多。但是如果倍增次数过多的话，算法需要积累的样本就越多，我们需要规定一个上界，这是一个神奇的数——$127$，将它作为上界可以更优。

最后我们结合 Miller Rabin 算法，先判断这个数是否为质数，若不是则用 Pollard Rho 找因子 $p$，然后把 $n$ 中所有的 $p$ 除掉，递归分解时同时记录最大值即可。

以下是基于路径倍增的代码。

/*
 * @Author: clorf 
 * @Date: 2020-09-21 10:19:30 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-21 11:51:39
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=100010;
const int mod=1e9+7;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
int t;
long long n,pr[6]={0,2,3,7,61,24251},maxx;
inline long long power(long long a,long long b,long long p)
{
    long long ans=1;
    for(;b;b>>=1)
    {
        if(b&1)
            ans=(__int128)ans*a%p;
        a=(__int128)a*a%p;
    }
    return ans%p;
}
inline bool Miller_Rabin(long long x)
{
    if(x==2)
        return true;
    if(x<2||x==46856248255981||(!(x&1)))
        return false;
    long long k=x-1,s=0;
    while(!(k&1))
    {
        k>>=1;
        s++;
    }
    for(int i=1;i<=5&&pr[i]<x;i++)
    {
        long long r=power(pr[i],k,x);
        for(int j=1;j<=s;j++)
        {
            long long q=(__int128)r*r%x;
            if(q==1&&r!=1&&r!=x-1)
                return false;
            r=q;
        }
        if(r!=1)
            return false;
    }
    return true;
}
long long gcd(long long a,long long b)
{
    if(!b)
        return a;
    return gcd(b,a%b);
}
inline long long f(long long x,long long c,long long p)
{
    return ((__int128)x*x+c)%p;
}
inline long long Pollard_Rho(long long x)
{
    long long s=0,t=0;
    long long c=1ll*rand()%(x-1)+1,val=1;
    for(int k=1;;k<<=1,s=t,val=1)
    {
        for(int step=1;step<=k;step++)
        {
            t=f(t,c,x);
            val=(__int128)val*abs(t-s)%x;
            if((step%127)==0)
            {
                long long d=gcd(val,x);
                if(d>1)
                    return d;
            }
        }
        long long d=gcd(val,x);
        if(d>1)
            return d;
    }
}
void solve(long long x)
{
    if(x<=maxx||x<2)
        return ;
    if(Miller_Rabin(x))
    {
        maxx=max(maxx,x);
        return ;
    }
    long long p=x;
    while(p>=x)//while防止PR(x)=x
        p=Pollard_Rho(x);
    while(!(x%p))
        x/=p;
    solve(p);
    solve(x);
}
int main()
{
    scanf("%d",&t);
    while(t--)
    {
        srand((unsigned)time(NULL));
        scanf("%lld",&n);
        maxx=0;
        solve(n);
        if(maxx==n)
            printf("Prime\n");
        else
            printf("%lld\n",maxx);
    }
    return 0;
}

阶乘质因数分解

Luogu2043 质因数分解

对 $n!$ 进行质因数分解。

若我们采用暴力做法，时间复杂度为 $O(n\sqrt{n})$，当 $n=400000$ 时就不能承受了。

我们考虑更优秀的算法。首先用线性筛筛出 $n$ 以内所有的质数，由于 $n!=1\times 2\times \cdots \times n$，所以对于 $n$ 以内的质数 $p$，它的出现次数为：

$$\left\lfloor\dfrac{n}{p}\right\rfloor+\left\lfloor\dfrac{n}{p ^{2}}\right\rfloor+\cdots +\left\lfloor\dfrac{n}{p ^{\log _{p} n}}\right\rfloor$$

为什么？首先 $p$ 的倍数有 $1$ 个质因子 $p$，总共有 $\left\lfloor\dfrac{n}{p}\right\rfloor$ 个 $p$ 的倍数；$p ^{2}$ 的倍数有 $2$ 个质因子 $p$，由于 $p ^{2}$ 的倍数一定是 $p$ 的倍数，所以先前已经计入过 $1$ 个因子 $p$，只需要再计入 $1$ 次即可，即 $\left\lfloor\dfrac{n}{p ^{2}}\right\rfloor$ 个因子，$p$ 的更高次方同理。

时间复杂度小于 $n\log \log n$，类比埃氏筛即可。

/*
 * @Author: clorf 
 * @Date: 2020-09-21 14:35:35 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-21 16:07:17
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=10010;
const int mod=1e9+7;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
int n,pr[maxn],cnt,c[maxn];
bool vis[maxn];
int main()
{
    scanf("%d",&n);
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])
        {
            pr[++cnt]=i;
            vis[i]=1;
        }
        for(int j=1;j<=cnt;j++)
        {
            if(i*pr[j]>n)
                break;
            vis[i*pr[j]]=1;
            if(i%pr[j]==0)
                break;
        }
    }
    for(int i=1;i<=cnt;i++)
    {
        int t=n;
        while(t)
        {
            t/=pr[i];
            c[i]+=t;
        }
    }
    for(int i=1;i<=cnt;i++)
        printf("%d %d\n",pr[i],c[i]);
    return 0;
}

欧拉函数

欧拉函数是基础数论函数之一，用 $\varphi$ 表示。

定义

$\varphi(n)$ 表示 $1\sim n$ 中与 $n$ 互质的数的数量，假设 $\displaystyle{n=\prod _{i=1} ^{m}p _{i} ^{c _{i}}}$，那么 $\displaystyle{\varphi(n)=n\times \prod _{i=1} ^{m}\left(1-\dfrac{1}{p _{i}}\right)}$。

求一个数的欧拉函数值

直接根据定义质因数分解求即可，可以用 Pollard Rho 算法优化。

inline long long eular(long long x)
{
    long long ans=x;
    long long p=x;
    for(int i=2;i*i<=x;i++)
        if(p%i==0)
        {
            ans=ans/i*(i-1);
            while(p%i==0)
                p/=i;
        }
    if(p>1)
        ans=ans/p*(p-1);
    return ans;
}

求多个数的欧拉函数值

求 $2\sim n$ 的欧拉函数值。用筛法来做，如果用埃氏筛，复杂度为 $O(n\log \log n)$。

inline void eular(long long n)
{
    for(int i=2;i<=n;i++)
        phi[i]=i;
    for(int i=2;i<=n;i++)
        if(phi[i]==i)
            for(int j=i;j<=n;j+=i)
                phi[j]=phi[j]/i*(i-1);
}

利用线性筛和下面的最后一个性质，可以在 $O(n)$ 内求出欧拉函数值。

inline void eular(long long n)
{
    for(int i=2;i<=n;i++)
	{
		if(!vis[i])
		{
			vis[i]=1;
			pr[++cnt]=i;
			phi[i]=i-1;
		}
		for(int j=1;j<=cnt;j++)
		{
			if(pr[j]*i>n)
				break;
			vis[i*pr[j]]=pr[j];
			phi[i*pr[j]]=phi[i]*(i%pr[j]?pr[j]-1:pr[j]);
            if(i%pr[j]==0)
                break;
		}
	}
}

性质

• 当 $n$ 是质数时，$\varphi(n)=n-1$；

• $\varphi(n)$ 是积性函数，即 $\gcd(a,b)=1$，那么 $\varphi(ab)=\varphi(a)\times \varphi(b)$；

• $n=\sum _{d|n} \varphi(d)$；

• 若 $n=p ^{k}$，$p$ 是质数，那么 $\varphi(n)=p ^{k}-p ^{k-1}$；

• $\forall n>1$，$1\sim n$ 中与 $n$ 互质数的和为 $\dfrac{n\times \varphi(n)}{2}$；

• 若 $p$ 是质数且 $p|n$，若 $p ^{2}|n$，则 $\varphi(n)=\varphi\left(\dfrac{n}{p}\right)\times p$；若 $p ^{2}\nmid n$，则 $\varphi(n)=\varphi\left(\dfrac{n}{p}\right)\times (p-1)$。

欧拉定理

若 $\gcd(a,n)=1$，则 $a ^{\varphi(n)}\equiv 1\pmod n$。当 $n$ 是质数时，由于 $\varphi(n)=n-1$，代入欧拉定理即可得到费马小定理。

扩展欧拉定理

扩展欧拉定理的主要作用是用于降幂。

$$a ^{b}=\begin{cases} a ^{b\bmod \varphi(p)}, & \gcd(a,p)=1\\ a ^{b}, & \gcd(a,p)\ne 1,b<\varphi(p)\\ a ^{b\bmod \varphi(p)+\varphi(p)}, & \gcd(a,p)\neq 1,b\ge \varphi(p) \end{cases} \pmod p$$

其中当 $\gcd(a,p)=1$ 时，也可以用 $a ^{b\bmod \varphi(p)+\varphi(p)}$ 计算。

Luogu5091 扩展欧拉定理

/*
 * @Author: clorf 
 * @Date: 2020-09-21 16:17:55 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-21 16:33:56
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=20000010;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
long long a,b,m,t,mod;
char s[maxn];
inline long long eular(long long x)
{
    long long ans=x;
    long long p=x;
    for(int i=2;i*i<=x;i++)
        if(p%i==0)
        {
            ans=ans/i*(i-1);
            while(p%i==0)
                p/=i;
        }
    if(p>1)
        ans=ans/p*(p-1);
    return ans;
}
inline long long power(long long a,long long b,long long p)
{
    if(b<=0)
        return 1;
    long long ans=1;
    for(;b;b>>=1)
    {
        if(b&1)
            ans=ans*a%p;
        a=a*a%p;
    }
    return ans%p;
}
int main()
{
    scanf("%lld%lld",&a,&m);
    mod=eular(m);
    bool flag=0;
    scanf("%s",s+1);
    int len=strlen(s+1);
    for(int i=1;i<=len;i++)
    {
        b=10*b+s[i]-'0';
        if(b>=mod)
        {
            b%=mod;
            flag=1;
        }
    }
    if(!flag)
        printf("%lld",power(a,b,m));
    else
        printf("%lld",power(a,b%mod+mod,m));
    return 0;
}

最大公约数

最大公约数即 $\text{Greatest Common Divisor}$，缩写为 $\gcd$。一组数的最大公约数是指它们的公约数中最大的那个。

求最大公约数

欧几里得算法

已知两个数 $a,b$，求 $\gcd(a,b)$。

欧几里得算法可以在 $O(\log \max(a,b))$ 的时间内求出 $\gcd$，由于它的实现方式是辗转相除，所以也被叫做辗转相除法。辗转相除法主要用到了下面这个式子：

$$\gcd(a,b)=\gcd(b,a\bmod b)$$

然后就可以递归处理了。

int gcd(int a,int b)
{
    if(!b)
        return a;
    return gcd(b,a%b);
}

更相减损法

其实不太常用。出自《九章算术》。式子如下：

$$\gcd(a,b)=\gcd(b,a-b)$$

int gcd(int a,int b)
{
    if(a==b)
        return a;
    else if(a>b)
        a-=b;
    else
        b-=a;
    return gcd(a,b);
}

二进制算法

避免了欧几里得算法的大量取模操作，可以有效减少时间消耗。

对于两个数 $a,b$，我们进行如下操作：

• 若 $a,b$ 都为偶数，那么 $\gcd(a,b)=2\gcd\left(\dfrac{a}{2},\dfrac{b}{2}\right)$；
• 若 $a$ 为偶数，$b$ 为奇数，那么 $\gcd(a,b)=\gcd\left(\dfrac{a}{2},b\right)$；
• 若 $a$ 为奇数，$b$ 为偶数，那么 $\gcd(a,b)=\gcd\left(a,\dfrac{b}{2}\right)$；
• 若 $a,b$ 都为奇数，那么 $\gcd(a,b)=\gcd(b,a-b)$。

递归即可。

int gcd(int a,int b)
{
    if(a<b)
    	return gcd(b,a);
    if((a&1)&&(b&1))
        return gcd(b,a-b);
    else if((a&1)&&(!(b&1)))
        return gcd(a,b/2);
    else if((!(a&1))&&(b&1))
        return gcd(a/2,b);
    else
        return gcd(a/2,b/2)*2;
}

最小公倍数

最小公倍数即 $\text{Least Common Multiple}$，缩写为 $\operatorname{lcm}$。

求最小公倍数先要求最大公约数，因为满足这样一个定理：

$$a\times b=\gcd(a,b)\times \operatorname{lcm}(a,b)$$

于是求出 $\gcd(a,b)$ 即可。

inline int lcm(int a,int b)
{
    return a/gcd(a,b)*b;
}

裴蜀定理

又称贝祖定理($\text{Bézout's lemma}$)。若 $a,b$ 是不全为 $0$ 的整数，则存在整数 $x,y$，使得 $ax+by=\gcd(a,b)$。它的另一个形式是关于 $x,y$ 的不定方程 $ax+by=c$ 有整数解的充要条件是 $\gcd(a,b)|c$。

扩展到多个数的情况，对于 $\sum _{i=1} ^{n}a _{i}x _{i}=c$ 有整数解的充要条件是 $\gcd(a _{1},a _{2},\cdots,a _{n})|c$。

Luogu4549 裴蜀定理

/*
 * @Author: clorf 
 * @Date: 2020-09-21 19:47:13 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-21 20:00:56
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=100010;
const int mod=1e9+7;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
int n,ans;
int gcd(int a,int b)
{
    if(!b)
        return a;
    return gcd(b,a%b);
}
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
    {
        int x;
        scanf("%d",&x);
        ans=gcd(ans,abs(x));
    }
    printf("%d",ans);
    return 0;
}

Luogu4571 瓶子和燃料

一道需要运用裴蜀定理的题，就是要从 $n$ 个数中找 $k$ 个数使得其最大公约数最大。我们对这 $n$ 个数分解因数(注意不是质因数！！)，然后找分解出的因数中个数大于 $k$ 的数，取最大值即可。注意为了防止 MLE 要用 map 存储。

/*
 * @Author: clorf 
 * @Date: 2020-09-29 10:00:29 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-29 10:11:33
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#include<unordered_map>
#define INF 1e9
using namespace std;
const int maxn=1010,maxv=10000010;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
long long n,k;
long long v[maxn],cnt,ans;
unordered_map<int, long long> fac, p;
int main()
{
    scanf("%lld%lld",&n,&k);
    for(int i=1;i<=n;i++)   
        scanf("%lld",&v[i]);
    for(int i=1;i<=n;i++)
        for(int j=1;j*j<=v[i];j++)
            if(v[i]%j==0)
            { 
                fac[++cnt]=j;
                p[j]++;
                int r=v[i]/j;
                if(j!=r)
                {
                    fac[++cnt]=r;
                    p[r]++;
                }
            }
    for(int i=1;i<=cnt;i++)
        if(p[fac[i]]>=k)
            ans=max(ans,fac[i]);
    printf("%lld",ans);
    return 0;
}

扩展欧几里得

扩展欧几里得算法用于求 $ax+by=\gcd(a,b)$ 的一组特解 $x,y$。

因为 $\gcd(a,b)=\gcd(b,a\bmod b)$，所以我们可以列出以下方程：

$$\begin{cases} ax _{1}+by _{1}=\gcd(a,b)\\ bx _{2}+(a\bmod b)y _{2}=\gcd(b,a\bmod b) \end{cases}$$

于是可以得到 $ax _{1}+by _{1}=bx _{2}+(a\bmod b)y _{2}$，又因为 $a\bmod b=a-\left(\left\lfloor\dfrac{a}{b}\right\rfloor\times b\right)$，代入化简后即可得到 $ax _{1}+by _{1}=ay _{2}+b\left(x _{2}-\left\lfloor\dfrac{a}{b}\right\rfloor y _{2}\right)$。于是得到 $x _{1}=y _{2},y _{1}=x _{2}-\left\lfloor\dfrac{a}{b}\right\rfloor y _{2}$。

于是可以递归求解。边界条件是 $b=0$ 的时候，这时候让 $x=1,y=0$ 即可。

inline int exgcd(int a,int b,int &x,int &y)
{
    if(!b)
    {
        x=1;
        y=0;
        return a;
    }
    int ans=exgcd(b,a%b,x,y);
    int t=x;
    x=y;
    y=t-a/b*y;
    return ans;
}

得到了一组特解 $x _{0},y _{0}$ 后，现在我们考虑求出通解 $x,y$。通解如下：
$$
\begin{aligned}
x & =x _{0}+k\dfrac{b}{\gcd(a,b)}\
y & =y _{0}-k\dfrac{a}{\gcd(a,b)}
\end

$$

然后我们就可以解普通的二元一次方程 $ax+by=c$ 了。

[Luogu5656 二元一次不定方程(exgcd)](https://www.luogu.com.cn/problem/P5656)

首先如果 $\gcd(a,b)\nmid c$，那么由于裴蜀定理，很明显无解。接着将用 ``exgcd`` 求出的 $ax _{0}+by _{0}=\gcd(a,b)$ 的特解 $x _{0},y _{0}$ 乘以 $\dfrac{c}{\gcd(a,b)}$ 即可。通解同理：
$$
\begin{aligned}
x & =\dfrac{c}{\gcd(a,b)}x _{0}+k\dfrac{b}{\gcd(a,b)}\\
y & =\dfrac{c}{\gcd(a,b)}y _{0}-k\dfrac{a}{\gcd(a,b)}
\end{aligned}
~~~~ ~ (k\in \mathbb{Z})
$$

最小正整数解 $x _{1},y _{1}$ 的式子如下：
$$
\begin{aligned}
x _{1} & =(x _{0}\bmod \dfrac{b}{\gcd(a,b)}+\dfrac{b}{\gcd(a,b)})\bmod \dfrac{b}{\gcd(a,b)}\\
y _{1} & =(y _{0}\bmod \dfrac{a}{\gcd(a,b)}+\dfrac{a}{\gcd(a,b)})\bmod \dfrac{a}{\gcd(a,b)}
\end{aligned}
$$

同时，当 $x$ 取最小正整数解 $x _{1}$ 时，得到的 $y$ 是最大整数解；$y$ 取得 $y _{1}$ 时同理。

/*
 * @Author: clorf 
 * @Date: 2020-09-21 20:01:22 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-21 22:49:50
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=100010;
const int mod=1e9+7;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
long long t,a,b,c;
long long exgcd(long long a,long long b,long long &x,long long &y)
{
    if(!b)
    {
        x=1;
        y=0;
        return a;
    }
    long long ans=exgcd(b,a%b,x,y);
    long long t=x;
    x=y;
    y=t-a/b*y;
    return ans;
}
int main()
{
    scanf("%lld",&t);
    while(t--)
    {
        long long x,y;
        scanf("%lld%lld%lld",&a,&b,&c);
        long long num=exgcd(a,b,x,y);
        if(c%num!=0)
        {
            printf("-1\n");
            continue;
        }
        long long k=c/num;
        long long x0=x*k,y0=y*k;
        long long p=a/num,q=b/num;
        long long miny=(y0%p+p)%p;
        long long minx=(x0%q+q)%q;
        if(!minx)
            minx+=q;
        if(!miny)
            miny+=p;
        long long maxy=(c-a*minx)/b;
        long long maxx=(c-b*miny)/a;
        if(maxy<=0)
            printf("%lld %lld\n",minx,miny);
        else
        {
            long long num1=(maxy-miny)/p+1;
            long long num2=(maxx-minx)/q+1;
            long long ans=min(num1,num2);
            printf("%lld %lld %lld %lld %lld\n",ans,minx,miny,maxx,maxy);
        }
    }   
    return 0;
}

# 类欧几里得

给定 $n,a,b,c$，分别求 $\displaystyle{\sum _{i=0} ^{n}\lfloor\dfrac{ai+b}{c}\rfloor,\sum _{i=0} ^{n}\lfloor\dfrac{ai+b}{c}\rfloor ^{2},\sum _{i=0} ^{n}i\lfloor\dfrac{ai+b}{c}\rfloor}$。

首先令：
$$
\begin{aligned}
f\left(a,b,c,n\right) & =\sum _{i=0} ^{n} \left\lfloor \dfrac{ai+b}{c} \right\rfloor\\
h\left(a,b,c,n\right) & =\sum _{i=0} ^{n} \left\lfloor \dfrac{ai+b}{c} \right\rfloor ^{2}\\
g\left(a,b,c,n\right) & =\sum _{i=0} ^{n}i \left\lfloor \dfrac{ai+b}{c} \right\rfloor 
\end{aligned}
$$
然后我们分情况开始推：

对于 $f$：
$$
\begin{aligned}
f\left(a,b,c,n\right) & =\sum _{i=0} ^{n}\left\lfloor \dfrac{ai+b}{c}\right\rfloor\\
 & = \begin{cases}
 x _{1}, & n=0\\
 x _{2}, & a=0\\
 x _{3}, & a\ge c~\text{or}~b\ge c\\
 x _{4}, & a\le c~\text{and}~a\le c\\
 \end{cases}
\end{aligned}
$$
对于 $x _{1}$：
$$
x _{1}=\left\lfloor \dfrac{b}{c}\right\rfloor
$$
对于 $x _{2}$：
$$
x _{2}=\left(n+1\right)\left\lfloor\dfrac{b}{c}\right\rfloor
$$
对于 $x _{3}$：
$$
\begin{aligned}
x _{3} & =\sum _{i=0} ^{n}\left\lfloor\dfrac{\left(\left\lfloor\dfrac{a}{c}\right\rfloor c+a\bmod c\right)i+\left(\left\lfloor\dfrac{b}{c}\right\rfloor c+b\bmod c\right)}{c}\right\rfloor\\
& =\sum _{i=0} ^{n}\left(\left\lfloor\dfrac{a}{c}\right\rfloor i+\left\lfloor\dfrac{b}{c}\right\rfloor +\left\lfloor\dfrac{\left(a\bmod c\right)i+\left(b\bmod c\right)}{c}\right\rfloor\right)\\
& =\dfrac{n\left(n+1\right)}{2}\left\lfloor\dfrac{a}{c}\right\rfloor+n\left\lfloor\dfrac{b}{c}\right\rfloor+\sum _{i=0} ^{n}\left\lfloor\dfrac{\left(a\bmod c\right)i+\left(b\bmod c\right)}{c}\right\rfloor\\
& =\dfrac{n\left(n+1\right)}{2}\left\lfloor\dfrac{a}{c}\right\rfloor+n\left\lfloor\dfrac{b}{c}\right\rfloor+f\left(a\bmod c,b\bmod c,c,n\right)
\end{aligned}
$$
对于 $x _{4}$：
$$
\begin{aligned}
m & =\left\lfloor\dfrac{an+b}{c}\right\rfloor-1 \\
x _{4} & =\sum _{j=1} ^{m+1}\sum _{i=0} ^{n}[j\le \left\lfloor\dfrac{ai+b}{c}\right\rfloor]\\
& =\sum _{j=0} ^{m}\sum _{i=0} ^{n}[j+1\le \dfrac{ai+b}{c}]\\
& =\sum _{j=0} ^{m}\sum _{i=0} ^{n}[ai> cj+c-b-1]\\
& =\sum _{j=0} ^{m}\sum _{i=0} ^{n}[i> \left\lfloor\dfrac{cj+c-b-1}{a}\right\rfloor]\\
& =\sum _{j=0} ^{m}\left(n-\left\lfloor\dfrac{cj+c-b-1}{a}\right\rfloor\right)\\
& =n\left(m+1\right)-f\left(c,c-b-1,a,m\right)
\end{aligned}
$$

对于 $h$：
$$
\begin{aligned}
h\left(a,b,c,n\right) & =\sum _{i=0} ^{n}\left\lfloor \dfrac{ai+b}{c}\right\rfloor ^{2}\\
 & = \begin{cases}
 y _{1}, & n=0\\
 y _{2}, & a=0\\
 y _{3}, & a\ge c~\text{or}~b\ge c\\
 y _{4}, & a\le c~\text{and}~a\le c\\
 \end{cases}
\end{aligned}
$$
对于 $y _{1}$：
$$
y _{1}=\left\lfloor \dfrac{b}{c}\right\rfloor ^{2}
$$
对于 $y _{2}$：
$$
y _{2}=\left(n+1\right)\left\lfloor\dfrac{b}{c}\right\rfloor ^{2}
$$
对于 $y _{3}$：
$$
\begin{aligned}
y _{3} & =\sum _{i=0} ^{n}\left(\left\lfloor\dfrac{a}{c}\right\rfloor i+\left\lfloor\dfrac{b}{c}\right\rfloor +\left\lfloor\dfrac{\left(a\bmod c\right)i+\left(b\bmod c\right)}{c}\right\rfloor\right) ^{2}\\
& =\sum _{i=0} ^{n}\bigg(\left\lfloor\dfrac{a}{c}\right\rfloor ^{2} i ^{2}+\left\lfloor\dfrac{b}{c}\right\rfloor ^{2}+\left\lfloor\dfrac{\left(a\bmod c\right)i+\left(b\bmod c\right)}{c}\right\rfloor ^{2} +2\left\lfloor\dfrac{a}{c}\right\rfloor \left\lfloor\dfrac{b}{c}\right\rfloor i\\
 &+2\left\lfloor\dfrac{a}{c}\right\rfloor i \left\lfloor\dfrac{\left(a\bmod c\right)i+\left(b\bmod c\right)}{c}\right\rfloor+2\left\lfloor\dfrac{b}{c}\right\rfloor \left\lfloor\dfrac{\left(a\bmod c\right)i+\left(b\bmod c\right)}{c}\right\rfloor\bigg) \\
& =\dfrac{n\left(n+1\right)\left(2n+1\right)}{6} \left\lfloor\dfrac{a}{c}\right\rfloor ^{2}+\left(n+1\right)\left\lfloor\dfrac{b}{c}\right\rfloor ^{2}+n\left(n+1\right)\left\lfloor\dfrac{a}{c}\right\rfloor \left\lfloor\dfrac{b}{c}\right\rfloor \\
& +h\left(a\bmod c,b\bmod c,c,n\right)+2\left\lfloor\dfrac{a}{c}\right\rfloor g\left(a\bmod c,b\bmod c,c,n\right)\\
& +2\left\lfloor\dfrac{b}{c}\right\rfloor f\left(a\bmod c,b\bmod c,c,n\right)
\end{aligned}
$$
对于 $y _{4}$：
$$
\begin{aligned}
m & =\left\lfloor\dfrac{an+b}{c}\right\rfloor-1 \\
y _{4} & =\sum _{i=0} ^{n} \left(\sum _{j=1} ^{m+1}[j\le \left\lfloor\dfrac{ai+b}{c}\right\rfloor] \right) ^{2}\\
& =\sum _{i=0} ^{n}  \left(\sum _{j=1} ^{m+1}[j\le \left\lfloor\dfrac{ai+b}{c}\right\rfloor] \right)\times  \left(\sum _{k=1} ^{m+1}[k\le \left\lfloor\dfrac{ai+b}{c}\right\rfloor] \right)\\
& =\sum _{j=0} ^{m}\sum _{k=0} ^{m}\sum _{i=0} ^{n}[i> \max\left(\left\lfloor\dfrac{cj+c-b-1}{a}\right\rfloor,\left\lfloor\dfrac{ck+c-b-1}{a}\right\rfloor\right)]\\
& = \sum _{j=0} ^{m}\sum _{k=0} ^{m}\left( n-\max\left(\left\lfloor\dfrac{cj+c-b-1}{a}\right\rfloor,\left\lfloor\dfrac{ck+c-b-1}{a}\right\rfloor\right) \right)\\
& = n(m+1) ^{2}-\sum _{j=0} ^{m}\sum _{k=0} ^{m}\max\left(\left\lfloor\dfrac{cj+c-b-1}{a}\right\rfloor,\left\lfloor\dfrac{ck+c-b-1}{a}\right\rfloor\right) \\
& = n(m+1) ^{2}-\left( 2 \sum _{j=0} ^{m}j\left\lfloor\dfrac{cj+c-b-1}{a}\right\rfloor+\sum _{j=0} ^{m}\left\lfloor\dfrac{cj+c-b-1}{a}\right\rfloor \right)\\
& =n(m+1) ^{2}-2g(c,c-b-1,a,m)-f(c,c-b-1,a,m)
\end{aligned}
$$
对于 $g$：
$$
\begin{aligned}
g\left(a,b,c,n\right) & =\sum _{i=0} ^{n}i\left\lfloor \dfrac{ai+b}{c}\right\rfloor\\
 & = \begin{cases}
 z _{1}, & n=0\\
 z _{2}, & a=0\\
 z _{3}, & a\ge c~\text{or}~b\ge c\\
z _{4}, & a\le c~\text{and}~a\le c\\
 \end{cases}
\end{aligned}
$$
对于 $z _{1}$：
$$
\begin{aligned}
z _{1} & =0
\end{aligned}
$$
对于 $z _{2}$：
$$
\begin{aligned}
z _{2} & =\left\lfloor\dfrac{b}{c}\right\rfloor\cdot\dfrac{n(n+1)}{2}
\end{aligned}
$$
对于 $z _{3}$：
$$
\begin{aligned}
z _{3} & =\sum _{i=0} ^{n}i\left\lfloor \dfrac{ai+b}{c}\right\rfloor \\
& =\sum _{i=0} ^{n}i\left(\left\lfloor\dfrac{a}{c}\right\rfloor i+\left\lfloor\dfrac{b}{c}\right\rfloor +\left\lfloor\dfrac{\left(a\bmod c\right)i+\left(b\bmod c\right)}{c}\right\rfloor\right) \\
& =\sum _{i=0} ^{n}\left\lfloor\dfrac{a}{c}\right\rfloor i ^{2}+\left\lfloor\dfrac{b}{c}\right\rfloor i+i\left\lfloor\dfrac{\left(a\bmod c\right)i+\left(b\bmod c\right)}{c}\right\rfloor \\
& =\left\lfloor\dfrac{a}{c}\right\rfloor \cdot \dfrac{n(n+1)(2n+1)}{6}+\left\lfloor\dfrac{b}{c}\right\rfloor\cdot \dfrac{n(n+1)}{2}+g(a\bmod c,b\bmod c,c,n)
\end{aligned}
$$
对于 $z _{4}$：
$$
\begin{aligned}
m & =\left\lfloor\dfrac{an+b}{c}\right\rfloor-1 \\
z _{4} & =\sum _{j=0} ^{m}\sum _{i=0} ^{n}i[i>\left\lfloor\dfrac{jc+c-b-1}{a}\right\rfloor]\\
& =\sum _{j=0} ^{m}\sum _{i=\lfloor\frac{jc+c-b-1}{a}\rfloor+1} ^{n}i\\
& =\sum _{j=0} ^{m}\left(\sum _{i=0} ^{n} i-\sum _{i=0} ^{\lfloor\frac{jc+c-b-1}{a}\rfloor}i\right)\\
& =\sum _{j=0} ^{m}\dfrac{n(n+1)}{2}-\dfrac{\lfloor\frac{jc+c-b-1}{a}\rfloor(\lfloor\frac{jc+c-b-1}{a}\rfloor+1)}{2}\\
& =\dfrac{1}{2}\left(n(n+1)(m+1)-\sum _{j=0} ^{m}\left\lfloor\dfrac{jc+c-b-1}{a}\right\rfloor ^{2}-\sum _{j=0} ^{m} \left\lfloor\dfrac{jc+c-b-1}{a}\right\rfloor\right)\\
& =\dfrac{1}{2}\left(n(n+1)(m+1)-h(c,c-b-1,a,m)-f(c,c-b-1,a,m)\right)
\end{aligned}
$$
由此我们得到了 $f,h,g$ 的求解式，发现三个函数在不同情况时用到的递归式的参数也相同，于是一起递归求解即可。

[Luogu5170 类欧几里得算法](https://www.luogu.com.cn/problem/P5170)

/*
 * @Author: clorf 
 * @Date: 2020-09-22 09:18:16 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-22 10:26:48
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=100010;
const int mod=998244353;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
struct euclid
{
    long long f;
    long long h;
    long long g;
};
long long t;
const long long inv2=499122177ll;
const long long inv6=166374059ll;
inline long long sum1(long long x)
{
    return x*(x+1)/2%mod;
}
inline long long sum2(long long x)
{
    return x*(x+1)%mod*(x+x+1)%mod*inv6%mod;
}
euclid solve(long long a,long long b,long long c,long long n)
{
    // printf("%lld %lld %lld %lld\n",a,b,c,n);
    long long x=0,y=0,z=0;
    euclid ans;
    ans.f=ans.g=ans.h=0;
    if(!n)
    {
        long long q=b/c;
        q%=mod;
        x=q;
        y=q*q%mod;
        ans=(euclid){x%mod,y%mod,0};
        return ans;
    }
    if(!a)
    {
        long long u=sum1(n);
        long long q=b/c;
        q%=mod;
        x=(n+1)*q%mod;
        y=(n+1)*q%mod*q%mod;
        z=q*u%mod;
        ans=(euclid){x%mod,y%mod,z%mod};
        return ans;
    }
    if((a>=c)||(b>=c))
    {
        long long r=a/c;
        r%=mod;
        long long q=b/c;
        q%=mod;
        euclid p=solve(a%c,b%c,c,n);
        long long u=sum1(n);
        long long v=sum2(n);
        x=r*u%mod;
        x=(x+(n+1)*q%mod)%mod;
        x=(x+p.f)%mod;
        y=r*r%mod*v%mod;
        y=(y+(n+1)*q%mod*q%mod)%mod;
        y=(y+p.h)%mod;
        y=(y+2*r%mod*q%mod*u%mod)%mod;
        y=(y+2*q%mod*p.f)%mod;
        y=(y+2*r%mod*p.g)%mod;
        z=r*v%mod;
        z=(z+q*u%mod)%mod;
        z=(z+p.g)%mod;
        ans=(euclid){x%mod,y%mod,z%mod};
        return ans;
    }
    long long m=(a*n+b)/c-1;
    euclid p=solve(c,c-b-1,a,m);
    x=n*(m+1)%mod;
    x=(x-p.f+mod)%mod;
    y=n*(m+1)%mod*(m+1)%mod;
    long long j=(2*p.g%mod+p.f)%mod;
    y=(y-j+mod)%mod;
    z=n*(n+1)%mod*(m+1)%mod;
    long long i=(p.h+p.f)%mod;
    z=(z-i+mod)%mod;
    z=z*inv2%mod;
    ans=(euclid){x%mod,y%mod,z%mod};
    return ans;
}
int main()
{
    scanf("%lld",&t);
    while(t--)
    {
        long long n,a,b,c;
        scanf("%lld%lld%lld%lld",&n,&a,&b,&c);
        euclid ans=solve(a,b,c,n);
        printf("%lld %lld %lld\n",ans.f%mod,ans.h%mod,ans.g%mod);
    }
    return 0;
}

[BZOJ2712 棒球](https://darkbzoj.tk/problem/2712)

给定分数 $\dfrac{p}{q}$ 四舍五入到第 $n$ 位的值，求 $q$ 的最小值。

我们的难点就是在于把小数转成 $q$ 最小的分数，即求出一个范围，设出 $a,b,c,d$ 使得 $\dfrac{a}{b}\le \dfrac{p}{q}<\dfrac{c}{d}$，且要让 $q$ 最小。这下就转为了求两个分数之间分母最小的分数。

首先当 $\dfrac{a}{b}\sim \dfrac{c}{d}$ 之间存在整数时，我们显然取最小的整数。当 $a=0$ 时，条件只有 $\dfrac{p}{q}<\dfrac{c}{d}$，即 $q>\dfrac{dp}{c}$，当 $p=1$ 时显然 $q$ 取得最小值 $\left\lfloor\dfrac{d}{c}\right\rfloor+1$。这些是边界条件。

然后分类讨论，当 $a<b$ 时，我们将分数取倒数，得到 $\dfrac{d}{c}<\dfrac{q}{p}\le \dfrac{b}{a}$，这样就转化为了下面的情况。当 $a\ge b$ 时，我们只保留小数部分，令 $t=\left\lfloor\dfrac{a}{b}\right\rfloor$，那么 $\dfrac{a\bmod b}{b}\le \dfrac{p-qt}{q}<\dfrac{c-dt}{d}$，继续递归即可。

/*
 * @Author: clorf 
 * @Date: 2020-09-27 15:54:05 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-27 16:16:16
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=100010;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
inline void read(long long &x)
{
    x=0;char ch;
    while((ch=getchar())^'.');
    while(isdigit(ch=getchar())){x=(x<<1)+(x<<3)+(ch&15);}
    return;
}
int n;
long long a,b,c,d,x,y,p;
long long v;
long long gcd(long long a,long long b)
{
    if(!b)
        return a;
    return gcd(b,a%b);
}
void solve(long long a,long long b,long long c,long long d,long long &x,long long &y,bool type)//求a/b和c/d中间的分数x/y
{
    long long minn=type?a/b+1:(a-1)/b+1;
    long long maxx=type?c/d:(c-1)/d;
    if(minn<=maxx)
    {
        x=minn;
        y=1;
        return ;
    }
    if(!a)
    {
        x=1;
        y=d/c+1;
        return ;
    }
    if(a<b)
    {
        solve(d,c,b,a,x,y,type^1);
        swap(x,y);
    }
    else
    {
        solve(a%b,b,c-(a/b)*d,d,x,y,type);
        x+=y*(a/b);
    }
}
int main()
{
    while(scanf("%d",&n)!=EOF)
    {
        read(v);
        p=1;
        for(int i=1;i<=n+1;i++)   
            p*=10;
        a=max(v*10-5,0ll);
        b=p;
        long long k=gcd(a,b);
        a/=k;
        b/=k;
        c=v*10+5;
        d=p;
        k=gcd(c,d);
        c/=k;
        d/=k;
        solve(a,b,c,d,x,y,0);//type=0左边取等右边不取等,=1右边取等左边不取等
        printf("%lld\n",y);
    }
    return 0;
}

# 线性同余方程

形如 $ax\equiv c\pmod p$ 的方程式为线性同余方程。

## 乘法逆元

对于一个线性同余方程 $ax\equiv 1\pmod p$，且 $\gcd(a,p)=1$，那么 $x$ 就称为 $a\bmod p$ 的逆元，记作 $x ^{-1}$。乘法逆元的作用是在模意义下将除法变为乘法，除以 $x$ 等同于乘以 $x ^{-1}$。

求解乘法逆元的方法很多，下面一一讲解。

### 扩展欧几里得法

用扩展欧几里得来求解乘法逆元。我们可以把同余方程展开，可以得到 $ax=kp+1$，其中 $k\in \mathbb{Z}$。移项后得到 $ax+p(-k)=1$，我们令 $y=-k$，就可得到一个二元一次方程 $ax+py=1$。由于 $\gcd(a,p)=1$，所以必然有解。然后用扩展欧几里得解出 $x$ 的最小正整数解即可。

inline int exgcd(int a,int b,int &x,int &y)
{
    if(!b)
    {
        x=1;
        y=0;
        return a;
    }
    int ans=exgcd(b,a%b,x,y);
    int t=x;
    x=y;
    y=t-a/b*y;
    return ans;
}
inline int inv(int a,int p)
{
    int x,y;
    x=y=0;
    exgcd(a,p,x,y);
    x=(x%p+p)%p;
    return x;
}

### 费马小定理法

只适用于 $p$ 是质数的情况。因为 $ax\equiv 1\pmod p$，且费马小定理 $a ^{p-1}\equiv 1\pmod p$，那么我们可以做出以下化简：
$$
a ^{p-1}\equiv a\times a^{p-2}\equiv 1\pmod p
$$
所以此时 $a$ 的逆元 $x$ 即是 $a ^{p-2}\bmod p$。直接快速幂求解即可。

inline long long power(long long a,long long b,long long p)
{
    long long ans=1;
    while(b)
    {
        if(b&1)
            ans=ans*a%p;
        a=a*a%p;
        b>>=1;
    }
    return ans%p;
}
inline long long inv(long long a,long long p)
{
    return power(a,p-2,p);
}

### 线性求逆元

求 $1\sim n$ 中每个数关于 $p$ 的逆元。

我们需要在 $O(n)$ 之内求得答案。首先很显然 $1 ^{-1}\equiv 1\pmod p$，即 $1 ^{-1}=1$。然后我们把 $p$ 用带余除法打开：$p=ki+j$，即可得到 $ki+j\equiv 0\pmod p$，两边同时乘以 $i ^{-1},j ^{-1}$，即可得到：
$$
\begin{aligned}
kj ^{-1}+i ^{-1} & \equiv 0\pmod p\\
i ^{-1} & \equiv -kj ^{-1}\pmod p\\
i ^{-1} & \equiv -\left\lfloor\dfrac{p}{i}\right\rfloor(p\bmod i) ^{-1}\pmod p
\end{aligned}
$$
于是我们就得到了 $i$ 的逆元的表达式，即可线性求逆元了。需要注意的是 $1\sim n$ 中的每个数都要与 $p$ 互质，不能有 $p$ 的倍数。同时式子需要加 $p$ 来防止负数。

inv[1]=1;
for(int i=2;i<=n;i++)
    inv[i]=(p-p/i)*inv[p%i]%p;

求任意 $n$ 个数 $a _{1},a _{2}\cdots a _{n}$ 的逆元。

我们首先计算 $n$ 个数的前缀积记为 $s _{i}$，然后求出 $s _{n}$ 的逆元 $invs _{n}$。因为 $invs _{n}$ 是  $n$ 个数的积的逆元，所以乘上 $a _{n}$ 后，就会和 $a _{n}$ 的逆元抵消，得到了 $a _{1}\sim a _{n-1}$ 的积的逆元，即 $inv s_{n-1}$。

然后我们可以算出所有的 $invs _{i}$，于是 $a _{i} ^{-1}$ 可以用 $s _{i-1}\times invs _{i}$ 求到，即将其他数的逆元都抵掉即可。

s[0]=0;
for(int i=1;i<=n;i++)
    s[i]=s[i-1]*a[i]%p;
invs[n]=power(invs[n],p-2,p);
for(int i=n;i>=1;i--)
    invs[i-1]=invs[i]*a[i]%p;
for(int i=1;i<=n;i++)
    inv[i]=invs[i]*s[i-1]%p;

求 $1!,2!\cdots n!$ 的逆元。

推导方式和上面大同小异。先将 $n!$ 的逆元 $inv _{n}$ 求出，然后倒着求出 $inv _{i}$。$inv _{i}=inv _{i+1}\times (i+1)$。

fac[0]=1;
for(int i=1;i<=n;i++)
    fac[i]=fac[i-1]*i%p;
inv[n]=power(fac[n],p-2,p);
for(int i=n-1;i>=1;i--)
    inv[i]=inv[i+1]*(i+1)%p;

[Luogu3811 乘法逆元](https://www.luogu.com.cn/problem/P3811)

/*
 * @Author: clorf 
 * @Date: 2020-09-21 16:08:35 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-21 16:13:20
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=3000010;
const int mod=1e9+7;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
long long n,p,inv[maxn];
int main()
{
    scanf("%lld%lld",&n,&p);
    inv[0]=inv[1]=1;
    for(int i=2;i<=n;i++)
        inv[i]=inv[p%i]*(p-p/i)%p;
    for(int i=1;i<=n;i++)
        printf("%lld\n",inv[i]);
    return 0;
}

[Luogu5431 乘法逆元2](https://www.luogu.com.cn/problem/P5431)

这道题时间卡的很紧，但是主要有两种解法。

首先容易想到通分求解。处理出 $a _{i}$ 的前缀积和后缀积即可。

#include<cstdio>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=5e6+1;
inline void read(int &x)
{
    x=0;char ch=getchar();
    while(!isdigit(ch)) {ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    return;
}
int n,p,k,a[maxn];
int M=1,l[maxn]={1},r[maxn],b[maxn]={1},ans;
int power(int a,int b,int p)
{
    int ans=1;
    while(b)
    {
        if(b&1)
            ans=(long long)ans*a%p;
        a=(long long)a*a%p;
        b>>=1;
    }
    return ans;
}
int main()
{
    read(n),read(p),read(k);
    for(register int i=1;i<=n;++i)
    {
        read(a[i]);
        l[i]=(long long)l[i-1]*a[i]%p;
        b[i]=(long long)b[i-1]*k%p;
    }
    r[n+1]=1;
    for(register int i=n;i;--i)
    {
        r[i]=(long long)r[i+1]*a[i]%p;
        ans=(ans+(long long)b[i]*l[i-1]%p*r[i+1])%p;
    }
    printf("%d",(long long)ans*power(l[n],p-2,p)%p);
    return 0;
}

第二种方法则比较奇妙，数组甚至都不用开。我们就对每一个新的项依次通分即可，总共需要通分 $n-1$ 次。

/*
 * @Author: clorf 
 * @Date: 2020-09-15 15:10:20 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-15 16:22:49
 */
#include<cstdio>
#include<cctype>
using namespace std;
const int maxn=5000001;
inline void read(int &x)
{
    x=0;char ch=getchar();
    while(!isdigit(ch)) {ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    return;
}
int n,p,k,a[maxn],ans,num=1,M=1,x;
inline int inv(int x)
{
	return ~-x?(long long)(p-p/x)*inv(p%x)%p:1;
}
int main()
{
    read(n),read(p),read(k);
    for(register int i=1;i<=n;++i)
    {
        read(x);
        num=(long long)num*k%p;
        ans=((long long)ans*x+(long long)M*num)%p;
        M=(long long)M*x%p;
    }
    printf("%d",(long long)ans*inv(M)%p);
    return 0;
}

## 解线性同余方程

具体方法和扩展欧几里得法求逆元一样。把它用带余除法转成二元一次方程求解即可。

[Luogu1082 同余方程](https://www.luogu.com.cn/problem/P1082)

#include<cstdio>
#include<cmath>
#include<cstring>
#include<cstdlib>
#include<algorithm>
using namespace std;
long long a,b,x,y,num;
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 ans=exgcd(b,a%b,x,y);
	long long t=x;
	x=y;
	y=t-a/b*y;
	return ans;
}
int main() 
{
	scanf("%lld%lld",&a,&b);
	num=exgcd(a,b,x,y);
	x/=num;
	printf("%lld",(x%b+b)%b);
    return 0;
}

## 中国剩余定理

### 中国剩余定理(CRT)

中国剩余定理主要用于求解如下形式的一元线性同余方程组($p _{1},p _{2}\cdots p _{k}$ 两两互质)。
$$
\begin{cases}
x \equiv a _{1}\pmod {p _{1}}\\
x \equiv a _{2}\pmod {p _{2}}\\
\vdots\\
x \equiv a _{n}\pmod {p _{n}}
\end{cases}
$$
首先我们要求出 $M=\prod _{i=1} ^{n} p _{i}$，接着对于第 $i$ 个方程算出 $m _{i}=\dfrac{M}{p _{i}}$，算出 $m _{i}$ 在模 $p _{i}$ 意义下的逆元 $t _{i}$，答案即是 $ans=\sum _{i=1} ^{n}a _{i}m _{i}t _{i}\bmod M$。

[Luogu1495 中国剩余定理(CRT)](https://www.luogu.com.cn/problem/P1495)

/*
 * @Author: clorf 
 * @Date: 2020-09-23 12:22:56 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-23 15:41:22
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=100010;
const int mod=1e9+7;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
long long n,a[maxn],m[maxn],x,y,ans;
long long t[maxn],M[maxn],all=1;
long long exgcd(long long a,long long b,long long &x,long long &y)
{
    if(!b)
    {
        x=1;
        y=0;
        return a;
    }
    long long ans=exgcd(b,a%b,x,y);
    long long t=x;
    x=y;
    y=t-a/b*y;
    return ans;
}
long long inv(long long a,long long p)
{
    x=y=0;
    long long num=exgcd(a,p,x,y);
    return (x%p+p)%p;
}
int main()
{
    scanf("%lld",&n);
    for(int i=1;i<=n;i++)
    {
        scanf("%lld%lld",&m[i],&a[i]);
        all*=m[i];
    }
    for(int i=1;i<=n;i++)
        M[i]=all/m[i];
    for(int i=1;i<=n;i++)
        t[i]=inv(M[i],m[i]);
    for(int i=1;i<=n;i++)
        ans=(ans+a[i]*M[i]%all*t[i]%all)%all;
    printf("%lld",(ans%all+all)%all);
    return 0;
}

### 扩展中国剩余定理(EXCRT)

扩展中国剩余定理适用于解如上同余方程时，模数 $p _{i}$ 不两两互质的情况。针对于这种情况有两种方法求解。

第一种方法是合并两个同余方程法。我们每次把两个同余方程合并成一个同余方程，这样合并 $n-1$ 后，就只用对剩下的那个同余方程求解即可。我们设两个方程分别如下所示：
$$
\begin{aligned}
x & \equiv a _{1}\pmod {p _{1}}\\
x & \equiv a _{2}\pmod {p _{2}}
\end{aligned}
$$
于是我们得到 $x=p _{1}i+a _{1}=p _{2}j+a _{2}$，其中 $i,j\in \mathbb{Z}$。继续可以得到 $p _{1}i-p _{2}j=a _{2}-a _{1}$。由裴蜀定理可知，当 $a _{2}-a _{1}$ 不能被 $\gcd(p _{1},p _{2})$ 整除时无解。否则可以用扩展欧几里得解出一组解 $(i _{0},j _{0})$，就得到了这两个方程 $x$ 的特解 $x _{0}=p _{1}i _{0}+a _{1}$。同时可以得出 $i$ 的通解 $i=i _{0}+k\dfrac{p _{2}}{\gcd(p _{1},p _{2})},k\in \mathbb{Z}$，带回去得到 $x$ 的通解 $x=k\dfrac{p _{1}p _{2}}{\gcd(p _{1},p _{2})}+p _{1}i _{0}+a _{1}=k\operatorname{lcm}(p _{1},p _{2})+x _{0}$。然后把它转成同余方程的形式，即：
$$
x\equiv x _{0}\pmod{\operatorname{lcm}(p _{1},p _{2})}
$$
这样便将两个方程合并成了一个。多个方程两两合并后会变成一个方程：
$$
x\equiv x _{n}\pmod{\operatorname{lcm}(p _{1},p _{2},\cdots,p _{n})}
$$
令 $m=\operatorname{lcm}(p _{1},p _{2},\cdots,p _{n})$，那么 $(x _{n}\bmod m+m)\bmod m$ 即是答案。

第二种方法更容易理解。我们假设我们已经求出了前 $i-1$ 个方程的特解 $x$，正在求解第 $i$ 个方程。设 $m=\operatorname{lcm}(p _{1},p _{2},\cdots,p _{n})$，那么 $x+km(k\in \mathbb{Z})$ 是前 $i$ 个方程的通解。代入第 $i$ 个方程，即可得到：
$$
x+km\equiv a _{i}\pmod {p _{i}}
$$
其中 $k$ 是未知量，我们可以用扩展欧几里得求出 $k$，于是我们就得到了前 $i$ 个方程的特解 $x'=x+km$。综上，使用 $n$ 次扩展欧几里得算法即可。

[Luogu4777 扩展中国剩余定理(EXCRT)](https://www.luogu.com.cn/problem/P4777)

/*
 * @Author: clorf 
 * @Date: 2020-09-23 15:45:37 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-23 17:19:44
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=100010;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
long long n,a[maxn],p[maxn];
long long m,ans,x,y;
long long exgcd(long long a,long long b,long long &x,long long &y)
{
    if(!b)
    {
        x=1;
        y=0;
        return a;
    }
    long long ans=exgcd(b,a%b,x,y);
    long long t=x;
    x=y;
    y=t-a/b*y;
    return ans;
}
inline long long multi(long long a,long long b,long long p)
{
    long long ans=0;
    for(;b;b>>=1)
    {
        if(b&1)
            ans=(ans+a)%p;
        a=(a+a)%p;
    }
    return ans%p;
}
int main()
{
    scanf("%lld",&n);
    for(int i=1;i<=n;i++)
        scanf("%lld%lld",&p[i],&a[i]);
    m=p[1];
    ans=a[1];
    for(int i=2;i<=n;i++)
    {
        long long num=exgcd(m,p[i],x,y);
        long long c=(a[i]-ans%p[i]+p[i])%p[i];
        long long mod=p[i]/num;
        if(c%num)
        {
            printf("-1");
            return 0;
        }
        x=multi(x,c/num,mod);
        x=(x+mod)%mod;
        ans+=x*m;
        m*=mod;
        ans=(ans%m+m)%m;
    }
    printf("%lld",ans);
    return 0;
}

[Luogu4774 屠龙勇士](https://www.luogu.com.cn/problem/P4774)

一道扩展中国剩余定理的套路题。先用 `set` 维护出每次的攻击力，然后就化成了求解一个同余方程组的问题了。如下所示：
$$
\begin{cases}
b _{1}x \equiv a _{1}\pmod{p _{1}}\\
b _{2}x \equiv a _{2}\pmod{p _{2}}\\
\vdots\\
b _{n}x\equiv a _{n}\pmod{p _{n}}
\end{cases}
$$

其实和扩展中国剩余定理的式子差不多，多了一个系数而已。我们依然按照推扩展中国剩余定理的式子那样去推答案即可。注意每次砍多少刀有个下限，至少一定要砍那么多刀(针对于 $\forall p _{i}=1$ 的情况)。

/*
 * @Author: clorf 
 * @Date: 2020-09-29 10:24:46 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-29 16:55:19
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#include<set>
#define INF 1e9
using namespace std;
const int maxn=200010;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
int t,n,m;
long long a[maxn],p[maxn];
long long r[maxn],mich[11];
multiset<long long> e;
long long exgcd(long long a,long long b,long long &x,long long &y)
{
    if(!b)
    {
        x=1;
        y=0;
        return a;
    }
    long long ans=exgcd(b,a%b,x,y);
    long long t=x;
    x=y;
    y=t-a/b*y;
    return ans;
}
int main()
{
    scanf("%d",&t);
    for(int _=1;_<=t;_++)
    {
        e.clear();
        bool flag=0;
        long long x,y;
        x=y=0;
        scanf("%d%d",&n,&m);
        for(int i=1;i<=n;i++)
            scanf("%lld",&a[i]);
        for(int i=1;i<=n;i++)
            scanf("%lld",&p[i]);
        for(int i=1;i<=n;i++)
            scanf("%lld",&r[i]);
        for(int i=1;i<=m;i++)
        {
            scanf("%lld",&x);
            e.insert(x);
        }
        long long atk,all=1,ans=0,maxx=0;
        for(int i=1;i<=n;i++)
        {
            if(e.upper_bound(a[i])==e.begin())
                atk=*e.begin();
            else
                atk=*(--e.upper_bound(a[i]));
            long long A=(__int128)atk*all%p[i];
            long long B=p[i];
            long long C=(a[i]-atk*ans%p[i]+p[i])%p[i];
            long long num=exgcd(A,B,x,y);
            x=(x%B+B)%B;
            if(C%num)
            {
                flag=1;
                break;
            }
            ans+=(__int128)(C/num)*x%(B/num)*all;
            all*=B/num;
            ans%=all;
            maxx=max(maxx,(a[i]-1)/atk+1);//每次至少要砍这么多刀
            e.erase(e.find(atk));
            e.insert(r[i]);
        }
        if(ans<maxx)    
            ans=maxx;
        if(flag)
        {
            mich[_]=-1;
            continue;
        }
        mich[_]=ans;
    }
    for(int i=1;i<=t;i++)
        printf("%lld\n",mich[i]);
    return 0;
}

# BSGS

## BSGS

BSGS(baby step giant step)，中文译作大步小步算法。它能在 $O(\sqrt{p})$ 的时间内求解以下高次同余方程：
$$
a ^{x}\equiv b\pmod{p}
$$
其中 $\gcd(a,p)=1$，解满足 $0\le x<p$。

我们令 $x=k\left\lceil\sqrt{p}\right\rceil-t$，其中 $0\le k,t\le \left\lceil\sqrt{p}\right\rceil$。则有 $a ^{k\sqrt{p}-t}\equiv b\pmod{p}$。变换以下即 $a ^{k\left\lceil\sqrt{p}\right\rceil}\equiv ba ^{t}\pmod{p}$。

我们枚举 $t$ 从 $0\sim \left\lceil\sqrt{p}\right\rceil$，把 $ba ^{t}$ 插入 `map` 或哈希表中，然后枚举 $k$ 从 $0\sim \left\lceil\sqrt{p}\right\rceil$，计算 $a ^{k\left\lceil\sqrt{p}\right\rceil}$，看看 `map` 或哈希表中是否有相等的 $ba ^{t}$，如果有答案就是 $k\left\lceil\sqrt{p}\right\rceil-t$。注意要特判底数 $a$ 为 $0$ 的情况。

[Luogu3846 可爱的质数](https://www.luogu.com.cn/problem/P3846)

/*
 * @Author: clorf 
 * @Date: 2020-09-24 16:33:30 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-24 19:02:09
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#include<unordered_map>
#define INF 1e9
using namespace std;
const int maxn=100010;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
namespace BSGS
{
    long long a,b,p;
    inline long long power(long long a,long long b,long long p)
    {
        long long ans=1;
        for(;b;b>>=1)
        {
            if(b&1)
                ans=ans*a%p;
            a=a*a%p;
        }
        return ans%p;
    }
    inline long long BSGS(long long a,long long b,long long p)
    {
        unordered_map<long long,int> hash;
        b%=p;
        int t=sqrt(p)+1;
        long long val=b;
        hash[val]=0;
        for(int j=1;j<t;j++)
        {
            val=val*a%p;
            hash[val]=j;
        }
        if(a==0)
        {
            if(b==0)
                return 1;
            else
                return -1;
        }
        long long num=power(a,t,p)%p;
        val=1;
        for(int i=1;i<=t;i++)
        {
            val=val*num%p;
            int j=hash.find(val)==hash.end()?-1:hash[val];
            if(j>=0&&i*t-j>=0)
                return i*t-j;
        }
        return -1;
    }
    inline void solve()
    {
        scanf("%lld%lld%lld",&p,&a,&b);
        long long ans=BSGS(a,b,p);
        if(ans>=0)
            printf("%lld",ans);
        else
            printf("no solution");
    }
}
int main()
{
    BSGS::solve();
    return 0;
}

## 扩展BSGS

即求解上题 $a,p$ 不一定互质的情况。

对于 $a,p$ 不互质的情况，我们设 $d _{1}=\gcd(a,p)$，如果 $d _{1}\nmid b$，那么方程无解。否则我们同时除以 $d _{1}$：
$$
\dfrac{a}{d _{1}}\cdot a ^{x-1}\equiv \dfrac{b}{d _{2}}\pmod{\dfrac{p}{d _{1}}}
$$
若 $a$ 和$\dfrac{p}{d _{1}}$ 还不互质就继续除，直到 $a$ 与 $\dfrac{p}{d _{1}d _{2}\cdots d _{k}}$ 互质，设 $m=\prod _{i=1} ^{k}d _{i}$，那么方程变成了这样：
$$
\dfrac{a ^{k}}{m}\cdot a ^{x-k}\equiv \dfrac{b}{m}\pmod{\dfrac{p}{m}}
$$
注意在这一步中如果某一步 $\dfrac{a ^{k'}}{m'}=\dfrac{b}{m'}$，那么此时 $k$ 就是解，直接返回即可。

此时由于 $\dfrac{a ^{k}}{m}$ 与 $\dfrac{p}{m}$ 互质，我们可以求出 $\dfrac{a ^{k}}{m}$ 关于 $\dfrac{p}{m}$ 的逆元移到右边，于是就变成了普通的 BSGS 问题，我们求出 $x-k$ 再加上 $k$ 即是答案。

[Luogu4195 扩展BSGS](https://www.luogu.com.cn/problem/P4195)

/*
 * @Author: clorf 
 * @Date: 2020-09-24 17:50:47 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-24 18:55:07
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#include<unordered_map>
#define INF 1e9
using namespace std;
const int maxn=100010;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
long long a,b,p;
inline long long power(long long a,long long b,long long p)
{
    long long ans=1;
    for(;b;b>>=1)
    {
        if(b&1)
            ans=ans*a%p;
        a=a*a%p;
    }
    return ans;
}
inline long long gcd(long long a,long long b)
{
    if(!b)
        return a;
    return gcd(b,a%b);
}
long long exgcd(long long a,long long b,long long &x,long long &y)
{
    if(!b)
    {
        x=1;
        y=0;
        return a;
    }
    long long ans=exgcd(b,a%b,x,y);
    long long t=x;
    x=y;
    y=t-a/b*y;
    return ans;
}
inline long long inv(long long a,long long p)
{
    long long x,y;
    x=y=0;
    exgcd(a,p,x,y);
    return (x%p+p)%p;
}
inline long long BSGS(long long a,long long b,long long p,long long k)
{
    unordered_map<long long,int> hash;
    hash.clear();
    b%=p;
    int t=sqrt(p)+1;
    long long val=b;
    hash[val]=0;
    for(int j=1;j<t;j++)
    {
        val=val*a%p;
        hash[val]=j;
    }
    if(a==0)
    {
        if(b==0)
            return 1;
        else
            return -1;
    }
    long long num=power(a,t,p)%p;
    val=1;
    for(int i=1;i<=t;i++)
    {
        val=val*num%p;
        int j=hash.find(val)==hash.end()?-1:hash[val];
        if(j>=0&&i*t-j+k>=0)
            return i*t-j+k;
    }
    return -1;
}
inline long long exBSGS(long long a,long long b,long long p)
{
    if(b==1)
        return 0;
    long long k=0,x=1,num=1;
    while(1)
    {
        long long d=gcd(a,p);
        if(d==1)
            break;
        if(b%d)
            return -1;
        p/=d;
        b/=d;
        k++;
        x=x*d%p;
        num=num*a/d%p;
        if(num==b)
            return k;
    }
    long long sum=power(a,k,p)*inv(x,p)%p;
    b=b*inv(sum,p)%p;
    return BSGS(a,b,p,k);
}
inline void solve()
{
    while(scanf("%lld%lld%lld",&a,&p,&b)!=EOF&&a+p+b)
    {
        long long ans=exBSGS(a,b,p);
        if(ans>=0)
            printf("%lld\n",ans);
        else
            printf("No Solution\n");
    }
}
int main()
{
    solve();
    return 0;
}

# 卢卡斯定理

卢卡斯($\text{Lucas}$)定理，用于求解大组合数取模的问题。

## 卢卡斯定理

对于质数 $p$，有：
$$
\dbinom{n}{m}\bmod p=\dbinom{\left\lfloor\dfrac{n}{p}\right\rfloor}{\left\lfloor\dfrac{m}{p}\right\rfloor}\cdot\binom{n\bmod p}{m\bmod p}\bmod p
$$
对于 $\dbinom{n\bmod p}{m\bmod p}$ 我们用组合数直接求解，对于 $\dbinom{\left\lfloor\dfrac{n}{p}\right\rfloor}{\left\lfloor\dfrac{m}{p}\right\rfloor}$ 可以继续用 Lucas 定理求解。

[Luogu3807 卢卡斯定理](https://www.luogu.com.cn/problem/P3807)

/*
 * @Author: clorf 
 * @Date: 2020-09-23 19:14:09 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-23 19:41:41
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=100010;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
int t;
long long fac[maxn],inv[maxn],n,m,p;
long long power(long long a,long long b,long long p)
{
    long long ans=1%p;
    for(;b;b>>=1)
    {
        if(b&1)
            ans=ans*a%p;
        a=a*a%p;
    }
    return ans%p;
}
long long C(long long n,long long m,long long p)
{
    return fac[n]*inv[m]%p*inv[n-m]%p;
}
long long lucas(long long n,long long m,long long p)
{
    if(!m)
        return 1;
    return lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}
int main()
{
    scanf("%d",&t);
    while(t--)
    {
        scanf("%lld%lld%lld",&n,&m,&p);
        fac[0]=inv[0]=1;
        for(int i=1;i<=n+m;i++)
        {
            fac[i]=fac[i-1]*i%p;
            inv[i]=power(fac[i],p-2,p);
        }
        printf("%lld\n",lucas(n+m,n,p));
    }
    return 0;
}
//有些阶乘是p的倍数，就不存在逆元，因此不能用阶乘逆元倒推。

[Luogu2480 古代猪文](https://www.luogu.com.cn/problem/P2480)

给定整数 $q,n(1\le q,n\le 10 ^{9})$，计算 $q ^{\sum _{d|n}\binom{n}{d}}\bmod 999911659$。

首先 $\sum _{d|n}\binom{n}{d}$ 太大了，由于 $999911659$ 是质数，我们用欧拉定理降幂，指数变为 $\sum _{d|n}\binom{n}{d}\bmod 999911658$。现在重点就在于如何求解这个。

首先已知 $999911658=2\times 3\times 4679\times 35617$。我们称这样所有质因数的指数都是 $1$ 的数为 square free number。

枚举 $n$ 的约数 $d$，运用卢卡斯定理求 $\dbinom{n}{d}$，分别计算出 $\sum _{d|n}\binom{n}{d}$ 对这四个数取模的结果，然后用中国剩余定理即可。最后用快速幂求出答案。

/*
 * @Author: clorf 
 * @Date: 2020-09-25 15:24:22 
 * @Last Modified by:   clorf 
 * @Last Modified time: 2020-09-25 15:24:22 
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=40010;
const int mod=999911658;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
long long n,g;
long long m[5]={0,2,3,4679,35617};
long long a[5],fac[maxn];
inline long long power(long long a,long long b,long long p)
{
    long long ans=1;
    for(;b;b>>=1)
    {
        if(b&1)
            ans=ans*a%p;
        a=a*a%p;
    }
    return ans%p;
}
inline long long inv(long long x,long long p)
{
    return power(x,p-2,p);
}
inline long long C(long long n,long long m,long long p)
{
    if(n<m)
        return 0;
    return fac[n]*inv(fac[m],p)%p*inv(fac[n-m],p)%p;
}
long long lucas(long long n,long long m,long long p)
{
    if(n<m)
        return 0;
    if(!n)
        return 1;
    return lucas(n/p,m/p,p)%p*C(n%p,m%p,p)%p;
}
inline void ready(int n)
{
    fac[0]=1;
    for(int i=1;i<=n;i++)   
        fac[i]=fac[i-1]*i%n;
}
inline long long CRT()
{
    long long ans=0;
    for(int i=1;i<=4;i++)
        ans=(ans+a[i]*(mod/m[i])%mod*power(mod/m[i],m[i]-2,m[i]))%mod;
    return ans;
}
int main()
{
    scanf("%lld%lld",&n,&g);
    if(g==mod+1)
    {
        printf("0");
        return 0;
    }
    for(int i=1;i<=4;i++)
    {
        ready(m[i]);
        for(int j=1;j*j<=n;j++)
            if(n%j==0)
            {
                a[i]=(a[i]+lucas(n,j,m[i])%m[i])%m[i];
                int k=n/j;
                if(j!=k)
                    a[i]=(a[i]+lucas(n,k,m[i])%m[i])%m[i];
            }
    }
    long long ans=CRT();
    printf("%lld",power(g,ans,mod+1));
    return 0;
}

## 扩展卢卡斯定理

扩展卢卡斯定理用于解决上题 $p$ 不是素数的情况。

首先对 $p$ 进行质因数分解，设 $p=p _{1} ^{c _{1}}p _{2} ^{c _{2}}\cdots p _{k} ^{c _{k}}$，如果可以求出每个 $\dbinom{n}{m}\equiv a _{i}\pmod{p _{i} ^{c _{i}}}$，于是就转移成了以下方程组：
$$
\begin{cases}
\dbinom{n}{m} \equiv a _{1}\pmod{p _{1} ^{c _{1}}} \\
\dbinom{n}{m} \equiv a _{2}\pmod{p _{2} ^{c _{2}}} \\
\vdots\\
\dbinom{n}{m} \equiv a _{k}\pmod{p _{k} ^{c _{k}}} \\
\end{cases}
$$
我们可以用中国剩余定理求出 $\dbinom{n}{m}$ 的值。

 接下来就是求 $a _{i}=\dbinom{n}{m}\bmod p _{i} ^{c _{i}}$ 了。我们已知组合数公式为 $\dbinom{n}{m}=\dfrac{n!}{m!(n-m)!}$，发现 $m!$ 和 $(n-m)!$ 可能包含质因子 $p _{i}$，所以不能直接求它们对于 $p _{i} ^{c _{i}}$ 的逆元。我们选择将 $n!,m!,(n-m)!$ 的质因子 $p$ 全部提出来，最后再乘回去，即转化为如下形式：
$$
\dfrac{\dfrac{n!}{p _{i} ^{a}}}{\dfrac{m!}{p _{i} ^{b}}\times \dfrac{(n-m)!}{p _{i} ^{c}}}\times p _{i} ^{a-b-c}
$$
其中 $\dfrac{m!}{p _{i} ^{b}},\dfrac{(n-m)!}{p _{i} ^{c}}$ 与 $p _{i} ^{c _{i}}$ 互质，可以直接求逆元。

现在问题就转化到了这里，如何求以下式子：
$$
\dfrac{n!}{p ^{a}}\bmod p ^{k}
$$
 我们先考虑如何算 $n!\bmod p ^{k}$。我们举个例子来理解：$n=22,p=3,k=2$。
$$
\begin{aligned}
22! & =1\times 2\times 3\times 4\times 5\times 6\times 7\times 8\times 9\times 10\times 11\times 12\\
& \times 13\times 14\times 15\times 16\times 17\times 18\times 19\times 20\times 21\times 22\\
& =3 ^{7}\times(1\times 2\times 3\times 4\times 5\times 6\times 7)\times(1\times 2\times 4\times 5\\
& \times 7\times 8\times 10\times 11\times 13\times 14\times 16\times 17\times 19\times 20\times 22)\\
& =3 ^{7}\times 7!\times (1\times 2\times 4\times 5 \times 7\times 8\times 10\times 11\times 13\times 14\\
& \times 16\times 17\times 19\times 20\times 22)
\end{aligned}
$$
通过把 $p$ 的倍数提出来之后，可以发现上面的式子变成了 $3$ 个部分，第一个部分是 $3$ 的幂，次数是小于等于 $22$ 的 $3$ 的倍数的个数，即 $\left\lfloor\dfrac{n}{p}\right\rfloor$。故第一部分是 $p ^{\lfloor\frac{n}{p}\rfloor}$；第二部分是 $7!$，即 $\left\lfloor\dfrac{n}{p}\right\rfloor !$；第 $3$ 部分是 $n!$ 中与 $p$ 互质的数的乘积，这一部分在模 $p ^{k}$ 意义下具有循环的性质，即：
$$
1\times 2\times 4\times 5\times 7\times 8\equiv 10\times 11\times 13\times 14\times 16\times 17\pmod {3 ^{2}}
$$
循环节共循环了 $\left\lfloor\dfrac{n}{p ^{k}}\right\rfloor$ 次，我们暴力求出 $p ^{k}$ 之内与其互质的数的乘积，然后用快速幂求出它的 $\left\lfloor\dfrac{n}{p ^{k}}\right\rfloor$ 次幂，最后乘上多出来的那一截 $19\times 20\times 22$，即 $\prod _{\gcd(i,p ^{k})=1} ^{n~\bmod~ p ^{k}}i$，暴力乘上去即可。

最后由于第一部分要除掉，我们将第二和第三部分乘起来即可，注意第二部分要递归求解，因为在这个阶乘里可能还有 $p$ 的倍数没有除完。最后还有一个优化，我们可以每次先把 $p ^{k}$ 内与 $p$ 互质的数先预处理出来，即用一个数组 $fac _{i}$ 代表 $1\sim i$ 中与 $p$ 互质的数的乘积模 $p ^{k}$，到时候直接调用即可，这样效率会提升很多(亲测)。

最后回到这个式子：
$$
\dfrac{\dfrac{n!}{p _{i} ^{a}}}{\dfrac{m!}{p _{i} ^{b}}\times \dfrac{(n-m)!}{p _{i} ^{c}}}\times p _{i} ^{a-b-c}
$$
借助上面的方法求出分子分母三个式子，求出逆元与 $p _{i} ^{a-b-c}$ 即可。

然后质因数分解对每个质数这样来一遍，中国剩余定理求解。

(https://www.luogu.com.cn/problem/P4720)[Luogu4720 扩展卢卡斯]

/*
 * @Author: clorf 
 * @Date: 2020-09-23 19:42:16 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-23 22:40:40
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=1000100;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
namespace EXlucas
{
    long long n,m,p,a[maxn],mod[maxn],cnt,fac[maxn];
    inline long long power(long long a,long long b,long long p)
    {
        long long ans=1;
        for(;b;b>>=1)
        {
            if(b&1)
                ans=ans*a%p;
            a=a*a%p;
        }
        return ans%p;
    }
    long long exgcd(long long a,long long b,long long &x,long long &y)
    {
        if(!b)
        {
            x=1;
            y=0;
            return a;
        }
        long long ans=exgcd(b,a%b,x,y);
        long long t=x;
        x=y;
        y=t-a/b*y;
        return ans;
    }
    inline long long inv(long long a,long long p)
    {
        long long x,y;
        x=y=0;
        long long num=exgcd(a,p,x,y);
        return (x%p+p)%p;
    }
    inline long long CRT()
    {
        long long all=1,ans=0;
        for(int i=1;i<=cnt;i++)
            all*=mod[i];
        for(int i=1;i<=cnt;i++)
            ans=(ans+a[i]*(all/mod[i])%all*inv(all/mod[i],mod[i])%all)%all;
        return ans;
    }
    long long calc(long long n,long long p,long long pk)
    {
        if(!n)
            return 1;
        if(n<p)
            return fac[n];
        long long ans;
        ans=power(fac[pk-1],n/pk,pk);
        ans=ans*fac[n%pk]%pk;
        return ans*calc(n/p,p,pk)%pk;
    }
    inline long long C(long long n,long long m,long long p,long long pk)
    {
        if(n<m)
            return 0;
        fac[0]=1;
        for(int i=1;i<=pk-1;i++)
        {
            if(i%p)
                fac[i]=fac[i-1]*i%pk;
            else
                fac[i]=fac[i-1];
        }//预处理！！！
        long long num1=calc(n,p,pk);
        long long num2=calc(m,p,pk);
        long long num3=calc(n-m,p,pk);
        long long num=0;
        for(long long i=n;i;i/=p)
            num+=i/p;
        for(long long i=m;i;i/=p)
            num-=i/p;
        for(long long i=n-m;i;i/=p)
            num-=i/p;
        long long ans=num1*inv(num2,pk)%pk*inv(num3,pk)%pk*power(p,num,pk)%pk;
        return ans;
    }
    long long exlucas(long long n,long long m,long long p)
    {
        cnt=0;
        long long q=sqrt(p);
        for(int i=2;p>1&&i<=q;i++)
        {
            long long num=1;
            while(p%i==0)
            {
                p/=i;
                num*=i;
            }
            if(num>1)
            {
                a[++cnt]=C(n,m,i,num);
                mod[cnt]=num;
            }
        }
        if(p>1)
        {
            a[++cnt]=C(n,m,p,p);
            mod[cnt]=p;
        }
        return CRT();
    }
    inline void solve()
    {
        scanf("%lld%lld%lld",&n,&m,&p);
        printf("%lld\n",exlucas(n,m,p));
    }
};
int main()
{
    EXlucas::solve();
    return 0;
}

[Luogu2183 礼物](https://www.luogu.com.cn/problem/P2183)

裸的扩展卢卡斯，推出答案的式子为如下所示：
$$
\prod _{i=1} ^{m}\dbinom{n-\sum _{j=1} ^{i-1}w _{j}}{w _{i}}\bmod p
$$
直接扩展卢卡斯。

/*
 * @Author: clorf 
 * @Date: 2020-09-29 19:44:20 
 * @Last Modified by:   clorf 
 * @Last Modified time: 2020-09-29 19:44:20 
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=7;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
namespace EXlucas
{
    long long a[maxn],mod[maxn],cnt;
    inline long long power(long long a,long long b,long long p)
    {
        long long ans=1;
        for(;b;b>>=1)
        {
            if(b&1)
                ans=ans*a%p;
            a=a*a%p;
        }
        return ans%p;
    }
    long long exgcd(long long a,long long b,long long &x,long long &y)
    {
        if(!b)
        {
            x=1;
            y=0;
            return a;
        }
        long long ans=exgcd(b,a%b,x,y);
        long long t=x;
        x=y;
        y=t-a/b*y;
        return ans;
    }
    inline long long inv(long long a,long long p)
    {
        long long x,y;
        x=y=0;
        long long num=exgcd(a,p,x,y);
        return (x%p+p)%p;
    }
    inline long long CRT()
    {
        long long all=1,ans=0;
        for(int i=1;i<=cnt;i++)
            all*=mod[i];
        for(int i=1;i<=cnt;i++)
            ans=(ans+a[i]*(all/mod[i])%all*inv(all/mod[i],mod[i])%all)%all;
        return ans;
    }
    long long fac(long long n,long long p,long long pk)
    {
        if(!n)
            return 1;
        long long ans=1;
        for(int i=1;i<pk;i++)
            if(i%p)
                ans=ans*i%pk;
        ans=power(ans,n/pk,pk);
        for(int i=1;i<=n%pk;i++)
            if(i%p)
                ans=ans*i%pk;
        return ans*fac(n/p,p,pk)%pk;
    }
    inline long long C(long long n,long long m,long long p,long long pk)
    {
        if(n<m)
            return 0;
        long long num1=fac(n,p,pk);
        long long num2=fac(m,p,pk);
        long long num3=fac(n-m,p,pk);
        long long num=0;
        for(long long i=n;i;i/=p)
            num+=i/p;
        for(long long i=m;i;i/=p)
            num-=i/p;
        for(long long i=n-m;i;i/=p)
            num-=i/p;
        long long ans=num1*inv(num2,pk)%pk*inv(num3,pk)%pk*power(p,num,pk)%pk;
        return ans;
    }
    long long exlucas(long long n,long long m,long long p)
    {
        cnt=0;
        long long q=sqrt(p);
        for(int i=2;p>1&&i<=q;i++)
        {
            long long num=1;
            while(p%i==0)
            {
                p/=i;
                num*=i;
            }
            if(num>1)
            {
                a[++cnt]=C(n,m,i,num);
                mod[cnt]=num;
            }
        }
        if(p>1)
        {
            a[++cnt]=C(n,m,p,p);
            mod[cnt]=p;
        }
        return CRT();
    }
}
int mod,n,m;
long long w[maxn],ans=1,sum;
int main()
{
    scanf("%d%d%d",&mod,&n,&m);
    for(int i=1;i<=m;i++)
    {
        scanf("%lld",&w[i]);
        sum+=w[i];
    }
    if(sum>n)
    {
        printf("Impossible\n");
        return 0;
    }
    for(int i=1;i<=m;i++)
    {
        ans=(__int128)(ans*EXlucas::exlucas(n,w[i],mod))%mod;
        n-=w[i];
    }
    printf("%lld\n",ans);
    return 0;
}

# 拉格朗日插值

~~其实严格来说应该不算数论知识，应该算多项式。~~

拉格朗日插值能将多项式的点值表达式转为系数表达式。

给出 $n$ 个点 $p _{i}(x _{i},y _{i})$ 和 $k$，将这 $n$ 个点的最多 $n-1$ 的多项式记为 $f(x)$，求 $f(k)$。

我们用待定系数法可以变成一个 $n$ 个 $n$ 元一次方程组成的方程组，然后用高斯消元求解，时间复杂度为 $O(n ^{3})$。

但是这种方法的复杂度太高了，我们需要更简单的方法。

拉格朗日插值能在 $O(n ^{2})$ 之内求解。它给出了 $f(x)$ 如何用点值表达式求，如下所示：
$$
f(x)=\sum _{i=1} ^{n}y _{i}\prod _{j\ne i}\dfrac{x-x _{j}}{x _{i}-x _{j}}
$$
我们直接将 $k$ 代入求解即可。

[Luogu4781 拉格朗日插值](https://www.luogu.com.cn/problem/P4781)

/*
 * @Author: clorf 
 * @Date: 2020-09-23 17:53:58 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-23 18:05:23
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=2010;
const double Pi=acos(-1.0);
const int mod=998244353;
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
long long n,k,ans;
long long x[maxn],y[maxn];
long long power(long long a,long long b,long long p)
{
    long long ans=1;
    for(;b;b>>=1)
    {
        if(b&1)
            ans=ans*a%p;
        a=a*a%p;
    }
    return ans%p;
}
long long inv(long long x)
{
    return power(x,mod-2,mod);
}
int main()
{
    scanf("%lld%lld",&n,&k);
    for(int i=1;i<=n;i++)
        scanf("%lld%lld",&x[i],&y[i]);
    for(int i=1;i<=n;i++)
    {
        long long sum1=y[i]%mod;
        long long sum2=1;
        for(int j=1;j<=n;j++)
            if(i!=j)
            {
                sum1=sum1*(k-x[j])%mod;
                sum2=sum2*((x[i]%mod-x[j]%mod)%mod)%mod;
            }
        long long sum=sum1*inv(sum2)%mod;
        ans=(ans+sum)%mod;
        ans=(ans+mod)%mod;
    }
    printf("%lld",ans);
    return 0;
}

[Luogu4593 教科书般的亵渎](https://www.luogu.com.cn/problem/P4593)

这道题值得一做。首先易知 $k=m+1$，因为每一张亵渎必会杀死一个区间的怪物。我们算出在每一个空位上使用亵渎的贡献，设空位为 $p$，那么有贡献的区间为 $p+1\sim n$，贡献为 $\sum _{i=1} ^{n-p}i ^{k}$，然后要减去后面空位多算的贡献，那么答案即为：
$$
\sum _{j=0} ^{m}\left(\sum _{i=1} ^{n-a _{j}}i ^{k}-\sum _{i=j+1} ^{m}(a _{i}-a _{j}) ^{m+1}\right)
$$
其中出现了 $f(x)=\sum _{i=1} ^{x} i ^{k}$ 的式子，这样的式子叫做自然数幂和。它是一个 $k+1$ 次多项式，证明自行查找。我们用拉格朗日插值来求即可。

同时对于这种 $x$ 取值连续的插值，我们能优化到 $O(n)$ 复杂度。我们首先把 $x _{i}$ 变成 $i$，新的式子为 :
$$
f(k)=\sum _{i=0} ^{n}y _{i}\prod _{i\ne j}\dfrac{k-j}{i-j}
$$
如何快速计算 $\prod _{i\ne j}\dfrac{k-j}{i-j}$？对于分子我们可以维护出关于 $k$ 的前缀积和后缀积，即：
$$
\begin{aligned}
pre _{i} & =\prod _{j=0} ^{i}k-j\\
suf _{i} & =\prod _{j=i} ^{n}k-j
\end{aligned}
$$
对于分母来说，可以发现这就是阶乘的形式，那么式子即是：
$$
f(k)=\sum _{i=0} ^{n}y _{i}\dfrac{pre _{i-1} \times suf _{i+1}}{fac _{i}\times fac _{n-i}}
$$
需要注意的是分母可能会出现符号问题，当 $n-i$ 为奇数时应该变号。

最后这道题有一个小细节，如果末尾有一段连续空位，应该去掉它，因为它不需要多余的亵渎。

/*
 * @Author: clorf 
 * @Date: 2020-09-24 19:27:06 
 * @Last Modified by: clorf
 * @Last Modified time: 2020-09-24 20:34:05
 */
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<ctime>
#include<cctype>
#define INF 1e9
using namespace std;
const int maxn=55;
const int mod=1e9+7;
const double Pi=acos(-1.0);
const double eps=1e-7;
//1.数组开小
//2.没用long long/爆了类型存储范围
//3.写之前式子没推清楚
//4.变量写错(想当然typo/没想清楚含义)
//5.写之前没想清楚复杂度
template<class T>void read(T &x)
{
    x=0;int f=0;char ch=getchar();
    while(!isdigit(ch)) {f|=(ch=='-');ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}
    x=f?-x:x;
    return;
}
int t;
long long n,m,a[maxn];
long long k,d,fac[maxn];
long long y[maxn],ans;
long long pre[maxn],suf[maxn];
inline long long power(long long a,long long b,long long p)
{
    long long ans=1;
    for(;b;b>>=1)
    {
        if(b&1)
            ans=ans*a%p;
        a=a*a%p;
    }
    return ans%p;
}
inline long long Lagrange(long long x)//处理连续区间的插值
{
    pre[0]=1;
    suf[d+1]=1;
    for(int i=1;i<=d;i++)
        pre[i]=(pre[i-1]*(x-i)%mod+mod)%mod;
    for(int i=d;i>=1;i--)
        suf[i]=(suf[i+1]*(x-i)%mod+mod)%mod;
    long long sum1,sum2,t1,t2;
    sum1=sum2=t2=0;//分子分母
    t1=1;
    for(int i=1;i<=d;i++)
    {
        sum2=pre[i-1]*suf[i+1]%mod*y[i]%mod;
        t2=fac[i-1]*fac[d-i]%mod;
        if((d-i)&1)//d-i奇数，变号
            t2=(mod-t2)%mod;
        sum1=(sum1*t2%mod+sum2*t1%mod)%mod;
        t1=t1*t2%mod;//通分
    }
    return sum1*power(t1,mod-2,mod);
}
int main()
{
    scanf("%d",&t);
    while(t--)
    {
        scanf("%lld%lld",&n,&m);
        for(int i=1;i<=m;i++)
            scanf("%lld",&a[i]);
        sort(a+1,a+m+1);
        if(a[m]==n)
        {
            n--;
            m--;
        }//把末尾空格去掉
        k=m+1;//需要用m+1张
        d=k+2;//k+1次多项式需要这么多点
        y[0]=0;
        for(int i=1;i<=d;i++)
            y[i]=(y[i-1]+power(i,k,mod))%mod;
        fac[0]=1;
        for(int i=1;i<=d;i++)
            fac[i]=fac[i-1]*i%mod;
        ans=0;
        for(int i=0;i<=m;i++)
        {
            ans=(ans+Lagrange(n-a[i]))%mod;
            for(int j=i+1;j<=m;j++)
                ans=(ans-power(a[j]-a[i],k,mod)+mod)%mod;
        }
        printf("%lld\n",ans);
    }
    return 0;
}


        
<link rel="stylesheet" href="//cdn.jsdelivr.net/npm/katex/dist/katex.min.css">
        
<link rel="stylesheet" href="//cdn.jsdelivr.net/npm/markdown-it-texmath/css/texmath.min.css">