数论

· · 算法·理论

观前提示：本文为个人对于数论一部分知识点的整理，由于作者能力有限，可能会有错误或者不给出严格证明的情况，敬请谅解，欢迎提出修正，我会尽量修改。文中给出的代码大多为示意核心代码，可能会出现不严谨（如不开 long long）的情况，不保证能通过模板题或相关题目。可以配合 OI Wiki 或者其他资料食用，效果更佳。感谢阅读！

数论函数（算术函数）

数论函数（也称算数函数）指定义域为正整数的函数。数论函数也可以视作一个数列。

积性函数（乘性函数）

定义

若函数 f(n) 满足 f(1)=1，且对于任意的正整数 x,y，若 x,y 互质且满足 f(xy)=f(x)f(y)，则称 f(n) 为积性函数。\ 特别地，若 x,y 不互质时，f(xy)=f(x)f(y) 仍然满足，则称 f(n) 为完全积性函数。

性质

若 f(x) 和 g(x) 均为积性函数，则以下的 h(x) 也为积性函数：

h(x)=f(x^p)\\ h(x)=f^p(x)\\ h(x)=f(x)g(x)\\ h(x)=\sum_{d\mid x}f(d)g(\frac{x}{d})

设正整数 x=\prod p_i^{k_i}，其中 p_i 为素数。

若 F(x) 为积性函数，则有 F(x)=\prod F(p_i^{k_i})。\ 特别地，若 F(x) 为完全积性函数，有 F(x)=\prod F(p_i)^{k_i}。

例子

常数函数：I(n)=1。（完全积性）

恒等函数：id(n)=n。（完全积性）

除数函数：\sigma_k(n)=\sum_{d\mid n} d^k。

幂函数：id_k(n)=n^k。（完全积性）

单位函数：\varepsilon(n)=\begin{cases} 1 & n=1\\0 & n>1 \end{cases}，即 \varepsilon(n)=[n=1]。（完全积性）

以上积性函数常在狄利克雷卷积中用到。

欧拉函数 \varphi(n)。

莫比乌斯函数 \mu(n)。

加性函数

定义

若函数 f(n) 满足 f(1)=1，且对于任意的正整数 x,y，若 x,y 互质且满足 f(xy)=f(x)+f(y)，则称 f(n) 为加性函数。\ 特别地，若 x,y 不互质时，f(xy)=f(x)+f(y) 仍然满足，则称 f(n) 为完全加性函数。

性质

设正整数 x=\prod p_i^{k_i}，其中 p_i 为素数。

若 F(x) 为加性函数，则有 F(x)=\sum F(p_i^{k_i})。\ 特别地，若 F(x) 为完全加性函数，有 F(x)=\sum F(p_i)\times{k_i}。

例子

所有质因子个数 \Omega(n)。（完全加性）

相异质因子个数 \omega(n)。

所有质因子之和 a_0(n)。（完全加性）

相异质因子之和 a_1(n)。

扩展欧几里得（exgcd）

扩展欧几里得算法可以用于求解形如 ax+by=\gcd(a,b) 的方程。

推导过程略，见 OI Wiki，此处仅给出代码。

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

函数的返回值为 \gcd(a,b)。

费马小定理

若 p 为素数，对于任意的 a \in \mathbb Z，若 \gcd(a,p)=1，则有 a^{p-1} \equiv 1\pmod p。

证明是不会的。

逆元

定义

如果 ax\equiv 1 \pmod b，则称 x 是 a\bmod b 的逆元，记作 a^{-1}。

用处

\frac{a}{b}\bmod p=a\times b^{-1}\bmod p

求解

扩展欧几里得法

原式可转化为：ax+by=1。

定理方程 ax+by=m 有解的必要条件是 m \bmod \gcd(a,b)=0。\ 证明： 由最大公因数的定义易得，若 x,y 是整数，则有 ax+by \bmod \gcd(a,b)=0，因此 m\bmod \gcd(a,b)=0。

因此 ax+by=1 有解的必要条件是 1\bmod \gcd(a,b)=0，因此 \gcd(a,b)=1。

所以我们可以直接使用 exgcd 求解。如果要求最小的逆元，需在求完后加上一行 x=(x%b+b)%b。

费马小定理法

由 ax\equiv 1 \pmod b，

若 b 为素数，由费马小定理可得 ax\equiv a^{b-1} \pmod b，

所以 x\equiv a^{b-2} \pmod b。

可以快速幂求解。

注意：费马小定理法只能在 b 为素数时使用！

线性求逆元

一个一个求太慢，考虑推式子然后递推。

（已知 1^{-1} \equiv 1 \pmod p）

如何求 i^{-1} 呢？

显然 p=\lfloor \frac{p}{i}\rfloor \times i+(p \bmod i)，因此有

\lfloor \frac{p}{i}\rfloor\times i+(p \bmod i)\equiv 0 \pmod p \\

式子两边同乘 i^{-1},(p \bmod i)^{-1}，则有

\lfloor \frac{p}{i}\rfloor\times (p \bmod i)^{-1}+i^{-1} \equiv 0 \pmod p\\ i^{-1}\equiv -\lfloor \frac{p}{i}\rfloor\times (p \bmod i)^{-1} \pmod p

这样我们就可以递推了。

核心代码：

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

素数筛

埃氏筛

每次把该数的倍数（不包括自己）标记为合数，这样没有被标记过的就是素数。

code:

bool isprime[N];
long long prime[N];
void Eratosthenes(int n){  //埃氏筛
    isprime[1]=1;
    for(int i=2;i<=n;i++){
        if(!isprime[i]){
            prime[++prime[0]]=i;
            for(long long j=1ll*i*i;j<=n;j+=i){
                isprime[j]=1;
            }
        }
    }
}

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

线性筛（欧拉筛）

考虑优化埃氏筛。显然每个合数会被它的多个质因数重复标记，这样是不优的。让每个数只被其最小的质因子标记一次就可以做到 O(n) 的复杂度了。

code:

bool isprime[N];
long long prime[N];
void Euler(int n){  //欧拉筛 
    isprime[1]=1;
    for(int i=2;i<=n;i++){
        if(!isprime[i]){
            prime[++prime[0]]=i;
        }
        for(int j=1;j<=prime[0]&&i*prime[j]<=n;j++){
            isprime[i*prime[j]]=1;
            if(i%prime[j]==0){
                break;
            }
        }
    }
}

两份代码中，prime[0] 用来记录已经筛出的素数个数。

欧拉函数

定义

定义 \varphi (n)=\sum^n_{i=1} [\gcd(i,n)=1]，即 1,2,3,...,n 中，与 n 互质的数的个数。

其中 [A] 的值：若命题 A 成立为 1，否则为 0。

则 [\gcd(i,n)=1] 的值：若 i,n 互质则为 1，否则为 0。

## 性质 ### 性质一 若 $p$ 为素数，则 $\varphi(p)=p-1$。 显然，$1,2,3,...,p-1$ 都与 $p$ 互质。 ### 性质二 若 $p$ 为素数，则 $\varphi(p^k)=p^k-p^{k-1}=p^k(1-\frac{1}{p})=p^{k-1}(p-1)=p^{k-1}\varphi(p)$。 证明： 显然 $[1,n]$ 中，$p$ 的倍数有 $\lfloor \frac{n}{p} \rfloor$ 个，因此 $[1,p^k]$ 中，与 $p^k$ 不互质的数有 $p^{k-1}$ 个。 ### 性质三 $\varphi(n)$ 为积性函数。 ### 性质四 设正整数 $n=\prod^k_{i=1} p_i^{a_i}$，则有： $\varphi(n)=n\times \prod_{i=1}^k (1-\frac{1}{1-p_i})

证明：

\varphi(n)=\prod_{i=1}^k \varphi(p_i^{a_i})=\prod_{i=1}^k p_i^{a_i}\prod_{i=1}^k(1-\frac{1}{1-p_i})=n\times \prod_{i=1}^k (1-\frac{1}{1-p_i})

解释一下：第一步变换是根据性质三，第二步根据性质二与性质三，第三步注意到式子的前半部分就是 n，证毕。

性质五

做题时经常用到的一个结论。疑似叫欧拉反演。

n=\sum_{d\mid n} \varphi(d)

证明：\ 假设我们有 n 个分数 \frac{1}{n},\frac{2}{n},...,\frac{n-1}{n},\frac{n}{n}，将它们能约分的全部约分，那么约分后的各个分数的分母 d 显然满足 d\mid n，而分母为 d 的分数共出现 \varphi(d) 次。因此 \sum_{d\mid n} \varphi(d)=n。

更常用的是把 n 替换为 \gcd(x,y)，那么我们就得到了下面的式子：

\gcd(x,y)=\sum_{d\mid x,d\mid y} \varphi(d)

求欧拉函数

单个

利用性质四，枚举 n 的质因数即可。

code:

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

时间复杂度 O(\sqrt n)。

多个

如何求解 [1,n] 中各个数的欧拉函数值？

已知我们已经求得了 \varphi(i) 和 \varphi(p) 的值，其中 p 为素数，则有：

\varphi(i) \times p & i \bmod p =0 \\ \varphi(i) \times \varphi(p) & i \bmod p \ne 0 \end{cases}

利用性质三不难证明。

由此我们可以在线性筛的基础上写出求欧拉函数的代码：

bool isprime[N];
long long prime[N],phi[N];
void Euler(int n){
    isprime[1]=1;
    phi[1]=1;
    for(int i=2;i<=n;i++){
        if(!isprime[i]){
            prime[++prime[0]]=i;
            phi[i]=i-1;
        }
        for(int j=1;j<=prime[0]&&i*prime[j]<=n;j++){
            isprime[i*prime[j]]=1;
            if(i%prime[j]==0){
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            else{
                phi[i*prime[j]]=phi[i]*phi[prime[j]];
            }
        }
    }
}

莫比乌斯反演

莫比乌斯函数

设 f 为算术函数，F 为 f 的和函数，F(n)=\sum _{d\mid n} f(d)，其中 d\mid n 表示 d 整除 n，即 d 是 n 的因数。

如何用 F 求出 f 的值？

打表不难发现，f(n) 等于一些 F(\frac{n}{d}) 的和，而这些项的系数为 -1 或 0 或 1。猜想有以下式子：

f(n)=\sum_{d\mid n} \mu(d)F(\frac{n}{d})

其中 \mu 是算数函数。

若式子成立，不难发现，若 p 为素数，因为 f(p)=F(p)-F(1)，可得 \mu(p)=-1。

因为 f(p^2)=F(p^2)-F(p)，可得 \mu(p^2)=0，同理可得 \mu(p^k)=0。

定义

1 & n=1 \\ 0 & n \text 含有平方因子 \\ (-1)^r & r \text 为 n \text 本质不同质因子个数 \end{cases}

也就是说，如果 n=1，则 \mu(n)=1；否则只要 n 有任意质因子次数 >1，\mu(n)=0（证明见上）；否则 \mu(n)=(-1)^{\text{n的质因子个数}}。

性质

性质一

莫比乌斯函数 \mu(n) 是积性函数。

性质二

1 & n=1 \\ 0 & n>1 \end{cases}

即 \sum_{d\mid n} \mu(d)=[n=1]。\ 有关 [n=1]，见欧拉函数定义。

证明：

当 n=1 时，显然。

当 n>1 时，设 n=\prod^k_{i=1} p_i^{a_i}，则有 F(n)=\prod^k_{i=1} F(p_i^{a_i})，设其中一项为 F(p^t)，则\

### 性质三 由莫比乌斯函数的由来，可以得出莫比乌斯反演公式： 若 $f$ 为算术函数，$F$ 为 $f$ 的和函数，对于任意 $n \in \mathbb N^ \times$，$F(n)=\sum _{d\mid n} f(d)$，则有\ $f(n)=\sum_{d\mid n} \mu(d)F(\frac{n}{d})

重要结论

由性质二，将 n 替换为 \gcd(i,j)，可得：

[\gcd(i,j)=1]=\sum_{d\mid \gcd(i,j)} \mu(d)

有关 [\gcd(i,j)=1]，见欧拉函数定义。

求莫比乌斯函数

在线性筛的基础上修改，与欧拉函数类似。原理见定义与性质。

bool isprime[N];
long long prime[N],mu[N];
void Euler(int n){
    isprime[1]=1;
    mu[1]=1;
    for(int i=2;i<=n;i++){
        if(!isprime[i]){
            prime[++prime[0]]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=prime[0]&&i*prime[j]<=n;j++){
            isprime[i*prime[j]]=1;
            if(i%prime[j]==0){
                mu[i*prime[j]]=0;
                break;
            }
            else{
                mu[i*prime[j]]=-mu[i];
            }
        }
    }
}

数论分块（整除分块）

用于计算含有下取整的和式，如 \sum_{i=1}^n f(i)\lfloor \frac{n}{i}\rfloor。

对于 \sum_{i=1}^n \lfloor \frac{n}{i}\rfloor，如何更快速地计算？

对于小数据打表，不难发现，\lfloor \frac{n}{i}\rfloor 的值在一段区间内相同。

以 n=18 为例，\lfloor \frac{n}{i}\rfloor 的值如下：

18 9 6 4 3 3 2 2 2 1 1 1 1 1 1 1 1 1

那么我们可以想办法得出每一块中的值和每一块的长度，进而在更快时间内计算。

事实上，有如下结论：

结论

值 \lfloor \frac{n}{i}\rfloor 所在块的右端点为 \Large \lfloor \frac{n}{\lfloor \frac{n}{i}\rfloor}\rfloor。

证明：

令 k=\lfloor \frac{n}{i}\rfloor，可知 k\le \frac{n}{i}。

因此 \lfloor \frac{n}{k}\rfloor \ge\lfloor \frac{n}{\frac{n}{i}}\rfloor=\lfloor i\rfloor=i。

所以满足条件的 i_{max}=\large \lfloor \frac{n}{k}\rfloor=\Large \lfloor \frac{n}{\lfloor \frac{n}{i}\rfloor}\rfloor。

知道了块的左右端点，也就是知道了块长，我们前缀和预处理 f(i)，就可以求出 \sum_{i=1}^n f(i)\lfloor \frac{n}{i}\rfloor 的值了。

代码实现

for(int l=1,r;l<=n;l=r+1){
    r=n/(n/l);
    ans+=(sum[r]-sum[l-1])*(n/l);  //sum[]是预处理的前缀和
}

时间复杂度是 O(\sqrt n) 的，但是我不会证。

用途

一个重要用途是和莫比乌斯反演结合以进一步降低复杂度。这里给出一道两者结合的经典题目：[POI 2007] ZAP-Queries，大致思路是利用莫比乌斯反演中的重要结论得出和式，再利用整除分块在最后运算时加速。

狄利克雷卷积（Dirichlet 卷积）

定义

对于两个数论函数 f(x) 和 g(x)，它们的狄利克雷卷积得到的结果 h(x) 定义为：

h(x)=(f \times g)(n)=\sum_{d\mid n}f(d)g(\frac{x}{d})

上式可简记为：

h=f \times g

性质

类似于普通乘法，有以下性质：

交换律 f \times g=g \times f。

结合律 (f \times g) \times h=f \times (g \times h)。

分配律 (f+g) \times h=f \times h+g \times h。

等式的性质 f=g 的充要条件是 f \times h=g \times h，其中 h(x) 需满足 h(1)\ne 0。

单位元单位函数 \varepsilon 是狄利克雷卷积运算中的单位元，即对于任何数论函数 f，都有 f \times \varepsilon=f。

逆元对于任一满足 f(x)\ne 0 的数论函数，如果有另一个数论函数 g(x) 满足 f \times g=\varepsilon，则称 g(x) 是 f(x) 的逆元。由等式的性质可知，逆元是唯一的。

几个等式

  1. \mu \times I=\varepsilon
  2. \varphi \times I=id
  3. \mu \times id=\varphi

重要结论

srds，我不知道这俩结论为啥重要，也不知道咋证，所以只给个结论很合理吧。

结论一

两个积性函数的 Dirichlet 卷积也是积性函数。

结论二

积性函数的逆元也是积性函数。

证明好像可以归纳法？？？

杜教筛

推导

已知一个积性函数 f，现在要求它的前缀和，且复杂度要小于线性。

我们设 S(n)=\sum_{i=1}^n f(i)，再找一个积性函数 g，考虑 (f \times g) 的前缀和，则有：

\sum_{i=1}^n (f \times g)(i) \begin{aligned} &= \sum_{i=1}^n \sum _{d\mid i} f(d)g(\frac{i}{d})\\ &= \sum_{d=1}^n g(d)\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor} f(i)\\ &= \sum_{d=1}^n g(d)S(\lfloor \frac{n}{d}\rfloor) \end{aligned}

注意到，

\begin{aligned} g(1)S(n)&=\sum_{i=1}^n g(i)S(\lfloor \frac{n}{d}\rfloor)-\sum_{i=2}^n g(i)S(\lfloor \frac{n}{i}\rfloor)\\ &=\sum_{i=1}^n (f \times g)(i)-\sum_{i=2}^n g(i)S(\lfloor \frac{n}{i}\rfloor) \end{aligned}

最后这个就是杜教筛的核心式子。

因此如果我们能找到一个合适的积性函数 g，能够快速算出 g 的前缀和以及 \sum_{i=1}^n (f \times g)(i) 的值，就可以利用数论分块加上递归更快地求解 S(n)。

这么做的时间复杂度是 O(n^{\frac{3}{4}}) 的，但如果我们能线性筛预处理出前 n^{\frac{2}{3}} 个答案，再利用杜教筛求解，就可以做到 O(n^{\frac{2}{3}}) 的复杂度。

应用

\varphi 的前缀和

我们知道 \varphi \times I=id，其中 id 的前缀和为 \frac{n\times(n+1)}{2}，因此显然有下列代码：

int get_S_phi(int n){
    if(n<=N){
        return S_phi_init[n];
    }
    if(S_phi[n]){
        return S_phi[n];
    }
    int ans=n*(n+1)/2;
    for(int l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        ans-=(r-l+1)*get_S_phi(n/l);
    }
    return S_phi[n]=ans;  //记忆化
}

\mu 的前缀和

同理，由 \mu \times I=\varepsilon，可得下列代码：

int get_S_mu(int n){
    if(n<=N){
        return S_mu_init[n];
    }
    if(S_mu[n]){
        return S_mu[n];
    }
    int ans=1;
    for(int l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        ans-=(r-l+1)*get_S_mu(n/l);
    }
    return S_mu[n]=ans;
}

这里一并给出 P4213 【模板】杜教筛 的全部代码：

#include<iostream>
#include<unordered_map>
using namespace std;
#define int long long
const int N=3e6;
bool isprime[N+10];
int prime[N+10],phi[N+10],mu[N+10];
void Euler(int n){
    isprime[1]=1;
    phi[1]=mu[1]=1;
    for(int i=2;i<=n;i++){
        if(!isprime[i]){
            prime[++prime[0]]=i;
            phi[i]=i-1;
            mu[i]=-1;
        }
        for(int j=1;j<=prime[0]&&i*prime[j]<=n;j++){
            isprime[i*prime[j]]=1;
            if(i%prime[j]==0){
                phi[i*prime[j]]=phi[i]*prime[j];
                mu[i*prime[j]]=0;
                break;
            }
            else{
                phi[i*prime[j]]=phi[i]*phi[prime[j]];
                mu[i*prime[j]]=-mu[i];
            }
        }
    }
}
int S_phi_init[N+10],S_mu_init[N+10];
void init(){
    Euler(N);
    for(int i=1;i<=N;i++){
        S_phi_init[i]=S_phi_init[i-1]+phi[i];
        S_mu_init[i]=S_mu_init[i-1]+mu[i];
    }
}
unordered_map<int,int> S_phi,S_mu;
int get_S_phi(int n){
    if(n<=N){
        return S_phi_init[n];
    }
    if(S_phi[n]){
        return S_phi[n];
    }
    int ans=n*(n+1)/2;
    for(int l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        ans-=(r-l+1)*get_S_phi(n/l);
    }
    return S_phi[n]=ans;
}
int get_S_mu(int n){
    if(n<=N){
        return S_mu_init[n];
    }
    if(S_mu[n]){
        return S_mu[n];
    }
    int ans=1;
    for(int l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        ans-=(r-l+1)*get_S_mu(n/l);
    }
    return S_mu[n]=ans;
}
signed main(){
    init();
    int T;
    cin>>T;
    while(T--){
        int n;
        cin>>n;
        cout<<get_S_phi(n)<<" "<<get_S_mu(n)<<"\n";
    }
    return 0;
}

卢卡斯定理 / Lucas 定理

定义

定义 Lucas(n,m,p)=C^n_m \bmod p，

则有 Lucas(n,m,p)=Lucas(n/p,m/p,p)\times C^{n\bmod p}_{m\bmod p}。

其中 p 为素数。

代码实现

递归即可。

int Lucas(int n,int m,int p){
    if(!m){
        return 1;
    }
    return Lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}

扩展欧拉定理

定义

a^b & b<\varphi(m)\\ a^{(b\bmod \varphi(m))+\varphi(m)} & b\ge \varphi(m) \end{cases} \pmod m

证明（不严谨）

不难想到，a^i \bmod m 的值会是一个循环，我们知道了循环节的长度就可以降幂，进而更快速地求解了。

代码实现

我直接略略略略略。

仅给出示意核心代码。

cout<<qpow(a,(b<phi(p))?b:(b%phi(p)+phi(p)),p);

中国剩余定理（CRT）

定义

中国剩余定理可以用于求解如下形式关于 x 的一元线性同余方程组（其中 m_1,m_2,...,m_k 两两互质）：

\begin{cases} x \equiv a_1\pmod {m_1}\\ x \equiv a_2\pmod {m_2}\\ ...\\ x \equiv a_k\pmod {m_k} \end{cases}

过程

  1. 计算所有模数的积 M；
  2. 遍历每个方程，对于第 i 个方程：\ a. 计算 m\_now=\frac{M}{m_i}；\ b. 计算 m\_now 在模 m_i 意义下的逆元 m\_now^{-1}（使用 exgcd 算法）；\ c. 计算 c_i=m\_now\times m\_now^{-1}；
  3. 方程组在模 M 意义下的唯一解为：x=\sum^k_{i=1} a_ic_i \pmod M。

代码实现

int CRT(){
    int ans=0;
    for(int i=1;i<=n;i++){
        int m_now=M/m[i],x,y;
        exgcd(m_now,m[i],x,y);
        ans=((ans+m_now*x*a[i])%M+M)%M;
    }
    return ans;
}