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