Skip to content

Interpolation

Introduction

Interpolation is a method of estimating new data points within the range of a discrete set of known data points. Interpolation is often used in function fitting.

For example, given the data points:

\(x\) \(0\) \(1\) \(2\) \(3\) \(4\) \(5\) \(6\)
\(f(x)\) \(0\) \(0.8415\) \(0.9093\) \(0.1411\) \(-0.7568\) \(-0.9589\) \(-0.2794\)

where \(f(x)\) is unknown, interpolation can estimate unknown data points by fitting \(f(x)\) in a certain form.

For instance, we can fit \(f(x)\) with a piecewise linear function:

This interpolation method is called linear interpolation.

We can also fit \(f(x)\) with a polynomial:

This interpolation method is called polynomial interpolation.

The general form of polynomial interpolation is as follows:

Polynomial Interpolation

Given \(n+1\) known points \((x_0,y_0),(x_1,y_1),\dots,(x_n,y_n)\), find a polynomial \(f(x)=\sum_{i=0}^n a_ix^i\) that satisfies

\[ f(x_i)=y_i,\qquad\forall i=0,1,\dots,n \]

Below, two methods for polynomial interpolation are introduced: Lagrange interpolation and Newton interpolation. It can be shown that the results obtained by these two methods are equivalent.

Lagrange Interpolation

We need to construct a function \(f(x)\) passing through points \(P_1(x_1, y_1), P_2(x_2,y_2),\cdots,P_n(x_n,y_n)\). First, let the projection of the \(i\)-th point onto the \(x\)-axis be \(P_i^{\prime}(x_i,0)\).

Consider constructing \(n\) functions \(f_1(x), f_2(x), \cdots, f_n(x)\), such that for the \(i\)-th function \(f_i(x)\), its graph passes through \(\begin{cases}P_j^{\prime}(x_j,0),(j\neq i)\\P_i(x_i,y_i)\end{cases}\). Then the desired function \(f(x)=\sum\limits_{i=1}^nf_i(x)\).

We can set \(f_i(x)=a\cdot\prod_{j\neq i}(x-x_j)\). Substituting point \(P_i(x_i,y_i)\) gives \(a=\dfrac{y_i}{\prod_{j\neq i} (x_i-x_j)}\), so

\[ f_i(x)=y_i\cdot\dfrac{\prod_{j\neq i} (x-x_j)}{\prod_{j\neq i} (x_i-x_j)}=y_i\cdot\prod_{j\neq i}\dfrac{x-x_j}{x_i-x_j} \]

Thus, the Lagrange interpolation form is:

\[ f(x)=\sum_{i=1}^ny_i\cdot\prod_{j\neq i}\dfrac{x-x_j}{x_i-x_j} \]

A naive implementation has time complexity \(O(n^2)\), which can be optimized to \(O(n\log^2 n)\). See Polynomial Fast Interpolation.

Luogu P4781【Template】Lagrange Interpolation

Given \(n\) point pairs \((x_i,y_i)\) and \(k\), with \(\forall i,j\) satisfying \(i\neq j \iff x_i\neq x_j\), \(f(x_i)\equiv y_i\pmod{998244353}\), and \(\deg(f(x)) < n\) (defining \(\deg(0)=-\infty\)), find \(f(k)\bmod{998244353}\).

Solution

In this problem, we only need to compute the value of \(f(k)\), so we can directly substitute \(k\) into the formula. Sometimes, multiple evaluations or more complex operations are required, in which case we need to find all coefficients of \(f\). The code provides an implementation for obtaining the coefficients.

\[ f(k)=\sum_{i=1}^{n}y_i\prod_{j\neq i }\frac{k-x_j}{x_i-x_j} \]

In this problem, modular inverses also need to be computed. If we first compute the numerator and denominator separately, then multiply the numerator by the modular inverse of the denominator and accumulate into the final answer, the bottleneck of time complexity will not be in computing inverses, resulting in \(O(n^2)\) time complexity.

Since operations are under the fixed modulus \(998244353\), we consider the time complexity of computing modular inverses as constant here.

Code Implementation 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
#include <iostream>
#include <vector>

constexpr int MOD = 998244353;
using LL = long long;

int inv(int k) {
  int res = 1;
  for (int e = MOD - 2; e; e /= 2) {
    if (e & 1) res = (LL)res * k % MOD;
    k = (LL)k * k % MOD;
  }
  return res;
}

// Returns f such that f(x_i) = y_i
// Without considering the time for modular inverse, complexity is O(n^2)
std::vector<int> lagrange_interpolation(const std::vector<int> &x,
                                        const std::vector<int> &y) {
  const int n = x.size();
  std::vector<int> M(n + 1), px(n, 1), f(n);
  M[0] = 1;
  // Compute M(x) = prod_(i=0..n-1)(x - x_i)
  for (int i = 0; i < n; ++i) {
    for (int j = i; j >= 0; --j) {
      M[j + 1] = (M[j] + M[j + 1]) % MOD;
      M[j] = (LL)M[j] * (MOD - x[i]) % MOD;
    }
  }
  // Compute px_i = prod_(j=0..n-1, j!=i) (x_i - x_j)
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < n; ++j)
      if (i != j) {
        px[i] = (LL)px[i] * (x[i] - x[j] + MOD) % MOD;
      }
  }
  // Combine to get f(x) = sum_(i=0..n-1)(y_i / px_i)(M(x) / (x - x_i))
  for (int i = 0; i < n; ++i) {
    LL t = (LL)y[i] * inv(px[i]) % MOD, k = M[n];
    for (int j = n - 1; j >= 0; --j) {
      f[j] = (f[j] + k * t) % MOD;
      k = (M[j] + k * x[i]) % MOD;
    }
  }
  return f;
}

int main() {
  std::ios::sync_with_stdio(false);
  std::cin.tie(nullptr);
  int n, k;
  std::cin >> n >> k;
  std::vector<int> x(n), y(n);
  for (int i = 0; i < n; ++i) std::cin >> x[i] >> y[i];
  const auto f = lagrange_interpolation(x, y);
  int v = 0;
  for (int i = n - 1; i >= 0; --i) v = ((LL)v * k + f[i]) % MOD;
  std::cout << v << '\n';
  return 0;
}

Lagrange Interpolation with Consecutive Integer x-coordinates

If the x-coordinates of the known points are consecutive integers, we can achieve \(O(n)\) interpolation.

Let the desired polynomial be \(f(x)\), and we know \(f(1),\cdots,f(n+1)\) (\(1\le i\le n+1\)). Substituting into the interpolation formula:

\[ \begin{aligned} f(x)&=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j\ne i}\frac{x-x_j}{x_i-x_j}\\ &=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j\ne i}\frac{x-j}{i-j} \end{aligned} \]

The product can be considered separately for numerator and denominator. The numerator is:

\[ \dfrac{\prod\limits_{j=1}^{n+1}(x-j)}{x-i} \]

The denominator's product of \(i-j\) can be split into two factorials:

\[ (-1)^{n+1-i}\cdot(i-1)!\cdot(n+1-i)! \]

Thus, the interpolation formula for x-coordinates \(1,\cdots,n+1\) is:

\[ f(x)=\sum\limits_{i=1}^{n+1}(-1)^{n+1-i}y_i\cdot\frac{\prod\limits_{j=1}^{n+1}(x-j)}{(i-1)!(n+1-i)!(x-i)} \]

Preprocess the prefix and suffix products of \((x-i)\), factorials and factorial inverses, then substitute into this formula. Complexity is \(O(n)\).

Example CF622F The Sum of the k-th Powers

Given \(n,k\), find \(\sum\limits_{i=1}^ni^k\) modulo \(10^9+7\).

Solution

In this problem, the answer is a polynomial of degree \(k+1\). Therefore, we can linearly sieve the values of \(1^i,\cdots,(k+2)^i\) and then perform \(O(n)\) interpolation.

Alternatively, the following formula can be derived from the difference method using combinatorial mathematics:

\[ f(x)=\sum_{i=1}^{n+1}\binom{x-1}{i-1}\sum_{j=1}^{i}(-1)^{i+j}\binom{i-1}{j-1}y_{j}=\sum\limits_{i=1}^{n+1}y_i\cdot\frac{\prod\limits_{j=1}^{n+1}(x-j)}{(x-i)\cdot(-1)^{n+1-i}\cdot(i-1)!\cdot(n+1-i)!} \]
Code Implementation 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
// By: Luogu@rui_er(122461)
#include <iostream>
using namespace std;
constexpr int N = 1e6 + 5, mod = 1e9 + 7;

int n, k, tab[N], p[N], pcnt, f[N], pre[N], suf[N], fac[N], inv[N], ans;

int qpow(int x, int y) {
  int ans = 1;
  for (; y; y >>= 1, x = 1LL * x * x % mod)
    if (y & 1) ans = 1LL * ans * x % mod;
  return ans;
}

void sieve(int lim) {
  f[1] = 1;
  for (int i = 2; i <= lim; i++) {
    if (!tab[i]) {
      p[++pcnt] = i;
      f[i] = qpow(i, k);
    }
    for (int j = 1; j <= pcnt && 1LL * i * p[j] <= lim; j++) {
      tab[i * p[j]] = 1;
      f[i * p[j]] = 1LL * f[i] * f[p[j]] % mod;
      if (!(i % p[j])) break;
    }
  }
  for (int i = 2; i <= lim; i++) f[i] = (f[i - 1] + f[i]) % mod;
}

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  cin >> n >> k;
  sieve(k + 2);
  if (n <= k + 2) return cout << f[n], 0;
  pre[0] = suf[k + 3] = 1;
  for (int i = 1; i <= k + 2; i++) pre[i] = 1LL * pre[i - 1] * (n - i) % mod;
  for (int i = k + 2; i >= 1; i--) suf[i] = 1LL * suf[i + 1] * (n - i) % mod;
  fac[0] = inv[0] = fac[1] = inv[1] = 1;
  for (int i = 2; i <= k + 2; i++) {
    fac[i] = 1LL * fac[i - 1] * i % mod;
    inv[i] = 1LL * (mod - mod / i) * inv[mod % i] % mod;
  }
  for (int i = 2; i <= k + 2; i++) inv[i] = 1LL * inv[i - 1] * inv[i] % mod;
  for (int i = 1; i <= k + 2; i++) {
    int P = 1LL * pre[i - 1] * suf[i + 1] % mod;
    int Q = 1LL * inv[i - 1] * inv[k + 2 - i] % mod;
    int mul = ((k + 2 - i) & 1) ? -1 : 1;
    ans = (ans + 1LL * (Q * mul + mod) % mod * P % mod * f[i] % mod) % mod;
  }
  cout << ans << '\n';
  return 0;
}

Newton Interpolation

Newton interpolation is a method based on higher-order differences, with the advantage of supporting \(O(n)\) insertion of new data points.

To achieve \(O(n)\) insertion of new data points, we let:

\[ f(x)=\sum_{j=0}^n a_jn_j(x) \]

where \(n_j(x):=\prod_{i=0}^{j-1}(x-x_i)\) is called the Newton basis.

If we solve for \(a_j\), we obtain the interpolation polynomial \(f(x)\). We define the forward divided differences as follows:

\[ \begin{aligned} \lbrack y_k\rbrack & := y_k, & k=0,\dots,n, \\ [y_k,\dots,y_{k+j}] & := \dfrac{[y_{k+1},\dots,y_{k+j}]-[y_k,\dots,y_{k+j-1}]}{x_{k+j}-x_k}, & k=0,\dots,n-j,~j=1,\dots,n. \end{aligned} \]

Then:

\[ \begin{aligned} f(x)&=[y_0]+[y_0,y_1](x-x_0)+\dots+[y_0,\dots,y_n](x-x_0)\dots(x-x_{n-1})\\ &=\sum_{j=0}^n [y_0,\dots,y_j]n_j(x) \end{aligned} \]

This is the Newton interpolation form. A naive implementation has time complexity \(O(n^2)\).

If the sample points are equally spaced (i.e., \(x_i=x_0+ih\), \(i=1,\dots,n\)), we can derive

\[ [y_k,\dots,y_{k+j}]=\frac{1}{j!h^j}\Delta^{(j)}y_k, \]

where \(\Delta^{(j)}y_k\) are the forward differences, defined as:

\[ \begin{aligned} \Delta^{(0)}y_k & := y_k, & k=0,\dots,n, \\ \Delta^{(j)}y_k & := \Delta^{(j-1)} y_{k+1}-\Delta^{(j-1)} y_k, & k=0,\dots,n-j,~j=1,\dots,n. \end{aligned} \]

Let \(x=x_0+sh\). Then the Newton interpolation formula becomes

\[ f(x)=\sum_{j=0}^n \binom{s}{j}j!h^j[y_0,\dots,y_j]=\sum_{j=0}^n \binom{s}{j}\Delta^{(j)}y_0. \]
Code Implementation (Luogu P4781【Template】Lagrange Interpolation) 
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
#include <cstdint>
#include <iostream>
#include <vector>
using namespace std;

constexpr uint32_t MOD = 998244353;

struct mint {
  uint32_t v_;

  mint() : v_(0) {}

  mint(int64_t v) {
    int64_t x = v % (int64_t)MOD;
    v_ = (uint32_t)(x + (x < 0 ? MOD : 0));
  }

  friend mint inv(mint const &x) {
    int64_t a = x.v_, b = MOD;
    if ((a %= b) == 0) return 0;
    int64_t s = b, m0 = 0;
    for (int64_t q = 0, _ = 0, m1 = 1; a;) {
      _ = s - a * (q = s / a);
      s = a;
      a = _;
      _ = m0 - m1 * q;
      m0 = m1;
      m1 = _;
    }
    return m0;
  }

  mint &operator+=(mint const &r) {
    if ((v_ += r.v_) >= MOD) v_ -= MOD;
    return *this;
  }

  mint &operator-=(mint const &r) {
    if ((v_ -= r.v_) >= MOD) v_ += MOD;
    return *this;
  }

  mint &operator*=(mint const &r) {
    v_ = (uint32_t)((uint64_t)v_ * r.v_ % MOD);
    return *this;
  }

  mint &operator/=(mint const &r) { return *this = *this * inv(r); }

  friend mint operator+(mint l, mint const &r) { return l += r; }

  friend mint operator-(mint l, mint const &r) { return l -= r; }

  friend mint operator*(mint l, mint const &r) { return l *= r; }

  friend mint operator/(mint l, mint const &r) { return l /= r; }
};

template <class T>
struct NewtonInterp {
  // {(x_0,y_0),...,(x_{n-1},y_{n-1})}
  vector<pair<T, T>> p;
  // dy[r][l] = [y_l,...,y_r]
  vector<vector<T>> dy;
  // (x-x_0)...(x-x_{n-1})
  vector<T> base;
  // [y_0]+...+[y_0,y_1,...,y_n](x-x_0)...(x-x_{n-1})
  vector<T> poly;

  void insert(T const &x, T const &y) {
    p.emplace_back(x, y);
    size_t n = p.size();
    if (n == 1) {
      base.push_back(1);
    } else {
      size_t m = base.size();
      base.push_back(0);
      for (size_t i = m; i; --i) base[i] = base[i - 1];
      base[0] = 0;
      for (size_t i = 0; i < m; ++i)
        base[i] = base[i] - p[n - 2].first * base[i + 1];
    }
    dy.emplace_back(p.size());
    dy[n - 1][n - 1] = y;
    if (n > 1) {
      for (size_t i = n - 2; ~i; --i)
        dy[n - 1][i] =
            (dy[n - 2][i] - dy[n - 1][i + 1]) / (p[i].first - p[n - 1].first);
    }
    poly.push_back(0);
    for (size_t i = 0; i < n; ++i) poly[i] = poly[i] + dy[n - 1][0] * base[i];
  }

  T eval(T const &x) {
    T ans{};
    for (auto it = poly.rbegin(); it != poly.rend(); ++it) ans = ans * x + *it;
    return ans;
  }
};

int main() {
  NewtonInterp<mint> ip;
  int n, k;
  cin >> n >> k;
  for (int i = 1, x, y; i <= n; ++i) {
    cin >> x >> y;
    ip.insert(x, y);
  }
  cout << ip.eval(k).v_;
  return 0;
}

Newton Interpolation with Consecutive Integer x-coordinates

For example: Find the coefficients of polynomial \(f(x)=\sum_{i=0}^{3} a_ix^i\), given that the values of \(f(1)\) to \(f(6)\) are \(1, 5, 14, 30, 55, 91\) respectively.

\[ \begin{array}{cccccccccccc} 1 & & 5 & & 14 & & 30 & & 55 & & 91 & \\ & 4 & & 9 & & 16 & & 25 & & 36 & \\ & & 5 & & 7 & & 9 & & 11 & \\ & & & 2 & & 2 & & 2 & \\ \end{array} \]

The first row contains the first \(n\) consecutive values of \(f(x)\); each subsequent row contains the differences between adjacent terms of the previous row. It can be observed that if this operation is repeated enough times (provided \(f(x)\) is a polynomial), it will eventually yield a constant value.

The first term of the \((i-1)\)-th order difference is \(\sum_{j=1}^{i}(-1)^{i+j}\binom{i-1}{j-1}f(j)\). The contribution of the first term of the \((i-1)\)-th order difference to \(f(k)\) is \(\binom{k-1}{i-1}\) times.

\[ f(k)=\sum_{i=1}^n\binom{k-1}{i-1}\sum_{j=1}^{i}(-1)^{i+j}\binom{i-1}{j-1}f(j) \]

Time complexity is \(O(n^2)\).

Implementation in C++

Since C++20, the standard library has added std::midpoint and std::lerp functions, for computing the midpoint and linear interpolation, respectively.

Exercises

References

  1. Interpolation - Wikipedia
  2. Newton polynomial - Wikipedia
本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：OI-wiki
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用