'*********************************************************
'* Calculate the coefficients of the characteristic poly-*
'* nomial P(l)=A-l*I of a real tridiagonal matrix A(i,j) *
'* ----------------------------------------------------- *
'* Ref.: "Algebre, Algorithmes et programmes en Pascal   *
'*        By Jean-Louis Jardrin, DUNOD Paris, 1988" [10].*
'* ----------------------------------------------------- *
'* SAMPLE RUN:                                           *
'* (Find the caracteristic polynomial of matrix:         *
'*           4 1 0 0 0 0                                 *
'*           2 4 1 0 0 0                                 *
'*      A =  0 2 4 1 0 0                                 *
'*           0 0 2 4 1 0                                 *
'*           0 0 0 2 4 1                                 *
'*           0 0 0 0 2 4   )                             *
'*                                                       *
'* Input matrix size (max 100): 6                        *
'* Input 5 elements of subdiagonal:                      *
'*   Element 1: 2                                        *
'*   Element 2: 2                                        *
'*   Element 3: 2                                        *
'*   Element 4: 2                                        *
'*   Element 5: 2                                        *
'* Input 6 elements of main diagonal:                    *
'*   Element 1: 4                                        *
'*   Element 2: 4                                        *
'*   Element 3: 4                                        *
'*   Element 4: 4                                        *
'*   Element 5: 4                                        *
'*   Element 6: 4                                        *
'* Input 5 elements of supradiagonal:                    *
'*   Element 1: 1                                        *
'*   Element 2: 1                                        *
'*   Element 3: 1                                        *
'*   Element 4: 1                                        *
'*   Element 5: 1                                        *
'*                                                       *
'* Chracteristic polynomial P(l):                        *
'* Coefficient for degree 6:      1.0000                 *
'* Coefficient for degree 5:    -24.0000                 *
'* Coefficient for degree 4:    230.0000                 *
'* Coefficient for degree 3:  -1120.0000                 *
'* Coefficient for degree 2:   2904.0000                 *
'* Coefficient for degree 1:  -3776.0000                 *
'* Coefficient for degree 0:   1912.0000                 *
'*                                                       *
'* (The characteristic polynomial of matrix A is:        *
'*   P(l)=l^6-24*l^5+230*l^4-1120*l^3+2904*l^2-3376*l    *
'*        +1912   )                                      *
'* ----------------------------------------------------- *
'* Subroutines used:                                     *
'*                                                       *
'* 1000 ReadData: Read elements of input matrix          *
'* 2000 PCTR    : Calculate the coefficients  of the     *
'*                characteristic polynomial P(l)=A-l*I   *
'*                of a real square tridiagonal matrix.   *
'*                                                       * 
'*                  Basic version by J-P Moreau, Paris.  *
'*                           (www.jpmoreau.fr)           *
'*********************************************************
'Program CARPOL
DEFDBL A-H, O-Z
DEFINT I-N

DIM B(25), C(25), D(25), E(25), P(25)

  CLS
  PRINT
  INPUT " Input size of matrix: ", n

  GOSUB 1000  'call Read_data(n,C,D,E)

' prepare B vector
  FOR k = 1 TO n - 1
    B(k) = C(k) * E(k)
  NEXT k
' call routine to calculate coefficients
  GOSUB 2000 'call PCTR(n,B,D,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 'of main program

1000 'Subroutine ReadData
'integer n, i
  PRINT
  PRINT " Input "; n - 1; " elements of subdiagonal:"
  FOR i = 1 TO n - 1
    PRINT "   Element "; i; ": "; : INPUT "", C(i)
  NEXT
  PRINT " Input "; n; " elements of main diagonal:"
  FOR i = 1 TO n
    PRINT "   Element "; i; ": "; : INPUT "", D(i)
  NEXT
  PRINT " Input "; n - 1; " elements of supradiagonal:"
  FOR i = 1 TO n - 1
    PRINT "   Element "; i; ": "; : INPUT "", E(i)
  NEXT
RETURN

'*****************************************************
'*  The subroutine PCTR calculates the coefficients  *
'*  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)              *
'*          B: vector of codiagonal terms products   *
'*             (see main program).                   *
'*          D: vector of main diagonal terms         *
'* OUTPUT:                                           *
'*          P: vector of coefficients of P(l)        *
'*                                                   *
'*****************************************************
'Subroutine PCTR(n,B,D,P)
'integer n,i,k, ialloc
2000 DIM T(25)
  P(2) = D(1): T(1) = 1#: P(1) = -T(1)
  FOR k = 2 TO n
    P(k + 1) = D(k) * P(k) - B(k - 1) * T(k - 1)
    FOR i = k TO 3 STEP -1
      T(i) = P(i)
      P(i) = -T(i) + D(k) * P(i - 1) - B(k - 1) * T(i - 2)
    NEXT i
    T(2) = P(2): P(2) = -T(2) + D(k) * P(1)
    T(1) = P(1): P(1) = -T(1)
  NEXT k
RETURN

'end of file carpol.bas