FPS composition f(g(x)) の実装について

01_guide.md

02.md は、 noshi91 さんによる発表 に始まるアルゴリズム群に対する、 Library Checker 提出 #199636 に至るまでの効率化のカジュアルなメモです。(当初は ScrapBox で公開、最終更新 : 2024-03-29)

ぜひ maspy さんによってもっと整理された投稿「 FPS 合成・逆関数の解説 (1) 逆関数と Power Projection」を参考にしてください。

02.md

FPS composition f(g(x)) の実装について

noshi91 さんがこの間発表した↓

https://noshi91.hatenablog.com/entry/2024/03/16/224034

のアルゴリズムを頑張って実装してみたが、意味不明ソースコードになってしまって読むのに適さない。写経する、とかの方法もあるが、きっと各々のライブラリの都合でうまく行かないだろう。

inverse https://judge.yosupo.jp/submission/197749

composition https://judge.yosupo.jp/submission/197722

まあ、（ multipoint eval とか Bostan--Mori のように）じきに定数倍まで研究されて発表されるんだろうけども、今実装したい人もいるでしょう、ということで

misc

$2$ 変数だが、 $1$ 次元配列で管理してもそんなに支障は出ないだろう、ただし FFT の削減のために、 $x$ の添え字を内側にとる。今後のさらなる FFT の削減のことを考えると、 $2$ 次元 FFT にしたほうがよさそうかな～ （ 2 次元配列はすごくメモリ食うのでだめ！（過激な主張））

$\lbrack x^k \rbrack \dfrac{g(x)}{1 - y f(x)}$ が計算できるということだから、 $g(x)=1$ の場合は代わりに $g(x)=x^c$ で計算することで $k$ に融通が効く。

$f(0)\neq 0$ でも計算できる ( https://twitter.com/maspy_stars/status/1769050170694779100 ) 。ただし、計算を進めて $\dfrac{p(y)}{q(y)}$ の形にしたとき $q(y)$ が定数ではないので注意。

FFT の削減

https://noshi91.hatenablog.com/entry/2023/12/10/163348

AtCoder Library 等の FFT の実装（よくある実装 ・「バタフライ演算」）では、 FFT した後の列が分かっているときにもとの列（を多項式としてみたとき）の偶数次の係数を取り出すのが非常に簡単。線形漸化式のための Bostan--Mori 法の高速化で重要だった方法である。一方で奇数次の係数を取り出したいときは、 FFT するときにずらしてやればいい。

上位桁の切り捨てがあるので線形漸化式のときほどうまくはいかないが、少なくとも、偶数次を先に取り出せば inverse FFT の長さが半分で済む。

$\lbrack x^k \rbrack \dfrac{g(x)}{1 - y f(x)}$ ではなく、 $y$ の次数について reverse して $\lbrack x^k \rbrack \dfrac{g(x)}{y - f(x)}$ の気持ちになってみる（この式変形は全然おかしいんだけど）。分母の計算で $y$ の次数が $2^n$ になって嫌な気持ちになるが、最高次だけ気を付けてやれば FFT のサイズがぴったりになる。つまり、巡回畳込みで $y^{2\cdot 2^n}$ と $1$ が同一視されてしまっても処理できる。分子は気を付けなくても普通に収まる。

ここで 2 次元 FFT の使用を仮定する。 $x$ , $y$ それぞれでゼロ埋めの領域があるので、それを利用して FFT を削減できる。特に、片方の変数について FFT するときはもう一方の変数のゼロ埋め領域は変換しなくてもよい。強いて言うなら、 $x$ と $y$ のどちらの軸で先に FFT するかで FFT の総量が変動して、それを考慮しようとすると大変。

         int fg = (d <= e ? 1 : 0);     // x 軸と y 軸のどちらが長いか？
         int m = n/2;
         for(int t=0; t<2; t++){
             if(fg^t){     // 長いほうをさきに処理する : 先にやるほうは長さが半分でよい
                 ibat(q, d, e, 1);
                 if(d == 1) break;
                 for(int i=0; i<d; i++) q1[i] = q[i];
                 for(int i=0; i<m; i++) q1[d+i] = Zero;
                 for(int i=d; i<m; i++) q1[m+i] = q[i];
                 bat(q1, d, e*2, 1);
                 e *= 2; m *= 2;   // 各軸の長さを更新
             } else {
                 for(int i=0; i<m*2; i+=d*2){
                     for(int j=0; j<d; j++) q1[i+j] = q[i/2+j];
                     for(int j=0; j<d; j++) q1[i+d+j] = Zero;
                 }
                 bat(q1, 1, d*2, e);
                 d *= 2; m *= 2;   // 各軸の長さを更新
             }
             std::swap(q, q1); std::swap(q1, q2); // バッファを切り替える
         }

転置原理

$g$ を事前情報としたとき $f(g(x))$ の計算を表す行列は $A_{i,j}=\lbrack x^i\rbrack (g(x))^j$ で、これを転置してコネコネするとべき乗係数列挙になる

FFT の転置は、「先頭以外を reverse して FFT 」である。（よくある実装の） FFT をした後でも、区間 $\lbrack 1,2),\lbrack 2,4),\lbrack 4,8) \cdots$ をそれぞれ reverse すれば得られる。

要素の並べ替えは普通に行列を書いて転置すればイメージしやすそう。

分母は、計算される順番と使う順番が逆だから、分母に FFT （の転置）を実行したものを全部スタックに積んでおく。（分母の FFT を転置しておけば、普通に畳みこんだときに middle product が計算できる。）