积性函数前缀和

最近突然做到一些求积性函数前缀和的题，用到了各种筛，有一题用到 min_25 筛法，于是好好学习了一波，运用极不熟练。后来又遇到一道杜教筛的题，结果发现自己连 phi(x) 前缀和都不会推了？吓得我赶紧复习 + 写博客。

 

前置技能 

常见（完全）积性函数；整除分块； dirichlet 卷积；埃氏筛

 

---

 

这里还是简单介绍一下 dirichlet 卷积：

我们用 * 表示两个函数的 dirichlet 卷积。

定义函数 f 和 g 的 dirichlet 卷积形式为 $(f*g)(n)=\sum_{d|n} f(d)g\left (\frac{n}{d}\right )$ 。

 

接下来，设 $f$ 为积性函数，我们要求 $\sum_{i=1}^{n} f(i)$ ， n 是 $10^9\sim 10^{12}$ 左右的级别。本文忽略复杂度分析（需用到积分相关知识，本人太弱不会）。

 

---

 

杜教筛

（先扔链接：https://www.cnblogs.com/zzqsblog/p/5461392.html）

适用范围：存在另两个积性函数 $g,h$ ，且满足 $f*g=h$ ， $g,h$ 函数的前缀和比较容易求。

先来举几个例子：

• 比如说要求 $\mu$ 的前缀和，我们知道 $\mu * 1=e$ （ $e(i)=[i=1]$ ）， $e$ 的前缀和显然恒为 $[n\ge 1]$ ， $1$ 的前缀和也很好求，所以可以使用杜教筛；
• 又比如说要求 $\varphi$ 的前缀和，我们知道 $\varphi * 1=id$ ， $1,id$ 的前缀和都很好求 ，所以也可以使用杜教筛。

那么具体怎么求呢？

 

设 $S(n)=\sum_{i=1}^{n} f(i)$ ，那么有

$$\begin{aligned} \sum_{i=1}^{n} (f*g)(i) &= \sum_{i=1}^{n} \sum_{d|i} f(d)g\left (\frac{i}{d} \right ) \\ S(n)=\sum_{i=1}^n f(i) &= \sum_{i=1}^{n} (f*g)(i) - \sum_{i=1}^{n} \sum_{d|i,d<i} f(d)g\left (\frac{i}{d}\right ) \\ &=\sum_{i=1}^n (f*g)(i) - \sum_{i=2}^n g(i) \sum_{d=1}^{\left \lfloor \frac{n}{i} \right\rfloor} f(d)\\ &=\sum_{i=1}^n (f*g)(i) - \sum_{i=2}^n g(i)S\left (\left \lfloor \frac{n}{i} \right \rfloor \right ) \end{aligned}$$

可以发现我们得到了一个 $S(n)$ 的子问题，可以递归下去求解。

至于边界情况，考虑如果 $n\leq B$ 的时候我们直接用线性筛线性求解 $S(n)$ ，当 B 取 $n^{\frac 23}$ 的时候复杂度最优，大概为 $\mathcal{O}(n^{\frac23})$ 。足以解决这个问题。

 

贴代码 （一百年前写的）

#include<bits/stdc++.h>
#define rep(i,x,y) for (int i=(x);i<=(y);i++)
#define ll long long
using namespace std;
const int N=1e6+5;
int p[N/10],cnt,vis[N],mu[N]; map<ll,ll> mp;
void init(int n){
    mu[1]=1;
    rep (i,2,n){
        if (!vis[i]) p[++cnt]=i,mu[i]=-1;
        for (int j=1;j<=cnt&&i*p[j]<=n;j++){
            vis[i*p[j]]=1;
            if (i%p[j]==0) break;
            mu[i*p[j]]=-mu[i];
        }
    }
    rep (i,2,n) mu[i]+=mu[i-1];
}
ll S(ll n){
    if (n<N) return mu[n];
    if (mp.count(n)) return mp[n];
    ll ret=1;
    for (ll i=2,last;i<=n;i=last+1){
        last=n/(n/i);
        ret-=(last-i+1)*S(n/i);
    }
    return mp[n]=ret;
}
int main(){
    init(N-1);
    ll a,b; cin>>a>>b; cout<<S(b)-S(a-1);
    return 0;
}

 

---

 

min_25 筛

（先扔链接：https://www.cnblogs.com/cjyyb/p/9185093.html）

一些吐槽：这个杜教筛真的是筛？？好吧其实我也不知道它为什么叫杜教筛，这确实只是个递归求和的方法。（ yhx : 应该叫杜氏求和法）

不管这些，下面介绍一种真正的筛—— min_25 筛。它本质和洲阁筛相同，但小一倍常数。众所周知数论题非常容易被卡常数，所以学好 min_25 筛还是很重要的 ~

 

适用范围：函数 f 需要满足 $f(p)$ 是一个低阶多项式， $f(p^c)$ 有较为容易的求法，最好是 $\mathcal{O}(1)$ 。这里 p 代表任意素数。

首先 min_25 筛和洲阁筛都是基于埃氏筛的。

要求 $\sum_{i=1}^{n} f(i)$ ，我们把 $1\sim n$ 中的素数和合数分开考虑。

假设素数集合为 $P$ ， $P_i$ 为第 i 个素数， $minp(i)$ 表示 i 的最小素因子。 

 

素数

以下假设 $f(p)$ 为完全积性函数（但 $f(i)$ 不一定。原因下面会说）。

要求 $\sum_{i=1}^n [i \in P]f(i)$ 。

为了方便计算，这里所有数的函数值都定义为 $f(p)$ ，即都和 $f(p)$ 使用一样的表达式。

设 $g(n,j)=\sum_{i=1}^n [i\in P\ or \ minp(i)>P_j]f(i)$ ，即筛到第 j 个素数为止剩下的数和所有素数的函数值之和。

考虑 $g(n,j)$ 的转移：

•  $P_{j}^2>n$ ：显然 $g(n,j)=g(n,j-1)$ 。
•  $P_{j}^2\leq n$ ：此时需要减去一些贡献。这个贡献怎么求呢？考虑筛 $P_j$ 的过程：找到剩下的数中所有 $P_j$ 的倍数删去，减去它们的函数值。由于 $f$ 为完全积性函数（上面说的原因在这里，如果不是完全积性就无法提出 $f(P_j)$ ，因而无法进行下面的变换），可以提出 $f(P_j)$ ，即 $f(P_j) g\left (\left \lfloor \frac{n}{P_j} \right \rfloor,j-1 \right )$ 。注意到若直接减去这个值，我们会把质数部分减去，而事实上我们需要把质数保留，因此这个贡献实际上是 $$f(P_j) [g\left (\left \lfloor \frac{n}{P_j} \right \rfloor,j-1\right ) - \sum_{i=1}^{j-1} f(P_i)]\\ =f(P_j) [g\left (\left \lfloor \frac{n}{P_j} \right \rfloor,j-1\right ) - g(P_j-1,j-1)]$$ 
• 故 $$g(n,j)=\begin{cases} g(n,j-1) & P_j^2>n \\ g(n,j-1)-f(P_j) [g\left (\lfloor \frac{n}{P_j} \rfloor,j-1 \right) - g(P_j-1,j-1)] & P_j^2\leq n \end{cases}$$

那么 $\sum_{i=1}^n [i \in P]f(i)=g(n,|P|)$ ，这里 $P$ 是前 $\sqrt{n}$ 个素数的集合。

 

【稍微解释一下这个过程】

由于我们最终用到的函数值只有 $g(n,|P|)$ ，即我们这步求的是素数的函数值，按照我们的筛法，实际上是先将所有数都按照 $f(p)$ 的式子计算函数值，然后把合数筛掉，最后只剩下素数，那么就能得到正确的函数值之和。但是中间筛的过程中，如果 $f(x)$ 的表达式和 $f(p)$ 不一样的话， $g$ 的值并不是真实的函数值。所以注意这里只求出了素数函数值之和。

 

合数

这里没有素数那栏里的若干限制。

设 $S(n,j)=\sum_{i=1}^n [minp(i)\geq P_j]f(i)$ 。和 $g$ 类似的定义，但注意这里取到了等号，且不包含素数。

$S$ 的值分 2 部分计算：质数，合数。 1 的函数值放到最后计算。

质数的函数值即 $g(n,|P|)-\sum_{i=1}^{j-1} f(P_i)$ 。至于合数，由于 $f(x)$ 是积性函数，那么可以把最小质因子全部提出来（它和剩下部分互质），所以枚举它的最小质因子 $P_k$ 及幂次 $e$ ，得 $$S(n,j)=g(n,|P|)-\sum_{i=1}^{j-1}f(P_i)+\sum_{k\geq j}\sum_{e>0}^{P_k^{e+1}\leq n}(f(P_k^e)S\left (\left \lfloor \frac{n}{P_k^e} \right \rfloor,k+1 \right )+f(P_k^{e+1}))$$

后面的 $f(P_k^{e+1})$ 是单独处理的特殊情况，由于前面无法统计到。

 

这样总复杂度为 $\mathcal{O}(\frac{n^{\frac34}}{\log n})$ 。

 

【为什么要求 $f(p)$ 为低阶多项式呢】

我们在处理素数部分的时候说过假设 $f(p)$ 为完全积性函数。那么低阶多项式一定能拆成 $\sum_{i} a_ix^i$ ，这里每一项都是完全积性的，所以只需拆项使用上述方法求解即可。复杂度 $\times 阶数$ 。

 

栗子

loj6053 简单的函数

这题的函数比较奇怪，但好在它是积性的； $f(p)$ 除 $p=2$ 时均为 $p-1$ ，是一个低阶多项式； $f(p^c)$ 可以 $\mathcal{O}(1)$ 求解。所以可以用 min_25 筛求解。

把 $f(p)$ 拆掉， $p=2$ 的时候特判一下。别的直接套用上面的公式即可。

code:

#include<bits/stdc++.h>
#define rep(i,x,y) for (int i=(x);i<=(y);i++)
#define ll long long

using namespace std;

const int N=2e5+10,mod=1e9+7;
ll n,w[N]; int m,sq,id1[N],id2[N],h[N],g[N],p[N],cnt,pre[N],vis[N];

inline int get_id(ll x){return x<=sq?id1[x]:id2[n/x];}

void sieve(int n){
    rep (i,2,n){
        if (!vis[i]) p[++cnt]=i,pre[cnt]=(pre[cnt-1]+i)%mod;
        for (int j=1;j<=cnt&&i*p[j]<=n;j++){
            vis[i*p[j]]=1;
            if (i%p[j]==0) break;
        }
    }
}

int S(ll n,int i){
    if (n<=1||p[i]>n) return 0;
    int k=get_id(n);
    int res=((ll)g[k]-h[k]-pre[i-1]+i-1)%mod; res=(res+mod)%mod;
    if (i==1) res+=2;
    for (int j=i;j<=cnt&&(ll)p[j]*p[j]<=n;j++){
        ll pw=p[j];
        for (int e=1;pw<=n/p[j];e++)
            res=(res+((ll)(p[j]^e)*S(n/pw,j+1)%mod+(p[j]^(e+1)))%mod)%mod,pw*=p[j];
    }
    return res;
}

int main(){
    scanf("%lld",&n); sq=sqrtl(n); sieve(sq);
    for (ll i=1,j;i<=n;i=j+1){
        ll x=n/i; j=n/x; w[++m]=x;
        if (x<=sq) id1[x]=m; else id2[j]=m;
        h[m]=(x-1)%mod,g[m]=(x+2)%mod*((x-1)%mod)%mod;
        (g[m]&1)?g[m]+=mod:0; g[m]>>=1;
    }
    rep (j,1,cnt)
        for (int i=1;i<=m&&(ll)p[j]*p[j]<=w[i];i++){
            int k=get_id(w[i]/p[j]);
            h[i]=((ll)h[i]-h[k]+j-1)%mod,h[i]=(h[i]+mod)%mod;
            g[i]=((ll)g[i]+mod-(ll)p[j]*(g[k]+mod-pre[j-1])%mod)%mod;
        }
    printf("%d\n",S(n,1)+1);
    return 0;
}