定义
现在我们要求：$$S_m(n)=sum_{k=0}^n k^m$$
其中(m>0)。
分析
[egin{align*}
S_m(n-1)&=sum_{k=0}^{n-1}k^m\
&=sum_{k=0}^{n-1}sum_{j=0}^m{mrace j}k^underline{j} ag{将普通幂展开为阶乘幂}\
&=sum_{j=0}^m{mrace j}sum_{k=0}^{n-1}k^underline{j} ag{交换求和次序}\
&=sum_{j=0}^m{mrace j}sideset{}{_0^n}sum x^underline{j}delta x ag{写成离散微积分形式}\
&=sum_{j=0}^m{mrace j}left.frac{x^underline{j+1}}{j+1}
ight|_0^n ag{逆差分}\
&=sum_{j=0}^m{mrace j}frac{n^underline{j+1}}{j+1} ag{化简}
end{align*}]
用(n+1)代替(n)：
[egin{align*}
S_m(n)&=sum_{j=0}^m{mrace j}frac{(n+1)^underline{j+1}}{j+1}\
&=(n+1)sum_{j=0}^m{mrace j}frac{n^underline{j}}{j+1}
end{align*}]
用递推式(O(m^2))或者NTT(O(mlog m))预处理第二类斯特林数，然后就可以直接(O(m))递推了。
推广
要求(S_m(n))的一个关于(n)的多项式。
展开为幂级数：
[egin{align*}
S_m(n)&=(n+1)sum_{k=0}^m{mrace k}frac{n^{underline k}}{k+1}\
&=(n+1)sum_{k=0}^m{mrace k}frac{1}{k+1}sum_{j=0}^m{krack j}(-1)^{k-j}n^j\
&=(n+1)sum_{k=0}^msum_{j=0}^m{mrace k}frac{1}{k+1}{krack j}(-1)^{k-j}n^j\
&=(n+1)sum_{j=0}^mn^jsum_{k=0}^m{mrace k}frac{1}{k+1}{krack j}(-1)^{k-j}\
end{align*}]
则我们令：
[egin{align*}
T(x)&=sum_{k=0}^m{mrace k}frac{1}{k+1}{krack x}(-1)^{k-x}\
&=(-1)^xsum_{k=0}^m{mrace k}{krack 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
",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]));
}