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