好虚。。不会做啊
Ans=∑ni=1∑mj=1gcd(i,j)k
,我们考虑枚举公约数
d
,然后令
Ans=∑d=1ndk×∑x=1⌊nd⌋μ(x)⌊ndx⌋⌊mdx⌋
我们令 T=dx ,
则
Ans=∑T=1n⌊nT⌋⌊mT⌋∑d∣Tdkμ(Td)
令 g(T)=∑d∣Tdkμ(Td) ,求出 g(T) 就可以分块辣。
可以发现 g(T) 显然是一个积性函数,然后我们可以进行线性筛。
若 is_prime(x) ,则 g(x)=xk−1
否则,若 gcd(x,p)=1 ,则 g(x×p)=g(x)×g(p) 。
否则 g(x×p)=g(x)×pk 。
前两个是显然的,我们来YY一下第三个。
我们加入一个已经存在的质数 p ,由于含
#include<iostream>
#include<cstdio>
#define P 1000000007
#define ll long long
#define MAXN 5000005
using namespace std;
int prime[MAXN];
ll f[MAXN],p[MAXN];
bool flag[MAXN];
int k;
inline int read()
{
int a=0,f=1; char c=getchar();
while (c<'0'||c>'9') {if (c=='-') f=-1; c=getchar();}
while (c>='0'&&c<='9') {a=a*10+c-'0'; c=getchar();}
return a*f;
}
ll quick_power(ll a,ll b)
{
ll ans=1;
for (ll i=b;i;i>>=1,a=a*a%P)
if (i&1) ans=ans*a%P;
return ans;
}
inline void prepare()
{
f[1]=1;
for (int i=2;i<MAXN;i++)
{
if (!flag[i]) prime[++prime[0]]=i,p[prime[0]]=quick_power(i,k),f[i]=p[prime[0]]-1;
for (int j=1;j<=prime[0]&&i*prime[j]<MAXN;j++)
{
flag[i*prime[j]]=1;
if (i%prime[j]==0)
{
f[i*prime[j]]=f[i]*p[j]%P;
break;
}
else f[i*prime[j]]=f[i]*f[prime[j]]%P;
}
}
for (int i=1;i<MAXN;i++) (f[i]+=f[i-1])%=P;
}
int main()
{
int testcase=read(); k=read();
prepare();
while (testcase--)
{
ll ans=0;
int n=read(),m=read();
if (n>m) swap(n,m);
for (int i=1,pos;i<=n;i=pos+1)
{
pos=min(n/(n/i),m/(m/i));
(ans+=(f[pos]-f[i-1]+P)%P*(n/i)%P*(m/i)%P)%=P;
}
printf("%lld\n",ans);
}
return 0;
}