path: root/src/algebra/numeigen.spad.pamphlet

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra numeigen.spad}
\author{Patrizia Gianni}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package INEP InnerNumericEigenPackage}
<<package INEP InnerNumericEigenPackage>>=
)abbrev package INEP InnerNumericEigenPackage
++ Author:P. Gianni
++ Date Created: Summer 1990
++ Date Last Updated:Spring 1991
++ Basic Functions:
++ Related Constructors: ModularField
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This package is the inner package to be used by NumericRealEigenPackage
++ and NumericComplexEigenPackage for the computation of numeric
++ eigenvalues and eigenvectors.
InnerNumericEigenPackage(K,F,Par) : C == T
 where
   F    :   Field  -- this is the field where the answer will be
                   -- for dealing with the complex case
   K    :   Field  -- type of the  input
   Par  :   Join(Field,OrderedRing)  -- it will be NF or RN

   SE    ==> Symbol()
   RN    ==> Fraction Integer
   I     ==> Integer
   NF    ==> Float
   CF    ==> Complex Float
   GRN   ==> Complex RN
   GI    ==> Complex Integer
   PI    ==> PositiveInteger
   NNI   ==> NonNegativeInteger
   MRN   ==> Matrix RN

   MK          ==> Matrix K
   PK          ==> Polynomial K
   MF          ==> Matrix F
   SUK         ==> SparseUnivariatePolynomial K
   SUF         ==> SparseUnivariatePolynomial F
   SUP         ==> SparseUnivariatePolynomial
   MSUK        ==> Matrix SUK

   PEigenForm  ==> Record(algpol:SUK,almult:Integer,poleigen:List(MSUK))

   outForm     ==> Record(outval:F,outmult:Integer,outvect:List MF)

   IntForm     ==> Union(outForm,PEigenForm)
   UFactor     ==> (SUK -> Factored SUK)
   C == with

     charpol  :  MK   ->  SUK
       ++ charpol(m) computes the characteristic polynomial of a matrix
       ++ m with entries in K.
       ++ This function returns a polynomial
       ++ over K, while the general one (that is in EiegenPackage) returns
       ++ Fraction P K

     solve1   : (SUK, Par) -> List F
       ++ solve1(pol, eps) finds the roots of the univariate polynomial
       ++ polynomial pol to precision eps. If K is \spad{Fraction Integer}
       ++ then only the real roots are returned, if K is
       ++ \spad{Complex Fraction Integer} then all roots are found.

     innerEigenvectors    : (MK,Par,UFactor)   ->  List(outForm)
       ++ innerEigenvectors(m,eps,factor) computes explicitly
       ++ the eigenvalues and the correspondent eigenvectors
       ++ of the matrix m. The parameter eps determines the type of
       ++ the output, factor is the univariate factorizer to br used
       ++ to reduce the characteristic polynomial into irreducible factors.

   T == add

     numeric(r:K):F ==
       K is RN =>
         F is NF => convert(r)$RN
         F is RN    => r
         F is CF    => r :: RN :: CF
         F is GRN   => r::RN::GRN
       K is GRN =>
         F is GRN => r
         F is CF  => convert(convert r)
       error "unsupported coefficient type"

    ---- next functions neeeded for defining  ModularField ----

     monicize(f:SUK) : SUK ==
       (a:=leadingCoefficient f) =1 => f
       inv(a)*f

     reduction(u:SUK,p:SUK):SUK == u rem p

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

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

         ----  eval a vector of F in a radical expression  ----
     evalvect(vect:MSUK,alg:F) : MF ==
       n:=nrows vect
       w:MF:=zero(n,1)$MF
       for i in 1..n repeat
         polf:=map(numeric,
           vect(i,1))$UnivariatePolynomialCategoryFunctions2(K,SUK,F,SUF)
         v:F:=elt(polf,alg)
         setelt(w,i,1,v)
       w

       ---- internal function for the computation of eigenvectors  ----
     inteigen(A:MK,p:SUK,fact:UFactor) : List(IntForm) ==
       dimA:NNI:=  nrows A
       MM:=ModularField(SUK,SUK,reduction,merge,exactquo)
       AM:=Matrix(MM)
       lff:=factors fact(p)
       res: List IntForm  :=[]
       lr : List MF:=[]
       for ff in lff repeat
         pol:SUK:= ff.factor
         if (degree pol)=1 then
           alpha:K:=-coefficient(pol,0)/leadingCoefficient pol
           -- compute the eigenvectors, rational case
           B1:MK := zero(dimA,dimA)$MK
           for i in 1..dimA repeat
             for j in 1..dimA repeat B1(i,j):=A(i,j)
             B1(i,i):= B1(i,i) - alpha
           lr:=[]
           for vecr in nullSpace B1 repeat
             wf:MF:=zero(dimA,1)
             for i in 1..dimA repeat wf(i,1):=numeric vecr.i
             lr:=cons(wf,lr)
           res:=cons([numeric alpha,ff.exponent,lr]$outForm,res)
         else
           ppol:=monicize pol
           alg:MM:= reduce(monomial(1,1),ppol)
           B:AM:= zero(dimA,dimA)$AM
           for i in 1..dimA  repeat
             for j in 1..dimA repeat B(i,j):=reduce(A(i,j) ::SUK,ppol)
             B(i,i):=B(i,i) - alg
           sln2:=nullSpace B
           soln:List MSUK :=[]
           for vec in sln2 repeat
             wk:MSUK:=zero(dimA,1)
             for i in 1..dimA repeat wk(i,1):=(vec.i)::SUK
             soln:=cons(wk,soln)
           res:=cons([ff.factor,ff.exponent,soln]$PEigenForm,
                            res)
       res

     if K is RN then
         solve1(up:SUK, eps:Par) : List(F) ==
           denom := "lcm"/[denom(c::RN) for c in coefficients up]
           up:=denom*up
           upi := map(numer,up)$UnivariatePolynomialCategoryFunctions2(RN,SUP RN,I,SUP I)
           innerSolve1(upi, eps)$InnerNumericFloatSolvePackage(I,F,Par)
     else if K is GRN then
         solve1(up:SUK, eps:Par) : List(F) ==
           denom := "lcm"/[lcm(denom real(c::GRN), denom imag(c::GRN))
                                for c in coefficients up]
           up:=denom*up
           upgi := map(complex(numer(real #1), numer(imag #1)),
                      up)$UnivariatePolynomialCategoryFunctions2(GRN,SUP GRN,GI,SUP GI)
           innerSolve1(upgi, eps)$InnerNumericFloatSolvePackage(GI,F,Par)
     else error "unsupported matrix type"

          ----  the real eigenvectors expressed as floats  ----

     innerEigenvectors(A:MK,eps:Par,fact:UFactor) : List outForm ==
       pol:= charpol A
       sln1:List(IntForm):=inteigen(A,pol,fact)
       n:=nrows A
       sln:List(outForm):=[]
       for lev in sln1 repeat
         lev case outForm => sln:=cons(lev,sln)
         leva:=lev::PEigenForm
         lval:List(F):= solve1(leva.algpol,eps)
         lvect:=leva.poleigen
         lmult:=leva.almult
         for alg in lval repeat
           nsl:=[alg,lmult,[evalvect(ep,alg) for ep in lvect]]$outForm
           sln:=cons(nsl,sln)
       sln

     charpol(A:MK) : SUK ==
       dimA :PI := (nrows A):PI
       dimA ~= ncols A => error " The matrix is not square"
       B:Matrix SUK :=zero(dimA,dimA)
       for i in 1..dimA repeat
         for j in 1..dimA repeat  B(i,j):=A(i,j)::SUK
         B(i,i) := B(i,i) - monomial(1,1)$SUK
       determinant B


@
\section{package NREP NumericRealEigenPackage}
<<package NREP NumericRealEigenPackage>>=
)abbrev package NREP NumericRealEigenPackage
++ Author:P. Gianni
++ Date Created:Summer 1990
++ Date Last Updated:Spring 1991
++ Basic Functions:
++ Related Constructors: FloatingRealPackage
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This package computes explicitly eigenvalues and eigenvectors of
++ matrices with entries over the Rational Numbers.
++ The results are expressed as floating numbers or as rational numbers
++ depending on the type of the parameter Par.
NumericRealEigenPackage(Par) : C == T
 where
   Par   :   Join(Field,OrderedRing) -- Float or RationalNumber

   SE    ==> Symbol()
   RN    ==> Fraction Integer
   I     ==> Integer
   NF    ==> Float
   CF    ==> Complex Float
   GRN   ==> Complex RN
   GI    ==> Complex Integer
   PI    ==> PositiveInteger
   NNI   ==> NonNegativeInteger
   MRN   ==> Matrix RN

   MPar        ==> Matrix Par
   outForm     ==> Record(outval:Par,outmult:Integer,outvect:List MPar)

   C == with
     characteristicPolynomial :   MRN    -> Polynomial RN
       ++ characteristicPolynomial(m) returns the characteristic polynomial
       ++ of the matrix m expressed as polynomial
       ++ over RN with a new symbol as variable.
       -- while the function in EigenPackage returns Fraction P RN.
     characteristicPolynomial : (MRN,SE) -> Polynomial RN
       ++ characteristicPolynomial(m,x) returns the characteristic polynomial
       ++ of the matrix m expressed as polynomial
       ++ over RN with variable x.
       -- while the function in EigenPackage returns
       ++ Fraction P RN.
     realEigenvalues  :   (MRN,Par)   ->  List Par
       ++ realEigenvalues(m,eps) computes the eigenvalues of the matrix
       ++ m to precision eps. The eigenvalues are expressed as floats or
       ++ rational numbers depending on the type of eps (float or rational).
     realEigenvectors    : (MRN,Par)   ->  List(outForm)
       ++ realEigenvectors(m,eps)  returns a list of
       ++ records each one containing
       ++ a real eigenvalue, its algebraic multiplicity, and a list of
       ++ associated eigenvectors. All these results
       ++ are computed to precision eps as floats or rational
       ++ numbers depending on the type of eps .


   T == add

     import InnerNumericEigenPackage(RN, Par, Par)

     characteristicPolynomial(m:MRN) : Polynomial RN ==
       x:SE:=new()$SE
       multivariate(charpol(m),x)

            ----  characteristic polynomial of a matrix A ----
     characteristicPolynomial(A:MRN,x:SE):Polynomial RN ==
       multivariate(charpol(A),x)

     realEigenvalues(m:MRN,eps:Par) : List Par  ==
       solve1(charpol m, eps)

     realEigenvectors(m:MRN,eps:Par) :List outForm ==
       innerEigenvectors(m,eps,factor$GenUFactorize(RN))

@
\section{package NCEP NumericComplexEigenPackage}
<<package NCEP NumericComplexEigenPackage>>=
)abbrev package NCEP NumericComplexEigenPackage
++ Author: P. Gianni
++ Date Created: Summer 1990
++ Date Last Updated: Spring 1991
++ Basic Functions:
++ Related Constructors: FloatingComplexPackage
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This package computes explicitly eigenvalues and eigenvectors of
++ matrices with entries over the complex rational numbers.
++ The results are expressed either as complex floating numbers or as
++ complex rational numbers
++ depending on the type of the precision parameter.
NumericComplexEigenPackage(Par) : C == T
 where
   Par   :   Join(Field,OrderedRing)   -- Float or RationalNumber

   SE    ==> Symbol()
   RN    ==> Fraction Integer
   I     ==> Integer
   NF    ==> Float
   CF    ==> Complex Float
   GRN   ==> Complex RN
   GI    ==> Complex Integer
   PI    ==> PositiveInteger
   NNI   ==> NonNegativeInteger
   MRN   ==> Matrix RN

   MCF         ==> Matrix CF
   MGRN        ==> Matrix GRN
   MCPar       ==> Matrix Complex Par
   SUPGRN      ==> SparseUnivariatePolynomial GRN
   outForm     ==> Record(outval:Complex Par,outmult:Integer,outvect:List MCPar)

   C == with
     characteristicPolynomial :   MGRN    -> Polynomial GRN
       ++ characteristicPolynomial(m) returns the characteristic polynomial
       ++ of the matrix m expressed as polynomial
       ++ over complex rationals with a new symbol as variable.
       -- while the function in EigenPackage returns Fraction P GRN.
     characteristicPolynomial : (MGRN,SE) -> Polynomial GRN
       ++ characteristicPolynomial(m,x) returns the characteristic polynomial
       ++ of the matrix m expressed as polynomial
       ++ over Complex Rationals with variable x.
       -- while the function in EigenPackage returns Fraction P GRN.
     complexEigenvalues  :   (MGRN,Par)   ->  List Complex Par
       ++ complexEigenvalues(m,eps) computes the eigenvalues of the matrix
       ++ m to precision eps. The eigenvalues are expressed as complex floats or
       ++ complex rational numbers depending on the type of eps (float or rational).
     complexEigenvectors    : (MGRN,Par)   ->  List(outForm)
       ++ complexEigenvectors(m,eps)  returns a list of
       ++ records each one containing
       ++ a complex eigenvalue, its algebraic multiplicity, and a list of
       ++ associated eigenvectors. All these results
       ++ are computed to precision eps and are expressed as complex floats
       ++ or complex rational numbers depending on the type of eps (float or rational).
   T == add

     import InnerNumericEigenPackage(GRN,Complex Par,Par)

     characteristicPolynomial(m:MGRN) : Polynomial GRN  ==
       x:SE:=new()$SE
       multivariate(charpol m, x)

            ----  characteristic polynomial of a matrix A ----
     characteristicPolynomial(A:MGRN,x:SE):Polynomial GRN ==
       multivariate(charpol A, x)

     complexEigenvalues(m:MGRN,eps:Par) : List Complex Par  ==
       solve1(charpol m, eps)

     complexEigenvectors(m:MGRN,eps:Par) :List outForm ==
       innerEigenvectors(m,eps,factor$ComplexFactorization(RN,SUPGRN))

@
\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 INEP InnerNumericEigenPackage>>
<<package NREP NumericRealEigenPackage>>
<<package NCEP NumericComplexEigenPackage>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}