# [BZOJ4407]于神之怒加强版

Description

$\sum\limits_{i=1}^n\sum\limits_{j=1}^m \gcd(i,j)^k \%(10^9+7)$

Input

Output

Sample Input
1 2
3 3

Sample Output
20

HINT
1<=N,M,K<=5000000,1<=T<=2000

$\sum\limits_{i=1}^n\sum\limits_{j=1}^m \gcd(i,j)^k$

$\sum\limits_{i=1}^n\sum\limits_{j=1}^m\sum\limits_{d=1}^n d^k [\gcd(i,j)=1]$

$\sum\limits_{d=1}^n d^k\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sum\limits_{x|i,x|j}\mu(x)$

$\sum\limits_{d=1}^n d^k\sum\limits_{x=1}^{\lfloor\frac{n}{d}\rfloor}\mu(x)\lfloor\dfrac{n}{dx}\rfloor\lfloor\dfrac{n}{dx}\rfloor$

$\sum\limits_{T=1}^n\lfloor\dfrac{n}{T}\rfloor\lfloor\dfrac{m}{T}\rfloor\sum\limits_{d|T}d^k\mu(\dfrac{T}{d})$

$g(T)=\sum\limits_{d|T}d^k\mu(\dfrac{T}{d})$，因此我们只要预处理出$g(T)$及其前缀和，我们便又可以开心地分块了

$g(T)=\prod\limits_{i=1}^t g(P_i^{x_i})$

$=\prod\limits_{i=1}^t(P_i^{k\times (x_i-1)}\times \mu(P_i)+P_i^{k\times x_i}\times \mu(1))$

$=\prod\limits_{i=1}^t P_i^{k\times (x_i-1)}\times(P_i^k-1)$

/*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;
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(){
prepare();
while (Data--){