BZOJ 4407 于神之怒加强版

题目链接：于神之怒加强版

　　这个式子还是很妙的，只是我已经思维僵化了

\begin{aligned}
 &\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k \\
=&\sum_{g=1}^ng^k\sum_{i=1}^{\lfloor \frac{n}{g} \rfloor}\sum_{j=1}^{\lfloor \frac{m}{g} \rfloor}\sum_{d|i,d|j}\mu(d) \\
=&\sum_{g=1}^ng^k\sum_{d=1}^{\lfloor \frac{n}{g} \rfloor}\mu(d)\lfloor \frac{n}{dg} \rfloor\lfloor \frac{m}{dg} \rfloor \\
=&\sum_{g=1}^ng^k\sum_{g|x}^{n}\mu(\frac{x}{g}) \lfloor \frac{n}{x} \rfloor\lfloor \frac{m}{x} \rfloor \\
=&\sum_{x=1}^n \lfloor \frac{n}{x} \rfloor\lfloor \frac{m}{x} \rfloor  \sum_{g|x} g^k \mu(\frac{x}{g})
\end{aligned}

　　然后我们可以令\(f(x)=\sum_{g|x} g^k \mu(\frac{x}{g})\)，显然\(f(x)\)是\(id^k\)与\(\mu\)的狄利克雷卷积，所以\(f(x)\)也是一个积性函数

　　由于\(n\)的范围有\(5\times 10^6\)，所以我们来考虑一下如何筛这个\(f\)函数，也就是考虑在\(i\)是\(p\)的倍数的情况下计算\(f(ip)\)(\(p\)为质数)

　　令\(i=p^kq(q\perp p)\)，然后分情况讨论一下：

　　当\(q\ne 1\)时，由于\(i\)的最小质因子为\(p\)，那么\(p^{k+1}<i\)，所以\(f(\frac{i}{p^k})\times f(p^{k+1})\)就是所求；

　　当\(q=1\)时，由于\(f(p^x)=p^{kx}-p^{kx-k}\)，所以\(f(i)\times p^k\)即为所求，\(p^k\)则可以通过预处理所有质数的\(k\)次方解决。

　　最后再分块计算就好了。

　　下面贴代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#define File(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
#define maxn 5000010
#define mod 1000000007

using namespace std;
typedef long long llg;

int T,k,n,m,pr[maxn],lp;
int g[maxn],c[maxn];
bool vis[maxn];
llg f[maxn];

void gi(llg &x){if(x>=mod) x%=mod;}
llg mi(llg a,int b){
	llg s=1;
	while(b){
		if(b&1) s=s*a,gi(s);
		a=a*a,gi(a); b>>=1;
	}
	return s;
}

int main(){
	File("a");
	scanf("%d %d",&T,&k); f[1]=c[1]=1;
	for(int i=2;i<maxn;i++){
		if(!vis[i]) pr[++lp]=i,g[i]=i,c[i]=mi(i,k),f[i]=c[i]-1;
		for(int j=1;i*pr[j]<maxn;j++){
			if(i*pr[j]==31823){
				int aa;
				aa++;
			}
			vis[i*pr[j]]=1;
			if(i%pr[j]) g[i*pr[j]]=pr[j],f[i*pr[j]]=f[i]*f[pr[j]]%mod;
			else{
				f[i*pr[j]]=g[i]!=i?f[i/g[i]]*f[g[i]*pr[j]]:f[i]*c[pr[j]];
				gi(f[i*pr[j]]); g[i*pr[j]]=g[i]*pr[j]; break;
			}
		}
	}
	for(int i=2;i<maxn;i++) f[i]+=f[i-1],gi(f[i]);
	while(T--){
		scanf("%d %d",&n,&m);
		if(n>m) n^=m^=n^=m;
		llg ans=0;
		for(int i=1,nt;i<=n;i=nt+1){
			nt=min(n/(n/i),m/(m/i));
			ans+=(f[nt]-f[i-1]+mod)*(n/i)%mod*(m/i)%mod;
			if(ans>=mod) ans%=mod;
		}
		printf("%lld\n",ans);
	}
	return 0;
}