'********************************************************
'*             Equations of degree 2 to 4               *
'* ---------------------------------------------------- *
'* This program calculates the real or complex roots of *
'* algebraic equations of degree 2, 3 and 4.            *
'* ---------------------------------------------------- *
'* Reference: Mathematiques et statitiques - Programmes *
'*            en BASIC. Editions du P.S.I., 1981        *
'*            [BIBLI 13].                               *
'* ---------------------------------------------------- *
'* SAMPLE RUN:                                          *
'*                                                      *
'* ROOTS OF EQUATIONS OF DEGREE 2, 3 OU 4               *
'* Example: X^4+4X^2-3X+3 = 0                           *
'*                                                      *
'* INSTRUCTIONS:                                        *
'*                                                      *
'* 1. Input degree of equation                          *
'* 2. Input coefficients beginning with highest exponent*
'* 3. Solutions are given under the form:               *
'*     (real part) + I (imaginary part)                 *
'*                                                      *
'* Degree of equation (2 to 4): 4                       *
'*                                                      *
'* Coefficient of X^4 = 1                               *
'* Coefficient of X^3 = 0                               *
'* Coefficient of X^2 = 4                               *
'* Coefficient of X^1 = -3                              *
'* Coefficient of X^0 = 3                               *
'*                                                      *
'* Root 1 = 0.449702650 + I 0.731070280                 *
'* Root 2 = 0.449702650 - I 0.731070280                 *
'* Root 1 = -0.449702650 + I 1.967231870                *
'* Root 2 = -0.449702650 - I 1.967231870                *
'*                                                      *
'********************************************************
CLS
DIM a(4, 4), r(4)
F$="##.#########"
PRINT " ROOTS OF EQUATIONS OF DEGREE 2, 3 OU 4"
PRINT " Example: X^4+4X^2-3X+3 = 0"
PRINT
PRINT " INSTRUCTIONS:"
PRINT
PRINT " 1. Input degree of equation"
PRINT " 2. Input coefficients beginning by highest exponent"
PRINT " 3. Solutions are given under the form:"
PRINT "    (real part) + I (imaginary part)"
PRINT
INPUT " Degree of equation (2 to 4): ", n
IF n < 2 THEN
  n = 2
END IF
IF n > 4 THEN
  n = 4
END IF
PRINT
FOR i = n TO 0 STEP -1
PRINT " Coefficient of X^"; i; : INPUT " = ", a(n, i)
NEXT i
PRINT : PRINT
' calling root search subroutine
GOSUB 1000
' printing results
ON n - 1 GOTO 400, 500, 600
' case n=2
400 IF im = 0 THEN
  PRINT " Root 1 = "; : print using F$; r(1)
  PRINT " Root 2 = "; : print using F$; r(2)
ELSE
  PRINT " Root 1 = "; : print using F$; r(1);
  PRINT "  +I  "; : print using F$; im
  PRINT " Root 2 = "; : print using F$; r(2);
  PRINT "  -I  "; : print using F$; im
END IF
END
' case n=3
500 PRINT " Root 1 = "; : print using F$; r(3)
IF im = 0 THEN
  PRINT " Root 2 = "; : print using F$; r(1)
  PRINT " Root 3 = "; : print using F$; r(2)
ELSE
  PRINT " Root 2 = "; : print using F$; r(1);
  PRINT "  +I  "; : print using F$; im
  PRINT " Root 3 = "; : print using F$; r(2);
  PRINT "  -I  "; : print using F$; im
END IF
END
' cas n=4
600 IF im = 0 THEN
  PRINT " Root 1 = "; : print using F$; r(1)
  PRINT " Root 2 = "; : print using F$; r(2)
ELSE
  PRINT " Root 1 = "; : print using F$; r(1);
  PRINT "  +I  "; : print using F$; im
  PRINT " Root 2 = "; : print using F$; r(2);
  PRINT "  -I  "; : print using F$; im
END iF

IF ii = 0 THEN
  PRINT " Root 3 = "; : print using F$; r(3)
  PRINT " Root 4 = "; : print using F$; r(4)
ELSE
  PRINT " Root 3 = "; : print using F$; r(3);
  PRINT "  +I  "; : print using F$; ii
  PRINT " Root 4 = "; : print using F$; r(4);
  PRINT "  -I  "; : print using F$; ii
END IF
END

'Roots search subroutine
1000 'Normalization of coefficients
FOR i = 0 TO n - 1
  a(n, i) = a(n, i) / a(n, n)
NEXT i
' Branch according to n
ON n - 1 GOTO 1400, 1500, 1600
' cas n=2
1400 GOSUB 3000
RETURN
' cas n=3
1500 GOSUB 2000
RETURN
' cas n=4
1600 a = a(4, 3)
b = a(4, 2)
c = a(4, 1)
d = a(4, 0)
q = b - (3 * a * a / 8)
r = c - (a * b / 2) + (a * a * a / 8)
s = d - (a * c / 4) + (a * a * b / 16) - (3 * a * a * a * a / 256)
' Definition of coefficients of cubic equation
a(3, 2) = q / 2
a(3, 1) = (q * q - 4 * s) / 16
a(3, 0) = -(r * r / 64)
' Calculates a real root of this equation
IF (r <> 0) OR (a(3, 1) >= 0) THEN GOTO 1650
' Case when this equation is reduced to 2nd degree
a(2, 1) = a(3, 2): a(2, 0) = a(3, 1)
GOSUB 3000
' One takes positive root
r(3) = r(1)
GOTO 1700
' Calling subroutine 2000 with sw=-1 to calculate only one root
1650 sw = -1
GOSUB 2000
' real root of above equation
1700 k = SQR(r(3))
' Calculates L and M when k=0
IF k = 0 THEN
  r = SQR(q * q - 4 * s): GOTO 1750
END IF
q = q + (4 * r(3))
r = r / (2 * k)
1750 l = (q - r) / 2: m = (q + r) / 2
' Solving two 2nd degree equations
a(2, 1) = 2 * k: a(2, 0) = l
' 1st equation
GOSUB 3000
' Solutions  transferred in r(3), r(4), ii
r(3) = r(1) - (a(4, 3) / 4): r(4) = r(2) - (a(4, 3) / 4): ii = im
a(2, 1) = -2 * k: a(2, 0) = m
' 2nd equation
GOSUB 3000
r(2) = r(2) - (a(4, 3) / 4): r(1) = r(1) - (a(4, 3) / 4)
RETURN

'********************************************
'* This subroutine calculates the roots of  *
'* X^3 + A(3,2)*X^2 + A(3,1)*X + A(3,0) = 0 *
'********************************************
' Case when a root equals zero 
2000 IF a(3, 0) = 0 THEN
  xc = 0: GOTO 2500
END IF
' Searching highest coefficient in absolute value
am = ABS(a(3, 0))
FOR i = 1 TO 3
  tt = ABS(a(3, i))
  IF (am < tt) THEN am = tt
NEXT i
' Define interval where a real root exists 
' according to sign of A(3,0)
IF a(3, 0) > 0 THEN
  a = -am - 1
  b = 0
  GOTO 2100
END IF
a = 0: b = am + 1

2100
' Searching for xc = real root in interval (a,b)
' (by Bisection method)
DEF fnf (x) = x * x * x + a(3, 2) * x * x + a(3, 1) * x + a(3, 0)
' Define segment (xa,xb) including the root
xa = a: xb = b: y1 = fnf(a): y2 = fnf(b)
' Desired precision 
er = .000001
IF y1 = 0 THEN
  xc = a: GOTO 2500
END IF
IF y2 = 0 THEN
  xc = b: GOTO 2500
END IF
' xc = middle of segment
2200 xc = (xa + xb) / 2: te = fnf(xc)
IF te = 0 THEN GOTO 2500
IF (xb - xa) < er THEN GOTO 2500
IF (fnf(xa) * te) > 0 THEN
  xa = xc: GOTO 2200
END IF
xb = xc: GOTO 2200
' r(3) is a real root
2500 r(3) = xc
IF sw = -1 THEN RETURN
' Calculates the roots of remaining quadratic equation
' Define the equation coefficients
a(2, 1) = a(3, 2) + xc
a(2, 0) = a(3, 1) + a(3, 2) * xc + xc * xc
GOSUB 3000
RETURN
'End subroutine 2000

'*******************************************
'* This Subroutine calculates the roots of *
'*    X^2 + A(2,1)*X + A(2,0) = 0          *
'*    W = Determinant of equation          *
'*******************************************
3000 w = a(2, 1) * a(2, 1) - 4 * a(2, 0)
q1 = -a(2, 1) / 2
ON (SGN(w) + 2) GOTO 3200, 3300, 3400
' Cas W < 0
3200 q2 = SQR(-w) / 2
im = q2
r(1) = q1: r(2) = q1
RETURN
' Cas W = 0
3300 r(1) = q1: r(2) = q1: im = 0
' Cas W > 0
3400 q2 = SQR(w) / 2: im = 0
r(1) = q1 + q2: r(2) = q1 - q2
RETURN

'End of file root4.bas