This mathematica notebook computes the coloured Jones polynomial following the formulas in Ohtsuki’s book

First we define the quantum integer [n], note that by definition [0]=1

quantum[n_]:=If[n==0,1,(((q^(n/2))-(q^(-(n/2))))/((q^(1/2))-(q^(-(1/2)))))
];

Then the quantum factorial

qfact[n_]:=Product[quantum[i],{i,0,n}]

For normalization purpose, we also compute the quantum dimension of an n dimensional representation

qdim[n_]:=Sum[(q^(((n-2i+1))/2)),{i,1,n}]

Now we define the action of the crossings, a.k.a. R matrix, following the formula in Ohtsuki’s book. Here Rp stands for positive R matrix, we visualize the crossing as going downwards, i, j are the labels of the top and k,l are the label of the bottom, and n is the dimension of the representation, a.k.a. the colour

Rp[i_,j_,k_,l_,n_]:=
Which[i+j≠k+l,0,i<l,0,True,
(m=i-l;
If[m>Min[k-1,n-l],0,(q^((((n-2k+1)(n-2l+1))/4)-((((n-1))^2)/4)-(((n-1))/2)))(q^((m(m-1))/4))((((q^(1/2))-(q^(-(1/2)))))^m)((qfact[n-k+m]qfact[l+m-1])/(qfact[m]qfact[n-k]qfact[l-1]))])
]

And the negative crossing is given by Rm. Here Rm stands for  the inverse of the R matrix, with the same convention as above

Rm[i_,j_,k_,l_,n_]:=
Which[i+j≠k+l,0,i>l,0,True,
(m=l-i;
If[m>Min[l-1,n-k],0,(q^(-(((n-2k+1)(n-2l+1))/4)+((((n-1))^2)/4)+((n-1)/2)))(q^(-((m(m-1))/4)))((((q^(-(1/2)))-(q^(1/2))))^m)((qfact[n-l+m]qfact[k+m-1])/(qfact[m]qfact[n-l]qfact[k-1]))])
]

Now we are ready to compute the coloured Jones polynomial. Here the input consists of 4 pieces of data: Knot, m, n, c where m is the maximum arc label and n is the colour and c is the number of components of the link. 

For the knot piece, we need to put our knot in a ‘’Morse position’’, i.e. there is a Morse projection to a vertical direction and we can slice our knot by horizontal slices so that the portion of the knot between two consecutive slices only consists of a crossing or a cup and cap (a.k.a. minimum and maximum).

For the cups and caps, note that we have two types of cups and two types of caps, namely the ones that go from left to right and the ones that go from right to left. The former give no contribution to the coloured Jones polynomial, therefore in what follows we only concern ourselves with the latter type.

We then label the arcs of our knot consecutively starting from 1 to m, where we change the label everytime we pass a crossing (note that this also applies to an arc that goes over). We now input the Knot as a product of crossings and cups and caps as we go from top to bottom. We input positive crossing as Xp[i,j,k,l] (same convention as the R matrix), negative crossing as Xm[i,j,k,l], a cap that goes from right to left of an arc labelled k as Ca[k], a cup that goes from right to left of an arc labelled k as Cu[k].

ColouredJones[Knot_,m_,n_,c_]:=Module[{t},
t=(Sum[Knot/.{
Xp[i_,j_,k_,l_]:> Rp[(Subscript[s, i]),(Subscript[s, j]),(Subscript[s, k]),(Subscript[s, l]),n],
Xm[i_,j_,k_,l_]:> Rm[(Subscript[s, i]),(Subscript[s, j]),(Subscript[s, k]),(Subscript[s, l]),n],
Cu[k_]:> (q^(-(((n-2(Subscript[s, k])+1))/2))),
Ca[k_]:> (q^(((n-2(Subscript[s, k])+1))/2))
},##]&@@({(Subscript[s, #]),1,n}&/@Range[m]))/qdim[n];
Simplify[(((-1))^(c-1))t/.q:>(q^(-1))]
]

As an example, let’s compute the coloured Jones with colour 4 of the trefoil

ColouredJones[Ca[1]Ca[4]Xp[1,4,5,2]Xp[5,2,3,6]Xp[3,6,1,4],6,4,1]

It is known that when n=2, we obtain the usual Jones polynomial, let’s check that this is indeed the case for several knots. For convenience, we include the program that computes the Jones polynomial

Jones[Knot_]:=Module[{t1,t2,t3},
t1=Expand[Knot /.{
Xp[i_,j_,k_,l_]:> -(q^(1/2))*P[i,j]*P[k,l]-q*P[i,l]*P[j,k],
Xm[i_,j_,k_,l_]:>  -(q^(-1))*P[i,j]*P[k,l]-(q^(-1/2))*P[i,l]*P[j,k]}];
 SetAttributes[P,Orderless];
  t2=t1//.P[a_,b_]*P[b_,c_]:> P[a,c];  t3=Simplify[t2/.{P[i_,i_]:> (-(q^(1/2))-(q^(-1/2))), P[i_,j_]^2:> (-(q^(1/2))-(q^(-1/2)))}];
If[FreeQ[t3,P],Simplify[t3/(-(q^(1/2))-(q^(-1/2)))],t3]];

Let’s first check the trefoil

ColouredJones[Xp[3,6,1,4]Xp[5,2,3,6]Xp[1,4,5,2]Ca[4]Ca[1],6,2,1]

Jones[Xp[1,5,2,4] Xp[3,1,4,6] Xp[5,3,6,2]]

Now the figure eight knot

ColouredJones[Xm[6,3,4,7]Xp[5,8,1,6]Xm[2,7,8,3]Xp[1,4,5,2]Ca[7]Ca[4]Ca[1],8,2,1]

Jones[Xm[6,1,7,2]Xm[2,5,3,6]Xp[8,4,1,3]Xp[4,8,5,7]]

The Hopf link

ColouredJones[Ca[3]Xp[1,3,4,2]Xp[4,2,1,3]Cu[1],4,2,2]

Jones[Xp[1,3,2,4]Xp[3,1,4,2]]

The Borromean ring, notice how slow the ColouredJones program is :(

Jones[Xp[1,5,2,8]Xm[9,2,10,3]Xm[3,6,4,7]Xm[7,12,8,9]Xp[5,11,6,10]Xp[11,1,12,4]]

ColouredJones[Ca[1]Ca[5]Ca[13]Xp[1,5,12,4]Xp[12,4,3,11]Xp[11,13,16,10]Xp[16,10,9,15]Xm[3,9,8,2]Xm[8,2,1,7]Xm[7,15,14,6]Xm[14,6,5,13],16,2,3]