実根分離

数学、特に数値解析とコンピュータ代数において、多項式の実根分離（じつこんぶんり、英語: Real-root isolation）とは、実数直線の互いに素な区間を生成することであり、これらの区間は多項式の実根をそれぞれ1つだけ含み、また全体としては多項式のすべての実根を含む。

実根分離が有用なのは、多項式の実根を計算する通常の求根アルゴリズムは、いくつかの実根（の探索結果）を生成することはあっても、一般にすべての実根が探索されたことを証明できないためである。特に、そのようなアルゴリズムで根が探索されない場合、それが実根が存在しないからなのかどうかは分からない。いくつかのアルゴリズムはすべての複素根を計算するが、一般に実根は複素根よりもはるかに少ないため、計算時間の大部分は非実根の計算に費やされる（平均して、次数 n の多項式には n 個の複素根があり、実根は log n 個しかない。多項式の根の性質の『実根』（英語版）を参照）。さらに、虚部が小さい非実根と実根を区別することが困難な場合がある（次のセクションのウィルキンソンの多項式（英語版）の例を参照）。

最初の完全な実根分離アルゴリズムは スツルムの定理（1829年）から得られた。しかし、実根分離アルゴリズムがコンピュータに実装され始めたとき、スツルムの定理から導かれたアルゴリズムはデカルトの符号法則（1637年）から導かれたアルゴリズムよりも効率が悪いことがわかった。20世紀初頭以来、デカルトの符号法則から導かれるアルゴリズムの改良と、きわめて効率的な実装の実現、そしてそれらの計算複雑度（英語版）の決定に関する研究が盛んに行われてきた。最良の実装は、1,000次を超える多項式の実根を問題なく分離することができる。^[1][2]

仕様と一般的な戦略

多項式の実根を求める一般的な戦略は、実数直線（または根を求める区間）を、各区間に根が最大で1つになるまで互いに素な区間に分割することである。このような手順は根分離と呼ばれ、結果として得られる区間にちょうど1つの根が含まれる場合、その根の分離区間となる。

ウィルキンソンの多項式（英語版）は、多項式の1つの係数をわずかに変更するだけで、根の値だけでなく、根の性質（実数か複素数か）も劇的に変化する可能性があることを示している。また、良好な近似値を用いても、近似根で多項式を評価すると、結果がゼロに非常に近くなる場合がある。例えば、20次（ウィルキンソン多項式の次数）の多項式の根が10に近い場合、その根における多項式の導関数は ${\displaystyle 10^{20}}$ のオーダーになる可能性がある。これは、根の値に ${\displaystyle 10^{-10}}$ の誤差があると、近似根における多項式の値が ${\displaystyle 10^{10}}$ のオーダーになる可能性があることを意味する。ごく低い次数の場合を除いて、根の分離手法では正確な演算を用いなければ信頼できる結果を得ることができない、ということになる。したがい、浮動小数点係数を持つ多項式の根を分離したい場合には、しばしば、それらを有理数に変換し、得られた多項式の原始成分を取得して、整数係数の多項式を得るのが適切である。

このため、以下に説明する方法は理論的には実数に対して有効であるが、実際には整数係数の多項式や有理数で終わる区間に対して一般的に用いられる。また、多項式は常に平方無縁（英語版）であると想定されている。これには2つの理由がある。第一に、平方無縁因数分解（英語版）に対するYunのアルゴリズム（英語版）の計算コストは、多項式とその導関数の最大公約数の計算コストの2倍未満である。これにより低次の因数が生成できるため、アルゴリズムで必要とされない場合でも、多重根を持たない多項式に対してのみ根分離アルゴリズムを適用するのが一般に有利である。平方無縁多項式のみを考慮する2つ目の理由は、最高速の根分離アルゴリズムが多重根の場合は機能しないからである。

根の分離には、区間内の多項式の実根を計算せずに数える手順、あるいは少なくとも区間に根が0個、1個、あるいはそれ以上含まれているかどうかを判定する手順が必要である。このような判定手順を用いれば、実根を含む可能性のある区間の作業リストを扱うことができる。リストには、対象となるすべての根（通常は実数直線全体またはその正の部分）を含む単一の区間が含まれる。次に、リストの各区間は2つの小さな区間に分割される。新しい区間の1つに根が含まれていない場合、その区間はリストから削除される。根が1個含まれている場合、その区間は分離区間の出力リストに格納される。そうでない場合、その区間はさらなる分割のために作業リストに保持され、すべての根が最終的に分離されるまでこの処理が続けられる。

スツルムの定理

最初の完全な根分離手順はスツルムの定理（1829年）である。これは、区間内の実根の数を、区間の両端における「スツルム列」と呼ばれる多項式列の値の符号変化の数で表す。スツルム列は、多項式とその導関数にユークリッドの互除法の変種を適用した際に発生する剰余の列である。コンピュータに実装した場合、スツルムの定理を用いた根分離は、以下に説明する他の方法よりも効率が悪いことがわかった^[3]。そのため、スツルムの定理は机上目的には有用であるが、高効率な計算にはほとんど使用されていない。

デカルトの符号法則とその一般化

デカルトの符号法則によれば、多項式の係数列における符号変化の数と、その正の実根の数との差は非負の偶数である。したがって、この符号変化の数が0の場合、その多項式には正の実根が存在せず、この数が1の場合、その多項式には唯一の正の実根、つまり単一の根が存在する。残念ながら、逆は成り立たない。つまり、正の実根を持たない多項式、または正の単根を1つだけ持つ多項式は、1を超える符号変化の数を持つ可能性がある。

これはブダンの定理（1807年）によって一般化され、半開区間 (a, b] 内の実根についても同様の結果が得られた。f(x) が多項式で、v が f(x + a) と f(x + b) の係数列の符号変化の数の差である場合、v から区間内の実根の数（重複度を含む）を引いた値は非負の偶数になる。これはデカルトの符号法則の一般化である。なぜなら、b が十分に大きい場合、f(x + b) の係数に符号変化はなく、すべての実根は b より小さくなる。

ブダンの定理は、平方無縁多項式（英語版）（多重根を持たない多項式）に対する実根分離アルゴリズムを提供する。多項式の係数から、根の絶対値の上限 M と、2つの根の差の絶対値の下限 m を計算できる（多項式の根の性質（英語版）を参照）。次に、区間 [–M, M] を m 未満の長さの区間に分割すると、すべての実根はいずれかの区間に含まれ、2つの根を含む区間はない。したがって、分離区間はブダンの定理によって奇数の根が存在するとされる区間である。

しかし、このアルゴリズムはきわめて非効率的である。区間 [–M, M] の粗い分割を使用できないためある。これは、ブダンの定理がより大きなサイズの区間に対して1より大きい結果を与える場合、複数の根が含まれていないことを保証する方法がないためである。

ヴァンサンの定理と関連定理

ヴァンサンの定理（1834年）^[4]は、実根分離のための手法を提供し、これは最も効率的な実根分離アルゴリズムの基礎となっている。これは、平方無縁多項式（英語版）（つまり、重根を持たない多項式）の正の実根に関するものである。${\displaystyle a_{1},a_{2},\ldots }$ が正の実数の列である場合、

${\displaystyle c_{k}=a_{1}+{\cfrac {1}{a_{2}+{\cfrac {1}{a_{3}+{\cfrac {1}{\ddots +{\cfrac {1}{a_{k}}}}}}}}}}$

を連分数

${\displaystyle a_{1}+{\cfrac {1}{a_{2}+{\cfrac {1}{a_{3}+{\cfrac {1}{\ddots }}}}}}}$

の k 番目の収束とする。

ヴァンサンの定理 ― ${\displaystyle p_{0}(x)}$ を n 次の平方無縁多項式とし、${\displaystyle a_{1},a_{2},\ldots }$ を実数列とする。 i = 1, 2,..., に対して多項式

${\displaystyle p_{i}(x)=x^{n}p_{i-1}\left(a_{i}+1/x\right)}$

を考える。すると、 ${\displaystyle p_{k}(0)=0}$ または ${\displaystyle p_{k}}$ の係数列のいずれかが、最大で1つの符号変化を持つような整数 k が存在する。

最初のケースでは、収束する c_k は ${\displaystyle p_{0}}$ の正の根である。それ以外の場合、この符号変化の数（0 または 1）は、${\displaystyle c_{k-1}}$ と ${\displaystyle c_{k}}$ で定義される区間における ${\displaystyle p_{0}}$ の実根の数である。

ヴァンサンは定理を証明するために、それ自体でも有用な結果を証明した^[4]。

ヴァンサンの補助定理 ― p(x) が n 次の平方無縁多項式であり、a, b, c, d が ${\displaystyle \left|{\frac {a}{c}}-{\frac {b}{d}}\right|}$ のような、十分に小さい（ただし0ではない）非負の実数である場合、多項式

${\displaystyle q(x)=(cx+d)^{n}p\left({\frac {ax+b}{cx+d}}\right)}$

の係数には最大で1つの符号変化があり、この符号変化の数は ${\displaystyle {\frac {a}{c}}}$ と ${\displaystyle {\frac {b}{d}}}$ で定義される開区間における p(x) の実根の数である。

実数を扱う場合、常に c = d = 1 を選択できるが、効果的な計算は有理数で行われるため、一般的に a、b、c、d は整数であると仮定すると便利である。

「十分に小さい」条件は、Nikola Obreshkov（英語版）^[5]と Alexander Ostrowski（英語版）^[6]によって独立に定量化されている。

Obreschkoff–Ostrowski の定理：青と黄色は、0または1の符号変化を持つべき複素平面上の根が存在しない領域。左側は p の根が除外される領域、右側は変換された多項式 q の根が除外される領域。青色で示されているのは、符号変化が1回であるため除外されている領域で、符号変化が0回の場合を許容している。

定理 (Obreschkoff–Ostrowski) ― ヴァンサンの補助定理の結論は、多項式 p(x) が

${\displaystyle \left(\alpha -{\frac {a}{c}}\right)\left(\alpha -{\frac {b}{d}}\right)+\beta ^{2}\leq {\frac {1}{\sqrt {3}}}\left|\beta \left({\frac {a}{c}}-{\frac {b}{d}}\right)\right|}$

となる根 α + iβ を最大で1つ持つ場合に成立する。 特に、この結論は

${\displaystyle \left|{\frac {a}{c}}-{\frac {b}{d}}\right|<{\frac {\operatorname {sep} (p)}{2{\sqrt {3}}}}}$

の場合に成立する。 ここで sep(p) は p の2つの根間の最小距離である。

整数係数を持つ多項式の場合、最小距離 sep(p) は、多項式の次数とその係数の絶対値の最大値によって下限が定められる（多項式の根の性質の『根の分離』（英語版）を参照）。これにより、ヴァンサンの定理に基づくアルゴリズムの最悪ケースの複雑度（英語版）を解析できる。しかし、Obreschkoff–Ostrowski の定理によれば、これらのアルゴリズムの反復回数は、動作区間の近傍における根間の距離に依存するため、同じ多項式でも異なる根に対して反復回数は大きく異なる可能性がある。

James V. Uspensky（英語版）は、連分数（ヴァンサンの定理における整数 k）の長さについて、符号変化が0または1となる上限を与えた^[1][7]。

定理 (Uspensky) ―  p(x) を n 次の多項式とし、 sep(p) を p の2つの根の間の最小距離とする。

${\displaystyle \varepsilon =\left(n+{\frac {1}{n}}\right)^{\frac {1}{n-1}}-1}$

とすると、ヴァンサンの定理で存在するとされている整数 k は、

${\displaystyle F_{h-1}\operatorname {sep} (p)>2\quad {\text{and}}\quad F_{h-1}F_{k}\operatorname {sep} (p)>{\frac {1}{\varepsilon }}}$

となる最小の整数 h より大きくない。 ここで ${\displaystyle F_{h}}$ は h 番目のフィボナッチ数である。

連分数法

実根分離における連分数の使用はヴァンサンによって導入されたが、彼はジョゼフ＝ルイ・ラグランジュの成果を引用しながらも、出典を示していない^[4]。ヴァンサンの定理のアルゴリズムを作成するには、彼の定理に現れる ${\displaystyle a_{i}}$ を選択するための基準を与える必要がある。ヴァンサン自身もいくつかの選択肢を提示している（下記参照）。他にも選択肢があり得るが、アルゴリズムの効率はこれらの選択肢に大きく依存する可能性がある。以下に、これらの選択肢が後述する補助関数によって決定されるアルゴリズムを示す。

このアルゴリズムを実行するには、特定のデータ構造で表される区間のリストを使用する必要がある。このアルゴリズムは、区間を選択し、リストから削除し、0、1、または2つのより小さな区間をリストに追加し、場合によっては分離区間を出力するという動作をする。

次数 n の多項式 p(x) の実根を分離するために、各区間は ${\displaystyle (A(x),M(x))}$ のペアで表される。ここで、A(x) は次数 n の多項式であり、${\displaystyle M(x)={\frac {px+r}{qx+s}}}$ は整数係数のメビウス変換である。ひとつは

${\displaystyle A(x)=p(M(x)),}$

で、このデータ構造で表される区間は、${\displaystyle M(\infty )={\frac {p}{q}}}$ と ${\displaystyle M(0)={\frac {r}{s}}}$ を端点とする区間である。メビウス変換は、この区間における p の根を、(0, +∞) における A の根に写像する。

このアルゴリズムは、区間のリストを用いて動作する。リストには、最初に、実数を正と負の数に分割する2つの区間 ${\displaystyle (A(x)=p(x),M(x)=x)}$ と ${\displaystyle (A(x)=p(-x),M(x)=-x)}$ が含まれる（ゼロは根ではないと仮定することもできる。もしゼロが根であるならば、このアルゴリズムを p(x)/x に適用すれば良い）。次に、リスト内の各区間 (A(x), M(x)) について、アルゴリズムはそれをリストから削除する。A の係数の符号変化の数が 0 の場合、区間に根は存在せず、次の区間に進む。符号変化の数が 1 の場合、${\displaystyle M(0)}$ と ${\displaystyle M(\infty )}$ によって定義される区間は分離区間である。それ以外の場合は、正の実数 b を使用して区間 (0, +∞) を (0, b) と (b, +∞) に分割し、各部分区間について、区間を (0, +∞) にマッピングするメビウス変換を使用して M を合成し、リストに追加する2つの新しい区間を取得する。擬似コードは次のようになる。ここで、var(A) は多項式 A の係数の符号変化の数を示す。

function continued fraction is
    input: P(x) 平方無縁多項式（英語版）
    output: 分離区間を定義する有理数のペアのリスト
    /* 初期化 */
        L := [(P(x), x), (P(–x), –x)]                /* 2つの開始区間 */
        Isol := [ ]
    /* 計算 */
    while L ≠ [ ] do
        Choose (A(x), M(x)) in L, and remove it from L
        v := var(A)
        if v = 0 then exit                /* 区間内に根がない */
        if v = 1 then                     /* 分離区間が見つかった */
            add (M(0), M(∞)) to Isol
            exit
        b := some positive integer
        B(x) = A(x + b)
        w := v – var(B)
        if B(0) = 0 then                         /* 有理根が見つかった */
            add (M(b), M(b)) to Isol
            B(x) := B(x)/x
        add (B(x),  M(b + x)) to L           /* (M(b), M(+∞)) の根 */
        if w = 0 then exit                  /* ブダンの定理 */ 
        if w = 1 then                       /* 再度ブダンの定理 */ 
            add (M(0), M(b)) to Isol
        if w > 1 then
            add ( A(b/(1 + x)),  M(b/(1 + x)) )to L      /* (M(0), M(b)) の根 */

アルゴリズムの様々なバリエーションは、基本的に b の選択に依存する。ヴァンサンの論文とUspenskyの著書では、常に b = 1 としているが、Uspenskyは(0, b) に関連付けられた区間の更なる二分を避けるため、ブダンの定理を用いていない、という違いがある。

常に b = 1 を選択することの欠点は、x → 1 + x という形式の変数変換を何度も連続して行わなければならないことである。これは単一の変数変換 x → n + x で置き換えることもできるが、それでもブダンの定理を適用するためには、中間の変数変換を行う必要がある。

アルゴリズムの効率を改善する方法は、b に対して、多項式の係数から計算された正の実根の下限値を取ることである（このような下限値については多項式の根の性質（英語版）を参照）。^[8][1]

二分法

大まかに言うと、二分法は多項式のすべての実根を含む区間から出発し、それを再帰的に2つの部分に分割し、最終的に0個または1個の根を含む区間を得るというものである。開始区間は (-B, B) という形式を取ることができる。ここで、B は根の絶対値の上限であり、例えば「多項式の根の性質の『（複素）多項式の根の境界』（英語版）」で与えられる。技術的な理由（変数変換の単純化、複雑度解析の単純化、コンピュータの2元解析を活用できる可能性など）により、アルゴリズムは一般的に区間 [0, 1] から開始するものとして提示される。変数 x = By と x = –By の変換により、それぞれ区間 [0, 1] 内の正の根と負の根が移動するため、一般性が失われることはない。(変数 x = (2By – B) の変換も使用できる)

この方法では、区間に0個、1個、あるいは複数の根が存在するかどうかをテストするアルゴリズムが必要である。また、終了を保証するために、このテストアルゴリズムは「複数の根が存在する可能性」が無限回出力される可能性を排除する必要がある。スツルムの定理 とヴァンサンの補助定理は、このような便利なテストを提供する。デカルトの符号法則とヴァンサンの補助定理の使用は、スツルムの定理の使用よりもはるかに計算効率が高いため、このセクションでは前者についてのみ説明する。

デカルトの符号法則とヴァンサンの補助定理に基づく二分法は、1976年にAkritasとCollins（英語版）によって「修正Uspenskyアルゴリズム」^[3]という名称で導入され、「Uspenskyアルゴリズム」、「Vincent–Akritas–Collins アルゴリズム」、あるいは「デカルト法」とも呼ばれてきたが、デカルト、ヴァンサン、Uspenskyが実際にこの方式を記述したことはない。

この方式は以下のように機能する。ある区間における根を求めるには、まずその区間を [0, 1] に写像するため変数を変換し、新しい多項式 q(x) を生成する。[0, 1] における q の根を求めるには、変数変換 ${\displaystyle x\to {\frac {1}{x+1}}}$ によって区間 [0, 1] を [0, +∞])^{[訳注 1]} に写像して、多項式 r(x) を得る。多項式 r にデカルトの符号法則を適用すると、区間 [0, 1] における q の実根の個数が示され、したがって、[0, 1] に写像された区間における初期多項式の根の個数が示される。r の係数のシーケンスに符号変化がない場合、考慮されている区間に実根は存在しない。符号変化が1つある場合は、孤立区間となる。それ以外の場合、区間 [0, 1] を [0, 1/2] と [1/2, 1] に分割し、変数変換 x = y/2 と x = (y + 1)/2 によって、これらを [0, 1] に写像する。ヴァンサンの補助定理により、この手順は確実に終了する。

初期化を除けば、これらの変数変換はすべて、最大で2つの非常に単純な変数変換、すなわち2倍のスケーリング x → x/2、平行移動 x → x + 1、および逆数 x → 1/x の組み合わせで構成される。逆数は、多項式の係数の順序を単純に反転するだけである。計算時間のほとんどは変数の変換に費やされるため、すべての区間を [0, 1] にマッピングする方法は、優れた効率を確保するための基本となる。

擬似コード

以下の擬似コードでは、次の表記法が使用されている。

• p(x) は、区間 [0, 1] 内の実根が孤立している必要がある多項式である。
• var(q(x)) は、多項式 q の係数列における 符号変化 の数を表す。
• 作業リストの要素は (c, k, q(x)) の形式を持つ。
  • c と k は、c < 2^k を満たす 2 つの非負整数で、区間 ${\displaystyle \left[{\frac {c}{2^{k}}},{\frac {c+1}{2^{k}}}\right]}$ を表す。
  • ${\displaystyle q(x)=2^{kn}p\left({\frac {x+c}{2^{k}}}\right)}$ 、ただし n は p の次数。（多項式 q は p, c および k から直接計算できるが、アルゴリズム内で行われるよう逐次的に計算する方がコストが低くなる。p に整数係数がある場合、q についても同様）

function bisection is
    input: p(x) p(0) p(1) ≠ 0 を満たす平方無縁多項式（英語版）、
                      区間 [0, 1] 内の根が探索される
    output: 3つの数の組 (c, k, h) のリスト、
                      形式 ${\displaystyle \left[{\frac {c}{2^{k}}},{\frac {c+h}{2^{k}}}\right]}$ の分離区間を表す
    /* 初期化 */
    L := [(0, 0, p(x))] /* 作業リスト L 内の単一要素 */
    Isol := [ ]
    n := degree(p)
    /* 計算 */
    while L ≠ [ ] do
        Choose (c, k, q(x)) in L, and remove it from L
        if q(0) = 0 then
            q(x) := q(x)/x
            n := n – 1                /* 有理根が見つかった */
            add (c, k, 0) to Isol
        v := ${\displaystyle \operatorname {var} ((x+1)^{n}q(1/(x+1)))}$
        if v = 1 then                /* 分離区間が見つかった */
            add (c, k, 1) to Isol
        if v > 1 then                /* 二分する */
            add (2c, k + 1, ${\displaystyle 2^{n}q(x/2)}$  to L
            add (2c + 1, k + 1, ${\displaystyle 2^{n}q((x+1)/2)}$  to L
    end

この手順は、CollinsとAkritasによって記述されたものと本質的に同じである^[3]。実行時間は主に、考慮しなければならない区間の数と変数の変化に依存する。効率を向上させる方法はいくつかあり、アルゴリズムの発表以来、特に21世紀初頭から活発に研究されてきた。

最近の改良

Akritas-Collins二分法アルゴリズムを改良するための様々な方式が提案されている。それらには、変数変換の単純さを損なうことなく長い多項式リストの保存を回避する方式^[9]、符号変化の数に対して適切な値が得られる場合に近似演算（浮動小数点および区間演算）を使用する方式^[9]、可能な場合にニュートン法を使用する方式^[9]、高速多項式演算を使用する方式^[10]、閉根のクラスターの場合に長い二分法の連鎖を短縮する方式^[10]、多項式評価における不安定性問題を限定するための不等部分二分法^[10]などが含まれる。

これらの改良により、整数係数を持つ多項式のすべての実根を分離するアルゴリズムが生まれる。このアルゴリズムの複雑度（英語版）は

${\displaystyle {\tilde {O}}(n^{2}(k+t))}$

（対数係数を省略するためにソフトO-記法（英語版） Õ を用いた場合）である。ここで、n は多項式の次数、k は非ゼロ項の数、t は係数の桁数 の最大値である。^[10]

このアルゴリズムの実装は、非常に近い根を持つ多項式（以前は二分法で最も困難だったケース）であっても、多項式の実根を計算する他のどの実装方式よりも効率的であると考えられる^[2]。