洲阁筛解决的问题主要是\(n\)范围较大的积性函数前缀和。 已知一积性函数\(f(i)\)，求\(\sum_{i=1}^nf(i)\)。 \(n\leq10^{12}\).
如果\(f(i)\)在质数处的取值比较简单，那么可以运用洲阁筛来求解。 我们需要两个辅助数组。
定义如下：\[ \begin{aligned} g_{i,j}&=\sum_{k=2}^i[k与p_1,p_2,...,p_j互质或就是其中某个质数]\; s(k)\\ &=\sum_{k=2}^i[k是\leq p_j的质数或k的最小质因子大于p_j]\; s(k) \end{aligned} \] 其中\(s(x)\)是一个完全积性函数，它可以是\(s(x)=x\)，或\(s(x)=1\)，等等。 我们一般要求第二维计算到\(m\)，其中\(m\)为满足\(p_m\le\sqrt n\)的最大正整数。 这时候，\(g_{x,m}\)就表示函数\(s\)在前缀范围\([1,x]\)内的质数处的取值之和，这也是\(g\)数组的主要作用。但这里\(x\)的定义域不是\([1,n]\)，下文会提到。 这个数组怎么求呢？ 首先边界条件比较简单：\[ g_{i,0}=\sum_{k=2}^is(k) \] 可是这里有一个限制：\(\sum_{k=2}^is(k)\)必须足够简单好算。如果函数前缀和比较复杂，如\(s=\mu\)，就不能使用在\(s=\mu\)意义下计算\(g\)的这种方法。但是如果函数在质数的取值比较简单，如\(\mu(p)=-1\)，我们可以考虑换一个角度：令\(s(x)=1\)，可以计算出\([1,n]\)的质数个数\(p_{sum}\)，那么我们可以用\(-1*p_{sum}\)表示所求的东西，一样达到了我们的目的（黑体字）。 考虑由\(g_{i,j-1}\)推出\(g_{i,j}\)。 若\(i<p_j^2\)，显然\(g_{i,j}=g_{i,j-1}\)。 否则当\(i\ge p_j^2\)时，如何考虑？\(g_{i,j-1}\)代表着一些在\(j-1\)时合法的数的函数值之和，而从\(j-1\)变成\(j\)时，若某些原本合法的数变得不合法，显然一定是因为触犯了第二个条件：最小质因子恰好为\(p_j\)。如何算出这一部分的数的函数值之和呢？\[ g_{i,j}=g_{i,j-1}-s(p_j)(g_{\lfloor\frac i {p_j}\rfloor,j-1}-g_{p_j-1,j-1}) \] 后面减去的部分就是这一部分数的函数值之和。首先它们都有一个最小质因子\(p_j\)，随后，它们除去\(p_j\)剩下的部分，不可以含有小于\(p_j\)的质因子，因此剩下的部分至少大于等于\(p_j\)，但要小于等于\(\frac{i}{p_j}\)。这恰好对应了\(g\)的定义！后面一部分算的就是\[ \sum_{k=p_j}^{\frac i {p_j}}[最小质因子>p_{j-1}]\; s(k) \] 但是\(n\)太大，存不下怎么办？ 记集合\(S=\{\lfloor\frac n x\rfloor |x\in[1,n]\}\;\;(\)\(|S|=2\sqrt n)\)，我们只需要关注\(g_{i,...}\;\;(i\in S)\)即可，那为什么只用考虑第一维取这些值的情况呢？可以证明，所有以后需要调用的第一维都属于\(S\)。从\(g\)本身的递推需要来看，我们需要\(g_{\lfloor\frac i {p_j}\rfloor,j-1}\)，它的第一维\(\lfloor \frac i {p_j}\rfloor\in S\)；还需要\(g_{p_j-1,j-1}\)，它的第一维\(p_j-1\le\sqrt n\)，一定属于\(S\)。所以，其他的第一维取值就不需要计算了。从下文的\(h\)计算的调用来看，也都只会调用到\(S\)内的第一维，下文会提到。 再之，我们发现每次递归时只用到第二维为\(j-1\)的数据，这给了我们滚动的思想：我们依次计算\(j=0...m\)的情况。将\(|S|\)个元素排成一排，分别表示在当前\(j\)意义下，\(g_{i,j},i\in S\)的取值。由于更新\(g_{i,j}\)的时候只需要用到诸如\(k\le i\)的\(g_{k,j-1}\)，我们从大到小枚举并更新那些需要更新的元素\(i\in S,i\ge p_j^2\)，这样可以保证调用的元素仍然是\(j-1\)意义下的；而\(i<p_j^2\)的情况，与\(j-1\)时完全没有改变，不需要操作，直接继承。 \(j\)越是大，我们需要更新的元素就越来越少。如此一层层地覆盖上去（很像一维背包的那种继承思想），我们，就可以计算到\(g_{n,m}\)了。这一步的复杂度是\(O(\frac {n^{\frac 3 4}}{\log n})\)。 离散\(S\)中元素到数组上的方法很简单，对于\(x\in S\)，如果\(x\le \sqrt n\)，用\(pos_1[x]\)表示\(x\)的离散位置；如果\(x>\sqrt n\)，用\(pos_2[\lfloor \frac n x\rfloor]\)表示\(x\)的离散位置。可以写一个简短的\(get()\)函数来处理询问离散位置的操作。
定义如下：\[ h_{i,j}=\sum_{k=2}^i[k的最小质因子\ge p_j]\;f(k) \] 其递推式为：\[ \begin{aligned} h_{i,j}&=\sum_{p_k\ge p_j}\sum_{e\ge 1且p_k^e\le i}f(p_k^e)h_{\lfloor \frac i{p_k^e}\rfloor,j+1}+f(p_k^e)\\ &=\sum_{p_k\ge p_j}\sum_{e\ge 1且p_k^e\le i}f(p_k^e)(h_{\lfloor \frac i{p_k^e}\rfloor,j+1}+1) \end{aligned} \] 它相当于枚举所有在\([2,i]\)范围内，最小质因子大于等于\(p_j\)的数，并对函数值求和。首先枚举它们的最小质因数\(p_k\)，再枚举最小质因数\(p_k\)的幂；最后统计有多少个数，满足除去\(p_k^e\)后剩余部分最小质因子大于\(p_j\)，也就是大于等于\(p_{j+1}\)，这和\(h\)本身的定义恰好符合。如此枚举，不会算重，但是会漏掉一种情况：\(f(p_k^e)\)没有被算入，因为\(h\)的\(\sum\)是从2开始枚举的。额外加上即可。 我们使用搜索计算\(h\)。这里的搜索不需要记忆化，因为有结论是：如果要求不同的\(h_{i,j}\)，它们在搜索时不会搜索到重复的地方。 第一个循环不可以完全枚举\(k=j...m\)，不然的话复杂度过高。当\(p_k^2>i\)时，余下的未计算的函数值之和恰好是函数在\(p\in(\sqrt i,i]\)处的取值之和，可以直接用\(g_{i,m}-g_{p_{m'},m}\)表示，其中\(m'\)为满足\(p_{m'}\le \sqrt i\)的最大正整数。我们发现\(p_m',i\in S\)，回应上文，只需要计算\(g\)的第一维属于\(S\)时的情况。
\(\sum_{i=2}^nf(i)=h_{n,0}\)。 那么\(\sum_{i=1}^nf(i)=h_{n,0}+f(1)\)，简洁明了。 当然，有时候题目甚至不需要调用\(h\)数组，这要依题目而灵活变化。
整体思路稍微有点复杂，但只要明白了洲阁筛对积性函数求和的关键步骤，就可以比较好地理解。 首先，我们所讨论的积性函数，最好在质数处有简明的表达式。我们可以将表达式写出后，对于每个和式，用\(g\)在不同\(s(x)\)的意义下逐个求解。 然后要理解洲阁筛对分配律的利用。它通过枚举最小质因数、枚举其作为最小质因数的指数、最后统计除去最小质因数后，剩余部分最小质因数大于自己的情况数这一种枚举方法，可以结合积性函数性质、运用\(h\)数组快速枚举遗漏的数，并得到其函数值之和。 题目所求的问题并不会总是中规中矩，我们需要把题目求的表达式拆成一个一个部分，尽量用g或h的定义来表示，依次求解。 总的时间复杂度为\(O(\frac{n^{\frac 3 4}}{\log n})\)。
以SPOJ的DivcntK为例，求的是\(\sum_{i=1}^n \sigma_0(i^k)\)，\(n,k\le 10^{10}\)。这份代码是单组数据的。
#include <cstdio> #include <cmath> #include <algorithm> #include <iostream> using namespace std; typedef long long ll; typedef unsigned long long ull; const int SQRTN=100005; bool vis[SQRTN]; ll p[SQRTN],pcnt; ll n,k,sqrtn; int m; ll a[SQRTN*2],cnt; int pos1[SQRTN],pos2[SQRTN]; ull g[SQRTN*2]; void prework(){ for(int i=2;i<SQRTN;i++){ if(!vis[i]) p[++pcnt]=i; for(int j=1;j<=pcnt&&i*p[j]<SQRTN;j++){ int x=i*p[j]; vis[x]=true; if(i%p[j]==0) break; } } } int gp(ll x){//getpos return x<=sqrtn?pos1[x]:pos2[n/x]; } void Discretization(){ for(ll i=1,j;i<=n;i=j+1){ a[++cnt]=n/i; j=n/(n/i); } reverse(a+1,a+1+cnt); for(int i=1;i<=cnt;i++) if(a[i]<=sqrtn) pos1[a[i]]=i; else pos2[n/a[i]]=i; } void calc_g(){ for(int i=1;i<=cnt;i++) g[i]=a[i]-1; for(int j=1;j<=m;j++) for(int i=cnt;i>=1&&a[i]>=p[j]*p[j];i--) g[i]-=g[gp(a[i]/p[j])]-g[gp(p[j]-1)]; } ull calc_h(ll i,ll j){ if(i<=1) return 0; ull res=0; int a; for(a=j;a<=m&&p[a]*p[a]<=i;a++) for(ll pe=p[a],e=1;pe<=i;pe*=p[a],e++) res+=(ull)(e*k+1)*(calc_h(i/pe,a+1)+1); if(p[a-1]<=i) res+=(ull)(k+1)*(g[gp(i)]-g[gp(p[a-1])]); return res; } int main(){ prework(); scanf("%lld%lld",&n,&k); sqrtn=(ll)sqrt(n); m=upper_bound(p+1,p+1+pcnt,sqrtn)-p-1; Discretization(); calc_g(); ull ans=(ull)(k+1)*(g[gp(n)]-m); for(int i=1;i<=m;i++) for(ll pe=p[i],e=1;pe<=n;pe*=p[i],e++) ans+=(ull)(e*k+1)*(calc_h(n/pe,i+1)+1); ans++; cout<<ans<<endl; return 0; }转载于:https://www.cnblogs.com/RogerDTZ/p/9184689.html
