P5850 calc

P5850 calc

别做P4463，来做P5850，可以少写一个MTT

题面

题目链接

解法

首先想到能构造生成函数。

$$F(x) = \prod_{t=1}^k (1 + tx)$$

答案就是$n! [x^n]F(x)$，正确性显然。

$k$巨大（$10^9$级别），套路地将$\prod$化为$\sum$：

$$F(x) = \exp \sum_{t=1}^k \ln(1+tx)$$

使用麦克劳林级数将$\ln$化简：

$$F(x) = \exp \sum_{t=1}^k {\sum_{j=1}^{+\infty}(-1)^{j-1}\frac{t^jx^j}{j}} \\ \ln F(x) = \sum_{t=1}^k {\sum_{j=1}^{+\infty}(-1)^{j-1}\frac{t^jx^j}{j}} \\ [x^n] \ln F(x) = (-1)^{n-1} \frac{1}{n} \sum_{t=1}^k t^n$$

将$\ln F(x)$每项系数求出再进行多项式$\exp$得到$F(x)$。$\sum$左边的部分可以$O(1)$求出，右边是个自然数幂和的形式。令自然数幂和的指数生成函数为$S(x)$，则有

$$[x^n] \ln F(x) = (-1)^{n-1} \frac{n!}{n} S(n) = (-1)^{n-1} (n-1)! S(n)$$ $$S(x) = \sum_{j=0}^{+\infty} \frac{x^j}{j!} \sum_{i=1}^k i^j = \sum_{i=1}^k \sum_{j=0}^{+\infty} \frac{i^jx^j}{j!}$$

发现右边一大串其实是$e^{ix}$的泰勒展开式。

$$S(x) = \sum_{i=1}^k e^{ix}$$

而且是个等比数列求和。

$$S(x) = \frac{e^{(k+1)x}-e^x}{e^x-1} = \frac{\exp ((k+1)x) - \exp x}{\exp x - 1}$$

会发现分子和分母的常数项都是$0$，而常数项为$0$的多项式没有逆，所以分子分母同时乘上$\frac{1}{x}$即可。

AC代码

/**
 * @file:           P5850.cpp
 * @author:         yaoxi-std
 * @url:            https://www.luogu.com.cn/problem/P5850
*/
// #pragma GCC optimize ("O2")
// #pragma GCC optimize ("Ofast", "inline", "-ffast-math")
// #pragma GCC target ("avx,sse2,sse3,sse4,mmx")
#include <bits/stdc++.h>
using namespace std;
#define resetIO(x) \
    freopen(#x ".in", "r", stdin), freopen(#x ".out", "w", stdout)
#define debug(fmt, ...) \
    fprintf(stderr, "[%s:%d] " fmt "\n", __FILE__, __LINE__, ##__VA_ARGS__)
namespace fastio {
const int BUFFSIZ = 1 << 21;
char pbuf[BUFFSIZ], *pp = pbuf;
inline void flush() {
    fwrite(pbuf, 1, BUFFSIZ, stdout), pp = pbuf;
}
inline void pc(char ch) {
    (pp == pbuf + BUFFSIZ) ? flush() : void();
    *pp++ = ch;
}
};  // namespace fastio
using fastio::pc;
using fastio::flush;
template <class _Tp>
inline _Tp& read(_Tp& x) {
    bool sign = false;
    char ch = getchar();
    long double tmp = 1;
    for (; !isdigit(ch); ch = getchar())
        sign |= (ch == '-');
    for (x = 0; isdigit(ch); ch = getchar())
        x = x * 10 + (ch ^ 48);
    if (ch == '.')
        for (ch = getchar(); isdigit(ch); ch = getchar())
            tmp /= 10.0, x += tmp * (ch ^ 48);
    return sign ? (x = -x) : x;
}
template <class _Tp>
inline void write(_Tp x) {
    if (x < 0)
        pc('-'), x = -x;
    if (x > 9)
        write(x / 10);
    pc((x % 10) ^ 48);
}
const int MAXN = 1 << 20;
const int INF = 0x3f3f3f3f;
const int MOD = 998244353;
namespace polynomial {
int inv[MAXN];
inline int add(int x, int y) {
    x += y;
    return x >= MOD ? x - MOD : x;
}
inline int sub(int x, int y) {
    x -= y;
    return x < 0 ? x + MOD : x;
}
inline int qpow(int x, int y, int p = MOD) {
    int ret = 1;
    for (; y; y >>= 1, x = 1ll * x * x % p)
        if (y & 1)
            ret = 1ll * ret * x % p;
    return ret;
}
template <class _Tp>
void change(_Tp* f, int len) {
    static int rev[MAXN];
    for (int i = rev[0] = 0; i < len; ++i) {
        rev[i] = rev[i >> 1] >> 1;
        if (i & 1)
            rev[i] |= len >> 1;
    }
    for (int i = 0; i < len; ++i)
        if (i < rev[i])
            swap(f[i], f[rev[i]]);
}
void ntt(int* f, int len, int on) {
    change(f, len);
    for (int h = 2; h <= len; h <<= 1) {
        int gn = qpow(3, (MOD - 1) / h);
        for (int j = 0; j < len; j += h) {
            int g = 1;
            for (int k = j; k < j + h / 2; ++k) {
                int u = f[k], t = 1ll * g * f[k + h / 2] % MOD;
                f[k] = add(u, t), f[k + h / 2] = sub(u, t);
                g = 1ll * g * gn % MOD;
            }
        }
    }
    if (on == -1) {
        reverse(f + 1, f + len);
        int inv = qpow(len, MOD - 2);
        for (int i = 0; i < len; ++i)
            f[i] = 1ll * f[i] * inv % MOD;
    }
}
int polymul(const int* f, int n, const int* g, int m, int* ans) {
    static int tf[MAXN], tg[MAXN];
    int len = 1;
    while (len < n + m - 1)
        len <<= 1;
    copy(f, f + n, tf);
    fill(tf + n, tf + len, 0);
    copy(g, g + m, tg);
    fill(tg + m, tg + len, 0);
    ntt(tf, len, 1);
    ntt(tg, len, 1);
    for (int i = 0; i < len; ++i)
        tf[i] = 1ll * tf[i] * tg[i] % MOD;
    ntt(tf, len, -1);
    copy(tf, tf + n + m - 1, ans);
    return n + m - 1;
}
int polyinv(const int* f, int n, int* ans) {
    static int tmp[MAXN];
    int len = 1;
    while (len < n)
        len <<= 1;
    fill(ans, ans + len + len, 0);
    ans[0] = qpow(f[0], MOD - 2);
    for (int h = 2; h <= len; h <<= 1) {
        copy(f, f + h, tmp);
        fill(tmp + h, tmp + h + h, 0);
        ntt(tmp, h + h, 1);
        ntt(ans, h + h, 1);
        for (int i = 0; i < h + h; ++i)
            ans[i] = 1ll * ans[i] * (2 - 1ll * ans[i] * tmp[i] % MOD + MOD) % MOD;
        ntt(ans, h + h, -1);
        fill(ans + h, ans + h + h, 0);
    }
    return n;
}
int derivation(const int* f, int n, int* ans) {
    for (int i = 0; i < n - 1; ++i)
        ans[i] = 1ll * f[i + 1] * (i + 1) % MOD;
    return ans[n - 1] = 0, n - 1;
}
int integral(const int* f, int n, int* ans) {
    for (int i = n; i >= 1; --i)
        ans[i] = 1ll * f[i - 1] * inv[i] % MOD;
    return ans[0] = 0, n + 1;
}
int ln(const int* f, int n, int* ans) {
    static int tf[MAXN], tg[MAXN];
    derivation(f, n, tf);
    polyinv(f, n, tg);
    polymul(tf, n - 1, tg, n, ans);
    integral(ans, n - 1, ans);
    fill(ans + n, ans + n + n, 0);
    return n;
}
int exp(const int* f, int n, int* ans) {
    static int tmp[MAXN];
    ans[0] = 1, ans[1] = 0;
    for (int h = 2; h <= (n << 1); h <<= 1) {
        ln(ans, h, tmp);
        for (int i = 0; i < h; ++i)
            tmp[i] = add(i == 0, sub(f[i], tmp[i]));
        polymul(ans, h, tmp, h, ans);
    }
    return n;
}
};  // namespace polynomial
using namespace polynomial;
int k, m, ta[MAXN], tb[MAXN], tc[MAXN], ts[MAXN], fac[MAXN], ans[MAXN];
signed main() {
    read(k), read(m);
    fac[0] = inv[1] = 1;
    for (int i = 1; i <= m; ++i)
        fac[i] = 1ll * fac[i - 1] * i % MOD;
    for (int i = 2; i < MAXN; ++i)
        inv[i] = 1ll * (MOD - MOD / i) * inv[MOD % i] % MOD;
    tc[1] = (k + 1) % MOD;
    exp(tc, m + 2, ta);
    tc[1] = 1;
    exp(tc, m + 2, tb);
    for (int i = 0; i < m + 2; ++i)
        ta[i] = sub(ta[i], tb[i]);
    tb[0] = sub(tb[0], 1);
    for (int i = 0; i <= m; ++i)
        ta[i] = ta[i + 1], tb[i] = tb[i + 1];
    polyinv(tb, m + 1, tc);
    polymul(ta, m + 1, tc, m + 1, ts);
    ts[0] = 0;
    for (int i = 1; i <= m; ++i) {
        if (i & 1)
            ts[i] = add(0, 1ll * fac[i - 1] * ts[i] % MOD);
        else
            ts[i] = sub(0, 1ll * fac[i - 1] * ts[i] % MOD);
    }
    exp(ts, m + 1, ans);
    for (int i = 1; i <= m; ++i)
        write(1ll * ans[i] * fac[i] % MOD), pc('\n');
    flush();
    return 0;
}

subtask1过了，subtask2一直卡常，好像常数卡到现在的$\frac{1}{4}$倍，所以不卡了，代码贴出来作为参考。