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

Posted 2020-12-28 Orion545

题面

传送门

思路

这题妙啊

先把式子摆出来

$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}{\\frac{n}{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}{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)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);
}