SPOJ DIVCNT2（莫比乌斯反演+杜教筛）

题目传送门：http://www.spoj.com/problems/DIVCNT2/

---

题目分析：
千古神题，这题想了我两天也没想出来……
好吧，其实这题并没有用杜教筛，它并没有使用记忆化递归，也没有用狄利克雷卷积，只不过它的时间复杂度证明和杜教筛类似，也是O(n23)（虽然我还是不知道具体是怎么证的）。该怎么说呢？至少加深了对积性函数前缀和的理解吧。
回到正文：

---

题目要我们求：

∑i=1nσ0(i2)

相当于求：

∑i=1n∑d|i21

如果我们改为先枚举d，接下来我们还要枚举一层什么呢？我们很容易发现，假设：

d=pk11pk22……pkxx

那么它能够贡献的最小的i就是：

h(d)=i=p⌊k1+12⌋1p⌊k2+12⌋2……p⌊kx+12⌋x

同时，h(d)的平方以及它的倍数的平方也会被d更新到。
于是，原式变为：

∑d=1n2⌊nh(d)⌋

然后我们再改一下枚举的顺序，我们令j=h(d)，那么当我们枚举j的时候，有多少个d会满足要求呢？我们分析一下，设：

j=pk11pk22……pkxx

则满足h(d)=j的d为：

d=p2k1(−1)1p2k2(−1)2……p2kx(−1)x

（-1）表示可以减一，也可以不减。
于是我们就发现了，满足条件的d有2w(j)个，其中w(j)表示j的不同质因数个数。
于是原式化成：

∑j=1n2w(j)⌊nj⌋

然而这还是不可做，我们又发现，2w(j)本质上就是指j所有约数中的无平方因数个数，而它们有一个共同的特点：μ不为0（即为正负1），于是μ2等于1。原式化为：

∑j=1n⌊nj⌋∑d|jμ2(d)

我们令i=jd，然后将枚举d拉到外层，内层枚举i，得到：

∑d=1nμ2(d)∑i=1⌊nd⌋⌊nid⌋

然后我们记右边的值为G(⌊nd⌋)，并记F(n)=∑nd=1μ2(d)。
我们同时还发现：F(n)本质上就是1~n中无平方因子的数的个数，我们要求单独的一个F(n)，可以通过一个n√的枚举来实现：

F(n)=∑i=1n√μ(i)⌊ni2⌋

至于还不知道为什么的同学可以去看看这题：bzoj2440完全平方数。
而我们知道，单独求一个G(n)也是可以用n√的时间做的。
而且我们可以用线性筛在O(n)的时间内求出1~n的所有F值（求个μ再平方再前缀和），G值（即约数个数的前缀和），并O(1)解决询问。
现在我们看回原式：

∑d=1nμ2(d)G(⌊nd⌋)

我们对G进行下底函数分块，假设⌊nd⌋相等的一段是[L,R]，那就对答案产生贡献：(F(R)-F(L-1))*G(⌊nd⌋)。
之前说过，我们可以通过预处理达到O(1)返回这三个值，也可以用n√的时间计算这些值，于是我们要在两者之间寻找一个平衡。由类似杜教筛时间复杂度证明可知，当我们预处理n23的时候，时间复杂度最优，达到O(n23)。

---

CODE：

#include<iostream>#include<string>#include<cstring>#include<cmath>#include<cstdio>#include<cstdlib>#include<stdio.h>#include<algorithm>using namespace std;const int maxn=1e8;const int maxt=10010;typedef long long LL;bool vis[maxn];int prime[maxn/10];int cur=0;int miu[maxn];LL f[maxn];int pk[maxn];LL g[maxn];LL que[maxt];int t;LL n=0,N;void Preparation(){    miu[1]=1,g[1]=1;    for (int i=2; i<N; i++)    {        if (!vis[i]) prime[++cur]=i,miu[i]=-1,pk[i]=1,g[i]=2;        for (int j=1; j<=cur && i*prime[j]<N; j++)        {            int k=i*prime[j];            vis[k]=true;            if (i%prime[j]) miu[k]=-miu[i],pk[k]=1,g[k]=2*g[i];            else            {                miu[k]=0;                pk[k]=pk[i]+1;                g[k]=g[i]/pk[k]*(pk[k]+1);                break;            }        }    }    for (int i=1; i<N; i++) f[i]=miu[i]*miu[i];    for (int i=2; i<N; i++) f[i]+=f[i-1],g[i]+=g[i-1];}LL G(LL x){    if (x<N) return g[x];    LL temp=0,last;    for (LL i=1; i<=x; i=last+1)    {        last=x/(x/i);        temp+=( x/i*(last-i+1LL) );    }    return temp;}LL F(LL x){    if (x<N) return f[x];    int sx=(int)floor( sqrt( (double)x )+0.0000001 );    LL temp=0;    for (int i=1; i<=sx; i++)        temp+=( (long long)miu[i]*( x/( (long long)i*(long long)i ) ) );    return temp;}int main(){    freopen("c.in","r",stdin);    freopen("c.out","w",stdout);    scanf("%d",&t);    for (int i=1; i<=t; i++)    {        cin>>que[i];        n=max(n,que[i]);    }    if (n<=10000) N=n+1;    else N=maxn;    Preparation();    for (int q=1; q<=t; q++)    {        n=que[q];        LL ans=0,last;        for (LL i=1; i<=n; i=last+1)        {            last=n/(n/i);            ans+=( G(n/i)*( F(last)-F(i-1) ) );        }        cout<<ans<<endl;    }    return 0;}