1
0
mirror of https://github.com/PDP-10/its.git synced 2026-04-18 00:57:09 +00:00
Files
PDP-10.its/src/paulw/mat.286
Eric Swenson 85994ed770 Added files to support building and running Macsyma.
Resolves #284.

Commented out uses of time-origin in maxtul; mcldmp (init) until we
can figure out why it gives arithmetic overflows under the emulators.

Updated the expect script statements in build_macsyma_portion to not
attempt to match expected strings, but simply sleep for some time
since in some cases the matching appears not to work.
2018-03-11 13:10:19 -07:00

453 lines
14 KiB
Common Lisp
Raw Blame History

This file contains invisible Unicode characters
This file contains invisible Unicode characters that are indistinguishable to humans but may be processed differently by a computer. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
;;;;;;;;;;;;;;;;;;; -*- Mode: Lisp; Package: Macsyma -*- ;;;;;;;;;;;;;;;;;;;
;;; (c) Copyright 1982 Massachusetts Institute of Technology ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(macsyma-module mat)
(COMMENT THIS IS THE MAT PACKAGE)
(DECLARE (SPECIAL PIVSIGN* *ECH* *TRI* LSOLVEFLAG $ALGEBRAIC
$MULTIPLICITIES EQUATIONS
MUL* FORMATFORM DOSIMP $DISPFLAG $RATFAC
*TB $NOLABELS ERRRJFFLAG *DET* GENVAR
XM* XN* VARLIST AX LINELABLE $LINECHAR $LINENUM SOL)
(*LEXPR $SOLVE $RAT)
(ARRAY* (NOTYPE XA* 2))
(FIXNUM TIM)
(GENPREFIX MAT))
;; The array declarations of ROW, COL, and COLINV aren't having any
;; effect on the Lisp Machine. Should be fixed somehow.
(DECLARE (SPECIAL *DET*)
(ARRAY* (FIXNUM ROW 1 COL 1 COLINV 1)))
(DEFMVAR $GLOBALSOLVE NIL)
(DEFMVAR $SPARSE NIL)
(DEFMVAR $BACKSUBST T)
(DEFVAR *RANK* NIL)
(DEFVAR *INV* NIL)
(DEFVAR SOLVEXP NIL)
(DEFUN SOLCOEF
(M *C VARL FLAG)
(PROG (CC ANSWER LEFTOVER)
(SETQ CC (CDR (RATREP* *C)))
(IF (OR (ATOM (CAR CC))
(NOT (EQUAL (CDAR CC) '(1 1)))
(NOT (EQUAL 1 (CDR CC))))
(MERROR "Unacceptable variable to SOLVE:~%~M" *C))
(SETQ ANSWER (RATREDUCE (PRODCOEF (CAR CC) (CAR M)) (CDR M)))
(IF (NOT FLAG) (RETURN ANSWER))
(SETQ LEFTOVER
(RDIS (RATPLUS M (RATTIMES (RATMINUS ANSWER) CC T))))
(IF (OR (NOT (FREEOF *C LEFTOVER))
(DEPENDSALL (RDIS ANSWER) VARL))
(ERRRJF "NON-LINEAR"))
(RETURN ANSWER)))
(DEFUN FORMX (FLAG NAM EQL VARL)
(PROG (B AX X IX J)
(SETQ VARLIST VARL)
(MAPC #'NEWVAR EQL)
(AND (NOT $ALGEBRAIC)
(ORMAPC #'ALGP VARLIST)
(SETQ $ALGEBRAIC T))
(*ARRAY NAM T (1+ (SETQ XN* (LENGTH EQL)))
(1+ (SETQ XM* (1+ (LENGTH VARL)))))
(SETQ IX 0)
LOOP1(COND ((NULL EQL) (RETURN VARLIST)))
(SETQ AX (CAR EQL))
(SETQ EQL (CDR EQL))
(SETQ IX (1+ IX))
(STORE (FUNCALL NAM IX XM*) (CONST AX VARL))
(SETQ J 0)
(SETQ B VARL) (SETQ AX (CDR (RATREP* AX)))
LOOP2(SETQ X (CAR B))
(SETQ B (CDR B))
(SETQ J (1+ J))
(STORE (FUNCALL NAM IX J) (SOLCOEF AX X VARL FLAG))
(COND (B (GO LOOP2)))
(GO LOOP1)))
(DEFUN DEPENDSALL
(EXP L)
(COND ((NULL L) NIL)
((OR (NOT (FREEOF (CAR L) EXP)) (DEPENDSALL EXP (CDR L))) T)
(T NIL)))
(SETQ *DET* NIL *ECH* NIL *RANK* NIL *INV* NIL *TRI* NIL)
(DEFUN PTORAT (AX M N)
(PROG (I J)
(SETQ AX (GET-ARRAY-POINTER AX))
(SETQ I (1+ M) N (1+ N))
LOOP1
(COND ((EQUAL I 1) (RETURN NIL)))
(SETQ I (1- I) J N)
LOOP2
(COND ((EQUAL J 1) (GO LOOP1)))
(SETQ J (1- J))
(STORE (ARRAYCALL T AX I J) (CONS (ARRAYCALL T AX I J) 1))
(GO LOOP2)
))
(DEFUN MEQHK (Z)
(COND ((AND (NOT (ATOM Z)) (EQ (CAAR Z) 'MEQUAL))
(SIMPLUS (LIST '(MPLUS) (CADR Z) (LIST '(MTIMES) -1 (CADDR Z))) 1 NIL))
(T Z)))
(DEFUN CONST (E VARL)
(PROG (ZL)
(SETQ VARL (MAPCAR (FUNCTION (LAMBDA(X) (CAADR (RATREP* X)))) VARL))
(SETQ E (CDR(RATREP* E)))
(SETQ ZL (NZEROS (LENGTH VARL) NIL))
(RETURN (RATREDUCE (PCTIMES -1 (PCSUBSTY ZL VARL (CAR E)))
(PCSUBSTY ZL VARL (CDR E))))))
(DEFVAR *MOSESFLAG NIL)
(DEFMVAR $%RNUM 0)
(DEFMFUN MAKE-PARAM ()
(LET ((PARAM (CONCAT '$%R (SETQ $%RNUM (1+ $%RNUM)))))
(TUCHUS $%RNUM_LIST PARAM)
PARAM))
(DEFMVAR $LINSOLVE_PARAMS T "LINSOLVE generates %Rnums")
(DECLARE (FIXNUM N))
(DEFUN NCDR (X N) (NTHCDR (1- N) X))
(DEFUN ITH (X N) (COND ((ATOM X) NIL) (T (CAR (NCDR X N)))))
(DECLARE (NOTYPE N))
(DEFUN POLYIZE (AX R M MUL)
(DECLARE (FIXNUM M C))
(DO ((C 1 (1+ C)) (D))
((> C M) NIL)
(SETQ D (ARRAYCALL T AX R C))
(SETQ D (COND ((EQUAL MUL 1) (CAR D))
(T (PTIMES (CAR D)
(PQUOTIENTCHK MUL (CDR D))))))
(STORE (ARRAYCALL T AX R C) (IF $SPARSE (CONS D 1) D))))
; TWO-STEP FRACTION-FREE GAUSSIAN ELIMINATION ROUTINE
(DEFUN TFGELI (AX N M &AUX ($SPARSE (AND $SPARSE (OR *DET* *INV*))))
;;$sparse is also controlling whether polyize stores polys or ratforms
(SETQ AX (GET-ARRAY-POINTER AX))
(SETQ MUL* 1)
(DO ((R 1 (1+ R)))
((> R N) (COND ((AND $SPARSE *DET*)(SPRDET AX N))
((AND *INV* $SPARSE)(NEWINV AX N M))
(T (TFGELI1 AX N M))))
(DO ((C 1 (1+ C))
(D)
(MUL 1))
((> C M)
(AND *DET* (SETQ MUL* (PTIMES MUL* MUL)))
(POLYIZE AX R M MUL))
(COND ((EQUAL 1 (SETQ D (CDR (ARRAYCALL T AX R C)))) NIL)
(T (SETQ MUL (PTIMES MUL (PQUOTIENT D (PGCD MUL D)))))))))
(SETQ LSOLVEFLAG NIL)
;; The author of the following programs is Tadatoshi Minamikawa (TM).
;; This program is one-step fraction-free Gaussian elimination with
;; optimal pivotting. DRB claims the hair in this program is not
;; necessary and that straightforward Gaussian elimination is sufficient,
;; for sake of future implementors.
;; To debug, delete the comments around PRINT and BREAK statements.
(DECLARE (SPECIAL PERMSIGN A RANK DELTA NROW NVAR N M VARIABLEORDER
DEPENDENTROWS INCONSISTENTROWS L K)
;; We could just use fortran, you know.
(FIXNUM NROW NVAR RANK I J K L M N))
(DEFUN TFGELI1 (AX N M)
(PROG (K L DELTA VARIABLEORDER INCONSISTENTROWS
DEPENDENTROWS NROW NVAR RANK PERMSIGN RESULT)
(*ARRAY 'ROW 'FIXNUM (1+ N))
(*ARRAY 'COL 'FIXNUM (1+ M))
(*ARRAY 'COLINV 'FIXNUM (1+ M))
#+LISPM (FILLARRAY (FUNCTION ROW) '(0))
#+LISPM (FILLARRAY (FUNCTION COL) '(0))
#+LISPM (FILLARRAY (FUNCTION COLINV) '(0))
(SETQ AX (GET-ARRAY-POINTER AX))
;; (PRINT 'ONESTEP-LIPSON-WITH-PIVOTTING)
(SETQ NROW N)
(SETQ NVAR (COND (*RANK* M) (*DET* M) (*INV* N) (*ECH* M) (*TRI* M) (T (1- M))))
(DO I 1 (1+ I) (> I N) (STORE (ROW I) I))
(DO I 1 (1+ I) (> I M)
(STORE (COL I) I) (STORE (COLINV I) I))
(SETQ RESULT
(COND
(*RANK* (FORWARD T) RANK)
(*DET* (FORWARD T)
(COND ((= NROW N) (COND (PERMSIGN (PMINUS DELTA))
(T DELTA)))
(T 0)))
(*INV* (FORWARD T) (BACKWARD) (RECOVERORDER1))
(*ECH* (FORWARD NIL) (RECOVERORDER2))
(*TRI* (FORWARD NIL) (RECOVERORDER2))
(T (FORWARD T) (COND ($BACKSUBST (BACKWARD)))
(RECOVERORDER2)
(LIST DEPENDENTROWS INCONSISTENTROWS VARIABLEORDER))))
(*REARRAY 'ROW) (*REARRAY 'COL) (*REARRAY 'COLINV)
(RETURN RESULT)))
;FORWARD ELIMINATION
;IF THE SWITCH *CPIVOT IS NIL, IT AVOIDS THE COLUMN PIVOTTING.
(DEFUN FORWARD (*CPIVOT)
(SETQ DELTA 1) ;DELTA HOLDS THE CURRENT DETERMINANT
(DO ((K 1 (1+ K))
(NVAR NVAR) ;PROTECTS AGAINST TEMPORARAY RESETS DONE IN PIVOT
(M M))
((OR (> K NROW) (> K NVAR)))
(COND ((PIVOT AX K *CPIVOT) (RETURN NIL)))
;PIVOT IS T IF THERE IS NO MORE NON-ZERO ROW LEFT. THEN GET OUT OF THE LOOP
(DO ((I (1+ K) (1+ I)))
((> I NROW))
(DO ((J (1+ K) (1+ J)))
((> J M))
(STORE ( ARRAYCALL T AX (ROW I) (COL J))
(PQUOTIENT (PDIFFERENCE (PTIMES ( ARRAYCALL T AX (ROW K) (COL K))
(ARRAYCALL T AX (ROW I) (COL J)))
(PTIMES ( ARRAYCALL T AX (ROW I) (COL K))
( ARRAYCALL T AX (ROW K) (COL J))))
DELTA))))
(DO ((I (1+ K) (1+ I)))
((> I NROW))
(STORE ( ARRAYCALL T AX (ROW I) (COL K))
0))
(SETQ DELTA ( ARRAYCALL T AX (ROW K) (COL K))))
; UNDOES COLUMN HACK IN PIVOT.
(OR *CPIVOT (DO I 1 (1+ I) (> I M) (STORE (COL I) I)))
(SETQ RANK (MIN NROW NVAR)))
; BACKWARD SUBSTITUTION
(DEFUN BACKWARD ()
(DO ((I (1- RANK) (1- I)))
((< I 1))
(DO ((L (1+ RANK) (1+ L)))
((> L M))
(STORE (ARRAYCALL T AX (ROW I) (COL L)) (PQUOTIENT (PDIFFERENCE (PTIMES (ARRAYCALL T AX (ROW I) (COL L)) (ARRAYCALL T AX (ROW RANK) (COL RANK)) ) (DO ((J (1+ I) (1+ J)) (SUM 0))
((> J RANK) SUM)
(SETQ SUM (PPLUS SUM (PTIMES (ARRAYCALL T AX (ROW I) (COL J)) (ARRAYCALL T AX (ROW J) (COL L))))))) (ARRAYCALL T AX (ROW I) (COL I)))))
(DO ((L (1+ I) (1+ L)))
((> L RANK))
(STORE (ARRAYCALL T AX (ROW I) (COL L)) 0)))
;PUT DELTA INTO THE DIAGONAL MATRIX
(SETQ DELTA (ARRAYCALL T AX (ROW RANK) (COL RANK)))
(DO ((I 1 (1+ I)))
((> I RANK))
(STORE (ARRAYCALL T AX (ROW I) (COL I)) DELTA)))
;RECOVER THE ORDER OF ROWS AND COLUMNS.
(DEFUN RECOVERORDER1 ()
;(COMMENT (PRINT 'REARRANGE))
(DO ((I NVAR (1- I)))
((= I 0))
(SETQ VARIABLEORDER (CONS I VARIABLEORDER)))
(DO ((I (1+ RANK) (1+ I)))
((> I N))
(COND ((EQUAL (ARRAYCALL T AX (ROW I) (COL M)) 0)
(SETQ DEPENDENTROWS (CONS (ROW I) DEPENDENTROWS)))
(T (SETQ INCONSISTENTROWS (CONS (ROW I) INCONSISTENTROWS)))))
(DO ((I 1 (1+ I)))
((> I N))
(COND ((NOT (= (ROW (COLINV I)) I)) (PROG () (MOVEROW AX N M I 0) (SETQ L I) LOOP (SETQ K (ROW (COLINV L))) (STORE (ROW (COLINV L)) L) (COND ((= K I) (MOVEROW AX N M 0 L)) (T (MOVEROW AX N M K L) (SETQ L K) (GO LOOP))))))))
(DEFUN RECOVERORDER2 ()
(DO ((I NVAR (1- I)))
((= I 0))
(SETQ VARIABLEORDER (CONS (COL I) VARIABLEORDER)))
(DO ((I (1+ RANK) (1+ I)))
((> I N))
(COND ((EQUAL (ARRAYCALL T AX (ROW I) (COL M)) 0) (SETQ DEPENDENTROWS (CONS (ROW I) DEPENDENTROWS)))
(T (SETQ INCONSISTENTROWS (CONS (ROW I) INCONSISTENTROWS)))))
(DO ((I 1 (1+ I)))
((> I N))
(COND ((NOT (= (ROW I) I))
(PROG ()
(MOVEROW AX N M I 0)
(SETQ L I)
LOOP (SETQ K (ROW L))
(STORE (ROW L) L)
(COND ((= K I) (MOVEROW AX N M 0 L))
(T (MOVEROW AX N M K L)
(SETQ L K)
(GO LOOP)))))))
(DO ((I 1 (1+ I)))
((> I NVAR))
(COND ((NOT (= (COL I) I))
(PROG ()
(MOVECOL AX N M I 0)
(SETQ L I)
LOOP2 (SETQ K (COL L))
(STORE (COL L) L)
(COND ((= K I) (MOVECOL AX N M 0 L))
(T (MOVECOL AX N M K L)
(SETQ L K)
(GO LOOP2))))))))
;THIS PROGRAM IS USED IN REARRANGEMENT
(DEFUN MOVEROW (AX N M I J)
(DO ((K 1 (1+ K))) ((> K M))
(STORE (ARRAYCALL T AX J K) (ARRAYCALL T AX I K))))
(DEFUN MOVECOL (AX N M I J)
(DO ((K 1 (1+ K))) ((> K N))
(STORE (ARRAYCALL T AX K J) (ARRAYCALL T AX K I))))
;COMPLEXITY IS DEFINED AS FOLLOWS
; COMPLEXITY(0)=0
; COMPLEXITY(CONSTANT)=1
; COMPLEXITY(POLYNOMIAL)=1+SUM(COMPLEXITY(C(N))+COMPLEXITY(E(N)), FOR N=0,1 ...M)
; WHERE POLYNOMIAL IS OF THE FORM
; SUM(C(N)*X^E(N), FOR N=0,1 ... M) X IS THE VARIABLE
(DEFUN COMPLEXITY (EXP)
(COND ((NULL EXP) 0) ((EQUAL EXP 0) 0) ((ATOM EXP) 1)
(T (PLUS (COMPLEXITY (CAR EXP)) (COMPLEXITY (CDR EXP))
))))
(DEFUN COMPLEXITY/ROW (AX I J1 J2)
(DO ((J J1 (1+ J)) (SUM 0))
((> J J2) SUM)
(SETQ SUM (PLUS SUM (COMPLEXITY (ARRAYCALL T AX (ROW I) (COL J)))))))
(DEFUN COMPLEXITY/COL (AX J I1 I2)
(DO ((I I1 (1+ I)) (SUM 0))
((> I I2) SUM)
(SETQ SUM (PLUS SUM (COMPLEXITY (ARRAYCALL T AX (ROW I) (COL J)))))))
(DEFUN ZEROP/ROW (AX I J1 J2)
(DO ((J J1 (1+ J)))
((> J J2) T)
(COND ((NOT (EQUAL (ARRAYCALL T AX (ROW I) (COL J)) 0)) (RETURN NIL)))))
;PIVOTTING ALGORITHM
(DEFUN PIVOT (AX K *CPIVOT)
(PROG (ROW/OPTIMAL COL/OPTIMAL COMPLEXITY/I/MIN COMPLEXITY/J/MIN
COMPLEXITY/I COMPLEXITY/J COMPLEXITY/DET COMPLEXITY/DET/MIN DUMMY)
(SETQ ROW/OPTIMAL K COMPLEXITY/I/MIN 1000000. COMPLEXITY/J/MIN 1000000.)
;TEST THE SINGULARITY
(COND ((DO ((I K (1+ I)) (ISALLZERO T))
((> I NROW) ISALLZERO)
LOOP (COND ((ZEROP/ROW AX I K NVAR)
(COND (*INV* (MERROR "Singular"))
(T (EXCHANGEROW I NROW)
(SETQ NROW (1- NROW))))
(COND ((NOT (> I NROW)) (GO LOOP))))
(T (SETQ ISALLZERO NIL))))
(RETURN T)))
;FIND AN OPTIMAL ROW
;IF *CPIVOT IS NIL , (AX I K) WHICH IS TO BE THE PIVOT MUST BE NONZERO. BUT IF *CPIVOT IS T, IT IS UNNECESSARY BECAUSE WE CAN DO THE COLUMN PIVOT.
FINDROW
(DO ((I K (1+ I)))
((> I NROW))
(COND ((OR *CPIVOT (NOT (EQUAL (ARRAYCALL T AX (ROW I) (COL K)) 0)))
(COND ((> COMPLEXITY/I/MIN
(SETQ COMPLEXITY/I (COMPLEXITY/ROW AX I K M)))
(SETQ ROW/OPTIMAL I COMPLEXITY/I/MIN COMPLEXITY/I))))))
;EXCHANGE THE ROWS K AND ROW/OPTIMAL
(EXCHANGEROW K ROW/OPTIMAL)
;IF THE FLAG *CPIVOT IS NIL, THE FOLLOWING STEPS ARE NOT EXECUTED. THIS TREATMENT WAS DONE FOR THE LSA AND ECHELONTHINGS WHICH ARE NOT HAPPY WITH THE COLUMN OPERATIONS.
(COND ((NULL *CPIVOT)
(COND ((NOT (EQUAL (ARRAYCALL T AX (ROW K) (COL K)) 0))
(RETURN NIL))
(T (DO I K (1+ I) (= I NVAR)
(STORE (COL I) (COL (1+ I))))
(SETQ NVAR (1- NVAR) M (1- M))
(GO FINDROW)))))
;STEP3 ...FIND THE OPTIMAL COLUMN
(SETQ COL/OPTIMAL 0
COMPLEXITY/DET/MIN 1000000.
COMPLEXITY/J/MIN 1000000.)
(DO ((J K (1+ J)))
((> J NVAR))
(COND ((NOT (EQUAL (ARRAYCALL T AX (ROW K) (COL J)) 0))
(COND ((> COMPLEXITY/DET/MIN
(SETQ COMPLEXITY/DET
(COMPLEXITY (ARRAYCALL T AX (ROW K) (COL J)))))
(SETQ COL/OPTIMAL J
COMPLEXITY/DET/MIN COMPLEXITY/DET
COMPLEXITY/J/MIN (COMPLEXITY/COL AX J (1+ K) N)))
((EQUAL COMPLEXITY/DET/MIN COMPLEXITY/DET)
(COND ((> COMPLEXITY/J/MIN
(SETQ COMPLEXITY/J
(COMPLEXITY/COL AX J (1+ K) N)))
(SETQ COL/OPTIMAL J
COMPLEXITY/DET/MIN COMPLEXITY/DET
COMPLEXITY/J/MIN COMPLEXITY/J))))))))
;(COND ((ZEROP COL/OPTIMAL) (COMMENT (PRINT '"SINGULAR!"))))
;EXCHANGE THE COLUMNS K AND COL/OPTIMAL
(EXCHANGECOL K COL/OPTIMAL)
(SETQ DUMMY (COLINV (COL K)))
(STORE (COLINV (COL K)) (COLINV (COL COL/OPTIMAL)))
(STORE (COLINV (COL COL/OPTIMAL)) DUMMY)
(RETURN NIL)))
(DEFUN EXCHANGEROW (I J)
(PROG (DUMMY)
(SETQ DUMMY (ROW I))
(STORE (ROW I) (ROW J))
(STORE (ROW J) DUMMY)
(COND ((= I J) (RETURN NIL))
(T (SETQ PERMSIGN (NOT PERMSIGN))))))
(DEFUN EXCHANGECOL (I J)
(PROG (DUMMY)
(SETQ DUMMY (COL I))
(STORE (COL I) (COL J))
(STORE (COL J) DUMMY)
(COND ((= I J) (RETURN NIL))
(T (SETQ PERMSIGN (NOT PERMSIGN))))))
;; Displays list of solutions.
(DEFUN SOLVE2 (List)
(SETQ $MULTIPLICITIES NIL)
(MAP2C #'(LAMBDA (EQUATN MULTIPL)
(SETQ EQUATIONS
(NCONC EQUATIONS (LIST (DISPLINE EQUATN))))
(PUSH MULTIPL $MULTIPLICITIES)
(IF (AND (> MULTIPL 1) $DISPFLAG)
(MTELL "Multiplicity ~A~%" MULTIPL)))
List)
(SETQ $MULTIPLICITIES (CONS '(MLIST SIMP) (NREVERSE $MULTIPLICITIES))))
;; Displays an expression and returns its linelabel.
(DEFMFUN DISPLINE (EXP)
(LET ($NOLABELS (TIM 0))
(ELABEL EXP)
(COND ($DISPFLAG (REMPROP LINELABLE 'NODISP)
(SETQ TIM (RUNTIME))
(MTERPRI)
(DISPLA (LIST '(MLABLE) LINELABLE EXP))
(TIMEORG TIM))
(T (PUTPROP LINELABLE T 'NODISP)))
LINELABLE))