from convex import HyperplaneSolid
from rings import *
import numpy as np
# from math import cos, sin, sqrt, pi
import mpmath as mp
from mpmath import cos, sin, sqrt, pi, log, fabs
from quaternions import *

mp.mp.dps = 40
# print mp.mp


phi = (1+sqrt(5))/2
phidot = (1-sqrt(5))/2
p = 7+5*phi # maybe implement differently for more precision?
pdot = 12-5*phi
p_exact = Zphi(7,5)

# print [phi, phidot, p, pdot]

# bi-invariant metric on PU(2)
def metric(a,b):
	return 1-0.5*abs(np.trace(np.dot(a,b)))

def searchX1(t, eps, m):
	m *= mp.mpf(1)
	term1 = fabs(cos(t)) * eps * sqrt(2-eps**2)
	term2 = (1-eps**2)*sin(t)
	X1 = HyperplaneSolid(2)
	# X1.enter([1,phi],p**(m/2)) ### redundant w/ 5th ineq
	X1.enter([-1,-phi],p**(m/2))
	X1.enter([1,phidot],pdot**(m/2))
	X1.enter([-1,-phidot],pdot**(m/2))
	X1.enter([1,phi],p**(m/2)*(term1+term2))
	X1.enter([-1,-phi],p**(m/2)*(term1-term2))
	X1.enter([sin(t),sin(t)*phi],p**(m/2)*(1-eps**2))
	# print X1
	pts = X1.Znpoints()
	return pts # [pt for pt in pts if X1.check(pt)]

"""
for m in range(10,20):
	print "==================== m =", m, "===================="
	n,pts = searchX1(pi/8, mp.mpf(10.0**(-10)), m)
	if len(pts) > 0:
		print "n was", n
		print pts
		if len(pts) > 3:
			break
"""

def searchX0(t, eps, m, pt):
	x1 = pt[0] + pt[1]*phi
	x1dot = pt[0] + pt[1]*phidot
	m *= mp.mpf(1)
	X0 = HyperplaneSolid(2)
	X0.enter([1,phi],sqrt(p**m - x1**2))
	X0.enter([-1,-phi],sqrt(p**m - x1**2))
	X0.enter([1,phidot],sqrt(pdot**m - x1dot**2))
	X0.enter([-1,-phidot],sqrt(pdot**m - x1dot**2))
	X0.enter([cos(t),cos(t)*phi],p**(m/2)-x1*sin(t))
	X0.enter([-cos(t),-cos(t)*phi],-p**(m/2)*(1-eps**2)+x1*sin(t))
	# print X0
	pts = X0.Znpoints()
	return pts # [pt for pt in pts if X0.check(pt)]

"""
X1points = [[-1246585, -2022220], [-1247195, -2021843], [-1248182, -2021233], [-1249169, -2020623]]
m = 12
for pt in X1points:
	print "==================== pt =", pt, "===================="
	n,pts = searchX0(pi/8, mp.mpf(10.0**(-10)), m, pt)
	if len(pts) > 0:
		print "n was", n
		print pts
"""


"""
t is an angle corresponding to the diagonal matrix
/  it	   \
| e	0	|
|	   -it |
\ 0	e	/
"""


def diag(t, eps, start_m = 1, stop_m = 0):
	if stop_m < start_m:
		stop_m = 2*log(1/eps**3)/log(59)
	m = start_m - 1
	stopQ = False
	ticks = 0 # number of times we give up on comp for having too large factors
	while not stopQ:
		m += 1
		print("==================== m =", m, "====================")
		points1 = searchX1(t, eps, m)
		print("there are", len(points1), "candidate point(s) for x1")
		for pt1 in points1:
			points0 = searchX0(t, eps, m, pt1)
			if len(points0) > 0:
				print(pt1, "gives rise to", len(points0), "candidate point(s) for x0")
				for pt0 in points0:
					x0,x1 = Zphi(pt0[0],pt0[1]),Zphi(pt1[0],pt1[1])
					comp = p_exact**m - x0**2 - x1**2
					if np.all(np.array(list(factorint(comp.norm()).keys())) < 50000):
						(x2,x3) = twoSquares(comp)
						if not x2.isNone:
							# print "comp =", comp
							# print "goal (numerical) =", p**m
							# print "goal (exact) =", p_exact**m
							# print "goal?", x0**2 + x1**2 + x2**2 + x3**2
							print("found it!")
							x = HZphi(x0,x1,x2,x3)
							print(ticks, "ticks")
							return TeXGamma(factorGamma(x))
					else:
						ticks += 1
		if m >= stop_m:
			print(ticks, "ticks")
			return ""

# print diag(pi/8, sqrt(2)*mp.mpf(10.0**(-8)))
# looks good!

# print diag(mp.mpf(0.1602149635879141605714577864586057804666), sqrt(2)*mp.mpf(10**(-6)), 3)

# print diag(pi/8, sqrt(2)*mp.mpf(10.0**(-10)))
### current goal 

"""
s = HyperplaneSolid(2)
s.enter([mp.mpf(1),mp.mpf(0)],2)
s.enter([mp.mpf(-1),mp.mpf(0)],0)
s.enter([mp.mpf(0),mp.mpf(1)],3)
s.enter([mp.mpf(0),mp.mpf(-1)],0)
print s.Znpoints()
# looks okay (other than rounding errors at corners (e17)
"""