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;
}