【数学】多项式插值 - 拉格朗日插值

阅读目录

多项式插值简介
原理
模板
对于任意点求解
- 直接求解
- 转换成多项式的系数表达
对横坐标为连续整数的点求解
对横坐标为连续整数的点求解（不使用额外空间）
拓展
- 选取等差数列进行插值
- 重心拉格朗日插值法
验证

多项式插值简介

这部分内容收录在这篇文章中：【数学】多项式插值。

原理

由这 $n + 1$ 个点，可以构造对应的 $n + 1$ 个多项式，其中第 $i$ 的多项式的形式为：

l_{i} (x) = \prod_{j = 0, j \neq i}^{n} \frac{(x - x_{j})}{(x_{i} - x_{j})}

$l_i(x) = \prod\limits_{j = 0, \; j \neq i}^{n} \frac{(x - x_j)} {(x_i - x_j)}$

把下标为 $i$ 的点 $x_i$ 代入，可以得到：

\begin{array}{rcl} l_{i} (x_{i}) & = & \prod_{j = 0, j \neq i}^{n} \frac{(x_{i} - x_{j})}{(x_{i} - x_{j})} \\ = & \prod_{j = 0, j \neq i}^{n} 1 \\ = & 1 \end{array}

$\begin{eqnarray} l_i(x_i) &=& \prod\limits_{j = 0, \; j \neq i}^{n} \frac{(x_i - x_j)} {(x_i - x_j)} \nonumber \\ &=& \prod\limits_{j = 0, \; j \neq i}^{n} 1 \nonumber \\ &=& 1 \nonumber \end{eqnarray}$

把给定的点中，下标不为 $i$ 的某个点 $x_{j0}$ 代入，可以得到：

\begin{array}{rcl} l_{i} (x_{j 0}) & = & \prod_{j = 0, j \neq i}^{n} \frac{(x_{j 0} - x_{j})}{(x_{i} - x_{j})} \\ = & 0 \end{array}

$\begin{eqnarray} l_i(x_{j0}) &=& \prod\limits_{j = 0, \; j \neq i}^{n} \frac{(x_{j0} - x_j)} {(x_i - x_j)} \nonumber \\ &=& 0 \nonumber \end{eqnarray}$

因为分子部分总存在 $(x_{j0} - x_{j0}) = 0$ 所以上式为 $0$ 。

把不在给定的点中的任意点 $x$ 代入，显然有 $n$ 项，所以，这是一个 $n$ 次的多项式。

于是，把上面的 $n + 1$ 个不同的 $l_i(x)$ 和 $y_i$ 相乘，再相加，即可构造下面的多项式：

L (x) = \sum_{i = 0}^{n} y_{i} \cdot l_{i} (x)

$L(x) = \sum\limits_{i = 0}^{n} y_i \cdot l_i(x)$

由上面的推理可知， $L(x_i) = y_i$ ，并且 $L(x_i)$ 是一个不超过 $n$ 次的多项式（由于最高次的系数可能会被抵消，导致次数降低）。

由唯一性定理可知，这个 $L(x)$ 就是我们要找的 $P(x)$ ，即

P (x) = L (x)

$P(x) = L(x)$

回到顶部

模板

在代码中，为了统一我所有的模板和我一直以来的习惯，也就是从1开始计数的写法以及用 $n$ 表示点数，代码部分中的下标和其他文章（包括前述内容）中或许有差异。

本文中所有的代码都用 $n$ 表示点数，他们分别为 $(x_1, y_1), (x_2, y_2), (x_3, y_3), \cdots, (x_{n}, y_{n})$ ，也就是多项式的最高次数不超过 $n - 1$ 次。对一个不超过 $n - 1$ 次的多项式插值，至少需要 $n$ 个点。

回到顶部

对于任意点求解

直接求解

复杂度

时间复杂度 $O(n^2)$ ，空间复杂度 $O(1)$

使用说明

检查数据规模MAXN，注意数据规模指的是**点数 $n$ **（或者多项式的次数 $n - 1$ ）和模数MOD。
输入点数 $n$ ，和两个1-index的数组 $x, y$ 表示给定的点。
无需初始化，直接调用 ll lagrange_interpolation (ll x0, ll *x, ll *y, int n)，返回值即 $L(x_0)$ 。

代码

copynamespace LagrangeInterpolation {

static const int MOD = 998244353;

ll qpow (ll x, ll n) {
    ll res = 1;
    while (n) {
        if (n & 1) {
            res = res * x % MOD;
        }
        x = x * x % MOD;
        n >>= 1;
    }
    return res;
}

ll inv (ll x) {
    return qpow (x, MOD - 2);
}

/**
 * 用 n 个点 (x[1], y[1]), (x[2], y[2]), ..., (x[n], y[n])
 * 插值得到 n - 1 次多项式 L，直接求解 L(x0) 的值
 *
 * 时间复杂度: O(n^2)
 * 空间复杂度：O(1)
*/
ll lagrange_interpolation (ll x0, ll *x, ll *y, int n) {
    for (int i = 1; i <= n; ++i) {
        if (x0 == x[i]) {
            return y[i];
        }
    }
    ll Lx0 = 0;
    ll P = 1;
    for (int i = 1; i <= n; ++i) {
        P = P * (x0 - x[i] + MOD) % MOD;
    }
    for (int i = 1; i <= n; ++i) {
        ll q = 1;
        for (int j = 1; j <= n; ++j) {
            if (j == i) {
                continue;
            }
            q = q * (x[i] - x[j] + MOD) % MOD;
        }
        ll p = P * inv (x0 - x[i] + MOD) % MOD;
        ll lix0 = p * inv (q) % MOD;
        Lx0 = (Lx0 + lix0 * y[i] % MOD) % MOD;
    }
    return Lx0;
}

};

转换成多项式的系数表达

原理

当 $x \neq x_i$ 恒成立时：

\begin{array}{rcl} l_{i} (x) & = & \prod_{j = 0, j \neq i}^{n} \frac{(x - x_{j})}{(x_{i} - x_{j})} \\ = & \frac{\prod_{j = 0, j \neq i}^{n} (x - x_{j})}{\prod_{j = 0, j \neq i}^{n} (x_{i} - x_{j})} \\ = & \frac{\prod_{j = 0}^{n} (x - x_{j})}{(x - x_{i}) \prod_{j = 0, j \neq i}^{n} (x_{i} - x_{j})} \end{array}

$\begin{eqnarray} l_i(x) &=& \prod\limits_{j = 0, \; j \neq i}^{n} \frac{(x - x_j)} {(x_i - x_j)} \nonumber \\ &=& \frac{\prod\limits_{j = 0, \; j \neq i}^{n} (x - x_j)} {\prod\limits_{j = 0, \; j \neq i}^{n} (x_i - x_j)} \nonumber \\ &=& \frac{\prod\limits_{j = 0}^{n} (x - x_j)} {(x - x_i) \prod\limits_{j = 0, \; j \neq i}^{n} (x_i - x_j)} \nonumber \end{eqnarray}$

记 $P = \prod\limits_{j = 0}^{n} (x - x_j)$ 这是一个和 $i$ 无关的量，结果是一个 $n$ 次多项式
记 $C = \prod\limits_{j = 0, \; j \neq i}^{n} (x_i - x_j)$ 这是一个常数
记 $Q_i = (x - x_i)$ 这是一个二项式

得到

l_{i} (x) = P \cdot \frac{1}{C_{i} \cdot Q_{i}}

$l_i(x) = P \cdot \frac{1}{C_i \cdot Q_i}$

预处理多项式 $P$ ，然后用多项式除法除以 $Q_i$ 再每一项除以常数 $C_i$ 即可得到所需的系数，最后全部加起来得到 $L(x)$ 的系数。

复杂度

预处理系数：时间复杂度 $O(n^2)$ ，空间复杂度 $O(n)$
单次求值：时间复杂度 $O(n)$ ，空间复杂度 $O(1)$

回到顶部

对横坐标为连续整数的点求解

根据 $x[] = 1, 2, 3, \cdots, k$ 时的值 $y[]$ ，
插值计算不超过 $k-1$ 次的多项式 $L$ 在 $x_0$ 处的取值 $L(x_0)$ 。

时间复杂度 $O(k)$ ，空间复杂度 $O(MAXK)$ 。

$P1[i]$ 是 $x - i$ 的前缀积，
$P2[i]$ 是 $x - i$ 的后缀积，
$Q[i]$ 是阶乘 $i!$ 的逆元。

注意：
1. $k$ 为点的数量， $MAXK$ 是否足够大？
2. $y[]$ 是否传入负数？
3. $int$ 是否可能溢出？
4. $MOD$ 必须是大于其他所有数字的质数，否则逆元会不存在。

卡常：
1. $q$ 只与 $k$ 有关，对于固定的 $k$ ，可以预处理后多次使用，但是真没什么必要。

使用：
1. 直接调用 int lagrange (int x0, int *y, int k) 即可，其会自动初始化需要的数组，并且可以多次调用。

copynamespace Lagrange1 {

const int MAXK = 1e6 + 5;
static int P1[MAXK], P2[MAXK], Q[MAXK];

ll qpow (ll x, ll n) {
    ll res = 1;
    for (; n; n >>= 1) {
        if (n & 1) res = res * x % MOD;
        x = x * x % MOD;
    }
    return res;
}

void init_P (int x0, int k) {
    P1[0] = 1;
    for (int i = 1; i <= k; ++i)
        P1[i] = 1LL * P1[i - 1] * (x0 - i) % MOD;
    P2[k + 1] = 1;
    for (int i = k; i >= 1; --i)
        P2[i] = 1LL * P2[i + 1] * (x0 - i) % MOD;
}

void init_Q (int k) {
    if (Q[k] != 0) return;
    Q[k] = 1;
    for (int i = 1; i <= k; ++i)
        Q[k] = 1LL * Q[k] * i % MOD;
    Q[k] = qpow (Q[k], MOD - 2);
    for (int i = k; i >= 1; --i)
        Q[i - 1] = 1LL * Q[i] * i % MOD;
}

int lagrange (int x0, int *y, int k) {
    if (x0 <= k) return y[x0];
    init_P (x0, k);
    init_Q (k);
    int Lx0 = 0;
    for (int i = 1; i <= k; ++i) {
        ll p = 1LL * P1[i - 1] * P2[i + 1] % MOD;
        ll q = 1LL * Q[i - 1] * Q[k - i] % MOD;
        if ( (k - i) & 1) q = MOD - q;
        Lx0 = (Lx0 + p * q % MOD * y[i] % MOD) % MOD;
    }
    return Lx0;
}

};

using Lagrange1::lagrange;

回到顶部

对横坐标为连续整数的点求解（不使用额外空间）

时间复杂度 $O(k\log{M})$ ，空间复杂度 $O(1)$ 。

使用：
1. 直接调用 int lagrange (int x0, int *y, int k) 即可，并且可以多次调用。

copynamespace Lagrange2 {

ll qpow (ll x, ll n) {
    ll res = 1;
    for (; n; n >>= 1) {
        if (n & 1) res = res * x % MOD;
        x = x * x % MOD;
    }
    return res;
}

ll inv (ll x) {
    return qpow (x, MOD - 2);
}

int lagrange (int x0, int *y, int k) {
    if (x0 <= k)
        return y[x0];
    ll Lx0 = 0, P = 1, Q1 = 1, Q2 = 1;
    for (int i = 1; i <= k; ++i) {
        P = P * (x0 - i) % MOD;
        if (i < k) Q2 = Q2 * (MOD - i) % MOD;
    }
    for (int i = 1; i <= k; ++i) {
        ll p = P * inv (x0 - i) % MOD;
        ll q = inv (Q1 * Q2 % MOD);
        Lx0 = (Lx0 + p * q % MOD * y[i]) % MOD;
        Q1 = Q1 * i % MOD;
        Q2 = Q2 * inv (MOD - k + i) % MOD;
    }
    return Lx0;
}

};

using Lagrange2::lagrange;

回到顶部

拓展

选取等差数列进行插值

例如：求 $\sum\limits_{i=1}^ni^{k}$ 的值，众所周知，这是一个 $k+1$ 次多项式，至少需要 $k+2$ 个点。

这个值有很多种求法，有 $O(n\log k)$ 的朴素解法（枚举+快速幂），有 $O(k^2)$ 的系数递推法，使用上述的适用范围广泛的拉格朗日插值法，也可以做到 $O(k^2)$ ，算法的瓶颈在于计算插值多项式的分母部分。

这里可以选取一些特殊点来加速计算插值多项式的分母部分。例如通过选取一个等差数列，那么这个分母的临项会变得非常有规律。简单起见可以取正整数。

观察分子部分，对每个 $i$ 来说，分子部分乘上 $(x-x_i)$ 后，都变为公共的 $\prod\limits_{j=0}^{k}(x-x_j)$ ，那么在使用的时候再把这个 $(x-x_i)$ 除掉就行，复杂度为 $O(k)$ 预处理， $O(\log M)$ 单次使用，复杂度为 $O(k\log M)$ 。算法的复杂度瓶颈在于计算分母部分。

观察分母部分，在选取 $x_i=i$ 后，分母变为两个自然数的前缀积之积，符号是正负交替。更准确地说，第 $i$ 项的值是： $[\prod\limits_{j=0}^{i-1}(i-j)][\prod\limits_{j=i+1}^{k}(i-j)]=(\prod\limits_{j=1}^{i}j)(\prod\limits_{j=1}^{k-i}j)(-1)^{k-i}$

由于算法竞赛一般是在模 $M$ 意义下进行，算上求解逆元时间复杂度应为 $O(k\log{M})$ 。

注： $\sum_{i=1}^ni^{k}$ 的暴力求解其实可以是 $O(n)$ 的，是使用线性筛去筛出每个数的 $k$ 次方，然后线性递推出前缀和，不过这实际上是在画蛇添足，因为快速幂的常数实际上很小，而且 $k$ 本身不大。

卡常：对同一个多项式多次求值时，插值完成后可以保存其分母（毒瘤）。

重心拉格朗日插值法

观察上面这个式子：

l_{i} (x_{0}) = \prod_{j = 1}^{k} [i \neq j] \frac{(x_{0} - x_{j})}{(x_{i} - x_{j})}

$l_i(x_0)=\prod\limits_{j=1}^{k} [i\neq j]\frac{(x_0-x_j)}{(x_i-x_j)}$

L (x_{0}) = \sum_{i = 1}^{k} y_{i} l_{i} (x_{0})

$L(x_0)=\sum\limits_{i=1}^{k} y_il_i(x_0)$

若记 $g=\prod\limits_{j=1}^{k} (x_0-x_j)$ ，代入上式得到：

L (x_{0}) = \sum_{i = 1}^{k} y_{i} \prod_{j = 1}^{k} [i \neq j] \frac{(x_{0} - x_{j})}{(x_{i} - x_{j})} = g \sum_{i = 1}^{k} \prod_{j = 1}^{k} [i \neq j] \frac{y_{i}}{(x_{0} - x_{i}) (x_{i} - x_{j})}

$L(x_0)=\sum\limits_{i=1}^{k} y_i\prod\limits_{j=1}^{k} [i\neq j]\frac{(x_0-x_j)}{(x_i-x_j)} \\ = g\sum\limits_{i=1}^{k} \prod\limits_{j=1}^{k} [i\neq j]\frac{y_i}{(x_0-x_i)(x_i-x_j)}$

若记 $t_i=\prod\limits_{j=1}^{k} [i\neq j]\frac{y_i}{(x_i-x_j)}$ ，代入上式得到