P5409 第一类斯特林数·列
第一类斯特林数\(\begin{bmatrix}n\\ m\end{bmatrix}\)表示将\(n\)个不同元素构成\(m\)个圆排列的数目。
给定\(n,k\)，对于所有的整数\(i\in[0,n]\)，你要求出\(\begin{bmatrix}i\\ k\end{bmatrix}\)。
由于答案会非常大，所以你的输出需要对\(167772161\)（\(2^{25}\times 5+1\)，是一个质数）取模。
\(1\le k,n< 131072\)
\(S(n,k)=S(n-1,k-1)+(n-1)S(n-1,k)\)
考虑利用指数型生成函数来做
当\(k=1\)时，有\(F(x)=\sum_i(i-1)!\frac{x^i}{i!}\)
则答案为\(\frac{F(x)^k}{k!}\)
利用快速幂可以实现\(nlog^2\)
更快的采用多项式快速幂\(ln+exp\)实现\(nlogn\)
code
//#include<bits/stdc++.h>
#include<iostream>
#include<cstdio>
#include<ctime>
#include<cctype>
#include<queue>
#include<deque>
#include<stack>
#include<iostream>
#include<iomanip>
#include<cstdio>
#include<cstring>
#include<string>
#include<ctime>
#include<cmath>
#include<cctype>
#include<cstdlib>
#include<queue>
#include<deque>
#include<stack>
#include<vector>
#include<algorithm>
#include<utility>
#include<bitset>
#include<set>
#include<map>
#define ll long long
#define db double
#define INF 1000000000
#define inf 100000000000000000ll
#define ldb long double
#define pb push_back
#define put_(x) printf("%d ",x)
#define putl_(x) printf("%lld ",x)
#define get(x) x=read()
#define putl(x) printf("%lld\n",x)
#define rep(p,n,i) for(int i=p;i<=n;i+=1)
#define fep(n,p,i) for(int i=n;i>=p;--i)
#define go(x) for(int i=lin[x],tn=ver[i];i;tn=ver[i=nex[i]])
#define pii pair<int,int>
#define mk make_pair
#define gf(x) scanf("%lf",&x)
#define pf(x) ((x)*(x))
#define uint unsigned long long
#define ui unsigned
#define sq sqrt
#define x(w) t[w].x
#define r(w) t[w].r
#define id(w) t[w].id
#define R(w) s[w].r
#define yy p<<1|1
#define zz p<<1
#define sum(w) t[w].sum
#define mod 167772161
#define sc(A) scanf("%d",&A)
#define scs(A) scanf("%s",A);
#define put(A) printf("%d\n",A)
#define min(x,y) (x>=y?y:x)
#define max(x,y) (x>=y?x:y)
#define sub(x,y) (x-y<0?x-y+mod:x-y)
using namespace std;
const int MAXN=1<<19|1,G=3;
int n,k,m,T,cnt,top,lim;
int f[MAXN],g[MAXN],h[MAXN],a[MAXN],b[MAXN],F[MAXN],w[MAXN],w1[MAXN];
int fac[MAXN],inv[MAXN],in[MAXN],rev[MAXN];
inline int ksm(int b,int p)
{
int cnt=1;
while(p)
{
if(p&1)cnt=(ll)cnt*b%mod;
b=(ll)b*b%mod;
p=p>>1;
}
return cnt;
}
inline void NTT(int *a,int op)
{
rep(0,lim-1,i)if(i<rev[i])swap(a[i],a[rev[i]]);
for(int len=2;len<=lim;len=len<<1)
{
int mid=len>>1;
int wn=ksm(G,op==1?(mod-1)/len:mod-1-(mod-1)/len);
for(int j=0;j<lim;j+=len)
{
int d=1;
rep(0,mid-1,i)
{
int x=a[i+j];
int y=(ll)a[i+j+mid]*d%mod;
a[i+j]=(x+y)%mod;
a[i+j+mid]=(x-y+mod)%mod;
d=(ll)d*wn%mod;
}
}
}
if(op==-1)
{
int IN=ksm(lim,mod-2);
rep(0,lim-1,i)a[i]=(ll)a[i]*IN%mod;
}
}
inline void qn(int *g,int *f,int n)//求逆
//g会更改 f不会。
{
if(n==1)
{
g[0]=1;//特殊了
return;
}
qn(g,f,(n+1)>>1);
rep(0,n-1,i)a[i]=g[i],b[i]=f[i];
lim=1;
while(lim<=n-1+n-1)lim=lim<<1;
rep(0,lim-1,i)rev[i]=rev[i>>1]>>1|((i&1)?lim>>1:0);
NTT(a,1);NTT(b,1);
rep(0,lim-1,i)a[i]=(ll)a[i]*b[i]%mod*a[i]%mod;
NTT(a,-1);
rep(0,lim-1,i)
{
if(i<=n-1)g[i]=((ll)2*g[i]-a[i]+mod)%mod;
a[i]=b[i]=0;
}
}
inline void qd(int *g,int *f,int n)
{
rep(1,n,i)g[i-1]=(ll)f[i]*i%mod;
}
inline void jf(int *g,int n)
{
fep(n,1,i)g[i]=(ll)g[i-1]*in[i]%mod;
g[0]=0;
}
inline void ql(int *g,int *f,int n)
{
qn(h,f,n+1);
qd(F,f,n);
lim=1;
while(lim<=n+n-1)lim=lim<<1;
rep(0,lim-1,i)rev[i]=rev[i>>1]>>1|((i&1)?lim>>1:0);
NTT(h,1);NTT(F,1);
rep(0,lim-1,i)h[i]=(ll)h[i]*F[i]%mod;
NTT(h,-1);
rep(0,lim-1,i)
{
if(i<=n)g[i]=h[i];
h[i]=F[i]=0;
}
jf(g,n);
}
inline void exp(int *w,int *g,int n)
{
if(n==1)
{
w[0]=1;
return;
}
exp(w,g,(n+1)>>1);
ql(w1,w,n-1);
rep(0,n-1,i)w1[i]=(mod-w1[i]+g[i])%mod;
++w1[0];
lim=1;
while(lim<=n-1+n-1)lim=lim<<1;
rep(0,lim-1,i)rev[i]=rev[i>>1]>>1|((i&1)?lim>>1:0);
NTT(w,1);NTT(w1,1);
rep(0,lim-1,i)w[i]=(ll)w[i]*w1[i]%mod;
NTT(w,-1);
rep(0,lim-1,i)
{
w1[i]=0;
if(i>=n)w[i]=0;
}
}
signed main()
{
//freopen("1.in","r",stdin);
sc(n);sc(k);
fac[0]=1;
rep(1,n,i)fac[i]=(ll)fac[i-1]*i%mod;
inv[n]=ksm(fac[n],mod-2);
for(int i=n-1;i>=0;--i)inv[i]=(ll)inv[i+1]*(i+1)%mod,in[i+1]=(ll)fac[i]*inv[i+1]%mod;
rep(1,n,i)f[i-1]=in[i];//求f^k
//g(x)^k=f(x)^kx^k
ql(g,f,n-1);
rep(0,n,i)g[i]=(ll)g[i]*k%mod;
exp(w,g,n+1);
fep(n,0,i)
{
if(i>=k)w[i]=w[i-k];
else w[i]=0;
}
rep(0,n,i)put_((ll)w[i]*fac[i]%mod*inv[k]%mod);
return 0;
}
代码有空再写 需要很复杂的东西。
