[日常摸鱼]积性函数求和——杜教筛
因为博主比较菜所以可能一些地方写的有问题或者不清楚，以及我的废话好像有点多…
本文（大概）会不定期更新…一些东西的证明这里可能暂时没有
以及在这里先感谢下小伙伴 @MoebiusMeow 的帮助~ww
参考资料：
[1]浅谈一类积性函数的前缀和（skywalkert）
[2]杜教筛——省选前的学习1（_rqy）
[3]杜教筛（zzq）
[4]具体数学（Concrete Mathematics）
前置技能（一些定义）
约定$[p]$表示满足条件$p$时这个值为1，否则为0
欧拉函数：用$\phi(n)$表示小于$n$且和$n$互质的整数的个数，若$n=p_1^{q_1}*p_2^{q_2}*...*p_k^{q_k}$那么有一种求法是$\phi(n)=n*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(1-\frac{1}{p_k})$
莫比乌斯函数：一种$\mu(n)$的定义是$\mu(1)=1$，对于无平方因子的$n=\prod_{i=1}^k p_i$，$\mu(n)=(-1)^k$，其他情况$\mu(n)=0$，另一种等价的定义是$\sum_{d|n}\mu(d)=[n=1]$
数论函数：若$f:Z^{+} \rightarrow C$，则称$f$为数论函数
积性函数：若一个数论函数$f(n)$对于所有$m_1 \bot m_2$有$f(m_1*m_2)=f(m_1)*f(m_2)$且$f(1)=1$那么称$f(n)$是积性函数，如果对于任意$m_1,m_2$都满足上面那个条件那么就称$f(n)$是完全积性函数，常见的比如莫比乌斯函数、欧拉函数、幂函数（$id_k(n)=n^{k}）$、除数函数（$\sigma_k(n)=\sum_{d|n}d^k$））（[1]）
如果$f$是一个积性函数，把$n$写成$p_1*p_2*...*p_k$的形式，那么可以根据积性就得到$f(n)=\prod_{i=1}^k f(p_i)$
狄利克雷卷积：对两个函数$f,g$记他们的狄利克雷卷积为：$$ (f*g)(n)=\sum_{d|n}f(d)*g(\frac{n}{d}) $$ 很明显$(f*g)=(g*f)$
狄利克雷卷积的单位元函数为$e(n)=[n=1]$，即$f*e=e*f=f$，如果$f,g$都是积性函数那么$f*g$也是积性函数，函数集和狄利克雷卷积能够构成群
莫比乌斯反演：$g(n)=\sum_{d|n}f(d) \Rightarrow f(n)=\sum_{d|n}\mu(d)g(\frac{n}{d})$，反过来也成立，写成卷积的形式就是$g=1*f \Rightarrow f=\mu*g$，证明就是给$g=1*f$两边卷上$\mu$，变成$g*\mu=1*f*\mu=1*\mu*f=e*f=f$
杜教筛
问题：快速计算一个$f$的前缀和：$S(n)=\sum_{i=1}^n f(i)$（比如欧拉函数，$n$可以到1e10的级别）
假设找到了一个合适的函数$g$，他们卷积的前缀和（QwQ我这里写的会比较繁琐很多…）：
$$ \begin{aligned} \sum_{i=1}^n(f*g)(n) & = \sum_{i=1}^n \sum_{d \mid i} g(d)f( \frac{i}{d}) = \sum_{i=1}^n \sum_{d=1}^n [d|i] g(d) f(\frac{i}{d}) \\ & = \sum_{d=1}^n g(d) \sum_{i=1}^n [d|i] f(i) = \sum_{d=1}^n g(d) \sum_{i=1}^n \sum_{t=1}^n [i=dt] f(i) \\ & = \sum_{d=1}^n g(d) \sum_{t=1}^n \sum_{i=1}^n [i=dt] f(i) = \sum_{d=1}^n g(d) \sum_{t=1}^{\lfloor\frac{n}{d} \rfloor } \sum_{i=1}^n [i=dt]f(i) \\ & = \sum_{d=1}^n g(d) \sum_{t=1}^{\lfloor \frac{n}{d} \rfloor } f(i) = \sum_{d=1}^n g(d)S(\lfloor \frac{n}{d} \rfloor)\end{aligned}$$
（最后三行是因为当$t>\lfloor \frac{n}{d} \rfloor$时后面的和式为0，所以这时候可以直接把上界缩小，以及对于每一对$d,t$页只有一个$i$与之对应，这时候后面的和式为1）
然后把$d=1$的提出来就得到了：
$$g(1)S(n)=\sum_{i=1}^n (f*g)(i)-\sum_{i=2}^n g(i)S(\lfloor \frac{n}{i} \rfloor)$$
这里的$\lfloor \frac{n}{i} \rfloor$一共只有$O(\sqrt{n})$种取值，如果能够快速处理掉$g(1)$以及算出卷积的前缀和跟每一段连续的$S(\lfloor \frac{n}{i} \rfloor)$的和那么根据[1]里面说的就可以在$O(n^{\frac{3}{4}})$的时间里算出$S(n)$，如果$f$是一个积性函数，那么利用积性可以筛出前$n^{\frac{2}{3}}$项让复杂度降到$O(n^{\frac{2}{3}})$级别。
实现上有个trick就是可以用两个数组把答案存下来，对$<\sqrt{N}$的$x$放到$p1[x]$里，对$>=\sqrt{N}$的放到$p2[n/x]$里[3]
类似这样
inline lint & get(lint x) { if(x<G)return p1[x]; return p2[n/x]; }
如果有预处理的话这里也可以省去这个p1
举栗子
1.莫比乌斯函数的前缀和（51nod1244）
根据$\mu$的定义卷上一个$1$就得到了$[n=1]$，代到上面的式子直接得到$S(n)=1-\sum_{i=2}^nS(\lfloor \frac{n}{i} \rfloor)$
2.欧拉函数的前缀和（51nod1239）
注意到$\sum_{d|n} \varphi(d)=n$，这就相当于$\varphi * 1=n$，同样直接代进去就可以得到$S(n)=\frac{n(n+1)}{2}-\sum_{i=2}^n S(\lfloor \frac{n}{i} \rfloor)$
（话说bzoj3944那个就是要同时求$\mu$跟$\phi$的前缀和，然后我一开始还想着先全部读进来然后让最大的那个等于$n$然后再去配合前面那个trick算它可以快一点…然后发现这样根本就是错的…根本就不是一个n…）
（bzoj那题的代码）
#include<iostream> #include<cstdio> #include<cstring> using namespace std; typedef long long lint; typedef pair<lint,lint> pa; const lint N=2000005; const lint G=100005; lint T,tot,n; lint pri[N],phi[N],mu[N]; bool p[N]; pa p1[G],p2[G]; inline pa & get(lint x) { if(x<G)return p1[x]; return p2[n/x]; } inline void init() { phi[1]=mu[1]=1;p[1]=1; for(register int i=2;i<N;i++) { if(!p[i]) { pri[++tot]=i;phi[i]=i-1;mu[i]=-1; } for(register int j=1;j<=tot&&i*pri[j]<N;j++) { p[i*pri[j]]=1; if(i%pri[j]==0) { phi[i*pri[j]]=phi[i]*pri[j];mu[i*pri[j]]=0;break; }phi[i*pri[j]]=phi[i]*(pri[j]-1);mu[i*pri[j]]=-mu[i]; } } for(register int i=1;i<N;i++)phi[i]+=phi[i-1],mu[i]+=mu[i-1]; } inline pa calc(lint x) { if(x<N)return (pa){phi[x],mu[x]}; pa & cur=get(x); if((cur.first)!=-1)return cur; lint res1,res2,pos;res1=x*(x+1)/2;res2=1; for(register lint i=2;i<=x;i=pos+1) { pos=x/(x/i);pa temp=calc(x/i); res1-=(pos-i+1)*temp.first; res2-=(pos-i+1)*temp.second; }return cur=(pa){res1,res2}; } int main() { scanf("%lld",&T);init();memset(p1,-1,sizeof(p1)); while(T--) { memset(p2,-1,sizeof(p2)); scanf("%lld",&n); if(n<N)printf("%lld %lld\n",phi[n],mu[n]); else { pa temp=calc(n); printf("%lld %lld\n",temp.first,temp.second); } } return 0; }
3.$f(n)=\varphi(n)*n$的前缀和$S(n)$
这时候可以构造出$g(n)=n$然后卷上$f$就得到：$(f*g)(n)=\sum_{d|n} \varphi(d)*d*\frac{n}{d}=n \sum_{d|n} \varphi(d)=n^2$，然后就得到了$S(n)=\frac{n(n+1)(2n+1)}{6}-\sum_{i=2}^n i*S(\lfloor \frac{n}{i} \rfloor)$，每一段相同的$S(\lfloor \frac{n}{i} \rfloor)$前面的系数也可以直接算出来
4.已知函数$f(n):\sum_{d|n}f(d)=n^2-3n+2$，求$\sum{i=1}^n f(i)$，多组数据，$n<=1e9$（HDU5608）
注意到已知的条件可以写成$1*f=n^2-3n+2$右边是一个很好求和的东西，那根据上面的结论就得到了$f$的前缀和$S(n)=\sum_{i=1}^n (i^2-3i+2)-\sum_{i=2}^nS(\lfloor \frac{n}{i}\rfloor)$，经过一番化简就得到了$S(n)=\frac{n(n-1)(n-2)}{3}-\sum_{i=2}^nS(\lfloor \frac{n}{i} \rfloor)$，对于前面小的数据用莫比乌斯反演预处理，注意数据范围~
（存个代码）
#include<cstdio> #include<cstring> typedef long long lint; const lint MOD=(1e9+7); const lint N=1000005; const lint G=40005; lint T,n,tot,inv3; lint p2[G],mu[N],pri[N],sum[N]; bool p[N]; inline lint pow_mod(lint a,lint b,lint p) { lint res=1; for(;b;b>>=1,a=(a*a)%p)if(b&1)res=(res*a)%p; return res; } inline void init() { mu[1]=1;p[1]=1; for(register int i=2;i<N;i++) { if(!p[i]) { pri[++tot]=i;mu[i]=-1; } for(register int j=1;j<=tot&&i*pri[j]<N;j++) { lint temp=i*pri[j];p[temp]=1; if(i%pri[j]==0){mu[temp]=0;break;} mu[temp]=-mu[i]; } } for(register lint i=1;i<N;i++) for(register lint j=i;j<N;j+=i)sum[j]=(sum[j]+(i*i%MOD-3*i%MOD+2)*mu[j/i])%MOD; for(register int i=1;i<N;i++)sum[i]+=sum[i-1],sum[i]=(sum[i]%MOD+MOD)%MOD; } inline lint calc(lint x) { if(x<N)return sum[x]; lint &cur=p2[n/x]; if(cur!=-1)return cur; lint res,pos; res=x*(x-1)%MOD*(x-2)%MOD*inv3%MOD; for(register lint i=2;i<=x;i=pos+1) { pos=x/(x/i); res-=(pos-i+1)*calc(x/i); while(res<0)res+=MOD; }res=(res%MOD+MOD)%MOD; return cur=res; } int main() { inv3=pow_mod(3,MOD-2,MOD); scanf("%lld",&T);init(); while(T--) { scanf("%lld",&n); if(n<N)printf("%lld\n",sum[n]); else { memset(p2,-1,sizeof(p2)); printf("%lld\n",calc(n)); } } return 0; }
待更...
