Dyd's Blog

He who has a strong enough why can bear almost any how.

luoguP5667 拉格朗日插值2

Posted on In
Symbols count in article: 3.4k Reading time ≈ 3 mins.

拉插 + 多项式卷积，难点在平移多项式

拉格朗日插值2

思路

直接线性插值是 $O(n^2)$ 的复杂度，考虑到询问比较特殊，可以进一步优化，看式子
$$
\begin{aligned}
f(m + i) = & \sum_{j = 0}^n y_j \prod_{k \ne j}\frac{(m + i - k)}{(j - k)} \\
= & \frac{(m + i)!}{(m - n + i - 1)!} \sum_{j = 0}^n \frac{f(j) * (-1)^{n - j}}{j! (n - j)!} \times \frac{1}{m + i - j}
\end{aligned}
$$
仔细看，后面的东西是个卷积（ $j$ 和 $i - j$ 贡献给 $i$ ），但问题是循环条件是 $j \in [0, n]$ 而不是 $j \in [0, i]$ ，考虑平移多项式

令 $g_i = \frac{f(i) * (-1)^{n - i}}{i! (n - i)!}, h_i = \frac{1}{m - n + i}$ ，注意，本来 $h_i$ 应该是 $\frac{1}{m + i}$ 的，但我们让它 $-n$ ，即是 $h_{x - n}$ ，那么相当于整个多项式向右平移 $n$ 位，原来 $i$ 的值现在应该在 $i + n$ 处，那么我们把 $G, H$ 乘起来得 $L$ ，考虑：
$$
\begin{aligned}
& l_x = \sum_{i = 0}^{x} g_i * h_{x - i} \\
& l_{x + n} = \sum_{i = 0}^{x + n} g_i * h_{x + n - i} \\
& l_{x + n} = \sum_{i = 0}^{x + n} \frac{f(i) * (-1)^{n - i}}{i! (n - i)!} * \frac{1}{m + x - i}
\end{aligned}
$$
不难发现当 $i > n$ 的多算了，所以我们对于 $i > n$ 让 $g_i = 0$ ，那么后面那个东西就是 $l_{x + n}$

现在问题是 $m \le 10^8$ ，我们要求 $(m - n + i)!$ 和其逆元，好像 $O(m)$ 都有点危险，发现其实很大一部分消去了，考虑预处理出 $mfac[i] = \prod_{j = 0}^{i} (m - n + j)$ 和其逆元，只算没消的部分

代码

注意 $0$ 处， $m - n + i - 1$ 的下标到了 $-1$ 会越界，特殊处理

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
#include <cstdio>
#include <vector>
using LL = long long;
using poly = std::vector<int>;
const int N = 1.6e5 + 100, P = 998244353, G = 3, iG = 332748118;
int n, m, f[N];
int rev[1 << 20 | 100];
int ifac[N], miv[N << 1], mfac[N << 1], imfac[N << 1];
void adj(int &x){ x += (x >> 31) & P; }
int qpow(int x, int y = P - 2)
{
    int res = 1;
    for (; y; y >>= 1, x = LL(x) * x % P) if (y & 1) res = LL(res) * x % P;
    return res;
}
void init()
{
    int fac = 1;
    for (int i = 2; i <= n; ++i) fac = LL(i) * fac % P;
    ifac[n] = qpow(fac);
    for (int i = n - 1; i > 1; --i) ifac[i] = LL(i + 1) * ifac[i + 1] % P;
    ifac[0] = ifac[1] = 1;
    mfac[0] = m - n;
    for (int i = 1; i <= (n << 1); ++i) mfac[i] = LL(m - n + i) * mfac[i - 1] % P;
    imfac[n << 1] = qpow(mfac[n << 1]);
    for (int i = (n << 1) - 1; ~i; --i) imfac[i] = LL(m - n + i + 1) * imfac[i + 1] % P;
    miv[0] = imfac[0];
    for (int i = 1; i <= (n << 1); ++i) miv[i] = LL(imfac[i]) * mfac[i - 1] % P;
}
void rdy(int &bit, int &tot, int len){ for (tot = 1, bit = 0; tot <= len; tot <<= 1, ++bit); }
void get_r(int bit, int tot){ for (int i = 1; i < tot; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bit - 1)); }
void ntt(int *x, int tot, int op)
{
    static int gn[1 << 20 | 100];
    for (int i = 1; i < tot; ++i) if (i < rev[i]) std::swap(x[i], x[rev[i]]);
    for (int mid = 1; mid < tot; mid <<= 1)
    {
        int g0 = qpow(op == 1 ? G : iG, (P - 1) / (mid << 1));
        gn[0] = 1;
        for (int i = 1; i < mid; ++i) gn[i] = LL(g0) * gn[i - 1] % P;
        for (int i = 0; i < tot; i += (mid << 1))
            for (int *x1 = x + i, *x2 = x + i + mid, *ed = x2, *g = gn; x1 != ed; ++x1, ++x2, ++g)
            {
                int p = *x1, q = LL(*x2) * *g % P;
                adj(*x1 = p + q - P), adj(*x2 = p - q);
            }
    }
    if (op == 1) return ;
    int t = qpow(tot);
    for (int i = 0; i < tot; ++i) x[i] = LL(t) * x[i] % P; 
} 
poly operator * (poly x, poly y)
{
    if (x.empty() || y.empty()) return {};
    int n = x.size(), m = y.size(), tot, bit;
    rdy(bit, tot, n + m), get_r(bit, tot);
    x.resize(tot), y.resize(tot);
    ntt(x.data(), tot, 1);
    if (x == y) y = x;
    else ntt(y.data(), tot, 1);
    for (int i = 0; i < tot; ++i) x[i] = LL(x[i]) * y[i] % P;
    ntt(x.data(), tot, -1);
    return x.resize(n + m - 1), x;
}
int main()
{
    scanf("%d %d", &n, &m);
    for (int i = 0; i <= n; ++i) scanf("%d", f + i), f[i] %= P;
    init();
    poly g(n + 1), h(2 * n + 1);
    for (int i = 0; i <= n; ++i)
    {
        g[i] = LL(f[i]) * ifac[i] % P * ifac[n - i] % P;
        if ((n - i) & 1) g[i] = P - g[i];
    }
    for (int i = 0; i <= (n << 1); ++i) h[i] = miv[i];
    poly l = g * h;
    printf("%d ", LL(mfac[n]) * l[n] % P);
    for (int i = 1; i <= n; ++i) printf("%d ", LL(mfac[i + n]) * imfac[i - 1] % P * l[i + n] % P);
    return 0;
}