実係数代数方程式の解（ベアストウ法）

#***********************************/
# 代数方程式の解（ベアストウ法）   */
#      例：(x+1)(x-2)(x-3)(x2+x+1) */
#           =x5-3x4-2x3+3x2+7x+6=0 */
#      coded by Y.Suganuma         */
#***********************************/

#************************************************/
# 実係数代数方程式の解（ベアストウ法）          */
#      n : 次数                                 */
#      ct : 最大繰り返し回数                    */
#      eps : 収束判定条件                       */
#      p0, q0 : x2+px+qにおけるp,qの初期値      */
#      a : 係数（最高次から与え，値は変化する） */
#      b,c : 作業域（(n+1)次の配列）            */
#      rl, im : 結果の実部と虚部                */
#      k : 結果の位置                           */
#      return : =0 : 正常                       */
#               =1 : 収束せず                   */
#      coded by Y.Suganuma                      */
#************************************************/
def Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, k)
			# 初期設定
	p1    = p0
	p2    = 0.0
	q1    = q0
	q2    = 0.0
	ind   = 0
	count = 0
#
#          １次の場合
#
	if n == 1
		if a[0].abs() < eps
			ind = 1
		else
			rl[k] = -a[1] / a[0]
			im[k] = 0.0
		end
#
#          ２次の場合
#
	elsif n == 2
					# １次式
		if a[0].abs() < eps
			if a[1].abs() < eps
				ind = 1
			else
				rl[k] = -a[2] / a[1]
				im[k] = 0.0
			end
					# ２次式
		else
			d = a[1] * a[1] - 4.0 * a[0] * a[2]
			if d < 0.0   # 虚数
				d        = Math.sqrt(-d)
				a[0]    *= 2.0
				rl[k]    = -a[1] / a[0]
				rl[k+1]  = -a[1] / a[0]
				im[k]    = d / a[0]
				im[k+1]  = -im[0]
			else           # 実数
				d       = Math.sqrt(d)
				a[0]    = 1.0 / (2.0 * a[0])
				rl[k]   = a[0] * (-a[1] + d)
				rl[k+1] = a[0] * (-a[1] - d)
				im[k]   = 0.0
				im[k+1] = 0.0
			end
		end
					# ３次以上の場合
	else
						# 因数分解
		ind = 1
		while ind > 0 && count <= ct
			for i1 in 0 ... n+1
				if i1 == 0
					b[i1] = a[i1]
				elsif i1 == 1
					b[i1] = a[i1] - p1 * b[i1-1]
				else
					b[i1] = a[i1] - p1 * b[i1-1] - q1 * b[i1-2]
				end
			end
			for i1 in 0 ... n+1
				if i1 == 0
					c[i1] = b[i1]
				elsif i1 == 1
					c[i1] = b[i1] - p1 * c[i1-1]
				else
					c[i1] = b[i1] - p1 * c[i1-1] - q1 * c[i1-2]
				end
			end
			d = c[n-2] * c[n-2] - c[n-3] * (c[n-1] - b[n-1])
			if d.abs() < eps
				return ind
			else
				dp = (b[n-1] * c[n-2] - b[n] * c[n-3]) / d
				dq = (b[n] * c[n-2] - b[n-1] * (c[n-1] - b[n-1])) / d
				p2 = p1 + dp
				q2 = q1 + dq
				if dp.abs() < eps && dq.abs() < eps
					ind = 0
				else
					count += 1
					p1 = p2
					q1 = q2
				end
			end
		end

		if ind == 0
						# ２次方程式を解く
			d = p2 * p2 - 4.0 * q2
			if d < 0.0   # 虚数
				d       = Math.sqrt(-d)
				rl[k]   = -0.5 * p2
				rl[k+1] = -0.5 * p2
				im[k]   = 0.5 * d
				im[k+1] = -im[k]
			else           # 実数
				d       = Math.sqrt(d)
				rl[k]   = 0.5 * (-p2 + d)
				rl[k+1] = 0.5 * (-p2 - d)
				im[k]   = 0.0
				im[k+1] = 0.0
			end
						# 残りの方程式を解く
			n -= 2
			for i1 in 0 ... n+1
				a[i1] = b[i1]
			end
			ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, k+2)
		end
	end

	return ind
end

				# データの設定
ct   = 1000
eps  = 1.0e-10
p0   = 0.0
q0   = 0.0
n    = 5
a    = [1.0, -3.0, -2.0, 3.0, 7.0, 6.0]
b    = Array.new(n+1)
c    = Array.new(n+1)
rl   = Array.new(n)
im   = Array.new(n)
				# 計算
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, rl, im, 0)
				# 出力
if ind > 0
	printf("収束しませんでした！\n")
else
	for i1 in 0 ... n
		printf("   %f  i %f\n", rl[i1], im[i1])
	end
end