(* Run cabling.nb before running this program *)

Clear[JW]

JW/:JW[a_,b_]:=Arc[a,b]

JW[pp__]:=(JW@@(z$/@Range[Length[{pp}]]))/.{z$[i_]:>{pp}[[i]]}
JW[pp__z$]:=JW[pp]=Module[{p={pp},b,i,n},
n=Length[{pp}]/2;
Collect[Sum[b=Drop[p,{2 n-i,2 n-i+1}];(((-1))^(i+n))((qI[i]/qI[n])/.qIRules) Expand[((JW@@(b/.z$->z$$))/.z$$->z$) Arc[p[[2 n-i]],p[[2 n-i+1]]]],{i,1,n}],_Arc,Together]]

JW@@Table[a[i],{i,1,16}]//AbsoluteTiming

TwistyJW/:TwistyJW[pp__]:=JW@@RotateLeft[{pp},Length[{pp}]/4]

δ=-(1/(q^2))-(q^2);

QI[0]=0;
QI[1]=1;
QI[n_]:=QI[n]=Simplify[δ QI[n-1]-QI[n-2]]

qIRules={qI[n_]:> QI[n]};

Simplify[qI[4]//.qIRules]
Simplify[qI[5]//.qIRules]
Simplify[qI[6]//.qIRules]

Cable[BoxConnectingComponents[a_,b_,c_,d_],strandReplacementRules_List]:=
Module[{n,ar,br,cr,dr,result},
(*Print["while cabling BoxConnectingComponents, received strandReplacementRules", strandReplacementRules];*)
DebugPrintHeld["Beginning ",Cable[BoxConnectingComponents[a,b,c,d],strandReplacementRules]];
{ar,br}={a,b}/.strandReplacementRules;
{cr,dr}=Reverse/@({c,d}/.strandReplacementRules);
n=Length[strandReplacementRules[[1]][[2]] ];
(* If you're expecting a qI[2n+1] in the denominator, notice that it's equal to 1! *)
result= Collect[(1/4)(Α (((-1))^n)(qI[n+1]/.qIRules) * TwistyJW@@Flatten[{ar,br,cr,dr}]+ JW@@Flatten[{ar,dr}] * JW@@Flatten[{br,cr}]),_Arc,Together];
DebugPrintHeld["Finishing ",Cable[BoxConnectingComponents[a,b,c,d],strandReplacementRules]];
result
]

Components[OPD[F_]]:=
Module[{glom},
glom/:glom[x___,y_,z___]glom[u___,y_,v___]:=Union[glom[x,y,z,u,v]];productToList[Times@@(Cases[productToList[F],_Arc|_Xp|_Xn|_X]/.{Arc[a_,b_]:>glom[a,b],(Xp|Xn|X)[a_,b_,c_,d_]:>glom[a,c]glom[b,d]})]
/.glom->List]

Components[OPD[X[a,b,c,d]]]

Components[OPD[Arc[a,b]Arc[c,d]Arc[e,f]]]

Components[OPD[X[a,b,c,d]X[b,a,d,c]]]

Components[OPD[Knot[3,1]]]

findLabels[α_ X_]/;NumericQ[α/.q->1]:=findLabels[X];
findLabels[X_+Y__]:=findLabels[X];
findLabels[X_]:=Union[Flatten[Cases[X,p:(_Arc|_Xp|_Xn):>List@@p,∞]]]

expandOneStep[OPD[F_Times],rules___]:=Module[{sumFactors,otherFactors,coefficients,remainingFactors,labels,bestPair,boundary,bestBoundary},
DebugPrint["Expanding one step, with ",Length[productToList[F]]," factors"];
sumFactors=Cases[productToList[F],f_Plus];
otherFactors=DeleteCases[productToList[F],f_Plus];
coefficients=Cases[otherFactors,f_/;NumericQ[f/.{rules}]];
remainingFactors=DeleteCases[otherFactors,f_/;NumericQ[f/.{rules}]];
If[Length[sumFactors]≤1,
Return[Expand[Times@@coefficients/.{rules}]OPD[Expand[Times@@(remainingFactors~Join~sumFactors)]]/.{rules}];
];
labels=findLabels/@sumFactors;

(* now find the pair of sumFactors which have the longest common boundary! *)
bestPair={1,2};
boundary[i_,j_]:=Length[Intersection[labels[[i]],labels[[j]]]];
bestBoundary:=boundary@@bestPair;
Do[If[boundary[i,j]>bestBoundary,bestPair={i,j}],{i,1,Length[sumFactors]},{j,i+1,Length[sumFactors]}];

DebugPrint["Found a pair of factors which share ",bestBoundary, " labels."];

Expand[Times@@coefficients/.{rules}]OPD[(Times@@remainingFactors) (Times@@Delete[sumFactors,{{bestPair[[1]]},{bestPair[[2]]}}]) (Expand[sumFactors[[bestPair[[1]]]] sumFactors[[bestPair[[2]]]]]/.{rules})]
]

expandOneStep[X_,rules___]:=Expand[X/.{rules}]

expandOneStep[OPD[X_Plus],rules___]:=expandOneStep[OPD[#],rules]&/@X
expandOneStep[X_Plus,rules___]:=expandOneStep[#,rules]&/@X
expandOneStep[α_ opd_OPD,rules___]:=α expandOneStep[opd,rules]

progressiveExpand[opd_OPD,rules___]:=(DebugPrint["beginning progressiveExpand"];FixedPoint[(expandOneStep[#,rules])&,opd])

(* At the moment, this is only defined for two-component links *)

D2nInvariant[L:OPD[F_],n_]/;(Length[Components[L]]==2):=Module[{components=Components[L], firstComponent,secondComponent,twoComponentXings,twoComponentXing,twoComponentXingReplacement,n1,n2,twoComponentLabels,cablingLabels,tempLink,a,b,cable,writheCorrection},

(* Find a crossing of L that has both components coming together *)
firstComponent=components[[1]];
secondComponent=components[[2]];
twoComponentXings=Select[productToList[F],(Intersection[List@@#,firstComponent]≠{}) ∧ (Intersection[List@@#,secondComponent]≠{})&,1];
If[Length[twoComponentXings]==0,
Print["It looks like you gave me a split presentation of the link, and I don't know how to deal with that."];
Return[$Failed]
];
twoComponentXing=twoComponentXings[[1]];

writheCorrection=((((((-1))^n)(q^((n^2)+2n))/.q-> I Exp[I Pi/(8 n -4)]))^(-Writhe[L]));

(* put a special box next to that crossing *)
twoComponentLabels=Flatten[List@@twoComponentXing];
(* ... we have to be careful to put the box in so that it's oriented correctly! *)
(* ... see p. 1573 of Scott's notes! *)
twoComponentXingReplacement=
Switch[Head[twoComponentXing],
Xn,Xn[twoComponentLabels[[1]],twoComponentLabels[[2]],n1,n2] *BoxConnectingComponents[n2,n1,twoComponentLabels[[3]],twoComponentLabels[[4]]],
Xp,Xp[n1,twoComponentLabels[[2]],twoComponentLabels[[3]],n2] *BoxConnectingComponents[twoComponentLabels[[4]],twoComponentLabels[[1]],n1,n2],
_,Print["I found what I thought was a crossing involving both components, but it doesn't look right: ",twoComponentXing];
Return[$Failed]
];

cablingLabels={twoComponentLabels[[1]]->Table[a[k],{k,1,2n-2}],twoComponentLabels[[2]]->Table[b[k],{k,1,2n-2}]};

tempLink=(L/.twoComponentXing->twoComponentXingReplacement);
DebugPrint["about to cable!"];
cable=Cable[tempLink,cablingLabels];
DebugPrint["finished cabling"];
Together[writheCorrection progressiveExpand[cable/.kauffmanRules,Sequence@@qIRules,δ-> -(q^2)-(q^(-2)),q-> I Exp[I Pi/(8 n -4)]]]
]

Simplify[δ/.{q-> I Exp[I Pi/(8 2 -4)]}]

Simplify[D2nInvariant[OPD[Xp[1,2,3,4]Xp[2,1,4,3]],1]//.qIRules]

The unlink!

Simplify[D2nInvariant[OPD[Xp[2,3,4,1]Xn[4,3,2,1]],2]]

The hopf link :

(W[2]=D2nInvariant[OPD[PD[BR[2,{1,1,1,-1}]]],2])//AbsoluteTiming

Simplify[(W[2]/.(Α->0))-(1/4)QuantumKnotInvariant[A1,Irrep[A1][{2}]][BR[2,{1,1}]][q]/.q->Exp[((2 π I)/12)]]

W[2]/.Α->1

FullSimplify[%]

W[4]=D2nInvariant[OPD[PD[BR[2,{1,1,1,1}]]],2]

Simplify[(W[4]/.(Α->0))-(1/4)QuantumKnotInvariant[A1,Irrep[A1][{2}]][BR[2,{1,1,1,1}]][q]/.q->Exp[((2 π I)/12)]]

FullSimplify[W[4]/.Α->1]

(W[6]=D2nInvariant[OPD[PD[BR[2,{1,1,1,1,1,1}]]],2])//AbsoluteTiming

Simplify[(W[6]/.(Α->0))-(1/4)QuantumKnotInvariant[A1,Irrep[A1][{2}]][BR[2,{1,1,1,1,1,1}]][q]/.q->Exp[((2 π I)/12)]]

FullSimplify[W[6]/.Α->1]

Clear[W]

W[k_]:=W[k]=D2nInvariant[OPD[PD[BR[2,Table[1,{2k}]]]],2]

checkW[k_]:=0==Simplify[(W[k]/.(Α->0))-(1/4)QuantumKnotInvariant[A1,Irrep[A1][{2}]][BR[2,Table[1,{2k}]]][q]/.q->Exp[((2 π I)/12)]]

Table[checkW[k],{k,2,5}]

Table[FullSimplify[W[k]/.Α->1],{k,2,6}]

Let's try D6!

DebugSet[True]

DebugPrint["hey"]

progressiveExpand[OPD[foo]]

D2nInvariant[OPD[Xp[2,4,3,1]Xp[4,2,1,3]],3]//AbsoluteTiming

WW=%[[2]]

FullSimplify[WW/.Α->0]

N[%]

(1/4)QuantumKnotInvariant[A1,Irrep[A1][{3}]][BR[2,{1,1}]][q]/.q->Exp[((2 π I)/20)]

FullSimplify[%]

FullSimplify[WW/.Α->1]

N[%]

D2nInvariant[OPD[PD[BR[2,{1,1,1,1}]]],3]

D2nInvariant[OPD[PD[BR[2,{1,1,1,1,1,1}]]],3]

FullSimplify[QuantumKnotInvariant[A1,Irrep[A1][{2}]][BR[2,{1,1}]][q]/.q->Exp[((2 π I)/12)]]

FullSimplify[QuantumKnotInvariant[A1,Irrep[A1][{2}]][PD[BR[2,{1,1,1}]]][q]/.q->Exp[((2 π I)/12)]]

FullSimplify[QuantumKnotInvariant[A1,Irrep[A1][{2}]][PD[BR[2,{1,1,1,1}]]][q]/.q->Exp[((2 π I)/12)]]

FullSimplify[QuantumKnotInvariant[A1,Irrep[A1][{2}]][PD[BR[2,{1,1,1,1,1}]]][q]/.q->Exp[((2 π I)/12)]]

FullSimplify[QuantumKnotInvariant[A1,Irrep[A1][{2}]][PD[BR[2,{1,1,1,1,1,1}]]][q]/.q->Exp[((2 π I)/12)]]

FullSimplify[QuantumKnotInvariant[A1,Irrep[A1][{2}]][PD[BR[2,{1,1,1,1,1,1,1,1}]]][q]/.q->Exp[((2 π I)/12)]]