【BZOJ4407】于神之怒加强版

题面

题目分析

\[\begin{split} \sum\limits_{i=1}^n\sum\limits_{j=1}^mgcd(i,j)^k&=\sum\limits_{d=1}^nd^k\sum\limits_{i=1}^n\sum\limits_{j=1}^m[gcd(i,j)==d]\\ \end{split} \]

设\(f(x)\)表示\(gcd(i,j)=x\)，\(g(x)\)表示\(gcd(i,j)==kx,k\in Z\)。

\[\begin{split} g(x)&=\sum\limits_{x|d}^nf(d)\\ &=\sum\limits_{i=1}^n\sum\limits_{j=1}^m[x|gcd(i,j)]\\ &=\sum\limits_{i=1}^{\lfloor\frac n x\rfloor}\sum\limits_{j=1}^{\lfloor\frac m x\rfloor}\lfloor\frac n x\rfloor\lfloor\frac m x\rfloor\\ f(x)&=\sum\limits_{x|d}^n\mu(\frac dx)g(d)=\sum\limits_{x|d}^n\mu(\frac dx)\lfloor\frac n d\rfloor\lfloor\frac m d\rfloor \end{split} \]

\[\begin{split} ans&=\sum\limits_{d=1}^nd^k\cdot f(d)\\ &=\sum\limits_{d=1}^nd^k\sum\limits_{d|T}^n\mu(\frac Td)\lfloor\frac n T\rfloor\lfloor\frac m T\rfloor\\ &=\sum\limits_{T=1}^n\lfloor\frac n T\rfloor\lfloor\frac m T\rfloor\sum\limits_{d|T}\mu(\frac Td)d^k \end{split} \]

由于\(\mu\)和\(d^k\)均为积性函数，所以\(\sum\limits_{d|T}\mu(\frac Td)d^k\)也为积性函数，可以在线性筛中\(O(n\log n)\)预处理。

前面部分用整除分块加速。

代码实现

#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdio>
#include<iomanip>
#include<cstdlib>
#define MAXN 0x7fffffff
typedef long long LL;
const int N=5000005,mod=1e9+7;
using namespace std;
inline int Getint(){register int x=0,f=1;register char ch=getchar();while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}return x*f;}
int g[N],mu[N],prime[N];
bool vis[N];
LL ksm(LL x,LL k){
	LL ret=1;
	while(k){
		if(k&1)ret=ret*x%mod;
		x=x*x%mod;
		k>>=1;
	}
	return ret;
}
int low[N];
int main(){
	int T=Getint(),K=Getint();
	
	mu[1]=g[1]=1;
	for(int i=2;i<=5e6;i++){
		if(!vis[i]){
			prime[++prime[0]]=i,mu[i]=-1;
			low[i]=i,g[i]=ksm(i,K)-1;
		}
		for(int j=1;j<=prime[0]&&1ll*prime[j]*i<=5e6;j++){
			vis[i*prime[j]]=1;
			if(i%prime[j]==0){
				low[i*prime[j]]=low[i]*prime[j];
				if(low[i*prime[j]]==i*prime[j])
					g[i*prime[j]]=g[i]*ksm(prime[j],K)%mod;
				else 
					g[i*prime[j]]=(1ll*g[low[i*prime[j]]]*g[i*prime[j]/low[i*prime[j]]])%mod;
				break;
			}
			low[i*prime[j]]=prime[j];
			g[i*prime[j]]=(1ll*g[i]*g[prime[j]])%mod;
			mu[i*prime[j]]=-mu[i];
		}
	}
	for(int i=1;i<=5e6;i++)g[i]=(g[i]+g[i-1])%mod;
	
	while(T--){
		int n=Getint(),m=Getint();
		if(n>m)swap(n,m);
		int ans=0;
		for(int l=1,r;l<=n;l=r+1){
			r=min(n/(n/l),m/(m/l));
			ans=(ans+1ll*(n/l)*(m/l)%mod*(g[r]-g[l-1])%mod+mod)%mod; 
		}
		cout<<ans<<'\n';
	}
	return 0;
}