$[x^{N}]\frac{P(x)}{Q(x)}$ の高速計算 (Bostan-Mori アルゴリズム)
(Src/FPS/BostanMori.hpp)

• View this file on GitHub
• Last update: 2026-01-03 22:37:32+09:00
• Include: #include "Src/FPS/BostanMori.hpp"

概要

ライブラリの説明

template <concepts::IndexedFPS FPS, class Conv = FPSMult>
requires concepts::Convolution<FPS, Conv>
typename FPS::value_type BostanMori(usize N, FPS P, FPS Q, Conv conv = {})

$[x^N]\frac{P(x)}{Q(x)}$を計算する

convはFPSNTTFriendlyのときは何も指定しなくても良い。それ以外のときは畳み込みをする関数オブジェクトを与える。

$\Theta O(NM)$ の畳み込みはFPS.hppのNaiveConvolutionを与えると良い。

計算量

• 畳み込みを $2\times \lceil\log_{2}(N + 1)\rceil$ 回行うのがボトルネック

Depends on

• Src/FPS/FPS.hpp
• 標準データ型のエイリアス (Src/Template/TypeAlias.hpp)

Required by

• 線形漸化式のK項目を計算する (Src/FPS/KthTerm.hpp)

Verified with

• ABC436-G Linear Inequation (Test/AtCoder/abc436_g.test.cpp)
• Test/AtCoder/tdpc_fibonacci.test.cpp
• Test/LC/kth_term_of_linearly_recurrent_sequence.test.cpp
• Test/yukicoder/3044.test.cpp

Code

#pragma once

#include "FPS.hpp"

namespace zawa {

template <concepts::IndexedFPS FPS, class Conv = FPSMult>
requires concepts::Convolution<FPS, Conv>
typename FPS::value_type BostanMori(usize N, FPS P, FPS Q, Conv conv = {}) {
    assert(P.size());
    assert(Q.size() and Q[0] != 0); 
    auto takeParity = [&](const FPS& f, usize p) {
        FPS res;
        res.reserve(f.size() / 2);
        for (usize i = p ; i < f.size() ; i += 2)
            res.push_back(f[i]);
        return res;
    };
    while (N) {
        FPS Qm(Q.size());
        for (usize i = 0 ; i < Q.size() ; i++)
            Qm[i] = i % 2 ? -Q[i] : Q[i];
        P = takeParity(conv(P, Qm), N % 2);
        Q = takeParity(conv(Q, Qm), 0);
        N >>= 1;
    }
    return P[0] / Q[0];
}

} // namespace zawa

#line 2 "Src/FPS/BostanMori.hpp"

#line 2 "Src/FPS/FPS.hpp"

#line 2 "Src/Template/TypeAlias.hpp"

#include <cstdint>
#include <cstddef>

namespace zawa {

using i16 = std::int16_t;
using i32 = std::int32_t;
using i64 = std::int64_t;
using i128 = __int128_t;

using u8 = std::uint8_t;
using u16 = std::uint16_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;

using usize = std::size_t;

} // namespace zawa
#line 4 "Src/FPS/FPS.hpp"

#include <concepts>

namespace zawa {

namespace concepts {

template <class FPS>
concept IndexedFPS = requires(FPS f, usize i) {
    typename FPS::value_type;
    { f.size() } -> std::convertible_to<usize>;
    { f[i] } -> std::convertible_to<typename FPS::value_type>;
    f.reserve(0);
    f.push_back(f[i]);
};

template <class FPS, class Conv>
concept Convolution = 
    std::regular_invocable<Conv, const FPS&, const FPS&> &&
    std::same_as<std::invoke_result_t<Conv, const FPS&, const FPS&>, FPS>;

} // namespace concepts

struct FPSMult {
    template <class FPS>
    requires requires(const FPS& a, const FPS& b) {
        { a * b } -> std::same_as<FPS>;
    }
    FPS operator()(const FPS& a, const FPS& b) const {
        return a * b;
    }
};

struct NaiveConvolution {
    template <class FPS>
    FPS operator()(const FPS& a, const FPS& b) const {
        if (a.empty())
            return b;
        if (b.empty())
            return a;
        FPS res(a.size() + b.size() - 1);
        for (usize i = 0 ; i < a.size() ; i++)
            for (usize j = 0 ; j < b.size() ; j++)
                res[i + j] += a[i] * b[j];
        return res;
    }
};

} // namespace zawa
#line 4 "Src/FPS/BostanMori.hpp"

namespace zawa {

template <concepts::IndexedFPS FPS, class Conv = FPSMult>
requires concepts::Convolution<FPS, Conv>
typename FPS::value_type BostanMori(usize N, FPS P, FPS Q, Conv conv = {}) {
    assert(P.size());
    assert(Q.size() and Q[0] != 0); 
    auto takeParity = [&](const FPS& f, usize p) {
        FPS res;
        res.reserve(f.size() / 2);
        for (usize i = p ; i < f.size() ; i += 2)
            res.push_back(f[i]);
        return res;
    };
    while (N) {
        FPS Qm(Q.size());
        for (usize i = 0 ; i < Q.size() ; i++)
            Qm[i] = i % 2 ? -Q[i] : Q[i];
        P = takeParity(conv(P, Qm), N % 2);
        Q = takeParity(conv(Q, Qm), 0);
        N >>= 1;
    }
    return P[0] / Q[0];
}

} // namespace zawa

Back to top page