library

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

View the Project on GitHub kk2a/library

:heavy_check_mark: fps/poly_interpolation.hpp

Depends on

Verified with

Code

#ifndef KK2_FPS_POLY_INTERPOLATION_HPP
#define KK2_FPS_POLY_INTERPOLATION_HPP 1

#include <cassert>
#include <functional>
#include <vector>

#include "chirp_Z.hpp"
#include "poly_multi_eval.hpp"

namespace kk2 {

// return FPS f s.t. f(x[i]) = y[i] for all i
template <class FPS, class mint = typename FPS::value_type>
FPS poly_interpolation(const std::vector<mint> &x, const std::vector<mint> &y) {
    assert(x.size() == y.size());
    SubProductTree<FPS> mpe(x);
    FPS gp = mpe.pr[1].diff();
    std::vector<mint> vs = mpe.query(gp);
    auto rec = [&](auto self, int idx) -> FPS {
        if (idx >= mpe.size) {
            if (idx - mpe.size < (int)y.size()) {
                return {y[idx - mpe.size] / vs[idx - mpe.size]};
            } else return {mint(1)};
        }
        if (mpe.pr[idx << 1 | 0].empty()) return {};
        else if (mpe.pr[idx << 1 | 1].empty()) return self(self, idx << 1 | 0);
        return self(self, idx << 1 | 0) * mpe.pr[idx << 1 | 1]
               + self(self, idx << 1 | 1) * mpe.pr[idx << 1 | 0];
    };
    return rec(rec, 1);
}

// return FPS f s.t. f(ar^i) = y[i] for all i
template <class FPS, class mint = typename FPS::value_type>
FPS poly_interpolation_geo(const mint &a, const mint &r, const std::vector<mint> &y) {
    // reference:
    // https://noshi91.github.io/algorithm-encyclopedia/polynomial-interpolation-geometric#fn:Bostan

    if (y.empty()) return {};
    if (y.size() == 1) return FPS{y[0]};
    assert(a != mint(0) && r != mint(0) && r != mint(1));

    int n = (int)y.size();
    // https://yukicoder.me/wiki/%E9%80%86%E5%85%83
    std::vector<mint> s(n + 1), invs(n), pre(n + 1);
    s[0] = pre[0] = invs[0] = 1;
    mint q = r;
    for (int i = 1; i < n + 1; i++) {
        s[i] = s[i - 1] * (-q + 1);
        pre[i] = pre[i - 1] * s[i];
        q *= r;
    }
    mint inv = pre[n - 1].inv();
    for (int i = n - 1; i >= 1; i--) {
        invs[i] = inv * pre[i - 1];
        inv *= s[i];
    }

    std::vector<mint> w(n);
    q = 1;
    int idx1 = n - 1, idx2 = 0;
    w[n - 1] = r.pow(1ll * (n - 1) * (n - 2) / 2).inv() * invs[idx1] * invs[idx2];
    if ((n - 1) & 1) w[n - 1] *= -1;
    for (int i = n - 1; i > 0; i--) {
        w[i - 1] = w[i] * q * (-1) * s[idx1] * invs[idx1 - 1] * s[idx2] * invs[idx2 + 1];
        q *= r;
        idx1--;
        idx2++;
    }
    for (int i = 0; i < n; i++) w[i] *= y[i];

    FPS g{begin(w), end(w)};
    std::vector<mint> tmp = ChirpZ(g, r, n);
    FPS gq{begin(tmp), end(tmp)};
    FPS prod(n);
    q = 1;
    mint plus = r;
    prod[0] = 1;
    for (int i = 1; i < n; i++) {
        prod[i] = s[n] * invs[i] * invs[n - i] * q;
        if (i & 1) prod[i] *= -1;
        q *= plus;
        plus *= r;
    }

    FPS ret = (prod * gq).pre(n).rev();
    if (a != mint(1)) {
        mint x = 1, inva = a.inv();
        for (int i = 0; i < n; i++) {
            ret[i] *= x;
            x *= inva;
        }
    }
    return ret;
}

} // namespace kk2

#endif // KK2_FPS_POLY_INTERPOLATION_HPP
#line 1 "fps/poly_interpolation.hpp"



#include <cassert>
#include <functional>
#include <vector>

#line 1 "fps/chirp_Z.hpp"



#include <algorithm>
#line 6 "fps/chirp_Z.hpp"

namespace kk2 {

// return f(aw^0), f(aw^1), ..., f(aw^(n - 1))
template <class FPS, class mint = typename FPS::value_type>
std::vector<mint> ChirpZ(const FPS &f_, mint w, int n = -1, mint a = 1) {
    FPS f(f_.begin(), f_.end());
    if (n == -1) n = f.size();
    if (f.empty() || n == 0) return std::vector<mint>(n, mint(0));
    int m = f.size();
    if (a != mint(1)) {
        mint x = 1;
        for (int i = 0; i < m; i++) {
            f[i] *= x;
            x *= a;
        }
    }
    if (w == mint(0)) {
        std::vector<mint> g(n, f[0]);
        for (int i = 1; i < m; i++) g[0] += f[i];
        return g;
    }
    FPS wc(n + m), iwc(std::max(n, m));
    mint ws = 1, iw = w.inv(), iws = 1;
    wc[0] = iwc[0] = 1;
    for (int i = 1; i < n + m; i++) {
        wc[i] = ws * wc[i - 1];
        ws *= w;
    }
    for (int i = 1; i < std::max(n, m); i++) {
        iwc[i] = iws * iwc[i - 1];
        iws *= iw;
    }
    for (int i = 0; i < m; i++) f[i] *= iwc[i];
    std::reverse(std::begin(f), std::end(f));
    FPS g = f * wc;
    std::vector<mint> ret{std::begin(g) + m - 1, std::begin(g) + m + n - 1};
    for (int i = 0; i < n; i++) ret[i] *= iwc[i];
    return ret;
}

} // namespace kk2


#line 1 "fps/poly_multi_eval.hpp"



#line 6 "fps/poly_multi_eval.hpp"

namespace kk2 {

template <class FPS, class mint = typename FPS::value_type> struct SubProductTree {
    int _n, size;
    std::vector<int> l, r;
    std::vector<FPS> pr;
    std::vector<mint> v;
    FPS f;

    SubProductTree(const std::vector<mint> &v_) : _n(int(v_.size())), v(v_) {
        size = 1;
        while (size < _n) size <<= 1;
        pr.resize(size << 1);
        l.resize(size << 1, _n);
        r.resize(size << 1, _n);
        build();
    }

    SubProductTree(const std::vector<mint> &v_, const FPS &f_) : SubProductTree(v_) {
        this->f = f_;
    }

    void set(const FPS &f_) { this->f = f_; }

    void build() {
        for (int i = 0; i < _n; i++) {
            l[i + size] = i;
            r[i + size] = i + 1;
            pr[i + size] = {-v[i], 1};
        }
        for (int i = size - 1; i > 0; i--) {
            l[i] = l[i << 1 | 0];
            r[i] = r[i << 1 | 1];
            if (pr[i << 1 | 0].empty()) continue;
            else if (pr[i << 1 | 1].empty()) pr[i] = pr[i << 1 | 0];
            else pr[i] = pr[i << 1 | 0] * pr[i << 1 | 1];
        }
    }

    std::vector<mint> query(const FPS &f) {
        this->f = f;
        return query();
    }

    std::vector<mint> query() {
        if (f.empty() || !_n) return FPS(_n, mint(0));
        std::vector<mint> ret;
        ret.reserve(_n);
        auto rec = [&](auto self, FPS a, int idx) -> void {
            if (l[idx] == r[idx]) return;
            a %= pr[idx];
            if (a.size() <= 64u) {
                for (int i = l[idx]; i < r[idx]; i++) { ret.push_back(a.eval(v[i])); }
                return;
            }
            self(self, a, idx << 1 | 0);
            self(self, a, idx << 1 | 1);
        };
        rec(rec, f, 1);
        return ret;
    }
};

template <class FPS, class mint = typename FPS::value_type>
std::vector<mint> MultiEval(const std::vector<mint> &v, const FPS &f) {
    SubProductTree<FPS> mpe(v, f);
    return mpe.query();
}

} // namespace kk2


#line 10 "fps/poly_interpolation.hpp"

namespace kk2 {

// return FPS f s.t. f(x[i]) = y[i] for all i
template <class FPS, class mint = typename FPS::value_type>
FPS poly_interpolation(const std::vector<mint> &x, const std::vector<mint> &y) {
    assert(x.size() == y.size());
    SubProductTree<FPS> mpe(x);
    FPS gp = mpe.pr[1].diff();
    std::vector<mint> vs = mpe.query(gp);
    auto rec = [&](auto self, int idx) -> FPS {
        if (idx >= mpe.size) {
            if (idx - mpe.size < (int)y.size()) {
                return {y[idx - mpe.size] / vs[idx - mpe.size]};
            } else return {mint(1)};
        }
        if (mpe.pr[idx << 1 | 0].empty()) return {};
        else if (mpe.pr[idx << 1 | 1].empty()) return self(self, idx << 1 | 0);
        return self(self, idx << 1 | 0) * mpe.pr[idx << 1 | 1]
               + self(self, idx << 1 | 1) * mpe.pr[idx << 1 | 0];
    };
    return rec(rec, 1);
}

// return FPS f s.t. f(ar^i) = y[i] for all i
template <class FPS, class mint = typename FPS::value_type>
FPS poly_interpolation_geo(const mint &a, const mint &r, const std::vector<mint> &y) {
    // reference:
    // https://noshi91.github.io/algorithm-encyclopedia/polynomial-interpolation-geometric#fn:Bostan

    if (y.empty()) return {};
    if (y.size() == 1) return FPS{y[0]};
    assert(a != mint(0) && r != mint(0) && r != mint(1));

    int n = (int)y.size();
    // https://yukicoder.me/wiki/%E9%80%86%E5%85%83
    std::vector<mint> s(n + 1), invs(n), pre(n + 1);
    s[0] = pre[0] = invs[0] = 1;
    mint q = r;
    for (int i = 1; i < n + 1; i++) {
        s[i] = s[i - 1] * (-q + 1);
        pre[i] = pre[i - 1] * s[i];
        q *= r;
    }
    mint inv = pre[n - 1].inv();
    for (int i = n - 1; i >= 1; i--) {
        invs[i] = inv * pre[i - 1];
        inv *= s[i];
    }

    std::vector<mint> w(n);
    q = 1;
    int idx1 = n - 1, idx2 = 0;
    w[n - 1] = r.pow(1ll * (n - 1) * (n - 2) / 2).inv() * invs[idx1] * invs[idx2];
    if ((n - 1) & 1) w[n - 1] *= -1;
    for (int i = n - 1; i > 0; i--) {
        w[i - 1] = w[i] * q * (-1) * s[idx1] * invs[idx1 - 1] * s[idx2] * invs[idx2 + 1];
        q *= r;
        idx1--;
        idx2++;
    }
    for (int i = 0; i < n; i++) w[i] *= y[i];

    FPS g{begin(w), end(w)};
    std::vector<mint> tmp = ChirpZ(g, r, n);
    FPS gq{begin(tmp), end(tmp)};
    FPS prod(n);
    q = 1;
    mint plus = r;
    prod[0] = 1;
    for (int i = 1; i < n; i++) {
        prod[i] = s[n] * invs[i] * invs[n - i] * q;
        if (i & 1) prod[i] *= -1;
        q *= plus;
        plus *= r;
    }

    FPS ret = (prod * gq).pre(n).rev();
    if (a != mint(1)) {
        mint x = 1, inva = a.inv();
        for (int i = 0; i < n; i++) {
            ret[i] *= x;
            x *= inva;
        }
    }
    return ret;
}

} // namespace kk2
Back to top page