From 799d9837551e311c70ecb65358bc6097f73bdd0c Mon Sep 17 00:00:00 2001
From: ganmodokix <1450602+ganmodokix@users.noreply.github.com>
Date: Sun, 14 Dec 2025 22:27:28 +0900
Subject: [PATCH] =?UTF-8?q?=E2=9C=A8=20Add=20Bostan-Mori=20algorithm?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 template/fps/bostan_mori.hpp | 52 ++++++++++++++++++++++++++++++++++++
 1 file changed, 52 insertions(+)
 create mode 100644 template/fps/bostan_mori.hpp

diff --git a/template/fps/bostan_mori.hpp b/template/fps/bostan_mori.hpp
new file mode 100644
index 0000000..e60b68d
--- /dev/null
+++ b/template/fps/bostan_mori.hpp
@@ -0,0 +1,52 @@
+#pragma once
+#include "base_template.hpp"
+#include "modint/modint.hpp"
+#include "conv/ntt.hpp"
+
+// Bostan-Mori [x^m] f(x) / g(x) を求める O(NlogNlogM) (N: deg g)
+// verified at https://judge.yosupo.jp/submission/336936
+// verified at https://atcoder.jp/contests/abc436/submissions/71731332
+template <mod_integral T>
+constexpr T bostan_mori(vector<T> f, vector<T> g, ll m) {
+    assert(m >= 0);
+    assert(!g.empty());
+    assert(g[0] != 0);
+    assert(ssize(f) < ssize(g));
+    if (m == 0) {
+        if (f.empty()) return 0;
+        return f[0] / g[0];
+    }
+    // P(x)Q(-x) = A(x^2) + x B(x^2),
+    // Q(x)Q(-x) = R(x^2) ∵ Q(x)Q(-x) は偶
+    // 以上のように分解すると
+    // P(x)/Q(x) = P(x)Q(-x) / Q(x)Q(-x) = A(x^2)/R(x^2) + x B(x^2)/R(x^2)
+    // FPS が x^2 のみからなる <=> 偶 <=> 奇数次数の係数がすべて 0
+    // なので
+    // [x^m] P(x)/Q(x)
+    // = [x^{m/2}] A(x)/R(x) if m は偶数
+    // = [x^{(m-1)/2}] B(x)/R(x)
+
+    // ステップごとにmが半減するのでダブリングの要領で O(NlogNlogM)
+    // P,Qの次数は畳み込んで半減なので±1を除いて変わらず
+    
+    // 分母
+    auto gnx = g;  // g(-x)
+    STEP(i, 1, ssize(gnx) - 1, 2) {
+        gnx[i] *= -1;
+    }
+    auto rx2 = convolve_p(move(g), gnx);
+    auto r = vector<T>{};
+    r.reserve((rx2.size() + 1) / 2);
+    STEP(i, 0, ssize(rx2) - 1, 2) {
+        r.emplace_back(rx2[i]);
+    }
+
+    auto ab = convolve_p(move(f), move(gnx));
+    auto ab2 = vector<T>{};
+    ab2.reserve((ssize(ab) + (m + 1) % 2) / 2);
+    STEP(i, m % 2, ssize(ab) - 1, 2) {
+        ab2.emplace_back(ab[i]);
+    }
+
+    return bostan_mori(move(ab2), move(r), m / 2);
+}