はじめに
多項式GCDは近年の形式的冪級数関連のアルゴリズムの隆盛に伴い競技プログラミングの文脈にて語られる事が徐々に増えてきた超高度典型の一つです。
応用には 多項式の多項式mod上での逆元 や、F(x)=P(x)/Q(x) となる F(x) から P, Q を推定する問題、 因数分解の乱択アルゴリズム 等があり、また、個人的には一般に部分分数分解が $\mathcal{O}(N(\log N)^3)$ で出来るのではないかと思ってます。 部分分数分解が O(N(log N)^2) で出来ます(37zigen さんありがとうございます！)。$10^{1000000}$ 位の大きな自然数同士の $\gcd$ の計算にも同様の手法が使われます。
$F(x)=P(x)/Q(x)$ となる $F(x)$ から $P, Q$ を推定する問題に関してもこの記事に書きたかったのですが、難航したため、お蔵入りを避けるためにもこれを省いて暫定版として出します。
前提知識
多項式除算 の $\mathrm{O}(N\log N)$ での計算を断りなく使用します
【競技プログラミング】形式的冪級数の応用テクニック(前編)
拡張ユークリッドの互除法の考え方を使用しています。
拡張ユークリッドの互除法 〜 一次不定方程式 ax + by = c の解き方 〜
多項式GCD とは
多項式GCDは２つの多項式を両方割り切る次数最大の多項式です（複数存在する場合はどれか一つ）。例えば、 $\gcd(1+2x+x^2,2+3x+x^2)=1+x$ となります。（正確には実数係数等の場合）。こういった物を高速に求めたいというのがこの記事の目標です。
ところで自然数同士の GCD は一般的には以下の様な実装になります。これを多項式について適用するのはどうでしょうか？
int gcd(int s, int t){
if(t==0)return s;
else return gcd(t,s%t);
}
結論から言うと正しい結果が得られます。以上でこの記事は終了です、とはいきません。正しい計算結果は得られるのですが、 $s\% t$ の次数は $t$ の次数未満の値になるため、単にこれを使うと $\gcd(s,t)$ にかかる計算量は $N=\max(\deg(s),\deg(t))$ として、 $\mathcal{O}(N^2\log N)$ となります。
ここで、 half GCD を用いてこの計算を $\mathcal{O}(N(\log N)^2)$ に落とすことが目標となります。
half GCD の主な方針
$30002$ と $20004$ の gcdを考えます。ここで、値が約 3:2 になっている事がわかると思います。
ここで、
$30010 \times 1 - 20004 \times 1 = 10006$
$30010 \times 2 - 20004 \times 3 = 8$
$\gcd(30010, 20004)=\gcd(10006, 8)= \gcd(8,6) = 2$ と求められる。
これを分割統治法でやっていこうというのが half GCD の主な方針です。
ところで、 $\gcd(30010, 20004)=\gcd(10006, 8)$ というのは本当でしょうか？実は $\gcd(3,2)$ の操作にそってやる（同じ変換行列を適用する）限りこれは本当です。
$\gcd$ のアルゴリズムはある数 $c$ を用いて $\gcd(x, y)=\gcd(y, x-cy)$ と変形する操作の繰り返しによって構成されますが、これは $c$が $x\%y$ と等しくなくても 常に成り立つためです。
補助関数
以下の様な補助関数 $\mathrm{hgcd}$ を作る事を考えます。
多項式 $f,g$ ( $\deg(f) \geq \deg(g)$ )に対し、 $c=af+bg$ を満たすような多項式 $a,b,c$ を返す。ただし、 $a,b,c$ は次数が $\frac{1}{2}\max(\deg(f)) + \mathcal{O}(1)$ 以下となるような多項式。
step 1. 補助関数により次数を約 3/4 倍に
$$
f=f_0+f_1 X^{\deg(f)/2},~g=g_0+g_1 X^{\deg(f)/2}
$$
とおき、 $\mathrm{hgcd}(f_1,g_1)$ を呼び出すことで $c=af_1+bg_1$ となる多項式 $a,b,c$ を得ます。
すると、
$$
af+bg=(af_0+bg_0) +(af_1+bg_1)X^{\deg(f)/2}=(af_0+bg_0)+cX^{\deg(f)/2}
$$
となりますが、よく見てみるとこれは
$$
\deg(af+bg) \leq \frac{3}{4}\deg(f) + \mathcal{O}(1)
$$
となることがわかります。ここで、
$$
\gcd(f,g)=\gcd(g,af+bg)= \gcd(af+bg, g\%(af+bg))
$$
と式変形することで両者の次数を $\frac{3}{4}\deg(f) +O(1)$ にすることが出来ました。
step 2. 補助関数により次数を約 1/2 倍に
step 1. とほぼ同じです。処理が重いので、この時点で 1/2倍になっていたら返す高速化を入れると良いです。
$p=af+bg,q=g\%(af+bg)$ と置き直し、
$p=p_0+p_1 X^{\deg(f)/4}$ , $q=q_0+q_1 X^{\deg(f)/4}$ として、もう一度 $\mathrm{hgcd}(p_1,q_1)$ を呼び出すことで $c=ap_1+bq_1$ となる多項式 $a,b,c$ を得ます。
すると、
$$
ap+bq=(ap_0+bq_0) +(ap_1+bq_1)X^{\deg(f)/4}=(ap_0+bq_0)+cX^{\deg(f)/4}
$$
となりますが、よく見てみるとこれは
$$
\deg(af+bg) \leq \frac{1}{2} \deg(f) + \mathcal{O}(1)
$$
となることがわかります。ここで、
$$
\gcd(f, g)=\gcd(p,q)=\gcd(g,ap+bq)= \gcd(ap+bq, q\%(ap+bq))
$$
と式変形することで両者の次数を $\frac{1}{2}\deg(f) +O(1)$ にすることが出来ました。
計算量解析
$N=\deg(f)$ となる $f,g$ に対し、 $\mathrm{hgcd}(f,g)$ の時間計算量を $T(N)$ とします。
すると、この関数からは次数半分の $\mathrm{hgcd}$ が二回、多項式除算が定数回呼ばれるので
$$
T(N) = 2T(N/2) + \mathcal{O}(N\log N)
$$
となります。これを展開して考えると、
$$
T(N) = \mathrm{O}(N\log N)+\mathrm{O}(N\log \frac{N}{2})+ \dots +\mathrm{O}(N\log(\frac{N}{2^K})) + \dots = \mathrm{O}(N(\log{N})^2)
$$
となりました。
多項式 GCD
$\mathrm{hgcd}$ を一度呼ぶごとに $\gcd(f,g)=\gcd(p,q)$ となる次数が半分以下の多項式になるため、 $\log(\deg(f))$ 回呼び出せば $\gcd(f,g)$ が求まります。この計算量 $T(N)$ は
$$
T(N) = T(N/2) + \mathcal{O}(N(\log N)^2)
$$
より、 $\mathrm{hgcd}(N)$ と同様 $\mathrm{O}(N(\log{N})^2)$ となりました。
TODO
Find Linear Recurrence は変換行列を適切な所まで計算する事で出来るらしいが、最小を飛び越えてしまう場合がない理由にいまいち自信がないので理解したら書きます...
参考文献
一番詳しい
https://codeforces.com/blog/entry/101850
殆どこの記事と同じことが書いてある(この記事を書く意味...)
https://scrapbox.io/37zigen-43465887/half-GCD
その他 half GCD でググると出てくる物
verify
half GCDで解ける問題
https://judge.yosupo.jp/problem/inv_of_polynomials
応用すると解けると思われる問題(ACしたが正当性に自信ない)
https://judge.yosupo.jp/problem/find_linear_recurrence