ベアストウ法による求解

別名ヒチッコック法（というか，ヒチコック・ベアストウ法として覚えていたが）により，実係数の高次代数方程式を解く。虚数解になる場合も含み，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)