问题描述
洲阁筛解决的问题主要是(n)范围较大的积性函数前缀和。
已知一积性函数(f(i))，求(sum_{i=1}^nf(i))。
(nleq10^{12}).
求解方法
如果(f(i))在质数处的取值比较简单，那么可以运用洲阁筛来求解。
我们需要两个辅助数组。
(g_{i,j})
定义如下：
其中(s(x))是一个完全积性函数，它可以是(s(x)=x)，或(s(x)=1)，等等。
我们一般要求第二维计算到(m)，其中(m)为满足(p_mlesqrt n)的最大正整数。
这时候，(g_{x,m})就表示函数(s)在前缀范围([1,x])内的质数处的取值之和，这也是(g)数组的主要作用。但这里(x)的定义域不是([1,n])，下文会提到。
这个数组怎么求呢？
首先边界条件比较简单：
可是这里有一个限制：(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})。
否则当(ige p_j^2)时，如何考虑？(g_{i,j-1})代表着一些在(j-1)时合法的数的函数值之和，而从(j-1)变成(j)时，若某些原本合法的数变得不合法，显然一定是因为触犯了第二个条件：最小质因子恰好为(p_j)。如何算出这一部分的数的函数值之和呢？
后面减去的部分就是这一部分数的函数值之和。首先它们都有一个最小质因子(p_j)，随后，它们除去(p_j)剩下的部分，不可以含有小于(p_j)的质因子，因此剩下的部分至少大于等于(p_j)，但要小于等于(frac{i}{p_j})。这恰好对应了(g)的定义！后面一部分算的就是
但是(n)太大，存不下怎么办？
记集合(S={lfloorfrac n x
floor |xin[1,n]};;()(|S|=2sqrt n))，我们只需要关注(g_{i,...};;(iin S))即可，那为什么只用考虑第一维取这些值的情况呢？可以证明，所有以后需要调用的第一维都属于(S)。从(g)本身的递推需要来看，我们需要(g_{lfloorfrac i {p_j}
floor,j-1})，它的第一维(lfloor frac i {p_j}
floorin S)；还需要(g_{p_j-1,j-1})，它的第一维(p_j-1lesqrt n)，一定属于(S)。所以，其他的第一维取值就不需要计算了。从下文的(h)计算的调用来看，也都只会调用到(S)内的第一维，下文会提到。
再之，我们发现每次递归时只用到第二维为(j-1)的数据，这给了我们滚动的思想：我们依次计算(j=0...m)的情况。将(|S|)个元素排成一排，分别表示在当前(j)意义下，(g_{i,j},iin S)的取值。由于更新(g_{i,j})的时候只需要用到诸如(kle i)的(g_{k,j-1})，我们从大到小枚举并更新那些需要更新的元素(iin S,ige p_j^2)，这样可以保证调用的元素仍然是(j-1)意义下的；而(i<p_j^2)的情况，与(j-1)时完全没有改变，不需要操作，直接继承。
(j)越是大，我们需要更新的元素就越来越少。如此一层层地覆盖上去（很像一维背包的那种继承思想），我们，就可以计算到(g_{n,m})了。这一步的复杂度是(O(frac {n^{frac 3 4}}{log n}))。
离散(S)中元素到数组上的方法很简单，对于(xin S)，如果(xle sqrt n)，用(pos_1[x])表示(x)的离散位置；如果(x>sqrt n)，用(pos_2[lfloor frac n x
floor])表示(x)的离散位置。可以写一个简短的(get())函数来处理询问离散位置的操作。
(h_{i,j})
定义如下：
其递推式为：
它相当于枚举所有在([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)时，余下的未计算的函数值之和恰好是函数在(pin(sqrt i,i])处的取值之和，可以直接用(g_{i,m}-g_{p_{m'},m})表示，其中(m')为满足(p_{m'}le sqrt i)的最大正整数。我们发现(p_m',iin 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,kle 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;
}