Rename *ll* and *ul* to ll and ul in in-interval
[maxima.git] / share / tensor / symtry.lisp
blob68a86b8688231255b7bcd9d4c0da505d17380649
1 ;;; -*- Mode:LISP; Package:MACSYMA -*-
2 ;;
3 ;; This program is free software; you can redistribute it and/or
4 ;; modify it under the terms of the GNU General Public License as
5 ;; published by the Free Software Foundation; either version 2 of
6 ;; the License, or (at your option) any later version.
7 ;;
8 ;; This program is distributed in the hope that it will be
9 ;; useful, but WITHOUT ANY WARRANTY; without even the implied
10 ;; warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 ;; PURPOSE. See the GNU General Public License for more details.
13 ;; Comments: Symmetrization module for itensor.lisp
15 ;;; symtry 100 Feb 12, 1982
17 (in-package :maxima)
19 ; ** (c) Copyright 1979 Massachusetts Institute of Technology **
21 (declare-top (special symtypes $symmetries $allsym csign smlist $idummyx))
23 (setq symtypes '($sym $anti $cyc) $symmetries '((mlist simp)))
25 ;$SYMMETRIES is a list of indexed objects with declared symmetries
26 ;$ALLSYM if non-nil means that all indexed objects are assumed symmetric
28 (defun $decsym (name ncov ncontr covl contrl) ;DEClare SYMmetries
29 (prog (tensor)
30 (cond ((not (symbolp name))
31 (merror "First argument must be a possible tensor name"))
32 ((not (and (fixnump ncov)
33 (fixnump ncontr)
34 (signp ge ncov)
35 (signp ge ncontr)))
36 (merror
37 "2nd and 3rd arguments must be non-negative integers"))
38 ((not (and (eq (caar covl) 'mlist)
39 (eq (caar contrl) 'mlist)))
40 (merror "4th and 5th arguments must be lists"))
41 ((and (< ncov 2) (< ncontr 2))
42 (merror "This object can have no symmetry properties"))
43 ((or (and (< ncov 2) (not (null (cdr covl))))
44 (and (< ncontr 2) (not (null (cdr contrl)))))
45 (merror
46 "Non-null list associated with zero or single index specification")))
47 (setq tensor (implode (nconc (exploden name) (ncons 45)
48 (exploden ncov) (ncons 45)
49 (exploden ncontr))))
50 (do ((covl (cdr covl) (cdr covl)) (carl) (arglist) (prop))
51 ((null covl))
52 (cond ((not (member (setq prop (caar (setq carl (car covl))))
53 symtypes :test #'equal))
54 (merror "Invalid symmetry operator: ~M" carl))
55 ((and (null (cddr carl)) (eq (cadr carl) '$all))
56 (setq arglist (interval 1 ncov)))
57 (t (setq arglist (check-symargs (cdr carl) (1+ ncov)))))
58 (setq carl (zl-get tensor prop))
59 (putprop tensor (cons (cons arglist (car carl)) (cdr carl))
60 prop))
61 (do ((contl (cdr contrl) (cdr contl)) (carl) (arglist) (prop))
62 ((null contl))
63 (cond ((not (member (setq prop (caar (setq carl (car contl))))
64 symtypes :test #'equal))
65 (merror "Invalid symmetry operator: ~M" carl))
66 ((and (null (cddr carl)) (eq (cadr carl) '$all))
67 (setq arglist (interval 1 ncontr)))
68 ((setq arglist (check-symargs (cdr carl) (1+ ncontr)))))
69 (setq carl (zl-get tensor prop))
70 (putprop tensor (cons (car carl) (cons arglist (cdr carl)))
71 prop))
72 (add2lnc tensor $symmetries)
73 (return '$done)))
75 ;(defun interval (i j) ;INTERVAL returns the list of integers from I thru J.
76 ; (do ((n i (1+ n)) (ans)) ;Thus (INTERVAL 3 5) yields (3 4 5)
77 ; ((> n j) (nreverse ans))
78 ; (setq ans (cons n ans))))
80 (defun check-symargs (ll n) ;Returns an ascending list of the unique
81 ;elements of LL and checks that they are
82 (do ((l ll (cdr l)) (c) (ans)) ;integers between 0 and N
83 ((null l) (cond ((null (cdr ans))
84 (merror "Only one distinct index in symmetry property declaration"))
85 (t (sort ans '<))))
86 (setq c (car l))
87 (cond ((not (and (fixnump c) (< 0 c n)))
88 (merror "Bad argument encountered for symmetry operator"))
89 ((not (member c ans :test #'equal)) (setq ans (cons c ans))))))
91 (defun $dispsym (name ncov ncontr) ;DISPlay SYMmetries
92 (prog (tensor)
93 (setq tensor (implode (nconc (exploden name) (ncons 45)
94 (exploden ncov) (ncons 45)
95 (exploden ncontr))))
96 (cond ((not (member tensor (cdr $symmetries) :test #'equal))
97 (return (ncons smlist))))
98 (return
99 (do ((q symtypes (cdr q)) (l) (prop))
100 ((null q) (consmlist l))
101 (cond ((not (null (setq prop (zl-get tensor (car q)))))
102 (setq l
103 (append l
104 (ncons
105 (consmlist
106 (list
107 (car q)
108 (consmlist (mapcar 'consmlist (car prop)))
109 (consmlist (mapcar 'consmlist (cdr prop))))
110 ))))))))))
112 (defun $remsym (name ncov ncontr)
113 (prog (tensor)
114 (setq tensor (implode (nconc (exploden name) (ncons 45)
115 (exploden ncov) (ncons 45)
116 (exploden ncontr))))
117 (cond ((not (member tensor (cdr $symmetries) :test #'equal))
118 (mtell "~&No symmetries have been declared for this tensor.~%"))
119 (t (setq $symmetries (delete tensor $symmetries :test #'equal))
120 (zl-remprop tensor '$sym)
121 (zl-remprop tensor '$anti)
122 (zl-remprop tensor '$cyc)))
123 (return '$done)))
125 (defun $canform (&rest args) ;Convert E into CANonical FORM
126 (prog (e f)
127 (cond
128 ( (equal (length args) 1) (setq f t))
129 ( (equal (length args) 2) (setq f (cadr args)))
130 (t (merror "CANFORM requires one or two arguments"))
132 (setq e (car args))
133 (return
134 (cond ((atom e) e)
135 ((eq (caar e) 'mequal)
136 (mysubst0 (list (car e) ($canform (cadr e) f) ($canform (caddr e) f))
138 ((eq (caar e) 'mplus)
139 (mysubst0 (simplus (cons '(mplus) (mapcar (lambda (ee) ($canform ee f)) (cdr e)))
140 1 nil) e))
141 ((eq (caar e) 'mtimes) (mysubst0 (simplifya (canprod e f) nil) e))
142 ((rpobj e) (canten e f))
143 (t (mysubst0 (simplifya (cons (ncons (caar e))
144 (mapcar (lambda (ee) ($canform ee f)) (cdr e))) t) e))))))
146 (defun canten (e nfprpobjs) ;CANonical TENsor
147 (prog (cov contr deriv tensor)
148 ((lambda (dummy) (and nfprpobjs dummy (setq e (rename1 e dummy))))
149 (nonumber (cdaddr ($indices e)))) ;NFPRPOBJS is Not From Product
150 (setq cov (copy-tree (cdadr e)) ;of RP (indexed) OBJects
151 contr (copy-tree (cdaddr e))
152 deriv (copy-tree (cdddr e))
153 tensor (implode (nconc (exploden (caar e)) (ncons 45)
154 (exploden (length cov)) (ncons 45)
155 (exploden (length contr))))
156 csign nil) ;Set when reordering antisymmetric indices.
157 ;Indicates whether overall sign of
158 ;expression needs changing.
159 (cond ((or (or (eq (caar e) '$levi_civita) (eq (caar e) '%levi_civita))
160 (or (eq (caar e) '$kdelta) (eq (caar e) '%kdelta)))
161 (setq cov (antisort cov) contr (antisort contr)))
162 ((or $allsym (eq (caar e) '$kdels) (eq (caar e) '%kdels))
163 (setq cov (itensor-sort cov) contr (itensor-sort contr)))
164 ((member ($verbify tensor) (cdr $symmetries) :test #'equal)
165 (do ((q symtypes (cdr q)) (type))
166 ((null q))
167 (setq type (car q))
168 (do ((props (car (zl-get ($verbify tensor) type)) (cdr props)) (p))
169 ((null props))
170 (setq p (car props)
171 cov (inserts (symsort (extract-elements p cov) type)
172 cov p)))
173 (do ((props (cdr (zl-get ($verbify tensor) type)) (cdr props)) (p))
174 ((null props))
175 (setq p (car props)
176 contr (inserts (symsort (extract-elements p contr)
177 type)
178 contr p))))))
179 (cond
181 (and ;; (rpobj e) ;; g([a,b],[],y) = -g([a,u],[])*g([b,v],[])*g([],[u,v],y)
182 (boundp '$imetric)
183 (eq (caar e) $imetric)
184 (eql (length cov) 2)
185 (eql (length contr) 0)
186 (= (length deriv) 1)
188 (return
189 (prog (d1 d2)
190 (setq d1 ($idummy) d2 ($idummy))
191 (return
192 ($canform
193 (list '(mtimes simp)
194 (cond (csign 1) (t -1))
195 (list (cons $imetric '(simp))
196 (list '(mlist simp) (nth 0 cov) d1)
197 '((mlist simp))
199 (list (cons $imetric '(simp))
200 (list '(mlist simp) (nth 1 cov) d2)
201 '((mlist simp))
203 (append (list (cons $imetric '(simp))
204 '((mlist simp))
205 (list '(mlist simp) d1 d2)
207 deriv
210 nfprpobjs
217 (and ;; (rpobj e) ;; g([a,b],[],y,d)
218 (boundp '$imetric)
219 (eq (caar e) $imetric)
220 (eql (length cov) 2)
221 (eql (length contr) 0)
222 (= (length deriv) 2)
224 (return
225 (prog (d1 d2)
226 (setq d1 ($idummy) d2 ($idummy))
227 (return
228 ($canform
229 (list '(mplus simp)
230 (list '(mtimes simp)
231 (cond (csign 1) (t -1))
232 (list (cons $imetric '(simp))
233 (list '(mlist simp) (nth 0 cov) d1)
234 '((mlist simp))
235 (nth 1 deriv)
237 (list (cons $imetric '(simp))
238 (list '(mlist simp) (nth 1 cov) d2)
239 '((mlist simp))
241 (list (cons $imetric '(simp))
242 '((mlist simp))
243 (list '(mlist simp) d1 d2)
244 (nth 0 deriv)
247 (list '(mtimes simp)
248 (cond (csign 1) (t -1))
249 (list (cons $imetric '(simp))
250 (list '(mlist simp) (nth 0 cov) d1)
251 '((mlist simp))
253 (list (cons $imetric '(simp))
254 (list '(mlist simp) (nth 1 cov) d2)
255 '((mlist simp))
256 (nth 1 deriv)
258 (list (cons $imetric '(simp))
259 '((mlist simp))
260 (list '(mlist simp) d1 d2)
261 (nth 0 deriv)
264 (list '(mtimes simp)
265 (cond (csign 1) (t -1))
266 (list (cons $imetric '(simp))
267 (list '(mlist simp) (nth 0 cov) d1)
268 '((mlist simp))
270 (list (cons $imetric '(simp))
271 (list '(mlist simp) (nth 1 cov) d2)
272 '((mlist simp))
274 (list (cons $imetric '(simp))
275 '((mlist simp))
276 (list '(mlist simp) d1 d2)
277 (nth 0 deriv)
278 (nth 1 deriv)
282 nfprpobjs
289 (setq tensor (mysubst0 (append (list (car e)
290 (consmlist cov)
291 (consmlist contr))
292 (cond ($iframe_flag deriv)
293 (t (itensor-sort deriv) ))
294 ) e))
295 (cond (csign (setq tensor (neg tensor))))
296 (return tensor)))
298 (defun rename1 (e dummy) ;Renames dummy indices in a consistent manner
299 (sublis (cleanup0 dummy) e))
301 (defun cleanup0 (a)
302 (do ((b a (cdr b))
303 (n 1 (1+ n))
305 (dumx))
306 ((null b) l)
307 (setq dumx (intern (format nil "~a~d" $idummyx n)))
308 (cond ((not (eq dumx (car b)))
309 (setq l (cons (cons (car b) dumx) l))))))
311 (defun extract-elements (a b) ;Extracts the elements from B specified by the indices in
312 ;i.e. (EXTRACT-elements '(2 5) '(A B C D E F)) yields (B E)
313 (do ((a a) (b b (cdr b)) (n 1 (1+ n)) (l))
314 ((null a) (nreverse l))
315 (cond ((equal (car a) n)
316 (setq l (cons (car b) l) a (cdr a))))))
318 (defun inserts (a b c) ;Substitutes A into B with respect to the index
319 (do ((a a) (b b (cdr b)) (c c) (n 1 (1+ n)) (l)) ;specification C
320 ((null a) (nreconc l b))
321 (cond ((equal (car c) n)
322 (setq l (cons (car a) l) a (cdr a) c (cdr c)))
323 (t (setq l (cons (car b) l))))))
325 (defun symsort (l type)
326 (cond ((eq type '$sym) (sort l #'less)) ;SORT SYMmetric indices
327 ((eq type '$anti) (antisort l))
328 (t (cycsort l))))
330 (defun antisort (l) ;SORT ANTIsymmetric indices and set CSIGN as needed
331 ((lambda (q) (cond ((equal ($levi_civita (consmlist (mapcar 'cdr q))) -1)
332 (setq csign (not csign))))
333 (mapcar 'car q))
334 (sort (index l) #'less :key #'car)))
336 (defun index (l) ;(INDEX '(A B C)) yields ((A . 1) (B . 2) (C . 3))
337 (do ((l l (cdr l)) (n 1 (1+ n)) (q))
338 ((null l) (nreverse q))
339 (setq q (cons (cons (car l) n) q))))
341 (defun cycsort (l) ;SORT CYClic indices
342 ((lambda (n) (cond ((equal n 0) l)
343 (t (append (nthcdr n l)
344 (reverse (nthcdr (f- (length l) n)
345 (reverse l)))))))
346 (1- (cdr (least l)))))
348 (defun least (l) ;Returns a dotted pair containing the alphanumerically least
349 ;element in L in the car and its index in L in the cdr
350 (do ((l (cdr l) (cdr l)) (a (cons (car l) 1)) (n 2 (1+ n)))
351 ((null l) a)
352 (cond ((less (car l) (car a)) (setq a (cons (car l) n))))))
354 (declare-top (special free-indices))
356 (defun canprod (e f)
357 (prog (scalars indexed)
358 (cond
360 (catch 'foo
362 ((f (cdr e) (cdr f)) (obj))
364 (null f)
365 (setq scalars (nreverse scalars)
366 indexed (nreverse indexed)
370 (setq obj (car f))
371 (cond
372 ((atom obj) (setq scalars (cons obj scalars)))
373 ((rpobj obj) (setq indexed (cons obj indexed)))
374 ((eq (caar obj) 'mplus) (throw 'foo t))
375 (t (setq scalars (cons obj scalars)))
379 (return ($canform ($expand e) f))
381 ((null indexed) (return e))
383 (null (cdr indexed))
384 (return
385 (nconc (ncons '(mtimes)) scalars (ncons (canten (car indexed) t)))
389 (return
390 (nconc
391 (ncons '(mtimes))
392 scalars
393 (mapcar
394 (function (lambda (z) (canten z nil)))
396 (lambda (q)
397 (cond
398 (f (rename1 q
399 (nonumber
400 (cdaddr ($indices2 (cons '(mtimes) (reverse q))))
404 (t q)
407 (mapcar
408 #'cdr
410 (lambda (x) (cond ($flipflag (reverse x)) (t x)))
411 (sort
412 (progn
413 (setq free-indices (nonumber (cdadr ($indices e))))
414 (mapcar 'describe-tensor indexed)
416 #'tensorpred :key #'car
429 (defun tensorpred (x y)
430 (do ((x x (cdr x)) (y y (cdr y)) (a) (b))
431 ((null x))
432 (setq a (car x) b (car y))
433 (and (not (equal a b)) (return
434 (cond ((fixnump a) (< a b))
435 ((and (listp a) (listp b)) (tensorpred a b))
436 ((null a) t)
437 ((null b) nil)
438 (t (alphalessp a b)))))))
440 (defun describe-tensor (f)
441 (cons (tdescript f) f))
443 (defun tdescript (f)
444 (prog (name indices lcov lcontr lderiv)
445 (setq name (caar f)
446 indices (append (cdadr f) (cdaddr f) (cdddr f))
447 lcov (length (cdadr f))
448 lcontr (length (cdaddr f))
449 lderiv (length (cdddr f)))
450 (return (list (car (least (intersect indices free-indices)))
451 (f+ lcov (f+ lcontr lderiv) )
452 lderiv lcov name))))
454 (declare-top (unspecial free-indices))