share/tensor/atensor.mac

   1 /* Copyright (C) 2004 Viktor T. Toth <http://www.vttoth.com/>
   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  * Algebraic tensor manipulation
  14  *
  15  */
  16
  17 if get('atensor,'version)#false then error("ATENSOR already loaded!")$
  18
  19 alg_type: 'universal;
  20 asymbol:v;
  21 adim:0;
  22 aform:diagmatrix(3,1);
  23
  24 dotscrules:true;
  25 dotdistrib:true;
  26 dotexptsimp:false;
  27
  28 abasep(v):=if not atom(v) and mapatom(v) and op(v)=asymbol and length(v)=1
  29               and part(v,1)>0 and part(v,1)<=adim then true else false;
  30
  31 /* declare(..,scalarp) doesn't work on function names, which is why
  32    we need some special rules here. */
  33 scalarfunp(x):=if not atom(x) and (op(x)=nounify(sf) or op(x)=nounify(af))
  34                then true else scalarp(x);
  35 nonscfunp(x):=if atom(x) or (op(x)#nounify(sf) and op(x)#nounify(af))
  36               then not scalarp(x) else false;
  37
  38 matchdeclare([nonsc1,nonsc2], nonscfunp());
  39 matchdeclare([scalarf],scalarfunp());
  40 matchdeclare([abasev1,abasev2],abasep());
  41
  42 defrule(scalarfun1,nonsc1.scalarf,nonsc1*scalarf);
  43 defrule(scalarfun2,scalarf.nonsc2,nonsc2*scalarf);
  44
  45 defrule(grassmann1,nonsc1.nonsc1,0);
  46 defrule(grassmann2,nonsc1.nonsc2,
  47         if ordergreatp(nonsc1,nonsc2) then -nonsc2.nonsc1 else nonsc1.nonsc2);
  48 defrule(clifford1,nonsc1.nonsc1,sf(nonsc1,nonsc1));
  49 defrule(clifford2,nonsc1.nonsc2,
  50         if ordergreatp(nonsc1,nonsc2) then -nonsc2.nonsc1+2*sf(nonsc2,nonsc1) else nonsc1.nonsc2);
  51 defrule(symmetric,nonsc1.nonsc2,
  52         if ordergreatp(nonsc1,nonsc2) then nonsc2.nonsc1 else nonsc1.nonsc2);
  53 defrule(symplectic,nonsc1.nonsc2,
  54         if ordergreatp(nonsc1,nonsc2) then nonsc2.nonsc1-2*af(nonsc2,nonsc1)
  55         else nonsc1.nonsc2);
  56 defrule(lie_envelop,nonsc1.nonsc2,
  57         if ordergreatp(nonsc1,nonsc2) then nonsc2.nonsc1-2*av(nonsc2,nonsc1)
  58         else nonsc1.nonsc2);
  59 defrule(complex1,abasev1.abasev1,-1);
  60
  61 atenrules:[
  62  [grassmann,[grassmann1,grassmann2]],
  63  [clifford,[clifford1,clifford2]],
  64  [symmetric,[symmetric]],
  65  [symplectic,[symplectic]],
  66  [lie_envelop,[lie_envelop]]
  67 ];
  68
  69 /* Recursively simplify an atensor expression, with special attention
  70    paid to the fact that mnctimes can have multiple arguments;
  71    furthermore, we don't want dot products to end up as arguments
  72    to av(), af(), or sf(). The counter j is to ensure that the function
  73    always terminates. Admittedly, this is not a very efficient
  74    implementation and can use improvements, but it does provide results
  75    that are compatible with commercial MACSYMA. */
  76 atensimp(exp):=block
  77 (
  78   [i,j,exp1,exp2,exp3,done],
  79   if not mapatom(exp) then exp:map(atensimp,exp),
  80
  81   for j thru 10 do
  82   (
  83     done:false,
  84     exp1:exp,
  85     while not done do
  86     (
  87       done:true,
  88       if not atom(exp1) and op(exp1)="." and length(args(exp1))>2 then block
  89       (
  90         [argv:args(exp1)],
  91         for i thru length(argv)-1 do
  92         (
  93           exp2:"."(argv[i],argv[i+1]),
  94           exp3:atensimp(exp2),
  95           if exp2#exp3 then
  96           (
  97             exp1:apply(".",append(rest(argv,-length(argv)+i-1),[exp3],
  98                                   rest(argv,i+1))),
  99             exp1:ev(exp1,simp),
 100             i:length(argv),
 101             done:false
 102           )
 103         )
 104       )
 105     ),
 106     if exp#exp1 then exp:atensimp(exp1) else
 107       for i thru length(atenrules) do
 108         if atenrules[i][1]=alg_type then
 109     (
 110       exp2:apply('applyb1,cons(exp,atenrules[i][2])),
 111       exp2:ev(exp2,simp),
 112       exp2:applyb1(exp2,scalarfun1,scalarfun2),
 113       exp2:ev(exp2,simp),
 114       if exp=exp2 then j:10
 115       else exp:(if mapatom(exp2) then exp2 else map(atensimp,exp2))
 116     )
 117   ),
 118   exp
 119 );
 120
 121 init_atensor(algname, [optdims]):=
 122 (
 123   if algname='universal or algname='grassmann or algname='clifford or
 124      algname='symmetric or algname='symplectic or algname='lie_envelop then
 125   (
 126     alg_type:algname,
 127     if alg_type='symplectic then
 128     (
 129       adim:0,
 130       kill(aform),
 131       if length(optdims)>0 then
 132       (
 133         adim:optdims[1],
 134         aform:zeromatrix(optdims[1],optdims[1]),
 135         for i thru adim do
 136           for j thru i-1 do
 137         (
 138           aform[i,j]:if evenp(i+j) then 1 else -1,
 139           aform[j,i]:-aform[i,j]
 140         ),
 141         if length(optdims)>1 then
 142         (
 143           for i thru optdims[2] do
 144             aform:addrow(aform,makelist(0,j,1,adim)),
 145           adim:adim+optdims[2],
 146           for i thru optdims[2] do
 147             aform:addcol(aform,makelist(0,j,1,adim))
 148         ),
 149         if length(optdims)>2 then error("Invalid optional dimensions"),
 150         optdims:[]
 151       )
 152     )
 153     else if alg_type='clifford then block
 154     (
 155       [p:0,z:0,n:0],
 156       if length(optdims)>0 then p:optdims[1],
 157       if length(optdims)>1 then z:optdims[2],
 158       if length(optdims)>2 then n:optdims[3],
 159       if length(optdims)>3 then error("Invalid optional dimensions"),
 160       adim:p+z+n,
 161       if adim>0 then
 162       (
 163         aform:ident(adim),
 164         for i from p+1 thru p+z do aform[i,i]:0,
 165         for i from p+z+1 thru adim do aform[i,i]:-1
 166       )
 167       else kill(aform),
 168       optdims:[]
 169     )
 170     else if length(optdims)>0 then
 171     (
 172       if length(optdims)>1 then error("Invalid optional dimensions"),
 173       adim:optdims[1],
 174       aform:ident(adim),
 175       if alg_type='lie_envelop then
 176       (
 177         for i thru adim do for j thru i do
 178           if i=j then aform[i,j]:0
 179           else aform[i,j]:-(aform[j,i]:(remainder(2*adim+2-i-j,adim)+1)*
 180                               signum(i-j)*(if oddp(i+j) then 1 else -1))
 181       ),
 182       optdims:[]
 183     )
 184     else
 185     (
 186       kill(aform),
 187       adim:0
 188     )
 189   )
 190   else if algname='complex then init_atensor(clifford,0,0,1)
 191   else if algname='quaternion then init_atensor(clifford,0,0,2)
 192   else if algname='pauli then init_atensor(clifford,3)
 193   else if algname='dirac then init_atensor(clifford,3,0,1)
 194   else error("Unknown algebra name"),
 195   if optdims#[] then error("No optional dimensions permitted for this algebra"),
 196   done
 197 );
 198
 199 af(u,v):=if abasep(u) and abasep(v) then aform[part(u,1),part(v,1)]
 200          else 'af(u,v);
 201 sf(u,v):=if abasep(u) and abasep(v) then aform[part(u,1),part(v,1)]
 202          else 'sf(u,v);
 203 av(%u,%v):=if abasep(%u) and abasep(%v) then block
 204            (
 205              [i:aform[part(%u,1),part(%v,1)]],
 206              if abs(i)>0 and abs(i)<=adim then asymbol[abs(i)]*signum(i) else 0
 207            )
 208            else 'av(%u,%v);
 209
 210
 211 put('atensor,'v20081223,'version)$