算法学习之：拉格朗日插值

最新推荐文章于 2025-03-27 19:40:38 发布

lvzelong2014

最新推荐文章于 2025-03-27 19:40:38 发布

阅读量5.1k

点赞数 1

分类专栏： # 数学相关

本文链接：https://blog.csdn.net/lvzelong2014/article/details/79159346

版权

数学相关专栏收录该内容

64 篇文章

订阅专栏

本文详细介绍了拉格朗日插值法，包括定义、伪代码和应用实例。通过对给定点进行插值，可以求得k次多项式函数。在OI竞赛中，可以通过预处理降低时间复杂度至O(k)。文章还提到了使用差分方法作为解决此类问题的小技巧。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

算法学习之：拉格朗日插值

定义

对某个k次多项式 $f^k(n)$ 函数，已知有给定的k + 1个点：
$(x_{0},y_{0}),\ldots ,(x_{k},y_{k})$
假设任意两个不同的 $x_j$ 都互不相同，那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：
$L(x)=\sum\limits _{{j=0}}^{{k}}y_{j}\ell _{j}(x)$
其中每个 $\ell _{j}(x)$ 为拉格朗日基本多项式（或称插值基函数），其表达式为：
$\ell _{j}(x):=\prod\limits _{{i=0,\,i\neq j}}^{{k}}{\frac {x-x_{i}}{x_{j}-x_{i}}}={\frac {(x-x_{0})}{(x_{j}-x_{0})}}\cdots {\frac {(x-x_{{j-1}})}{(x_{j}-x_{{j-1}})}}{\frac {(x-x_{{j+1}})}{(x_{j}-x_{{j+1}})}}\cdots {\frac {(x-x_{{k}})}{(x_{j}-x_{{k}})}}.$
拉格朗日基本多项式 $\ell _{j}(x)$ 的特点是在 $x_{j}$ 上取值为1，在其它的点 $x_{i},\,i\neq j$ 上取值为0。
这样对于每个L(x)，它必定过所有 $(x_{0},y_{0}),\ldots ,(x_{k},y_{k})$

伪代码

double Lagrange(double *ax, double *ay, int n, int x) {
    double ret = 0;
    for(int i = 1;i <= n; ++i) {
        double t = 1;
        for(int j = 0;j <= n; ++j)
            if(i != j)
            t = t * (x- ax[j]) / (ax[i] - ax[j]);
       }
       ret += t;
   }
   return ret;
}

应用

在实际的OI竞赛中，选手常常要自己带值插值，所以我们一般选择 $(0,f^k(0)),\ldots ,(k,f^k(k))$ ，带入L(x)有
$f^k(n)=L(n)=\sum _{j=0}^{k+1}\frac{f^k(j)n(n−1)\dots(n−j+1)(n−j−1)\dots(n−k−1)(−1)^{k+1−j}}{j!(k+1−j)!}$
于是我们可以预处理阶乘等来使得时间复杂度变成O(k)的
缺点就是对模数有要求。

伪代码

long long Lagrange(long long *a, int n, long long pos) {
    if(pos <= n) return a[pos];
    s1[0] = s2[n + 1] = 1;
    for(int i = 1;i <= n; ++i) s1[i] = s1[i - 1] * (pos - i) % p;
    for(int i = n; i; --i) s2[i] = s2[i + 1] * (pos - i) % p;
    long long ans = 0;
    for(int i = 1;i <= n; ++i) 
        up(ans, a[i] * s1[i - 1] % p * s2[i + 1] % p * inv[i - 1] % p * 
        inv[n - i] % p * ((n - i) & 1 ? -1 : 1));
    return ans;
}

例题

题目链接

bzoj3453: tyvj 1858 XLkxc

分析

$f(n)=\sum\limits_{i=0}^n\sum\limits_{j=1}^{a+id}\sum\limits_{l=1}^jl^k$
我们令 $g(n)=\sum\limits_{j=1}^{n}\sum\limits_{l=1}^jl^k$
那么有 $f(n)=\sum\limits_{i=0}^ng(a+id)$
有一个神奇的引理

引理： $f^k(n)=\sum\limits_{i=1}^ni^k$ 必定能用k+1次多项式函数表示

那么g(n)就是k+2次多项式，f(n)就是k+4次多项式
然而我不会。

小技巧

有一种神奇的方法，叫做差分。
首先把值带入，然后差分，差分第k次为零，那么它就是k-1次多项式。

for(int i = 1;i <= k + 10; ++i) g[i] = mpower(i, k);
        for(int i = 2;i <= k + 10; ++i) up(g[i], g[i - 1]);
        for(int i = 2;i <= k + 10; ++i) up(g[i], g[i - 1]);
        for(int i = 1;i <= k + 10; ++i) cout<<g[i]<<" "; cout<<endl;
        for(int i = 1;i <= k + 10; ++i) {
            for(int j = 1;j <= k + 10 - i; ++j) g[j] = g[j + 1] - g[j];
            for(int j = 1;j <= k + 10 - i; ++j) cout<<g[j]<<" "; cout<<endl;
        }
/*运行结果
3 3 3 3
1 10 46 146 371 812 1596 2892 4917 7942 12298 18382 26663
9 36 100 225 441 784 1296 2025 3025 4356 6084 8281
27 64 125 216 343 512 729 1000 1331 1728 2197
37 61 91 127 169 217 271 331 397 469
24 30 36 42 48 54 60 66 72
6 6 6 6 6 6 6 6
0 0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0
0 0 0 0
0 0 0
0 0
0
*/

对于f(n)也是一样的操作
然后我们就可以g(n)插值一次，f(n)插值一次就可以啦。

代码

/**************************************************************
    Problem: 3453
    User: 2014lvzelong
    Language: C++
    Result: Accepted
    Time:32 ms
    Memory:1300 kb
****************************************************************/

#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<map>
#include<cmath>
using namespace std;
const int N = 155, p = 1234567891;
long long read() {
    char ch = getchar(); long long x = 0, f = 1;
    while(ch < '0' || ch > '9') {if(ch == '-') f = -1; ch = getchar();}
    while(ch >= '0' && ch <= '9') {x = x * 10 + ch - '0'; ch = getchar();}
    return x * f;
}
int T, k;
long long inv[N << 1], s1[N], s2[N], g[N], f[N], a, d, n, m;
void init(int n) {
    inv[0] = inv[1] = 1;
    for(int i = 2;i <= n; ++i) inv[i] = inv[p % i] * (p - p / i) % p;
    for(int i = 2;i <= n; ++i) inv[i] = inv[i] * inv[i - 1]  % p;
}
void up(long long &a, long long b) {(a += b) %= p;}
long long Lagrange(long long *a, int n, long long pos) {
    if(pos <= n) return a[pos];
    s1[0] = s2[n + 1] = 1;
    for(int i = 1;i <= n; ++i) s1[i] = s1[i - 1] * (pos - i) % p;
    for(int i = n; i; --i) s2[i] = s2[i + 1] * (pos - i) % p;
    long long ans = 0;
    for(int i = 1;i <= n; ++i) 
        up(ans, a[i] * s1[i - 1] % p * s2[i + 1] % p * inv[i - 1] % p * 
        inv[n - i] % p * ((n - i) & 1 ? -1 : 1));
    return ans;
}
long long mpower(long long x, long long y) {
    long long ans = 1;
    for(;y; y >>= 1, (x *= x) %= p)
        if(y & 1) (ans *= x) %= p;
    return ans;
}

int main() {
    T = read();
    init(300);
    while(T--) {
        k = read(); a = read(); n = read(); d = read();
        for(int i = 1;i <= k + 3; ++i) g[i] = mpower(i, k);
        for(int i = 2;i <= k + 3; ++i) up(g[i], g[i - 1]);
        for(int i = 2;i <= k + 3; ++i) up(g[i], g[i - 1]);
        f[0] = Lagrange(g, k + 3, a);
        for(int i = 1;i <= k + 5; ++i) 
            f[i] = Lagrange(g, k + 3, (i * d + a) % p), up(f[i], f[i - 1]);
        printf("%lld\n", (Lagrange(f, k + 5, n) + p) % p);
    }
    return 0;
}