洛谷P3327：[SDOI2015]约数个数和 （莫比乌斯反演）

题目传送门：https://www.luogu.org/problem/show?pid=3327

---

题目分析：这题我又没有自己想出来……
主要是本题要用到一个很神的结论：

d(ij)=∑x|i∑y|j[(x,y)=1]

这个是怎么推出来的呢？
我们考虑质数p对d(ij)的贡献，假设i的质因数分解中有k个p，j的质因数分解中q个p，那么d(ij)中就会有因数k+q+1，而：

k+q+1=∑x=0k∑y=0q[(px,py)=1]

这就相当于一个(k+1)*(q+1)的矩形，只取了左下角的”L”字型的一块。
现在我们假设i=pk11∗pk22……pknn,j=pq11∗pq22……pqnn，那么我们逐个考虑每一个质数pi，得到一条式子：

d(ij)=∑x1=1k1∑y1=1q1[(px11,py11)=1]……∑xn=1kn∑yn=1qn[(pxnn,pynn)=1]

现在我们记x=px11px22……pxnn,y=py11py22……pynn，则变成：

d(ij)=∑x|i∑y|j[(x,y)=1]

是不是很神奇？

---

现在回到原题，经上面变形，我们要求：

∑i=1n∑j=1m∑x|i∑y|j[(x,y)=1]

我们将x,y的枚举拉到i,j的前面，得到：

∑x=1n∑y=1m[(x,y)=1]⌊nx⌋⌊my⌋

我们用莫比乌斯函数将中括号内的东西替代掉，得到：

∑x=1n∑y=1m∑k|x,k|yμ(k)⌊nx⌋⌊my⌋

将k的枚举拉到前面：

∑k=1min(n,m)μ(k)∑x=1⌊nk⌋⌊nkx⌋∑y=1⌊mk⌋⌊mky⌋

我们设g(n)=∑ni=1⌊ni⌋，那么右边的部分变成g(⌊nk⌋)g(⌊mk⌋)，于是我们左边用下底函数分块，右边预处理g函数即可。
我们可以用O(nn√)通过下底函数分块求出每一个g，但我们发现其实g函数就是约数个数的前缀和，直接线性筛即可。

---

CODE:

#include<iostream>#include<string>#include<cstring>#include<cmath>#include<cstdio>#include<cstdlib>#include<stdio.h>#include<algorithm>using namespace std;const int maxn=50005;typedef long long LL;int sum[maxn];bool vis[maxn];int prime[maxn];int cur=0;LL g[maxn];int t,n,m;LL ans;void Preparation(){    for (int i=1; i<maxn; i++)    {        g[i]=0;        int last;        for (int j=1; j<=i; j=last+1)        {            last=i/(i/j);            g[i]=g[i]+(long long)( i/j*(last-j+1) );        }    }    sum[1]=1;    for (int i=2; i<maxn; i++)    {        if (!vis[i]) sum[i]=-1,prime[++cur]=i;        for (int j=1; j<=cur && i*prime[j]<maxn; j++)        {            int k=i*prime[j];            vis[k]=true;            if (i%prime[j]) sum[k]=-sum[i];            else            {                sum[k]=0;                break;            }        }    }    for (int i=2; i<maxn; i++) sum[i]+=sum[i-1];}int main(){    freopen("c.in","r",stdin);    freopen("c.out","w",stdout);    Preparation();    scanf("%d",&t);    while (t--)    {        scanf("%d%d",&n,&m);        if (n>m) swap(n,m);        ans=0;        int last;        for (int i=1; i<=n; i=last+1)        {            last=min( n/(n/i) , m/(m/i) );            ans=ans+( (long long)(sum[last]-sum[i-1])*g[n/i]*g[m/i] );        }        cout<<ans<<endl;    }    return 0;}