(* CliffSymNil10 Add-on for CliffMath10 package *) 
(* Defines "symmetric" Clifford product and commutative "null-square" product ("zeon" algebra) *)
(* G. Stacey Staples *)

BeginPackage["CliffSymNil10`"];

(* Define symmetric product for commutative algebra ((Subscript[Cl, p,n-p])^sym). Generators satisfy the following rules *)
(* (Subscript[e, {j}])⊕(Subscript[e, {k}])=(Subscript[e, {k}])⊕(Subscript[e, {y}])   k ≠ j  *)
(* (Subscript[e, {j}])⊕(Subscript[e, {j}]) = 1 if 1 ≤ j ≤ p *)
(* (Subscript[e, {j}])⊕(Subscript[e, {j}]) = -1 if p+1 ≤ j ≤ n  *)
  
Unprotect[CirclePlus];
ClearAll[CirclePlus];
SetAttributes[CirclePlus,{Flat, OneIdentity, Listable}];
CirclePlus[x_?NumericQ ,y_?NumericQ]:=x y ;                                         
CirclePlus[arg1_,y_?NumericQ]:=y  arg1;                                                 
CirclePlus[x_?NumericQ, arg2_]:=x arg2;
CirclePlus[x_?NumericQ  arg1_,y_?NumericQ]:=x y  arg1;
CirclePlus[x_?NumericQ,y_?NumericQ  arg2_]:=x y  arg2;
CirclePlus[x_?NumericQ  arg1_,y_?NumericQ  arg2_]:= x y  CirclePlus[arg1,arg2];   
CirclePlus[x_?NumericQ  arg1_,  arg2_]:=x CirclePlus[arg1,arg2];
CirclePlus[arg1_,y_?NumericQ  arg2_]:= y  CirclePlus[arg1,arg2];

CirclePlus[(Subscript[e, a_]),(Subscript[e, b_])]:=(((-1))^(Length[Complement[a⋂b,pSet]]))*(Subscript[e, (a⋃b)⋂(Complement[Union[a,b],Intersection[a,b]])])/.{(Subscript[e, {}])->1};

CirclePlus=Symbol["CirclePlus"];
Protect[CirclePlus];
Unprotect[ClSymExpand];
ClearAll[ClSymExpand];
SetAttributes[ClSymExpand,Listable];
ClSymExpand[x_]:= x/;FreeQ[x, _⊕_];
ClSymExpand[x_ + y_]:= Simplify[ClSymExpand[x] + ClSymExpand[y]];
ClSymExpand[x_?NumericQ  arg1_]:=x ClSymExpand[arg1];
ClSymExpand[arg1_]:= Module[{arg2},
arg2 = Distribute[ExpandAll[arg1],Plus, CirclePlus];
If[!FreeQ[arg2,_⊕_] && arg2  =!=  arg1,
arg2 =  ClSymExpand[arg2];
];
Return[arg2];
];
Protect[ClSymExpand];

(* Compute powers of (Cl^sym) elements *)
ClSymPower[x_,n_Integer]:=Module[{y},y=ClSymExpand[x];Switch[EvenQ[n],
True,If[n==0,Return[1],Return[ClSymPower[ClSymExpand[y⊕y],(n/2)]]],
False,If[n==1,Return[y],Return[ClSymExpand[y⊕ClSymPower[ClSymExpand[y⊕y],((n-1)/2)]]]]];];

(* Procedure to multiply matrices with (Cl^sym) entries *)
ClSymMatrixProduct[A_,B_]:=If[Dimensions[A][[2]]≠Dimensions[B][[1]],Abort[],
Table[Total[ClSymExpand[CirclePlus[A[[i]],Transpose[B][[j]]]]],{i,1,Length[A]},{j,1,Dimensions[B][[2]]}]];

(* A procedure for computing powers of matrices with (Cl^sym) entries. In this method, A^m is computed by recursive squaring ((A^2)^2...)A *)
ClSymMatrixPower[A_,m_]:=Module[{y},y=ClSymExpand[A];Switch[EvenQ[m],
True,If[m==0,Return[IdentityMatrix[Length[y]]],Return[ClSymExpand[ClSymMatrixPower[ClSymExpand[ClSymMatrixProduct[y,y]],(m/2)]]]],
False,If[m==1,Return[y],Return[ClSymExpand[ClSymMatrixProduct[ClSymMatrixPower[ClSymExpand[ClSymMatrixProduct[y,y]],((m-1)/2)],y]]]]]];

(* Define product for commutative algebra ((Subscript[Cl, n])^nil).  Generators satisfy the following conditions: *) 
(*(Subscript[e, {j}]) ⊖ (Subscript[e, {k}]) = (Subscript[e, {k}]) ⊖ (Subscript[e, {j}])   k≠j *)
(* (Subscript[e, {j}]) ⊖ (Subscript[e, {j}])=0                *)
  
Unprotect[CircleMinus];
Unprotect[CircleMinus];
ClearAll[CircleMinus];
SetAttributes[CircleMinus,{Flat, OneIdentity, Listable}];
CircleMinus[x_?NumericQ ,y_?NumericQ]:=x y ;                                         (* One 'no arguments' case *)
CircleMinus[arg1_,y_?NumericQ]:=y  arg1;                                                  (* Four 'one argument' cases *)
CircleMinus[x_?NumericQ, arg2_]:=x arg2;
CircleMinus[x_?NumericQ  arg1_,y_?NumericQ]:=x y  arg1;
CircleMinus[x_?NumericQ,y_?NumericQ  arg2_]:=x y  arg2;
CircleMinus[x_?NumericQ  arg1_,y_?NumericQ  arg2_]:= x y  CircleMinus[arg1,arg2];    (* Three 'two arguments' cases *)
CircleMinus[x_?NumericQ  arg1_,  arg2_]:=x CircleMinus[arg1,arg2];
CircleMinus[arg1_,y_?NumericQ  arg2_]:= y  CircleMinus[arg1,arg2];
CircleMinus[(Subscript[e, a_]),(Subscript[e, b_])]:= If[Length[a⋂b]>0,0,(Subscript[e, a⋃b])];
CircleMinus=Symbol["CircleMinus"];
Protect[CircleMinus];
Unprotect[ClNilExpand];
ClearAll[ClNilExpand];
SetAttributes[ClNilExpand,Listable];
ClNilExpand[x_]:= x/;FreeQ[x, _⊖_];
ClNilExpand[x_ + y_]:= Simplify[ClNilExpand[x] + ClNilExpand[y]];
ClNilExpand[x_?NumericQ  arg1_]:=x ClNilExpand[arg1];
ClNilExpand[arg1_]:= Module[{arg2},
arg2 = Distribute[ExpandAll[arg1],Plus, CircleMinus];
If[!FreeQ[arg2,_⊖_] && arg2  =!=  arg1,
arg2 =  ClNilExpand[arg2];
];
Return[arg2];
];
Protect[ClNilExpand];

(* Compute powers of (Cl^nil) elements *)
ClNilPower[x_,n_Integer]:=Module[{y},y=ClNilExpand[x];Switch[EvenQ[n],
True,If[n==0,Return[1],Return[ClNilPower[ClNilExpand[y⊖y],(n/2)]]],
False,If[n==1,Return[y],Return[ClNilExpand[y⊖ClNilPower[ClNilExpand[y⊖y],((n-1)/2)]]]]];];

(* Procedure to multiply matrices with (Cl^nil) entries *)
ClNilMatrixProduct[A_,B_]:=If[Dimensions[A][[2]]≠Dimensions[B][[1]],Abort[];,
Table[Total[ClNilExpand[CircleMinus[A[[i]],Transpose[B][[j]]]]],{i,1,Length[A]},{j,1,Dimensions[B][[2]]}]];

(* A procedure for computing powers of (Cl^(nil ))matrices. In this method, A^m is computed by recursive squaring ((A^2)^2...)A *)
ClNilMatrixPower[A_,m_]:=Module[{y},y=ClNilExpand[A];Switch[EvenQ[m],
True,If[m==0,Return[IdentityMatrix[Length[y]]],Return[ClNilExpand[ClNilMatrixPower[ClNilExpand[ClNilMatrixProduct[y,y]],(m/2)]]]],
False,If[m==1,Return[y],Return[ClNilExpand[ClNilMatrixProduct[ClNilMatrixPower[ClNilExpand[ClNilMatrixProduct[y,y]],((m-1)/2)],y]]]]]];

(* Define idempotent product for commutative algebra ((Subscript[Cl, p+q])^idem). Generators satisfy the following rules *)
(* (Subscript[e, {j}])⋄(Subscript[e, {k}])=(Subscript[e, {k}])⋄(Subscript[e, {y}])   k ≠ j     *)
(* (Subscript[e, {j}])⋄(Subscript[e, {j}]) = (Subscript[e, {j}])   1 ≤ j ≤ p+q *)
(* x⋄y is entered by typing x[ESC]dia[ESC]y *)  

Unprotect[Diamond];
ClearAll[Diamond];
SetAttributes[Diamond,{Flat, OneIdentity, Listable}];
Diamond[x_?NumericQ ,y_?NumericQ]:=x y ;                                         (* One 'no arguments' case *)
Diamond[arg1_,y_?NumericQ]:=y  arg1;                                                  (* Four 'one argument' cases *)
Diamond[x_?NumericQ, arg2_]:=x arg2;
Diamond[x_?NumericQ  arg1_,y_?NumericQ]:=x y  arg1;
Diamond[x_?NumericQ,y_?NumericQ  arg2_]:=x y  arg2;
Diamond[x_?NumericQ  arg1_,y_?NumericQ  arg2_]:= x y  Diamond[arg1,arg2];    (* Three 'two arguments' cases *)
Diamond[x_?NumericQ  arg1_,  arg2_]:=x Diamond[arg1,arg2];
Diamond[arg1_,y_?NumericQ  arg2_]:= y  Diamond[arg1,arg2];
Diamond[(Subscript[e, a_]),(Subscript[e, b_])]:= (Subscript[e, a⋃b])/.{(Subscript[e, {}])->1};
Diamond=Symbol["Diamond"];
Protect[Diamond];
Unprotect[ClIdExpand];
ClearAll[ClIdExpand];
SetAttributes[ClIdExpand,Listable];
ClIdExpand[x_]:= x/;FreeQ[x, _⋄_];
ClIdExpand[x_ + y_]:= Simplify[ClIdExpand[x] + ClIdExpand[y]];
ClIdExpand[x_?NumericQ  arg1_]:=x ClIdExpand[arg1];
ClIdExpand[arg1_]:= Module[{arg2},
arg2 = Distribute[ExpandAll[arg1],Plus, Diamond];
If[!FreeQ[arg2,_⋄_] && arg2  =!=  arg1,
arg2 =  ClIdExpand[arg2];
];
Return[arg2];
];
Protect[ClIdExpand];

(* Compute powers of (Cl^sym) elements *)
ClIdPower[x_,n_Integer]:=Module[{y},y=ClIdExpand[x];Switch[EvenQ[n],
True,If[n==0,Return[1],Return[ClIdPower[ClIdExpand[y⋄y],(n/2)]]],
False,If[n==1,Return[y],Return[ClIdExpand[y⋄ClIdPower[ClIdExpand[y⋄y],((n-1)/2)]]]]];];

(* Procedure to multiply matrices with (Cl^sym) entries *)
ClIdMatrixProduct[A_,B_]:=If[Dimensions[A][[2]]≠Dimensions[B][[1]],Abort[],
Table[Total[ClIdExpand[Diamond[A[[i]],Transpose[B][[j]]]]],{i,1,Length[A]},{j,1,Dimensions[B][[2]]}]];

(* A procedure for computing powers of matrices with (Cl^sym) entries. In this method, A^m is computed by recursive squaring ((A^2)^2...)A *)
ClIdMatrixPower[A_,m_]:=Module[{y},y=ClIdExpand[A];Switch[EvenQ[m],
True,If[m==0,Return[IdentityMatrix[Length[y]]],Return[ClIdExpand[ClIdMatrixPower[ClIdExpand[ClIdMatrixProduct[y,y]],(m/2)]]]],
False,If[m==1,Return[y],Return[ClIdExpand[ClIdMatrixProduct[ClIdMatrixPower[ClIdExpand[ClIdMatrixProduct[y,y]],((m-1)/2)],y]]]]]];

EndPackage[];