#!/usr/bin/python -i

########################################################################
#
#  Software for "Semisimple tunnels"
#
#     by Sangbum Cho and Darryl McCullough
#
#  Version of June 21, 2010
#
#  Contact:  dmccullough@math.ou.edu
#
#  This is a Python script, written for Python 2.3
#
########################################################################
#  
#  To find the cabling slope sequence for the upper tunnel of the (1,1)-knot
#  determined by the braid word  m^3 s^{-2} l^4 s^{-4} m^{-1} s^{-4} l^3 m s l^2
#  use
#
#  Semisimple> upperSlopes( 'm 3 s -2 l 4 s -4 m -1 s -4 l 3 m 1 s 1 l 2' )
#  [ 5/9 ], 131/78, -11, -11
#
#  and for the lower slopes,
#
#  Semisimple> lowerSlopes( 'm 3 s -2 l 4 s -4 m -1 s -4 l 3 m 1 s 1 l 2' )
#  [ 16/19 ], -7, -7, -7, -195/31
#
#  You can find a braid word given the slope sequence:
#  Semisimple> braidWord( [ 5, 9, 131, 78, -11, 1, -11, 1 ] )
#  <__main__.Word instance at 0x2b364bee7a28>
#  
#  Semisimple> print braidWord( [ 5, 9, 131, 78, -11, 1, -11, 1 ] )
#  m 3 s -3 m -1 l -3 m 1 l -1 s -4 m 1 s -4 l -1
#
#  Semisimple> upperSlopes( 'm 3 s -3 m -1 l -3 m 1 l -1 s -4 m 1 s -4 l -1' )
#  [ 5/9 ], 131/78, -11, -11
#  
#  These functions also accept the actual word type:
#
#  Semisimple> upperSlopes( braidWord( [ 5, 9, 131, 78, -11, 1, -11, 1 ] ) )
#  [ 5/9 ], 131/78, -11, -11
#
#  For the slopes of the dual tunnel of a tunnel with a given slope sequence:
#
#  Semisimple> dualSlopes( [ 5, 9, 131, 78, -11, 1, -11, 1 ] )
#  [ 16/19 ], -7, -7, -7, -195/31
#
########################################################################
#
#  To find the slope sequences for the tunnels of the two-bridge knot
#  with rational invariant a/b, use twoBridge(b, a)  (b  must be odd):
#
#  Semisimple> twoBridge( 413, 290 )
#  Upper simple tunnel:     [ 47/413 ]
#  Upper semisimple tunnel: [ 1/3 ], 3, 3, 3, 7/4, -1, -3
#  Lower simple tunnel:     [ 290/413 ]
#  Lower semisimple tunnel: [ 2/3 ], -7/3, -1, -3/2, 3, 3, 3
#
#  Semisimple> print dualInvariant( 413, 290 )
#  413/47
#
#  Semisimple> twoBridge( 413, 47 )
#  Upper simple tunnel:     [ 290/413 ]
#  Upper semisimple tunnel: [ 2/3 ], -7/3, -1, -3/2, 3, 3, 3
#  Lower simple tunnel:     [ 47/413 ]
#  Lower semisimple tunnel: [ 1/3 ], 3, 3, 3, 7/4, -1, -3
#
#  There are also functions to find the slope sequence of any one of the
#  tunnels:
#
#  Semisimple> upperSemisimpleSlopes( 413, 290 )
#  [ 1/3 ], 3, 3, 3, 7/4, -1, -3
#
#  The actual braid words used for the calculations can be obtained.
#  In each case, the given tunnel is the upper tunnel for the
#  (1,1)-position:
#
#  Semisimple> print upperSemisimpleBraidWord( 413, 290 )
#  m -1 l 1 m -1 l -1 m -1 s -4 l 1 m 4 l -2
#  Semisimple> print upperSimpleBraidWord( 413, 290 )
#  m -1 s -1 l 2 m -1 l -1 m 1 l 1 m -1 l 2 m 1 l -5
#  Semisimple> print lowerSimpleBraidWord( 413, 290 )
#  m -2 l 4 m 1 s -4 l -1 m -1 l -1 m 1 l -1
#  Semisimple> print lowerSemisimpleBraidWord( 413, 290 )
#  m -5 l 1 m 2 l -1 m 1 l 1 m -1 l -1 m 2 s -1 l -1
#
#########################################################################
#
#  To find the slope sequences of the upper and lower tunnels
#  of a torus knot, use
#
#  Semisimple> torusUpperSlopes( 24, 5 )
#  [ 1/9 ], 19, 29, 39
#  Semisimple> torusLowerSlopes( 24, 5 )
#  [ 1/3 ], 3, 3, 3, 3, 5, 5, 5, 5, 5, 7, 7, 7, 7, 7, 9, 9, 9, 9
#
#  And the braid word is given by
#
#  Semisimple> print torusBraidWord( 27, 5 )
#  m 1 l -5 m 1 l -6 m 1 l -5 m 1 l -6
#
#  or for the unreduced braid word
#
#  Semisimple> print fullTorusBraidWord( 27, 5 )
#  l -5 m 1 l -5 m 1 l -6 m 1 l -5 m 1 l -6 m 1
#
########################################################################

import sys
sys.ps1 = 'Semisimple> '
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)

    def __cmp__( self, r ):
        if self.num * r.denom < r.num * self.denom:
            return -1
        elif self.num * r.denom == r.num * self.denom:
            return 0
        else:
            return 1
        
    # 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 ) )

    def isPositive( self ):
        return self.num * self.denom > 0

#  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()

    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:
            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 isInfinite( self ):
        return self.value().denom == 0

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

class Letter:
    """Class of letters.
    A letter has two variables, the lettertype--- one of 'm', 'l',
    's', or 'ms'--- and the exponent, an integer."""
    
    def __init__( self, lettertype = 'm', exponent = 0 ):
        self.lettertype = lettertype
        self.exponent = int( exponent )
        if not self.lettertype in [ 'm', 'l', 's', 'ms' ]:
           print "Error creating Letter: lettertype not 'm', 'l', 's', or 'ms'"
           sys.exit()

    def __mul__( self, y ):
        if not self.lettertype == y.lettertype:
            print 'Error: undefined multiplication of letters ', self, ' * ', y
            sys.exit()
        return Letter( self.lettertype, self.exponent + y.exponent )

    def __str__( self ):
        return self.lettertype + ' ' + str( self.exponent )

class Word:

    def __init__( self, letterList = [ ] ):
        self.letterList = letterList
    
    def stripZeros( self ):
        for letter in self.letterList:
            if letter.exponent == 0:
                self.letterList.remove( letter )

    def stripLeftCoset( self ):
        while len( self.letterList ) > 0 and \
                  self.letterList[0].lettertype in ['l', 's']:
            self.letterList = self.letterList[1:]

    def stripRightCoset( self ):
        while len( self.letterList ) > 0 and \
                  self.letterList[-1].lettertype in ['m', 's', 'ms' ]:
            self.letterList = self.letterList[:-1]

    def freelyReduce( self ):
        # Letters of the form  ms 1  are unaffected.
        self.stripZeros()
        if len( self.letterList ) > 1:
            newLetterList = [ ]
            for letter in self.letterList[0:]:
                if len( newLetterList ) == 0:
                    newLetterList.append( letter )
                elif letter.lettertype == "ms":
                    newLetterList.append( letter )
                elif newLetterList[-1].lettertype == letter.lettertype:
                    newLetterList[-1] = newLetterList[-1] * letter
                    if newLetterList[-1].exponent == 0:
                        del newLetterList[-1]
                else:
                    newLetterList.append( letter )
            self.letterList = newLetterList

    def reduce( self ):
        self.freelyReduce()
        self.stripLeftCoset()
        self.stripRightCoset()

    def __mul__( self, y ):
        return Word( self.letterList + y.letterList )

    def __str__( self ):
        returnString = ''
        return ' '.join( [ str( letter ) for letter in self.letterList ] )

    def dual( self ):
        newletterList = [ ]
        for aLetter in self.letterList:
            if aLetter.lettertype == 'l':
                newletterList.append( Letter( 'm', aLetter.exponent ) )
            elif aLetter.lettertype == 'm':
                newletterList.append( Letter( 'l', aLetter.exponent ) )
            elif aLetter.lettertype == 'ms':
                # since (ms)^2 = 1, ignore if exponent is even, else treat as 1
                if aLetter.exponent % 2 == 1:
                    newletterList.append( Letter( 'l', 1 ) )
                    newletterList.append( Letter( 's', 1 ) )
            elif aLetter.lettertype == 's':
                newletterList.append( Letter( 's', aLetter.exponent ) )
        newletterList.reverse()
        returnWord = Word( newletterList )
        returnWord.reduce()
        return returnWord

def prettyWord( aWord ):
    # replace each \sigma^2 or \sigma^P{-2} by a commutator to
    # get a prettier word
    newWordList = [ ]
    commutator = [ Letter( 'l', 1), Letter( 'm', -1 ), \
                   Letter( 'l', -1), Letter( 'm', 1 ) ]
    inverseCommutator = [ Letter( 'm', -1), Letter( 'l', 1 ), \
                         Letter( 'm', 1), Letter( 'l', -1 ) ]

    for aLetter in aWord.letterList:
        if aLetter.lettertype == 's' and aLetter.exponent == -2:
            newWordList = newWordList + commutator
        elif aLetter.lettertype == 's' and aLetter.exponent == 2:
            newWordList = newWordList + inverseCommutator
        else:
            newWordList.append( aLetter )
    returnWord = Word( newWordList )
    returnWord.reduce()
    return returnWord

def listToWord( tokenList ):
    """Make a word from a list
    [ 'x_1', n_1, 'x_2', ..., 'x_m', n_m], where each x_i is one of
    the letters m, l, or s, and each n_i is an integer.
    [ 'm', 2, 's', -2, 'l', 3 ] represents the word m^2 s^{-2} l^3.
    """
    if len( tokenList ) % 2 == 1:
        print 'Error forming word: input list must have even length'
        return None
    for position in range(0, len( tokenList ) - 1 , 2):
        if not tokenList[position] in [ 'm', 'l', 's', 'ms' ]:
            print "Error forming word: odd terms must be 'm', 'l', 's', or 'ms'"
            return None
        try:
            exponent = int( tokenList[ position + 1 ] )
        except:
            print 'Error forming word: even terms must be integers.'
            return None
    letterList = [ ]
    for count in range(0, len( tokenList )/2 ):
        letterList.append( Letter( tokenList[2 * count],
                                tokenList[1 + 2 * count ] ) )
    return Word( letterList )

def stringToWord( inputString ):
    return listToWord( inputString.split() )

def wordToString( word ):
    returnString = ''
    for letter in word.letterList:
        returnString = returnString + letter.lettertype + ' '
        returnString = returnString + str(letter.exponent) + ' '
    return returnString[:-1]

def slopeForm( word ):
    word.reduce()

    # make all appearances of m have positive exponent
    mPositiveLetterList = [ ]
    for thisLetter in word.letterList:
        if thisLetter.lettertype == 'm' and thisLetter.exponent < 0:
            for count in range(0, -thisLetter.exponent):
                mPositiveLetterList.append( Letter( 's', 1 ) )
                mPositiveLetterList.append( Letter( 'm', 1 ) )
                mPositiveLetterList.append( Letter( 's', 1 ) )
        else:
            mPositiveLetterList.append( thisLetter )

    # insert ms letters corresponding to cablings
    msLetterList = [ ]
    for thisLetter in mPositiveLetterList:
        if thisLetter.lettertype == 'm':
            for count in range(0, thisLetter.exponent ):
                msLetterList.append( Letter( 'ms', 1 ) )
                msLetterList.append( Letter( 's', -1 ) )
        else:
            msLetterList.append( thisLetter )

    msWord = Word( msLetterList )
    msWord.stripLeftCoset()
    msWord.freelyReduce()

    # cancel ms ms pairs 
    returnLetterList = [ ]
    for thisLetter in msWord.letterList:
        if len( returnLetterList ) == 0:
            returnLetterList.append( thisLetter )
        elif returnLetterList[-1].lettertype == "ms" and thisLetter.lettertype == "ms":
            del returnLetterList[-1]
        else:
            returnLetterList.append( thisLetter )
    returnWord = Word( returnLetterList )
    returnWord.reduce()

    return returnWord
    
def t( word ):
    letterList = word.letterList
    total = 0
    parity = 1
    for letter in letterList:
        if letter.lettertype == 'ms' or letter.lettertype == 's':
            parity = (parity + letter.exponent) % 2
        if letter.lettertype == 'l':
            total = total + (-1)**parity * letter.exponent
    return total

def padSlopeWord( word ):
    """Operates on a word having only 'l' and 's' letters.
    Assuming that it is reduced, tack on 's 0' in front if the word
    starts with an 'l' letter, and tack on 'l 0' at the end if it
    ends in an 's' letter."""
    
    if len( word.letterList ) == 0:
        return word
    newLetterList = [ ]
    if word.letterList[0].lettertype == 'l':
        newLetterList.append( Letter( 's', 0 ) )
    newLetterList = newLetterList + word.letterList
    if word.letterList[-1].lettertype == 's':
        newLetterList.append( Letter( 'l', 0 ) )
    return Word( newLetterList )

def findSlopeWords( word ):
    slopeFormWord = slopeForm( word )
    if len( slopeFormWord.letterList ) < 2:
        return [ ]

    cablingWordList = [ ]
    thisCabling = [ ]
    for letter in slopeFormWord.letterList[1:]:
        if not letter.lettertype == "ms":
            thisCabling.append( letter )
        else:
            thisCablingWord = Word( thisCabling )
            thisCablingWord.freelyReduce()
            thisCablingWord = padSlopeWord( thisCablingWord )
            cablingWordList.append( thisCablingWord )
            thisCabling = [ ]

    # when the end of the word has been reached, make the
    # last block into a slope word
    if len( thisCabling ) > 0:
        thisCablingWord = Word( thisCabling )
        thisCablingWord.freelyReduce()
        thisCablingWord = padSlopeWord( thisCablingWord )
        cablingWordList.append( thisCablingWord ) 

    return cablingWordList

def slopeWordToCFrac( word ):
    exponentList = [ letter.exponent for letter in word.letterList ]
    exponentList.reverse()
    for position in range(0, len(exponentList), 2):
        exponentList[ position ] = -exponentList[ position ] * 2
    return CFrac( exponentList )

def omitInfiniteSlopes( slopeWords ):
    cleanPassMade = False
    while not cleanPassMade:
        cleanPassMade = True

        finiteSlopes = [ ]
        infiniteSlopeWord = None
        remainingSlopes = [ ]

        # Look for an infinite slope. If one is found, split the list
        # at that point
        for s in slopeWords:
            if slopeWordToCFrac( s ).isInfinite():
                cleanPassMade = False
                finiteSlopes = slopeWords[ :slopeWords.index(s) ]
                infiniteSlopeWord = s
                remainingSlopes = slopeWords[ slopeWords.index(s)+1: ]
                break

        if not cleanPassMade:
            # if there was only one slope, it is deleted
            # and the process finishes
            if len(slopeWords) == 1:
                slopeWords = [ ]
                cleanPassMade = True
                
            # if the first slope is infinite, delete the first two
            elif len( finiteSlopes ) == 0:
                slopeWords = slopeWords[2:]
                
            else:
                # first, move the powers of l to the previous word
                # an infinite slope word w is equivalent to l^{-t(w)}
                # pushing it past ms makes it l^{t(w)}
                previousWord = finiteSlopes[-1]
                exp = t(infiniteSlopeWord )
                previousWord = previousWord * listToWord( ['l', exp] )

                # if w was not the last word, must combine the last
                # finite word with the next word, since
                # the ms * ms disappears
                if len( remainingSlopes ) > 0:
                    previousWord = previousWord * remainingSlopes[0]

                previousWord.freelyReduce()
                finiteSlopes[-1] = padSlopeWord( previousWord )
                slopeWords = finiteSlopes + remainingSlopes[1:]

    return slopeWords

def calculateSlopes( word ):
    slopeWords = findSlopeWords( word )
    slopeWords = omitInfiniteSlopes( slopeWords )
    slopes = [ ]

    # accumulate words in w2, to calculate algebraic winding numbers
    w2 = Word( [ ] )

    slopeWords.reverse()
    for thisWord in slopeWords:

        exponentList = [ letter.exponent for letter in thisWord.letterList ]
        exponentList.reverse()
        for position in range(0, len(exponentList), 2):
            exponentList[ position ] = -exponentList[ position ] * 2
        exponentList[0] += 2 * t( w2 )

        f = CFrac( exponentList )
        slopes.append( f.value() )

        w2 = (listToWord( [ "ms", 1 ] ) * thisWord) * w2

    while len(slopes) > 0 and abs(slopes[0].num) == 1:
        del slopes[0]
    if len(slopes) > 0:
        slopes[0] = slopes[0].reciprocal() - intPart( slopes[0].reciprocal() )
    return slopes

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

def upperSlopes( input ):
    if isinstance( input, Word ):
        inputWord = input
    else:
        inputWord = stringToWord( input )
    print prettyPrintSlopes( calculateSlopes( inputWord ) )

def lowerSlopes( input ):
    if isinstance( input, Word ):
        inputWord = input
    else:
        inputWord = stringToWord( input )
    inputWord = inputWord.dual()
    print prettyPrintSlopes( calculateSlopes( inputWord ) )

def dualSlopes( input ):

    if isinstance( input, Word ):
        aBraidWord = input
    elif isinstance( input, str ):
        aBraidWord = stringToWord( input )
    elif isinstance( input, list ):
        aBraidWord = braidWord( input )
    else:
        print 'Input type not recognized.'
        return None

    dualWord = aBraidWord.dual()
    print prettyPrintSlopes( calculateSlopes( dualWord ) )

########################################################################
#
#  Now we carry out the reverse construction: start with the slope sequence
#  and work out a braid word description of the tunnel.
#

#  functions to compute continued fraction expansions of rational numbers

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

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 )

def slopeToWord( cFracTerms ):
    """Auxiliary function for braidWord()."""

    mSigma = listToWord( ['m', 1 ] ) * listToWord( ['s', 1 ] )

    cFracTerms.reverse()
    returnWord = Word()
    for place in range(0, len(cFracTerms)/2 ):
        returnWord = returnWord * listToWord( ['s', cFracTerms[ 2 * place ]] )
        returnWord = returnWord * listToWord( ['l', -cFracTerms[ 2 * place + 1]/2] )  
    return mSigma * returnWord    

def braidWord( inputList ):
    
    if len( inputList ) == 0:
        return Word()
    if len( inputList ) % 2 == 1:
        print 'Input list must have an even number of terms.'
        return None
    for term in inputList:
        try:
            thisTerm = int( term )
        except:
            print 'Input list must have integer terms.'
            return None
    if inputList[1] % 2 == 0:
        print 'Second entry must be odd.'
        return None
    for place in range(2, len( inputList ), 2):
        if inputList[ place ] % 2 == 0:
            print 'Third, fifth, seventh, ... entries must be odd.'

    slopeList = [ ]
    for place in range(0, len(inputList)/2 ):
        slopeList.append( Rational( inputList[ 2 * place ],\
                                    inputList[ 2 * place + 1 ] ) )
    slopeList[0] = slopeList[0].reciprocal()
    
    expansions = [ findEvenCFrac( r ).terms for r in slopeList ]
    wordList = [ slopeToWord( exp ) for exp in expansions ]

    returnWord = wordList[0]
    for word in wordList[1:]:
        windingNumber = t( returnWord )
        returnWord = word * listToWord( [ 'l', windingNumber ] ) * returnWord

    returnWord = prettyWord( returnWord )
    return returnWord

##################################################################################
#
#  To compute slope sequence of tunnels of 2-bridge knots, we work out the
#  braid words for the two (1,1)-positions, then read off the slopes using
#  calculateSlopes()
#

def dualInvariant( p, q ):
    """Return the rational invariant for the inverted position."""
    if p % 2 == 0:
        print 'p must be odd for a 2-bridge knot.'
        return None

    q = q % p
    r = Rational( p, q )
    fracTerms = findEvenCFrac( r ).terms
    fracTerms = [ -1 * n for n in fracTerms ]
    fracTerms.reverse()
    invariant = CFrac( fracTerms ).value()
    if invariant.num < 0:
        r, s = -invariant.num, -invariant.denom
        s = s % p
        invariant = Rational( r, s )
    return invariant

def upperSemisimpleBraidWord( p, q ):
    """Find the braid word whose upper tunnel is the upper 
    semisimple tunnel of a 2-bridge knot."""
    if p % 2 == 0:
        print 'p must be odd for a 2-bridge knot.'
        return None
    q = q % p
    r = Rational( p, q )
    frac = findEvenCFrac( r ).terms
    letterList = [ ]
    for position in range(0, len( frac )/2):
        letterList.append( Letter( 'm', -frac[ 2 * position ]/2 ) )
        letterList.append( Letter( 's', frac[ 2 * position + 1 ] ) )
    letterList.append( Letter( 'l', -1 ) )
    return prettyWord ( Word( letterList ) )

def upperSimpleBraidWord( p, q ):
    """Find the braid word whose upper tunnel is the lower
    simple tunnel of a 2-bridge knot."""
    if p % 2 == 0:
        print 'p must be odd for a 2-bridge knot.'
        return None
    dualSemisimpleBraidWord = upperSemisimpleBraidWord( \
        dualInvariant(p,q).num, dualInvariant(p,q).denom )
    return dualSemisimpleBraidWord.dual()

def lowerSemisimpleBraidWord( p, q ):
    return upperSimpleBraidWord( p, q ).dual()

def lowerSimpleBraidWord( p, q ):
    return upperSemisimpleBraidWord( p, q ).dual()

def upperSemisimpleSlopes( p, q ):
    print prettyPrintSlopes( calculateSlopes( upperSemisimpleBraidWord( p, q ) ) )

def upperSimpleSlopes( p, q ):
    print prettyPrintSlopes( calculateSlopes( upperSimpleBraidWord( p, q ) ) )    

def lowerSemisimpleSlopes( p, q ):
    print prettyPrintSlopes( calculateSlopes( lowerSemisimpleBraidWord( p, q ) ) )    

def lowerSimpleSlopes( p, q ):
    print prettyPrintSlopes( calculateSlopes( lowerSimpleBraidWord( p, q ) ) )    

def twoBridge( p, q ):
    if p % 2 == 0:
        print 'p must be odd for a 2-bridge knot.'
        return None
    print 'Upper simple tunnel:    ', prettyPrintSlopes( \
        calculateSlopes( upperSimpleBraidWord( p, q ) ) )
    print 'Upper semisimple tunnel:', prettyPrintSlopes( \
        calculateSlopes( upperSemisimpleBraidWord( p, q ) ) )
    print 'Lower simple tunnel:    ', prettyPrintSlopes( \
        calculateSlopes( lowerSimpleBraidWord( p, q ) ) )
    print 'Lower semisimple tunnel:', prettyPrintSlopes( \
        calculateSlopes( lowerSemisimpleBraidWord( p, q ) ) )

##################################################################################
#
#  For semisimple tunnels of torus knots, we compute a braid word,
#  then apply the upper and lower slope functions.
#

def fullTorusBraidWord( p, q ):
    """The unreduced braid word for 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 None

    letterList = [ ]

    if p * q > 0:
        p, q = pp, qq
        pkList = [ 0 ] + [ 1 + (k * p)/q for k in range(1,q) ] + [ p ]

        deltaM = Letter( 'm', 1 )
        for index in range(q-1, -1, -1):
            letterList.append( Letter( 'l', pkList[ index ] - pkList[ index + 1 ] ) )
            letterList.append( deltaM )

    if p * q < 0:
        p, q = pp, -qq
        qkList = [ 0 ] + [ 1 + (k * q)/p for k in range(1, p) ] + [ q ]

        deltaL = Letter( 'l', 1 )
        for index in range(0, p):
            letterList.append( deltaL )
            letterList.append( Letter( 'm', qkList[ index ] - qkList[ index + 1 ] ) )

    returnWord = Word( letterList )
    returnWord.freelyReduce()
    return returnWord

def torusBraidWord( p, q ):
    """The reduced braid word for the (p,q) torus knot."""
    returnWord = fullTorusBraidWord( p, q )
    returnWord.reduce()
    return returnWord

def torusUpperSlopes( p, q ):
    # check for bad inputs
    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) )
        p, q = p/gcd( pp, qq ), q/gcd( pp, qq )
    if pp < 2 or qq < 2:
        print 'The knot is trivial.'
        return None

    return upperSlopes( torusBraidWord( p, q ) )

def torusLowerSlopes( p, q ):
    # check for bad inputs
    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) )
        p, q = p/gcd( pp, qq ), q/gcd( pp, qq )
    if pp < 2 or qq < 2:
        print 'The knot is trivial.'
        return None

    return lowerSlopes( torusBraidWord( p, q ) )

########################################################################
#
#   Find the 2-bridge knot producing a given 2-bridge slope sequence
#

def kPart( ratl ):
    # finds k in expressions 2 + 1/k or -2 + 1/k
    if ratl.isPositive():
        reciprocalK = ratl - Rational( 2 )
    else:
        reciprocalK = ratl - Rational( -2 )

    if reciprocalK.num == 1:
        return reciprocalK.denom
    elif reciprocalK.num == -1:
        return -reciprocalK.denom
    else:
        return False

def constructSlopelist( inputList ):
    # convert an even-length integer list into a list of rational numbers
    returnList = [ ]
    for index in range( 0, len(inputList), 2):
        returnList.append( Rational( inputList[index], inputList[index + 1] ) )
    return returnList

def findKValues( slopeList ):
    returnList = [ ]
    for x in slopeList[1:]:
        returnList.append( kPart( x ) )
    return returnList

def computeCFrac( inputList ):
    # reconstruct a continued fraction expression for the rational invariant
    # a/b, given the list of slopes (as an even-length list of integers)
    
    slopeList = constructSlopelist( inputList )

    firstSlope = slopeList[0]
    firstSlope = Rational( firstSlope.num % firstSlope.denom, firstSlope.denom )
    if firstSlope.denom == 2 * firstSlope.num + 1:
        n0 = firstSlope.num
    else:
        n0 = -firstSlope.num
    
    kValues = findKValues( slopeList )

    # the list of a-values will be built in reverse order
    aList = [ ]

    if n0 % 2 == 0:
        aList.append( 2 )
    else:
        aList.append( -2 )

    for index in range(0, len(kValues)):
        if kValues[index] % 2 == 1:
            aList.append( aList[-1] )
        else:
            aList.append( (-1) * aList[-1] )

    # the continued-fraction list will also be built in reverse order    
    returnList = [ ]
    if n0 % 2 == 1:
        returnList.append( n0 + 1)
    else:
        returnList.append( n0 )
    returnList.append( aList[0] )

    for index in range(0, len( kValues )):
        if aList[index + 1] == aList[index]:
            if aList[index + 1] == 2:
                returnList.append( kValues[index] - 1 )
            else:
                returnList.append( kValues[index] + 1 )
        else:
            returnList.append( kValues[index] )
        returnList.append( aList[index + 1] )

    returnList.reverse()
    return CFrac( returnList )

def find2BridgeInvariant( inputList ):
    cFrac = computeCFrac( inputList )
    invariant = cFrac.value()
    if not invariant.isPositive():
        invariant = Rational( invariant.num, invariant.denom + invariant.num )
    return invariant

#  for the main calling program, we will need to test input data

def test2BridgeSlopeList( inputList ):
    if len( inputList ) == 0:
        print 'Empty input list.'
        return False
    if len( inputList ) % 2 == 1:
        print 'Input list must have even length.'
        return False
    for x in inputList:
        if not isinstance( x, int ):
            print 'Input list must have integer entries.'
            return False
    if inputList[1] % 2 == 0:
        print 'Second entry must be odd.'
        return False

    if len( inputList ) == 2:
        return True
    else:
        for index in range(2, len(inputList), 2):
            if inputList[ index ] % 2 == 0:
                print 'Third, fifth, ... entries must be odd.'
                return False

    slopeList = constructSlopelist( inputList )

    firstSlope = slopeList[0]
    firstSlope = Rational( firstSlope.num % firstSlope.denom, firstSlope.denom )
    if firstSlope.denom == 2 * firstSlope.num + 1:
        n0 = firstSlope.num
    elif firstSlope.denom == 2 * firstSlope.num - 1:
        n0 = -firstSlope.num
    else:
        print 'First slope must be equivalent in Q/Z to [n0/(2 * n0 + 1)],\n' \
              'with n0 not -1 or 0.'
        return False
    if firstSlope.denom == 1 or n0 == 0:
        print 'First slope must be equivalent in Q/Z to [n0/(2 * n0 + 1)],\n' \
              'with n0 not -1 or 0.'
        return False

    if len( slopeList ) > 1:
        if slopeList[1].num > 0 and n0 % 2 == 0 or \
           slopeList[1].num < 0 and n0 % 2 == 1:
            print 'm1 must be positive or negative according as n0 is odd or even.'
            return False
        for index in range(1, len( slopeList )):
            if not kPart( slopeList[index] ):
                print 'Slopes other than first must be of the form 2 + 1/k or 2 - 1/k.'
                return False
        for index in range(1, len(slopeList) - 1):
            if kPart( slopeList[index] ) % 2 == 0:
                if (slopeList[index] * slopeList[ index + 1 ]).isPositive():
                    print 'The ith and (i+1)st slopes must have opposite signs'
                    print 'when k_i is even.'
                    return False
            elif kPart( slopeList[index] ) % 2 == 1:
                if not (slopeList[index] * slopeList[ index + 1 ]).isPositive():
                    print 'The ith and (i+1)st slopes must have the same sign'
                    print 'when k_i is odd.'
                    return False

    return True

def find2BridgeKnot( inputList ):
    if test2BridgeSlopeList( inputList ):
        if len( inputList ) == 2:
            b, a = inputList[0], inputList[1]
            if a < 0:
                a, b = -a, -b
            b = b % a
            invariant = Rational( a, b )
            if not invariant.isPositive():
                invariant = Rational( invariant.num, invariant.denom + invariant.num )
            dual = dualInvariant( invariant.num, invariant.denom )
            print 'The slope sequence has length 1.'
            print 'The tunnel is the lower simple tunnel of K( %d, %d) , or' \
                % ( invariant.num, invariant.denom )
            print 'equivalently the upper simple tunnel of K( %d, %d ).' \
                  % ( dual.num, dual.denom )
        else:
            a, b = find2BridgeInvariant( inputList ).num, \
                   find2BridgeInvariant( inputList ).denom
            dual = dualInvariant( a, b )
            print 'The tunnel is the upper semisimple tunnel of K( %d, %d ), or' \
                  % (a, b)
            print 'equivalently the lower semisimple tunnel of K( %d, %d).' \
                  % ( dual.num, dual.denom )

########################################################################

testlist = [ 6, 13, -11, 5, -11, 6, 25, 13 ]
slopeList = [ 2, 3, -3, 1, -3, 2, 3, 1 ]

def slopeList2IntegerList( slopeList ):
    returnList = [ ]
    for r in slopeList:
        returnList.append( r.num )
        returnList.append( r.denom )
    return returnList

def iSlopes( p, q ):
    return slopeList2IntegerList( calculateSlopes( \
        upperSemisimpleBraidWord( p, q ) ) )

def testTest( p, q ):
    if test2BridgeSlopeList( iSlopes( p, q ) ):
        print 'K( %d, %d ) is OK.' % ( p, q )
    else:
        print 'K( %d, %d ) fails test.' % ( p, q )
        
def testTestRange( p, q, r ):
    for k in range( q, r, 2):
        testTest( p, k )

def testAndFind2BridgeInvariant( inputList ):
    terms = computeCFrac( inputList ).terms
    if test2BridgeSlopeList( terms ):
        invariant = cFrac.value()
        if not invariant.isPositive():
            invariant = Rational( invariant.num, invariant.denom + invariant.num )
        return invariant
    return 'No invariant found.'
    
def test( p, q ):
    if Rational( p, q) == find2BridgeInvariant( iSlopes( p, q ) ):
        print 'K( %d, %d ) OK.' % ( p, q )
    else:
        print 'K( %d, %d ) fails test.' % ( p, q )

def testRange( p, q, r ):
    for k in range( q, r, 2):
        test( p, k )
