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(in-package :maxima) 11 12;; Non-commutative product and exponentiation simplifier 13;; Written: July 1978 by CWH 14 15;; Flags to control simplification: 16 17(macsyma-module mdot) 18 19(defmvar $dotconstrules t 20 "Causes a non-commutative product of a constant and 21another term to be simplified to a commutative product. Turning on this 22flag effectively turns on DOT0SIMP, DOT0NSCSIMP, and DOT1SIMP as well.") 23 24(defmvar $dot0simp t 25 "Causes a non-commutative product of zero and a scalar term to 26be simplified to a commutative product.") 27 28(defmvar $dot0nscsimp t 29 "Causes a non-commutative product of zero and a nonscalar term 30to be simplified to a commutative product.") 31 32(defmvar $dot1simp t 33 "Causes a non-commutative product of one and another term to be 34simplified to a commutative product.") 35 36(defmvar $dotscrules nil 37 "Causes a non-commutative product of a scalar and another term to 38be simplified to a commutative product. Scalars and constants are carried 39to the front of the expression.") 40 41(defmvar $dotdistrib nil 42 "Causes every non-commutative product to be expanded each time it 43is simplified, i.e. A . (B + C) will simplify to A . B + A . C.") 44 45(defmvar $dotexptsimp t "Causes A . A to be simplified to A ^^ 2.") 46 47(defmvar $dotassoc t 48 "Causes a non-commutative product to be considered associative, so 49that A . (B . C) is simplified to A . B . C. If this flag is off, dot is 50taken to be right associative, i.e. A . B . C is simplified to A . (B . C).") 51 52(defmvar $doallmxops t 53 "Causes all operations relating to matrices (and lists) to be 54carried out. For example, the product of two matrices will actually be 55computed rather than simply being returned. Turning on this switch 56effectively turns on the following three.") 57 58(defmvar $domxmxops t "Causes matrix-matrix operations to be carried out.") 59 60(defmvar $doscmxops nil "Causes scalar-matrix operations to be carried out.") 61 62(defmvar $domxnctimes nil 63 "Causes non-commutative products of matrices to be carried out.") 64 65(defmvar $scalarmatrixp t 66 "Causes a square matrix of dimension one to be converted to a 67scalar, i.e. its only element.") 68 69(defmvar $dotident 1 "The value to be returned by X^^0.") 70 71(defmvar $assumescalar t 72 "This governs whether unknown expressions 'exp' are assumed to behave 73like scalars for combinations of the form 'exp op matrix' where op is one of 74{+, *, ^, .}. It has three settings: 75 76FALSE -- such expressions behave like non-scalars. 77TRUE -- such expressions behave like scalars only for the commutative 78 operators but not for non-commutative multiplication. 79ALL -- such expressions will behave like scalars for all operators 80 listed above. 81 82Note: This switch is primarily for the benefit of old code. If possible, 83you should declare your variables to be SCALAR or NONSCALAR so that there 84is no need to rely on the setting of this switch.") 85 86;; Specials defined elsewhere. 87 88(declare-top (special $expop $expon ; Controls behavior of EXPAND 89 errorsw)) 90 91;; The operators "." and "^^" distribute over equations. 92 93(defprop mnctimes (mequal) distribute_over) 94(defprop mncexpt (mequal) distribute_over) 95 96(defun simpnct (exp vestigial simp-flag) 97 (declare (ignore vestigial)) 98 (let ((check exp) 99 (first-factor (simpcheck (cadr exp) simp-flag)) 100 (remainder (if (cdddr exp) 101 (ncmuln (cddr exp) simp-flag) 102 (simpcheck (caddr exp) simp-flag)))) 103 (cond ((null (cdr exp)) $dotident) 104 ((null (cddr exp)) first-factor) 105 106 ;; This does (. sc m) --> (f* sc m) and (. (f* sc m1) m2) --> (f* sc (. m1 m2)) 107 ;; and (. m1 (f* sc m2)) --> (f* sc (. m1 m2)) where sc can be a scalar 108 ;; or constant, and m1 and m2 are non-constant, non-scalar expressions. 109 110 ((commutative-productp first-factor remainder) 111 (mul2 first-factor remainder)) 112 ((product-with-inner-scalarp first-factor) 113 (let ((p-p (partition-product first-factor))) 114 (outer-constant (car p-p) (cdr p-p) remainder))) 115 ((product-with-inner-scalarp remainder) 116 (let ((p-p (partition-product remainder))) 117 (outer-constant (car p-p) first-factor (cdr p-p)))) 118 119 ;; This code does distribution when flags are set and when called by 120 ;; $EXPAND. The way we recognize if we are called by $EXPAND is to look at 121 ;; the value of $EXPOP, but this is a kludge since $EXPOP has nothing to do 122 ;; with expanding (. A (f+ B C)) --> (f+ (. A B) (. A C)). I think that 123 ;; $EXPAND wants to have two flags: one which says to convert 124 ;; exponentiations to repeated products, and another which says to 125 ;; distribute products over sums. 126 127 ((and (mplusp first-factor) (or $dotdistrib (not (zerop $expop)))) 128 (addn (mapcar #'(lambda (x) (ncmul x remainder)) 129 (cdr first-factor)) 130 t)) 131 ((and (mplusp remainder) (or $dotdistrib (not (zerop $expop)))) 132 (addn (mapcar #'(lambda (x) (ncmul first-factor x)) 133 (cdr remainder)) 134 t)) 135 136 ;; This code carries out matrix operations when flags are set. 137 138 ((matrix-matrix-productp first-factor remainder) 139 (timex first-factor remainder)) 140 ((or (scalar-matrix-productp first-factor remainder) 141 (scalar-matrix-productp remainder first-factor)) 142 (simplifya (outermap1 'mnctimes first-factor remainder) t)) 143 144 ;; (. (^^ x n) (^^ x m)) --> (^^ x (f+ n m)) 145 146 ((and (simpnct-alike first-factor remainder) $dotexptsimp) 147 (simpnct-merge-factors first-factor remainder)) 148 149 ;; (. (. x y) z) --> (. x y z) 150 151 ((and (mnctimesp first-factor) $dotassoc) 152 (ncmuln (append (cdr first-factor) 153 (if (mnctimesp remainder) 154 (cdr remainder) 155 (ncons remainder))) 156 t)) 157 158 ;; (. (^^ (. x y) m) (^^ (. x y) n) z) --> (. (^^ (. x y) m+n) z) 159 ;; (. (^^ (. x y) m) x y z) --> (. (^^ (. x y) m+1) z) 160 ;; (. x y (^^ (. x y) m) z) --> (. (^^ (. x y) m+1) z) 161 ;; (. x y x y z) --> (. (^^ (. x y) 2) z) 162 163 ((and (mnctimesp remainder) $dotassoc $dotexptsimp) 164 (setq exp (simpnct-merge-product first-factor (cdr remainder))) 165 (if (and (mnctimesp exp) $dotassoc) 166 (simpnct-antisym-check (cdr exp) check) 167 (eqtest exp check))) 168 169 ;; (. x (. y z)) --> (. x y z) 170 171 ((and (mnctimesp remainder) $dotassoc) 172 (simpnct-antisym-check (cons first-factor (cdr remainder)) check)) 173 174 (t (eqtest (list '(mnctimes) first-factor remainder) check))))) 175 176;; Predicate functions for simplifying a non-commutative product to a 177;; commutative one. SIMPNCT-CONSTANTP actually determines if a term is a 178;; constant and is not a nonscalar, i.e. not declared nonscalar and not a 179;; constant list or matrix. The function CONSTANTP determines if its argument 180;; is a number or a variable declared constant. 181 182(defun commutative-productp (first-factor remainder) 183 (or (simpnct-sc-or-const-p first-factor) 184 (simpnct-sc-or-const-p remainder) 185 (simpnct-onep first-factor) 186 (simpnct-onep remainder) 187 (zero-productp first-factor remainder) 188 (zero-productp remainder first-factor))) 189 190(defun simpnct-sc-or-const-p (term) 191 (or (simpnct-constantp term) (simpnct-assumescalarp term))) 192 193(defun simpnct-constantp (term) 194 (and $dotconstrules 195 (or (mnump term) 196 (and ($constantp term) (not ($nonscalarp term)))))) 197 198(defun simpnct-assumescalarp (term) 199 (and $dotscrules (scalar-or-constant-p term (eq $assumescalar '$all)))) 200 201(defun simpnct-onep (term) (and $dot1simp (onep1 term))) 202 203(defun zero-productp (one-term other-term) 204 (and (zerop1 one-term) 205 $dot0simp 206 (or $dot0nscsimp (not ($nonscalarp other-term))))) 207 208;; This function takes a form and determines if it is a product 209;; containing a constant or a declared scalar. Note that in the 210;; next three functions, the word "scalar" is used to refer to a constant 211;; or a declared scalar. This is a bad way of doing things since we have 212;; to cdr down an expression twice: once to determine if a scalar is there 213;; and once again to pull it out. 214 215(defun product-with-inner-scalarp (product) 216 (and (mtimesp product) 217 (or $dotconstrules $dotscrules) 218 (do ((factor-list (cdr product) (cdr factor-list))) 219 ((null factor-list) nil) 220 (if (simpnct-sc-or-const-p (car factor-list)) 221 (return t))))) 222 223;; This function takes a commutative product and separates it into a scalar 224;; part and a non-scalar part. 225 226(defun partition-product (product) 227 (do ((factor-list (cdr product) (cdr factor-list)) 228 (scalar-list nil) 229 (nonscalar-list nil)) 230 ((null factor-list) (cons (nreverse scalar-list) 231 (muln (nreverse nonscalar-list) t))) 232 (if (simpnct-sc-or-const-p (car factor-list)) 233 (push (car factor-list) scalar-list) 234 (push (car factor-list) nonscalar-list)))) 235 236;; This function takes a list of constants and scalars, and two nonscalar 237;; expressions and forms a non-commutative product of the nonscalar 238;; expressions, and a commutative product of the constants and scalars and 239;; the non-commutative product. 240 241(defun outer-constant (constant nonscalar1 nonscalar2) 242 (muln (nconc constant (ncons (ncmul nonscalar1 nonscalar2))) t)) 243 244(defun simpnct-base (term) (if (mncexptp term) (cadr term) term)) 245 246(defun simpnct-power (term) (if (mncexptp term) (caddr term) 1)) 247 248(defun simpnct-alike (term1 term2) 249 (alike1 (simpnct-base term1) (simpnct-base term2))) 250 251(defun simpnct-merge-factors (term1 term2) 252 (ncpower (simpnct-base term1) 253 (add2 (simpnct-power term1) (simpnct-power term2)))) 254 255(defun matrix-matrix-productp (term1 term2) 256 (and (or $doallmxops $domxmxops $domxnctimes) 257 (mxorlistp1 term1) 258 (mxorlistp1 term2))) 259 260(defun scalar-matrix-productp (term1 term2) 261 (and (or $doallmxops $doscmxops) 262 (mxorlistp1 term1) 263 (scalar-or-constant-p term2 (eq $assumescalar '$all)))) 264 265 266(defun simpncexpt (exp vestigial simp-flag) 267 (declare (ignore vestigial)) 268 (let ((factor (simpcheck (cadr exp) simp-flag)) 269 (power (simpcheck (caddr exp) simp-flag)) 270 (check exp)) 271 (twoargcheck exp) 272 (cond ((zerop1 power) 273 (if (zerop1 factor) 274 (if (not errorsw) 275 (merror (intl:gettext "noncommutative exponent: ~M is undefined.") 276 (list '(mncexpt) factor power)) 277 (throw 'errorsw t))) 278 (if (mxorlistp1 factor) (identitymx factor) $dotident)) 279 ((onep1 power) factor) 280 ((and (simpnct-sc-or-const-p factor) 281 (simpnct-sc-or-const-p power)) (power factor power)) 282 ((and (zerop1 factor) $dot0simp) factor) 283 ((and (onep1 factor) $dot1simp) factor) 284 ((and (or $doallmxops $domxmxops) 285 (mxorlistp1 factor) 286 (fixnump power)) 287 (let (($scalarmatrixp (or ($listp factor) $scalarmatrixp))) 288 (simplify (powerx factor power)))) 289 290 ;; This does (A+B)^^2 --> A^^2 + A.B + B.A + B^^2 291 ;; and (A.B)^^2 --> A.B.A.B 292 293 ((and (or (mplusp factor) (not $dotexptsimp)) 294 (fixnump power) 295 (not (> power $expop)) 296 (plusp power)) 297 (ncmul factor (ncpower factor (1- power)))) 298 299 ;; This does the same thing as above for (A+B)^^(-2) 300 ;; and (A.B)^^(-2). Here the "-" operator does the trick 301 ;; for us. 302 303 ((and (or (mplusp factor) (not $dotexptsimp)) 304 (fixnump power) 305 (not (> (- power) $expon)) 306 (< power -1)) 307 (let (($expop $expon)) 308 (ncpower (ncpower factor (- power)) -1))) 309 310 ((product-with-inner-scalarp factor) 311 (let ((p-p (partition-product factor))) 312 (mul2 (power (muln (car p-p) t) power) 313 (ncpower (cdr p-p) power)))) 314 ((and $dotassoc (mncexptp factor)) 315 (ncpower (cadr factor) (mul2 (caddr factor) power))) 316 (t (eqtest (list '(mncexpt) factor power) check))))) 317 318 319(defun simpnct-invert (exp) 320 (cond ((mnctimesp exp) 321 (ncmuln (nreverse (mapcar #'simpnct-invert (cdr exp))) t)) 322 ((and (mncexptp exp) (integerp (caddr exp))) 323 (ncpower (cadr exp) (- (caddr exp)))) 324 (t (list '(mncexpt simp) exp -1)))) 325 326(defun identitymx (x) 327 (if (and ($listp (cadr x)) (= (length (cdr x)) (length (cdadr x)))) 328 (simplifya (cons (car x) (cdr ($ident (length (cdr x))))) t) 329 $dotident)) 330 331;; This function incorporates the hairy search which enables such 332;; simplifications as (. a b a b) --> (^^ (. a b) 2). It assumes 333;; that FIRST-FACTOR is not a dot product and that REMAINDER is. 334;; For the product (. a b c d e), three basic types of comparisons 335;; are done: 336;; 337;; 1) a <---> b first-factor <---> inner-product 338;; a <---> (. b c) 339;; a <---> (. b c d) 340;; a <---> (. b c d e) (this case handled in SIMPNCT) 341;; 342;; 2) (. a b) <---> c outer-product <---> (car rest) 343;; (. a b c) <---> d 344;; (. a b c d) <---> e 345;; 346;; 3) (. a b) <---> (. c d) outer-product <---> (firstn rest) 347;; 348;; Note that INNER-PRODUCT and OUTER-PRODUCT share list structure which 349;; is clobbered as new terms are added. 350 351(defun simpnct-merge-product (first-factor remainder) 352 (let ((half-product-length (ash (1+ (length remainder)) -1)) 353 (inner-product (car remainder)) 354 (outer-product (list '(mnctimes) first-factor (car remainder)))) 355 (do ((merge-length 2 (1+ merge-length)) 356 (rest (cdr remainder) (cdr rest))) 357 ((null rest) outer-product) 358 (cond ((simpnct-alike first-factor inner-product) 359 (return 360 (ncmuln 361 (cons (simpnct-merge-factors first-factor inner-product) 362 rest) 363 t))) 364 ((simpnct-alike outer-product (car rest)) 365 (return 366 (ncmuln 367 (cons (simpnct-merge-factors outer-product (car rest)) 368 (cdr rest)) 369 t))) 370 ((and (not (> merge-length half-product-length)) 371 (alike1 outer-product 372 (cons '(mnctimes) 373 (subseq rest 0 merge-length)))) 374 (return 375 (ncmuln (cons (ncpower outer-product 2) 376 (nthcdr merge-length rest)) 377 t))) 378 ((= merge-length 2) 379 (setq inner-product 380 (cons '(mnctimes) (cddr outer-product))))) 381 (rplacd (last inner-product) (ncons (car rest)))))) 382 383(defun simpnct-antisym-check (l check) 384 (cond ((and (get 'mnctimes '$antisymmetric) (cddr l)) 385 (multiple-value-bind (l antisym-sign) (bbsort1 l) 386 (cond ((equal l 0) 0) 387 ((prog1 (null antisym-sign) 388 (setq l (eqtest (cons '(mnctimes) l) check))) 389 l) 390 (t (neg l))))) 391 (t (eqtest (cons '(mnctimes) l) check)))) 392