多項式求根

概要

多項式 $p(x) = c_0 + c_1 x + \cdots + c_n x^n$ の全根 (複素数含む) を求める問題は、汎用の 1D 求根とは別の専用手法群が発展している。 次数 $n$ に応じて閉形式解・反復法・固有値問題への帰着が使い分けられる。

関連 API: solveLinear, solveQuadratic, solveCubic, solveQuartic, jenkinsTraub, laguerre, durandKernerAberth, bairstow。

低次の閉形式解 (1〜4 次)

$n \leq 4$ では公式解が存在する。 Galois 理論により $n \geq 5$ では一般に閉形式解は存在しない (Abel-Ruffini)。

2 次方程式の桁落ちは古典的な数値解析の例で、sangi は $a \approx 0$ や $b^2 \gg 4ac$ で逆 Vieta 公式を使う実装を採用する。

Companion 行列 + QR (Francis 法)

monic 多項式 $p(x) = x^n + c_{n-1} x^{n-1} + \cdots + c_0$ の零点は、次のcompanion 行列の固有値と一致する:

$$C(p) = \begin{pmatrix} 0 & 0 & \cdots & 0 & -c_0 \\ 1 & 0 & \cdots & 0 & -c_1 \\ 0 & 1 & \cdots & 0 & -c_2 \\ \vdots & & \ddots & & \vdots \\ 0 & 0 & \cdots & 1 & -c_{n-1} \end{pmatrix}$$

$\det(xI - C(p)) = p(x)$ なので $C(p)$ の固有値 = $p$ の零点。 この帰着で多項式求根を固有値問題として扱える。

Francis 法 (QR 反復)

QR 反復はシフト戦略 (Wilkinson, Francis double shift) と組み合わせて二次・三次収束し、$O(n^3)$ で全固有値を計算する。 NumPy の np.roots や MATLAB の roots() はこの経路を採用する標準実装である。

長所: 数値的に成熟、複素根も一括

短所: $O(n^3)$、係数のスケーリングが必要、近接根で精度低下

sangi は companion + QR を直接公開していないが、Eigen 等の線形代数バックエンドと組み合わせれば呼び出せる。 汎用経路は Jenkins-Traub と Aberth-Ehrlich が中心。

Durand-Kerner (Weierstrass)

全 $n$ 根を同時反復する手法。初期推定 $z_1^{(0)}, \ldots, z_n^{(0)}$ から:

$$z_i^{(k+1)} = z_i^{(k)} - \frac{p(z_i^{(k)})}{\prod_{j \neq i} (z_i^{(k)} - z_j^{(k)})}$$

分母は他の根との差の積で、$z_i = \alpha_i$ (真の根) のとき右辺 = 0 が固定点になる。 Weierstrass (1891) が発見し、Durand (1960) と Kerner (1966) が独立に再発見した。

収束次数: 単根に対して 2 次

初期値の標準的な選び方は単位円上の等分点 $z_i^{(0)} = \beta \cdot e^{2\pi i (i-1)/n + \theta}$ で $\beta$ は係数から決まるスケール。

Aberth-Ehrlich

Aberth (1973) と Ehrlich (1967) による Durand-Kerner の改良。 Newton 法と Durand-Kerner を組み合わせ、より高い収束次数と安定性を得る:

$$z_i^{(k+1)} = z_i^{(k)} - \frac{N_i}{1 - N_i \cdot \sum_{j \neq i} (z_i^{(k)} - z_j^{(k)})^{-1}}, \quad N_i = \frac{p(z_i^{(k)})}{p'(z_i^{(k)})}$$

$N_i$ は標準的な Newton ステップで、補正項 $\sum_{j \neq i} (z_i - z_j)^{-1}$ が他の根からの「斥力」を表す。

収束次数: 単根に対して 3 次

全根同時反復の利点

Jenkins-Traub のような縮約 (deflation) 系手法は、1 根を取り出すたびに $p(x) / (x - \alpha)$ を計算して次数を下げる。 浮動小数では縮約が元の多項式と数値的にずれていくため、後半の根ほど精度が落ちる。 Aberth は元の $p(x)$ で全反復評価するため、縮約起因の誤差累積がない。

並列化との相性も良く、近年は MPSolve のような実装で大規模多項式 (次数 $10^4$ 以上) の標準手法になっている。 sangi の durandKernerAberth はこの手法を採用する。

Jenkins-Traub

Jenkins と Traub (1970) による 3 段階反復法。実係数版 (RPOLY) と複素係数版 (CPOLY) がある。 数値的安定性が極めて高く、ACM Algorithm 419/493 として配布されている。

3 段階の構成

1. Stage 1 (no-shift): $H$-多項式の無シフト反復で大局的な情報を集める。$M$ ステップ ($M \approx 5$)。
2. Stage 2 (fixed-shift): 1 つの根の近似値 $s$ を固定して反復し、$H$ の方向を決定する。$L$ ステップ ($L \approx 9 + n/2$)。
3. Stage 3 (variable-shift): シフト $s$ を毎ステップ更新し、二次収束に乗せる。

$H$-多項式は $H^{(j+1)}(x) = (H^{(j)}(x) - H^{(j)}(s) p(x)/p(s)) / (x - s)$ で更新され、徐々に「不要な根」を取り除いて目標根に集中する。 Stage 3 で目標根の収束を確認後、その根を縮約で除いて次の根に進む。

収束次数: Stage 3 で約 1.62 (3 段階全体の漸近)

実装は複雑だが数値的に成熟しており、sangi の jenkinsTraub はsolvePolynomial の 5 次以上のデフォルト経路。

Laguerre 法

Laguerre 法は多項式専用の三次収束法で、Newton-Halley 族の発展形:

$$x_{n+1} = x_n - \frac{n \cdot p(x_n)}{G(x_n) \pm \sqrt{(n-1)(n \cdot H(x_n) - G(x_n)^2)}}$$

ここで $G(x) = p'(x)/p(x)$, $H(x) = G(x)^2 - p''(x)/p(x)$。 平方根の符号は分母の絶対値を大きくする側を選ぶ。

大域収束性

Laguerre 法は単根に対して三次収束し、ほぼすべての初期値から収束する大域収束性を持つ (Adams-Smith の結果)。 Newton や Halley が局所的にしか収束しないのと対照的。 実係数多項式の複素根を求めるとき、初期値の選択に悩む必要が少ない。

sangi の laguerre はsolvePolynomial の選択肢のひとつ。

Bairstow 法

実係数多項式の複素根を実数演算のみで求める手法。 多項式を 2 次因子 $(x^2 - rx - s)$ で割り、剰余 $b_1 x + b_0$ が 0 になるよう $(r, s)$ を Newton 反復で更新する。

$$\begin{pmatrix} \partial b_0/\partial r & \partial b_0/\partial s \\ \partial b_1/\partial r & \partial b_1/\partial s \end{pmatrix} \begin{pmatrix} \Delta r \\ \Delta s \end{pmatrix} = - \begin{pmatrix} b_0 \\ b_1 \end{pmatrix}$$

2 次因子から複素共役な根のペアが直接得られ、複素数演算を一切しなくて済む。 ハードウェアが複素演算を持たない・複素実装したくない場合に有用 (組込み等)。 現代の汎用環境では Jenkins-Traub や Aberth に劣るが、sangi は bairstow として用意する。

Sturm 列と実根計数

多項式 $f$ から Sturm 列 $f_0 = f, f_1 = f', f_{i+1} = -(f_{i-1} \bmod f_i)$ を構成する。 $a$ での Sturm 列の符号変化数 $v(a)$ を使うと、区間 $[a, b]$ 内の異なる実根の個数が:

$$\#\{\alpha \in [a, b] : f(\alpha) = 0\} = v(a) - v(b)$$

で厳密に得られる (Sturm の定理)。 重根も 1 回だけ数える。

用途

• 実根の個数の厳密な事前カウント
• 区間二分による実根の分離 (各小区間に 1 根だけ)
• 分離後に 1 次元囲い込み法 (Brent 等) で精密化

sangi の signVariations(p) は Descartes の符号則による上界を返し、 vincentRealRootIsolation は Vincent-Akritas-Strzeboński の連分数アプローチで実根を分離する (Sturm より高速)。

sangi のディスパッチ判断

solvePolynomial(p) は次のように経路を選ぶ:

• 次数 1〜4: 閉形式解 (solveLinear/Quadratic/Cubic/Quartic)
• 次数 5 以上: jenkinsTraub

Jenkins-Traub をデフォルトにする理由:

• 実装が成熟し、係数のスケーリング・劣条件係数への耐性が確認されている
• RPOLY/CPOLY として 50 年以上の運用実績
• 5 次〜数百次の多項式で精度・速度のバランスが安定

大規模次数 (数千以上) や並列環境では Aberth-Ehrlich が有利になる。 重根や近接根が多いケース、または初期値の選択に悩む場合は Laguerre 法の大域収束性が有用。 実係数で複素演算を避けたい場合は Bairstow。 利用者が手動で経路を選ぶ場合は API で jenkinsTraub, laguerre, durandKernerAberth, bairstow を直接呼び出せる。

係数型 (Int / Rational) の場合

厳密な根が代数的に表現できる場合 (有理根定理が適用できる場合) は、因数分解パイプラインで 1 次因子を取り出してから求根を行う。 浮動小数で全根を求めたい場合は double または Complex<double> 経由の数値求根に委譲する。

参考文献

• Jenkins, M. A. and Traub, J. F. (1970). "A three-stage algorithm for real polynomials using quadratic iteration". SIAM Journal on Numerical Analysis, 7(4), 545–566.
• Aberth, O. (1973). "Iteration methods for finding all zeros of a polynomial simultaneously". Mathematics of Computation, 27(122), 339–344.
• Bini, D. A. and Robol, L. (2014). "Solving secular and polynomial equations: A multiprecision algorithm". Journal of Computational and Applied Mathematics, 272, 276–292.
• Adams, D. A. (1967). "A stopping criterion for polynomial root finding". Communications of the ACM, 10(10), 655–658.