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



(* CliffMath08.nb  Procedures for Clifford algbebra computations in Mathematica *)
(* June 2008 G. Stacey Staples *)

(* Thanks to Marco Budinich for helpful comments and suggestions. *)
 
(* 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["CliffMath08`"];


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

SetSignature[p_,q_]:=Module[{i,j},
Unprotect[MaxIndex];
MaxIndex=p+q;
Protect[MaxIndex];
Unprotect[PseudoScalar];
PseudoScalar:=Subscript[e, Table[i,{i,1,MaxIndex}]];
Protect[PseudoScalar];
Unprotect[pPart, pSet,qPart];   
pPart=p; qPart=q;
pSet = Table[i,{i,1,pPart}];
Protect[pPart, pSet,qPart];
Unprotect[CircleDot];
ClearAll[CircleDot];   
SetAttributes[CircleDot,{Flat, OneIdentity, Listable}];

CircleDot[x_?NumericQ ,y_?NumericQ]:=x y ;                                         (* One 'no arguments' case *)
CircleDot[arg1_,y_?NumericQ]:=y  arg1;                                                  (* Four 'one argument' cases *)
CircleDot[x_?NumericQ, arg2_]:=x arg2;
CircleDot[x_?NumericQ  arg1_,y_?NumericQ]:=x y  arg1;
CircleDot[x_?NumericQ,y_?NumericQ  arg2_]:=x y  arg2;
CircleDot[x_?NumericQ  arg1_,y_?NumericQ  arg2_]:= x y  CircleDot[arg1,arg2];    (* Three 'two arguments' cases *)
CircleDot[x_?NumericQ  arg1_,  arg2_]:=x CircleDot[arg1,arg2];
CircleDot[arg1_,y_?NumericQ  arg2_]:= y  CircleDot[arg1,arg2];CircleDot[Subscript[e, a_],Subscript[e, b_]]:=(-1)^(Length[Complement[a\[Intersection]b,pSet]]+Ceiling[FractionalPart[1/2 Sum[Sum[If[a[[i]]>b[[j]],1,0],{i,1,Length[a]}],{j,1,Length[b]}]]])  Subscript[e, (a\[Union]b)\[Intersection](Complement[Union[a,b],Intersection[a,b]])]/. {
\!\(\*SubscriptBox["e", 
RowBox[{"{", "}"}]]\)->1};
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];

(* Define wedge product *)

Unprotect[Wedge];
ClearAll[Wedge];
SetAttributes[Wedge, {Flat, OneIdentity, Listable}];
Wedge[x_?NumericQ ,y_?NumericQ]:=x y ;                                  
Wedge[arg1_,y_?NumericQ]:=y  arg1;                                            
Wedge[x_?NumericQ, arg2_]:=x arg2;
Wedge[x_?NumericQ  arg1_,y_?NumericQ]:=x y  arg1;
Wedge[x_?NumericQ,y_?NumericQ  arg2_]:=x y  arg2;
Wedge[x_?NumericQ  arg1_,y_?NumericQ  arg2_]:= x y Wedge[arg1,arg2];    
Wedge[x_?NumericQ  arg1_,  arg2_]:=x Wedge[arg1,arg2];
Wedge[arg1_,y_?NumericQ  arg2_]:= y  Wedge[arg1,arg2];
Wedge[Subscript[e, a_],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]}]]]  Subscript[e, (a\[Union]b)]/. {
\!\(\*SubscriptBox["e", 
RowBox[{"{", "}"}]]\)->1}];
Wedge=Symbol["Wedge"];
Protect[Wedge];
Unprotect[WExpand];
ClearAll[WExpand];
SetAttributes[WExpand,Listable];
WExpand[x_]:= x/;FreeQ[x, _\[Wedge]_];
WExpand[x_ + y_]:= Simplify[WExpand[x] + WExpand[y]];
WExpand[x_?NumericQ  arg1_]:=x WExpand[arg1];
WExpand[arg1_]:= Module[{arg2},
arg2 = Distribute[ExpandAll[arg1],Plus, Wedge];
If[!FreeQ[arg2,_\[Wedge]_] && arg2  =!=  arg1,
arg2 =  WExpand[arg2];
];
Return[arg2];
];
Protect[WExpand];];

(* 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_?NumericQ ,y_?NumericQ]:=x y ;                                         
CircleDot[arg1_,y_?NumericQ]:=y  arg1;                                                  
CircleDot[x_?NumericQ, arg2_]:=x arg2;
CircleDot[x_?NumericQ  arg1_,y_?NumericQ]:=x y  arg1;
CircleDot[x_?NumericQ,y_?NumericQ  arg2_]:=x y  arg2;
CircleDot[x_?NumericQ  arg1_,y_?NumericQ  arg2_]:= x y  CircleDot[arg1,arg2]; 
CircleDot[x_?NumericQ  arg1_,  arg2_]:=x CircleDot[arg1,arg2];
CircleDot[arg1_,y_?NumericQ  arg2_]:= y  CircleDot[arg1,arg2];CircleDot[Subscript[e, a_],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]}]]])  Subscript[e, (a\[Union]b)\[Intersection](Complement[Union[a,b],Intersection[a,b]])]/. {
\!\(\*SubscriptBox["e", 
RowBox[{"{", "}"}]]\)->1};
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];
];

(* 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[WExpand[y\[Wedge]y],n/2]]],
False,If[n==1,Return[y],Return[WExpand[y\[Wedge]WedgePower[WExpand[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];

(* 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]]]]]];

(* 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", 
RowBox[{"{", "}"}]]\)->1})];

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



EndPackage[];