fps/Multipoint_Evaluation.hpp

View this file on GitHub
Last update: 2024-11-23 03:09:08+09:00
Include: #include "fps/Multipoint_Evaluation.hpp"

Depends on

Verified with

test/fps/multpoint_evalution.test.cpp

Code

#pragma once
#include <atcoder/convolution>
#include "FPS_division.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 2 "fps/Multipoint_Evaluation.hpp"
#include <atcoder/convolution>
#line 2 "fps/FPS_division.hpp"
#include <vector>
#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;
}
}