如何快速求解第一类斯特林数--nlog^2n + nlogn
- 参考资料
- 前言
- 暴力
- nlog^2n的做法
- nlogn的做法
- 代码
参考资料
百度百科
前言
首先是因为这道题,才去研究了这个玩意:
感觉这个东西非常的...巧妙。
暴力
第一类斯特林树S(n,k)就是将n个数字划分为k个不相区分的圆排列的方案数(即忽略顺序)。
首先,第一类斯特林数有一个人尽皆知的\(O(n^2)\)递推式:
\[S(n,k)=S(n-1,k-1)+(n-1)*S(n-1,k) \]理解起来也是比较容易的。就是考虑新来的一个元素,可以自成一个圆排列,也可以放在前面已经有的圆排列的空位置,而一共有n-1个空位。所以说上式成立。
nlog^2n的做法
这个还是比较好理解的。
类比于二项式定理的形式,其实也有一个关于第一类斯特林数的多项式的形式:
其中上面(1)式的左半部分称作上升幂,(2)式的左半部分称作下降幂。
其实不知道推理也没关系,因为这和二项式定理一样,本身就是一个结论,规定
推导过程:
\[=>f(x,n+1)=(x-n)f(x,n) \]\[=>f(x,n+1)=\sum_{i=0}^{n}S(n,i)(1-)^{n-i}x^{i+1}-n\sum_{i=0}^{n}(-1)^{n-i}S(n,i)x^i \]数学归纳法:
先令\(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):
再合并:
\[=>f(x,n++1)=\sum_{i=1}^{n+1}(S(n,i-1)+n*S(n,i))(-1)^{n-i+1}x^i \]\[=>f(x,n+1)=\sum_{i=1}^{n+1}S(n+1,i)(-1)^{n-i+1}x^i \]对于上升幂来说,也是如此。
为了免去(-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>=1,~j&s;);
if(i=1;i--)
B[i]=(1LL*B[i]*(1LL*deg-1LL)%MO+1LL*B[i-1])%MO;//可以列个式子看一下
}