\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra radeigen.spad}
\author{Patrizia Gianni}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package REP RadicalEigenPackage}
<<package REP RadicalEigenPackage>>=
)abbrev package REP RadicalEigenPackage
++ Author: P.Gianni
++ Date Created: Summer 1987
++ Date Last Updated: October 1992
++ Basic Functions:
++ Related Constructors: EigenPackage, RadicalSolve
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ Package for the computation of eigenvalues and eigenvectors.
++ This package works for matrices with coefficients which are
++ rational functions over the integers.
++ (see \spadtype{Fraction Polynomial Integer}).
++ The eigenvalues and eigenvectors are expressed in terms of radicals.
RadicalEigenPackage() : C == T
 where
   R     ==> Integer
   P     ==> Polynomial R
   F     ==> Fraction P
   RE    ==> Expression R
   SE    ==> Symbol()
   M     ==> Matrix(F)
   MRE   ==> Matrix(RE)
   ST    ==> SuchThat(SE,P)
   NNI   ==> NonNegativeInteger

   EigenForm    ==> Record(eigval:Union(F,ST),eigmult:NNI,eigvec:List(M))
   RadicalForm  ==> Record(radval:RE,radmult:Integer,radvect:List(MRE))



   C == with
     radicalEigenvectors    :  M     ->  List(RadicalForm)
       ++ radicalEigenvectors(m) computes
       ++ the eigenvalues and the corresponding eigenvectors of the
       ++ matrix m;
       ++ when possible, values are expressed in terms of radicals.

     radicalEigenvector       :  (RE,M)     ->  List(MRE)
       ++ radicalEigenvector(c,m) computes the eigenvector(s) of the
       ++ matrix m corresponding to the eigenvalue c;
       ++ when possible, values are
       ++ expressed in terms of radicals.

     radicalEigenvalues     :  M     ->  List RE
       ++ radicalEigenvalues(m) computes the eigenvalues of the matrix m;
       ++ when possible, the eigenvalues are expressed in terms of radicals.

     eigenMatrix       :    M        ->  Union(MRE,"failed")
       ++ eigenMatrix(m) returns the matrix b
       ++ such that \spad{b*m*(inverse b)} is diagonal,
       ++ or "failed" if no such b exists.

     normalise         :    MRE      ->  MRE
       ++ normalise(v) returns the column
       ++ vector v
       ++ divided by its euclidean norm;
       ++ when possible, the vector v is expressed in terms of radicals.

     gramschmidt      : List(MRE)   ->  List(MRE)
       ++ gramschmidt(lv) converts the list of column vectors lv into
       ++ a set of orthogonal column vectors
       ++ of euclidean length 1 using the Gram-Schmidt algorithm.

     orthonormalBasis  :    M        ->  List(MRE)
       ++ orthonormalBasis(m)  returns the orthogonal matrix b such that
       ++ \spad{b*m*(inverse b)} is diagonal.
       ++ Error: if m is not a symmetric matrix.

   T == add
     PI       ==> PositiveInteger
     RSP := RadicalSolvePackage R
     import EigenPackage R

                 ----  Local  Functions  ----
     evalvect         :  (M,RE,SE)  ->  MRE
     innerprod        :   (MRE,MRE)   ->  RE

         ----  eval a vector of F in a radical expression  ----
     evalvect(vect:M,alg:RE,x:SE) : MRE ==
       n:=nrows vect
       xx:=kernel(x)$Kernel(RE)
       w:MRE:=zero(n,1)$MRE
       for i in 1..n repeat
         v:=eval(vect(i,1) :: RE,xx,alg)
         setelt(w,i,1,v)
       w
                      ---- inner product ----
     innerprod(v1:MRE,v2:MRE): RE == (((transpose v1)* v2)::MRE)(1,1)

                 ----  normalization of a vector  ----
     normalise(v:MRE) : MRE ==
       normv:RE := sqrt(innerprod(v,v))
       normv = 0$RE => v
       (1/normv)*v

                ----  Eigenvalues of the matrix A  ----
     radicalEigenvalues(A:M): List(RE) ==
       x:SE :=new()$SE
       pol:= characteristicPolynomial(A,x) :: F
       radicalRoots(pol,x)$RSP

      ----  Eigenvectors belonging to a given eigenvalue  ----
            ----  expressed in terms of radicals ----
     radicalEigenvector(alpha:RE,A:M) : List(MRE) ==
       n:=nrows A
       B:MRE := zero(n,n)$MRE
       for i in 1..n repeat
         for j in 1..n repeat B(i,j):=(A(i,j))::RE
         B(i,i):= B(i,i) - alpha
       [v::MRE  for v in nullSpace B]

             ----  eigenvectors and eigenvalues  ----
     radicalEigenvectors(A:M) : List(RadicalForm) ==
       leig:List EigenForm := eigenvectors A
       n:=nrows A
       sln:List RadicalForm := empty()
       veclist: List MRE
       for eig in leig repeat
         eig.eigval case F =>
           veclist := empty()
           for ll in eig.eigvec repeat
             m:MRE:=zero(n,1)
             for i in 1..n repeat m(i,1):=(ll(i,1))::RE
             veclist:=cons(m,veclist)
           sln:=cons([(eig.eigval)::F::RE,eig.eigmult,veclist]$RadicalForm,sln)
         sym := eig.eigval :: ST
         xx:= lhs sym
         lval : List RE := radicalRoots((rhs sym) :: F ,xx)$RSP
         for alg in lval repeat
           nsl:=[alg,eig.eigmult,
                 [evalvect(ep,alg,xx) for ep in eig.eigvec]]$RadicalForm
           sln:=cons(nsl,sln)
       sln

            ----  orthonormalization of a list of vectors  ----
                  ----  Grahm - Schmidt process  ----

     gramschmidt(lvect:List(MRE)) : List(MRE) ==
       lvect=[]  => []
       v:=lvect.first
       n := nrows v
       RMR:=RectangularMatrix(n:PI,1,RE)
       orth:List(MRE):=[(normalise v)]
       for v in lvect.rest repeat
         pol:=((v:RMR)-(+/[(innerprod(w,v)*w):RMR for w in orth])):MRE
         orth:=cons(normalise pol,orth)
       orth


              ----  The matrix of eigenvectors  ----

     eigenMatrix(A:M) : Union(MRE,"failed") ==
       lef:List(MRE):=[:eiv.radvect  for eiv in radicalEigenvectors(A)]
       n:=nrows A
       #lef <n => "failed"
       d:MRE:=copy(lef.first)
       for v in lef.rest repeat d:=(horizConcat(d,v))::MRE
       d

         ----  orthogonal basis for a symmetric matrix  ----

     orthonormalBasis(A:M):List(MRE) ==
       not symmetric?(A) => error "the matrix is not symmetric"
       basis:List(MRE):=[]
       lvec:List(MRE) := []
       alglist:List(RadicalForm):=radicalEigenvectors(A)
       n:=nrows A
       for alterm in alglist repeat
         if (lvec:=alterm.radvect)=[] then error "sorry "
         if #(lvec)>1  then
           lvec:= gramschmidt(lvec)
           basis:=[:lvec,:basis]
         else basis:=[normalise(lvec.first),:basis]
       basis

@
\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 REP RadicalEigenPackage>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}