$ n $次多項式
f(z) = c_0 z^n + c_1 z^{n-1} + \cdots + c_{n-1} z + c_n
の係数$ c_0, c_1, \dots, c_n$を引数として与えた場合に、その$ n $個の根――すなわち
f(z) = c_0 \prod_{k=1}^n (z-\alpha_k) = 0
を満たす$ \alpha_0, \alpha_1, \dots, \alpha_n$を得るVBAプログラムを次のとおり自作した。
Classモジュール
複素数を表すClassモジュール
complex.cls
Option Explicit
Private Const strClassName As String = "Complex"
'============================
'プロパティ
'============================
Public re As Double '実部
Public im As Double '虚部
'=============================
'コンストラクタ（初期化処理）
'=============================
Private Sub Class_Initialize()
re = 0
im = 0
End Sub
'=============================
'複素数の定義(実部と虚部) a+bi
'=============================
Public Function SetC(a As Double, b As Double)
re = a
im = b
End Function
'=============================
'複素数の定義(絶対値と偏角) re^{iΘ}
'=============================
Public Function SetC_RT(r As Double, theta As Double)
re = r * Cos(theta)
im = r * Sin(theta)
End Function
'=============================
'複素数の定義(実部と虚部の文字列)
'=============================
Public Function SetC_str(strComplex As String)
re = WorksheetFunction.ImReal(strComplex)
im = WorksheetFunction.Imaginary(strComplex)
End Function
'=============================
'絶対値 |a| = a.AbsC()
'=============================
Public Function AbsC() As Double
AbsC = Sqr(re ^ 2 + im ^ 2)
End Function
'=============================
'絶対値 |a| = a.AbsCSq()
'=============================
Public Function AbsCSq() As Double
AbsCSq = re ^ 2 + im ^ 2
End Function
'=============================
'偏角 ∠a = a.Angle()
'=============================
Public Function Angle() As Double
If re ^ 2 + im ^ 2 = 0 Then
Angle = 0
Else
Angle = WorksheetFunction.Atan2(re, im)
End If
End Function
'=============================
'足し算 a+b = a.Plus(b)
'=============================
Public Function Plus(ByVal b As Variant) As Variant
If (2 <= VarType(b) And VarType(b) <= 6) Then
Set Plus = New Complex
Plus.re = re + b
Plus.im = im
ElseIf (VarType(b) = 9 And TypeName(b) = strClassName) Then
Set Plus = New Complex
Plus.re = re + b.re
Plus.im = im + b.im
Else
Plus = CVErr(xlErrNum)
End If
End Function
'=============================
'引き算 a-b = a.Minus(b)
'=============================
Public Function Minus(ByVal b As Variant) As Variant
If (2 <= VarType(b) And VarType(b) <= 6) Then
Set Minus = New Complex
Minus.re = re - b
Minus.im = im
ElseIf (VarType(b) = 9 And TypeName(b) = strClassName) Then
Set Minus = New Complex
Minus.re = re - b.re
Minus.im = im - b.im
Else
Minus = CVErr(xlErrNum)
End If
End Function
'=============================
'掛け算 a*b = a.Times(b)
'=============================
Public Function Times(ByVal b As Variant) As Variant
If (2 <= VarType(b) And VarType(b) <= 6) Then
Set Times = New Complex
Times.re = re * b
Times.im = im * b
ElseIf (VarType(b) = 9 And TypeName(b) = strClassName) Then
Set Times = New Complex
Times.re = re * b.re - im * b.im
Times.im = re * b.im + im * b.re
Else
Times = CVErr(xlErrNum)
End If
End Function
'============================
'割り算 a/b = a.Div(b)
'============================
Public Function Div(ByVal b As Variant) As Variant
If (2 <= VarType(b) And VarType(b) <= 6) Then
If b = 0 Then
Div = CVErr(xlErrNum)
Else
Set Div = New Complex
Div.re = re / b
Div.im = im / b
End If
ElseIf (VarType(b) = 9 And TypeName(b) = strClassName) Then
If b.AbsC() = 0 Then
Div = CVErr(xlErrNum)
Else
Set Div = New Complex
Div.re = (re * b.re + im * b.im) / (b.AbsC() ^ 2)
Div.im = (-re * b.im + im * b.re) / (b.AbsC() ^ 2)
End If
Else
Div = CVErr(xlErrNum)
End If
End Function
'============================
'共役複素数
'============================
Public Function CompCnj() As Complex
Set CompCnj = New Complex
CompCnj.re = re
CompCnj.im = (-1) * im
End Function
'============================
'文字列(直交座標)
'============================
Public Function CompStr() As String
If re = 0 Then
If im = 0 Then
CompStr = "0"
ElseIf im > 0 Then
If im = 1 Then
CompStr = "i"
Else
CompStr = Str(im) & "i"
End If
Else
If im = -1 Then
CompStr = "-i"
Else
CompStr = "-" & Str(Abs(im)) & "i"
End If
End If
ElseIf re > 0 Then
If im = 0 Then
CompStr = Str(re)
ElseIf im > 0 Then
If im = 1 Then
CompStr = Str(re) & "+i"
Else
CompStr = Str(re) & "+" & Str(im) & "i"
End If
Else
If im = -1 Then
CompStr = Str(re) & "-i"
Else
CompStr = Str(re) & Str(im) & "i"
End If
End If
Else
If im = 0 Then
CompStr = "-" & Str(Abs(re))
ElseIf im > 0 Then
If im = 1 Then
CompStr = "-" & Str(Abs(re)) & "+i"
Else
CompStr = "-" & Str(Abs(re)) & "+" & Str(im) & "i"
End If
Else
If im = -1 Then
CompStr = "-" & Str(Abs(re)) & "-i"
Else
CompStr = "-" & Str(Abs(re)) & "-" & Str(Abs(im)) & "i"
End If
End If
End If
CompStr = Replace(CompStr, " ", "")
End Function
'============================
'文字列(直交座標,極座標cossin,極座標rExpΘ)
'============================
Public Function CompStrPolar() As String
Dim dblAbs As Double
Dim dblAngle As Double
dblAbs = Sqr(re ^ 2 + im ^ 2)
If re ^ 2 + im ^ 2 = 0 Then
dblAngle = 0
Else
dblAngle = WorksheetFunction.Atan2(re, im)
End If
CompStrPolar = Str(dblAbs) & " Exp(" & Str(dblAngle) & "i)"
End Function
標準モジュール
各種計算用の函数を定義する標準モジュール
ComplexFunc.bas
Option Explicit
Private Const strTypeComplex As String = "Complex"
Private Const cnstDblEpsilon As Double = 10 ^ (-12) 'ニュートン法の収束基準
Private Const PI As Double = 3.14159265358979 '円周率
'=============================
'Sum函数
'=============================
Public Function SumComp(ParamArray varArray()) As Variant
Dim cnt As Long
Dim cmpSum As Complex: Set cmpSum = New Complex
For cnt = LBound(varArray) To UBound(varArray)
If (2 <= VarType(varArray(cnt)) And VarType(varArray(cnt)) <= 6) Then
cmpSum.re = cmpSum.re + varArray(cnt)
cmpSum.im = cmpSum.im + 0
ElseIf (VarType(varArray(cnt)) = 9 And TypeName(varArray(cnt)) = strTypeComplex) Then
cmpSum.re = cmpSum.re + varArray(cnt).re
cmpSum.im = cmpSum.im + varArray(cnt).im
Else
SumComp = CVErr(xlErrNum)
Exit Function
End If
Next cnt
Set SumComp = cmpSum
End Function
'=============================
'Sub函数（第1引数を被減数として、第2引数以降で順次減じていく）
'=============================
Public Function SubComp(ParamArray varArray()) As Variant
Dim cnt As Long
Dim dblTmp As Double
Dim cmpSub As Complex: Set cmpSub = New Complex
Dim dblSq As Double
For cnt = LBound(varArray) To UBound(varArray)
If (2 <= VarType(varArray(cnt)) And VarType(varArray(cnt)) <= 6) Then
If cnt = LBound(varArray) Then
cmpSub.re = varArray(cnt)
cmpSub.im = 0
Else
cmpSub.re = cmpSub.re - varArray(cnt)
End If
ElseIf (VarType(varArray(cnt)) = 9 And TypeName(varArray(cnt)) = strTypeComplex) Then
If cnt = LBound(varArray) Then
cmpSub.re = varArray(cnt).re
cmpSub.im = varArray(cnt).im
Else
cmpSub.re = cmpSub.re - varArray(cnt).re
cmpSub.im = cmpSub.im - varArray(cnt).im
End If
Else
SubComp = CVErr(xlErrNum)
Exit Function
End If
Next cnt
Set SubComp = cmpSub
End Function
'=============================
'Product函数
'=============================
Public Function ProductComp(ParamArray varArray()) As Variant
Dim cnt As Long
Dim dblTmp As Double
Dim cmpProduct As Complex: Set cmpProduct = New Complex
For cnt = LBound(varArray) To UBound(varArray)
If (2 <= VarType(varArray(cnt)) And VarType(varArray(cnt)) <= 6) Then
If cnt = LBound(varArray) Then
cmpProduct.re = varArray(cnt)
cmpProduct.im = 0
Else
dblTmp = cmpProduct.re * varArray(cnt) - cmpProduct.im * 0
cmpProduct.im = cmpProduct.re * 0 + cmpProduct.im * varArray(cnt)
cmpProduct.re = dblTmp
End If
ElseIf (VarType(varArray(cnt)) = 9 And TypeName(varArray(cnt)) = strTypeComplex) Then
If cnt = LBound(varArray) Then
cmpProduct.re = varArray(cnt).re
cmpProduct.im = varArray(cnt).im
Else
dblTmp = cmpProduct.re * varArray(cnt).re - cmpProduct.im * varArray(cnt).im
cmpProduct.im = cmpProduct.re * varArray(cnt).im + cmpProduct.im * varArray(cnt).re
cmpProduct.re = dblTmp
End If
Else
ProductComp = CVErr(xlErrNum)
Exit Function
End If
Next cnt
Set ProductComp = cmpProduct
End Function
'=============================
'Div函数（第1引数を被除数として、第2引数以降で順次除していく）
'=============================
Public Function DivComp(ParamArray varArray()) As Variant
Dim cnt As Long
Dim dblTmp As Double
Dim cmpDiv As Complex: Set cmpDiv = New Complex
Dim dblSq As Double
For cnt = LBound(varArray) To UBound(varArray)
If (2 <= VarType(varArray(cnt)) And VarType(varArray(cnt)) <= 6) Then
If cnt = LBound(varArray) Then
cmpDiv.re = varArray(cnt)
cmpDiv.im = 0
Else
If varArray(cnt) = 0 Then
DivComp = CVErr(xlErrNum)
Exit Function
Else
dblTmp = cmpDiv.re / varArray(cnt)
cmpDiv.im = cmpDiv.im / varArray(cnt)
cmpDiv.re = dblTmp
End If
End If
ElseIf (VarType(varArray(cnt)) = 9 And TypeName(varArray(cnt)) = strTypeComplex) Then
If cnt = LBound(varArray) Then
cmpDiv.re = varArray(cnt).re
cmpDiv.im = varArray(cnt).im
Else
dblSq = varArray(cnt).re ^ 2 + varArray(cnt).im ^ 2
If dblSq = 0 Then
DivComp = CVErr(xlErrNum)
Exit Function
Else
dblTmp = cmpDiv.re * varArray(cnt).re + cmpDiv.im * varArray(cnt).im
cmpDiv.im = cmpDiv.im * varArray(cnt).re - cmpDiv.re * varArray(cnt).im
cmpDiv.re = dblTmp / dblSq
cmpDiv.im = cmpDiv.im / dblSq
End If
End If
Else
DivComp = CVErr(xlErrNum)
Exit Function
End If
Next cnt
Set DivComp = cmpDiv
End Function
'=============================
'Power函数。ただし指数は整数に限る
'=============================
Public Function PowerComp(ByVal varBase As Variant, ByVal varExp As Variant) As Variant
Dim cnt As Long
Dim dblTmp As Double
Dim cmpPower As Complex: Set cmpPower = New Complex
If (2 <= VarType(varExp) And VarType(varExp) <= 6) And (varExp = Int(varExp)) Then
If varExp = 0 Then
Call cmpPower.SetC(1, 0)
Else
Call cmpPower.SetC(1, 0)
For cnt = 1 To Abs(varExp)
Set cmpPower = ProductComp(cmpPower, varBase)
Next cnt
If varExp < 0 Then
Set cmpPower = DivComp(1, cmpPower)
End If
End If
Set PowerComp = cmpPower
Else
Debug.Print "error"
Exit Function
End If
End Function
'=============================
'絶対値
'=============================
Public Function AbsComp(varParam As Variant) As Variant
If (2 <= VarType(varParam) And VarType(varParam) <= 6) Then
AbsComp = Abs(varParam)
ElseIf (VarType(varParam) = 9 And TypeName(varParam) = strTypeComplex) Then
AbsComp = varParam.AbsC
Else
End If
End Function
'=============================
'多項式 f(z)=c_0 z^n + c_1 z^(n-1) + … + c_{n-2} z^2 + c_{n-1} z + c_n
'の係数c_0, c_1, ..., c_n を与えると
'その根α_1, α_2, ..., α_n を返す函数
'=============================
Public Function GetRoots(ByRef varRoots As Variant, ParamArray varCoeff() As Variant) As Boolean
Dim cnt As Long
Dim varTmp As Variant
Dim compCoeff As Variant
If UBound(varCoeff) = -1 Then
GetRoots = False
Exit Function
ElseIf UBound(varCoeff) = 0 Then
If UBound(varCoeff(0)) <= 0 Then
GetRoots = False
Exit Function
Else
If GetRoots2(varRoots, varCoeff(0)) Then
GetRoots = True
End If
End If
Else
ReDim compCoeff(UBound(varCoeff))
For cnt = LBound(varCoeff) To UBound(varCoeff)
Set compCoeff(cnt) = New Complex
Set compCoeff(cnt) = SumComp(varCoeff(cnt), 0)
Next cnt
GetRoots = True
If GetRoots2(varRoots, compCoeff) Then
GetRoots = True
End If
End If
End Function
'=============================
'多項式 f(z)=c_0 z^n + c_1 z^(n-1) + … + c_{n-2} z^2 + c_{n-1} z + c_n
'の係数c_0, c_1, ..., c_n を与えると
'その根α_1, α_2, ..., α_n をByRef引数に代入する
'根を得られればTrueを返し、得られなければFalseを返す
'=============================
Public Function GetRoots2(ByRef varRoots As Variant, ByVal varCoeff As Variant) As Boolean
Dim cnt As Long, cnt2 As Long, cnt3 As Long, cnt4 As Long
Dim varTmp As Variant
Dim compCoeff As Variant '与えられたvarCoeffがComplex型とは限らないので、引数をComplex型として格納する変数
Dim compCoeff2 As Variant, varRoots2 As Variant
Dim compRoots01 As Variant 'ニュートン法を適用する根 α(k)_i
Dim compRoots02 As Variant 'ニュートン法を適用する根 α(k+1)_i
Dim compRoots11 As Variant 'f(z) = c_0(z-α_1)(z-α_2)・・・(z-α_n)
Dim compRoots12 As Variant 'f'(z) = c_0 Π(α_i-α_j) (i≠j)
Dim compTemp As Complex: Set compTemp = New Complex
Dim varError As Variant 'ニュートン法を適用する根の差分の絶対値 ε_i = |α(k)_i - α(k+1)_i|
Dim dblError As Double 'ニュートン法を適用する根の差分の合計 Σε_i
Dim dblEpsilon As Double: dblEpsilon = cnstDblEpsilon
' On Error GoTo ErrLabal
' varCoeffの型を確認
' 配列型でなければ
If VarType(varCoeff) < vbArray Then
'引数として与えられた係数の配列が不適切なので
'処理結果をFlaseとする
GetRoots2 = False
Else
'配列の要素数が0から数えて0なら
If UBound(varCoeff) < 1 Then
'1次式ですらないので
'処理結果をFlaseとする
GetRoots2 = False
Else
ReDim compCoeff(UBound(varCoeff))
'配列の要素の型に想定外のものがないか確認
For cnt = LBound(varCoeff) To UBound(varCoeff)
Select Case VarType(varCoeff(cnt))
'実数型なら
Case vbByte, vbLong, vbSingle, vbDouble, vbCurrency
Set compCoeff(cnt) = New Complex
Call compCoeff(cnt).SetC(varCoeff(cnt), 0)
'文字列型なら
Case vbString
Set compCoeff(cnt) = New Complex
Call compCoeff(cnt).SetC_str(varCoeff(cnt))
'オブジェクト型なら
Case vbObject
'Complex型なら
If TypeName(varCoeff(cnt)) = strTypeComplex Then
Set compCoeff(cnt) = New Complex
Set compCoeff(cnt) = varCoeff(cnt)
End If
'それ以外なら
Case Else
'想定外の配列が存在したので
'処理結果をFlaseとする
GetRoots2 = False
Exit Function
End Select
Next cnt
End If
'根を格納する変数の配列数を再定義し、各要素をComplex型とする
ReDim compRoots01(UBound(compCoeff) - 1)
For cnt = LBound(compCoeff) To UBound(compCoeff) - 1
Set compRoots01(cnt) = New Complex
Next cnt
Select Case UBound(compCoeff)
'係数配列の要素数が0から数えて1なら
Case 1
'1次方程式なのでニュートン法を使うまでもなく
Set compRoots01(0) = ProductComp(DivComp(compCoeff(1), compCoeff(0)), -1)
GetRoots2 = compRoots01
'係数配列の要素数が0から数えて2以上なら
Case Else
'まずはn次方程式がn重根かどうか確認
If IsNFoldedRoots(compCoeff, compTemp) Then
For cnt = LBound(compCoeff) To UBound(compCoeff) - 1
Set compRoots01(cnt) = compTemp
Next cnt
varRoots = compRoots01
'n重根でなければ
Else
'Aberthの初期値を代入
compRoots01 = GetAberthInitial(compCoeff)
'compRoots01以外の配列数を再定義（compRoots01はCase 1（1次方程式）の場合も使用するので、Selectに入る前に再定義済み）
ReDim compRoots02(UBound(compCoeff) - 1)
ReDim compRoots11(UBound(compCoeff) - 1)
ReDim compRoots12(UBound(compCoeff) - 1)
ReDim varError(UBound(compCoeff) - 1)
For cnt = LBound(compCoeff) To UBound(compCoeff) - 1
Set compRoots02(cnt) = New Complex
Set compRoots11(cnt) = New Complex
Set compRoots12(cnt) = New Complex
varError(cnt) = CDbl(0)
Next cnt
'ループ処理
Do
'ループ回数を数える
cnt3 = cnt3 + 1
dblError = 0
For cnt = LBound(compCoeff) To UBound(compCoeff) - 1
Call compRoots11(cnt).SetC(0, 0)
For cnt2 = LBound(compCoeff) To UBound(compCoeff)
Set compRoots11(cnt) = SumComp(compRoots11(cnt), ProductComp(compCoeff(cnt2), PowerComp(compRoots01(cnt), UBound(compCoeff) - cnt2)))
Next cnt2
Next cnt
For cnt = LBound(compCoeff) To UBound(compCoeff) - 1
Set compRoots12(cnt) = compCoeff(0)
For cnt2 = LBound(compCoeff) To UBound(compCoeff) - 1
If cnt2 <> cnt Then Set compRoots12(cnt) = ProductComp(compRoots12(cnt), SumComp(compRoots01(cnt), ProductComp(-1, compRoots01(cnt2))))
Next cnt2
Set compRoots02(cnt) = SumComp(compRoots01(cnt), ProductComp(-1, DivComp(compRoots11(cnt), compRoots12(cnt))))
Next cnt
'ループ処理の繰り返しによる根の収束の確認
For cnt = LBound(compCoeff) To UBound(compCoeff) - 1
'ループ処理cnt3-1回目の根とループ処理cnt3回目の根との差分を取る
Set compTemp = SubComp(compRoots01(cnt), compRoots02(cnt))
'差分の絶対値
varError(cnt) = compTemp.AbsC
'n個の根について差分を合計
dblError = dblError + varError(cnt)
'ループ処理cnt3-1回目の根を削除して、そこへループ処理cnt3回目の根を代入
Set compRoots01(cnt) = compRoots02(cnt)
Next cnt
'反復処理回数が多い場合の処理
If cnt3 Mod 100 = 0 Then
For cnt = LBound(varError) To UBound(varError)
If varError(cnt) < dblEpsilon Then
compRoots01(cnt).re = Round(compRoots01(cnt).re, 12)
compRoots01(cnt).im = Round(compRoots01(cnt).im, 12)
Call PolyN2Nminus1(compRoots01(cnt), compCoeff, compCoeff2)
If GetRoots2(varRoots2, compCoeff2) Then
For cnt2 = LBound(varError) To UBound(varError)
If cnt2 <> cnt Then
Set compRoots01(cnt2) = varRoots2(cnt4)
cnt4 = cnt4 + 1
End If
Next cnt2
varRoots = compRoots01
GetRoots2 = True
Exit Function
End If
End If
Next cnt
If cnt3 > 500 Then
Debug.Print "Error: The number of loops exceeded the specified number."
Call RootsRounding(compRoots01, compCoeff)
varRoots = compRoots01
GetRoots2 = False
Exit Function
End If
dblEpsilon = dblEpsilon * 10
End If
Loop While dblError > dblEpsilon
Call RootsRounding(compRoots01, compCoeff)
varRoots = compRoots01
End If
End Select
'処理結果をTrueとする
GetRoots2 = True
End If
'エラー処理に進入しないようExit
Exit Function
'
'ErrLabal:
' Debug.Print "An error occurred."
' GetRoots2 = False
End Function
'=============================
'Aberthの初期値の半径
'http://nalab.mind.meiji.ac.jp/~mk/lecture/ouyoukaiseki4/notebook/nonlinear-algebra/node40.html
'=============================
Private Function GetAberthRadius(ByVal varCoeff As Variant) As Variant
Dim cnt As Long
Dim varTmp As Variant
Dim dblTmp As Double
ReDim varTmp(UBound(varCoeff))
For cnt = LBound(varCoeff) To UBound(varCoeff)
Set varTmp(cnt) = DivComp(varCoeff(cnt), varCoeff(LBound(varCoeff)))
If dblTmp < varTmp(cnt).AbsC Then dblTmp = varTmp(cnt).AbsC
Next cnt
GetAberthRadius = varTmp(1).AbsC / (UBound(varCoeff) - 1) + 1 + dblTmp
End Function
'=============================
'Aberthの初期値
'http://nalab.mind.meiji.ac.jp/~mk/lecture/ouyoukaiseki4/notebook/nonlinear-algebra/node40.html
'=============================
Private Function GetAberthInitial(ByVal varCoeff As Variant) As Variant
Dim cnt As Long
Dim varTmp As Variant
Dim compTmp As Complex: Set compTmp = New Complex
Dim dblRadius As Double
dblRadius = GetAberthRadius(varCoeff)
ReDim varTmp(UBound(varCoeff) - 1)
For cnt = LBound(varCoeff) To UBound(varCoeff) - 1
Call compTmp.SetC_RT(dblRadius, (2 * cnt + 0.5) * PI / UBound(varCoeff))
' ↑Aberthの初期値、↓伊理の初期値
' Call compTmp.SetC_RT(dblRadius, (2 * PI * cnt + 1.5) / UBound(varCoeff))
Set compTmp = SumComp(compTmp, DivComp(ProductComp(-1, varCoeff(1)), ProductComp(varCoeff(0), UBound(varCoeff))))
Set varTmp(cnt) = compTmp
Set compTmp = New Complex
Next cnt
GetAberthInitial = varTmp
End Function
'=============================
' 複素数Classの配列を、実数のものだけ実数型に置き換える
'=============================
Private Sub IfRealThenComp2Real(ByRef varBf, ByRef varAf)
If Not (LBound(varBf) = LBound(varAf) And UBound(varBf) = UBound(varAf)) Then
Exit Sub
End If
Dim cnt As Long
For cnt = LBound(varBf) To UBound(varBf)
If (2 <= VarType(varBf(cnt)) And VarType(varBf(cnt)) <= 6) Then
varAf(cnt) = varBf(cnt)
ElseIf (VarType(varBf(cnt)) = 9 And TypeName(varBf(cnt)) = strTypeComplex) Then
If varBf(cnt).im = 0 Then
varAf(cnt) = varBf(cnt).re
Else
Set varAf(cnt) = varBf(cnt)
End If
Else
Debug.Print "error"
Exit Sub
End If
Next cnt
End Sub
'=============================
'複素数Classの配列を、文字列に置き換える
'=============================
Private Sub RealAndComp2String(ByRef varBf, ByRef varAf)
If Not (LBound(varBf) = LBound(varAf) And UBound(varBf) = UBound(varAf)) Then
Exit Sub
End If
Dim cnt As Long
For cnt = LBound(varBf) To UBound(varBf)
If (2 <= VarType(varBf(cnt)) And VarType(varBf(cnt)) <= 6) Then
varAf(cnt) = Str(varBf(cnt))
ElseIf (VarType(varBf(cnt)) = 9 And TypeName(varBf(cnt)) = strTypeComplex) Then
If varBf(cnt).im = 0 Then
varAf(cnt) = Str(varBf(cnt).re)
Else
varAf(cnt) = varBf(cnt).CompStr
End If
Else
Debug.Print "error"
Exit Sub
End If
Next cnt
End Sub
'=============================
'-1.999999のような根が-2でよいか確認するSub Procedure
'=============================
Private Sub RootsRounding(ByRef varRoots As Variant, ByVal varCoeff As Variant, Optional ByVal byteDigit As Byte = 7)
Dim varRootsRounded As Variant
Dim cnt As Long, cnt2 As Long
Dim compTemp As Complex
ReDim varRootsRounded(UBound(varRoots))
For cnt = LBound(varRoots) To UBound(varRoots)
Set varRootsRounded(cnt) = New Complex
varRootsRounded(cnt).re = Round(varRoots(cnt).re, byteDigit)
varRootsRounded(cnt).im = Round(varRoots(cnt).im, byteDigit)
Next cnt
Set compTemp = New Complex
For cnt2 = LBound(varRoots) To UBound(varRoots)
Call compTemp.SetC(0, 0)
For cnt = LBound(varCoeff) To UBound(varCoeff)
Set compTemp = SumComp(compTemp, ProductComp(varCoeff(cnt), PowerComp(varRootsRounded(cnt2), UBound(varCoeff) - cnt)))
Next cnt
If Round(compTemp.re, byteDigit) = 0 And Round(compTemp.im, byteDigit) = 0 Then
Set varRoots(cnt2) = varRootsRounded(cnt2)
End If
Next cnt2
End Sub
'=============================
'n重根の確認
'=============================
Private Function IsNFoldedRoots(ByVal varCoeff As Variant, ByRef compRoot As Complex, Optional ByVal byteDigit As Byte = 7) As Boolean
Dim compNFoldedRoots As Complex: Set compNFoldedRoots = New Complex
Dim cnt As Long
Dim compTemp As Complex: Set compTemp = New Complex
'係数配列の要素数が0から数えて2（2次方程式）未満であれば
If UBound(varCoeff) < 2 Then
'重根判定にFalseを返す
IsNFoldedRoots = False
'係数配列の要素数が0から数えて2以上であれば
Else
'パスカルの三角形の左から2つめ（n次多項式ならn-1。配列の要素はゼロから数えるので配列の要素数）に基づき
'c_0(x-(a+bi))^n=c_0 x^n+c_0 (n-1) (-(a+bi)) x^(n-1)+…
'c_1=c_0 (n-1) (-(a+bi)) ⇔ (a+bi) = c_1 / c_0 / (-(n-1))
'としてn重根を推定
Set compNFoldedRoots = DivComp(varCoeff(1), varCoeff(0), -1 * UBound(varCoeff))
'推定したn重根を多項式に当てはめて数値計算
For cnt = LBound(varCoeff) To UBound(varCoeff)
Set compTemp = SumComp(compTemp, ProductComp(varCoeff(cnt), PowerComp(compNFoldedRoots, UBound(varCoeff) - cnt)))
Next cnt
'推定したn重根を多項式に当てはめた結果がゼロとなっているか確認
If Round(compTemp.re, byteDigit) = 0 And Round(compTemp.im, byteDigit) = 0 Then
'与えられているByRef引数にn重根を代入
Set compRoot = New Complex
Set compRoot = compNFoldedRoots
'重根判定にTrueを返す
IsNFoldedRoots = True
Else
'重根判定にFalseを返す
IsNFoldedRoots = False
End If
End If
End Function
'=============================
'n次多項式の根の1つが判明している場合に、残りのn-1個の根を持つn-1次多項式を返す
'=============================
Public Function PolyN2Nminus1(ByVal compRoot As Variant, ByVal varCoeffIn, ByRef varCoeffOut, Optional ByVal byteDigit As Byte = 7) As Boolean
Dim cnt As Long
Dim compTemp As Complex: Set compTemp = New Complex
'n次多項式の係数配列の要素数が0から数えて2（2次方程式）未満であれば
If UBound(varCoeffIn) < 2 Then
PolyN2Nminus1 = False
Else
'与えられた根を多項式に当てはめて数値計算
For cnt = LBound(varCoeffIn) To UBound(varCoeffIn)
Set compTemp = SumComp(compTemp, ProductComp(varCoeffIn(cnt), PowerComp(compRoot, UBound(varCoeffIn) - cnt)))
Next cnt
'数値計算結果がゼロなら（与えられた根が真に根なら）
If Round(compTemp.re, byteDigit) = 0 And Round(compTemp.im, byteDigit) = 0 Then
'残りのn-1個の根を持つn-1次多項式の係数をByRef引数に返す
Call compTemp.SetC(0, 0)
ReDim varCoeffOut(UBound(varCoeffIn) - 1)
For cnt = LBound(varCoeffIn) To UBound(varCoeffIn) - 1
Set compTemp = ProductComp(compTemp, compRoot)
Set compTemp = SumComp(compTemp, varCoeffIn(cnt))
Set varCoeffOut(cnt) = New Complex
Set varCoeffOut(cnt) = compTemp
Next cnt
PolyN2Nminus1 = True
'数値計算結果がゼロでないなら（与えられた根が真に根でないなら）
Else
PolyN2Nminus1 = False
End If
End If
End Function
'=============================
'WorkSheet函数用
'=============================
Public Function WSFuncGetRoots(rngCoeff As Range) As Variant
Dim strCoeff() As String
Dim compCoeff() As Complex
Dim cnt As Long
Dim varRoots As Variant
If rngCoeff.Rows.Count > 1 And rngCoeff.Columns.Count > 1 Then
WSFuncGetRoots = CVErr(xlErrNum)
ElseIf rngCoeff.Rows.Count = 1 Then
ReDim strCoeff(rngCoeff.Columns.Count - 1)
For cnt = 0 To rngCoeff.Columns.Count - 1
strCoeff(cnt) = rngCoeff(1, cnt + 1)
Next cnt
ElseIf rngCoeff.Columns.Count = 1 Then
ReDim strCoeff(rngCoeff.Rows.Count - 1)
For cnt = 0 To rngCoeff.Rows.Count - 1
strCoeff(cnt) = rngCoeff(cnt + 1, 1)
Next cnt
End If
ReDim compCoeff(UBound(strCoeff))
For cnt = LBound(compCoeff) To UBound(compCoeff)
Set compCoeff(cnt) = New Complex
Call compCoeff(cnt).SetC_str(strCoeff(cnt))
Next cnt
Dim rt As Variant
Call GetRoots(rt, compCoeff)
ReDim varRoots(1 To UBound(rt) + 1)
For cnt = LBound(rt) To UBound(rt)
varRoots(cnt + 1) = rt(cnt).CompStr
Next cnt
WSFuncGetRoots = varRoots
End Function
計算例・使い方
テスト用モジュール
Module1.bas
Option Explicit
Public Sub test001()
Dim cnt As Long
Dim a0 As Complex: Set a0 = New Complex
Dim a1 As Complex: Set a1 = New Complex
Dim a2 As Complex: Set a2 = New Complex
Dim a3 As Complex: Set a3 = New Complex
Dim a4 As Complex: Set a4 = New Complex
Dim a5 As Complex: Set a5 = New Complex
Dim a6 As Complex: Set a6 = New Complex
Dim a7 As Complex: Set a7 = New Complex
Dim a8 As Complex: Set a8 = New Complex
Dim rt As Variant
Set rt = Nothing
Debug.Print Chr(10) & "x^2+1=0"
If GetRoots(rt, 1, 0, 1) Then
Debug.Print "収束しました"
Else
Debug.Print "収束しませんでした"
End If
For cnt = LBound(rt) To UBound(rt)
Debug.Print rt(cnt).CompStr
Next cnt
Call a1.SetC(3, 2)
Call a2.SetC(5, 0)
Call a3.SetC(2, 1)
Set rt = Nothing
Debug.Print Chr(10) & "(3+2i)x^2+5x+(2+i)=0"
If GetRoots(rt, a1, a2, a3) Then
Debug.Print "収束しました"
Else
Debug.Print "収束しませんでした"
End If
For cnt = LBound(rt) To UBound(rt)
Debug.Print rt(cnt).CompStr
Next cnt
Call a1.SetC(0, -8)
Set rt = Nothing
Debug.Print Chr(10) & "(x^3-8i=0"
If GetRoots(rt, 1, 0, 0, a1) Then
Debug.Print "収束しました"
Else
Debug.Print "収束しませんでした"
End If
For cnt = LBound(rt) To UBound(rt)
Debug.Print rt(cnt).CompStr
Next cnt
Set rt = Nothing
If GetRoots(rt, 33, 140, 266, 140, 33) Then
Debug.Print "収束しました"
Else
Debug.Print "収束しませんでした"
End If
Debug.Print Chr(10) & "33x^4+140x^3+266x^2+140x+33=0"
For cnt = LBound(rt) To UBound(rt)
Debug.Print rt(cnt).CompStr
Next cnt
Call a0.SetC(1, 0)
Call a1.SetC(0, -16)
Call a2.SetC(-112, 0)
Call a3.SetC(0, 448)
Call a4.SetC(1120, 0)
Call a5.SetC(0, -1792)
Call a6.SetC(-1792, 0)
Call a7.SetC(0, 1024)
Call a8.SetC(256, 0)
Set rt = Nothing
Debug.Print Chr(10) & "(x-2i)^8=0"
If GetRoots(rt, a0, a1, a2, a3, a4, a5, a6, a7, a8) Then
Debug.Print "収束しました"
Else
Debug.Print "収束しませんでした"
End If
For cnt = LBound(rt) To UBound(rt)
Debug.Print rt(cnt).CompStr
Next cnt
Call a0.SetC(1, 0)
Call a1.SetC(-16, -24)
Call a2.SetC(-140, 336)
Call a3.SetC(2576, -504)
Call a4.SetC(-8330, -8400)
Call a5.SetC(-6832, 33432)
Call a6.SetC(56980, -23184)
Call a7.SetC(-52432, -35592)
Call a8.SetC(-239, 28560)
Set rt = Nothing
Debug.Print Chr(10) & "{x-(2+3i)}^8=0"
If GetRoots(rt, a0, a1, a2, a3, a4, a5, a6, a7, a8) Then
Debug.Print "収束しました"
Else
Debug.Print "収束しませんでした"
End If
For cnt = LBound(rt) To UBound(rt)
Debug.Print rt(cnt).CompStr
Next cnt
Call a0.SetC(1, 0)
Call a1.SetC(-15, -21)
Call a2.SetC(-91, 273)
Call a3.SetC(1715, -567)
Call a4.SetC(-5775, -3885)
Call a5.SetC(1603, 16737)
Call a6.SetC(16807, -18333)
Call a7.SetC(-20799, 1347)
Call a8.SetC(6554, 4449)
Set rt = Nothing
Debug.Print Chr(10) & "{x-(2+3i)}^7 (x-1)=0"
If GetRoots(rt, a0, a1, a2, a3, a4, a5, a6, a7, a8) Then
Debug.Print "収束しました"
Else
Debug.Print "収束しませんでした"
End If
For cnt = LBound(rt) To UBound(rt)
Debug.Print rt(cnt).CompStr
Next cnt
Call a0.SetC(1, 0)
Call a1.SetC(0, -8 * Sqr(3))
Call a2.SetC(-84, 0)
Call a3.SetC(0, 168 * Sqr(3))
Call a4.SetC(630, 0)
Call a5.SetC(0, -504 * Sqr(3))
Call a6.SetC(-756, 0)
Call a7.SetC(0, 216 * Sqr(3))
Call a8.SetC(81, 0)
Set rt = Nothing
Debug.Print Chr(10) & "{x-sqrt(3)i}^8=0"
If GetRoots(rt, a0, a1, a2, a3, a4, a5, a6, a7, a8) Then
Debug.Print "収束しました"
Else
Debug.Print "収束しませんでした"
End If
For cnt = LBound(rt) To UBound(rt)
Debug.Print rt(cnt).CompStr
Next cnt
End Sub
テスト用モジュール実行結果
Module1_Result.bas
'x^2+1=0
'収束しました
'i
'-i
'(3+2i)x^2+5x+(2+i)=0
'収束しました
'-.836846641040605+1.09033036388437i
'-.316999512805549-.3210995946536i
'(x^3-8i=0
'収束しました
'1.73205080756888+i
'-1.73205080756888+i
'-2i
'収束しました
'33x^4+140x^3+266x^2+140x+33=0
'-.323693680587964+.274409342564836i
'-1.79751844062416+1.5238352897217i
'-1.79751844062416-1.5238352897217i
'-.323693680587964-.274409342564836i
'(x-2i)^8=0
'収束しました
'2i
'2i
'2i
'2i
'2i
'2i
'2i
'2i
'{x-(2+3i)}^8=0
'収束しました
'2+3i
'2+3i
'2+3i
'2+3i
'2+3i
'2+3i
'2+3i
'2+3i
'{x-(2+3i)}^7 (x-1)=0
'収束しました
'2+3i
'2+3i
'2+3i
'2+3i
'2+3i
'1
'2+3i
'2+3i
'{x-sqrt(3)i}^8=0
'収束しました
'1.73205080756888i
'1.73205080756888i
'1.73205080756888i
'1.73205080756888i
'1.73205080756888i
'1.73205080756888i
'1.73205080756888i
'1.73205080756888i