reduce-historical/r33/matr.red

module matr;

% Author: Anthony C. Hearn;

% This module is rife with essential references to RPLAC-based
% functions.

fluid '(!*sub2);

global '(nxtsym!* subfg!*);

symbolic procedure matrix u;
   %declares list U as matrices;
   begin scalar v,w,x;
        for each j in u do
           if atom j then if null (x := gettype j)
                            then put(j,'rtype,'matrix)
                           else if x eq 'matrix             
                            then <<lprim list(x,j,"redefined");
                                   put(j,'rtype,'matrix)>>
                           else typerr(list(x,j),"matrix")
            else if not idp car j
                   or length (v := revlis cdr j) neq 2
                   or not natnumlis v
             then errpri2(j,'hold)
            else if not (x := gettype car j) or x eq 'matrix
             then <<w := nil;
                    for n := 1:car v do w := nzero cadr v . w;
                    put(car j,'rtype,'matrix);
                    put(car j,'rvalue,'mat . w)>>
            else typerr(list(x,car j),"matrix")
   end;

symbolic procedure natnumlis u;
   % True if U is a list of natural numbers.
   null u
      or numberp car u and fixp car u and car u>0 and natnumlis cdr u;

rlistat '(matrix);

symbolic procedure nzero n;
   % Returns a list of N zeros.
   if n=0 then nil else 0 . nzero(n-1);

% Parsing interface.

symbolic procedure matstat;
   % Read a matrix.
   begin scalar x,y;
   a: scan();
      scan();
      y := xread 'paren;
      if not eqcar(y,'!*comma!*) then y := list y else y := remcomma y;
      x := y . x;
      if nxtsym!* eq '!)
        then return <<scan(); scan(); 'mat . reversip x>>
       else if not(nxtsym!* eq '!,) then symerr("Syntax error",nil);
      go to a
   end;

put('mat,'stat,'matstat);

symbolic procedure formmat(u,vars,mode);
   'list . mkquote 'mat
     . for each x in cdr u collect('list . formlis(x,vars,mode));

put('mat,'formfn,'formmat);

put('mat,'i2d,'mkscalmat);

put('mat,'inversefn,'matinverse);

put('mat,'lnrsolvefn,'lnrsolve);

put('mat,'rtypefn,'(lambda (x) 'matrix));

flag('(mat tp),'matflg);

flag('(mat),'noncommuting);

put('mat,'prifn,'matpri);

flag('(mat),'struct);      % for parsing

put('matrix,'fn,'matflg);

put('matrix,'evfn,'matsm!*);

flag('(matrix),'sprifn);

put('matrix,'tag,'mat);

put('matrix,'lengthfn,'matlength);

put('matrix,'getelemfn,'getmatelem);

put('matrix,'setelemfn,'setmatelem);

symbolic procedure mkscalmat u;
   % Converts id u to 1 by 1 matrix.
   list('mat,list u);

symbolic procedure getmatelem u;
   begin scalar x;
      x := get(car u,'rvalue);
      if not eqcar(x,'mat) then rederr list("Matrix",car u,"not set")
       else if not numlis (u := revlis cdr u) or length u neq 2
        then errpri2(x . u,t);
      return nth(nth(cdr x,car u),cadr u)
   end;

symbolic procedure setmatelem(u,v); letmtr(u,v,get(car u,'rvalue));

symbolic procedure matlength u;
   if not eqcar(u,'mat) then rederr list("Matrix",u,"not set")
    else list('list,length cdr u,length cadr u);

symbolic procedure matsm!*(u,v);
   % Matrix expression simplification function.
   begin
        u := 'mat . for each j in matsm u collect
                        for each k in j collect mk!*sq2 k;
        !*sub2 := nil;   %since all substitutions done;
        return u
   end;

symbolic procedure mk!*sq2 u;
   begin scalar x;
        x := !*sub2;   %since we need value for each element;
        u := subs2 u;
        !*sub2 := x;
        return mk!*sq u
   end;

symbolic procedure matsm u;
   begin scalar x,y;
      for each j in nssimp(u,'matrix) do 
         <<y := multsm(car j,matsm1 cdr j);
           x := if null x then y else addm(x,y)>>;
      return x
   end;

symbolic procedure matsm1 u;
   %returns matrix canonical form for matrix symbol product U;
   begin scalar x,y,z; integer n;
    a:  if null u then return z
         else if eqcar(car u,'!*div) then go to d
         else if atom car u then go to er
         else if caar u eq 'mat then go to c1
         else x := apply(caar u,cdar u);
    b:  z := if null z then x
              else if null cdr z and null cdar z then multsm(caar z,x)
              else multm(x,z);
    c:  u := cdr u;
        go to a;
    c1: if not lchk cdar u then rederr "Matrix mismatch";
        x := for each j in cdar u collect
                for each k in j collect xsimp k;
        go to b;
    d:  y := matsm cadar u;
        if (n := length car y) neq length y
          then rederr "Non square matrix"
         else if (z and n neq length z) then rederr "Matrix mismatch"
         else if cddar u then go to h
         else if null cdr y and null cdar y then go to e;
        x := subfg!*;
        subfg!* := nil;
        if null z then z := apply1(get('mat,'inversefn),y)
         else if null(x := get('mat,'lnrsolvefn))
          then z := multm(apply1(get('mat,'inversefn),y),z)
         else z := apply2(get('mat,'lnrsolvefn),y,z);
        subfg!* := x;
        % Make sure there are no power substitutions.
        z := for each j in z collect for each k in j collect
                 <<!*sub2 := t; subs2 k>>;
        go to c;
    e:  if null caaar y then rederr "Zero denominator";
        y := revpr caar y;
        z := if null z then list list y else multsm(y,z);
        go to c;
     h: if null z then z := generateident n;
        go  to c;
    er: rederr list("Matrix",car u,"not set")
   end;

symbolic procedure lchk u;
   begin integer n;
        if null u or atom car u then return nil;
        n := length car u;
        repeat u := cdr u
           until null u or atom car u or length car u neq n;
        return null u
   end;

symbolic procedure addm(u,v);
   %returns sum of two matrix canonical forms U and V;
   for each j in addm1(u,v,function cons)
      collect addm1(car j,cdr j,function addsq);

symbolic procedure addm1(u,v,w);
   if null u and null v then nil
    else if null u or null v then rederr "Matrix mismatch"
    else apply(w,list(car u,car v)) . addm1(cdr u,cdr v,w);

symbolic procedure tp u; tp1 matsm u;

put('tp,'rtypefn,'getrtypecar);

symbolic procedure tp1 u;
   %returns transpose of the matrix canonical form U;
   %U is destroyed in the process;
   begin scalar v,w,x,y,z;
        v := w := list nil;
        while car u do
         <<x := u;
           y := z := list nil;
           while x do
             <<z := cdr rplacd(z,list caar x);
               x := cdr rplaca(x,cdar x)>>;
           w := cdr rplacd(w,list cdr y)>>;
        return cdr v
   end;

symbolic procedure scalprod(u,v);
   %returns scalar product of two lists (vectors) U and V;
   if null u and null v then nil ./ 1
    else if null u or null v then rederr "Matrix mismatch"
    else addsq(multsq(car u,car v),scalprod(cdr u,cdr v));

symbolic procedure multm(u,v);
   %returns matrix product of two matrix canonical forms U and V;
    (lambda x;
        for each y in u collect for each k in x collect scalprod(y,k))
     tp1 v;

symbolic procedure multsm(u,v);
   %returns product of standard quotient U and matrix standard form V;
   if u = (1 ./ 1) then v
    else for each j in v collect for each k in j collect multsq(u,k);

symbolic procedure letmtr(u,v,y);
   %substitution for matrix elements;
   begin scalar z;
        if not eqcar(y,'mat) then rederr list("Matrix",car u,"not set")
         else if not numlis (z := revlis cdr u) or length z neq 2
          then return errpri2(u,'hold);
        rplaca(pnth(nth(cdr y,car z),cadr z),v);
   end;

endmodule;


module matpri;   % Matrix printing routines.

% Author: Anthony C. Hearn;

global '(!*nat);

symbolic procedure setmatpri(u,v);
   matpri1(cdr v,u);

put('mat,'setprifn,'setmatpri);

symbolic procedure matpri u;
   matpri1(cdr u,"MAT");

symbolic procedure matpri1(u,x);
   %prints a matrix canonical form U with name X;
   begin scalar m,n;
        m := 1;
        for each y in u do
         <<n := 1;
           for each z in y do
              <<varpri(z,list('setq,list(x,m,n),z),'only); n := n+1>>;
        m := m+1>>
   end;

endmodule;


module bareiss; % Inversion routines using the Bareiss 2-step method.

% Author: Anthony C. Hearn;

% This module is rife with essential references to RPLAC-based
% functions.

fluid '(!*exp asymplis!*);

global '(wtl!*);

symbolic procedure matinverse u;
   lnrsolve(u,generateident length u);

symbolic procedure lnrsolve(u,v);
   %U is a matrix standard form, V a compatible matrix form.
   %Value is U**(-1)*V.
   begin integer n; scalar !*exp,temp;
        !*exp := t; n := length u;
        if asymplis!* or wtl!*
          then <<temp := asymplis!* . wtl!*;
                 asymplis!* := wtl!* := nil>>;
        v := backsub(bareiss car normmat augment(u,v),n);
        if temp then <<asymplis!* := car temp; wtl!* := cdr temp>>;
        u := rhside(car v,n);
        v := cdr v;
        return if temp
                 then for each j in u collect
                          for each k in j collect resimp(k ./ v)
                else for each j in u collect
                          for each k in j collect cancel(k ./ v)
   end;

symbolic procedure augment(u,v);
   if null u then nil else append(car u,car v) . augment(cdr u,cdr v);

symbolic procedure generateident n;
  %returns matrix canonical form of identity matrix of order N.
   begin scalar u,v;
        for i := 1:n do
         <<u := nil;
           for j := 1:n do u := ((if i=j then 1 else nil) . 1) . u;
           v := u . v>>;
        return v
   end;

symbolic procedure rhside(u,m);
   if null u then nil else pnth(car u,m+1) . rhside(cdr u,m);

symbolic procedure bareiss u;
  %The 2-step integer preserving elimination method of Bareiss
  %based on the implementation of Lipson.
  %If the value of procedure is NIL then U is singular, otherwise the
  %value is the triangularized form of U (in matrix polynomial form).
  begin scalar aa,c0,ci1,ci2,ik1,ij,kk1,kj,k1j,k1k1,ui,u1,x;
        integer k,k1;
        %U1 points to K-1th row of U
        %UI points to Ith row of U
        %IJ points to U(I,J)
        %K1J points to U(K-1,J)
        %KJ points to U(K,J)
        %IK1 points to U(I,K-1)
        %KK1 points to U(K,K-1)
        %K1K1 points to U(K-1,K-1)
        %M in comments is number of rows in U
        %N in comments is number of columns in U.
        aa:= 1;
        k:= 2;
        k1:=1;
        u1:=u;
        go to pivot;
   agn: u1 := cdr u1;
        if null cdr u1 or null cddr u1 then return u;
        aa:=nth(car u1,k);              %aa := u(k,k).
        k:=k+2;
        k1:=k-1;
        u1:=cdr u1;
   pivot:  %pivot algorithm.
        k1j:= k1k1 := pnth(car u1,k1);
        if car k1k1 then go to l2;
        ui:= cdr u1;                    %i := k.
   l:   if null ui then return nil
         else if null car(ij := pnth(car ui,k1))
          then go to l1;
   l0:  if null ij then go to l2;
        x:= car ij;
        rplaca(ij,negf car k1j);
        rplaca(k1j,x);
        ij:= cdr ij;
        k1j:= cdr k1j;
        go to l0;
   l1:  ui:= cdr ui;
        go to l;
   l2:  ui:= cdr u1;                    %i:= k;
   l21: if null ui then return; %if i>m then return;
        ij:= pnth(car ui,k1);
        c0:= addf(multf(car k1k1,cadr ij),
                    multf(cadr k1k1,negf car ij));
        if c0 then go to l3;
        ui:= cdr ui;                    %i:= i+1;
        go to l21;
   l3:  c0:= quotf!*(c0,aa);
        kk1 := kj := pnth(cadr u1,k1);  %kk1 := u(k,k-1);
        if cdr u1 and null cddr u1 then go to ev0
         else if ui eq cdr u1 then go to comp;
   l31: if null ij then go to comp;     %if i>n then go to comp;
        x:= car ij;
        rplaca(ij,negf car kj);
        rplaca(kj,x);
        ij:= cdr ij;
        kj:= cdr kj;
        go to l31;
        %pivoting complete.
    comp:
        if null cdr u1 then go to ev;
        ui:= cddr u1;                   %i:= k+1;
    comp1:
        if null ui then go to ev;       %if i>m then go to ev;
        ik1:= pnth(car ui,k1);
        ci1:= quotf!*(addf(multf(cadr k1k1,car ik1),
                           multf(car k1k1,negf cadr ik1)),
                     aa);
        ci2:= quotf!*(addf(multf(car kk1,cadr ik1),
                           multf(cadr kk1,negf car ik1)),
                     aa);
        if null cddr k1k1 then go to comp3;%if j>n then go to comp3;
        ij:= cddr ik1;                     %j:= k+1;
        kj:= cddr kk1;
        k1j:= cddr k1k1;
    comp2:
        if null ij then go to comp3;
        rplaca(ij,quotf!*(addf(multf(car ij,c0),
                               addf(multf(car kj,ci1),
                                  multf(car k1j,ci2))),
                     aa));
        ij:= cdr ij;
        kj:= cdr kj;
        k1j:= cdr k1j;
        go to comp2;
    comp3:
        ui:= cdr ui;
        go to comp1;
    ev0:if null c0 then return;
    ev: kj := cdr kk1;
        x := cddr k1k1;                 %x := u(k-1,k+1);
        rplaca(kj,c0);
    ev1:kj:= cdr kj;
        if null kj then go to agn;
        rplaca(kj,quotf!*(addf(multf(car k1k1,car kj),
                               multf(car kk1,negf car x)),
                     aa));
        x := cdr x;
        go to ev1
   end;

symbolic procedure backsub(u,m);
   begin scalar detm,det1,ij,ijj,ri,summ,uj,ur; integer i,jj;
   %N in comments is number of columns in U.
        if null u then rederr "Singular matrix";
        ur := reverse u;
        detm := car pnth(car ur,m);             %detm := u(i,j).
        if null detm then rederr "Singular matrix";
        i := m;
    rows:
        i := i-1;
        ur := cdr ur;
        if null ur then return u . detm;
                %if i=0 then return u . detm.
        ri := car ur;
        jj := m+1;
        ijj:=pnth(ri,jj);
    r2: if null ijj then go to rows;    %if jj>n then go to rows;
        ij := pnth(ri,i);               %j := i;
        det1 := car ij;                 %det1 := u(i,i);
        uj := pnth(u,i);
        summ := nil;                    %summ := 0;
    r3: uj := cdr uj;                   %j := j+1;
        if null uj then go to r4;       %if j>m then go to r4;
        ij := cdr ij;
        summ := addf(summ,multf(car ij,nth(car uj,jj)));
                %summ:=summ+u(i,j)*u(j,jj);
        go to r3;
    r4: rplaca(ijj,quotf!*(addf(multf(detm,car ijj),negf summ),det1));
                %u(i,j):=(detm*u(i,j)-summ)/det1;
        jj := jj+1;
        ijj := cdr ijj;
        go to r2
   end;

symbolic procedure normmat u; 
   %U is a matrix standard form.
   %Value is dotted pair of matrix polynomial form and factor.
   begin scalar x,y,z; 
      x := 1; 
      for each v in u do
         <<y := 1; 
           for each w in v do y := lcm(y,denr w);
           z := (for each w in v
                    collect multf(numr w,quotf(y,denr w)))
              . z; 
           x := multf(y,x)>>; 
      return reverse z . x
   end;

endmodule;


module det;   % Determinant and trace routines.

% Author: Anthony C. Hearn;

symbolic procedure simpdet u; detq matsm carx(u,'det);

% The hashing and determinant routines below are due to M. L. Griss.

comment Some general purpose hashing functions;

flag('(array),'eval);      %declared again for bootstrapping purposes;

array !$hash 64;  %general array for hashing;

symbolic procedure gethash key;
   %access previously saved element;
   assoc(key,!$hash(remainder(key,64)));

symbolic procedure puthash(key,valu);
   begin integer k; scalar buk;
      k := remainder(key,64);
      buk := (key . valu) . !$hash k;
      !$hash k := buk;
      return car buk
   end;

symbolic procedure clrhash;
   for i := 0:64 do !$hash i := nil;

comment Determinant Routines;

symbolic procedure detq u;
   %top level determinant function;
   begin integer len;
      len := length u;   %number of rows;
      for each x in u do
        if length x neq len then rederr "NON SQUARE MATRIX";
      if len=1 then return caar u;
      clrhash();
      u := detq1(u,len,0);
      clrhash();
      return u
   end;

symbolic procedure detq1(u,len,ignnum);
   %U is a square matrix of order LEN. Value is the determinant of U;
   %Algorithm is expansion by minors of first row;
   %IGNNUM is packed set of column indices to avoid;
   begin integer n2; scalar row,sign,z;
      row := car u;   %current row;
      n2 := 1;
      if len=1
        then return <<while twomem(n2,ignnum)
                         do <<n2 := 2*n2; row := cdr row>>;
                      car row>>;   %last row, single element;
      if z := gethash ignnum then return cdr z;
      len := len-1;
      u := cdr u;
      z := nil ./ 1;
      for each x in row do
        <<if not twomem(n2,ignnum)
            then <<if numr x
                     then <<if sign then x := negsq x;
                            z:= addsq(multsq(x,detq1(u,len,n2+ignnum)),
                                        z)>>;
                   sign := not sign>>;
          n2 := 2*n2>>;
      puthash(ignnum,z);
      return z
   end;

symbolic procedure twomem(n1,n2);
   %for efficiency reasons, this procedure should be coded in assembly
   %language;
   not evenp(n2/n1);

put('det,'simpfn,'simpdet);

symbolic procedure simptrace u;
   begin integer n; scalar z;
        u := matsm carx(u,'trace);
        if length u neq length car u then rederr "NON SQUARE MATRIX";
        n := 1;
        z := nil ./ 1;
        for each x in u do <<z := addsq(nth(x,n),z); n := n+1>>;
        return z
   end;

put('trace,'simpfn,'simptrace);

endmodule;


module glmat; % Routines for inverting matrices and finding eigen-values
              % and vectors. Methods are the same as in glsolve module.
 
% Author: Eberhard Schruefer.
 
fluid '(!*cramer !*gcd 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,'(lambda (x) 'matrix));

flag('(minv),'matflg);
 
%put('mateigen,'rtypefn,'(lambda (x) 'matrix));
remprop('mateigen,'rtypefn);
%flag('(mateigen),'matflg);
remflag('(mateigen),'matflg);
 
put('detex,'simpfn,'detex);
 
symbolic procedure matinv u;
   %u is a matrix form. Result is the inverse of matrix u.
   begin scalar sgn,x,y,z;
     integer l,m,lm;
     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 rederr "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 reverse 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).
%*** needs still checking; seems to do fairly well.
  begin scalar arbvars,exu,sgn,q,r,s,x,y,z,eivec;
     integer l,m,lm;
     z := 1;
     if not(getrtype u eq 'matrix) then typerr(u,"matrix");
     u := matsm u;
     lm := length car u;
     exu := 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 l=m then j := addsq(j,negsq !*k2q !*a2k eival);
            if numr j then
            x := list m .* multf(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>>;
     r := if minusf lc z then negf lc z else lc z;
     r := numr subs2(r ./ 1);
     kord!* := eival . kord!*;
     if domainp r then s := 1 else
     s := comfac!-to!-poly comfac(r := reorder r);
     r := quotf1(r,s);
     r := if domainp r then nil else sqfrf r;
     if null domainp s and (mvar s eq eival) then
     if red s then r := (s . 1) . r
     else r := (!*k2f eival . ldeg s) . r;
     for each j in q do r := (absf reorder car j . cdr j) . r;
     kord!* := cdr kord!*;
     r := for each j in r collect reorder car j . cdr j;
     l := length r;
     return 'list .
       for each j in r 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)>>
   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);
 
symbolic procedure reduce!-eival!-powers(v,u);
   if domainp u then u ./ 1
    else if mvar u eq caar v then reduce!-eival!-powers1(v,u ./ 1)
    else if ordop(caar v,mvar u) then u ./ 1
    else addsq(multpq(lpow u,reduce!-eival!-powers(v,lc u)),
           reduce!-eival!-powers(v,red u));
 
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));
 
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;
 
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
 
 
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 multf(lc u,lc v)
                   else multf(lc u,lc v))
              .+ 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!: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;