Zeon Essentials

Rev. 20250205

G. Stacey Staples

Department of Mathematics & Statistics
Southern Illinois University Edwardsville

(*
Author: George Stacey Staples
© 2025. This work is openly licensed via CC BY 4.0.

===============
Disclaimer: This software package is provided "AS IS" without warranty of any kind. 
Should anything in the program prove defective, the end-user assumes the cost of any necessary service, repair, or correction.
===============

*)

Unprotect[ζ];
ζ=Symbol["ζ"];
Protect[ζ];

Unprotect["ℨ"];
ℨ=Symbol["ℨ"];
Protect["ℨ"];


Unprotect[etaMatrix,zeonMinGradeProjection,zeonMaxGradeProjection,zeonMaxGrade,zeonCharacteristicPolynomial];
Unprotect[generalZeonElement,generalZeonAntiDerivative,generalZeonDerivative,zeonIndices,zeonGrade];
Unprotect[zeonGradeBasis,zeonGradeInvolution,zeonGradeProjection,matrixKroneckerPower,zeonEigenvalues];
Unprotect[nilpotentZeonQ,zeonParityQ,matrixListProduct,scalarSum,zeonMatrixInverse,zeonPutzer,zeonMatrixExp];
Unprotect[zDu,zSc,zeonAntiDerivative,zeonArcCos,zeonArcSin,zeonArcTan,zeonBasis,zeonCos,nilpotentMonomialQ];
Unprotect[zeonCosh,zeonCubeRoot,zeonDerivative,zeonDiscriminant,zeonDivisionAlgorithm,zeonHodgeDual,zeonElementaryFactorization];
Unprotect[zeonExp,zeonExtensionEval,zeonFrobeniusNorm,zeonGeneratorRep,zeonGeometricSeries,zeonCubicDiscriminant,zeonCubicFormula];
Unprotect[zeonInnerProduct,zeonInverse,zeonKthRoot,zeonLeftContraction,zeonLog,zeonMaxGenDecomp,zeonMaxGenerator];
Unprotect[zeonMinGrade,zeonNilpotency,zeonNorm,zeonPower,zeonPowerZeon,zeonQuadraticFormula];
Unprotect[zeonRep,zeonRightContraction,zeonSin,zeonSinh,zeonSpectralSeminorm,zeonSqRoot,zeonTan,zeonTanh,zeonZero];
Unprotect[nilMonoSqRoots,nilMonoCubeRoots,zInvQ,zeonMatrixStarDual,zListProduct,zeonVectorIP,zeonVectorHIP,nilMatrixPwrSum];
Unprotect[lambdaFactor,laplacianEigenvalue,reduceTLaplacian,laplacianEigenvector,rowSwap,triangularizeLaplacian,nonunitTReduceRow,nonunitRevReduceRow];
Unprotect[genRandGraph,randomZeonPolynomial,lowerreduceLaplacian,quickLaplacianEigenvector];

ClearAll[etaMatrix,zeonMinGradeProjection,zeonMaxGradeProjection,zeonMaxGrade,zeonCharacteristicPolynomial];
ClearAll[generalZeonElement,generalZeonAntiDerivative,generalZeonDerivative,zeonIndices,zeonGrade];
ClearAll[zeonGradeBasis,zeonGradeInvolution,zeonGradeProjection,matrixKroneckerPower,zeonEigenvalues];
ClearAll[nilpotentZeonQ,nilpotentMonomialQ,zeonParityQ,matrixListProduct,scalarSum,zeonMatrixInverse,zeonPutzer,zeonMatrixExp];
ClearAll[zDu,zSc,zeonAntiDerivative,zeonArcCos,zeonArcSin,zeonArcTan,zeonBasis,zeonCos];
ClearAll[zeonCosh,zeonDerivative,zeonDiscriminant,zeonDivisionAlgorithm,zeonHodgeDual,zeonElementaryFactorization];
ClearAll[zeonExp,zeonExtensionEval,zeonFrobeniusNorm,zeonGeneratorRep,zeonGeometricSeries,zeonCubicDiscriminant,zeonCubicFormula];
ClearAll[zeonInnerProduct,zeonInverse,zeonKthRoot,zeonLeftContraction,zeonLog,zeonMaxGenDecomp,zeonMaxGenerator];
ClearAll[zeonMinGrade,zeonNilpotency,zeonNorm,zeonPower,zeonPowerZeon,zeonQuadraticFormula];
ClearAll[zeonRep,zeonRightContraction,zeonSin,zeonSinh,zeonSpectralSeminorm,zeonSqRoot,zeonTan,zeonTanh,zeonZero];
ClearAll[nilMonoSqRoots,nilMonoCubeRoots,zInvQ,zeonMatrixStarDual,zListProduct,zeonVectorIP,zeonVectorHIP,nilMatrixPwrSum];
ClearAll[lambdaFactor,laplacianEigenvalue,reduceTLaplacian,laplacianEigenvector,rowSwap,triangularizeLaplacian,nonunitTReduceRow,nonunitRevReduceRow];
ClearAll[genRandGraph,randomZeonPolynomial,lowerreduceLaplacian,quickLaplacianEigenvector];

zeonKeywords::usage="zeonKeywords[] - Returns list of zeons package reserved symbols.";
zeonKeywords[]:=Print[MatrixForm[{
"etaMatrix",
generalZeonAntiDerivative,
generalZeonDerivative,
generalZeonElement,
genRandGraph,
laplacianEigenvector,
lowerreduceLaplacian
matrixKroneckerPower,
matrixListProduct,
nilMatrixPwrSum,
nilMonoSqRoots,
nilpotentMonomialQ,
nilpotentZeonQ,
nonunitRevReduceRow,
nonunitTReduceRow,
quickLaplacianEigenvector,
randomZeonPolynomial,
reduceTLaplacian
scalarSum,
triangularizeLaplacian,
zDu,
zeonAntiDerivative,
zeonArcCos,
zeonArcSin,
zeonArcTan,
zeonBasis,
zeonCharacteristicPolynomial,
zeonCos,
zeonCosh,
zeonCubeRoot,
zeonCubicDiscriminant,
zeonCubicFormula,
zeonDerivative,
zeonDiscriminant,
zeonDivisionAlgorithm,
zeonEigenvalues,
zeonElementaryFactorization,
zeonExp,
zeonExtensionEval,
zeonFrobeniusNorm,
zeonGeneratorRep,
zeonGeometricSeries,
zeonGrade,
zeonGradeBasis,
zeonGradeInvolution,
zeonGradeProjection,
zeonHodgeDual,
zeonIndices,
zeonInnerProduct,
zeonInverse,
zeonKthRoot,
zeonLeftContraction,
zeonLog,
zeonMatrixExp,
zeonMatrixInverse,
zeonMatrixStarDual,
zeonMaxGenDecomp,
zeonMaxGenerator,
zeonMaxGrade,
zeonMaxGradeProjection,
zeonMinGrade,
zeonMinGradeProjection,
zeonNilpotency,
zeonNorm,
zeonParityQ,
zeonPower,
zeonPowerZeon,
zeonPutzer,
zeonQuadraticFormula,
zeonRep,
zeonRightContraction,
zeonSin,
zeonSinh,
zeonSpectralSeminorm,
zeonSqRoot,
zeonTan,
zeonTanh,
zeonVectorHIP,
zeonVectorIP,
zeonZero,
zInvQ,
zListProduct,
zSc}]];

(* Initialize zeon product *)  

Unprotect[Times,Power];
Times[(Subscript[ζ, a_]),(Subscript[ζ, b_])]:=If[IntersectingQ[a,b],0,1](Subscript[ζ, a⋃b]); (* zeon multiplication *)
Power[(Subscript[ζ, a_]),p_]:=0/;p>1; (* zeon powers *)
Protect[Times,Power];

Unprotect[Conjugate];
Conjugate[(Subscript[ζ, a_])]:=(Subscript[ζ, a]);
Protect[Conjugate];

(* Scalar part of a zeon element *)
zSc::usage="zSc[u] returns scalar part, ℭu, of u ∈ (Complexes)(Subscript[ℨ, ∞]).";
zSc[u_]:=u/.{(Subscript[ζ, _])->0}; 

(* Dual part of zeon element *)
zDu::usage="zDu[u] returns dual part, 𝔇u, of u ∈ (Complexes)(Subscript[ℨ, ∞]).";
zDu[u_]:=u-zSc[u]; 

zeonGrade::usage="zeonGrade[u] returns grades of all nonzero terms in canonical expansion of u ∈ (Complexes)(Subscript[ℨ, ∞]).";
zeonGrade[x_]:=0/;(FreeQ[x,ζ]);
zeonGrade[(Subscript[ζ, a_])]:=Length[a];
zeonGrade[x_ (Subscript[ζ, a_])]:=Length[a];
zeonGrade[x_+y_]:=Flatten[{zeonGrade[x],zeonGrade[y]}];
SetAttributes[zeonGrade,Listable];

(* 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 minimal grade of 𝔇u, provided 𝔇u≠0.  Returns 0 otherwise.";
zeonMinGrade[x_]:=Min[DeleteCases[zeonGrade[x],0]];

(* Maximal grade of zeon u.   Returns 0 if u is trivial, otherwise returns minimal grade of the dual part of u. *)
zeonMaxGrade::usage="zeonMaxGrade[u] returns maximal grade of 𝔇u, provided 𝔇u≠0.  Returns 0 otherwise.";
zeonMaxGrade[x_]:=Max[DeleteCases[zeonGrade[x],0]];

(* Compute zeon inverses *)
zeonInverse::usage="zeonInverse[u] - Multiplicative inverse of u ∈ (Complexes)(Subscript[ℨ, ∞]), provided ℭu≠0.";
zeonInverse[u_]:=Module[{a,b,β,inv},
a=u/.{(Subscript[ζ, _])->0};
If[PossibleZeroQ[a],Message[zeonInverse::noninv,u];,
b=Chop[u-a];
β=Expand[-(a^(-2))b];
inv=(a^(-1));
While[zeonMaxGrade[β]>0,
inv=inv+β;
β=Expand[(-(a^(-1))b)β];
];
Return[inv];];
];
zeonInverse::noninv="ℭu=0. The argument `1` is not invertible.";
SetAttributes[zeonInverse,Listable];

(* zeon grade involution *)
zeonGradeInvolution::usage="zeonGradeInvolution[u] acts by linear extension of (Subscript[ζ, I]) |-> (-1)^( | I|)(Subscript[ζ, I]).";
zeonGradeInvolution[x_]:=x/;(FreeQ[x,ζ]);
zeonGradeInvolution[(Subscript[ζ, a_])]:=(((-1))^(Length[a]))(Subscript[ζ, a]);
zeonGradeInvolution[x_ (Subscript[ζ, a_])]:=(((-1))^(Length[a]))x (Subscript[ζ, a])/;(FreeQ[x,ζ]);
zeonGradeInvolution[x_ +y_]:=zeonGradeInvolution[x]+zeonGradeInvolution[y];

(* Determine nilpotency of dual part of x.  If x is a scalar, return ∞ *)
zeonNilpotency::usage="zeonNilpotency[u] - Returns index of nilpotency of 𝔇u, the dual part of u ∈ (Complexes)(Subscript[ℨ, ∞]). Returns ∞ if u is a scalar.";
zeonNilpotency[x_]:=Module[{dx,j},dx=zDu[x];
If[PossibleZeroQ[dx],j=∞,
j=2;
While[!PossibleZeroQ[dx=Expand[dx zDu[x]]],j++];];
Return[j]];

(* Minimal grade projection *)
zeonMinGradeProjection::usage="zeonMinGradeProjection[u] returns the minimal grade part of u ∈ (Complexes)(Subscript[ℨ, ∞]).";
zeonMinGradeProjection[u_]:=zeonGradeProjection[u,zeonMinGrade[u]];

(* Maximal grade projection *)
zeonMaxGradeProjection::usage="zeonMaxGradeProjection[u] returns the maximal grade part of u ∈ (Complexes)(Subscript[ℨ, ∞]).";
zeonMaxGradeProjection[u_]:=zeonGradeProjection[u,zeonMaxGrade[u]];

zeonCharacteristicPolynomial::usage="zeonCharacteristicPolynomial[M,u] returns characteristic polynomial of matrix M in the zeon variable u.";
zeonCharacteristicPolynomial[mat_,var_]:=Collect[Expand[Det[var IdentityMatrix[Dimensions[mat][[1]]]-mat]],var];

(* Return spectrally simple zeon zero with specified scalar part *)
zeonZero::usage="zeonZero[ϕ,u,(Subscript[u, 0])]: Given zeon polynomial ϕ of zeon variable u,  returns spectrally simple zeon zero λ satisfying ℭλ = (Subscript[u, 0]).  If no value is given for (Subscript[u, 0]), then the first zero of complex polynomial ℭϕ is used.";
SetAttributes[zeonZero, Listable];
zeonZero[zpoly_,var_,firstroot_:"x"]:=Module[{λ,g,μ,i,t,κ},
If[firstroot=="x",λ=N[var/.Solve[zSc[zpoly]==0,var][[1]]],λ=firstroot];
g=ExpandAll[PolynomialQuotient[zSc[zpoly],var-λ,var]];
μ=ExpandAll[g/.{var->λ}];
κ=Chop[N[Expand[zpoly/.{var->λ}]]];
t=zeonMinGrade[κ];
While[t≠0,
λ=Chop[N[λ-Expand[κ/μ]]];
κ=zeonMinGradeProjection[Chop[N[ExpandAll[zpoly/.{var->λ}]],(10^(-1))]];
t=zeonMinGrade[κ];
];
Return[λ];
];


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,mxi,proX},
ra=zSc[A];
da=zDu[A];
mxi=Max[Flatten[Variables[A]/.((Subscript[ζ, a_])->a)]];
If[!SquareMatrixQ[A],
Print["Matrix is not square!"];,
If[Det[ra]≠0,
ira=Inverse[ra];
X=Expand[da.ira];
U=ira;
proX=IdentityMatrix[Length[A]];
While[!FreeQ[PossibleZeroQ[proX],False],
proX=ExpandAll[proX.(-X)];
U=U+ira.proX;];
Return[U];,Print["Matrix is not invertible!"];]]];


zeonEigenvalues::usage="zeonEigenvalues[A] returns zeon eigenvalues of square matrix A, provided all eigenvalues of A are spectrally simple.";
zeonEigenvalues[A_]:=Module[{eigs,zeigs,zcp},
If[SquareMatrixQ[A],
eigs=N[Eigenvalues[zSc[A]]];
If[CountDistinct[eigs]==Length[eigs],
zcp=zeonCharacteristicPolynomial[A,u];
zeigs=zeonZero[zcp,u,eigs];
Return[zeigs],
Print["Eigenvalues are not spectrally simple!"];]
,Print["Matrix is not square!"];];
];


zeonInnerProduct::usage="zeonInnerProduct[u,v] - returns <u,v>=(∑_I)(Subscript[u, I])(OverBar[(Subscript[v, I])]) inner product <u,v> of zeon elements u,v ∈ (Complexes)(Subscript[ℨ, ∞])";
SetAttributes[zeonInnerProduct,Listable];
zeonInnerProduct[u_,v_]:=Module[{u0,v0,udotv,mxI},
u0=zSc[u];
v0=zSc[v];
mxI=Max[{zeonMaxGenerator[u],zeonMaxGenerator[v]}];
udotv=Prepend[Coefficient[(u-u0),zeonBasis[mxI][[2;;]]],u0].Conjugate[Prepend[Coefficient[(v-v0),zeonBasis[mxI][[2;;]]],v0]];
Return[udotv];];


zeonNorm::usage="zeonNorm[u,p] - returns the p-norm \[LeftDoubleBracketingBar]u(Subscript[\[RightDoubleBracketingBar], p]) of u ∈ (Complexes)(Subscript[ℨ, ∞]).  (Default is p=2.)";
zeonNorm[u_,p_:2]:=Norm[Prepend[Coefficient[u,zeonBasis[zeonMaxGenerator[u]][[2;;]]],zSc[u]],p];

zeonSpectralSeminorm::usage="zeonSpectralSeminorm[u] - returns |ℭu| for u ∈ (Complexes)(Subscript[ℨ, ∞]).";
zeonSpectralSeminorm[u_]:=Abs[zSc[u]];

zListProduct::usage="zListProduct[{(Subscript[u, 1]), (Subscript[u, 2]), ..., (Subscript[u, k])}] returns the product (Subscript[u, 1]) (Subscript[u, 2]) ... (Subscript[u, k]).";
zListProduct[L_]:=ExpandAll[Product[L[[j]],{j,1,Length[L]}]];

(* Returns sum of scalar coefficients from canonical expansion *)
scalarSum::usage="scalarSum[u] returns sum of scalar coefficients in canonical expansion of u ∈ (Complexes)(Subscript[ℨ, ∞]).";
scalarSum[u_]:=Expand[u]/.{(Subscript[ζ, _])->1};

(* Returns the grade k part of an expression *)
zeonGradeProjection::usage="zeonGradeProjection[u,k] - Returns grade-k part <u(Subscript[>, k]) of u.";
zeonGradeProjection[x_,k_Integer]:=
DeleteCases[DeleteCases[If[k==0,Expand[x]/.{(Subscript[ζ, _])->0},Expand[x-(x/.{(Subscript[ζ, Table[_,{k}]])->0})]],0._],0.`];

(* Return canonical basis for (Complexes)(Subscript[ℨ, n]) *)

zeonBasis::usage="zeonBasis[n] - Returns set of basis blades {(Subscript[ζ, I]): I∈(2^([n]))} for (Complexes)(Subscript[ℨ, n]).";
zeonBasis[n_]:=Module[{i,table,pwrset},
pwrset=Subsets[Range[n]][[2;;]];
table=Prepend[Sort[pwrset,Total[(2^(#1))]<Total[(2^(#2))]&],{}];
Table[(Subscript[ζ, table[[i]]]),{i,2^n}]/.{(Subscript[ζ, {}])->1}];


zeonGradeBasis::usage="zeonGradeBasis[k,n] - Returns basis zeon k-blades in (Complexes)(Subscript[ℨ, n]) for positive integer k.";
zeonGradeBasis[k_,n_]:=Module[{i,table,indices},
indices=Subsets[Range[n],{k}];
table=Sort[indices,Total[(2^(#1))]<Total[(2^(#2))]&];
Table[(Subscript[ζ, table[[i]]]),{i,Length[indices]}]];


(* True if all terms of u are of even grade or all terms are of odd grade *)
zeonParityQ::usage="zeonParityQ[u] - True if all terms of u are of even grade or all terms are of odd grade.";
zeonParityQ[u_]:=Xor[AllTrue[zeonGrade[DeleteCases[Variables[u],Except[(Subscript[ζ, _])]]],EvenQ],((u/.{(Subscript[ζ, _])->0})==0∧NoneTrue[zeonGrade[DeleteCases[Variables[u],Except[(Subscript[ζ, _])]]],EvenQ])];


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

zeonIndices::usage="zeonIndices[u] - Extract all indices appearing among multi-indices of u.";
zeonIndices[x_]:={}/;(FreeQ[x,(Subscript[ζ, _])]);
zeonIndices[_ (Subscript[ζ, a_])]:=a;
zeonIndices[(Subscript[ζ, a_])]:=a;
zeonIndices[x_+y_]:=Union[zeonIndices[x],zeonIndices[y]];


zeonMaxGenerator::usage="zeonMaxGenerator[u] - Maximum generator index appearing in canonical expansion of u ∈ (Complexes)(Subscript[ℨ, ∞]).";
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[ζ, _])]);
SetAttributes[zeonMaxGenerator,Listable];

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

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

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-Expand[(Subscript[ζ, {m}])ψ];
Return[{ω,ψ}]];

zeonSqRoot::usage="zeonSqRoot[u] Recursive computation of principal square root of u ∈ (Complexes)(Subscript[ℨ, n]), provided ℭu ≠ 0.";
zeonSqRoot[u_]:=If[zeonMaxGenerator[u]==0,(u^(1/2)),Expand[zeonSqRoot[zeonMaxGenDecomp[u][[1]]]+Expand[(Subscript[ζ, {zeonMaxGenerator[u]}]) Expand[(1/2)zeonInverse[zeonSqRoot[zeonMaxGenDecomp[u][[1]]]]] zeonMaxGenDecomp[u][[2]]]]];

nilMonoSqRoots::usage="nilMonoSqRoots[u] returns `principal' square roots of nilpotent monomial u.";
nilMonoSqRoots[u_]:=Module[{coeff,aye,parts,roots,s,k},
If[nilpotentMonomialQ[u],
coeff=u/.{(Subscript[ζ, _])->1};
aye=zeonIndices[u];
parts=IntegerPartitions[Length[aye],{2}];
roots={};
For[j=1,j<=Length[parts],j++,
s=DeleteCases[Subsets[aye,parts[[j,1]]],{}];
roots=Join[roots,Table[Sqrt[coeff/2]((Subscript[ζ, s[[k]]])+(Subscript[ζ, Complement[aye,s[[k]]]])),{k,Length[s]}]];
];
Return[DeleteDuplicates[roots]];
,Return["Not a nilpotent monomial."];];];

nilMonoCubeRoots::usage="nilMonoCubeRoots[u] returns `principal' cube roots of nilpotent monomial u, provided zeonGrade[u] is at least 3.";
nilMonoCubeRoots[u_]:=Module[{coeff,aye,parts,roots,s,k},
If[nilpotentMonomialQ[u]∧zeonGrade[u]>=3,
coeff=u/.{(Subscript[ζ, _])->1};
aye=zeonIndices[u];
parts=IntegerPartitions[Length[aye],{3}];
roots={};
For[j=1,j<=Length[parts],j++,
s=DeleteCases[Subsets[aye,parts[[j,1]]],{}];
roots=Flatten[Join[roots,Table[If[RealValuedNumberQ[coeff],((coeff/6)^(1/3)),(((coeff/6))^(1/3))]((Subscript[ζ, s[[k]]])+(Sqrt[2])nilMonoSqRoots[(Subscript[ζ, Complement[aye,s[[k]]]])]),{k,Length[s]}]]];
];
Return[DeleteDuplicates[roots]];
,Return["Not a nilpotent monomial or insufficient grade."];];];

SetAttributes[zeonCubeRoot,Listable];
zeonCubeRoot::usage="zeonCubeRoot[u] Computation of principal cube root of u ∈ (Complexes)(Subscript[ℨ, ∞]), provided ℭu ≠ 0.";
zeonCubeRoot[u_]:=If[zeonMaxGenerator[u]==0,If[RealValuedNumberQ[u],Surd[u,3],Power[u,1/3]],ExpandAll[Expand[zeonCubeRoot[zeonMaxGenDecomp[u][[1]]]]+Expand[(Subscript[ζ, {zeonMaxGenerator[u]}]) ExpandAll[(1/3)Power[zeonInverse[zeonCubeRoot[zeonMaxGenDecomp[u][[1]]]],2]] zeonMaxGenDecomp[u][[2]]]]];
  
(* Zeon kth root - often slower than zeonPower[u,1/k]*)
zeonKthRoot::usage="zeonKthRoot[u,k] Recursive computation of principal kth root of u ∈ (Complexes)(Subscript[ℨ, ∞]), provided ℭu ≠ 0.";
zeonKthRoot[u_,k_]:=If[zeonMaxGenerator[u]==0,If[RealValuedNumberQ[u],Surd[u,k],Power[u,1/k]],ExpandAll[Expand[zeonKthRoot[zeonMaxGenDecomp[u][[1]],k]]+Expand[(Subscript[ζ, {zeonMaxGenerator[u]}]) ExpandAll[(1/k)Power[zeonInverse[zeonKthRoot[zeonMaxGenDecomp[u][[1]],k]],k-1]] zeonMaxGenDecomp[u][[2]]]]];

zeonPower::usage="zeonPower[u,r] returns (u^r) (rth power of zeon u) for r ∈ (Reals) and u ∈ (Complexes)(Subscript[ℨ, ∞]).  ℭu≠0 required if r<0.";
zeonPower[u_,r_]:=Module[{x0,Du,x,mI},
x0=zSc[u];
Du=u-x0;
mI=zeonNilpotency[Du];
Return[ExpandAll[Sum[(1/(k!))(D[(x^r),{x,k}]/.{x->x0})Expand[Power[Du,k]],{k,0,mI}]]]];


zeonExp::usage="zeonExp[u] returns the exponential of u ∈ (Complexes)(Subscript[ℨ, ∞]).";
zeonExp[u_]:=Module[{x0,Du,mxI},
x0=zSc[u];
Du=u-x0;
mxI=zeonNilpotency[Du];
Return[If[mxI<∞,ExpandAll[Sum[((Exp[x0])/(k!))Expand[Power[Du,k]],{k,0,mxI}]],Exp[x0]]]];


zeonLog::usage="zeonLog[u] Evaluate logarithm of invertible u ∈  (Complexes)(Subscript[ℨ, ∞]); i.e., provided ℭu≠0.";
zeonLog[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Log[x],x,z]]];

zeonSin::usage="zeonSin[u] evaluates the sine function at u ∈ (Complexes)(Subscript[ℨ, ∞]).";
zeonSin[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Sin[x],x,z]]];

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

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

zeonArcSin::usage="zeonArcSin[u] Evaluates arcsine at u ∈ (Complexes)(Subscript[ℨ, ∞]).";
zeonArcSin[z_]:=Module[{x},zeonExtensionEval[ArcSin[x],x,z]];

zeonArcCos::usage="zeonArcCos[u] Evaluates arccosine at u ∈ (Complexes)(Subscript[ℨ, ∞]).";
zeonArcCos[z_]:=Module[{x},zeonExtensionEval[ArcCos[x],x,z]];

zeonArcTan::usage="zeonArcTan[u] Evaluates arctangent at u ∈ (Complexes)(Subscript[ℨ, ∞]).";
zeonArcTan[z_]:=Module[{x},zeonExtensionEval[ArcTan[x],x,z]];

zeonSinh::usage="zeonSinh[u] - Hyperbolic sine evaluated at u ∈  (Complexes)(Subscript[ℨ, ∞]).";
zeonSinh[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Sinh[x],x,z]]];
zeonCosh::usage="zeonCosh[u] - Hyperbolic cosine evaluated at u ∈  (Complexes)(Subscript[ℨ, ∞]).";
zeonCosh[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Cosh[x],x,z]]];
zeonTanh::usage="zeonTanh[u] - Hyperbolic tangent evaluated at u ∈ (Complexes)(Subscript[ℨ, ∞]).";
zeonTanh[z_]:=Module[{x},ExpandAll[zeonExtensionEval[Tanh[x],x,z]]];



(*zeonDivisionAlgorithm[u,v] returns q,r ∈ (Complexes)(Subscript[ℨ, ∞]) 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} ⊂ (Complexes)(Subscript[ℨ, ∞]) 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,mxI},
a=zSc[u];
m=zeonMinGrade[u];
mxI=Max[{zeonMaxGenerator[u],zeonMaxGenerator[v]}];
monomials=DeleteCases[zeonGradeProjection[zeonBasis[mxI],m],0];
normalized=ExpandAll[u/a];
factormonomials=Coefficient[zDu[normalized],monomials]monomials;
q=a Expand[zeonInverse[v]zListProduct[1+factormonomials]];
qv=Expand[q v];
If[Variables[Chop[u-qv]]=={},qr={q,0},
q=Expand[(a/2)Expand[Times[zeonInverse[v],zListProduct[1+2factormonomials]]]];
r=Expand[u-Expand[Times[q,v]]];
qr={q,r};];
Return[{q,r}]];

(*ZeonElementaryFactorization[u] returns factors a, {(1+(Subscript[a, I])(Subscript[ζ, I])): 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⊂(PositiveIntegers), I≠∅} such that u = a∏ (1+(Subscript[a, I])(Subscript[ζ, I])).";
zeonElementaryFactorization[u_]:=Module[{a,factormonomials,factors,m,monomials,normalized,mxI,q},
a=zSc[u];
mxI=zeonMaxGenerator[u];
m=zeonMinGrade[u];
factors={a};
monomials=DeleteCases[zeonGradeProjection[zeonBasis[mxI],m],0];
normalized=ExpandAll[u/a];
factormonomials=Coefficient[zDu[normalized],monomials]monomials;
factors=Union[factors,1+factormonomials];
q=Expand[Times[normalized,zListProduct[1-factormonomials]]];
m=zeonMinGrade[q];
While[m>0,
monomials=DeleteCases[zeonGradeProjection[zeonBasis[mxI],m],0];
factormonomials=Coefficient[zDu[q],monomials]monomials;
factors=Union[factors,1+factormonomials];
q=Expand[Times[q,zListProduct[1-factormonomials]]];
m=zeonMinGrade[q];
];
Return[factors]];

(* Zeon powers (u^v), where u,v are zeons and u is invertible. *)
zeonPowerZeon::usage="zeonPowerZeon[u,v] returns zeon to zeon power (u^v), where u,v ∈ (Complexes)(Subscript[ℨ, ∞]) and u is invertible.";
zeonPowerZeon[u_,v_]:=Expand[zeonExp[Expand[v zeonLog[u]]]];

(* Zeon series representation of suitably differentiable function f evaluated at u *) 
zeonExtensionEval::usage="zeonExtensionEval[f[x],x,(Subscript[u, 0])] evaluates the zeon extension ϕ (of suitably differentiable) complex function f at (Subscript[u, 0]) ∈ (Complexes)(Subscript[ℨ, ∞]).";
zeonExtensionEval[f_,var_,u_]:=Module[{x0,Du,mxI},
x0=zSc[u];
Du=u-x0;
mxI=zeonNilpotency[Du];
Return[If[mxI<∞,ExpandAll[Sum[((D[f,{var,k}]/.{var->x0})/(k!))Expand[Power[Du,k]],{k,0,mxI}]],f/.{var->x0}]]];

(* Algebraic derivative of zeon-valued function ϕ of a zeon variable, induced by scalar function f, evaluated at u *)
zeonDerivative::usage="zeonDerivative[f[x],x,(Subscript[u, 0]),k] returns (ϕ^((k)))((Subscript[u, 0])): the kth order derivative of zeon extension ϕ of real or complex function f, evaluated at (Subscript[u, 0]) ∈ (Complexes)(Subscript[ℨ, ∞]).  If no order is specified, k=1 is assumed.";
zeonDerivative[f_,var_,u_,ord_:1]:=Module[{x0,Du,mxI},
x0=zSc[u];
Du=u-x0;
mxI=zeonNilpotency[Du];
Return[If[!PossibleZeroQ[Du],ExpandAll[Sum[((D[f,{var,k+ord}]/.{var->x0})/(k!))Expand[Power[Du,k]],{k,0,mxI}]],D[f,{var,ord}]/.{var->x0}]]];

(* Algebraic antiderivative of zeon-valued function ϕ of a zeon variable, induced by scalar function f, evaluated at u *)
zeonAntiDerivative::usage="zeonAntiDerivative[f[x],x,(Subscript[u, 0])] returns Φ[(Subscript[u, 0])], where Φ[x] is the zeon extension of real or complex function F[x], 
the antiderivative of f[x].  Evaluated at (Subscript[u, 0]) ∈ (Complexes)(Subscript[ℨ, ∞])";
zeonAntiDerivative[f_,var_,u_]:=Module[{x0,Du,mxI,F},
x0=zSc[u];
Du=u-x0;
mxI=zeonNilpotency[Du];
F=Integrate[f,var];
Return[ExpandAll[Sum[((D[F,{var,k}]/.{var->x0})/(k!))Expand[Power[Du,k]],{k,0,mxI}]]]];

(* Algebraic derivative of zeon-valued function ϕ of a zeon variable, induced by complex function f as a function of u ∈ (Complexes)(Subscript[ℨ, n]).*)
generalZeonDerivative::usage="generalZeonDerivative[f[x],x,k,n] returns (ϕ^((k)))(u): the kth order derivative of zeon extension ϕ of complex function f 
as function of u ∈ (Complexes)(Subscript[ℨ, n]).  If no order is specified,  k=1 is assumed.";
generalZeonDerivative[f_,var_,ord_:1,n_]:=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,n}]]]];

(* General antiderivative of zeon-valued function ϕ of a zeon variable, induced by real function f, evaluated at u *)
generalZeonAntiDerivative::usage="generalZeonAntiDerivative[f[x],x,n] returns Φ[u]+k, where Φ[u] is the zeon extension of real function F[x], 
the antiderivative of f[x] and u ∈ (Complexes)(Subscript[ℨ, n]).";
generalZeonAntiDerivative[f_,var_,n_]:=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,n}]]+κ]];


(* 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 α,β,γ ε (Complexes)(Subscript[ℨ, ∞]).";
zeonDiscriminant[f_,x_]:=Module[{α,β,γ},
α=Coefficient[f,(x^2)];
β=Coefficient[f,x];
γ=f/.{x->0};
Return[Expand[Expand[β β]-4Expand[α γ]]];
];

(* Apply the zeon quadratic formula to ϕ[x] *)
zeonQuadraticFormula::usage="zeonQuadraticFormula[ϕ[x],x] applies the zeon quadratic formula to ϕ[x]=((αx)^2)+βx+γ, 
where α,β,γ ε (Complexes)(Subscript[ℨ, ∞]). When Δ=0, includes arbitrary κ assumed to be nilpotent of index 2.";
zeonQuadraticFormula[ϕ_,x_]:=Module[{α,β,γ,Δ},
α=Coefficient[ϕ,(x^2)];
β=Coefficient[ϕ,x];
γ=ϕ/.{x->0};
Δ=zeonDiscriminant[ϕ,x];
If[zSc[α]==0,Return[{}],
If[zeonNorm[Δ,1]==0,Return[-Expand[β Expand[zeonInverse[2α]]]+ κ],
If[zSc[Δ]!=0,Return[{Expand[(-β+zeonSqRoot[Δ]) zeonInverse[Expand[2α]]],
Expand[(-β-zeonSqRoot[Δ]) zeonInverse[Expand[2α]]]}],
Return[Indeterminate]]]]];


(* 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[zSc[u]],True,False];

nilpotentMonomialQ::usage="nilpotentMonomialQ[u] returns True if u is a *nonzero* nilpotent monomial; false otherwise.";
nilpotentMonomialQ[u_]:=nilpotentZeonQ[u]∧(Count[u,ζ,Infinity]==1);

nilMatrixPwrSum::usage="nilMatrixPwrSum[A,d] returns <d|∑ (A^k)|d> where A is a square matrix over (Complexes)(ℨ^°)";
nilMatrixPwrSum[A_,d_]:=Module[{Mx,Mxp},
Mx=A[[d]];
Mxp=Mx;
While[(Mx/.{(Subscript[ζ, _])->1}).(Mx/.{(Subscript[ζ, _])->1})!=0,
Mx=Expand[Mx.A];
Mxp=Mxp+Mx;
];
Return[Mxp[[d]]]
];

(* Code to generate a random polynomial *)
(* deg is polynomial degree, var is the variable, gens is number of zeon generators (default 3) and coeffs is an interval from which coefficients are chosen (default [-3,3])  *)
randomZeonPolynomial::usage="randomZeonPolynomial[n,var,gens,coeffs] returns a degree-n monic zeon polynomial in variable var whose coefficients 
are zeon elements of (Complexes)(Subscript[ℨ, gens]) having integer coefficients randomly selected from the set coeffs={min,max}.  By default, gens=3 and coeffs={-3,3}.";
randomZeonPolynomial[deg_,var_, gens_:3,coeffs_:{-3,3}]:=Module[{roots},
roots=Table[Table[Random[Integer,coeffs],{j,1,(2^gens)}].zeonBasis[gens],{k,deg}];
Return[{Expand[Product[var-roots[[k]],{k,1,deg}]],roots}]];

zeonCubicDiscriminant::usage="zeonCubicDiscriminant[ϕ[x],x] computes the cubic discriminant Δ = 18αβγ-4(α^3)γ+(α^2)(β^2)-4(β^3)-27(γ^2) of the monic zeon cubic polynomial ϕ(x)=(x^3)+((αx)^2)+βx+γ.";
zeonCubicDiscriminant[f_,var_]:=Module[{a0,a1,a2,Δ},
a0=f/.{var->0};
a1=Coefficient[f,var];
a2=Coefficient[f,(var^2)];
Δ=ExpandAll[18 a1 a2 a0 - 4(a2^3)a0+(a2^2)(a1^2)-4(a1^3)-27(a0^2) ];
Return[Δ]];


(* Zeon Cubic Formula *)
zeonCubicFormula::usage="zeonCubicFormula[ϕ[x],x] applies the monic zeon cubic formula to the *monic* zeon polynomial ϕ(x)=(x^3)+((αx)^2)+βx+γ, 
where α,β,γ ∈ \!\(\*TemplateBox[{},\"Complexes\"]\)(Subscript[ℨ, ∞]).  Returns all zeon zeros when both the cubic discriminant and β-3(α^2) are invertible.";
zeonCubicFormula[f_,var_]:=Module[{q,r,δ,δrt,s1,s2,z1,z2,z3,a1,a2,a0,dups},
a0=f/.{var->0};
a1=Coefficient[f,var];
a2=Coefficient[f,(var^2)];
If[3zSc[a1]==((zSc[a2])^2), Return["Indeterminate."],
q=ExpandAll[(1/3)a1-(1/9)(a2^( 2))];
r=ExpandAll[(1/6)(a1 a2-3 a0)-(1/27)(a2^( 3))];
δ=ExpandAll[(q^3)+(r^2)];
If[!nilpotentZeonQ[δ]∨PossibleZeroQ[δ],δrt=zeonSqRoot[δ];
dups=0;,
If[(!nilpotentMonomialQ[δ]∨zeonMinGrade[δ]<2)∧!PossibleZeroQ[δ],
Return["Indeterminate"],
δrt=nilMonoSqRoots[δ];
dups=1;
]];
If[(dups==1∧zSc[r]==0)∨(dups==0∧nilpotentMonomialQ[r]),Return["Indeterminate"],
s1=FullSimplify[zeonCubeRoot[Expand[r+δrt]]];
s2=FullSimplify[zeonCubeRoot[Expand[r-δrt]]];
z1=FullSimplify[(s1+s2)-(a2/3)];
z2=FullSimplify[((-1)/2)(s1+s2)-(1/3)a2+I((Sqrt[3])/2)(s1-s2)];
z3=FullSimplify[((-1)/2)(s1+s2)-(1/3)a2-I((Sqrt[3])/2)(s1-s2)];
Return[If[dups==1,DeleteDuplicates[Flatten[{z1,z2,z3}]],{z1,z2,z3}]]]]];


(* 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)^(κ (𝔇u)))(1-ℭu)^(-j)(𝔇u)^j = (a/(1 - u)).";
zeonGeometricSeries[a_,u_,k_:∞]:=If[k<∞, Expand[a Sum[Expand[zeonPower[u,j]],{j,0,k}]],If[Abs[zSc[u]]<1,
Expand[a/(1-zSc[u])Sum[Expand[(((1-zSc[u]))^(-j))zeonPower[zDu[u],j]],{j,0,zeonNilpotency[u]}]]]];

(* Matrix representations *)

etaMatrix::usage="etaMatrix - Returns 2\[Cross]2 matrix η = (Subscript[σ, x]) + (Subscript[Iσ, y]) used to construct representations of zeon generators.";
etaMatrix=(1/2)(PauliMatrix[1]+I PauliMatrix[2]);

matrixListProduct::usage="matrixListProduct[M,q] - Compute ordered product of matrices in list M containing q matrices.";
matrixListProduct[M_,q_]:=If[q==1,M[[1]],M[[1]].matrixListProduct[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]]]];

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


zeonHodgeDual::usage="zeonHodgeDual[u,n] - Returns Zeon-Hodge dual (u^⋆) given by linear extension of ((Subscript[ζ, I])^⋆)=(Subscript[ζ, [n]∖I]) in (Complexes)(Subscript[ℨ, n]).
Satisfies <α   (SuperStar[β])(Subscript[>, n]) = <α, β> (Subscript[ζ, [n]]). Value must be specified for n.";

SetAttributes[zeonHodgeDual,Listable];
zeonHodgeDual[x_,n_]:=x (Subscript[ζ, Range[n]])/;FreeQ[x,(Subscript[ζ, _])];
zeonHodgeDual[(Subscript[ζ, a_]),n_]:=If[Range[n]!=a,(Subscript[ζ, Complement[Range[n],a]]),1];
zeonHodgeDual[x_ (Subscript[ζ, a_]),n_]:=x zeonHodgeDual[(Subscript[ζ, a]),n]/;FreeQ[x,(Subscript[ζ, _])];
zeonHodgeDual[x_+y_,n_]:=zeonHodgeDual[x,n]+zeonHodgeDual[y,n];


zeonVectorHIP::usage="zeonVectorHIP[u,v,n] - returns ℭ[((SuperDagger[v])u)^( * )] in (Complexes)(Subscript[ℨ, n]).";
zeonVectorHIP[u_,v_,n_]:=zSc[zeonHodgeDual[Expand[zeonHodgeDual[Conjugate[v],n].u],n]];

zeonVectorIP::usage="zeonVectorHP[u,v] - returns (SuperDagger[v])u ∈ (Complexes)(Subscript[ℨ, n]).";
zeonVectorIP[u_,v_]:=Expand[Conjugate[v].u];


zeonGeneratorRep::usage="zeonGeneratorRep[j,n] - (2^n)\[Cross](2^n) matrix representation of zeon generator in (Complexes)(Subscript[ℨ, n]): (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_,n_]:=KroneckerProduct[KroneckerProduct[IdentityMatrix[(2^(j-1))],etaMatrix],IdentityMatrix[(2^(n-j))]];

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

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

zeonPutzer::usage="zeonPutzer[A] returns (via Putzer's algorithm) the matrix exponential (E^A), where A is an m×m ``spectrally simple'' zeon matrix.";
zeonPutzer[A_]:=Module[{n,eigs,zeigs,zcp,vars,X,Xprime,v0,g,λsubs,denoms,nums,rfns,P},
If[SquareMatrixQ[A],
eigs=Eigenvalues[zSc[A]];
If[CountDistinct[eigs]==Length[eigs],
n=Length[A];
zcp=zeonCharacteristicPolynomial[A,u];
zeigs=zeonZero[zcp,u,eigs];
vars=Table[(Subscript[r, j])[x],{j,1,n}];
X=DiagonalMatrix[Table[(Subscript[λ, i]),{i,n}]]+Table[If[i==j+1,1,0],{i,n},{j,n}];
Xprime=Table[(Subscript[r, i])'[x],{i,n}];
v0=vars/.{x->0};
g=Flatten[DSolve[{X.vars==Xprime,v0==UnitVector[n,1]},vars,x]];
λsubs=Table[(Subscript[λ, j])->zeigs[[j]],{j,Length[g-1]}];
denoms=zeonInverse[Expand[Table[Denominator[(Subscript[r, j])[x]/.g],{j,Length[g]}]/.λsubs]];
nums=Expand[Table[Numerator[(Subscript[r, j])[x]/.g],{j,Length[g]}]/.{x->1}/.λsubs/.{Exp[a_]:>zeonExp[a]}];
rfns=Expand[denoms nums];
P[k_]:=If[k==0,IdentityMatrix[n],(A-zeigs[[k]]IdentityMatrix[n]).P[k-1]];
Return[Chop[Sum[Expand[rfns[[i]]P[i-1]],{i,1,n}]]];]];];

zeonMatrixExp::usage="zeonMatrix[A,k] approximates the exponential (E^A) of square zeon matrix A via finite sum (∑_(n = 0)^k)((A^n)/(n!)).  Default value k=10.";
zeonMatrixExp[M_,k_:10]:=Module[{Mx,Mxp,it},
Mx=IdentityMatrix[Length[M]];
Mxp=Mx;
it=1;
While[it<=k,
Mxp=Expand[Mxp.M];
Mx=Mx+(1/(it!))Mxp;
it++;];
Return[Mx];];

reduceTLaplacian::usage="reduceTLaplacian[L] returns (by elementary row operations) matrix L' which is lower triangular with exception of last column. ";
reduceTLaplacian[L_]:=Module[{rows,cols,row,oldrow,newL,col,newA,z},
{rows,cols}=Dimensions[L];
row=rows-1;
col=rows-1;
newL=L;
While[col>1,
newL=nonunitRevReduceRow[newL,col,col];
col--;
];
Return[newL];
];

zInvQ::usage="zInvQ[u] returns True if u is invertible and False otherwise.";
zInvQ[u_]:=zSc[u]!=0;
SetAttributes[zInvQ,Listable];

lambdaFactor::usage="lambdaFactor[(Subscript[ζ, I]),Λ,exp,d,pos] returns zeon coefficient of (Subscript[ζ, I]) in Laplacian eigenvalue associated with vertex of degree d.";
lambdaFactor[(Subscript[ζ, a_]),laplacian_,exp_,δ_,pos_]:=Module[{a0},
a0=DeleteCases[a,pos];
Return[(((-1))^(Length[a])) Product[(((δ-laplacian[[a0[[j]],a0[[j]]]]))^(-1)),{j,1,Length[a0]}]Coefficient[exp,(Subscript[ζ, a])]];
];

nilMatrixPowerSum::usage="nilMatrixPwrSum[A,i,j] returns <i|∑ (A^k)|j> where A is a square matrix over \!\(\*TemplateBox[{},\"Complexes\"]\)(ℨ^°)";
nilMatrixPowerSum[A_,row_,col_]:=Module[{Mx,Mxp},
Mx=A[[row]];
Mxp=Mx;
While[(Mx/.{(Subscript[ζ, _])->1}).(Mx/.{(Subscript[ζ, _])->1})!=0,
Mx=Expand[Mx.A];
Mxp=Mxp+Mx;
];
Return[Mxp[[col]]]
];

lowerreduceLaplacian::usage="reduceLaplacian[L,p] returns (by elementary row operations) matrix L' which is lower triangular with exception of pth column. ";
lowerreduceLaplacian[L_,pos_]:=Module[{rows,cols,row,oldrow,newL,col,newA,z},
{rows,cols}=Dimensions[L];
row=rows;
col=rows;
newL=L;
While[col>=1,
newL=If[col!=pos,nonunitRevReduceRow[newL,col,col],newL];
col--;
];
Return[newL];
];

uppertriangularizeLaplacian::usage="uppertriangularizeLaplacian[Λ,p] returns the upper triangular form of the zeon combinatorial Laplacian Λ.  Assumes invertible elements along the main diagonal with 
exception of pth element.";
uppertriangularizeLaplacian[L_,pos_]:=Module[{rows,cols,row,oldrow,newL,col,newA,z},
{rows,cols}=Dimensions[L];
row=1;
col=1;
newL=L;
While[col<=cols,
newL=If[col!=pos,nonunitTReduceRow[newL,col,col],newL];
col++;
];
Return[newL];
];


laplacianEigenvalue::usage="laplacianEigenvalue[Λ,d] returns unique zeon eigenvalue of the zeon combinatorial Λ of graph G associated with 
unique vertex of degree d in G.";
laplacianEigenvalue[laplacian_,δ_]:=Module[{exp,last,dL,pos,lapvars},
pos=Flatten[Position[zSc[laplacian],δ]][[1]];
exp=nilMatrixPwrSum[-zDu[laplacian],pos];
lapvars=Variables[exp];
dL=Table[lambdaFactor[lapvars[[j]],laplacian,exp,δ,pos],{j,1,Length[lapvars]}].lapvars;
Return[δ+dL]];


laplacianEigenvectorComponent::usage="laplacianEigenvectorComponent[Λ,δ,comp] returns \!\(\*StyleBox[\"comp\",FontSlant->\"Italic\"]\)th component of the unit eigenvector ν associated 
with the unique eigenvalue ω of Λ having scalar part zSc[ω]=δ.";
laplacianEigenvectorComponent[laplacian_,δ_,comp_]:=Module[{exp,last,dL,pos,lapvars},
pos=Flatten[Position[zSc[laplacian],δ]][[1]];
exp=nilMatrixPowerSum[-zDu[laplacian],comp,pos];
lapvars=Variables[exp];
dL=Table[lambdaFactor[lapvars[[j]],laplacian,exp,δ,pos],{j,1,Length[lapvars]}].lapvars;
Return[If[δ==laplacian[[comp,comp]],1,Expand[dL/(δ-laplacian[[comp,comp]])]]]];

quickLaplacianEigenvector::usage="quickLaplacianEigenvector[Λ,δ] returns eigenvector associated with the unique eigenvalue ω of 
zeon combinatorial Laplacian Λ satisfying zSc[ω]=δ.";
quickLaplacianEigenvector[laplacian_,δ_]:=Module[{exp,dmat,dL,pos,lapvars},
dmat=DiagonalMatrix[Table[If[δ==laplacian[[j,j]],0,(((δ-laplacian[[j,j]]))^(-1))],{j,1,Length[laplacian]}]];
dL=Expand[Table[laplacianEigenvectorComponent[laplacian,δ,j],{j,1,Length[laplacian]}]];
Return[dL]];

laplacianEigenvector::usage="laplacianEigenvector[Λ,ω] returns zeon eigevector of zeon Laplacian Λ associated with eigenvalue ω, assuming the vertex associated
with ω is represented by the last row & column of Λ.";
laplacianEigenvector[Λ_,ω_]:=Module[{reducedform,position, matrix,dim,inverses,finalform,j},
dim=Length[Λ];
position=Flatten[Position[zSc[Λ],zSc[ω]]][[1]];
matrix=Expand[ω IdentityMatrix[dim]-Λ];
reducedform=lowerreduceLaplacian[uppertriangularizeLaplacian[matrix,position],position];
inverses=DiagonalMatrix[Table[If[j!=position,zeonInverse[reducedform[[j,j]]],1],{j,1,dim}]];
finalform=Expand[inverses.reducedform];
Return[UnitVector[dim,position]-(Transpose[finalform][[position]]Table[If[j!=position,1,0],{j,1,dim}])]];


rowSwap::usage="rowSwap[A,r1,r2] returns matrix A' obtained from A by swapping rows r1 and r2.";
rowSwap[A_,row1_,row2_]:=Module[{Anew,ord,j},
ord=Table[j,{j,Length[A]}]/.{row1->row2,row2->row1};
Return[A[[ord]]]];

nonunitTReduceRow::usage="nonunitTReduceRow[A,row,col] returns matrix A' obtained from A by Gaussian elimination establishing zeros below (Subscript[A, row, col]), 
assuming (Subscript[A, row, col]) is invertible.";
nonunitTReduceRow[A_,row_,col_]:=Module[{zI,other,ell,oldrow,newA},
newA=A;zI=Expand[zeonInverse[A[[row,col]]]];oldrow=A[[row]];other=Table[j,{j,row+1,Length[A]}];For[ell=1,ell<=Length[other],ell++,oldrow=newA[[other[[ell]]]];
newA[[other[[ell]]]]=Expand[oldrow-Expand[oldrow[[col]]Expand[zI newA[[row]]]]];];Return[newA];];

nonunitRevReduceRow::usage="nonunitRevReduceRow[A,row,col] returns matrix A' obtained from A by Gaussian elimination establishing zeros above (Subscript[A, row, col]), 
assuming (Subscript[A, row, col]) is invertible.";
nonunitRevReduceRow[A_,row_,col_]:=Module[{zI,other,ell,oldrow,newA},
newA=A;zI=Expand[zeonInverse[A[[row,col]]]];oldrow=A[[row]];other=Table[j,{j,1,row-1}];For[ell=Length[other],ell>=1,ell--,oldrow=newA[[other[[ell]]]];
newA[[other[[ell]]]]=Expand[oldrow-Expand[oldrow[[col]]Expand[zI newA[[row]]]]];];Return[newA];];

triangularizeLaplacian::usage="triangularizeLaplacian[Λ] returns the upper triangular form of the zeon combinatorial Laplacian Λ.  Assumes invertible elements along the main diagonal with possible exception of last element.";
triangularizeLaplacian[L_]:=Module[{rows,cols,row,oldrow,newL,col,newA,z},
{rows,cols}=Dimensions[L];
row=1;
col=1;
newL=L;
While[col<cols,
newL=nonunitTReduceRow[newL,col,col];
col++;
];
Return[newL];
];

genRandGraph::usage="genRandGraph[n,p] returns the adjacency matrix of a randomly generated simple graph on n vertices in which each pair of vertices has probability p of being adjacent.";
genRandGraph[n_,p_]:=Module[{i,j,A,rw,co},
A=Table[0,{i,1,n},{j,1,n}];For[rw=1,rw≤n,rw++,
For[co=1,co<rw,co++,
If[RandomReal[]≤p,A[[rw,co]]=1;A[[co,rw]]=1;,Null];]];Return[A];];

hopMinimal::usage="hopMinimal[Ψ,(Subscript[v, 0]),(Subscript[v, t])] returns zeon representation of all hop-minimal paths (Subscript[v, 0])->(Subscript[v, t]) in a simple graph with nilpotent adjacency matrix Ψ.";
hopMinimal[NAM_,init_,term_]:=Module[{x},
x=Expand[(Subscript[ζ, {init}])NAM[[init]]];
While[(x/.{(Subscript[ζ, _])->1}).(x/.{(Subscript[ζ, _])->1})≠0∧(x[[term]]/.{(Subscript[ζ, _])->1})==0,
x=Expand[x.NAM];];
Return[x[[term]]]];

relabelGraph::usage="relabelGraph[A,v] Given adjacency matrix A of finite graph on m vertices, returns adjacency matrix A' of relabeled graph in which vertex v is the last vertex.";
relabelGraph[m_,pos_]:=Module[{p},
p=rowSwap[IdentityMatrix[Length[m]],pos,Length[m]];
Return[p.m.p];];

zeonMatrixStarDual::usage="zeonMatrixStarDual[A] returns conjugate transpose Hodge dual (A^‡)=(A^(†^⋆))of matrix A over (Complexes)(Subscript[ℨ, n]), where n defaults to the maximum generator appearing in terms of A.";
zeonMatrixStarDual[A_,n_:0]:=If[n==0,Conjugate[Transpose[zeonHodgeDual[A,Max[zeonMaxGenerator[A]]]]],Conjugate[Transpose[zeonHodgeDual[A,n]]]];

Protect[etaMatrix,zeonMinGradeProjection,zeonMaxGradeProjection,zeonMaxGrade,zeonCharacteristicPolynomial];
Protect[generalZeonElement,generalZeonAntiDerivative,generalZeonDerivative,zeonIndices,zeonGrade];
Protect[zeonGradeBasis,zeonGradeInvolution,zeonGradeProjection,matrixKroneckerPower,zeonEigenvalues];
Protect[nilpotentZeonQ,zeonParityQ,matrixListProduct,scalarSum,zeonMatrixInverse,zeonPutzer,zeonMatrixExp];
Protect[zDu,zSc,zeonAntiDerivative,zeonArcCos,zeonArcSin,zeonArcTan,zeonBasis,zeonCos,nilpotentMonomialQ];
Protect[zeonCosh,zeonCubeRoot,zeonDerivative,zeonDiscriminant,zeonDivisionAlgorithm,zeonHodgeDual,zeonElementaryFactorization];
Protect[zeonExp,zeonExtensionEval,zeonFrobeniusNorm,zeonGeneratorRep,zeonGeometricSeries,zeonCubicDiscriminant,zeonCubicFormula];
Protect[zeonInnerProduct,zeonInverse,zeonKthRoot,zeonLeftContraction,zeonLog,zeonMaxGenDecomp,zeonMaxGenerator];
Protect[zeonMinGrade,zeonNilpotency,zeonNorm,zeonPower,zeonPowerZeon,zeonQuadraticFormula];
Protect[zeonRep,zeonRightContraction,zeonSin,zeonSinh,zeonSpectralSeminorm,zeonSqRoot,zeonTan,zeonTanh,zeonZero];
Protect[nilMonoSqRoots,nilMonoCubeRoots,zInvQ,zeonMatrixStarDual,zListProduct,zeonVectorIP,zeonVectorHIP,nilMatrixPwrSum];
Protect[lambdaFactor,laplacianEigenvalue,reduceTLaplacian,laplacianEigenvector,rowSwap,triangularizeLaplacian,nonunitTReduceRow,nonunitRevReduceRow];
Protect[genRandGraph,randomZeonPolynomial,lowerreduceLaplacian,quickLaplacianEigenvector];

(* All source code is in the (collapsed) cell above.  Highlight and choose cell-properties-open to see it.*)

(* Allowing execution of the initialization cell when prompted loads all procedures. *)

All available commands...

zeonKeywords[]

Help is available...

?zeonSqRoot

?zeonInverse

Zeon Polynomials

Cubic Polynomials

?randomZeonPolynomial

f=randomZeonPolynomial[3,z,2,{-5,5}];
ϕ[u_]:=Collect[f[[1]]/.{z->u},u];
Print["ϕ[u] = ",ϕ[u]];
Print["zSc[ϕ[u]] = ",zSc[ϕ[u]]];
Print["Roots: ",f[[2]]];
Print["---"];
Print["Δ = ",zeonCubicDiscriminant[ϕ[u],u]//Expand];
Print["α = ",α=Coefficient[ϕ[u],(u^2)]];
Print["β = ",β=Coefficient[ϕ[u],u]];
Print["γ = ",γ=ϕ[0]];
Print["q =  ",q=Expand[(3β-(α^2))/9]];
Print["r = ",r=Expand[(β α-3γ)/6-(α^3)/27]];
Print["-108 δ = ",Expand[-108((q^3)+(r^2))]];
Print["3β-(α^2) = ",Expand[3β-(α^2)]];
Print["Depressed cubic: ϕ[u-α/3] = ",Collect[Expand[(u^3)+3q u-2r],u]];
Print["Solutions : ",zeonCubicFormula[ϕ[u],u]];

Graph Theory

Finding all hop-minimal paths in a randomly generated graph

Generate a graph

vertices=1000;
m=genRandGraph[vertices,.01];
A=Transpose[Table[m[[i]](Subscript[ζ, {i}]),{i,1,vertices}]];

Find all vertex sets (with multiplicity) on which hop-minimal paths from random initial vertex to random terminal vertex exist...

initial=Random[Integer,{1,vertices}];
terminal=Random[Integer,{1,vertices}];
g=AdjacencyGraph[m,VertexSize->Larger,VertexStyle->{initial->Green,terminal->Red},EdgeStyle->{Opacity[.3,RGBColor[.3,.7,.8]]},GraphLayout->"BalloonEmbedding"];
Print[g];
Print["Initial: ",initial];
Print["Terminal: ",terminal];

{time,paths}=Timing[hopMinimal[A,initial,terminal]];
Print["Hop-minimal paths ",initial," -> ",terminal,": "];
Print[paths];
Print["---"];
Print["Computed in ",time," seconds."];