/*
FileName: Clifford.mac
Contents: Maxima program to manupulate the Clifford algebra
Creation: 2022.12.8  by Kodama Hideo
2022/12/14: bug on the handling of "-" is corrected in CLsimplify
2022/12/16: maxima bug related to the confusion of global and local variables is coped with.
2022/12/16: version 1.00
Last modified: 2022/12/16
*/

/******>Initialization */

scalarmatrixp: false$
stringdisp: true$


/*****>Local setting*/
LibDir: pathname_directory(load_pathname)$

/*****> Library */
load(concat(LibDir,"SymbolicTensorCalculus_v1_00_pub.mac"))$
load(concat(LibDir, "mylib_dict_v1_00_pub.mac"))$


/******>Declaration */
declare( eta, symmetric)$  /* metric  eta[a,b] */
declare(Gamm, antisymmetric)$  /* Gamm[a b c ...] */
dim: 'dim$
infix( "&&" )$

/******>Commond defs */
Pauli[0]: matrix([1,0],[0,1])$
Pauli[1]: matrix([0,1],[1,0])$
Pauli[2]: matrix([0, -I],[I,0])$
Pauli[3]: matrix([1,0],[0,-1])$

/******>Job controle */
Job_Cliff: true$
  Job_S3_Fierz: false$
  Job_11DSugra: false$
    Job_Psi2DPsi: false$
      Job_Psi2DPsi_Z: false$
    Job_Psi3G: false$
      Job_Psi3G_dXdYdZ: false$
    Job_SC: false$       /* assuming Job_Psi2DPsi & Job_Psi3G are true */
Job_TPR: false$
  Job_10D: false$
  Job_10D_nullKilling: false$
  Job_10D_Spin7: false$
  Job_10D_Majorana: false$
 Job_4D: false$
  Job_4D_Majorana: false$
 Job_8DL: false$
  Job_8DL_Majorana: false$
 Job_8DE: false$
  Job_8DE_Real: false$
    Job_8DE_N8Sugra: false$

 
/******>Symbolic Clifford Calculus */
print("Job1,Symbolic Clifford Calculus");
print("Usage, e.g. CLsimplify(Gamm[a,b] && Gamm[c,d,e])");

/*****>Procedures */
/****>Defining the binary operator && */
/***>proc: isAbelian */
/*  isAbelian( eq ) => false  if  eq is monomial and contains non-Abelian name
    The name of a Clifford var or expression starts with G !!
*/
isAbelian(x) := block(
  if (not atom(x) ) and op(x)="-" then (
      x: -x
   ),  /*end if*/
   if ( symbolp(x) and substring(string(x), 1,2)="G" )
      or (not atom(x))  and ( substring(string(op(x)),1,2)="G" or op(x)="&&") then (
        return (false)
   ) else (
        return (true)
   ) /*end if*/
)$  /*end proc: isAbelian*/


/***>proc: && */
/* "&&"(X, Y) => non-Abelian product of eqs containg Clifford objects
*/
"&&"(X, Y) := block(
    [ x, y,  A1, NA1, A2, NA2],

    x: expand(X), y: expand(Y),
    if ( not atom(x) ) and op(x)="+" then (
        return(map(lambda([u], "&&"(u,  y)), x))
    ) elseif  (not atom(y)) and op(y)="+" then (
        return( map(lambda([u], "&&"(x,  u)), y) )
    ) elseif (not atom(x) ) and op(x)="-"  then (
        return (- ((-x) && y)) 
    ) elseif (not atom(y)) and op(y)="-"  then (
        return (-( x && (-y)) )
    ) else (
        X:  selectremove(isAbelian, x),
        A1: X[1], NA1: X[2], 
        Y:  selectremove(isAbelian, y),
        A2: Y[1], NA2: Y[2], 
/* display([NA1,NA2]),*/
        if NA1=1 or NA2=1 then (
            return( A1*A2*NA1*NA2)
        ) else (
             return( A1*A2*ev( "&&"(NA1, NA2), noeval) )
        ) /*end if*/
    )  /*end if*/
)$  /*end proc: &&*/


/****>simplifications  */
/***>Clifford simplifications */
/**>proc: isContractable */
/* isContractable(_ex) => true/false
*/
isContractable (_ex) := block( [_ex1, _result],
  
    _ex1: expand(_ex),
    if atom(_ex1) then (
        return (false)
    ),  /*end if*/
    
    _result: false,
    if  op(_ex1)="+"  then ( 
        for _term in _ex1 do (
            if  isContractable(_term)  then (
                _result : true, 
                break
            )  /*end if*/
        )  /*end do*/
    ) elseif  op(_ex1)="&&"  then (
        if op(args(_ex1)[1])=Gamm and op(args(_ex1)[2])=Gamm then (
            _result: true
        ) elseif  op(args(_ex1)[1])="&&" or op(args(_ex1)[2])="&&" then (
            _result: ev( isContractable(args(_ex1)[1]) or isContractable(args(_ex1)[2]), eval)
        )  /*end if*/
    ) elseif op(_ex1)="/" then (
        _result: isContractable(args(_ex1)[1])
    ) elseif  op(_ex1)="*" then (
        for _term in _ex1 do (
            if (not atom(_term)) and op(_term)="&&"  and isContractable(_term) then (
                _result: true
            )  /*end if*/
        )  /*end do*/
    ),  /*end if*/
    return ( _result )
)$  /*end proc: isContractable*/   
  
/**>proc: CLsimplify */
/*  CLsimplify(expr [, dim]): Gamm[a,b,...]&&Gamm[c,d,..] => Gamm[a,b,c,..]+ ...
     dim = dimension of the CL algebra
*/
CLsimplify( _ex0, [opt] ) :=  block(
    [ /*list*/
      /*misc*/
      _ex, _aa, _bb, _lena, _lenb, _inda, _indb,  _subind,_val,_ml,_gamf,
        _xx,_yy,__val1,_repl, _X
    ],
    local(_aa,_bb,_xx,_yy),
    
    _ex: expand(_ex0),
    if length(opt)>0 and integerp(opt[1]) and opt[1]>0 then (
        _ml: opt[1]
    ) else (
        _ml: 11  /* default */
    ),
    
    if atom(_ex) or atom(-_ex) then (
/* print("Case atom"),*/
        _val: _ex
    ) else (
        if  op(_ex)="+"  then (
/*print("Case +"),*/
            _val: 0,
            for term in _ex do (
                _val: _val+CLsimplify(term,_ml)
            )  /*end do*/
        ) elseif op(_ex)="-" then (
            _val : -CLsimplify(-_ex, _ml)
 /*       ) elseif op(_ex)="/" then (
            _val: CLsimplify(args(_ex)[1],_ml)/args(_ex)[2]
 */
        ) elseif  op(_ex)="*" then (
/*print("Case *"), */
            _val: 1,
            for term in _ex do (
                _val: _val*CLsimplify(term,_ml)
            )  /*end do*/
        ) elseif op(_ex)="&&"  then (
/*print("Case &&"),*/
            [_aa, _bb]: args(_ex),
 /*display([_aa,_bb]),*/
            if isContractable(_aa) then (
/*print("Case &&: _aa"),*/
                 _yy: CLsimplify(_aa, _ml),
/*print("_yy"=_yy),*/
                _xx:  "&&" (_yy,  _bb),
/* print("_xx"=_xx),*/
                _val:  ev(CLsimplify( _xx, _ml), e_val)
 /*,print("_val"=_val)*/
            ) elseif  isContractable(_bb) then (
                _xx: "&&"(_aa, CLsimplify(_bb,_ml)),
/* print("_xx"=_xx), */
                _val: ev( CLsimplify( _xx, _ml), e_val)
/* ,print("_val"=_val) */
            ) elseif  (not atom(_aa)) and op(_aa)=Gamm and
                        (not atom(_bb)) and op(_bb)=Gamm then (
/*print("Case &&: Gamm && Gamm"),*/
                _inda: args(_aa), _lena : length(_inda),
                _indb: args(_bb), _lenb : length(_indb),
                if _lena>_ml or _lenb >_ml then (
                    return( 0)
                ),  /*end if*/
                if _lena=1 then (
                    _val: arraymake(Gamm, append(_inda, _indb) ),
                    if  _lena+_lenb>_ml then (
                        _val: 0
                    ),  /*end if*/
 /* print([_lena, _lenb]=_val),*/
                    for i: 1 thru _lenb do (
                        _subind: append( slist(_indb, 1, (i-1)), slist(_indb, (i+1), _lenb) ),
                        if _subind=[] then (
                            _X : 1
                         ) else (
                            _X: arraymake(Gamm, _subind)
                         ), /*end if*/
                        _val: _val+(-1)^(i-1)*eta[_inda[1], _indb[i]]*_X
                    )  /*end do*/
                ) elseif  _lenb=1 then (
                    _val: arraymake(Gamm, append(_inda, _indb) ),
                    if _lena+_lenb>_ml then (
                        _val: 0
                    ),  /*end if*/
/*print([_lena, _lenb]=_val),*/
                    for  i: 1 thru  _lena do (
                        _subind: append( slist(_inda, 1, (i-1)), slist(_inda, (i+1), _lena) ),
                        if _subind=[] then (
                            _X : 1
                         ) else (
                            _X: arraymake(Gamm, _subind)
                         ), /*end if*/
                        _val: _val+(-1)^(_lena - i)*eta[_indb[1], _inda[i]]*_X
                    )  /*end do*/
                ) else (
                    if _lena <=_lenb  then (
                        _val: _bb,
                        for i: 1 thru _lena do (
                            _val: CLsimplify( (Gamm[_xx[_lena -i+1]]  && _val), _ml)
                        ),  /*end do*/
                        if _val=0 then ( 
                            return (0)
                        ),  /*end if*/
                        _yy: permlist1(_inda),
                        __val1: 0,
                        for j: 1 thru length(_yy) do (
                            _repl: makelist( _xx[k]=_yy[j][1][k], k, 1, _lena),
                            __val1: __val1+_yy[j][2]*subst(_repl, _val)
                        ),  /*end do*/
                        return ( expand(__val1/factorial(_lena)))
                    ) else (
                         _val: a, 
                         for i: 1 thru _lenb do (
                             _val: CLsimplify( (_val && Gamm[_xx[i]]),_ml)
                         ),  /*end do*/
                         if _val=0 then ( 
                             return( 0)
                         ),  /*end if*/
                         _yy: permlist1(_indb),
                         __val1: 0,
                         for j: 1 thru length(_yy) do (
                             _repl:  makelist( _xx[k]=_yy[j][1][k], k, 1, _lenb ),
                             __val1: __val1+_yy[j][2]*subst(_repl, _val) 
                         ),  /*end do*/
                         return (expand(__val1/factorial(_lenb)) )
                    )  /*end if*/ 
                )  /*end if*/
            ) else (
/* print("Case &&: impossible"),*/
                _val: _ex  /* should not happen!*/
            )  /*end if*/
        ) else (
 /*print("Case noncontractible"),*/
            _gamf: selecremove_has(Gamm, _ex),
            if _gamf[1]#1  and length(args(_gamf[1]))>_ml then (
                _val: 0
            ) else (
                _val: _ex
            )  /*end if*/
        )  /*end if*/
    ), /*end if*/
    return ( _val)
)$  /*end proc:CLsymplify*/


/***>sum simplifications */
/**>proc: issumindex */
/*  issumindex (x) => true/false
*/      
sumindexset: {"s"}$
issumindex(x) := ( member(substring(string(x), 1,2), sumindexset) 
)$  /*end proc: issumindex*/

/**>proc: gContraction */
/* gContraction(expr[, dim])
    eta[i,s]*eta[j,s]=eta[i,j], eta[s,s]=dim, eta[s1,s2]^2=dim
    T[..., s]*eta[s,a]=T[....,a]
*/
gContraction( _ex, [opt]) := block(
    [ _ex1, _dim0, _ex2, _aa,_bb,_X, _Y, _Y1, _Ylist, _Ylen, _yy, _zz],

    _ex1: expand(_ex),
    if  length(opt)>0 and integerp(opt[1]) and opt[1]>0 then (
        _dim0: opt[1]
    ) else ( 
        _dim0: 'dim
    ),  /*end if*/
    if (not atom(_ex1)) and op(_ex1)="+" then (
        return( map(lambda([x], gContraction(x,_dim0)), _ex1))
    ),  /*end if*/
    if (not atom(_ex1)) and op(_ex1)="-" then (
        return ( - gContraction(-_ex1, _dim0))
    ), /*end if*/
    if (not atom(_ex1)) and op(_ex1)="/" then (
        return( gContraction(args(_ex1)[1])/args(_ex1)[2] )
    ), /*end if*/
    
    /* _ex1= _Y * _X,  _Y=eta[*,*] eta[*,*] ... 
         => _Y1*_X1, _Y1 contains no sumindex */
    [_Y, _X]:  selectremove_has(eta, _ex1),
/* display([_Y, _X]),*/
    if _Y=1 then (
          return (_ex1)
    ),  /*end if*/    
    if op(_Y)="*" then (
          _Ylist: args(_Y)
    ) else (
          _Ylist: [_Y]
    ),  /*end if*/
/*   display(_ex1, _X, _Y, _Ylist), */

    _Y1: 1,  _Ylen: length(_Ylist),
    for ii: 1 thru _Ylen do ( 
          _yy: _Ylist[ii],
          if   op(_yy)="^" and args(_yy)[2]=2 and op(args(_yy)[1])='eta  then (      /* eta[a,b]^2*/
              [_aa, _bb] : args(args(_yy)[1]),
              if issumindex(_aa) or issumindex(_bb) then (
                  if _aa=_bb then (
                      _Y1: _Y1*_dim0^2
                  ) elseif  issumindex(_aa) and issumindex(_bb) then (
                      _Y1: _Y1*_dim0
                  ) elseif  issumindex(_aa) then (
                      _Y1: _Y1*eta[_bb,_bb]
                  ) else ( 
                      _Y1: _Y1*eta[_aa,_aa]
                  )  /*end if*/
              ) else (
                  _Y1: _Y1*_yy
              )  /*end if*/   
          ) elseif  op(_yy)='eta then (              /* eta[a,b] */
              [_aa, _bb] : args(_yy),
              if issumindex(_aa) or issumindex(_bb) then (
                  if _aa=_bb then (
                      _Y1: _Y1*_dim0
                  ) elseif  issumindex(_aa) then (
                      _zz: selectremove_has(_aa, _X*product(_Ylist[j], j, (ii+1), length(_Ylist))),
                      if _zz[1]#1 then (
                          _X: subst(_aa=_bb,_X),
                          _Ylist: subst(_aa=_bb,_Ylist)
                      ) else (
                          _Y1: _Y1*_yy
                      )  /*end if*/
                  ) elseif  issumindex(_bb) then (
                      _zz: selectremove_has(_bb, _X*product(_Ylist[j], j, (ii+1), length(_Ylist))) ,
                      if  _zz[1]#1 then (
                           _X: subst(_bb=_aa,_X),
                           _Ylist: subst(_bb=_aa,_Ylist)
                      ) else (
                           _Y1: _Y1*_yy
                      )  /*end if*/
                  ) else (
                          _Y1: _Y1*_yy
                  )  /*end if*/
              ) else (
                  _Y1: _Y1*_yy
              )  /*end if*/
          )  /*end if*/
     ),  /*end do*/
    return (_Y1*_X)
)$  /*end proc:gContraction*/


/**>proc: sumsimp */
/* sumsimp(_ex [, _sx]) : replace the sum indeces by the canonical symbols sx1, ..., sxn
    sx is a name with s as its first character
*/
sumsimp(_ex, [opt]) := block( 
    [ _sx, _ex1, _termlist, _indset1, _indset2,  _sumindlist] ,
  
    if length(opt)=0 then (
        _sx: 'ss
    ) else (
        _sx: opt[1]
    ),  /*end if*/
    _ex1: expand(_ex),
    if atom(_ex1) or atom(-_ex1) then (
        return( _ex1)
    ), /*end if*/
    if op(_ex1)="+"  then (
        map(sumsimp,_ex1)
    ) else (
        _indset1: [], _indset2: [],
        if op(_ex1)="*" then (
            _termlist: args(_ex1)
        ) else (
            _termlist: [ _ex1]
        ), /*end if*/
        for term in _termlist do (
            if not atom(term) and op(term)= "-" then (
                term: -term
            ),  /*end if*/
            if (not atom(term)) then (
                 if op(term)="^"  then (
                     if   (not atom(args(term)[1]) ) and op(args(term)[1])#"-"  then (
                         if op(args(term)[1])= 'Gamm  then (
                              _indset1: append(_indset1,  args(args(term)[1]) )
                          ) else (
                              _indset2:  append( _indset2,  args(args(term)[1]) )
                          )  /*end if*/
                     ) elseif (not atom(-args(term)[1])) and (not atom(args(term)[1])) then (  
                         if op(-args(term)[1])= 'Gamm   then (
                              _indset1: append(_indset1,  args(-args(term)[1]) )
                          ) else (
                              _indset2:  append( _indset2,  args(-args(term)[1]) )
                          )  /*end if*/
                     )  /*end if*/
                 ) elseif  (not atom(term))  and (not atom(-term) ) then (
                     if op(term)="-" then (
                          term: -term
                     ),  /*end if*/
                     if op(term)='Gamm  then (
                         _indset1: append( _indset1, args(term) )
                     ) else (
                         _indset2: append( _indset2, args(term) )
                     )  /*end if*/
                )  /*end if*/
            )  /*end if*/
        ),  /*end do*/
/*display(_indset1, _indset2),*/
        _indset1: setify(sublist(_indset1, issumindex)),
        _indset2: setify(sublist(_indset2, issumindex)),
        _sumindlist: union(_indset1, _indset2),
        _sumindlist: listify(_sumindlist),
        if length(_sumindlist)>0 then (
            ev(_ex1, makelist( _sumindlist[i]=concat(_sx, i) ,i, 1, length(_sumindlist)))
        ) else ( 
            _ex1
        )  /*end if*/    
    )  /*end if*/
)$  /*end proc: sumsimp*/ 

/****>expression helpers */
/***>proc: altsum */
/* f(a_1,...,a_n) -> 1/n! sum_{s} sign(s) f(a_s(1),...,a_s(n))
*/
altsum(_ex, indlist) := block( [_ex1,_ex2,_ex3,_permlist,_indnum,_tmplist],

    _ex1: expand(_ex),
    if (not atom(_ex1)) and op(_ex1)="+"  then (
        return ( map(lambda([x], altsum(x, indlist)), _ex1))
    ),  /*end if*/
    
    _permlist: listpermutations(indlist),
    _indnum: length(indlist),
    _ex2: ev(_ex1, makelist( indlist[i]=_tmplist[i], i, 1, _indnum)),
    _ex3: 0,
    for i: 1 thru _indnum! do (
        _ex3: _ex3
                +_permlist[i,2]*ev(_ex2,  makelist(_tmplist[j]=_permlist[i,1][j], j, 1, _indnum))
    ),  /*end do*/
    return (1/(_indnum!)*_ex3)
)$  /*end proc: altsum*/
  
/***>proc: Gsort */
/* Gsort(expr [,dim]) -> [expr0, expr1, .., exprn, exprtotal] 
     dim = dimension of  the Clifford algebra
     exprk =  the partial sums of the terms with the degree k
     exprtotal = a sorted  linear equation of Gamm[*] with respect to the rank of Gamma
*/
Gsort(ex0, [opt]) := block( 
    [/*list */
       _ranklist,
     /* array*/
       _exs,
     /* _dict */
       _dict,
     /* misc*/
      __dm, _rank,  _maxrank, _gg, _len, _jj
    ],
    local(_exs),

    _dm: 11,
    if length(opt)>0 then (
        _dm: opt[1],
        if not integerp(_dm) or _dm<0 then (
            error("dimension must be a non-negative integer!")
        )  /*end if*/
    ),  /*end if*/
    ex: expand(CLsimplify(ex0,_dm)),

    if not ( (not atom(ex)) and op(ex)="+" )  then (
        print(ex),
        return (ex)
    ) else (
         _ranklist: [],
         for term in args(ex) do (
/*display(term),*/
             _gg: selectremove_has(Gamm, term),
/*display(_gg),*/
             if  _gg[1]#1 then (
                 _rank: length(args(_gg[1])),
                 if member(_rank,_ranklist) then (
                     _exs[_rank]: _exs[_rank]+term
                 ) else (
                     _ranklist: append(_ranklist, [_rank]),
                     _exs[_rank]: term
                 )  /*end if*/
             ) else (
                 if member(0,_ranklist) then (
                     _exs[0]: _exs[0]+term
                 ) else (
                     _ranklist: append(_ranklist, [0]),
                     _exs[0]: term
                 )  /*end if*/
             )  /*end if*/
         ),  /*end do*/
         _ranklist: sort(_ranklist),
         _len: length(_ranklist),
          _maxrank: _ranklist[_len],
          if _maxrank > _dm then (
              _jj:1, 
              while _jj<=_len and _ranklist[_jj]<= _dm do ( 
                  _jj: _jj+1
              ),  /*end do*/
              _jj: _jj-1
          ) else (
              _jj: _len
          ),  /*end if*/
          _ranklist: slist(_ranklist, 1, _jj),
          print("_ranklist"=append(_ranklist, [all])),
          _dict: newdict([all]=[collectterms(ex,'Gamm, 'eta)]),
          for i in _ranklist do (
              adddict(_dict, [i]=[collectterms(_exs[i],'Gamm, 'eta)])
          ), /*end do*/
          return (copy(_dict))
    )  /*end if*/
)$  /*end proc:Gsort*/


/***>proc: swaptrf */
/*  swaptrf(ex,a,b) => ex1 = ex(a<->b)
*/
swaptrf( ex, a, b) := block( [c],

    if not symbolp(a) or not symbolp(b) then (
        error("swaptrf(ex,a,b), a and b must be names")
    ),  /*end if*/
    ex: subst(a=c,ex),
    ex: subst(b=a,ex),
    ex: subst(c=b,ex),
    return (ex)
)$  /*end proc*/


/******> End */
print("loading Clifford.mac --- done")$
print("loaded main functions: \"&&\",  CLsimplify, gContraction, sumsimp, Gsort")$