147

第 13 章 代数方程式の解法
最近は代数学の講義がすっかりモダンになってしまって，群とか束とかイデヤルだとか同型定理だとか代数的整数論だとか，そういうものばかり教えているようであるが，昔は方程式論というものがあって，数値計算に密着した定理をいろいろ教えてくれたものである。
戸川隼人「マトリクスの数値計算」(オーム社)

代数方程式は前章で取り扱った非線型方程式の部分集合をなすものである。前述のように，非線型方程式は多種多様でありそれを統一的に取り扱うことは難しいが，代数方程式，特に一次元のものについては，解が存在することが証明されている上，多項式に対する演算は取り扱いが容易であるため，それらの性質を生かした多様な解法が提案されている。本章ではその一部を表面的に取り上げてみたい。

13.1 代数方程式の特徴
本章で取り扱う代数方程式をn∑ai xi = 0 (an , 0, ai ∈ R)

(13.1)

i=0とする。場合に応じて左辺の多項式を pn (x) =∑nii=0 ai x と置き換えて

pn (x) = 0(13.2)とも書くことにする。
行列の固有値・固有ベクトル計算で述べたように，5 次未満の 1 次元代数方程式には解の公式が存在する。即ち，有限回の代数的演算 (四則演算と有理数のべき乗) で真の解を得ることが可能である。また，5 次以上の代数方程式については一般の非線型方程式と同様，Newton 法などの無限回反復法を適用しなければならない。
一方，一般の非線型方程式と決定的に違うのは，解の存在が保証されているということである。
定理 13.1.1 (代数学の基本定理)

n 次代数方程式 (13.1) は複素数体 C 内に，重複も込めて n 個の解 α1 , α2 , ..., αn が必ず存在する。
(13.2) の解が αi (i = 1, 2, ..., n) であれば，左辺の多項式の因数分解
pn (x) = ann∏i=1(x − αi )(13.3)

第 13 章 代数方程式の解法

148を求めたことと同じ意味を持つ。

また，一般の非線型方程式と異なり， pn (x) の任意階導関数は簡単に求めることが出来，Newton
法の反復で必要となる，k 回反復時の近似値 xk における関数値 pn (xk ) と微分係数 p′n (xk ) の値は次のようにして同時に求めることが可能である。
アルゴリズム 30 (多項式の値と微分係数の計算)

1. bn := an
2. cn−1 := bn
3. i = n − 1, ..., 1(a) bi := ai + bi+1 xk(b) ci−1 := bi + ci xk
4. b0 := a0 + b1 · xk以上の計算により，
b0 = pn (xk )(13.4)

c0 = p′n (xk )(13.5)

を得る。
問題 13.1.1
アルゴリズム 30 によって pn (xk ) と p′n (xk ) が求められることを示せ。

13.2 4 次以下の代数方程式の解公式
1 次方程式
a1 x + a0 = 0 (a1 , 0)(13.6)はx = −a0 /a1とすればよい。2 次方程式

a2 x2 + a1 x + a0 = 0
は第 3 章で述べた通り

x=−a1 ±√

a21 − 4a2 a02a2を桁落ちを避けるように計算すればよい。
以下 3 次, 4 次代数方程式の解公式について述べる。以下の記述は「数値計算ハンドブック」(オーム社) による。

13.2. 4 次以下の代数方程式の解公式

13.2.1

149

Cardano 法 — 3 次代数方程式の解公式

3 次方程式
a3 x3 + a2 x2 + a1 x + a0 = 0 (a3 , 0)(13.7)

に対する解公式は Cardano 法と呼ばれる。この解公式を導く手順を大雑把に示す。
まず

x=y−a23a3とし，式 (13.7) に代入して

y3 + 3py + q = 0(13.8)と式変形する。この時

p =q =−a229a232a3227a33+a13a3−a1 a2 a0+a33a23である。更に y = u + v とし，かつ u と v は uv = −p を満足するものとすると，式 (13.8) の y に代入してu3 + v3 = −qを得る。また u3 v3 = −p3 なので，u3 , v3 は 2 次方程式

z2 + qz − p3 = 0の解である。従って
3u=v3=√q2 + 4p32√
−q − q2 + 4p32−q +である。これを満足するもののうち一つを û, v̂ とし，1 の複素数 3 乗根のうちどちらか一つを ω とすれば，uv = −p を満足するものは

β1=û + v̂β2=

ωû + ω2 v̂

β3=

ω2 û + ωv̂

である。従って式 (13.8) の解 β1 , β2 , β3 を具体的に書くと

β1=β2=β3=û + v̂

√√31(û − v̂)− (û + v̂) − −122√√
13− (û + v̂) + −1(û − v̂)22(13.9)

第 13 章 代数方程式の解法

150となる。ここで√−1 は虚数単位である。よって，これらを用いて元の 3 次代数方程式 (13.7) の解αi (i = 1, 2, 3) は
αi = βi −a2(i = 1, 2, 3)3a3となる。
問題 13.2.1

1. 1 の n 乗根を ωi (i = 1, 2, ..., n) とする時，
ωi = cos√
2π2π(i − 1) + −1 sin (i − 1)nnであることを確認せよ。

2. 解の公式 (13.9) を計算する際に留意する点を考えよ (Hint: 2 次方程式の解を求める時の留意事項を当てはめればよい)。

13.2.2 Ferrari 法 — 4 次代数方程式の解公式
4 次方程式
a4 x4 + a3 x3 + a2 x2 + a1 x + a0 = 0 (a4 , 0)(13.10)に対する解公式は Ferrari 法と呼ばれる。これも Cardano 法と同様に

x=y−a34a4

とおいて式 (13.10) に代入することによりy4 + py2 + qy + r = 0となる。ここでp =q =r=−3a238a24a33
8a34−+−

3a43a2a4a2 a3 a1+a42a24+
4256a4a2 a23
16a34−a1 a3 a0+a44a24である。

q = 0 の時は直ちに因数分解でき√√ 2 −p − p2 − 4r   2 −p + p2 − 4r  y − = 0y −22を y について解けばよい。

(13.11)

13.3. 減次を用いた Newton 法

151

q , 0 の時は，式 (13.11) をy4 = −py2 − qy − rとし，この両辺に y2 z + z2 /4 を加えて(()2q1z )2y += (z − p) y −+(z3 − pz2 − 4rz + 4pr − q2 )22(z − p)4(z − p)2(13.12)と式変形する。さすれば右辺の 3 次式が z3 − pz2 − 4rz + 4pr − q2 = 0 となる z を一つ見つければ，

q = 0 の時と同様，直ちに因数分解でき((()) ())
z √qqz √22y + − z− p y−y + + z− p y−=022(z − p)22(z − p)(13.13)を y について解けばよいことになる。

13.3 減次を用いた Newton 法
前述のように，代数方程式は非線型方程式の一種であるから，単純に Newton 法を適用すれば近似根が得られる。しかし，複数個の根を求めるには，一つの根 α が求められる毎にpn (x)=0x−α(13.14)

として，次数が一つ減った代数方程式に対して Newton 法を適用する必要がある。さもなければ同じ解を何度も求めることになりかねない。 このように代数方程式の次数を下げることを減次(deflation) と呼ぶ。
一般に減次は次のようにして行うことが推奨されている。
1. 絶対値最小の根から求める際には，高次項の係数から計算する。
2. 絶対値最大の根から求める際には，低次項の係数から計算する。
減次した代数方程式 (13.14) の左辺の係数を高次項から求めるには，アルゴリズム 30 を使えばよい。
アルゴリズム 31 (高次項の係数から求める減次)

1. bn−1 := an
2. i = n − 2, ..., 0(a) bi := ai+1 + bi+1 α
3. pn (x)/(x − α) =∑n−1i=0 bi xi

また，低次項の係数から求める際には次のようにする。
アルゴリズム 32 (低次項の係数から求める減次)

1. b0 := −a0 /α

第 13 章 代数方程式の解法

152

表 13.1: 高次項から係数を求めた場合○絶対値最小の根から×絶対値最大の根から

11.00000000000000089e+00

9.99999999998291855e+00

21.99999999999990452e+00

9.00000000016882851e+00

33.00000000000209166e+00

7.99999999927011718e+00

43.99999999997718714e+00

7.00000000182702831e+00

55.00000000011708146e+00

5.99999999709054155e+00

65.99999999968412556e+00

5.00000000305080494e+00

77.00000000047785509e+00

3.99999999789091065e+00

87.99999999959239183e+00

3.00000000092989616e+00

99.00000000018273916e+00

1.99999999976192555e+00

109.99999999996662403e+00

1.00000000002703437e+00

2. i = 1, 2, ..., n − 1(a) bi := −(ai − bi−1 )/α
3. pn (x)/(x − α) =∑n−1

i=0 bi xi

例題 13.3.1

10 次の実係数代数方程式(x − 10)(x − 9) · · · (x − 1) = 0(13.15)

に対して，減次を伴った Newton 法を適用した結果を表 13.1 と表 13.2 に示す。
減次の原則が示す通り，絶対値最小の根から求める際には高次項の係数から，絶対値最大の根から求める際には低次項の係数から減次した方が近似根の精度に与える影響は少なくなることが分かる。
問題 13.3.1
次の 3 次方程式の解をすべて求めよ。

x3 − 7.2741x2 + 16.8267x − 12.077 = 0

13.4 Durand-Kerner-Aberth 法
ここでは pn (x) を

qn (x) = xn + cn−1 xn−1 + · · · + c1 x + c0と変形しqn (x) = 0(13.16)

13.4. Durand-Kerner-Aberth 法

153

表 13.2: 低次項から係数を求めた場合
×絶対値最小の根から○絶対値最大の根から

11.00000000000000089e+00

9.99999999998291855e+00

21.99999999999668643e+00

9.00000000007256773e+00

33.00000000076505957e+00

7.99999999986114130e+00

43.99999996825865445e+00

7.00000000015056667e+00

55.00000044368798147e+00

5.99999999990402166e+00

65.99999725179189802e+00

5.00000000003449774e+00

77.00000856141343686e+00

3.99999999999397948e+00

87.99998605011867259e+00

3.00000000000018829e+00

99.00001132664817405e+00

2.00000000000007105e+00

109.99999639010785124e+00

9.99999999999995337e-01という代数方程式について考えることにする。ここでci =ai(i = 0, 1, ..., n − 1)anである。
代数方程式 (13.16) に対して n 個の解 α1 , α2 , ..., αn と係数の関係式を書き下してみると
∑(−1)1 ni=1 αi − cn−1 = 0∑(−1)2 ni1 <i2 αi1 αi2 − cn−2 = 0
∑


(−1)3 ni1 <i2 <i3 αi1 αi2 αi3 − cn−3 = 0
(13.17)
..


.


∑

(−1)n−1 ni1 <i2 <···<in−1 αi1 αi2 · · · αin−1 − c1 = 0(−1)n α α · · · α − c = 0
1 2n0となる。この n 次元非線型方程式 (13.17) の解は代数方程式 (13.16) の解でもある。これに Newton法を適用する方法を DKA 法 (Durand-Kerner-Aberth 法)1 と呼ぶ。
アルゴリズム 33 (DKA 法)

1. 初期値 z(0)i (i = 1, 2, ..., n) を設定する。
2. k = 0, 1, 2, ... に対して次の反復を繰り返す。
z(k+1):= z(k)ii − ∏nqn (z(k)
i )(k)(k)
j=1, j,i (zi − z j )(13.18)

この際，初期値としては Aberth が提案したz(0)i := −{)}(√3cn−12(i − 1)π+ r exp+−1nn2n(13.19)1 一松 [9] によれば，Durand が代数方程式用に Newton 法を改良し，Kerner が対称式と係数の関係式からの解釈を与え，Aberth が初期値を提案した，という流れになるらしい。

第 13 章 代数方程式の解法

154が良いとされている。ここで r は全ての根が複素平面上の円盤z+cn−1≤rnに存在するように取る。例えば伊理 [13] は実係数代数方程式

xn − |cn−2 |xn−2 − |cn−3 |xn−3 − · · · − |c1 |x − |c0 | = 0(13.20)の正の実数根を r とすることを提案している。
例題 13.4.1

p5 (x) = x5 − 15x4 + 85x3 − 225x2 + 274x − 120とすると，これはp5 (x) = (x − 1)(x − 2)(x − 3)(x − 4)(x − 5)と因数分解できる。これを IEEE754 倍精度の DKA 法で計算すると図 13.1 のような収束状況になる。これは初期値の r をかなり大きめに取ってある。
5虚数軸

虚数軸45032100-1-2-3-50-4-500

実軸50-5-2-10123実軸45678

図 13.1: DKA 法の収束 (左図の中心部を拡大したのが右図)

問題 13.4.1

1. DKA 法の反復式が非線型方程式 (13.17) に対する Newton 法であることを，n = 2 の場合について確認せよ。

2. 実係数代数方程式 (13.20) が正の実数解を持つことを示せ。

13.5. 代数方程式の数値解の誤差解析

155

表 13.3: (x − 10)(x − 9) · · · (x − 1) = 0 の数値解と相対誤差の評価値
i
数値解
相対誤差の評価値

19.99999999999999112e-01

2.2 × 10−15

22.00000000000017719e+00

2.0 × 10−14

33.00000000000170752e+00

1.2 × 10−1343.99999999999103784e+00

5.0 × 10−13

55.00000000000905143e+00

1.5 × 10−12

66.00000000014250023e+00

3.0 × 10−1276.99999999954060925e+00

4.4 × 10−128

7.99999999972644638e+00

4.2 × 10−1298.99999999990030908e+00

2.4 × 10−12109.99999999997828226e+00

6.1 × 10−13

13.5 代数方程式の数値解の誤差解析
代数方程式の解は近接しているほど悪条件になる。特に m 重根の場合， pn (x) = (x − α)m qn−m (x)であるから，近似解 z が|z − α| ≤ ε1/mMの範囲に入ってしまえば|pn (z)| ≤ ε M |qn−m (z)|となる。|qn−m (z)| ≈ 1 程度であれば，この近似値 z の精度をこれ以上上げることは難しい。一般に「m 重解は 1/m 桁程度しか求めることが出来ない」と言われている理由はここにある。
しかし，重解でなくとも，それに近い状態，即ち解が近接していればそれも悪条件と言える。それを表す指標として，次の評価式が提案されている [26]。

ei の相対誤差
代数方程式の (13.3) 式のように表せる時，解は α1 , α2 , ..., αn であるが，その近似解 αrE(eαi ) はn∏max(|αi |, |α j |)rE(eαi ) ≈εM(13.21)
αi − α jj=1, j,iで評価できるというものである。Smith の評価式と呼ばれているものの簡易形と解釈出来る。
次の数値例でこの評価式の持つ意味を考えてみたい。
問題 13.5.1 (Wilkinson の数値例)

Winkinson が最初に取り上げた例として(x − 10)(x − 9) · · · (x − 1) = 0(13.22)というものがある。αi = i (i = 1, 2, .., 10) として，これを上記評価式に当て嵌めて，数値解を比較みると表 13.3 のようになる。以下全て IEEE754 倍精度で DKA 法を適用した結果である。根が「近第 13 章 代数方程式の解法156接」してかどうかは根の絶対値によって相対的に決まるものである。このまま次数を高くすると段々とその距離が「相対的に近く」なってくる。
更に次数を高くした数値例を示す。同様にして 20 次まで次数を増やし(x − 20)(x − 19) · · · (x − 1) = 0(13.23)とすると，多項式の係数すら表 13.4 の左に示すような値となり，IEEE754 倍精度では収まりきれずに丸め誤差が混入してしまう。これに DKA 法を適用してみると表 13.4 の右に示すような誤差の多い値となり，特に 12〜16 までの数値解の相対誤差が増大していることが分かる。

演習問題
1. monic な代数方程式 (13.16) は，次のコンパニオン行列 C ∈ Mn (R) 0
 0 .
C =  ..
 0−c010..
.

01..
.

00..
.

······..
.

···

00100..
.
0···−cn−4−cn−3−cn−2−cn−1
の固有値多項式と同一視できることを確認せよ。また，この性質を利用して，代数方程式を行列 C の固有値問題として解き，DKA 法を用いて解いた結果と比較せよ。

13.5. 代数方程式の数値解の誤差解析

157

表 13.4: (x − 20)(x − 19) · · · (x − 1) = 0 の係数 (左) と数値解 (右)a0=2432902008176640000

a1=−8752948036761600000

=i

数値解

a213803759753640704000

11.00000000000000178e+00

a3=−12870931245150988800

22.00000000000251710e+00

a4=8037811822645051776

32.99999999997963007e+00

44.00000000275738898e+00

54.99999989747896922e+00

66.00000035498705131e+00

76.99999671931183531e+00

a5=−3599979517947607200

a6=

1206647803780373360

a7=−311333643161390640

a8=63030812099294896

88.00004935667591077e+00

a9= −10142299865511450

98.99982516210797456e+00

a10= 1307535010540395

10

9.99978049671140212e+00

a11= −135585182899530

111.10009569379923811e+01

a12= 11310276995381

12

1.20021141240172842e+01

a13= −756111184500

131.29979541385837347e+01= 40171771630

14

1.39970718883922469e+01

a14

15

1.49928298862719380e+01

a15= −1672280820

161.59948146111038749e+01

a16= 53327946

171.70011031622414208e+01

a17= −1256850

18

1.79999214964136165e+01

a18= 2061519

1.89998301783606749e+01

a19= −21020

2.00000269515716731e+01

a20= 1