Bostan-Mori Algorithm
(math/bostan_mori_algorithm.hpp)

• View this file on GitHub
• Last update: 2023-06-04 13:46:26+09:00
• Include: #include "math/bostan_mori_algorithm.hpp"

説明

$P(x)$ を高々 $d-1$ 次の多項式、 $Q(x)$ を $d$ 次の多項式とし、 $[x^0]Q(x) = 1$ とする。

$[x^N] \frac{P(x)}{Q(x)}$

を求める。$d$ 次の多項式積を $M(d)$ でできるとすると、計算量は $O(M(d) \log N)$ となる。

bostan_mori_algorithm(std::int64_t n, std::vector p, std::vector q)

$[x^N] \frac{P(x)}{Q(x)}$ を求める。

kitamasa(std::int64_t n, std::vector a, std::vector c)

$a$ を 数列の先頭 $d$ 項として、線形漸化式 $a_i = \sum_{j = 1}^d c_j a_{i-j}$ を満たす数列の第 $n$ 項を求める。 $O(M(d) \log N)$

Verified with

• test/math/Kth_term_of_Linearly_Recurrent_Sequence.test.cpp

Code

#pragma once

#include <cstdint>
#include <vector>

namespace ebi {

template <class T, std::vector<T> (*convolution)(const std::vector<T> &,
                                                 const std::vector<T> &)>
T bostan_mori_algorithm(std::int64_t n, std::vector<T> p, std::vector<T> q) {
    while (n > 0) {
        auto q_neg = q;
        for (int i = 1; i < (int)q_neg.size(); i += 2) q_neg[i] = -q_neg[i];
        p = convolution(p, q_neg);
        q = convolution(q, q_neg);
        for (int i = (n & 1LL); i < (int)p.size(); i += 2) p[i >> 1] = p[i];
        p.resize(((int)p.size() + 1 - (n & 1LL)) / 2);
        for (int i = 0; i < (int)q.size(); i += 2) q[i >> 1] = q[i];
        q.resize(((int)q.size() + 1) / 2);
        n >>= 1;
    }
    return p[0] / q[0];
}

template <class T, std::vector<T> (*convolution)(const std::vector<T> &,
                                                 const std::vector<T> &)>
T kitamasa(std::int64_t n, std::vector<T> a, std::vector<T> c) {
    if (n < (int)a.size()) return a[n];
    const int d = c.size();
    for (auto &val : c) val = -val;
    c.insert(c.begin(), 1);
    auto p = convolution(a, c);
    p.resize(d);
    return bostan_mori_algorithm<T, convolution>(n, p, c);
}

}  // namespace ebi

#line 2 "math/bostan_mori_algorithm.hpp"

#include <cstdint>
#include <vector>

namespace ebi {

template <class T, std::vector<T> (*convolution)(const std::vector<T> &,
                                                 const std::vector<T> &)>
T bostan_mori_algorithm(std::int64_t n, std::vector<T> p, std::vector<T> q) {
    while (n > 0) {
        auto q_neg = q;
        for (int i = 1; i < (int)q_neg.size(); i += 2) q_neg[i] = -q_neg[i];
        p = convolution(p, q_neg);
        q = convolution(q, q_neg);
        for (int i = (n & 1LL); i < (int)p.size(); i += 2) p[i >> 1] = p[i];
        p.resize(((int)p.size() + 1 - (n & 1LL)) / 2);
        for (int i = 0; i < (int)q.size(); i += 2) q[i >> 1] = q[i];
        q.resize(((int)q.size() + 1) / 2);
        n >>= 1;
    }
    return p[0] / q[0];
}

template <class T, std::vector<T> (*convolution)(const std::vector<T> &,
                                                 const std::vector<T> &)>
T kitamasa(std::int64_t n, std::vector<T> a, std::vector<T> c) {
    if (n < (int)a.size()) return a[n];
    const int d = c.size();
    for (auto &val : c) val = -val;
    c.insert(c.begin(), 1);
    auto p = convolution(a, c);
    p.resize(d);
    return bostan_mori_algorithm<T, convolution>(n, p, c);
}

}  // namespace ebi

Back to top page