Artifact be96d4b0f11d530ff46e238ea930a7639dc4d7cebccda3807ca5832a388a4bd1:
- File
r33/matr.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 27388) [annotate] [blame] [check-ins using] [more...]
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;