题目
题目链接：https://www.luogu.com.cn/problem/P4449
给定 (n,m,k)，计算
[sum_{i=1}^n sum_{j=1}^m gcd(i,j)^k
]
对 (10^9+7) 取模的结果。
思路
菜到真的只会模板题了 /kk。
[sum_{i=1}^n sum_{j=1}^m gcd(i,j)^k
]
[=sum^{n}_{d=1}(sum^n_{d|i}mu(frac{i}{d})lfloor frac{n}{i}
floor lfloor frac{m}{i}
floor imes d^k)
]
[=sum^{n}_{i=1}lfloor frac{n}{i}
floor lfloor frac{m}{i}
floor(sum_{d|i}mu(frac{i}{d}) imes d^k)
]
前面整除分块，后面预处理即可。
时间复杂度 (O(nlog n+Tsqrt{n}))。
代码
#include <bits/stdc++.h>
#define reg register
using namespace std;
typedef long long ll;
const int N=5000010,MOD=1e9+7;
int Q,cnt,n,m,t,mu[N],prm[N],sum[N];
ll ans;
bool v[N];
ll fpow(ll x,ll k)
{
ll ans=1;
for (;k;k>>=1,x=x*x%MOD)
if (k&1) ans=ans*x%MOD;
return ans;
}
void findprm(int n)
{
mu[1]=1;
for (reg int i=2;i<=n;i++)
{
if (!v[i]) prm[++cnt]=i,mu[i]=-1;
for (reg int j=1;j<=cnt;j++)
{
if (i>n/prm[j]) break;
v[i*prm[j]]=1; mu[i*prm[j]]=-mu[i];
if (i%prm[j]==0)
{
mu[i*prm[j]]=0;
break;
}
}
}
}
int main()
{
scanf("%d%d",&Q,&t);
findprm(N-10);
for (reg int i=1;i<=N-10;i++)
{
ll power=fpow(i,t);
for (reg int j=i;j<=N-10;j+=i)
sum[j]=(sum[j]+power*mu[j/i])%MOD;
}
for (reg int i=1;i<=N-10;i++)
sum[i]=(sum[i]+sum[i-1])%MOD;
while (Q--)
{
ans=0;
scanf("%d%d",&n,&m);
for (reg int l=1,r;l<=min(n,m);l=r+1)
{
r=min(n/(n/l),m/(m/l));
ans=(ans+1LL*(n/l)*(m/l)%MOD*(sum[r]-sum[l-1]))%MOD;
}
printf("%lld
",(ans%MOD+MOD)%MOD);
}
return 0;
}