Bead[x_] := (
  Deg[x] = 2;
  BeadQ[x^n_.] = True
);

BareStates[td_TD] := Flatten[Outer @@ Prepend[
          BareStates /@ (Head /@ td),
          BS
          ]];

Parity[x_.*S_BS] := Mod[Plus @@ (Parity /@ S), 2];

Deg[S_BS] := Plus @@ (Deg /@ S);
Deg[p_Times] := Plus @@ (Deg /@ List @@ p);
Deg[c_Integer] := 0;
Deg[x_^n_Integer] := n*Deg[x];
Deg[x_Plus] := Union[Deg /@ (List @@ x)] /. {d_Integer} :> d;

SprinkleBeads[_, k_] /; k < 0 := {};
SprinkleBeads[_, 0] = {1};
SprinkleBeads[vars_, k_Integer] /; k > 0 :=
    List @@ Expand[(Plus @@ vars)^k] /. c_Integer*vs_ :> vs;
SprinkleBeads[_, _] = {};

States[td_TD, deg_Integer, opts___] := Module[
      {
        parity = Parity /. {opts} /. Parity -> {0, 1},
        bs, vars
        },
      bs = BareStates[td];
      If[Head[parity] == Integer, parity = {parity}];
      bs = Select[bs, MemberQ[parity, Parity[#]] &];
      vars = Union @@ (List @@@ td);
      Flatten[
        (#*SprinkleBeads[vars, (deg - Deg[#])/2]) & /@ bs
        ]
      ];

d[pd_][sum_Plus] := d[pd] /@ sum;
d[pd_][c_Integer*x_] := Expand[c*d[pd][x]];
d[pd_][b_?BeadQ*x_] := Expand[b*d[pd][x]];
d[td_TD][S_BS] := Sum[
        Times[
          (-1)^Parity[Take[S, k - 1]],
          Distribute[MapAt[d, S, k] /. x[i_] :> td[[k, i]]] 
          ],
        {k, Length[S]}
        ] //. {
        BS[l___, c_Integer*s_, r___] :> c*BS[l, s, r],
        BS[l___, b_?BeadQ*s_, r___] :> b*BS[l, s, r]
        };

TDC[td_TD, 
          deg_Integer, opts___Rule] := TDC[td, deg, opts] =
                 Module[   (* Tagged Doodles Chains *)
      {
        parity = Parity /. {opts} /. Parity -> {0, 1},
        s, l
        },
      l = Length[s = States[td, deg, Parity -> parity]];
      SQ[Thread[Rule[s, Array[e, l]]], l, s]
      ]

TDH[td_TD, deg_Integer, opts___Rule] := TDH[td, deg, 
        opts] = Module[  (* Tagged Doodles Homology *)
      {
        parity = Parity /. {opts} /. Parity -> {0, 1},
        s0, sq0, s1, sq1, sq2
        },
      sq0 = FreeSpan[s0 = States[td, deg - 1 - n, Parity -> 1 - parity]];
      sq1 = FreeSpan[s1 = States[td, deg, Parity -> parity]];
      sq2 = FreeSpan[States[td, deg + 1 + n, Parity -> 1 - parity]];
      CoKernel[
        sq0,
        Kernel[sq1, sq2, (# -> d[td][#]) & /@ s1],
        (# -> d[td][#]) & /@ s0
        ]
      ]

pi[n_, x_, y_] := pi[n, x, y] = Cancel[(x^(n + 1) - y^(n + 1))/(x - y)];
pi[x_, y_] := pi[n, x, y];
v1 := x[1] + x[2] - x[3] - x[4];
v2 := x[1]*x[2] - x[3]*x[4];
g[n_, s_, p_] := g[n, s, p] = s^(n + 1) + (n + 1)Sum[
            (-1)^i*Binomial[n - i, i - 1]s^(n + 1 - 2i)p^i/i,
            {i, 1, (n + 1)/2}
            ];
g[s_, p_] := g[n, s, p];
u1 := Cancel[(g[x[1] + x[2], x[1]*x[2]] - g[x[3] + x[4], x[1]*x[2]])/v1];
u2 := Cancel[(g[x[3] + x[4], x[1]*x[2]] - g[x[3] + x[4], x[3]*x[4]])/v2];
U21  := Cancel[(u1 + x[4]*u2 - pi[x[2], x[3]])/(x[1] - x[4])];
V21  := Cancel[(u1 + x[1]*u2 - pi[x[2], x[3]])/(x[4] - x[1])];

BareStates[LL] = {L0, L1};
Parity[L0] = 0; Deg[L0] = 0;
Parity[L1] = 1; Deg[L1] := 1 - n;
d[L0] := Expand[pi[x[1], x[2]]L1];
d[L1] := Expand[(x[1] - x[2])L0];

BareStates[PP] = {P00, P11, P10, P01};
Parity[P00] = 0; Deg[P00] = 0;
Parity[P11] = 0; Deg[P11] := 2 - 2n;
Parity[P10] = 1; Deg[P10] := 1 - n;
Parity[P01] = 1; Deg[P01] := 1 - n;
d[P00] := Expand[pi[x[1], x[4]]P10 + pi[x[2], x[3]]P01];
d[P11] := Expand[(x[2] - x[3])P10 + (x[4] - x[1])P01];
d[P10] := Expand[(x[1] - x[4])P00 + pi[x[2], x[3]]P11];
d[P01] := Expand[(x[2] - x[3])P00 - pi[x[1], x[4]]P11];

BareStates[QQ] = {Q00, Q11, Q10, Q01};
Parity[Q00] = 0; Deg[Q00] = -1;
Parity[Q11] = 0; Deg[Q11] := 3 - 2n;
Parity[Q10] = 1; Deg[Q10] := -n;
Parity[Q01] = 1; Deg[Q01] := 2 - n;
d[Q00] := Expand[u1 Q10 + u2 Q01];
d[Q11] := Expand[v2 Q10 - v1 Q01];
d[Q10] := Expand[v1 Q00 + u2 Q11];
d[Q01] := Expand[v2 Q00 - u1 Q11];

KRC[pd_PD, etc___] := Module[
      {x, xi},
      x[i_] := (Bead[xi = ToExpression["x" <> ToString[i]]]; xi);
      KRC[KRPD @@ pd /. {
            X[i_, j_, k_, l_] /; j - l == 1 || l - j > 1 :> Xp[x@k, x@j, x@
          i, x@l],
            X[
      i_, j_, k_, l_] /; l - j == 1 || j - l > 1 :> Xm[x@l, x@k, x@j, x@i]
            }, etc
        ]
      ];
KRC[krpd_KRPD] := krpd;

KRC[krpd_KRPD, s_String, etc___] := KRC[krpd, Characters[s] /. {"0" -> 0, "1
            " -> 1}, etc];
KRC[krpd_KRPD, a : {(0 | 1) ...}] := Module[
      {x},
      BareStates[
        TD @@ (Thread[{List @@ krpd, a}] /. {
                {x_Xp, 0} :> (QQ @@ x), {x_Xp, 1} :> (PP @@ x),
                {x_Xm, 0} :> (PP @@ x), {x_Xm, 1} :> (QQ @@ x)
                })
        ]
      ];
KRC[krpd_KRPD, a : {(0 | 1) ...}, deg_Integer, opts___] := Module[
      {np, nm, x},
      np = Count[krpd, _Xp]; nm = Count[krpd, _Xm];
      TDH[
        TD @@ (Thread[{List @@ krpd, a}] /. {
                {x_Xp, 0} :> (QQ @@ x), {x_Xp, 1} :> (PP @@ x),
                {x_Xm, 0} :> (PP @@ x), {x_Xm, 1} :> (QQ @@ x)
                }),
        deg + n(nm - np) - nm + (Plus @@ a), opts
        ]
      ];
KRC[krpd_KRPD, h_Integer, deg_Integer, opts___] := Module[
    {np, nm, verts},
    np = Count[krpd, _Xp]; nm = Count[krpd, _Xm];
    If[nm - h < 0 || np + h < 0, SQ[{}, 0, {}],
      verts = Permutations[Join[
            Table[0, {nm - h}],
            Table[1, {np + h}]
            ]];
      SQDirectSum[
        KRC[krpd, #, deg, opts] & /@ verts,
        v @@@ verts
        ]
      ]
    ]

KRC[krpd_, a_][sum_Plus] := KRC[krpd, a] /@ sum;
KRC[krpd_, a_][c_*S_BS] := Expand[c*KRC[krpd, a][S]]; KRC[krpd_KRPD, a : {(
      0 | 1 | "*") ...}][S_BS] /; (Count[a, "*"] == 1) := 
  With[{p = Position[a, "*"][[1, 1]]},
    Head[krpd[[p]]] /. {
          Xp :> (S[[p]] /. {
                  Q00 -> ReplacePart[
                  S, P00, p] + Expand[V21*ReplacePart[S, P11, p]],
                  Q11 -> Expand[(x[1] - x[3])ReplacePart[S, P11, p]],
                  Q10 -> ReplacePart[S, P10, p] + ReplacePart[S, P01, p],
                  Q01 -> x[3]*ReplacePart[S, P10, p] + x[1]*ReplacePart[S, 
                  P01, p]
                  }),
          Xm :> (S[[p]] /. {
                  P00 -> Expand[(x[4] - x[
                        2])ReplacePart[S, Q00, 
                          p] + U21*ReplacePart[S, Q11, p]],
                  P11 -> ReplacePart[S, Q11, p],
                  
            P10 -> x[4]*ReplacePart[S, Q10, p] - ReplacePart[S, Q01, p],
                
                      P01 -> -x[2]*ReplacePart[S, 
                              Q10, p] + ReplacePart[S, Q01, p]
                  })
          } /. x[i_] :> krpd[[p, i]]
    ];

d[pd_PD] := d[KRC[pd]];
d[krpd_KRPD][a_v*S_BS] := Expand[
      sign = 1;
      Sum[If[a[[i]] == 1, 
          sign *= -1; 0, 
          sign*ReplacePart[a, 1, i]*KRC[krpd, List @@ ReplacePart[
      a, "*", i]][S]
          ],
        {i, Length[a]}]
];

KRH[pd_PD, etc___] := KRH[KRC[pd], etc];
KRH[krpd_KRPD, h_Integer,
                     deg_Integer, opts___Rule] := KRH[krpd, h,
                           deg, opts] = Module[  (* Khovanov Rozansky 
                        Homology *)
      {
        parity = Parity /. {opts} /. Parity -> {0, 1},
        s0, sq0, s1, sq1, sq2
        },
      s0 = First /@ First[sq0 = KRC[krpd, h - 1, deg, opts]];
      s1 = First /@ First[sq1 = 
      KRC[krpd, h, deg, opts]];
      sq2 = KRC[krpd, h + 1, deg, opts];
      CoKernel[
        sq0,
        Kernel[sq1, 
      sq2, (# -> d[krpd][#]) & /@ s1],
        (# -> d[krpd][#]) & /@ s0
        ]
      ]

