Artifact 3a4079024acf18be0346c36f48f904e87ad2cc5aebc0e9e67d9edba29eb6a155:
- File
r30/matr.red
— part of check-in
[eb17ceb7f6]
at
2020-04-21 19:40:01
on branch master
— Add Reduce 3.0 to the historical section of the archive, and some more
files relating to version sof PSL from the early 1980s. Thanks are due to
Paul McJones and Nelson Beebe for these, as well as to all the original
authors.git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/historical@5328 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 13054) [annotate] [blame] [check-ins using] [more...]
%********************************************************************* %********************************************************************* % MATRIX PACKAGE %********************************************************************* %********************************************************************; %Copyright (c) 1983 The Rand Corporation; SYMBOLIC; %********************************************************************* % REQUIRES SIMPLIFICATION FUNCTIONS FOR NON-SCALAR QUANTITIES %********************************************************************; FLUID '(!*EXP !*S!*); %Used in this module; GLOBAL '(SUBFG!* !*SUB2 !*NAT); SYMBOLIC PROCEDURE MATSM!* U; %matrix expression simplification function; BEGIN U := MATSM U; U := IF NULL CDR U AND NULL CDAR U THEN MK!*SQ2 CAAR U ELSE 'MAT . MAPC2(U,FUNCTION MK!*SQ2); !*SUB2 := NIL; %since all substitutions done; RETURN U END; SYMBOLIC PROCEDURE MAPC2(U,V); %this very conservative definition is to allow for systems with %poor handling of functional arguments, and because of bootstrap- %ping difficulties, which are no longer really relevant; BEGIN SCALAR X,Y,Z; A: IF NULL U THEN RETURN REVERSIP Z; X := CAR U; Y := NIL; B: IF NULL X THEN GO TO C; Y := APPLY(V,LIST CAR X) . Y; X := CDR X; GO TO B; C: U := CDR U; Z := REVERSIP Y . Z: GO TO A 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; U := NSSIMP(U,'MATP); A: IF NULL U THEN RETURN X; Y := MULTSM(CAAR U,MATRIXTIMES CDAR U); X := IF NULL X THEN Y ELSE ADDM(X,Y); U := CDR U; GO TO A END; SYMBOLIC PROCEDURE MATRIXTIMES 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 IF (X := GET(CAAR U,'MSIMPFN)) THEN X := APPLY(X,CDAR U) ELSE GO TO ER; 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 := MAPC2(CDAR U,FUNCTION XSIMP); 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:= GENERATEIDENT N; Z := LNRSOLVE(Y,Z); SUBFG!* := X; 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; 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(!*S!*,U); %returns product of standard quotient !*S!* and matrix standard %form U; IF !*S!* = (1 ./ 1) THEN U ELSE MAPC2(U,FUNCTION (LAMBDA J; MULTSQ(!*S!*,J))); 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; SYMBOLIC PROCEDURE MATPRI!*(U,V,W); %symbolic interface to VARPRI; MATPRI(CDR U,IF V THEN EVAL CAR V ELSE NIL); SYMBOLIC PROCEDURE MATPRI(U,X); %prints a matrix canonical form U with name X; BEGIN SCALAR M,N; M := 1; IF NULL X THEN X := 'MAT; FOR EACH Y IN U DO <<N := 1; FOR EACH Z IN Y DO <<VARPRI(Z,LIST MKQUOTE LIST(X,M,N),T); IF !*NAT THEN TERPRI!* T; N := N+1>>; M := M+1>> END; %********************************************************************* % MATRIX INVERSION ROUTINES %********************************************************************; 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 X,!*S!*; X := !*EXP; !*EXP := T; N := LENGTH U; !*S!* := BACKSUB(BAREISS CAR NORMMAT AUGMENT(U,V),N); U := MAPC2(RHSIDE(CAR !*S!*,N), FUNCTION (LAMBDA J; CANCEL(J . CDR !*S!*))); !*EXP := X; RETURN U 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; %********************************************************************* % DETERMINANT AND TRACE ROUTINES %********************************************************************; SYMBOLIC PROCEDURE SIMPDET U; DETQ MATSM CARX(U,'DET); COMMENT 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; REMAINDER(N2/N1,2)=1; 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"; Z := NIL ./ 1; N := 1; A: IF NULL U THEN RETURN Z; Z := ADDSQ(NTH(CAR U,N),Z); U := CDR U; N := N+1; GO TO A END; PUT('TRACE,'SIMPFN,'SIMPTRACE); END;