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

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

---

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

---

题目要我们求：

$$\sum_{i=1}^{n}\sigma_{0}(i^2)$$
相当于求：
$$\sum_{i=1}^{n}\sum_{d|i^2}1$$
如果我们改为先枚举d，接下来我们还要枚举一层什么呢？我们很容易发现，假设：
$$d=p_1^{k_1}p_2^{k_2}……p_x^{k_x}$$
那么它能够贡献的最小的i就是：
$$h(d)=i=p_1^{\lfloor\frac{k_1+1}{2}\rfloor}p_2^{\lfloor\frac{k_2+1}{2}\rfloor}……p_x^{\lfloor\frac{k_x+1}{2}\rfloor}$$
同时，h(d)的平方以及它的倍数的平方也会被d更新到。
于是，原式变为：
$$\sum_{d=1}^{n^2}\lfloor\frac{n}{h(d)}\rfloor$$
然后我们再改一下枚举的顺序，我们令j=h(d)，那么当我们枚举j的时候，有多少个d会满足要求呢？我们分析一下，设：
$$j=p_1^{k_1}p_2^{k_2}……p_x^{k_x}$$
则满足h(d)=j的d为：
$$d=p_1^{2k_1(-1)}p_2^{2k_2(-1)}……p_x^{2k_x(-1)}$$
（-1）表示可以减一，也可以不减。
于是我们就发现了，满足条件的d有 $2^{w(j)}$个，其中w(j)表示j的不同质因数个数。
于是原式化成：
$$\sum_{j=1}^{n}2^{w(j)}\lfloor\frac{n}{j}\rfloor$$
然而这还是不可做，我们又发现， $2^{w(j)}$本质上就是指j所有约数中的无平方因数个数，而它们有一个共同的特点： $\mu$不为0（即为正负1），于是 $\mu^2$等于1。原式化为：
$$\sum_{j=1}^{n}\lfloor\frac{n}{j}\rfloor\sum_{d|j}\mu^2(d)$$
我们令 $i=\frac{j}{d}$，然后将枚举d拉到外层，内层枚举i，得到：
$$\sum_{d=1}^{n}\mu^2(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{id}\rfloor$$
然后我们记右边的值为 $G(\lfloor\frac{n}{d}\rfloor)$，并记 $F(n)=\sum_{d=1}^{n}\mu^2(d)$。
我们同时还发现：F(n)本质上就是1~n中无平方因子的数的个数，我们要求单独的一个F(n)，可以通过一个 $\sqrt{n}$的枚举来实现：
$$F(n)=\sum_{i=1}^{\sqrt{n}}\mu(i)\lfloor\frac{n}{i^2}\rfloor$$
至于还不知道为什么的同学可以去看看这题： bzoj2440完全平方数。
而我们知道，单独求一个G(n)也是可以用 $\sqrt{n}$的时间做的。
而且我们可以用线性筛在O(n)的时间内求出1~n的所有F值（求个 $\mu$再平方再前缀和），G值（即约数个数的前缀和），并O(1)解决询问。
现在我们看回原式：
$$\sum_{d=1}^{n}\mu^2(d)G(\lfloor\frac{n}{d}\rfloor)$$
我们对G进行下底函数分块，假设 $\lfloor\frac{n}{d}\rfloor$相等的一段是[L,R]，那就对答案产生贡献：(F(R)-F(L-1))*G( $\lfloor\frac{n}{d}\rfloor$)。
之前说过，我们可以通过预处理达到O(1)返回这三个值，也可以用 $\sqrt{n}$的时间计算这些值，于是我们要在两者之间寻找一个平衡。由类似杜教筛时间复杂度证明可知，当我们预处理 $n^{\frac{2}{3}}$的时候，时间复杂度最优，达到 $O(n^{\frac{2}{3}})$。

---

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;
}