【BZOJ3994】约数个数和（SDOI2015）-莫比乌斯反演+数论分块

测试地址：约数个数和
做法：本题需要用到莫比乌斯反演+数论分块。
首先，因为$ij$的每个因数都可以唯一表示成$\frac{pj}{q}$的形式，其中$p|i,q|j,\gcd(p,q)=1$，所以有：
$d(ij)=\sum_{x|i}\sum_{y|j}[\gcd(x,y)=1]$
把这个式子带进原式中去，有：
$ans=\sum_{i=1}^n\sum_{j=1}^m\sum_{x|i}\sum_{y|j}[\gcd(x,y)=1]$
把$x,y$换出来，得到：
$ans=\sum_{x=1}^n\sum_{y=1}^m[\gcd(x,y)=1]\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor$
由莫比乌斯反演定理，有$\sum_{d|n}\mu(d)=[n=1]$（可以从狄利克雷卷积的形式来理解，因为$I=I*e$，所以有$e=\mu*I$），所以：
$ans=\sum_{x=1}^n\lfloor\frac{n}{x}\rfloor\sum_{y=1}^m\lfloor\frac{m}{y}\rfloor\sum_{d|\gcd(x,y)}\mu(d)$
把$d$换出来，得到：
$ans=\sum_{d=1}^{\min(n,m)}\mu(d)\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{xd}\rfloor\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\lfloor\frac{m}{yd}\rfloor$
我们知道$\lfloor\frac{n}{xy}\rfloor=\Big\lfloor\frac{\lfloor\frac{n}{x}\rfloor}{y}\Big\rfloor$，简单证明如下：
令$\lfloor\frac{n}{xy}\rfloor=k$，则有$n=k(xy)+r$，其中$0\le r<xy$，又令$r=r_0x+r_1$，其中$0\le r_1<x$，可以知道$0\le r_0<y$，那么有$\lfloor\frac{n}{x}\rfloor=ky+r_0$，所以$\Big\lfloor\frac{\lfloor\frac{n}{x}\rfloor}{y}\Big\rfloor=k$。
所以上式可以写成：
$ans=\sum_{d=1}^{\min(n,m)}\mu(d)\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\Big\lfloor\frac{\lfloor\frac{n}{d}\rfloor}{x}\Big\rfloor\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\Big\lfloor\frac{\lfloor\frac{m}{d}\rfloor}{y}\Big\rfloor$
那么我们只需预处理出$\mu$的前缀和，以及找到计算后面那一串东西的快速方法，就可以数论分块了。令$S(n)=\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor$，那么上式可以写成：
$ans=\sum_{d=1}^{\min(n,m)}\mu(d)S(\lfloor\frac{n}{d}\rfloor)S(\lfloor\frac{m}{d}\rfloor)$
我们又知道$S(n)=\sum_{i=1}^nd(i)$，为什么呢？因为看$S$的第一种表示，实际上是枚举一个数，这个数对答案做出它的倍数个数的贡献，那么就可以换乘枚举每个数，这个数对答案做出它的因数个数的贡献，所以两种表示是等价的。那么显然我们可以线性筛出$d$来预处理$S$，那么这个问题就解决了。
于是采用数论分块算$ans$，时间复杂度为$O(T\sqrt n)$。
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
int T;
ll n[50010],m[50010],maxn=0;
ll prime[50010],mu[50010],summu[50010];
ll id[50010],d[50010],sumd[50010];
bool vis[50010]={0};

void calc(ll limit)
{
    mu[1]=d[1]=1;
    prime[0]=0;
    for(ll i=2;i<=limit;i++)
    {
        if (!vis[i])
        {
            prime[++prime[0]]=i;
            mu[i]=-1;
            id[i]=2;
            d[i]=2;
        }
        for(ll j=1;j<=prime[0]&&i*prime[j]<=limit;j++)
        {
            vis[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                id[i*prime[j]]=id[i]+1;
                d[i*prime[j]]=d[i]*(id[i]+1ll)/id[i];
                break;
            }
            mu[i*prime[j]]=-mu[i];
            id[i*prime[j]]=2;
            d[i*prime[j]]=d[i]*id[i*prime[j]];
        }
    }
    summu[0]=sumd[0]=0;
    for(ll i=1;i<=limit;i++)
    {
        summu[i]=summu[i-1]+mu[i];
        sumd[i]=sumd[i-1]+d[i];
    }
}

ll solve(ll n,ll m)
{
    ll ans=0;
    for(ll i=min(n,m);i>=1;i=max(n/(n/i+1),m/(m/i+1)))
    {
        ll l=max(n/(n/i+1),m/(m/i+1))+1,r=i;
        ans+=(summu[r]-summu[l-1])*sumd[n/i]*sumd[m/i];
    }
    return ans;
}

int main()
{
    scanf("%d",&T);
    for(int i=1;i<=T;i++)
    {
        scanf("%lld%lld",&n[i],&m[i]);
        maxn=max(maxn,max(n[i],m[i]));
    }

    calc(maxn);
    for(int i=1;i<=T;i++)
        printf("%lld\n",solve(n[i],m[i]));

    return 0;
}