(* ::Package:: *)

(* ::Text:: *)
(*Pensieve Header: The U(I2D) program, extending U(I2D)-PolyTime.nb*)


I2DRules = {
u[lft___,1,rgt___]:>u[lft,rgt],
u[lft___,s_^n_.,s_^m_.,rgt___]:>u[lft,s^(n+m),rgt],
u[lft___,y^n_.,x^m_.,rgt___]:>\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(k = 0\), \(m\)]\(
\*SuperscriptBox[\((\(-n\))\), \(k\)]\ Binomial[m, k]\ u[lft, 
\*SuperscriptBox[\(x\), \(\(-k\) + m\)], 
\*SuperscriptBox[\(y\), \(n\)], rgt]\)\),
u[lft___,py^n_.,px^m_.,rgt___]:>u[lft,px^m,py^n,rgt],u[lft___,x^n_.,px^m_.,rgt___]:>u[lft,px^m,x^n,rgt],
u[lft___,x^n_.,py^m_.,rgt___]:>\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(k = 0\), \(n\)]\(
\*SuperscriptBox[\((\(-m\))\), \(k\)]\ Binomial[n, k]\ u[lft, 
\*SuperscriptBox[\(py\), \(m\)], 
\*SuperscriptBox[\(x\), \(\(-k\) + n\)], rgt]\)\),
u[lft___,y^n_.,px^m_.,rgt___]:>u[lft,px^m,y^n,rgt],
u[lft___,y^n_.,py^m_.,rgt___]:>\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(k = 0\), \(Min[m, n]\)]\(Binomial[n, k]\ Binomial[m, k]\ \(k!\)\ u[lft, 
\*SuperscriptBox[\(px\), \(k\)], 
\*SuperscriptBox[\(py\), \(m - k\)], 
\*SuperscriptBox[\(y\), \(n - k\)], rgt]\)\)
};


If[Head[$DegreeStack] =!= List,$DegreeStack={Infinity}];
$ModDegree=First[$DegreeStack];
SetAttributes[ModDegree,HoldRest];
ModDegree[m_,expr_] := Module[{res},
PrependTo[$DegreeStack,$ModDegree=m];
res=expr;
$DegreeStack=Rest[$DegreeStack];
$ModDegree=First[$DegreeStack];
res
];
Deg[t_T] := Plus @@ Take[Plus @@ t, 2];


(* ::Input:: *)
(*Outer[u[#1,#2]-u[#2,#1]&, {px,py,x,y}, {px,py,x,y}]//. I2DRules// Expand // MatrixForm*)


(* ::Input:: *)
(*u[y^3, x,x,py^2, px^4] //. I2DRules// Expand*)


Format[t_T, StandardForm] := Module[
{k, out},
out = {};
Do[
If[t[[k]] === {0,0,0,0},
AppendTo[out, 1],
(* Else *)
If[t[[k,1]] == 1, AppendTo[out, "\[Xi]"]];
If[t[[k,1]] > 1, AppendTo[out, SuperscriptBox["\[Xi]", t[[k,1]]]]];
If[t[[k,2]] == 1, AppendTo[out, "\[Eta]"]];
If[t[[k,2]] > 1, AppendTo[out, SuperscriptBox["\[Eta]", t[[k,2]]]]];
If[t[[k,3]] == 1, AppendTo[out, "x"]];
If[t[[k,3]] > 1, AppendTo[out, SuperscriptBox["x", t[[k,3]]]]];
If[t[[k,4]] == 1, AppendTo[out, "y"]];
If[t[[k,4]] > 1, AppendTo[out, SuperscriptBox["y", t[[k,4]]]]]
];
If[k<Length[t], AppendTo[out, "\[CircleTimes]"]],
{k, Length[t]}
];
DisplayForm[RowBox[out]]
]


Unprotect[NonCommutativeMultiply];
0**_ = 0;
_**0 = 0;
(c_?NumberQ * a_) ** b_ := Expand[c*(a**b)];
a_ ** (c_?NumberQ * b_) := Expand[c*(a**b)];
a_Plus ** b_ := (# ** b)& /@ a;
a_ ** b_Plus := (a ** #)& /@ b;
T[{m1_, m2_, m3_, m4_}]**T[{n1_, n2_, n3_, n4_}] /; m1+m2+n1+n2 >= $ModDegree := 0;
T[{m1_, m2_, m3_, m4_}]**T[{n1_, n2_, n3_, n4_}] := Expand[
u[px^m1, py^m2, x^m3, y^m4, px^n1, py^n2, x^n3, y^n4] //. I2DRules
] /. ut_u :> T[Exponent[Times @@ ut, #]& /@ {px, py, x, y}];
T[t1__] ** T[t2__] /; Deg[T[t1]]+Deg[T[t2]] >= $ModDegree := 0;
T[t1__] ** T[t2__] := Module[
{cc,tt},
Distribute[
OT @@ MapThread[T[#1]**T[#2]&, {{t1}, {t2}}]
] /. ot_OT :> (
cc=1; tt=T[];
Replace[ot, c_.*t_T :> (cc*=c;  AppendTo[tt, t];), {1}];
cc*Flatten[tt]
)
]


(* ::Input:: *)
(*T[{2,1,3,4}]**T[{1,1,0,1}]*)


(* ::Input:: *)
(*T[{2,1,3,4},{0,0,0,1}]**T[{1,1,0,1}, {0,0,1,0}]*)


TPower[expr_, p_Integer] /; p>0 := NonCommutativeMultiply @@ Table[expr, {p}];
TExp[n_Integer, expr_] := Module[
{total, term, k},
k=0;
total = term = T@@Table[{0,0,0,0}, {n}];
While[term =!= 0,
++k;
total += (term=Expand[term**expr/k])
];
total
];
TExp[expr_] := Module[
{n},
{n}=Cases[expr, t_T :> Length[t], Infinity, 1];
TExp[n, expr]
];


(* ::Input:: *)
(*ModDegree[7, TExp[1, T[{0,1,1,0}]]]*)


PutOn[s_Integer, ind_List, T[{0,0,0,0}...]] := T @@ Table[{0,0,0,0}, {s}];
PutOn[s_Integer, ind_List, t_T] := Module[
{k, indk, j, base, rest},
base = T @@ Table[{0,0,0,0}, {s}];
{{k,j}} = Position[t, p_Integer /; p>0, {2}, 1];
indk=Replace[ind[[k]], i_Integer :> {i}];
rest = t;
--rest[[k,j]];
rest = PutOn[s, ind, rest];
Sum[
ReplacePart[base, {indk[[i]], j} -> 1]**rest,
{i, Length[indk]}
]
];
PutOn[s_Integer, ind_List, expr_] := Expand[expr /. t_T :> PutOn[s, ind, t]];


(* ::Input:: *)
(*PutOn[3, {{1,3}}, T[{0,1,1,0}]]*)


(* ::Input:: *)
(*OrderedPartitions[n_,1] := {{n}};*)
(*OrderedPartitions[n_Integer, p_Integer] /; p>1 := OrderedPartitions[n, p] = Join @@ Table[*)
(*Append[#, k]& /@ OrderedPartitions[n-k, p-1],*)
(*{k,0,n}*)
(*];*)
(*Basis[n_,d_] := Flatten[Outer[*)
(*Function[{p1,p2},*)
(*T @@ MapThread[*)
(*Join,*)
(*{Partition[p1,2], Partition[p2,2]}*)
(*]*)
(*],*)
(*OrderedPartitions[d, 2n],*)
(*OrderedPartitions[d, 2n+1],*)
(*1*)
(*]]*)


(* ::Input:: *)
(*OrderedPartitions[5,3]*)


(* ::Input:: *)
(*Partition[Range[9],2]*)


(* ::Input:: *)
(*Basis[2,2] *)


(* ::Input:: *)
(*Basis[2,5] // Length*)
