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



(* ::Code::Initialization:: *)
BeginPackage["CliffMath2018`"];

(*

===============
Disclaimer: These software packages are provided "AS IS" without warranty of any kind.   
Should any program prove defective, the end-user assumes the cost of any necessary service, repair, or correction. 
===============

*)

(* Set aside notation for canonical generators *)

Unprotect[e];
e=Symbol["e"];
Protect[e];
Unprotect[\[Zeta]];
\[Zeta]=Symbol["\[Zeta]"];
Protect[\[Zeta]];
Unprotect[\[Stigma]];
\[Stigma]=Symbol["\[Stigma]"];
Protect[\[Stigma]];
Unprotect[\[CurlyEpsilon]];
\[CurlyEpsilon]=Symbol["\[CurlyEpsilon]"];
Protect[\[CurlyEpsilon]];
Unprotect["Cl"];
Cl=Symbol["Cl"];
Protect["Cl"];

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

Unprotect[setSignature];
setSignature::usage="setSignature[p,q,r] - Initializes Clifford multiplication in \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\), as well as `commutative Clifford' multiplication in \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(p, q, r\)], \(sym\)]\).  For n=p+q+r, initializes exterior (Grassmann) multiplication in \!\(\*SuperscriptBox[\(\[CapitalLambda]R\), \(n\)]\), `idempotent-Clifford' multiplication in \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(idem\)]\)  and `zeon multiplication' in \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\)";
setSignature[P_,Q_:0,R_:0]:=Module[{},
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];

(* Default recursion limit to be used by particular CliffMath procedures. Can be "tuned" for performance. *)
Unprotect[recLimit];
recLimit=1024;
Protect[recLimit];

ClearAll[CircleDot,clGrade,cliffordConjugate,clIdemNorm,cliffordDecomp,cliffordDecomposableQ,cliffordGeneratorRep];
ClearAll[cliffordInnerProduct,cliffordInverse,cliffordInvertibleQ,cliffordListProduct,cliffordMatrixPower,cliffordMatrixProduct];
ClearAll[cliffordNorm,cliffordPower,conjugatePseudoInverse,cliffordRep,clRecursionLimit,cliffordReversion,clMu];
ClearAll[clSymNorm, etaMatrix,fastBladeFactor,generalCliffordElement,generalIdemElement,generalSymElement];
ClearAll[generalZeonElement,generalZeonAntiDerivative,generalZeonDerivative,getIndex,getSignature,clGrade];
ClearAll[gradeBasis,gradeInvolution,gradeKPart,idemBasis,idemInnerProduct,leftContract,matrixKroneckerPower];
ClearAll[nilpotentZeonQ,OverBar,OverHat,OverTilde,parityCheck,repMatrixProduct,rightContract,scalarSum];
ClearAll[setSymmetricDifference,signedSeminormSquare,symBasis,symInnerProduct,versorFactor,wedge,wedgeListProduct];
ClearAll[wedgePower,zDu,zRe,zeonAntiDerivative,zeonAbs,zeonArcCos,zeonArcSin,zeonArcTan,zeonBasis,zeonCos];
ClearAll[zeonCosh,zeonDerivative,zeonDiscriminant,zeonDivisionAlgorithm,zeonDual,zeonElementaryFactorization];
ClearAll[zeonExp,zeonExtension,zeonExtensionEval,zeonFrobeniusNorm,zeonGeneratorRep,zeonGeometricSeries];
ClearAll[zeonInnerProduct,zeonInverse,zeonKthRoot,zeonLeftContraction,zeonLog,zeonMaxGenDecomp,zeonMaxGenerator];
ClearAll[zeonMinGrade,zeonNilpotency,zeonNorm,zeonPower,zeonPowerZeon,zeonQuadraticEval,zeonQuadraticFormula];
ClearAll[zeonRep,zeonRightContraction,zeonSin,zeonSinh,zeonSpectralSeminorm,zeonSqRoot,zeonTan,zeonTanh];

ClearAll[clExpand,productSignatureMap,symProductSignatureMap,zeonProductSignatureMap];

(* mu_i[J]:=|{j\[Element]J: j>i}|  
J: subset of natural numbers, i: a nonnegative integer
Count elements of set J that are greater than nonnegative integer i *)

clMu::usage="clMu[i,J] - Counts number of elements in multi-index J that are greater than i \[Element] \!\(\*SubscriptBox[\(\[DoubleStruckCapitalN]\), \(0\)]\).";
clMu[i_,J_]:=clMu[i,J]=Module[{j},Sum[If[j>i,1,0],{j,J}]];

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

antisymmetricProductSignatureMap[I_,J_]:=antisymmetricProductSignatureMap[I,J]=
Module[{j},If[clMu[0,Intersection[I,J]]!=0,0,(-1)^(Sum[clMu[j,I],{j,J}])]];

(* Compute the set-symmetric difference of sets A and B *)
setSymmetricDifference::usage="setSymmetricDifference[A,B] - returns (A\[Union]B))\[Intersection](A\[Intersection]B)', the symmetric difference of sets A and B.";
setSymmetricDifference[A_,B_]:=setSymmetricDifference[A,B]=
Complement[Union[A,B],Intersection[A,B]];

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

symProductSignatureMap[I_,J_]:=symProductSignatureMap[I,J]=
If[clMu[pPart+qPart,Intersection[I,J]]!=0,0,(-1)^(clMu[pPart,Intersection[I,J]])];

(*Protect[productSignatureMap];*)

zeonProductSignatureMap[a_,b_]:=zeonProductSignatureMap[a,b]=If[Length[a\[Intersection]b]>0,0,1];

(*Unprotect[CircleDot];(* Clifford product *)*)
SetAttributes[CircleDot,{Flat, OneIdentity, Listable}];
CircleDot::usage="CircleDot[u,v] - returns product u \[CircleDot] v in Clifford algebra or subalgebra determined by canonical basis elements.";
CircleDot[x_, y_]:=Expand[x y] /;((FreeQ[x,Subscript[\[Zeta], _]]\[And]FreeQ[x,Subscript[e, _]]\[And]FreeQ[x,Subscript[\[Stigma], _]]\[And]FreeQ[x,Subscript[\[CurlyEpsilon], _]])\[Or](FreeQ[y,Subscript[\[Zeta], _]]\[And]FreeQ[y,Subscript[e, _]]\[And]FreeQ[y,Subscript[\[Stigma], _]]\[And]FreeQ[y,Subscript[\[CurlyEpsilon], _]]));
CircleDot[x_, y_]:=0/;{x==0\[Or]y==0};

CircleDot[x_. Subscript[e, a_],y_. Subscript[e, b_]]:=Expand[CircleDot[x,y] productSignatureMap[a,b]]Subscript[e, setSymmetricDifference[a,b]]/. {\!\(\*SubscriptBox[\(e\), \({}\)]\)->1};
CircleDot[x_. Subscript[\[Zeta], a_],y_. Subscript[\[Zeta], b_]]:=zeonProductSignatureMap[a,b]Expand[CircleDot[x,y] Subscript[\[Zeta], a\[Union]b]];
CircleDot[x_. Subscript[\[Stigma], a_],y_. Subscript[\[Stigma], b_]]:=Expand[CircleDot[x,y] symProductSignatureMap[a,b]]Subscript[\[Stigma], setSymmetricDifference[a,b]]/. {\!\(\*SubscriptBox[\(\[Stigma]\), \({}\)]\)->1};
CircleDot[x_. Subscript[\[CurlyEpsilon], a_],y_. Subscript[\[CurlyEpsilon], b_]]:=Expand[CircleDot[x,y]]Subscript[\[CurlyEpsilon], a\[Union]b]/. {\!\(\*SubscriptBox[\(\[CurlyEpsilon]\), \({}\)]\)->1};


(*Protect[CircleDot];*)

clRecursionLimit::usage="clRecursionLimit[m_] - sets recursion limit for particular CliffMath procedures.  
Mathematica's default value is 1024.  In some cases, higher limits are desirable, but 1024 is minimum allowed.";
clRecursionLimit[m_]:=Module[{},
Unprotect[recLimit];
Print[recLimit=Max[1024, m]];
Protect[recLimit];];

(*Unprotect[clExpand];*)
SetAttributes[clExpand,Listable];
clExpand[x_ + y_]:= Block[{$RecursionLimit=recLimit},ExpandAll[clExpand[x] + clExpand[y]]];
clExpand[x_?NumericQ  arg1_]:=Expand[x clExpand[arg1]];
clExpand[arg1_]:= Module[{arg2,arg3},
arg2 = Distribute[ExpandAll[arg1],Plus, CircleDot];
If[!FreeQ[arg2,_\[CircleDot]_] && arg2  !=  arg1,
arg2 =  clExpand[arg2];
];
arg3=Distribute[ExpandAll[arg2],Plus,wedge];
If[!FreeQ[arg3,_\[Wedge]_] && arg3  !=  arg2,
arg3=  clExpand[arg3];
];
arg2=arg3;
Return[If[!FreeQ[e][arg2],Collect[arg2,Subscript[e, _]],If[!FreeQ[\[Zeta]][arg2],Collect[arg2,Subscript[\[Zeta], _]],If[!FreeQ[\[Stigma]][arg2],Collect[arg2,Subscript[\[Stigma], _]],If[!FreeQ[\[CurlyEpsilon]][arg2],Collect[arg2,Subscript[\[CurlyEpsilon], _]],arg2]]]
]];
(*Protect[clExpand];*)];
(*Protect[setSignature];*)

(* Define exterior (wedge) product *)
SetAttributes[wedge, {Flat, OneIdentity, Listable}];
wedge::usage="wedge[u,v] - returns exterior (wedge) product u\[Wedge]v.";
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, _]]);

(* Accommodate symbolic wedge product. *)
SetAttributes[Wedge, {Flat, OneIdentity, Listable}];
Wedge[a_,b_]:=wedge[a,b];

(* Return the (maximal) grade of an element u \[Element] Subscript[Cl, p,q,r] or u \[Element] Subscript[Cl, *]^* *)
clGrade::usage="clGrade[u] - Maximal clGrade of nonzero terms in canonical expansion of u.";
clGrade[B_]:=Module[{m,X},
If[(FreeQ[B,e]\[And]FreeQ[B,\[Zeta]]\[And]FreeQ[B,\[Stigma]]\[And]FreeQ[B,\[CurlyEpsilon]]),0,
X=Variables[B]/.{Subscript[_, a_]->a};
Return[Max[Table[Length[X[[i]]],{i,Length[X]}]]]]];
SetAttributes[clGrade,Listable];


(* Clifford algebra contraction operators *)
leftContract::usage="leftContract[u,v] - Clifford left contraction u\[RightFloor]v.";
leftContract[x_+y_,u_]:=Block[{$RecursionLimit=recLimit},Expand[leftContract[x,u]+leftContract[y,u]]];
leftContract[x_,y_+u_]:=Block[{$RecursionLimit=recLimit},Expand[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]],(clGrade[Subscript[e, b]]-clGrade[Subscript[e, a]])]],0];

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

rightContract::usage="rightContract[u,v] - Clifford right contraction u\[LeftFloor]v.";
rightContract[x_+y_,u_]:=Block[{$RecursionLimit=recLimit},rightContract[x,u]+rightContract[y,u]];
rightContract[x_,y_+u_]:=Block[{$RecursionLimit=recLimit},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]],(clGrade[Subscript[e, a]]-clGrade[Subscript[e, b]])]],0];

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


(* Compute powers of Clifford elements *)
cliffordPower::usage="cliffordPower[u,m] - Compute \!\(\*SuperscriptBox[\(u\), \(\(\\\ \)\(\[CircleDot]m\)\)]\) for nonnegative integer m and u \[Element] \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\) or other \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(*\)], \(*\)]\).";
cliffordPower[x_,n_Integer]:=cliffordPower[x,n]=Module[{y},y=clExpand[x];Switch[EvenQ[n],
True,If[n==0,Return[1],Return[clExpand[CircleDot[cliffordPower[y,n/2],cliffordPower[y,n/2]]]]],
False,If[n==1,Return[y],Return[clExpand[y \[CircleDot] clExpand[CircleDot[cliffordPower[y,(n-1)/2],cliffordPower[y,(n-1)/2]]]]]]];];

(* Compute wedge powers of Clifford elements *)
wedgePower::usage="wedgePower[u,m] - Compute \!\(\*SuperscriptBox[\(u\), \(\[Wedge]m\)]\) for nonegative integer m.";
wedgePower[x_,n_Integer]:=Module[{y},y=clExpand[x];Switch[EvenQ[n],
True,If[n==0,Return[1],Return[clExpand[wedge[wedgePower[y,n/2],wedgePower[y,n/2]]]]],
False,If[n==1,Return[y],Return[clExpand[y\[Wedge]clExpand[wedge[wedgePower[y,(n-1)/2],wedgePower[y,(n-1)/2]]]]]]];];

(* Procedure to multiply matrices with Clifford entries *)
cliffordMatrixProduct::usage="cliffordMatrixProduct[A,B] Compute matrix product AB using multiplication in appropriate algebra \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\), \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\), \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(idem\)]\), or \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(p, q, r\)], \(sym\)]\).";
cliffordMatrixProduct[A_,B_]:=Inner[CircleDot, A, B, Plus];

(* Procedure to multiply matrices using wedge (exterior) product *)
wedgeMatrixProduct::usage="wedgeMatrixProduct[A,B] Compute wedge(exterior) matrix product AB.";
wedgeMatrixProduct[A_,B_]:=Inner[wedge, A, B, Plus];

(* Returns the previously set signature *)
getSignature::usage="getSignature - Returns \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\) , where p,q,r represent the currently set signature.";
getSignature:=If[rPart>0,Subscript[Cl, pPart,qPart,rPart],Subscript[Cl, pPart,qPart]];

(* Inner products in Subscript[Cl, p,q,r], Subscript[Cl, n]^nil, Subscript[Cl, p,q,r]^sym, and Subscript[Cl, n]^idem *)

cliffordInnerProduct::usage="cliffordInnerProduct[u,v] - returns <u,v>=\!\(\*UnderscriptBox[\(\[Sum]\), \(I\)]\)\!\(\*SubscriptBox[\(u\), \(I\)]\)\!\(\*SubscriptBox[\(v\), \(I\)]\) inner product of Clifford elements u,v \[Element] \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\).";
cliffordInnerProduct[u_,v_]:=Module[{u0,v0,udotv},
u0=gradeKPart[u,0];
v0=gradeKPart[v,0];
udotv=Prepend[Coefficient[(u-u0),clBasis[[2;;]]],u0].Prepend[Coefficient[(v-v0),clBasis[[2;;]]],v0];
Return[udotv];];

zeonInnerProduct::usage="zeonInnerProduct[u,v] - returns <u,v>=\!\(\*UnderscriptBox[\(\[Sum]\), \(I\)]\)\!\(\*SubscriptBox[\(u\), \(I\)]\)\!\(\*SubscriptBox[\(v\), \(I\)]\) inner product <u,v> of zeon elements u,v \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
SetAttributes[zeonInnerProduct,Listable];
zeonInnerProduct[u_,v_]:=Module[{u0,v0,udotv},
u0=gradeKPart[u,0];
v0=gradeKPart[v,0];
udotv=Prepend[Coefficient[(u-u0),zeonBasis[[2;;]]],u0].Prepend[Coefficient[(v-v0),zeonBasis[[2;;]]],v0];
Return[udotv];];

idemInnerProduct::usage="idemInnerProduct[u,v] - returns <u,v>=\!\(\*UnderscriptBox[\(\[Sum]\), \(I\)]\)\!\(\*SubscriptBox[\(u\), \(I\)]\)\!\(\*SubscriptBox[\(v\), \(I\)]\) inner product <u,v> of elements u,v \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(idem\)]\).";
idemInnerProduct[u_,v_]:=Module[{u0,v0,udotv},
u0=gradeKPart[u,0];
v0=gradeKPart[v,0];
udotv=Prepend[Coefficient[(u-u0),idemBasis[[2;;]]],u0].Prepend[Coefficient[(v-v0),idemBasis[[2;;]]],v0];
Return[udotv];];

symInnerProduct::usage="symInnerProduct[u,v] - returns <u,v>=\!\(\*UnderscriptBox[\(\[Sum]\), \(I\)]\)\!\(\*SubscriptBox[\(u\), \(I\)]\)\!\(\*SubscriptBox[\(v\), \(I\)]\) inner product <u,v> of elements u,v \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(p, q, r\)], \(sym\)]\).";
symInnerProduct[u_,v_]:=Module[{u0,v0,udotv},
u0=gradeKPart[u,0];
v0=gradeKPart[v,0];
udotv=Prepend[Coefficient[(u-u0),symBasis[[2;;]]],u0].Prepend[Coefficient[(v-v0),symBasis[[2;;]]],v0];
Return[udotv];];


(* Returns the p-norm in the underlying 2^n dimensional linear space *)
cliffordNorm::usage="cliffordNorm[u,m] - returns \[LeftDoubleBracketingBar]u\!\(\*SubscriptBox[\(\[RightDoubleBracketingBar]\), \(m\)]\) m-norm of u \[Element] \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\).";
cliffordNorm[u_,p_]:=Norm[Prepend[Coefficient[u,clBasis[[2;;]]],gradeKPart[u,0]],p];

clSymNorm::usage="clSymNorm[u,m] - returns \[LeftDoubleBracketingBar]u\!\(\*SubscriptBox[\(\[RightDoubleBracketingBar]\), \(m\)]\) m-norm of u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(p, q, r\)], \(sym\)]\).";
clSymNorm[u_,p_]:=Norm[Prepend[Coefficient[u,symBasis[[2;;]]],gradeKPart[u,0]],p];
clIdemNorm::usage="clIdemNorm[u,m] - returns \[LeftDoubleBracketingBar]u\!\(\*SubscriptBox[\(\[RightDoubleBracketingBar]\), \(m\)]\) m-norm of u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(idem\)]\).";
clIdemNorm[u_,p_]:=Norm[Prepend[Coefficient[u,idemBasis[[2;;]]],gradeKPart[u,0]],p];

zeonNorm::usage="zeonNorm[u,m] - returns the m-norm \[LeftDoubleBracketingBar]u\!\(\*SubscriptBox[\(\[RightDoubleBracketingBar]\), \(m\)]\) of u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
zeonNorm[u_,p_]:=Norm[Prepend[Coefficient[u,zeonBasis[[2;;]]],zRe[u]],p];

zeonSpectralSeminorm::usage="zeonSpectralSeminorm[u] - returns |\[GothicCapitalR]u| for u\[Element]\!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
zeonSpectralSeminorm[u_]:=Abs[zRe[u]];


(* Returns sum of scalar coefficients from canonical expansion *)
scalarSum::usage="scalarSum[u] returns sum of scalar coefficients in canonical expansion of u \[Element] \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\) or u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(*\)], \(*\)]\).";
scalarSum[u_]:=Expand[u]/. {Subscript[e, _]->1,Subscript[\[Zeta], _]->1,Subscript[\[CurlyEpsilon], _]->1,Subscript[\[Stigma], _]->1};

(* A procedure for computing powers of matrices. In this method, A^m is computed by recursive squaring ((A^2)^2...)^2 *)
cliffordMatrixPower::usage="cliffordMatrixPower[A,m] - Computes \!\(\*SuperscriptBox[
StyleBox[\"A\",\nFontSlant->\"Italic\"], \(\[CircleDot]\*
StyleBox[\"m\",\nFontSlant->\"Italic\"]\)]\) for square matrix \!\(\*
StyleBox[\"A\",\nFontSlant->\"Italic\"]\) having entries from \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\) or \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(*\)], \(*\)]\).";
cliffordMatrixPower[A_,m_]:=cliffordMatrixPower[A,m]=Module[{y},y=clExpand[A];Switch[EvenQ[m],
True,If[m==0,Return[IdentityMatrix[Length[A]]],Return[clExpand[cliffordMatrixProduct[cliffordMatrixPower[y,m/2],cliffordMatrixPower[y,m/2]]]]],
False,If[m==1,Return[y],Return[clExpand[cliffordMatrixProduct[y,clExpand[cliffordMatrixProduct[cliffordMatrixPower[y,(m-1)/2],cliffordMatrixPower[y,(m-1)/2]]]]]]]];];

(* A procedure for computing wedge powers of matrices. In this method, A^(^m) is computed by recursive squaring ((A^(^2))^(^2)...) *)
wedgeMatrixPower::usage="wedgeMatrixPower[A,m] - Computes \!\(\*SuperscriptBox[
StyleBox[\"A\",\nFontSlant->\"Italic\"], \(\[Wedge]m\)]\) for square matrix \!\(\*
StyleBox[\"A\",\nFontSlant->\"Italic\"]\) having entries from \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\).";
wedgeMatrixPower[A_,m_]:=wedgeMatrixPower[A,m]=Module[{y},y=clExpand[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[clExpand[wedgeMatrixProduct[wedgeMatrixPower[wedgeMatrixProduct[y,y],(m-1)/2],y]]]]]];


(* Returns the grade k part of an expression *)
gradeKPart::usage="gradeKPart[u,k] - Returns grade-k part <u\!\(\*SubscriptBox[\(>\), \(k\)]\) of u \[Element] \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\) or \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(*\)], \(*\)]\).";
gradeKPart[x_,k_Integer]:=
DeleteCases[DeleteCases[If[k==0,clExpand[x]/. {Subscript[e, _]->0,Subscript[\[Zeta], _]->0,Subscript[\[Stigma], _]->0,Subscript[\[CurlyEpsilon], _]->0},clExpand[x-(x/. {Subscript[e, Table[_,{k}]]->0,Subscript[\[Zeta], Table[_,{k}]]->0,Subscript[\[Stigma], Table[_,{k}]]->0,Subscript[\[CurlyEpsilon], Table[_,{k}]]->0})]],0._],0.`];

(* The involutions are self-explanatory *)
SetAttributes[gradeInvolution, Listable];
SetAttributes[cliffordReversion, Listable];
SetAttributes[cliffordConjugate, Listable];

gradeInvolution::usage="gradeInvolution[u] - Clifford grade involution of u \[Element] \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\).   u \[ControlKey]& ^";
gradeInvolution[x_]:=Module[{k},Sum[(-1)^k gradeKPart[x,k],{k,0,maxIndex}]];

cliffordReversion::usage="cliffordReversion[u] - Clifford reversion of u \[Element] \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\).   u \[ControlKey]& ~";
cliffordReversion[x_]:=Module[{k},Sum[(-1)^(k*(k-1)/2)* gradeKPart[x,k],{k,0,maxIndex}]];

cliffordConjugate::usage="cliffordConjugate[u] - Clifford conjugate of u \[Element] \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\).   u \[ControlKey]& _";
cliffordConjugate[x_]:=Module[{k},Sum[(-1)^(k*(k+1)/2)* gradeKPart[x,k],{k,0,maxIndex}]];

(* Common notation for involutions is also supported. *)

OverTilde::usage="OverTilde[u] - Reversion of u \[Element] \!\(\*SubscriptBox[\(\[ScriptCapitalC]\[ScriptL]\), 
StyleBox[
RowBox[{\"p\", \",\", \"q\", \",\", \"r\"}],\nFontSlant->\"Italic\"]]\)";
OverTilde[x_]:=cliffordReversion[x];

OverHat::usage="OverHat[u] - grade involution of u \[Element] \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\)";
OverHat[x_]:=gradeInvolution[x];

OverBar::usage="OverBar[u] - Clifford conjugate of u \[Element] \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\)";
OverBar[x_]:=cliffordConjugate[x];

(* Returns the collection of all basis blades (Grassmann monomials) associated with a signature 
Ordering is InvLex *)

clBasis::usage="clBasis - Returns ordered basis blades {\!\(\*SubscriptBox[\(e\), \(\(I\)\(\\\ \)\)]\): I\[Element]\!\(\*SuperscriptBox[\(2\), \([p + q + r]\)]\)} \[SubsetEqual] \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\).";
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}];

zeonBasis::usage="zeonBasis - Returns ordered basis blades {\!\(\*SubscriptBox[\(\[Zeta]\), \(\(I\)\(\\\ \)\)]\): I\[Element]\!\(\*SuperscriptBox[\(2\), \([n]\)]\)} \[SubsetEqual] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
zeonBasis:=clBasis/.{e->\[Zeta]};

symBasis::usage="symBasis - Returns ordered basis blades {\!\(\*SubscriptBox[\(\[Stigma]\), \(\(I\)\(\\\ \)\)]\): I\[Element]\!\(\*SuperscriptBox[\(2\), \([n]\)]\)} \[SubsetEqual] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
symBasis:=clBasis/.{e->\[Stigma]};

idemBasis::usage="idemBasis - Returns ordered basis blades {\!\(\*SubscriptBox[\(\[CurlyEpsilon]\), \(\(I\)\(\\\ \)\)]\): I\[Element]\!\(\*SuperscriptBox[\(2\), \([n]\)]\)} \[SubsetEqual] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
idemBasis:=clBasis/.{e->\[CurlyEpsilon]};

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

gradeBasis::usage="gradeBasis[k,symb] - Returns basis k-blades for positive integer k. 
Parameter symb specifies the algebra: e (default): \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\);  \[Zeta]: \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\); \[CurlyEpsilon]: \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(idem\)]\); \[Stigma]: \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(p, q, r\)], \(sym\)]\).";
gradeBasis[k_,symb_:e]:=Module[{i,table,indices},
indices=Subsets[Range[maxIndex],{k}];
table=Sort[indices,Total[2^#1]<Total[2^#2]&];
Table[Subscript[symb, table[[i]]],{i,Length[indices]}]];

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

cliffordListProduct::usage="cliffordListProduct[{\!\(\*SubscriptBox[\(u\), \(1\)]\), \!\(\*SubscriptBox[\(u\), \(2\)]\), ..., \!\(\*SubscriptBox[\(u\), \(k\)]\)}] returns the ordered product \!\(\*SubscriptBox[\(u\), \(1\)]\) \[CircleDot] \!\(\*SubscriptBox[\(u\), \(2\)]\) \[CircleDot] ... \[CircleDot] \!\(\*SubscriptBox[\(u\), \(k\)]\).";
cliffordListProduct[L_,chop_:10^-12]:=
If[Length[L]==1,L[[1]],If[Length[L]==2,Chop[clExpand[CircleDot[L[[1]],L[[2]]]],chop],Chop[clExpand[CircleDot[cliffordListProduct[L[[1;;Floor[Length[L]/2]]]],cliffordListProduct[L[[Floor[Length[L]/2]+1;;Length[L]]]]]],chop]]];

(* wedge List Product: Computes ordered wedge product of a list of Clifford elements  *)
wedgeListProduct::usage="wedgeListProduct[{\!\(\*SubscriptBox[\(u\), \(1\)]\), \!\(\*SubscriptBox[\(u\), \(2\)]\), ..., \!\(\*SubscriptBox[\(u\), \(k\)]\)}] computes wedge product \!\(\*SubscriptBox[\(u\), \(1\)]\) \[Wedge] \!\(\*SubscriptBox[\(u\), \(2\)]\) \[Wedge] ... \[Wedge] \!\(\*SubscriptBox[\(u\), \(k\)]\)";
wedgeListProduct[L_]:=
If[Length[L]==1,L[[1]],If[Length[L]==2,Chop[clExpand[wedge[L[[1]],L[[2]]]]],Chop[clExpand[wedge[wedgeListProduct[L[[1;;Floor[Length[L]/2]]]],wedgeListProduct[L[[Floor[Length[L]/2]+1;;]]]]]]]];

(* True if all terms of u are of even grade or all terms are of odd grade *)
parityCheck::usage="parityCheck[u] - True if all terms of u are of even grade or all terms are of odd grade.";
parityCheck[u_]:=Module[{v},
v=u/.{\[Zeta]->e,\[CurlyEpsilon]->e,\[Stigma]->e};
Return[Xor[AllTrue[clGrade[DeleteCases[Variables[v],Except[Subscript[e, _]]]],EvenQ],((v/.{Subscript[e, _]->0})==0\[And]NoneTrue[clGrade[DeleteCases[Variables[v],Except[Subscript[e, _]]]],EvenQ])]];];

(* Works in definite signatures.  Decomposition by geometric considerations. *)
versorFactor::usage="versorFactor[u] - Factor versor u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\) (Euclidean Clifford algebra) via geometric considerations.  
Returns vectors {\!\(\*SubscriptBox[\(v\), \(i\)]\):1\[LessEqual]i\[LessEqual]clGrade[u]} and scalar \[Alpha] such that u= \[Alpha]\!\(\*UnderoverscriptBox[\(\[Product]\), \(i = 1\), \(clGrade[u]\)]\)\!\(\*SubscriptBox[\(v\), \(i\)]\).";
versorFactor[B_]:=Module[{ell,j,X,gb,bi,Bprime,btemp,k,x,z,BNorm},
BNorm=N[cliffordNorm[B,2]];
ell=1;
Bprime=Expand[B/BNorm];
While[clGrade[Bprime]>1,
bi=getIndex[Bprime];
x=N[\!\(\*SubscriptBox[\(e\), \({bi[\([\)\(1  Random[Integer, {1, Length[bi]}]\)\(]\)]}\)]\)];

X=Chop[N[clExpand[CircleDot[Bprime,x,cliffordConjugate[Bprime]]]]];
gb=(-1)^clGrade[Bprime];
If[Chop[leftContract[x,X]!=1],
btemp=Chop[Expand[x-X]];
Subscript[b, ell]=Chop[Expand[btemp/cliffordNorm[btemp,2]]];

(* If odd grade, pull out one reflection first *)
If[gb<0,
Bprime=Chop[clExpand[CircleDot[Bprime,Subscript[b, ell]]]];
ell=ell+1;,
(* For even-grade Bprime, compute bivector *)
z=Chop[Expand[N[clExpand[CircleDot[Bprime,Subscript[b, ell],cliffordReversion[Bprime]]]]]];
If[Chop[leftContract[Subscript[b, ell],z]]!=-1,
Subscript[b, ell+1]=N[Chop[Expand[(Subscript[b, ell]+z)/cliffordNorm[(Subscript[b, ell]+z),2]]]];
Bprime=Expand[Chop[clExpand[CircleDot[Bprime,Subscript[b, ell],Subscript[b, ell+1]]]]];
ell=ell+2;,
Bprime=Expand[Chop[clExpand[CircleDot[Bprime,Subscript[b, ell]]]]];
ell=ell+1;];];];];
Return[Append[Table[Subscript[b, j],{j,ell-1}],Expand[BNorm Bprime]]];];

(* A very fast algorithm dealing strictly with blades *)
fastBladeFactor::usage="fastBladeFactor[B] - Returns scalar \[Alpha] and vectors {b_i: 1\[LessEqual]i\[LessEqual]clGrade(B)} such that B=\[Alpha]\!\(\*UnderoverscriptBox[\(\[Product]\), \(i = 1\), \(clGrade \((B)\)\)]\)\!\(\*SubscriptBox[\(b\), \(i\)]\), provided B \[Element] \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\) is a blade.";
fastBladeFactor[B_]:=Module[{ell,\[Alpha],gb,F,Fi,Bprime,btemp,u},
\[Alpha]=B[[1]]/.{Subscript[e, _]->1};
ell=1;
Bprime=Expand[B/\[Alpha]];
F=Expand[B[[1]]/\[Alpha]];
Fi=getIndex[F];
gb=clGrade[F];
While[ell<=gb,
btemp=Subscript[e, Complement[Fi,{Fi[[ell]]}]];
Subscript[b, ell]=gradeKPart[Chop[clExpand[CircleDot[cliffordInverse[btemp],Bprime]]],1];
ell=ell+1];
Return[Reverse[Prepend[Table[Subscript[b, j],{j,ell-1}], \[Alpha]]]];];

(* Clifford algebraic Gram-Schmidt orthogonalization *)
clGramSchmidt::usage="clGramSchmidt[L] - Given set L of linearly independent vectors, returns set B of vectors pairwise-orthogonal 
with respect to the Clifford inner product <u,v> = \!\(\*FractionBox[\(1\), \(2\)]\)(uv+vu).";
clGramSchmidt[L_]:=Module[{B,lp,ell},
B={L[[1]]};
lp=2;
While[lp<=Length[L],B=Append[B,clExpand[L[[lp]]-Sum[(clExpand[L[[lp]]\[CircleDot]B[[ell]]]+clExpand[B[[ell]]\[CircleDot]L[[lp]]])/(2clExpand[B[[ell]]\[CircleDot]B[[ell]]])B[[ell]],{ell,1,lp-1}]]];
lp++;];
Return[B];];


(* Decomposability Test. Necessary in indefinite signatures *)
cliffordDecomposableQ::usage="cliffordDecomposableQ[u] - Decomposability test for Clifford algebra element. Necessary in indefinite signatures";
cliffordDecomposableQ[u_]:=parityCheck[u]\[And]FreeQ[clExpand[CircleDot[u,cliffordReversion[u]]],e]\[And]signedSeminormSquare[gradeKPart[u,clGrade[u]]]!=0;

(* Decompose elements in arbitrary signatures *)
cliffordDecomp::usage="cliffordDecomp[u] - Returns scalar vectors {\!\(\*SubscriptBox[\(v\), \(i\)]\) : 1 \[LessEqual] i \[LessEqual]clGrade[u]} such that u= \!\(\*UnderoverscriptBox[\(\[Product]\), \(i = 1\), \(clGrade[u]\)]\)\!\(\*SubscriptBox[\(v\), \(i\)]\).
Note: Assumes u is decomposable \!\(\*
StyleBox[\"a\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"priori\",\nFontSlant->\"Italic\"]\). \!\(\*
StyleBox[\"If\",\nFontSlant->\"Plain\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Plain\"]\)\!\(\*
StyleBox[\"in\",\nFontSlant->\"Plain\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Plain\"]\)\!\(\*
StyleBox[\"doubt\",\nFontSlant->\"Plain\"]\)\!\(\*
StyleBox[\",\",\nFontSlant->\"Plain\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Plain\"]\)\!\(\*
StyleBox[\"apply\",\nFontSlant->\"Plain\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Plain\"]\)\!\(\*
StyleBox[\"cliffordDecomposableQ\",\nFontSlant->\"Plain\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Plain\"]\)\!\(\*
StyleBox[\"first\",\nFontSlant->\"Plain\"]\)\!\(\*
StyleBox[\".\",\nFontSlant->\"Plain\"]\)";
cliffordDecomp[B_]:=Module[{X,bi,gb,Bprime,btemp,facts},
Bprime=B;
facts={};
While[clGrade[Bprime]>1,
gb=clGrade[Bprime];
bi=getIndex[Bprime];
x=\!\(\*SubscriptBox[\(e\), \({bi[\([\)\(Random[Integer, {1, Length[bi]}]\)\(]\)]}\)]\);X=Chop[clExpand[CircleDot[Bprime,x,conjugatePseudoInverse[Bprime]]]];
btemp=Expand[x- X];
If[Chop[signedSeminormSquare[btemp]]!=0,
Bprime=Chop[clExpand[CircleDot[Bprime, cliffordInverse[btemp]]]];
facts=Append[facts, btemp];
If[clGrade[Bprime]>gb,Bprime=B; facts={};];
];];
Return[Reverse[Append[facts,Bprime]]];];

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

getIndex::usage="getIndex[u] - Extract all indices appearing among multi-indices of u.";
getIndex[x_]:={}/;(FreeQ[x,Subscript[e, _]]\[And]FreeQ[x,Subscript[\[Stigma], _]]\[And]FreeQ[x,Subscript[\[Zeta], _]]\[And]FreeQ[x,Subscript[\[CurlyEpsilon], _]]);
getIndex[_ Subscript[e, a_]]:=a;
getIndex[Subscript[e, a_]]:=a;
getIndex[_ Subscript[\[Zeta], a_]]:=a;
getIndex[Subscript[\[Zeta], a_]]:=a;
getIndex[_ Subscript[\[Stigma], a_]]:=a;
getIndex[Subscript[\[Stigma], a_]]:=a;
getIndex[_ Subscript[\[CurlyEpsilon], a_]]:=a;
getIndex[Subscript[\[CurlyEpsilon], a_]]:=a;
getIndex[x_+y_]:=Block[{$RecursionLimit=recLimit},Union[getIndex[x],getIndex[y]]];

(* Signed seminorm square <Overscript[\!\(\(uu\)
\*SubscriptBox[\(>\), \(0\)]\), ~] *)
signedSeminormSquare::usage="signedSeminormSquare[u] - returns <u\!\(\*OverscriptBox[\(u\), \(~\)]\)\!\(\*SubscriptBox[\(>\), \(0\)]\).";
signedSeminormSquare[u_]:=gradeKPart[clExpand[CircleDot[u,cliffordReversion[u]]],0];

(* Invertibility Test. Necessary in indefinite signatures *)
cliffordInvertibleQ::usage="cliffordInvertibleQ[u] - Invertibility test for Clifford algebra element. Necessary in indefinite signatures";
cliffordInvertibleQ[u_]:=PossibleZeroQ[clExpand[u \[CircleDot] cliffordReversion[u]]-signedSeminormSquare[u]];

(* Clifford inverse of an element (assuming it's invertible). *)
cliffordInverse::usage="cliffordInverse[u] - returns \!\(\*FractionBox[OverscriptBox[\(u\), \(~\)], \(u \*OverscriptBox[\(u\), \(~\)]\)]\), the multiplicative inverse of u \[Element] \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\) when u is invertible.";
cliffordInverse[B_]:=Chop[Expand[cliffordReversion[B]/signedSeminormSquare[B]]];


(* grade involution of inverse: Overscript[u, _]/(uOverscript[u), ~]. *)
conjugatePseudoInverse::usage="conjugatePseudoInverse[u] - returns grade involution of \!\(\*SuperscriptBox[\(u\), \(-1\)]\), i.e.,  \!\(\*FractionBox[OverscriptBox[\(u\), \(_\)], \(u \*OverscriptBox[\(u\), \(~\)]\)]\)   (u \[Element] \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\)).";
conjugatePseudoInverse[B_]:=Chop[Expand[cliffordConjugate[B]/signedSeminormSquare[B]]];

(* ALGEBRAIC PROPERIES OF ZEONS *)

(*Zeon real and dual parts *)
zRe::usage="zRe[u] - Returns real part, \[GothicCapitalR]u (`body'), of u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\)";
zRe[u_]:=gradeKPart[u,0];

zDu::usage="zDu[u] - Returns dual part, \[GothicCapitalD]u (`soul'), of u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\)";
zDu[u_]:=Chop[u-zRe[u]];

(* The zeon absolute value is defined as the principal square root of (u^2). *)

zeonAbs::usage="zeonAbs[u] - returns `zeon absolute value' of u: = \!\(\*SqrtBox[SuperscriptBox[\(u\), \(2\)]]\).";
zeonAbs[u_]:=If[nilpotentZeonQ[u],u,Expand[(zRe[u]/Abs[zRe[u]])u]];

(* Compute zeon inverses *)
zeonInverse::usage="zeonInverse[u] - Multiplicative inverse of u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\), provided Re[u]\[NotEqual]0.";
zeonInverse[u_]:=Module[{a,b,\[Beta],inv},
a=gradeKPart[u,0];
If[PossibleZeroQ[a],Message[zeonInverse::noninv,u];,
b=Chop[u-a];
\[Beta]=Expand[-a^-2b];
inv=a^-1;
While[clGrade[\[Beta]]>0,
inv=inv+\[Beta];
\[Beta]=clExpand[(-a^-1b)\[CircleDot]\[Beta]];
];
Return[inv];];
];
zeonInverse::noninv="\[GothicCapitalR]u=0. The argument `1` is not invertible.";

(* Obtain maximum generator index appearing in an arbitrary zeon *)
zeonMaxGenerator::usage="zeonMaxGenerator[u] - Maximum generator index appearing in u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
zeonMaxGenerator[u_+v_]:=Max[{zeonMaxGenerator[u],zeonMaxGenerator[v]}];
zeonMaxGenerator[Subscript[\[Zeta], a_]]:=Max[a];
zeonMaxGenerator[x_]:=0/;(FreeQ[x,Subscript[\[Zeta], _]]);
zeonMaxGenerator[x_ Subscript[\[Zeta], a_]]:=Max[a]/;(FreeQ[x,Subscript[\[Zeta], _]]);

(* Zeon left contraction operator *)
zeonLeftContraction::usage="zeonLeftContraction[u,v] - Zeon contraction u\[RightCeiling]v defined by linear extension of \!\(\*SubscriptBox[\(\[Zeta]\), \(I\)]\)\[RightCeiling]\!\(\*SubscriptBox[\(\[Zeta]\), \(J\)]\) = \!\(\*SubscriptBox[\(\[Zeta]\), \(J\[Backslash]I\)]\) if I \[SubsetEqual] J, 0 otherwise.";
zeonLeftContraction[x_ Subscript[\[Zeta], a_], y_ Subscript[\[Zeta], J_]]:=Expand[x y zeonLeftContraction[Subscript[\[Zeta], a],Subscript[\[Zeta], J]]]/;(FreeQ[x,Subscript[\[Zeta], _]]\[And]FreeQ[y,Subscript[\[Zeta], _]]);
zeonLeftContraction[x_, Subscript[\[Zeta], J_]]:=0/;(FreeQ[x,Subscript[\[Zeta], _]]);
zeonLeftContraction[Subscript[\[Zeta], J_], y_]:=0/;(FreeQ[y,Subscript[\[Zeta], _]]);
zeonLeftContraction[Subscript[\[Zeta], a_], y_ Subscript[\[Zeta], J_]]:=y zeonLeftContraction[Subscript[\[Zeta], a],Subscript[\[Zeta], J]]/;(FreeQ[y,Subscript[\[Zeta], _]]);
zeonLeftContraction[x_ Subscript[\[Zeta], J_], Subscript[\[Zeta], a_]]:=x zeonLeftContraction[Subscript[\[Zeta], J],Subscript[\[Zeta], a]]/;(FreeQ[x,Subscript[\[Zeta], _]]);
zeonLeftContraction[Subscript[\[Zeta], a_],Subscript[\[Zeta], J_]]:=If[a==J,1,If[SubsetQ[J,a],Subscript[\[Zeta], Complement[J,a]],0]];
zeonLeftContraction[u_+v_, w_]:=zeonLeftContraction[u,w]+zeonLeftContraction[v,w];
zeonLeftContraction[u_,v_+w_]:=zeonLeftContraction[u,v]+zeonLeftContraction[u,w];

(* Zeon right contraction operator *)
zeonRightContraction::usage="zeonRightContraction[u,v] - Zeon contraction u\[LeftCeiling]v defined by linear extension of \!\(\*SubscriptBox[\(\[Zeta]\), \(I\)]\)\[LeftCeiling]\!\(\*SubscriptBox[\(\[Zeta]\), \(J\)]\) = \!\(\*SubscriptBox[\(\[Zeta]\), \(I\[Backslash]J\)]\) if J \[SubsetEqual] I, 0 otherwise.";
zeonRightContraction[x_ Subscript[\[Zeta], a_], y_ Subscript[\[Zeta], J_]]:=Expand[x y zeonRightContraction[Subscript[\[Zeta], a],Subscript[\[Zeta], J]]]/;(FreeQ[x,Subscript[\[Zeta], _]]\[And]FreeQ[y,Subscript[\[Zeta], _]]);
zeonRightContraction[x_, Subscript[\[Zeta], J_]]:=0/;(FreeQ[x,Subscript[\[Zeta], _]]);
zeonRightContraction[Subscript[\[Zeta], J_], y_]:=0/;(FreeQ[y,Subscript[\[Zeta], _]]);
zeonRightContraction[Subscript[\[Zeta], a_], y_ Subscript[\[Zeta], J_]]:=y zeonRightContraction[Subscript[\[Zeta], a],Subscript[\[Zeta], J]]/;(FreeQ[y,Subscript[\[Zeta], _]]);
zeonRightContraction[x_ Subscript[\[Zeta], J_], Subscript[\[Zeta], a_]]:=x zeonRightContraction[Subscript[\[Zeta], J],Subscript[\[Zeta], a]]/;(FreeQ[x,Subscript[\[Zeta], _]]);
zeonRightContraction[Subscript[\[Zeta], a_],Subscript[\[Zeta], J_]]:=If[a==J,1,If[SubsetQ[a,J],Subscript[\[Zeta], Complement[a,J]],0]];
zeonRightContraction[u_+v_, w_]:=zeonRightContraction[u,w]+zeonRightContraction[v,w];
zeonRightContraction[u_,v_+w_]:=zeonRightContraction[u,v]+zeonRightContraction[u,w];

(* Given zeon u with maximum generator m, find \[Omega], \[Psi] such that u = \[Omega] + Subscript[\[Psi]\[Zeta], {m}] *)
zeonMaxGenDecomp::usage="zeonMaxGenDecomp[u] - Given zeon u with maximum generator index m, return {\[Omega],\[Psi]} such that u = \[Omega] + \!\(\*SubscriptBox[\(\[Psi]\[Zeta]\), \({m}\)]\).";
zeonMaxGenDecomp[u_]:=Module[{\[Psi],\[Omega],m},
m=zeonMaxGenerator[u];
\[Psi]=zeonRightContraction[u,\!\(\*SubscriptBox[\(\[Zeta]\), \({m}\)]\)];
\[Omega]=u-clExpand[\!\(\*SubscriptBox[\(\[Zeta]\), \({m}\)]\)\[CircleDot]\[Psi]];
Return[{\[Omega],\[Psi]}]];

(*Zeon Square Root *)
zeonSqRoot::usage="zeonSqRoot[u] Recursive computation of principal square root of u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\), provided Re[u] \[GreaterEqual] 0.";
zeonSqRoot[u_]:=If[zeonMaxGenerator[u]==0,u^(1/2),clExpand[zeonSqRoot[zeonMaxGenDecomp[u][[1]]]+clExpand[\!\(\*SubscriptBox[\(\[Zeta]\), \({zeonMaxGenerator[u]}\)]\)\[CircleDot]Expand[1/2 zeonInverse[zeonSqRoot[zeonMaxGenDecomp[u][[1]]]]]\[CircleDot]zeonMaxGenDecomp[u][[2]]]]];
 
(* Zeon kth root - often slower than zeonPower[u,1/k]*)
zeonKthRoot::usage="zeonKthRoot[u] Recursive computation of principal kth root of u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\), provided Re[u] \[GreaterEqual] 0.";
zeonKthRoot[u_,k_]:=If[zeonMaxGenerator[u]==0,u^(1/k),ExpandAll[clExpand[zeonKthRoot[zeonMaxGenDecomp[u][[1]],k]]+clExpand[\!\(\*SubscriptBox[\(\[Zeta]\), \({zeonMaxGenerator[u]}\)]\)\[CircleDot]ExpandAll[1/k cliffordPower[zeonInverse[zeonKthRoot[zeonMaxGenDecomp[u][[1]],k]],k-1]]\[CircleDot]zeonMaxGenDecomp[u][[2]]]]];

(* Zeon rth power: r can be any real number *)
zeonPower::usage="zeonPower[u,r] returns \!\(\*SuperscriptBox[\(u\), \(r\)]\) (rth power of zeon u) for r\[Element]\!\(\*TemplateBox[{},\n\"Reals\"]\) and u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).  zRe[u]\[NotEqual]0 required if r<0.";
zeonPower[u_,r_]:=Module[{x0,Du,x},
x0=zRe[u];
Du=u-x0;
Return[ExpandAll[Sum[1/k! (D[x^r,{x,k}]/.{x->x0})clExpand[cliffordPower[Du,k]],{k,0,maxIndex}]]]];


(* General zeon exponential, sine, cosine, and tangent *)
zeonExp::usage="zeonExp[u] returns the exponential of u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
zeonExp[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Exp[x],x,z]]];

zeonLog::usage="zeonLog[u] Evaluate logarithm of u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\), provided zRe[u]>0.";
zeonLog[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Log[x],x,z]]];

zeonSin::usage="zeonSin[u] evaluates the sine function at u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
zeonSin[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Sin[x],x,z]]];

zeonCos::usage="zeonCos[u] evaluates the cosine function at u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
zeonCos[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Cos[x],x,z]]];

zeonTan::usage="zeonTan[u] evaluates the tangent function at u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\), provided cos(u) is invertible; i.e., Re[cos(u)] \[NotEqual] 0.";
zeonTan[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Tan[x],x,z]]];

zeonArcSin::usage="zeonArcSin[u] Evaluates arcsine at u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\), provided -1 \[LessEqual] zRe[u] \[LessEqual] 1.";
zeonArcSin[z_]:=Module[{x},zeonArcSin[z]=zeonExtensionEval[ArcSin[x],x,z]];

zeonArcCos::usage="zeonArcCos[u] Evaluates arccosine at u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\), provided -1 \[LessEqual] zRe[u] \[LessEqual] 1.";
zeonArcCos[z_]:=Module[{x},zeonArcCos[z]=zeonExtensionEval[ArcCos[x],x,z]];

zeonArcTan::usage="zeonArcTan[u] Evaluates arctangent at u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
zeonArcTan[z_]:=Module[{x},zeonArcTan[z]=zeonExtensionEval[ArcTan[x],x,z]];

(* Zeon hyperbolic functions *)
zeonSinh::usage="zeonSinh[u] - Hyperbolic sine evaluated at u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
zeonSinh[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Sinh[x],x,z]]];
zeonCosh::usage="zeonCosh[u] - Hyperbolic cosine evaluated at u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
zeonCosh[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Cosh[x],x,z]]];
zeonTanh::usage="zeonTanh[u] - Hyperbolic tangent evaluated at u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
zeonTanh[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Tanh[x],x,z]]];

(* Zeon-Hodge dual operator.  Multiindices are replaced with their relative complements. *)
zeonDual::usage="zeonDual[u] - Returns Zeon-Hodge dual \!\(\*SuperscriptBox[\(u\), \(\[Star]\)]\) given by linear extension of \!\(\*SuperscriptBox[SubscriptBox[\(\[Zeta]\), \(I\)], \(\[Star]\)]\)=\!\(\*SubscriptBox[\(\[Zeta]\), \(\([n]\)\[Backslash]I\)]\) in \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).
Satisfies <\[Alpha] \[CircleDot] \!\(\*SuperscriptBox[\(\[Beta]\), \(*\)]\)\!\(\*SubscriptBox[\(>\), \(n\)]\) = <\[Alpha], \[Beta]> \!\(\*SubscriptBox[\(\[Zeta]\), \([n]\)]\).";
SetAttributes[zeonDual,Listable];
zeonDual[u_+v_]:=Block[{$RecursionLimit=recLimit},zeonDual[u]+zeonDual[v]];
zeonDual[x_]:=x Subscript[\[Zeta], Range[maxIndex]]/;(FreeQ[x,Subscript[\[Zeta], _]]);
zeonDual[x_ Subscript[\[Zeta], a_]]:=x zeonDual[Subscript[\[Zeta], a]]/;(FreeQ[x,Subscript[\[Zeta], _]]);
zeonDual[Subscript[\[Zeta], a_]]:=Subscript[\[Zeta], Complement[Range[maxIndex],a]]/;(Length[a]<maxIndex);
zeonDual[Subscript[\[Zeta], a_]]:=1/;(Length[a]>=maxIndex);

(* Minimal grade of zeon u.   Returns 0 if u is trivial, otherwise returns minimal grade of the dual part of u. *)
zeonMinGrade::usage="zeonMinGrade[u] Returns 0 if u is trivial (i.e, real), otherwise returns minimal grade of the dual part of u.";
zeonMinGrade[u_]:=If[FreeQ[u,Subscript[\[Zeta], _]],0,maxIndex-clGrade[zeonDual[zDu[u]]]];

(*zeonDivisionAlgorithm[u,v] returns q,r \[Element] Subscript[Cl, n]^nil such that u = qv+r, where (i.) qv is homogeneously decomposable and (ii.) r=0 or zeonMinGrade[r] > zeonMinGrade[u].*)
zeonDivisionAlgorithm::usage="zeonDivisionAlgorithm[u,v] returns {q,r} \[Subset] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\) such that u = qv+r, where (i.) qv is homogeneously decomposable and (ii.) r=0 or zeonMinGrade[r] > zeonMinGrade[u].";
zeonDivisionAlgorithm[u_,v_]:=Module[{q,r,factormonomials,a,m,monomials,normalized,qv,qr},
a=zRe[u];
m=zeonMinGrade[u];
monomials=DeleteCases[gradeKPart[zeonBasis,m],0];
normalized=ExpandAll[u/a];
factormonomials=Coefficient[zDu[normalized],monomials]monomials;
q=a clExpand[CircleDot[zeonInverse[v],cliffordListProduct[1+factormonomials]]];
qv=clExpand[CircleDot[q,v]];
If[Variables[Chop[u-qv]]=={},qr={q,0},
q=Expand[a/2 clExpand[ CircleDot[zeonInverse[v],cliffordListProduct[1+2factormonomials]]]];
r=Expand[u-clExpand[CircleDot[q,v]]];
qr={q,r};];
Return[{q,r}]];

(*ZeonElementaryFactorization[u] returns factors a, {(1+Subscript[a, I]Subscript[\[Zeta], I]): I\[Element]2^[n], I\[NotEqual]\[EmptySet]} such that u = a\[Product] (1+Subscript[a, I]Subscript[\[Zeta], I])*)
zeonElementaryFactorization::usage="zeonElementaryFactorization[u] returns factors a, {(1+\!\(\*SubscriptBox[\(a\), \(I\)]\)\!\(\*SubscriptBox[\(\[Zeta]\), \(I\)]\)): I\[Element]\!\(\*SuperscriptBox[\(2\), \([n]\)]\), I\[NotEqual]\[EmptySet]} such that u = a\[Product] (1+\!\(\*SubscriptBox[\(a\), \(I\)]\)\!\(\*SubscriptBox[\(\[Zeta]\), \(I\)]\))";
zeonElementaryFactorization[u_]:=Module[{a,factormonomials,factors,m,monomials,normalized},
a=zRe[u];
m=zeonMinGrade[u];
factors={a};
monomials=DeleteCases[gradeKPart[zeonBasis,m],0];
normalized=ExpandAll[u/a];
factormonomials=Coefficient[zDu[normalized],monomials]monomials;
factors=Union[factors,1+factormonomials];
q=clExpand[CircleDot[normalized,cliffordListProduct[1-factormonomials]]];
m=zeonMinGrade[q];
While[m>0,
monomials=DeleteCases[gradeKPart[zeonBasis,m],0];
factormonomials=Coefficient[zDu[q],monomials]monomials;
factors=Union[factors,1+factormonomials];
q=clExpand[CircleDot[q,cliffordListProduct[1-factormonomials]]];
m=zeonMinGrade[q];
];
Return[factors]];

(* Zeon powers u^v, where u,v are both zeons and Re u > 0. *)
zeonPowerZeon::usage="zeonPowerZeon[u,v] returns zeon to zeon power \!\(\*SuperscriptBox[\(u\), \(v\)]\), where u,v \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\) and Re u > 0.";
zeonPowerZeon[u_,v_]:=clExpand[zeonExp[clExpand[v\[CircleDot]zeonLog[u]]]];

(* Zeon series representation of suitably differentiable function f evaluated at u *) 
zeonExtensionEval::usage="zeonExtensionEval[f[x],x,u] evaluates the zeon extension \[Phi] (of suitably differentiable) real function f at u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
zeonExtensionEval[f_,var_,u_]:=Module[{x0,Du},
x0=zRe[u];
Du=u-x0;
Return[ExpandAll[Sum[(D[f,{var,k}]/.{var->x0})/k! clExpand[cliffordPower[Du,k]],{k,0,maxIndex}]]]];

(* Zeon series representation of suitably differentiable function f evaluated at u *) 
zeonExtension::usage="zeonExtension[f[x],x] returns \[Phi][u], the series representation of zeon extension \[Phi] of (suitably differentiable) real function f.";
zeonExtension[f_,var_]:=Module[{Ru,Du},
Ru=Symbol["Ru"];
Du=Symbol["Du"];
Return[ExpandAll[Sum[(D[f,{var,k}]/.{var->Ru})/k! (Du)^k,{k,0,maxIndex}]]]];


(* Algebraic derivative of zeon-valued function \[Phi] of a zeon variable, induced by real function f, evaluated at u *)
zeonDerivative::usage="zeonDerivative[f[x],x,u,k] returns \!\(\*SuperscriptBox[\(\[Phi]\), \((k)\)]\)(u): the kth order derivative of zeon extension \[Phi] of real function f, evaluated at u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).  If no order is specified, k=1 is assumed.";
zeonDerivative[f_,var_,u_,ord_:1]:=Module[{x0,Du},
x0=zRe[u];
Du=u-x0;
Return[ExpandAll[Sum[(D[f,{var,k+ord}]/.{var->x0})/k! clExpand[cliffordPower[Du,k]],{k,0,maxIndex}]]]];

(* Algebraic antiderivative of zeon-valued function \[Phi] of a zeon variable, induced by real function f, evaluated at u *)
zeonAntiDerivative::usage="zeonAntiDerivative[f[x],x,u] returns \[CapitalPhi][u], where \[CapitalPhi][x] is the zeon extension of real function F[x], the antiderivative of f[x].";
zeonAntiDerivative[f_,var_,u_]:=Module[{x0,Du,F},
x0=zRe[u];
Du=u-x0;
F=Integrate[f,var];
Return[ExpandAll[Sum[(D[F,{var,k}]/.{var->x0})/k! clExpand[cliffordPower[Du,k]],{k,0,maxIndex}]]]];

(* Algebraic derivative of zeon-valued function \[Phi] of a zeon variable, induced by real function f, evaluated at u *)
generalZeonDerivative::usage="generalZeonDerivative[f[x],x,k] returns \!\(\*SuperscriptBox[\(\[Phi]\), \((k)\)]\)(u): the kth order derivative of zeon extension \[Phi] of real function f.  If no order is specified,  k=1 is assumed.";
generalZeonDerivative[f_,var_,ord_:1]:=Module[{\[GothicCapitalR]u,\[GothicCapitalD]u,k},
\[GothicCapitalR]u=Symbol["\[GothicCapitalR]u"];
\[GothicCapitalD]u=Symbol["\[GothicCapitalD]u"];
Return[ExpandAll[Sum[(D[f,{var,k+ord}]/.{var->\[GothicCapitalR]u})/k! (\[GothicCapitalD]u)^k,{k,0,maxIndex}]]]];

(* General antiderivative of zeon-valued function \[Phi] of a zeon variable, induced by real function f, evaluated at u *)
generalZeonAntiDerivative::usage="generalZeonAntiDerivative[f[x],x] returns \[CapitalPhi][u]+k, where \[CapitalPhi][u] is the zeon extension of real function F[x], the antiderivative of f[x].";
generalZeonAntiDerivative[f_,var_]:=Module[{\[GothicCapitalR]u,\[GothicCapitalD]u,g,\[Kappa],k},
g=Integrate[f,var];
\[GothicCapitalR]u=Symbol["\[GothicCapitalR]u"];
\[GothicCapitalD]u=Symbol["\[GothicCapitalD]u"];
\[Kappa]=Symbol["\[Kappa]"];
Return[ExpandAll[Sum[(D[g,{var,k}]/.{var->\[GothicCapitalR]u})/k! (\[GothicCapitalD]u)^k,{k,0,maxIndex}]+\[Kappa]]]];

(* Evaluate a quadratic function \[Phi][x] of variable x at zeon u *)
zeonQuadraticEval::usage="zeonQuadraticEval[\[Phi][x],x,u] evaluates the quadratic function \[Phi][x]=\!\(\*SuperscriptBox[\(\[Alpha]x\), \(2\)]\)+\[Beta]x+\[Gamma] at u, where \[Alpha],\[Beta],\[Gamma],u \[CurlyEpsilon] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
zeonQuadraticEval[\[Phi]_,x_,u_]:=Module[{\[Alpha],\[Beta],\[Gamma]},
\[Alpha]=Coefficient[\[Phi],x^2];
\[Beta]=Coefficient[\[Phi],x];
\[Gamma]=\[Phi]/.{x->0};
Return[clExpand[\[Alpha]\[CircleDot]u\[CircleDot]u+\[Beta]\[CircleDot]u+\[Gamma]]]];

(* Compute the discriminant of quadratic polynomial f[x] having zeon coefficients *)
zeonDiscriminant::usage="zeonDiscriminant[f[x],x] computes the zeon discriminant of \[Phi][x]=\!\(\*SuperscriptBox[\(\[Alpha]x\), \(2\)]\)+\[Beta]x+\[Gamma], where \[Alpha],\[Beta],\[Gamma] \[CurlyEpsilon] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
zeonDiscriminant[f_,x_]:=Module[{\[Alpha],\[Beta],\[Gamma]},
\[Alpha]=Coefficient[f,x^2];
\[Beta]=Coefficient[f,x];
\[Gamma]=f/.{x->0};
Return[Expand[clExpand[\[Beta]\[CircleDot]\[Beta]]-4clExpand[\[Alpha]\[CircleDot]\[Gamma]]]];
];

(* Apply the zeon quadratic formula to \[Phi][x] *)
zeonQuadraticFormula::usage="zeonQuadraticFormula[\[Phi][x],x] applies the zeon quadratic formula to \[Phi][x]=\!\(\*SuperscriptBox[\(\[Alpha]x\), \(2\)]\)+\[Beta]x+\[Gamma], 
where \[Alpha],\[Beta],\[Gamma] \[CurlyEpsilon] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\). When \[CapitalDelta]=0, includes arbitrary \[Kappa] assumed to be nilpotent of index 2.";
zeonQuadraticFormula[\[Phi]_,x_]:=Module[{\[Alpha],\[Beta],\[Gamma],\[CapitalDelta]},
\[Alpha]=Coefficient[\[Phi],x^2];
\[Beta]=Coefficient[\[Phi],x];
\[Gamma]=\[Phi]/.{x->0};
\[CapitalDelta]=zeonDiscriminant[\[Phi],x];
If[zRe[\[Alpha]]==0,Return[{}],
If[zeonNorm[\[CapitalDelta],1]==0,Return[-clExpand[\[Beta]\[CircleDot]Expand[zeonInverse[2\[Alpha]]]]+ \[Kappa]],
If[zRe[\[CapitalDelta]]>0,Return[{clExpand[(-\[Beta]+zeonSqRoot[\[CapitalDelta]])\[CircleDot]zeonInverse[Expand[2\[Alpha]]]],
clExpand[(-\[Beta]-zeonSqRoot[\[CapitalDelta]])\[CircleDot]zeonInverse[Expand[2\[Alpha]]]]}],
If[zRe[\[CapitalDelta]]<0,Return[{}],Return[Indeterminate]]]]]];

(* Inverses of zeon matrices *)

zeonMatrixInverse::usage="zeonMatrixInverse[A] - Given square zeon matrix A, returns zeon matrix U such that AU=UA=I, provided |\[GothicCapitalR]A|\[NotEqual]0.";
zeonMatrixInverse[A_]:=Module[{ra,da,X,ira,U},
ra=zRe[A];
da=zDu[A];
If[!SquareMatrixQ[A],
Print["Matrix is not square!"];,
If[Det[ra]!=0,
ira=Inverse[ra];
X=clExpand[cliffordMatrixProduct[da,Inverse[ra]]];
clExpand[cliffordMatrixProduct[Inverse[ra],Sum[clExpand[cliffordMatrixPower[-X,ell]],{ell,0,maxIndex}]]],Print["Matrix is not invertible!"];]]];

(* Determine whether or not a zeon element is nilpotent. *)
nilpotentZeonQ::usage="nilpotentZeonQ[u] returns True if u is nilpotent; False if u is invertible.";
nilpotentZeonQ[u_]:=If[PossibleZeroQ[zRe[u]],True,False];

(* Determine the index of nilpotency of a zeon element (or its dual part, if invertible). *)
zeonNilpotency::usage="zeonNilpotency[u] returns least positive integer \[Kappa] such that 
 (\[GothicCapitalD]u\!\(\*SuperscriptBox[\()\), \(\[Kappa]\)]\)=0.";
zeonNilpotency[u_]:= 1+Exponent[zeonExp[t zDu[u]],t];

(* Zeon geometric series *)
zeonGeometricSeries::usage="zeonGeometricSeries[a,u,k] - returns partial sum \!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 0\), \(k\)]\)\!\(\*SuperscriptBox[\(au\), \(j\)]\) of geometric series \!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 0\), \(\[Infinity]\)]\)\!\(\*SuperscriptBox[\(au\), \(j\)]\).  
If k is omitted and |\[GothicCapitalR]u|<1, returns limit of geometric series: \!\(\*FractionBox[\(a\), \(1 - \[GothicCapitalR]u\)]\)\!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 0\), \(maxIndex\)]\)(1-\[GothicCapitalR]u\!\(\*SuperscriptBox[\()\), \(-j\)]\)(\[GothicCapitalD]u\!\(\*SuperscriptBox[\()\), \(j\)]\).";
zeonGeometricSeries[a_,u_,k_:\[Infinity]]:=If[k<\[Infinity], Expand[a Sum[clExpand[zeonPower[u,j]],{j,0,k}]],If[Abs[zRe[u]]<1,
Expand[a/(1-zRe[u])Sum[clExpand[(1-zRe[u])^-j zeonPower[zDu[u],j]],{j,0,maxIndex}]]]];

(* Matrix representations *)

etaMatrix::usage="etaMatrix - Returns 2\[Cross]2 matrix \[Eta] = \!\(\*SubscriptBox[\(\[Sigma]\), \(z\)]\) + \!\(\*SubscriptBox[\(\[ImaginaryI]\[Sigma]\), \(y\)]\) used to construct representations of zeon generators.";
etaMatrix=PauliMatrix[3]+I PauliMatrix[2];

repMatrixProduct::usage="repMatrixProduct[M,q] - Compute ordered product of matrices in list M containing q matrices.";
repMatrixProduct[M_,q_]:=If[q==1,M[[1]],M[[1]].repMatrixProduct[M[[2;;q]],q-1]];

matrixKroneckerPower::usage="matrixKroneckerPower[M,k] - returns \!\(\*
StyleBox[\"k\",\nFontSlant->\"Italic\"]\)th Kronecker power of Matrix M, i.e., \!\(\*SuperscriptBox[\(M\), \(\[CircleTimes]k\)]\).";
matrixKroneckerPower[M_,pwr_]:=If[pwr==0,{1},Expand[KroneckerProduct[M,matrixKroneckerPower[M,pwr-1]]]];

cliffordGeneratorRep::usage="cliffordGeneratorRep[j] - \!\(\*SuperscriptBox[\(2\), \(n\)]\)\[Cross]\!\(\*SuperscriptBox[\(2\), \(n\)]\) real matrix representation of generator \!\(\*SubscriptBox[\(e\), \({j}\)]\)\[Element]\!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\):  
\!\(\*SubscriptBox[\(e\), \({j}\)]\)\[Rule] \!\(\*SuperscriptBox[SubscriptBox[\(\[Sigma]\), \(x\)], \(\[CircleTimes]\((j - 1)\)\)]\) \[CircleTimes] \!\(\*SubscriptBox[\(\[Sigma]\), \(z\)]\) \[CircleTimes] \!\(\*SuperscriptBox[SubscriptBox[\(\[Sigma]\), \(0\)], \(\[CircleTimes]\((n - j)\)\)]\) (1 \[LessEqual] j \[LessEqual]p)
\!\(\*SubscriptBox[\(e\), \({j}\)]\)\[Rule] \!\(\*SuperscriptBox[SubscriptBox[\(\[Sigma]\), \(x\)], \(\[CircleTimes]\((j - 1)\)\)]\) \[CircleTimes] \!\(\*SubscriptBox[\(\[ImaginaryI]\[Sigma]\), \(y\)]\) \[CircleTimes] \!\(\*SuperscriptBox[SubscriptBox[\(\[Sigma]\), \(0\)], \(\[CircleTimes]\((n - j)\)\)]\) (p < j \[LessEqual] p+q)
\!\(\*SubscriptBox[\(e\), \({j}\)]\)\[Rule] \!\(\*SuperscriptBox[SubscriptBox[\(\[Sigma]\), \(x\)], \(\[CircleTimes]\((j - 1)\)\)]\) \[CircleTimes] \[Eta] \[CircleTimes] \!\(\*SuperscriptBox[SubscriptBox[\(\[Sigma]\), \(0\)], \(\[CircleTimes]\((n - j)\)\)]\) (p+q < j \[LessEqual]n)
where \[Eta]=\!\(\*SubscriptBox[\(\[Sigma]\), \(z\)]\)+\!\(\*SubscriptBox[\(\[ImaginaryI]\[Sigma]\), \(y\)]\), and \!\(\*SubscriptBox[\(\[Sigma]\), \(0\)]\), \!\(\*SubscriptBox[\(\[Sigma]\), \(y\)]\), and \!\(\*SubscriptBox[\(\[Sigma]\), \(z\)]\) are Pauli matrices.";
cliffordGeneratorRep[j_]:=If[j<=pPart,KroneckerProduct[KroneckerProduct[matrixKroneckerPower[PauliMatrix[1],j-1],PauliMatrix[3]],IdentityMatrix[2^(maxIndex-j)]],
If[j<=pPart+qPart,KroneckerProduct[KroneckerProduct[matrixKroneckerPower[PauliMatrix[1],j-1],I PauliMatrix[2]],IdentityMatrix[2^(maxIndex-j)]],
KroneckerProduct[KroneckerProduct[matrixKroneckerPower[PauliMatrix[1],j-1],etaMatrix],IdentityMatrix[2^(maxIndex-j)]]]];

cliffordRep::usage="cliffordRep[u] - Returns \!\(\*SuperscriptBox[\(2\), \(n\)]\)\[Cross]\!\(\*SuperscriptBox[\(2\), \(\(n\)\(\\\ \)\)]\)(real) matrix representation of u \[Element] \!\(\*SubscriptBox[\(Cl\), \(p, q, r\)]\).";
cliffordRep[x_]:=Expand[x IdentityMatrix[2^maxIndex]]/;FreeQ[x,e];
cliffordRep[Subscript[e, I_]]:=Module[{M},
M=Table[cliffordGeneratorRep[j],{j,I}];
Return[repMatrixProduct[M,Length[I]]];
];
cliffordRep[x_ Subscript[e, I_]]:=Module[{M},
M=Table[cliffordGeneratorRep[j],{j,I}];
Return[x repMatrixProduct[M,Length[I]]];
]/;FreeQ[x,e];
cliffordRep[x_+y_]:=cliffordRep[x]+cliffordRep[y];


zeonFrobeniusNorm::usage="zeonFrobeniusNorm[A] - returns \!\(\*SqrtBox[\(\[GothicCapitalR] \((\*SuperscriptBox[\(tr[\*SuperscriptBox[\(A\), \(\[Dagger]\)] A]\), \(*\)])\)\)]\) = \!\(\*SqrtBox[\(\[Sum]\\\ \*SuperscriptBox[SubscriptBox[\(\[LeftDoubleBracketingBar]\*SubscriptBox[\(a\), \(ij\)]\[RightDoubleBracketingBar]\), \(2\)], \(2\)]\)]\).";
zeonFrobeniusNorm[A_]:=Sqrt[Total[Flatten[zeonInnerProduct[A,A]]]];


zeonGeneratorRep::usage="zeonGeneratorRep[j] - \!\(\*SuperscriptBox[\(2\), \(n\)]\)\[Cross]\!\(\*SuperscriptBox[\(2\), \(n\)]\) matrix representation of zeon generator \!\(\*SubscriptBox[\(\[Zeta]\), \({j}\)]\)\[Rule] \!\(\*SuperscriptBox[SubscriptBox[\(\[Sigma]\), \(0\)], \(\[CircleTimes]\((j - 1)\)\)]\) \[CircleTimes] \[Eta] \[CircleTimes] \!\(\*SuperscriptBox[SubscriptBox[\(\[Sigma]\), \(0\)], \(\[CircleTimes]\((n - j)\)\)]\),
where \[Eta]=\!\(\*SubscriptBox[\(\[Sigma]\), \(z\)]\)+\!\(\*SubscriptBox[\(\[ImaginaryI]\[Sigma]\), \(y\)]\), and \!\(\*SubscriptBox[\(\[Sigma]\), \(0\)]\), \!\(\*SubscriptBox[\(\[Sigma]\), \(y\)]\), and \!\(\*SubscriptBox[\(\[Sigma]\), \(z\)]\) are Pauli matrices.";
zeonGeneratorRep[j_]:=KroneckerProduct[KroneckerProduct[IdentityMatrix[2^(j-1)],etaMatrix],IdentityMatrix[2^(maxIndex-j)]];

zeonRep::usage="zeonRep[u] - Returns \!\(\*SuperscriptBox[\(2\), \(n\)]\)\[Cross]\!\(\*SuperscriptBox[\(2\), \(\(n\)\(\\\ \)\)]\)real matrix representation of u \[Element] \!\(\*SuperscriptBox[SubscriptBox[\(Cl\), \(n\)], \(nil\)]\).";
zeonRep[x_]:=Expand[x IdentityMatrix[2^maxIndex]]/;FreeQ[x,\[Zeta]];
zeonRep[Subscript[\[Zeta], I_]]:=Module[{M},
M=Table[zeonGeneratorRep[j],{j,I}];
Return[repMatrixProduct[M,Length[I]]];
];
zeonRep[x_ Subscript[\[Zeta], I_]]:=Module[{M},
M=Table[zeonGeneratorRep[j],{j,I}];
Return[x repMatrixProduct[M,Length[I]]];
]/;FreeQ[x,\[Zeta]];
zeonRep[x_+y_]:=zeonRep[x]+zeonRep[y];

generalZeonElement::usage="generalZeonElement[\[Alpha]] - Canonical expansion\!\(\*UnderscriptBox[\(\[Sum]\), \(I \[Element] \*SuperscriptBox[\(2\), \([n]\)]\)]\)\!\(\*SubscriptBox[\(\[Alpha]\), \(I\)]\)\!\(\*SubscriptBox[\(\[Zeta]\), \(I\)]\), where \[Alpha] is a symbol.";
generalZeonElement[var_]:=Prepend[(clBasis[[2;;]]/.{e->var}),Subscript[var, 0]].zeonBasis;

generalCliffordElement::usage="generalCliffordElement[\[Alpha]] - Canonical expansion\!\(\*UnderscriptBox[\(\[Sum]\), \(I \[Element] \*SuperscriptBox[\(2\), \([n]\)]\)]\)\!\(\*SubscriptBox[\(\[Alpha]\), \(I\)]\)\!\(\*SubscriptBox[\(e\), \(I\)]\), where \[Alpha] is a symbol.";
generalCliffordElement[var_]:=Prepend[(clBasis[[2;;]]/.{e->var}),Subscript[var, 0]].clBasis;

generalIdemElement::usage="generalIdemElement[\[Alpha]] - Canonical expansion\!\(\*UnderscriptBox[\(\[Sum]\), \(I \[Element] \*SuperscriptBox[\(2\), \([n]\)]\)]\)\!\(\*SubscriptBox[\(\[Alpha]\), \(I\)]\)\!\(\*SubscriptBox[\(\[CurlyEpsilon]\), \(I\)]\), where \[Alpha] is a symbol.";
generalIdemElement[var_]:=Prepend[(clBasis[[2;;]]/.{e->var}),Subscript[var, 0]].idemBasis;

generalSymElement::usage="generalSymElement[\[Alpha]] - Canonical expansion\!\(\*UnderscriptBox[\(\[Sum]\), \(I \[Element] \*SuperscriptBox[\(2\), \([n]\)]\)]\)\!\(\*SubscriptBox[\(\[Alpha]\), \(I\)]\)\!\(\*SubscriptBox[\(\[Stigma]\), \(I\)]\), where \[Alpha] is a symbol.";
generalSymElement[var_]:=Prepend[(clBasis[[2;;]]/.{e->var}),Subscript[var, 0]].symBasis;

];

EndPackage[];





