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

题意

设d(x)为x的约数个数，给定N、M，求∑ni=1∑mj=1d(ij)

n,m,T<=50000

分析

一开始推了半天还是停留在原来的式子，于是无奈地点开了题解。。。

首先模一下PoPoQQQ大佬%%%

要把这个讨厌的d(ij)给去掉的话，怎么换元或交换主题都没用，得先知道一个神奇的结论：

d(nm)=∑i|n∑j|m[gcd(i,j)==1]

我们来试着用归纳法证明一下：

显然当n=1,m=1时原式成立

考虑每个素数p的贡献

设m=m′∗pk1,n=n′∗pk2

若当m=m′,n=n′时等式成立

显然p对左边式子的贡献为(k1+k2+1)

p对右边式子的贡献为：满足条件的数对有(pk1,1),(pk1−1,1)...(p,1),(1,1),(1,p)...(1,pk2)

恰好有(k1+k2+1)对，那么对右边的式子的贡献也为(k1+k2+1)，那么等式成立，证毕。

然后我们就可以开始愉快的推式子啦！

那么∑ni=1∑mj=1d(ij)=∑ni=1∑mj=1⌊ni⌋⌊mj⌋[gcd(i,j)==1]

然后一波推之后就可以得到ans=∑nd=1μ(d)∗f(⌊nd⌋)∗f(⌊md⌋)

f(d)=∑di=1⌊di⌋

上板子即可。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#define N 50005
#define LL long long
using namespace std;

int mu
,tot,prime
,f
;
bool not_prime
;

void get_prime(int n)
{
mu[1]=1;
for (int i=2;i<=n;i++)
{
if (!not_prime[i]) prime[++tot]=i,mu[i]=-1;
for (int j=1;j<=tot&&i*prime[j]<=n;j++)
{
not_prime[i*prime[j]]=1;
if (i%prime[j]==0)
{
mu[i*prime[j]]=0;
break;
}
mu[i*prime[j]]=-mu[i];
}
}
for (int i=1;i<=n;i++) mu[i]+=mu[i-1];
for (int d=1;d<=n;d++)
for (int i=1,last;i<=d;i=last+1)
{
last=d/(d/i);
f[d]+=(last-i+1)*(d/i);
}
}

LL solve(int n,int m)
{
if (n>m) swap(n,m);
LL ans=0;
for (int i=1,last;i<=n;i=last+1)
{
last=min(n/(n/i),m/(m/i));
ans+=(LL)(mu[last]-mu[i-1])*f[n/i]*f[m/i];
}
return ans;
}

int main()
{
get_prime(50000);
int T;
scanf("%d",&T);
while (T--)
{
int n,m;
scanf("%d%d",&n,&m);
printf("%lld\n",solve(n,m));
}
return 0;
}