$ 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