(*******************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which
typically will have the same name as this file except ending in
".nb" instead of ".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent
to using the Evaluate Initialization Cells menu command in the front
end.

DO NOT EDIT THIS FILE.  This entire file is regenerated
automatically each time the parent Notebook file is saved in the
Mathematica front end.  Any changes you make to this file will be
overwritten.
***********************************************************************)







If[MemberQ[$Packages,"MSTA2`"],
    Unprotect["MSTA2`*"];
    Clear["MSTA2`*"];
    Unprotect[Power,Exp,Times,Expand,Simplify,FullSimplify];
    Clear[Power,Exp,Times,Expand,CenterDot,Wedge,Simplify,FullSimplify];
    Protect[Power,Exp,Times,Expand,Simplify,FullSimplify];
    DeleteCases[$Packages,"MSTA2`"];
    ];

BeginPackage["MSTA2`"];




\[Gamma]::usage="\[Gamma] is the wrapper symbol for blade indices.";
setMSTA::usage="setMSTA[n_Integer, d_Integer] sets the number of particles to n and the number of spatial dimensions to d.";\

rev::usage="rev[x_] takes an expression x and reverses the order of blade indices.";\

dotprod::usage="dotprod[x_, y_] takes any two expressions x and y, which can be linear combinations of blades, and returns the geometric dot product x\[CenterDot]y. The infix operator [ESC].[ESC] is an alias.";\

wedgeprod::usage="wedgeprod[x_, y_] takes any two expressions x and y, which can be linear combinations of blades, and returns the geometric wedge product x\[Wedge]y. The infix operator [ESC]^[ESC] is an alias.";\

gather::usage="gather[x_] takes any expression x and groups together the coefficients of distinct blades.";\

comprod::usage="comprod[x_, y_] computes the commutator product of x and y (with a factor of 1/2).";\

listgrades::usage="listgrades[x_] lists the spectrum of grades present in the expression x.";\

project::usage="project[x_, n_Integer] projects out the n-grade part of the expression x.";\

shiftind::usage="shiftind[x_, n_Integer] shifts all indices of every blade occuring in x by adding an integer n.";\

Exp::usage="Exp[x_] is the exponential function.\n Exp[x_] has been overloaded so that exponentials of monoblades \"a \[Gamma][__]\" can be computed.";\

order::usage="order[x_] transforms blades by anticommuting indices until they are in increasing order.";\

red::usage="red[x_] cancels out duplicate blade indices of an expression x.";
Times::usage="\"x*y*z\" or \"x y z\" represents a product of terms.\n Times[x__] has been overloaded such that products between blades \[Gamma][__] are correctly executed.";\

whatMSTA::usage="whatMSTA[] returns the number of particles and space dimensions of the currently defined algebra.";\

listblades::usage="listblades[x_] lists the spectrum of blades present in x.";\

monobladeQ::usage=
    "monobladeQ[x_] returns True if x is a sum of terms each containing only one blade.";\

listcoeff::usage=
    "listcoeff[x_] returns the scalar part and blade coefficients in canonical order, together with the corresponding list of blades.";\

Expand::usage="Expand[x_] has been overloaded to expand out expressions involving blade elements.\n
    Expand[x_] expands out products and positive integer powers in x. Expand[x_, patt] leaves unexpanded any parts of expr that are free of the pattern patt.";\

Power::usage="\"x^y\" gives x to the power y.\nPower has been overloaded to deal with GA expressions correctly.";\

setMetric::usage=
    "setMetric[g_List] takes a list of the diagonal elements of the metric tensor. One can specify the signature for arbitrary dimensions {\[Gamma][0],\[Gamma][1],...,\[Gamma][N-1]}.";



Begin["`Private`"];





Unprotect@Times;
ClearAttributes[Times,Orderless];
Times[a__/;Not@FreeQ[{a},\[Gamma]],b__/;FreeQ[{b},\[Gamma]],c___]:=
    Times[b,a,c];
Times[a__/;(FreeQ[{a},\[Gamma]]&& Not@OrderedQ[{a}]),b___]:=
    Times[Sequence@@Sort[{a}],b];
Protect@Times;



Unprotect@Expand;
Expand[x_/;Not@FreeQ[x,\[Gamma]]]:=Distribute[x,Plus,Times];
Protect@Expand;





Unprotect[Power];
Power[x_/;Not@FreeQ[x,\[Gamma][__]],n_Integer]:=
  Module[{y,f},
    y=Expand[x];
    Switch[EvenQ[n],
      True,
      If[n\[Equal]0,Return[1],
        Power[Composition[gather,order,red][
            Distribute[f[y,y]]/.f\[Rule]Times],n/2]],
      False,
      If[n\[Equal]1,Return[x],
        Power[Composition[gather,order,red][
              Distribute[f[y,y]]/.f\[Rule]Times],(n-1)/2]x]
      ]
    ]
Protect[Power];

\!\(\(Unprotect@Power;\)\[IndentingNewLine]
  \(Power[E, \[Gamma][x__]] := \[IndentingNewLine]Switch[\[Gamma][x]\^2, \[IndentingNewLine]\[Gamma][], Cosh[1] \[Gamma][] + Sinh[1] \[Gamma][x], \[IndentingNewLine]\(-\[Gamma][]\), Cos[1] + Sin[1] \[Gamma][x]\[IndentingNewLine]];\)\[IndentingNewLine]
  Power[E, Times[\[Lambda]__ /; FreeQ[{\[Lambda]}, \[Gamma]], \[Gamma][x__]]] := Switch[\[Gamma][x]\^2, \[Gamma][], Cosh[Times[\[Lambda]]] \[Gamma][] + Sinh[Times[\[Lambda]]] \[Gamma][x], \(-\[Gamma][]\), Cos[Times[\[Lambda]]] + Sin[Times[\[Lambda]]] \[Gamma][x]]\[IndentingNewLine]
  \(Protect@Power;\)\)



Unprotect[Simplify,FullSimplify];
Simplify::blades=
    FullSimplify::blades="Simplify and FullSimplify do not operate correctly on expressions containing blades.";\

Simplify[x_/;Not@FreeQ[x,\[Gamma]],___]:=(Message[Simplify::blades];
      Abort[];);
FullSimplify[x_/;Not@FreeQ[x,\[Gamma]],___]:=(Message[FullSimplify::blades];
      Abort[];);
Protect[Simplify,FullSimplify];



setMSTA[n_Integer,dim_Integer]:=
  (
    numParticles=n;
    spacedim=dim;
    maxIndex=(spacedim+1)numParticles-1;
    metricsign=Table[If[Mod[k,spacedim+1]\[Equal]0,1,-1],{k,0,maxIndex}];
    Print["setMSTA:\n\t",n," particle\n\t",dim," space dimension"];
    )

whatMSTA[]:=
    Print["MSTA:\n\t",numParticles," particle\n\t",spacedim,
      " space dimension"];





setMetric[diag_List]:=
  (
    maxIndex=Length[diag]-1;
    metricsign=diag;
    Print["setMetric:\n\t metric set to ",diag];
    )



reduce::index="Error: index is out of range. Maximum allowed is `1`.";
reduce::noninteger="Error: the index `1` is a non integer.";



reduce[blade_\[Gamma]]:=
  (
    If[MemberQ[blade,n_/;n>maxIndex],Message[reduce::index,maxIndex];
      Abort[]];
    If[Length@#\[NotEqual]0,Message[reduce::noninteger,First@#];Abort[]]&@
      Select[blade,Not@IntegerQ[#]&];
    If[Signature[blade]\[NotEqual]0,(*no duplicates*) Return@blade];
    
    Times@@Fold[
        CancelDup,
        {+1,blade},
        Union@Flatten@DeleteCases[Split@Sort@(List@@blade),{_}]
        ]
    )



CancelDup[blade_,n_]:=
  NestWhile[CancelPair,{n,blade},(Count[#[[2,2]],n]>1)&][[2]]



CancelPair[{n_,{sign_,blade_}}]:=
  {n,{sign*metricsign[[n+1]]*(-1)^First[#[[2]]-#[[1]]-1],
          Delete[blade,Take[#,2]]}&@Position[blade,n]}

red=#/.\[Gamma][x__]\[RuleDelayed] reduce@\[Gamma][x] &;





\[Gamma]/:\[Gamma][x___]*\[Gamma][y___]:=reduce[Join[\[Gamma][x],\[Gamma][y]]]



rev[y_]:=
  Module[{HAHA},
    y/.HoldPattern[Times[x__]]\[RuleDelayed]HAHA[x]/.HAHA[x__]\[RuleDelayed]
            Reverse@HAHA[x]/.\[Gamma][x__]\[RuleDelayed]
          Reverse@\[Gamma][x]/.HAHA\[Rule]Times
    ]





orderblade::duplicate="The input blade `1` contains duplicate indices which first need to be eliminated.";

orderblade[blade_\[Gamma]]:=
  (
    If[Length@Union@blade\[NotEqual]Length@blade,
      Message[orderblade::duplicate,blade];Abort[];];(Signature*Sort)@blade//
      Through
    )



order=(#/.\[Gamma][x__]\[RuleDelayed]orderblade@\[Gamma][x])&;





gather=(#/.\[Gamma][__]\[Rule]0)+Coefficient[#1,#2].#2&[#,
        Union@Cases[#,\[Gamma][__],\[Infinity]]]&;



\!\(\(comprod = \(#1\ #2 - #2  #1\)\/2 &;\)\)





listgrades[x_]:=
    Union[Length/@Cases[{x//reduce},\[Gamma][__],Infinity]]//
      If[ToString[x/.\[Gamma][__]\[Rule]0]\[Equal]"0",#,Prepend[#,0]]&;



listblades[x_]:=
    If[(x/.\[Gamma][__]\[Rule]0)===0,
      Cases[x,\[Gamma][__],\[Infinity]]//Union,
      Prepend[Union@Cases[x,\[Gamma][__],\[Infinity]],1]
      ];



listcoeff::monoblade="The input expression needs to be put into monoblade form.";

listcoeff[x_]:=
    Module[{scalar,blades,coeffs,x$},
      If[!monobladeQ[x],Message[listcoeff::monoblade];Abort[];];
      x$=x//order;
      scalar=x$/.\[Gamma][__]\[Rule]0;
      blades=Cases[x$,\[Gamma][__],\[Infinity]]//Union;
      coeffs=Coefficient[x$,blades];
      If[scalar===0,
        {coeffs,blades},
        MapThread[Prepend,{{coeffs,blades},{scalar,\[Gamma][]}}]
        ]
      ];







project[x_,n_Integer]:=
    If[n\[Equal]0,
          #/.\[Gamma][__]\[Rule]0,
          #-(#/.\[Gamma]@@Table[_,{n}]\[Rule]0)
          ]&@red@x;





shiftind::index="Error: index will go out of range. Maximum allowed is `1`.";

shiftind[x_,n_Integer]:=
  (
    If[
          MemberQ[#,p_Integer/;Not[0\[LessEqual]p+n\[LessEqual]maxIndex]],
          Message[shiftind::index,maxIndex];Abort[];
          ]&/@Cases[x,\[Gamma][__],\[Infinity]];
    x/.\[Gamma][y__]\[RuleDelayed]\[Gamma][Sequence@@({y}+n)]
    )





dotprod[x_,y_]:=project[x*y,Abs[Subtract@@Flatten[listgrades/@{x,y}]]];



dotprod[x_,y_+z_]:=dotprod[x,y]+dotprod[x,z]
dotprod[x_+y_,z_]:=dotprod[x,z]+dotprod[y,z]



wedgeprod[x_,y_]:=project[x*y,Plus@@Flatten[listgrades/@{x,y}]]

wedgeprod[x_,y_+z_]:=wedgeprod[x,y]+wedgeprod[x,z]
wedgeprod[x_+y_,z_]:=wedgeprod[x,z]+wedgeprod[y,z]



monobladeQ[x_]:=
  If[Head[x]===Plus, And@@((Count[#,\[Gamma][__],\[Infinity]]<2)&/@List@@x),
    Count[x,\[Gamma][__],\[Infinity]]<2
    ]



package$defined=True;
End[];

EndPackage[];





Print["
MSTA2:\n\tCreated by S. Furuta.\n\tA package for manipulating multiparticle spacetime algebras.\n\t{+,-,-,-} metric set by default.
"];





setMSTA[1,1];



CenterDot=dotprod;
Wedge=wedgeprod;





SetAttributes[\[Gamma],NHoldAll];

\[Gamma][]=1;
Protect[\[Gamma]];



Times::Orderless="Warning: Times has spuriously turned orderless. Autofix: Attribute \"Orderless\" has been removed.";

$Post=(If[MemberQ[Attributes[Times],Orderless],Message[Times::Orderless];
          Unprotect@Times;ClearAttributes[Times,Orderless];Protect@Times];#)&;