蒟蒻TJY的博客

[第二类斯特林数]自然数幂求和

定义

现在我们要求:$$S_m(n)=\sum_{k=0}^n k^m$$

其中\(m>0\)

分析

\[\begin{align*} S_m(n-1)&=\sum_{k=0}^{n-1}k^m\\ &=\sum_{k=0}^{n-1}\sum_{j=0}^m{m\brace j}k^\underline{j}\tag{将普通幂展开为阶乘幂}\\ &=\sum_{j=0}^m{m\brace j}\sum_{k=0}^{n-1}k^\underline{j}\tag{交换求和次序}\\ &=\sum_{j=0}^m{m\brace j}\sideset{}{_0^n}\sum x^\underline{j}\delta x\tag{写成离散微积分形式}\\ &=\sum_{j=0}^m{m\brace j}\left.\frac{x^\underline{j+1}}{j+1}\right|_0^n\tag{逆差分}\\ &=\sum_{j=0}^m{m\brace j}\frac{n^\underline{j+1}}{j+1}\tag{化简} \end{align*}\]

\(n+1\)代替\(n\)

\[\begin{align*} S_m(n)&=\sum_{j=0}^m{m\brace j}\frac{(n+1)^\underline{j+1}}{j+1}\\ &=(n+1)\sum_{j=0}^m{m\brace j}\frac{n^\underline{j}}{j+1} \end{align*}\]

用递推式\(O(m^2)\)或者NTT\(O(m\log m)\)预处理第二类斯特林数,然后就可以直接\(O(m)\)递推了。

推广

要求\(S_m(n)\)的一个关于\(n\)的多项式。

展开为幂级数:

\[\begin{align*} S_m(n)&=(n+1)\sum_{k=0}^m{m\brace k}\frac{n^{\underline k}}{k+1}\\ &=(n+1)\sum_{k=0}^m{m\brace k}\frac{1}{k+1}\sum_{j=0}^m{k\brack j}(-1)^{k-j}n^j\\ &=(n+1)\sum_{k=0}^m\sum_{j=0}^m{m\brace k}\frac{1}{k+1}{k\brack j}(-1)^{k-j}n^j\\ &=(n+1)\sum_{j=0}^mn^j\sum_{k=0}^m{m\brace k}\frac{1}{k+1}{k\brack j}(-1)^{k-j}\\ \end{align*}\]

则我们令:

\[\begin{align*} T(x)&=\sum_{k=0}^m{m\brace k}\frac{1}{k+1}{k\brack x}(-1)^{k-x}\\ &=(-1)^x\sum_{k=0}^m{m\brace k}{k\brack x}\frac{(-1)^k}{k+1} \end{align*}\]

则有:$$S_m(n)=\sum_{j=1}^{m+1} n^j[T(j)+T(j-1)]$$

这样,我们就得到了一个\(O(m^2)\)的算法 。

代码

给出\(n\)\(m\),求\(S_m(n)\)

所有结果对\(998244353\)取模。

#include<iostream>
#include<cstdio>
using namespace std;
typedef long long ll;
const ll p=998244353;
ll add(ll a,ll b){return a+b>=p?a+b-p:a+b;}
ll cut(ll a,ll b){return a-b<0?a-b+p:a-b;}
ll mul(ll a,ll b){return a*b%p;}
ll pow(ll a,ll b){
	ll ans=1;
	while(b){
		if(b&1)ans=mul(ans,a);
		a=mul(a,a);
		b>>=1;
	}
	return ans;
}
ll div(ll a,ll b){return mul(a,pow(b,p-2));}
int n,m;
ll ans,S2[1001][1001];
int main(){
	scanf("%d%d",&n,&m);
	S2[0][0]=1;
	for(int i=1;i<=m;i++)for(int j=1;j<=m;j++)S2[i][j]=add(S2[i-1][j-1],mul(S2[i-1][j],j));
	ll facpw=n;
	for(int i=1;i<=m;i++){
		ans=add(ans,mul(S2[m][i],div(facpw,i+1)));
		facpw=mul(facpw,n-i);
	}
	ans=mul(ans,n+1);
	printf("%lld\n",ans);
}

给出\(m\),求多项式\(S_m(n)\)

结果为\(x^0,x^1,x^2,\cdots,x^{m+1}\)的系数。

#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
typedef long long ll;
const ll p=998244353;
inline ll add(ll a,ll b){return a+b>=p?a+b-p:a+b;}
inline ll cut(ll a,ll b){return a-b<0?a-b+p:a-b;}
inline ll mul(ll a,ll b){return a*b%p;}
inline ll pow(ll a,ll b){
	ll ans=1;
	while(b){
		if(b&1)ans=mul(ans,a);
		a=mul(a,a);
		b>>=1;
	}
	return ans;
}
inline ll div(ll a,ll b){return mul(a,pow(b,p-2));}
int m;
ll S1[1001][1001],S2[1001][1001],t[1001];
bool fir;
int main(){
    scanf("%d",&m);
    S1[0][0]=S2[0][0]=1;
    for(int i=1;i<=m;i++){
        for(int j=1;j<=m;j++){
            S1[i][j]=add(S1[i-1][j-1],mul(S1[i-1][j],(i-1)));
            S2[i][j]=add(S2[i-1][j-1],mul(S2[i-1][j],j));
        }
    }
    for(int i=1;i<=m;i++){
        for(int j=1;j<=m;j++){
            ll pw=(j%2?p-1:1);
            t[i]=add(t[i],mul(mul(S2[m][j],S1[j][i]),div(pw,j+1)));
        }
        t[i]=mul(t[i],i%2?p-1:1);
    }
    for(int i=0;i<=m+1;i++)printf("%lld ",add(t[i],t[i-1]));
}
posted @ 2018-08-21 15:57  蒟蒻TJY  阅读(1215)  评论(0编辑  收藏  举报