【BZOJ3601】一个人的数论（莫比乌斯反演+高斯消元）

点此看题面

大致题意： 定义\(f_d(n)\)为所有小于\(n\)且与\(n\)互质的正整数的\(d\)次幂和。现给定\(n=\sum_{i=1}^wp_i^{\alpha_i}\)（\(p_i\)为互不相同的质数），求\(f_d(n)\)。

莫比乌斯反演

显然，根据定义，有：

\[f_d(n)=\sum_{i=1}^ni^d[gcd(i,n)=1] \]

显然\([gcd(i,n)=1]\)可以反演掉，变成：

\[f_d(n)=\sum_{i=1}^ni^d\sum_{k|i,k|n}\mu(k) \]

调整枚举顺序，就得到：

\[f_d(n)=\sum_{k|n}\mu(k)\sum_{i=1}^{\frac nk}(ik)^d \]

从后一项的式子中提出\(k^d\)，就得到：

\[f_d(n)=\sum_{k|n}\mu(k)k^d\sum_{i=1}^{\frac nk}i^d \]

进一步推式子

我推到这一步之后就做不下去了，想了半天都没什么想法，于是就去参考了一下闪指导\(hl666\)的题解，发现其实后面的也不难推，只是我智（yǎn）障（xiā）了。

由于眼瞎，我竟然没有看到数据范围的表格中第一列\(d\le100\)这个条件......

捕捉到这一关键信息之后，赶紧又埋头继续推起了式子。

注意到\(d\)这么小，显然可以利用拉格朗日插值的知识，发现\(\sum_{i=1}^{\frac nk}i^d\)就是一个\(d+1\)项的多项式。

因此，我们设\(\sum_{i=1}^xi^d=\sum_{i=1}^{d+1}a_ix^i\)。

代入到原式之中，得到：

\[f_d(n)=\sum_{k|n}\mu(k)k^d\sum_{i=1}^{d+1}a_i(\frac nk)^i=\sum_{i=1}^{d+1}a_i\sum_{k|n}\mu(k)k^d(\frac nk)^i \]

显然，如果我们提出\(\sum_{i=1}^{d+1}a_i\)这一项，剩下的\(\sum_{k|n}\mu(k)k^d(\frac nk)^i\)是一个积性函数。

看到积性函数，再想到题目中给我们的是\(n\)的质因数表示法，就该顿时明白出题人的用意了吧。

我们设\(g_i(n)=\sum_{k|n}\mu(k)k^d(\frac nk)^i\)，就有：

\[f_d(n)=\sum_{i=1}^{d+1}a_ig_i(n)=\sum_{i=1}^{d+1}a_i\prod_{v=1}^wg_i(p_v^{\alpha_v}) \]

其中，\(g_i(p_v^{\alpha_v})\)的值，如果直接枚举因数，复杂度同样爆炸。

但是，由于式子中有一个\(\mu(k)\)，根据莫比乌斯函数的定义，只要含有平方因子，这个式子值就为\(0\)。

所以，只有\(k=1\)或\(k=p_v\)时原式有值，因此：

\[g_i(p_v^{\alpha_v})=(p_v^{\alpha_v})^i-p_v^d(p_v^{\alpha_v-1})^i=p_v^{\alpha_v\times i}-p_v^{(\alpha_v-1)\times i+d} \]

代回原式就得到：

\[f_d(n)=\sum_{i=1}^{d+1}a_i\prod_{v=1}^w(p_v^{\alpha_v\times i}-p_v^{(\alpha_v-1)\times i+d}) \]

至于\(a_i\)怎么求，数学老师告诉过我们，对于一个这样的多项式，只要知道\(d+1\)个点值，然后用待定系数法列出\(d+1\)个方程，就能解出所有的系数。

因此，我们把\(x=1\sim d+1\)代入，然后高斯消元，就能解出\(a_i\)了。

代码

#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 1000
#define X 1000000007
#define swap(x,y) (x^=y^=x^=y)
#define Qinv(x) Qpow(x,X-2)
#define Inc(x,y) ((x+=(y))>=X&&(x-=X))
using namespace std;
int n,d,p[N+5],q[N+5];
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;}
namespace Gauss//高斯消元
{
	#define SZ 100
	int a[SZ+5][SZ+5],v[SZ+5];
	I void Find(CI x) {RI i=x;W(!a[i][x]) ++i;if(i^x) for(RI j=x;j<=n;++j) swap(a[x][j],a[i][j]);}
	I void Solve(CI n)
	{
		RI i,j,k,t;for(i=1;i<=n;++i) for(Find(i),j=i+1;j<=n;++j)
		{
			t=X-1LL*a[j][i]*Qinv(a[i][i])%X,v[j]=(1LL*t*v[i]+v[j])%X;
			for(k=i;k<=n;++k) a[j][k]=(1LL*t*a[i][k]+a[j][k])%X;
		}
		for(i=n;i;--i)
		{
			v[i]=1LL*v[i]*Qinv(a[i][i])%X;
			for(j=i-1;j;--j) v[j]=(v[j]-1LL*v[i]*a[j][i]%X+X)%X;
		}
	}
}
int main()
{
	using namespace Gauss;
	RI i,j;for(scanf("%d%d",&d,&n),i=1;i<=n;++i) scanf("%d%d",p+i,q+i);
	RI t=0;for(i=1;i<=d+1;++i) for(Inc(t,Qpow(i,d)),v[i]=t,j=1;j<=d+1;++j) a[i][j]=Qpow(i,j);//建立方程
	RI ans=0,res;for(Solve(d+1),i=1;i<=d+1;++i)//解出系数，然后枚举每一个系数去求答案
	{
		for(res=j=1;j<=n;++j) res=1LL*res*//枚举每个质因子，统计函数值
			(Qpow(p[j],1LL*q[j]*i%(X-1))-Qpow(p[j],(1LL*(q[j]-1)*i+d)%(X-1))+X)%X;
		ans=(1LL*v[i]*res+ans)%X;//统计答案
	}return printf("%d",ans),0;
}