{-----------------------------------------------------------------------

Software for "The tree of knot tunnels"

    by Sangbum Cho and Darryl McCullough

Version of April 19, 2007
     
Contact:  dmccullough@math.ou.edu
    
This is a Haskell script, written for the hugs Haskell 98 interpreter.

------------------------------------------------------------------------

These functions calculate the cabling slope sequences for tunnels of
two-bridge knots.

The algorithm is described in the section on two-bridge knots in the paper.

------------------------------------------------------------------------

To find the slope sequence for a two-bridge knot with rational invariant  
b/a, use slopes (b/a)  (b  must be odd)

TwoBridge> slopes (9357/2434)
[ 1/3 ], 3, 3, 3, 1, 3, 5/3, 11/6, -1

To find the slope sequence for all the two-bridge knot tunnels  with invariant
b/a for a from m to n, use  slopesRange b m n

TwoBridge> slopesRange 9357 2433 2440

9357/2433 -> [ 2/3 ], -43/21, -13/6, 3
9357/2434 -> [ 1/3 ], 3, 3, 3, 1, 3, 5/3, 11/6, -1
9357/2435 -> [ 1/3 ], 19/9, 1, 3/2, -13/6, 3
9357/2436 -> [ 2/3 ], -3, -3, -5/2, 3, 1, 11/6, -1
9357/2437 -> [ 4/7 ], -5/2, 7/4, -3, -13/6, 3
9357/2438 -> [ 3/7 ], 3, 11/5, 1, 1, 11/6, -1
9357/2439 -> [ 3/5 ], 1, 1, 3/2, -3, -3, -3, -13/6, 3
9357/2440 -> [ 2/3 ], -3, -3, -9/4, 1, 1, 1, 1, 1, 1, 1, 1, 11/6, -1

-----------------------------------------------------------------------}

module TwoBridge where

import Prelude
import Ratio

--------------------------------------------------------------------------
-- continued fraction methods

unCFrac :: [Integer] -> Rational
-- convert from continued fraction back to number
unCFrac [n] = toRational n
unCFrac (n:ns) = (toRational n) + 1/unCFrac(ns)

-- functions to calculate the continued fraction with
-- all terms even, except possibly the last term

integerPart :: Rational -> Integer
integerPart x = quot (numerator ( x ) ) (denominator ( x ))

nearestEvenInteger :: Rational -> Integer
nearestEvenInteger r
  | r > 0 = 2 * (integerPart ( ( r + 1 )/2 ))
  | otherwise = (-2) * (integerPart ( ( (-r) + 1 )/2 ))

evenCFrac :: Rational -> [Integer]
evenCFrac r 
  | even ( length eCFrac ) = eCFrac
  | even ( last eCFrac )   = eCFrac
  | otherwise              = take (length eCFrac - 1) eCFrac ++
    ( if lastTerm > 1 then [ lastTerm -1, 1] else [ lastTerm + 1, -1 ] )
  where
  lastTerm = last eCFrac
  eCFrac = evenCFrac' r

evenCFrac' :: Rational -> [Integer]
-- auxiliary function for evenCFrac, allows last term odd
evenCFrac' r
  | denominator( r ) == 1 = [numerator( r )]
  | otherwise  = (nearestEvenInteger r ) : 
    ( evenCFrac' ( 1/(evenRemainder r ) ) )
  where
  evenRemainder x = x - toRational( nearestEvenInteger x )

-----------------------------------------------------------------------
-- calculation of slope coefficients for a 2-bridge knot tunnel

longCFrac :: Rational -> [Integer]
-- replace  a_i  with [2, 0, 2, ... , 0, 2 ] or [-2, 0, -2, ... , 0, -2 ]
longCFrac r = longCFrac' ( evenCFrac r )
  where
  longCFrac' [ ] = [ ]
  longCFrac' (x:y:rest) = (expandEven x) ++ [y] ++ (longCFrac' rest)

expandEven :: Integer -> [Integer]
expandEven n
  | odd n     = error "cannot expand odd integer in form [2, 0, 2, ... , 2]"
  | n == 0    = [0]
  | n > 0     = concat ( replicate (halfN-1) [2,0] ) ++ [2]
  | otherwise = map (\x -> -x ) (expandEven (-n))
  where
  halfN = fromInteger ( quot n 2 )

firstSlope :: Integer -> Integer -> Rational
-- calculates  m_0, the slope of the initial cabling
firstSlope b a
  | a ==  2  = b % ( 2 * b + 1 )
  | a == -2  = ( b - 1 ) % ( 2 * b - 1 )

kValue :: Integer -> Integer -> Integer -> Integer
-- kValue prevA b a  is the  k-value for a cabling, where
-- prevA  is  a_{i+1},  b  is  b_i, and  a  is a_i
kValue prevA b a
  | prevA ==  2 && a ==  2 = b + 1
  | prevA == -2 && a == -2 = b - 1
  | otherwise              = b

laterSlope :: Integer -> Integer -> Integer -> Integer -> Rational
-- calculates the slope  m_i  for a cabling
laterSlope parity prevA b a 
  | even parity =  2 + reciprocalK
  | odd parity  = -2 + reciprocalK
  where
  reciprocalK = 1 % ( kValue prevA b a )

slopeCoefficients :: Rational -> [Rational]
slopeCoefficients r
  | length reverseLongExpansion == 2 = [ firstSlope bn an ]
  | otherwise = ( firstSlope bn an ) : ( slopeCoefficients' bn reverseLongExpansion )
  where
  reverseLongExpansion = reverse ( longCFrac ( normalizedR ) )
  normalizedR = 1 / ( (1/r) - ( toRational ( integerPart (1/r) ) ) )
  bn = head reverseLongExpansion
  an = head (tail reverseLongExpansion )

slopeCoefficients' :: Integer -> [Integer] -> [Rational]
-- auxiliary function for slopeCoefficients
slopeCoefficients' bn (b2 : a2 : b1 : a1 : rest)
  | rest == [ ] = [ laterSlope (bn + a2ParityCorrection) a2 b1 a1 ]
  | otherwise   = ( laterSlope (bn + a2ParityCorrection) a2 b1 a1 ) :
    ( slopeCoefficients' bn ( b1 : a1 : rest ) )
  where
  a2ParityCorrection = if a2 > 0 then 1 else 0

-- output routines

prettyPrint :: Rational -> String
prettyPrint r 
  | denominator(r) == 1 = (show (numerator r))
  | otherwise           = (show (numerator r)) ++ "/" ++ (show (denominator r))

prettyPrintList :: [Rational] -> String
prettyPrintList list
  | length list == 0 = "Knot is trivial, has empty cabling sequence."
  | length list == 1 = "[ " ++ (prettyPrint (head list)) ++ " ]"
  | otherwise        = "[ " ++ (prettyPrint (head list)) ++ " ], " ++
      drop 2 ( concat ( [ ", " ++ (prettyPrint r) | r <- (tail list) ] ) )

slopes :: Rational -> IO()
slopes r
  | even (numerator r) = error "Numerator of  b/a  must be odd."
  | otherwise = ( putStr . prettyPrintList . slopeCoefficients ) r

showSlopes :: Integer -> IO()
showSlopes n 
  | even n    = error "n must be odd."
  | otherwise = putStr ( concat ( map 
    (\x -> "\n" ++ show n ++ "/" ++ show x ++ " -> " ++
      prettyPrintList ( slopeCoefficients (n % x) ) ) 
      [1..(quot (n-1) 2)] ) )

slopesRange :: Integer -> Integer -> Integer -> IO()
slopesRange n a b
  | even n    = error "n must be odd."
  | otherwise = putStr ( concat ( map 
    (\x -> "\n" ++ show n ++ "/" ++ show x ++ " -> " ++
      prettyPrintList ( slopeCoefficients (n % x) ) ) 
      [a..b] ) )

--------------------------------------------------------------------------