別名ヒチッコック法（というか，ヒチコック・ベアストウ法として覚えていたが）により，実係数の高次代数方程式を解く。虚数解になる場合も含み，n 次式の場合は n 個の解を全て求める。
インターネット上にもプログラム例が見られるが，数値計算的，プログラミング的に見た場合，ちょっとひどいなあというのも（箇所も）見られるので，お手本とはいわないまでも，まずまずのプログラムを示しておこう。
まずは R によるプログラム例，その後に Python によるプログラム例を掲示する。
quadratic.root = function(p, q) {
d = p ^ 2 - 4 * q
if (d >= 0) {
x1 = (-p - ifelse(p >= 0, 1,-1) * sqrt(d)) / 2 # 絶対値の大きいほうの解を求める
x2 = q / x1 # もう一方は解と係数の関係から得る
} else {
x1 = (-p + sqrt(d + 0i)) / 2 # 大抵のプログラミング言語は虚数を取り扱える
x2 = Conj(x1) # 共役複素数を求める関数もある
}
return(c(x1, x2))
}
Bairstow = function(a, EPSILON = 1e-14) {
a = a / a[1]
ans = NULL
repeat {
n = length(a)
if (n == 3) { # n = 3, n = 2 の場合の処理は，ループの最初に書かねばならない
ans = c(ans, quadratic.root(a[2], a[3]))
break
} else if (n == 2) {
ans = c(ans, -a[2] / a[1])
break
}
p = q = 1
b = c = numeric(n)
c[1] = b[1] = 1 # a[1], b[1], c[1] は常に 1 である
repeat {
b[2] = a[2] - p * b[1]
c[2] = b[2] - p * c[1]
for (i in 3:n) {
b[i] = a[i] - p * b[i - 1] - q * b[i - 2]
c[i] = b[i] - p * c[i - 1] - q * c[i - 2]
}
d = c[n - 2] ^ 2 - c[n - 3] * (c[n - 1] - b[n - 1])
delta.p = (b[n - 1] * c[n - 2] - b[n] * c[n - 3]) / d
delta.q = (b[n] * c[n - 2] - b[n - 1] * (c[n - 1] - b[n - 1])) / d
if (abs(delta.p) < EPSILON && abs(delta.q) < EPSILON) {
break
}
p = p + delta.p
q = q + delta.q
}
ans = c(ans, quadratic.root(p, q)) # x^2 + px + q = 0 を解く
a = b[1:(n - 2)] # 元の多項式を x^2 + px + q = 0 で割った多項式の係数
}
return(sort(ans))
}
> options(digits = 14)
> Bairstow(c(1.0, -15.0, 85.0, -225.0, 274.0, -120.0))
[1] 1 2 3 4 5
> Bairstow(c(1, 10, 35, 50, 24))
[1] -4 -3 -2 -1
> Bairstow(c(1, 10, 25, 50, 24))
[1] -7.49826796187678+0.0000000000000i
[2] -0.93451222322734-2.0458454872479i
[3] -0.93451222322734+2.0458454872479i
[4] -0.63270759166854+0.0000000000000i
> Bairstow(c(1, 1, 1, 1))
[1] -1+0i 0-1i 0+1i
> Bairstow(c(1, 1, 1))
[1] -0.5-0.86602540378444i -0.5+0.86602540378444i
> Bairstow(c(3, 5))
[1] -1.6666666666667
> Bairstow(c(1,-2, 1))
[1] 1 1
> quadratic.root(-1.000000001, 0.000000001)
[1] 1e+00 1e-09
> quadratic.root(0, 1)
[1] 0+1i 0-1i
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
Python による場合
import scipy as sp
def quadratic_root(p, q):
d = p**2 - 4 * q
if d >= 0:
x1 = (-p - sp.copysign(sp.sqrt(d), p)) / 2
x2 = q / x1
else:
x1 = (-p + sp.sqrt(d + 0j)) / 2
x2 = x1.conjugate()
return x1, x2
def Bairstow(a, EPSILON = 1e-14):
a = sp.ravel(a)
a = a / a[0]
ans = []
while True:
n = len(a) - 1
if n == 2:
ans.extend(quadratic_root(a[1], a[2]))
break
elif n == 1:
ans.append(-a[1] / a[0])
break
p = q = 1
b = sp.ones_like(a)
c = sp.ones_like(a)
while True:
b[1] = a[1] - p * b[0]
c[1] = b[1] - p * c[0]
for i in range(2, n + 1):
b[i] = a[i] - p * b[i - 1] - q * b[i - 2]
c[i] = b[i] - p * c[i - 1] - q * c[i - 2]
d = c[n - 2]**2 - c[n - 3] * (c[n - 1] - b[n - 1])
delta_p = (b[n - 1] * c[n - 2] - b[n] * c[n - 3]) / d
delta_q = (b[n] * c[n - 2] - b[n - 1] * (c[n - 1] - b[n - 1])) / d
if abs(delta_p) < EPSILON and abs(delta_q) < EPSILON:
break
p += delta_p
q += delta_q
ans.extend(quadratic_root(p, q))
a = b[:n - 1]
return sp.sort(ans)
>>> Bairstow(sp.array([1, -15, 85, -225, 274, -120.0]))
array([1., 2., 3., 4., 5.])
>>> Bairstow([1, 10, 35, 50, 24])
array([-4., -3., -2., -1.])
>>> Bairstow([1, 10, 25, 50, 24])
array([-7.49826796+0.j , -0.93451222-2.04584549j,
-0.93451222+2.04584549j, -0.63270759+0.j ])
>>> Bairstow([1, 1, 1, 1])
array([-1.+0.j, 0.-1.j, 0.+1.j])
>>> Bairstow([1, 1, 1])
array([-0.5-0.8660254j, -0.5+0.8660254j])
>>> Bairstow([3, 5])
array([-1.66666667])
>>> Bairstow([1, -2, 1])
array([1., 1.])
>>> quadratic_root(-1.000000001, 0.000000001)
(1.0, 1e-09)
>>> quadratic_root(0, 1)
(1j, -1j)
