[BZOJ 4407] 于神之怒加强版

题意

求下式的值:

\[\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k \]

\(n,m\le 5\times 10^6, q\le 2000\)

题解

\[\begin{align} f(x)&= \sum_i^N\sum_j^M[\gcd(i,j)=x] \\ F(x)&= \sum_{x|m}f(m) \\ &=\left \lfloor \frac N x \right \rfloor\left \lfloor \frac M x \right \rfloor \\ \text{Ans}&=\sum_d d^k f(d) \\ &=\sum_dd^k\sum_{d|m}F(m)\mu\left(\frac m d\right) \\ &= \sum_d\sum_{d|m}d^k F(m)\mu \left(\frac m d \right) \\ &=\sum_m\sum_{d|m} d^kF(m)\mu\left(\frac m d \right) \\ &=\sum_m F(m)\sum_{d|m} d^k \mu \left (\frac m d \right) \\ \\ \end{align} \]

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

筛最后的 \(g(x)=\sum\limits_{d|x}d^k\mu\left (\frac x d \right)\) 的时候, 考虑 \(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\n",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;
}