import math
from sympy.ntheory import factorint
from rings import Zphi

class HZphi:
	"""
	Z[phi].
	Functions implemented are:
		 Arithmetic functions: +,*,//,%,**(exponentiation)
		 a.gcd(b) - Compute the greatest common divisor of a and b.
		 a.xgcd(b) - Extended gcd - return gcd(a,b) and x,y such that gcd(a,b)=xa+yb.
		 n.isprime() - Is n prime (pseudoprime tests)
		 n.factor() - Return a factor of n.
		 n.factors() - Return list of the factors of n.
	Can be created by:
		 n = HZphi(5,7)  # Create (5 + 7i)
		 n = HZphi(13)  # Create (5 + 0i)
	"""
	# Adapted from Robert Campell and Hubert Holin, https://github.com/Robert-Campbell-256/Number-Theory-Python/blob/master/gaussint.py
	
	def __init__(self, a = 0, b = 0, c = 0, d = 0):
		self.one, self.i, self.j, self.k = Zphi(a), Zphi(b), Zphi(c), Zphi(d)
	
	def __str__(self): # Overload string conversion used by print
		return "(" + str(self.one) + ") + (" + str(self.i) + ")i + (" + str(self.j) + ")j + (" + str(self.k) + ")k"
		
	def __format__(self, spec): # Overload string conversion used by format
		return "(" + str(self.one) + ") + (" + str(self.i) + ")i + (" + str(self.j) + ")j + (" + str(self.k) + ")k"
	
	def __repr__(self): # Overload conversion used for output # return "HZphi(" + str(self.one) + ", " + str(self.phi) + ")"
		return "(" + str(self.one) + ") + (" + str(self.i) + ")i + (" + str(self.j) + ")j + (" + str(self.k) + ")k"
	
	def tuple(self):
		return (self.one,self.i,self.j,self.k)

	def __eq__(self,other):  # Overload the "==" test operator
		if (type(other) is not HZphi):
			return False
		else:
			return (self.one == other.one) and (self.i == other.i) and (self.j == other.j) and (self.k == other.k)
	
	def __ne__(self,other):  # Overload the "!=" test operator
		return not (self == other)
	
	"""
	def neutral_element_for_multiplication():
		return __class__(1)
	"""
	
	def conjugate(self):
		return HZphi(self.one, -self.i, -self.j, -self.k)
	
	def norm(self):
		return self.one*self.one + self.i*self.i + self.j*self.j + self.k*self.k
	
	def __pos__(self):	  # Overload the "+" unary operator
		return self
	
	def add(self,summand):
		sum_one = self.one + summand.one
		sum_i = self.i + summand.i
		sum_j = self.j + summand.j
		sum_k = self.k + summand.k
		return HZphi(sum_one, sum_i, sum_j, sum_k)
	
	def __add__(self,summand):  # Overload the "+" binary operator
		return self.add(summand)
	
	def __radd__(self,summand):  # Overload the "+" binary operator
		return self.add(summand)
	
	def __iadd__(self,summand): # Overload the "+=" operator
		self = self + summand
		return self
	
	def __neg__(self):  # Overload the "-" unary operator
		return HZphi(-self.one,-self.i,-self.j,-self.k)
	
	def __sub__(self,summand):  # Overload the "-" binary operator
		return self.__add__(-summand)
	
	def __rsub__(self,summand):  # Overload the "-" binary operator
		return summand-self
	
	def __isub__(self,summand): # Overload the "-=" operator
		self = self - summand
		return self
	
	def mult(self,multip):
		m = multip
		if type(multip) is int or type(multip) is Zphi:
			m = HZphi(multip)
		prod_one = self.one*m.one - self.i*m.i - self.j*m.j - self.k*m.k
		prod_i = self.one*m.i + self.i*m.one + self.j*m.k - self.k*m.j
		prod_j = self.one*m.j + self.j*m.one - self.i*m.k + self.k*m.i
		prod_k = self.one*m.k + self.k*m.one - self.j*m.i + self.i*m.j
		return HZphi(prod_one, prod_i, prod_j, prod_k)
	
	### double check that these do the right thing (order)
	def __mul__(self,multip):  # Overload the "*" operator
		return self.mult(multip)
	
	def __rmul__(self,multip):  # Overload the "*" operator
		return self.mult(multip)
	
	def projectivelyEqual(self,other): # do self and other differ by a Qphi multiple?
		if self == 0:
			return other == 0
		if self.one != 0:
			return self.one*other.i == self.i*other.one and self.one*other.j == self.j*other.one and self.one*other.k == self.k*other.one
		if self.i != 0:
			return self.i*other.one == self.one*other.i and self.i*other.j == self.j*other.i and self.i*other.k == self.k*other.i
		if self.j != 0:
			return self.j*other.one == self.one*other.j and self.j*other.i == self.i*other.j and self.j*other.k == self.k*other.j
		return self.k*other.one == self.one*other.k and self.k*other.i == self.i*other.k and self.k*other.j == self.j*other.k
	
	###is this needed?
	#"""
	def floordiv(self,divisor):
		if type(divisor) is int:
			numerator = (-self if (divisor < 0) else self)
			denominator = (-divisor if (divisor < 0) else divisor)
			if denominator == 0:
				raise ZeroDivisionError("{0:s} is null!".format(divisor))
		else:
			numerator = self*divisor.conjugate()
			denominator = divisor.norm()	# Recall that denominator >= 0
			if denominator == 0:
				raise ZeroDivisionError("{0:s} is null!".format(divisor))
		candidate_r = numerator.one//denominator
		candidate_i = numerator.phi//denominator

		# i.e. (candidate_r+1)*denominator-numerator.one < numerator.one-candidate_r*denominator
		if (2*candidate_r+1)*denominator < 2*numerator.one:
			candidate_r += 1

		# i.e. (candidate_i+1)*denominator-numerator.phi < numerator.phi-candidate_i*denominator
		if (2*candidate_i+1)*denominator < 2*numerator.phi:
			candidate_i += 1
		return HZphi(candidate_r,candidate_i)
	
	def __floordiv__(self,divisor):  # Overload the "//" operator
		# return self.floordiv(divisor)
		a = self.one//divisor
		b = self.i//divisor
		c = self.j//divisor
		d = self.k//divisor
		return HZphi(a,b,c,d)
	
	def __ifloordiv__(self,divisor): # Overload the "//=" operator
		self = self//divisor
		return self
	
	"""
	def mod(self,divisor):
		return self - divisor * (self//divisor)
	
	def __mod__(self,divisor):  # Overload the "%" operator
		return self.mod(divisor)
	
	def __imod__(self,divisor): # Overload the "%=" operator
		self = self % divisor
		return self
	"""
	
	### is this needed?
	"""
	def divmod(self,divisor):
		q = self//divisor
		return q, self - divisor * q
	"""
	
	"""
	def xgcd(self,other):
		quot = HZphi()
		a1 = HZphi(1,0)
		b1 = HZphi(0,0)
		a2 = HZphi(0,0)
		b2 = HZphi(1,0)
		a = self
		b = other
		if(b.norm() > a.norm()):	# Need to start with a>b
			a,b = b,a					# Swap a and b
			a1,b1,a2,b2 = a2,b2,a1,b1	# Swap (a1,b1) with (a2,b2)
		while (True):
			quot = a // b
			a %= b
			a1 -= quot*a2
			b1 -= quot*b2
			if (a == HZphi(0,0)):
				return b, a2, b2
			quot = b // a
			b %= a
			a2 -= quot*a1
			b2 -= quot*b1
			if (b == HZphi()):
				return a, a1, b1
	"""
	"""
	def Bezout(self, other):
		a = self
		b = other
		if a.norm() < b.norm():
			(u, v, pgcd) = b.Bezout(a)
			return (v, u, pgcd)
		if b == 0:
			return (1, 0, a)
		u_n, u_n_moins_1, v_n, v_n_moins_1 = 0, 1, 1, 0
		while b.norm() > 0:
			q,r = a.divmod(b)
			u_n_plus_1 = u_n_moins_1 - q*u_n
			v_n_plus_1 = v_n_moins_1 - q*v_n
			a, b = b, r
			u_n_moins_1, u_n, v_n_moins_1, v_n = u_n, u_n_plus_1, v_n, v_n_plus_1
		return (u_n_moins_1, v_n_moins_1, a)
	"""
	
	### is this needed?
	"""
	def gcd(self,other):
		a = self
		b = other
		if abs(a.norm()) < abs(b.norm()):
			return b.gcd(a)
		while abs(b.norm()) > 0:
			q,r = a.divmod(b)
			a,b = b,r
		return a
  """
	
	### is this needed?
	"""
	def maxPower(self,div):
		x = self*div.conjugate()
		n = div.norm()
  		if self == 0 or div == 0 or abs(n) == 1:
			return 0
		elif x.one%n == 0 and x.phi%n == 0:
			return 1 + HZphi(x.one/n,x.phi/n).maxPower(div)
		else:
			return 0
   """
	
	### is this needed?
	"""
	def factor(self):
		n = self.norm()
		fact = factorint(n)
		d = {}
		for p in fact.keys():
			if p == 5:
				d[(-1,2)] = self.maxPower(HZphi(-1,2))
			elif p%5 == 2 or p%3 == 3:
				d[(p,0)] = self.maxPower(HZphi(p))
			else:
				y = tonelli(5,p)
				x = (y+1)*(p+1)/2
				a = x - HZphi(0,1)
				q = a.gcd(HZphi(p))
				d[q.pair()] = self.maxPower(q)
				r = q.conjugate()
				d[r.pair()] = d[q.pair()]
		return d
	"""
	
	"""
	def powmod(self, a_power, a_modulus):
	# We adapt the Binary Exponentiation algorithm with modulo
		result = HZphi(1)
		auxilliary = HZphi(self.one, self.phi)
		while a_power:
			if a_power % 2:	# If power is odd
				result = (result * auxilliary) % a_modulus
			# Divide the power by 2
			a_power >>= 1
			# Multiply base to itself
			auxilliary = (auxilliary * auxilliary) % a_modulus
		return result
	"""
	
	""" stopping here to see if the stuff above actually works
	def __pow__(self, a_power):  # Overload the "**" operator
	# We adapt the Binary Exponentiation algorithm (without modulo!)
		
		result = HZphi(1)
		
		auxilliary = HZphi(self.one, self.phi)
		
		while a_power:
			
			if a_power % 2:	# If power is odd
				
				result = result * auxilliary
			
			# Divide the power by 2
			a_power >>= 1
			
			# Multiply base to itself
			auxilliary = auxilliary * auxilliary
		
		return result
	
	
	def isprime(self):
		# n.isprime() - Test whether the HZphi n is prime using a variety of pseudoprime tests.
		# Multiply by (1,i,-1,-i) to rotate to first quadrant (similar to abs)
		if (self.one < 0): self *= (-1)
		if (self.phi < 0): self *= HZphi(0,1)
		# Check some small non-primes
		if (self in [HZphi(0,0), HZphi(1,0), HZphi(0,1)]): return False
		# Check some small primes
		if (self in [HZphi(1,1), HZphi(2,1), HZphi(1,2), HZphi(3,0), HZphi(0,3), HZphi(3,2), HZphi(2,3)]):
			return True
		return self.isprimeF(2) and self.isprimeF(3) and self.isprimeF(5)
	
	
	def isprimeF(self,base):
		# n.isprimeF(base) - Test whether the HZphi n is prime using the Gaussian Integer analogue of the Fermat pseudoprime test.
		if type(base) is not HZphi:
			 base = HZphi(base)  # Coerce if base not HZphi (works for int or complex)
		return base.powmod(self.norm()-1,self) == HZphi(1,0)
	# Note: Possibly more effective would be to use the characterization of primes
	# in the Gaussian Integers based on the primality of their norm and reducing mod 4.
	# This depends on the characterization of the ideal class group, and only works for
	# simple rings of algebraic integers.
	
	
	def factor(self):
		# n.factor() - Find a prime factor of Gaussian Integer n using a variety of methods.
		if (self.isprime()): return n
		for fact in [HZphi(1,1), HZphi(2,1), HZphi(1,2),
				 HZphi(3,0), HZphi(3,2), HZphi(2,3)]:
			if self%fact == 0: return fact
		return self.factorPR()  # Needs work - no guarantee that a prime factor will be returned
	
	
	def factors(self):
		# n.factors() - Return a sorted list of the prime factors of Gaussian Integer n.
		if (self.isprime()):
			  return [self]
		fact = self.factor()
		if (fact == 1): return "Unable to factor "+str(n)
		facts = (self/fact).factors() + fact.factors()
		return facts
	
	
	def factorPR(self):	# TODO: learn and test
		# n.factorPR() - Find a factor of Gaussian Integer n using the analogue of the Pollard Rho method.
		# Note: This method will occasionally fail.
		for slow in [2,3,4,6]:
			numsteps=2*math.floor(math.sqrt(math.sqrt(self.norm()))); fast=slow; i=1
			while i<numsteps:
				slow = (slow*slow + 1) % self
				i = i + 1
				fast = (fast*fast + 1) % self
				fast = (fast*fast + 1) % self
				g = gcd(fast-slow,self)
				if (g != 1):
					if (g == self):
						break
					else:
						return g
			return 1
	# Note: Possibly more effective would be to factor the norm and then use the known
	# splitting properties of primes over the Gaussian Integers. This depends on the
	# characterization of the ideal class group, and only works for simple rings of algebraic integers.
	"""

"""
x = HZphi(3,4)
y = HZphi(5)
print("x = " + str(x))
print("y = " + str(y))
print("x-dot = " + str(x.conjugate()))
print("y-dot = " + str(y.conjugate()))
print("-x = " + str(-x))
print("-y = " + str(-y))
print("x+y = " + str(x+y))
print("x-y = " + str(x-y))
print("x*y = " + str(x*y))
print("N(x) = " + str(x.norm()))
print("N(y) = " + str(y.norm()))
print("x =?= y is " + str((x == y)))
print("x =/?= y is " + str((x != y)))
print("gcd(x,y) = " + str(x.gcd(y)))
print("x * x-dot = " + str(x*x.conjugate()))
ele = HZphi(11)
nin = HZphi(19)
twe = HZphi(29)
thi = HZphi(31)
print ele.factor(), nin.factor(), twe.factor(), thi.factor()
"""

"""
x = HZphi(1,2,0,0)
y = HZphi(2,4,0,0)
t = HZphi(0,Zphi(2,1),1,1)
print(t.norm())
r = HZphi(0,1,Zphi(-1,1),Zphi(0,1))
print r*r
print(r.norm())
s = HZphi(1,1,1,1)
print s*s
print(s.norm())

print(t*t)
print(r*s)
print(s*r)

rd = HZphi(2,2,2,2)

print r.projectivelyEqual(r), s.projectivelyEqual(t), x.projectivelyEqual(y), y.projectivelyEqual(x), r.projectivelyEqual(rd)

p = Zphi(7,5)
t = HZphi(0,Zphi(2,1),1,1)

Corig = [HZphi(1,0,0,0), r, s, r*s, s*r, s*s, r*s*r, r*s*s, s*r*s, s*s*r, r*s*r*s, r*s*s*r, s*r*s*r, s*r*s*s, s*s*r*s, r*s*r*s*r, r*s*r*s*s, r*s*s*r*s, s*r*s*r*s, s*r*s*s*r, s*s*r*s*r, s*s*r*s*s, r*s*r*s*r*s, r*s*r*s*s*r, r*s*s*r*s*r, r*s*s*r*s*s, s*r*s*r*s*s, s*r*s*s*r*s, s*s*r*s*r*s, r*s*r*s*r*s*s, r*s*r*s*s*r*s, r*s*s*r*s*r*s, s*r*s*r*s*s*r, s*r*s*s*r*s*r, s*r*s*s*r*s*s, s*s*r*s*r*s*s, r*s*r*s*r*s*s*r, r*s*r*s*s*r*s*r, r*s*r*s*s*r*s*s, r*s*s*r*s*r*s*s, s*r*s*r*s*s*r*s, s*r*s*s*r*s*r*s, s*s*r*s*r*s*s*r, r*s*r*s*r*s*s*r*s, r*s*r*s*s*r*s*r*s, r*s*s*r*s*r*s*s*r, s*r*s*r*s*s*r*s*r, s*r*s*r*s*s*r*s*s, s*r*s*s*r*s*r*s*s, s*s*r*s*r*s*s*r*s, r*s*r*s*r*s*s*r*s*r, r*s*r*s*r*s*s*r*s*s, r*s*r*s*s*r*s*r*s*s, r*s*s*r*s*r*s*s*r*s, s*r*s*s*r*s*r*s*s*r, s*s*r*s*r*s*s*r*s*r, r*s*r*s*s*r*s*r*s*s*r, r*s*s*r*s*r*s*s*r*s*r, s*r*s*s*r*s*r*s*s*r*s, r*s*r*s*s*r*s*r*s*s*r*s]
# looks good!
"""

p = Zphi(7,5)
t = HZphi(0,Zphi(2,1),1,1)
r = HZphi(0,1,Zphi(-1,1),Zphi(0,1))
s = HZphi(1,1,1,1)

C = [HZphi(1,0,0,0), r, s, r*s//2, s*r//2, s*s//2, r*s*r//4, r*s*s//4, s*r*s//4, s*s*r//4, r*s*r*s//8, r*s*s*r//8, s*r*s*r//8, s*r*s*s//8, s*s*r*s//8, r*s*r*s*r//16, r*s*r*s*s//16, r*s*s*r*s//16, s*r*s*r*s//16, s*r*s*s*r//16, s*s*r*s*r//16, s*s*r*s*s//16, r*s*r*s*r*s//32, r*s*r*s*s*r//64, r*s*s*r*s*r//32, r*s*s*r*s*s//32, s*r*s*r*s*s//32, s*r*s*s*r*s//32, s*s*r*s*r*s//32, r*s*r*s*r*s*s//64, r*s*r*s*s*r*s//64, r*s*s*r*s*r*s//64, s*r*s*r*s*s*r//64, s*r*s*s*r*s*r//64, s*r*s*s*r*s*s//64, s*s*r*s*r*s*s//64, r*s*r*s*r*s*s*r//128, r*s*r*s*s*r*s*r//128, r*s*r*s*s*r*s*s//128, r*s*s*r*s*r*s*s//128, s*r*s*r*s*s*r*s//128, s*r*s*s*r*s*r*s//128, s*s*r*s*r*s*s*r//128, r*s*r*s*r*s*s*r*s//256, r*s*r*s*s*r*s*r*s//256, r*s*s*r*s*r*s*s*r//256, s*r*s*r*s*s*r*s*r//256, s*r*s*r*s*s*r*s*s//512, s*r*s*s*r*s*r*s*s//256, s*s*r*s*r*s*s*r*s//512, r*s*r*s*r*s*s*r*s*r//512, r*s*r*s*r*s*s*r*s*s//512, r*s*r*s*s*r*s*r*s*s//512, r*s*s*r*s*r*s*s*r*s//512, s*r*s*s*r*s*r*s*s*r//512, s*s*r*s*r*s*s*r*s*r//512, r*s*r*s*s*r*s*r*s*s*r//1024, r*s*s*r*s*r*s*s*r*s*r//1024, s*r*s*s*r*s*r*s*s*r*s//1024, r*s*r*s*s*r*s*r*s*s*r*s//2048]
# print C

def isReduced(p,quat):
	a, b, c, d = quat.one, quat.i, quat.j, quat.k
	return a.maxPower(p) <= 0 or b.maxPower(p) <= 0 or c.maxPower(p) <= 0 or d.maxPower(p) <= 0

# print isReduced(p,HZphi(p,0,0,0)), isReduced(p,t), isReduced(p,t*t)

def inC60(quat): # returns the position of quat (projectively) in the list of C's elements
	for c in range(60):
		if quat.projectivelyEqual(C[c]):
			return c
	return 60

# print inC60(s), inC60(HZphi(2,2,2,2)), inC60(t)

# print "======== TESTING factorGamma ========"
def factorGamma(g):
	i = inC60(g)
	# print g, i
	if i < 60:
		return [i]
	for i in range(60):
		h = g*C[i]*t
		if not isReduced(p,h):
			# print "unreduced h is", h
			a, b, c, d = h.one, h.i, h.j, h.k
			a2 = (a*p.conjugate())//59
			b2 = (b*p.conjugate())//59
			c2 = (c*p.conjugate())//59
			d2 = (d*p.conjugate())//59
			h2 = HZphi(a2,b2,c2,d2)
			# print "reduced h is", h2
			cn = C[i].conjugate()
			j = inC60(cn)
			# print i, j
			L = factorGamma(h2)
			# print L
			return L + [j]

# print factorGamma(s), factorGamma(t), factorGamma(C[59]*t*C[58]*t*C[40]), factorGamma(s*t), factorGamma(t*r)

TeX = { 0 : " ", 1 : "\\rho ", 2 : "\\sigma ", 3 : "\\rho\\sigma ", 4 : "\\sigma\\rho ", 5 : "\\sigma\\sigma ", 6 : "\\rho\\sigma\\rho ", 7 : "\\rho\\sigma\\sigma ", 8 : "\\sigma\\rho\\sigma ", 9 : "\\sigma\\sigma\\rho ", 10 : "\\rho\\sigma\\rho\\sigma ", 11 : "\\rho\\sigma\\sigma\\rho ", 12 : "\\sigma\\rho\\sigma\\rho ", 13 : "\\sigma\\rho\\sigma\\sigma ", 14 : "\\sigma\\sigma\\rho\\sigma ", 15 : "\\rho\\sigma\\rho\\sigma\\rho ", 16 : "\\rho\\sigma\\rho\\sigma\\sigma ", 17 : "\\rho\\sigma\\sigma\\rho\\sigma ", 18 : "\\sigma\\rho\\sigma\\rho\\sigma ", 19 : "\\sigma\\rho\\sigma\\sigma\\rho ", 20 : "\\sigma\\sigma\\rho\\sigma\\rho ", 21 : "\\sigma\\sigma\\rho\\sigma\\sigma ", 22 : "\\rho\\sigma\\rho\\sigma\\rho\\sigma ", 23 : "\\rho\\sigma\\rho\\sigma\\sigma\\rho ", 24 : "\\rho\\sigma\\sigma\\rho\\sigma\\rho ", 25 : "\\rho\\sigma\\sigma\\rho\\sigma\\sigma ", 26 : "\\sigma\\rho\\sigma\\rho\\sigma\\sigma ", 27 : "\\sigma\\rho\\sigma\\sigma\\rho\\sigma ", 28 : "\\sigma\\sigma\\rho\\sigma\\rho\\sigma ", 29 : "\\rho\\sigma\\rho\\sigma\\rho\\sigma\\sigma ", 30 : "\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma ", 31 : "\\rho\\sigma\\sigma\\rho\\sigma\\rho\\sigma ", 32 : "\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho ", 33 : "\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\rho ", 34 : "\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\sigma ", 35 : "\\sigma\\sigma\\rho\\sigma\\rho\\sigma\\sigma ", 36 : "\\rho\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho ", 37 : "\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\rho ", 38 : "\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\sigma ", 39 : "\\rho\\sigma\\sigma\\rho\\sigma\\rho\\sigma\\sigma ", 40 : "\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma ", 41 : "\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\rho\\sigma ", 42 : "\\sigma\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho ", 43 : "\\rho\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma ", 44 : "\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\rho\\sigma ", 45 : "\\rho\\sigma\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho ", 46 : "\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\rho ", 47 : "\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\sigma ", 48 : "\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\rho\\sigma\\sigma ", 49 : "\\sigma\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma ", 50 : "\\rho\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\rho ", 51 : "\\rho\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\sigma ", 52 : "\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\rho\\sigma\\sigma ", 53 : "\\rho\\sigma\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma ", 54 : "\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho ", 55 : "\\sigma\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\rho ", 56 : "\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho ", 57 : "\\rho\\sigma\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\rho ", 58 : "\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma ", 59 : "\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma\\rho\\sigma\\sigma\\rho\\sigma " }

def TeXGamma(L):
	s = ""
	for i in L[:-1]:
		s += TeX[i] + "\\tau "
	return s + TeX[L[-1]]

"""
print TeXGamma(factorGamma(s*t))
print TeXGamma(factorGamma(t*r))
print TeXGamma(factorGamma(C[59]*t*C[58]*t*C[40]))
# looks good!
"""