Luogu6115

引言

这题我之前一直只是口胡着，没有来写，现在有空了来写一下。

做这题之前建议先去做一下快速阶乘算法。

容易发现快速阶乘算法其实就是整式递推的一个特例而已……

所以如果你已经会了快速阶乘算法，这道题目就很容易实现了。

以下默认你已经学会了快速阶乘算法。

---

思路

好了，你现在已经会了快速阶乘算法了。

快速阶乘算法的思想在于让机器帮你快速实现分段打表。

具体的实现方法上，其依赖于连续点值平移与值域倍增。

为了实现整式递推，我们应当注意到快速阶乘算法过程中，我们维护的点值其实并不是 \(n!\)，而是 \(\prod_{i=0}^{T-1}(aT+i)\)，即一对点值之间的倍数关系。

由于整式递推阶数 \(m\) 不止是 \(1\) 了，我们就不能直接维护一对数之间的倍数关系了；而是维护出一对 \(m\) 维向量之间的线性变换，即 \(m\times m\) 的一个矩阵，矩阵的每一项对应于某个多项式的一个点值。

（学过高等代数的同学应该知道，由于一元多项式环是整环，所以在多项式环上的矩阵有定义，通常称为 \(\lambda\)-矩阵，其可以把 \(\lambda\) 提出来写成多项式的形式。不过这不是重点。）

这个矩阵的形式不好构造成完全的多项式形式，我们不妨先令 \(P_0(n)=-1\)，那么相邻的向量之间的关系可以记为：

\[\begin{bmatrix}P_1(n)&P_2(n)&P_3(n)&\cdots&P_{m-1}(n)&P_m(n)\\1\\&1\\&&1\\&&&\ddots\\&&&&1\\\end{bmatrix} \begin{bmatrix}a_{n-1}\\a_{n-2}\\a_{n-3}\\\vdots\\a_{n-m+1}\\a_{n-m}\end{bmatrix} =\begin{bmatrix}a_n\\a_{n-1}\\a_{n-2}\\\vdots\\a_{n-m+2}\\a_{n-m+1}\end{bmatrix} \]

然后我们现在要考虑一般的 \(P_0\)，则为

\[-{1\over P_0(n)}\begin{bmatrix}P_1(n)&P_2(n)&P_3(n)&\cdots&P_{m-1}(n)&P_m(n)\\-P_0(n)\\&-P_0(n)\\&&-P_0(n)\\&&&\ddots\\&&&&-P_0(n)\\\end{bmatrix} \begin{bmatrix}a_{n-1}\\a_{n-2}\\a_{n-3}\\\vdots\\a_{n-m+1}\\a_{n-m}\end{bmatrix} =\begin{bmatrix}a_n\\a_{n-1}\\a_{n-2}\\\vdots\\a_{n-m+2}\\a_{n-m+1}\end{bmatrix} \]

设 \(B(\lambda)=\begin{bmatrix}P_1(\lambda)&P_2(\lambda)&P_3(\lambda)&\cdots&P_{m-1}(\lambda)&P_m(\lambda)\\-P_0(\lambda)\\&-P_0(\lambda)\\&&-P_0(\lambda)\\&&&\ddots\\&&&&-P_0(\lambda)\\\end{bmatrix} \)，我们先撇开前面的 \(-\frac1{P_0(n)}\) 因子不论，我们现在要维护 \(\prod_{i=0}^{T-1}B(aT+m+i)\) 这种形式的量，其中乘法自右往左。

容易发现 \(B_T(\lambda)=\prod_{i=0}^{T-1}B(\lambda+i)\) 是一个各项次数不高于 \(dT\) 的 \(\lambda\)-矩阵，只用 \(dT+1\) 个值即可维护。

于是我们维护出 \(B_T(m)\)，\(B_T(m+T)\)，\(B_T(m+2T)\)，\(\dots\)，\(B_T(m+(dT-1)T)\)，\(B_T(m+dT^2)\) 这几个 \(\lambda\)-矩阵点值，然后用类似于快速阶乘算法的方式暴力进行多项式点值平移和倍增就好了。

初值 \(T=1\)，\(B_1(\lambda)=B(\lambda)\)，要维护 \(d+1\) 处点值。

最后把原向量逐个乘过去即可。

完了？

别忘记还有 \(-\frac1{P_0(n)}\) 因子！

现在要求 \(\prod_{i=m}^n-P_0(i)\)，直接再同上的方式做一次快速阶乘算法式的倍增即可。

---

Code

常数太大，没过，就不贴了。

---

Update on 2022.10.10：

用循环卷积卡常卡过去了。

下面是核心代码。

const ullt Mod=998244353,g=3;
typedef ConstMod::mod_ullt<Mod>modint;
typedef std::vector<modint>modvec;
typedef NTT_POLY::poly_NTT<Mod,g>poly;
typedef NTT_POLY::poly_eval<Mod,g>eval;
typedef NTT_POLY::poly_inter<Mod,g>inter;
typedef NTT_POLY::poly_cpd<Mod,g>cpd;
typedef NTT_POLY::poly_nums<Mod,g>nums;
typedef FWT_MODINT::FWT_Mod<Mod>FWT;
modint P[2000005],Q[2000005];
poly::DIFDIT s;
modvec shift2(modvec A,uint n,modint v){
    static modint User[4000005],UserInv[4000005];
    User[0]=1;
    for(uint i=0;i<=n*2;i++)User[i+1]=User[i]*(v-n+i);
    UserInv[n*2+1]=User[n*2+1].inv();
    for(uint i=n*2;~i;i--)UserInv[i]=UserInv[i+1]*(v-n+i);
    modvec G;for(uint i=0;i<=n*2;i++)G.push_back(UserInv[i+1]*User[i]);
    s.dif(G);for(uint i=0;i<s.size();i++)G[i]*=A[i];
    s.dit(G);
    modvec ans(n+1);
    for(uint i=0;i<=n;i++)ans[i]=G[i+n]*UserInv[i]*User[i+n+1];
    return ans;
}
modint A[10];poly W[10];
modvec B[10][10],B2;
int main()
{
#ifdef MYEE
    freopen("QAQ.in","r",stdin);
    // freopen("QAQ.out","w",stdout);
#endif
    P[0]=1;for(uint i=1;i<=2000000;i++)P[i]=P[i-1]*i;
    Q[2000000]=P[2000000].inv();for(uint i=2000000;i;i--)Q[i-1]=Q[i]*i;
    uint n,m,d;scanf("%u%u%u",&n,&m,&d);
    uint T=2;while(m+T*T*d<=n)T<<=1; // T \sim \sqrt n/d
    for(uint i=0;i<m;i++)A[i].read();
    std::reverse(A,A+m);
    for(uint i=0;i<=m;i++){W[i].chg_deg(d);for(uint j=0;j<=d;j++)W[i][j].read();}
    for(uint i=0;i<m;i++)for(uint j=0;j<m;j++)B[i][j].resize(d+1);
    B2.resize(d+1);
    for(uint i=0;i<=d;i++)
    {
        for(uint j=0;j<m;j++)B[j][0][i]=cpd().point_eval(W[j+1],m+i);
        modint v=B2[i]=-cpd().point_eval(W[0],m+i);
        for(uint j=1;j<m;j++)B[j-1][j][i]=v;
    }
    for(uint t=2;t<=T;t<<=1){
        static modvec X,Y,Z,W;s.bzr(d*t+1); // 循环卷积卡常
        for(uint i=0;i<m;i++)for(uint j=0;j<m;j++){
            W=B[i][j];
            for(uint i=0;i<=d*(t/2);i++)W[i]*=Q[i]*Q[d*(t/2)-i]*((d*(t/2)-i)&1?-modint(1):1);
            s.dif(W),X=shift2(W,d*(t/2),d*(t/2)+1),Y=shift2(W,d*(t/2),d*t+2),Z=shift2(W,d*(t/2),d*(t/2)*3+3);
            for(auto v:X)B[i][j].push_back(v);
            for(auto v:Y)B[i][j].push_back(v);
            for(auto v:Z)B[i][j].push_back(v);
        }
        static modint QwQ[10][10];
        for(uint i=0;i<m;i++)for(uint j=0;j<m;j++)QwQ[i][j]=0;
        for(uint i=0;i<m;i++)for(uint j=0;j<m;j++)for(uint k=0;k<m;k++)QwQ[i][k]+=B[i][j][0]*B[j][k][1];
        for(uint i=0;i<m;i++)for(uint j=0;j<m;j++)B[i][j][0]=QwQ[i][j];
        for(uint dep=1;dep<=d*t;dep++){
            for(uint i=0;i<m;i++)for(uint j=0;j<m;j++)B[i][j][dep]=0;
            for(uint i=0;i<m;i++)for(uint j=0;j<m;j++)for(uint k=0;k<m;k++)
                B[i][k][dep]+=B[i][j][dep<<1]*B[j][k][dep<<1|1];
        }
        for(uint i=0;i<m;i++)for(uint j=0;j<m;j++)B[i][j].resize(t*d+1);
        W=B2;
        for(uint i=0;i<=d*(t/2);i++)W[i]*=Q[i]*Q[d*(t/2)-i]*((d*(t/2)-i)&1?-modint(1):1);
        s.dif(W),X=shift2(W,d*(t/2),d*(t/2)+1),Y=shift2(W,d*(t/2),d*t+2),Z=shift2(W,d*(t/2),d*(t/2)*3+3);
        for(auto v:X)B2.push_back(v);
        for(auto v:Y)B2.push_back(v);
        for(auto v:Z)B2.push_back(v);
        for(uint dep=0;dep<=d*t;dep++)B2[dep]=B2[dep<<1]*B2[dep<<1|1];
        B2.resize(d*t+1);
    }
    modint ans(1);
    for(uint i=0;i<(n-m)/T;i++){
        modvec User(m);
        for(uint j=0;j<m;j++)for(uint k=0;k<m;k++)User[k]+=B[j][k][i]*A[j];
        for(uint j=0;j<m;j++)A[j]=User[j];
        ans*=B2[i];
    }
    for(uint i=(n-m)/T*T+m;i<=n;i++){
        modint user;
        for(uint j=0;j<m;j++)user+=cpd().point_eval(W[j+1],i)*A[j];
        modint v=-cpd().point_eval(W[0],i);
        for(uint j=m-1;j;j--)A[j]=A[j-1]*v;
        A[0]=user,ans*=v;
    }
    ans=A[0]/ans;
    ans.println();
    return 0;
}