3994: [SDOI2015]约数个数和 莫比乌斯反演

去年省选不会做，过了一年又要省选了还是不会做，滚粗既视感。

---

一维的时候我们知道，$\sum_{i=1}^nd(i)=\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor$，然后这个可以拓展到二维，然而我并不会证明

$$ans=\sum_{i=1}^n\sum_{j=1}^m\lfloor\frac{n}{i}\rfloor\lfloor\frac{m}{j}\rfloor,gcd(i,j)=1$$
然后莫比乌斯反演。
$$ans=\sum_{i=1}^n\sum_{j=1}^m\lfloor\frac{n}{i}\rfloor\lfloor\frac{m}{j}\rfloor\sum_{d\mid gcd(i,j)}\mu(d)$$
$$=\sum_{d=1}^n\mu(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\lfloor\frac{n}{id}\rfloor\lfloor\frac{m}{jd}\rfloor$$
$$=\sum_{d=1}^n\mu(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{id}\rfloor\sum_{i=1}^{\lfloor\frac{m}{d}\rfloor}\lfloor\frac{n}{jd}\rfloor$$
$$=\sum_{d=1}^n\mu(d)f(\lfloor\frac{n}{d}\rfloor)f(\lfloor\frac{m}{d}\rfloor)$$
这里的$f(x)=\sum_{i=1}^x\lfloor\frac{x}{i}\rfloor$，可以发现就是$d(x)$的前缀和。
我们可以线性筛除$d(x)$然后就有了$f(x)$然后就分块就好辣。

#include<iostream>
#include<cstdio>
#define MAXN 50005
#define ll long long 
using namespace std;
int mobius[MAXN],prime[MAXN],d[MAXN],f[MAXN];
bool flag[MAXN];
inline ll read()
{
    ll a=0,f=1; char c=getchar();
    while (c<'0'||c>'9') {if (c=='-') f=-1; c=getchar();}
    while (c>='0'&&c<='9') {a=a*10+c-'0'; c=getchar();}
    return a*f;
}
inline void prepare()
{
    mobius[1]=d[1]=1;
    for (int i=2;i<MAXN;i++)
    {
        if (!flag[i]) prime[++prime[0]]=i,mobius[i]=-1,d[i]=2,f[i]=1;
        for (int j=1;j<=prime[0]&&i*prime[j]<MAXN;j++)
        {
            flag[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                mobius[i*prime[j]]=0;
                f[i*prime[j]]=f[i]+1;
                d[i*prime[j]]=d[i]/(f[i]+1)*(f[i]+2);
                break;
            }
            mobius[i*prime[j]]=-mobius[i];
            f[i*prime[j]]=1;
            d[i*prime[j]]=d[i]*2;
        }
    }
    for (int i=1;i<MAXN;i++) d[i]+=d[i-1],mobius[i]+=mobius[i-1];
}
int main()
{
    prepare();
    int testcase=read();
    while (testcase--)
    {
        ll n=read(),m=read();
        if (n>m) swap(n,m);
        ll ans=0;
        for (int i=1,pos;i<=n;i=pos+1)
        {
            pos=min(n/(n/i),m/(m/i));
            ans+=(ll)(mobius[pos]-mobius[i-1])*d[n/i]*d[m/i];
        }
        printf("%lld\n",ans);
    }
    return 0;
}