[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;
}
作者:Wolfycz
本文版权归作者和博客园共有,欢迎转载,但必须在文章开头注明原文出处,否则保留追究法律责任的权利
posted @ 2018-08-16 09:23  Wolfycz  阅读(362)  评论(0编辑  收藏  举报