Excel VBAでDKA法（Durand-Kerner-Aberth法）

$ 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

'=============================
'複素数の定義(実部と虚部) z=a+bi = z.SetC(a, b)
'=============================
Public Function SetC(a As Double, b As Double)
  re = a
  im = b
End Function

'=============================
' 複素数の定義(絶対値と偏角) z=re^{iΘ} = z.SetC_RT(r, Θ)
'=============================
Public Function SetC_RT(r As Double, theta As Double)
  re = r * Cos(theta)
  im = r * Sin(theta)
End Function

'=============================
'絶対値 |z| = sqrt(a^2+b^2) = z.AbsC
'=============================
Public Function AbsC() As Double
  AbsC = Sqr(re ^ 2 + im ^ 2)
End Function

'=============================
'絶対値の二乗 |z|^2 = a^2+b^2 = z.AbsCSq
'=============================
Public Function AbsCSq() As Double
  AbsCSq = re ^ 2 + im ^ 2
End Function
 
'=============================
'　偏角 ∠z = arctan(b/a) = z.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

'=============================
'　足し算 z_1+z_2 = z_1.Plus(z_2)
'=============================
Public Function Plus(ByVal z_2 As Variant) As Variant
  If (2 <= VarType(z_2) And VarType(z_2) <= 6) Then
    Set Plus = New Complex
    Plus.re = re + z_2
    Plus.im = im
  ElseIf (VarType(z_2) = 9 And TypeName(z_2) = strClassName) Then
    Set Plus = New Complex
    Plus.re = re + z_2.re
    Plus.im = im + z_2.im
  Else
    Plus = CVErr(xlErrNum)
  End If
End Function

'=============================
'　引き算 z_1-z_2 = z_1.Minus(z_2)
'=============================
Public Function Minus(ByVal z_2 As Variant) As Variant
  If (2 <= VarType(z_2) And VarType(z_2) <= 6) Then
    Set Minus = New Complex
    Minus.re = re - z_2
    Minus.im = im
  ElseIf (VarType(z_2) = 9 And TypeName(z_2) = strClassName) Then
    Set Minus = New Complex
    Minus.re = re - z_2.re
    Minus.im = im - z_2.im
  Else
    Minus = CVErr(xlErrNum)
  End If
End Function

'=============================
'　掛け算 z_1×z_2 = z_1.Times(z_2)
'=============================
Public Function Times(ByVal b As Variant) As Variant
  If (2 <= VarType(z_2) And VarType(z_2) <= 6) Then
    Set Times = New Complex
    Times.re = re * z_2
    Times.im = im * z_2
  ElseIf (VarType(z_2) = 9 And TypeName(z_2) = strClassName) Then
    Set Times = New Complex
    Times.re = re * z_2.re - im * z_2.im
    Times.im = re * z_2.im + im * z_2.re
  Else
    Times = CVErr(xlErrNum)
  End If
End Function

'============================
'　割り算 z_1 / z_2 = z_1.Div(z_2)
'============================
Public Function Div(ByVal z_2 As Variant) As Variant
  If (2 <= VarType(z_2) And VarType(b) <= 6) Then
    If z_2 = 0 Then
      Div = CVErr(xlErrNum)
    Else
      Set Div = New Complex
      Div.re = re / z_2
      Div.im = im / z_2
    End If
  ElseIf (VarType(z_2) = 9 And TypeName(z_2) = strClassName) Then
    If z_2.AbsC() = 0 Then
      Div = CVErr(xlErrNum)
    Else
      Set Div = New Complex
      Div.re = (re * z_2.re + im * z_2.im) / (z_2.AbsC() ^ 2)
      Div.im = (-re * z_2.im + im * z_2.re) / (z_2.AbsC() ^ 2)
    End If
  Else
    Div = CVErr(xlErrNum)
  End If
End Function

'============================
'共役複素数 z = a+bi → z.CompCnj = a-bi
'============================
Public Function CompCnj() As Complex
  Set CompCnj = New Complex
  CompCnj.re = re
  CompCnj.im = (-1) * im
End Function

'============================
'文字列(直交座標) z = a + bi → z.CompStr = "a+bi"
'============================
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Θ) z = r Exp(Θi) → z.CompStr = "r Exp(Θi)"
'============================
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
  
  If dblAbs = 1 Then
    CompStrPolar = "Exp(" & Str(dblAngle) & "i)"
  Else
    CompStrPolar = Str(dblAbs) & " Exp(" & Str(dblAngle) & "i)"
  End If
End Function

標準モジュール

各種計算用の函数を定義する標準モジュール

ComplexFunc.bas

Option Explicit

Private Const strTypeComplex As String = "Complex"
Private Const dblEpsilon 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

'=============================
'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(ParamArray varCoeff() As Variant) As Variant
  Dim cnt As Long, cnt2 As Long, cnt3 As Long
  Dim varTmp 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 dblError As Double

  If UBound(varCoeff) = -1 Then
    GetRoots = CVErr(xlErrNum)
    Exit Function
  ElseIf UBound(varCoeff) = 0 Then
    If UBound(varCoeff(0)) <= 0 Then
      GetRoots = CVErr(xlErrNum)
      Exit Function
    Else
      For cnt = LBound(varCoeff(0)) To UBound(varCoeff(0))
        If Not ((2 <= VarType(varCoeff(0)(cnt)) And VarType(varCoeff(0)(cnt)) <= 6) Or (VarType(varCoeff(0)(cnt)) = 9 And TypeName(varCoeff(0)(cnt)) = strTypeComplex)) Then
          GetRoots = CVErr(xlErrNum)
          Exit Function
        End If
      Next cnt
      ReDim compRoots01(UBound(varCoeff(0)) - 1)
      ReDim varTmp(UBound(varCoeff(0)) - 1)
      For cnt = LBound(varCoeff(0)) To UBound(varCoeff(0)) - 1
        Set compRoots01(cnt) = New Complex
      Next cnt
      Select Case UBound(varCoeff(0))
        Case 1
          Set compRoots01(0) = ProductComp(DivComp(varCoeff(0)(1), varCoeff(0)(0)), -1)
'          Call RealAndComp2String(compRoots01, varTmp)
'          GetRoots = varTmp
          GetRoots = compRoots01
        Case Else
          'https://www.chart.co.jp/subject/sugaku/suken_tsushin/84/84-5.pdf
          compRoots01 = GetAberthInitial(varCoeff(0))
          ReDim compRoots02(UBound(varCoeff(0)) - 1) 'compRoots01は Case1でも使うので、これよりも前に再定義済み
          ReDim compRoots11(UBound(varCoeff(0)) - 1)
          ReDim compRoots12(UBound(varCoeff(0)) - 1)
          For cnt = LBound(varCoeff(0)) To UBound(varCoeff(0)) - 1
            Set compRoots02(cnt) = New Complex
            Set compRoots11(cnt) = New Complex
            Set compRoots12(cnt) = New Complex
          Next cnt
          Do
            cnt3 = cnt3 + 1
            dblError = 0
            For cnt = LBound(varCoeff(0)) To UBound(varCoeff(0)) - 1
              Call compRoots11(cnt).SetC(0, 0)
              For cnt2 = LBound(varCoeff(0)) To UBound(varCoeff(0))
                Set compRoots11(cnt) = SumComp(compRoots11(cnt), ProductComp(varCoeff(0)(cnt2), PowerComp(compRoots01(cnt), UBound(varCoeff(0)) - cnt2)))
              Next cnt2
            Next cnt
            For cnt = LBound(varCoeff(0)) To UBound(varCoeff(0)) - 1
              If (2 <= VarType(varCoeff(0)(0)) And VarType(varCoeff(0)(0)) <= 6) Then
                Call compTemp.SetC(Val(varCoeff(0)(0)), 0)
                Set compRoots12(cnt) = compTemp
              ElseIf (VarType(varCoeff(0)(0)) = 9 And TypeName(varCoeff(0)(0)) = strTypeComplex) Then
                Set compRoots12(cnt) = varCoeff(0)(0)
              Else
                Debug.Print "error"
                Exit Function
              End If
                
              For cnt2 = LBound(varCoeff(0)) To UBound(varCoeff(0)) - 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(varCoeff(0)) To UBound(varCoeff(0)) - 1
              Set compTemp = compRoots01(cnt).Minus(compRoots02(cnt))
              dblError = dblError + compTemp.AbsC
              Set compRoots01(cnt) = compRoots02(cnt)
            Next cnt
            If cnt3 Mod 1000 = 0 Then
              Stop
            End If
          Loop While dblError > dblEpsilon
          GetRoots = compRoots01
      End Select
    End If
  Else
    For cnt = LBound(varCoeff) To UBound(varCoeff)
      If Not ((2 <= VarType(varCoeff(cnt)) And VarType(varCoeff(cnt)) <= 6) Or (VarType(varCoeff(cnt)) = 9 And TypeName(varCoeff(cnt)) = strTypeComplex)) Then
        GetRoots = CVErr(xlErrNum)
        Exit Function
      End If
    Next cnt
    ReDim compRoots01(UBound(varCoeff) - 1)
    ReDim varTmp(UBound(varCoeff) - 1)
    For cnt = LBound(varCoeff) To UBound(varCoeff) - 1
      Set compRoots01(cnt) = New Complex
    Next cnt
    Select Case UBound(varCoeff)
      Case 1
        Set compRoots01(0) = ProductComp(DivComp(varCoeff(1), varCoeff(0)), -1)
        GetRoots = compRoots01
      Case Else
        compRoots01 = GetAberthInitial(varCoeff)
        ReDim compRoots02(UBound(varCoeff) - 1) 'compRoots01は Case1でも使うので、これよりも前に再定義済み
        ReDim compRoots11(UBound(varCoeff) - 1)
        ReDim compRoots12(UBound(varCoeff) - 1)
        For cnt = LBound(varCoeff) To UBound(varCoeff) - 1
          Set compRoots02(cnt) = New Complex
          Set compRoots11(cnt) = New Complex
          Set compRoots12(cnt) = New Complex
        Next cnt
        Do
          cnt3 = cnt3 + 1
          dblError = 0
          For cnt = LBound(varCoeff) To UBound(varCoeff) - 1
            Call compRoots11(cnt).SetC(0, 0)
            For cnt2 = LBound(varCoeff) To UBound(varCoeff)
              Set compRoots11(cnt) = SumComp(compRoots11(cnt), ProductComp(varCoeff(cnt2), PowerComp(compRoots01(cnt), UBound(varCoeff) - cnt2)))
            Next cnt2
          Next cnt
          For cnt = LBound(varCoeff) To UBound(varCoeff) - 1
            If (2 <= VarType(varCoeff(0)) And VarType(varCoeff(0)) <= 6) Then
              Call compTemp.SetC(Val(varCoeff(0)), 0)
              Set compRoots12(cnt) = compTemp
            ElseIf (VarType(varCoeff(0)) = 9 And TypeName(varCoeff(0)) = strTypeComplex) Then
              Set compRoots12(cnt) = varCoeff(0)
            Else
              Debug.Print "error"
              Exit Function
            End If
              
            For cnt2 = LBound(varCoeff) To UBound(varCoeff) - 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(varCoeff) To UBound(varCoeff) - 1
            Set compTemp = compRoots01(cnt).Minus(compRoots02(cnt))
            dblError = dblError + compTemp.AbsC
            Set compRoots01(cnt) = compRoots02(cnt)
          Next cnt
          If cnt3 Mod 1000 = 0 Then
            Stop
          End If
        Loop While dblError > dblEpsilon
        GetRoots = compRoots01
    End Select
  End If
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))
    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

'=============================
'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
  rt = GetRoots(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 a1 As Complex:  Set a1 = New Complex
  Dim a2 As Complex:  Set a2 = New Complex
  Dim a3 As Complex:  Set a3 = New Complex

  Dim rt As Variant

  Set rt = Nothing
  rt = GetRoots(1, 0, 1)
  For cnt = LBound(rt) To UBound(rt)
    Debug.Print rt(cnt).CompStr
  Next cnt
  '+i
  '-i

  Call a1.SetC(3, 2)
  Call a2.SetC(5, 0)
  Call a3.SetC(2, 1)
  Set rt = Nothing
  rt = GetRoots(a1, a2, a3)
  For cnt = LBound(rt) To UBound(rt)
    Debug.Print rt(cnt).CompStr
  Next cnt
  '-.836846641040605+1.09033036388437i
  '-.316999512805549-.3210995946536i

  Call a1.SetC(0, -8)
  Set rt = Nothing
  rt = GetRoots(1, 0, 0, a1)
  For cnt = LBound(rt) To UBound(rt)
    Debug.Print rt(cnt).CompStr
  Next cnt
  '1.73205080756888+i
  '-1.73205080756888+i
  '-2i

  Set rt = Nothing
  rt = GetRoots(33, 140, 266, 140, 33)
  For cnt = LBound(rt) To UBound(rt)
    Debug.Print rt(cnt).CompStr
  Next cnt
  '-.323693680587964+.274409342564836i
  '-1.79751844062416+1.5238352897217i
  '-1.79751844062416-1.5238352897217i
  '-.323693680587964-.274409342564836i

End Sub