Lagrange interpolation
(math/lagrange_interpolation.hpp)

• View this file on GitHub
• Last update: 2025-01-28 23:59:05+09:00
• Include: #include "math/lagrange_interpolation.hpp"

Depends on

• math/binomial.hpp

Code

#pragma once

#include"binomial.hpp"

namespace noya2 {

/**
 * @brief Lagrange interpolation
 * @note y is at most n-degree polynomial of x
 * 
 * @tparam mint (use noya2::binomial<mint>::ifact(int))
 * @param y value of y(0), y(1), ... y(n)
 * @param x specific value of x
 * @return mint y(x)
 */
template<typename mint>
mint lagrange_interpolation(const std::vector<mint> &y, mint x){
    if (x.val() < y.size()){
        return y[x.val()];
    }
    int n = y.size() - 1;
    std::vector<mint> lui(n+1,1), rui(n+1,1);
    mint a = x;
    for (int i = 0; i < n; i++){
        lui[i+1] = lui[i] * a;
        a -= 1;
    }
    for (int i = n-1; i >= 0; i--){
        rui[i] = rui[i+1] * a;
        a += 1;
    }
    mint ans = 0;
    binomial<mint> bnm;
    for (int i = 0; i <= n; i++){
        mint tmp = y[i] * lui[i] * rui[i] * bnm.ifact(i) * bnm.ifact(n-i);
        ans += ((n-i) & 1) ? -tmp : tmp;
    }
    return ans;
}

} // namespace noya2

#line 2 "math/lagrange_interpolation.hpp"

#line 2 "math/binomial.hpp"

#include<vector>
namespace noya2 {

template<typename mint>
struct binomial {
    binomial(int len = 300000){ extend(len); }
    static mint fact(int n){
        if (n < 0) return 0;
        while (n >= (int)_fact.size()) extend();
        return _fact[n];
    }
    static mint ifact(int n){
        if (n < 0) return 0;
        while (n >= (int)_fact.size()) extend();
        return _ifact[n];
    }
    static mint inv(int n){
        return ifact(n) * fact(n-1);
    }
    static mint C(int n, int r){
        if (!(0 <= r && r <= n)) return 0;
        return fact(n) * ifact(r) * ifact(n-r);
    }
    static mint P(int n, int r){
        if (!(0 <= r && r <= n)) return 0;
        return fact(n) * ifact(n-r);
    }
    static mint catalan(int n){
        return C(n * 2, n) * inv(n + 1);
    }
    inline mint operator()(int n, int r) { return C(n, r); }
    template<class... Cnts>
    static mint M(const Cnts&... cnts){
        return multinomial(0,1,cnts...);
    }
    static void initialize(int len = 2){
        _fact.clear();
        _ifact.clear();
        extend(len);
    }
  private:
    static mint multinomial(const int& sum, const mint& div_prod){
        if (sum < 0) return 0;
        return fact(sum) * div_prod;
    }
    template<class... Tail>
    static mint multinomial(const int& sum, const mint& div_prod, const int& n1, const Tail&... tail){
        if (n1 < 0) return 0;
        return multinomial(sum+n1,div_prod*ifact(n1),tail...);
    }
    static inline std::vector<mint> _fact, _ifact;
    static void extend(int len = -1){
        if (_fact.empty()){
            _fact = _ifact = {1,1};
        }
        int siz = _fact.size();
        if (len == -1) len = siz * 2;
        len = (int)min<long long>(len, mint::mod() - 1);
        if (len < siz) return ;
        _fact.resize(len+1), _ifact.resize(len+1);
        for (int i = siz; i <= len; i++) _fact[i] = _fact[i-1] * i;
        _ifact[len] = _fact[len].inv();
        for (int i = len; i > siz; i--) _ifact[i-1] = _ifact[i] * i;
    }
};

} // namespace noya2
#line 4 "math/lagrange_interpolation.hpp"

namespace noya2 {

/**
 * @brief Lagrange interpolation
 * @note y is at most n-degree polynomial of x
 * 
 * @tparam mint (use noya2::binomial<mint>::ifact(int))
 * @param y value of y(0), y(1), ... y(n)
 * @param x specific value of x
 * @return mint y(x)
 */
template<typename mint>
mint lagrange_interpolation(const std::vector<mint> &y, mint x){
    if (x.val() < y.size()){
        return y[x.val()];
    }
    int n = y.size() - 1;
    std::vector<mint> lui(n+1,1), rui(n+1,1);
    mint a = x;
    for (int i = 0; i < n; i++){
        lui[i+1] = lui[i] * a;
        a -= 1;
    }
    for (int i = n-1; i >= 0; i--){
        rui[i] = rui[i+1] * a;
        a += 1;
    }
    mint ans = 0;
    binomial<mint> bnm;
    for (int i = 0; i <= n; i++){
        mint tmp = y[i] * lui[i] * rui[i] * bnm.ifact(i) * bnm.ifact(n-i);
        ans += ((n-i) & 1) ? -tmp : tmp;
    }
    return ans;
}

} // namespace noya2

Back to top page