[BZOJ4407]于神之怒加强版

Description
给下N,M,K.求

$$\sum\limits_{i=1}^n\sum\limits_{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，然后继续推柿子

$$\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$$

然后我们令$T=dx$，柿子变为

$$\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)$及其前缀和，我们便又可以开心地分块了

由于$d^k$是积性函数，$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)$$

然后就可以线筛了$\sim\sim\sim$

/*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\n",Ans);
	}
	return 0;
}