BZOJ 4407 于神之怒加强版

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

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

egin{aligned}
 &sum_{i=1}^nsum_{j=1}^mgcd(i,j)^k \
=&sum_{g=1}^ng^ksum_{i=1}^{lfloor frac{n}{g} floor}sum_{j=1}^{lfloor frac{m}{g} floor}sum_{d|i,d|j}mu(d) \
=&sum_{g=1}^ng^ksum_{d=1}^{lfloor frac{n}{g} floor}mu(d)lfloor frac{n}{dg} floorlfloor frac{m}{dg} floor \
=&sum_{g=1}^ng^ksum_{g|x}^{n}mu(frac{x}{g}) lfloor frac{n}{x} floorlfloor frac{m}{x} floor \
=&sum_{x=1}^n lfloor frac{n}{x} floorlfloor frac{m}{x} floor  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 imes 10^6)，所以我们来考虑一下如何筛这个(f)函数，也就是考虑在(i)是(p)的倍数的情况下计算(f(ip))((p)为质数)

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

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

　　当(q=1)时，由于(f(p^x)=p^{kx}-p^{kx-k})，所以(f(i) imes 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
",ans);
	}
	return 0;
}