share/tensor/canten.lisp

   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
 335 a   (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)))))))