BZOJ.4407.于神之怒加强版(莫比乌斯反演)

题目链接

---

Description

　　求$$\sum_{i=1}^n\sum_{j=1}m\gcd(i,j)^{K \mod 10}9+7$$

Solution

前面部分依旧套路。

\[\begin{aligned}\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^K&=\sum_{d=1}^{min(n,m)}d^K\sum_{i=1}^n\sum_{j=1}^m\left[(i,j)=d\right]\\&=\sum_{d=1}^{min(n,m)}d^K\sum_{i=1}^{min(\lfloor\frac{n}{d}\rfloor),\lfloor\frac{m}{d}\rfloor)}\mu(i)\lfloor\frac{n}{id}\rfloor\lfloor\frac{m}{id}\rfloor\\&=\sum_{T=1}^{min(n,m)}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\sum_{d\mid T}d^K\mu(\frac{T}{d})\end{aligned} \]

　　这数据范围比较大，后面没法枚举约数求。因为是个狄利克雷卷积的形式，而且\(d^K,\mu\)都是积性函数，所以\(F(T)=\sum_{d\mid T}d^K*\mu(\frac{T}{d})\)也是积性函数，可以线性筛。
　　若\(a,b\)互质，则\(F(ab)=F(a)*F(b)\)；
　　若\(a,b\)不互质，考虑同一个质因数，即\(F(p_i^{k_i})=\mu(1)\times\left(p_i^{k_i}\right)^K+\mu(p_i)\times\left(p_i^{k_i-1}\right)^K\)，后面项的\(\mu()\)都是0了。
　　所以\(F(p_i^{k_i}\times p_i)=F(p_i^{k_i})\times p_i^{K}\)。而当\(p_i,p_j\)不同时，它们互质，直接乘起来就好了。

挂个链接：https://www.cnblogs.com/zwfymqz/p/9337898.html。

---

//49648kb	16248ms
#include <cstdio>
#include <cctype>
#include <algorithm>
//#define gc() getchar()
#define MAXIN 200000
#define gc() (SS==TT&&(TT=(SS=IN)+fread(IN,1,MAXIN,stdin),SS==TT)?EOF:*SS++)
#define MAX 5000000
#define mod (1000000007)
const int N=5e6+5;

int cnt,P[N>>3],mu[N],F[N],pw[N>>3];
bool not_P[N];
char IN[MAXIN],*SS=IN,*TT=IN;

inline int read()
{
	int now=0;register char c=gc();
	for(;!isdigit(c);c=gc());
	for(;isdigit(c);now=now*10+c-'0',c=gc());
	return now;
}
inline int FP(int x,int k)
{
	int t=1;
	for(; k; k>>=1,x=1ll*x*x%mod)
		if(k&1) t=1ll*t*x%mod;
	return t;
}
void Pre(int K)
{
	F[1]=mu[1]=1;
	for(int i=2; i<=MAX; ++i)
	{
		if(!not_P[i])
			P[++cnt]=i, mu[i]=-1, pw[cnt]=FP(i,K), F[i]=(pw[cnt]-1)%mod;
		for(int v,j=1; j<=cnt&&(v=i*P[j])<=MAX; ++j)
		{
			not_P[v]=1;
			if(i%P[j]) mu[v]=-mu[i], F[v]=1ll*F[i]*F[P[j]]%mod;
			else {mu[v]=0, F[v]=1ll*F[i]*pw[j]%mod; break;}
		}
	}
	for(int i=2; i<=MAX; ++i) F[i]+=F[i-1], F[i]>=mod&&(F[i]-=mod);
}

int main()
{
	int T=read(),K=read(); Pre(K);
	for(int n,m; T--; )
	{
		n=read(), m=read();
		long long res=0;
		for(int i=1,nxt,lim=std::min(n,m); i<=lim; i=nxt+1)
		{
			nxt=std::min(n/(n/i),m/(m/i));
			res+=1ll*(F[nxt]-F[i-1])*(n/i)%mod*(m/i)%mod;
		}
		printf("%lld\n",(res%mod+mod)%mod);
	}
	return 0;
}