(* ::Package:: *)

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


BeginPackage["CliffSymNil08`"];


(* Define symmetric product for commutative algebra Subscript[Cl, p,q]^sym. Generators satisfy the following rules *)
(* Subscript[e, {j}]\[CirclePlus]Subscript[e, {k}]=Subscript[e, {k}]\[CirclePlus]Subscript[e, {y}]   k != j  *)
(* Subscript[e, {j}]\[CirclePlus]Subscript[e, {j}] = 1 if 1 <= j <= p *)
(* Subscript[e, {j}]\[CirclePlus]Subscript[e, {j}] = -1 if p+1 <= j <= p+q  *)
  
Needs["CliffMath08`"];
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\[Intersection]b,pSet]]*Subscript[e, (a\[Union]b)\[Intersection](Complement[Union[a,b],Intersection[a,b]])]/.{
\!\(\*SubscriptBox["e", 
RowBox[{"{", "}"}]]\)->1};

CirclePlus=Symbol["CirclePlus"];
Protect[CirclePlus];
Unprotect[ClSymExpand];
ClearAll[ClSymExpand];
SetAttributes[ClSymExpand,Listable];
ClSymExpand[x_]:= x/;FreeQ[x, _\[CirclePlus]_];
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,_\[CirclePlus]_] && 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\[CirclePlus]y],n/2]]],
False,If[n==1,Return[y],Return[ClSymExpand[y\[CirclePlus]ClSymPower[ClSymExpand[y\[CirclePlus]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}] \[CircleMinus] Subscript[e, {k}] = Subscript[e, {k}] \[CircleMinus] Subscript[e, {j}]   k!=j *)
(* Subscript[e, {j}] \[CircleMinus] 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\[Intersection]b]>0,0,Subscript[e, a\[Union]b]];
CircleMinus=Symbol["CircleMinus"];
Protect[CircleMinus];
Unprotect[ClNilExpand];
ClearAll[ClNilExpand];
SetAttributes[ClNilExpand,Listable];
ClNilExpand[x_]:= x/;FreeQ[x, _\[CircleMinus]_];
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,_\[CircleMinus]_] && 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\[CircleMinus]y],n/2]]],
False,If[n==1,Return[y],Return[ClNilExpand[y\[CircleMinus]ClNilPower[ClNilExpand[y\[CircleMinus]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}]\[Diamond]Subscript[e, {k}]=Subscript[e, {k}]\[Diamond]Subscript[e, {y}]   k != j     *)
(* Subscript[e, {j}]\[Diamond]Subscript[e, {j}] = Subscript[e, {j}]   1 <= j <= p+q *)
(* x\[Diamond]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\[Union]b]/.{
\!\(\*SubscriptBox["e", 
RowBox[{"{", "}"}]]\)->1};
Diamond=Symbol["Diamond"];
Protect[Diamond];
Unprotect[ClIdExpand];
ClearAll[ClIdExpand];
SetAttributes[ClIdExpand,Listable];
ClIdExpand[x_]:= x/;FreeQ[x, _\[Diamond]_];
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,_\[Diamond]_] && 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\[Diamond]y],n/2]]],
False,If[n==1,Return[y],Return[ClIdExpand[y\[Diamond]ClIdPower[ClIdExpand[y\[Diamond]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[];