积性函数前缀和

最近突然做到一些求积性函数前缀和的题,用到了各种筛,有一题用到 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})$ 。足以解决这个问题。

 

贴代码 (一百年前写的)

 1 #include<bits/stdc++.h>
 2 #define rep(i,x,y) for (int i=(x);i<=(y);i++)
 3 #define ll long long
 4 using namespace std;
 5 const int N=1e6+5;
 6 int p[N/10],cnt,vis[N],mu[N]; map<ll,ll> mp;
 7 void init(int n){
 8     mu[1]=1;
 9     rep (i,2,n){
10         if (!vis[i]) p[++cnt]=i,mu[i]=-1;
11         for (int j=1;j<=cnt&&i*p[j]<=n;j++){
12             vis[i*p[j]]=1;
13             if (i%p[j]==0) break;
14             mu[i*p[j]]=-mu[i];
15         }
16     }
17     rep (i,2,n) mu[i]+=mu[i-1];
18 }
19 ll S(ll n){
20     if (n<N) return mu[n];
21     if (mp.count(n)) return mp[n];
22     ll ret=1;
23     for (ll i=2,last;i<=n;i=last+1){
24         last=n/(n/i);
25         ret-=(last-i+1)*S(n/i);
26     }
27     return mp[n]=ret;
28 }
29 int main(){
30     init(N-1);
31     ll a,b; cin>>a>>b; cout<<S(b)-S(a-1);
32     return 0;
33 }
51nod1244 莫比乌斯函数之和

 


 

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:

 1 #include<bits/stdc++.h>
 2 #define rep(i,x,y) for (int i=(x);i<=(y);i++)
 3 #define ll long long
 4 
 5 using namespace std;
 6 
 7 const int N=2e5+10,mod=1e9+7;
 8 ll n,w[N]; int m,sq,id1[N],id2[N],h[N],g[N],p[N],cnt,pre[N],vis[N];
 9 
10 inline int get_id(ll x){return x<=sq?id1[x]:id2[n/x];}
11 
12 void sieve(int n){
13     rep (i,2,n){
14         if (!vis[i]) p[++cnt]=i,pre[cnt]=(pre[cnt-1]+i)%mod;
15         for (int j=1;j<=cnt&&i*p[j]<=n;j++){
16             vis[i*p[j]]=1;
17             if (i%p[j]==0) break;
18         }
19     }
20 }
21 
22 int S(ll n,int i){
23     if (n<=1||p[i]>n) return 0;
24     int k=get_id(n);
25     int res=((ll)g[k]-h[k]-pre[i-1]+i-1)%mod; res=(res+mod)%mod;
26     if (i==1) res+=2;
27     for (int j=i;j<=cnt&&(ll)p[j]*p[j]<=n;j++){
28         ll pw=p[j];
29         for (int e=1;pw<=n/p[j];e++)
30             res=(res+((ll)(p[j]^e)*S(n/pw,j+1)%mod+(p[j]^(e+1)))%mod)%mod,pw*=p[j];
31     }
32     return res;
33 }
34 
35 int main(){
36     scanf("%lld",&n); sq=sqrtl(n); sieve(sq);
37     for (ll i=1,j;i<=n;i=j+1){
38         ll x=n/i; j=n/x; w[++m]=x;
39         if (x<=sq) id1[x]=m; else id2[j]=m;
40         h[m]=(x-1)%mod,g[m]=(x+2)%mod*((x-1)%mod)%mod;
41         (g[m]&1)?g[m]+=mod:0; g[m]>>=1;
42     }
43     rep (j,1,cnt)
44         for (int i=1;i<=m&&(ll)p[j]*p[j]<=w[i];i++){
45             int k=get_id(w[i]/p[j]);
46             h[i]=((ll)h[i]-h[k]+j-1)%mod,h[i]=(h[i]+mod)%mod;
47             g[i]=((ll)g[i]+mod-(ll)p[j]*(g[k]+mod-pre[j-1])%mod)%mod;
48         }
49     printf("%d\n",S(n,1)+1);
50     return 0;
51 }
loj6053 - 简单的函数

 

posted @ 2019-04-22 10:10  bestfy  阅读(895)  评论(0编辑  收藏  举报