(* ::Package:: *)

(************************************************************************)
(* 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.                                                         *)
(************************************************************************)



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

(* Use Wedge for exterior product *)
(* x\[Wedge]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\[Intersection]b,rSet]]>=1,0,(-1)^(Length[Intersection[a\[Intersection]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\[Union]b)\[Intersection](Complement[Union[a,b],Intersection[a,b]])]/. {
\!\(\*SubscriptBox[\(e\), \({}\)]\)->1}];

CircleDot[x_, y_]:=x y /;(FreeQ[x,Subscript[e, _]]\[Or]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,_\[CircleDot]_] && 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\[Intersection]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\[Union]b)]/. {
\!\(\*SubscriptBox[\(e\), \({}\)]\)->1}];
Wedge[x_, y_]:=x y/;(FreeQ[x,Subscript[e, _]]\[Or]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,_\[Wedge]_] && 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,_\[CircleDot]_] && 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\[Intersection]b\[Intersection]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\[Union]b)\[Intersection](Complement[Union[a,b],Intersection[a,b]])]/. {
\!\(\*SubscriptBox[\(e\), \({}\)]\)->1};

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

CircleDot=Symbol["CircleDot"];
Protect[CircleDot];
Unprotect[ClExpand];
ClearAll[ClExpand];
SetAttributes[ClExpand,Listable];
ClExpand[x_]:= x/;FreeQ[x, _\[CircleDot]_];
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,_\[CircleDot]_] && 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]\[CircleDot]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]\[CircleDot]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\[CircleDot]y],n/2]]],
False,If[n==1,Return[y],Return[ClExpand[y\[CircleDot]CliffordPower[ClExpand[y\[CircleDot]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\[Wedge]y,n/2]]],
False,If[n==1,Return[y],Return[WExpand[y\[Wedge]WedgePower[y\[Wedge]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}]/. {
\!\(\*SubscriptBox[\(e\), \({}\)]\)->1})];

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



EndPackage[];