1;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;; 2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3;;; The data in this file contains enhancments. ;;;;; 4;;; ;;;;; 5;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;; 6;;; All rights reserved ;;;;; 7;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 8;;; (c) Copyright 1982 Massachusetts Institute of Technology ;;; 9;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 10 11(in-package :maxima) 12 13(macsyma-module mat) 14 15;; this is the mat package 16 17(declare-top (special *ech* *tri* $algebraic $multiplicities equations 18 mul* $dispflag $nolabels *det* 19 xm* xn* varlist ax *linelabel*)) 20 21;;these are arrays. 22(defvar *row*) 23(defvar *col*) 24(defvar *colinv*) 25 26(defmvar $globalsolve nil) 27(defmvar $sparse nil) 28(defmvar $backsubst t) 29 30(defmvar *rank* nil) 31(defmvar *inv* nil) 32 33(defun solcoef (m *c varl flag) 34 (prog (cc answer leftover) 35 (setq cc (cdr (ratrep* *c))) 36 (if (or (atom (car cc)) 37 (not (equal (cdar cc) '(1 1))) 38 (not (equal 1 (cdr cc)))) 39 ;; NOTE TO TRANSLATORS: NOT CLEAR WHAT IS "UNACCEPTABLE" HERE 40 (merror (intl:gettext "solve: unacceptable variable: ~M") *c)) 41 (setq answer (ratreduce (prodcoef (car cc) (car m)) (cdr m))) 42 (if (not flag) (return answer)) 43 (setq leftover 44 (rdis (ratplus m (rattimes (ratminus answer) cc t)))) 45 (if (or (not (freeof *c leftover)) 46 (dependsall (rdis answer) varl)) 47 (rat-error "`non-linear'")) 48 (return answer))) 49 50(defun formx (flag nam eql varl) 51 (prog (b ax x ix j) 52 (setq varlist varl) 53 (mapc #'newvar eql) 54 (and (not $algebraic) 55 (some #'algp varlist) 56 (setq $algebraic t)) 57 (setf (symbol-value nam) (make-array (list (1+ (setq xn* (length eql))) 58 (1+ (setq xm* (1+ (length varl))))))) 59 (setq nam (get-array-pointer nam)) 60 (setq ix 0) 61 loop1 62 (cond ((null eql) (return varlist))) 63 (setq ax (car eql)) 64 (setq eql (cdr eql)) 65 (incf ix) 66 (setf (aref nam ix xm*) (const ax varl)) 67 (setq j 0) 68 (setq b varl) (setq ax (cdr (ratrep* ax))) 69 loop2 70 (setq x (car b)) 71 (setq b (cdr b)) 72 (incf j) 73 (setf (aref nam ix j) (solcoef ax x varl flag)) 74 (cond (b (go loop2))) 75 (go loop1))) 76 77(defun dependsall (exp l) 78 (cond ((null l) nil) 79 ((or (not (freeof (car l) exp)) (dependsall exp (cdr l))) t) 80 (t nil))) 81 82(setq *det* nil *ech* nil *tri* nil) 83 84(defun ptorat (ax m n) 85 (prog (i j) 86 (setq ax (get-array-pointer ax)) 87 (setq i (1+ m)) 88 (incf n) 89 loop1 90 (when (equal i 1) (return nil)) 91 (decf i) 92 (setq j n) 93 loop2 94 (when (equal j 1) (go loop1)) 95 (decf j) 96 (setf (aref ax i j) (cons (aref ax i j) 1)) 97 (go loop2))) 98 99(defun meqhk (z) 100 "If Z is of the form lhs = rhs, return the expression lhs - rhs. 101 Otherwise, return Z unchanged." 102 (if (and (not (atom z)) (eq (caar z) 'mequal)) 103 (simplus (list '(mplus) (cadr z) (list '(mtimes) -1 (caddr z))) 1 nil) 104 z)) 105 106(defun const (e varl) 107 (setq varl (mapcar #'(lambda(x) (caadr (ratrep* x))) varl)) 108 (setq e (cdr (ratrep* e))) 109 (let ((zl (make-list (length varl) :initial-element 0))) 110 (ratreduce (pctimes -1 (pcsubsty zl varl (car e))) 111 (pcsubsty zl varl (cdr e))))) 112 113(defvar *mosesflag nil) 114 115(defmvar $%rnum 0) 116 117(defun make-param () 118 (let ((param (intern (format nil "~A~D" '$%r (incf $%rnum))))) 119 (tuchus $%rnum_list param) 120 param)) 121 122(defmvar $linsolve_params t "`linsolve' generates %Rnums") 123 124(defun ith (x n) 125 (if (atom x) nil (nth (1- n) x))) 126 127(defun polyize (ax r m mul) 128 (declare (fixnum m)) 129 (do ((c 1 (1+ c)) (d)) 130 ((> c m) nil) 131 (declare (fixnum c)) 132 (setq d (aref ax r c)) 133 (setq d (cond ((equal mul 1) (car d)) 134 (t (ptimes (car d) 135 (pquotientchk mul (cdr d)))))) 136 (setf (aref ax r c) (if $sparse (cons d 1) d)))) 137 138;; TWO-STEP FRACTION-FREE GAUSSIAN ELIMINATION ROUTINE 139 140(defun tfgeli (ax n m &aux ($sparse (and $sparse (or *det* *inv*)))) 141 ;;$sparse is also controlling whether polyize stores polys or ratforms 142 (setq ax (get-array-pointer ax)) 143 (setq mul* 1) 144 (do ((r 1 (1+ r))) 145 ((> r n) (cond ((and $sparse *det*)(sprdet ax n)) 146 ((and *inv* $sparse)(newinv ax n m)) 147 (t (tfgeli1 ax n m)))) 148 (do ((c 1 (1+ c)) 149 (d) 150 (mul 1)) 151 ((> c m) 152 (and *det* (setq mul* (ptimes mul* mul))) 153 (polyize ax r m mul)) 154 (cond ((equal 1 (setq d (cdr (aref ax r c)))) nil) 155 (t (setq mul (ptimes mul (pquotient d (pgcd mul d))))))))) 156 157;; The author of the following programs is Tadatoshi Minamikawa (TM). 158;; This program is one-step fraction-free Gaussian elimination with 159;; optimal pivotting. DRB claims the hair in this program is not 160;; necessary and that straightforward Gaussian elimination is sufficient, 161;; for sake of future implementors. 162 163;; To debug, delete the comments around PRINT and BREAK statements. 164 165(declare-top (special permsign a rank delta nrow nvar n m variableorder 166 dependentrows inconsistentrows l k)) 167 168(defun tfgeli1 (ax n m) 169 (prog (k l delta variableorder inconsistentrows 170 dependentrows nrow nvar rank permsign result) 171 (setq ax (get-array-pointer ax)) 172 (setq *col* (make-array (1+ m) :initial-element 0)) 173 (setq *row* (make-array (1+ n) :initial-element 0)) 174 (setq *colinv* (make-array (1+ m) :initial-element 0)) 175 ;; (PRINT 'ONESTEP-LIPSON-WITH-PIVOTTING) 176 (setq nrow n) 177 (setq nvar (cond (*rank* m) (*det* m) (*inv* n) (*ech* m) (*tri* m) (t (1- m)))) 178 (do ((i 1 (1+ i))) 179 ((> i n)) 180 (setf (aref *row* i) i)) 181 (do ((i 1 (1+ i))) 182 ((> i m)) 183 (setf (aref *col* i) i) (setf (aref *colinv* i) i)) 184 (setq result 185 (cond (*rank* (forward t) rank) 186 (*det* (forward t) 187 (cond ((= nrow n) (cond (permsign (pminus delta)) 188 (t delta))) 189 (t 0))) 190 (*inv* (forward t) (backward) (recoverorder1)) 191 (*ech* (forward nil) (recoverorder2)) 192 (*tri* (forward nil) (recoverorder2)) 193 (t (forward t) (cond ($backsubst (backward))) 194 (recoverorder2) 195 (list dependentrows inconsistentrows variableorder)))) 196 (return result))) 197 198;;FORWARD ELIMINATION 199;;IF THE SWITCH *CPIVOT IS NIL, IT AVOIDS THE COLUMN PIVOTTING. 200(defun forward (*cpivot) 201 (setq delta 1) ;DELTA HOLDS THE CURRENT DETERMINANT 202 (do ((k 1 (1+ k)) 203 (nvar nvar) ;PROTECTS AGAINST TEMPORARAY RESETS DONE IN PIVOT 204 (m m)) 205 ((or (> k nrow) (> k nvar))) 206 (cond ((pivot ax k *cpivot) (return nil))) 207 ;; PIVOT IS T IF THERE IS NO MORE NON-ZERO ROW LEFT. THEN GET OUT OF THE LOOP 208 (do ((i (1+ k) (1+ i))) 209 ((> i nrow)) 210 (do ((j (1+ k) (1+ j))) 211 ((> j m)) 212 (setf (aref ax (aref *row* i) (aref *col* j)) 213 (pquotient (pdifference (ptimes (aref ax (aref *row* k) (aref *col* k)) 214 (aref ax (aref *row* i) (aref *col* j))) 215 (ptimes (aref ax (aref *row* i) (aref *col* k)) 216 (aref ax (aref *row* k) (aref *col* j)))) 217 delta)))) 218 (do ((i (1+ k) (1+ i))) 219 ((> i nrow)) 220 (setf (aref ax (aref *row* i) (aref *col* k)) 0)) 221 (setq delta (aref ax (aref *row* k) (aref *col* k)))) 222 ;; UNDOES COLUMN HACK IN PIVOT. 223 (or *cpivot (do ((i 1 (1+ i))) ((> i m)) (setf (aref *col* i) i))) 224 (setq rank (min nrow nvar))) 225 226;; BACKWARD SUBSTITUTION 227(defun backward () 228 (declare(special ax delta m rank)) 229 (do ((i (1- rank) (1- i))) 230 ((< i 1)) 231 (do ((l (1+ rank) (1+ l))) 232 ((> l m)) 233 (setf (aref ax (aref *row* i) (aref *col* l)) 234 (let ((mess1 (pdifference 235 (ptimes (aref ax (aref *row* i) (aref *col* l)) 236 (aref ax (aref *row* rank) (aref *col* rank))) 237 (do ((j (1+ i) (1+ j)) (sum 0)) 238 ((> j rank) sum) 239 (setq sum (pplus sum (ptimes (aref ax (aref *row* i) (aref *col* j)) 240 (aref ax (aref *row* j) (aref *col* l))))))) ) 241 (mess2 (aref ax (aref *row* i) (aref *col* i)) )) 242 (cond ((equal 0 mess1) 0) 243 ((equal 0 mess2) 0) 244 (t ;; (pquotient mess1 mess2) ; fixed by line below. RJF 1/12/2017 245 246 (car (ratreduce mess1 mess2)) 247 ) 248 )))) 249 (do ((l (1+ i) (1+ l))) 250 ((> l rank)) 251 (setf (aref ax (aref *row* i) (aref *col* l)) 0))) 252 ;; PUT DELTA INTO THE DIAGONAL MATRIX 253 (setq delta (aref ax (aref *row* rank) (aref *col* rank))) 254 (do ((i 1 (1+ i))) 255 ((> i rank)) 256 (setf (aref ax (aref *row* i) (aref *col* i)) delta)) ) 257 258;;RECOVER THE ORDER OF ROWS AND COLUMNS. 259 260(defun recoverorder1 () 261 ;;(PRINT 'REARRANGE) 262 (do ((i nvar (1- i))) 263 ((= i 0)) 264 (setq variableorder (cons i variableorder))) 265 (do ((i (1+ rank) (1+ i))) 266 ((> i n)) 267 (cond ((equal (aref ax (aref *row* i) (aref *col* m)) 0) 268 (setq dependentrows (cons (aref *row* i) dependentrows))) 269 (t (setq inconsistentrows (cons (aref *row* i) inconsistentrows))))) 270 (do ((i 1 (1+ i))) 271 ((> i n)) 272 (cond ((not (= (aref *row* (aref *colinv* i)) i)) 273 (prog () 274 (moverow ax n m i 0) 275 (setq l i) 276 loop 277 (setq k (aref *row* (aref *colinv* l))) 278 (setf (aref *row* (aref *colinv* l)) l) 279 (cond ((= k i) (moverow ax n m 0 l)) 280 (t (moverow ax n m k l) 281 (setq l k) 282 (go loop)))))))) 283 284(defun recoverorder2 () 285 (do ((i nvar (1- i))) 286 ((= i 0)) 287 (setq variableorder (cons (aref *col* i) variableorder))) 288 (do ((i (1+ rank) (1+ i))) 289 ((> i n)) 290 (cond ((equal (aref ax (aref *row* i) (aref *col* m)) 0) 291 (setq dependentrows (cons (aref *row* i) dependentrows))) 292 (t (setq inconsistentrows (cons (aref *row* i) inconsistentrows))))) 293 (do ((i 1 (1+ i))) 294 ((> i n)) 295 (cond ((not (= (aref *row* i) i)) 296 (prog () 297 (moverow ax n m i 0) 298 (setq l i) 299 loop 300 (setq k (aref *row* l)) 301 (setf (aref *row* l) l) 302 (cond ((= k i) (moverow ax n m 0 l)) 303 (t (moverow ax n m k l) 304 (setq l k) 305 (go loop))))))) 306 (do ((i 1 (1+ i))) 307 ((> i nvar)) 308 (cond ((not (= (aref *col* i) i)) 309 (prog () 310 (movecol ax n m i 0) 311 (setq l i) 312 loop2 313 (setq k (aref *col* l)) 314 (setf (aref *col* l) l) 315 (cond ((= k i) (movecol ax n m 0 l)) 316 (t (movecol ax n m k l) 317 (setq l k) 318 (go loop2)))))))) 319 320;;THIS PROGRAM IS USED IN REARRANGEMENT 321(defun moverow (ax n m i j) 322 (do ((k 1 (1+ k))) ((> k m)) 323 (setf (aref ax j k) (aref ax i k)))) 324 325(defun movecol (ax n m i j) 326 (do ((k 1 (1+ k))) ((> k n)) 327 (setf (aref ax k j) (aref ax k i)))) 328 329;;COMPLEXITY IS DEFINED AS FOLLOWS 330;; COMPLEXITY(0)=0 331;; COMPLEXITY(CONSTANT)=1 332;; COMPLEXITY(POLYNOMIAL)=1+SUM(COMPLEXITY(C(N))+COMPLEXITY(E(N)), FOR N=0,1 ...M) 333;; WHERE POLYNOMIAL IS OF THE FORM 334;; SUM(C(N)*X^E(N), FOR N=0,1 ... M) X IS THE VARIABLE 335 336(defun complexity (exp) 337 (cond ((null exp) 0) 338 ((equal exp 0) 0) 339 ((atom exp) 1) 340 (t (+ (complexity (car exp)) (complexity (cdr exp)))))) 341 342(defun complexity/row (ax i j1 j2) 343 (do ((j j1 (1+ j)) (sum 0)) 344 ((> j j2) sum) 345 (incf sum (complexity (aref ax (aref *row* i) (aref *col* j)))))) 346 347(defun complexity/col (ax j i1 i2) 348 (do ((i i1 (1+ i)) (sum 0)) 349 ((> i i2) sum) 350 (incf sum (complexity (aref ax (aref *row* i) (aref *col* j)))))) 351 352(defun zerop/row (ax i j1 j2) 353 (do ((j j1 (1+ j))) 354 ((> j j2) t) 355 (cond ((not (equal (aref ax (aref *row* i) (aref *col* j)) 0)) (return nil))))) 356 357;;PIVOTTING ALGORITHM 358(defun pivot (ax k *cpivot) 359 (prog (row/optimal col/optimal complexity/i/min complexity/j/min 360 complexity/i complexity/j complexity/det complexity/det/min dummy) 361 (setq row/optimal k complexity/i/min 1000000. complexity/j/min 1000000.) 362 ;;TEST THE SINGULARITY 363 (cond ((do ((i k (1+ i)) (isallzero t)) 364 ((> i nrow) isallzero) 365 loop (cond ((zerop/row ax i k nvar) 366 (cond (*inv* (merror (intl:gettext "solve: singular matrix."))) 367 (t (exchangerow i nrow) 368 (decf nrow))) 369 (unless (> i nrow) (go loop))) 370 (t (setq isallzero nil)))) 371 (return t))) 372 373 ;; FIND AN OPTIMAL ROW 374 ;; IF *CPIVOT IS NIL, (AX I K) WHICH IS TO BE THE PIVOT MUST BE NONZERO. 375 ;; BUT IF *CPIVOT IS T, IT IS UNNECESSARY BECAUSE WE CAN DO THE COLUMN PIVOT. 376 findrow 377 (do ((i k (1+ i))) 378 ((> i nrow)) 379 (cond ((or *cpivot (not (equal (aref ax (aref *row* i) (aref *col* k)) 0))) 380 (cond ((> complexity/i/min 381 (setq complexity/i (complexity/row ax i k m))) 382 (setq row/optimal i complexity/i/min complexity/i)))))) 383 ;; EXCHANGE THE ROWS K AND ROW/OPTIMAL 384 (exchangerow k row/optimal) 385 386 ;; IF THE FLAG *CPIVOT IS NIL, THE FOLLOWING STEPS ARE NOT EXECUTED. 387 ;; THIS TREATMENT WAS DONE FOR THE LSA AND ECHELONTHINGS WHICH ARE NOT 388 ;; HAPPY WITH THE COLUMN OPERATIONS. 389 (cond ((null *cpivot) 390 (cond ((not (equal (aref ax (aref *row* k) (aref *col* k)) 0)) 391 (return nil)) 392 (t (do ((i k (1+ i))) ((= i nvar)) 393 (setf (aref *col* i) (aref *col* (1+ i)))) 394 (setq nvar (1- nvar) m (1- m)) 395 (go findrow))))) 396 397 ;;STEP3 ... FIND THE OPTIMAL COLUMN 398 (setq col/optimal 0 399 complexity/det/min 1000000. 400 complexity/j/min 1000000.) 401 402 (do ((j k (1+ j))) 403 ((> j nvar)) 404 (cond ((not (equal (aref ax (aref *row* k) (aref *col* j)) 0)) 405 (cond ((> complexity/det/min 406 (setq complexity/det 407 (complexity (aref ax (aref *row* k) (aref *col* j))))) 408 (setq col/optimal j 409 complexity/det/min complexity/det 410 complexity/j/min (complexity/col ax j (1+ k) n))) 411 ((equal complexity/det/min complexity/det) 412 (cond ((> complexity/j/min 413 (setq complexity/j 414 (complexity/col ax j (1+ k) n))) 415 (setq col/optimal j 416 complexity/det/min complexity/det 417 complexity/j/min complexity/j)))))))) 418 419 ;; EXCHANGE THE COLUMNS K AND COL/OPTIMAL 420 (exchangecol k col/optimal) 421 (setq dummy (aref *colinv* (aref *col* k))) 422 (setf (aref *colinv* (aref *col* k)) (aref *colinv* (aref *col* col/optimal))) 423 (setf (aref *colinv* (aref *col* col/optimal)) dummy) 424 (return nil))) 425 426(defun exchangerow (i j) 427 (prog (dummy) 428 (setq dummy (aref *row* i)) 429 (setf (aref *row* i) (aref *row* j)) 430 (setf (aref *row* j) dummy) 431 (cond ((= i j) (return nil)) 432 (t (setq permsign (not permsign)))))) 433 434(defun exchangecol (i j) 435 (prog (dummy) 436 (setq dummy (aref *col* i)) 437 (setf (aref *col* i) (aref *col* j)) 438 (setf (aref *col* j) dummy) 439 (cond ((= i j) (return nil)) 440 (t (setq permsign (not permsign)))))) 441 442;; Displays list of solutions. 443 444(defun solve2 (llist) 445 (setq $multiplicities nil) 446 (map2c #'(lambda (equatn multipl) 447 (setq equations 448 (nconc equations (list (displine equatn)))) 449 (push multipl $multiplicities) 450 (if (and (> multipl 1) $dispflag) 451 (mtell (intl:gettext "solve: multiplicity ~A~%") multipl))) 452 llist) 453 (setq $multiplicities (cons '(mlist simp) (nreverse $multiplicities)))) 454 455;; Displays an expression and returns its linelabel. 456 457(defun displine (exp) 458 (let ($nolabels (tim 0)) 459 (elabel exp) 460 (cond ($dispflag (remprop *linelabel* 'nodisp) 461 (setq tim (get-internal-run-time)) 462 (mterpri) 463 (displa (list '(mlabel) *linelabel* exp)) 464 (timeorg tim)) 465 (t (putprop *linelabel* t 'nodisp))) 466 *linelabel*)) 467 468(declare-top (unspecial permsign a rank delta nrow nvar n m variableorder 469 dependentrows inconsistentrows l k)) 470