SetDirectory[NotebookDirectory[]];
<<RISC`LinearSystemSolver`
<<"BimonicRecurrences.m";
temp=""; (* directory for storing temporary results *)
pretz=""; (* path to pretzel knot recurrences *)
twist=""; (* path to twist knot recurrences *)

The family of (-2,3,2p+3)-pretzel knots (Subscript[K, p]^pretz)

p=-2

(* We load the precomputed (inhomogeneous) q-recurrence for this knot *)
rec=Get[pretz<>"CJK.-2.1.rec.m"]/.f[a_]:>f[a]/(q^a-1)

(* We make some experiments to find out the optimal order and degree setting. *)
TableForm[Table[BimonicRecurrence[rec,ord,deg,Power->2*qn1+qn2,Return->"nsol"],{ord,4,8},{deg,3,10},{qn1,0,1},{qn2,0,1}]]

(* In case there are old temporary files, delete them. *)
Run["rm vec_*.m"];

(* The computation is done via evaluation-interpolation and rational reconstruction. *)
(* This computes the solution for certain evaluation points for q and modulo prime. *)
(* The results are stored in files. *)
BimonicRecurrence[rec,5,4,Power->1,Evaluate->Range[3,40],Modulus->{NextPrime[2^31,-2]}]

(* Next we do interpolation and rational function reconstruction. *)
(* Again, for each prime the result is stored in a file. *)
InterpolateAndReconstruct["vec",{NextPrime[2^31,-2]}]

(* Finally, chinese remaindering and rational number reconstruction. *)
(* The result is stored in a file named "vec_reconst.m". *)
RationalCoefficients["vec"]

(* From the previous result, we obtain the coefficients u_i for the linear combination. *)
reconst=Get["vec_reconst.m"];
ansatz=BimonicRecurrence[rec,5,4,Power->1,Return->"ansatz"];
lincomb=ansatz/.C[i_]:>reconst[[i]]

(* This is how the ansatz looks like. *)
ansatz

(* Here we demonstrate how to compute the C[i] "by hand". *)
(* (1) Plug in the recurrence into the ansatz (and replace q^n by a new variable qn) *)
eqns=ansatz/.f[a_]:>(rec/.n->a)/.q^a_:>q^Expand[a]/.q^(a_.*n+b_.)->qn^a*q^b;
(* (2) Coefficient comparison w.r.t. f[_] *)
eqns=Prepend[Coefficient[eqns,Table[f[n+i],{i,0,5}]],eqns/.f[_]->0];
(* (3) The common factor that we want to have (and that we want to cancel in the end) *)
fac=(-1+q^3*qn)*(1+q^3*qn)*(-1+q^4*qn^2)*(-1+q^5*qn^2);
(* (4) Polynomial division by fac. We impose the condition that the remainder is zero. *)
eqns=DeleteCases[PolynomialRemainder[#,fac,qn]&/@eqns,0];
(* (5) Coefficient comparison w.r.t. qn *)
eqns=Flatten[CoefficientList[#,qn]&/@eqns];
(* (6) Solve the linear system *)
Solve[Thread[eqns==0],Array[C,7]]

(* Performing this linear combination to rec yields the desired bimonic recurrence... *)
bmrec=Collect[lincomb/.f[a_]:>(rec/.n->a),f[_]];
(* ... in which we still have to divide by the common factor. *)
bmrec=Collect[bmrec/(Coefficient[bmrec,f[n+5]]/q^n),f[_],Together]

(* This result is stored in the following file: *)
Expand[bmrec-Get["CJK.-2.1.bmrec.m"]]

ListPlot[Monomials/@{bmrec,rec},PlotMarkers->{MyDisk[15,Blue],MyDisk[8,Red]},PlotRange->{{-2,6},{-1,12}}]

Palindromicity

(* The bimonic recurrence is somewhat palindromic. *)
Expand[bmrec]-(Expand[bmrec]/.q^a_*f[n+b_.]:>q^Expand[(a/.n->-n-5)+9n+20+b]*f[n+5-b])
inhom=Expand[bmrec/.f[_]->0];
inhom+(inhom/.q^a_:>q^Expand[(a/.n->-n-5)+10n+25])

Table[Together[Coefficient[Coefficient[bmrec,f[n+i]],q^n,j]-q^(5j-i-20)*Coefficient[Coefficient[bmrec,f[n+5-i]],q^n,9-j]],{i,0,5},{j,0,9}]
Table[Together[Coefficient[bmrec/.f[_]->0,q^n,j]+q^(5j-25)*Coefficient[bmrec/.f[_]->0,q^n,10-j]],{j,0,9}]

bmrecCoeffs=Prepend[Coefficient[bmrec,Table[f[n+i],{i,0,5}]],bmrec/.f[_]->0];
qnCoeffs=Transpose[PadRight[Factor[CoefficientList[#,q^n]&/@bmrecCoeffs]]];
TableForm[qnCoeffs]

TableForm[Flatten/@Map[If[#===0,{"",""},With[{m=Exponent[#,q,Min]},{(#/q^m)/.q->2,q^m}]]&,qnCoeffs,{2}]]

(* Load the first few values of the colored Jones polynomial of this knot and apply the normalization. *)
data=Get["CJK.-2.1.data.m"]*Table[q^n-1,{n,31}];

(* Consistency check *)
Expand[Table[{rec,bmrec},{n,1,25}]/.f[n_]:>data[[n]]]

(* Use the bimonic recurrence to compute f[n] at negative integers. *)
data1={};
Do[
AppendTo[data1,Expand[Together[Solve[(bmrec/.n->nn/.f[n_/;n>0]:>data[[n]]/.f[n_/;n≠nn]:>data1[[1-n]])==0,f[nn]][[1,1,2]]]]];
,{nn,0,-31,-1}]

(* We get an almost palindromic sequence *)
Table[Together[data[[i]]/data1[[i+1]]],{i,31}]

p=2

rec=Get[pretz<>"CJK.2.1.rec.m"]/.f[a_]:>f[a]/(q^a-1);

TableForm[Table[BimonicRecurrence[rec,ord,deg,Power->2*qn1+qn2,Return->"nsol"],{ord,9,12},{deg,0,5},{qn1,0,1},{qn2,0,1}]]

Run["rm vec_*.m"];
BimonicRecurrence[rec,10,1,Evaluate->Range[3,120],Modulus->{NextPrime[2^31,-2]}];
InterpolateAndReconstruct["vec",{NextPrime[2^31,-2]}];
RationalCoefficients["vec"];
reconst=Get["vec_reconst.m"]

ansatz=BimonicRecurrence[rec,10,1,Return->"ansatz"];
lincomb=ansatz/.C[i_]:>reconst[[i]];
bmrec=Collect[lincomb/.f[a_]:>(rec/.n->a),f[_]];
bmrec=Collect[bmrec/Coefficient[bmrec,f[n+10]],f[_],Together];

(* This result is stored in the following file: *)
Collect[bmrec+Get["CJK.2.1.bmrec.m"],f[_],Together]

ListPlot[Monomials/@{lincomb,bmrec,rec},PlotMarkers->{MyDisk[8,Green],MyDisk[8,Blue],MyDisk[5,Red]},ImageSize->600]

(* CJK.2.1.bmrec is somewhat palindromic. *)
inhom=Expand[bmrec/.f[_]->0];
hom=Expand[q*(bmrec-inhom)]; (* multiply by q to make sure that all terms are matched *)hom+(hom/.q^a_.*f[n+b_.]:>q^Expand[(a/.n->-n-10)+55n+270+b]*f[n+10-b])
inhom-(inhom/.q^a_.:>q^Expand[(a/.n->-n-10)+56n+280])

p=-3

rec=Get[pretz<>"CJK.-3.1.rec.m"]/.f[a_]:>f[a]/(q^a-1);
Max[Cases[rec,f[n+a_.]:>a,Infinity]]

TableForm[Table[BimonicRecurrence[rec,ord,deg,Return->"nsol"],{ord,8,10},{deg,10,40,5}]]

TableForm[Table[BimonicRecurrence[rec,ord,deg,Return->"nsol"],{ord,10,11},{deg,30,35}]]

Table[BimonicRecurrence[rec,10,30,Power->n0,Return->"nsol"],{n0,0,30}]

TableForm[Table[BimonicRecurrence[rec,10,30,Power->11,Return->"nsol",Begin->i,End->i+j],{i,0,10},{j,0,10-i}]]

(* Try to minimize the size of the linear system. *)
BimonicRecurrence[rec,10,30,Power->11,Begin->#1,End->#2,Return->"dims",Method->1]&@@@{{0,5},{1,5},{2,5},{3,7},{4,8},{5,8}}

Timing[BimonicRecurrence[rec,10,30,Power->11,Begin->2,End->5,Evaluate->Range[3,415],Modulus->Table[NextPrime[2^31,-i],{i,2,5}],Method->1,FileName->temp<>"cjk-31"]]

BimonicRecurrence[rec,10,30,Power->11,Begin->#1,End->#2,Return->"dims",Method->2]&@@@{{0,5},{1,5},{2,5},{3,7},{4,8},{5,8}}

Timing[BimonicRecurrence[rec,10,30,Power->11,Begin->2,End->5,Evaluate->Range[3,415],Modulus->Table[NextPrime[2^31,-i],{i,2,5}],Method->2,FileName->temp<>"cjk-31"]]

Timing[InterpolateAndReconstruct[temp<>"cjk-31",Table[NextPrime[2^31,-i],{i,2,5}]]]

Timing[RationalCoefficients[temp<>"cjk-31"]]

reconst=Get[temp<>"cjk-31_reconst.m"];
Max[Abs[Cases[reconst,_Integer,Infinity]]]

ansatz=BimonicRecurrence[rec,10,30,Power->11,Begin->2,End->5,Return->"ansatz"];
lincomb=ansatz/.C[i_]:>reconst[[i]];
Collect[lincomb+Get["CJK.-3.1.lincomb.m"],f[_],Together]

(* Delete temporary files *)
Run["rm "<>temp<>"cjk-31_*"];

bmrec=Collect[lincomb/.f[a_]:>(rec/.n->a),f[_]];
cf=Select[Factor[Coefficient[bmrec,f[n+10]]]/q^(11n),Not[FreeQ[#,n]]&];
bmrec=Collect[bmrec/cf,f[_],Together];
(* Put[Collect[bmrec,f[_],Collect[#,q^n,Expand]&],"CJK.-3.1.bmrec.m"]; *)

(* This is the content that we can remove. *)
cf

(* Indeed, after dividing all coefficients are still polynomials. *)
Denominator[Together[Join[{bmrec/.f[_]->0},Coefficient[bmrec,Table[f[n+i],{i,0,10}]]]]]

(* There is are large factors in Q[q], that's why we look only at an image of the leading/trailing coefficient. *)
PolynomialMod[Coefficient[bmrec,{f[n],f[n+10]}]/.q->2,2^13-1]

ListPlot[Monomials/@{lincomb,bmrec,rec},PlotMarkers->{MyDisk[10,Green],MyDisk[8,Blue],MyDisk[5,Red]}]

Length[Monomials[#]]&/@{rec,lincomb,bmrec}

(* Palindromicity *)
Union[Flatten[Table[Together[Coefficient[Coefficient[bmrec,f[n+i]],q^n,j]+q^(10j-i-110)*Coefficient[Coefficient[bmrec,f[n+10-i]],q^n,23-j]],{j,0,23},{i,0,10}]]]
Union[Table[Together[Coefficient[bmrec/.f[_]->0,q^n,j]-q^(10j-120)*Coefficient[bmrec/.f[_]->0,q^n,24-j]],{j,0,23}]]

Try to exploit palindromicity

mat=Mod[BimonicRecurrence[rec,10,30,Power->11,Begin->2,End->5,Return->"mat"],NextPrime[2^31,-2]];
pows=Palindromicity[mat,95,-1,NextPrime[2^31,-2]]

(* By exploiting palindromicity, we can reduce the number of unknowns. *)
(* Unfortunately, this doesn't reduce the computation time considerably. *)
Timing[BimonicRecurrence[rec,10,30,Power->11,Begin->2,End->5,Evaluate->Range[3,415],Modulus->Table[NextPrime[2^31,-i],{i,2,5}],FileName->temp<>"cjk-31",MatrixForm->((Join[Take[#,2],Take[#,{3,48}]-PowerMod["q",ps,"mod"]*Reverse[Take[#,{50,95}]],{#[[49]]}]&)/.ps->pows)]]

Timing[InterpolateAndReconstruct[temp<>"cjk-31",Table[NextPrime[2^31,-i],{i,2,5}]]]
Timing[RationalCoefficients[temp<>"cjk-31"]]

(* We get the same result as before. *)
reconst=Get[temp<>"cjk-31_reconst.m"];
reconst=Join[Take[#,2],Take[#,{3,48}],{#[[49]]},Reverse[-Take[#,{3,48}]*q^pows]]&[reconst];
reconst=(PolynomialLCM@@Denominator[reconst])*reconst;
lincomb=ansatz/.C[i_]:>reconst[[i]];
bmrec=Collect[lincomb/.f[a_]:>(rec/.n->a),f[_]];
cf=Select[Factor[Coefficient[bmrec,f[n+10]]]/q^(11n),Not[FreeQ[#,n]]&];
bmrec=Collect[bmrec/cf,f[_],Together];
Collect[lincomb+Get["CJK.-3.1.lincomb.m"],f[_],Together]
Collect[bmrec-Get["CJK.-3.1.bmrec.m"],f[_],Together]

(* Delete temporary files *)
Run["rm "<>temp<>"cjk-31_*"];

p=3

rec=Get[pretz<>"CJK.3.1.rec.m"]/.f[a_]:>f[a]/(q^a-1);
Max[Cases[rec,f[n+a_.]:>a,Infinity]]

TableForm[Table[BimonicRecurrence[rec,ord,deg,Return->"nsol"],{ord,10,15},{deg,5,40,5}]]

TableForm[Table[BimonicRecurrence[rec,ord,deg,Return->"nsol"],{ord,12,20,2},{deg,10,40,5}]]

TableForm[Table[BimonicRecurrence[rec,ord,deg,Return->"nsol"],{ord,15,16},{deg,25,30}]]

Table[BimonicRecurrence[rec,16,25,Power->qn,Return->"nsol"],{qn,0,50}]

TableForm[Table[BimonicRecurrence[rec,16,25,Power->6,Return->"nsol",Begin->i,End->i+j],{i,0,16},{j,0,16-i}]]

BimonicRecurrence[rec,16,25,Power->6,Return->"nsol",Begin->#1,End->#2]&@@@{{0,9},{1,9},{2,9},{3,10},{4,10},{5,11},{6,12},{7,14}}

BimonicRecurrence[rec,16,25,Power->6,Return->"dims",Begin->#1,End->#2]&@@@{{0,9},{1,9},{2,9},{3,10},{4,10},{5,11},{6,12},{7,14}}

BimonicRecurrence[rec,16,25,Power->6,Return->"dims",Begin->#1,End->#2,Method->1]&@@@{{0,9},{1,9},{2,9},{3,10},{4,10},{5,11},{6,12},{7,14}}

BimonicRecurrence[rec,16,{23,24,24,25,25,25},Power->6,Return->"nsol"]

BimonicRecurrence[rec,16,{25,25,25,24,24,23},Power->6,Return->"nsol"]

BimonicRecurrence[rec,16,{{2,25},{1,25},{0,25},{0,25},{0,24},{0,23}},Power->6,Return->"nsol"]

BimonicRecurrence[rec,16,{{0,23},{0,24},{0,25},{0,25},{1,25},{2,25}},Power->6,Return->"nsol"]

Timing[
Run["rm "<>temp<>"cjk31*"];
BimonicRecurrence[rec,16,{23,24,24,25,25,25},Power->6,Begin->5,End->11,Evaluate->Range[3,500],Modulus->Table[NextPrime[2^31,-i],{i,2,6}],FileName->temp<>"cjk31"];
InterpolateAndReconstruct[temp<>"cjk31",Table[NextPrime[2^31,-i],{i,2,6}]];
RationalCoefficients[temp<>"cjk31"];
]

ByteCount[reconst=Get[temp<>"cjk31_reconst.m"]]
Max[Abs[Cases[reconst,_Integer,Infinity]]]
Table[Flatten[Position[Together[reconst/reconst[[i]]],q^_.,{1}]],{i,3,77}]

Timing[
Run["rm "<>temp<>"cjk31*"];
BimonicRecurrence[rec,16,{{2,25},{1,25},{0,25},{0,25},{0,24},{0,23}},Power->6,Begin->5,End->11,Evaluate->Range[3,600],Modulus->Table[NextPrime[2^31,-i],{i,2,6}],FileName->temp<>"cjk31"];
InterpolateAndReconstruct[temp<>"cjk31",Table[NextPrime[2^31,-i],{i,2,6}]];
RationalCoefficients[temp<>"cjk31"];
]

ByteCount[reconst=Get[temp<>"cjk31_reconst.m"]]
Max[Abs[Cases[reconst,_Integer,Infinity]]]
Table[Flatten[Position[Together[reconst/reconst[[i]]],q^_.,{1}]],{i,3,77}]

Timing[
Run["rm "<>temp<>"cjk31*"];
BimonicRecurrence[rec,16,{{0,23},{0,24},{0,25},{0,25},{1,25},{2,25}},Power->6,Begin->5,End->11,Evaluate->Range[3,500],Modulus->Table[NextPrime[2^31,-i],{i,2,6}],FileName->temp<>"cjk31"];
InterpolateAndReconstruct[temp<>"cjk31",Table[NextPrime[2^31,-i],{i,2,6}]];
RationalCoefficients[temp<>"cjk31"];
]

ByteCount[reconst=Get[temp<>"cjk31_reconst.m"]]
Max[Abs[Cases[reconst,_Integer,Infinity]]]
Table[Flatten[Position[Together[reconst/reconst[[i]]],q^_.,{1}]],{i,3,77}]

Timing[
ansatz=BimonicRecurrence[rec,16,{{0,23},{0,24},{0,25},{0,25},{1,25},{2,25}},Power->6,Return->"ansatz"];
lincomb=ansatz/.C[i_]:>reconst[[i]];
Put[lincomb,"CJK.3.1.lincomb.m"];
bmrec=Collect[lincomb/.f[a_]:>(rec/.n->a)/.q^a_:>q^Expand[a],f[_]];
cf=Select[Factor[Coefficient[bmrec,f[n+16]]]/q^(6n),Not[FreeQ[#,n]]&];
bmrec=Collect[bmrec/cf,f[_],Together];
bmrec=Collect[bmrec,f[_],Function[c,Factor[c]/.a_Plus:>Collect[a,q^n,With[{d=Exponent[#,q,Min]},q^d*Expand[#/q^d]]&]]];
Put[bmrec,"CJK.3.1.bmrec.m"];
]

mons=Monomials/@{lincomb,bmrec,rec};
ListPlot[mons,PlotMarkers->{MyDisk[6,Green],MyDisk[6,Blue],MyDisk[4,Red]}]

MyEll[r_,c_]:=Graphics[{c,Disk[{0,0},{2,1}]},ImageSize->{2r,r}];
gr=ListPlot[mons,PlotMarkers->{MyEll[2,RGBColor[0.2,0.7,0]],MyEll[2,Blue],MyDisk[2,Red]}];
Export["support_pretz_3.pdf",gr];

p=-4

rec=Get[pretz<>"CJK.-4.1.rec.m"]/.f[a_]:>f[a]/(q^a-1);
rec=Collect[rec,f[_],Together];
Max[Cases[rec,f[n+a_.]:>a,Infinity]]

TableForm[Table[BimonicRecurrence[rec,ord,deg,Return->"nsol"],{ord,12,16},{deg,60,80,4}]]

Table[BimonicRecurrence[rec,16,deg,Return->"nsol"],{deg,68,72}]

Table[BimonicRecurrence[rec,16,71,Power->qn,Return->"nsol"],{qn,0,100}]

Position[%,1]-1

TableForm[Table[BimonicRecurrence[rec,16,71,Power->22,Return->"nsol",Begin->i,End->i+j],{i,0,16},{j,0,16-i}]]

BimonicRecurrence[rec,16,71,Power->22,Return->"nsol",Begin->#1,End->#2]&@@@{{0,8},{1,8},{2,9},{3,9},{4,9},{5,10},{6,11},{7,12},{8,15}}

TableForm[(T@@Timing[BimonicRecurrence[rec,16,71,Power->22,Return->"dims",Begin->#1,End->#2]])&@@@{{0,8},{1,8},{2,9},{3,9},{4,9},{5,10},{6,11},{7,12},{8,15}}]

TableForm[(T@@Timing[BimonicRecurrence[rec,16,71,Power->22,Method->1,Return->"dims",Begin->#1,End->#2]])&@@@{{0,8},{1,8},{2,9},{3,9},{4,9},{5,10},{6,11},{7,12},{8,15}}]

TableForm[Table[First[Timing[BimonicRecurrence[rec,16,71,Method->i,Power->22,Begin->4,End->9,Evaluate->Range[3,r],FileName->temp<>"cjk-41"]]],{i,2},{r,12,22,10}]]

Timing[BimonicRecurrence[rec,16,71,Power->22,Begin->4,End->9,Evaluate->Range[3,1000],FileName->temp<>"cjk-41"]]

Timing[BimonicRecurrence[rec,16,71,Power->22,Begin->4,End->9,Evaluate->Range[1001,1500],FileName->temp<>"cjk-41"]]

Timing[BimonicRecurrence[rec,16,71,Power->22,Begin->4,End->9,Evaluate->Range[1501,2000],FileName->temp<>"cjk-41"]]

Timing[InterpolateAndReconstruct[temp<>"cjk-41",Table[NextPrime[2^31,-i],{i,2,2}],Print->True]]

Max[Get[temp<>"cjk-41_qdegs.m"]]

(* 24798s on thebrain *)
InterpolateAndReconstruct[temp<>"cjk-41",Table[NextPrime[2^31,-i],{i,3,6}],Print->True];
(* In the end we used primes {i, 2, 15}. *)

(* 335s on thebrain *)
RationalCoefficients[temp<>"cjk-41"];

reconst=Get[temp<>"cjk-41_reconst.m"];
{Dimensions[#],ByteCount[#],N[Max[Abs[Cases[#,_Integer,Infinity]]]]}&[reconst]

Timing[
ansatz=BimonicRecurrence[rec,16,71,Power->22,Return->"ansatz"];
lincomb=ansatz/.C[i_]:>reconst[[i]];
Put[lincomb,"CJK.-4.1.lincomb.m"];
bmrec=Collect[lincomb/.f[a_]:>(rec/.n->a)/.q^a_:>q^Expand[a],f[_]];
cf=(Coefficient[ansatz,f[n]]/.C[1]->1)*Coefficient[rec,f[n]];
cf=Together[cf/q^(Exponent[cf,q^n,Min]n)];
cf=Together[cf/(PolynomialGCD@@CoefficientList[cf,q^n])];
]

(* 167736s on thebrain *)
Timing[
bmrec=Collect[bmrec,f[_],(quo=Together[#/cf];Print[PolynomialQ[quo,{q,q^n}],"   ",DateString[]];quo)&];
]

bmrec=Collect[bmrec,f[_],Function[c,Factor[c]/.a_Plus:>Collect[a,q^n,With[{d=Exponent[#,q,Min]},q^d*Expand[#/q^d]]&]]];

bmrec=Get["CJK.-4.1.bmrec.m"];
ByteCount[bmrec]

Timing[mons=Monomials/@{lincomb,bmrec,rec};]

ListPlot[mons+{-.1,0,.1},PlotMarkers->{MyDisk[5,Green],MyDisk[5,Blue],MyDisk[5,Red]}]

The family of twist knots (Subscript[K, p]^twist)

(* Formulas from Garoufalidis+Sun paper *)
Clear[J,Jhat,myc];
Jhat[p_,n_]:=Jhat[p,n]=
Together[Sum[q^(n(n+3)/2+p k(k+1)+k(k-1)/2)*(-1)^(n+k+1)*(q^(2k+1)-1)*FunctionExpand[QPochhammer[q^(n-k+1),q,k]/QPochhammer[q,q,n+k+1]],{k,0,n}]];
myc[n_,k_]:=myc[n,k]=(-1)^k*q^(-k(k+1)/2)*FunctionExpand[QPochhammer[q^(1-n),q,k]*QPochhammer[q^(1+n),q,k]];
J[p_,n_]:=J[p,n]=Expand[Together[Sum[myc[n,k]*Jhat[p,k],{k,0,n}]]];
Table[J[-1,n],{n,4}]

(* Simplified formula (putting everything together). *)
Clear[J];
J[p_,n_]:=J[p,n]=
Expand[Together[Sum[
(-1)^(j+1)*q^(k+p j(j+1)+j(j-1)/2)*(q^(2j+1)-1)*FunctionExpand[QPochhammer[q^(1-n),q,k]*QPochhammer[q^(1+n),q,k]*QPochhammer[q^(k-j+1),q,j]/QPochhammer[q,q,k+j+1]],
{k,0,n},{j,0,k}]]];

(* Consistency check *)
Table[
test=Get[twist<>"rec.twist.knot."<>ToString[p]<>".m"];
Expand[Table[test,{n,5}]/.f[n_]:>J[p,n]],{p,-5,5}]

(* For the calculations, see twistknots.calc.m and twistknots.out. *)

TwistPlot[tw_Integer,r_:{6,6,4}]:=
Module[{s=ToString[tw],r1,r2,r3},
{r1,r2,r3}=r;
Return[ListPlot[Monomials/@{
Get[twist<>"lincomb.twist."<>s<>".m"],
Get[twist<>"bmrec.twist."<>s<>".m"],
Get[twist<>"rec.twist.knot."<>s<>".m"]},
PlotMarkers->{MyDisk[r1,Green],MyDisk[r2,Blue],MyDisk[r3,Red]}]];
]

TwistPlot[-1,{25,15,10}]

TwistPlot[2,{20,12,7}]

TwistPlot[-2,{15,10,5}]

TwistPlot[3,{15,10,5}]

TwistPlot[-3,{12,8,5}]

TwistPlot[4,{12,8,5}]