Day 042-Q3 — Bostan-Mori アルゴリズム（線形漸化式 N 項目）

2026-05-25 赤色 Master / Phase 8+ ★★★★★★★★★ Bostan-Mori / 有理母関数係数抽出

問題

数列 $a_n$ が線形漸化式 $a_n = c_1 a_{n-1} + \cdots + c_d a_{n-d}$ で定義される。初期値と係数、$N$ が与えられたとき $a_N \bmod (10^9+7)$ を求めよ。

制約

$1 \le d \le 500$
$0 \le N \le 10^{18}$
$0 \le a_i, c_i < 10^9+7$
時間制限: 2sec / メモリ: 256MB

入出力例

入力例 1

3 10
0 1 1
0 1 1

出力例 1

55

フィボナッチ数列の $a_{10} = 55$

概念図: Bostan-Mori の半分化

$[x^N]\ P/Q$ deg P < deg Q = d × Q* $U = P \cdot Q^*$, $V = Q \cdot Q^*$ V の奇数次 = 0（$Q \cdot Q(-x)$ の性質） N 偶奇で U の偶/奇次抽出 $[x^{\lfloor N/2 \rfloor}]\ P'/Q'$ N が半分に → $O(\log N)$ 回の再帰で終了 Kitamasa: $O(d^2 \log N)$ 行列累乗: $O(d^3 \log N)$ Bostan-Mori: $O(d \log d \log N)$ ✨ （NTT使用時; 本実装は $O(d^2 \log N)$）

ヒント（段階的開示）

ヒント1: 方向性
線形漸化式の解は有理母関数 $F(x) = P(x)/Q(x)$ で表される。$[x^N] F(x)$ を求める Bostan-Mori アルゴリズムは Kitamasa 法より高速（NTT使用時）。
ヒント2: アプローチ
  1. 特性多項式 $Q(x) = 1 - c_1 x - \cdots - c_d x^d$ を構成
  2. 分子 $P(x)$ を初期値から計算（$P = Q \cdot F \bmod x^d$）
  3. Bostan-Mori ループ: $Q^*(x) = Q(-x)$ を掛けて $N$ を半分に
  4. $N = 0$ になったら $P[0]$ を返す
ヒント3: 実装骨格 
def bostan_mori(P, Q, N):
    while N > 0:
        Qs = [c * pow(-1, i, MOD) % MOD for i, c in enumerate(Q)]
        U = poly_mul(P, Qs)  # P * Q*
        V = poly_mul(Q, Qs)  # Q * Q* (奇数次 = 0)
        if N & 1:
            P = U[1::2]  # 奇数次のみ
        else:
            P = U[0::2]  # 偶数次のみ
        Q = V[0::2]      # 偶数次のみ (= Q' の係数)
        N >>= 1
    return P[0] % MOD

模範解答 (Python)

import sys
input = sys.stdin.readline
MOD = 10**9 + 7

def poly_mul(A, B):
    n = len(A) + len(B) - 1
    C = [0] * n
    for i, a in enumerate(A):
        if a == 0: continue
        for j, b in enumerate(B):
            C[i+j] = (C[i+j] + a * b) % MOD
    return C

def bostan_mori(P, Q, N):
    """[x^N] P(x)/Q(x) mod MOD"""
    P = [x % MOD for x in P]
    Q = [x % MOD for x in Q]
    while N > 0:
        # Q*(x) = Q(-x): 奇数次の係数を符号反転
        Qs = [c * pow(-1, i, MOD) % MOD for i, c in enumerate(Q)]
        U = poly_mul(P, Qs)
        V = poly_mul(Q, Qs)
        # N の偶奇で P を選択
        if N & 1:
            P = U[1::2]
        else:
            P = U[0::2]
        Q = V[0::2]  # V の偶数次のみ
        N >>= 1
    return P[0] % MOD if P else 0

def solve():
    d, N = map(int, input().split())
    a = list(map(int, input().split()))
    c = list(map(int, input().split()))

    # 特性多項式: Q(x) = 1 - c_1*x - c_2*x^2 - ... - c_d*x^d
    Q = [1] + [(-ci) % MOD for ci in c]

    # 分子 P(x): P[n] = a[n] - sum_{k=1}^{n} c[k-1]*a[n-k]  (n < d)
    P = [0] * d
    for n in range(d):
        P[n] = a[n] % MOD
        for k in range(1, n + 1):
            P[n] = (P[n] - c[k-1] * a[n-k]) % MOD

    print(bostan_mori(P, Q, N))

solve()

Step-by-Step 解説

1母関数の構成
線形漸化式の解 $F(x) = \sum_{n \ge 0} a_n x^n$ は $F(x) = P(x)/Q(x)$ の形（有理関数）。分母 $Q$ は特性多項式、分子 $P$ は $Q \cdot F \bmod x^d$。
2Bostan-Mori の核心
$[x^N] P/Q$ を求める。両辺に $Q^*(x) = Q(-x)$ を掛けると分母 $Q \cdot Q^*$ の奇数次が消え（実際には消えるよう設計）、$N$ が半分の同型問題に帰着。
3N の偶奇での分岐
$[x^N] U/V$（$V$ の奇数次 = 0）を計算するとき、$N$ が偶数なら $U$ の偶数次を $x^{N/2}$ に対応させ、奇数なら奇数次を使う。
4停止条件
$N = 0$ のとき $[x^0] P/Q = P[0]$（$Q[0] = 1$ なので $P[0]$ が答え）。

計算量

多項式乗算 (ナイーブ): $O(d^2)$ / ステップ
ステップ数: $O(\log N)$
合計: $O(d^2 \log N)$（本実装）
NTT を使えば: $O(d \log d \log N)$（最適版）

よくあるミス

ミス原因正しい書き方
$Q^*(x)$ の符号反転ミス$Q(-x)$ の定義の誤解奇数インデックスの係数を $-1$ 倍
N 偶奇での P の選択ミスoff-by-oneN&1 で U[1::2] vs U[0::2]
分子 P の計算誤り$Q \cdot F \bmod x^d$ の構成漸化式を逆用して各 $P[n]$ を計算
MOD での符号扱いPython の % は正値(-c) % MOD で正規化

次のステップ

  • 発展問題: NTT で多項式乗算を $O(d \log d)$ に高速化し $d \le 10^5$ 対応
  • 類題: Library Checker "Find Linear Recurrence" + "Sparse Matrix Determinant"
  • 応用: Berlekamp-Massey で漸化式を自動発見 → Bostan-Mori で $N$ 項目計算

自己評価