--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

)abbrev package EP EigenPackage
++ Author: P. Gianni
++ Date Created: summer 1986
++ Date Last Updated: October 1992
++ Basic Functions:
++ Related Constructors: NumericRealEigenPackage,  NumericComplexEigenPackage,
++ RadicalEigenPackage
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This is a package for the exact computation of eigenvalues and eigenvectors.
++ This package can be made to work for matrices with coefficients which are
++ rational functions over a ring where we can factor polynomials.
++ Rational eigenvalues are always explicitly computed while the
++ non-rational ones are expressed in terms of their minimal
++ polynomial.
-- Functions for the numeric computation of eigenvalues and eigenvectors
-- are in numeigen spad.
EigenPackage(R) : C == T
 where
   R     : GcdDomain
   P     ==> Polynomial R
   F     ==> Fraction P
   SE    ==> Symbol()
   SUP   ==> SparseUnivariatePolynomial(P)
   SUF   ==> SparseUnivariatePolynomial(F)
   M     ==> Matrix(F)
   NNI   ==> NonNegativeInteger
   ST    ==> SuchThat(SE,P)

   Eigenvalue  ==> Union(F,ST)
   EigenForm   ==> Record(eigval:Eigenvalue,eigmult:NNI,eigvec : List M)
   GenEigen    ==> Record(eigval:Eigenvalue,geneigvec:List M)

   C == with
     characteristicPolynomial :  (M,Symbol)  ->  P
       ++ characteristicPolynomial(m,var) returns the
       ++ characteristicPolynomial of the matrix m using
       ++ the symbol var as the main variable.

     characteristicPolynomial :      M       ->  P
       ++ characteristicPolynomial(m) returns the
       ++ characteristicPolynomial of the matrix m using
       ++ a new generated symbol symbol as the main variable.

     eigenvalues       :    M        ->  List Eigenvalue
       ++ eigenvalues(m) returns the
       ++ eigenvalues of the matrix m which are expressible
       ++ as rational functions over the rational numbers.

     eigenvector       :   (Eigenvalue,M)  ->  List M
       ++ eigenvector(eigval,m) returns the
       ++ eigenvectors belonging to the eigenvalue eigval
       ++ for the matrix m.

     generalizedEigenvector  : (Eigenvalue,M,NNI,NNI) -> List M
       ++ generalizedEigenvector(alpha,m,k,g)
       ++ returns the generalized eigenvectors
       ++ of the matrix relative to the eigenvalue alpha.
       ++ The integers k and g are respectively the algebraic and the
       ++ geometric multiplicity of tye eigenvalue alpha.
       ++ alpha can be either rational or not.
       ++ In the seconda case apha is the minimal polynomial of the
       ++ eigenvalue.

     generalizedEigenvector  : (EigenForm,M) -> List M
       ++ generalizedEigenvector(eigen,m)
       ++ returns the generalized eigenvectors
       ++ of the matrix relative to the eigenvalue eigen, as 
       ++ returned by the function eigenvectors.

     generalizedEigenvectors  : M -> List GenEigen
       ++ generalizedEigenvectors(m)
       ++ returns the generalized eigenvectors
       ++ of the matrix m.

     eigenvectors      :    M        ->  List(EigenForm)
       ++ eigenvectors(m) returns the eigenvalues and eigenvectors
       ++ for the matrix m.
       ++ The rational eigenvalues and the correspondent eigenvectors
       ++ are explicitely computed, while the non rational ones
       ++ are given via their minimal polynomial and the corresponding
       ++ eigenvectors are expressed in terms of a "generic" root of
       ++ such a polynomial.

   T == add
     PI       ==> PositiveInteger


     MF  := GeneralizedMultivariateFactorize(SE,IndexedExponents SE,R,R,P)
     UPCF2:= UnivariatePolynomialCategoryFunctions2(P,SUP,F,SUF)
    

                 ----  Local  Functions  ----
     tff              :  (SUF,SE)      ->  F
     fft              :  (SUF,SE)      ->  F
     charpol          :   (M,SE)       ->   F
     intRatEig        :  (F,M,NNI)    ->   List M
     intAlgEig        :  (ST,M,NNI)    ->   List M 
     genEigForm       : (EigenForm,M)  ->   GenEigen

    ---- next functions needed for defining  ModularField ----
     reduction(u:SUF,p:SUF):SUF == u rem p

     merge(p:SUF,q:SUF):Union(SUF,"failed") ==
         p = q => p
         p = 0 => q
         q = 0 => p
         "failed"

     exactquo(u:SUF,v:SUF,p:SUF):Union(SUF,"failed") ==
        val:=extendedEuclidean(v,p,u)
        val case "failed" => "failed"
        val.coef1

               ----  functions for conversions  ----
     fft(t:SUF,x:SE):F ==
       n:=degree(t)
       cf:=monomial(1,x,n)$P :: F
       cf * leadingCoefficient t

     tff(p:SUF,x:SE) : F ==
       degree p=0 => leadingCoefficient p
       r:F:=0$F
       while not zero? p repeat
         r:=r+fft(p,x)
         p := reductum p
       r

      ---- generalized eigenvectors associated to a given eigenvalue ---       
     genEigForm(eigen : EigenForm,A:M) : GenEigen ==
       alpha:=eigen.eigval
       k:=eigen.eigmult
       g:=#(eigen.eigvec)
       k = g  => [alpha,eigen.eigvec]
       [alpha,generalizedEigenvector(alpha,A,k,g)]

           ---- characteristic polynomial  ----
     charpol(A:M,x:SE) : F ==
       dimA :PI := (nrows A):PI
       dimA ~= ncols A => error " The matrix is not square"
       B:M:=zero(dimA,dimA)
       for i in 1..dimA repeat
         for j in 1..dimA repeat  B(i,j):=A(i,j)
         B(i,i) := B(i,i) - monomial(1$P,x,1)::F
       determinant B

          --------  EXPORTED  FUNCTIONS  --------
   
            ----  characteristic polynomial of a matrix A ----
     characteristicPolynomial(A:M):P ==
       x:SE:=new()$SE
       numer charpol(A,x)

            ----  characteristic polynomial of a matrix A ----
     characteristicPolynomial(A:M,x:SE) : P == numer charpol(A,x)
     
                ----  Eigenvalues of the matrix A  ----
     eigenvalues(A:M): List Eigenvalue  ==
       x:=new()$SE
       pol:= charpol(A,x)
       lrat:List F :=empty()
       lsym:List ST :=empty()
       for eq in solve(pol,x)$SystemSolvePackage(R) repeat
         alg:=numer lhs eq
         degree(alg, x)=1 => lrat:=cons(rhs eq,lrat)
         lsym:=cons([x,alg],lsym)
       append([lr::Eigenvalue for lr in lrat],
              [ls::Eigenvalue for ls in lsym])

          ----  Eigenvectors belonging to a given eigenvalue  ----
                ----  the eigenvalue must be exact  ----
     eigenvector(alpha:Eigenvalue,A:M) : List M  ==
       alpha case F => intRatEig(alpha::F,A,1$NNI)
       intAlgEig(alpha::ST,A,1$NNI)

   ----  Eigenvectors belonging to a given rational eigenvalue  ----
                ---- Internal function -----
     intRatEig(alpha:F,A:M,m:NNI) : List M  ==
       n:=nrows A
       B:M := zero(n,n)$M
       for i in 1..n repeat
         for j in 1..n repeat B(i,j):=A(i,j)
         B(i,i):= B(i,i) - alpha
       [v::M for v in nullSpace(B**m)]
   
   ----  Eigenvectors belonging to a given algebraic eigenvalue  ----
         ------   Internal  Function  -----
     intAlgEig(alpha:ST,A:M,m:NNI) : List M  ==
       n:=nrows A
       MM := ModularField(SUF,SUF,reduction,merge,exactquo)
       AM:=Matrix MM
       x:SE:=lhs alpha
       pol:SUF:=unitCanonical map(coerce,univariate(rhs alpha,x))$UPCF2
       alg:MM:=reduce(monomial(1,1),pol)
       B:AM := zero(n,n)
       for i in 1..n repeat
         for j in 1..n repeat B(i,j):=reduce(A(i,j)::SUF,pol)
         B(i,i):= B(i,i) - alg
       sol: List M :=empty()
       for vec in nullSpace(B**m) repeat
         w:M:=zero(n,1)
         for i in 1..n repeat w(i,1):=tff((vec.i)::SUF,x)
         sol:=cons(w,sol)
       sol

     ----  Generalized Eigenvectors belonging to a given eigenvalue  ----
     generalizedEigenvector(alpha:Eigenvalue,A:M,k:NNI,g:NNI) : List M  ==
       alpha case F => intRatEig(alpha::F,A,(1+k-g)::NNI)
       intAlgEig(alpha::ST,A,(1+k-g)::NNI)

     ----  Generalized Eigenvectors belonging to a given eigenvalue  ----
     generalizedEigenvector(eigen :EigenForm,A:M) : List M  ==
       generalizedEigenvector(eigen.eigval,A,eigen.eigmult,# eigen.eigvec)

          ----  Generalized Eigenvectors -----
     generalizedEigenvectors(A:M) : List GenEigen  ==
       n:= nrows A
       leig:=eigenvectors A
       [genEigForm(leg,A) for leg in leig]
         
                 ----  eigenvectors and eigenvalues  ----
     eigenvectors(A:M):List(EigenForm) ==
       n:=nrows A
       x:=new()$SE
       p:=numer charpol(A,x)
       MM := ModularField(SUF,SUF,reduction,merge,exactquo)
       AM:=Matrix(MM)
       ratSol : List EigenForm := empty()
       algSol : List EigenForm := empty()
       lff:=factors factor  p
       for fact in lff repeat
         pol:=fact.factor   
         degree(pol,x)=1 =>
           vec:F :=-coefficient(pol,x,0)/coefficient(pol,x,degree(pol,x))
           ratSol:=cons([vec,fact.exponent :: NNI,
                         intRatEig(vec,A,1$NNI)]$EigenForm,ratSol)
         alpha:ST:=[x,pol]     
         algSol:=cons([alpha,fact.exponent :: NNI,
                       intAlgEig(alpha,A,1$NNI)]$EigenForm,algSol)
       append(ratSol,algSol)

)abbrev package CHARPOL CharacteristicPolynomialPackage
++ Author: Barry Trager
++ Date Created:
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This package provides a characteristicPolynomial function
++ for any matrix over a commutative ring.

CharacteristicPolynomialPackage(R:CommutativeRing):C == T where
   PI ==> PositiveInteger
   M ==> Matrix R
   C == with
      characteristicPolynomial: (M, R) -> R
        ++ characteristicPolynomial(m,r) computes the characteristic
        ++ polynomial of the matrix m evaluated at the point r.
        ++ In particular, if r is the polynomial 'x, then it returns
        ++ the characteristic polynomial expressed as a polynomial in 'x.
   T == add

           ---- characteristic polynomial  ----
     characteristicPolynomial(A:M,v:R) : R ==
       dimA :PI := (nrows A):PI
       dimA ~= ncols A => error " The matrix is not square"
       B:M:=zero(dimA,dimA)
       for i in 1..dimA repeat
         for j in 1..dimA repeat  B(i,j):=A(i,j)
         B(i,i) := B(i,i) - v
       determinant B

