BZOJ 4407: 于神之怒加强版 莫比乌斯反演 + 线筛积性函数

#include<bits/stdc++.h>
#define ll long long 
#define maxn 5000003 
#define N 5000001  
using namespace std;  
const ll mod=1000000007; 
namespace IO
{
	void setIO(string s)
	{
		string in=s+".in"; 
		freopen(in.c_str(),"r",stdin); 
	}
	void shut()
	{
		fclose(stdin); 
		fclose(stdout); 
	}
};  
ll qpow(ll base,ll o)
{  
	ll tmp=1;          
	while(o) 
	{
		if(o&1)tmp=1ll*tmp*base%mod; 
		base=base*base%mod;            
		o>>=1; 
	}     
	return tmp;  
}
int k,cnt; 
int mu[maxn],vis[maxn],prime[maxn]; 
ll h[maxn],g[maxn];      
void prepare()
{ 
	mu[1]=h[1]=1;   
	for(int i=2;i<=N;++i) 
	{
		if(!vis[i]) g[i]=qpow(i,k), prime[++cnt]=i,mu[i]=-1,h[i]=(g[i]-1+mod)%mod;  
		for(int j=1;j<=cnt&&1ll*prime[j]*i<=N;++j) 
		{
			vis[i*prime[j]]=1; 
			if(i%prime[j]==0) 
			{
				mu[i*prime[j]]=0; 
				h[i*prime[j]]=(h[i]*g[prime[j]])%mod;   
				break; 
			}
			mu[i*prime[j]]=-mu[i]; 
			h[i*prime[j]]=(h[i]*h[prime[j]])%mod; 
		}
	}
	for(int i=1;i<=N;++i) h[i]=(h[i-1]+h[i]+mod)%mod;   
}
ll work(int n,int m)
{
	ll re=0,tmp=0; 
	int i,j; 
	for(i=1;i<=n;i=j+1) 
	{
		j=min(n/(n/i), m/(m/i));       
		tmp=(h[j]-h[i-1]+mod)%mod*(n/i)%mod*(m/i)%mod; 
		re=(re+tmp+mod)%mod; 
	}
	return re;   
}
int main()
{ 
	//IO::setIO("input");   
	int T,n,m; 
	scanf("%d%d",&T,&k); 
	prepare(); 
	while(T--)
	{
		scanf("%d%d",&n,&m); 
		if(n>m) swap(n,m); 
		printf("%lld\n",work(n,m));        
	} 
	//IO::shut();
	return 0; 
}