[BZOJ 4407] 于神之怒加强版

题意

求下式的值:

[sum_{i=1}^nsum_{j=1}^mgcd(i,j)^k ]

(n,mle 5 imes 10^6, qle 2000)

题解

[egin{align} f(x)&= sum_i^Nsum_j^M[gcd(i,j)=x] \ F(x)&= sum_{x|m}f(m) \ &=left lfloor frac N x ight floorleft lfloor frac M x ight floor \ ext{Ans}&=sum_d d^k f(d) \ &=sum_dd^ksum_{d|m}F(m)muleft(frac m d ight) \ &= sum_dsum_{d|m}d^k F(m)mu left(frac m d ight) \ &=sum_msum_{d|m} d^kF(m)muleft(frac m d ight) \ &=sum_m F(m)sum_{d|m} d^k mu left (frac m d ight) \ \ end{align} ]

这题目可以有一个推式子的思路: 把所有东西扔进最里层求和号里, 然后就可以对外层为所欲为然后再把一部分东西拽回来.

筛最后的 (g(x)=sumlimits_{d|x}d^kmuleft (frac x d ight)) 的时候, 考虑 (x=p^k) 时, 只有 (p^k) 和 (-p^{k-1}) 会作出贡献, 于是每次加入一个新的质因子的时候直接乘以 (p-1), 如果是旧质因子的话只需要给幂次 (+1) 即可, 乘上一个 (p) 就星了

代码实现

#include <bits/stdc++.h>

const int MOD=1e9+7;
const int MAXN=5e6+10;

int k;
int cnt;
int g[MAXN];
int pr[MAXN];
int pw[MAXN];
int mu[MAXN];
bool npr[MAXN];

void EulerSieve(int);
int Pow(int,int,int);

int main(){
	int T;
	scanf("%d%d",&T,&k);
	EulerSieve(5e6);
	while(T--){
		int n,m;
		scanf("%d%d",&n,&m);
		if(n>m)
			std::swap(n,m);
		int ans=0;
		for(int i=1,j;i<=n;i=j+1){
			j=std::min(n/(n/i),m/(m/i));
			(ans+=1ll*(n/i)*(m/i)%MOD*(g[j]-g[i-1]+MOD)%MOD)%=MOD;
		}
		printf("%d
",ans);
	}
	return 0;
}

void EulerSieve(int n){
	npr[0]=npr[1]=true;
	mu[1]=1;
	g[1]=1;
	for(int i=2;i<=n;i++){
		if(!npr[i]){
			pr[cnt++]=i;
			pw[i]=Pow(i,k,MOD);
			g[i]=(pw[i]-1+MOD)%MOD;
			mu[i]=-1;
		}
		for(int j=0,t;j<cnt&&(t=pr[j]*i)<=n;j++){
			npr[t]=true;
			if(i%pr[j]){
				mu[t]=-mu[i];
				g[t]=1ll*g[i]*g[pr[j]]%MOD;
			}
			else{
				mu[t]=0;
				g[t]=1ll*g[i]*pw[pr[j]]%MOD;
				break;
			}
		}
	}
	for(int i=1;i<=n;i++)
		(g[i]+=g[i-1])%=MOD;
}

inline int Pow(int a,int n,int p){
	int ans=1;
	while(n>0){
		if(n&1)
			ans=1ll*a*ans%p;
		a=1ll*a*a%p;
		n>>=1;
	}
	return ans;
}