【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^dsum_{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^dsum_{i=1}^{frac nk}i^d ]

进一步推式子

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

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

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

注意到(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^dsum_{i=1}^{d+1}a_i(frac nk)^i=sum_{i=1}^{d+1}a_isum_{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_iprod_{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 imes i}-p_v^{(alpha_v-1) imes i+d} ]

代回原式就得到：

[f_d(n)=sum_{i=1}^{d+1}a_iprod_{v=1}^w(p_v^{alpha_v imes i}-p_v^{(alpha_v-1) imes i+d}) ]

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

因此，我们把(x=1sim 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;
}