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




BeginPackage["CliffMath11`"];

(* Set aside notation for canonical generators *)

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


(* Mu_i[J]:=|{j\[Epsilon]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]]/. {
\!\(\*SubscriptBox[\(e\), \({}\)]\)->1};

CircleDot[x_, y_]:=Expand[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_]:= Block[{$RecursionLimit=\[Infinity]},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[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\[Intersection]b]!=0,0,Expand[AntisymmetricProductSignatureMap[a,b]x y ] Subscript[e, (a\[Union]b)]];
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_]:= Block[{$RecursionLimit=\[Infinity]},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[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]]/. {
\!\(\*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_]:=  Block[{$RecursionLimit=\[Infinity]},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[Collect[arg2,Subscript[e, _]]];
];
Protect[ClExpand];
];
Protect[SetAlternatingSignature];


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

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

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

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



(* Vector representation of Clifford element *)
VecRep[x_+y_]:=Block[{$RecursionLimit=\[Infinity]},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\[CircleDot]y),n/2]]],
False,If[n==1,Return[y],Return[ClExpand[y\[CircleDot]CliffordPower[(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_]:=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=\[Infinity]},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}]/.{
\!\(\*SubscriptBox[\(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]}]/.{
\!\(\*SubscriptBox[\(e\), \({}\)]\)->1}];

EndPackage[];