src/mdot.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancements.                   ;;;;;
   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 $domxnctimes nil
  37   "Causes non-commutative products of matrices to be carried out.")
  38
  39 (defmvar $dotident 1 "The value to be returned by X^^0.")
  40
  41 ;; The operators "." and "^^" distribute over equations.
  42
  43 (defprop mnctimes (mequal) distribute_over)
  44 (defprop mncexpt (mequal) distribute_over)
  45
  46 (defun simpnct (exp vestigial simp-flag)
  47   (declare (ignore vestigial))
  48   (let ((check exp)
  49         (first-factor (simpcheck (cadr exp) simp-flag))
  50         (remainder (if (cdddr exp)
  51                        (ncmuln (cddr exp) simp-flag)
  52                        (simpcheck (caddr exp) simp-flag))))
  53     (cond ((null (cdr exp)) $dotident)
  54           ((null (cddr exp)) first-factor)
  55
  56           ;;  This does (. sc m) --> (f* sc m)  and  (. (f* sc m1) m2) --> (f* sc (. m1 m2))
  57           ;;  and (. m1 (f* sc m2)) --> (f* sc (. m1 m2)) where sc can be a scalar
  58           ;;  or constant, and m1 and m2 are non-constant, non-scalar expressions.
  59
  60           ((commutative-productp first-factor remainder)
  61            (mul2 first-factor remainder))
  62           ((product-with-inner-scalarp first-factor)
  63            (let ((p-p (partition-product first-factor)))
  64              (outer-constant (car p-p) (cdr p-p) remainder)))
  65           ((product-with-inner-scalarp remainder)
  66            (let ((p-p (partition-product remainder)))
  67              (outer-constant (car p-p) first-factor (cdr p-p))))
  68
  69           ;;  This code does distribution when flags are set and when called by
  70           ;;  $EXPAND.  The way we recognize if we are called by $EXPAND is to look at
  71           ;;  the value of $EXPOP, but this is a kludge since $EXPOP has nothing to do
  72           ;;  with expanding (. A (f+ B C)) --> (f+ (. A B) (. A C)).  I think that
  73           ;;  $EXPAND wants to have two flags:  one which says to convert
  74           ;;  exponentiations to repeated products, and another which says to
  75           ;;  distribute products over sums.
  76
  77           ((and (mplusp first-factor) (or $dotdistrib (not (zerop $expop))))
  78            (addn (mapcar #'(lambda (x) (ncmul x remainder))
  79                          (cdr first-factor))
  80                  t))
  81           ((and (mplusp remainder) (or $dotdistrib (not (zerop $expop))))
  82            (addn (mapcar #'(lambda (x) (ncmul first-factor x))
  83                          (cdr remainder))
  84                  t))
  85
  86           ;;  This code carries out matrix operations when flags are set.
  87
  88           ((matrix-matrix-productp first-factor remainder)
  89            (timex first-factor remainder))
  90           ((or (scalar-matrix-productp first-factor remainder)
  91                (scalar-matrix-productp remainder first-factor))
  92            (simplifya (outermap1 'mnctimes first-factor remainder) t))
  93
  94           ;;  (. (^^ x n) (^^ x m)) --> (^^ x (f+ n m))
  95
  96           ((and (simpnct-alike first-factor remainder) $dotexptsimp)
  97            (simpnct-merge-factors first-factor remainder))
  98
  99           ;;  (. (. x y) z) --> (. x y z)
 100
 101           ((and (mnctimesp first-factor) $dotassoc)
 102            (ncmuln (append (cdr first-factor)
 103                            (if (mnctimesp remainder)
 104                                (cdr remainder)
 105                                (ncons remainder)))
 106                    t))
 107
 108           ;;  (. (^^ (. x y) m) (^^ (. x y) n) z) --> (. (^^ (. x y) m+n) z)
 109           ;;  (. (^^ (. x y) m) x y z) --> (. (^^ (. x y) m+1) z)
 110           ;;  (. x y (^^ (. x y) m) z) --> (. (^^ (. x y) m+1) z)
 111           ;;  (. x y x y z) --> (. (^^ (. x y) 2) z)
 112
 113           ((and (mnctimesp remainder) $dotassoc $dotexptsimp)
 114            (setq exp (simpnct-merge-product first-factor (cdr remainder)))
 115            (if (and (mnctimesp exp) $dotassoc)
 116                (simpnct-antisym-check (cdr exp) check)
 117                (eqtest exp check)))
 118
 119           ;;  (. x (. y z)) --> (. x y z)
 120
 121           ((and (mnctimesp remainder) $dotassoc)
 122            (simpnct-antisym-check (cons first-factor (cdr remainder)) check))
 123
 124           (t (eqtest (list '(mnctimes) first-factor remainder) check)))))
 125
 126 ;;  Predicate functions for simplifying a non-commutative product to a
 127 ;;  commutative one.  SIMPNCT-CONSTANTP actually determines if a term is a
 128 ;;  constant and is not a nonscalar, i.e. not declared nonscalar and not a
 129 ;;  constant list or matrix.  The function CONSTANTP determines if its argument
 130 ;;  is a number or a variable declared constant.
 131
 132 (defun commutative-productp (first-factor remainder)
 133   (or (simpnct-sc-or-const-p first-factor)
 134       (simpnct-sc-or-const-p remainder)
 135       (simpnct-onep first-factor)
 136       (simpnct-onep remainder)
 137       (zero-productp first-factor remainder)
 138       (zero-productp remainder first-factor)))
 139
 140 (defun simpnct-sc-or-const-p (term)
 141   (or (simpnct-constantp term) (simpnct-assumescalarp term)))
 142
 143 (defun simpnct-constantp (term)
 144   (and $dotconstrules
 145        (or (mnump term)
 146            (and ($constantp term) (not ($nonscalarp term))))))
 147
 148 (defun simpnct-assumescalarp (term)
 149   (and $dotscrules (scalar-or-constant-p term (eq $assumescalar '$all))))
 150
 151 (defun simpnct-onep (term) (and $dot1simp (onep1 term)))
 152
 153 (defun zero-productp (one-term other-term)
 154   (and (zerop1 one-term)
 155        $dot0simp
 156        (or $dot0nscsimp (not ($nonscalarp other-term)))))
 157
 158 ;;  This function takes a form and determines if it is a product
 159 ;;  containing a constant or a declared scalar.  Note that in the
 160 ;;  next three functions, the word "scalar" is used to refer to a constant
 161 ;;  or a declared scalar.  This is a bad way of doing things since we have
 162 ;;  to cdr down an expression twice: once to determine if a scalar is there
 163 ;;  and once again to pull it out.
 164
 165 (defun product-with-inner-scalarp (product)
 166   (and (mtimesp product)
 167        (or $dotconstrules $dotscrules)
 168        (do ((factor-list (cdr product) (cdr factor-list)))
 169            ((null factor-list) nil)
 170          (if (simpnct-sc-or-const-p (car factor-list))
 171              (return t)))))
 172
 173 ;;  This function takes a commutative product and separates it into a scalar
 174 ;;  part and a non-scalar part.
 175
 176 (defun partition-product (product)
 177   (do ((factor-list (cdr product) (cdr factor-list))
 178        (scalar-list nil)
 179        (nonscalar-list nil))
 180       ((null factor-list) (cons (nreverse scalar-list)
 181                                 (muln (nreverse nonscalar-list) t)))
 182     (if (simpnct-sc-or-const-p (car factor-list))
 183         (push (car factor-list) scalar-list)
 184         (push (car factor-list) nonscalar-list))))
 185
 186 ;;  This function takes a list of constants and scalars, and two nonscalar
 187 ;;  expressions and forms a non-commutative product of the nonscalar
 188 ;;  expressions, and a commutative product of the constants and scalars and
 189 ;;  the non-commutative product.
 190
 191 (defun outer-constant (constant nonscalar1 nonscalar2)
 192   (muln (nconc constant (ncons (ncmul nonscalar1 nonscalar2))) t))
 193
 194 (defun simpnct-base (term) (if (mncexptp term) (cadr term) term))
 195
 196 (defun simpnct-power (term) (if (mncexptp term) (caddr term) 1))
 197
 198 (defun simpnct-alike (term1 term2)
 199   (alike1 (simpnct-base term1) (simpnct-base term2)))
 200
 201 (defun simpnct-merge-factors (term1 term2)
 202   (ncpower (simpnct-base term1)
 203            (add2 (simpnct-power term1) (simpnct-power term2))))
 204
 205 (defun matrix-matrix-productp (term1 term2)
 206   (and (or $doallmxops $domxmxops $domxnctimes)
 207        (mxorlistp1 term1)
 208        (mxorlistp1 term2)))
 209
 210 (defun scalar-matrix-productp (term1 term2)
 211   (and (or $doallmxops $doscmxops)
 212        (mxorlistp1 term1)
 213        (scalar-or-constant-p term2 (eq $assumescalar '$all))))
 214
 215
 216 (defun simpncexpt (exp vestigial simp-flag)
 217   (declare (ignore vestigial))
 218   (let ((factor (simpcheck (cadr exp) simp-flag))
 219         (power (simpcheck (caddr exp) simp-flag))
 220         (check exp))
 221     (twoargcheck exp)
 222     (cond ((zerop1 power)
 223            (if (zerop1 factor)
 224                (if (not errorsw)
 225                    (merror (intl:gettext "noncommutative exponent: ~M is undefined.")
 226                            (list '(mncexpt) factor power))
 227                    (throw 'errorsw t)))
 228            (if (mxorlistp1 factor) (identitymx factor) $dotident))
 229           ((onep1 power) factor)
 230           ((and (simpnct-sc-or-const-p factor)
 231                 (simpnct-sc-or-const-p power)) (power factor power))
 232           ((and (zerop1 factor) $dot0simp) factor)
 233           ((and (onep1 factor) $dot1simp) factor)
 234           ((and (or $doallmxops $domxmxops)
 235                 (mxorlistp1 factor)
 236                 (fixnump power))
 237            (let (($scalarmatrixp (or ($listp factor) $scalarmatrixp)))
 238              (simplify (powerx factor power))))
 239
 240           ;; This does (A+B)^^2 --> A^^2 + A.B + B.A + B^^2
 241           ;; and (A.B)^^2 --> A.B.A.B
 242
 243           ((and (or (mplusp factor) (not $dotexptsimp))
 244                 (fixnump power)
 245                 (not (> power $expop))
 246                 (plusp power))
 247            (ncmul factor (ncpower factor (1- power))))
 248
 249           ;; This does the same thing as above for (A+B)^^(-2)
 250           ;; and (A.B)^^(-2).  Here the "-" operator does the trick
 251           ;; for us.
 252
 253           ((and (or (mplusp factor) (not $dotexptsimp))
 254                 (fixnump power)
 255                 (not (> (- power) $expon))
 256                 (< power -1))
 257            (let (($expop $expon))
 258              (ncpower (ncpower factor (- power)) -1)))
 259
 260           ((product-with-inner-scalarp factor)
 261            (let ((p-p (partition-product factor)))
 262              (mul2 (power (muln (car p-p) t) power)
 263                    (ncpower (cdr p-p) power))))
 264           ((and $dotassoc (mncexptp factor))
 265            (ncpower (cadr factor) (mul2 (caddr factor) power)))
 266           (t (eqtest (list '(mncexpt) factor power) check)))))
 267
 268
 269 (defun simpnct-invert (exp)
 270   (cond ((mnctimesp exp)
 271          (ncmuln (nreverse (mapcar #'simpnct-invert (cdr exp))) t))
 272         ((and (mncexptp exp) (integerp (caddr exp)))
 273          (ncpower (cadr exp) (- (caddr exp))))
 274         (t (list '(mncexpt simp) exp -1))))
 275
 276 (defun identitymx (x)
 277   (if (and ($listp (cadr x)) (= (length (cdr x)) (length (cdadr x))))
 278       (simplifya (cons (car x) (cdr ($ident (length (cdr x))))) t)
 279       $dotident))
 280
 281 ;;  This function incorporates the hairy search which enables such
 282 ;;  simplifications as (. a b a b) --> (^^ (. a b) 2).  It assumes
 283 ;;  that FIRST-FACTOR is not a dot product and that REMAINDER is.
 284 ;;  For the product (. a b c d e), three basic types of comparisons
 285 ;;  are done:
 286 ;;
 287 ;;  1)  a <---> b               first-factor <---> inner-product
 288 ;;      a <---> (. b c)
 289 ;;      a <---> (. b c d)
 290 ;;      a <---> (. b c d e)     (this case handled in SIMPNCT)
 291 ;;
 292 ;;  2) (. a b) <---> c          outer-product <---> (car rest)
 293 ;;     (. a b c) <---> d
 294 ;;     (. a b c d) <---> e
 295 ;;
 296 ;;  3) (. a b) <---> (. c d)    outer-product <---> (firstn rest)
 297 ;;
 298 ;;  Note that INNER-PRODUCT and OUTER-PRODUCT share list structure which
 299 ;;  is clobbered as new terms are added.
 300
 301 (defun simpnct-merge-product (first-factor remainder)
 302   (let ((half-product-length (ash (1+ (length remainder)) -1))
 303         (inner-product (car remainder))
 304         (outer-product (list '(mnctimes) first-factor (car remainder))))
 305     (do ((merge-length 2 (1+ merge-length))
 306          (rest (cdr remainder) (cdr rest)))
 307         ((null rest) outer-product)
 308       (cond ((simpnct-alike first-factor inner-product)
 309              (return
 310                (ncmuln
 311                 (cons (simpnct-merge-factors first-factor inner-product)
 312                       rest)
 313                 t)))
 314             ((simpnct-alike outer-product (car rest))
 315              (return
 316                (ncmuln
 317                 (cons (simpnct-merge-factors outer-product (car rest))
 318                       (cdr rest))
 319                 t)))
 320             ((and (not (> merge-length half-product-length))
 321                   (alike1 outer-product
 322                           (cons '(mnctimes)
 323                                 (subseq rest 0 merge-length))))
 324              (return
 325                (ncmuln (cons (ncpower outer-product 2)
 326                              (nthcdr merge-length rest))
 327                        t)))
 328             ((= merge-length 2)
 329              (setq inner-product
 330                    (cons '(mnctimes) (cddr outer-product)))))
 331       (rplacd (last inner-product) (ncons (car rest))))))
 332
 333 (defun simpnct-antisym-check (l check)
 334   (cond ((and (get 'mnctimes '$antisymmetric) (cddr l))
 335          (multiple-value-bind (l antisym-sign) (bbsort1 l)
 336            (cond ((equal l 0) 0)
 337                  ((prog1 (null antisym-sign)
 338                     (setq l (eqtest (cons '(mnctimes) l) check)))
 339                   l)
 340                  (t (neg l)))))
 341         (t (eqtest (cons '(mnctimes) l) check))))