'*********************************************************
'*   Calculate the coefficients of the characteristic    *
'*    polynomial P(l)=A-l*I of a real matrix A(i,j)      *
'*                 by Souriau's method                   *
'* ----------------------------------------------------- *
'* Ref.: "Algebre, Algorithmes et programmes en Pascal   *
'*        By Jean-Louis Jardrin, DUNOD Paris, 1988" [10].*
'* ----------------------------------------------------- *
'* SAMPLE RUN:                                           *
'* (Find the caracteristic polynomial of matrix:         *
'*           1  0  0  0  0  1                            *
'*           1  1  0  0  0 -1                            *
'*      A = -1  1  1  0  0  1                            *
'*           1 -1  1  1  0 -1                            *
'*          -1  1 -1  1  1  1                            *
'*           1 -1  1 -1  1 -1  )                         *
'*                                                       *
'* Input matrix size (max 100): 6                        *
'*  Line 1:                                              *
'*   Element 1: 1                                        *
'*   Element 2: 0                                        *
'*   Element 3: 0                                        *
'*   Element 4: 0                                        *
'*   Element 5: 0                                        *
'*   Element 6: 1                                        *
'*  Line 2:                                              *
'*   Element 1: 1                                        *
'*   Element 2: 1                                        *
'*   Element 3: 0                                        *
'*   Element 4: 0                                        *
'*   Element 5: 0                                        *
'*   Element 6: -1                                       *
'*  Line 3:                                              *
'*   Element 1: -1                                       *
'*   Element 2: 1                                        *
'*   Element 3: 1                                        *
'*   Element 4: 0                                        *
'*   Element 5: 0                                        *
'*   Element 6: 1                                        *
'*  Line 4:                                              *
'*   Element 1: 1                                        *
'*   Element 2: -1                                       *
'*   Element 3: 1                                        *
'*   Element 4: 1                                        *
'*   Element 5: 0                                        *
'*   Element 6: -1                                       *
'*  Line 5:                                              *
'*   Element 1: -1                                       *
'*   Element 2: 1                                        *
'*   Element 3: -1                                       *
'*   Element 4: 1                                        *
'*   Element 5: 1                                        *
'*   Element 6: 1                                        *
'*  Line 6:                                              *
'*   Element 1: 1                                        *
'*   Element 2: -1                                       *
'*   Element 3: 1                                        *
'*   Element 4: -1                                       *
'*   Element 5: 1                                        *
'*   Element 6: -1                                       *
'*                                                       *
'* Chracteristic polynomial P(l):                        *
'* Coefficient for degree 6:    1.000000                 *
'* Coefficient for degree 5:   -4.000000                 *
'* Coefficient for degree 4:    0.000000                 *
'* Coefficient for degree 3:   30.000000                 *
'* Coefficient for degree 2:  -75.000000                 *
'* Coefficient for degree 1:   79.000000                 *
'* Coefficient for degree 0:  -32.000000                 *
'*                                                       *
'* (The characteristic polynomial of matrix A is:        *
'*   P(l)=l^6-4*l^5+30*l^3-75*l^2+79*l-32 )              *
'*                                                       *
'*                 Basic version by J-P Moreau, Paris    *
'*                           (www.jpmoreau.fr)           *
'*********************************************************
'Program CARPOL2
defdbl a-h,o-z
defint i-n

NMAX  = 20
NMAXP = 21

' k,n: Integer
dim  A(NMAX,NMAX)
dim  P(NMAXP)

  cls
  gosub 1000  'Read data
' call routine to calculate coefficients
  gosub 2000  'PCMS(n,A,P)
' print results
  print
  print " Characteristic polynomial P(l):"
  For k=1 to n+1
    print " Coefficient for degree "; n-k+1; ": ";
    print using "####.######"; P(k)
  next k
  print

END

'read matrix A from screen
1000 print
  print " Input matrix size (max ";NMAX;"): "; : input "", n
  for i=1 to n
    print : print " Line "; i
    for j=1 to n
      print "  Element "; j; ": "; : input "", A(i,j)
   next j
  next i
return

'calculate the trace of real matrix C
'Function TRM(n:integer; C:matrix):Double
1500 t0=0.0
  for i=1 to n
    t0=t0+C(i,i)
  next i
return

'*****************************************************
'* The procedure PCMS calculates the coefficients of *
'* the characteristic polynomial P(l)=A-l*I of a real*
'* square 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: real matrix (n,n)                     *
'*             (see main program).                   *
'* OUTPUT:                                           *
'*          P: vector of coefficients of P(l)        *
'*                                                   *
'*****************************************************
2000 dim B(NMAX,NMAX), C(NMAX,NMAX)
  if INT(n/2) <> n/2 then
    P(1)=-1#
  else
    P(1)= 1#
  end if
  for l=1 to n
    if (l=1) then
      for i=1 to n
        for j=1 to n
          C(i,j)=A(i,j)
        next j
      next i
    else
      for i=1 to n
        for j=1 to n
          s=0#
          for k=1 to n
            s=s+B(i,k)*A(k,j)
          next k
          C(i,j)=s
        next j
      next i
    end if
    gosub 1500  't0=TRM(n,C)/l
    t0=t0/L
    P(l+1)=-t0*P(1)
    if l<n then
      for i=1 to n
        for j=1 to n
          if j=i then
            B(i,j)=C(i,j)-t0
          else
            B(i,j)=C(i,j)
          end if
        next j
      next i
    end if
  next l
return

'end of file carpol2.bas