'*********************************************************
'* Calculate the coefficients of the characteristic poly-*
'* nomial P(l)=A-l*I of a square complex matrix A(i,j)   *
'* ----------------------------------------------------- *
'* Ref.: "Algèbre, Algorithmes et programmes en Pascal   *
'*        by Jean-Louis Jardrin, DUNOD Paris, 1988" [10].*
'* ----------------------------------------------------- *
'* SAMPLE RUN:                                           *
'* (Find the characteristic polynomial of matrix:        *
'*           1 -1  3                                     *
'*      A = -1  i -1-2i                                  *
'*           i  1 -2+i   )                               *
'*                                                       *
'* Input matrix size (max 20): 3                         *
'*                                                       *
'* Line 1                                                *
'*  Column 1                                             *
'*   Real  part: 1                                       *
'*   Imag. part: 0                                       *
'*  Column 2                                             *
'*   Real  part: -1                                      *
'*   Imag. part: 0                                       *
'*  Column 3                                             *
'*   Real  part: 3                                       *
'*   Imag. part: 0                                       *
'* Line 2                                                *
'*  Column 1                                             *
'*   Real  part: -1                                      *
'*   Imag. part: 0                                       *
'*  Column 2                                             *
'*   Real  part: 0                                       *
'*   Imag. part: 1                                       *
'*  Column 3                                             *
'*   Real  part: -1                                      *
'*   Imag. part: -2                                      *
'* Line 3                                                *
'*  Column 1                                             *
'*   Real  part: 0                                       *
'*   Imag. part: 1                                       *
'*  Column 2                                             *
'*   Real  part: 1                                       *
'*   Imag. part: 0                                       *
'*  Column 3                                             *
'*   Real  part: -2                                      *
'*   Imag. part: 1                                       *
'*                                                       *
'* Characteristic polynomial P(l):                       *
'* Coefficient for degree 3:                             *
'* -1.000000 + 0.000000 i                                *
'* Coefficient for degree 2:                             *
'* -1.000000 + 2.000000 i                                *
'* Coefficient for degree 1:                             *
'*  3.000000 + 1.000000 i                                *
'* Coefficient for degree 0:                             *
'*  0.000000 + 0.000000 i                                *
'*                                                       *
'* (The characteristic polynomial of matrix A is:        *
'*   P(l)=-l^3+(-1+2i)l^2+(3+i)l                         *
'*                                                       *
'*                 Basic version by J-P Moreau, Paris    *
'*                          (www.jpmoreau.fr)            *
'*********************************************************
'Program CARPOL1
'Note: a complex number is represented by a couple of two
       real numbers (Ex.:z0 by z01 and z02).
DEFDBL A-H, O-Z
DEFINT I-N

NMAX = 20
NMAXP = 21

DIM A1(NMAX, NMAX), A2(NMAX, NMAX)
DIM P1(NMAXP), P2(NMAXP)

  GOSUB 1000  'Read data
  GOSUB 2000  'Call routine to calculate coefficients

' print results
  F$="###.######"
  PRINT
  PRINT " Characteristic polynomial P(l):"
  FOR k = 1 TO n + 1
    PRINT " Coefficient for degree "; n - k + 1; ": "
    PRINT USING F$; P1(k);
    IF P2(k) >= 0 THEN
      PRINT " +";
    ELSE
      PRINT " -";
    END IF
    PRINT USING F$; ABS(P2(k)); : PRINT " i"
  NEXT
END

1000 'Read data
  CLS : PRINT
  PRINT " Input matrix size (max "; NMAX; "): ";
  INPUT "", n
  FOR i = 1 TO n
    PRINT : PRINT " Line "; i
    FOR j = 1 TO n
      PRINT "  Column "; j
      PRINT "   real  part: "; : INPUT "", A1(i, j)
      PRINT "   imag. part: "; : INPUT "", A2(i, j)
    NEXT j
  NEXT i
RETURN

1500 'calculate complex trace of matrix C
't01,t02:Complex
  c01 = 0#: c02 = 0#
  FOR i = 1 TO n
    t01 = c01: t02 = c02
    'CADD(t0,C(i,i),c0)
    c01 = t01 + C1(i, i)
    c02 = t02 + C2(i, i)
  NEXT
RETURN

'******************************************************
'* The procedure PCCS calculates the complex coeffi-  *
'* cients of the characteristic polynomial P(l)=A-l*I *
'* of a real square tridiagonal matrix A(i,j).        *
'*                                                    *
'* Note: the roots of P(l) are the eigenvalues of     *
'*       matrix A(i,j).                               *
'* -------------------------------------------------- *
'* INPUTS:                                            *
'*          n: zize of matrix (integer)               *
'*          A: complex square matrix (n,n)            *
'*             (see main program).                    *
'* OUTPUT:                                            *
'*          P: vector of complex coefficients of P(l) *
'*                                                    *
'******************************************************
2000 'Procedure PCCS
' c0,c1,s,t0,z0,z1: Complex
DIM B(NMAX, NMAX), C(NMAX, NMAX)
  z01 = 0#: z02 = 0#
  z11 = 1#: z12 = 0#
  IF n / 2 <> INT(n / 2) THEN
    'CPRO(-1.0,z1,P^[1])
    P1(1) = -1# * z11: P2(1) = -1# * z12
  ELSE
    P1(1) = z11: P2(1) = z12
  END IF
  FOR l = 1 TO n
    IF l = 1 THEN
      FOR i = 1 TO n
        FOR j = 1 TO n
          C1(i, j) = A1(i, j)
          C2(i, j) = A2(i, j)
        NEXT j
      NEXT i
    ELSE
      FOR i = 1 TO n
        FOR j = 1 TO n
          s1 = z01: s2 = z02
          FOR k = 1 TO n
            c01 = s1: c02 = s2
            'CMUL(B(i,k),A(k,j),c1)
            c11 = B1(i, k) * A1(k, j) - B2(i, k) * A2(k, j)
            c12 = B1(i, k) * A2(k, j) + B2(i, k) * A1(k, j)
            'CADD(c0,c1,s)
            s1 = c01 + c11: s2 = c02 + c12
          NEXT k
          C1(i, j) = s1: C2(i, j) = s2
        NEXT j
      NEXT i
    END IF
    GOSUB 1500  'TRC(n,C,c0)
    'CPRO(1.0/l,c0,t0)
    t01 = 1# / l * c01: t02 = 1# / l * c02
    'CMUL(t0,P(1),c0)
    c01 = t01 * P1(1) - t02 * P2(1)
    c02 = t01 * P2(1) + t02 * P1(1)
    'CPRO(-1.0,c0,P(l+1))
    P1(l + 1) = -1# * c01: P2(l + 1) = -1# * c02
    IF l < n THEN
      FOR i = 1 TO n
        FOR j = 1 TO n
          IF j = i THEN
            'CDIF(C(i,j),t0,B(i,j))
            B1(i, j) = C1(i, j) - t01
            B2(i, j) = C2(i, j) - t02
          ELSE
            B1(i, j) = C1(i, j): B2(i, j) = C2(i, j)
          END IF
        NEXT j
      NEXT i
    END IF
  NEXT l
RETURN

'end of file carpol1.bas