# 题解 P5282 【模板】快速阶乘算法

## 【核心代码】

poly c, d;
inline int get_fac(int n) {
int pos=curStk, s=sqrt(n)+1e-6;
vir invs=Inv[s];
for(int i=s;i>1;i>>=1) Stk[++curStk]=i;//手动压栈

c.resize(2); c[0]=1; c[1]=s+1;//初始化 h_1 的状态
for(int l=Stk[curStk]; curStk>pos; l=Stk[--curStk]) {
ValueTrans(c, d, l>>1, (l>>1)+1);//点值平移出额外状态
c.resize(2*sz(c));
for(int i=0; i<sz(d); ++i) c[sz(d)+i]=d[i];
ValueTrans(c, d, sz(c)-1, invs*vir(l>>1));//点值平移出点乘状态
for(int i=0; i<sz(c); ++i) c[i]=c[i]*d[i];
if(l&1) {
for(int i=0; i<=l; ++i) c[i]=c[i]*vir(i*s+l);//h_l 转为 h_{l+1}
}
else c.resize(l+1);
}

vir res=1;
for(int i=0; i<s; ++i) res=res*c[i];//累乘
for(int i=s*s+1; i<=n; ++i) res=res*vir(i);//处理尾巴
return res;
}


## 【拓展】

\begin{aligned} \begin{pmatrix} (n+1)! \\S_n \end{pmatrix}&= \begin{pmatrix} n+1&0 \\1&1 \end{pmatrix}\cdot \begin{pmatrix} n! \\S_{n-1} \end{pmatrix} \\&=\prod_{i=1}^n \begin{pmatrix} i+1&0 \\1&1 \end{pmatrix}\cdot \begin{pmatrix} 1 \\0 \end{pmatrix} \end{aligned}

$$\begin{pmatrix} a_1&\ \\b_1&c_1 \end{pmatrix}\cdot\begin{pmatrix} a_2&\ \\b_2&c_2 \end{pmatrix}=\begin{pmatrix} a_1a_2&\ \\b_1a_2+c_1b_2&c_1c_2 \end{pmatrix}$$

\begin{aligned} &h_{d+1}(i) \\=&\begin{pmatrix}mc_{d+1,0}(i)&\ \\mc_{d+1,1}(i)&1\end{pmatrix} \\=&\begin{pmatrix}i+d+2&0\\1&1\end{pmatrix}h_d(i) \\=&\begin{pmatrix}i+d+2&0\\1&1\end{pmatrix}\begin{pmatrix}mc_{d,0}(x)&\ \\mc_{d,1}(x)&1\end{pmatrix} \\=&\begin{pmatrix}(i+d+2)mc_{d,0}(x)&\ \\mc_{d,0}+mc_{d,1}(x)&1\end{pmatrix} \end{aligned}

poly mc[2], md[2];
inline int get_sumfac(int n) {
int pos=curStk, s=sqrt(n)+1e-6;
vir invs=Inv[s];
for(int i=s;i>1;i>>=1) Stk[++curStk]=i;

mc[0].resize(2); mc[0][0]=2; mc[0][1]=s+2;
mc[1].resize(2); mc[1][0]=mc[1][1]=1;
for(int l=Stk[curStk]; curStk>pos; l=Stk[--curStk]) {
for(int t=0; t<2; ++t){
poly &c = mc[t], &d = md[t];
ValueTrans(c, d, l>>1, (l>>1)+1);
c.resize(2*sz(c));
for(int i=0; i<sz(d); ++i) c[sz(d)+i]=d[i];
ValueTrans(c, d, sz(c)-1, invs*vir(l>>1));
}
for(int i=0; i<sz(mc[0]); ++i) {
mc[1][i]=mc[0][i]*md[1][i]+mc[1][i];
mc[0][i]=mc[0][i]*md[0][i];
}
if(l&1) {
for(int i=0; i<=l; ++i) {
mc[1][i]=mc[0][i]+mc[1][i];
mc[0][i]=mc[0][i]*vir(i*s+l+1);
}
}
else {
mc[0].resize(l+1);
mc[1].resize(l+1);
}
}

vir res0=1, res1=0;
for(int i=0; i<s; ++i){
res1=res0*mc[1][i]+res1;
res0=res0*mc[0][i];
}
for(int i=s*s+1; i<=n; ++i){
res1=res1+res0;
res0=res0*vir(i+1);
}
return res1;
}


posted @ 2022-01-26 14:59  JustinRochester  阅读(200)  评论(0编辑  收藏  举报