「算法笔记」杜教筛

2021 年写的（已折叠）

一、前置概念

具体在 「算法笔记」莫比乌斯反演 写过，所以「前置概念」就简单写写。积性函数和完全积性函数就不写了。

狄利克雷卷积：对于两个数论函数 \(f,g\)，定义它们的狄利克雷卷积 \(h=f*g\) 为

\[\displaystyle h(x)=\sum_{d\mid x}f(d)g\left(\frac{x}{d}\right) \]

狄利克雷卷积满足交换律、结合律、对加法的分配律，有单位元 \(\varepsilon\)（其中 \(\varepsilon\) 为单位函数 \(\varepsilon(x)=[x=1]\)）。若 \(f,g\) 是积性函数，则 \(f*g\) 也是积性函数。

在狄利克雷卷积的意义下，\(\mu*1=\varepsilon\)，即 \(\mu\) 和 \(1\) 互为逆元。

常用卷积：

• \(\varepsilon=\mu*1\)，即 \(\varepsilon(n)=\sum_{d\mid n}\mu(d)\)。
• \(\varphi*1=\text{Id}\)，即 \(\sum_{d\mid n}\varphi(d)=n\)。两边同时卷一个 \(\mu\) 可得 \(\varphi=\text{Id}*\mu\)，即 \(\varphi(n)=\sum_{d\mid n}d\cdot \mu(\frac n d)\)。
• \(d=1*1 \Leftrightarrow d(n)=\sum_{d\mid n}1\)，\(\sigma_k=\text{Id}_k*1 \Leftrightarrow \sigma_k(n)=\sum_{d|n} d^k\)

提共因式：对于任意数论函数 \(f,g\) 和 完全积性函数 \(h\) 均有 \((f*g)\cdot h=(f\cdot h)*(g\cdot h)\)（注意中间是点乘，\((f\cdot g)(i)=f(i)g(i)\)），展开即可证。比如 \((\varphi\cdot \text{Id}_k)* (1\cdot \text{Id}_k)=(\varphi*1)\cdot \text{Id}_k=\text{Id}\cdot \text{Id}_k=\text{Id}_{k+1}\)。

二、杜教筛

1. 基本思想

对于数论函数 \(f\)，求 \(S(n)=\sum_{i=1}^nf(i)\)。

假如我们能快速计算 \(g\) 和 \(f*g\) 的前缀和，考虑怎么求 \(S\)。

\[\sum_{i=1}^n (f*g)(i)=\sum_{i=1}^n\sum_{d\mid i}f(\frac i d)g(d)=\sum_{d=1}^n g(d)\sum_{i=1}^{\lfloor\frac n i\rfloor}f(\frac{id}{d})=\sum_{d=1}^n g(d)\sum_{i=1}^{\lfloor\frac n i\rfloor}f(i)=\sum_{i=1}^n g(i)S(\lfloor\frac n i\rfloor) \]

则可以得到 \(S(n)\) 关于 \(S(\lfloor \frac{n}{i}\rfloor)\) 的递推式：

\[g(1)S(n)=\sum_{i=1}^n(f*g)(i)-\sum_{i=2}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \]

如果我们可以快速求出 \(\sum_{i=1}^n(f*g)(i)\)，再用整除分块来求 \(\sum_{i=2}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\)，就能求得 \(g(1)S(n)\)。也就是说，我们要找到一个合适的 \(g\)，使得 \(f*g\) 和 \(g\) 都能快速求前缀和。

\(S(\lfloor\frac{n}{i}\rfloor)\) 是一个子问题，对 \(S(n)\) 递归求解并记忆化即可。

线性筛预处理前 \(n^{\frac{2}{3}}\) 个 \(f(x)\) 的值，则时间复杂度为 \(\mathcal O(n^{\frac{2}{3}})\)。

（当然如果我们能快速求出 \(f,g\) 的前缀和，也可以求 \(f*g\) 的前缀和，根据最初那个式子整除分块右边推左边就好了）

2. 复杂度证明

先不考虑预处理。首先我们每次记忆化搜索用到的 \(S\) 可以写成 \(S(\lfloor\large\frac n k\normalsize\rfloor)\) 的形式，再加上每次整除分块用的 \(\sqrt x\)，总复杂度就是 \(\sum_{\large x,\exists k, \lfloor\frac n k\rfloor=x}\sqrt x\)。

• Lemma：\(\lfloor\sqrt x\rfloor=p\)，那么 \(x\in[p^2,(p+1)^2)\)，即 \(x\) 有 \((p+1)^2-p^2=\mathcal O(p)\) 个不同的取值（下面的根号为向下取整）。

• 若 \(x=\lfloor\large\frac n k\normalsize\rfloor\leq \sqrt n\)，对和式的贡献为 \(\sum\limits_{x=1}^{\sqrt n}\sqrt x\)，枚举 \(i=\sqrt x\)，根据 Lemma 有 \(\sum\limits_{i=1}^{n^{\frac 1 4}}i^2\)（有 \(\mathcal O(i)\) 个 \(\sqrt x=i\) 的 \(x\)，每个贡献 \(i\)），平方求和式个三次方的东西所以有 \(\mathcal O((n^\frac 1 4)^3)=\mathcal O(n^{\frac 3 4})\)。

• 若 \(x=\lfloor\large\frac n k\normalsize\rfloor>\sqrt n\)，那么对应的 \(k\leq \sqrt n\)，对和式的贡献为 \(\sum\limits_{k=1}^{\sqrt n}\sqrt{\large\frac n k\normalsize }=\sqrt n\sum\limits_{k=1}^{\sqrt n}\frac{1}{\sqrt k}\)，枚举 \(i=\sqrt k\)，根据 Lemma 有 \(\sqrt n\sum\limits_{i=1}^{n^{\frac 1 4}}\frac{i}{i}=n^{\frac 3 4}\)。

因此总复杂度也是 \(\mathcal O(n^{\frac 3 4})\) 级别的。还可以进一步优化杜教筛，即先线性筛出前 \(m\) 个答案，之后再用杜教筛。假设 \(m>\sqrt n\)，优化后复杂度变为 \(\mathcal O(m+\sum\limits_{k=1}^{\large\lfloor\frac n m\rfloor}\sqrt{\large\frac n k})=\mathcal O(m+\sqrt n\cdot \sqrt{\large\frac n m})\)，根据均值不等式当 \(m=\sqrt n\cdot \sqrt{\large\frac n m}\) 时复杂度最优，此时 \(m=n^{\frac 2 3}\)，时间复杂度达到最小值 \(\mathcal O(n^{\frac 2 3})\)。

实测 1s 杜教筛能跑 \(n=10^{10}\)。

upd on 2022.1.14：

pwj 指正：并不能按 \(m=\sqrt n\cdot \sqrt{\large\frac n m}=\large\frac{n}{\sqrt m}\) 算，因为 \(m+\large\frac{n}{\sqrt m}\normalsize\geq 2\sqrt{m\cdot \large\frac{n}{\sqrt m}}\) 会发现右边的 \(m\) 消不掉，而 \(m\) 是我们自己定的。但是数量级没错。

这里有个拆项凑系数用来证不等式的 trick：

\[m+\frac{n}{\sqrt m}=m+\frac{\frac n 2}{\sqrt m}+\frac{\frac n 2}{\sqrt m}\geq 3\sqrt[3]{m\cdot \frac{\frac n 2}{\sqrt m}\cdot \frac{\frac n 2}{\sqrt m}}=3\sqrt[3]{\frac{n^2}{4}}=\mathcal O(n^{\frac 2 3}) \]

当 \(m=\large\frac{\frac n 2}{\sqrt m}\) 时复杂度最优。

3. 例子

求 \(S(n)=\sum_{i=1}^n\varphi(i)\)。

\(\varphi*1=\text{Id}\)。显然 \(1\)、\(\text{Id}(x)=x\) 都可以快速求前缀和。取 \(g(x)=1\) 即可。

\(S(n)=\sum_{i=1}^n i-\sum_{i=2}^nS\left(\left\lfloor\frac{n}{i}\right\rfloor\right)=\frac{n(n+1)}{2}-\sum_{i=2}^nS\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\)。

求 \(S(n)=\sum_{i=1}^n\mu(i)\)。

\(\mu*1=\varepsilon\)。显然 \(1\) 和 \(\varepsilon(x)=[x=1]\) 可以快速求前缀和。取 \(g(x)=1\) 即可。

\(S(n)=1-\sum_{i=2}^nS\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\)。

//Luogu P4213
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=2e6+5;
int t,n=2e6,cnt,p[N],phi[N],mu[N];
bool vis[N];
map<int,int>Smu,Sphi;
int S_mu(int n){
    if(n<=2e6) return mu[n];
    if(Smu[n]) return Smu[n];
    int ans=0;
    for(int l=2,r=0;l<=n;l=r+1)    //注意从 2 开始
        r=n/(n/l),ans+=(r-l+1)*S_mu(n/l); 
    return Smu[n]=1-ans;
}
int S_phi(int n){
    if(n<=2e6) return phi[n];
    if(Sphi[n]) return Sphi[n];
    int ans=0;
    for(int l=2,r=0;l<=n;l=r+1)
        r=n/(n/l),ans+=(r-l+1)*S_phi(n/l);
    return Sphi[n]=n*(n+1)/2-ans; 
}
signed main(){
    scanf("%lld",&t);
    vis[0]=vis[1]=1,phi[1]=mu[1]=1;
    for(int i=2;i<=n;i++){
        if(!vis[i]) p[++cnt]=i,phi[i]=i-1,mu[i]=-1; 
        for(int j=1;j<=cnt&&i*p[j]<=n;j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j],mu[i*p[j]]=0;break;}
            phi[i*p[j]]=phi[i]*phi[p[j]],mu[i*p[j]]=-mu[i]; 
        }
    }
    for(int i=1;i<=n;i++) mu[i]+=mu[i-1],phi[i]+=phi[i-1];
    while(t--){
        scanf("%lld",&n);
        printf("%lld %lld\n",S_phi(n),S_mu(n));
    }
    return 0;
} 

再比如说求 \(S(n)=\sum_{i=1}^n i\cdot \varphi(i)\)。

很明显 \(f(x)=x\cdot\varphi(x)\) 是积性函数。我们尝试构造出合适的 \(g(x)\)。

\[(f*g)(x)=\sum_{d\mid x}d\cdot \varphi(d)\cdot g(\frac x d) \]

考虑将 \(d\) 消掉，这样可以利用 \(\sum_{d\mid x}\varphi(d)=x\) 推出一些东西。尝试 \(g=\text{Id}\)。

\[(f*g)(x)=\sum_{d\mid x}d\cdot \varphi(d)\cdot \frac{x}{d}=x\sum_{d\mid x}\varphi(d)=x^2 \]

取 \(g=\text{Id}\) 即可。

求 \(S(n)=\sum_{i=1}^n\varphi(i)\cdot i^k\)。

按照提公因式卷上一个 \(1\cdot \text{Id}_k\)，那么 \((\varphi\cdot \text{Id}_k)*(1\cdot\text{Id}_k)=(\varphi*1)\cdot \text{Id}_k=\text{Id}\cdot\text{Id}_k=\text{Id}_{k+1}\)。

取 \(g=\text{Id}_k\) 即可。

求 \(S(n)=\sum_{i=1}^n \mu(i)\cdot i^k\)。

按照提公因式卷上一个 \(1\cdot \text{Id}_k\)，那么 \((\mu\cdot \text{Id}_k)*(1\cdot \text{Id}_k)=(\mu*1)\cdot \text{Id}_k=\varepsilon\cdot \text{Id}_k=\varepsilon\)。取 \(g=\text{Id}_k\) 即可。

三、狄利克雷生成函数

有些时候不容易看出 \(g(x)\) 取什么，比如：求 \(S(n)=\sum_{i=1}^n i\cdot \varphi(i)\)。凑是一种方法；而使用 狄利克雷生成函数（DGF） 可以从另一个角度求出 \(g(x)\)。

对于无穷序列 \(f_1,f_2,\cdots\)，定义其狄利克雷生成函数为

\[\tilde F(x)=\sum_{i\geq 1}\frac{f_i}{i^x} \]

与普通生成函数以及指数生成函数的对比：

普通生成函数：\(F(x)=\sum_{i\geq 0}f_ix^i\)。

指数生成函数：\(\hat F(x)=\sum_{i\geq 0}f_i\frac{x^i}{i!}\)。

对于两个序列 \(f,g\)，其 DGF 的乘积对应 \(f,g\) 的狄利克雷卷积序列：

\[\tilde F(x)\tilde G(x)=\sum_{i\geq 1}\sum_{j\geq 1}\frac{f_i}{i^x}\frac{g_j}{j^x}=\sum_{i\geq 1}\sum_{j\geq 1}\frac{f_ig_j}{(ij)^x}=\sum_d\sum_{i\mid d}\frac{f_ig_{\frac{d}{i}}}{d^x} \]

如果序列 \(f\) 满足积性（积性函数），其 DGF 可以由质数幂处的取值表示：

\[\tilde F=\prod_{p\in\text{prime}}(1+\frac{f(p)}{p^x}+\frac{f(p^2)}{p^{2x}}+\frac{f(p^3)}{p^{3x}}+\cdots) \]

可以考虑 \(i=\prod p_j^{e_j}\)，\(f(i)=\prod f(p_j^{e_j})\)。那么有：

\[\frac{f_i}{i^x}=\prod \frac{f(p_j^{e_j})}{(p_j^{e_j})^x} \]

咕咕咕……

四、例题

Luogu P3768 简单的数学题

Problem：给定 \(n,p\)，求：

\[\left(\sum_{i=1}^n\sum_{j=1}^n i\cdot j\cdot \gcd(i,j)\right)\bmod p \]

\(n\leq 10^{10},5\times 10^8\leq p\leq 1.1\times 10^9\)，\(p\) 是质数。

Solution：

可以用 \(\varphi*1=\text{Id}\) 推式子，即 \(\sum_{d\mid x}\varphi(d)=x\)。

\[ans=\sum_{i=1}^n\sum_{j=1}^n i\cdot j\cdot \sum_{d\mid \gcd(i,j)}\varphi(d)=\sum_{d=1}^n \varphi(d)d^2(\sum_{i=1}^{\lfloor\frac n d\rfloor}i)^2 \]

设 \(s(n)=(\sum_{i=1}^n i)^2=(\frac{(1+n)n}{2})^2\)，则 \(ans=\sum_{d=1}^n\varphi(d)d^2s(\lfloor\frac{n}{d}\rfloor)\)，如果能快速求出 \(f(x)=\varphi(x)x^2\) 的前缀和，那么整除分块就做完了。

考虑用杜教筛来筛 \(f(x)=\varphi(x)x^2\)，要找到一个合适的 \(g\)，使得 \(f*g\) 和 \(g\) 都可以快速求前缀和。

\[(f*g)(x)=\sum_{d\mid x}\varphi(d)d^2g(\frac x d) \]

考虑把这里的 \(d^2\) 消掉，我们尝试 \(g=\text{Id}_2\)（幂函数 \(\text{Id}_k(n)=n^k\)），那么：

\[(f*g)(x)=\sum_{d\mid x}\varphi(d)d^2(\frac x d)^2=\left(\sum_{d\mid x}\varphi(d)\right)x^2=x\cdot x^2=x^3 \]

取 \(g=\text{Id}_2\) 即可，\(f*g=\text{Id}_3\)。当然也可以用 狄利克雷生成函数（DGF） 推出。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=2e6+5;
int n,m=2e6,mod,cnt,p[N],phi[N],inv2,inv6,ans;
bool vis[N];
map<int,int>s;
int mul(int x,int n,int mod){
    int ans=mod!=1;
    for(x%=mod;n;n>>=1,x=x*x%mod)
        if(n&1) ans=ans*x%mod;
    return ans;
}
int s2(int n){    //平方和 
    n%=mod;
    return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
}
int s3(int n){    //立方和 
    n%=mod; 
    int x=n%mod*(n+1)%mod*inv2%mod;
    return x*x%mod;
}
int S(int n){
    if(n<=m) return phi[n];
    if(s[n]) return s[n];
    int ans=0;
    for(int l=2,r=0;l<=n;l=r+1)
        r=n/(n/l),ans=(ans+(s2(r)-s2(l-1))%mod*S(n/l)%mod)%mod;
    return s[n]=((s3(n)-ans)%mod+mod)%mod;
}
signed main(){
    scanf("%lld%lld",&mod,&n);
    vis[0]=vis[1]=1,phi[1]=1;
    for(int i=2;i<=m;i++){
        if(!vis[i]) p[++cnt]=i,phi[i]=i-1;
        for(int j=1;j<=cnt&&i*p[j]<=m;j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j];break;}
            phi[i*p[j]]=phi[i]*phi[p[j]]; 
        }
    }
    for(int i=1;i<=m;i++) phi[i]=(phi[i-1]+i*i%mod*phi[i]%mod)%mod;
    inv2=mul(2,mod-2,mod),inv6=mul(6,mod-2,mod);
    for(int l=1,r=0;l<=n;l=r+1)
        r=n/(n/l),ans=(ans+s3(n/l)*(S(r)-S(l-1))%mod)%mod;
    printf("%lld\n",(ans+mod)%mod);
    return 0;
} 

一、杜教筛

1. 基本思想

对于数论函数 \(f\)，求 \(S(n)=\sum_{i=1}^nf(i)\)。

\[\sum_{i=1}^n (f*g)(i)=\sum_{i=1}^n\sum_{d\mid i}f(\frac i d)g(d)=\sum_{d=1}^n g(d)\sum_{i=1}^{\lfloor\frac n i\rfloor}f(\frac{id}{d})=\sum_{i=1}^n g(i)S(\left\lfloor\frac n i\right\rfloor)\\ g(1)S(n)=\sum_{i=1}^n(f*g)(i)-\sum_{i=2}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \]

找到一个合适的 \(g\)，使得 \(f*g\) 和 \(g\) 都能快速求前缀和。整除分块求 \(\sum_{i=2}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\)。

对 \(S\) 递归求解并记忆化。线性筛预处理前 \(n^{\frac{2}{3}}\) 个 \(f(x)\) 的值，时间复杂度 \(\mathcal O(n^{\frac{2}{3}})\)。

2. 例子

求 \(S(n)=\sum_{i=1}^n\varphi(i)\)。

\(\varphi*1=\text{Id}\)。显然 \(1\)、\(\text{Id}(x)=x\) 都可以快速求前缀和。取 \(g(x)=1\) 即可。

\(S(n)=\sum_{i=1}^n i-\sum_{i=2}^nS\left(\left\lfloor\frac{n}{i}\right\rfloor\right)=\frac{n(n+1)}{2}-\sum_{i=2}^nS\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\)。

求 \(S(n)=\sum_{i=1}^n\mu(i)\)。

\(\mu*1=\varepsilon\)。显然 \(1\) 和 \(\varepsilon(x)=[x=1]\) 可以快速求前缀和。取 \(g(x)=1\) 即可。

\(S(n)=1-\sum_{i=2}^nS\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\)。

//Luogu P4213
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=2e6+5;
int t,n=2e6,cnt,p[N],phi[N],mu[N];
bool vis[N];
map<int,int>Smu,Sphi;
int S_mu(int n){
	if(n<=2e6) return mu[n];
	if(Smu[n]) return Smu[n];
	int ans=0;
	for(int l=2,r=0;l<=n;l=r+1)	//注意从 2 开始
		r=n/(n/l),ans+=(r-l+1)*S_mu(n/l); 
	return Smu[n]=1-ans;
}
int S_phi(int n){
	if(n<=2e6) return phi[n];
	if(Sphi[n]) return Sphi[n];
	int ans=0;
	for(int l=2,r=0;l<=n;l=r+1)
		r=n/(n/l),ans+=(r-l+1)*S_phi(n/l);
	return Sphi[n]=n*(n+1)/2-ans; 
}
signed main(){
	scanf("%lld",&t);
	vis[0]=vis[1]=1,phi[1]=mu[1]=1;
	for(int i=2;i<=n;i++){
		if(!vis[i]) p[++cnt]=i,phi[i]=i-1,mu[i]=-1; 
		for(int j=1;j<=cnt&&i*p[j]<=n;j++){
			vis[i*p[j]]=1;
			if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j],mu[i*p[j]]=0;break;}
			phi[i*p[j]]=phi[i]*phi[p[j]],mu[i*p[j]]=-mu[i]; 
		}
	}
	for(int i=1;i<=n;i++) mu[i]+=mu[i-1],phi[i]+=phi[i-1];
	while(t--){
		scanf("%lld",&n);
		printf("%lld %lld\n",S_phi(n),S_mu(n));
	}
	return 0;
} 

再比如说求 \(S(n)=\sum_{i=1}^n i\cdot \varphi(i)\)。

很明显 \(f(x)=x\cdot\varphi(x)\) 是积性函数。我们尝试构造出合适的 \(g(x)\)。

\[(f*g)(x)=\sum_{d\mid x}d\cdot \varphi(d)\cdot g(\frac x d) \]

考虑将 \(d\) 消掉，这样可以利用 \(\sum_{d\mid x}\varphi(d)=x\) 推出一些东西。尝试 \(g=\text{Id}\)。

\[(f*g)(x)=\sum_{d\mid x}d\cdot \varphi(d)\cdot \frac{x}{d}=x\sum_{d\mid x}\varphi(d)=x^2 \]

取 \(g=\text{Id}\) 即可。

求 \(S(n)=\sum_{i=1}^n \varphi(i)\cdot i^2\)：

\[(f*g)(x)=\sum_{d\mid x}\varphi(d)d^2g(\frac x d) \]

考虑把这里的 \(d^2\) 消掉，我们尝试 \(g=\text{Id}_2\)，那么：

\[(f*g)(x)=\sum_{d\mid x}\varphi(d)d^2(\frac x d)^2=\left(\sum_{d\mid x}\varphi(d)\right)x^2=x\cdot x^2=x^3 \]

取 \(g=\text{Id}_2\) 即可，\(f*g=\text{Id}_3\)。

求 \(S(n)=\sum_{i=1}^n\varphi(i)\cdot i^k\)。

按照提公因式卷上一个 \(1\cdot \text{Id}_k\)，那么 \((\varphi\cdot \text{Id}_k)*(1\cdot\text{Id}_k)=(\varphi*1)\cdot \text{Id}_k=\text{Id}\cdot\text{Id}_k=\text{Id}_{k+1}\)。

取 \(g=\text{Id}_k\) 即可。

求 \(S(n)=\sum_{i=1}^n \mu(i)\cdot i^k\)。

按照提公因式卷上一个 \(1\cdot \text{Id}_k\)，那么 \((\mu\cdot \text{Id}_k)*(1\cdot \text{Id}_k)=(\mu*1)\cdot \text{Id}_k=\varepsilon\cdot \text{Id}_k=\varepsilon\)。

取 \(g=\text{Id}_k\) 即可。

\(f=\mu^2*(\mu\cdot \text{Id})\)，求 \(S_f(n)\)。

卷 \(1\cdot\text{Id}\)，那么 \(\mu^2*(\mu\cdot \text{Id})*(1\cdot \text{Id})=\mu^2*((\mu*1)\cdot \text{Id})=\mu^2*\varepsilon=\mu^2\)。其中 \(\sum_{i=1}^n \mu^2(i)=\sum_{i=1}^{\sqrt n}\lfloor\large\frac{n}{i^2}\normalsize\rfloor\mu(i)\)。

取 \(g=\text{Id}\) 即可。

3. 补充

如果 \(g,f*g\) 是 \(1,\text{Id}_k,\varepsilon\)，我们可以 \(\mathcal O(1)/\mathcal O(\log n)\) 求出前缀和（由于复杂度瓶颈在于整除分块的根号，即便这里是 \(\log n\) 也不影响复杂度），那么如果 \(g,f*g\) 的前缀和不好求怎么办？不难发现，在整除分块的过程中，我们需要对 \(g,h\) 求前缀和的部分可以写成 \(\lfloor\large\frac n k\normalsize\rfloor\) 的形式，如果 \(g,h\) 可以杜教筛，那 \(f\) 也可以杜教筛。

如果我们能快速求出 \(f,g\) 的前缀和，要求 \(f*g\) 的前缀和呢？还是上面那个式子：

\[\sum_{i=1}^n (f*g)(i)=\sum_{i=1}^n g(i)S(\lfloor\frac n i\rfloor) \]

整除分块右边推左边就好了。

常见杜教筛筛法：

• \(f=\mu\)，根据 \(\mu*1=\varepsilon\)，取 \(g=1\)。

• \(f=\varphi\)，根据 \(\varphi*1=\text{Id}\)，取 \(g=\text{Id}\)。

• \(h=\sigma_k\)，根据 \(1*\text{Id}_k=\sigma_k\)，取 \(f=1,g=\text{Id}_k\)。

• \(f=\mu\cdot \text{Id}_k\)，按照提公因式卷上一个 \(1\cdot \text{Id}_k\)，那么 \((\mu\cdot \text{Id}_k)*(1\cdot \text{Id}_k)=(\mu*1)\cdot \text{Id}_k=\varepsilon\cdot \text{Id}_k=\varepsilon\)。取 \(g=\text{Id}_k\)。

  提公因式：对于任意数论函数 \(f,g\) 和完全积性函数 \(h\) 均有 \((f*g)\cdot h=(f\cdot h)*(g\cdot h)\)（注意中间是点积）。

• \(f=\varphi\cdot \text{Id}_k\)，卷 \(1\cdot \text{Id}_k\)，那么 \((\varphi\cdot \text{Id}_k)* (1\cdot \text{Id}_k)=(\varphi*1)\cdot \text{Id}_k=\text{Id}\cdot \text{Id}_k=\text{Id}_{k+1}\)。取 \(g=\text{Id}_k\)。

• \(f=\mu^2*(\mu\cdot \text{Id})\)，卷 \(1\cdot\text{Id}\)，那么 \(\mu^2*(\mu\cdot \text{Id})*(1\cdot \text{Id})=\mu^2*((\mu*1)\cdot \text{Id})=\mu^2*\varepsilon=\mu^2\)，\(\mu^2\) 的前缀和见例题 P4318，\(\sum_{i=1}^n \mu^2(i)=\sum_{i=1}^{\sqrt n}\lfloor\large\frac{n}{i^2}\normalsize\rfloor\mu(i)\)。

二、贝尔级数

对于积性函数 \(f(n)\)，定义其在质数 \(p\) 意义下的贝尔级数为：

\[F_p(x)=\sum_{i=0}^{\infty}f(p^i)x^i \]

可以发现 \(f*g=h\Leftrightarrow \forall p,F_p(x)\times G_p(x)=H_p(x)\)。

如果 \(f\cdot h=h\) 就是对应项系数相乘。

栗子：

• \(\varepsilon_p(x)=1\)。

• \(1_p(x)=\sum_{i\geq 0}1\cdot x^i=\frac{1}{1-x}\)。

• \((\text{Id}_k)_p(x)=\sum_{i\geq 0}(p^i)^k\cdot x^i=\sum_{i\geq 0}(p^kx)^i=\frac{1}{1-p^kx}\)。

• \(\mu_p(x)=1\cdot x^0+(-1)\cdot x^1=1-x\)。

  另一种理解是，\(\mu*1=\varepsilon\Leftrightarrow \mu_p(x)=\frac{\varepsilon_p(x)}{1_p(x)}=\frac{1}{\frac{1}{1-x}}=1-x\)。

• \((\mu^2)_p(x)=1\cdot x^0+1\cdot x^1=1+x\)。

• \(d_p(x)=\sum_{i\geq 0}(i+1)x^i=\frac{1}{(1-x)^2}\)。

  另一种理解是，\(d=1*1\Leftrightarrow d_p(x)=1_p(x)\times 1_p(x)\)。

• \(f(x)=d(x^2)\)。那么 \(f(p^k)=d(p^{2k})=2k+1\)。

  \(f_p(x)=\sum_{i\geq 0}(2i+1)x^i=\frac{1+x}{(1-x)^2}\)。说明 \(f=\mu^2*d\)。

• \((\sigma_k)_p(x)=\sum_{i\geq 0}(\sum_{j=0}^i (p^j)^k)x^i=\frac{1}{(1-x)(1-p^kx)}\)。

  另一种理解是，\(\sigma_k=\text{Id}_k*1\Leftrightarrow (\sigma_k)_p=(\text{Id}_k)_p\times 1_p(x)=\frac{1}{1-p^kx}\times \frac{1}{1-x}=\frac{1}{(1-x)(1-p^kx)}\)。

• \(\varphi_p(x)=1+\sum_{i\geq 1}(p^i-p^{i-1})x^i=1+(1-\frac{1}{p})\sum_{i\geq 1}p^ix^i=1+\frac{p-1}{p}(\frac{1}{1-px}-1)=\frac{1-x}{1-px}\)。

  另一种理解是，\(\varphi=\text{Id}*\mu\Leftrightarrow \varphi_p(x)=\frac{1}{1-px}\times (1-x)=\frac{1-x}{1-px}\)。

• \(w(x)=\sum_{d\mid x}\mu^2(d)=2^{x\ 的不同质因子数}\)。显然 \(w(1)=1\)，\(w(p^k)=2\)（\(d=1,p\)）。

  \(w_p(x)=1+\sum_{i\geq 1}2x^i=1+2(\frac{1}{1-x}-1)=\frac{1+x}{1-x}\)。

  另一种理解是，\(w=\mu^2*1\Leftrightarrow w_p(x)=(1+x)\times \frac{1}{1-x}=\frac{1+x}{1-x}\)

点积 \(\text{Id}_k\) 在贝尔级数中有很好的性质：\((f\cdot \text{Id}_k)_p(x)=\sum_{i\geq 0}f(p^i)(p^i)^kx^i=\sum_{i\geq 0}f(p^i)(p^kx)^i=f_p(p^kx)\)。即将 \(x\) 整体替换为 \(p^kx\)。

• \((\mu\cdot \text{Id}_k)_p(x)=1-p^kx\)。
• \((\varphi\cdot \text{Id}_k)_p(x)=\frac{1-p^kx}{1-p^{k+1}x}\)。

应用到杜教筛找 \(g\) 中。一般都形如 \(\frac{1}{1-p^kx}\)（\(1\cdot \text{Id}_k\)）、\(1-p^kx\)（\(\mu\cdot\text{Id}_k\)）、\(1+p^kx\)（\(\mu^2\cdot \text{Id}_k\)）若干个相乘（\(k\) 可以为 \(0\)）。比如：

• 筛 \(f=\mu\cdot\text{Id}_k\)。

  \(f_p(x)=1-p^kx\)，考虑取 \(g_p(x)=\frac{1}{1-p^kx}\)，此时 \((f*g)_p(x)=1\)。

  \(f=\text{Id}_k\)，\(f*g=\varepsilon\)。

• 筛 \(f=\varphi\cdot \text{Id}_k\)。

  \(f_p(x)=\frac{1-p^kx}{1-p^{k+1}x}\)，考虑取 \(g_p(x)=\frac{1}{1-p^kx}\)，此时 \((f*g)_p(x)=\frac{1}{1-p^{k+1}x}\)。

  \(g=\text{Id}_k\)，\(f*g=\text{Id}_{k+1}\)。

• 筛 \(f=\mu^2*(\mu\cdot \text{Id})\)。

  \(f_p(x)=(1+x)(1-px)\)，考虑取 \(g_p(x)=\frac{1}{1-px}\)，此时 \((f*g)_p(x)=1+x\)。

  \(g=\text{Id}\)，\(f*g=\mu^2\)。

• 筛 \(f(p^k)=\begin{cases}1&(k=0)\\p^k+(-1)^k&(k\geq 1 )\end{cases}\)。

  \(f_p(x)=1+\sum_{i\geq 1}(p^i+(-1)^i)x^i=1+(\frac{1}{1-px}+\frac{1}{1+x}-2)=\frac{1}{1-px}+\frac{1}{1+x}-1\)。

  考虑取 \(g_p(x)=(1-px)(1+x)\)，此时 \(f*g=1+px^2\)。

  \(g=(\mu\cdot \text{Id})*\mu^2\)，是可以套个杜教筛来筛的。

  \(f*g=1+px^2\) 说明 \((f*g)(n)=\prod_{p_i^{e_i}}[e_i=2]p_i=\mu^2(\sqrt n)\cdot \sqrt n\)。所以 \(S_{f*g}(n)=\sum_{i=1}^{\sqrt n}\mu^2(i)\cdot i\)。

  \(S_g,S_{f*g}\) 都能求了，所以 \(S_f\) 也能求。

三、例题

P3768 简单的数学题

给出 \(n,p\)，求：

\[\sum_{i=1}^n\sum_{j=1}^n i\cdot j\cdot \gcd(i,j)\pmod p \]

\(n\leq 10^{10},5\times 10^8\leq p\leq 1.1\times 10^9\)，\(p\) 是质数。

根据 \(\text{Id}=\varphi*1\)：

\[ans=\sum_{i=1}^n\sum_{j=1}^n i\cdot j\cdot \sum_{d\mid \gcd(i,j)}\varphi(d)=\sum_{d=1}^n \varphi(d)d^2(\sum_{i=1}^{\lfloor\frac n d\rfloor}i)^2 \]

设 \(s(n)=(\sum_{i=1}^n i)^2=(\frac{(1+n)n}{2})^2\)，则 \(ans=\sum_{d=1}^n\varphi(d)d^2s(\lfloor\frac{n}{d}\rfloor)\)，如果能快速求出 \(f(x)=\varphi(x)x^2\) 的前缀和，那么整除分块就做完了。

用杜教筛来筛 \(f(x)=\varphi(x)\cdot \text{Id}_2\)：取 \(g=1\cdot \text{Id}_2\)，那么 \(f*g=\text{Id}_3\)。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=2e6+5;
int n,m=2e6,mod,cnt,p[N],phi[N],inv2,inv6,ans;
bool vis[N];
map<int,int>s;
int mul(int x,int n,int mod){
	int ans=mod!=1;
	for(x%=mod;n;n>>=1,x=x*x%mod)
		if(n&1) ans=ans*x%mod;
	return ans;
}
int s2(int n){	//平方和 
	n%=mod;
	return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
}
int s3(int n){	//立方和 
	n%=mod; 
	int x=n%mod*(n+1)%mod*inv2%mod;
	return x*x%mod;
}
int S(int n){
	if(n<=m) return phi[n];
	if(s[n]) return s[n];
	int ans=0;
	for(int l=2,r=0;l<=n;l=r+1)
		r=n/(n/l),ans=(ans+(s2(r)-s2(l-1))%mod*S(n/l)%mod)%mod;
	return s[n]=((s3(n)-ans)%mod+mod)%mod;
}
signed main(){
	scanf("%lld%lld",&mod,&n);
	vis[0]=vis[1]=1,phi[1]=1;
	for(int i=2;i<=m;i++){
		if(!vis[i]) p[++cnt]=i,phi[i]=i-1;
		for(int j=1;j<=cnt&&i*p[j]<=m;j++){
			vis[i*p[j]]=1;
			if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j];break;}
			phi[i*p[j]]=phi[i]*phi[p[j]]; 
		}
	}
	for(int i=1;i<=m;i++) phi[i]=(phi[i-1]+i*i%mod*phi[i]%mod)%mod;
	inv2=mul(2,mod-2,mod),inv6=mul(6,mod-2,mod);
	for(int l=1,r=0;l<=n;l=r+1)
		r=n/(n/l),ans=(ans+s3(n/l)*(S(r)-S(l-1))%mod)%mod;
	printf("%lld\n",(ans+mod)%mod);
	return 0;
} 

P1587 [NOI2016] 循环之美

给出 \(n,m,k\)，求对于 \(1\leq x\leq n,1\leq y\leq m\)，\(\large\frac x y\) 满足数值上互不相等且在 \(k\) 进制下为纯循环小数的个数。

\(1\leq n\leq 10^9\)，\(1\leq m\leq 10^9\)，\(2\leq k\leq 2\times 10^3\)。

取最简分数，即 \(\gcd(i,j)=1\)。

要求是纯循环小数，设 \(\large\frac x y\) 的循环节长度为 \(l\)，那么 \(\large \frac x y\) 的小数部分和 \(\large\frac x y\normalsize\cdot k^l\) 的小数部分相同（也就是说小数点往后挪 \(l\) 位小数部分还是相同，比如 \(1.123123\cdots\) 小数点往后挪 \(3\) 位变成 \(1123.123\cdots\)）。那么 \(\large\frac{x}{y}\normalsize-\lfloor\large\frac{x}{y}\normalsize\rfloor=\large\frac{xk^l}{y}\normalsize-\lfloor\large\frac{xk^l}{y}\normalsize\rfloor\)，\(x-\lfloor\large\frac x y\normalsize\rfloor y=xk^l-\lfloor\large\frac{xk^l}{y}\normalsize\rfloor y\)，\(x\equiv xk^l \pmod y\)。

由于 \(\gcd(x,y)=1\)，所以 \(k^l\equiv 1\pmod y\)，\(\gcd(k,y)=1\)。因此：

\[\begin{aligned} ans=&\sum_{x=1}^n\sum_{y=1}^m[\gcd(x,y)=1][\gcd(y,k)=1]\\ =&\sum_{d\mid k}\mu(d)\sum_{y=1}^{m/d}\sum_{x=1}^n[\gcd(x,yd)=1]\\ =&\sum_{d\mid k}\mu(d)\sum_{y=1}^{m/d}\sum_{x=1}^n[\gcd(x,y)=1][\gcd(x,d)=1] \end{aligned} \]

即 \(f(n,m,k)=\sum_{x=1}^n\sum_{y=1}^m[\gcd(x,y)=1][\gcd(y,k)=1]\)。那么

\[ans=\sum_{d\mid k}\mu(d)f(m/d,n,d) \]

递归求解即可，\(n=0\) 或 \(m=0\) 时 \(f(n,m,k)=0\)。注意特判 \(k=1\) 的情况，此时 \(f(n,m,k)=\sum_{x=1}^n\sum_{y=1}^m[\gcd(x,y)=1]=\sum_{d=1}^{\min(n,m)}\mu(d)\lfloor\frac n d\rfloor\lfloor\frac n d\rfloor\)，整除分块 + 杜教筛即可。

由于每次杜教筛求前缀和的位置都是 \(n,m\) 的关键点，总复杂度是 \(\mathcal O(n^{\frac 2 3}+\text{递归次数})\) 的。

#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=1e6+5;
int n,m,k,w=1e6,cnt,p[N],mu[N],smu[N];
bool vis[N];
map<int,int>s;
int S(int n){
	if(n<=w) return smu[n];
	if(s[n]) return s[n];
	int ans=0;
	for(int l=2,r;l<=n;l=r+1)
		r=n/(n/l),ans+=(r-l+1)*S(n/l); 
	return s[n]=1-ans;
}
LL f(int n,int m,int k){
	if(!n||!m) return 0;
	LL ans=0;
	if(k==1){
		for(int l=1,r;l<=min(n,m);l=r+1)
			r=min(n/(n/l),m/(m/l)),ans+=1ll*(S(r)-S(l-1))*(n/l)*(m/l);
		return ans;
	}
	for(int i=1;i*i<=k;i++) if(k%i==0){
		if(mu[i]) ans+=mu[i]*f(m/i,n,i);
		if(i*i!=k&&mu[k/i]) ans+=mu[k/i]*f(m/(k/i),n,k/i);
	}
	return ans;
}
signed main(){
	scanf("%d%d%d",&n,&m,&k);
	vis[0]=vis[1]=1,mu[1]=1;
	for(int i=2;i<=w;i++){
		if(!vis[i]) p[++cnt]=i,mu[i]=-1; 
		for(int j=1;j<=cnt&&i*p[j]<=w;j++){
			vis[i*p[j]]=1;
			if(i%p[j]==0){mu[i*p[j]]=0;break;}
			mu[i*p[j]]=-mu[i]; 
		}
	}
	for(int i=1;i<=w;i++) smu[i]=smu[i-1]+mu[i];
	printf("%lld\n",f(n,m,k));
	return 0;
}