# Boundary Slopes of 2-bridge links
# Jim Hoste,  jhoste@pitzer.edu
# and 
# Patrick D.  Shanahan,  pshanahan@lmu.edu
# May, 2005

# This program will compute the boundary slopes of 2-bridge links.
# All incompressible, boundary incompressible, merdionally incompressible surfaces in the
# complement of a 2-bridge link are classified, up to isotopy, but Floyd and Hatcher in their paper,
# "The space of incompressible surfaces in a 2-bridge link complement", TAMS Vol 305, No 2, 1988 575-599. 
# The  algorithm used here for computing the boundary slopes of these surfaces is explained in the paper 
# "Computing Boundary Slopes of 2-bridge Links"
# by Jim Hoste and Patrick D. Shanahan, to appear in Math. Comp.

# A 2-bridge link is determined by the pair of relatively prime integers p,q
# where 0<p<q, and q is even.

# Each essential surface in the link complement corresponds to a minimal path in a graph Dt, which specilizes
# in different ways to the graphs D0 and D1. See the above papers for defintions of these graphs.
# Given p/q, all minimal paths from 1/0 to p/q in the two graphs D1 and D0 must be found.
# Each minimal path in Dt from 1/0 to p/q corresponds to a pair of compatible minimalpath in D1 and D0.
# A branched surface corresponds to each minimal path in Dt. Each branched surface carries 
# essential surfaces whose boundary slopes we then compute.

# The routine bdyslopes(p,q) will take p and q as input and return a list of boundary slopes, each of the
# form [x,y] or [x,y,z]. It may be the case that no pairs are returned, but
# at least one triple will always be returned. Pairs precede triples, and both the set of pairs
# and the set of triples are sorted with respect to lexicograhpic order

# How to interpret the output:
# (See the above mentioned paper for more details.)
# Let X be the exterior of the link with boundary tori N1 and N2.
#
# Doubles [x,y]:
# This means there are incompressible, boundary incompressible surfaces in the
# link exterior X that have an equal number of boundary components on N1 and N2,
# and have slopes of x+ys and x-ys respectively on each boundary component of X
# where s is any rational number in the interval [-1,1]. Slopes are given with respect 
# to a preferred longitude and meridian.
#
# Triples [x,y,z]:
# This means there are incompressible, boundary incompressible survaces in X with 
# alpha boundary components on N1, beta boundary components
# on N2, and slopes of x+y beta/alpha and z+ y alpha/beta on N1 and N2 respectively, 
# where alpha>beta. Furthermore, by symmetry of the link, there are surfaces with slopes
# z+y beta/alpha and x+y alpha/beta  on N1 and N2 respectively, where alpha<beta.
# We may derive limiting cases by letting beta=0 or alpha=0, but the cases with alpha=beta
# are described by the doubles [x,y] already mentioned.


import sys, re

out = sys.stdout.write
get_input = sys.stdin.readline
##################################################
def euclid(p,q):
	""" Returns [m, n] where mp+nq=gcd(p,q)
	"""
	# performs Euclidean algorithm
	pp=p
	qq=q
	A=[1,0,0,1]
	while pp>0:
		(k,r)=divmod(qq,pp)
		qq=pp
		pp=r
		A=[A[2], A[3], A[0]-k* A[2], A[1]-k*A[3]]
	
	return [A[1], A[0]]
##################################################
def parents1(p,q):
	"""finds the parents of p/q in the graph D_1
	assumes p and q are relatively prime, and nonnegative
	"""
	[d,b]= euclid(p,q)
	d=abs(d)
	b=abs(b)
	a=p-b
	c=q-d
	return[[a,c],[b,d]]
	
##################################################
def minpaths1(p,q):
	"""recursively  finds all min paths in D_1 from 1/0 to p/q
	assumes p and q are relatively prime and nonnegative
	"""
	if [p,q]==[1,0]:
		return [[[1,0]]]
	elif q==1:
		return [[[1,0], [p,1]]]
	else:
		allpaths=[]
		for [r,s] in parents1(p,q):
			shorterpaths=minpaths1(r,s)
			# decide if each minimal path to the parent remains 
			# minimal when extended to the child
			for path in shorterpaths:
				[a,b]=path[-2]
				if abs(a*q-b*p)!=1:
					path.append ([p,q])
					allpaths.append(path)
		return allpaths
##################################################
def rot(p,q):
	""" rotates the Poincare disk 90 degrees counterclockwise, 
	taking the graph D1 to the graph D0 and vice vera
	expects p and q relatively prime, nonnegative and p<q
	unless (p,q)=(1,0)
	"""
	if (p,q)==(1,0):
		return [0,1]
	elif  q%2==0:
		return [q/2, q-p]
	else:
		return [q, 2*(q-p)]
##################################################
def Rot(p,q):
	""" opposite rotation of rot
	expect p and q relatively prime, nonnegative and p>q
	"""
	if (p,q)==(0,1):
		return [1,0]
	elif q%2==0:
		return[p-q/2, p]
	else:
		return [2*p-q, 2*p]
##################################################
def minpaths0(p,q):
	""" finds all min paths in D_0 from 1/0 to p/q
	"""
	# the idea is to take the D0 graph to D1, find minimal paths there
	# using minpath1(p,q), then rotate back. This gives minimal
	# paths in D0 that start at 1/1 and end at p/q. So we then need to see
	# how to modify the beginnings of these paths to get minimal paths
	# from 1/0 to p/q
	if q==1 or q==2:
		return[[[1,0],[p,q]]]
	
	else:
		if p>q:
			#use Rot(p,q) so that r and s will be nonnegative
			[r,s]=Rot(p,q)
			paths1=minpaths1(r,s)
			paths0=[]
			# rotate back with rot
			for path in paths1:
				paths0.append([rot(a[0], a[1]) for a in path])
		else:
			#use rot(p,q) so that r and s will be nonnegative
			[r,s]=rot(p,q)
			paths1=minpaths1(r,s)
			paths0=[]
			#rotate back with Rot
			for path in paths1:
				paths0.append([Rot(a[0], a[1]) for a in path])
		finalpaths=[]
		for path in paths0:
			#adjust the beginning of each path
			if path[1]==[1,0]:
				finalpaths.append(path[1:])
			elif path[1]==[1,2]:
				if  path[2]==[0,1]:
					new=[[1,0]]
					new.extend(path[2:])
					finalpaths.append(new)
				else:	
					new=[[1,0]]
					new.extend(path[1:])
					finalpaths.append(new)
			else:
				new=[[1,0]]
				new.extend(path)
				finalpaths.append(new)
		return finalpaths			
##################################################
def M1(path):
	""" returns value of M(path) for path a minimal path in D_1
	we know M=((x+ys)beta, (x-ys)beta) with -1<=s<=1
	and x=k-P+N, y=P+N. This function returns [x,y]
	"""
	k=P=N=0
	for i in range(len(path)-2):
		[p,q]=path[i+1]
		[r,s]=path[i+2]
		det=p*s-r*q
		if (q+s)%2==1:
			#this is an A-type edge
			k+=det
		else:
			#this is a C-type edge
			if  p%2==0:
				k+=det+1
				P+=1
			else:
				k+=det-1
				N+=1
	return[k-P+N, P+N]			
##################################################
def lk(p,q):
	"""computes linking number between K_i and the longitude lambda_i
	"""
	sum=0
	for j in range((q-2)/2):
		sum+=-pow(-1, (2*(j+1)*p)/q)
	return sum
	
##################################################
def bdyslopes(p,q):
	""" returns the set of boundary slopes for the p/q 2-bridge link
	""" 						

	linkingnumber=lk(p,q) 
	#linking number between K_1 and the longitude lambda_1
	# with respect to which all initial calculations are done.
	
	#first find the slopes for t=alpha/beta=1 that are derived from minimal paths in D_1
	#these are of the form (x+ys, x-ys) where x and y are integers of opposite parity and 
	# s is a rational parameter with -1<=s<=1. The slope x+ys is the boundary slope
	# (preferred long/meridian) on the first component and the slope x-ys is the slope
	#on the second component. The program returns the set of all such pairs [x,y] for
	# which y!=0
	
	slopes1=[]
	pathset1=minpaths1(p,q)
	for path in pathset1:
		[x,y]=M1(path)
		if y!=0:
			slopes1.append([x+linkingnumber, y]) #adjusts to preferred longitude
			

	
	# now find boundary slopes for t>1. These correpsond to minimal paths in D_t
	# Each such path specializes to a minimal path in D_0 and D_1.
	# Thus, we must examine each pairing of a minimal path in D_1 with a minimal
	#path in D_0 to see if they give rise to a minimal path in D_t. Then for that
	#path we must compute M.
	
	pathset0=minpaths0(p,q)
	
	slopest=[]
	for path1 in pathset1:
		for path0 in pathset0:
			compatible=True
			i=x=y=z=w=0
			while compatible==True and i<len(path1)-1:
				[p,q]=path1[i]
				[r,s]=path1[i+1]
				det=p*s-q*r
				if (q*s)%2==0:
					# this is an A-type edge
					if s%2==0: #trade the vertices so that q is even
						temp=[r,s]
						[r,s]=[p,q]
						[p,q]=temp
					tempdet=p*s-q*r
					if [p,q] not  in path0 or [p+tempdet*r, q+tempdet*s] in path0:
						compatible=False
					elif i>0:
						x+=det
						w+=det
													
						
				else:
					#this is a C-type edge
					if [p+r, q+s] in path0 and [r-p, s-q] in path0:
						compatible=False
					elif [p+r, q+s] not in path0 and [r-p, s-q] not in path0:
						compatible=False
					elif [p+r, q+s] in path0 :
						x+=2*det
					else:
						w+=2*det
				i+=1
		
			if compatible==True:
				# ajdust for the D-type edges of path0
				for i in range(len(path0)-1):
					[p,q]=path0[i]
					[r,s]=path0[i+1]
					det=p*s-q*r
					if abs(det)==2:
						if  [(p+r)/2, (s+q)/2] in path1:
							x+=-det/2;y+=det/2;z+=det/2;w+=-det/2
						else:
							x+=det/2;y+=-det/2;z+=-det/2;w+=det/2
						
				slopest.append([x+linkingnumber,y,w+linkingnumber])#adjust to preferred longitude
				
# choose one of following two lines for different ways to return the output				
#	return sortandtrim(slopes1)+sortandtrim(slopest)
	return [sortandtrim(slopes1), sortandtrim(slopest)]
##################################################
def sortandtrim(data):
	""" sorts a list and removes all repeated elements.
	returns the new list
	"""
	
	data.sort()
	i=0
	while i<len(data):
		item=data[i]
		j=i
		while j<len(data) and item==data[j]:
			j+=1
		del data[i+1:j]
		i+=1
	return data
##################################################
# a simple routine to compute gcd, needed to test input in UI
def gcd(a,b):
	while b:
		a,b = b, a % b
	return a
#########################################################
# ask user Question and get back a one character response which must be from Answers
# used in UI

def ask_question(Question, Answers):
	while 1:
		out(Question)
		inp = get_input()[0]    
		if re.search(inp, Answers):
	   		break
	   	out("Sorry, input not understood.\n")
	return inp
##########################################################
# This is a simple front-end that asks the user for 
# p and q and then computes the boundary slope of the
# p/q 2-bridge link. Output is reformatted to match
# Table 5 in "Computing boundary slopes of 2-bridge links"
# by Hoste and Shanahan

if __name__ == "__main__":
# main program
	while 1: # get p/q  and compute boundary slopes until user chooses to quit
		while 1: #get p/q until valid input
			while 1: # get p until valid input
				try:
					p=int(raw_input("Enter p."))
					break
				except ValueError:
					print "Not a valid integer. Try again."
			while 1: # get q until valid input
				try:
					q=int(raw_input("Enter q."))
					break
				except ValueError:
					print "Not a valid integer. Try again."
			#see if p and q are valid input
			if 0<p and p<q and q%2==0 and gcd(p,q)==1:
				break
			else:
				print "You must enter 0<p<q, q even, and gcd(p,q)=1"
		# compute boundary slopes using these values of p and q
		slopes=bdyslopes(p,q)
		print "\n"
		print "Boundary slopes of surfaces with equal number of sheets on each boundary component."
		print "Here s is any rational number in [-1,1].\n" 

		for [x,y] in slopes[0]:
#			print [x,y]
			if x==0:
				print  "(%ds, %ds)" % (y,-y)
			else:
				if y<-1:
					print "(%d%ds,%d+%ds)" %(x,y,x,-y)
				elif y==-1:
					print "(%d-s,%d+s)" %(x,x)
				elif y==0:
					print "(%d,%d)" %(x,x)
				elif y==1:
					print "(%d+s,%d-s)" %(x,x)
				else:
					print "(%d+%ds,%d-%ds)" %(x,y,x,y)
		print "\n"
		print "Boundary slopes of surfaces with alpha sheets on component 1 and beta sheets"
		print "on component 2 where t=alpha/beta>1.\n"
		for [x,y,z] in slopes[1]:
#			print [x,y,z]
			if x==0:
				if y==0:
					s1="0"
				else:
					s1="%d/t" % y
			else:
				if y<0:
					s1="%d%d/t" % (x,y)
				elif y==0:
					s1="%d" % x
				else:
					s1="%d+%d/t" % (x,y)
	
			if z==0:
				if y<-1:
					s2="%dt" % y
				elif y==-1:
					s2="-t"
				elif y==0:
					s2="0"
				elif y==1:
					s2="t"
				else:
					s2="%dt" % y
			else:
				if y<-1:
					s2="%d%dt" % (z,y)
				elif y==-1:
					s2="%d-t" % z
				elif y==0:
					s2="%d" % z
				elif y==1:
					s2="%d+t" % z
				else:
					s2="%d+%dt" % (z,y)
	
	
			print "".join(["(", s1, ",", s2, ")"])

		if ask_question("Another link? (y/n): ", "yn") == "n":
			break