【BZOJ4407】于神之怒加强版（莫比乌斯反演）

点此看题面

大致题意： 求\(\sum_{i=1}^n\sum_{j=1}^mgcd(i,j)^k\)。

前言

\(Jan\ 28th\)刷题计划（3/6），算法标签：莫比乌斯反演&杜教筛。

学了一个新技巧：线性筛筛自定义积性函数。

莫比乌斯反演

首先按照套路去枚举\(gcd\)得到：

\[\sum_{d=1}^{min(n,m)}d^k\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}[gcd(i,j)==1] \]

莫比乌斯反演，得到：

\[\sum_{d=1}^{min(n,m)}d^k\sum_{p=1}^{min(\lfloor\frac nd\rfloor,\lfloor\frac md\rfloor)}\mu(p)\lfloor\frac n{dp}\rfloor\lfloor\frac m{dp}\rfloor \]

考虑枚举\(D=dp\)，得到：

\[\sum_{D=1}^{min(n,m)}\lfloor\frac nD\rfloor\lfloor\frac mD\rfloor\sum_{d|D}d^k\mu(\frac Dd) \]

设\(f(D)=\sum_{d|D}d^k\mu(\frac Dd)\)，则原式等同于：

\[\sum_{D=1}^{min(n,m)}\lfloor\frac nD\rfloor\lfloor\frac mD\rfloor f(D) \]

则只要求出\(f(D)\)，就可以通过除法分块轻松解决此题。

于是关键问题来了：如何求出\(f(D)\)。

线性筛筛积性函数

显然\(f(n)\)是一个积性函数。

如果我们设\(n=\prod_{i=1}^tP_i^{a_i}\)（\(P\)为质数），那么：

\[f(n)=\prod_{i=1}^tf(P_i^{a_i}) \]

考虑化简\(f(P_i^{a_i})\)：

\[f(P_i^{a_i})=\sum_{d|P_i^{a_i}}d^k\mu(\frac Dd)=1^k\times\mu(P_i^{a_i})+P_i^k\times\mu(P_i^{a_i-1})+...+P_i^{a_i\times k}\times\mu(1) \]

根据\(\mu\)的定义可知，\(\mu(P_i^{a_i}),\mu(P_i^{a_i-1}),...,\mu(P_i^2)\)都等于\(0\)，也就是说，我们只需要考虑这个式子的最后两项。即：

\[f(P_i^{a_i})=P_i^{(a_i-1)\times k}\times\mu(P_i)+P_i^{a_i\times k}\times\mu(1)=-P_i^{(a_i-1)\times k}+P_i^{a_i\times k}=P_i^{(a_i-1)\times k}(P_i^k-1) \]

也就是说，有：

\[f(n)=\prod_{i=1}^tP_i^{(a_i-1)\times k}(P_i^k-1) \]

然后就很好做了。

以线性筛为例，我们可以得到边界状态为：

\[f(n)=\begin{cases}1&(n=1)\\p^k-1&(n\ is\ prime)\end{cases} \]

然后考虑筛的过程中的筛法为：

\[f(i\times P_j)=\begin{cases}f(i)\times f(P_j)&(P_j\not{|}i)\\f(i)\times P_i^k&(P_j|i)\end{cases} \]

这可以根据前面推出的式子自行理解。

代码

#include<bits/stdc++.h>
#define Tp template<typename Ty>
#define Ts template<typename Ty,typename... Ar>
#define Reg register
#define RI Reg int
#define Con const
#define CI Con int&
#define I inline
#define W while
#define N 10000000
#define X 1000000007
using namespace std;
int n,m,k;
I int Qpow(RI x,RI y) {RI t=1;W(y) y&1&&(t=1LL*t*x%X),x=1LL*x*x%X,y>>=1;return t;}
template<int SZ> class LinearSieve//线性筛
{
	private:
		int Pt,P[SZ+5],f[SZ+5],p[SZ+5];
	public:
		I int operator [] (CI x) Con {return f[x];}
		I void Init()
		{
			RI i,j;for(f[1]=1,i=2;i<=SZ;++i)
			{
				!P[i]&&(P[++Pt]=i,p[i]=Qpow(i,k),f[i]=p[i]-1);
				for(j=1;j<=Pt&&i*P[j]<=SZ;++j)
					if(P[i*P[j]]=1,i%P[j]) f[i*P[j]]=1LL*f[i]*f[P[j]]%X;
					else {f[i*P[j]]=1LL*f[i]*p[P[j]]%X;break;}
			}
			for(i=2;i<=SZ;++i) f[i]=(f[i]+f[i-1])%X;
		}
};LinearSieve<N> F;
int main()
{
	RI Tt,i,j,l,r,t;scanf("%d%d",&Tt,&k),F.Init();W(Tt--)
	{
		for(scanf("%d%d",&n,&m),t=0,l=1;l<=min(n,m);l=r+1)//除法分块
			r=min(n/(n/l),m/(m/l)),t=(1LL*(n/l)*(m/l)%X*(F[r]-F[l-1]+X)+t)%X;
		printf("%d\n",t);
	}return 0;
}