Dyd's Blog

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

拉格朗日插值

Posted on In
Symbols count in article: 4.3k Reading time ≈ 4 mins.

已知点值求多项式

拉格朗日插值

as we all know ， $n$ 个点可以唯一确定一个 $n - 1$ 次多项式，那么现在，给你 $n$ 个点，你该如何确定这个多项式呢？当然，你可以 gauss 消元解方程，但过于麻烦，且复杂度较高，这里就有一个新方法——拉格朗日插值法

式子

给定 $n$ 个点 $(x_1, y_1), (x_2, y_2), …, (x_n, y_n)$ 它们确定的 $n - 1$ 次多项式为：
$$
f(x) = \sum_{i = 1}^n y_i \prod_{j \ne i}\frac{(x - x_j)}{(x_i - x_j)}
$$
显然，当 $x = x_i$ 的时候， $y_i$ 后面的系数为 $1$ ，而 $\forall j \ne i, y_j$ 后面的系数都为 $0$ （因为有 $(x - x_i)$ 这一项），那么 $f(x)$ 就过点 $(x_i, y_i)$ ，这就是我们要求的多项式

一般拉插

【模板】拉格朗日插值

就是上面的式子，预处理是 $O(n^2)$ （有些题要预处理逆元但一般瓶颈不在求逆元上），每次带入求点是 $O(n^2)$ ，想得到原多项式也是 $O(n^2)$ 的

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
#include <cstdio>
#include <algorithm>
#define fi first
#define se second
using LL = long long;
const int N = 2e3 + 100, P = 998244353;
int n, b, la[N];
std::pair<int, int> a[N];
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()
{
    for (int i = 1, t; i <= n; ++i)
    {
        la[i] = 1;
        for (int j = 1; j <= n; ++j) if (j != i)
        {
            adj(t = a[i].fi - a[j].fi);
            la[i] = LL(t) * la[i] % P;
        }
        la[i] = qpow(la[i]);
    }   
}
int lag(int x)
{
    int res = 0;
    for (int i = 1, t, tt; i <= n; ++i)
    {
        tt = la[i];
        for (int j = 1; j <= n; ++j) if (j != i)
        {
            adj(t = x - a[j].fi);
            tt = LL(t) * tt % P;
        }
        adj(res += LL(tt) * a[i].se % P - P);
    }
    return res;
}
int main()
{
    scanf("%d %d", &n, &b);
    for (int i = 1; i <= n; ++i) scanf("%d %d", &a[i].fi, &a[i].se);
    init();
    printf("%d\n", lag(b));
    return 0;
}

线性拉插

在很多题目中，取的 $n$ 个点是连续的（或者让我们自己取，那就取连续的），此时我们一般可以做到询问线性

考虑我们求的点是 $(1, f(1)), (2, f(2)), …, (n, f(n))$ ，那么带入式子就是：
$$
\begin{aligned}
f(x) = & \sum_{i = 1}^n f(i) \prod_{j \ne i}\frac{(x - j)}{(i - j)} \\
= & \sum_{i = 1}^n f(i) * (-1)^{n - i} * \frac{(x - 1)(x - 2)…(x - (i - 1))(x - (i + 1))(x - (i + 2))…(x - n)}{(i - 1)!(n - i)!} \\
\end{aligned}
$$
预处理 $f(i)$ 和分母，复杂度的关键在 $f(i)$ 的求取和逆元的处理，每次查询的时候维护前、

后缀积，可以做到线性查询

值得一提的是，如果 $x - x_i$ 的逆元好求或者求它不在瓶颈上，查询时可以先算出 $\prod (x - x_i)$ ，然后每项用逆元除去

The Sum of the k-th Powers

经典的求自然数幂和，一个结论是： $\sum_{i = 1}^n i^k$ 是一个关于 $n$ 的 $k + 1$ 次多项式，复杂度瓶颈在预处理求 $f(i)$ 的地方，总复杂度 $O(k \log k)$ 

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
#include <cstdio>
#include <algorithm>
#define fi first
#define se second
using LL = long long;
const int N = 1e6 + 100, P = 1e9 + 7;
int n, k;
std::pair<int, int> a[N];
int la[N], la2[N];
int ifac[N];
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 n)
{
    int fac = 1, fi = 0;
    ifac[0] = ifac[1] = 1;
    for (int i = 2; i <= n; ++i) fac = LL(fac) * i % P;
    ifac[n] = qpow(fac);
    for (int i = n - 1; i > 1; --i) ifac[i] = LL(i + 1) * ifac[i + 1] % P;
    for (int i = 1; i <= n; ++i)
    {
        adj(fi += qpow(i, k) - P);
        la[i] = LL(fi) * ifac[i - 1] % P * ifac[n - i] % P;
        if ((n - i) & 1) la[i] = P - la[i];
    }
}
int lag(int x, int n)
{
    int res = 0;
    for (int i = 1, t = 1; i <= n; ++i)
    {
        la2[i] = t;
        adj(t = LL(x - i) * t % P);    
    }
    for (int i = n, t = 1; i; --i)
    {
        la2[i] = LL(t) * la2[i] % P;
        adj(t = LL(x - i) * t % P);    
    }
    for (int i = 1; i <= n; ++i) adj(res += LL(la2[i]) * la[i] % P - P); 
    return res;
}
int main()
{
    scanf("%d %d", &n, &k);
    init(k + 2);
    printf("%d\n", lag(n, k + 2));
    return 0;
}

重心拉插

设 $g = \prod (x - x_i), w_i = \frac{y_i}{\prod_{j \ne i} (x_i - x_j)}$ ，则式子为：
$$
\begin{aligned}
f(x) = g \sum_{i = 1}^n \frac{w_i}{x - x_i}
\end{aligned}
$$
预处理是 $O(n^2)$ 不变，每次计算的复杂度变为 $O(n * t)$ ，其中 $t$ 是计算逆元的复杂度，但这个式子最好的地方在于，可以用 $O(n * t)$ 的时间插入一个新的点，只要计算出 $w_i$ 即可

注意处理 $x$ 相同的情况

拉格朗日插值

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
#include <cstdio>
#include <vector>
#define eb emplace_back
using LL = long long;
const int N = 3000 + 100, P = 998244353;
int n;
std::vector<int> xn, yn, w;
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 adj(int &x){ x += (x >> 31) & P; }
void ins(int x, int y)
{
    yn.eb(y);
    int t;
    for (int i = xn.size() - 1; ~i; --i)
    {
        adj(t = x - xn[i]);
        t = qpow(t);
        y = LL(y) * t % P;
        w[i] = LL(P - t) * w[i] % P;
    }
    xn.eb(x), w.eb(y);
}
int ask(int x)
{
    int res = 0, t, g = 1;
    for (int i = xn.size() - 1; ~i; --i)
        if (xn[i] == x) return yn[i];
        else
        {
            adj(t = x - xn[i]);
            g = LL(t) * g % P;
            t = LL(w[i]) * qpow(t) % P;
            adj(res += t - P);
        }
    return LL(g) * res % P;
}
int main()
{
    scanf("%d", &n);
    for (int op, x, y; n--; )
    {
        scanf("%d %d", &op, &x);
        if (op == 1)
        {
            scanf("%d", &y);
            ins(x, y);
        }
        else printf("%d\n", ask(x));
    }
    return 0;
}

后话

其实还有一种 拉插 + 分治 ntt + 多项式多点求值 的 $O(n \log^2 n)$ 快速插值方法，后面有时间补充吧