參考資料
百度百科
斯特林數 學習筆記-by zhouzhendong
前言
首先是因為這道題，才去研究了這個玩意：【2019雅禮集訓】【第一類斯特林數】【NTT&多項式】permutation
感覺這個東西非常的...巧妙。
暴力
第一類斯特林樹S(n,k)就是將n個數字划分為k個不相區分的圓排列的方案數（即忽略順序）。
首先，第一類斯特林數有一個人盡皆知的\(O(n^2)\)遞推式：
理解起來也是比較容易的。就是考慮新來的一個元素，可以自成一個圓排列，也可以放在前面已經有的圓排列的空位置，而一共有n-1個空位。所以說上式成立。
nlog^2n的做法
這個還是比較好理解的。
類比於二項式定理的形式，其實也有一個關於第一類斯特林數的多項式的形式：
其中上面（1）式的左半部分稱作上升冪，（2）式的左半部分稱作下降冪。
其實不知道推理也沒關系，因為這和二項式定理一樣，本身就是一個結論，規定
推導過程：
數學歸納法：
先令\(f(x,n)=x(x-1)(x-2)...(x-n+1)\)
左邊：令i=i-1；
右邊：因為S(n,0)和S(n,n+1)都等於0，所以可以將S(n,0)替換為S(n,n+1):
再合並：
對於上升冪來說，也是如此。
為了免去(-1)這個系數的影響，下面統一使用上升冪進行討論。
觀察上述式子，我們發現，好像可以直接暴力分治FFT：假設當前是計算\((x+L)(x+L+1)...(x+R)\)，然后我們考慮將L ~ mid和mid+1 ~ R分別計算后，再將兩邊的式子乘起來，繼續往上遞回溯就行了。復雜度是\(O(n\log ^2n)\)的。
nlogn的做法
然而我們還可以進一步進行優化。大致的思路是在之前的基礎上，將左半部分，也就是L ~ mid的部分算出來之后，直接使用卷積計算出右半部分(mid+1 ~ R)回溯回來之后的多項式，然后將兩個式子相乘之后直接回溯。如果能夠實現的話，時間復雜度就降到了\(O(nlogn)\)（雖然常數很大）。
假設當前的最高項的次數為2n(先不管奇數的情況)。
現在我們具體來考慮如何通過將L ~ mid進行遞歸計算后返回的多項式(即\(\sum_{i=0}^{n-1}(x+i)\))的系數直接推出遞歸mid+1 ~ R得到的多項式(\(\sum_{i=n}^{2*n-1}(x+i)\))的系數：
左邊的多項式是\(a_0+a_1x+a_2x^2+...+a_nx^n\)，即\(\sum_{i=0}^{n}a_ix^i\)。
那么就可以直接得到右邊一半的式子是：\(\sum_{i=0}^{n}a_i(x+n)^i\)
然后就可以得到:
其中(1)到(2)就是簡單的二項式定理暴力展開。(2)到(3)則是交換了i和j... (3)到(4)是把\(x^i\)單獨提出來了。(3)到(4)則是將inv[j!]單獨提出來，並將j-i和j兩個參量分開，形成了一個標准的減法卷積的形式。可以在最后的時候每一項單獨乘上inv[j!]就可以了。
減法卷積怎么搞呢？
當然是選擇將其中的一個多項式進行反轉，最后再反轉/平移回來啦(巨惡心)！然后我的做法是令\(p_i=\frac{n^{i}}{i!}\),\(q_i=i!a_{i}\)，然后構造多項式\(f_1(x)=\sum_{i=0}^{n}p_ix^i\)和多項式\(f_2(x)=\sum_{i=0}^{n}q_ix^i\)，然后我們將第一個多項式進行翻轉，然后將兩個多項式卷積起來，並設最終式(其中p,q,k的定義都是翻轉前的定義)為：\(\sum_{i=0}^{n}k_ix^i\)就有：
這樣子我們就強行將這個式子轉化為了一個加法卷積式。注意最后只需要將求出來左移n格就是答案了。但是感覺將翻轉后的多項式進行卷積后，答案竟然需要平移回來，有點怪？？？
據說，將\(f_2(x)\)進行翻轉的話，最后就是將卷積得到的數組再翻轉回來了，有興趣的讀者可以自行嘗試。
對了剛剛還沒提次數為奇數的情況怎么處理。假設最高項為2n+1，那么你可以先將其看做是2n次的多項式，按照套路計算完后，再將最后一項(x+2n)暴力O(n)乘上去就可以了。
寫起來總體思路並不是那么惡心，只是細節較多（寫着寫着就卡住了）
代碼
int PowMod(int x,int y)
{
int ret=1;
while(y)
{
if(y&1)
ret=1LL*ret*x%MO;
x=1LL*x*x%MO;
y>>=1;
}
return ret;
}
void Prepare()
{
fact[0]=1;
for(int i=1;i<=MAXN;i++)
fact[i]=1LL*fact[i-1]*i%MO;
inv[MAXN]=PowMod(fact[MAXN],MO-2);//求階乘的逆元
for(int i=MAXN-1;i>=0;i--)
inv[i]=1LL*inv[i+1]*(1LL*i+1LL)%MO;
}
void Reverse(int A[],int deg)
{
for(int i=0;i<deg/2;i++)
swap(A[i],A[deg-i-1]);
}
void NTT(int P[],int len,int oper)
{
for(int i=1,j=0;i<len-1;i++)
{
for(int s=len;j^=s>>=1,~j&s;);
if(i<j) swap(P[i],P[j]);
}
int unit,unit_p0;
for(int d=0;(1<<d)<len;d++)
{
int m=(1<<d),m2=m*2;
unit_p0=PowMod(G,(MO-1)/m2);
if(oper==-1)
unit_p0=PowMod(unit_p0,MO-2);
for(int i=0;i<len;i+=m2)
{
unit=1;
for(int j=0;j<m;j++)
{
int &P1=P[i+j+m],&P2=P[i+j];
int t=1LL*unit*P1%MO;
P1=((1LL*P2-1LL*t)%MO+MO)%MO;
P2=(1LL*P2+1LL*t)%MO;
unit=1LL*unit*unit_p0%MO;
}
}
}
if(oper==-1)
{
int inv=PowMod(len,MO-2);
for(int i=0;i<len;i++)
P[i]=1LL*P[i]*inv%MO;
}
}
void Mul(int ret[],int _x[],int l1,int _y[],int l2)//多項式乘法
{
static int RET[MAXN+5],X[MAXN+5],Y[MAXN+5];
int len=1;
while(len<l1+l2) len<<=1;
copy(_x,_x+l1,X);copy(_y,_y+l2,Y);
fill(X+l1,X+len,0);fill(Y+l2,Y+len,0);
NTT(X,len,1);NTT(Y,len,1);
for(int i=0;i<len;i++)
RET[i]=1LL*X[i]*Y[i]%MO;
NTT(RET,len,-1);
copy(RET,RET+l1+l2,ret);
}
void Get(int deg,int A[],int B[])
{
static int tmpA[MAXN+5],tmpB[MAXN+5];
int len=deg/2;
for(int i=0;i<len+1;i++)
tmpA[i]=1LL*PowMod(len,i)*inv[i]%MO;//先預處理f1(x)
fill(tmpA+len+1,tmpA+deg+1,0);
for(int i=0;i<len+1;i++)
tmpB[i]=1LL*fact[i]*A[i]%MO;//再預處理f2(x)
fill(tmpB+len+1,tmpB+deg+1,0);
Reverse(tmpA,len+1);//將f1進行翻轉
Mul(tmpA,tmpA,len+1,tmpB,len+1);//相乘
for(int i=0;i<=len;i++)
tmpA[i]=1LL*tmpA[i+len]*inv[i]%MO;//最后將結果移位回來，同時乘上inv[j!]
copy(tmpA,tmpA+len+1,B);
}
void Solve(int deg,int B[])//快速求解第一類斯特林數
{//deg為最高次數項的次數（和以前的代碼風格不一樣，不舒服）
static int tmpB[MAXN+5];
if(deg==1)
{
B[1]=1;//終止狀態為x
return;
}
Solve(deg/2,B);//先遞歸求左半部分
int hf=deg/2;
copy(B,B+hf+1,tmpB);//注意+1
fill(tmpB+hf+1,tmpB+deg+1,0);
Get(deg-deg%2,tmpB,tmpB+hf+1);//通過左邊的返回多項式直接求右邊的返回多項式
Mul(B,tmpB,hf+1,tmpB+hf+1,hf+1);//將左右兩邊相乘
if(deg%2==1)//奇數的情況特殊處理
for(int i=deg;i>=1;i--)
B[i]=(1LL*B[i]*(1LL*deg-1LL)%MO+1LL*B[i-1])%MO;//可以列個式子看一下
}