#!/usr/bin/python -i ######################################################################## # # Tunnels of torus knots # # by Sangbum Cho and Darryl McCullough # # Version of July 7, 2008 # # Contact: dmccullough@math.ou.edu # # This is a Python script that implements algorithms from the paper # # "Cabling sequences of torus knot tunnels". # # Also, some functions are included to calculate the depth invariant. # ######################################################################## # # To find the cabling slope sequence for the middle tunnel of # a (p,q) torus knot ( p > q ), use middleSlopes( p, q ) # # TorusKnots> middleSlopes( 41, 29 ) # [ 1/3 ], 5, 17, 29, 99, 169, 577 # # To find the sequence of binary invariants for the middle tunnel, # use binaries( p, q ) # # TorusKnots> binaries( 41, 29 ) # [1, 0, 1, 0, 1] # # To find the intermediate torus knots for the cabling sequence for the # middle tunnel, use intermediates( p, q ) # # TorusKnots> intermediates( 41, 29 ) # (3,2), (4,3), (7,5), (10,7), (17,12), (24,17), (41,29) # # To find the depth of the middle tunnel, use depth( p, q ) # # TorusKnots> depth( 41, 29 ) # 4 # # To find the depths of the middle tunnels for the (p, q) torus knots with # 1 < q < p, use depths p # # TorusKnots> depths( 41 ) # [1,1,1,1,1,1,1,2,1,2,3,2,1,2,3,3,3,2,1,1,2,3,3,3,3,2,3,4,3,2,3,2,2,2,2,2,2,2,1] # TorusKnots> depths( 42 ) # [-,-,-,2,-,-,-,-,-,2,-,2,-,-,-,2,-,2,-,-,-,2,-,2,-,-,-,3,-,3,-,-,-,-,-,3,-,-,-,1] # # (the "-" entries indicate q not relatively prime to 42) # # To find the slope sequence of the upper tunnel, use upperSlopes( p, q ) # # TorusKnots> upperSlopes( 18, 7 ) # [ 1/ 5 ], 11, 15, 21, 25, 31 # # TorusKnots> upperSlopes( 7, 18 ) # [ 1/ 3 ], 3, 3, 5, 5, 7, 7, 7, 9, 9, 11, 11, 11, 13, 13 # # and for the slopes of the lower tunnel, use lowerSlopes( p, q ) # # TorusKnots> lowerSlopes( 18, 7 ) # [ 1/ 3 ], 3, 3, 5, 5, 7, 7, 7, 9, 9, 11, 11, 11, 13, 13 # ######################################################################## import re, sys sys.ps1 = 'TorusKnots> ' 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." def __init__( self, num=0 , denom=1 ): self.num = num self.denom = denom d = gcd( num, denom ) if d != 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 if self.denom == 0: self.num = 1 def __str__( self ): if self.denom == 1: return "%d" % self.num else: return "%d/%d" % (self.num, self.denom) # define operations for r + s, where r is rational and s is # 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 reciprocal( self ): return( Rational( self.denom, self.num ) ) 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 ): returnList = [] for aTerm in lst: if f( aTerm ): returnList.append( aTerm ) else: return returnList return returnList def dropWhile( f, lst ): 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() if len(self.terms) is 0: print "Warning: Continued fraction initialized with empty terms list." 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} ].""" while len( self.terms ) > 1 and self.terms[-1] is 0: self.terms = self.terms[:-2] if len(self.terms) < 2: return True if self.terms[0] == 0: initialZero=[0] self.terms = dropWhile( isZero, self.terms ) else: initialZero = [] # stripping out zeros can create new zeros, as in # [ 2, 0, -2] --> [ 0 ], # so we need to use a while loop for the next step while len( self.terms ) > 1 and len( filter( isZero, self.terms )) > 0: 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 Rational( 1, 0 ) elif len( self.terms ) is 1: return Rational( self.terms[ 0 ] ) 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 ] ) + ' ]' # define functions to compute the continued fraction expansions of # rational numbers def intPart( r ): return r.num / r.denom def findCFrac( r ): 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 ) class Matrix: def __init__( self, a=0, b=0, c=0, d=0 ): self.a, self.b, self.c, self.d = a, b, c, d def __str__(self): return "[ [ %d, %d ], [ %d, %d ] ]" % (self.a, self.b, self.c, self.d) def __mul__(self,n): # defines self * n return Matrix( self.a * n.a + self.b * n.c, self.a * n.b + self.b * n.d, self.c * n.a + self.d * n.c, self.c * n.b + self.d * n.d ) # begin calculation of slope and depth invariants for the middle tunnels U = Matrix(1, 1, 0, 1) L = Matrix(1, 0, 1, 1) def prettyPrintSlopes( slopeList ): slopeString = '[ ' + str( slopeList[0] ) + ' ]' if len( slopeList ) > 1: slopeString += ', ' + ', '.join( [ str(r) for r in slopeList[1:] ] ) return slopeString def diagonals( m ): return m.a * m.d + m.b * m.c def middleSlopes( p, q ): pp, qq = abs(p), abs(q) # check for bad inputs if gcd( pp, qq ) is not 1: print 'p and q are not relatively prime.', \ 'Their greatest common divisor is d = %d.' % gcd( pp, qq ) print "Using (p/d, q/d) = ( %d, %d )." % ( p/gcd( pp, qq), q/gcd( pp, qq) ) pp, qq = pp/gcd( pp, qq ), qq/gcd( pp, qq ) if pp < 2 or qq < 2: print "The knot is trivial." return True if pp < qq: pp, qq = qq, pp expansion = findCFrac( Rational( pp, qq ) ).terms slopeList = [] n1 = expansion[0] slopeList.append( Rational( 1, 2 * n1 + 1 ) ) matrix = Matrix( n1 + 1, 1, n1 , 1 ) reducedExpansion = expansion[1:] reducedExpansion[0] = reducedExpansion[0] - 1 reducedExpansion[-1] = reducedExpansion[-1] - 1 parity = 0 for n in reducedExpansion: for k in range(0,n): if parity % 2 is 0: matrix = U * matrix else: matrix = L * matrix slopeList.append( diagonals( matrix ) ) parity = parity + 1 if p * q < 0: slopeList = [ slopeList[0] * (-1) + 1 ]\ + [ r * (-1) for r in slopeList[1:] ] print prettyPrintSlopes( slopeList ),'\n' def intermediates( p, q ): pp, qq = abs( p ), abs( q ) # check for bad inputs if gcd( pp, qq ) is not 1: print 'p and q are not relatively prime.', \ 'Their greatest common divisor is d = %d.' % gcd( pp, qq ) print "Using (p/d, q/d) = ( %d, %d )." % ( p/gcd( pp, qq), q/gcd( pp, qq) ) pp, qq = pp/gcd( pp, qq ), qq/gcd( pp, qq ) if pp < 2 or qq < 2: print "The knot is trivial." return True if pp < qq: pp, qq = qq, pp expansion = findCFrac( Rational( pp, qq ) ).terms n1 = expansion[0] matrix = Matrix( n1 + 1, 1, n1 , 1 ) intermediateList = [] intermediateList.append( ( 2 * n1 + 1, 2 ) ) reducedExpansion = expansion[1:] reducedExpansion[0] = reducedExpansion[0] - 1 reducedExpansion[-1] = reducedExpansion[-1] - 1 parity = 0 for n in reducedExpansion: for k in range(0,n): if parity % 2 is 0: matrix = U * matrix else: matrix = L * matrix intermediateList.append( ( matrix.a + matrix.c, matrix.b + matrix.d ) ) parity = parity + 1 if p * q < 0: intermediateList = [ (a, -b) for (a, b) in intermediateList ] print ', '.join( [ '('+str(a)+','+str(b)+')' for (a,b) in intermediateList ] ) def binaries( p, q ): pp, qq = abs( p ), abs( q ) # check for bad inputs if gcd( pp, qq ) is not 1: print 'p and q are not relatively prime.', \ 'Their greatest common divisor is d = %d.' % gcd( pp, qq ) print "Using (p/d, q/d) = ( %d, %d )." % ( p/gcd( pp, qq), q/gcd( pp, qq) ) pp, qq = pp/gcd( pp, qq ), qq/gcd( pp, qq ) if pp < 2 or qq < 2: print "The knot is trivial." return True if pp < qq: pp, qq = qq, pp expansion = findCFrac( Rational( pp, qq ) ).terms nsum = 0 for n in expansion[1:]: nsum += n; binariesList = [0] binariesList*= nsum + 1 place = 0 for n in expansion[1:]: place += n binariesList[place] = 1 return binariesList[2:-2] def depthCalc( p, q ): expansion = findCFrac( Rational( p, q ) ).terms depth = 0 tail = expansion[1:] place = 0 while place < len( tail ): if tail[place] is 1: if place < len(tail) - 1: depth += 1 place += 2 else: depth += 1 place += 1 return depth def depth( p, q ): pp, qq = abs( p ), abs( q ) if gcd( pp, qq ) is not 1: print 'p and q are not relatively prime.', \ 'Their greatest common divisor is d = %d.' % gcd( pp, qq ) print "Using (p/d, q/d) = ( %d, %d )." % ( p/gcd( pp, qq), q/gcd( pp, qq) ) pp, qq = pp/gcd( pp, qq ), qq/gcd( pp, qq ) if pp < qq: pp, qq = qq, pp if pp < 2 or qq < 2: print 0 else: print depthCalc( pp, qq ) def depths( n ): depthsList = [] for k in range(2,n): if gcd( n, k ) is not 1: depthsList.append( '-' ) else: depthsList.append( str(depthCalc( n, k )) ) print '[' + ','.join( depthsList ) + ']' # begin calculation of slope invariants for the upper and lower tunnels def pk( p, q , k): if ( k * p ) % q == 0: return k * p / q else: return 1 + k * p / q def upperSlopes( p, q ): """Compute and print the slopes of the upper tunnel of the (p, q) torus knot.""" pp, qq = abs(p), abs(q) # check for bad inputs if gcd( pp, qq ) is not 1: print 'p and q are not relatively prime.', \ 'Their greatest common divisor is d = %d.' % gcd( pp, qq ) print "Using (p/d, q/d) = ( %d, %d )." % ( p/gcd( pp, qq), q/gcd( pp, qq) ) pp, qq = pp/gcd( pp, qq ), qq/gcd( pp, qq ) if pp < 2 or qq < 2: print "The knot is trivial." return True # inputs are good preSlopes = [ 2 * pk(pp, qq, k) - 1 for k in \ range(1,qq) if pk( pp, qq, k) > 1 ] slopes = [ Rational( 1, preSlopes[0] ) ] + \ [ Rational( k ) for k in preSlopes[1:] ] if p*q < 0: slopes = (-1) * slopes print prettyPrintSlopes( slopes ) def lowerSlopes( p, q ): pp, qq = abs(p), abs(q) # check for bad inputs if gcd( pp, qq ) is not 1: print 'p and q are not relatively prime.', \ 'Their greatest common divisor is d = %d.' % gcd( pp, qq ) print "Using (p/d, q/d) = ( %d, %d )." % ( p/gcd( pp, qq), q/gcd( pp, qq) ) pp, qq = pp/gcd( pp, qq ), qq/gcd( pp, qq ) if pp < 2 or qq < 2: print "The knot is trivial." return True upperSlopes( qq, pp)