[莫比乌斯反演 积性函数前缀和] BZOJ 4407 于神之怒加强版

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;

inline char nc()
{
	static char buf[100000],*p1=buf,*p2=buf;
	if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }
	return *p1++;
}

inline void read(ll &x)
{
	char c=nc(),b=1;
	for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;
	for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;
}

const int P=1e9+7;
const int N=5000000;

inline ll Pow(ll a,ll b){
	ll ret=1;
	for (;b;b>>=1,a=a*a%P)
		if (b&1)
			ret=ret*a%P;
	return ret;
}

ll T,K;
ll n[2005],m[2005],maxn;

int prime[N+5],vst[N+5],tot;
ll sum[N+5];

inline void Init(){
	sum[1]=1;
	for (int i=2;i<=maxn;i++)
	{
		if (!vst[i]){
			prime[++tot]=i; sum[i]=Pow(i,K)-1; sum[i]=(sum[i]+P)%P;
		}
		for (int j=1;(ll)prime[j]*i<=maxn && j<=tot;j++)
		{
			vst[i*prime[j]]=1;
			if (i%prime[j]==0){
				sum[i*prime[j]]=sum[i]*Pow(prime[j],K)%P;
				break;
			}
			sum[i*prime[j]]=sum[i]*sum[prime[j]]%P;
		}
	}
	for (int i=2;i<=maxn;i++)
		(sum[i]+=sum[i-1])%=P;
}

inline void Solve(ll n,ll m){
	if (n>m) swap(n,m);
	int last=0; ll ans=0;
	for (int i=1;i<=n;i=last+1)
	{
		last=min(n/(n/i),m/(m/i));
		ans+=(n/i)*(m/i)%P*(sum[last]-sum[i-1]+P)%P;
		ans=(ans+P)%P;
	}
	printf("%lld\n",ans);
}

int main()
{
	freopen("t.in","r",stdin);
	freopen("t.out","w",stdout);
	read(T); read(K);
	for (int i=1;i<=T;i++) read(n[i]),read(m[i]),maxn=max(maxn,min(n[i],m[i]));
	Init();
	for (int i=1;i<=T;i++)
		Solve(n[i],m[i]);
	return 0;
}