DECLARE SUB CDIV (z1#(), z2#(), z#())
DECLARE SUB CMUL (z1#(), z2#(), z#())
DECLARE SUB CSQRT (z#(), z1#())
DECLARE SUB CUBRT (z#(), z1#(), z2#(), z3#())
DECLARE SUB F (z#(), z1#())
DECLARE FUNCTION CABS# (z#())
DECLARE SUB ATAN (Numerator AS DOUBLE, denominator#, Phase#)
'***********************************************************
'*          Solve complex equation of degree 3:            *
'*             a z^3 + b z^2 + c z + d = 0                 *
'* ------------------------------------------------------- *
'* SAMPLE RUN:                                             *
'*             a = (1,0) b = (2,-10) c = (-4,1) d = (5,2)  *
'*                                                         *
'* Equation az^3 + bz^2 + cz + d = 0                       *
'* Solutions z:                                            *
'* z = ( -0.47030,  0.66849)                               *
'* F(z) = ( -0.00000,  0.00000)                            *
'* z = (  0.60178, -0.29111)                               *
'* F(z) = (  0.00000, -0.00000)                            *
'* z = ( -2.13148,  9.62262)                               *
'* F(z) = ( -0.00000,  0.00000)                            *
'*                                                         *
'*                                     J-P Moreau, Paris.  *
'*                                     (www.jpmoreau.fr)   *
'***********************************************************
' NOTE:
' z = Z - (b/3a)
' p = (c/a) - (b²/3a²)
' q = (2b^3/27a^3) + (d/a) - (bc)/3a²)
'
' Now the equation is reduced to Z^3 + pZ + q = 0
'
' with the roots:
'
' Z[1..3] = cubroot(-q/2 + sqrt(q²/4 + p^3/27)) + cubroot(-q/2 - sqrt(q²/4 + p^3/27))
'------------------------------------------------------------------------------------
'Program CRoot3
DEFDBL A-H, O-Z
DEFINT I-N

'Complex numbers
DIM SHARED a(2), b(2), c(2), d(2)
DIM SHARED u(2), v(2), w(2)
DIM det(2), p(2), q(2), sdet(2)
DIM u1(2), u2(2), v1(2), v2(2), z1(2), zz(2)
DIM z(9, 2)

' seek three complex roots of complex equation:
'  a z^3 + b z^2 + c z + d = 0

' define complex coefficients a, b, c, d
  a(1) = 1#: a(2) = 0#
  b(1) = 2#: b(2) = -10#
  c(1) = -4#: c(2) = 1#
  d(1) = 5#: d(2) = 2#

' calculate p = c/a - b²/3a²
  CDIV c(), a(), u(): CMUL b(), b(), v(): CMUL a(), a(), w()
  w(1) = 3# * w(1): w(2) = 3# * w(2): CDIV v(), w(), z1()
  p(1) = u(1) - z1(1): p(2) = u(2) - z1(2)

' calculate q = 2b^3/27a^3 + d/a - bc/3a²
  CMUL b(), v(), w(): w(1) = 2# * w(1): w(2) = 2# * w(2)
  CMUL a(), a(), u(): CMUL u(), a(), v()
  v(1) = 27# * v(1): v(2) = 27# * v(2)
  CDIV w(), v(), q()
  CDIV d(), a(), u()
  q(1) = q(1) + u(1): q(2) = q(2) + u(2)
  CMUL b(), c(), u(): CMUL a(), a(), v()
  v(1) = 3# * v(1): v(2) = 3# * v(2)
  CDIV u(), v(), w()
  q(1) = q(1) - w(1): q(2) = q(2) - w(2)

' calculate det = q²/4 + p^3/27
  CMUL q(), q(), u(): u(1) = u(1) / 4#: u(2) = u(2) / 4#
  CMUL p(), p(), v(): CMUL p(), v(), w(): w(1) = w(1) / 27#: w(2) = w(2) / 27#
  det(1) = u(1) + w(1): det(2) = u(2) + w(2)

  CSQRT det(), sdet()  'now sdet contains sqrt(q²/4 + p^3/27)

  v(1) = -q(1) / 2# + sdet(1): v(2) = -q(2) / 2# + sdet(2)

  CUBRT v(), u(), u1(), u2()'3 cubic roots

  w(1) = -q(1) / 2# - sdet(1): w(2) = -q(2) / 2# - sdet(2)

  CUBRT w(), v(), v1(), v2()'3 cubic roots

  z(1, 1) = u(1) + v(1): z(1, 2) = u(2) + v(2)
  z(2, 1) = u(1) + v1(1): z(2, 2) = u(2) + v1(2)
  z(3, 1) = u(1) + v2(1): z(3, 2) = u(2) + v2(2)
  z(4, 1) = u1(1) + v(1): z(4, 2) = u1(2) + v(2)
  z(5, 1) = u1(1) + v1(1): z(5, 2) = u1(2) + v1(2)
  z(6, 1) = u1(1) + v2(1): z(6, 2) = u1(2) + v2(2)
  z(7, 1) = u2(1) + v(1): z(7, 2) = u2(2) + v(2)
  z(8, 1) = u2(1) + v1(1): z(8, 2) = u2(2) + v1(2)
  z(9, 1) = u2(1) + v2(1): z(9, 2) = u2(2) + v2(2)

' Note: only 3 z[i] are correct solutions

  u(1) = 3# * a(1): u(2) = 3# * a(2)
  CDIV b(), u(), v()

  FOR i = 1 TO 9
    z(i, 1) = z(i, 1) - v(1)
    z(i, 2) = z(i, 2) - v(2)
  NEXT i

  F1$ = "####.#####"

  CLS
  PRINT
  PRINT " Equation az^3 + bz^2 + cz + d = 0"
  PRINT " Solutions z:"

  FOR i = 1 TO 9
    zz(1) = z(i, 1): zz(2) = z(i, 2)
    F zz(), z1()
    IF CABS(z1()) < .001 THEN   'detect if solution is correct
      PRINT " z = (";
      PRINT USING F1$; z(i, 1);
      PRINT ",";
      PRINT USING F1$; z(i, 2);
      PRINT ")"
      PRINT " F(z) = (";
      PRINT USING F1$; z1(1);
      PRINT ",";
      PRINT USING F1$; z1(2);
      PRINT ")"
    END IF
  NEXT i
END 'of main program

'end of file croot3.bas

SUB ATAN (Numerator AS DOUBLE, denominator, Phase)
'Return a phase between -PI and +PI
  PI = 4# * ATN(1#)
  IF ABS(denominator) < 1E-15 THEN
    IF (Numerator < 0) THEN
      Phase = -PI / 2
    ELSE
      Phase = PI / 2
    END IF
  ELSE
    Phase = ATN(Numerator / denominator)
    IF denominator < 0 THEN Phase = Phase + PI
  END IF
END SUB 'ATAN

'return ABS(z)
FUNCTION CABS (z())
  CABS = SQR(z(1) * z(1) + z(2) * z(2))
END FUNCTION

' Z = Z1 / Z2
SUB CDIV (z1(), z2(), z())
DIM cc(2)
  d = z2(1) * z2(1) + z2(2) * z2(2)
  IF d < 1E-15 THEN
    PRINT " Complex Divide by zero!"
  ELSE
    cc(1) = z2(1): cc(2) = -z2(2)
    CMUL z1(), cc(), z()
    z(1) = z(1) / d: z(2) = z(2) / d
  END IF
END SUB

' z = z1 * z2
SUB CMUL (z1(), z2(), z())
  z(1) = z1(1) * z2(1) - z1(2) * z2(2)
  z(2) = z1(1) * z2(2) + z1(2) * z2(1)
END SUB

' z1 = SQRT(z)
SUB CSQRT (z(), z1())
  r = SQR(z(1) * z(1) + z(2) * z(2))
  z1(1) = SQR((r + z(1)) / 2#)
  z1(2) = SQR((r - z(1)) / 2#)
  IF z(2) < 0 THEN z1(2) = -z1(2)
END SUB

' z1, z2, z3 = z^(1/3)
SUB CUBRT (z(), z1(), z2(), z3())
  PI = 4# * ATN(1#)
  r = SQR(z(1) * z(1) + z(2) * z(2))
  r = r ^ (1# / 3#)
  ATAN z(2), z(1), t
  t = t / 3#
  z1(1) = r * COS(t)
  z1(2) = r * SIN(t)
  tt = t - (2# * PI / 3#)
  z2(1) = r * COS(tt)
  z2(2) = r * SIN(tt)
  tt = t + (2# * PI / 3#)
  z3(1) = r * COS(tt)
  z3(2) = r * SIN(tt)
END SUB

' z1 = a*z^3 + b*z2 + c*z + d
SUB F (z(), z1())
  CMUL a(), z(), u(): CMUL z(), u(), v(): CMUL z(), v(), w()
  z1(1) = w(1): z1(2) = w(2)
  CMUL b(), z(), u(): CMUL u(), z(), v()
  z1(1) = z1(1) + v(1): z1(2) = z1(2) + v(2)
  CMUL c(), z(), u()
  z1(1) = z1(1) + u(1) + d(1): z1(2) = z1(2) + u(2) + d(2)
END SUB