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

时间：2021-08-09 00:49:21

题面

传送门

思路

这题妙啊

先把式子摆出来

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

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

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

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

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

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

我们对左边的$x$，取$1...m$求和，得到$\sum_{i=1}^m in=\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 ik=\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)jd\sum_{i=1}{d+1}c_i(\frac{n}{j})i=\sum_{i=1}{d+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)kd(\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);

}

相关文章