CliffMath11`

Clifford algebra computations in Mathematica

G. Stacey Staples

August 2011

BeginPackage["CliffMath11`"];

(* Set aside notation for canonical generators *)

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


(* Mu_i[J]:=|{jϵJ: j>i}|  
J: subset of natural numbers, i: a nonnegative integer

Count elements of set J that are greater than nonnegative integer i *)

Mu[i_,J_]:=Module[{j},Sum[If[j>i,1,0],{j,J}]];
Protect[Mu];

(* ProductSignatureMap[I,J]  Multi-index product signature map 
I,J: Subsets of natural numbers *)
Unprotect[ProductSignatureMap];
ProductSignatureMap[I_,J_]:=Module[{j},If[Mu[pPart+qPart,Intersection[I,J]]≠0,0,(-1)^(Sum[Mu[j,I],{j,J}]+Mu[pPart,Intersection[I,J]])]];
Protect[ProductSignatureMap];

(* AltProductSignatureMap[I,J]  Alternating multi-index product signature map 
even indices square to 1 
odd indices square to -1
I,J: Subsets of natural numbers *)

Unprotect[AltProductSignatureMap];
AltProductSignatureMap[I_,J_]:=Module[{j},(-1)^(Sum[Mu[j,I],{j,J}]+Length[Select[Intersection[I,J],OddQ]])];
Protect[AltProductSignatureMap];


(* AntisymmetricSigMap[I,J]  Multi-index product signature map 
I,J: Subsets of natural numbers *)

Unprotect[AntisymmetricProductSignatureMap];
AntisymmetricProductSignatureMap[I_,J_]:=Module[{j},If[Mu[0,Intersection[I,J]]≠0,0,(-1)^(Sum[Mu[j,I],{j,J}])]];
Protect[AntisymmetricProductSignatureMap];


(* Compute the set-symmetric difference of sets A and B *)

SetSymmetricDifference[A_,B_]:=Complement[Union[A,B],Intersection[A,B]];


(* Initialize Clifford product and wedge product 
Creates symbols pPart, qPart, rPart holding details of signature
MaxIndex indicates number of generators *)


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


Unprotect[CircleDot];(* Clifford product *)
ClearAll[CircleDot];   
SetAttributes[CircleDot,{Flat, OneIdentity, Listable}];

CircleDot[x_. (Subscript[e, a_]),y_. (Subscript[e, b_])]:=Expand[x y ProductSignatureMap[a,b]](Subscript[e, SetSymmetricDifference[a,b]])/.{(Subscript[e, {}])->1};

CircleDot[x_, y_]:=Expand[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_]:= Block[{$RecursionLimit=∞},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[Collect[arg2,(Subscript[e, _])]];
];
Protect[ClExpand];
Protect[SetSignature];


(* Define exterior (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,Expand[AntisymmetricProductSignatureMap[a,b]x y ] (Subscript[e, (a⋃b)])];
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_]:= Block[{$RecursionLimit=∞},Simplify[WExpand[x] + WExpand[y]]];
WExpand[arg1_]:= Module[{arg2},
arg2 = Distribute[ExpandAll[arg1],Plus, Wedge];
If[!FreeQ[arg2,_⋀_] && arg2  !=  arg1,
arg2 =  WExpand[arg2];
];
Return[Collect[arg2,(Subscript[e, _])]];
];

Protect[WExpand];

];



SetAlternatingSignature[n_]:=Module[{},
Unprotect[MaxIndex];
MaxIndex=2n;
Protect[MaxIndex];
Unprotect[PseudoScalar];
PseudoScalar:=(Subscript[e, Range[MaxIndex]]);
Protect[PseudoScalar];
Unprotect[pPart, pSet,qPart, rPart];   
pPart="Alt"; qPart=2n;rPart=0;
Protect[pPart,qPart,rPart];
Unprotect[CircleDot];
ClearAll[CircleDot];  
SetAttributes[CircleDot,{Flat, OneIdentity, Listable}];
CircleDot[x_. (Subscript[e, a_]),y_. (Subscript[e, b_])]:=Expand[AltProductSignatureMap[a,b]x y ](Subscript[e, SetSymmetricDifference[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_]:=  Block[{$RecursionLimit=∞},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[Collect[arg2,(Subscript[e, _])]];
];
Protect[ClExpand];
];
Protect[SetAlternatingSignature];


(* Contraction operators *)
LeftContract[x_+y_,z_]:=Block[{$RecursionLimit=∞},LeftContract[x,z]+LeftContract[y,z]];
LeftContract[x_,y_+z_]:=Block[{$RecursionLimit=∞},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];

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

RightContract[x_+y_,z_]:=Block[{$RecursionLimit=∞},RightContract[x,z]+RightContract[y,z]];
RightContract[x_,y_+z_]:=Block[{$RecursionLimit=∞},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];

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



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

(* Convert vector in (R^(2^n)) to Clifford multivector *)

VectorToMultivector[x_]:=Module[{i},Sum[x[[i]] ClBasis[[i]],{i,(2^MaxIndex)}]];
Protect[VectorToMultivector];


(* 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[(y⊙y),(n/2)]]],
False,If[n==1,Return[y],Return[ClExpand[y⊙CliffordPower[(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_]:=Inner[CircleDot, A, B, Plus];

(* Procedure to multiply matrices using wedge product *)
WedgeMatrixProduct[A_,B_]:=Inner[Wedge, A, B, Plus];

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

(* 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},y=ClExpand[A];Switch[EvenQ[m],
True,If[m==0,Return[IdentityMatrix[Length[y]]],Return[CliffordMatrixPower[CliffordMatrixProduct[y,y],(m/2)]]],
False,If[m==1,Return[y],Return[ClExpand[CliffordMatrixProduct[CliffordMatrixPower[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},y=WExpand[A];Switch[EvenQ[m],
True,If[m==0,Return[IdentityMatrix[Length[y]]],Return[WedgeMatrixPower[WedgeMatrixProduct[y,y],(m/2)]]],
False,If[m==1,Return[y],Return[WExpand[WedgeMatrixProduct[WedgeMatrixPower[WedgeMatrixProduct[y,y],((m-1)/2)],y]]]]]];

(* Return max grade of arbitrary element, or grade of Grassmann monomial *)
Grade[x_+y_]:=Block[{$RecursionLimit=∞},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 blades (Grassmann monomials) associated with a signature 
Ordering is InvLex *)


ClBasis:=Module[{i,table,pwrset},
pwrset=Subsets[Range[MaxIndex]][[2;;]];
table=Prepend[Sort[pwrset,Total[(2^(#1))]<Total[(2^(#2))]&],{}];
Table[(Subscript[e, table[[i]]]),{i,2^MaxIndex}]/.{(Subscript[e, {}])->1}];

(* Returns the basis k-vectors associated with a signature
k is assumed to be a positive integer less than or equal to MaxIndex *)

ClKVectors[k_]:=Module[{i},
pwrset=Subsets[Range[MaxIndex],{k}];
table=Sort[pwrset,Total[(2^(#1))]<Total[(2^(#2))]&];
Table[(Subscript[e, table[[i]]]),{i,Binomial[MaxIndex,k]}]/.{(Subscript[e, {}])->1}];

EndPackage[];