########################################################################
#
#  Software for "The tree of knot tunnels"
#
#     by Sangbum Cho and Darryl McCullough
#
#  Version of April 19, 2007
#
#  Contact:  dmccullough@math.ou.edu
#
#  Written for GAP 4.4.5
#
########################################################################

########################################################################
#
#  These routines convert between the Scharlemann-Thompson
#  invariant and the principal slope invariant.
#  First, a rational number  r = q/p  with  q  odd is written
#  as a continued fraction of the form 
#           [2a_1, 2b_1, 2a_2, 2b_2, ..., 2a_n, b_n ]
#  (i. e. with an even number of terms, all of which are even
#  except possibly the last one.
#  Putting  a  equal to the sum of the a_i, the other invariant
#  is then [(-1)^p 2a, -b_n, -2a_{n-1}, ..., -2a_2, -2b_1 ].
#
########################################################################
#
#  To convert from either invariant  p/q  to the other, use  convert (p/q)
#  
#  gap> convert (54723/13363);
#  -199299/13363
#  
#  gap> convertRange(436,135,155);   
#  (p/q, p'/q')
#  
#  135/436, -4599/436
#  137/436, -8249/436
#  139/436, -4291/436
#  141/436, -1509/436
#  143/436, -9903/436
#  145/436, -63217/436
#  147/436, 14213/436
#  149/436, 4003/436
#  151/436, -231/436
#  153/436, 1687/436
#  155/436, 3533/436
#  
########################################################################
#
#  Note: Haskell scripts for all calculations discussed in the paper
#  are available for download at
#
#  http://www.math.ou.edu/~dmccullough/research/software.html
#
########################################################################

unCFrac := function( integerList )
#  Converts a continued fraction back to the rational number
#  that it represents.

   local
      length, value, termNumber;

   length := Length( integerList );
   value := integerList[ length ];
   if ( length > 1 ) then
      for termNumber in [1..length - 1] do
         value := integerList[ length - termNumber ] + 1/value;
      od;
   fi;

   return value;
   end;
      
integerPart := function( r )
#  Returns the greatest integer less than or equal to  r.
#  Assumes that r > 0.
   local
      p, q, integerQuotient;
   p := NumeratorRat( r );
   q := DenominatorRat( r );
   return QuoInt( p, q);
   end;

nearestEvenInt := function( r )
#  Returns the nearest even integer; for an odd integer, it
#  returns  r + 1 if  r > 0 and  r - 1 if  r < 0.
   if ( r > 0 ) then
      return 2 * integerPart( (r + 1)/2 );
   else
      return -2 * integerPart( (-r + 1)/2 );
   fi;
   end;

evenRemainder := function( r )
   return  r - nearestEvenInt( r );
   end;

isInteger := function( r )
   return DenominatorRat( r ) = 1;
   end;

evenCFrac := function( r )
   local   
      expansion, remainder, reciprocal, length;

   expansion := [];
   Add( expansion, nearestEvenInt( r ) );
   remainder := evenRemainder(r);

   # If the remainder is nonzero, i. e. the original rational number 
   # was not an even integer, compute the continued fraction terms 
   # until the reciprocal of the remainder is an integer.
   if ( not remainder = 0 ) then
      reciprocal := 1/remainder;
      while( not isInteger( reciprocal ) ) do
         Add( expansion, nearestEvenInt( reciprocal ) );
         remainder := evenRemainder( reciprocal );      
         reciprocal := 1/remainder;
      od;
 
      # The reciprocal is now an integer.
      # The reciprocal cannot be  1  or  -1, since if so then the 
      # previous reciprocal would have been integral.
      
      # If the current length is even, and the reciprocal is odd, then
      # the last term requires special handling to ensure that there
      # are an even number of terms:

      length := Length( expansion );
      if ( IsOddInt( length ) or IsEvenInt( reciprocal ) ) then
         Add( expansion, reciprocal );
      elif ( reciprocal > 1 ) then
         Add( expansion, reciprocal - 1); 
         Add( expansion, 1); 
      elif ( reciprocal < -1 ) then
         Add( expansion, reciprocal + 1); 
         Add( expansion, -1); 
      fi;
   fi;

   return expansion;
   end;

sumOddTerms := function( integerList )
   local
      term, sum;
   sum := 0;
   term := 1;
   while( 2 * term - 1 <= Length( integerList ) ) do
      sum := sum + integerList[ 2 * term - 1 ];
      term := term + 1;
   od;
   return sum;
   end;

convertIntegerList := function( integerList, p )
   local
      sumOdd, newExpansion;

   sumOdd := sumOddTerms( integerList );
   newExpansion := [];
   Add( newExpansion, (-1)^p * sumOdd );
   Remove( integerList, 1);
   Append( newExpansion, -Reversed( integerList ) );
   return newExpansion;
   end;

convert := function( r )
   return unCFrac( convertIntegerList( evenCFrac( r ), DenominatorRat( r ) ) );
   end;

convertRange := function (q, m, n)
   local
      p, step;

   Print("\n(p/q, p'/q')\n\n");

   if m = n then
      Print(p,"/",q,", ", convert( p/q ),"\n\n");
      return;
   fi;

   if m < n then
      step := 2;
   else
      step := -2;
   fi;

   if m mod 2 = 0 then
      if m < n then
         m := m + 1;
      else
         m := m - 1;
      fi;
   fi;

   if n mod 2 = 0 then
      if m < n then
         n := n - 1;
      else
         n := n + 1;
      fi;
   fi;

   for p in [m, m + step..n] do
      if GcdInt (p, q) = 1 then
         Print(p,"/",q,", ", convert( p/q ),"\n");
      fi;
   od;
   Print("\n");
   end;