多項式 GCD と終結式

概要

多項式の GCD (最大公約因子) は因数分解・無平方分解・有理関数の約分・代数曲線の交点計算など、計算機代数の多くの操作の核となる。 体上 ($\mathbb{Q}[x]$, $\mathrm{GF}(p)[x]$ など) では Euclidean アルゴリズムが直接動くが、 整数環上 ($\mathbb{Z}[x]$) では中間係数の爆発を抑える特別な扱いが必要になる。

関連 API: gcd(p, q), intPolyGcd, resultant, discriminant, sylvesterMatrix。

Euclidean アルゴリズムと係数の爆発

体上の Euclidean アルゴリズムは整数の場合とまったく同じ形である:

$$r_0 = f,\; r_1 = g,\; r_{i+1} = r_{i-1} \bmod r_i, \;\gcd(f, g) = r_k \;(r_{k+1} = 0)$$

$\mathbb{Q}[x]$ なら問題なく動くが、$\mathbb{Z}[x]$ で素朴に擬剰余列を走らせると中間係数が急激に巨大化する。 Knuth は古典的な例として $f = x^8 + x^6 - 3x^4 - 3x^3 + 8x^2 + 2x - 5$, $g = 3x^6 + 5x^4 - 4x^2 - 9x + 21$ を挙げ、 最終 GCD は $\gcd(f, g) = 1$ (定数) であるにも関わらず中間係数が 5 桁程度まで膨らむことを示している。 一般に擬剰余列で中間係数の bit 長は次数に対して指数的に増える。

このため $\mathbb{Z}[x]$ の GCD には素朴な擬剰余列ではなく、後述するサブレザルタント PRS かモジュラー GCD を使う。

サブレザルタント PRS

サブレザルタント擬剰余列 (subresultant pseudo-remainder sequence, SR-PRS) は、 Collins (1967) と Brown (1971) によって整備された手法で、$\mathbb{Z}[x]$ 上の GCD を中間係数の指数爆発を伴わずに計算する。

アイデア

擬剰余 $r_{i+1} = \mathrm{prem}(r_{i-1}, r_i)$ は主係数の冪を含むため大きすぎる。 SR-PRS は各ステップで適切な「共通因子」$\beta_i$ で除算する:

$$r_{i+1} = \frac{\mathrm{prem}(r_{i-1}, r_i)}{\beta_i}$$

$\beta_i$ は前ステップの主係数 $\ell_i$ と次数差から閉形式で書ける ($\beta_i = -\ell_{i-1} \cdot \psi_i^{\delta_i}$ など)。 この共通因子で割っても結果は必ず整数係数に留まり、 中間係数の bit 長は入力サイズに対して多項式オーダーに抑えられる。

sangi の実装

intPolyGcd(a, b) はサブレザルタント PRS で実装されている。 pseudoRemainder(f, g) は SR-PRS の構成要素として公開されている (擬剰余そのもの)。 最終ステップで内容 (content) の GCD を取り、原始部 (primitive part) と合わせて $\mathbb{Z}[x]$ の GCD を得る。

計算量: 古典 SR-PRS は $O(n^2 M(n L))$ 程度 ($n$ は次数、$L$ は入力係数 bit 長、$M$ は多倍長乗算のコスト)

関連記事: 多項式の GCD 計算 — Euclidean / SR-PRS / モジュラー法の体系的解説

Sylvester 行列と終結式

2 つの多項式 $f = \sum_{i=0}^n a_i x^i$ と $g = \sum_{i=0}^m b_i x^i$ のSylvester 行列は $(m + n) \times (m + n)$ の正方行列で、上半部に $f$ の係数を $m$ 行ずらしながら詰め、下半部に $g$ の係数を $n$ 行詰めたものである:

$$\mathrm{Syl}(f, g) = \begin{pmatrix} a_n & a_{n-1} & \cdots & a_0 & & \\ & a_n & a_{n-1} & \cdots & a_0 & \\ & & \ddots & & & \ddots \\ b_m & b_{m-1} & \cdots & b_0 & & \\ & b_m & b_{m-1} & \cdots & b_0 & \\ & & \ddots & & & \ddots \\ \end{pmatrix}$$

終結式 $\mathrm{Res}(f, g)$ はこの行列の行列式として定義される:

$$\mathrm{Res}(f, g) = \det \mathrm{Syl}(f, g)$$

性質

  • $\mathrm{Res}(f, g) = 0 \;\Leftrightarrow\; f$ と $g$ が代数閉包の中で共通根を持つ $\;\Leftrightarrow\; \deg \gcd(f, g) \geq 1$
  • $\mathrm{Res}(f, g) = a_n^m \prod_i g(\alpha_i) = (-1)^{nm} b_m^n \prod_j f(\beta_j)$ ($\alpha_i, \beta_j$ はそれぞれ $f, g$ の根)
  • $f, g \in \mathbb{Z}[x]$ なら $\mathrm{Res}(f, g) \in \mathbb{Z}$ (係数のみで定まる)

計算

$\det \mathrm{Syl}$ の直接計算は $O((n+m)^3)$ で、係数が大きいとさらに重い。 漸近的に有利な経路としては SR-PRS の副産物として終結式を読み取る方法があり、 最終 PRS の主係数を、各ステップで掛けた係数の積で割ると終結式が得られる ($O((n+m)^2 \cdot$ 多倍長係数コスト$)$)。 sangi の現行 resultant(f, g) は Sylvester 行列を直接構築して Gauss 消去で行列式を計算する経路 ($O((n+m)^3)$) を採用しており、SR-PRS 経路への切替は今後の改善候補である。

関連記事: 終結式

判別式

$f$ の判別式は導関数との終結式から定義される:

$$\mathrm{Disc}(f) = \frac{(-1)^{n(n-1)/2}}{a_n} \mathrm{Res}(f, f')$$

2 次・3 次の馴染みの判別式 ($b^2 - 4ac$ など) の一般化である。

性質と用途

  • $\mathrm{Disc}(f) = 0 \;\Leftrightarrow\; f$ が代数閉包で重根を持つ
  • $\mathrm{Disc}(f)$ の正負で実根の個数 (の偶奇) が分かる場合がある
  • 因数分解の前処理で「無平方かどうか」を 1 回の終結式計算で判定できる

sangi の discriminant(f) は内部で resultant(f, derivative(f)) を呼び、符号と $a_n$ で割る。

// f(x) = x^2 - 2 の判別式 = -1·Res(f, 2x)/1 = ... = 8
Polynomial<int> f = {-2, 0, 1};
std::cout << discriminant(f) << std::endl;  // 8

Half-GCD (HGCD)

古典的な Euclidean / SR-PRS は次数 $n$ に対して各反復で次数が 1 ずつ落ちる場合に $O(n)$ 回反復し、 毎反復の擬除算が $O(n M(n L))$ なので全体は $O(n^2 M(n L))$ になる。 Half-GCD はこの $O(n)$ 反復を分割統治で $O(\log n)$ レベルに圧縮する手法である。

アイデア

2 つの多項式 $(r_0, r_1)$ の上半分の係数だけ見て、剰余列を「半分の次数」まで進めるユニモジュラー変換行列 $\begin{pmatrix} a & b \\ c & d \end{pmatrix}$ を再帰的に計算する。 変換行列は $(r_0, r_1) \to (r_k, r_{k+1})$ ($k \approx n/2$) の関係を保持し、 最終的に $O(M(n) \log n)$ で GCD と Bezout 係数が得られる。

計算量: $O(M(n) \log n)$ ここで $M(n)$ は多項式乗算のコスト

次数が数千以上で FFT 多項式乗算が効く領域で素朴 GCD を大きく上回る。 sangi の現行 GCD は SR-PRS 系を主軸とし、Half-GCD はオプション経路として将来拡張で用意する。

モジュラー GCD (Brown-Collins)

整数係数多項式の GCD を、複数の素数 $p_1, p_2, \ldots$ で有限体上の GCDとして計算し、 中国剰余定理 (CRT) で $\mathbb{Z}[x]$ に持ち上げる手法。Brown (1971) と Collins により整備された。

手順

  1. 適当な素数 $p$ を選ぶ ($p$ は $\mathrm{lc}(f), \mathrm{lc}(g)$ を割らない)
  2. $\mathrm{GF}(p)[x]$ 上で Euclidean アルゴリズムで $\gcd_p(f, g)$ を計算 ($O(n^2)$ in $\mathrm{GF}(p)$)
  3. 複数の素数で繰り返し、CRT で係数を統合
  4. 結果係数の絶対値が事前境界 (Mignotte bound) を超えなくなったら終了

利点と注意点

  • 係数爆発が起きない: 各 $\mathrm{GF}(p)$ 上の計算は係数が $< p$ に収まる
  • 並列化に向く: 異なる素数の計算は独立
  • 悪い素数: 一部の $p$ で次数が下がる場合があり、これらは捨てて再試行する (確率は低い)

係数 bit 長 $L$ と次数 $n$ が両方大きい場合に SR-PRS より速く、CAS の標準実装で用いられる。 sangi の intPolyGcd は現状 SR-PRS を主経路とし、Brown-Collins モジュラー GCD は計画中の拡張。

近似 GCD (浮動小数点向け)

係数が浮動小数 (double) の場合、厳密な共通因子は存在しないことが多い (係数の摂動で簡単に互いに素になる)。 近似 GCD はある許容誤差 $\epsilon$ の下で「実質的な共通因子」を抽出する手法。

sangi の approximateGCD(p, q, tol) は次の戦略を取る:

  • Sylvester 行列の数値階数 (SVD で特異値 $< \epsilon$ を切り捨て) から $\deg \gcd$ を推定
  • 推定された次数の多項式を最小二乗で再構成

本格的な信号処理応用では Karmarkar-Lakshman 法や ApaTools の手法が知られるが、 sangi は実装の単純さを優先して階数推定 + LS 再構成を採用している。

設計判断

SR-PRS を主経路にする理由

sangi の典型的ユースケース (因数分解の前処理、Yun 無平方分解、終結式) では次数 $n \lesssim 100$、係数 bit 長 $L \lesssim 64$ が大半である。 この領域では SR-PRS の定数倍が小さく、Half-GCD やモジュラー GCD の前処理コスト (FFT・複数素数選択・CRT) を上回るほどの利得はない。 SR-PRS 単独で実装の複雑度と性能のバランスが最良となる。

終結式の現行経路と将来の SR-PRS 化

Sylvester 行列式の直接展開 (Gauss 消去) は $O((n+m)^3)$ で、$\mathbb{Z}$ の中間係数も巨大化しやすい。 漸近的には SR-PRS の副産物として終結式を取り出すほうが安く、GCD 計算と終結式計算を 1 パスにまとめられる。 現行の sangi resultant は実装の単純さを優先して Sylvester+det 経路を採用しているが、 将来は SR-PRS 経路への切替で多倍長領域の性能を改善する予定である。

参考文献

  • Brown, W. S. (1971). "On Euclid's algorithm and the computation of polynomial greatest common divisors". Journal of the ACM, 18(4), 478–504.
  • Collins, G. E. (1967). "Subresultants and reduced polynomial remainder sequences". Journal of the ACM, 14(1), 128–142.
  • Knuth, D. E. The Art of Computer Programming, Volume 2: Seminumerical Algorithms. 3rd ed. Addison-Wesley, 1997.
  • von zur Gathen, J. and Gerhard, J. Modern Computer Algebra. 3rd ed. Cambridge University Press, 2013.