Luogu 4705 玩游戏

看见这个题依稀想起了$5$月月赛时候的事情，到现在仍然它感觉非常神仙。

游戏$k$次价值的期望答案

$$ans_k = frac{1}{nm}sum_{i = 1}^{n}sum_{j = 1}^{n}(a_i + b_j)^k$$

二项式定理展开

$$ans_k=frac{1}{nm}sum_{i = 1}^{n}sum_{j = 1}^{m}sum_{t = 0}^{k}inom{k}{t}a_i^tb_j^{k - t}$$

$$ = frac{1}{nm}sum_{t = 0}^{k}inom{k}{t}sum_{i = 1}^{n}a_i^tsum_{j = 1}^{m}b_j^{k - t}$$

$$= frac{k!}{nm}frac{sum_{i = 1}^{n}a_i^t}{t!}frac{sum_{j = 1}^{m}b_j^{k - t}}{(k - t)!}$$

容易发现这后面是一个卷积的形式。

我们设$A(x) = frac{sum_{i = 1}^{n}a_i^x}{x!}$，$B(x) = frac{sum_{i = 1}^{n}b_i^x}{x!}$，答案就变成了

$$frac{k!}{nm}(A * B)(k)$$

现在考虑怎么算这个$sum_{i = 1}^{n}a_i^k$。

不会了，点开题解

我们构造一个多项式$F(x) = prod_{i = 1}^{n}(1 - a_ix)$，这个$F(x)$可以通过分治算出来，时间复杂度$T(n) = 2T(frac{n}{2}) + O(nlogn)$趋向$O(nlogn)$？？？。

接下来要开始变魔术了，

设$G(x) = lnF(x)$，

$$G(x) = sum_{i = 1}^{n}ln(1 - a_ix)$$

对里面的$ln$在$x_0 = 0$处泰勒展开。

我们知道

$$ln(1 - x) = 0 - frac{x}{1} - frac{x^2}{2} - frac{x^3}{3} - cdots = -sum_{i = 1}^{infty}frac{x^i}{i}$$

所以

$$G(x) = sum_{i = 1}^{n}sum_{j = 1}^{k}-frac{a_i^j}{j}x^j$$

因为这个题对$x^{k + 1}$次取模了，所以后面的项可以不用再写了。

变形一下

$$G(x) = sum_{j = 1}^{k}(-frac{1}{j}sum_{i = 1}^{n}a_i^j)x^j$$

发现$1$到$k$的$sum_{i = 1}^{n}a_i^k$直接算出来了，而$i = 0$时候这东西显然为$n$。

鼓掌~~~

现在回过头来考虑一下这个魔术是怎么变出来的。

考虑构造每一个$a_i$的生成函数

$$1 + a_ix + a_i^2x^2 + a_i^3x^3 + cdots$$

把每一个$a_i$的生成函数都构造出来然后加起来的第$i$项系数就是$k = i$时候的答案了。

$$F(x) = sum_{i = 1}^{n}frac{1}{1 - a_ix}$$

发现这个$F(x)$并不是很好算，而

$$ln'(1 - a_ix) = frac{1}{1 - a_ix}$$

$$(ln(1 - a_ix))' = frac{-a_i}{1 - a_ix}$$

所以先计算

$$G(x) = sum_{i = 1}^{n}frac{-a_i}{1 - a_ix} = (lnprod_{i = 1}^{n}(1 - a_ix))'$$

把这两个式子写回到生成函数的形式，发现

$$F(x) = -x * G(x) + n$$

于是大功告成。

用$vector$写多项式真爽。

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
using namespace std;
typedef long long ll;
typedef vector <ll> poly;

const int N = 1e5 + 5;
const int Maxn = 1e5;

int n, m, K;
ll a[N], b[N], fac[N], ifac[N], inv[N];

template <typename T>
inline void read(T &X) {
    X = 0; char ch = 0; T op = 1;
    for (; ch > '9' || ch < '0'; ch = getchar())
        if (ch == '-') op = -1;
    for (; ch >= '0' && ch <= '9'; ch = getchar())
        X = (X << 3) + (X << 1) + ch - 48;
    X *= op;
}

inline void deb(poly c) {
    for (int i = 0; i < (int)c.size(); i++)
        printf("%lld%c", c[i], " 
"[i == (int)c.size() - 1]);
}

namespace Poly {
    const int L = 1 << 20;
    const ll P = 998244353LL;

    int lim, pos[L];

    inline ll fpow(ll x, ll y) {
        ll res = 1;
        for (; y > 0; y >>= 1) {
            if (y & 1) res = res * x % P;
            x = x * x % P;
        }
        return res;
    }

    template <typename T>
    inline void inc(T &x, T y) {
        x += y;
        if (x >= P) x -= P;
    }

    template <typename T>
    inline void sub(T &x, T y) {
        x -= y;
        if (x < 0) x += P;
    }

    inline void prework(int len) {
        int l = 0;
        for (lim = 1; lim < len; lim <<= 1, ++l);
        for (int i = 0; i < lim; i++)
            pos[i] = (pos[i >> 1] >> 1) | ((i & 1) << (l - 1));
    }

    inline void ntt(poly &c, int opt) {
        c.resize(lim, 0);
        for (int i = 0; i < lim; i++)
            if (i < pos[i]) swap(c[i], c[pos[i]]);
        for (int i = 1; i < lim; i <<= 1) {
            ll wn = fpow(3, (P - 1) / (i << 1));
            if (opt == -1) wn = fpow(wn, P - 2);
            for (int len = i << 1, j = 0; j < lim; j += len) {
                ll w = 1;
                for (int k = 0; k < i; k++, w = w * wn % P) {
                    ll x = c[j + k], y = w * c[j + k + i] % P;
                    c[j + k] = (x + y) % P, c[j + k + i] = (x - y + P) % P;
                }
            }
        }

        if (opt == -1) {
            ll inv = fpow(lim, P - 2);
            for (int i = 0; i < lim; i++) c[i] = c[i] * inv % P;
        }
    }

    inline poly operator * (const poly &x, const poly &y) {
        poly res, u = x, v = y;
        prework(u.size() + v.size() - 1);
        ntt(u, 1), ntt(v, 1);
        for (int i = 0; i < lim; i++) res.push_back(v[i] * u[i] % P);
        ntt(res, -1);
        res.resize(u.size() + v.size() - 1);
        return res;
    }

    poly getInv(poly x, int len) {
        x.resize(len);
        if (len == 1) {
            poly res;
            res.push_back(x[0]);
            return res;
        }
        poly y = getInv(x, (len + 1) >> 1);
        prework(len << 1);

        poly u = x, v = y, res;
        ntt(u, 1), ntt(v, 1);
        for (int i = 0; i < lim; i++) res.push_back(v[i] * (2LL - u[i] * v[i] % P + P) % P);
        ntt(res, -1);

        res.resize(len);
        return res;
    }

    inline void direv(poly &c) {
        for (int i = 0; i < (int)c.size() - 1; i++)
            c[i] = c[i + 1] * (i + 1) % P;
        c[c.size() - 1] = 0;
    }

    inline void integ(poly &c) {
        for (int i = (int)c.size() - 1; i > 0; i--)
            c[i] = c[i - 1] * inv[i] % P;
        c[0] = 0;
    }

    inline poly getLn(poly c) {
        poly a = getInv(c, (int)c.size());
        poly b = c;
        direv(b);

        poly res = b * a;
        res.resize(c.size());
        integ(res);
        return res;
    }

    inline poly solve(int l, int r, ll *c) {
        if (l == r) {
            poly res;
            res.push_back(1), res.push_back((P - c[l]) % P);
            return res;
        }

        int mid = ((l + r) >> 1);
        poly u = solve(l, mid, c), v = solve(mid + 1, r, c);
        poly res = u * v;
        res.resize(u.size() + v.size() - 1);
        return res;
    }

    inline poly mul(const poly &x, const poly &y) {
        return x * y;
    }

}

using Poly::P;
using Poly::fpow;
using Poly::inc;
using Poly::sub;
using Poly::solve;
using Poly::getLn;

inline void prework() {
    fac[0] = inv[0] = inv[1] = 1;
    for (int i = 1; i <= Maxn; i++) {
        fac[i] = fac[i - 1] * i % P;
        if (i > 1) inv[i] = (P - P / i) * inv[P % i] % P;
    }
    ifac[Maxn] = fpow(fac[Maxn], P - 2);
    for (int i = Maxn - 1; i >= 0; i--) ifac[i] = ifac[i + 1] * (i + 1) % P;
}

int main() {
/*    #ifndef ONLINE_JUDGE
        freopen("Sample.txt", "r", stdin);
    #endif    */
    
    freopen("input.txt", "r", stdin);
    freopen("my.out", "w", stdout);

    prework();

    read(n), read(m);
    for (int i = 1; i <= n; i++) read(a[i]);
    for (int i = 1; i <= m; i++) read(b[i]);
    read(K);
    ++K;

    poly Fa = solve(1, n, a), Fb = solve(1, m, b);

//    deb(Fa), deb(Fb);

    Fa.resize(K, 0), Fb.resize(K, 0);
    poly Ga = getLn(Fa), Gb = getLn(Fb);

    Ga.resize(K, 0), Gb.resize(K, 0);

//    deb(Ga), deb(Gb);

    poly A, B, C;

    A.push_back(n), B.push_back(m);
    for (int i = 1; i < K; i++) {
        A.push_back((P - Ga[i]) % P * i % P * ifac[i] % P);
        B.push_back((P - Gb[i]) % P * i % P * ifac[i] % P);
    }

//    deb(A), deb(B);

    C = Poly::mul(A, B);
    C.resize(K);

//    deb(C);

//    printf("
");
    ll invn = fpow(n, P - 2), invm = fpow(m, P - 2);
    for (int i = 1; i < K; i++)
        printf("%lld
", invn * invm % P * fac[i] % P * C[i] % P);

    return 0;
}