# -*- coding: UTF-8 -*-
from math import *
import numpy as np
############################################
# 実係数代数方程式の解（ベアストウ法）
# n : 次数
# ct : 最大繰り返し回数
# eps : 収束判定条件
# p0, q0 : x2+px+qにおけるp,qの初期値
# a : 係数（最高次から与え，値は変化する）
# b,c : 作業域（(n+1)次の配列）
# r : 結果
# k : 結果の位置
# return : =0 : 正常
# =1 : 収束せず
# coded by Y.Suganuma
############################################
def Bairstow(n, ct, eps, p0, q0, a, b, c, r, k) :
p1 = p0
p2 = 0.0
q1 = q0
q2 = 0.0
ind = 0
count = 0
# １次の場合
if n == 1 :
if abs(a[0]) < eps :
ind = 1
else :
r[k] = complex(-a[1] / a[0], 0)
# ２次の場合
elif n == 2 :
# １次式
if abs(a[0]) < eps :
if abs(a[1]) < eps :
ind = 1
else :
r[k] = complex(-a[2] / a[1], 0)
# ２次式
else :
D = a[1] * a[1] - 4.0 * a[0] * a[2]
if D < 0.0 : # 虚数
D = sqrt(-D)
a[0] *= 2.0
r[k] = complex(-a[1] / a[0], D / a[0])
r[k+1] = complex(-a[1] / a[0], -D / a[0])
else : # 実数
D = sqrt(D)
a[0] = 1.0 / (2.0 * a[0])
r[k] = complex(a[0] * (-a[1] + D), 0)
r[k+1] = complex(a[0] * (-a[1] - D), 0)
# ３次以上の場合
else :
# 因数分解
ind = 1
while ind > 0 and count <= ct :
for i1 in range(0, n+1) :
if i1 == 0 :
b[i1] = a[i1]
elif i1 == 1 :
b[i1] = a[i1] - p1 * b[i1-1]
else :
b[i1] = a[i1] - p1 * b[i1-1] - q1 * b[i1-2]
for i1 in range(0, n+1) :
if i1 == 0 :
c[i1] = b[i1]
elif i1 == 1 :
c[i1] = b[i1] - p1 * c[i1-1]
else :
c[i1] = b[i1] - p1 * c[i1-1] - q1 * c[i1-2]
D = c[n-2] * c[n-2] - c[n-3] * (c[n-1] - b[n-1])
if fabs(D) < 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 abs(dp) < eps and fabs(dq) < eps :
ind = 0
else :
count += 1
p1 = p2
q1 = q2
if ind == 0 :
# ２次方程式を解く
D = p2 * p2 - 4.0 * q2
if D < 0.0 : # 虚数
D = sqrt(-D)
r[k] = complex(-0.5 * p2, 0.5 * D)
r[k+1] = complex(-0.5 * p2, -0.5 * D)
else : # 実数
D = sqrt(D)
r[k] = complex(0.5 * (-p2 + D), 0)
r[k+1] = complex(0.5 * (-p2 - D), 0)
# 残りの方程式を解く
n -= 2
for i1 in range(0, n+1) :
a[i1] = b[i1]
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, r, k+2)
return ind
----------------------------------
# -*- coding: UTF-8 -*-
import numpy as np
from math import *
from function import Bairstow
############################################
# 代数方程式の解（ベアストウ法）
# 例：(x+1)(x-2)(x-3)(x2+x+1)
# =x5-3x4-2x3+3x2+7x+6=0
# coded by Y.Suganuma
############################################
# データの設定
ct = 1000
eps = 1.0e-10
p0 = 0.0
q0 = 0.0
n = 5
a = np.array([1.0, -3.0, -2.0, 3.0, 7.0, 6.0])
b = np.empty(n+1, np.float)
c = np.empty(n+1, np.float)
r = np.empty(n, np.complex)
# 計算
ind = Bairstow(n, ct, eps, p0, q0, a, b, c, r, 0)
# 出力
if ind > 0 :
print("収束しませんでした！")
else :
for i1 in range(0, n) :
print(" " + str(r[i1]))