Attribute VB_Name = "Module1"
Public Function LambertW(x As Double)
  ' Lambert W function for x (real number)
  ' real result, domain: [-1/e , infinity]
  ' complex result, domain: [infinity , infinity]
  ' return: Complex number (in the form a+bi)
    
  LambertW = LambertWc(x, 0)
End Function
Public Function LambertWc(a As Double, b As Double)
  ' Lambert W function for a+bi
  ' real result, domain: [-1/e , infinity]
  ' complex result, domain: [infinity , infinity]
  ' return: Complex number (in the form a+bi)
  Dim z, x, z1, z2, z3, q1
  If a = 0 And b = 0 Then
    z = 0
  Else
    ' complex parameter
    x = WorksheetFunction.Complex(a, b)

    ' approx initial value = sqtr(2*(1+e*x))-1
    z = WorksheetFunction.ImProduct(x, WorksheetFunction.ImExp(1))
    z = WorksheetFunction.ImProduct(2, WorksheetFunction.ImSum(1, z))
    If WorksheetFunction.ImReal(z) = 0 And WorksheetFunction.Imaginary(z) = 0 Then
      z = -1
    Else
      z = WorksheetFunction.ImSub(WorksheetFunction.ImSqrt(z), 1)
      For i = 1 To 100
        ' the comments are right to left evaluation
        ' z1 = 0 sub z sub ln x div z
        z1 = WorksheetFunction.ImSub(z, WorksheetFunction.ImLn(WorksheetFunction.ImDiv(x, z)))
        z1 = WorksheetFunction.ImSub(0, z1)
        ' q1 = 2 mul (1 add z) mul 1 add z add 2 mul z1 div 3
        z2 = WorksheetFunction.ImProduct(2, WorksheetFunction.ImDiv(z1, 3))
        z2 = WorksheetFunction.ImSum(1, WorksheetFunction.ImSum(z, z2))
        z2 = WorksheetFunction.ImProduct(WorksheetFunction.ImSum(z, 1), z2)
        q1 = WorksheetFunction.ImProduct(2, z2)
        ' z1 = z mul 1 add (z1 div 1 add z) mul (q1 sub z1) div q1 sub 2 mul z1
        z2 = WorksheetFunction.ImSub(q1, WorksheetFunction.ImProduct(2, z1))
        z2 = WorksheetFunction.ImDiv(WorksheetFunction.ImSub(q1, z1), z2)
        z3 = WorksheetFunction.ImDiv(z1, WorksheetFunction.ImSum(1, z))
        z2 = WorksheetFunction.ImSum(1, WorksheetFunction.ImProduct(z3, z2))
        z1 = WorksheetFunction.ImProduct(z, z2)
        If z = z1 Then
          Exit For
        Else
          z = z1
        End If
      Next
    End If
  End If
  LambertWc = z
End Function