' ******************************************************************
' This program solves a polynomial equation
'
' Analytical solutions are used up to degree 4, then the polynomial
' is solved by the method of the companion matrix.
'
' The example polynomial is:
'
' 720 - 1764 * X + 1624 * X^2 - 735 * X^3 +
'       175 * X^4 - 21 * X^5 + X^6
'
' The roots are: X = 1, 2 ... 6
' ******************************************************************

#INCLUDE "fbmath.bi"

' ------------------------------------------------------------------
' Define polynomial here, in the form:
' Coef(0) + Coef(1) * X + Coef(2) * X^2 + ... + Coef(Deg) * X^Deg
' ------------------------------------------------------------------

CONST Deg = 6

DIM AS DOUBLE Coef(0 TO Deg) = {720, -1764, 1624, -735, 175, -21, 1}

' ------------------------------------------------------------------

DIM AS Complex Z(1 TO Deg)  ' Roots of polynomial

DIM AS INTEGER Nr, Nc       ' Number of real and complex roots
DIM AS INTEGER I, J         ' Loop variables
DIM AS DOUBLE  A            ' Auxiliary variable

PRINT : PRINT "Polynomial:" : PRINT

FOR I = 0 TO Deg
  IF Coef(I) <> 0 THEN
    IF Coef(I) > 0 THEN PRINT " + ";
    IF Coef(I) < 0 THEN PRINT " - ";
    PRINT USING "####.######"; ABS(Coef(I));
    IF I > 0 THEN PRINT " X";
    IF I > 1 THEN PRINT "^"; STR$(I);
    PRINT
  END IF
NEXT I

PRINT : PRINT

' Solve polynomial. Nr is the number of real roots

SELECT CASE Deg
  CASE 1
    Nr = RootPol1(Coef(0), Coef(1), Z(1).X)
  CASE 2
    Nr = RootPol2(Coef(), Z())
  CASE 3
    Nr = RootPol3(Coef(), Z())
  CASE 4
    Nr = RootPol4(Coef(), Z())
  CASE ELSE
    Nr = RootPol(Coef(), Deg, Z())
END SELECT

' Case when an error occurs

IF Nr < 0 THEN
  PRINT "Error during root evaluation !"
  STOP
END IF

' Set the small imaginary parts to zero (optional)

Nr = SetRealRoots(Deg, Z(), 1D-6)

' Sort roots: first real roots, in ascending order,
' then complex roots (unordered)

SortRoots Deg, Z()

' Print real roots

IF Nr > 0 THEN
  PRINT Nr; " real root(s):"
  PRINT

  FOR I = 1 TO Nr
    PRINT USING "Z(#) = ####.######"; I; Z(I).X
  NEXT I

  PRINT
END IF

' Print complex roots

Nc = Deg - Nr

IF Nc > 0 THEN
  PRINT Nc; " complex roots:"
  PRINT

  FOR I = 1 TO Nc
    J = I + Nr
    PRINT USING "Z(#) = ####.######"; J; Z(J).X;
    IF Z(J).Y > 0 THEN PRINT " + ";  ELSE PRINT " - ";
    PRINT USING "####.######"; ABS(Z(J).Y); " * i"
  NEXT I
END IF

SLEEP
