/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 0.8.2 ] */

/* [wxMaxima: comment start ]
This is a routine for computing the A-polynomials of the first few two-bridge knots. 
Along the way it computes the equations for the spaces of representations of the fundamental groups 
to SL(2,C) (Riley's equation or Le's equation). The method is described in the original paper 
where the A-polynomial was defined (see e.g. my papers for references and details). 
By rewriting the routine in the free computer algebra system Maxima 
(http://maxima.sourceforge.net) I hope it might become more accessible.

July 2009
Vu Quang Huynh
http://www.math.hcmus.edu.vn/~hqvu
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
trigreduce(%);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
reset;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
omega(p,q):=block(
/*computing the word w in the presentation of the fundamental group*/
local(word,i,epsilon),
word:[],
for i:1 thru (p-1) do
(epsilon:(-1)^(floor(i*q/p)),
if (mod(i,2)=1) then  word:endcons(a^(epsilon),word)
else word:endcons(b^(epsilon),word)),
return(word));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
omega_exponent(p,q):=block(
/* computing the sum of the exponents of the word w in the presentation of the fundamental group */
local(value,i),
 
value:0,
for i:1 thru (p-1) do (epsilon:(-1)^(floor(i*q/p)),
value:value+epsilon),
return(value));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
lamda(p,q):=block(
/* description `the longitude` */
local(omega1,value1,i),
 
omega1:omega(p,q),
value1:omega1,
for i:1 thru (p-1) do value1:endcons(omega1[p-i],value1),
value1:endcons(b^(-2*omega_exponent(p,q)),value1),
return(value1));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
rho(x):=block(
local(value),
if x=1 then value:ident(2) else
if x=a then value:matrix([M,1],[0,M^(-1)]) else
if x=1/a then value:matrix([M^(-1),-1],[0,M]) else
if x=b then value:matrix([M,0],[-z,M^(-1)]) else
if x=1/b then value:matrix([M^(-1),0],[z,M]),
return(value));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
Le_polynomial(p,q):=block(
/* compute the defining polynomial of the non-abelian part of representations to SL(2,C),*/
/* according to Le93, p. 198 */
local(omega1,omega1_prime_matrix,n,phi,i),

omega1:omega(p,q),
n:(p-1)/2,
phi:(-1)^n,
omega1_prime_matrix:ident(2),
for i:1 thru n do
  (omega1_prime_matrix:expand(rho(omega1[n-i+1]).omega1_prime_matrix.rho(omega1[n+i])),
   phi:phi+(-1)^(n-i)*(omega1_prime_matrix[1,1]+omega1_prime_matrix[2,2])),
return(ratnumer(phi)));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
Riley_polynomial(p,q):=block(
/* Riley polynomial describing the nonabelian part of the character variety of */
/* a 2-bridge knot (p,q) */
local(C,D,w,epsilon,i),
 
C:matrix([M,1],[0,M^(-1)]),
D:matrix([M,0],[-z,M^(-1)]), 
w:ident(2),
for i:1 thru (p-1) do
  (epsilon:(-1)^(floor(i*q/p)),
   if (mod(i,2)=1) then  w:expand(w.C^^(epsilon))
   else w:expand(w.D^^(epsilon))),
return(ratnumer(w[1,1]+(M^(-1)-M)*w[1,2])));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
L(p,q):=block(
/* find the variable L, thereby giving the A-polynomial */
local(i,value,omega1, omega_exponent, epsilon),

value:ident(2),
omega1:omega(p,q),
for i:1 thru (p-1) do
  value:expand(value.rho(omega1[i])),

for i:1 thru (p-1) do
  value:expand(value.rho(omega1[p-i])),

omega_exponent:0,
for i:1 thru (p-1) do
  (epsilon:(-1)^(floor(i*q/p)),
   omega_exponent:omega_exponent+epsilon),

value:expand(value.(rho(b)^^(-2*omega_exponent))),
value:value[1,1],
return(value));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
A_polynomial(p,q):=block(
local(value),

LL(z):=ratnumer(L(p,q)-L),
phi(z):=Riley_polynomial(p,q),
value:resultant(phi(z),LL(z),z),
return(factor(expand(value))));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
trace_lambda(p,q):=block(
/* find the trace of the word lambda - the longitude */
local(i,value,omega1, omega_exponent, epsilon),

/* M:(x+(x^2-4)^(1/2))/2,*/
value:ident(2),
omega1:omega(p,q),
for i:1 thru (p-1) do
  value:expand(value.rho(omega1[i])),

for i:1 thru (p-1) do
  value:expand(value.rho(omega1[p-i])),

omega_exponent:0,
for i:1 thru (p-1) do
  (epsilon:(-1)^(floor(i*q/p)),
   omega_exponent:omega_exponent+epsilon),

value:expand(value.(rho(b)^^(-2*omega_exponent))),
value:value[1,1]+value[2,2],
return(ratsymp(value));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
omega(7,3);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
Riley_polynomial(7,3);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
Le_polynomial(7,3);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
lamda(7,3);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
L(7,3);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
A_polynomial(7,3);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
trace_lambda(7,3);
/* [wxMaxima: input   end   ] */

/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$