#!/usr/bin/python -i

########################################################################
#
#  Software for "The tree of knot tunnels"
#
#     by Sangbum Cho and Darryl McCullough
#
#  Version of July 31, 2007
#
#  Contact:  dmccullough@math.ou.edu
#
#  This is a Python script, written for Python 2.3
#
########################################################################

########################################################################
#
#  The first set of 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
#
########################################################################
#
#
#  The next set of 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):
#
#  Tree> 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 ):
#  
#  Tree> 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
#  
#  
########################################################################
#
#  Note: Haskell scripts for all calculations discussed in the paper
#  are also available for download at
#
#  http://www.math.ou.edu/~dmccullough/research/software.html
#
########################################################################

import re, sys
sys.ps1 = 'Tree> '
sys.ps2 = '   ... '

def gcd ( a, b ):
    """Greatest common divisor of two integers"""
    while b != 0:
        a, b = b, a%b
    return abs( a )

class Rational:
    """Class of rational numbers. Zero denominators are allowed. Division
    by the zero rational number yields zero."""
    def __init__( self, num=0 , denom=1 ):
        self.num = num
        self.denom = denom
        d = gcd( num, denom )
        if d is not 0:
            self.num = self.num/d
            self.denom = self.denom/d
        if self.denom < 0:
            self.num = -self.num
            self.denom = -self.denom
        if self.num == 0:
            self.denom = 1

    def __str__( self ):
        if self.denom == 1:
            return "%d" % self.num
        else:
            return "%d/%d" % (self.num, self.denom)
    
    def reciprocal( self ):
        return( Rational( self.denom, self.num ) )

    # define operations with first operand rational and second
    # either rational or integer
    
    def __add__( self, r ):
        if isinstance( r, Rational):
            return( Rational( self.num * r.denom + self.denom * r.num,
                              self.denom * r.denom ) )
        elif isinstance( r, int):
            return( Rational( self.num + self.denom * r, self.denom ) )
        
    def __sub__( self, r ):
        if isinstance( r, Rational):
            return( Rational( self.num * r.denom - self.denom * r.num,
                              self.denom * r.denom ) )
        elif isinstance( r, int):
            return( Rational( self.num - self.denom * r, self.denom ) )            
 
    def __mul__( self, r ):
        if isinstance( r, Rational):
            return( Rational( self.num * r.num, self.denom * r.denom ) )
        elif isinstance( r, int):
            return( Rational( self.num * r , self.denom ) )

    def __div__( self, r ):
        if isinstance( r, Rational):
            return( Rational( self.num * r.denom, self.denom * r.num ) )
        elif isinstance( r, int):
            return( Rational( self.num, self.denom * r ) )

# define a couple of handy list-manipulation functions

def takeWhile( f, lst ):
    """Return the initial segment of elements satisfying condition f."""
    returnList = []
    for aTerm in lst:
        if f( aTerm ):
            returnList.append( aTerm )
        else:
            return returnList
    return returnList

def dropWhile( f, lst ):
    """Return the list with the initial segment of elements that
    satisfy condition f dropped."""
    return lst[ len( takeWhile( f, lst ) ): ]

def isZero( n ):
    return n is 0

def isNonzero( n ):
    return n is not 0

class CFrac:
    """Class of continued fractions.
    The constructor function creates a continued fraction from a list
    of integers. Its terms list is rewritten without zero entries."""

    def __init__( self, terms = [] ):
        self.terms = terms
        self.stripZeros()

    def stripZeros( self ):
        """Rewrite the sequence of terms to eliminate 0 terms, except possibly
        one initial 0, using the formulas
        [ ..., a_{i-1}, 0, a_{i+1}, ... ] = [ ..., a_{i-1} + a_{i+1}, ... ] and
        [ ..., a_{n-2}, a_{n-1}, 0 ] = [ ..., a_{n-2} ]."""
        if len(self.terms) is 0:
            print "Warning: Continued fraction with empty terms list."
            return True
        if len(self.terms) is 1:
            return True
        while len( self.terms ) > 2 and self.terms[-1] is 0:
            self.terms = self.terms[:-2]
        if len( self.terms ) is 2 and self.terms[1] is 0:
            print "Warning: Continued fraction has infinite value."
            return True
        if len( takeWhile( isZero, self.terms ) ) > 0:
            initialZero=[0]
            self.terms = dropWhile( isZero, self.terms )
        else:
            initialZero = []
        if len( self.terms ) > 1:
            nonzeroBlocks = []
            while self.terms != []:
                nonzeroBlocks.append( takeWhile( isNonzero, self.terms ) )
                self.terms = dropWhile( isNonzero, self.terms )
                self.terms = dropWhile( isZero, self.terms )
            self.terms += nonzeroBlocks[0]
            for aBlock in nonzeroBlocks[1:]:
                self.terms[-1] += aBlock[0]
                self.terms += aBlock[1:]
        self.terms = initialZero + self.terms

    def value( self ):
        if len( self.terms ) is 0:
            print "Warning: evaluation of continued fraction with empty list of terms.\n"
            return None
        elif len( self.terms ) is 1:
            return Rational( self.terms[ 0 ] )
        elif self.terms[1] is 0:
            print "Error: evaluation of continued fraction gives infinite value."
            return None
        else:
            eval = Rational( self.terms[ -1 ], 1 )
            for aTerm in self.terms[:-1:][::-1]:
                eval = eval.reciprocal() + aTerm
            return eval

    def __str__( self ):
        return '[ '  + ', '.join( [ str( n ) for n in self.terms ] ) + ' ]'

# functions to compute the continued fraction expansions
# of rational numbers

def intPart( r ):
    return r.num / r.denom

def findCFrac( r ):
    """Expand the rational number r as a continued fraction."""
    if r.denom == 1:
        return [r.num]
    else:
        expansion = []
        expansion.append( intPart ( r ) )
        remainder = r - intPart( r )
        while remainder.num != 1 and remainder.num != -1:
            remainder = remainder.reciprocal()
            expansion.append( intPart( remainder ) )
            remainder = remainder - intPart( remainder )
        expansion.append( remainder.num * remainder.denom )    
        return CFrac( expansion )

def evenPart( r ):
    return(  intPart( ( r+1 ) / 2  ) * 2 )
    
def findEvenCFrac( r ):
    """Expands the rational number r as a continued fraction of the forma
    [ 2a_1, 2b_1, ..., 2a_n, b_n ],
    where if b_n is 1 or -1 then it has the same sign as a_n."""

    expansion = []
    if r.denom == 1:
        expansion.append( r.num )
    else:
        expansion.append( evenPart ( r ) )
        remainder = r - evenPart( r )
        while remainder.num != 1 and remainder.num != -1:
            remainder = remainder.reciprocal()
            expansion.append( evenPart( remainder ) )
            remainder = remainder - evenPart( remainder )
        expansion.append( remainder.num * remainder.denom )

    if len(expansion) % 2 == 1 and expansion[-1] % 2 == 1:
        if expansion[-1] > 0:
            expansion = expansion[:-1] + [ expansion[-1] - 1,  1]
        else:
            expansion = expansion[:-1] + [ expansion[-1] + 1, -1]

    return CFrac( expansion )

# conversion formula between the Scharlemann-Thompson and principal slope
# invariants

def convertInvariants( p, q=1 ):
    """Compute the converted invariant."""
    r = Rational( p, q )
    expansion = findEvenCFrac( r ).terms
    aSum = sum( expansion[::2] )
    sign = 1
    if r.denom % 2 is 1:
        sign = -1
    return CFrac( [ sign * aSum ] + [ -k for k in expansion[1:][::-1] ] ).value()

def convert( p, q=1 ):
    """Print the converted invariant."""
    print convertInvariants( p, q )

def convertRange( q, m, n ):
    """Convert all invariants in a given range."""
    if m % 2 is 0:
        mm = m + 1
    else:
        mm = m
    for k in range( mm, n+1 )[::2]:
        print( str( Rational( k, q )) + ', ' + str( convertInvariants( k, q ) ) )

# now begin the slope calculations for semisimple tunnels of
# 2-bridge knots

def expandABlocks( frac ):
    """An auxiliary function for slopes( p, q ). It replaces each 2a_i
    term with [ 2, 0, 2, 0, 2, 0, ..., 0, 2 ] (or its negative, when
    a_i < 0), where there are a_i 2's."""
    newExp = []
    place = 0
    while place < len( frac ):
        if frac[place] > 2:
            newExp = newExp + [2, 0] * (frac[place]/2 - 1) + [2]
        elif frac[place] == 2 or frac[place] == -2:
            newExp.append( frac[place] )
        elif frac[place] < -2:
            newExp = newExp + [-2, 0] * ( abs(frac[place])/2 - 1) + [-2]
        newExp.append( frac[place + 1] )
        place += 2
    return newExp

def slopeExpression( k, parity ):
    """An auxiliary function for slopes( p, q )."""
    if parity % 2 is 0:
        return Rational( (2 * k) + 1, k )
    else:
        return Rational( (-2 * k) + 1, k )

def prettyPrintSlopes( slopeList ):
    slopeString = '[ ' + str( slopeList[0] ) + ' ]'
    if len( slopeList ) > 1:
        slopeString += ', ' + ', '.join( [ str(r) for r in slopeList[1:] ] )
    return slopeString

def slopeList( p, q ):
    q = q % p
    r = Rational( p, q )
    expansion = expandABlocks( findEvenCFrac( r ).terms )
    slopeList = []

    # calculate the first slope parameter using a_n and b_n
    if expansion[-2] is 2:
        slopeList.append( Rational( expansion[-1], (expansion[-1] * 2) + 1 ) )
    else:
        slopeList.append( Rational( expansion[-1] - 1, (expansion[-1] * 2) - 1 ) )

    # save b_n for parity calculations
    bn = expansion[-1]

    # for the remaining slopes, first calculate the k_i and the corresponding
    # parities, as a list [ (k_1, p_1), ..., (k_{n-1}, p_{n-1}) ]
    # these will be fed to the slopeExpression( k, p ) function
    
    kvaluesWithParities = []
    place = 0
    while place < len( expansion ) - 2:
        ai = expansion[place]
        bi = expansion[place+1]
        aip1 = expansion[place+2]

        if ai is 2 and aip1 is 2:
            kvaluesWithParities.append( ( bi + 1, bn + 1 ) )
        elif ai is 2 and aip1 is -2:
            kvaluesWithParities.append( ( bi, bn ) )
        elif ai is -2 and aip1 is 2:
            kvaluesWithParities.append( ( bi, bn + 1 ) )
        elif ai is -2 and aip1 is -2:
            kvaluesWithParities.append( ( bi - 1, bn ) )
        place += 2

    # compute the slopes using the reverse-ordered list of k-parity pairs
    # slopeList already contains the first slope
    slopeList += [ slopeExpression( k, parity )\
                   for (k, parity) in kvaluesWithParities[::-1] ]
    return slopeList

def slopes( p, q ):
    print prettyPrintSlopes( slopeList( p, q ) )

def slopesRange( b, m, n ):
    for a in range( m, n+1 ):
        print str(b) + '/' + str(a) + ' -> ' + prettyPrintSlopes( slopeList( b, a ) )