题目描述
传送门
题目大意:
∑ni=1∑mj=1gcd(i,j)k mod 109+7
题解
一道非常不错的莫比乌斯反演,想了一晚啊。。。。
首先引入反演的两个公式
(1)如果
F(n)=∑d|nf(n)
, 那么
f(n)=∑d|nμ(d)F(nd)
(2)如果
F(n)=∑n|df(d)
,那么
f(n)=∑n|dμ(dn)∗F(d)
那么反演有什么用呢?可以实现F(n),f(n)之间的转换求解,使不易求的式子转化成易求的形式。
最常见的应用就是
[n=1]=∑d|nμ(d)
那么对于这道题来说我们先对式子进行改动,设
n<=m
∑d=1ndk∑i=1n∑j=1m[gcd(i,j)=d]
设 f(d)=∑ni=1∑mj=1[gcd(i,j)=d] 表示满足 gcd(i,j)=d 且 1<=i<=n,1<=j<=m 的 (i,j) 的对数
F(d)=∑d|nf(n) 表示满足 d|gcd(i,j) 且 1<=i<=n,1<=j<=m 的 (i,j) 的对数
显然 F(d)=⌊nd⌋∗⌊md⌋
利用反演公式(2),可以得到
∑d=1ndk∑i=1nμ(i)∗F(d∗i)
注意 F(d)=∑d|nf(n) 中n的范围是无穷,但是如果n超过给出的min(n,m),那么F(d)就会是0
所以对于枚举的d来说,最大的倍数只能是min(n,m).
∑d=1ndk∑i=1nμ(i)∗⌊nd∗i⌋∗⌊md∗i⌋
设 T=d∗i ,然后我们调整枚举的顺序
∑T=1n⌊nT⌋∗⌊mT⌋∑d|Tμ(Td)∗dk
⌊nT⌋∗⌊mT⌋ 可以在 O(n√+m−−√) 的时间内求解,那么问题的关键就是如果快速的计算 h(d)=∑d|Tμ(Td)∗dk 发现函数 h 是狄利克雷卷积的形式。
所以当a,b互质的时候,我们可以利用积性函数的性质 h(ab)=h(a)∗h(b)
关键是a,b不互质的时候怎么求解。
h(pixi)=μ(1)∗(pixi)k+μ(pi)∗(pixi−1)k 剩下的项 μ 都是0
那么对于一个质因数来书, h(pixi) 变成 h(pixi+1) 的影响就是在 h(pixi) 的基础上乘 pik ,其他质因数与其互质,直接乘起来就好啦。
那么对于函数h我们就可用线性筛进行预处理,利用的时候直接 O(1)
代码
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define N 5000000
#define p 1000000007
#define LL long long
using namespace std;
LL f[N+3],g[N+3];
int n,m,k,T,pd[N+3],prime[N+3];
LL quickpow(LL num,int x)
{
LL ans=1; LL base=num%p;
while (x) {
if (x&1) ans=base*ans%p;
x>>=1;
base=base*base%p;
}
return ans;
}
void init()
{
f[1]=1;
for (int i=2;i<=N;i++) {
if (!pd[i]) {
prime[++prime[0]]=i;
g[prime[0]]=quickpow(i,k);
f[i]=(g[prime[0]]-1+p)%p;
}
for (int j=1;j<=prime[0];j++) {
if (i*prime[j]>N) break;
pd[i*prime[j]]=1;
if (i%prime[j]==0) {
f[i*prime[j]]=f[i]*g[j]%p;
break;
}
else f[i*prime[j]]=f[i]*f[prime[j]]%p;
}
}
for (int i=1;i<=N;i++) f[i]=(f[i-1]+f[i])%p;
//for (int i=1;i<=10;i++) cout<<f[i]<<" ";
//cout<<endl;
}
int main()
{
freopen("a.in","r",stdin);
scanf("%d%d",&T,&k);
init();
while (T--) {
scanf("%d%d",&n,&m);
if (n>m) swap(n,m);
int j=0; LL ans=0;
for (int i=1;i<=n;i=j+1) {
j=min(n/(n/i),m/(m/i));
j=min(j,n);
LL t=(LL)(n/i)*(m/i)%p;
ans=(ans+(LL)t*(f[j]-f[i-1]+p)%p)%p;
}
printf("%I64d\n",(ans%p+p)%p);
}
}