(*The following worksheet computes the polynomial for a hyperbolic 2-bridge link that is described in http://arxiv.org/abs/1211.5089 *)





Clear["Global`*"]     (*memory cleaning*)




(*Please set up the following parameters manually, by looking a reduced alternating link diagram.*)

k=3; (*k is the number of twists in a reduced alternating diagram of a hyperbolic 2-bridge link*)

n[1]=3; n[2]=2; n[3]=3; (*n[i] is the number of crossings in the i-th twist*)

ka[1]=1; ka[2]=-1; ka[3]=1;(* Choose an orientation of the strands of the link. Then ka[i] is 1 if two strands in the i-th twist are oriented coherently, and -1 otherwise*)

(*The end of manual input*)




w[1]=x; (*variable for the polynomial, for simplicity only*)




(*Computing the Fibonacci-like polynomials for every twist*)

g[0, j_]=0; g[1, j_]=1;

g[i_, j_]:=g[i-1, j]-ka[j]*w[j]*g[i-2, j]; (*g[i, j] is the i-th polynomial for the j-th twist,w[j] is the crossing label for the j-th twist*)

(*The end of computing the Fibonacci-like polynomials*)




(*Computing labels for the j-th region with n_j twists*)

u[1, 1]:=(ka[1]*w[1]*g[n[1]-1,1])/g[n[1],1];  (*the first index of u[i,j] corresponds to the lower index of u in the ``Exact volume ..." paper*)

u[3,1]:=u[1,1];

w[2]:=-((ka[1]*w[1])^n[1])/(ka[2](g[n[1],1])^2);

u[1,j_]:=u[3,(j-1)]+(-1)^(j-1);

u[3,j_]:=If[n[j]>1, (w[j](ka[j]u[1,j]g[n[j]-1,j](-1)^(j+1)-w[j]g[n[j]-2,j]))/(u[1,j]g[n[j],j]+(-1)^(j)*ka[j]w[j]g[n[j]-1,j]), (ka[j]w[j])/u[1,j]]; (*if there are more than 1 crossing in twist goes first, if there is just one crossing in the twist goes second*)

w[j_]:=If[n[j]>1, (ka[j-2](ka[j-1]w[j-1])^(n[j-1])*w[j-2])/(ka[j]*(u[1,j-1]g[n[j-1],j-1]+(-1)^(j-1)*ka[j-1]w[j-1]g[n[j-1]-1,j-1])^2),(ka[j-2]k[j-1]w[j-2]w[j-1])/(ka[j](u[1,j])^2)];

u[3,k]:=((-1)^(k+1)*ka[k]w[k]g[n[k]-1,k])/(g[n[k],k]);

(*The end of computing labels*)




(*Computing the polynomial*)

Q:=u[1,k]*g[n[k],k]+(-1)^k*ka[k]w[k]g[n[k]-1,k]; (*Q is the resulting rational exression, P is the resultng, possibly reducible, polynomial*)

P:=Numerator[Together[Expand[Q]]]

(*The end of computing the polynomial*)




a=If[ka[1]==1,1/2+I/2, 1/2-I/2]; (*initial values for the Newton's method below*)





Print["The polynomial is  ",P]
Print ["It factors as  ", Factor[P]]   (*Wolfram Mthematica claims that the Factor function fully factors the polynomial*)
(*Print["A decimal approximation to the geometric root is  ", FindRoot[P,{x,a}, Method->"Newton"]  Solving polynomial using Newton's method and the initial values for a regular region polygon*)
Print["The decimal approximations of the roots are  ", N[Solve[P==0]]]