Polynomial Interpolation (Lagrange Interpolation)
(poly_interpolation.hpp)

• View this file on GitHub
• Last update: 2025-01-19 15:28:01+08:00
• Include: #include "poly_interpolation.hpp"

Lagrange Interpolation

Given sample points $f(0), f(1), \dots ,f(n - 1)$ of polynomial $f \in \mathbb{C}\left\lbrack x\right\rbrack$ with $\deg f \lt n$, we want to compute $f(c)$ for $c \in \mathbb{Z}$.

Recall the Lagrange interpolation formula,

\[f(x) = \sum _ {0 \leq i \lt n}\left(f(i) \prod _ {0 \leq j \lt n \atop j \neq i} \frac{x - j}{i - j}\right)\]

and let $k_i := \frac{f(i)}{i!} \prod _ {0 \leq j \lt i}(x - j)$ for some $i$, we have

\[f(x) = \sum _ {0 \leq i \lt n}\left(k_i \frac{(-1)^{n - 1 - i}}{(n - 1 - i)!} \prod _ {i \lt j \lt n}(x - j)\right).\]

C-Finite sequence

The O.g.f. of $\left(f(j)\right) _ {j \geq 0}$ is

\[F(x) = \sum _ {j \geq 0} f(j) x^j\]

Let’s introduce the difference operator $\Delta : f(x) \mapsto f(x+1) - f(x)$, note that $\deg \left(\Delta f(x)\right) \lt n-1$, so $\Delta^n f(x)=0$. Now we have

\[\begin{aligned} (1-x)F(x)-f(0) &= x\sum _ {j\geq 0}\left(\Delta f(j)\right)x^j \\ (1-x)\left((1-x)F(x)-f(0)\right)-x\Delta f(0) &= x^2 \sum _ {j\geq 0}\left(\Delta^2 f(j)\right)x^j \\ \vdots \end{aligned}\]

since $\Delta^n f(x)=0$, we know $\deg\left(\left(1-x\right)^n F(x)\right) \lt n$, which is saying that

\[F(x)=\frac{P(x)}{(1-x)^n},\qquad P\in\mathbb{C}\left\lbrack x\right\rbrack _ {\lt n}.\]

Depends on

• batch_inv.hpp
• binomial.hpp

Verified with

• test/enumerative_combinatorics/factorial.0.test.cpp

Code

#pragma once

#include "batch_inv.hpp"
#include "binomial.hpp"
#include <cassert>
#include <vector>

// returns f
// x: x_0,...
// y: f(x_0),...
template<typename Tp> inline std::vector<Tp>
lagrange_interpolation_naive(const std::vector<Tp> &x, const std::vector<Tp> &y) {
    assert(x.size() == y.size());
    const int n = x.size();
    std::vector<Tp> M(n + 1), xx(n), f(n);
    M[0] = 1;
    for (int i = 0; i < n; ++i)
        for (int j = i; j >= 0; --j) M[j + 1] += M[j], M[j] *= -x[i];
    for (int i = n - 1; i >= 0; --i)
        for (int j = 0; j < n; ++j) xx[j] = xx[j] * x[j] + M[i + 1] * (i + 1);
    xx = batch_inv(xx);
    for (int i = 0; i < n; ++i) {
        Tp t = y[i] * xx[i], k = M[n];
        for (int j = n - 1; j >= 0; --j) f[j] += k * t, k = M[j] + k * x[i];
    }
    return f;
}

// returns f(c)
// f: polynomial f -> f[0]=f(0),...
template<typename Tp> inline Tp lagrange_interpolation_iota(const std::vector<Tp> &f, Tp c) {
    if (f.empty()) return 0;
    const int n = f.size();
    auto &&bin  = Binomial<Tp>::get(n);
    std::vector<Tp> k(n);
    k[0] = f[0];
    Tp v = 1;
    for (int i = 1; i < n; ++i) k[i] = f[i] * (v *= (c - (i - 1)) * bin.inv(i));
    Tp res = k[n - 1];
    v      = 1;
    for (int i = n - 2; i >= 0; --i) res += k[i] * (v *= -(c - (i + 1)) * bin.inv(n - 1 - i));
    return res;
}

#line 2 "poly_interpolation.hpp"

#line 2 "batch_inv.hpp"

#include <cassert>
#include <vector>

template<typename Tp> inline std::vector<Tp> batch_inv(const std::vector<Tp> &a) {
    if (a.empty()) return {};
    const int n = a.size();
    std::vector<Tp> b(n);
    Tp v = 1;
    for (int i = 0; i < n; ++i) b[i] = v, v *= a[i];
    assert(v != 0);
    v = v.inv();
    for (int i = n - 1; i >= 0; --i) b[i] *= v, v *= a[i];
    return b;
}
#line 2 "binomial.hpp"

#include <algorithm>
#line 5 "binomial.hpp"

template<typename Tp> class Binomial {
    std::vector<Tp> factorial_, invfactorial_;

    Binomial() : factorial_{Tp(1)}, invfactorial_{Tp(1)} {}

    void preprocess(int n) {
        if (const int nn = factorial_.size(); nn < n) {
            int k = nn;
            while (k < n) k *= 2;
            k = std::min<long long>(k, Tp::mod());
            factorial_.resize(k);
            invfactorial_.resize(k);
            for (int i = nn; i < k; ++i) factorial_[i] = factorial_[i - 1] * i;
            invfactorial_.back() = factorial_.back().inv();
            for (int i = k - 2; i >= nn; --i) invfactorial_[i] = invfactorial_[i + 1] * (i + 1);
        }
    }

public:
    static const Binomial &get(int n) {
        static Binomial bin;
        bin.preprocess(n);
        return bin;
    }

    Tp binom(int n, int m) const {
        return n < m ? Tp() : factorial_[n] * invfactorial_[m] * invfactorial_[n - m];
    }
    Tp inv(int n) const { return factorial_[n - 1] * invfactorial_[n]; }
    Tp factorial(int n) const { return factorial_[n]; }
    Tp inv_factorial(int n) const { return invfactorial_[n]; }
};
#line 7 "poly_interpolation.hpp"

// returns f
// x: x_0,...
// y: f(x_0),...
template<typename Tp> inline std::vector<Tp>
lagrange_interpolation_naive(const std::vector<Tp> &x, const std::vector<Tp> &y) {
    assert(x.size() == y.size());
    const int n = x.size();
    std::vector<Tp> M(n + 1), xx(n), f(n);
    M[0] = 1;
    for (int i = 0; i < n; ++i)
        for (int j = i; j >= 0; --j) M[j + 1] += M[j], M[j] *= -x[i];
    for (int i = n - 1; i >= 0; --i)
        for (int j = 0; j < n; ++j) xx[j] = xx[j] * x[j] + M[i + 1] * (i + 1);
    xx = batch_inv(xx);
    for (int i = 0; i < n; ++i) {
        Tp t = y[i] * xx[i], k = M[n];
        for (int j = n - 1; j >= 0; --j) f[j] += k * t, k = M[j] + k * x[i];
    }
    return f;
}

// returns f(c)
// f: polynomial f -> f[0]=f(0),...
template<typename Tp> inline Tp lagrange_interpolation_iota(const std::vector<Tp> &f, Tp c) {
    if (f.empty()) return 0;
    const int n = f.size();
    auto &&bin  = Binomial<Tp>::get(n);
    std::vector<Tp> k(n);
    k[0] = f[0];
    Tp v = 1;
    for (int i = 1; i < n; ++i) k[i] = f[i] * (v *= (c - (i - 1)) * bin.inv(i));
    Tp res = k[n - 1];
    v      = 1;
    for (int i = n - 2; i >= 0; --i) res += k[i] * (v *= -(c - (i + 1)) * bin.inv(n - 1 - i));
    return res;
}

Back to top page