from math import sqrt, log
from sympy.ntheory import factorint
from numpy import inf
from Ziphi_dict import D
import mpmath as mp
from mpmath import sqrt

phi = (1+sqrt(5))/2

def legendre(a, p):
	return pow(a, (p - 1) // 2, p)

def tonelli(n, p): # from https://rosettacode.org/wiki/Tonelli-Shanks_algorithm#Python
	assert legendre(n, p) == 1, "not a square mod p"
	q = p - 1
	s = 0
	while q % 2 == 0:
		q //= 2
		s += 1
	if s == 1:
		return pow(n, (p + 1) // 4, p)
	for z in range(2, p):
		if p - 1 == legendre(z, p):
			break
	c = pow(z, q, p)
	r = pow(n, (q + 1) // 2, p)
	t = pow(n, q, p)
	m = s
	t2 = 0
	while (t - 1) % p != 0:
		t2 = (t * t) % p
		for i in range(1, m):
			if (t2 - 1) % p == 0:
				break
			t2 = (t2 * t2) % p
		b = pow(c, 1 << (m - i - 1), p)
		r = (r * b) % p
		c = (b * b) % p
		t = (t * c) % p
		m = i
	return r

class Zphi:
	"""
	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 = Zphi(5,7)  # Create (5 + 7i)
		 n = Zphi(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):
		if isinstance(a,Zphi):
			self.one, self.phi = a.one, a.phi
			self.isNone = a.isNone
		else:
			self.one, self.phi = a, b
			if isinstance(a, int):
				self.isNone = False
			else: # a had better be None
				self.isNone = True
	
	def __str__(self): # Overload string conversion used by print
		return str(self.one) + ((" + "+str(self.phi)) if (self.phi >= 0) else (" - "+str(-self.phi))) + " phi"
		
	def __format__(self, spec): # Overload string conversion used by format
		return str(self.one) + ((" + "+str(self.phi)) if (self.phi >= 0) else (" - "+str(-self.phi))) + " phi"
	
	def __repr__(self): # Overload conversion used for output # return "Zphi(" + str(self.one) + ", " + str(self.phi) + ")"
		if self.isNone:
			return "None"
		elif (self.phi >= 0):
			return str(self.one) + " + " + str(self.phi) + " phi"
		else:
			return str(self.one) + " - " + str(-self.phi) + " phi" # return str(self.one) + ((" + "+str(self.phi)) if (self.phi >= 0) else (" - "+str(-self.phi))) + " phi"
	
	def pair(self):
		return (self.one,self.phi)
	
	def float(self):
		if self.isNone:
			return 0
		return self.one+self.phi*phi

	def __eq__(self,other):  # Overload the "==" test operator
		if isinstance(other, int):
			return self.one == other and self.phi == 0
		else:
			return self.one == other.one and self.phi == other.phi
	
	def __ne__(self,other):  # Overload the "!=" test operator
		return not (self == other)
	
	"""
	def neutral_element_for_multiplication():
		return __class__(1)
	"""
	
	def conjugate(self):
		return Zphi(self.one + self.phi, -self.phi)
	
	def norm(self):
		return self.one*self.one + self.one*self.phi - self.phi*self.phi
	
	def __pos__(self):	  # Overload the "+" unary operator
		return self
	
	def add(self,summand):
		sum_one = self.one + summand.one
		sum_phi = self.phi + summand.phi
		return Zphi(sum_one, sum_phi)
	
	def __add__(self,summand):  # Overload the "+" binary operator
		if isinstance(summand, int):
			return Zphi(self.one+summand, self.phi)
		else:
			return self.add(summand)
	
	def __radd__(self,summand):  # Overload the "+" binary operator
		if isinstance(summand, int):
			return Zphi(self.one+summand, self.phi)
		else:
			return self.add(summand)
	
	def __iadd__(self,summand): # Overload the "+=" operator
		self = self + summand
		return self
	
	def __neg__(self):  # Overload the "-" unary operator
		return Zphi(-self.one,-self.phi)
	
	def __sub__(self,summand):  # Overload the "-" binary operator
		return self.__add__(-summand)
	
	def __rsub__(self,summand):  # Overload the "-" binary operator
		if isinstance(summand, int):
			return Zphi(summand-self.one, -self.phi)
		elif isinstance(summand, Zphi): # Handle subtraction with another Zphi object
			return Zphi(summand.one - self.one, summand.phi - self.phi)
		else: # For any other type, return NotImplemented
			return NotImplemented
	
	def __isub__(self,summand): # Overload the "-=" operator
		self = self - summand
		return self
	
	def mult(self,multip):
		prod_r = (self.one * multip.one) + (self.phi * multip.phi)
		prod_i = (self.phi * multip.one) + (self.one * multip.phi) + (self.phi * multip.phi)
		return Zphi(prod_r, prod_i)
	
	def __mul__(self,multip):  # Overload the "*" operator
		if isinstance(multip, int):
			return Zphi(self.one*multip, self.phi*multip)
		else:
			return self.mult(multip)
	
	def __rmul__(self,multip):  # Overload the "*" operator
		if isinstance(multip, int):
			return Zphi(self.one*multip, self.phi*multip)
		else:
			return self.mult(multip)
	
	def __imul__(self,multip): # Overload the "*=" operator
		self = self * multip
		return self
	
	#"""
	def floordiv(self,divisor):
		if isinstance(divisor, 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 Zphi(candidate_r,candidate_i)
	
	def __floordiv__(self,divisor):  # Overload the "//" operator
		return self.floordiv(divisor)
	
	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
	"""
	
	def divmod(self,divisor):
		q = self//divisor
		return q, self - divisor * q
	#"""
	
	"""
	def xgcd(self,other):
		quot = Zphi()
		a1 = Zphi(1,0)
		b1 = Zphi(0,0)
		a2 = Zphi(0,0)
		b2 = Zphi(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 == Zphi(0,0)):
				return b, a2, b2
			quot = b // a
			b %= a
			a2 -= quot*a1
			b2 -= quot*b1
			if (b == Zphi()):
				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)
	"""
	
	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
	
	def maxPower(self,div):
		x = self*div.conjugate()
		n = div.norm()
		if self == 0 or div == 0 or abs(n) == 1:
			return inf
		elif x.one%n == 0 and x.phi%n == 0:
			return 1 + Zphi(x.one/n,x.phi/n).maxPower(div)
		else:
			return 0
	
	def factor(self):
		n = abs(self.norm())
		fact = factorint(n)
		d = {}
		for p in list(fact.keys()):
			if p == 5:
				d[(-1,2)] = self.maxPower(Zphi(-1,2))
			elif p%5 == 2 or p%5 == 3:
				d[(p,0)] = self.maxPower(Zphi(p))
			else:
				y = tonelli(5,p)
				x = (y+1)*(p+1)//2
				a = x - Zphi(0,1)
				q = a.gcd(Zphi(p))
				r = q.conjugate()
				pow_1, pow_2 = self.maxPower(q), self.maxPower(r)
				if pow_1 > 0:
					d[q.pair()] = pow_1
				if pow_2 > 0:
					d[r.pair()] = pow_2
		return d

	"""
	def powmod(self, a_power, a_modulus):
	# We adapt the Binary Exponentiation algorithm with modulo
		result = Zphi(1)
		auxilliary = Zphi(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
	"""
	
	def __pow__(self, a_power):  # Overload the "**" operator
	# We adapt the Binary Exponentiation algorithm (without modulo!)
		
		result = Zphi(1)
		auxilliary = Zphi(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
	
	""" stopping here to see if the stuff above actually works
	def isprime(self):
		# n.isprime() - Test whether the Zphi 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 *= Zphi(0,1)
		# Check some small non-primes
		if (self in [Zphi(0,0), Zphi(1,0), Zphi(0,1)]): return False
		# Check some small primes
		if (self in [Zphi(1,1), Zphi(2,1), Zphi(1,2), Zphi(3,0), Zphi(0,3), Zphi(3,2), Zphi(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 Zphi n is prime using the Gaussian Integer analogue of the Fermat pseudoprime test.
		if isinstance(base, Zphi):
			 base = Zphi(base)  # Coerce if base not Zphi (works for int or complex)
		return base.powmod(self.norm()-1,self) == Zphi(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 [Zphi(1,1), Zphi(2,1), Zphi(1,2),
				 Zphi(3,0), Zphi(3,2), Zphi(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 = Zphi(3,4)
y = Zphi(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 = Zphi(11)
nin = Zphi(19)
twe = Zphi(29)
thi = Zphi(31)
print ele.factor(), nin.factor(), twe.factor(), thi.factor()
"""

"""
for n in range(3,10):
	x = Zphi(n)
	print x.factor(),
	check = Zphi(1)
	for y in x.factor().keys():
		z = Zphi(y[0],y[1])
		check *= z**x.factor()[y]
	print check
"""

"""
zero = Zphi(0,0)
print zero == 0
negone = Zphi(-1,0)
print abs(negone.norm())
x = Zphi(3,4)
print x.maxPower(zero)
print x.maxPower(negone)
print zero.maxPower(x)
"""
def unitDegree(x): # if x is a unit in Z[phi] of the form phi**n then unitDegree finds n. note that this only works for "positive" units
	if x == Zphi(1,0):
		return 0
	if x == Zphi(0,1):
		return 1
	if x == Zphi(-1,1):
		return -1
	if abs(x.norm()) != 1:
		return None
	(a,b) = x.pair()
	phi = (1+sqrt(5))/2
	m = int(round(log(sqrt(5)*abs(b))/log(phi)))
	if a < 0 and b < 0:
		return None
	if a > 0 and b > 0:
		return m
	else:
		if b == -(-1)**m * int(round(phi**m/sqrt(5))):
			return -m
		else:
			return None

"""
for n in range(6):
	print n, Zphi(0,1)**n, unitDegree(Zphi(0,1)**n)
for n in range(1,6):
	print -n, Zphi(-1,1)**n, unitDegree(Zphi(-1,1)**n)
"""

class Ziphi:
	"""
	Z[i,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.
	"""
	# 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):
		if isinstance(a,Zphi):
			if isinstance(b,Zphi):
				self.one, self.phi, self.i, self.iphi = a.one, a.phi, b.one, b.phi
			else:
				self.one, self.phi, self.i, self.iphi = a.one, a.phi, 0, 0
		elif isinstance(a,Zphi):
			a.one, a.phi, a.i, a.iphi
		else:
			self.one, self.phi, self.i, self.iphi = int(a), int(b), int(c), int(d)
	
	def __str__(self): # Overload string conversion used by print
		return str(self.one) + ((" + "+str(self.phi)) if (self.phi >= 0) else (" - "+str(-self.phi))) + " phi" + ((" + "+str(self.i)) if (self.i >= 0) else (" - "+str(-self.i))) + " i" + ((" + "+str(self.iphi)) if (self.iphi >= 0) else (" - "+str(-self.iphi))) + " i phi"
		
	def __format__(self, spec): # Overload string conversion used by format
		return str(self.one) + ((" + "+str(self.phi)) if (self.phi >= 0) else (" - "+str(-self.phi))) + " phi" + ((" + "+str(self.i)) if (self.i >= 0) else (" - "+str(-self.i))) + " i" + ((" + "+str(self.iphi)) if (self.iphi >= 0) else (" - "+str(-self.iphi))) + " i phi"
	
	def __repr__(self): # Overload conversion used for output # return "Zphi(" + str(self.one) + ", " + str(self.phi) + ")"
		return str(self.one) + ((" + "+str(self.phi)) if (self.phi >= 0) else (" - "+str(-self.phi))) + " phi" + ((" + "+str(self.i)) if (self.i >= 0) else (" - "+str(-self.i))) + " i" + ((" + "+str(self.iphi)) if (self.iphi >= 0) else (" - "+str(-self.iphi))) + " i phi"
	
	def tuple(self):
		return (self.one,self.phi,self.i,self.iphi)
	
	def pair(self):
		return (Zphi(self.one,self.phi),Zphi(self.i,self.iphi))

	def __eq__(self,other):  # Overload the "==" test operator
		if isinstance(other, int):
			return self.one == other and self.phi == 0 and self.i == 0 and self.iphi == 0
		else:
			return self.one == other.one and self.phi == other.phi and self.i == other.i and self.iphi == other.iphi
	
	def __ne__(self,other):  # Overload the "!=" test operator
		return not (self == other)
	
	"""
	def neutral_element_for_multiplication():
		return __class__(1)
	"""
	
	def rational_conjugate(self):
		return Ziphi(self.one + self.phi, -self.phi, self.i + self.iphi, -self.iphi)
	
	def complex_conjugate(self):
		return Ziphi(self.one, self.phi, -self.i, -self.iphi)
	
	def conjugate(self):
		return self.rational_conjugate() * self.complex_conjugate() * (self.rational_conjugate()).complex_conjugate()
	
	def norm(self):
		return self.i**4 + 2 * self.i**3 * self.iphi - self.i**2 * self.iphi**2 - 2 * self.i * self.iphi**3 + self.iphi**4 + 2 * self.i**2 * self.one**2 + 2 * self.i * self.iphi * self.one**2 + 3 * self.iphi**2 * self.one**2 + self.one**4 + 2 * self.i**2 * self.one * self.phi - 8 * self.i * self.iphi * self.one * self.phi - 2 * self.iphi**2 * self.one * self.phi + 2 * self.one**3 * self.phi + 3 * self.i**2 * self.phi**2 - 2 * self.i * self.iphi * self.phi**2 + 2 * self.iphi**2 * self.phi**2 - self.one**2 * self.phi**2 - 2 * self.one * self.phi**3 + self.phi**4
	
	def __pos__(self):	  # Overload the "+" unary operator
		return self
	
	def add(self,summand):
		sum_one = self.one + summand.one
		sum_phi = self.phi + summand.phi
		sum_i = self.i + summand.i
		sum_iphi = self.iphi + summand.iphi
		return Ziphi(sum_one, sum_phi, sum_i, sum_iphi)
	
	def __add__(self,summand):  # Overload the "+" binary operator
		if isinstance(summand, int):
			return Ziphi(self.one+summand, self.phi, self.i, self.iphi)
		else:
			return self.add(summand)
	
	def __radd__(self,summand):  # Overload the "+" binary operator
		if isinstance(summand, int):
			return Ziphi(self.one+summand, self.phi, self.i, self.iphi)
		else:
			return self.add(summand)
	
	def __iadd__(self,summand): # Overload the "+=" operator
		self = self + summand
		return self
	
	def __neg__(self):  # Overload the "-" unary operator
		return Ziphi(-self.one,-self.phi,-self.i,-self.iphi)
	
	def __sub__(self,summand):  # Overload the "-" binary operator
		return self.__add__(-summand)
	
	def __rsub__(self,summand):  # Overload the "-" binary operator
		if isinstance(summand, int):
			return Ziphi(summand-self.one, -self.phi, -self.i, -self.iphi)
		else:
			return summand-self
	
	def __isub__(self,summand): # Overload the "-=" operator
		self = self - summand
		return self
	
	def mult(self,multip): ###
		prod_one = - multip.i*self.i - multip.iphi*self.iphi + multip.one*self.one + multip.phi*self.phi
		prod_phi = - multip.iphi*self.i - multip.i*self.iphi - multip.iphi*self.iphi + multip.phi*self.one + multip.one*self.phi + multip.phi*self.phi
		prod_i = multip.one*self.i + multip.phi*self.iphi + multip.i*self.one + multip.iphi*self.phi
		prod_iphi = multip.phi*self.i + multip.one*self.iphi + multip.phi*self.iphi + multip.iphi*self.one + multip.i*self.phi + multip.iphi*self.phi
		return Ziphi(prod_one, prod_phi, prod_i, prod_iphi)
	
	def __mul__(self,multip):  # Overload the "*" operator
		if isinstance(multip, int):
			return Ziphi(self.one*multip, self.phi*multip, self.i*multip, self.iphi*multip)
		else:
			return self.mult(multip)
	
	def __rmul__(self,multip):  # Overload the "*" operator
		if isinstance(multip, int):
			return Ziphi(self.one*multip, self.phi*multip, self.i*multip, self.iphi*multip)
		else:
			return self.mult(multip)
	
	def __imul__(self,multip): # Overload the "*=" operator
		self = self * multip
		return self
	
	#"""
	def floordiv(self,divisor):
		if isinstance(divisor, 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 Zphi(candidate_r,candidate_i)
	
	def __floordiv__(self,divisor):  # Overload the "//" operator
		return self.floordiv(divisor)
	
	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
	"""
	
	def divmod(self,divisor):
		q = self//divisor
		return q, self - divisor * q
	#"""
	
	"""
	def xgcd(self,other):
		quot = Zphi()
		a1 = Zphi(1,0)
		b1 = Zphi(0,0)
		a2 = Zphi(0,0)
		b2 = Zphi(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 == Zphi(0,0)):
				return b, a2, b2
			quot = b // a
			b %= a
			a2 -= quot*a1
			b2 -= quot*b1
			if (b == Zphi()):
				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)
	"""
	
	"""
	def fancy_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
	"""
	def gcd(self,other):
		a = self
		b = other
		bn = b.norm()
		if a.norm() < bn:
			return b.gcd(a)
		if bn > 0:
			abbb = a*b.conjugate()
			w,x,y,z = abbb.tuple()
			q1,q2,q3,q4 = int(round(w*1.0/bn)), int(round(x*1.0/bn)), int(round(y*1.0/bn)), int(round(z*1.0/bn))
			s1,s2,s3,s4 = w*1.0/bn - q1, x*1.0/bn - q2, y*1.0/bn - q3, z*1.0/bn - q4
			i1,i2,i3,i4 = 3+int(round(6*s1)), 3+int(round(6*s2)), 3+int(round(6*s3)), 3+int(round(6*s4))
			delta = D[i1][i2][i3][i4]
			q1,q2,q3,q4 = q1-delta[0],q2-delta[1],q3-delta[2],q4-delta[3]
			q = Ziphi(q1,q2,q3,q4)
			a,b = a-q*b, b
			return a.gcd(b)
		return a
"""
x = Ziphi(1,1,1,1)
y = Ziphi(2,3,4,5)
z = Ziphi(2,3,4,6)
print x+y
print x*y
print x-y
print x.rational_conjugate()
print y.complex_conjugate()
print x.norm()
print x*x.conjugate()
print x.gcd(y)
print y.gcd(x)
print x.gcd(z)
print y.gcd(z)
"""

def twoSquareProduct(s1,s2,t1,t2):
	return s1*t1 - s2*t2, s1*t2 + s2*t1

def product(a):
	if len(a) == 0:
		return 1
	return a[-1]*product(a[:-1])

def twoSquares(self):
	if self == 0:
		return Zphi(0),Zphi(0)
	f = self.factor()
	# print f
	s1 = 1
	s2 = 0
	for x in list(f.keys()):
		# print s1, "and", s2
		y = Zphi(x[0],x[1])
		# print " y =", y
		p = abs(y.norm())
		if int(sqrt(p))**2 == p:
			p = int(sqrt(p))
		# print " p =", p
		m = f[x]
		
		if p == 2:
			(t1,t2) = (1,1)
			# print "  (t1,t2) =", (t1,t2)
			for j in range(m):
				s1,s2 = twoSquareProduct(s1,s2,t1,t2)
			# u1 = s1*t1 - s2*t2
			# u2 = s1*t2 + s2*t1
			# s1,s2 = u1,u2
			""" # i think this case should fall under "else" since y = sqrt(5) which can't be written as a sum of squares
		elif p == 5:
			t1 = Zphi(-1,2)
			t2 = 0
			print "  (t1,t2) =", (t1,t2)
			for j in range(m):
				s1,s2 = twoSquareProduct(s1,s2,t1,t2)
			# u1 = s1*t1 - s2*t2
			# u2 = s1*t2 + s2*t1
			# s1,s2 = u1,u2
			"""
		elif p%4 == 1:
			u = tonelli(p-1,p)
			# print u
			y_ext = Ziphi(y)
			# print y_ext
			u_mod = Ziphi(u,0,1,0)
			# print u_mod
			z = y_ext.gcd(u_mod)
			(t1,t2) = z.pair()
			# print "  (t1,t2) =", (t1,t2)
			# print "  t1**2 + t2**2 =", t1**2 + t2**2
			# print "  N(t1**2 + t2**2) =", (t1**2 + t2**2).norm()
			for j in range(m):
				s1,s2 = twoSquareProduct(s1,s2,t1,t2)
			# u1 = s1*t1 - s2*t2
			# u2 = s1*t2 + s2*t1
			# s1,s2 = u1,u2
		elif p%5 == 2 or p%5 == 3:
			u = tonelli(p-5,p)
			y_ext = Ziphi(y)
			u_mod = Ziphi(u,0,-1,2)
			z = y_ext.gcd(u_mod)
			(t1,t2) = z.pair()
			# print "  (t1,t2) =", (t1,t2)
			# print "  t1**2 + t2**2 =", t1**2 + t2**2
			# print "  N(t1**2 + t2**2) =", (t1**2 + t2**2).norm()
			for j in range(m):
				s1,s2 = twoSquareProduct(s1,s2,t1,t2)
			# u1 = s1*t1 - s2*t2
			# u2 = s1*t2 + s2*t1
			# s1,s2 = u1,u2
		else:
			# print "  y can't be written as a sum of squares"
			if m%2 == 0:
				s1 *= y**(m//2)
				s2 *= y**(m//2)
			else:
				s1,s2 = Zphi(None),Zphi(None)
				break
	# print "current s1 =", s1, "and s2 =", s2
	if not Zphi(s1).isNone:
		s = Zphi(s1*s1 + s2*s2)
		# print "s =", s
		prod = self*s.conjugate()
		# print "prod =", prod
		norm = s.norm()
		# print "norm =", norm
		quot = Zphi(prod.one/norm, prod.phi/norm)
		# print "quot =", quot
		n = unitDegree(quot)
		# print "n =", n
		if n == None or n%2 == 1:
			s1,s2 = Zphi(None),Zphi(None)
		elif n < 0:
			s1 *= Zphi(-1,1)**(-n//2)
			s2 *= Zphi(-1,1)**(-n//2)
		else:
			s1 *= Zphi(0,1)**(n//2)
			s2 *= Zphi(0,1)**(n//2)
	return s1, s2

"""
for n in range(0,50):
	x = Zphi(n)
	print x,
	# print x.factor()
	ts = twoSquares(x)
	print ts
	if not ts[0].isNone:
		print "  ", ts[0]**2+ts[1]**2
# looks good!
"""

"""
x = Zphi(74247,-45881)
print twoSquares(x)
"""