module glmat; % Routines for inverting matrices and finding eigen-values
% and vectors. Methods are the same as in glsolve module.
% Author: Eberhard Schruefer.
% Modification: James Davenport and Fran Burstall.
fluid '(!*cramer !*factor !*gcd !*sqfree !*sub2 kord!*);
global '(!!arbint);
if null !!arbint then !!arbint := 0;
switch cramer;
put('cramer,'simpfg,
'((t (put 'mat 'lnrsolvefn 'clnrsolve)
(put 'mat 'inversefn 'matinv))
(nil (put 'mat 'lnrsolvefn 'lnrsolve)
(put 'mat 'inversefn 'matinverse))));
% algebraic operator arbcomplex;
% Done this way since it's also defined in the solve1 module.
deflist('((arbcomplex simpiden)),'simpfn);
symbolic procedure clnrsolve(u,v);
% Interface to matrix package.
multm(matinv u,v);
symbolic procedure minv u;
matinv matsm u;
put('minv,'rtypefn,'quotematrix);
%put('mateigen,'rtypefn,'quotematrix);
remprop('mateigen,'rtypefn);
symbolic procedure matinv u;
% U is a matrix form. Result is the inverse of matrix u.
begin scalar sgn,x,y,z,!*exp;
integer l,m,lm;
!*exp := t;
z := 1;
lm := length car u;
for each v in u do
<<y := 1;
for each w in v do y := lcm(y,denr w);
m := lm;
x := list(nil . (l := l + 1)) .* negf y .+ nil;
for each j in reverse v do
<<if numr j then
x := list m .* multf(numr j,quotf(y,denr j)) .+ x;
m := m - 1>>;
z := c!:extmult(x,z)>>;
if singularchk lpow z then rerror(matrix,13,"Singular matrix");
sgn := evenp length lpow z;
return for each k in lpow z collect
<<sgn := not sgn;
for each j in lpow z collect mkglimat(k,z,sgn,j)>>
end;
symbolic procedure singularchk u; pairp car lastpair u;
flag('(mateigen),'opfn);
flag('(mateigen),'noval);
symbolic procedure mateigen(u,eival);
% U is a matrix form, eival an indeterminate naming the eigenvalues.
% Result is a list of lists:
% {{eival-eq1,multiplicity1,eigenvector1},....},
% where eival-eq is a polynomial and eigenvector is a matrix.
% How much should we attempt to solve the eigenvalue eq.? sqfr?
% Sqfr is necessary if we want to have the full eigenspace. If there
% are multiple roots another pass through eigenvector calculation
% is needed(done).
% We should actually perform the calculations in the extension
% field generated by the eigenvalue equation(done inside).
begin scalar arbvars,exu,sgn,q,r,s,x,y,z,eivec,!*factor,!*sqfree,
!*exp;
integer l;
!*exp := t;
if not(getrtype u eq 'matrix) then typerr(u,"matrix");
eival := !*a2k eival;
kord!* := eival . kord!*;
exu := mateigen1(matsm u,eival);
q := car exu;
y := cadr exu;
z := caddr exu;
exu := cdddr exu;
!*sqfree := t;
for each j in cdr fctrf numr subs2(lc z ./ 1)
do if null domainp car j and mvar car j eq eival
then s := (if null red car j
then !*k2f mvar car j . (ldeg car j*cdr j)
else j) . s;
for each j in q
do (if x then rplacd(x,cdr x + cdr j)
else s := (y . cdr j) . s)
where x := assoc(y,s) where y := absf reorder car j;
l := length s;
r := 'list .
for each j in s collect
<<if null((cdr j = 1) and (l = 1)) then
<<y := 1;
for each k in exu do
if x := reduce!-mod!-eig(car j,c!:extmult(k,y))
then y := x>>;
arbvars := nil;
for each k in lpow z do
if (y=1) or null(k member lpow y) then
arbvars := (k . makearbcomplex()) . arbvars;
sgn := (y=1) or evenp length lpow y;
eivec := 'mat . for each k in lpow z collect list
if x := assoc(k,arbvars)
then mvar cdr x
else prepsq!* mkgleig(k,y,
sgn := not sgn,arbvars);
list('list,prepsq!*(car j ./ 1),cdr j,eivec)>>;
kord!* := cdr kord!*;
return r
end;
symbolic procedure mateigen1(u,eival);
begin scalar q,x,y,z; integer l,lm,m;
lm := length car u;
z := 1;
u := for each v in u collect
<<y := 1;
for each w in v do y := lcm(y,denr w);
m := lm;
l := l + 1;
x := nil;
for each j in reverse v do
<<if numr j or l = m then
x := list m .* multf(if l = m then
addf(numr j,
negf multf(!*k2f eival,
denr j)) else numr j,
quotf(y,denr j)) .+ x;
m := m - 1>>;
y := z;
z := c!:extmult(if null red x then <<
q := (if p then (car p . (cdr p + 1)) . delete(p,q)
else (lc x . 1) . q) where p = assoc(lc x,q);
!*p2f lpow x>> else x,z);
x>>;
return q . y . z . u
end;
symbolic procedure reduce!-mod!-eig(u,v);
% Reduces exterior product v wrt eigenvalue equation u.
begin scalar x,y;
for each j on v do
if numr(y := reduce!-mod!-eigf(u,lc j)) then
x := lpow j .* y .+ x;
y := 1;
for each j on x do y := lcm(y,denr lc j);
return for each j on reverse x collect
lpow j .* multf(numr lc j,quotf(y,denr lc j))
end;
symbolic procedure reduce!-mod!-eigf(u,v);
(subs2 reduce!-eival!-powers(lpow u . negsq cancel(red u ./ lc u),v))
where !*sub2 = !*sub2;
symbolic procedure reduce!-eival!-powers(v,u);
if domainp u or null(mvar u eq caar v) then u ./ 1
else reduce!-eival!-powers1(v,u ./ 1);
symbolic procedure reduce!-eival!-powers1(v,u);
% Reduces powers with the help of the eigenvalue polynomial.
if domainp numr u or (ldeg numr u<pdeg car v) then u
else if ldeg numr u=pdeg car v then
addsq(multsq(cdr v,lc numr u ./ denr u),
red numr u ./ denr u)
else reduce!-eival!-powers1(v,
addsq(multsq(multpf(mvar numr u .** (ldeg numr u-pdeg car v),
lc numr u) ./ denr u,
cdr v),red numr u ./ denr u));
% Determinant calculation using exterior multiplication.
symbolic procedure detex u;
% U is a matrix form, result is determinant of u.
begin scalar f,x,y,z;
integer m,lm;
z := 1;
u := matsm car u;
lm := length car u;
f := 1;
for each v in u do
<<y := 1;
for each w in v do y := lcm(y,denr w);
f := multf(y,f);
m := lm;
x := nil;
for each j in v do
<<if numr j then
x := list m .* multf(numr j,quotf(y,denr j)) .+ x;
m := m - 1>>;
z := c!:extmult(x,z)>>;
return cancel(lc z ./ f)
end;
% Not supported at algebraic user level since it is in general slower
% than other methods.
% put('detex,'simpfn,'detex);
symbolic procedure mkglimat(u,v,sgn,k);
begin scalar s,x,y;
x := nil ./ 1;
y := lpow v;
for each j on red v do
if s := glmatterm(u,y,j,k)
then x := addsq(cancel(s ./ lc v),x);
return if sgn then negsq x else x
end;
symbolic procedure glmatterm(u,v,w,k);
begin scalar x,y,sgn;
x := lpow w;
a: if null x then return
if pairp car y and (cdar y = k) then lc w else nil;
if car x = u then return nil
else if car x member v then <<x := cdr x;
if y then sgn := not sgn>>
else if y then return nil
else <<y := list car x; x := cdr x>>;
go to a
end;
symbolic procedure mkgleig(u,v,sgn,arbvars);
begin scalar s,x,y,!*gcd;
x := nil ./ 1;
y := lpow v;
!*gcd := t;
for each j on red v do
if s := glsoleig(u,y,j,arbvars)
then x := addsq(cancel(s ./ lc v),x);
return if sgn then negsq x else x
end;
symbolic procedure glsoleig(u,v,w,arbvars);
begin scalar x,y,sgn;
x := lpow w;
a: if null x then return
if null car y then lc w
else multf(cdr assoc(car y,arbvars),
if sgn then negf lc w else lc w);
if car x = u then return nil
else if car x member v then <<x := cdr x;
if y then sgn := not sgn>>
else if y then return nil
else <<y := list car x; x := cdr x>>;
go to a
end;
%**** Support for exterior multiplication ****
% Data structure is lpow ::= list of col.-ind. in exterior product
% | nil . number of eq. for inhomog. terms.
% lc ::= standard form
% Exterior multiplication and p-forms:
% Let V be a vector space of dimension n.
% We call the elements of V 1-forms and build new objects called
% p-forms as follows: define a multiplication on 1-forms ^ such that
% v^w=-w^v
% then the linear span of such objects is the space of 2-forms and has
% dimension n(n-1)/2. Indeed, if v_1,...,v_n is a basis of V then
% v_i^v_j for i<j is a basis for the 2-forms.
% We extend this multiplication (called exterior multiplication)
% to get products of p vectors, linear combinations of which are
% called p-forms: this extension is defined by the rule that v_1^...^v_p
% vanishes whenever some v_i=v_j (for i not j). Thus the effect of
% permuting the order of the vectors in such a product is to multiply
% the product by the sign of the permutation.
% Using this it is not difficult to show:
% Theorem: Vectors v_1,...,v_p are linear independent iff their exterior
% product v_1^...^v_p is non-zero.
%
% For more information see F. Warner "Foundations of Differentiable
% Manifolds and Lie Groups" (Springer) Chapter 2. (Or any other book
% on differential geometry or multilinear algebra
symbolic procedure c!:extmult(u,v);
% Special exterior multiplication routine. Degree of form v is
% arbitrary, u is a one-form.
if null u or null v then nil
else if v = 1 then u %unity
% else (if x then cdr x .* (if car x
% then negf c!:subs2multf(lc u,lc v)
% else c!:subs2multf(lc u,lc v))
% .+ c!:extadd(c!:extmult(!*t2f lt u,red v),
% ^^ this is bogus: .+ may not be valid in this circumstance
% c!:extmult(red u,v))
else (if x then
c!:extadd(cdr x .* (if car x
then negf c!:subs2multf(lc u,lc v)
else c!:subs2multf(lc u,lc v)) .+ nil,
c!:extadd(c!:extmult(!*t2f lt u,red v),
c!:extmult(red u,v)))
else c!:extadd(c!:extmult(!*t2f lt u,red v),
c!:extmult(red u,v)))
where x = c!:ordexn(car lpow u,lpow v);
symbolic procedure c!:subs2multf(u,v);
(if denr x neq 1 then rerror(matrix,14,"Sub error in glnrsolve")
else numr x)
where x = subs2(multf(u,v) ./ 1) where !*sub2 = !*sub2;
symbolic procedure c!:extadd(u,v);
if null u then v
else if null v then u
else if lpow u = lpow v then
(lambda x,y; if null x then y else lpow u .* x .+ y)
(addf(lc u,lc v),c!:extadd(red u,red v))
else if c!:ordexp(lpow u,lpow v) then lt u .+ c!:extadd(red u,v)
else lt v .+ c!:extadd(u,red v);
symbolic procedure c!:ordexp(u,v);
if null u then t
else if car u = car v then c!:ordexp(cdr u,cdr v)
else c!:ordxp(car u,car v);
symbolic procedure c!:ordexn(u,v);
% U is a single index, v a list. Returns nil if u is a member
% of v or a dotted pair of a permutation indicator and the ordered
% list of u merged into v.
begin scalar s,x;
a: if null v then return(s . reverse(u . x))
else if (u = car v) or (pairp u and pairp car v)
then return nil
else if c!:ordxp(u,car v) then
return(s . append(reverse(u . x),v))
else <<x := car v . x;
v := cdr v;
s := not s>>;
go to a
end;
symbolic procedure c!:ordxp(u,v);
if pairp u then if pairp v then cdr u < cdr v
else nil
else if pairp v then t
else u < v;
endmodule;
end;