洛谷 7440 「KrOI2021」Feux Follets

首先把给定的多项式转成牛顿级数，即转写成 F(x)=k−1∑i=0ai(xi)

让我们考虑一个错排，显然它是由非自环的循环置换为基本单位构成的，即 e−x−ln(1−x)

然后考虑从中选 k 个。不难发现这只需要对一个循环置换附加一个因子 (1+y) 即可做到： e(1+y)(−x−ln(1−x))

令其为 G，考虑求算 k−1∑i=0ai[yi]G

直接计算是困难的；但是我们不妨考虑应用转置原理，首先考虑 k−1∑i=0ai[xi]G

设 Fi=[xi]G，对 x 求导可以得到 ∂G∂x=(1+y)x1−xG

即 iFi=(i−1)Fi−1+(1+y)Fi−2

考虑构造矩阵满足 [FiFi−1]=[Fi−1Fi−2]Ai

不妨假定 F−1=1。

那么我们要求的其实就是 n−1∑i=0[ai0]A1A2⋯Ai

用线性算法描述就是 A[[a00][a10]⋮[an−10]]

其中 Aij=[yi]A0A1⋯Aj

我们在计算此问题时，使用分治，同时维护当前区间乘出来的列向量 res，合并时 res=resL+ΠLresR

而在递归到长度为 1 的区间时，让 res 返回其左乘 ai 的向量的结果。

其中 ΠL 表示左区间的 Ai 之积。

那么我们的转置其实就已经明晰了。分治时传入一个 2 维行向量的列向量 res，分治时 resL←res,resR←res×TΠTL

然后递归进入子问题计算。

时间复杂度瓶颈为转牛顿级数，总复杂度 O(klog2k+nlog2n)。

代码：

#include <cstdio>
#include <vector>
#include <cstring>
#include <algorithm>
#define add(a,b) (a + b >= mod ? a + b - mod : a + b)
#define dec(a,b) (a < b ? a - b + mod : a - b)
#define ls (p << 1)
#define rs (ls | 1)
using namespace std;
const int N = 1e5;
const int mod = 998244353;
inline int fpow(int a,int b)
{
    int ret = 1;
    for(;b;b >>= 1)
        (b & 1) && (ret = (long long)ret * a % mod),a = (long long)a * a % mod;
    return ret;
}
int n,k;
namespace Poly
{
    const int LG = 17;
    const int N = 1 << LG + 1;
    const int G = 3;
    int lg2[N + 5];
    int rev[N + 5],fac[N + 5],ifac[N + 5],inv[N + 5];
    int rt[N + 5];
    inline void init()
    {
        for(register int i = 2;i <= N;++i)
            lg2[i] = lg2[i >> 1] + 1;
        rt[0] = 1,rt[1 << LG] = fpow(G,(mod - 1) >> LG + 2);
        for(register int i = LG;i;--i)
            rt[1 << i - 1] = (long long)rt[1 << i] * rt[1 << i] % mod;
        for(register int i = 1;i < N;++i)
            rt[i] = (long long)rt[i & i - 1] * rt[i & -i] % mod;
        fac[0] = 1;
        for(register int i = 1;i <= N;++i)
            fac[i] = (long long)fac[i - 1] * i % mod;
        ifac[N] = fpow(fac[N],mod - 2);
        for(register int i = N;i;--i)
            ifac[i - 1] = (long long)ifac[i] * i % mod;
        for(register int i = 1;i <= N;++i)
            inv[i] = (long long)ifac[i] * fac[i - 1] % mod;
    }
    struct poly
    {
        vector<int> a;
        inline poly() {}
        inline poly(int x)
        {
            a.push_back(x);
        }
        inline poly(int x,int y)
        {
            a.push_back(x),a.push_back(y);
        }
        inline poly(const vector<int> &o)
        {
            a = o;
        }
        inline poly(const poly &o)
        {
            a = o.a;
        }
        inline int size() const
        {
            return a.size();
        }
        inline void resize(int x)
        {
            a.resize(x);
        }
        inline int operator[](int x) const
        {
            if(x < 0 || x >= size())
                return 0;
            return a[x];
        }
        inline void clear()
        {
            vector<int>().swap(a);
        }
        inline poly rever() const
        {
            return poly(vector<int>(a.rbegin(),a.rend()));
        }
        inline void dif()
        {
            int n = size();
            for(register int i = 0,len = n >> 1;len;++i,len >>= 1)
                for(register int j = 0,*w = rt;j < n;j += len << 1,++w)
                    for(register int k = j,R;k < j + len;++k)
                        R = (long long)*w * a[k + len] % mod,
                        a[k + len] = dec(a[k],R),
                        a[k] = add(a[k],R);
        }
        inline void dit()
        {
            int n = size();
            for(register int i = 0,len = 1;len < n;++i,len <<= 1)
                for(register int j = 0,*w = rt;j < n;j += len << 1,++w)
                    for(register int k = j,R;k < j + len;++k)
                        R = add(a[k],a[k + len]),
                        a[k + len] = (long long)(a[k] - a[k + len] + mod) * *w % mod,
                        a[k] = R;
            reverse(a.begin() + 1,a.end());
            for(register int i = 0;i < n;++i)
                a[i] = (long long)a[i] * inv[n] % mod;
        }
        inline void ntt(int type = 1)
        {
            type == 1 ? dif() : dit();
        }
        friend inline poly operator+(const poly &a,const poly &b)
        {
            vector<int> ret(max(a.size(),b.size()));
            for(register int i = 0;i < ret.size();++i)
                ret[i] = add(a[i],b[i]);
            return poly(ret);
        }
        friend inline poly operator-(const poly &a,const poly &b)
        {
            vector<int> ret(max(a.size(),b.size()));
            for(register int i = 0;i < ret.size();++i)
                ret[i] = dec(a[i],b[i]);
            return poly(ret);
        }
        friend inline poly operator*(poly a,poly b)
        {
            if(a.a.empty() || b.a.empty())
                return poly();
            if(a.size() < 40 || b.size() < 40)
            {
                if(a.size() > b.size())
                    swap(a,b);
                poly ret;
                ret.resize(a.size() + b.size() - 1);
                for(register int i = 0;i < ret.size();++i)
                    for(register int j = 0;j <= i && j < a.size();++j)
                        ret.a[i] = (ret[i] + (long long)a[j] * b[i - j]) % mod;
                return ret;
            }
            int lim = 1,tot = a.size() + b.size() - 1;
            for(;lim < tot;lim <<= 1);
            a.resize(lim),b.resize(lim);
            a.ntt(),b.ntt();
            for(register int i = 0;i < lim;++i)
                a.a[i] = (long long)a[i] * b[i] % mod;
            a.ntt(-1),a.resize(tot);
            return a;
        }
        poly &operator+=(const poly &o)
        {
            resize(max(size(),o.size()));
            for(register int i = 0;i < o.size();++i)
                a[i] = add(a[i],o[i]);
            return *this;
        }
        poly &operator-=(const poly &o)
        {
            resize(max(size(),o.size()));
            for(register int i = 0;i < o.size();++i)
                a[i] = dec(a[i],o[i]);
            return *this;
        }
        poly &operator*=(poly o)
        {
            return (*this) = (*this) * o;
        }
        poly deriv() const
        {
            if(a.empty())
                return poly();
            vector<int> ret(size() - 1);
            for(register int i = 0;i < size() - 1;++i)
                ret[i] = (long long)(i + 1) * a[i + 1] % mod;
            return poly(ret);
        }
        poly integ() const
        {
            if(a.empty())
                return poly();
            vector<int> ret(size() + 1);
            for(register int i = 0;i < size();++i)
                ret[i + 1] = (long long)a[i] * inv[i + 1] % mod;
            return poly(ret);
        }
        inline poly modxn(int n) const
        {
            if(a.empty())
                return poly();
            n = min(n,size());
            return poly(vector<int>(a.begin(),a.begin() + n));
        }
        inline poly inver(int m) const
        {
            poly ret(fpow(a[0],mod - 2)),f,g;
            for(register int k = 1;k < m;)
            {
                k <<= 1,f.resize(k),g.resize(k);
                for(register int i = 0;i < k;++i)
                    f.a[i] = (*this)[i],g.a[i] = ret[i];
                f.ntt(),g.ntt();
                for(register int i = 0;i < k;++i)
                    f.a[i] = (long long)f[i] * g[i] % mod;
                f.ntt(-1);
                for(register int i = 0;i < (k >> 1);++i)
                    f.a[i] = 0;
                f.ntt();
                for(register int i = 0;i < k;++i)
                    f.a[i] = (long long)f[i] * g[i] % mod;
                f.ntt(-1);
                ret.resize(k);
                for(register int i = (k >> 1);i < k;++i)
                    ret.a[i] = dec(0,f[i]);
            }
            return ret.modxn(m);
        }
        inline pair<poly,poly> div(poly o) const
        {
            if(size() < o.size())
                return make_pair(poly(),*this);
            poly f,g;
            f = (rever().modxn(size() - o.size() + 1) * o.rever().inver(size() - o.size() + 1)).modxn(size() - o.size() + 1).rever();
            g = (modxn(o.size() - 1) - o.modxn(o.size() - 1) * f.modxn(o.size() - 1)).modxn(o.size() - 1);
            return make_pair(f,g);
        }
        inline poly log(int m) const
        {
            return (deriv() * inver(m)).integ().modxn(m);
        }
        inline poly exp(int m) const
        {
            poly ret(1),iv,it,d = deriv(),itd,itd0,t1;
            if(m < 70)
            {
                ret.resize(m);
                for(register int i = 1;i < m;++i)
                {
                    for(register int j = 1;j <= i;++j)
                        ret.a[i] = (ret[i] + (long long)j * operator[](j) % mod * ret[i - j]) % mod;
                    ret.a[i] = (long long)ret[i] * inv[i] % mod;
                }
                return ret;
            }
            for(register int k = 1;k < m;)
            {
                k <<= 1;
                it.resize(k >> 1);
                for(register int i = 0;i < (k >> 1);++i)
                    it.a[i] = ret[i];
                itd = it.deriv(),itd.resize(k >> 1);
                iv = ret.inver(k >> 1),iv.resize(k >> 1);
                it.ntt(),itd.ntt(),iv.ntt();
                for(register int i = 0;i < (k >> 1);++i)
                    it.a[i] = (long long)it[i] * iv[i] % mod,
                    itd.a[i] = (long long)itd[i] * iv[i] % mod;
                it.ntt(-1),itd.ntt(-1),it.a[0] = dec(it[0],1);
                for(register int i = 0;i < k - 1;++i)
                    itd.a[i % (k >> 1)] = dec(itd[i % (k >> 1)],d[i]);
                itd0.resize((k >> 1) - 1);
                for(register int i = 0;i < (k >> 1) - 1;++i)
                    itd0.a[i] = d[i];
                itd0 = (itd0 * it).modxn((k >> 1) - 1);
                t1.resize(k - 1);
                for(register int i = (k >> 1) - 1;i < k - 1;++i)
                    t1.a[i] = itd[(i + (k >> 1)) % (k >> 1)];
                for(register int i = k >> 1;i < k - 1;++i)
                    t1.a[i] = dec(t1[i],itd0[i - (k >> 1)]);
                t1 = t1.integ();
                for(register int i = 0;i < (k >> 1);++i)
                    t1.a[i] = t1[i + (k >> 1)];
                for(register int i = (k >> 1);i < k;++i)
                    t1.a[i] = 0;
                t1.resize(k >> 1),t1 = (t1 * ret).modxn(k >> 1),t1.resize(k);
                for(register int i = (k >> 1);i < k;++i)
                    t1.a[i] = t1[i - (k >> 1)];
                for(register int i = 0;i < (k >> 1);++i)
                    t1.a[i] = 0;
                ret -= t1;
            }
            return ret.modxn(m);
        }
        inline poly sqrt(int m) const
        {
            poly ret(1),f,g;
            for(register int k = 1;k < m;)
            {
                k <<= 1;
                f = ret,f.resize(k >> 1);
                f.ntt();
                for(register int i = 0;i < (k >> 1);++i)
                    f.a[i] = (long long)f[i] * f[i] % mod;
                f.ntt(-1);
                for(register int i = 0;i < k;++i)
                    f.a[i % (k >> 1)] = dec(f[i % (k >> 1)],(*this)[i]);
                g = (2 * ret).inver(k >> 1),f = (f * g).modxn(k >> 1),f.resize(k);
                for(register int i = (k >> 1);i < k;++i)
                    f.a[i] = f[i - (k >> 1)];
                for(register int i = 0;i < (k >> 1);++i)
                    f.a[i] = 0;
                ret -= f;
            }
            return ret.modxn(m);
        }
        inline poly pow(int m,int k1,int k2 = -1) const
        {
            if(a.empty())
                return poly();
            if(k2 == -1)
                k2 = k1;
            int t = 0;
            for(;t < size() && !a[t];++t);
            if((long long)t * k1 >= m)
                return poly();
            poly ret;
            ret.resize(m);
            int u = fpow(a[t],mod - 2),v = fpow(a[t],k2);
            for(register int i = 0;i < m - t * k1;++i)
                ret.a[i] = (long long)operator[](i + t) * u % mod;
            ret = ret.log(m - t * k1);
            for(register int i = 0;i < ret.size();++i)
                ret.a[i] = (long long)ret[i] * k1 % mod;
            ret = ret.exp(m - t * k1),t *= k1,ret.resize(m);
            for(register int i = m - 1;i >= t;--i)
                ret.a[i] = (long long)ret[i - t] * v % mod;
            for(register int i = 0;i < t;++i)
                ret.a[i] = 0;
            return ret;
        }
    };
}
using Poly::init;
using Poly::poly;
inline int C(int n,int m)
{
    return n < m ? 0 : (long long)Poly::fac[n] * Poly::ifac[m] % mod * Poly::ifac[n - m] % mod;
}
namespace Newton_Series
{
    poly seg[(N << 2) + 5];
    inline poly comp_xplusa(const poly &f,int a)
    {
        int n = f.size();
        poly t1,t2,ret;
        t1.resize(n),t2.resize(n);
        for(register int i = 0,pw = 1;i < n;++i,pw = (long long)pw * a % mod)
            t1.a[n - 1 - i] = (long long)Poly::fac[i] * f[i] % mod,
            t2.a[i] = (long long)Poly::ifac[i] * pw % mod;
        t1 *= t2;
        ret.resize(n);
        for(register int i = 0;i < n;++i)
            ret.a[i] = (long long)Poly::ifac[i] * t1[n - 1 - i] % mod;
        return ret;
    }
    void build(int n,int p)
    {
        if(n == 1)
        {
            seg[p].resize(2),seg[p].a[1] = 1;
            return ;
        }
        int mid = n + 1 >> 1;
        build(mid,ls),build(n - mid,rs);
        seg[p] = seg[ls] * comp_xplusa(seg[rs],dec(0,mid));
    }
    poly solve(const poly &f,int n,int p)
    {
        if(n == 1)
            return f;
        int mid = n + 1 >> 1;
        poly ret;
        ret.resize(n);
        pair<poly,poly> res = f.div(seg[ls]);
        poly t1 = solve(res.second,mid,ls);
        poly t2 = solve(comp_xplusa(res.first,mid),n - mid,rs);
        for(register int i = 0;i < mid;++i)
            ret.a[i] = t1[i];
        for(register int i = mid;i < n;++i)
            ret.a[i] = t2[i - mid];
        return ret;
    }
    poly calc(const poly &f)
    {
        int n = f.size();
        poly ret;
        build(n,1),ret = solve(f,n,1);
        for(register int i = 0;i < n;++i)
            ret.a[i] = (long long)ret[i] * Poly::fac[i] % mod;
        return ret;
    }
}
inline poly mulT(poly a,poly b,int k = -1)
{
    if(a.a.empty() || b.a.empty())
        return poly();
    int n = a.size(),m = b.size();
    if(k == -1)
        k = n - m + 1;
    if(k <= 0)
        return poly();
    if(k < 40 || b.size() < 40)
    {
        poly ret;
        ret.resize(k);
        for(register int i = 0;i < k;++i)
            for(register int j = 0;j < b.size();++j)
                ret.a[i] = (ret[i] + (long long)a[i + j] * b[j]) % mod;
        return ret;
    }
    reverse(b.a.begin(),b.a.end());
    if(k == n - m + 1)
    {
        int lim = 1;
        for(;lim < n;lim <<= 1);
        a.resize(lim),b.resize(lim);
        a.ntt(),b.ntt();
        for(register int i = 0;i < lim;++i)
            a.a[i] = (long long)a[i] * b[i] % mod;
        a.ntt(-1);
        for(register int i = 0;i < k;++i)
            a.a[i] = a[m - 1 + i];
        for(register int i = k;i < lim;++i)
            a.a[i] = 0;
        a.resize(n - m + 1);
    }
    else
    {
        a *= b,a.resize(n + m - 1);
        for(register int i = 0;i < k;++i)
            a.a[i] = a[m - 1 + i];
        for(register int i = k;i < n + m - 1;++i)
            a.a[i] = 0;
        a.resize(k);
    }
    return a;
}
poly f,seg[(N << 2) + 5][2][2];
int ans[N + 5];
void build(int p,int l,int r)
{
    if(l == r)
    {
        seg[p][1][0] = 1,seg[p][1][1] = 0,
        seg[p][0][0] = l < 2 ? !l : (long long)(l - 1) * Poly::inv[l] % mod;
        if(l < 2)
            seg[p][0][1] = 0;
        else
            seg[p][0][1] = poly(Poly::inv[l],Poly::inv[l]);
        return ;
    }
    int mid = l + r >> 1;
    build(ls,l,mid),build(rs,mid + 1,r);
    for(register int k = 0;k < 2;++k)
        for(register int i = 0;i < 2;++i)
            for(register int j = 0;j < 2;++j)
                seg[p][i][j] += seg[rs][i][k] * seg[ls][k][j];
}
void solve(int p,int l,int r,poly f[2])
{
    if(l == r)
    {
        f[0].size()==0?f[0].resize(1):(void)1;
        f[1].size()<=1?f[1].resize(2):(void)1;
        if(l < 2)
            ans[l] = l ? 0 : f[0][0];
        else
            ans[l] = ((long long)f[0][0] * (l - 1) + f[1][0] + f[1][1]) % mod * Poly::inv[l] % mod;
        return ;
    }
    int mid = l + r >> 1;
    poly g[2];
    g[0] = f[0].modxn((mid - l >> 1) + 2),
    g[1] = f[1].modxn((mid - l >> 1) + 2);
    solve(ls,l,mid,g);
    g[0] = mulT(f[0],seg[ls][0][0]) + mulT(f[1],seg[ls][0][1]),
    g[1] = mulT(f[0],seg[ls][1][0]) + mulT(f[1],seg[ls][1][1]);
    solve(rs,mid + 1,r,g);
}
int main()
{
    init();
    scanf("%d%d",&n,&k),f.resize(k);
    for(register int i = 0;i < k;++i)
        scanf("%d",&f.a[i]);
    f = Newton_Series::calc(f);
    build(1,0,n);
    poly g[2];
    g[0].resize((n + 1 >> 1) + 1),g[1].resize((n + 1 >> 1) + 1);
    for(register int i = 0;i < g[0].size();++i)
        g[0].a[i] = f[i];
    solve(1,0,n,g);
    for(register int i = 1;i <= n;++i)
        ans[i] = (long long)ans[i] * Poly::fac[i] % mod;
    for(register int i = 1;i <= n;++i)
        printf("%d%c",ans[i]," \n"[i == n]);
}