cp-library

This documentation is automatically generated by online-judge-tools/verification-helper

View the Project on GitHub rniya/cp-library

:heavy_check_mark: Lagrange Interpolation
(src/polynomial/lagrange_interpolation.hpp)

入力

$N$ 次多項式 $f$ の評価値 $f(0), f(1), \dots , f(N)$ 及び整数 $T$．

出力

$f(T)$．

計算量

時間計算量 $\mathrm{O}(N)$

概要

$x _ i$ が相異なりかつ $f(x _ i) = y _ i$ を満たすような組 $(x _ 0, y _ 0), (x _ 1, y _ 1), \dots , (x _ N, y _ N)$ が与えられるとき，これを満たす多項式 $f$ は一意に定まる． 具体的には，

\[f _ i(x) = \prod _ {j \neq i} (x - x _ j) = \frac{\prod _ j (x - x _ j)}{x - x _ i}\]

として，

\[f(x) = \sum _ {i = 0} ^ n y _ i \frac{f _ i(x)}{f _ i(x _ i)}\]

とすることで復元可能である． 一般の場合には $\mathrm{O}(N ^ 2)$ で $f$ の各項の係数を計算可能であるが，特に $\lbrace x _ i \rbrace$ が等差数列を成す場合には，

\[f _ {i + 1}(x _ {i + 1}) = \prod _ {j \neq i} (x _ {i + 1} - x _ j) = f _ i(x _ i) \frac{x _ {i + 1} - x _ 0}{x _ i - x _ N}\]

を利用して $f _ i(x _ i)$ から $f _ {i + 1}(x _ {i + 1})$ を適切な前計算の元 $\mathrm{O}(1)$ で計算可能なので，任意の $T$ に対する $f(T)$ を全体 $\mathrm{O}(N)$ で計算することができる． $x _ i = i$ の場合には階乗前計算のもとでより簡単に計算することもできる．

出題例

Depends on

Verified with

Code

#pragma once
#include <vector>
#include "../math/binomial.hpp"

template <typename T> T lagrange_interpolation(const std::vector<T>& y, long long x, Binomial<T>& BINOM) {
    int n = y.size() - 1;
    if (0 <= x and x <= n) return y[x];
    std::vector<T> left(n + 1, 1), right(n + 1, 1);
    for (int i = 0; i < n; i++) left[i + 1] = left[i] * (x - i);
    for (int i = n; i > 0; i--) right[i - 1] = right[i] * (x - i);
    T res = 0;
    for (int i = 0; i <= n; i++) {
        T add = y[i] * left[i] * right[i] * BINOM.finv(i) * BINOM.finv(n - i);
        res += ((n - i) & 1) ? -add : add;
    }
    return res;
}
#line 2 "src/polynomial/lagrange_interpolation.hpp"
#include <vector>
#line 2 "src/math/binomial.hpp"
#include <algorithm>
#include <cassert>
#line 5 "src/math/binomial.hpp"

template <typename T> struct Binomial {
    Binomial(int MAX = 0) : n(1), facs(1, T(1)), finvs(1, T(1)), invs(1, T(1)) {
        assert(T::mod() != 0);
        if (MAX > 0) extend(MAX + 1);
    }

    T fac(int i) {
        assert(i >= 0);
        while (n <= i) extend();
        return facs[i];
    }

    T finv(int i) {
        assert(i >= 0);
        while (n <= i) extend();
        return finvs[i];
    }

    T inv(int i) {
        assert(i >= 0);
        while (n <= i) extend();
        return invs[i];
    }

    T P(int n, int r) {
        if (n < r or r < 0) return 0;
        return fac(n) * finv(n - r);
    }

    T C(int n, int r) {
        if (n < r or r < 0) return 0;
        return fac(n) * finv(n - r) * finv(r);
    }

    T H(int n, int r) {
        if (n < 0 or r < 0) return 0;
        return r == 0 ? 1 : C(n + r - 1, r);
    }

    T negative_binom(int n, int k) { return H(k, n); }

    T C_naive(int n, int r) {
        if (n < r or r < 0) return 0;
        T res = 1;
        r = std::min(r, n - r);
        for (int i = 1; i <= r; i++) res *= inv(i) * (n--);
        return res;
    }

    T catalan(int n) {
        if (n < 0) return 0;
        return fac(2 * n) * finv(n + 1) * finv(n);
    }

    T catalan_pow(int n, int k) {
        if (n < 0 or k < 0) return 0;
        if (k == 0) return n == 0 ? 1 : 0;
        return inv(n + k) * k * C(2 * n + k - 1, n);
    }

    T calatan1(int n, int m) { return C(n + m, m) - C(n + m, m - 1); }

    T catalan2(int n, int m, int k) { return n - m <= -k ? 0 : C(n + m, m) - C(n + m, m - k); }

    T narayana(int n, int k) {
        if (n < k or k <= 0) return 0;
        return C(n, k) * C(n, k - 1) * inv(n);
    }

    T grid_sum(int x, int y) {
        if (x < 0 or y < 0) return 0;
        return C(x + y + 2, x + 1) - 1;
    }

    T grid_sum2(int xl, int xr, int yl, int yr) {
        if (xl >= xr or yl >= yr) return 0;
        xl--, xr--, yl--, yr--;
        return grid_sum(xr, yr) - grid_sum(xl, yr) - grid_sum(xr, yl) + grid_sum(xl, yl);
    }

  private:
    int n;
    std::vector<T> facs, finvs, invs;

    inline void extend(int m = -1) {
        if (m == -1) m = n * 2;
        m = std::min(m, T::mod());
        if (n >= m) return;
        facs.resize(m);
        finvs.resize(m);
        invs.resize(m);
        for (int i = n; i < m; i++) facs[i] = facs[i - 1] * i;
        finvs[m - 1] = T(1) / facs[m - 1];
        invs[m - 1] = finvs[m - 1] * facs[m - 2];
        for (int i = m - 2; i >= n; i--) {
            finvs[i] = finvs[i + 1] * (i + 1);
            invs[i] = finvs[i] * facs[i - 1];
        }
        n = m;
    }
};
#line 4 "src/polynomial/lagrange_interpolation.hpp"

template <typename T> T lagrange_interpolation(const std::vector<T>& y, long long x, Binomial<T>& BINOM) {
    int n = y.size() - 1;
    if (0 <= x and x <= n) return y[x];
    std::vector<T> left(n + 1, 1), right(n + 1, 1);
    for (int i = 0; i < n; i++) left[i + 1] = left[i] * (x - i);
    for (int i = n; i > 0; i--) right[i - 1] = right[i] * (x - i);
    T res = 0;
    for (int i = 0; i <= n; i++) {
        T add = y[i] * left[i] * right[i] * BINOM.finv(i) * BINOM.finv(n - i);
        res += ((n - i) & 1) ? -add : add;
    }
    return res;
}
Back to top page