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. 12;; 13;; Comments: Canonical simplification for itensor.lisp 14;; 15 16; ** (c) Copyright 1979 Massachusetts Institute of Technology ** 17 18(in-package :maxima) 19 20(declare-top (special frei bouni $canterm breaklist smlist $idummyx)) 21 22(setq nodown '($ichr2 $ichr1 %ichr2 %ichr1 $kdelta %kdelta)) 23 24(defun ndown (x) (putprop x t 'nodown)) 25 26(defun ndownl (l) (mapcar 'ndown l)) 27 28(ndownl nodown) 29 30(setq breaklist nil $canterm nil) 31 32(defun brek (i) 33 (cond ((member i breaklist :test #'equal) t) )) 34 35; This package blindly rearranges the indices of RPOBJs, even those for 36; which such index mangling is forbidden, like our special tensors. To 37; avoid this, we use a private version of RPOBJ that excludes special 38; tensors like the Levi-Civita or Christoffel symbols 39 40(defun specialrpobj (e) 41 (cond ((or (atom e) (atom (car e))) nil) 42 (t (or (member (caar e) christoffels) 43 (eq (caar e) '$kdelta) (eq (caar e) '%kdelta) 44 (eq (caar e) '$levi_civita) (eq (caar e) '%levi_civita) 45 ) 46 ) 47 ) 48) 49 50(defun reallyrpobj (e) (cond ((specialrpobj e) nil) (t (rpobj e)))) 51 52;L IS A LIST OF FACTORS WHICH RPOBS SEPARATES INTO A LIST OF TWO LISTS. THE 53;FIRST IS A LIST OF THE RPOBECTS IN L. THE SECOND IS A LIST OF NON-RP OBJECTS 54 55(defun rpobs (l) 56 (do ( (x l (cdr x)) 57 (y nil (cond ((reallyrpobj (car x)) (append (list (car x)) y) ) 58 (t y) ) ) 59 (z nil (cond ((reallyrpobj (car x)) z) 60 (t (append (list (car x)) z)) ) ) ) 61 62 ( (null x) (cons y (list z))) )) 63 64;(defun name (rp) (cond ((rpobj rp) (caar rp) ) 65; (t (merror "not rpobject")))) 66;(defun conti (rp) (cond ((rpobj rp) (cdaddr rp)) 67; (t (merror "not rpobject")))) 68; 69;(defun covi (rp) (cond ((rpobj rp) (cdadr rp)) 70; (t (merror "not rpobject")))) 71; 72;(defun deri (rp) (cond ((rpobj rp) (cdddr rp)) 73; (t (merror "not rpobject")))) 74; 75;(defmfun $name (rp) (cond (($tenpr rp) (caar rp) ) ;test the name of tensor 76; (t (merror "not tenprect")))) 77;(defmfun $conti (rp) (cond (($tenpr rp) (cons smlist (cdaddr rp))) 78; (t (merror "not rpobject")))) ;test the contravariant indices 79; 80;(defmfun $covi (rp) (cond (($tenpr rp) (cons smlist (cdadr rp))) 81; (t (merror "not rpobject")))) ;test the contravariant indices 82; 83;(defun $deri (rp) (cond (($tenpr rp) (cons smlist (cdddr rp))) 84; (t (merror "not rpobject")))) 85 86(defun fre (l) (intersect l frei)) 87 88(defun boun (l) (intersect l bouni)) 89 90 91(defun grade (rp) 92 (+ (length (covi rp)) 93 (length (conti rp)) 94 (length (deri rp)))) 95 96;MON IS A MONOMIAL WHOSE "APPARENT" RANK IS ARANK 97 98(defun arank (mon) 99 (do ( (i 0 (+ (length (allind (car l))) i)) 100 (l (cdr mon) (cdr l) ) ) 101 ((null l) i) )) 102 103(defun eqarank (m1 m2) (= (arank m1) (arank m2))) 104 105;RP1 AND RP2 ARE RPOBJECTS WITH THE SAME GRADE 106 107(defun alphadi (rp1 rp2) 108 (alphalessp 109 (catenate (indlis rp1)) 110 (catenate (indlis rp2)))) 111 112 113(defun catenate (lis) (implode (exploden lis))) 114 115(defun indlis (rp) 116 (append (sort (append (fre (covi rp)) 117 (fre (conti rp)) 118 (fre (deri rp))) 'alphalessp) 119 (list (caar rp)) 120 (sort (append (boun (covi rp)) 121 (boun (conti rp)) 122 (boun (deri rp))) 'alphalessp))) 123 124;HOW TO USE ARRAY NAME AS PROG VARIABLE? 125 126(defun asort (l p) 127 (cond ((< (length l) 2) l) 128 (t 129 (prog (i j k az) 130 (setq i 0 j 0 k (length l) az (make-array k)) 131 (fillarray az l) 132 a (cond ((equal j (1- k)) 133 (return (listarray az))) 134 ((equal i (- k 1)) (setq i 0) (go a)) 135 ((apply p (list (aref az i) (aref az (1+ i)))) 136 (setq i (1+ i) j (1+ j)) (go a)) 137 ((apply p (list (aref az (1+ i)) (aref az i))) 138 (permute az i (1+ i)) 139 (setq i (1+ i) j 0) (go a) )) )))) 140 141(defun permute (arra i j) 142 (prog (x) (setq x (aref arra i)) 143 (setf (aref arra i) (aref arra j)) 144 (setf (aref arra j) x) )) 145 146(defun prank (rp1 rp2) (cond 147 ((> (grade rp1) (grade rp2)) t) 148 ((equal (grade rp1) (grade rp2)) (alphadi rp1 rp2)) )) 149 150 151(defun sa (x) 152 (sort (copy-list x) 'alphalessp)) 153 154(defun top (rp) (cdaddr rp)) 155(defun bot (rp) (append (cdadr rp) (cdddr rp))) 156(defun allind (rp) (cond ((not (reallyrpobj rp)) nil) 157 (t (append (cdadr rp) (cdaddr rp) (cdddr rp))))) 158 159;MON IS A MONOMIAL WHOSE FACTORS ARE ANYBODY 160;$IRPMON AND IRPMON RETURN LIST IS FREE AND DUMMY INDICES 161 162(defun $irpmon (mon) (prog (l cf dum free cl ci) 163 (setq l (cdr mon) 164 cf (car l) 165 dum nil 166 free nil 167 cl (allind cf) 168 ci nil ) 169 a (cond ((null l) (when (eq t (brek 19)) (break "19")) 170 (return (append (list smlist) 171 (la (list smlist) free) 172 (la (list smlist) dum) ) )) 173 ((null cl) (when (eq t (brek 18)) (break "18")) 174 (setq l (cdr l) cf (car l) cl (allind cf)) 175 (go a) ) 176 (t (setq ci (car cl)) 177 (cond ((not (member ci free :test #'eq)) (when (eq t (brek 17)) (break "17")) 178 (setq free (endcons ci free) 179 cl (cdr cl)) (go a) ) 180 (t (when (eq t (brek 16)) (break "16")) 181 (setq free (comp free (list ci)) 182 dum (endcons ci dum) 183 cl (cdr cl)) (go a) ) ) )))) 184 185(defun irpmon (mon) (prog (l cf dum free unio cl ci) 186 (setq l (cdr mon) 187 cf (car l) 188 dum nil 189 free nil 190 unio nil 191 cl (allind cf) 192 ci nil ) 193 a (cond ((null l) (when (eq t (brek 15)) (break "15")) 194 (setq free (comp unio dum) 195 dum (comp unio free)) 196 (return (append (list free) (list dum)) )) 197 198 ((null cl) (when (eq t (brek 14)) (break "14")) 199 (setq l (cdr l) cf (car l) cl (allind cf)) 200 (go a) ) 201 (t (setq ci (car cl)) 202 (cond ((not (member ci unio :test #'eq)) (when (eq t (brek 13)) (break "13")) 203 (setq unio (endcons ci unio) 204 cl (cdr cl)) (go a) ) 205 (t (when (eq t (brek 12)) (break "12")) 206 (setq dum (endcons ci dum) 207 cl (cdr cl)) (go a) ) ) )))) 208 209;THE ARGUMENT E IS A PRODUCT OF FACTORS SOME OF WHICH ARE NOT RPOBJECTS. THE 210;FUNCTION RPOBS SEPARATES THESE AND PUTS THEM INTO NRPFACT. THE RPFACTORS ARE 211;SORTED AND PUT IN A 212 213(defun redcan (e) 214 (prog (a b c d l nrpfact cci coi ct cil ocil) (when (eq t (brek 6)) (break "6")) 215 (setq nrpfact (cadr (rpobs (cdr e))) 216 d (irpmon e) 217 frei (car d) bouni (cadr d) 218 a (sort (append (car (rpobs (cdr e))) nil) 'prank) 219 l (length a) 220 b (make-array l) 221 c (make-array (list l 4))) 222 (fillarray b a) (when (eq t (brek 7)) (break "7")) 223 (do ((i 0 (1+ i))) 224 ((equal i l)) 225 (setf (aref c i 0) (name (aref b i))) 226 (setf (aref c i 1) (conti (aref b i))) 227 (setf (aref c i 2) (covi (aref b i))) 228 (setf (aref c i 3) (deri (aref b i)))) 229 230 231 (setq ocil (do ((i 0 (1+ i)) 232 (m nil (append (aref c i 3) m) ) ) 233 ((equal i l) m)) 234 ocil (append ocil (aref c 0 2) (car d)) 235 ct 0 236 cil (aref c ct 1) 237 cci (car cil) ) 238 (setf (aref c ct 1) nil) 239 240 a (cond 241 ((equal ct (1- l)) 242 (when (eq t (brek 1)) (break "1")) 243 (setf (aref c ct 1) cil) 244 (return 245 (append nrpfact 246 (do ((i 0 (1+ i)) 247 (lis (apdval c 0) (append lis (apdval c (1+ i))) ) ) 248 ((equal i (1- l)) lis ) ))) ) 249 250 ((zl-get (aref c ct 0) 'nodown) 251 (setf (aref c ct 1) cil) 252 (incf ct) 253 (setq cil (aref c ct 1) 254 ocil (append (aref c ct 2) ocil)) 255 (setf (aref c ct 1) nil) (go a)) 256 257 ((null cil) 258 (when (eq t (brek 2)) (break "2")) 259 (incf ct) 260 (setq cil (aref c ct 1) 261 ocil (append (aref c ct 2) ocil) ) 262 (setf (aref c ct 1) nil) (go a) ) 263 264 (t (setq cci (car cil)) (when (eq t (brek 5)) (break "5")) 265 (cond ((not (member cci ocil :test #'eq)) 266 (when (eq t (brek 3)) (break "3")) 267 (setq coi (do ((i (1+ ct) (1+ i) ) ) 268 ((member cci (aref c i 2) :test #'eq) i))) 269 270 (setf (aref c ct 2) 271 (cons cci (aref c ct 2))) 272 (setf (aref c coi 1) 273 (cons cci (aref c coi 1))) 274 (setf (aref c coi 2) 275 (comp (aref c coi 2) (list cci))) 276 (setq ocil (cons cci ocil) 277 cil (cdr cil) ) (go a) ) 278 (t (when (eq t (brek 4)) (break "4")) 279 (setf (aref c ct 1) (cons cci (aref c ct 1)) ) 280 (setq cil (cdr cil) ) (go a) ) )) ) ) ) 281 282(defun la (x y) (list (append x y))) 283 284(defun apdval (c i) (list (append (list (cons (aref c i 0) 285 (list 'simp))) 286 (la (list smlist) 287 (sa (aref c i 2))) 288 (la (list smlist) 289 (sa (aref c i 1))) 290 (sa (aref c i 3) )))) 291(defun canform (e) 292 (cond ((atom e) e) 293 ((rpobj e) e) 294 ((and (eq (caar e) 'mtimes) 295 (= 0 (length (car (rpobs (cdr e))))) ) e) 296 ((eq (caar e) 'mtimes) 297 (cons '(mtimes) (redcan e)) ) 298 ((eq (caar e) 'mplus) 299 (mysubst0 300 (simplus (cons '(mplus) 301 (mapcar #'(lambda (v) (cons '(mplus) (canarank v))) 302 (strata (cdr e) 'eqarank) )) 303 1. nil) e) ) 304 (t e) )) 305 306 307(defun endcons (x l) (reverse (cons x (reverse l)))) 308 309(defun comp (x y) 310 (do ((z (cond ((atom y) (ncons y)) (y)));patch for case when Y is not a list 311 (l x (cdr l)) 312 (a nil (cond ((member (car l) z :test #'equal) a) 313 (t (endcons (car l) a)) ))) 314 ((null l) a) ) ) 315 316(defun apdunion (x y) 317(do ((l y (cdr l)) 318 (a x (cond ((member (car l) a :test #'equal) a) 319 (t (endcons (car l) a)) ))) 320 ((null l) a) )) 321 322;LIST IS A LIST OF CANFORMED MONOMIALS OF THE SAME ARANK 323;CANARANK FINDS THE EQUIVALENT ONES 324 325(defun canarank (lis) (prog (a b c d ct cnt i) 326 (setq a lis 327 b nil 328 c (cdr a) 329 d (make-array (length a)) 330 ct (canform (car a)) 331 cnt (canform (car c)) 332 i 1) 333 (fillarray d a) 334 335a (cond ((null a) 336 (return b)) 337 338 ((and (null (cdr a)) (null c)) 339 (setq b (cons ct b)) 340 (return b) ) 341 342 343 ((null c) (when (eq t (brek 9)) (break "9")) 344 (setq b (cons ct b)) 345 (setf (aref d 0) nil) 346 (setq a (comp (listarray d) (list nil)) 347 c (cdr a) 348 i 1 349 ct (canform (car a)) 350 cnt (canform (car c)) ) 351 (cond ((null a) (return b)) 352 (t (setq d (make-array (length a))) 353 (fillarray d a))) 354 (go a)) 355 356 ((samestruc ct cnt) (when (eq t (brek 10)) (break "10")) 357 (setq b (cons (canform (transform cnt ct)) b)) 358 (setf (aref d i) nil) 359 (setq c (cdr c) 360 cnt (canform (car c)) 361 i (1+ i) ) (go a) ) 362 363 (t (when (eq t (brek 11)) (break "11")) 364 (setq c (cdr c) 365 cnt (canform (car c)) 366 i (1+ i)) (go a)) ))) 367 368;M1,M2 ARE (CANFORMED) MONOMIALS WHOSE INDEX STRUCTURE MAY BE THE SAME 369 370(defun samestruc (m1 m2) (equal (indstruc m1) (indstruc m2))) 371 372;MON IS (MTIMES) A LIST OF RP AND NON-RP FACTORS IN A MONOMIAL. THE NEXT 373;FUNCTION RETURNS A LIST WHOSE ELEMENTS ARE 4-ELEMENT LISTS GIVING THE NAME 374;(MON) AND THE LENGTHS OF THE COVARIANT,CONTRAVARIANT,DERIVATIVE INDICES. 375 376(defun indstruc (mon) 377 (do ( (l (cdr mon) (cdr l)) 378 (b nil (cond ((reallyrpobj (car l)) 379 (append b (list (findstruc (car l))) )) 380 (t b) )) ) 381 ( (null l) b) ) ) 382 383 384 385;FACT IS AN RP FACTOR IN MON. HERE WE FIND ITS INDEX STRUCTURE 386 387(defun findstruc (fact) 388 (append (list (caar fact) ) 389 (list (length (cdadr fact))) 390 (list (length (cdaddr fact))) 391 (list (length (cdddr fact))) )) 392 393;M1,M2 ARE MONOMIALS WITH THE SAMESTRUC TURE. THE NEXT FUNCTION TRANSFORMS 394;(PERMUTES) THE FORM OF M1 INTO M2. 395 396(defun transform (m1 m2) 397 (sublis (findperm m1 m2) m1)) 398 399(defun strata (lis p) 400 (cond ((or (null lis) (null (cdr lis))) (list lis)) 401 (t 402 403 (prog (l bl cs) (setq l lis cs nil bl nil) 404 405 a (cond ((null l) (when (eq t (brek 1)) (break "1")) 406 (return (cond ((null cs) bl) 407 (t (endcons cs bl))))) 408 409 ((and (null (cdr l)) (null cs)) (when (eq t (brek 2)) (break "2")) 410 (setq bl (endcons (list (car l)) bl)) 411 (return bl) ) 412 ((and (null (cdr l)) (not (null cs))) (when (eq t (brek 3)) (break "3")) 413 (return (cond ((funcall p (car l) (car cs)) 414 (setq cs (cons (car l) cs) 415 bl (endcons cs bl))) 416 (t (setq bl (endcons (list (car l)) (endcons cs bl))))))) 417 418 ((null cs) (when (eq t (brek 4)) (break "4")) 419 (setq cs (list (car l)) l (cdr l)) (go a)) 420 ((funcall p (car l) (car cs)) (when (eq t (brek 5)) (break "5")) 421 (setq cs (cons (car l) cs) 422 l (cdr l)) (go a) ) 423 424 (t (when (eq t (brek 6)) (break "6")) 425 (setq bl (endcons cs bl) 426 cs (list (car l)) 427 l (cdr l) ) (go a) ) ) )))) 428 429 430 431(defun tindstruc (mon) 432 (do ( (l (cdr mon) (cdr l)) 433 (b nil (cond ((reallyrpobj (car l)) 434 (append b (tfindstruc (car l)) )) 435 (t b) ))) 436 ((null l) b))) 437 438(defun tfindstruc (fact) 439 (append (cdadr fact) (cdaddr fact) (cdddr fact) )) 440 441(defun dumm (x) (equal (cadr (explodec x)) $idummyx)) 442 443 444(defun findpermut (i1 i2) 445 (comp (mapcar 'pij i1 i2) (list nil))) 446 447(defun pij (x y) 448 (cond ((and (dumm x) (dumm y) (not (eq x y))) (cons x y)))) 449 450 451;(SAMESTRUC M1 M2) IS T FOR THE ARGUMENTS BELOW 452;THE RESULTING PERMUTATION IS GIVEN TO TRANSFORM 453 454(defun findperm (m1 m2) 455 (do ((d1 (cadr (irpmon m1)) ) 456 (d2 (cadr (irpmon m2)) ) 457 (i1 (tindstruc m1) (cdr i1) ) 458 (i2 (tindstruc m2) (cdr i2) ) 459 (l nil (cond ((and (member (car i1) d1 :test #'eq) (member (car i2) d2 :test #'eq) 460 (not (eq (car i1) (car i2))) 461 (not (member (car i1) (car l) :test #'eq)) 462 (not (member (car i2) (cadr l) :test #'eq)) ) 463 (cons (endcons (car i1) (car l)) 464 (list (endcons (car i2) (cadr l))) ) ) 465 (t l)))) 466 467 ((null i1) (mapcar 'cons 468 (append (car l) (comp d1 (car l))) 469 (append (cadr l) (comp d2 (cadr l))))))) 470 471(defun $canten (x) 472 (cond ((not $allsym) 473 (merror "canten works only if allsym:true has been set")) 474 (t 475 (do ((i ($nterms x) ($nterms l)) 476 (l (canform x) (canform l))) 477 ((= i ($nterms l)) l) 478 (cond ((eq $canterm t) (print i))))))) 479 480(defun $concan (x) 481 (cond ((not $allsym) 482 (merror "concan works only if allsym:true has been set")) 483 (t 484 (do ((i ($nterms x) ($nterms l)) 485 (l (canform x) ($contract (canform l)))) 486 ((= i ($nterms l)) l) 487 (cond ((eq $canterm t) (print i))))))) 488