如何快速求解第一类斯特林数--nlog^2n + nlogn
目录
- 参考资料
- 前言
- 暴力
- nlog^2n的做法
- nlogn的做法
- 代码
参考资料
百度百科
斯特林数 学习笔记-by zhouzhendong
前言
首先是因为这道题,才去研究了这个玩意:【2019雅礼集训】【第一类斯特林数】【NTT&多项式】permutation
感觉这个东西非常的...巧妙。
暴力
第一类斯特林树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的做法
这个还是比较好理解的。
类比于二项式定理的形式,其实也有一个关于第一类斯特林数的多项式的形式:
\[x(x+1)(x+2)(x+3)...(x+n-1)=\sum_{i=0}^{n}S(n,i)x^i---(1)\]
\[x(x-1)(x-2)(x-3)...x(x-n+1)=\sum_{i=0}^{n}(-1)^{n-i}S(n,i)x^i---(2)\]
其中上面(1)式的左半部分称作上升幂,(2)式的左半部分称作下降幂。
其实不知道推理也没关系,因为这和二项式定理一样,本身就是一个结论,规定
推导过程:
数学归纳法:
先令\(f(x,n)=x(x-1)(x-2)...(x-n+1)\)
\[=>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\]
左边:令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)(-1)^{n-i+1}x^i-\sum_{i=1}^{n+1}S(n,i-1)(-1)^{n-i}x^i\]
再合并:
\[=>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\)
然后就可以得到:
\[ \begin{align} \sum_{i=1}^{n}a_i(x+n)^i &=\sum_{i=1}^{n}a_i\sum_{j=0}^{i}C_{i}^{j}x^jn^{i-j} \tag{1}\\ &=\sum_{j=0}^{n}x^{j}\sum_{i=j}^{n}C_{i}^{j}n^{i-j}a_{i} \tag{2}\\ &=\sum_{i=0}^{n}x^{i}\sum_{j=i}^{n}C_{j}^{i}n^{j-i}a_{j} \tag{3}\\ &=\sum_{i=0}^{n}x^{i}\sum_{j=i}^{n}\frac{1}{i!}*(\frac{n^{j-i}}{(j-i)!})*(j!a_j) \tag{4} \end{align} \]
其中(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\)就有:
\[ \begin{align} k_{i-j}&=\sum_{i,j} p_{j}q_i\\ k_{i+j}&=\sum_{i,j} p_{-j}q_i\\ k_{i+j-n}&=\sum_{i,j} p_{n-j}q_i\\ k'_{i+j}&=\sum_{i,j} p'_{j}q_i\\ \end{align} \]
这样子我们就强行将这个式子转化为了一个加法卷积式。注意最后只需要将求出来左移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;//终止状态为xreturn;}Solve(deg/2,B);//先递归求左半部分int hf=deg/2;copy(B,B+hf+1,tmpB);//注意+1fill(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;//可以列个式子看一下
}
转载于:https://www.cnblogs.com/T-Y-P-E/p/10274219.html
如何快速求解第一类斯特林数--nlog^2n + nlogn相关推荐
- 【2019雅礼集训】【CF 960G】【第一类斯特林数】【NTT多项式】permutation
目录 题意 输入格式 输出格式 思路 代码 题意 找有多少个长度为n的排列,使得从左往右数,有a个元素比之前的所有数字都大,从右往左数,有b个元素比之后的所有数字都大. n<=2*10^5,a, ...
- 第一类斯特林数 / 第二类斯特林数 / 贝尔数 小结
第一类斯特林数 有 nnn 个不同的小球,将它们串成 mmm 条项链,有多少种不同的方案? 第一类斯特林数的表示方法为 [nm]\left[\begin{matrix}n\\m\end{matrix} ...
- 第一类斯特林数学习记录
最近做题有时会碰到斯特林数(Stirling数),就觉得好好的学习一番,于是呢,写下这篇博客,来记录一些知识 简单介绍 第一类斯特林数表示表示将 n 个不同元素构成m个圆排列的数目.--百度百科 第一 ...
- CF960G-Bandit Blues【第一类斯特林数,分治,NTT】
正题 题目链接:https://www.luogu.com.cn/problem/CF960G 题目大意 求有多少个长度为nnn的排列,使得有AAA个前缀最大值和BBB个后缀最大值. 0≤n,A,B≤ ...
- 组合基础2 第一类斯特林数 第二类斯特林数 基础部分
记xn‾=x(x+1)(x+2)⋯(x+n−1)x^{\overline{n}}=x(x+1)(x+2)\cdots(x+n-1)xn=x(x+1)(x+2)⋯(x+n−1),xn‾=x(x−1)(x ...
- 洛谷P4609 [FJOI2016]建筑师 【第一类斯特林数】
题目链接 洛谷P4609 题解 感性理解一下: 一神带\(n\)坑 所以我们只需将除了\(n\)外的\(n - 1\)个元素分成\(A + B - 2\)个集合,每个集合选出最大的在一端,剩余进行排列 ...
- 建筑师 第一类斯特林数
文章目录 目录 题意: 思路: 目录 题意: 给你一个nnn的排列,排列中的数代表他的高度,问你有多少个排列能使得从左边能看到aaa个建筑,从右边能看到bbb个建筑. 如果建筑iii左边没有任何比他高 ...
- [HDU 3625] Examining the Rooms(第一类斯特林数)
Examining the Rooms problem solution code problem hdu 3625 solution 之前考试有一道题:最多砸开 kkk 扇门,采取最有操作,求把 n ...
- [数学最安逸][UVa1638改编][第一类斯特林数+组合数]杆子的排列
有高为1,2,3,...,n的杆子各一根排成一行.从左边能看到l根,从右边能看到r根,求有多少种可能. (l,r <= 200,n <= 200000) 给出T 组数据 (T <= ...
- 【BJOI2019】勘破神机(下降幂转自然幂)(第一类斯特林数)(特征方程)
传送门 题解: 完全自己推出来的第一道数学神题. 首先我们知道宽度为222的部分方案数是斐波那契数列. 设fnf_nfn表示长度为nnn的时候方案数,题目要求的实际上是这个东西: ∑n=lr(fnk ...
最新文章
- C语言比较好的风格梳理
- HBase scan setBatch和setCaching的区别
- 高德推出查岗功能_新型「查岗」工具?高德推出「家人地图」新功能
- 使用Instant Client配置PL/SQL Developer
- 使用Arquillian和LocalStack脱机测试AWS云堆栈
- Hibernate第十一篇【配置C3P0数据库连接池、线程Session】
- 在ASP.NET使用javascript的一点小技巧(转www.chinacs.net 中文C#技术站 )
- 数组的内存理解(图示)
- 阿里巴巴发布招聘微博:新财年新增超过1800岗位需求
- linux ps命令查看当前线程正在执行的程序
- z变换与s变换之间的转换(一些零碎且不严谨的想法)
- dropbox与public
- zb服务器连接不稳定,绝对惊人!全球服务器处理9.57ZB数据
- win10双显示屏,分屏显示内容
- ff14拆区后哪个服务器人最多,《最终幻想14》官宣拆区!国服大区调整计划公布...
- 淘管家一键铺货怎么弄?和分销下单有什么区别?
- bp神经网络推导以及物理意义
- 盛世昊通董车长2.0,数字化整合行业产业链变革
- html笔记(完整版)
- oracle10g SGA
热门文章
- Aibaba Dubbo 的前世今生以及黑历史 主程序员梁飞 阿里P9(2016年查看)
- 极客大学架构师训练营 微服务架构 领域驱动设计DDD 中台架构、组件设计原则 第十次作业
- 算法:Reverse String(反转字符串)
- linux挂载曙光存储,曙光I1620G30获取设备的cpu、内存、存储等参数信息。
- numpy.squeeze()的用法
- mysql join explain_详解 MySQL 中的 explain
- 计算机班英语试卷考法,计算机专业英语期末考试试卷A
- pwm 正弦波_增强型PWM抑制功能对于直列式电机控制的五大优势
- 知识点收录04:MAVEN相关的知识点
- Java Foundation serial ( 一 )