(* ::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["CliffMath2014`"];

(* Set aside notation for canonical generators *)

Unprotect[e];
e=Symbol["e"];
Protect[e];
Unprotect["Cl"];
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 *)

Unprotect[Mu];
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 *)
Unprotect[SetSymmetricDifference];
SetSymmetricDifference[A_,B_]:=Complement[Union[A,B],Intersection[A,B]];
Protect[SetSymmetricDifference];


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

Unprotect[SetSignature];
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];

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

(* Return the (maximal) grade of an element *)
Unprotect[Grade];

Grade[B_]:=Module[{m},
If[FreeQ[B,Subscript[e, _]],0,
m=MaxIndex;
While[(m>0\[And]Norm[Coefficient[B,ClKVectors[m]]]==0),
m--];
m]
];


SetAttributes[Grade,Listable];
Protect[Grade];
]

Unprotect[SetAlternatingSignature];
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 *)
Unprotect[{LeftContract, RightContract}];
LeftContract[x_+y_,u_]:=Block[{$RecursionLimit=\[Infinity]},LeftContract[x,u]+LeftContract[y,u]];
LeftContract[x_,y_+u_]:=Block[{$RecursionLimit=\[Infinity]},LeftContract[x,y]+LeftContract[x,u]];

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

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, _]]);
Protect[{LeftContract, RightContract}];


(* Vector representation of Clifford element *)
Unprotect[VecRep];
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, _]];
Protect[VecRep];

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

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


(* Compute powers of Clifford elements *)
Unprotect[CliffordPower];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]]]]];];
Protect[CliffordPower];

(* Compute wedge powers of Clifford elements *)
Unprotect[WedgePower];
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]]]]];];
Protect[WedgePower];

(* Procedure to multiply matrices with Clifford entries *)
Unprotect[CliffordMatrixProduct];
CliffordMatrixProduct[A_,B_]:=Inner[CircleDot, A, B, Plus];
Protect[CliffordMatrixProduct];

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

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

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

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

(* A procedure for computing powers of matrices. In this method, A^m is computed by recursive squaring ((A^2)^2...)^2 *)
Unprotect[CliffordMatrixPower];
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]]]]]];
Protect[CliffordMatrixPower];

(* A procedure for computing wedge powers of matrices. In this method, A^(^m) is computed by recursive squaring ((A^(^2))^(^2)...) *)
Unprotect[WedgeMatrixPower];
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]]]]]];
Protect[WedgeMatrixPower];


(* Returns the grade k part of an expression *)
Unprotect[GradeKPart];
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.`];
Protect[GradeKPart];

(* The involutions are self-explanatory *)
Unprotect[{GradeInvolution, Reversion, CliffordConjugate}];
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]];
Protect[{GradeInvolution, Reversion, CliffordConjugate}];

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


Unprotect[ClBasis];
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}];
Protect[ClBasis];

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

Unprotect[ClKVectors];
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}];
Protect[ClKVectors];

(* Clifford List Product:  Computes ordered Clifford product of a list of Clifford elements *)

Unprotect[ClListProd];
ClListProd[L_,chop_:10^-12]:=
If[Length[L]==1,L[[1]],If[Length[L]==2,Chop[ClExpand[CircleDot[L[[1]],L[[2]]]],chop],Chop[Expand[ClExpand[CircleDot[ClListProd[L[[1;;Floor[Length[L]/2]]]],ClListProd[L[[Floor[Length[L]/2]+1;;Length[L]]]]]]],chop]]];
Protect[ClListProd];

(* Wedge List Product: Computes ordered wedge product of a list of Clifford elements  *)
Unprotect[WListProd];
WListProd[L_]:=
If[Length[L]==1,L[[1]],If[Length[L]==2,Chop[WExpand[Wedge[L[[1]],L[[2]]]]],Chop[Expand[WExpand[Wedge[WListProd[L[[1;;Floor[Length[L]/2]]]],WListProd[L[[Floor[Length[L]/2]+1;;Length[L]]]]]]]]]];
Protect[WListProd];

(* Computes Overscript[B, ~]/(BOverscript[B, ~]), the inverse of B when B is invertible. *) 
Unprotect[ClPseudoInverse];
ClPseudoInverse[B_]:=Chop[Expand[Reversion[B]/Chop[ClExpand[CircleDot[B,Reversion[B]]]]]];
Protect[ClPseudoInverse];

(* Signed 2-seminorm of invertible B *)
Unprotect[ClSignedSeminorm];
ClSignedSeminorm[B_]:=Module[{Bnt},Bnt=Chop[N[ClExpand[CircleDot[B,Reversion[B]]]]];
Sign[Bnt]Sqrt[Abs[Bnt]]];
Protect[ClSignedSeminorm];

(* Extract all indices appearing in sum of multi-indexed elements. *)

Unprotect[GetIndex];
GetIndex[x_]:={}/;FreeQ[x,Subscript[e, _]];
GetIndex[_ Subscript[e, a_]]:=a;
GetIndex[Subscript[e, a_]]:=a;
GetIndex[x_+y_]:=Block[{$RecursionLimit=\[Infinity]},Union[GetIndex[x],GetIndex[y]]];
Protect[GetIndex];


EndPackage[];


