题解 Feux Follets

问 joke3579 杂题选讲有什么推荐的，主题库里他推了这个东西。一开始不知道是啥，点进去一看 Feux Follets。惊吓。

P7438

首先第一步把我卡死。先把这个 \(F\) 写成牛顿级数

\[F(x)=\sum_{i=0}^ka_i\binom xi \]

这个可以直接牛顿插值 \(O(k^2)\)。剩下的就平凡了。

\[\sum_{\pi}\sum_{i=0}^ka_i\binom {\text{cyc}_{\pi}}i=\sum_{i=0}^ka_i\sum_{\pi}\binom{\text{cyc}_{\pi}}i \]

算后边这个东西。首先我们知道错排的 EGF 是

\[e^{-\ln(1-x)-x} \]

推一下这玩意怎么来的。首先是环排列的 EGF

\[\sum_{i=1}\frac{(i-1)!x^i}{i!}=\sum_{i=1}\frac{x^i}i=-\ln(1-x) \]

然后不能有自环，减个 \(x\)。然后把所有环团起来，就是 \(\exp\) 一下。

然后我们是从所有的错排的环里边选 \(i\) 个，所以环要乘一个 \(1+y\) 来标记选没选。于是就变成了：

\[f(x,y)=e^{(-\ln(1-x)-x)(1+y)}=\left(\frac 1{(1-x)e^x}\right)^{1+y} \]

我们就要求所有 \(i\le n,j<k\) 的 \(i![x^iy^j]f(x,y)\)。

对于这种小结构，直接求导并表示自己通常是极为有效的（有时我们管它叫 ODE 手动机）。于是求导并表示自己：

\[\begin{aligned} &\frac{\partial f}{\partial x}\\ =&\frac{xe^x(1+y)}{((1-x)e^x)^2}\left(\frac 1{(1-x)e^x}\right)^y\\ =&\frac{x(1+y)}{1-x}f \end{aligned} \]

提取系数可以得到：

\[nf_{n,k}=\sum_{i=0}^{n-2}f_{i,k}+f_{i,k-1}=(n-1)f_{n-1,k}+f_{n-2,k}+f_{n-2,k-1} \]

边界条件 \(f_{2,0}=f_{2,1}=\dfrac 1{2!}\)。递推即可，复杂度 \(O(nk)\)。

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
using namespace std;
const int mod=998244353;
int n,k,f[600010][110],jc[6000010],inv[6000010],a[110],g[110][110];
int qpow(int a,int b){
    int ans=1;
    while(b){
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;
        b>>=1;
    }
    return ans;
}
int F(int x){
    int ret=1,ans=0;
    for(int i=0;i<k;i++)ans=(ans+1ll*ret*a[i])%mod,ret=1ll*ret*x%mod;
    return ans;
}
int main(){
    scanf("%d%d",&n,&k);jc[0]=inv[0]=1;
    for(int i=1;i<=n;i++)jc[i]=1ll*jc[i-1]*i%mod;
    inv[n]=qpow(jc[n],mod-2);
    for(int i=n-1;i>=1;i--)inv[i]=1ll*inv[i+1]*(i+1)%mod,inv[i+1]=1ll*inv[i+1]*jc[i]%mod;
    for(int i=0;i<k;i++)scanf("%d",&a[i]);
    for(int i=0;i<k;i++)g[0][i]=F(i);
    for(int i=1;i<k;i++){
        for(int j=0;j<k-i;j++){
            g[i][j]=(g[i-1][j+1]-g[i-1][j]+mod)%mod;
        }
    }
    for(int i=0;i<k;i++)a[i]=g[i][0];
    f[2][0]=f[2][1]=inv[2];
    for(int i=3;i<=n;i++){
        f[i][0]=(1ll*(i-1)*f[i-1][0]%mod+f[i-2][0])%mod*inv[i]%mod;
        for(int j=1;j<k;j++)f[i][j]=(1ll*(i-1)*f[i-1][j]%mod+f[i-2][j]+f[i-2][j-1])%mod*inv[i]%mod;
    }
    for(int i=1;i<=n;i++){
        int ans=0;
        for(int j=0;j<k;j++)ans=(ans+1ll*f[i][j]*a[j])%mod;
        ans=1ll*ans*jc[i]%mod;
        printf("%d ",ans);
    }
    return 0;
}

P7439

仍然考虑上边的思路。然而我们推一下可以发现其实不用什么牛顿级数。直接枚举 \(\text{cyc}_{\pi}\)，然后每个地方乘上多点求值得到的系数即可。上边那个不能这么搞是因为 \(n\) 太大了。

那么我们只需要考虑 \(\text{cyc}_{\pi}\) 。我们考虑它的生成函数

\[F(y)=[x^n]e^{(-\ln(1-x)-x)y} \]

后边的东西相对比较套路。看中间的一团 \(-\ln(1-x)-x\)，它露着一个 \(x\) 看起来就很像拉反。然而这玩意最低次项是 \(x^2\) 且系数是 \(\frac 12\)。那么设

\[\frac{f^2}2=-\ln(1-x)-x \]

显然 \(f\) 有复合逆，设为 \(g\)。那么

\[[x^n]e^{y\frac{f^2}2}=\frac 1n[x^{n-1}]\frac{\partial e^{y\frac{x^2}2}}{\partial x}\left(\frac xg\right)^n \]

观察中间

\[\begin{aligned} &\frac{\partial e^{y\frac{x^2}2}}{\partial x}\\ =&xye^{y\frac{x^2}2}\\ =&\sum_{i=0}\frac{x^{2i+1}y^{i+1}}{2^ii!} \end{aligned} \]

于是

\[\begin{aligned} &[x^ny^k]e^{y\frac{f^2}2}\\ =&\frac 1n[x^{n-1}y^k]\frac{\partial e^{y\frac{x^2}2}}{\partial x}\left(\frac xg\right)^n\\ =&\frac 1{n2^{k-1}(k-1)!}[x^{n-2k}]\left(\frac xg\right)^n \end{aligned} \]

看看怎么算 \(g\)。观察定义式

\[\frac{f^2}2=-\ln(1-x)-x \]

将 \(x\) 用 \(g\) 替换：

\[\frac{x^2}2=-\ln(1-g)-g \]

没写，因为懒。

P7440

我们要对所有的 \(n\) 计算

\[\sum_{i=1}^nF(i)[x^ky^i]e^{y(-\ln(1-x)-x)} \]

那设 \(f(x,y)=e^{y(-\ln(1-x)-x)}\)。看到这种东西感觉无从下手，试试转置原理。

一开始的转移矩阵显然是 \(a_{i,j}=[x^iy^j]f(x,y)\)。转置之后答案的生成函数为：

\[G=\sum_{i=1}^nF(i)[x^i]f(x,y) \]

考虑后边的 \([x^i]f(x,y)\)。设 \(f_i(y)=[x^i]f(x,y)\)，则求导可以得到递推式

\[f_i=\frac{i-1}if_{i-1}+\frac yif_{i-2} \]

那么这是个矩阵递推。转移矩阵为

\[A_i=\begin{bmatrix}\frac{i-1}i & \frac yi\\1 & 0\end{bmatrix} \]

那么答案的生成函数可以变成

\[G=\sum_{i=1}^nF(i)\prod_{j=1}^iA_j \]

把左上角拿出来。

这个实际上可以分治 FFT。设 \(P_{l,r}=\prod_{i=l}^rA_i,Q_{l,r}=\sum_{i=l}^rF(i)\prod_{j=l}^iA_j\)。转移显然：

\[P_{l,r}=P_{l,mid}P_{mid+1,r} \]

\[Q_{l,r}=Q_{l,mid}+Q_{mid+1,r}P_{l,mid} \]

然后将上述算法转置即可。转置后的矩阵乘就是把矩阵转置一下然后把输出当输入，乘法变转置乘法。

cnm，为什么卡常。

cnm，卡过了。

问 joke3579 杂题选讲有什么推荐的，主题库里他推了这个东西。一开始不知道是啥，点进去一看 Feux Follets。惊吓。

## P7438

首先第一步把我卡死。先把这个 $F$ 写成牛顿级数
$$F(x)=\sum_{i=0}^ka_i\binom xi$$
这个可以直接牛顿插值 $O(k^2)$。剩下的就平凡了。
$$\sum_{\pi}\sum_{i=0}^ka_i\binom {\text{cyc}_{\pi}}i=\sum_{i=0}^ka_i\sum_{\pi}\binom{\text{cyc}_{\pi}}i$$
算后边这个东西。首先我们知道错排的 EGF 是
$$e^{-\ln(1-x)-x}$$
推一下这玩意怎么来的。首先是环排列的 EGF
$$\sum_{i=1}\frac{(i-1)!x^i}{i!}=\sum_{i=1}\frac{x^i}i=-\ln(1-x)$$
然后不能有自环，减个 $x$。然后把所有环团起来，就是 $\exp$ 一下。

然后我们是从所有的错排的环里边选 $i$ 个，所以环要乘一个 $1+y$ 来标记选没选。于是就变成了：
$$f(x,y)=e^{(-\ln(1-x)-x)(1+y)}=\left(\frac 1{(1-x)e^x}\right)^{1+y}$$
我们就要求所有 $i\le n,j<k$ 的 $i![x^iy^j]f(x,y)$。

对于这种小结构，直接求导并表示自己通常是极为有效的（有时我们管它叫 ODE 手动机）。于是求导并表示自己：
$$
\begin{aligned}
&\frac{\partial f}{\partial x}\\
=&\frac{xe^x(1+y)}{((1-x)e^x)^2}\left(\frac 1{(1-x)e^x}\right)^y\\
=&\frac{x(1+y)}{1-x}f
\end{aligned}
$$
提取系数可以得到：
$$nf_{n,k}=\sum_{i=0}^{n-2}f_{i,k}+f_{i,k-1}=(n-1)f_{n-1,k}+f_{n-2,k}+f_{n-2,k-1}$$
边界条件 $f_{2,0}=f_{2,1}=\dfrac 1{2!}$。递推即可，复杂度 $O(nk)$。
```cpp
#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
using namespace std;
const int mod=998244353;
int n,k,f[600010][110],jc[6000010],inv[6000010],a[110],g[110][110];
int qpow(int a,int b){
    int ans=1;
    while(b){
        if(b&1)ans=1ll*ans*a%mod;
        a=1ll*a*a%mod;
        b>>=1;
    }
    return ans;
}
int F(int x){
    int ret=1,ans=0;
    for(int i=0;i<k;i++)ans=(ans+1ll*ret*a[i])%mod,ret=1ll*ret*x%mod;
    return ans;
}
int main(){
    scanf("%d%d",&n,&k);jc[0]=inv[0]=1;
    for(int i=1;i<=n;i++)jc[i]=1ll*jc[i-1]*i%mod;
    inv[n]=qpow(jc[n],mod-2);
    for(int i=n-1;i>=1;i--)inv[i]=1ll*inv[i+1]*(i+1)%mod,inv[i+1]=1ll*inv[i+1]*jc[i]%mod;
    for(int i=0;i<k;i++)scanf("%d",&a[i]);
    for(int i=0;i<k;i++)g[0][i]=F(i);
    for(int i=1;i<k;i++){
        for(int j=0;j<k-i;j++){
            g[i][j]=(g[i-1][j+1]-g[i-1][j]+mod)%mod;
        }
    }
    for(int i=0;i<k;i++)a[i]=g[i][0];
    f[2][0]=f[2][1]=inv[2];
    for(int i=3;i<=n;i++){
        f[i][0]=(1ll*(i-1)*f[i-1][0]%mod+f[i-2][0])%mod*inv[i]%mod;
        for(int j=1;j<k;j++)f[i][j]=(1ll*(i-1)*f[i-1][j]%mod+f[i-2][j]+f[i-2][j-1])%mod*inv[i]%mod;
    }
    for(int i=1;i<=n;i++){
        int ans=0;
        for(int j=0;j<k;j++)ans=(ans+1ll*f[i][j]*a[j])%mod;
        ans=1ll*ans*jc[i]%mod;
        printf("%d ",ans);
    }
    return 0;
}

P7439

仍然考虑上边的思路。然而我们推一下可以发现其实不用什么牛顿级数。直接枚举 \(\text{cyc}_{\pi}\)，然后每个地方乘上多点求值得到的系数即可。上边那个不能这么搞是因为 \(n\) 太大了。

那么我们只需要考虑 \(\text{cyc}_{\pi}\) 。我们考虑它的生成函数

\[F(y)=[x^n]e^{(-\ln(1-x)-x)y} \]

后边的东西相对比较套路。看中间的一团 \(-\ln(1-x)-x\)，它露着一个 \(x\) 看起来就很像拉反。然而这玩意最低次项是 \(x^2\) 且系数是 \(\frac 12\)。那么设

\[\frac{f^2}2=-\ln(1-x)-x \]

显然 \(f\) 有复合逆，设为 \(g\)。那么

\[[x^n]e^{y\frac{f^2}2}=\frac 1n[x^{n-1}]\frac{\partial e^{y\frac{x^2}2}}{\partial x}\left(\frac xg\right)^n \]

观察中间

\[\begin{aligned} &\frac{\partial e^{y\frac{x^2}2}}{\partial x}\\ =&xye^{y\frac{x^2}2}\\ =&\sum_{i=0}\frac{x^{2i+1}y^{i+1}}{2^ii!} \end{aligned} \]

于是

\[\begin{aligned} &[x^ny^k]e^{y\frac{f^2}2}\\ =&\frac 1n[x^{n-1}y^k]\frac{\partial e^{y\frac{x^2}2}}{\partial x}\left(\frac xg\right)^n\\ =&\frac 1{n2^{k-1}(k-1)!}[x^{n-2k}]\left(\frac xg\right)^n \end{aligned} \]

看看怎么算 \(g\)。观察定义式

\[\frac{f^2}2=-\ln(1-x)-x \]

将 \(x\) 用 \(g\) 替换：

\[\frac{x^2}2=-\ln(1-g)-g \]

没写，因为懒。

P7440

我们要对所有的 \(n\) 计算

\[\sum_{i=1}^nF(i)[x^ky^i]e^{y(-\ln(1-x)-x)} \]

那设 \(f(x,y)=e^{y(-\ln(1-x)-x)}\)。看到这种东西感觉无从下手，试试转置原理。

一开始的转移矩阵显然是 \(a_{i,j}=[x^iy^j]f(x,y)\)。转置之后答案的生成函数为：

\[G=\sum_{i=1}^nF(i)[x^i]f(x,y) \]

考虑后边的 \([x^i]f(x,y)\)。设 \(f_i(y)=[x^i]f(x,y)\)，则求导可以得到递推式

\[f_i=\frac{i-1}if_{i-1}+\frac yif_{i-2} \]

那么这是个矩阵递推。转移矩阵为

\[A_i=\begin{bmatrix}\frac{i-1}i & \frac yi\\1 & 0\end{bmatrix} \]

那么答案的生成函数可以变成

\[G=\sum_{i=1}^nF(i)\prod_{j=1}^iA_j \]

把左上角拿出来。

这个实际上可以分治 FFT。设 \(P_{l,r}=\prod_{i=l}^rA_i,Q_{l,r}=\sum_{i=l}^rF(i)\prod_{j=l}^iA_j\)。转移显然：

\[P_{l,r}=P_{l,mid}P_{mid+1,r} \]

\[Q_{l,r}=Q_{l,mid}+Q_{mid+1,r}P_{l,mid} \]

然后将上述算法转置即可。转置后的矩阵乘就是把矩阵转置一下然后把输出当输入，乘法变转置乘法。

cnm，为什么卡常。

cnm，卡过了。

int n,k;
poly f,val;
struct node{
    poly a[2][2];
    void resize(int n){
        a[0][0].resize(n);a[0][1].resize(n);
        a[1][0].resize(n);a[1][1].resize(n);
    }
    node operator*(const node&s)const{
        node ans;
        ans.a[0][0]=a[0][0]*s.a[0][0]+a[0][1]*s.a[1][0];
        ans.a[0][1]=a[0][0]*s.a[0][1]+a[0][1]*s.a[1][1];
        ans.a[1][0]=a[1][0]*s.a[0][0]+a[1][1]*s.a[1][0];
        ans.a[1][1]=a[1][0]*s.a[0][1]+a[1][1]*s.a[1][1];
        return ans;
    }
    node operator^(const node&s)const{
        node ans;
        ans.a[0][0]=(a[0][0]^s.a[0][0])+(a[0][1]^s.a[0][1]);
        ans.a[0][1]=(a[0][0]^s.a[1][0])+(a[0][1]^s.a[1][1]);
        ans.a[1][0]=(a[1][0]^s.a[0][0])+(a[1][1]^s.a[0][1]);
        ans.a[1][1]=(a[1][0]^s.a[1][0])+(a[1][1]^s.a[1][1]);
        return ans;
    }
}Q[400010];
#define lson rt<<1
#define rson rt<<1|1
void solve1(int rt,int l,int r){
    if(l==r){
        int inv=qpow(l,mod-2);
        Q[rt].a[0][0].f.emplace_back(1ll*(l-1)*inv%mod);
        Q[rt].a[0][1].f.emplace_back(0);Q[rt].a[0][1].f.emplace_back(inv);
        Q[rt].a[1][0].f.emplace_back(1);
        Q[rt].a[1][1].f.emplace_back(0);
        return;
    }
    int mid=(l+r)>>1;
    solve1(lson,l,mid);solve1(rson,mid+1,r);
    if(r<n)Q[rt]=Q[rson]*Q[lson];
}
int ans[100010];
void solve2(int rt,int l,int r,node&P){
    P.resize(r-l+2);
    if(l==r){
        P=P^Q[rt];
        ans[l]=P.a[0][0][0];
        ans[l]=1ll*ans[l]*jc[l]%mod;
        return;
    }
    int mid=(l+r)>>1;
    node PP=P;
    solve2(lson,l,mid,PP);
    PP=P^Q[lson];
    solve2(rson,mid+1,r,PP);
}
int main(){
    n=read();k=read();
    init(max(n,k)+1<<1);
    f.resize(k);
    for(int i=0;i<k;i++)f[i]=read();
    val.resize(n+2);
    for(int i=1;i<=n+1;i++)val[i]=i-1;
    val=multipoint(f,val,k,n+1);val=val>>1;
    solve1(1,1,n);
    node P;
    P.a[0][0]=val;
    P.a[0][1].f.emplace_back(0);
    P.a[1][0].f.emplace_back(0);
    P.a[1][1].f.emplace_back(0);
    solve2(1,1,n,P);
    for(int i=1;i<=n;i++)print(ans[i]),putchar(' ');puts("");
    return 0;
}