[数学记录]P7438~7440 「KrOI2021」Feux Follets 三连

三题三科技，好！（拍手）大棒子，学到许多！

## 二次弱化版

**题意** ： 设 $\text{cyc}_p$ 将长为 $n$ 的排列 $p$ 当成置换时所能分解成的循环个数。给定两个整数 $n,k$ 和一个 $k-1$ 次多项式，对 $1\leq m\leq n$ 分别求：

$$
\sum\limits_{p}F(\text{cyc}_{p})
$$

其中 $p$ 是长度为 $m$ 且不存在位置 $i$ 使得 $p_i=i$ 的排列。

$n\leq 6\times 10^5,k\leq100$ ，时限$\texttt{0.8s}$。

------------

考虑对 $1\leq i\leq n$ 分别求出循环个数为 $i$ 的错排的个数。然后多点求值得到贡献系数即可算出答案。

考虑如何用生成函数刻画错排。错排可以看做大小 $>1$ 的循环组成的无序集合。

大小 $>1$ 的循环的 $\rm EGF$ ： $-x-\ln(1-x)$

错排的 $\rm EGF$ : $e^{-x-\ln(1-x)}$

我们要根据循环个数分类，于是改写为二元函数 $G(x,y)=e^{\small y(-x-\ln(1-x))}$。

我们的目标即是求出 $[x^n]e^{\small y(-x-\ln(1-x))}$ ，然而 $y$ 的次数上界和 $n$ 有关，没有利用题目中 $k$ 较小的性质。

将多项式 $F$ 改写为牛顿级数，这部分容易 $O(k^2)$ 暴力完成。

然后分项求和，即对于 $0\leq t\leq k$ ，求

$$
\sum\limits_{p}\dbinom{\text{cyc}_{p}}{t}
$$

我们只需将生成函数改写为 $G(x,y)=e^{\small (1+y)(-x-\ln(1-x))}$ 。其中的 $1+y$ 表示对于每个循环可选可不选。

此时，求出 $[x^n]e^{\small (1+y)(-x-\ln(1-x))}$ 再与 $F$ 的系数（只有 $k$ 个）点乘即可得到最终答案。

考虑系数递推。

$$
\begin{aligned}
\dfrac{\partial}{\partial x}G(x,y)&=(1+y)\frac{x}{1-x}G(x,y)\\
[x^n]\dfrac{\partial}{\partial x}G(x,y)&=[x^n](1+y)\frac{x}{1-x}G(x,y)\\
(n+1)[x^{n+1}]G(x,y)&=[x^{n-1}](1+y)\frac{1}{1-x}G(x,y)\\
(n+1)[x^{n+1}y^m]G(x,y)&=[x^{n-1}y^{m-1}]\frac{1}{1-x}G(x,y)+[x^{n-1}y^m]\frac{1}{1-x}G(x,y)\\
(n+1)[x^{n+1}y^m]G(x,y)&=\sum\limits_{i=0}^{n-1}\Big([x^iy^{m-1}]G(x,y)+[x^iy^m]G(x,y)\Big)
\end{aligned}
$$

边界为 $[x^0]G=1,[x^1]G=0$。

按照上式递推即可，使用前缀和优化并滚动数组，复杂度为 $O(nk)$。

```cpp
#include<algorithm>
#include<cstdio>
#define ll long long
#define MaxN 600500
#define MaxK 105
using namespace std;
const int mod=998244353;
ll powM(ll a,int t=mod-2){
  ll ret=1;
  while(t){
    if (t&1)ret=ret*a%mod;
    a=a*a%mod;t>>=1;
  }return ret;
}
ll fac[MaxN],ifac[MaxN];
ll C(int n,int m)
{return fac[n]*ifac[m]%mod*ifac[n-m]%mod;}
void Init(int n)
{
  fac[0]=1;
  for (int i=1;i<=n;i++)
    fac[i]=fac[i-1]*i%mod;
  ifac[n]=powM(fac[n]);
  for (int i=n;i;i--)
    ifac[i-1]=ifac[i]*i%mod;
}
void CompT(ll *W,int n)
{
  static ll Y[MaxK];
  for (int k=0;k<n;k++){
    int buf=1;
    for (int i=0;i<n;i++){
      Y[k]=(Y[k]+W[i]*buf);
      buf=1ll*buf*k%mod;
    }
  }
  for (int k=0;k<n;k++){
    W[k]=0;
    for (int i=0;i<=k;i++){
      if ((k-i)&1)W[k]=(W[k]-C(k,i)*Y[i])%mod;
      else W[k]=(W[k]+C(k,i)*Y[i])%mod;
    }
  }
}
ll W[MaxK],f[MaxK],g[MaxK],s[MaxK];
int n,k;
int main()
{
  scanf("%d%d",&n,&k);Init(max(n,k));
  for (int i=0;i<k;i++)scanf("%lld",&W[i]);
  CompT(W,k);
  printf("0 ");g[0]=1;
  for (int i=2;i<=n;i++){
    ll ans=0,buf=ifac[i]*fac[i-1]%mod;
    f[0]=g[0]*buf%mod;
    for (int j=1;j<k;j++)f[j]=(g[j]+g[j-1])*buf%mod;
    for (int j=0;j<k;j++)g[j]=(g[j]+s[j])%mod;
    for (int j=0;j<k;j++)s[j]=f[j];
    for (int j=0;j<k;j++)ans=(ans+W[j]*f[j])%mod;
    printf("%lld ",ans*fac[i]%mod);
  }return 0;
}
```

## 弱化版

**题意** ： 给定两个整数 $n,k$ 和一个 $k-1$ 次多项式，求

$$
\sum\limits_{p}F(\text{cyc}_{p})
$$

其中 $p$ 是长度为 $n$ 的错排。

$n,k\leq 10^5$，时限$\texttt{2s}$。

------------

此时 $k$ 较大，故不一定需要使用牛顿级数。下文默认 $n,k$ 同阶。

首先用多点求值得到 $F(1)...F(n)$ ，然后求出 $[x^n]e^{\small y(-x-\ln(1-x))}$ 的系数，点乘即可得到答案。

> 因此，在加强后的问题中，多项式 $F$ 的存在仅仅是提供了一个贡献系数序列……

考虑拉格朗日反演来快速求出 $[x^n]e^{\small y(-x-\ln(1-x))}$。

首先用多项式复合描述问题。设 $G(x)=-x-\ln(1-x),H(x)=e^{xy}$ ，则 $F(x,y)=e^{\small y(-x-\ln(1-x))}=e^{yG(x)}=H(G(x))$。

此处 $G(x)$ 并不存在复合逆，但有 $G[2]=\frac{1}{2}$，于是设 $T(x)=\sqrt{2G(x)}$ ，此时 $T[1]=1$ ，故 $T(x)$ 显然有复合逆，设为 $T_2(x)$。

并改写 $H(x)=e^{x^2y/2}$ ，于是仍有 $F(x,y)=H(T(x))$。

套用扩展拉格朗日反演，可得 ：

$$[x^n]H(T(x))=\frac{1}{n}[x^{n-1}]H'(x)\Big(\dfrac{x}{T_2(x)}\Big)^{n}$$

算得 $H'(x)=xye^{x^2y/2}=\sum\limits_{i=0}\dfrac{x^{2i+1}y^{i+1}}{2^ii!}$ 。记 $P(x)=\Big(\dfrac{x}{T_2(x)}\Big)^{n}$。

$$
\begin{aligned}
[x^n]H(T(x))&=\frac{1}{n}[x^{n-1}]\sum\limits_{i=0}\dfrac{x^{2i+1}y^{i+1}}{2^ii!}P(x)\\
[x^ny^m]H(T(x))&=\frac{1}{n}[x^{n-1}y^m]\sum\limits_{i=0}\dfrac{x^{2i+1}y^{i+1}}{2^ii!}P(x)\\
[x^ny^m]H(T(x))&=\frac{1}{n}[x^{n-1}]\dfrac{x^{2(m-1)+1}}{2^{m-1}(m-1)!}P(x)\\
[x^ny^m]H(T(x))&=\dfrac{1}{n2^{m-1}(m-1)!}[x^{n-2m}]P(x)\\
\end{aligned}
$$

于是，在求出 $T_2(x)$ 之后可以 $O(n\log n)$ 得到 $[x^n]F(x,y)$。

注意到 $T$ 有一个较为简单的封闭形式，考虑利用牛顿迭代求复合逆 $T_2$。

$$
\begin{aligned}
T^2/2&=-x-\ln(1-x)\\
x^2/2&=-T_2-\ln(1-T_2)\\
2T_2+2\ln(1-T_2)+x^2&=0
\end{aligned}
$$

于是可以牛顿迭代，倍增式为 ：

$$T_2=T_2^*-\dfrac{2T_2^*+2\ln(1-T_2^*)+x^2}{2-2(1-T_2^*)^{-1}}$$

注意分母 $2-2(1-T_2^*)^{-1}$ 常数项为 $0$ ，故做除法时需要分子分母先除去一个 $x$。这会导致次数并非严格倍增。

复杂度同为 $O(n\log n)$。

多点求值的复杂度为 $O(n\log^2n)$ ，注意常数优化。

```cpp
#include<algorithm>
#include<cstring>
#include<cstdio>
#include<vector>
#define ll long long
#define ull unsigned ll
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n))
#define pb push_back
const int _G=3,mod=998244353,inv2=(mod+1)/2,Maxn=1<<17|500;
using namespace std; 
ll powM(ll a,ll t=mod-2){
  ll ans=1;
  while(t){
    if(t&1)ans=ans*a%mod;
    a=a*a%mod;t>>=1;
  }return ans;
}
const int invG=powM(_G);
int tr[Maxn<<1],tf;
void tpre(int n){
  if (tf==n)return ;tf=n;
  for(int i=0;i<n;i++)
    tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
}
void NTT(int *g,bool op,int n)
{
  tpre(n);
  static ull f[Maxn<<1],w[Maxn];w[0]=1;
  for (int i=0;i<n;i++)f[i]=(((ll)mod<<4)+g[tr[i]])%mod;
  for(int l=1;l<n;l<<=1){
    ull tG=powM(op?_G:invG,(mod-1)/(l+l));
    for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
    for(int k=0;k<n;k+=l+l)
      for(int p=0;p<l;p++){
        int tt=w[p]*f[k|l|p]%mod;
        f[k|l|p]=f[k|p]+mod-tt;
        f[k|p]+=tt;
      }
  }if (!op){
    ull invn=powM(n);
    for(int i=0;i<n;++i)
      g[i]=f[i]%mod*invn%mod;
  }else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
void px(int *f,int *g,int n)
{for(int i=0;i<n;++i)f[i]=1ll*f[i]*g[i]%mod;}
#define Poly vector<int>
Poly operator + (const Poly &A,const Poly &B)
{
  Poly C=A;C.resize(max(A.size(),B.size()));
  for (int i=0;i<B.size();i++)C[i]=(C[i]+B[i])%mod;
  return C;
}
Poly operator - (const Poly &A,const Poly &B)
{
  Poly C=A;C.resize(max(A.size(),B.size()));
  for (int i=0;i<B.size();i++)C[i]=(C[i]+mod-B[i])%mod;
  return C;
}
Poly operator * (int c,const Poly &A)
{
  Poly C;C.resize(A.size());
  for (int i=0;i<A.size();i++)C[i]=1ll*c*A[i]%mod;
  return C;
}
Poly operator * (const Poly &A,const Poly &B)
{
  Poly C;C.resize(A.size()+B.size()-1);
  if (min(A.size(),B.size())<=500){
    for (int i=0;i<A.size();i++)
      for (int j=0;j<B.size();j++)
        C[i+j]=(C[i+j]+1ll*A[i]*B[j])%mod;
  }else {
    static int a[Maxn<<1],b[Maxn<<1];
    cpy(a,&A[0],A.size());
    cpy(b,&B[0],B.size());
    int n=1;for(;n<C.size();n<<=1);
    NTT(a,1,n);NTT(b,1,n);
    px(a,b,n);NTT(a,0,n);
    cpy(&C[0],a,C.size());
    clr(a,n);clr(b,n);
  }return C;
}
Poly mulT(Poly &A,const Poly &B)
{
  reverse(A.begin(),A.end());
  Poly C=A*B;C.resize(A.size());
  reverse(C.begin(),C.end());reverse(A.begin(),A.end());
  return C;
}  
void pinv(const Poly &A,Poly &B,int n)
{
  if (n==1)B.pb(powM(A[0]));
  else {
    pinv(A,B,(n+1)/2);
    Poly sA;sA.resize(n);
    cpy(&sA[0],&A[0],n);
    B=2*B-B*B*sA;
    B.resize(n);
  }
}
Poly pinv(const Poly &A)
{Poly C;pinv(A,C,A.size());return C;}
int inv[Maxn];
void Init(int n)
{
  inv[1]=1;
  for (int i=2;i<=n;i++)
    inv[i]=1ll*inv[mod%i]*(mod-mod/i)%mod;
}
Poly dao(const Poly &A)
{
  Poly C=A;
  for (int i=1;i<C.size();i++)
    C[i-1]=1ll*C[i]*i%mod;
  C.pop_back();
  return C;
}
Poly ints(const Poly &A)
{
  Poly C=A;
  C.pb(0);
  for (int i=C.size()-1;i;i--)
    C[i]=1ll*C[i-1]*inv[i]%mod;
  C[0]=0;
  return C;
}
Poly ln(const Poly &A)
{
  Poly B=dao(A)*pinv(A);
  B.resize(A.size());
  return ints(B);
}
void pexp(const Poly &A,Poly &B,int n)
{
  if (n==1)B.pb(1);
  else {
    pexp(A,B,(n+1)/2);
    Poly lnB=B;
    lnB.resize(n);lnB=ln(lnB);
    for (int i=0;i<lnB.size();i++)
      lnB[i]=(mod+A[i]-lnB[i])%mod;
    lnB[0]++;
    B=B*lnB;
    B.resize(n);
  }
}
Poly pexp(const Poly &A)
{Poly C;pexp(A,C,A.size());return C;}
void strT(Poly &B,int n)
{
  if (n==2){B.pb(0);B.pb(1);}
  else {
    strT(B,n/2+1);
    Poly Q=B;Q[0]=(Q[0]+mod-1)%mod;
    Q=2*pinv(Q);Q.erase(Q.begin());Q=pinv(Q);
    Poly P=(mod-1)*B;P[0]++;
    P.resize(n+1);P=2*ln(P);
    P=P+2*B;P[2]++;
    P.erase(P.begin());
    B=B-P*Q;B.resize(n);
  }
}
Poly Q[Maxn],P[22];
void solve1(int l,int r,int u,Poly &X)
{
  if (l==r){
    Q[u].pb(1);
    Q[u].pb(mod-X[l]);
    return ;
  }int mid=(l+r)>>1;
  solve1(l,mid,u<<1,X);
  solve1(mid+1,r,u<<1|1,X);
  Q[u]=Q[u<<1]*Q[u<<1|1];
}
void solve2(int l,int r,int u,int dep,Poly &Y)
{
  if (l==r){Y[l]=P[dep][0];return ;}
  int mid=(l+r)>>1;
  P[dep].resize(r-l+1);
  P[dep+1]=mulT(P[dep],Q[u<<1|1]);
  solve2(l,mid,u<<1,dep+1,Y);
  P[dep+1]=mulT(P[dep],Q[u<<1]);
  solve2(mid+1,r,u<<1|1,dep+1,Y);
}
void calcP(Poly &F,Poly &X)
{
  solve1(0,X.size()-1,1,X);
  Q[1].resize(max(F.size(),X.size()));
  P[0]=mulT(F,pinv(Q[1]));
  solve2(0,X.size()-1,1,0,X);
}
int n,k;Poly F,W,T;
int main()
{
  scanf("%d%d",&n,&k);
  F.resize(k);Init(n);
  for (int i=0;i<k;i++)scanf("%d",&F[i]);
  W.resize(n/2+1);
  for (int i=0;i<=n/2;i++)W[i]=i;
  calcP(F,W);
  strT(T,n);
  T.erase(T.begin());
  T=pexp((mod-n)*ln(T));
  ll ans=0,ifac=1,pw2=1,fac=1;
  for (int i=1;i<n;i++)fac=fac*i%mod;
  for (int i=1;i<=n/2;i++){
    ans=(ans+1ll*W[i]*T[n-2*i]%mod*ifac%mod*pw2)%mod;
    ifac=ifac*inv[i]%mod;
    pw2=pw2*inv2%mod;
  }printf("%lld",ans*fac%mod);
	return 0;
}
```

## 本体

**题意** ： 给定两个整数 $n,k$ 和一个 $k-1$ 次多项式，对 $1\leq m\leq n$ 分别求

$$
\sum\limits_{p}F(\text{cyc}_{p})
$$

其中 $p$ 是长度为 $m$ 的错排。

$n,k\leq 10^5$，时限$\texttt{2s}$。

------------

前置芝士 ：[转置原理小记](https://www.luogu.com.cn/blog/command-block/zhuai-zhi-yuan-li-xiao-ji)

仍然首先搞一手多点求值，则有

$${\rm Ans}[k]=\sum\limits_{i}F(i)[x^ky^i]e^{\small y(-x-\ln(1-x))}$$

记 $s_i=F(i),\ G(x,y)=e^{\small y(-x-\ln(1-x))}$。

以线性变换的角度来看，将 ${\rm Ans}$ 视为输出向量，$s$ 视为输入向量，则贡献矩阵 $A_{k,i}=[x^ky^i]G(x,y)$。

考虑该变换的转置，即 $t=A^Ts$ ，写出和式 ：

$$t_k=\sum\limits_{i}s_i[x^iy^k]G(x,y)$$

考虑 $t_k$ 的生成函数 $T(y)$ ，则可写作 $T(y)=\sum\limits_{i}s_i[x^i]G(x,y)$

将 $G(x,y)$ 按 $x$ 切分，设 $G_i(y)=[x^i]G(x,y)$。

求导可得
  
$$\dfrac{\partial}{\partial x}G(x,y)=y\frac{x}{1-x}G(x,y)$$

由此不难得到递推式

$$G_{i}=\dfrac{i-1}{i}G_{i-1}+\dfrac{y}{i}G_{i-2}$$

边界为 $G_0=1$ ，为了方便，令 $G_{-1}=0$ ，这样就能一直套用递推式。

将转移写成矩阵 $W_i=\begin{bmatrix}0&1\\\frac{y}{i}&\frac{i-1}{i}\end{bmatrix}$ ，满足 $\begin{bmatrix}G_{i-1}\\G_i\end{bmatrix}=W_i[G_{i-2},G_{i-1}]$

则有 $\begin{bmatrix}G_{i-1}\\G_i\end{bmatrix}=W_iW_{i-1}...W_1[0,1]$ ，于是 $W_iW_{i-1}...W_1=\begin{bmatrix}?&G_{i-1}\\?&G_{i}\end{bmatrix}$

则答案可以改写为 ：
$$\sum\limits_{i=1}^ns_iW_iW_{i-1}..W_1$$

最终取出矩阵右下角即为答案。

这可以利用分治 FFT 计算。

设 $P_{[l,r)}=W_{r-1}W_{r-2}...W_l,\ Q_{[l,r)}=\sum\limits_{i\in[l,r)}s_iW_{i}W_{i-1}...W_{l}$

边界为 $P_{[i,i+1)}=s_iW_i$。

有转移 ：

$$
\begin{aligned}
P_{[l,r)}&=P_{[l,m)}P_{[m,r)}\\
Q_{[l,r)}&=Q_{[l,m)}+Q_{[m,r)}P_{[l,m)}
\end{aligned}
$$

注意矩阵乘法不满足交换律，故转移中的乘法顺序必须严格考虑。

将该算法转置即可。

矩阵乘法的转置相当于用转置乘法计算转置后矩阵的乘法。

将输入塞在原来答案的位置，输出就会出现在原来输入的位置。

具体计算方法可见代码。

复杂度为 $O(k\log^2k+n\log^2n)$ ，常数较大。

```cpp
#include<algorithm>
#include<cstring>
#include<cstdio>
#include<vector>
#define ll long long
#define ull unsigned ll
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n))
#define pb push_back
const int _G=3,mod=998244353,inv2=(mod+1)/2,Maxn=1<<17|500;
using namespace std; 
ll powM(ll a,ll t=mod-2){
  ll ans=1;
  while(t){
    if(t&1)ans=ans*a%mod;
    a=a*a%mod;t>>=1;
  }return ans;
}
const int invG=powM(_G);
int tr[Maxn<<1],tf;
void tpre(int n){
  if (tf==n)return ;tf=n;
  for(int i=0;i<n;i++)
    tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
}
void NTT(int *g,bool op,int n)
{
  tpre(n);
  static ull f[Maxn<<1],w[Maxn];w[0]=1;
  for (int i=0;i<n;i++)f[i]=(((ll)mod<<4)+g[tr[i]])%mod;
  for(int l=1;l<n;l<<=1){
    ull tG=powM(op?_G:invG,(mod-1)/(l+l));
    for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
    for(int k=0;k<n;k+=l+l)
      for(int p=0;p<l;p++){
        int tt=w[p]*f[k|l|p]%mod;
        f[k|l|p]=f[k|p]+mod-tt;
        f[k|p]+=tt;
      }
  }if (!op){
    ull invn=powM(n);
    for(int i=0;i<n;++i)
      g[i]=f[i]%mod*invn%mod;
  }else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
void px(int *f,int *g,int n)
{for(int i=0;i<n;++i)f[i]=1ll*f[i]*g[i]%mod;}
#define Poly vector<int>
Poly operator + (const Poly &A,const Poly &B)
{
  Poly C=A;C.resize(max(A.size(),B.size()));
  for (int i=0;i<B.size();i++)C[i]=(C[i]+B[i])%mod;
  return C;
}
Poly operator - (const Poly &A,const Poly &B)
{
  Poly C=A;C.resize(max(A.size(),B.size()));
  for (int i=0;i<B.size();i++)C[i]=(C[i]+mod-B[i])%mod;
  return C;
}
Poly operator * (int c,const Poly &A)
{
  Poly C;C.resize(A.size());
  for (int i=0;i<A.size();i++)C[i]=1ll*c*A[i]%mod;
  return C;
}
Poly operator * (const Poly &A,const Poly &B)
{
  Poly C;C.resize(A.size()+B.size()-1);
  if (min(A.size(),B.size())<=50){
    for (int i=0;i<A.size();i++)
      for (int j=0;j<B.size();j++)
        C[i+j]=(C[i+j]+1ll*A[i]*B[j])%mod;
  }else {
    static int a[Maxn<<1],b[Maxn<<1];
    cpy(a,&A[0],A.size());
    cpy(b,&B[0],B.size());
    int n=1;for(;n<C.size();n<<=1);
    NTT(a,1,n);NTT(b,1,n);
    px(a,b,n);NTT(a,0,n);
    cpy(&C[0],a,C.size());
    clr(a,n);clr(b,n);
  }return C;
}
Poly mulT(Poly &A,const Poly &B)
{
  reverse(A.begin(),A.end());
  Poly C=A*B;C.resize(A.size());
  reverse(C.begin(),C.end());reverse(A.begin(),A.end());
  return C;
}  
void pinv(const Poly &A,Poly &B,int n)
{
  if (n==1)B.pb(powM(A[0]));
  else {
    pinv(A,B,(n+1)/2);
    Poly sA;sA.resize(n);
    cpy(&sA[0],&A[0],n);
    B=2*B-B*B*sA;
    B.resize(n);
  }
}
Poly pinv(const Poly &A)
{Poly C;pinv(A,C,A.size());return C;}
namespace Eva
{ 
Poly Q[Maxn<<1],P[22];
void solve1(int l,int r,int u,Poly &X)
{
  if (l==r){
    Q[u].pb(1);
    Q[u].pb(mod-X[l]);
    return ;
  }int mid=(l+r)>>1;
  solve1(l,mid,u<<1,X);
  solve1(mid+1,r,u<<1|1,X);
  Q[u]=Q[u<<1]*Q[u<<1|1];
}
void solve2(int l,int r,int u,int dep,Poly &Y)
{
  if (l==r){Y[l]=P[dep][0];return ;}
  int mid=(l+r)>>1;
  P[dep].resize(r-l+1);
  P[dep+1]=mulT(P[dep],Q[u<<1|1]);
  solve2(l,mid,u<<1,dep+1,Y);
  P[dep+1]=mulT(P[dep],Q[u<<1]);
  solve2(mid+1,r,u<<1|1,dep+1,Y);
}
void calcP(Poly &F,Poly &X)
{
  solve1(0,X.size()-1,1,X);
  Q[1].resize(max(F.size(),X.size()));
  P[0]=mulT(F,pinv(Q[1]));
  solve2(0,X.size()-1,1,0,X);
}
}; 
struct PolyMat{
  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);}
}P[22];
PolyMat operator * (const PolyMat &A,const PolyMat &B)
{
  return (PolyMat){{{
    A.a[0][0]*B.a[0][0]+A.a[0][1]*B.a[1][0],
    A.a[0][0]*B.a[0][1]+A.a[0][1]*B.a[1][1]},{
    A.a[1][0]*B.a[0][0]+A.a[1][1]*B.a[1][0],
    A.a[1][0]*B.a[0][1]+A.a[1][1]*B.a[1][1]}
  }};
}
PolyMat mulT(PolyMat &A,const PolyMat &B)
{
  return (PolyMat){{{
    mulT(A.a[0][0],B.a[0][0])+mulT(A.a[0][1],B.a[0][1]),
    mulT(A.a[0][0],B.a[1][0])+mulT(A.a[0][1],B.a[1][1])},{
    mulT(A.a[1][0],B.a[0][0])+mulT(A.a[1][1],B.a[0][1]),
    mulT(A.a[1][0],B.a[1][0])+mulT(A.a[1][1],B.a[1][1])}
  }};
}
PolyMat Q[Maxn<<1];
int n;
void solve1(int l,int r,int u)
{
  if (l==r){
    int sav=powM(l);
    Q[u].a[0][0].pb(0);
    Q[u].a[0][1].pb(1);
    Q[u].a[1][0].pb(0);Q[u].a[1][0].pb(sav);
    Q[u].a[1][1].pb(1ll*(l-1)*sav%mod);
    return ;
  }int mid=(l+r)>>1;
  solve1(l,mid,u<<1);
  solve1(mid+1,r,u<<1|1);
  if (r<n)Q[u]=Q[u<<1|1]*Q[u<<1];
}
int ans[Maxn],fac[Maxn];
void Init(int n){
  fac[0]=1;
  for (int i=1;i<=n;i++)
    fac[i]=1ll*fac[i-1]*i%mod; 
}
void solve2(int l,int r,int u,int dep)
{
  P[dep].resize(r-l+2);
  if (l==r){
    P[dep]=mulT(P[dep],Q[u]);
    ans[l]=(P[dep].a[0][0][0]+P[dep].a[1][1][0])%mod;
    return ;
  }int mid=(l+r)>>1;
  P[dep+1]=P[dep];
  solve2(l,mid,u<<1,dep+1);
  P[dep+1]=mulT(P[dep],Q[u<<1]);
  solve2(mid+1,r,u<<1|1,dep+1);
}
int k;Poly F,W,T;
int main()
{
  scanf("%d%d",&n,&k);Init(n);
  F.resize(k);
  for (int i=0;i<k;i++)scanf("%d",&F[i]);
  W.resize(n+1);
  for (int i=0;i<=n;i++)W[i]=i;
  Eva::calcP(F,W);
  solve1(1,n,1);
  P[0].a[0][0].pb(0);
  P[0].a[0][1].pb(0);
  P[0].a[1][0].pb(0);
  P[0].a[1][1]=W;
  solve2(1,n,1,0);
  for (int i=1;i<=n;i++)printf("%d ",1ll*ans[i]*fac[i]%mod);
	return 0;
}
```