「笔记」拉格朗日插值


简介

对于 \(k\) 次多项式函数 \(F(x)\)
若已知 \(k+1\) 个点值,则可构造出多项式。
有:

\[F(x) = \sum_{i=1}^{k+1}y_i\prod_{i\not ={j}}\dfrac{x-x_j}{x_i-x_j} \]

正确性详见 拉格朗日插值


正确性

\(k + 1\) 个点值代入即可检验。
\(x=x_k\) 时,对 \(i\) 进行讨论:

  1. \(i\not ={k}\) 时,存在一个 \(j\) 满足 \(j=k\)
    对于乘积项的分子 \(\prod\limits_{i\not ={j}}x-x_j\),当 \(j=k\) 时,\(x-x_j = 0\)
    对答案无贡献。
  2. \(i=k\) 时,乘积项变为 \(\prod\limits_{i\not ={j}}\dfrac{x_i-x_j}{x_i-x_j}\)
    其值恒等于 \(1\)

模板题

P4781 【模板】拉格朗日插值

//
/*
By:Luckyblock
*/
#include <cstdio>
#include <cctype>
#include <cstdlib>
#include <algorithm>
#define ll long long
const int MARX = 2010;
const ll Mod = 998244353;
//=============================================================
ll ans, N, K, X[MARX], Y[MARX];
//=============================================================
inline int read()
{
    int f = 1, w = 0; char ch = getchar();
    for(; !isdigit(ch); ch = getchar()) if(ch == '-') f = -1;
    for(; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
    return f * w;
}
ll qpow(ll x, ll y, ll mod)
{
    ll ret = 1;
    for(; y; x = x * x % mod, y >>= 1)
      if(y & 1) ret = ret * x % mod;
    return ret;
}
//=============================================================
int main()
{
    N = read(), K = read();
    for(int i = 1; i <= N; i ++) X[i] = read(), Y[i] = read();
    for(int i = 1; i <= N; i ++)
    {
      ll tmp = 1ll;
      for(int j = 1; j <= N; j ++) 
        if(i != j) tmp = tmp * (X[i] + Mod - X[j]) % Mod;
      tmp = qpow(tmp, Mod - 2, Mod);
      for(int j = 1; j <= N; j ++) 
        if(i != j) tmp = tmp * (K + Mod - X[j]) % Mod;
      tmp = tmp * Y[i] % Mod, ans = (ans + tmp) % Mod;
    }
    printf("%lld", ans);
    system("pause");
    return 0;
}

自然幂数之和

CF622F The Sum of the k-th Powers

定义前 \(n\) 个自然数 \(k\) 次幂的和为:

\[S_k(n) = \sum_{i=1}^{n}i^k \]

给定 \(n,k\),求 \(S_k(n)\)
\(1\le n\le 10^9,\ 1\le k\le 10^6\)

性质

\(S_k(n)\) 为关于 \(n\)\(k+1\) 次多项式。

证明考虑对 \(k\) 进行归纳。
\(k=0\) 时,\(S_k(n) = n\),结论成立。
\(k=d\) 时:

通过二项式定理化简,提出一项相消,有:

\[\begin{aligned} &(i+1)^{d+1} - i^{d+1}\\ = &\sum_{j=0}^{d+1}\left(\begin{matrix}d+1\\j\end{matrix}\right)i^j - i^{d+1}\\ = &\sum_{j=0}^{d}\left(\begin{matrix}d+1\\j\end{matrix}\right)i^j \end{aligned}\]

\((i+1)^{d+1} - i^{d+1}\) 求和,有:

\[\begin{aligned}&\sum_{i=1}^{n}\left\{(i+1)^{d+1} - i^{d+1}\right\}\\=&\sum_{i=1}^{n}\sum_{j=0}^{d}\left(\begin{matrix}d+1\\j\end{matrix}\right)i^j\end{aligned} \]

发现 \(\left(\begin{matrix}d+1\\j\end{matrix}\right)\) 只与 \(j\) 有关,调换求和顺序,有:

\[\begin{aligned}&\sum_{i=1}^{n}\sum_{j=0}^{d}\left(\begin{matrix}d+1\\j\end{matrix}\right)i^j\\=&\sum_{j=0}^{d}\left(\begin{matrix}d+1\\j\end{matrix}\right)\sum_{i=1}^{n}i^j\\=&\sum_{j=0}^{d}\left(\begin{matrix}d+1\\j\end{matrix}\right)S_j(n)\end{aligned} \]

化出了有趣的玩意。
展开左侧 \(i^{d+1}\) 的求和式相消,则有:

\[(n+1)^{d+1} - 1 = \sum_{j=0}^{d}\left(\begin{matrix}d+1\\j\end{matrix}\right)S_j(n) \]

提出右侧 \(j=d\) 时的 \(S_d(n)\)

\[\left(\begin{matrix}d+1\\d\end{matrix}\right)S_d(n) = (n+1)^{d+1} -\sum_{j=0}^{d-1}\left\{\left(\begin{matrix}d+1\\j\end{matrix}\right)S_j(n)\right\} - 1 \]

二项式定理展开右侧,并略做处理:

\[S_d(n) = \dfrac{1}{d+1}\left\{\sum_{j=0}^{d+1}n^j -\sum_{j=0}^{d-1}\left\{\left(\begin{matrix}d+1\\j\end{matrix}\right)S_j(n)\right\} - 1\right\} \]

对于右侧,通过数学归纳得到 \(S_j(n)\) 为关于 \(n\)\(j\) 次多项式。
则整个式子最高次项出现在 \(\sum\limits_{j=0}^{d+1}n^j\) 中,次数为 \(d+1\)
得证。

具体实现

由性质,\(S_k(n)\) 为关于 \(n\)\(k + 1\) 次多项式,需要 \(k+2\) 个点值进行构造。
可以直接取 \(k+1\) 个点,按照定义计算出 \(S_k(n)\),构造多项式。
复杂度 \(O(k^2)\)?我觉得不行。

题目对选择的点并没有要求,考虑选取一段连续的点进行优化。
使 \(x_i = i\),则选择的点集为 \(\{(i, S_k(i))\, \mid \, i \in \N^+, i\le k+2\}\)
按照定义递推出来即可,复杂度 \(O(k \log k)\)

考虑要求的点 \((n, S_k(n))\),将其代入插值公式,有:

\[S_k(n)= \sum_{i=1}^{k+2}S_k(i)\prod_{i\not ={j}}\dfrac{n-j}{i-j} \]

发现乘积项的分母与 \(n\) 无关,提出来:

\[S_k(n)= \sum_{i=1}^{k+2}S_k(i)\dfrac{\prod\limits_{i\not ={j}}{(n-j)}}{\prod\limits_{i\not ={j}}(i-j)} \]

手玩一下乘积项的变化规律。

对于分母,\(j\le k +2\),在已知 \(i\) 时有:

\[\begin{aligned} &\dfrac{1}{\prod\limits_{i\not ={j}}(i-j)}\\ = &\dfrac{1}{i(i-1)(i-2) \dots 1 \times(-1) \dots (k+2-i-1)(k+2-i)}\\=&(-1)^{k+2-i}\dfrac{1}{i! (k+2-i)!} \end{aligned}\]

可先预处理阶乘再求逆元,快速得到分母的值。

对于分子,也暴力拆一波:

\[\begin{aligned} &\prod\limits_{i\not ={j}}{(n-j)}\\ =& n(n-1) \dots (n-(i-1))(n-(i+1))\dots (n-(k+2))\\ =&(\prod_{j=1}^{i-1}(n-j)) (\prod_{j=i+1}^{k+2}(n-j))& \end{aligned}\]

考虑维护 \((n-j)\) 的 前/后 缀积,可快速得到分子的值。

则有:

\[S_k(n)= \sum_{i=1}^{k+2}(-1)^{k+2-i}S_k(i)\dfrac{(\prod\limits_{j=1}^{i-1}(n-j)) (\prod\limits_{j=i+1}^{k+2}(n-j))}{i! (k+2-i)!} \]

复杂度

\(O(k \log k)\) 处理 \(k+2\) 个连续点值。
\(O(n)\) 预处理前/后 缀积,阶乘。
求逆元时可以 \(O(n)\) 预处理,也可以单次 \(O(\log k)\)

总复杂度 \(O(k\log k)\)

优缺点

优点:简单好写,复杂度优秀,就感觉到快。

缺点:出现了除法,需要求逆元,需保证模数为质数。

若模数不为质数怎么办? 可利用 斯特林数/伯努利数 求解。
伯努利数求解方法 详见 皎月半洒花的题解

斯特林数求解方法 复杂度 \(O(k^2)\),无法通过本题,详见 自然数幂之和 - Luckyblock


写在最后

参考资料:
拉格朗日插值 - OI Wiki

posted @ 2020-08-02 15:00  Luckyblock  阅读(336)  评论(0编辑  收藏  举报