公式解説 の内容を前提とし、\(x, y\) についての \(2\) 変数 FPS \(f\) を以下のように定義します。

\[ f \coloneqq \sum_{T \geq 0}y^T (VW)^T_{1,1}\]

すると、問題は \([x^N]f\) の \(U\) 次までの係数の計算と表現することができます。

\(f = (I - yVW)^{-1}_{1,1}\) より、\(f\) は低次の \(2\) 変数有理式であることが分かります。このような場合、\([x^N]f\) は \(N\) をパラメータとする低次の微分方程式を 持つことが知られています [1]。[追記: 今回のケースは参照した論文の適用条件を満たしていませんでした。ただし、発見した漸化式の正当性は数式処理で確認できます]。数式処理ライブラリなどで探索を行うと、

\[ (T-3)(3T-5)(2T - N)(2T-N-1) [x^N y^T] f = \sum_{1 \leq i \leq 5} p_i(N, T) [x^N y^{T-i}]f \]

という漸化式を得ることができます。ここで \(p_i\) は低次の \(2\) 変数多項式です。よって、左辺の係数が \(0\) にならない範囲ではこの式によって答えを順に計算することができます。

\(T = 3\) のケースは、\(N \gt 6\) なら答えが \(0\) になるため \(N \le 6\) の範囲を埋め込めば十分です。\(N = 2T, 2T-1\) の場合は単純な数え上げで答えが求まるため、これで全ての答えを計算することができます。

時間計算量は \(\Theta(U)\) です。

実装例

https://atcoder.jp/contests/fps-24/submissions/70671886

• \([x^N]f\) を \(N!\) 倍された状態で直接計算しています。
• 実装で手を抜いているため時間計算量が \(\Theta(U \log P)\) になっています

参考文献

[1] Bostan, A., Neiger, V., & Yurkevich, S. (2023, July). Beating binary powering for polynomial matrices. In Proceedings of the 2023 International Symposium on Symbolic and Algebraic Computation (pp. 70-79).