拉格朗日插值的一种做法

· · 算法·理论

背景

个人认为短除法插多项式比较麻烦，因此自创方法进行同样严格 O(n^2+nm) 的不连续多点插值。

问题

有 n 个已知的点，要求 m 个 x 坐标在经过这 n 个点的 n-1 次多项式上对应的值。

分析

那怎样快速得到多个点的点值呢？下面记已知的点数为 n，要插出的点数为 m。

我们知道，已知 n 个点 (a_i,b_i)，一个点 k 的对应值 f(k) 可以这样插值求：

f(k) = \sum^n_{i = 1} b_i \prod^n_{j \neq i} \dfrac{k - a_j}{a_i - a_j}

:::info[求单点]

int n,k;
cin>>n>>k;
for(int i = 1; i <= n; i ++) cin>>a[i]>>b[i];
LL res = 0;
for(int i = 1; i <= n; i ++)
{
    LL x = 1,y = 1;
    for(int j = 1; j <= n; j ++)
        if(i != j) x = x * (k - a[j]) % mod,y = y * (a[i] - a[j]) % mod;
    res = (res + b[i] * x % mod * qmi(y,mod - 2) % mod) % mod;
}
cout<<(res % mod + mod) % mod<<'\n';

::: 这样单点插值复杂度是 O(n^2) 的。如果直接这样算 m 个点的值，O(n^2m) 显然不可行。

可以先把分母的常数提出来，像下面这样： :::info[分母]

for(int i = 1; i <= n; i ++) f[i] = 1; 
for(int i = 1; i <= n; i ++)
{
    for(int j = 1; j <= n; j ++)
        if(i != j) f[i] = f[i] * (a[i] - a[j]) % mod;
    f[i] = inv(f[i]);
}

::: 然后因为 \prod\limits^n_{j\neq i}(k-a_j)=\dfrac{\prod\limits^n_{j=1} (k-a_j)}{k-a_i}，所以上面的也是可以先算出来的。 :::info[分子的一部分]

for(int k = 1; k <= m; k ++)
    for(int i = 1; i <= n; i ++) s[k] = s[k] * (q[k] - a[i]) % mod;

::: 然后就可以计算答案了。逆元离线预处理后复杂度 O(n^2+nm)。 :::info[计算答案]

for(int k = 1; k <= m; k ++)
    for(int i = 1; i <= n; i ++) res[k] = (res[k] + s[k] * inv(q[k] - a[i]) % mod * f[i] % mod * b[i] % mod) % mod;

:::

注意这种方法计算时不要对已知的点插值，否则会计算 0 的逆元。

实现

#include<bits/stdc++.h>
using namespace std;
const int mod = 1e9 + 7;
typedef long long LL;
typedef vector<LL> vec;
int n;
LL qmi(LL a,LL k)
{
    LL res = 1;
    while(k)
    {
        if(k & 1) res = res * a % mod;
        k >>= 1;
        a = a * a % mod;
    }
    return res;
}
vec get_inv(vec a)
{
    vec b,c = a,res;
    b.resize(a.size());
    res.resize(a.size());
    LL s = 1;
    for(auto x : a) s = s * x;
    for(int i = 1; i < a.size(); i ++) a[i] = a[i - 1] * a[i] % mod;
    b.back() = qmi(a.back(),mod - 2);
    for(int i = b.size() - 2; i >= 0; i --) b[i] = b[i + 1] * c[i + 1] % mod;
    res[0] = 1;
    for(int i = 1; i < a.size(); i ++) res[i] = a[i - 1] * b[i] % mod;
    for(int i = 0; i < c.size(); i ++) assert(c[i] * res[i] % mod == 1);
    return res;
}
vec get(const vec &a,const vec &b,vec q)
{
    vec res,s,f,q0,ans,inv;
    unordered_map<LL,int> st;
    for(int i = 0; i < a.size(); i ++) st[a[i]] = i;
    for(auto x : q)
        if(!st.count(x)) q0.push_back(x);
    swap(q,q0);
    int n = a.size(),m = q.size(),idx = 1;
    f.resize(n,1);
    res.resize(m,0);
    s.resize(m,1);
    inv.push_back(1);
    for(int i = 0; i < n; i ++)
        for(int j = 0; j < n; j ++)
            if(i != j) inv.push_back((a[i] - a[j] + mod) % mod);
    for(int k = 0; k < m; k ++)
        for(int i = 0; i < n; i ++) inv.push_back((q[k] - a[i] + mod) % mod);
    inv = get_inv(inv);
    for(int i = 0; i < n; i ++)
        for(int j = 0; j < n; j ++)
            if(i != j) f[i] = f[i] * inv[idx ++] % mod;
    for(int k = 0; k < m; k ++)
        for(int i = 0; i < n; i ++) s[k] = s[k] * (q[k] - a[i]) % mod;
    for(int k = 0; k < m; k ++)
        for(int i = 0; i < n; i ++) res[k] = (res[k] + s[k] * inv[idx ++] % mod * f[i] % mod * b[i] % mod) % mod;
    for(auto &x : res) x = (x + mod) % mod;
    for(int i = 0,j = 0; i < q0.size(); i ++)
    {
        if(st.count(q0[i])) ans.push_back(b[st[q0[i]]]);
        else ans.push_back(res[j ++]);
    }
    return ans;
}

例题

P9165 「INOH」Round 1 - 意外

P10832 [COTS 2023] 传 Mapa