' ******************************************************************
' Roots of a polynomial from the companion matrix
' ******************************************************************

' ------------------------------------------------------------------
' Type definition
' ------------------------------------------------------------------

TYPE Complex
  X AS DOUBLE
  Y AS DOUBLE
END TYPE

' ------------------------------------------------------------------
' Global variable
' ------------------------------------------------------------------

COMMON SHARED ErrCode AS INTEGER
' Error code from the last function evaluation

' ------------------------------------------------------------------
' External functions
' ------------------------------------------------------------------

DECLARE SUB Balance(      A()     AS DOUBLE,_ 
                    BYREF Lo      AS INTEGER,_ 
                    BYREF Hi      AS INTEGER,_ 
                          Scale() AS DOUBLE)
' Balances a matrix and tries to isolate eigenvalues

DECLARE SUB Hqr(H()      AS DOUBLE,_ 
                Lo       AS INTEGER,_ 
                Hi       AS INTEGER,_ 
                Lambda() AS Complex) 
' Eigenvalues of a real upper Hessenberg matrix by QR method

' ******************************************************************

FUNCTION RootPol(Coef() AS DOUBLE,_ 
                 Deg    AS INTEGER,_ 
                 Z()    AS Complex) AS INTEGER
' ------------------------------------------------------------------
' Solves the polynomial equation:
' Coef(0) + Coef(1) * Z + Coef(2) * Z^2 + ... 
'                       + Coef(Deg) * Z^Deg = 0
' ------------------------------------------------------------------

  DIM AS INTEGER Lo, Hi  ' Used by Balance
  DIM AS INTEGER I, J    ' Loop variables
  DIM AS INTEGER Nr      ' Number of real roots

  DIM AS DOUBLE A(1 TO Deg, 1 TO Deg)  ' Companion matrix
  DIM AS DOUBLE Scale(1 TO Deg)        ' Used by Balance

  ' Set up the companion matrix
  FOR J = 1 TO Deg
    A(1, J) = -Coef(Deg - J) / Coef(Deg)
  NEXT J

  FOR I = 2 TO Deg
    FOR J = 1 TO Deg
      IF I - 1 = J THEN A(I, J) = 1 ELSE A(I, J) = 0
    NEXT J
  NEXT I

  ' The roots of the polynomial are the
  ' eigenvalues of the companion matrix
  Balance A(), Lo, Hi, Scale()
  Hqr A(), Lo, Hi, Z()

  IF ErrCode <> 0 THEN RootPol = ErrCode: EXIT FUNCTION

  ' Count real roots
  Nr = 0
  FOR I = 1 TO Deg
    IF Z(I).Y = 0 THEN Nr = Nr + 1
  NEXT I

  RootPol = Nr
END FUNCTION

