### Ancillary SageMath code for the paper "Geometric transition from hyperbolic to Anti-de Sitter structures in dimension four" by S. Riolo and A. Seppi, to appear in Ann. Sc. Norm. Super. Pisa Cl. Sci.


# We will work in the affine chart A^4 with its affine coordinates y1,...,y4:

var('y1, y2, y3, y4, t')
Variables = [y1, y2, y3, y4]
x = vector([1, y1, y2, y3, y4])


# The  polytope \bar{Q_t} is defined in Lemma 6.6. The bounding hyperplanes of r_{|t|}(\bar{Q_t}) are obtained from Table 3 and Equation (11).  Here they are:

Vectors = [
vector([ -sqrt(2), +1, +1, +1, +1 ]),
vector([ -sqrt(2), +1, +1, -1, -1 ]),
vector([ -sqrt(2), +1, +1, +1, +t*t ]),
vector([ -sqrt(2), +1, +1, -1, -t*t ]),
vector([ -1, +sqrt(2), 0, 0, 0 ]),
vector([0,-1,1,0,0]),
vector([0,0,-1,1,0]),
vector([0,0,-1,-1,0]) ]


# The following is the main part of the computation, i.e. parts (1)--(4) in the proof of Lemma 6.8. We use Rouché--Capelli. FlagSanity is needed to check that we are avoiding a Sage bug.

import itertools
Solutions = []
ell = len(Vectors)
SanityFlag = 0
for l in range(0,ell-3):
	Solutions_l = []
	for combination in itertools.combinations(range(ell),ell-l):
		Equations = [x*Vectors[i] for i in combination]
		Rows = [Vectors[i] for i in combination]
		MatCompl = Matrix(Rows)
		rankC = MatCompl.rank()
		MatInc = MatCompl.delete_columns([0])
		rankI = MatInc.rank()
		if rankC == rankI and rankC == 4:
			point = x.substitute(solve(Equations,Variables)[0])
			flag = 0
			with assuming(t<0, t>-1):
				for v in Vectors:
					prod = point*v
					if prod == 0:
						prod = 0
					else:
						prod = prod.factor()
					if{bool(prod > 0), bool(prod <= 0)}!={True, False}:
						SanityFlag = 1
					if prod > 0:
						flag = 1
			if flag == 0:
				Solutions_l.append( [list(combination), point] )
	if Solutions_l != []:
		Solutions.append(Solutions_l)
if SanityFlag != 0:
	print('There are problems')
else:
	for list in Solutions:
		l = len(list[0][0])
		print(str(len(list)) + ' points of ' + str(l) + '-ple intersection')


# We eliminate the reundant solutions:
	
ideal_vertex = Solutions[0][0]
Vertices = [ideal_vertex]
for Sol_l in Solutions:
	for Sol in Sol_l:
		if Sol[1] != ideal_vertex[1]:
			Vertices.append(Sol)
print('The polytope has ' + str(len(Vertices)) + ' vertices')


# We check that all found vertices are pairwise distinct:

Vertices2 = []
for Sol_l in Solutions:
	for Sol in Sol_l:
		Vertices2.append(Sol)
ActualVertices = []
for i in range(len(Vertices2)):
	flag = 0
	for j in range(i):
		if Vertices2[j][1] == Vertices2[i][1]:
			flag = 1
	if flag == 0:
		ActualVertices.append(Vertices2[i])
Vertices2 = ActualVertices
if Vertices2 == Vertices:
	print('Check ok')
else:
	print('There are repetitions')
	

# We rescale back to AdS^4 and check that all the 12 finite vertices are in AdS^4:
	
Scale = diagonal_matrix([1,1,1,1,t])
Form = diagonal_matrix([-1,1,1,1,-1])
flag_sanity = 0
flag = 0
ideal_vertices_counter = 0
for v in Vertices:
	prod = Scale*v[1]*Form*Scale*v[1]
	if prod == 0:
		prod = 0
		ideal_vertices_counter += 1
	else:
		prod = prod.factor()
		with assuming(t<0, t>-1):
			if{bool(prod > 0), bool(prod <= 0)}!={True, False}:
				flag_sanity = 1
			if prod > 0:
				flag = 1
if flag_sanity == 0 and flag == 0 and ideal_vertices_counter == 1:
	print('Check ok')
else:
	print('There are problems')


# To conclude, we print the vertices with their names:

def Name(v):
	if v == Vectors[0]:
		return('0+')
	if v == Vectors[1]:
		return('3+')
	if v == Vectors[2]:
		return('0-')
	if v == Vectors[3]:
		return('3-')
	if v == Vectors[4]:
		return('A')
	if v == Vectors[5]:
		return('L')
	if v == Vectors[6]:
		return('M')
	if v == Vectors[7]:
		return('N')
def SolName(sol):
	name = ' '
	for i in sol[0]:
		name += Name(Vectors[i]) + ' '
	return(name)
def PrintSol(sol_list):
	for p in sol_list:
		print(SolName(p) + ':  '  + str(p[1]))

PrintSol(Vertices)