/usr/share/axiom-20170501/src/algebra/INEP.spad is in axiom-source 20170501-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 | )abbrev package INEP InnerNumericEigenPackage
++ Author:P. Gianni
++ Date Created: Summer 1990
++ Date Last Updated:Spring 1991
++ 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) : SIG == CODE 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)
SIG ==> 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.
CODE ==> 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((c:GRN):GI+->complex(numer(real c), numer(imag c)),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
|