Initial population.

1 files changed, 340 insertions, 0 deletions

diff --git a/src/algebra/eigen.spad.pamphlet b/src/algebra/eigen.spad.pamphlet
new file mode 100644
index 00000000..193b58f9
--- /dev/null
+++ b/src/algebra/eigen.spad.pamphlet
@@ -0,0 +1,340 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/algebra eigen.spad}
+\author{Patrizia Gianni, Barry Trager}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{package EP EigenPackage}
+<<package EP EigenPackage>>=
+)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 p^=0 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)
+
+@
+\section{package CHARPOL CharacteristicPolynomialPackage}
+<<package CHARPOL CharacteristicPolynomialPackage>>=
+)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
+
+@
+\section{License}
+<<license>>=
+--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.
+@
+<<*>>=
+<<license>>
+
+<<package EP EigenPackage>>
+<<package CHARPOL CharacteristicPolynomialPackage>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}