形式的冪級数(FPS)の inv,log,exp,pow の定数倍の軽いアルゴリズム

mod が NTT-friendly な場合について,形式的冪級数の inv,log,exp,pow を実装しました.

yosupo さんの Library Checker に投げてみたところ,それなりに定数倍が軽そうだったので,実装したアルゴリズムを紹介します.

より高速なアルゴリズムをご存じの方は教えてくださると嬉しいです!

乗法の逆元(inv)

問題

$a_0 \neq 0$ なる多項式 $f(x) = \sum_{i=0}^{D-1} a_i x^i$ と正の整数 $N$ が与えられる.$f(x)g(x) = 1 \pmod{x^N}$ を満たす $(N-1)$ 次多項式 $g$ を求めよ.

参考: Inv of Formal Power Series - Library Checker

アルゴリズム

$f(x)$ などを単に $f$ と書くことにします.

今,$f g_m = 1 \pmod{x^m}$ なる $(m-1)$ 次多項式 $g_m$ が求まっているとします.このとき,$(f g_m - 1)^2 = 0 \pmod{x^{2m}}$ が成り立ちます.左辺を展開することで,\[ f(2g_m - fg_m^2) = 1 \pmod{x^{2m}} \]を得ます.つまり,\begin{equation} \label{def-g2m} g_{2m} := ((2g_m - fg_m^2) \bmod{x^{2m}}) \end{equation}と定めれば,$f g_{2m} = 1 \pmod{x^{2m}}$ が成り立ちます.

一方,$g_1(x) := a_0^{-1}$ と定めると,$f g_1 = 1 \pmod{x}$ が成り立ちます.

これらより,帰納的に $g_1, g_2, g_4,\dots$ が計算でき,問題の答えが求まります.

高速化の工夫

式 \eqref{def-g2m} の計算の定数倍を軽くするために,以下の $3$ 点に注意します:

  • $g_{2m}$ の $[0, m)$ 次の項は既に求まっている($g_m$ に等しい)ので,再計算を避ける.$fg_m^2$ の $[m, 2m)$ 次の項だけ計算できればよい.
  • $f g_m = 1 \pmod{x^m}$ より,$f g_m$ には $[1, m)$ 次の項は存在しない.
  • $f g_m$ は $3m$ 次未満なので,$f g_m$ と $(f g_m) \bmod{x^{2m}}$ の $[m, 2m)$ 次の項は一致する.

すると,$g_m$ から $g_{2m}$ を求める部分は以下のように書けます:

  1. $f$ の $[0, 2m)$ 次の項だけ取り出したものを改めて $f$ とする.
  2. $h := ((f g_m) \bmod{x^{2m}})$ を計算する.
  3. $h$ の $[m, 2m)$ 次の項だけ取り出したものを改めて $h$ とする.
  4. $-hg_m$ を計算する.これと $g_{2m}$ の $[m, 2m)$ 次の項は一致する.

これは,長さ $2m$ の配列の FFT / IFTT(NTT / INTT)$5$ 回で実現できます.

対数関数との合成(log)

問題

$a_0 = 1$ なる多項式 $f(x) = \sum_{i=0}^{D-1} a_i x^i$ と正の整数 $N$ が与えられる.冪級数を \[\log f(x) := - \sum_{k=1}^\infty \frac{1}{k} (1 - f(x))^k \] と定める.$(\log f(x)) \bmod{x^N}$ を求めよ.

$(1 - f(x))^k$ には $k$ 次以上の項しかないため,$\log f(x)$ の定義式の無限和を $k = N-1$ までで打ち切り,冪級数ではなく多項式と見なしても答えは変わりません.

参考: Log of Formal Power Series - Library Checker

アルゴリズム

$f(x)g(x) = 1 \pmod{x^N}$ を満たす $(N-1)$ 次多項式 $g(x)$ を $\frac{1}{f(x)}$ と書きます.また,多項式 $f$ の微分を $f'$ と書きます.

手計算により,$\log f(x)$ と \[ \int f'(x) \frac{1}{f(x)}\ dx \] の $N$ 次未満の項は一致することが示せます(不定積分の定数項は $0$ とします).

この式に従って,inv,乗算,不定積分を組み合わせれば計算できます.不定積分の計算の際に,$1,\dots,N-1$ の(乗法の)逆元を $O(N)$ で前計算することに注意してください.

log の高速化のためには,計算量(の定数倍)が大きい inv を高速化するとよいです.

指数関数との合成(exp)

問題

$a_0 = 0$ なる多項式 $f(x) = \sum_{i=0}^{D-1} a_i x^i$ と正の整数 $N$ が与えられる.冪級数を \[\exp f(x) := \sum_{k=0}^\infty \frac{1}{k!} (f(x))^k \] と定める.$(\exp f(x)) \bmod{x^N}$ を求めよ.

log のときと同様に,$\exp f(x)$ の定義式の無限和を $k = N-1$ までで打ち切って,多項式と見なしてもよいです.

参考: Exp of Formal Power Series - Library Checker

アルゴリズム

以下の論文のアルゴリズムが,定数倍の軽さと実装の簡単さを両立していて良さそうでした:

Alin Bostan, Eric Schost. A simple and fast algorithm for computing exponentials of power series. Information Processing Letters, 109 (13), pp.754-756, Elsevier, 2009.
arXiv: https://arxiv.org/pdf/1301.5804.pdf

Figure 1. の右側が,提案されているアルゴリズムです.

高速化の工夫

上の論文では,Step 2.a' で

  • 長さ $2m$ の配列の FFT / IFTT を $2$ 回
  • 長さ $m$ の配列の FFT / IFTT を $1$ 回

実行すると書かれています.

しかし,Step 2.a' が inv 計算に対応していることに注意して,上で inv を高速化したのと同様にすれば,

  • 長さ $2m$ の配列の FFT / IFTT を $1$ 回
  • 長さ $m$ の配列の FFT / IFTT を $3$ 回

にできます.こちらの方が少しだけ高速だと思います(何か勘違いしていたら教えてください).

冪乗(pow)

問題

多項式 $f(x) = \sum_{i=0}^{D-1} a_i x^i$ と正の整数 $N, M$ が与えられる.$(f(x))^M \bmod{x^N}$ を求めよ.

参考: Pow of Formal Power Series - Library Checker

アルゴリズム

基本的には,$(f(x))^M$ と \[ \exp(M (\log f(x) \bmod{x^N}) ) \] の $N$ 次未満の項が一致することを使います.ただし,これは $a_0 = 1$ の場合しか使えません.

$a_0 \neq 1$ の場合は,定数項が $1$ の多項式 $g$ を用いて,$f(x) = c x^l g(x)$ と書き直し,\[(f(x))^M = c^M x^{lM} (g(x))^M \] を使えばよいです.

このように,pow は log と pow を組み合わせれば計算できます.

pow の高速化のためには,定数倍が大きい exp を高速化するとよいです.

実装例

この記事で紹介したものが全部入っています.ご自由にお使いください:
https://judge.yosupo.jp/submission/42413

タイトルとURLをコピーしました