(* CliffMath10.nb  Procedures for Clifford algbebra computations in Mathematica *)
(* May 2010 G. Stacey Staples *)
 
(* Use CircleDot for Clifford multiplication *)
(* x⊙y is entered by typing x[ESC]c.[ESC]y *)

(* Use Wedge for exterior product *)
(* x⋀y is entered by typing x[ESC]^[ESC]y *)

BeginPackage["CliffMath10`"];

e=Symbol["e"];
Protect[e];

SetSignature[p_,q_,r_:0]:=Module[{i,j},
Unprotect[MaxIndex];
MaxIndex=p+q+r;
Protect[MaxIndex];
Unprotect[PseudoScalar];
PseudoScalar:=(Subscript[e, Table[i,{i,1,MaxIndex}]]);
Protect[PseudoScalar];
Unprotect[pPart, pSet,qPart,qSet, rPart,rSet];   
pPart=p; qPart=q;rPart=r;
pSet = Table[i,{i,1,pPart}];
qSet=Table[i,{i,pPart+1,pPart+qPart}];
rSet=Table[i,{i,pPart+qPart+1,MaxIndex}];
Protect[pPart, pSet,qPart,qSet,rPart,rSet];

Unprotect[CircleDot];
ClearAll[CircleDot];   
SetAttributes[CircleDot,{Flat, OneIdentity, Listable}];

CircleDot[x_. (Subscript[e, a_]),y_. (Subscript[e, b_])]:=If[Length[Intersection[a⋂b,rSet]]≥1,0,(((-1))^(Length[Intersection[a⋂b,qSet]]+Ceiling[FractionalPart[(1/2) Sum[Sum[If[a[[i]]>b[[j]],1,0],{i,1,Length[a]}],{j,1,Length[b]}]]]))  x y (Subscript[e, (a⋃b)⋂(Complement[Union[a,b],Intersection[a,b]])])/.{(Subscript[e, {}])->1}];

CircleDot[x_, y_]:=x y /;(FreeQ[x,(Subscript[e, _])]∨FreeQ[y,(Subscript[e, _])]);

CircleDot=Symbol["CircleDot"];
Protect[CircleDot];

Unprotect[ClExpand];
ClearAll[ClExpand];
SetAttributes[ClExpand,Listable];
ClExpand[x_ + y_]:= Simplify[ClExpand[x] + ClExpand[y]];
ClExpand[x_?NumericQ  arg1_]:=x ClExpand[arg1];
ClExpand[arg1_]:= Module[{arg2},
arg2 = Distribute[ExpandAll[arg1],Plus, CircleDot];
If[!FreeQ[arg2,_⊙_] && arg2  !=  arg1,
arg2 =  ClExpand[arg2];
];
Return[arg2];
];
Protect[ClExpand];


(* Define wedge product *)

Unprotect[Wedge];
ClearAll[Wedge];
SetAttributes[Wedge, {Flat, OneIdentity, Listable}];
Wedge[x_. (Subscript[e, a_]),y_ ]:=x y (Subscript[e, a])/;FreeQ[y,_. (Subscript[e, _])];
Wedge[x_,y_. (Subscript[e, b_])]:=x y (Subscript[e, b])/;FreeQ[x,_. (Subscript[e, _])];
Wedge[x_,y_]:=x y /;FreeQ[x,_. (Subscript[e, _])]&&FreeQ[y,_. (Subscript[e, _])];
Wedge[x_.(Subscript[e, a_]),y_. (Subscript[e, b_])]:=If[Length[a⋂b]≠0,0,(((-1))^(Ceiling[FractionalPart[(1/2) Sum[Sum[If[a[[i]]>b[[j]],1,0],{i,1,Length[a]}],{j,1,Length[b]}]]])) x y  (Subscript[e, (a⋃b)])/.{(Subscript[e, {}])->1}];
Wedge[x_, y_]:=x y/;(FreeQ[x,(Subscript[e, _])]∨FreeQ[y,(Subscript[e, _])]);

Wedge=Symbol["Wedge"];
Protect[Wedge];
Unprotect[WExpand];
ClearAll[WExpand];
SetAttributes[WExpand,Listable];
WExpand[x_ + y_]:= Simplify[WExpand[x] + WExpand[y]];
WExpand[arg1_]:= Module[{arg2},
arg2 = Distribute[ExpandAll[arg1],Plus, Wedge];
If[!FreeQ[arg2,_⋀_] && arg2  !=  arg1,
arg2 =  WExpand[arg2];
];
Return[arg2];
];

Protect[WExpand];


Unprotect[ClExpand];
ClearAll[ClExpand];
SetAttributes[ClExpand,Listable];
ClExpand[x_ + y_]:= Simplify[ClExpand[x] + ClExpand[y]];
ClExpand[x_?NumericQ  arg1_]:=x ClExpand[arg1];
ClExpand[arg1_]:= Module[{arg2},
arg2 = Distribute[ExpandAll[arg1],Plus, CircleDot];
If[!FreeQ[arg2,_⊙_] && arg2  !=  arg1,
arg2 =  ClExpand[arg2];
];
Return[arg2];
];
Protect[ClExpand];];



(* Defines Clifford algebra with 2n generators.  Even index -> square to 1.  Odd index -> square to -1. *)
SetAlternatingSignature[n_]:=Module[{i,j,negtable},
Unprotect[MaxIndex];
MaxIndex=2n;
Protect[MaxIndex];
Unprotect[PseudoScalar];
PseudoScalar:=(Subscript[e, Table[i,{i,1,MaxIndex}]]);
Protect[PseudoScalar];
Unprotect[pPart,qPart,negtable];
pPart=n; qPart=n;
negtable=Table[2j-1,{j,1,n}];
Protect[pPart,qPart,negtable];
Unprotect[CircleDot];
ClearAll[CircleDot];  
SetAttributes[CircleDot,{Flat, OneIdentity, Listable}];

CircleDot[x_. (Subscript[e, a_]),y_. (Subscript[e, b_])]:=(((-1))^(Length[a⋂b⋂negtable]+Ceiling[FractionalPart[(1/2) Sum[Sum[If[a[[i]]>b[[j]],1,0],{i,1,Length[a]}],{j,1,Length[b]}]]]))  x y (Subscript[e, (a⋃b)⋂(Complement[Union[a,b],Intersection[a,b]])])/.{(Subscript[e, {}])->1};

CircleDot[x_, y_]:=x y/;(FreeQ[x,(Subscript[e, _])]∨FreeQ[y,(Subscript[e, _])]);

CircleDot=Symbol["CircleDot"];
Protect[CircleDot];
Unprotect[ClExpand];
ClearAll[ClExpand];
SetAttributes[ClExpand,Listable];
ClExpand[x_]:= x/;FreeQ[x, _⊙_];
ClExpand[x_ + y_]:= Simplify[ClExpand[x] + ClExpand[y]];
ClExpand[x_?NumericQ  arg1_]:=x ClExpand[arg1];
ClExpand[arg1_]:= Module[{arg2},
arg2 = Distribute[ExpandAll[arg1],Plus, CircleDot];
If[!FreeQ[arg2,_⊙_] && arg2  =!=  arg1,
arg2 =  ClExpand[arg2];
];
Return[arg2];
];
Protect[ClExpand];
];



(*New!  Contraction *)
LeftContract[x_+y_,z_]:=LeftContract[x,z]+LeftContract[y,z];
LeftContract[x_,y_+z_]:=LeftContract[x,y]+LeftContract[x,z];

LeftContract[c_. (Subscript[e, a_]),x_.(Subscript[e, b_])]:=
If[Length[Intersection[b,a]]==Length[a],ClExpand[c x GradeKPart[ClExpand[(Subscript[e, a])⊙(Subscript[e, b])],(Grade[(Subscript[e, b])]-Grade[(Subscript[e, a])])]],0];

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

RightContract[c_. (Subscript[e, a_]),x_.(Subscript[e, b_])]:=If[Length[Intersection[b,a]]==Length[b],ClExpand[c x GradeKPart[ClExpand[(Subscript[e, a])⊙(Subscript[e, b])],(Grade[(Subscript[e, a])]-Grade[(Subscript[e, b])])]],0];

(* Vector representation of Clifford element *)
VecRep[x_+y_]:=VecRep[x]+VecRep[y];
VecRep[x_.(Subscript[e, a_])]:=(x Table[If[j-1==Total[(2^(a-1))],1,0],{j,1,(2^MaxIndex)}])/;FreeQ[x,_.(Subscript[e, _])];
VecRep[x_]:=(x Table[If[j==1,1,0],{j,1,(2^MaxIndex)}])/;FreeQ[x,_.(Subscript[e, _])]

(* Compute powers of Clifford elements *)
CliffordPower[x_,n_Integer]:=Module[{y},y=ClExpand[x];Switch[EvenQ[n],
True,If[n==0,Return[1],Return[CliffordPower[ClExpand[y⊙y],(n/2)]]],
False,If[n==1,Return[y],Return[ClExpand[y⊙CliffordPower[ClExpand[y⊙y],((n-1)/2)]]]]];];

(* Compute wedge powers of Clifford elements *)
WedgePower[x_,n_Integer]:=Module[{y},y=WExpand[x];Switch[EvenQ[n],
True,If[n==0,Return[1],Return[WedgePower[y⋀y,(n/2)]]],
False,If[n==1,Return[y],Return[WExpand[y⋀WedgePower[y⋀y,((n-1)/2)]]]]];];

(* Procedure to multiply matrices with Clifford entries *)
CliffordMatrixProduct[A_,B_]:=If[Dimensions[A][[2]]≠Dimensions[B][[1]],Abort[],
Table[Total[ClExpand[CircleDot[A[[i]],Transpose[B][[j]]]]],{i,1,Length[A]},{j,1,Dimensions[B][[2]]}]];

(* Procedure to multiply matrices using wedge product *)
WedgeMatrixProduct[A_,B_]:=If[Dimensions[A][[2]]≠Dimensions[B][[1]],Abort[],
Table[Total[WExpand[Wedge[A[[i]],Transpose[B][[j]]]]],{i,1,Length[A]},{j,1,Dimensions[B][[2]]}]];

(* Returns the previously set signature *)
GetSignature:=(Subscript[Cl, pPart,qPart,rPart]);

(* Returns the p-norm in the underlying (2^n) dimensional linear space *)
CliffordNorm[u_,p_]:=Module[{x,V},V=Join[Coefficient[u,Variables[u]],{u/.{(Subscript[e, _])->0}}];Return[Norm[V,p]];];

(* Returns sum of scalar coefficients from canonical expansion *)
ScalarSum[u_]:=Expand[u]/.{(Subscript[e, _])->1};

(* A procedure for computing powers of matrices. In this method, A^m is computed by recursive squaring (((((((A^2)))^2)...))^2) *)
CliffordMatrixPower[A_,m_]:=Module[{y,f},y=ClExpand[A];Switch[EvenQ[m],
True,If[m==0,Return[IdentityMatrix[Length[y]]],Return[ClExpand[CliffordMatrixPower[ClExpand[CliffordMatrixProduct[y,y]],(m/2)]]]],
False,If[m==1,Return[y],Return[ClExpand[CliffordMatrixProduct[CliffordMatrixPower[ClExpand[CliffordMatrixProduct[y,y]],((m-1)/2)],y]]]]]];

(* A procedure for computing wedge powers of matrices. In this method, (A^(^m)) is computed by recursive squaring (((((A^(^2))))^(^2))...) *)
WedgeMatrixPower[A_,m_]:=Module[{y,f},y=ClExpand[A];Switch[EvenQ[m],
True,If[m==0,Return[IdentityMatrix[Length[y]]],Return[WExpand[WedgeMatrixPower[WExpand[WedgeMatrixProduct[y,y]],(m/2)]]]],
False,If[m==1,Return[y],Return[WExpand[WedgeMatrixProduct[WedgeMatrixPower[WExpand[WedgeMatrixProduct[y,y]],((m-1)/2)],y]]]]]];

(* Return max grade *)
Grade[x_+y_]:=Max[Grade[x],Grade[y]];
Grade[_. (Subscript[e, a_])]:=Length[a];

(* Returns the grade k part of an expression *)
GradeKPart[x_,k_Integer]:=
DeleteCases[DeleteCases[If[k==0,ClExpand[x]/.{(Subscript[e, _])->0},ClExpand[x-(x/.{(Subscript[e, Table[_,{k}]])->0})]],0. (Subscript[e, _])],0.`];

(* The involutions are self-explanatory *)
GradeInvolution[x_]:=Module[{k},Sum[(-1)^k GradeKPart[x,k],{k,0,MaxIndex}]];

Reversion[x_]:=Module[{k},Sum[(-1)^(k*(k-1)/2)* GradeKPart[x,k],{k,0,MaxIndex}]];

CliffordConjugate[x_]:=Reversion[GradeInvolution[x]];

(* Common notation for involutions is also supported. *)

Unprotect[OverTilde];
OverTilde[x_]:=Reversion[x];
Protect[OverTilde];

Unprotect[OverHat];
OverHat[x_]:=GradeInvolution[x];
Protect[OverHat];

Unprotect[OverBar];
OverBar[x_]:=CliffordConjugate[x];
Protect[OverBar];

(* Returns the collection of all basis multivectors associated with a signature *)
ClBasis:=Module[{T,j},(T=Table[j,{j,1,MaxIndex}];Table[(Subscript[e, DeleteCases[Reverse[IntegerDigits[i-1,2,MaxIndex]] T,0]]),{i,1,(2^MaxIndex)}]/.{(Subscript[e, {}])->1})];

(* Returns the basis k-vectors associated with a signature *)
ClKVectors[k_]:=DeleteCases[GradeKPart[ClBasis,k],0];

EndPackage[];