fps/Polynomial_Interpolation.hpp

• View this file on GitHub
• Last update: 2025-10-18 21:09:10+09:00
• Include: #include "fps/Polynomial_Interpolation.hpp"

Depends on

• fps/FPS_Product_Sequence.hpp
• fps/FPS_add.hpp
• fps/FPS_division.hpp
• fps/FPS_inv.hpp
• fps/Multipoint_Evaluation.hpp

Verified with

• test/fps/polynomial_interpolation.test.cpp

Code

#pragma once
#include "FPS_Product_Sequence.hpp"
#include "Multipoint_Evaluation.hpp"
#include "FPS_add.hpp"

namespace po167{
template<class T>
// ラグランジュの多項式補完
// f(X[i]) = Y[i] である f を返す
std::vector<T> Polynomial_Interpolation(std::vector<int> X, std::vector<T> Y){
    int N = X.size();
    assert(Y.size() == X.size());
    if (N == 0) return {};
    {
        auto Z = X;
        std::sort(Z.begin(), Z.end());
        for (int i = 0; i < N - 1; i++){
            assert(Z[i] != Z[i + 1]);
        }
    }
    std::vector<std::vector<T>> p(N);
    for (int i = 0; i < N; i++) p[i] = {-X[i], 1};
    auto g = FPS_Product_Sequence(p);
    for (int i = 0; i < N; i++){
        g[i] = g[i + 1] * (i + 1);
    }
    g.pop_back();
    std::vector<T> xt(N);
    for (int i = 0; i < N; i++) xt[i] = X[i];
    auto Z = Multipoint_Evaluation(g, xt);
    auto rec = [&](auto self, int l, int r) -> std::pair<std::vector<T>, std::vector<T>> {
        if (l + 1 == r){
            return {{Y[l] / Z[l]}, {-X[l], 1}};
        }
        int m = (l + r) / 2;
        auto L = self(self, l, m);
        auto R = self(self, m, r);
        auto D = atcoder::convolution(L.second, R.second);
        auto U = atcoder::convolution(L.second, R.first);
        FPS_add(U, atcoder::convolution(L.first, R.second));
        return {U, D};
    };
    return rec(rec, 0, N).first;
}
}

#line 2 "fps/FPS_Product_Sequence.hpp"
#include <vector>
#include <atcoder/convolution>

namespace po167{
template<class T>
std::vector<T> FPS_Product_Sequence(std::vector<std::vector<T>> f){
    if (f.empty()) return {1};
    auto op = [&](auto self,int l, int r) -> std::vector<T> {
        if (l + 1 == r) return f[l];
        int m = (l + r) / 2;
        return atcoder::convolution(self(self, l, m), self(self, m, r));
    };
    return op(op, 0, f.size());
}
}
#line 4 "fps/FPS_division.hpp"

#line 4 "fps/FPS_inv.hpp"

namespace po167{
// return 1 / f
template <class T>
std::vector<T> FPS_inv(std::vector<T> f, int len = -1){
    if (len == -1) len = f.size();
    assert(f[0] != 0);
    std::vector<T> g = {1 / f[0]};
    int s = 1;
    while(s < len){
        // g = 2g_s - f(g_s)^2 (mod x ^ (2 * s))
        // g = g - (fg - 1)g
        // (fg - 1) = 0 (mod x ^ (s))
        std::vector<T> n_g(s * 2, 0);
        std::vector<T> f_s(s * 2, 0);
        g.resize(s * 2);
        for (int i = 0; i < s * 2; i++){
            if (int(f.size()) > i) f_s[i] = f[i];
            n_g[i] = g[i];
        }
        atcoder::internal::butterfly(g);
        atcoder::internal::butterfly(f_s);
        for (int i = 0; i < s * 2; i++){
            f_s[i] *= g[i];
        }
        atcoder::internal::butterfly_inv(f_s);
        T iz = 1 / (T)(s * 2);
        for (int i = s; i < s * 2; i++){
            f_s[i] *= iz;
        }
        for (int i = 0; i < s; i++){
            f_s[i] = 0;
        }
        atcoder::internal::butterfly(f_s);
        for (int i = 0; i < s * 2; i++){
            f_s[i] *= g[i];
        }
        atcoder::internal::butterfly_inv(f_s);
        for (int i = s; i < s * 2; i++){
            n_g[i] -= f_s[i] * iz;
        }
        std::swap(n_g, g);
        s *= 2;
    }
    g.resize(len);
    return g;
}
}
#line 6 "fps/FPS_division.hpp"
namespace po167{
template<class T>
// f = g * res.first + res.second
// |res.first| <= |f| - |g| + 1
// |res.second| <= |g| - 1
std::pair<std::vector<T>, std::vector<T>>
FPS_division(std::vector<T> f, std::vector<T> g){
    while (!f.empty() && f.back() == 0) f.pop_back();
    assert(!g.empty() && g.back() != 0);
    if (f.size() < g.size()){
        return {{}, f};
    }
    // rev(f) / rev(g) = rev(q) (mod x ^ {|f| - |g| + 1})
    std::vector<T> r = f;
    std::reverse(f.begin(), f.end());
    std::reverse(g.begin(), g.end());
    int z = (int)f.size() - (int)g.size() + 1;
    f.resize(z);
    std::vector<T> q = atcoder::convolution(f, FPS_inv(g, z));
    q.resize(z);
    std::reverse(g.begin(), g.end());
    std::reverse(q.begin(), q.end());
    f = atcoder::convolution(q, g);
    for (int i = 0; i < (int)f.size(); i++) r[i] -= f[i];
    while (!q.empty() && q.back() == 0) q.pop_back();
    while (!r.empty() && r.back() == 0) r.pop_back();
    return {q, r};
}
}
#line 4 "fps/Multipoint_Evaluation.hpp"

namespace po167{
// return {f(p[0]), f(p[1]), f(p[2]), ... }
template <class T>
std::vector<T> Multipoint_Evaluation(
    std::vector<T> f,
    std::vector<T> p
){
    int m = p.size();
    if (m == 0) return {};
    if (m == 1){
        T res = 0;
        T tmp = 1;
        for (auto x : f) res += tmp * x, tmp *= p[0];
        return {res};
    }
    int size = 1;
    while (size < m) size *= 2;
    std::vector<std::vector<T>> prod(size * 2);
    for (int i = 0; i < size; i++){
        if (i < m) prod[i + size] = {(T)(-1) * p[i], 1};
        else prod[i + size] = {1};
    }
    for (int i = size - 1; i > 0; i--){
        prod[i] = atcoder::convolution(prod[i * 2], prod[i * 2 + 1]);
    }
    std::vector<T> res(m);
    auto calc = [&](auto self, int l, int r, int ind, std::vector<T> tmp) -> void {
        if (m <= l) return;
        if (l + 1 == r){
            res[l] = (tmp.empty() ? T(0) : tmp[0]);
            return;
        }
        int mid = (l + r) / 2;
        self(self, l, mid, ind * 2, po167::FPS_division(tmp, prod[ind * 2]).second);
        self(self, mid, r, ind * 2 + 1, po167::FPS_division(tmp, prod[ind * 2 + 1]).second);
    };calc(calc, 0, size, 1, f);
    return res;
}
}
#line 3 "fps/FPS_add.hpp"

namespace po167{
template <class T>
// a(x) += b(x) * c * x^d
void FPS_add(std::vector<T> &a, std::vector<T> b, T c = 1, int d = 0){
    for (int i = 0; i < (int)(b.size()); i++){
        while ((int)a.size() <= i + d) a.push_back((T)0);
        a[i + d] += b[i] * c;
    }
}
}
#line 5 "fps/Polynomial_Interpolation.hpp"

namespace po167{
template<class T>
// ラグランジュの多項式補完
// f(X[i]) = Y[i] である f を返す
std::vector<T> Polynomial_Interpolation(std::vector<int> X, std::vector<T> Y){
    int N = X.size();
    assert(Y.size() == X.size());
    if (N == 0) return {};
    {
        auto Z = X;
        std::sort(Z.begin(), Z.end());
        for (int i = 0; i < N - 1; i++){
            assert(Z[i] != Z[i + 1]);
        }
    }
    std::vector<std::vector<T>> p(N);
    for (int i = 0; i < N; i++) p[i] = {-X[i], 1};
    auto g = FPS_Product_Sequence(p);
    for (int i = 0; i < N; i++){
        g[i] = g[i + 1] * (i + 1);
    }
    g.pop_back();
    std::vector<T> xt(N);
    for (int i = 0; i < N; i++) xt[i] = X[i];
    auto Z = Multipoint_Evaluation(g, xt);
    auto rec = [&](auto self, int l, int r) -> std::pair<std::vector<T>, std::vector<T>> {
        if (l + 1 == r){
            return {{Y[l] / Z[l]}, {-X[l], 1}};
        }
        int m = (l + r) / 2;
        auto L = self(self, l, m);
        auto R = self(self, m, r);
        auto D = atcoder::convolution(L.second, R.second);
        auto U = atcoder::convolution(L.second, R.first);
        FPS_add(U, atcoder::convolution(L.first, R.second));
        return {U, D};
    };
    return rec(rec, 0, N).first;
}
}

Back to top page