Description
给下N,M,K.求
[sumlimits_{i=1}^nsumlimits_{j=1}^m gcd(i,j)^k \%(10^9+7)
]
Input
输入有多组数据，输入数据的第一行两个正整数T,K，代表有T组数据，K的意义如上所示，下面第二行到第T+1行，每行为两个正整数N,M，其意义如上式所示。
Output
如题
Sample Input
1 2
3 3
Sample Output
20
HINT
1<=N,M,K<=5000000,1<=T<=2000
我们令n<m，然后继续推柿子
[sumlimits_{i=1}^nsumlimits_{j=1}^m gcd(i,j)^k
]
[sumlimits_{i=1}^nsumlimits_{j=1}^msumlimits_{d=1}^n d^k [gcd(i,j)=1]
]
[sumlimits_{d=1}^n d^ksumlimits_{i=1}^{lfloorfrac{n}{d}
floor}sumlimits_{j=1}^{lfloorfrac{m}{d}
floor}sumlimits_{x|i,x|j}mu(x)
]
[sumlimits_{d=1}^n d^ksumlimits_{x=1}^{lfloorfrac{n}{d}
floor}mu(x)lfloordfrac{n}{dx}
floorlfloordfrac{n}{dx}
floor
]
然后我们令(T=dx)，柿子变为
[sumlimits_{T=1}^nlfloordfrac{n}{T}
floorlfloordfrac{m}{T}
floorsumlimits_{d|T}d^kmu(dfrac{T}{d})
]
令(g(T)=sumlimits_{d|T}d^kmu(dfrac{T}{d}))，因此我们只要预处理出(g(T))及其前缀和，我们便又可以开心地分块了
由于(d^k)是积性函数，(g(T))也是积性函数，所以有
[g(T)=prodlimits_{i=1}^t g(P_i^{x_i})
]
[=prodlimits_{i=1}^t(P_i^{k imes (x_i-1)} imes mu(P_i)+P_i^{k imes x_i} imes mu(1))
]
[=prodlimits_{i=1}^t P_i^{k imes (x_i-1)} imes(P_i^k-1)
]
然后就可以线筛了(simsimsim)
/*program from Wolfycz*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define inf 0x7f7f7f7f
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
inline int read(){
int x=0,f=1;char ch=getchar();
for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;
for (;ch>='0'&&ch<='9';ch=getchar()) x=(x<<1)+(x<<3)+ch-'0';
return x*f;
}
inline void print(int x){
if (x>=10) print(x/10);
putchar(x%10+'0');
}
const int N=5e6,p=1e9+7;
int prime[N+10],miu[N+10],g[N+10],f[N+10];
bool inprime[N+10];
int Data,k,tot;
int mlt(int a,int b){
int res=1;
for (;b;b>>=1,a=1ll*a*a%p) if (b&1) res=1ll*res*a%p;
return res;
}
void prepare(){
f[1]=1;
for (int i=2;i<=N;i++){
if (!inprime[i]){
prime[++tot]=i;
g[tot]=mlt(i,k);
f[i]=(g[tot]-1+p)%p;
}
for (int j=1;j<=tot&&i*prime[j]<=N;j++){
inprime[i*prime[j]]=1;
if (i%prime[j]==0){
f[i*prime[j]]=1ll*f[i]*g[j]%p;
break;
}
f[i*prime[j]]=1ll*f[i]*f[prime[j]]%p;
}
}
for (int i=1;i<=N;i++) f[i]=(f[i]+f[i-1])%p;
}
int main(){
Data=read(),k=read();
prepare();
while (Data--){
int n=read(),m=read(),T=min(n,m),pos,Ans=0;
for (int i=1;i<=T;i=pos+1){
pos=min(n/(n/i),m/(m/i));
Ans=(Ans+1ll*(f[pos]-f[i-1]+p)*(n/i)%p*(m/i)%p)%p;
}
printf("%d
",Ans);
}
return 0;
}