【模板】【洛谷5408】第一类斯特林数·行(倍增)(NTT)

传送门


题解:

首先我们设第一类斯特林数生成函数为 F n ( x ) = ∑ i = 0 n s n , i x i = x n ˉ F_n(x)=\sum_{i=0}^n s_{n,i}x^i=x^{\bar n} Fn(x)=i=0nsn,ixi=xnˉ

则对于奇数的 n n n,可以直接递归处理 F n − 1 F_{n-1} Fn1,然后 O ( n ) O(n) O(n)乘上一个一次多项式。

对于偶数的 n n n,可以处理出 F n / 2 F_{n/2} Fn/2,现在假设已经求出了 F n F_{n} Fn,考虑怎么得到 F 2 n F_{2n} F2n

F n = x n ˉ F 2 n = x 2 n ˉ = x n ˉ ⋅ ( x + n ) n ˉ F_{n}=x^{\bar{n}}\\F_{2n}=x^{\bar{2n}}=x^{\bar{n}}\cdot (x+n)^{\bar{n}} Fn=xnˉF2n=x2nˉ=xnˉ(x+n)nˉ

现在需要支持快速求出 f ( x + k ) f(x+k) f(x+k)的系数表达。

考虑暴力二项式展开:
f ( x + k ) = ∑ i = 0 n a i ( x + k ) i = ∑ i = 0 n a i ∑ j = 0 i ( i j ) x j k i − j = ∑ j = 0 n x j j ! ∑ i = j n a i ⋅ i ! ⋅ k i − j ( i − j ) ! \begin{aligned} f(x+k)&=&&\sum_{i=0}^na_i(x+k)^i\\ &=&&\sum_{i=0}^na_i\sum_{j=0}^i{i\choose j}x^jk^{i-j}\\ &=&&\sum_{j=0}^n \frac{x^j}{j!}\sum_{i=j}^n{a_i}\cdot {i!}\cdot \frac{k^{i-j}}{(i-j)!} \end{aligned} f(x+k)===i=0nai(x+k)ii=0naij=0i(ji)xjkijj=0nj!xji=jnaii!(ij)!kij

发现后面的系数是一个卷积形式,直接做就行了。

复杂度 T ( n ) = T ( n 2 ) + O ( n log ⁡ n ) = O ( n log ⁡ n ) T(n)=T(\dfrac{n}{2})+O(n\log n)=O(n\log n) T(n)=T(2n)+O(nlogn)=O(nlogn)


代码:

#include<bits/stdc++.h>
#define ll long long
#define re register
#define cs const

using std::cerr;
using std::cout;

cs int mod=167772161;
inline int add(int a,int b){return (a+=b)>=mod?a-mod:a;}
inline int dec(int a,int b){return (a-=b)<0?a+mod:a;}
inline int mul(int a,int b){static ll r;r=(ll)a*b;return r>=mod?r%mod:r;}
inline int power(int a,int b,int res=1){
	for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));
	return res;
}
inline void Inc(int &a,int b){(a+=b)>=mod?a-=mod:0;}
inline void Dec(int &a,int b){(a-=b)<0?a+=mod:0;}
inline void Mul(int &a,int b){a=mul(a,b);}

cs int bit=20,SIZE=1<<20|1;

int r[SIZE],*w[bit+1];
int fac[SIZE],ifac[SIZE];
inline void init_NTT(){
	for(int re i=1;i<=bit;++i)w[i]=new int[1<<i-1];
	int wn=power(3,mod-1>>bit);w[bit][0]=1;
	for(int re i=1;i<(1<<bit-1);++i)w[bit][i]=mul(w[bit][i-1],wn);
	for(int re i=bit-1;i;--i)
	for(int re j=0;j<(1<<i-1);++j)w[i][j]=w[i+1][j<<1];
	fac[0]=ifac[0]=fac[1]=1;
	for(int re i=2;i<SIZE;++i)fac[i]=mul(fac[i-1],i);
	ifac[SIZE-1]=power(fac[SIZE-1],mod-2);
	for(int re i=SIZE-2;i;--i)ifac[i]=mul(ifac[i+1],i+1);
}
inline void NTT(int *A,int len,int typ){
	for(int re i=0;i<len;++i)if(i<r[i])std::swap(A[i],A[r[i]]);
	for(int re i=1,d=1;i<len;i<<=1,++d)
	for(int re j=0;j<len;j+=i<<1)
	for(int re k=0;k<i;++k){
		int &t1=A[j+k],&t2=A[j+k+i],t=mul(t2,w[d][k]);
		t2=dec(t1,t),Inc(t1,t);
	}
	if(typ==-1){
		std::reverse(A+1,A+len);
		for(int re i=0,inv=power(len,mod-2);i<len;++i)
		Mul(A[i],inv);
	} 
}
inline void init_rev(int l){
	for(int re i=0;i<l;++i)r[i]=r[i>>1]>>1|((i&1)?l>>1:0);
}

int s[SIZE],b[SIZE],c[SIZE],pw[SIZE];
inline void solve(int len,int *a){
	if(len==1){a[1]=1;return ;}
	if(len&1){
		solve(len-1,a);
		for(int re i=len;i>=1;--i)a[i]=add(a[i-1],mul(a[i],len-1));
		return ;
	}
	solve(len>>1,a);
	int l2=len>>1,l=1;pw[0]=1;
	while(l<=len)l<<=1;init_rev(l);
	for(int re i=1;i<=l2;++i)pw[i]=mul(pw[i-1],l2);
	for(int re i=0;i<=l2;++i)c[i]=mul(a[i],fac[i]);
	for(int re i=0;i<=l2;++i)b[l2-i]=mul(pw[i],ifac[i]);
	memset(c+l2+1,0,sizeof(int)*(l-l2));NTT(c,l,1);
	memset(b+l2+1,0,sizeof(int)*(l-l2));NTT(b,l,1);
	for(int re i=0;i<l;++i)Mul(b[i],c[i]);NTT(b,l,-1);
	for(int re i=0;i<=l2;++i)b[i]=mul(b[i+l2],ifac[i]);
	memset(b+l2+1,0,sizeof(int)*(l-l2));NTT(b,l,1);NTT(a,l,1);
	for(int re i=0;i<l;++i)Mul(a[i],b[i]);NTT(a,l,-1);
}

signed main(){
#ifdef zxyoi
	freopen("first_stirling_row.in","r",stdin);
#endif
	init_NTT();
	int n;scanf("%d",&n);solve(n,s);
	for(int re i=0;i<=n;++i)std::cout<<s[i]<<" ";
	return 0;
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值