[LOJ6569] 仙人掌计数

Statement

带标号仙人掌计数问题.（\(n \le 131071\)）

Solution

设仙人掌个数的生成函数为\(C(x)\)

• 对于与根相邻的块, 还是仙人掌, 生成函数为\(C(x)\)

• 包含根的环, 生成函数为\(\sum_{i \ge 2}\frac{C(x)^i}{2}\)

组合起来:

\[C(x) = x \exp{\frac{2C(x)-C(x)^2}{2-2C(x)}} \]

设\(G(C(x)) = x\exp{\frac{2C(x)-C(x)^2}{2-2C(x)}}-C(x)\), 那么:

\[\small{ \begin{aligned} G'(C(x)) &= x\left(\exp{\frac{2C(x)-C(x)^2}{2-2C(x)}}\right)'-1 \\ &= x \exp\left(\frac{2C(x)-C(x)^2}{2-2C(x)}\right)\left(\frac{2C(x)-C(x)^2}{2-2C(x)}\right)' - 1 \\ &= x \exp\left(\frac{2C(x)-C(x)^2}{2-2C(x)}\right) \left(\frac{\left(2-2C(x)\right)^2-\left(2C(x) - C(x)^2\right)(-2)}{(2-2C(x))^2}\right) - 1\\ &= x \exp\left(\frac{2C(x)-C(x)^2}{2-2C(x)}\right) \left(1+\frac{4C(x) - 2C(x)^2}{(2-2C(x))^2}\right) - 1 \end{aligned} } \]

牛顿迭代:

\[\begin{aligned} C_1(x) &= C(x) - \frac{G(C(x))}{G'(C(x))} \\ &= C(x) - \frac{2x\exp\left(\frac{2C(x)-C(x)^2}{2-2C(x)}\right)-2C(x)} {x \exp\left(\frac{2C(x)-C(x)^2}{2-2C(x)}\right) \left(1+\frac{1}{(C(x)-1)^2}\right) - 2} \end{aligned} \]

Code

/************************************************
 * Au: Hany01
 * Prob: loj161
 * Email: hany01dxx@gmail.com & hany01@foxmail.com
 * Inst: Yali High School
************************************************/

#include <bits/stdc++.h>

using namespace std;

typedef long long LL;
typedef long double LD;
typedef pair<int, int> PII;
#define Rep(i, j) for (register int i = 0, i##_end_ = (j); i < i##_end_; ++ i)
#define For(i, j, k) for (register int i = (j), i##_end_ = (k); i <= i##_end_; ++ i)
#define Fordown(i, j, k) for (register int i = (j), i##_end_ = (k); i >= i##_end_; -- i)
#define Set(a, b) memset(a, b, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define X first
#define Y second
#define PB(a) push_back(a)
#define MP(a, b) make_pair(a, b)
#define SZ(a) ((int)(a).size())
#define ALL(a) a.begin(), a.end()
#define INF (0x3f3f3f3f)
#define INF1 (2139062143)
#define debug(...) fprintf(stderr, __VA_ARGS__)
#define y1 wozenmezhemecaia

template <typename T> inline bool chkmax(T &a, T b) {
    return a < b ? a = b, 1 : 0;
}
template <typename T> inline bool chkmin(T &a, T b) {
    return b < a ? a = b, 1 : 0;
}

template <typename T> inline T read() {
    static T _, __;
    static char c_;

    for (_ = 0, __ = 1, c_ = getchar(); c_ < '0' || c_ > '9'; c_ = getchar())
        if (c_ == '-')
            __ = -1;

    for (; c_ >= '0' && c_ <= '9'; c_ = getchar())
        _ = (_ << 1) + (_ << 3) + (c_ ^ 48);

    return _ * __;
}
//EOT


const int MAXN = 1 << 19, MOD = 998244353, g0 = 3;
int ig0;
int pw[MAXN], pw_[MAXN];
int fac[MAXN], ifac[MAXN];

inline int fpm(int a, int b = MOD - 2) {
    register int ans = 1;

    for (; b; b >>= 1, a = (LL)a * a % MOD)
        if (b & 1)
            ans = (LL)ans * a % MOD;

    return ans;
}
inline int ad(int x, int y) {
    return (x += y) >= MOD ? x - MOD : x;
}
inline void inc(int &x, int y) {
    if ((x += y) >= MOD)
        x -= MOD;
}
inline int times2(int x) {
    return (x += x) >= MOD ? x - MOD : x;
}

int Init(int n) {
    int pt, N;

    for (pt = 0, N = 1; N <= n; N <<= 1, ++ pt);

    ig0 = fpm(g0, MOD - 2);
    For(i, 1, pt + 1)
    pw[1 << i] = fpm(g0, (MOD - 1) / (1 << i)),
                 pw_[1 << i] = fpm(ig0, (MOD - 1) / (1 << i));
    fac[0] = 1;
    For(i, 1, N - 1) fac[i] = (LL)fac[i - 1] * i % MOD;
    ifac[N - 1] = fpm(fac[N - 1]);
    Fordown(i, N - 1, 1) ifac[i - 1] = (LL)ifac[i] * i % MOD;
    return N;
}

inline void NTT(int *a, int n, int ty) {
    static int rev[MAXN];
    static int W[MAXN];
    register int pt = __builtin_ctz(n);
    Rep(i, n)

    if (i < (rev[i] = ((rev[i >> 1] >> 1) | ((i & 1) << (pt - 1)))))
        swap(a[i], a[rev[i]]);

    for (register int i = 2, i2 = 1; i <= n; i2 = i, i <<= 1) {
        W[0] = 1, W[1] = ty > 0 ? pw[i] : pw_[i];
        For(j, 2, i2 - 1) W[j] = (LL)W[j - 1] * W[1] % MOD;

        for (register int j = 0; j < n; j += i) {
            Rep(k, i2) {
                register int x = a[j + k], y = (LL)a[j + k + i2] * W[k] % MOD;
                a[j + k] = ad(x, y), a[j + k + i2] = ad(x, MOD - y);
            }
        }
    }

    if (ty < 1) {
        register int inv = fpm(n);
        Rep(i, n) a[i] = (LL)a[i] * inv % MOD;
    }
}

inline void Mult(int *f, int *g, int n, int *h) {
    static int f_[MAXN], g_[MAXN];
    Rep(i, n) f_[i] = f[i], g_[i] = g[i];
    For(i, n, n * 2 - 1) f_[i] = g_[i] = 0;
    NTT(f_, n << 1, 1), NTT(g_, n << 1, 1);
    Rep(i, n << 1) h[i] = (LL)f_[i] * g_[i] % MOD;
    NTT(h, n << 1, -1);
}

inline void Mult(int *f1, int *f2, int *f3, int n, int *h) {
    static int f1_[MAXN], f2_[MAXN], f3_[MAXN];
    Rep(i, n) f1_[i] = f1[i], f2_[i] = f2[i], f3_[i] = f3[i];
    For(i, n, n * 2 - 1) f1_[i] = f2_[i] = f3_[i] = 0;
    NTT(f1_, n << 1, 1), NTT(f2_, n << 1, 1), NTT(f3_, n << 1, 1);
    Rep(i, n << 1) h[i] = (LL)f1_[i] * f2_[i] % MOD * f3_[i] % MOD;
    NTT(h, n << 1, -1);
}

namespace Inv {
static int f[MAXN];
inline void Inv_(int *g, int n) {
    static int h[MAXN];

    if (n == 1) {
        g[0] = fpm(f[0]);
        return;
    }

    Inv_(g, n >> 1);
    Mult(g, g, f, n, h);
    Rep(i, n) g[i] = ad(ad(g[i], g[i]), MOD - h[i]);
}
inline void Inv(int *A, int n, int *ans) {
    Rep(i, n) f[i] = A[i], ans[i] = 0;
    Inv_(ans, n);
}
}

inline void Int(int *f, int n, int *g) {
    Fordown(i, n - 1, 1) g[i] = (LL)f[i - 1] * fpm(i) % MOD;
    g[0] = 0;
}
inline void Der(int *f, int n, int *g) {
    For(i, 1, n - 1) g[i - 1] = (LL)f[i] * i % MOD;
    g[n - 1] = 0;
}

inline void Ln(int *f, int n, int *g) {
    static int h[MAXN];
    Der(f, n, h), Inv:: Inv(f, n, g);
    Mult(h, g, n, g), Int(g, n, g);
}

namespace Exp {
static int G[MAXN];
inline void Exp_(int *F, int n) {
    static int H[MAXN];

    if (n == 1) {
        F[0] = 1;
        return;
    }

    Exp_(F, n >> 1);
    Ln(F, n, H);
    Rep(i, n) H[i] = ad(G[i], MOD - H[i]);
    H[0] = ad(H[0], 1);
    Mult(H, F, n, F);
}
inline void Exp(int *g, int n, int *ans) {
    Rep(i, n) G[i] = g[i], ans[i] = 0;
    Exp_(ans, n);
}
}

inline void Pow(int *f, int n, int k, int *g) {
    static int h[MAXN];
    Ln(f, n, h);
    Rep(i, n) h[i] = (LL)h[i] * k % MOD;
    Exp:: Exp(h, n, g);
}

namespace Sqrt {
static int A[MAXN], B[MAXN], a[MAXN];
void Sqrt_(int *b, int n) {
    if (n == 1) {
        b[0] = sqrt(a[0]);
        return;
    }

    Sqrt_(b, n >> 1);
    Rep(i, n) A[i] = b[i];
    Mult(A, A, n, A);
    Rep(i, n) A[i] = ad(A[i], a[i]), B[i] = ad(b[i], b[i]);
    Inv:: Inv(B, n, B);
    Mult(A, B, n, b);
}
void Sqrt(int *x, int *y, int n) {
    Rep(i, n) a[i] = x[i], y[i] = 0;
    Sqrt_(y, n);
}
}

void Newton(int *C, int n) {
    if (n == 2) {
        C[1] = 1;
        return;
    }

    Newton(C, n >> 1);
    static int F1[MAXN], F2[MAXN], F[MAXN], G[MAXN];
    Mult(C, C, n, F1);
    Rep(i, n) F1[i] = ad(times2(C[i]), MOD - F1[i]);
    Rep(i, n) F2[i] = ad(0, MOD - times2(C[i]));
    inc(F2[0], 2);
    Inv:: Inv(F2, n, F2);
    Mult(F1, F2, n, F1);
    Exp:: Exp(F1, n, F1);
    Fordown(i, n - 1, 1) F1[i] = F1[i - 1];
    F1[0] = 0;
    Rep(i, n) F[i] = ad(times2(F1[i]), MOD - times2(C[i]));
    Rep(i, n) F2[i] = C[i];
    inc(F2[0], MOD - 1);
    Mult(F2, F2, n, F2);
    Inv:: Inv(F2, n, F2);
    inc(F2[0], 1);
    Mult(F1, F2, n, G);
    inc(G[0], MOD - 2);
    Inv:: Inv(G, n, G);
    Mult(F, G, n, F);
    Rep(i, n) C[i] = ad(C[i], MOD - F[i]);
}

int main() {
#ifdef hany01
    freopen("loj161.in", "r", stdin);
    freopen("loj161.out", "w", stdout);
#endif

    static int C[MAXN];
    int N = Init(131071);

    Newton(C, N);

    for (int Q = read<int>(), x; Q --;) {
        x = read<int>();
        printf("%lld\n", (LL)C[x] * fac[x - 1] % MOD);
    }

    return 0;
}