[SPOJ] DIVCNT2 - Counting Divisors (square) 积性函数前缀和

题面
题目要我们求

$$\sum_{i=1}^n\sigma_0(i^2)$$
把 $i$分解质因数$p_1^{\alpha_1}*p_2^{\alpha_2}...p_k^{\alpha_k}$
就是求 $$\sum_{i=1}^n\prod_{j=1}^k(2\alpha_j+1)$$
考虑每个乘式选 $2\alpha_j$还是 $1$，枚举集合就是 $$\sum_{i=1}^n\sum_{S\subseteq\{1,2,...,k\}}2^{|S|}*\prod_{j}\alpha_j$$
把 $S$集合看成选$\alpha_j>0$的，把 $\prod_{j}\alpha_j$看成枚举每个 $>0$的 $\alpha_j$，事实上就是枚举 $i$的约数，每个约数的贡献是$2^{\omega(d)}$，其中 $\omega(i)$表示 $i$的质因数个数，所以就是 $$\sum_{i=1}^n\sum_{d|n}2^{\omega(d)}$$
其实把 $\omega(d)$看成每个质因数选或不选，就是 $d$的无平方因子约数的个数，所以式子变成
$$\sum_{i=1}^n\sum_{d|n}\sum_{p|d}\mu^2(p)$$
改成枚举 $\mu^2(i)$，就是
$$\sum_{i=1}^n\mu^2(i)\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\sigma_0(j)$$
于是利用整除分块，我们只要知道 $\mu^2(i)$和 $\sigma_0(i)$的前缀和，按照杜教筛的复杂度分析，就能在预处理前 $n^{\frac{2}{3}}$项的情况下， $O(n^{\frac{2}{3}})$求出。
前者就是 $1~n$的无平方因子数个数，考虑容斥 $O(\sqrt n)$
$$\sum_{i=1}^n\mu^2(i)=\sum_{i=1}^{\sqrt n}\mu(i)*\lfloor\frac{n}{i^2}\rfloor$$
后者就是sb套路，整除分块 $O(\sqrt n)$
$$\sum_{i=1}^n\sigma_0(i)=\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor$$
代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#define ll long long
using namespace std;
const int maxn=50000000;
ll n;
int pri[10000000],num,mu[maxn+5],d[maxn+5],sm[maxn+5],sd[maxn+5],ci[maxn+5];
bool flag[maxn+5];
void getpri(int n)
{
    memset(flag,1,sizeof(flag));
    flag[1]=0;mu[1]=1;d[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(flag[i]) pri[++num]=i,mu[i]=-1,ci[i]=1,d[i]=2;
        for(int j=1;j<=num&&i*pri[j]<=n;j++)
        {
            int v=i*pri[j];
            flag[v]=0;ci[v]=1;
            if(i%pri[j]==0) {ci[v]+=ci[i];mu[v]=0;d[v]=d[i]*(ci[v]+1)/ci[v];break;}
            mu[v]=-mu[i];d[v]=(d[i]<<1);
        }
    }
}
ll qsd(ll m)
{
    if(m<=maxn) return sd[m];
    ll re=0;
    for(ll l=1,r=1,k;l<=m;l=r+1)
    {
        k=m/l;r=m/k;
        re+=(ll)(r-l+1)*k;
    }
    return re;
}
ll qsm(ll m)
{
    if(m<=maxn) return sm[m];
    ll gen=round(sqrt(m)+1),re=0;
    for(ll i=1;i<=gen;i++)
        re+=mu[i]*(m/(i*i));    
    return re;  
}
ll solve()
{
    ll re=0;
    for(ll l=1,r=1,k;l<=n;l=r+1)
    {
        k=n/l;r=n/k; 
        re+=qsd(k)*(qsm(r)-qsm(l-1));
    }
    return re;
}
int main()
{
    int ca,lim=maxn;
    scanf("%d",&ca);
    if(ca>800) lim=10000;
        getpri(lim);    
        for(int i=1;i<=lim;i++)
            sm[i]=sm[i-1]+mu[i]*mu[i],sd[i]=sd[i-1]+d[i];
    while(ca--)
    {
        scanf("%lld",&n);
        printf("%lld\n",solve());
    }
    return 0;
}