[bzoj3601] 一个人的数论 [莫比乌斯反演+高斯消元]

题面

传送门

思路

这题妙啊

先把式子摆出来

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

这个$gcd$看着碍眼，我们把它反演掉

$f_n(d)=\sum_{i=1}^{n\sum_{j|i,j|n}\mu(j)i}d=\sum_{j|n}\mu(j)\sum_{i=1}^{{\frac{n}{j}}(ij)}d=\sum_{j|n}\mu(j)j^d\sum_{i=1}{j}}i^d$

那么最后面这个东西就是个自然数幂求和了

在这篇关于斯特林数的blog最后，我给出了自然数幂求和转斯特林数的公式：

$x^n=\sum_{i=1}n \begin{Bmatrix} n\\i \end{Bmatrix} \frac{x!}{(x-i)!}$

我们对左边的$x$，取$1...m$求和，得到$\sum_{i=1}^m i^n=\sum_{j=1}m \sum_{i=1}^j \begin{Bmatrix} j\\i \end{Bmatrix} \frac{j!}{(j-i)!}$

由此可得，右边这个东西显然是一个关于$i$（也就是原来那个式子里面的$x$）的，在$1...n+1$项上有系数的多项式

（其实还有另外一个公式：$Sum_k(n)=\sum_{i=1}^n i^k=\sum_{j=1}k\begin{Bmatrix}k\\j\end{Bmatrix}\frac{{(n+1)}^\underline{j+1}}{j+1}$）

（好像这个简单易懂一点= =）

不管怎么样，我们可以设$\sum_{i=1}^x i^d =\sum_{i=1}^{d+1}c_i x^i$

然后我们对于$x=1...d+1$分别求出$c_i$那一项的系数，我们实际上得到了一个$d+1$元线性方程组

可以高斯消元之，得到$c$数组

再把$c$放进式子里面，得到：

$f_n(d)=\sum_{j|n}\mu(j)j^d\sum_{i=1}c_i(\frac{n}{j})^i=\sum_{i=1} c_i \sum_{j|n}\mu(j)j^d (\frac{n}{j})^i$

显然后面那个$\sum$里面的一坨东西是个积性函数（因为是两个积性函数的狄利克雷卷积）对吧

我们设$H(i)=\sum_{j|i}\mu(j)j^d (\frac{n}{j})^{i$，那么$H(n)=\prod_{i=1}}w H(p_i^{a_i})$

代回原式：

$f_n(d)=\sum_{i=1}^{d+1} c_i \prod_{j=1}^w H(p_j{a_j})=\sum_{i=1}^{d+1} c_i \prod_{j=1}^w \sum_{k|p_j^{{a_j}}\mu(k)k}d(\frac{p_j^{a_j}}{k})i$

后面这个式子，显然当且仅当$k=1$和$k=p_j$的时候有值（因为其他时候$\mu(k)=0$），那么把这个两项代入，可以得到：

$f_n(d)=\sum_{i+1}^{d+1} c_i \prod_{j=1}^w (p_j^{a_j\ast i}-p_j^{d+a_j\ast i -i})$

那么就解决了，总复杂度是$O(d^3+dw)$

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cassert>
#include<cmath>
#define ll long long
#define MOD 1000000007
using namespace std;
inline ll read(){
	ll re=0,flag=1;char ch=getchar();
	while(ch>'9'||ch<'0'){
		if(ch=='-') flag=-1;
		ch=getchar();
	}
	while(ch>='0'&&ch<='9') re=(re<<1)+(re<<3)+ch-'0',ch=getchar();
	return re*flag;
}
ll d,w,p[1010],a[1010],c[110][110],x[110];
ll qpow(ll a,ll b){
	ll re=1;
	while(b){
		if(b&1) re=re*a%MOD;
		a=a*a%MOD;b>>=1;
	}
	return re;
}
void Gauss(){
	ll i,j,k,num;ll tmp,tt;
	for(i=1;i<=d+1;i++){
		num=i;
		for(j=i+1;j<=d+1;j++) if(abs(c[j][i])>abs(c[num][i])) num=j;
		for(j=1;j<=d+2;j++) swap(c[i][j],c[num][j]);
		tmp=qpow(c[i][i],MOD-2);
		for(j=i+1;j<=d+1;j++){
			tt=c[j][i]*tmp%MOD;
			for(k=1;k<=d+2;k++) c[j][k]=(c[j][k]-tt*c[i][k]%MOD+MOD)%MOD;
		}
	}
	for(i=d+1;i>=1;i--){
		x[i]=c[i][d+2]=c[i][d+2]*qpow(c[i][i],MOD-2)%MOD;
		for(j=i-1;j>=1;j--) c[j][d+2]=(c[j][d+2]-c[j][i]*c[i][d+2]%MOD+MOD)%MOD;
	}
}
int main(){
	d=read();w=read();ll i,j,tmp,sum=0;
	for(i=1;i<=w;i++) p[i]=read(),a[i]=read();
	for(i=1;i<=d+1;i++){
		tmp=1;sum+=qpow(i,d);sum%=MOD;
		for(j=1;j<=d+1;j++){
			tmp=tmp*i%MOD;
			c[i][j]=tmp;
		}
		c[i][d+2]=sum;
	}
	Gauss();tmp=0;
	for(i=1;i<=d+1;i++){
		sum=1;
		for(j=1;j<=w;j++) sum*=(qpow(p[j],a[j]*i)-qpow(p[j],d+a[j]*i-i)+MOD),sum%=MOD;
		tmp=(tmp+x[i]*sum)%MOD;
	}
	printf("%lld\n",tmp);
}