Codeforces 923E Perpetual Subtraction

设 P′i 表示一次操作后的 P(x=i)。
再设 F,F′ 分别为 P,P′ 的生成函数。

注意到 P′i=n∑j=iPjj+1。
故 F′(x)=n∑i=0P′ixi=n∑i=0xin∑j=iPjj+1=n∑j=0Pjj+1j∑i=0xj=n∑j=0Pjj+1⋅xj+1−1x−1=1x−1n∑j=0Pjxj+1−Pjj+1=1x−1∫x1(n∑j=0Pjtj)dt=1x−1∫x1F(t)dt

然而这要怎么算啊……
考虑设 G(x)=F(x+1),G′(x)=F′(x+1)。
于是 G′(x)=F′(x+1)=1x∫x0G(t)dt

若 gi=[xi]G(x),g′i=[xi]G′(x)，则 g′i=gii+1。
做 m 次的话就是 gi(i+1)m。

但是如何在 P,g 间转化呢？
首先考虑从 P 求 g。
有 G(x)=F(x+1)=n∑i=0Pi(x+1)i=n∑i=0Pii∑j=0(ij)xj=n∑j=0xjn∑i=jPi(ij)

即 gi=n∑j=iPj(ji)。
再推一推 gi=n∑j=iPj(ji)=n∑j=ij!Pji!(j−i)!i!gi=n∑j=ij!Pj(j−i)!=n−i∑j=0(i+j)!Pi+jj!

翻转一下就是卷积了。
从 g′ 转回 P′ 类似。

代码：

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <utility>
#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)
using namespace std;
const int N = 1 << 18;
const int mod = 998244353;
const int G = 3;
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;
}
struct poly
{
    int a[N + 5];
    inline const int &operator[](int x) const
    {
        return a[x];
    }
    inline int &operator[](int x)
    {
        return a[x];
    }
    inline void clear(int x = 0)
    {
        memset(a + x,0,(N - x + 1) << 2);
    }
} f,g;
int len,k,n,lg2[N + 5];
long long m;
int rev[N + 5],fac[N + 5],ifac[N + 5],inv[N + 5];
int rt[N + 5],irt[N + 5];
inline void init(int len)
{
    for(n = 1;n < len;n <<= 1);
    for(register int i = 2;i <= n;++i)
        lg2[i] = lg2[i >> 1] + 1;
    int w = fpow(G,(mod - 1) / n);
    rt[n >> 1] = 1;
    for(register int i = (n >> 1) + 1;i <= n;++i)
        rt[i] = (long long)rt[i - 1] * w % mod;
    for(register int i = (n >> 1) - 1;i;--i)
        rt[i] = rt[i << 1];
    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;
}
inline void ntt(poly &a,int type,int n)
{
    type == -1 && (reverse(a.a + 1,a.a + n),1);
    int lg = lg2[n] - 1;
    for(register int i = 0;i < n;++i)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << lg),
        i < rev[i] && (swap(a[i],a[rev[i]]),1);
    for(register int w = 2,m = 1;w <= n;w <<= 1,m <<= 1)
        for(register int i = 0;i < n;i += w)
            for(register int j = 0;j < m;++j)
            {
                int t = (long long)rt[m | j] * a[i | j | m] % mod;
                a[i | j | m] = dec(a[i | j],t),a[i | j] = add(a[i | j],t);
            }
    if(type == -1)
        for(register int i = 0;i < n;++i)
            a[i] = (long long)a[i] * inv[n] % mod;
}
inline void mul(poly &a,const poly &b,int n)
{
    static poly x,y;
    int lim = 1;
    x.clear(),y.clear();
    for(;lim < (n << 1);lim <<= 1);
    x = a,y = b;
    x.clear(n),y.clear(n);
    ntt(x,1,lim),ntt(y,1,lim);
    for(register int i = 0;i < lim;++i)
        x[i] = (long long)x[i] * y[i] % mod;
    ntt(x,-1,lim);
    x.clear(n),a = x;
}
inline poly inverse(const poly &f,int n)
{
    static int s[30];
    static poly g,h,q;
    int lim = 1,top = 0;
    g.clear();
    for(;n > 1;s[++top] = n,n = (n + 1) >> 1);
    g[0] = fpow(f[0],mod - 2);
    for(;top;--top)
    {
        n = s[top];
        for(;lim < (n << 1);lim <<= 1);
        q = g,h = f,h.clear(n);
        ntt(g,1,lim),ntt(h,1,lim);
        for(register int i = 0;i < lim;++i)
            g[i] = (long long)g[i] * g[i] % mod * h[i] % mod;
        ntt(g,-1,lim);
        for(register int i = 0;i < n;++i)
            g[i] = dec(add(q[i],q[i]),g[i]);
        g.clear(n);
    }
    return g;
}
inline void derivative(poly &f,int n)
{
    for(register int i = 1;i < n;++i)
        f[i - 1] = (long long)f[i] * i % mod;
    f[n - 1] = 0;
}
inline void integral(poly &f,int n)
{
    for(register int i = n - 1;~i;--i)
        f[i + 1] = (long long)f[i] * inv[i + 1] % mod;
    f[0] = 0;
}
inline poly ln(const poly &f,int n)
{
    static poly g;
    g = f,derivative(g,n),mul(g,inverse(f,n),n),integral(g,n);
    return g;
}
inline poly exp(const poly &f,int n)
{
    static int s[30];
    static poly g,h;
    int lim = 1,top = 0;
    g.clear();
    for(;n > 1;s[++top] = n,n = (n + 1) >> 1);
    g[0] = 1;
    for(;top;--top)
    {
        n = s[top];
        for(;lim < (n << 1);lim <<= 1);
        h = g,g = ln(g,n);
        for(register int i = 0;i < n;++i)
            g[i] = dec(f[i],g[i]);
        g[0] = add(g[0],1);
        ntt(g,1,lim),ntt(h,1,lim);
        for(register int i = 0;i < lim;++i) 
            g[i] = (long long)g[i] * h[i] % mod;
        ntt(g,-1,lim);
        g.clear(n);
    }
    return g;
}
inline poly power(const poly &f,int k,int n)
{
    static poly g;
    g = ln(f,n);
    for(register int i = 0;i < n;++i)
        g[i] = (long long)g[i] * k % mod;
    g = exp(g,n);
    return g;
}
namespace Mod_sqrt
{
    typedef pair<int,int> cp;
    int w;
    inline cp operator*(const cp &a,const cp &b)
    {
        return cp(((long long)a.first * b.first % mod + (long long)a.second * b.second % mod * w % mod) % mod,((long long)a.first * b.second % mod + (long long)a.second * b.first % mod) % mod);
    }
    inline cp pow(cp a,int b)
    {
        cp ret(1,0);
        for(;b;b >>= 1)
            (b & 1) && (ret = ret * a,1),a = a * a;
        return ret;
    }
    inline int mod_sqrt(int x)
    {
        int y = rand() % mod;
        for(;fpow(w = ((long long)y * y % mod - x + mod) % mod,mod - 1 >> 1) <= 1;y = rand() % mod);
        cp ret = pow(cp(y,1),mod + 1 >> 1);
        return min(ret.first,mod - ret.first);
    }
}
using Mod_sqrt::mod_sqrt;
inline poly sqrt(const poly &f,int n)
{
    static int s[30];
    static poly g,h;
    int top = 0;
    g.clear();
    for(;n > 1;s[++top] = n,n = (n + 1) >> 1);
    g[0] = mod_sqrt(f[0]);
    for(;top;--top)
    {
        n = s[top];
        for(register int i = 0;i < n;++i)
            h[i] = add(g[i],g[i]);
        h = inverse(h,n),mul(g,g,n);
        for(register int i = 0;i < n;++i)
            g[i] = add(g[i],f[i]);
        mul(g,h,n);
    }
    return g;
}
int main()
{
    scanf("%d%lld",&len,&m),m %= mod - 1,init(++len << 1);
    for(register int i = 0;i < len;++i)
        scanf("%d",f.a + len - i - 1),f[len - i - 1] = (long long)f[len - i - 1] * fac[i] % mod,g[i] = ifac[i];
    mul(g,f,len);
    for(register int i = 0;i < len;++i)
        g[len - i - 1] = (long long)g[len - i - 1] * fpow(fpow(i + 1,m),mod - 2) % mod;
    for(register int i = 0;i < len;++i)
        f[i] = (long long)ifac[i] * (i & 1 ? mod - 1 : 1) % mod;
    f.clear(len);
    mul(f,g,len);
    for(register int i = 0;i < len;++i)
        printf("%d ",(int)((long long)f[len - i - 1] * ifac[i] % mod));
}