#!/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 ) )