AIOJ

Introduction

Interpolation is a method for estimating new data points within a certain range from known, discrete data points. It is commonly used for function fitting.

For example, consider 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$

$Figure: interp 1$

Here $f(x)$ is unknown, and interpolation can estimate the unknown data points by fitting $f(x)$ in a prescribed form.

For example, we can fit $f(x)$ with a piecewise linear function:

$Figure: interp 2$

This interpolation method is called linear interpolation.

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

$Figure: interp 3$

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$ satisfying

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

of a polynomial $f(x)$ .

Below we introduce two methods of polynomial interpolation: Lagrange interpolation and Newton interpolation. It is not difficult to prove that these two methods produce the same result.

Lagrange interpolation

We need to construct a function $f(x)$ passing through the 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}$ . It follows that the function required by the problem is $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 the 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}

Therefore, the Lagrange interpolation formula 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 Fast multipoint polynomial 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$ , and $f(x_i)\equiv y_i\pmod{998244353}$ and $\deg(f(x)) < n$ (define $\deg(0)=-\infty$ ), find $f(k)\bmod{998244353}$ .

Solution

In this problem, we only need the value of $f(k)$ , so we can directly substitute $k$ into the formula above during the computation. Sometimes we need to evaluate the polynomial at multiple points or perform other more complicated operations; in that case, we need to find all coefficients of $f$ . The code provides an implementation that computes the coefficients.

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

This problem also requires modular inverses. If we first compute the numerator and denominator separately, then multiply the numerator by the inverse of the denominator and add it to the final answer, the bottleneck in the time complexity will not be the modular inverse computation; the time complexity is $O(n^2)$ .

Because all computations are performed modulo $998244353$ , we temporarily regard the time complexity of computing a modular inverse as constant.

Code implementation

#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;
}

// Return f satisfying f(x_i) = y_i
// Clearly O(n^2), excluding the time for modular inverses
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 into 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 abscissas

If the abscissas of the known points are consecutive integers, interpolation can be done in $O(n)$ time.

Let the polynomial to be found be $f(x)$ , and suppose that we know $f(1),\cdots,f(n+1)$ ( $1\le i\le n+1$ ). Substituting into the interpolation formula above gives:

\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 products can be considered separately in the numerator and denominator. The numerator is:

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

The product of $i-j$ in the denominator can be split into two factorials:

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

Thus, the interpolation formula for abscissas $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)}

Precompute the prefix and suffix products of $(x-i)$ , factorials, and inverse factorials, then substitute into this formula; the complexity is $O(n)$ .

Example CF622F The Sum of the k-th Powers

Given $n,k$ , find the value of $\sum\limits_{i=1}^ni^k$ modulo $10^9+7$ .

Solution

The answer is a polynomial of degree $k+1$ , so we can use a linear sieve to compute the values of $1^i,\cdots,(k+2)^i$ , and then perform $O(n)$ interpolation.

Alternatively, the following formula can be derived from relevant combinatorial knowledge using finite differences:

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

// 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 of interpolation based on higher-order differences. Its advantage is that it supports inserting a new data point in $O(n)$ time.

To achieve $O(n)$ insertion of a new data point, 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.

Once $a_j$ is determined, the interpolation polynomial $f(x)$ is obtained. 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 formula. A naive implementation has time complexity $O(n^2)$ .

If the sample points are equally spaced (that is, $x_i=x_0+ih$ for $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$ is the forward difference, defined as follows:

\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)

#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 abscissas

For example, find the coefficients of the polynomial $f(x)=\sum_{i=0}^{3} a_ix^i$ , given that the values of $f(1)$ through $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 consists of the first consecutive terms of the function $f(x)$ , $n$ of them; each subsequent row contains the differences between adjacent terms in the preceding row. Observe that if this operation is repeated sufficiently many times (provided that $f(x)$ is a polynomial), it eventually always produces a constant.

The leading term of the $i-1$ -st difference is $\sum_{j=1}^{i}(-1)^{i+j}\binom{i-1}{j-1}f(j)$ , and the contribution of the leading term of the $i-1$ -st difference to $f(k)$ is multiplied by $\binom{k-1}{i-1}$ .

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)

The time complexity is $O(n^2)$ .

Implementation in C++

Since C++ 20, the standard library has provided std::midpoint and std::lerp, which are used to compute midpoints and perform linear interpolation, respectively.

Interpolation