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