DECLARE SUB Croot4 (a#(), b#(), c#(), d#(), e#(), z1#(), z2#(), z3#(), z4#())
DECLARE SUB F4 (a#(), b#(), c#(), d#(), e#(), z#(), z1#())
DECLARE FUNCTION CABS# (z#())
DECLARE SUB CMUL (z1#(), z2#(), z#())
DECLARE SUB CDIV (z1#(), z2#(), z#())
DECLARE SUB CSQRT (z#(), z1#())
DECLARE SUB ATAN (Xnumerator#, denominator#, phase#)
DECLARE SUB CUBRT (z#(), z1#(), z2#(), z3#())
DECLARE SUB F (a#(), b#(), c#(), d#(), z#(), z1#())
DECLARE SUB F1 (p#(), q#(), z#(), z1#())
DECLARE SUB Equa2c (a#(), b#(), c#(), s1#(), s2#())
DECLARE SUB Croot3 (a#(), b#(), c#(), d#(), z1#(), z2#(), z3#())
'***************************************************************
'*            Solve complex equation of degree 4:              *
'*            a z^4 + b z^3 + c z^2 + d z + e = 0              *
'* ----------------------------------------------------------- *
'* SAMPLE RUN:                                                 *
'*    a = (1,0) b = (2,-10) c = (-4,1) d = (5,2) e = (3,-7.5)  *
'*                                                             *
'* Equation az^4 + bz^3 + cz^2 + d z + e = 0                   *
'* Solutions z1..z4:                                           *
'*                                                             *
'* z1 = (  -2.139849,   9.625142)                              *
'* F(z1) = (0.00000000,0.00000000)                             *
'* z2 = (   0.690192,  -0.706380)                              *
'* F(z2) = (0.00000000,-.00000000)                             *
'* z3 = (  -0.875877,   0.223149)                              *
'* F(z3) = (0.00000000,-.00000000)                             *
'* z4 = (   0.325534,   0.858088)                              *
'* F(z4) = (0.00000000,-.00000000)                             *
'*                                                             *
'* ----------------------------------------------------------- *
'* Reference: Mathematiques et statitiques - Programmes en     *
'*            BASIC. Editions du P.S.I., 1981.                 *
'*                                                             *
'*            Adapted to complex Domain By J-P Moreau, Paris.  *
'*                           (www.jpmoreau.fr)                 *
'***************************************************************
' Explanations:
' ------------
' Let us solve complex equation z^4 + a*z^3 + b*z^2 + c*z + d = 0  (1)
'
' Let us set z = z1 - a/4 and equation (1) becomes:
'
'     (z1²+2ckz1+l)(z1²-2ckz1+m) = 0       (2)
'
' where:   l + m -4ck² = q  |            q = 3a²/8
'          2 ck (m-1)  = r  | (3) with:  r = c-(ab/2)+(a^3/8)
'                 lm   = s  |            s = d-(ac/4)+(a²b/16)-(3a^4/256)
'
' ck² = zz is a solution of complex cubic equation:
'
'          zz^3 + aa zz^2 + bb z + cc = 0   (4)
'
' with:    aa = q/2, bb = (q²-4s)/16, cc = -r²/64
'
' So we solve equation (4) using Croot3 (3 roots zz1, zz2, zz3), then
' we calculate complex l and m by solving system (3), finally, we solve
' both 2nd degree equations of product (2). We obtain three sets of four
' complex roots, only one of which is valid, i.e. annulates equation (1).
'
'------------------------------------------------------------------------
'Program Test_CRoot4
DEFDBL A-H, O-Z
DEFINT I-N

'Complex numbers
DIM a(2), b(2), c(2), d(2), e(2), z1(2), z2(2), z3(2), z4(2), zz(2)

  CLS
  PRINT
  PRINT " Equation az^4 + bz^3 + cz^2 + d z + e = 0"
  PRINT " Solutions z1..z4:"

  a(1) = 1#: a(2) = 0#
  b(1) = 2#: b(2) = -10#
  c(1) = -4#: c(2) = 1#
  d(1) = 5#: d(2) = 2#
  e(1) = 3#: e(2) = -7.5

  Croot4 a(), b(), c(), d(), e(), z1(), z2(), z3(), z4()

' restore coefficients a, b, c, d, e
  a(1) = 1#: a(2) = 0#
  b(1) = 2#: b(2) = -10#
  c(1) = -4#: c(2) = 1#
  d(1) = 5#: d(2) = 2#
  e(1) = 3#: e(2) = -7.5

  FR$ = "####.######"
  FR1$ = "#.########"

' print results using formats FR$ and FR1$

  PRINT
  PRINT " z1 = ("; : PRINT USING FR$; z1(1);
  PRINT ","; : PRINT USING FR$; z1(2); : PRINT ")"
  F4 a(), b(), c(), d(), e(), z1(), zz()
  PRINT " F(z1) = ("; : PRINT USING FR1$; zz(1);
  PRINT ","; : PRINT USING FR1$; zz(2); : PRINT ")"
  PRINT " z2 = ("; : PRINT USING FR$; z2(1);
  PRINT ","; : PRINT USING FR$; z2(2); : PRINT ")"
  F4 a(), b(), c(), d(), e(), z2(), zz()
  PRINT " F(z2) = ("; : PRINT USING FR1$; zz(1);
  PRINT ","; : PRINT USING FR1$; zz(2); : PRINT ")"
  PRINT " z3 = ("; : PRINT USING FR$; z3(1);
  PRINT ","; : PRINT USING FR$; z3(2); : PRINT ")"
  F4 a(), b(), c(), d(), e(), z3(), zz()
  PRINT " F(z3) = ("; : PRINT USING FR1$; zz(1);
  PRINT ","; : PRINT USING FR1$; zz(2); : PRINT ")"
  PRINT " z4 = ("; : PRINT USING FR$; z4(1);
  PRINT ","; : PRINT USING FR$; z4(2); : PRINT ")"
  F4 a(), b(), c(), d(), e(), z4(), zz()
  PRINT " F(z4) = ("; : PRINT USING FR1$; zz(1);
  PRINT ","; : PRINT USING FR1$; zz(2); : PRINT ")"
  PRINT

END 'of main program

SUB ATAN (Xnumerator, denominator, phase)
'Returna a phase between -PI and +PI
  PI = 4# * ATN(1#)
  IF ABS(denominator) < .000000000000001# THEN
    IF Xnumerator < 0# THEN
      phase = -PI / 2#
    ELSE
      phase = PI / 2#
    END IF
    EXIT SUB
  ELSE
    phase = ATN(Xnumerator / 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 c(2)
  d = z2(1) * z2(1) + z2(2) * z2(2)
  IF d < .000000000000001# THEN
    PRINT " Complex Divide by zero #"
  ELSE
    c(1) = z2(1): c(2) = -z2(2)
    CMUL z1(), c(), 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

'solve complex equation a*z^3+b*z^2+c*z+d=0 and return
' 3 solutions z1, z2, z3.
SUB Croot3 (a(), b(), c(), d(), z1(), z2(), z3())
' -----------------------------------------------------
' 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))
'
' -----------------------------------------------------------------------------------
DIM p(2), q(2), det(2), sdet(2), u(2), u1(2), u2(2), v(2), v1(2), v2(2), w(2), zz(2)
DIM z(9, 2)

' calculate p
  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
  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)

  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^2/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.

  index = 1
  FOR i = 1 TO 9
    u(1) = z(i, 1): u(2) = z(i, 2)
    F1 p(), q(), u(), zz()
    IF CABS(zz()) < .001 THEN   'detect if solution is correct
      IF index = 1 THEN
        z1(1) = z(i, 1): z1(2) = z(i, 2): index = index + 1
      ELSEIF index = 2 THEN
        z2(1) = z(i, 1): z2(2) = z(i, 2): index = index + 1
      ELSE
        z3(1) = z(i, 1): z3(2) = z(i, 2)
      END IF
    END IF
  NEXT i

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

  z1(1) = z1(1) - v(1): z1(2) = z1(2) - v(2)
  z2(1) = z2(1) - v(1): z2(2) = z2(2) - v(2)
  z3(1) = z3(1) - v(1): z3(2) = z3(2) - v(2)

END SUB 'croot3

'-----------------------------------------------------------
'solve complex equation a*z^4+b*z^3+c*z^2+d*z+e=0 and return
' 4 solutions z1,z2,z3,z4
'-----------------------------------------------------------
SUB Croot4 (a(), b(), c(), d(), e(), z1(), z2(), z3(), z4())
'Labels: 100,200,300,400
DIM a1(2), aa(2), b1(2), bb(2), c1(2), cc(2), q(2), r(2), s(2), t(2), u(2), v(2), w(2)
DIM ck(2), xl(2), xm(2)
DIM z(3, 2), zz(12, 2)

  t(1) = a(1): t(2) = a(2)
  CDIV b(), t(), a(): CDIV c(), t(), b(): CDIV d(), t(), c(): CDIV e(), t(), d()
' Now equation is: z^4+a*z^3+b*z^2+c*z+d = 0

' q = b - (3 * a * a / 8)
  CMUL a(), a(), u(): u(1) = u(1) * 3# / 8#: u(2) = u(2) * 3# / 8#
  q(1) = b(1) - u(1): q(2) = b(2) - u(2)

' r = c - (a * b / 2) + (a * a * a / 8)
  CMUL a(), b(), u(): r(1) = -u(1) / 2#: r(2) = -u(2) / 2#
  CMUL a(), a(), u(): CMUL a(), u(), v(): v(1) = v(1) / 8#: v(2) = v(2) / 8#
  r(1) = r(1) + c(1) + v(1): r(2) = r(2) + c(2) + v(2)

' s = d - (a * c / 4) + (a * a * b / 16) - (3 * a * a * a * a / 256
  CMUL a(), c(), u(): u(1) = u(1) / 4#: u(2) = u(2) / 4#
  CMUL a(), a(), v(): CMUL b(), v(), w(): v(1) = w(1) / 16#: v(2) = w(2) / 16#
  CMUL a(), a(), t(): CMUL t(), t(), w(): w(1) = w(1) * 3# / 256#: w(2) = w(2) * 3# / 256#
  s(1) = d(1) - u(1) + v(1) - w(1): s(2) = d(2) - u(2) + v(2) - w(2)

' Défine coefficients of cubic equation z^3+aa*z^2+bb*z+cc=0 (1)
' aa = q / 2
  aa(1) = q(1) / 2#: aa(2) = q(2) / 2#
' bb = (q * q - 4 * s) / 16
  CMUL q(), q(), u()
  bb(1) = (u(1) - 4# * s(1)) / 16#: bb(2) = (u(2) - 4# * s(2)) / 16#
' cc = -(r * r / 64
  CMUL r(), r(), u()
  cc(1) = -u(1) / 64#: cc(2) = -u(2) / 64#

' Calculate complex roots of equation (1)
  a1(1) = 1#: a1(2) = 0#
  IF CABS(r()) > .000000000000001# OR CABS(bb()) > .000000000000001# THEN
    GOTO 100
  END IF

' Particular case when equation (1) is of 2nd order
  b1(1) = aa(1): b1(2) = aa(2): c1(1) = bb(1): c1(2) = bb(2)
  Equa2c a1(), b1(), c1(), u(), v()
  z(1, 1) = u(1): z(1, 2) = u(2)
  z(2, 1) = v(1): z(2, 2) = v(2)
  z(3, 1) = 0#: z(3, 2) = 0#
  GOTO 200

100 Croot3 a1(), aa(), bb(), cc(), u(), v(), w()
  z(1, 1) = u(1): z(1, 2) = u(2)
  z(2, 1) = v(1): z(2, 2) = v(2)
  z(3, 1) = w(1): z(3, 2) = w(2)

200 FOR i = 1 TO 3
    u(1) = z(i, 1): u(2) = z(i, 2)
    CSQRT u(), ck()

'  Calculate xl and xm if ck=0
    IF CABS(ck()) = 0# THEN
      'r = SQR(q * q - 4 * s)
      CMUL q(), q(), u(): v(1) = 4# * s(1): v(2) = 4# * s(2)
      r(1) = u(1) - v(1): r(2) = u(2) - v(2)
      GOTO 300
    END IF
    'q = q + (4 * z(i))
    u(1) = 4# * z(i, 1): u(2) = 4# * z(i, 2)
    q(1) = q(1) + u(1): q(2) = q(2) + u(2)
    'r = r / (2 * ck)
    u(1) = 2# * ck(1): u(2) = 2# * ck(2)
    CDIV r(), u(), v()
    r(1) = v(1): r(2) = v(2)
300 'xl = (q - r) / 2, xm = (q + r) / 2
    xl(1) = (q(1) - r(1)) / 2#: xl(2) = (q(2) - r(2)) / 2#
    xm(1) = (q(1) + r(1)) / 2#: xm(2) = (q(2) + r(2)) / 2#

    'Solving two equations of degree 2
    b1(1) = 2# * ck(1): b1(2) = 2# * ck(2)
    '1st equation
    Equa2c a1(), b1(), xl(), u(), v()

    b1(1) = -b1(1): b1(2) = -b1(2)
    '2nd equation
    Equa2c a1(), b1(), xm(), w(), t()

    'Transferring solutions in zz(i)
    IF i = 1 THEN
      zz(1, 1) = u(1): zz(1, 2) = u(2)
      zz(2, 1) = v(1): zz(2, 2) = v(2)
      zz(3, 1) = w(1): zz(3, 2) = w(2)
      zz(4, 1) = t(1): zz(4, 2) = t(2)
    ELSEIF i = 2 THEN
      zz(5, 1) = u(1): zz(5, 2) = u(2)
      zz(6, 1) = v(1): zz(6, 2) = v(2)
      zz(7, 1) = w(1): zz(7, 2) = w(2)
      zz(8, 1) = t(1): zz(8, 2) = t(2)
    ELSE
      zz(9, 1) = u(1): zz(9, 2) = u(2)
      zz(10, 1) = v(1): zz(10, 2) = v(2)
      zz(11, 1) = w(1): zz(11, 2) = w(2)
      zz(12, 1) = t(1): zz(12, 2) = t(2)
    END IF  'Note: only 4 solutions are correct
  NEXT i

  u(1) = a(1) / 4#: u(2) = a(2) / 4#
  FOR i = 1 TO 12'shift zz by -a/4
    zz(i, 1) = zz(i, 1) - u(1): zz(i, 2) = zz(i, 2) - u(2)
  NEXT i
  'select 4 correct solutions
  i = 1
400 u(1) = zz(i, 1): u(2) = zz(i, 2)
  F a(), b(), c(), d(), u(), t()
  IF CABS(t()) < .0001 THEN
    IF i = 1 THEN
      z1(1) = zz(1, 1): z1(2) = zz(1, 2)
      z2(1) = zz(2, 1): z2(2) = zz(2, 2)
      z3(1) = zz(3, 1): z3(2) = zz(3, 2)
      z4(1) = zz(4, 1): z4(2) = zz(4, 2)
    ELSEIF i = 5 THEN
      z1(1) = zz(5, 1): z1(2) = zz(5, 2)
      z2(1) = zz(6, 1): z2(2) = zz(6, 2)
      z3(1) = zz(7, 1): z3(2) = zz(7, 2)
      z4(1) = zz(8, 1): z4(2) = zz(8, 2)
    ELSE
      z1(1) = zz(9, 1): z1(2) = zz(9, 2)
      z2(1) = zz(10, 1): z2(2) = zz(10, 2)
      z3(1) = zz(11, 1): z3(2) = zz(11, 2)
      z4(1) = zz(12, 1): z4(2) = zz(12, 2)
    END IF
  END IF
  i = i + 4
  IF i <= 9 THEN GOTO 400

  'now complex roots are in z1,z2,z3,z4.

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 = 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

'return complex roots of equation a*z^2+b*z+c = 0
SUB Equa2c (a(), b(), c(), s1(), s2())
DIM u(2), v(2), w(2)
  u(1) = b(1) * b(1) - b(2) * b(2) - 4 * a(1) * c(1) + 4 * a(2) * c(2)
  u(2) = 2 * b(1) * b(2) - 4 * a(1) * c(2) - 4 * a(2) * c(1)
  r = SQR(u(1) * u(1) + u(2) * u(2))
  v(1) = SQR((r + u(1)) / 2#)
  v(2) = SQR((r - u(1)) / 2#)
  IF u(2) < 0 THEN v(2) = -v(2)
  w(1) = (-b(1) - v(1)) / 2#
  w(2) = (-b(2) - v(2)) / 2#
  u(1) = (-b(1) + v(1)) / 2#
  u(2) = (-b(2) + v(2)) / 2#
  r = a(1) * a(1) + a(2) * a(2)
  s1(1) = (a(1) * w(1) + a(2) * w(2)) / r
  s1(2) = (a(1) * w(2) - a(2) * w(1)) / r
  s2(1) = (a(1) * u(1) + a(2) * u(2)) / r
  s2(2) = (a(1) * u(2) - a(2) * u(1)) / r
END SUB

'return z^4+a*z^3+b*z^2+c^z+d (called by Croot4)
SUB F (a(), b(), c(), d(), z(), z1())
DIM u(2), v(2), w(2)
  CMUL z(), z(), u(): CMUL u(), u(), v(): z1(1) = v(1): z1(2) = v(2)
  CMUL u(), z(), v(): CMUL a(), v(), w()
  z1(1) = z1(1) + w(1): z1(2) = z1(2) + w(2)
  CMUL b(), u(), 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

' z1 = z^3 + p*z + q (called by Croot3)
SUB F1 (p(), q(), z(), z1())
DIM u(2), v(2), w(2)
  CMUL z(), z(), u(): CMUL z(), u(), v()
  CMUL p(), z(), w()
  z1(1) = v(1) + w(1) + q(1)
  z1(2) = v(2) + w(2) + q(2)
END SUB

'return a*z^4+b*z^3+c*z^2+d^z+e (called by main)
SUB F4 (a(), b(), c(), d(), e(), z(), z1())
DIM u(2), v(2), w(2)
  CMUL z(), z(), u(): CMUL u(), u(), v(): CMUL a(), v(), z1()
  CMUL u(), z(), v(): CMUL b(), v(), w()
  z1(1) = z1(1) + w(1): z1(2) = z1(2) + w(2)
  CMUL c(), u(), v()
  z1(1) = z1(1) + v(1): z1(2) = z1(2) + v(2)
  CMUL d(), z(), u()
  z1(1) = z1(1) + u(1) + e(1): z1(2) = z1(2) + u(2) + e(2)
END SUB