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