import numpy as np
import mpmath as mp
from mpmath import cos, sin, atan2, sqrt, pi, sign
from diagonal import *
from rings import *
from quaternions import *

mp.mp.dps = 40
print(mp.mp)

phi = (1+sqrt(5))/2
phidot = (1-sqrt(5))/2
p = 7+5*phi
pdot = 7+5*phidot
p_exact = Zphi(7,5)

print(phi, phidot, p, pdot, p_exact)

# inputs will be characterized by the top row of the matrix (four real numbers); equivalently, the canonical quaternion. not necessarily normalized
def normalize(a,b,c,d):
	det = a**2 + b**2 + c**2 + d**2
	assert det != 0, "the 0 matrix is not in the group"
	s = sqrt(det)
	return a/s,b/s,c/s,d/s

def angle(re, im):
	if isinstance(re, Zphi):
		if re == 0:
			if im.float() > 0:
				return pi/2
			return -pi/2
		return atan2(im.float(), re.float())
	if re == 0:
		if im > 0:
			return pi/2
		return -pi/2
	return atan2(im, re)

def approx(a,b,c,d, eps, start_m = 1, stop_m = 0):
	if stop_m < start_m:
		stop_m = log(1/eps**3)/log(59)
	# p = 7+5*phi
	a,b,c,d = normalize(a,b,c,d)
	print(a, b, c, d)
	ga = sqrt(a**2 + b**2)
	m = start_m - 1
	stopQ = False
	while not stopQ:
		m += 1
		print("==================== m =", m, "====================")
		P = HyperplaneSolid(2)
		P.enter([1,phi], eps**1 * ga * p**m + ga**2 * p**m)
		P.enter([-1,-phi], eps**1 * ga * p**m - ga**2 * p**m)
		P.enter([1,phidot], pdot**m)
		P.enter([-1,-phidot], pdot**m)
		# print P
		pts = P.Znpoints()
		print("considering", len(pts), "points")
		for pt in pts:
			# print pt
			z = Zphi(pt[0],pt[1])
			if np.all(np.array(list(factorint(z.norm()).keys())) < 1000000):
				(x0,x1) = twoSquares(z)
				# print pt
				# print "  (", x0, ",", x1, ")"
				if not x0.isNone:
					comp = p_exact**m - z
					if np.all(np.array(list(factorint(comp.norm()).keys())) < 1000000):
						(x2,x3) = twoSquares(comp)
						print(pt)
						print("  (", x0, ",", x1, ")")
						print("  (", x2, ",", x3, ")")
						if not x2.isNone:
							a1, b1, a2, b2 = angle(a,b), angle(c,d), angle(x0,x1), angle(x2,x3)
							print(a1, b1, a2, b2)
							t1, t2 = 0.5*(a1 - a2 + b1 - b2), 0.5*(a1 - a2 - b1 + b2)
							print(t1, t2)
							# return x0, x1, x2, x3
							middle = TeXGamma(factorGamma(HZphi(x0, x1, x2, x3)))
							print("middle =", middle)
							front = diag(t1, sqrt(2)*eps, m//2)
							print("front =", front)
							end = diag(t2, sqrt(2)*eps, m//2)
							print("end =", end)
							return front + "\n" + middle + "\n" + end
		if m >= stop_m:
			break

print("\nFACTORIZATION \n" + approx(0, 1, 0, 1, mp.mpf(10**(-6)), 3))