拉格朗日插值学习笔记

拉格朗日插值学习笔记

简介

在 2022 年联合省选 D1T2 吃了这个的亏，来补一补。虽说我压根就没看出来跟插值有啥关系。

我们都知道，\(n\) 个点的点值 \((x_i,y_i)\) 可以唯一确定一个 \(n-1\) 次多项式（为了叙述方便，本文中所有“\(k\) 次多项式”“\(k\) 次函数”的最高次项系数可以为 \(0\)）。

拉格朗日插值就是用来求这个多项式的。

例子

例如，我们已知四个点值 \((-1,1)(0,2)(0.5,1.375)(1,1)\)，要求过这四个点的三次函数 \(f\)。

当然，你可以直接待定系数用高斯消元解方程，但那是 \(\mathcal{O}(n^3)\) 的，拉格朗日插值可以在 \(\mathcal{O}(n^2)\) 内求解。

约瑟夫·拉格朗日认为这个函数可以用四个三次函数线性组合出来。

首先构造一个三次函数 \(f_1\)，在 \(x=-1\) 的取值为 \(1\)，但在其他三个点的取值为 \(0\)。

类似地，构造 \(f_{2,3,4}\) 依次在每个点取值为 \(1\)，在其他三个点取值为 \(0\)。

画到一张图里就是这样：

那么这几个函数有啥用呢？

容易发现，\(f(x)=y_1f_1(x)+y_2f_2(x)+y_3f_3(x)+y_4f_4(x)\)，把那四个点点值代入进去就可以知道。

现在问题就转化为了怎么求 \(f_{1,2,3,4}\)。

类似于一次函数的两点式：

\[\frac{y-y_1}{y_2-y_1}=\frac{x-x_1}{x_2-x_1} \]

我们可以构造函数：

\[\begin{aligned} f_1(x)&=\frac{(x-x_2)(x-x_3)(x-x_4)}{(x_1-x_2)(x_1-x_3)(x_1-x_4)}\\ f_2(x)&=\frac{(x-x_1)(x-x_3)(x-x_4)}{(x_2-x_1)(x_2-x_3)(x_2-x_4)}\\ f_3(x)&=\frac{(x-x_1)(x-x_2)(x-x_4)}{(x_3-x_1)(x_3-x_2)(x_3-x_4)}\\ f_4(x)&=\frac{(x-x_1)(x-x_2)(x-x_3)}{(x_4-x_1)(x_4-x_2)(x_4-x_3)}\\ \end{aligned} \]

于是就做完了。

推导

是时候写一下推导了。

对于 \(n\) 个点值求 \(n-1\) 次多项式的问题，我们先有点值 \((x_j,y_j)\)（\(1\le j\le n\)），设函数 \(f_i\)（\(1\le i\le n\)），它们是 \(n-1\) 次函数，且满足：

\[f_i(x_j)= \begin{cases} 1,&i=j\\ 0,&i\ne j\\ \end{cases} \]

则 \(f_i\) 可以写成：

\[f_i(x)=\prod\limits_{j=1\\j\ne i}^n\frac{x-x_j}{x_i-x_j} \]

最终得到：

\[f(x)=\sum\limits_{i=1}^ny_if_i(x) \]

代码

题目是 洛谷 P4781 拉格朗日插值

套上面公式直接算即可。

如果你在计算每个 \(\dfrac{x-x_j}{x_i-x_j}\) 的时候都求一遍逆元，会导致复杂度变为 \(\mathcal{O}(n^2\log n)\)，多带一个逆元的 \(\log\)。为了让复杂度瓶颈不在逆元上，我们通常分开算分子分母，在每个函数算完后再进行有理数取模，这样的复杂度为 \(\mathcal{O}(n^2)\)。

//By: Luogu@rui_er(122461)
#include <bits/stdc++.h>
#define rep(x,y,z) for(int x=y;x<=z;x++)
#define per(x,y,z) for(int x=y;x>=z;x--)
#define debug printf("Running %s on line %d...\n",__FUNCTION__,__LINE__)
#define fileIO(s) do{freopen(s".in","r",stdin);freopen(s".out","w",stdout);}while(false)
using namespace std;
typedef long long ll;
const int N = 2e3+5, mod = 998244353; 

int n, k, x[N], y[N], ans;
template<typename T> void chkmin(T& x, T y) {if(x > y) x = y;}
template<typename T> void chkmax(T& x, T y) {if(x < y) x = y;}
int qpow(int x, int y) {
    int ans = 1;
    for(;y;y>>=1,x=1LL*x*x%mod) if(y & 1) ans = 1LL * ans * x % mod;
    return ans;
}
int inv(int x) {return qpow(x, mod-2);}

int main() {
    scanf("%d%d", &n, &k);
    rep(i, 1, n) scanf("%d%d", &x[i], &y[i]);
    rep(i, 1, n) {
        int p = y[i], q = 1;
        rep(j, 1, n) {
            if(i != j) {
                p = 1LL * p * (k - x[j]) % mod;
                q = 1LL * q * (x[i] - x[j]) % mod;
            }
        }
        ans = (ans + 1LL * p * inv(q) % mod + mod) % mod;
    }
    printf("%d\n", ans);
    return 0;
}

连续点值的插值

如果已知的点值是连续点的点值，我们可以做到 \(\mathcal{O}(n)\) 的插值。

有时候发现了一个函数是 \(n\) 次多项式，就求 \(n+1\) 个点值进去插值。为了省事，这里我们令 \(x_i=i\)（\(1\le i\le n+1\)）。注意这里 \(n\) 与上面意义不一样，是次数而不是点数。

我们有拉插公式：

\[f(x)=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j=1\\j\ne i}^{n+1}\frac{x-x_j}{x_i-x_j} \]

代入 \(x_i=i\)：

\[f(x)=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j=1\\j\ne i}^{n+1}\frac{x-j}{i-j} \]

考虑怎么快速求 \(\displaystyle\prod\limits_{j=1\\j\ne i}^{n+1}\frac{x-j}{i-j}\)。

这个东西的分子是:

\[\frac{\prod\limits_{j=1}^{n+1}(x-j)}{x-i} \]

分母的话把 \(i-j\) 累乘拆成两段阶乘，就是：

\[(-1)^{n+1-i}\cdot(i-1)!\cdot(n+1-i)! \]

于是连续点值的插值公式：

\[f(x)=\sum\limits_{i=1}^{n+1}y_i\cdot\frac{\prod\limits_{j=1}^{n+1}(x-j)}{(x-i)\cdot(-1)^{n+1-i}\cdot(i-1)!\cdot(n+1-i)!} \]

预处理前后缀积、阶乘阶乘逆然后代这个式子的复杂度为 \(\mathcal{O}(n)\)。

例子是 CodeForces 622F The Sum of the k-th Powers。

可以发现答案是 \(k+1\) 次多项式，因此代 \(k+2\) 个点进去拉插。

本题的 \(i^k\) 可以线性筛，因此复杂度可以做到 \(\mathcal{O}(n)\)。

//By: Luogu@rui_er(122461)
#include <bits/stdc++.h>
#define rep(x,y,z) for(int x=y;x<=z;x++)
#define per(x,y,z) for(int x=y;x>=z;x--)
#define debug printf("Running %s on line %d...\n",__FUNCTION__,__LINE__)
#define fileIO(s) do{freopen(s".in","r",stdin);freopen(s".out","w",stdout);}while(false)
using namespace std;
typedef long long ll;
const int N = 1e6+5, mod = 1e9+7; 

int n, k, tab[N], p[N], pcnt, f[N], pre[N], suf[N], fac[N], inv[N], ans;
template<typename T> void chkmin(T& x, T y) {if(x > y) x = y;}
template<typename T> void chkmax(T& x, T y) {if(x < y) x = y;}
int qpow(int x, int y) {
	int ans = 1;
	for(;y;y>>=1,x=1LL*x*x%mod) if(y & 1) ans = 1LL * ans * x % mod;
	return ans;
}
void sieve(int lim) {
    f[1] = 1;
    rep(i, 2, lim) {
        if(!tab[i]) {
            p[++pcnt] = i;
            f[i] = qpow(i, k);
        }
        for(int j=1;j<=pcnt&&1LL*i*p[j]<=lim;j++) {
            tab[i*p[j]] = 1;
            f[i*p[j]] = 1LL * f[i] * f[p[j]] % mod;
            if(!(i % p[j])) break;
        }
    }
    rep(i, 2, lim) f[i] = (f[i-1] + f[i]) % mod;
}

int main() {
    scanf("%d%d", &n, &k);
    sieve(k+2);
    if(n <= k + 2) return printf("%d\n", f[n])&0;
    pre[0] = suf[k+3] = 1;
    rep(i, 1, k+2) pre[i] = 1LL * pre[i-1] * (n - i) % mod;
    per(i, k+2, 1) suf[i] = 1LL * suf[i+1] * (n - i) % mod;
    fac[0] = inv[0] = fac[1] = inv[1] = 1;
    rep(i, 2, k+2) {
        fac[i] = 1LL * fac[i-1] * i % mod;
        inv[i] = 1LL * (mod - mod / i) * inv[mod%i] % mod;
    }
    rep(i, 2, k+2) inv[i] = 1LL * inv[i-1] * inv[i] % mod;
    rep(i, 1, k+2) {
        int P = 1LL * pre[i-1] * suf[i+1] % mod;
        int Q = 1LL * inv[i-1] * inv[k+2-i] % mod;
        int mul = ((k + 2 - i) & 1) ? -1 : 1;
        ans = (ans + 1LL * (Q * mul + mod) % mod * P % mod * f[i] % mod) % mod;
    }
    printf("%d\n", ans);
    return 0;
}

重心拉格朗日插值

不会。