src/mdot.lisp

   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
  21 another term to be simplified to a commutative product.  Turning on this
  22 flag 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
  26 be simplified to a commutative product.")
  27
  28 (defmvar $dot0nscsimp t
  29   "Causes a non-commutative product of zero and a nonscalar term
  30 to be simplified to a commutative product.")
  31
  32 (defmvar $dot1simp t
  33   "Causes a non-commutative product of one and another term to be
  34 simplified to  a commutative product.")
  35
  36 (defmvar $dotscrules nil
  37   "Causes a non-commutative product of a scalar and another term to
  38 be simplified to a commutative product.  Scalars and constants are carried
  39 to the front of the expression.")
  40
  41 (defmvar $dotdistrib nil
  42   "Causes every non-commutative product to be expanded each time it
  43 is 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
  49 that A . (B . C) is simplified to A . B . C.  If this flag is off, dot is
  50 taken 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
  54 carried out.  For example, the product of two matrices will actually be
  55 computed rather than simply being returned.  Turning on this switch
  56 effectively 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
  67 scalar, 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
  73 like scalars for combinations of the form 'exp op matrix' where op is one of
  74 {+, *, ^, .}.  It has three settings:
  75
  76 FALSE -- such expressions behave like non-scalars.
  77 TRUE  -- such expressions behave like scalars only for the commutative
  78          operators but not for non-commutative multiplication.
  79 ALL   -- such expressions will behave like scalars for all operators
  80          listed above.
  81
  82 Note:  This switch is primarily for the benefit of old code.  If possible,
  83 you should declare your variables to be SCALAR or NONSCALAR so that there
  84 is 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))))