#定義ios
如今咱們要求:$$S_m(n)=\sum_{k=0}^n k^m$$算法
其中$m>0$。ide
#分析spa
$$\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*}$$code
用$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*}$$io
用遞推式$O(m^2)$或者NTT$O(m\log m)$預處理第二類斯特林數,而後就能夠直接$O(m)$遞推了。class
#推廣stream
要求$S_m(n)$的一個關於$n$的多項式。di
展開爲冪級數:while
$$\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])); }