CliffMath2018

Clifford algebra computations in Mathematica

G. Stacey Staples

Southern Illinois University Edwardsville

June  2019
v. 20190605

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[ζ];
ζ=Symbol["ζ"];
Protect[ζ];
Unprotect[ϛ];
ϛ=Symbol["ϛ"];
Protect[ϛ];
Unprotect[ε];
ε=Symbol["ε"];
Protect[ε];
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 (Subscript[Cl, p, q, r]), as well as `commutative Clifford' multiplication in ((Subscript[Cl, p, q, r])^sym).  For n=p+q+r, initializes exterior (Grassmann) multiplication in ((ΛR)^n), `idempotent-Clifford' multiplication in ((Subscript[Cl, n])^idem)  and `zeon multiplication' in ((Subscript[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∈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 ∈ (Subscript[ℕ, 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⋃B))⋂(A⋂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⋂b]>0,0,1];

(*Unprotect[CircleDot];(* Clifford product *)*)
SetAttributes[CircleDot,{Flat, OneIdentity, Listable}];
CircleDot::usage="CircleDot[u,v] - returns product u ⊙ v in Clifford algebra or subalgebra determined by canonical basis elements.";
CircleDot[x_, y_]:=Expand[x y] /;((FreeQ[x,(Subscript[ζ, _])]∧FreeQ[x,(Subscript[e, _])]∧FreeQ[x,(Subscript[ϛ, _])]∧FreeQ[x,(Subscript[ε, _])])∨(FreeQ[y,(Subscript[ζ, _])]∧FreeQ[y,(Subscript[e, _])]∧FreeQ[y,(Subscript[ϛ, _])]∧FreeQ[y,(Subscript[ε, _])]));
CircleDot[x_, y_]:=0/;{x==0∨y==0};

CircleDot[x_. (Subscript[e, a_]),y_. (Subscript[e, b_])]:=Expand[CircleDot[x,y] productSignatureMap[a,b]](Subscript[e, setSymmetricDifference[a,b]])/.{(Subscript[e, {}])->1};
CircleDot[x_. (Subscript[ζ, a_]),y_. (Subscript[ζ, b_])]:=zeonProductSignatureMap[a,b]Expand[CircleDot[x,y] (Subscript[ζ, a⋃b])];
CircleDot[x_. (Subscript[ϛ, a_]),y_. (Subscript[ϛ, b_])]:=Expand[CircleDot[x,y] symProductSignatureMap[a,b]](Subscript[ϛ, setSymmetricDifference[a,b]])/.{(Subscript[ϛ, {}])->1};
CircleDot[x_. (Subscript[ε, a_]),y_. (Subscript[ε, b_])]:=Expand[CircleDot[x,y]](Subscript[ε, a⋃b])/.{(Subscript[ε, {}])->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,_⊙_] && arg2  !=  arg1,
arg2 =  clExpand[arg2];
];
arg3=Distribute[ExpandAll[arg2],Plus,wedge];
If[!FreeQ[arg3,_⋀_] && arg3  !=  arg2,
arg3=  clExpand[arg3];
];
arg2=arg3;
Return[If[!FreeQ[e][arg2],Collect[arg2,(Subscript[e, _])],If[!FreeQ[ζ][arg2],Collect[arg2,(Subscript[ζ, _])],If[!FreeQ[ϛ][arg2],Collect[arg2,(Subscript[ϛ, _])],If[!FreeQ[ε][arg2],Collect[arg2,(Subscript[ε, _])],arg2]]]
]];
(*Protect[clExpand];*)];
(*Protect[setSignature];*)

(* Define exterior (wedge) product *)
SetAttributes[wedge, {Flat, OneIdentity, Listable}];
wedge::usage="wedge[u,v] - returns exterior (wedge) product u⋀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⋂b]≠0,0,Expand[antisymmetricProductSignatureMap[a,b]x y ] (Subscript[e, (a⋃b)])];
wedge[x_, y_]:=x y/;(FreeQ[x,(Subscript[e, _])]∨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 ∈ (Subscript[Cl, p,q,r]) or u ∈ (SuperStar[(SubStar[Cl])]) *)
clGrade::usage="clGrade[u] - Maximal clGrade of nonzero terms in canonical expansion of u.";
clGrade[B_]:=Module[{m,X},
If[(FreeQ[B,e]∧FreeQ[B,ζ]∧FreeQ[B,ϛ]∧FreeQ[B,ε]),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⌋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])⊙(Subscript[e, b])],(clGrade[(Subscript[e, b])]-clGrade[(Subscript[e, a])])]],0];

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

rightContract::usage="rightContract[u,v] - Clifford right contraction u⌊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])⊙(Subscript[e, b])],(clGrade[(Subscript[e, a])]-clGrade[(Subscript[e, b])])]],0];

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


(* Compute powers of Clifford elements *)
cliffordPower::usage="cliffordPower[u,m] - Compute \!\(\*SuperscriptBox[\(u\), \(\(\\ \)\(⊙m\)\)]\) for nonnegative integer m and u ∈ (Subscript[Cl, p, q, r]) or other (SuperStar[(SubStar[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 ⊙ clExpand[CircleDot[cliffordPower[y,((n-1)/2)],cliffordPower[y,((n-1)/2)]]]]]]];];

(* Compute wedge powers of Clifford elements *)
wedgePower::usage="wedgePower[u,m] - Compute (u^(⋀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⋀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 (Subscript[Cl, p, q, r]), ((Subscript[Cl, n])^nil), ((Subscript[Cl, n])^idem), or ((Subscript[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 (Subscript[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>=(∑_I)(Subscript[u, I])(Subscript[v, I]) inner product of Clifford elements u,v ∈ (Subscript[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>=(∑_I)(Subscript[u, I])(Subscript[v, I]) inner product <u,v> of zeon elements u,v ∈ ((Subscript[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>=(∑_I)(Subscript[u, I])(Subscript[v, I]) inner product <u,v> of elements u,v ∈ ((Subscript[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>=(∑_I)(Subscript[u, I])(Subscript[v, I]) inner product <u,v> of elements u,v ∈ ((Subscript[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(Subscript[\[RightDoubleBracketingBar], m]) m-norm of u ∈ (Subscript[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(Subscript[\[RightDoubleBracketingBar], m]) m-norm of u ∈ ((Subscript[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(Subscript[\[RightDoubleBracketingBar], m]) m-norm of u ∈ ((Subscript[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\!(Subscript[\[RightDoubleBracketingBar], m]) of u ∈ \!((Subscript[Cl, n])^nil).";
zeonNorm[u_,p_]:=Norm[Prepend[Coefficient[u,zeonBasis[[2;;]]],zRe[u]],p];

zeonSpectralSeminorm::usage="zeonSpectralSeminorm[u] - returns |ℜu| for u∈((Subscript[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 ∈ (Subscript[Cl, p, q, r]) or u ∈ (SuperStar[(SubStar[Cl])]).";
scalarSum[u_]:=Expand[u]/.{(Subscript[e, _])->1,(Subscript[ζ, _])->1,(Subscript[ε, _])->1,(Subscript[ϛ, _])->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 (A^(⊙m)) for square matrix A having entries from (Subscript[Cl, p, q, r]) or (SuperStar[(SubStar[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 (A^(⋀m)) for square matrix A having entries from (Subscript[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(Subscript[>, k]) of u ∈ (Subscript[Cl, p, q, r]) or (SuperStar[(SubStar[Cl])]).";
gradeKPart[x_,k_Integer]:=
DeleteCases[DeleteCases[If[k==0,clExpand[x]/.{(Subscript[e, _])->0,(Subscript[ζ, _])->0,(Subscript[ϛ, _])->0,(Subscript[ε, _])->0},clExpand[x-(x/.{(Subscript[e, Table[_,{k}]])->0,(Subscript[ζ, Table[_,{k}]])->0,(Subscript[ϛ, Table[_,{k}]])->0,(Subscript[ε, 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 ∈ (Subscript[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 ∈ (Subscript[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 ∈ (Subscript[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 ∈ (Subscript[𝒞ℓ, p,q,r])";
OverTilde[x_]:=cliffordReversion[x];

OverHat::usage="OverHat[u] - grade involution of u ∈ (Subscript[Cl, p, q, r])";
OverHat[x_]:=gradeInvolution[x];

OverBar::usage="OverBar[u] - Clifford conjugate of u ∈ (Subscript[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∈(2^([p + q + r]))} ⊆ (Subscript[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}]/.{(Subscript[e, {}])->1}];

zeonBasis::usage="zeonBasis - Returns ordered basis blades {\!\(\*SubscriptBox[\(ζ\), \(\(I\)\(\\ \)\)]\): I∈(2^([n]))} ⊆ ((Subscript[Cl, n])^nil).";
zeonBasis:=clBasis/.{e->ζ};

symBasis::usage="symBasis - Returns ordered basis blades {\!\(\*SubscriptBox[\(ϛ\), \(\(I\)\(\\ \)\)]\): I∈(2^([n]))} ⊆ ((Subscript[Cl, n])^nil).";
symBasis:=clBasis/.{e->ϛ};

idemBasis::usage="idemBasis - Returns ordered basis blades {\!\(\*SubscriptBox[\(ε\), \(\(I\)\(\\ \)\)]\): I∈(2^([n]))} ⊆ ((Subscript[Cl, n])^nil).";
idemBasis:=clBasis/.{e->ε};

(* 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): (Subscript[Cl, p, q, r]);  ζ: ((Subscript[Cl, n])^nil); ε: ((Subscript[Cl, n])^idem); ϛ: ((Subscript[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[{(Subscript[u, 1]), (Subscript[u, 2]), ..., (Subscript[u, k])}] returns the ordered product (Subscript[u, 1]) ⊙ (Subscript[u, 2]) ⊙ ... ⊙ (Subscript[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[{(Subscript[u, 1]), (Subscript[u, 2]), ..., (Subscript[u, k])}] computes wedge product (Subscript[u, 1]) ⋀ (Subscript[u, 2]) ⋀ ... ⋀ (Subscript[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/.{ζ->e,ε->e,ϛ->e};
Return[Xor[AllTrue[clGrade[DeleteCases[Variables[v],Except[(Subscript[e, _])]]],EvenQ],((v/.{(Subscript[e, _])->0})==0∧NoneTrue[clGrade[DeleteCases[Variables[v],Except[(Subscript[e, _])]]],EvenQ])]];];

(* Works in definite signatures.  Decomposition by geometric considerations. *)
versorFactor::usage="versorFactor[u] - Factor versor u ∈ ((Subscript[Cl, n])^nil) (Euclidean Clifford algebra) via geometric considerations.  
Returns vectors {(Subscript[v, i]):1≤i≤clGrade[u]} and scalar α such that u= α(∏_(i = 1)^(clGrade[u]))(Subscript[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[(Subscript[e, {bi[[1Random[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 α and vectors {b_i: 1≤i≤clGrade(B)} such that B=α(∏_(i = 1)^(clGrade (B)))(Subscript[b, i]), provided B ∈ (Subscript[Cl, p, q, r]) is a blade.";
fastBladeFactor[B_]:=Module[{ell,α,gb,F,Fi,Bprime,btemp,u},
α=B[[1]]/.{(Subscript[e, _])->1};
ell=1;
Bprime=Expand[B/α];
F=Expand[B[[1]]/α];
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}], α]]];];

(* 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> = (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]]⊙B[[ell]]]+clExpand[B[[ell]]⊙L[[lp]]])/(2clExpand[B[[ell]]⊙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]∧FreeQ[clExpand[CircleDot[u,cliffordReversion[u]]],e]∧signedSeminormSquare[gradeKPart[u,clGrade[u]]]≠0;

(* Decompose elements in arbitrary signatures *)
cliffordDecomp::usage="cliffordDecomp[u] - Returns scalar vectors {(Subscript[v, i]) : 1 ≤ i ≤clGrade[u]} such that u= (∏_(i = 1)^(clGrade[u]))(Subscript[v, i]).
Note: Assumes u is decomposable a priori. If in doubt, apply cliffordDecomposableQ first.";
cliffordDecomp[B_]:=Module[{X,bi,gb,Bprime,btemp,facts},
Bprime=B;
facts={};
While[clGrade[Bprime]>1,
gb=clGrade[Bprime];
bi=getIndex[Bprime];
x=(Subscript[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, _])]∧FreeQ[x,(Subscript[ϛ, _])]∧FreeQ[x,(Subscript[ζ, _])]∧FreeQ[x,(Subscript[ε, _])]);
getIndex[_ (Subscript[e, a_])]:=a;
getIndex[(Subscript[e, a_])]:=a;
getIndex[_ (Subscript[ζ, a_])]:=a;
getIndex[(Subscript[ζ, a_])]:=a;
getIndex[_ (Subscript[ϛ, a_])]:=a;
getIndex[(Subscript[ϛ, a_])]:=a;
getIndex[_ (Subscript[ε, a_])]:=a;
getIndex[(Subscript[ε, a_])]:=a;
getIndex[x_+y_]:=Block[{$RecursionLimit=recLimit},Union[getIndex[x],getIndex[y]]];

(* Signed seminorm square <(OverTilde[uu(Subscript[>, 0])]) *)
signedSeminormSquare::usage="signedSeminormSquare[u] - returns <u(OverTilde[u])(Subscript[>, 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 ⊙ cliffordReversion[u]]-signedSeminormSquare[u]];

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


(* grade involution of inverse: (OverBar[u])/(u(OverTilde[u)]). *)
conjugatePseudoInverse::usage="conjugatePseudoInverse[u] - returns grade involution of (u^(-1)), i.e.,  ((OverBar[u])/(u (OverTilde[u])))   (u ∈ (Subscript[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, ℜu (`body'), of u ∈ ((Subscript[Cl, n])^nil)";
zRe[u_]:=gradeKPart[u,0];

zDu::usage="zDu[u] - Returns dual part, 𝔇u (`soul'), of u ∈ ((Subscript[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: = (Sqrt[(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 ∈ ((Subscript[Cl, n])^nil), provided Re[u]≠0.";
zeonInverse[u_]:=Module[{a,b,β,inv},
a=gradeKPart[u,0];
If[PossibleZeroQ[a],Message[zeonInverse::noninv,u];,
b=Chop[u-a];
β=Expand[-(a^(-2))b];
inv=(a^(-1));
While[clGrade[β]>0,
inv=inv+β;
β=clExpand[(-(a^(-1))b)⊙β];
];
Return[inv];];
];
zeonInverse::noninv="ℜ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 ∈ ((Subscript[Cl, n])^nil).";
zeonMaxGenerator[u_+v_]:=Max[{zeonMaxGenerator[u],zeonMaxGenerator[v]}];
zeonMaxGenerator[(Subscript[ζ, a_])]:=Max[a];
zeonMaxGenerator[x_]:=0/;(FreeQ[x,(Subscript[ζ, _])]);
zeonMaxGenerator[x_ (Subscript[ζ, a_])]:=Max[a]/;(FreeQ[x,(Subscript[ζ, _])]);

(* Zeon left contraction operator *)
zeonLeftContraction::usage="zeonLeftContraction[u,v] - Zeon contraction u⌉v defined by linear extension of (Subscript[ζ, I])⌉(Subscript[ζ, J]) = (Subscript[ζ, J∖I]) if I ⊆ J, 0 otherwise.";
zeonLeftContraction[x_ (Subscript[ζ, a_]), y_ (Subscript[ζ, J_])]:=Expand[x y zeonLeftContraction[(Subscript[ζ, a]),(Subscript[ζ, J])]]/;(FreeQ[x,(Subscript[ζ, _])]∧FreeQ[y,(Subscript[ζ, _])]);
zeonLeftContraction[x_, (Subscript[ζ, J_])]:=0/;(FreeQ[x,(Subscript[ζ, _])]);
zeonLeftContraction[(Subscript[ζ, J_]), y_]:=0/;(FreeQ[y,(Subscript[ζ, _])]);
zeonLeftContraction[(Subscript[ζ, a_]), y_ (Subscript[ζ, J_])]:=y zeonLeftContraction[(Subscript[ζ, a]),(Subscript[ζ, J])]/;(FreeQ[y,(Subscript[ζ, _])]);
zeonLeftContraction[x_ (Subscript[ζ, J_]), (Subscript[ζ, a_])]:=x zeonLeftContraction[(Subscript[ζ, J]),(Subscript[ζ, a])]/;(FreeQ[x,(Subscript[ζ, _])]);
zeonLeftContraction[(Subscript[ζ, a_]),(Subscript[ζ, J_])]:=If[a==J,1,If[SubsetQ[J,a],(Subscript[ζ, 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⌈v defined by linear extension of (Subscript[ζ, I])⌈(Subscript[ζ, J]) = (Subscript[ζ, I∖J]) if J ⊆ I, 0 otherwise.";
zeonRightContraction[x_ (Subscript[ζ, a_]), y_ (Subscript[ζ, J_])]:=Expand[x y zeonRightContraction[(Subscript[ζ, a]),(Subscript[ζ, J])]]/;(FreeQ[x,(Subscript[ζ, _])]∧FreeQ[y,(Subscript[ζ, _])]);
zeonRightContraction[x_, (Subscript[ζ, J_])]:=0/;(FreeQ[x,(Subscript[ζ, _])]);
zeonRightContraction[(Subscript[ζ, J_]), y_]:=0/;(FreeQ[y,(Subscript[ζ, _])]);
zeonRightContraction[(Subscript[ζ, a_]), y_ (Subscript[ζ, J_])]:=y zeonRightContraction[(Subscript[ζ, a]),(Subscript[ζ, J])]/;(FreeQ[y,(Subscript[ζ, _])]);
zeonRightContraction[x_ (Subscript[ζ, J_]), (Subscript[ζ, a_])]:=x zeonRightContraction[(Subscript[ζ, J]),(Subscript[ζ, a])]/;(FreeQ[x,(Subscript[ζ, _])]);
zeonRightContraction[(Subscript[ζ, a_]),(Subscript[ζ, J_])]:=If[a==J,1,If[SubsetQ[a,J],(Subscript[ζ, 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 ω, ψ such that u = ω + (Subscript[ψζ, {m}]) *)
zeonMaxGenDecomp::usage="zeonMaxGenDecomp[u] - Given zeon u with maximum generator index m, return {ω,ψ} such that u = ω + (Subscript[ψζ, {m}]).";
zeonMaxGenDecomp[u_]:=Module[{ψ,ω,m},
m=zeonMaxGenerator[u];
ψ=zeonRightContraction[u,(Subscript[ζ, {m}])];
ω=u-clExpand[(Subscript[ζ, {m}])⊙ψ];
Return[{ω,ψ}]];

(*Zeon Square Root *)
zeonSqRoot::usage="zeonSqRoot[u] Recursive computation of principal square root of u ∈ ((Subscript[Cl, n])^nil), provided Re[u] ≥ 0.";
zeonSqRoot[u_]:=If[zeonMaxGenerator[u]==0,(u^(1/2)),clExpand[zeonSqRoot[zeonMaxGenDecomp[u][[1]]]+clExpand[(Subscript[ζ, {zeonMaxGenerator[u]}])⊙Expand[(1/2)zeonInverse[zeonSqRoot[zeonMaxGenDecomp[u][[1]]]]]⊙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 ∈ ((Subscript[Cl, n])^nil), provided Re[u] ≥ 0.";
zeonKthRoot[u_,k_]:=If[zeonMaxGenerator[u]==0,(u^(1/k)),ExpandAll[clExpand[zeonKthRoot[zeonMaxGenDecomp[u][[1]],k]]+clExpand[(Subscript[ζ, {zeonMaxGenerator[u]}])⊙ExpandAll[(1/k)cliffordPower[zeonInverse[zeonKthRoot[zeonMaxGenDecomp[u][[1]],k]],k-1]]⊙zeonMaxGenDecomp[u][[2]]]]];

(* Zeon rth power: r can be any real number *)
zeonPower::usage="zeonPower[u,r] returns (u^r) (rth power of zeon u) for r∈(Reals) and u ∈ ((Subscript[Cl, n])^nil).  zRe[u]≠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 ∈ ((Subscript[Cl, n])^nil).";
zeonExp[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Exp[x],x,z]]];

zeonLog::usage="zeonLog[u] Evaluate logarithm of u ∈ ((Subscript[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 ∈ ((Subscript[Cl, n])^nil).";
zeonSin[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Sin[x],x,z]]];

zeonCos::usage="zeonCos[u] evaluates the cosine function at u ∈ ((Subscript[Cl, n])^nil).";
zeonCos[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Cos[x],x,z]]];

zeonTan::usage="zeonTan[u] evaluates the tangent function at u ∈ ((Subscript[Cl, n])^nil), provided cos(u) is invertible; i.e., Re[cos(u)] ≠ 0.";
zeonTan[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Tan[x],x,z]]];

zeonArcSin::usage="zeonArcSin[u] Evaluates arcsine at u ∈ ((Subscript[Cl, n])^nil), provided -1 ≤ zRe[u] ≤ 1.";
zeonArcSin[z_]:=Module[{x},zeonArcSin[z]=zeonExtensionEval[ArcSin[x],x,z]];

zeonArcCos::usage="zeonArcCos[u] Evaluates arccosine at u ∈ ((Subscript[Cl, n])^nil), provided -1 ≤ zRe[u] ≤ 1.";
zeonArcCos[z_]:=Module[{x},zeonArcCos[z]=zeonExtensionEval[ArcCos[x],x,z]];

zeonArcTan::usage="zeonArcTan[u] Evaluates arctangent at u ∈ ((Subscript[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 ∈ ((Subscript[Cl, n])^nil).";
zeonSinh[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Sinh[x],x,z]]];
zeonCosh::usage="zeonCosh[u] - Hyperbolic cosine evaluated at u ∈ ((Subscript[Cl, n])^nil).";
zeonCosh[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Cosh[x],x,z]]];
zeonTanh::usage="zeonTanh[u] - Hyperbolic tangent evaluated at u ∈ ((Subscript[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 (u^⋆) given by linear extension of ((Subscript[ζ, I])^⋆)=(Subscript[ζ, [n]∖I]) in ((Subscript[Cl, n])^nil).
Satisfies <α ⊙ (SuperStar[β])(Subscript[>, n]) = <α, β> (Subscript[ζ, [n]]).";
SetAttributes[zeonDual,Listable];
zeonDual[u_+v_]:=Block[{$RecursionLimit=recLimit},zeonDual[u]+zeonDual[v]];
zeonDual[x_]:=x (Subscript[ζ, Range[maxIndex]])/;(FreeQ[x,(Subscript[ζ, _])]);
zeonDual[x_ (Subscript[ζ, a_])]:=x zeonDual[(Subscript[ζ, a])]/;(FreeQ[x,(Subscript[ζ, _])]);
zeonDual[(Subscript[ζ, a_])]:=(Subscript[ζ, Complement[Range[maxIndex],a]])/;(Length[a]<maxIndex);
zeonDual[(Subscript[ζ, 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[ζ, _])],0,maxIndex-clGrade[zeonDual[zDu[u]]]];

(*zeonDivisionAlgorithm[u,v] returns q,r ∈ ((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} ⊂ ((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[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[ζ, I])): I∈(2^([n])), I≠∅} such that u = a∏ (1+(Subscript[a, I])(Subscript[ζ, I]))*)
zeonElementaryFactorization::usage="zeonElementaryFactorization[u] returns factors a, {(1+(Subscript[a, I])(Subscript[ζ, I])): I∈(2^([n])), I≠∅} such that u = a∏ (1+(Subscript[a, I])(Subscript[ζ, 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 (u^v), where u,v ∈ ((Subscript[Cl, n])^nil) and Re u > 0.";
zeonPowerZeon[u_,v_]:=clExpand[zeonExp[clExpand[v⊙zeonLog[u]]]];

(* Zeon series representation of suitably differentiable function f evaluated at u *) 
zeonExtensionEval::usage="zeonExtensionEval[f[x],x,u] evaluates the zeon extension ϕ (of suitably differentiable) real function f at u ∈ ((Subscript[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 ϕ[u], the series representation of zeon extension ϕ 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 ϕ of a zeon variable, induced by real function f, evaluated at u *)
zeonDerivative::usage="zeonDerivative[f[x],x,u,k] returns (ϕ^((k)))(u): the kth order derivative of zeon extension ϕ of real function f, evaluated at u ∈ ((Subscript[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 ϕ of a zeon variable, induced by real function f, evaluated at u *)
zeonAntiDerivative::usage="zeonAntiDerivative[f[x],x,u] returns Φ[u], where Φ[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 ϕ of a zeon variable, induced by real function f, evaluated at u *)
generalZeonDerivative::usage="generalZeonDerivative[f[x],x,k] returns (ϕ^((k)))(u): the kth order derivative of zeon extension ϕ of real function f.  If no order is specified,  k=1 is assumed.";
generalZeonDerivative[f_,var_,ord_:1]:=Module[{ℜu,𝔇u,k},
ℜu=Symbol["ℜu"];
𝔇u=Symbol["𝔇u"];
Return[ExpandAll[Sum[((D[f,{var,k+ord}]/.{var->ℜu})/(k!))(((𝔇u))^k),{k,0,maxIndex}]]]];

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

(* Evaluate a quadratic function ϕ[x] of variable x at zeon u *)
zeonQuadraticEval::usage="zeonQuadraticEval[ϕ[x],x,u] evaluates the quadratic function ϕ[x]=((αx)^2)+βx+γ at u, where α,β,γ,u ε ((Subscript[Cl, n])^nil).";
zeonQuadraticEval[ϕ_,x_,u_]:=Module[{α,β,γ},
α=Coefficient[ϕ,(x^2)];
β=Coefficient[ϕ,x];
γ=ϕ/.{x->0};
Return[clExpand[α⊙u⊙u+β⊙u+γ]]];

(* Compute the discriminant of quadratic polynomial f[x] having zeon coefficients *)
zeonDiscriminant::usage="zeonDiscriminant[f[x],x] computes the zeon discriminant of ϕ[x]=((αx)^2)+βx+γ, where α,β,γ ε ((Subscript[Cl, n])^nil).";
zeonDiscriminant[f_,x_]:=Module[{α,β,γ},
α=Coefficient[f,(x^2)];
β=Coefficient[f,x];
γ=f/.{x->0};
Return[Expand[clExpand[β⊙β]-4clExpand[α⊙γ]]];
];

(* Apply the zeon quadratic formula to ϕ[x] *)
zeonQuadraticFormula::usage="zeonQuadraticFormula[ϕ[x],x] applies the zeon quadratic formula to ϕ[x]=((αx)^2)+βx+γ, 
where α,β,γ ε ((Subscript[Cl, n])^nil). When Δ=0, includes arbitrary κ assumed to be nilpotent of index 2.";
zeonQuadraticFormula[ϕ_,x_]:=Module[{α,β,γ,Δ},
α=Coefficient[ϕ,(x^2)];
β=Coefficient[ϕ,x];
γ=ϕ/.{x->0};
Δ=zeonDiscriminant[ϕ,x];
If[zRe[α]==0,Return[{}],
If[zeonNorm[Δ,1]==0,Return[-clExpand[β⊙Expand[zeonInverse[2α]]]+ κ],
If[zRe[Δ]>0,Return[{clExpand[(-β+zeonSqRoot[Δ])⊙zeonInverse[Expand[2α]]],
clExpand[(-β-zeonSqRoot[Δ])⊙zeonInverse[Expand[2α]]]}],
If[zRe[Δ]<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 |ℜA|≠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 κ such that 
 (𝔇u)^κ=0.";
zeonNilpotency[u_]:= 1+Exponent[zeonExp[t zDu[u]],t];

(* Zeon geometric series *)
zeonGeometricSeries::usage="zeonGeometricSeries[a,u,k] - returns partial sum (∑_(j = 0)^k)(au^j) of geometric series (∑_(j = 0)^∞)(au^j).  
If k is omitted and |ℜu|<1, returns limit of geometric series: (a/(1 - ℜu))(∑_(j = 0)^maxIndex)(1-ℜu)^(-j)(𝔇u)^j.";
zeonGeometricSeries[a_,u_,k_:∞]:=If[k<∞, 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 η = (Subscript[σ, z]) + (Subscript[Iσ, 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 kth Kronecker power of Matrix M, i.e., (M^(⊗k)).";
matrixKroneckerPower[M_,pwr_]:=If[pwr==0,{1},Expand[KroneckerProduct[M,matrixKroneckerPower[M,pwr-1]]]];

cliffordGeneratorRep::usage="cliffordGeneratorRep[j] - (2^n)\[Cross](2^n) real matrix representation of generator (Subscript[e, {j}])∈(Subscript[Cl, p, q, r]):  
(Subscript[e, {j}])-> ((Subscript[σ, x])^(⊗(j - 1))) ⊗ (Subscript[σ, z]) ⊗ ((Subscript[σ, 0])^(⊗(n - j))) (1 ≤ j ≤p)
(Subscript[e, {j}])-> ((Subscript[σ, x])^(⊗(j - 1))) ⊗ (Subscript[Iσ, y]) ⊗ ((Subscript[σ, 0])^(⊗(n - j))) (p < j ≤ p+q)
(Subscript[e, {j}])-> ((Subscript[σ, x])^(⊗(j - 1))) ⊗ η ⊗ ((Subscript[σ, 0])^(⊗(n - j))) (p+q < j ≤n)
where η=(Subscript[σ, z])+(Subscript[Iσ, y]), and (Subscript[σ, 0]), (Subscript[σ, y]), and (Subscript[σ, 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 (2^n)\[Cross]\!\(\*SuperscriptBox[\(2\), \(\(n\)\(\\ \)\)]\)(real) matrix representation of u ∈ (Subscript[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 (Sqrt[ℜ ((SuperStar[tr[(SuperDagger[A]) A]]))]) = \!\(\*SqrtBox[\(∑\\ \*SuperscriptBox[SubscriptBox[\(\[LeftDoubleBracketingBar]\*SubscriptBox[\(a\), \(ij\)]\[RightDoubleBracketingBar]\), \(2\)], \(2\)]\)]\).";
zeonFrobeniusNorm[A_]:=Sqrt[Total[Flatten[zeonInnerProduct[A,A]]]];


zeonGeneratorRep::usage="zeonGeneratorRep[j] - (2^n)\[Cross](2^n) matrix representation of zeon generator (Subscript[ζ, {j}])-> ((Subscript[σ, 0])^(⊗(j - 1))) ⊗ η ⊗ ((Subscript[σ, 0])^(⊗(n - j))),
where η=(Subscript[σ, z])+(Subscript[Iσ, y]), and (Subscript[σ, 0]), (Subscript[σ, y]), and (Subscript[σ, z]) are Pauli matrices.";
zeonGeneratorRep[j_]:=KroneckerProduct[KroneckerProduct[IdentityMatrix[(2^(j-1))],etaMatrix],IdentityMatrix[(2^(maxIndex-j))]];

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

generalZeonElement::usage="generalZeonElement[α] - Canonical expansion(∑_(I ∈ (2^([n]))))(Subscript[α, I])(Subscript[ζ, I]), where α is a symbol.";
generalZeonElement[var_]:=Prepend[(clBasis[[2;;]]/.{e->var}),(Subscript[var, 0])].zeonBasis;

generalCliffordElement::usage="generalCliffordElement[α] - Canonical expansion(∑_(I ∈ (2^([n]))))(Subscript[α, I])(Subscript[e, I]), where α is a symbol.";
generalCliffordElement[var_]:=Prepend[(clBasis[[2;;]]/.{e->var}),(Subscript[var, 0])].clBasis;

generalIdemElement::usage="generalIdemElement[α] - Canonical expansion(∑_(I ∈ (2^([n]))))(Subscript[α, I])(Subscript[ε, I]), where α is a symbol.";
generalIdemElement[var_]:=Prepend[(clBasis[[2;;]]/.{e->var}),(Subscript[var, 0])].idemBasis;

generalSymElement::usage="generalSymElement[α] - Canonical expansion(∑_(I ∈ (2^([n]))))(Subscript[α, I])(Subscript[ϛ, I]), where α is a symbol.";
generalSymElement[var_]:=Prepend[(clBasis[[2;;]]/.{e->var}),(Subscript[var, 0])].symBasis;

];

EndPackage[];

Acknowledgements

Development of this package has been ongoing since 2008.  I would like to thank the following for helpful comments, suggestions, and general discussions.

Rafal Ablamowicz (Tennessee Technological University)
Arturas Acus (Vilnius University, Lithuania)
Marco Budinich (University of Trieste, Italy)
Stephen Sangwine (University of Essex, UK)

References

C. Cassiday, G.S. Staples. On representations of semigroups having hypercube-like Cayley graphs, Clifford Analysis, Clifford Algebras and Their Applications, 4 (2015), 111-130. 

L. Dollar, G.S. Staples. Zeon roots, Advances in Applied Clifford Algebras, 27 (2017), 1133-1145. 

L. Dorst, D. Fontijne. Efficient Algorithms for Factorization and Join of Blades, Geometric Algebra Computing, E.Bayro-Corrochano, G. Scheuermann, Eds., Springer, London, 2010, pp. 457--476.

D. Fontijne.  Efficient Implementation of Geometric Algebra, Ph.D. thesis, University of Amsterdam, 2007. 

E. Haake, G.S. Staples. Zeros of zeon polynomials, (2018) Submitted.

T. Lindell, G.S. Staples. Norm inequalities in zeon algebras, (2018) Submitted.

P. Lounesto, E. Latvamaa. Conformal transformations and Clifford algebras, Proc. Amer. Math. Soc., 79 (1980),  533--538.

C. Perwass. Geometric Algebra with Applications in Engineering, Springer-Verlag, Berlin, 2009.

G.S. Staples. Differential calculus of zeon functions,  (2018) Submitted. 

G.S. Staples, D. Wylie. Clifford algebra decompositions of conformal orthogonal group elements, Clifford Analysis, Clifford Algebras and Their Applications, 4 (2015), 223-240. 

G.S. Staples, A. Weygandt. Elementary functions and factorizations of zeons, Advances in Applied Clifford Algebras (2018) 28: 12.