前言
话说这是让我来学莫比乌斯反演的入坑题呢,然而学了莫比乌斯反演还是不会做=_=,连题解都看不懂QAQ
注意事项:
1.数恐症患者慎入!
2.所有除号如果没有作特殊说明,都是向下取整。
3.写的又臭又长,每一步都有一堆啰嗦的解释,因为本蒟蒻每一步都理解了半天(看神犇题解,一个等号看半个小时系列)
4.本蒟蒻发现自己学了假的懵逼钨丝繁衍,所以我重构了一遍这个博客
题目分析
第一步:预备的结论
首先,我们要证明一条结论:
d(nm)=∑i|n∑j|m[gcd(i,j)==1]
怎么证明呢?对于 nm 的每一个质因数 pi ,我们设 n=paii∗n′ , m=pbii∗m′ ,那么
d(nm)=∏i(ai+bi+1)
( ∏ 是累乘的意思)
这条结论意会一下就可以了吧?对了,我们记 pi 对 d(nm) 的贡献为 ai+bi+1 ,“贡献”这个词第三步要用
然后,我们看着要证的那个结论的等号右边部分,设 l 和
啊哈,得证了!
第二步:公式要变形
有了上面的结论,我们就可以玩公式变形了!
∑i=1n∑j=1md(ij)=∑i=1n∑j=1m∑k|i∑l|j[gcd(k,l)==1]
现在我们枚举每一个数对 (k,l) ,那么它们会被算多少遍呢?每一个 k 的倍数都会算一次
∑k=1n∑l=1mnk∗ml[gcd(k,l)==1]
然后大家应该知道莫比乌斯函数有一个神奇的性质(详见: 莫比乌斯反演的讲解):
∑d|nμ(d)=1(n=1)
∑d|nμ(d)=0(n>1)
所以呢:
∑k=1n∑l=1mnk∗ml∑d|gcd(k,l)μ(d)
这样子如果 gcd(i,j)=1 的话后面那团就是1,否则为0.
然后这一步变形非常显然对吧:
∑k=1n∑l=1mnk∗ml∑d|k且d|lμ(d)
现在我们先枚举 μ(d) ,而哪些数会和 μ(d) 一起乘呢?根据上一个式子,显然是 d 的倍数啊,那么是几倍呢?枚举啊!
第三步:转化成代码
已经变形的差不多了,再小小变一下得到最终公式:
∑d=1min(n,m)μ(d)(∑k=1ndndk)(∑l=1mdmdl)
现在我们令 g(x)=∑ii=1xi ,原式就是:
∑d=1min(n,m)μ(d)∗g(nd)∗g(md)
那么我们只要预处理 g(x) 就可以了!
现在观察一下 g(x) ,就会发现其实式子可以这么理解: x 以内有约数
要求 d(x) 的前缀和,可以求出每一个 d(x) ,也就是用线性筛求啦,要记录一下 c(x) :将 x 分解质因数后,最小质因数的指数。
那么参见第一步,
代码
经过刚才含辛茹苦地淬炼,蒟蒻终于写出一段精简的代码:
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define ll long long
int read(){
int q=0;char ch=' ';
while(ch<'0'||ch>'9')ch=getchar();
while(ch>='0'&&ch<='9')q=q*10+ch-'0',ch=getchar();
return q;
}
const int maxn=50001;
int tot,lim=50000,T,n,m;
int pri[30001],mu[maxn],c[maxn];
ll d[maxn];bool is[maxn];
void init(){
int i,j,k;
mu[1]=c[1]=d[1]=1;
for(i=2;i<=lim;i++){
if(!is[i]){
pri[++tot]=i;mu[i]=-1;
c[i]=1;d[i]=2;
}
for(j=1;j<=tot;j++){
k=pri[j]*i;if(k>lim)break;
is[k]=1;
if(i%pri[j]){ d[k]=d[i]*d[pri[j]];c[k]=1;mu[k]=-mu[i];}//有新的最小质因子了
else {d[k]=d[i]/(c[i]+1)*(c[i]+2);c[k]=c[i]+1;break;}
//d:删去原来的贡献,添加新的贡献
}
}
for(i=2;i<=lim;i++)d[i]=d[i]+d[i-1],mu[i]=mu[i]+mu[i-1];//前缀和
}
int main()
{
int i,j;
init();T=read();
while(T--){
n=read();m=read();
if(n>m)swap(n,m);ll ans=0;
for(i=1;i<=n;i=j+1){
j=min(n/(n/i),m/(m/i));//不懂这句请看我的莫比乌斯反演讲解
ans+=(mu[j]-mu[i-1])*d[n/i]*d[m/i];
}
printf("%lld\n",ans);
}
return 0;
}