[bzoj4407] 于神之怒加强版

Description

给出N,M,K.求

Input

输入有多组数据，输入数据的第一行两个正整数T,K，代表有T组数据，K的意义如上所示，下面第二行到第T+1行，每行为两个正整数N,M，其意义如上式所示。

Output

如题

Sample Input

1 2
3 3 

Sample Output

20

HINT

1<=N,M,K<=5000000,1<=T<=2000

前置知识：莫比乌斯反演

莫比乌斯反演，推下式子：

[egin{align} ans&=sum _{d=1}^{min(n,m)}d^k sum _{i=1}^{lfloorfrac{n}{d} floor} sum _{j=1}^{lfloorfrac{m}{d} floor} [gcd(i,j)=1]\ &=sum _{d=1}^{min(n,m)}d^k sum _{i=1}^{lfloorfrac{n}{d} floor} sum _{j=1}^{lfloorfrac{m}{d} floor} sum _{d^prime|i&d^prime|j} mu(d^prime)\ &=sum _{d=1}^{min(n,m)}d^k sum _{d^prime} mu(d^prime)lfloorfrac{n}{dd^prime} floor lfloorfrac{m}{dd^prime} floor \ &=sum _{d=1}^{min(n,m)}d^k sum _{d|T} mu(frac{T}{d})lfloorfrac{n}{T} floor lfloorfrac{m}{T} floor \ &=sum _{T} lfloorfrac{n}{T} floor lfloorfrac{m}{T} floor sum _{d|T}d^k mu(frac{T}{d})\ end {align} ]

设(f)为：

[f(n)=sum_{d|n}d^kmu(frac{n}{d}) ]

然后我们可以线筛出(mu)然后大力算(f)然后数论分块，由于调和级数，复杂度为(O(nlog(n)+qsqrt{n}))。

然而交一发T掉了，，，可以考虑线筛出(f)。

注意到(f)是(mu)和(g(x)=x^k)的狄利克雷卷积，可得(f)是个积性函数。

[f(p^a)=sum _{i=0}^{a}(p^i)^kmu(p^{a-i})=p^{ak}-p^{(a-1)k} ]

然后若(p|n)，可得(f(np)=f(n)cdot p^k)，否则可得(f(np)=f(n)cdot(p^k-1))。

然后线筛即可，复杂度(O(n+qsqrt{n}))。

暴力算(f)的代码（TLE）：

#include<bits/stdc++.h>
using namespace std;

#define int long long 

void read(int &x) {
	x=0;int f=1;char ch=getchar();
	for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
	for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}

void print(int x) {
	if(x<0) putchar('-'),x=-x;
	if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('
');}

const int maxn = 5e6+1;
const int mod = 1e9+7;

int pri[maxn],vis[maxn],mu[maxn],f[maxn],tot,n,m,k;

int qpow(int a,int x) {
	int res=1;
	for(;x;x>>=1,a=a*a%mod) if(x&1) res=res*a%mod;
	return res;
}

void sieve() {
	mu[1]=1;
	for(int i=2;i<=n;i++) {
		if(!vis[i]) pri[++tot]=i,mu[i]=-1;
		for(int j=1;j<=tot&&i*pri[j]<=n;j++) {
			vis[i*pri[j]]=1;
			if(!(i%pri[j])) {mu[i*pri[j]]=0;break;}
			mu[i*pri[j]]=-mu[i];
		}
	}
	for(int d=1;d<=n;d++) {
		int res=qpow(d,k);
		for(int i=1;i*d<=n;i++) (f[i*d]+=res*mu[i])%=mod;
	}
	for(int i=1;i<=n;i++) f[i]=(f[i]+f[i-1])%mod;
}

signed main() {
	int t;n=maxn-1;read(t),read(k);sieve();
	while(t--) {
		read(n),read(m);
		int T=1,ans=0;
		while(T<=n&&T<=m) {
			int pre=T;T=min(n/(n/T),m/(m/T));
			ans=(ans+(n/T)*(m/T)%mod*(f[T]-f[pre-1])%mod)%mod;
			T++;
		}
		write((ans%mod+mod)%mod);
	}
	return 0;
}

线筛(f)：

#include<bits/stdc++.h>
using namespace std;

#define int long long 

void read(int &x) {
	x=0;int f=1;char ch=getchar();
	for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
	for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}

void print(int x) {
	if(x<0) putchar('-'),x=-x;
	if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('
');}

const int maxn = 5e6+1;
const int mod = 1e9+7;

int pri[maxn],vis[maxn],mu[maxn],f[maxn],tot,n,m,k,p[maxn];

int qpow(int a,int x) {
	int res=1;
	for(;x;x>>=1,a=a*a%mod) if(x&1) res=res*a%mod;
	return res;
}

void sieve() {
	f[1]=1;
	for(int i=2;i<=n;i++) {
		if(!vis[i]) pri[++tot]=i,p[tot]=qpow(i,k),f[i]=p[tot]-1;
		for(int j=1;j<=tot&&i*pri[j]<=n;j++) {
			vis[i*pri[j]]=1;
			if(!(i%pri[j])) {f[i*pri[j]]=f[i]*p[j]%mod;break;}
			f[i*pri[j]]=f[i]*(p[j]-1)%mod;
		}
	}for(int i=1;i<=n;i++) f[i]=(f[i]+f[i-1])%mod;
}

signed main() {
	int t;n=maxn-1;read(t),read(k);
	//int PRE=clock();
	sieve();
	//cerr << (double) (clock()-PRE)/CLOCKS_PER_SEC << endl;
	while(t--) {
		read(n),read(m);
		int T=1,ans=0;
		while(T<=n&&T<=m) {
			int pre=T;T=min(n/(n/T),m/(m/T));
			ans=(ans+(n/T)*(m/T)%mod*(f[T]-f[pre-1])%mod)%mod;
			T++;
		}
		write((ans%mod+mod)%mod);
	}
	return 0;
}

注意下bzoj不能用clock函数，否则狂RE不止，别问我怎么知道的QAQ