Skip to content

题解-P5667

Posted on:2023年7月11日 at 14:54

本文章遵守知识共享协议 CC-BY-NC-SA,同步发表于洛谷题解区,转载时须在文章的任一位置附上原文链接和作者署名(rickyxrc)。推荐在我的个人博客阅读。

这是我的第一篇 NTT 题解,思路和其它题解几乎一样,不过推式子的过程和讲解更加详细,希望可以成为新人学习 NTT 的有力基础,当然会有讲的不周到的地方,欢迎各位补充与指正。

阅读本篇题解,您需要掌握拉格朗日插值的基本知识并 AC 拉格朗日插值

题面大意

给定一个不超过 nn 次的多项式的 n+1n+1 个点值 f(0),f(1)f(n)f(0),f(1) \dots f(n),和一个正整数 mm,求 f(m),f(m+1)f(m+n)f(m),f(m+1) \dots f(m+n)

解题思路

首先我们写出拉格朗日插值的式子:

li(x)=yijixxjxixjl_i(x) = y_i \prod_{j\neq i}^{}\frac{x-x_j}{x_i-x_j}
f(x)=i=0nli(x)f(x) = \sum_{i=0}^{n} l_i(x)

这个式子我不会证明,不过简单地来想想,将 xix_i 带入该式子,就会将除了 li(x)l_i(x) 之外的值全部消掉(等于 0),得到 yiy_i,所以可以认为它是对的(当然这显然不是严谨的证明)。

把所有的 li(x)l_i(x) 带入 f(x)f(x),我们就得到了拉格朗日插值的另一种表达形式:

f(x)=i=0nyijixxjxixjf(x) = \sum_{i=0}^{n}y_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j}

我们现在需要求出的是 f(m),f(m+1),...,f(m+n)f(m),f(m+1),...,f(m+n),我们做如下推导:

f(m+x)=i=0nyiji(m+x)xjxixjf(m+x) = \sum_{i=0}^{n}y_i\prod_{j\neq i}\frac{(m+x)-x_j}{x_i-x_j}

因为给入的 xx 是连续的且是 [0,n][0,n],所以式子可以简化如下:

f(m+x)=i=0nf(i)ji(m+x)jijf(m+x) = \sum_{i=0}^{n}f(i)\prod_{j\neq i}\frac{(m+x)-j}{i-j}

这一段连乘可以变成阶乘的形式,但是我没一步推出来所以我们写详细一点。

我们先下放连乘符号,于是式子可以写成这样:

f(m+x)=i=0nf(i)ji(m+xj)ji(ij)f(m+x) = \sum_{i=0}^{n}f(i)\frac{\prod_{j\neq i}(m+x-j)}{\prod_{j\neq i}(i-j)}

然后我们先看分子:

ji(m+xj)\prod_{j\neq i}(m+x-j)
=(m+x1)×(m+x2)×(m+x3)××(m+xn)m+xi= \frac{(m+x-1) \times (m+x-2) \times(m+x-3) \times \cdots \times (m+x-n)}{m+x-i}

明显可以简化成阶乘的式子:

=(mx1)!(m+xi)×(m+xn1)!= \frac{(m-x-1)!}{(m+x-i)\times(m+x-n-1)!}

然后看分母:

ji(ij)\prod_{j\neq i}(i-j)

这个式子可以化成:

j<i(ij)j=i+1n(ij)\prod_{j < i}(i-j) \prod_{j=i+1}^{n}(i-j)
=i!×j=i+1n(ji)=i!\times\prod_{j=i+1}^{n}-(j-i)
=i!×(1)ni×j=i+1n(ji)=i!\times(-1)^{n-i}\times\prod_{j=i+1}^{n}(j-i)
=i!×(1)ni×(ni)!=i!\times(-1)^{n-i}\times(n-i)!

所以整段分式就可以化成:

(mx1)!(m+xi)×(m+xn1)!i!×(1)ni×(ni)!\frac{\frac{(m-x-1)!}{(m+x-i)\times(m+x-n-1)!}}{i!\times(-1)^{n-i}\times(n-i)!}
=(mx1)!(m+xi)×(m+xn1)!×i!×(1)ni×(ni)!=\frac{(m-x-1)!}{(m+x-i)\times(m+x-n-1)!\times i!\times(-1)^{n-i}\times(n-i)!}

原式就变成了:

f(m+x)=i=0nf(i)×(mx1)!(m+xi)×(m+xn1)!×i!×(1)ni×(ni)!f(m+x) = \sum_{i=0}^{n}\frac{f(i)\times(m-x-1)!}{(m+x-i)\times(m+x-n-1)!\times i!\times(-1)^{n-i}\times(n-i)!}

提出公因式:

=(mx1)!(m+xn1)!i=0nf(i)(m+xi)×i!×(1)ni×(ni)!= \frac{(m-x-1)!}{(m+x-n-1)!} \sum_{i=0}^{n}\frac{f(i)}{(m+x-i)\times i!\times(-1)^{n-i}\times(n-i)!}

因为我们需要 NTT 求解,所以我们选择将其凑成卷积的形式:

=(mx1)!(m+xn1)!i=0n1(m+xi)×f(i)(ni)!×i!×(1)ni= \frac{(m-x-1)!}{(m+x-n-1)!} \sum_{i=0}^{n}\frac{1}{(m+x-i)}\times\frac{f(i)}{(n-i)!\times i!\times(-1)^{n-i}}

然后答案就比较显然了,为两个数列卷积之后再略微处理的形式。

形式化地说,设 Ai=f(i)(ni)!×i!×(1)niA_i=\frac{f(i)}{(n-i)!\times i!\times(-1)^{n-i}}Bi=1mn+iB_i = \frac{1}{m-n+i},则令 F=ABF = A * BFi×(mx1)!(m+xn1)!\frac{F_i\times(m-x-1)!}{(m+x-n-1)!} 即为答案。

然后需要求逆元和阶乘逆元,最好能线性求。

因为阶乘可能很大,所以我们不能预处理阶乘逆元,而是预处理上升幂,参考代码中 jinv 和 mjcsinv 的部分。

代码如下:

#include <stdio.h>

#define maxn 1000007

typedef long long i64;

i64 g = 3, gi;
const i64 mod = 998244353;

i64 t[maxn], a[maxn], b[maxn], n, m, len = 1, l, r[maxn], jcs[maxn], jinv[maxn], mjcs[maxn], mjcsinv[maxn], minv[maxn], ninv[maxn];

inline i64 pow(i64 x, i64 p);
inline i64 inv(i64 x);
void ntt(i64 *c, i64 op);

void getinv(i64 n, i64 m)
{
    jcs[0] = mjcs[0] = 1;
    for (i64 i = 1; i <= 2 * n + 1; i++)
        jcs[i] = jcs[i - 1] * i % mod,
        mjcs[i] = mjcs[i - 1] * (m - n + i - 1) % mod;
    jinv[2 * n + 1] = inv(jcs[2 * n + 1]);
    mjcsinv[2 * n + 1] = inv(mjcs[2 * n + 1]);
    for (int i = 2 * n + 1; i; i--)
        jinv[i - 1] = jinv[i] * i % mod,
                 mjcsinv[i - 1] = mjcsinv[i] * (m - n + i - 1) % mod,
                 ninv[i] = jinv[i] * jcs[i - 1] % mod,
                 minv[i] = mjcsinv[i] * mjcs[i - 1] % mod;
    minv[0] = 1;

    //	for(int i=1;i<=2*n+1;i++)
    //		printf("[%lld %lld %lld %lld] ",
    //			i*ninv[i]%mod,(m-n+i-1)*minv[i]%mod,jcs[i]*jinv[i]%mod,mjcs[i]*mjcsinv[i]%mod);
}

int main()
{
    gi = inv(g);

    scanf("%lld%lld", &n, &m);

    for (i64 i = 0; i <= n; i++)
        scanf("%lld", t + i);

    getinv(n, m);

    for (i64 i = 0; i <= n; i++)
    {
        a[i] = t[i] * jinv[i] % mod * jinv[n - i] % mod;
        if ((n - i) & 1)
            a[i] = mod - a[i];
    }
    for (i64 i = 0; i <= 2 * n; i++)
        b[i] = minv[i + 1];

    while (len <= 2 * n)
        len <<= 1, l++;

    for (i64 i = 0; i < len; i++)
        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));

    ntt(a, 1);
    ntt(b, 1);
    for (i64 i = 0; i < len; i++)
        a[i] = a[i] * b[i] % mod;
    ntt(a, -1);

    i64 linv = pow(len, mod - 2);
    for (i64 i = n; i <= 2 * n; i++)
        printf("%lld ", mjcs[i + 1] * a[i] % mod * mjcsinv[i - n] % mod * linv % mod);
    return 0;
}


在 Rickyxrc's blog 出现的文章,若无特殊注明,均采用 CC BY-NC-SA 4.0 协议共享,也就是转载时需要注明本文章的地址,并且引用本文章的文章也要使用相同的方式共享。