From 3e7361c0d5ace9129fd44d93d651f2c8f63afb73 Mon Sep 17 00:00:00 2001 From: Sam Freedman Date: Thu, 15 Dec 2022 16:57:29 -0500 Subject: [PATCH 1/4] create periodic_points.py --- flatsurf/geometry/periodic_points.py | 72 ++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) create mode 100644 flatsurf/geometry/periodic_points.py diff --git a/flatsurf/geometry/periodic_points.py b/flatsurf/geometry/periodic_points.py new file mode 100644 index 000000000..cdfe05991 --- /dev/null +++ b/flatsurf/geometry/periodic_points.py @@ -0,0 +1,72 @@ +###################################################################### +# This file is part of sage-flatsurf. +# +# Copyright (C) 2022 Julian Rüth +# 2022 Sam Freedman +# +# sage-flatsurf is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# +# sage-flatsurf is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with sage-flatsurf. If not, see . +###################################################################### + +from thurston_veech import ThurstonVeech + +class PeriodicPoints(): + def __init__(self, T): + # we will need hp, vp, hmult, vmult + pass + + ''' + given a rectangle (H cap V), write down the constraints in that rectangle + assuming that c(H)/h(V) is irrational + we'd need: + all the rectangles in the horizontal cylinders + global horizontal multitwist + multiplicities of the horizontal cylinder + + y / ht(H) in Q + x / ht(V) in Q + (x + ay - d) / ht(V') in Q, + where d = sum ht(V_k) adding from V --> V_k + and T = [1 a | 0 1] + Note: the left-hand sides are elements of K[x, y] + + pi_j : (K = Q^n) --> Q + pi_j (y / ht(H)) = 0 for all 1 <= j <= n + pi_j (x / ht) ... + pi_j (...) + + each constraint --> k1 x + k2 y + k3 in Q + and now project k1, k2, k3 to get coefficient for each number field basis elt. + + let alpha_1, .. alpha_n be a Q-basis for K + then + x = x^i alpha_i + k = a^i alpha_i + (a^i alpha_i) * (x^j alpha_j) + (b^i alpha_i) * (y^j alpha_j) + (c^i alpha_i) in Q + + y / ht(H) in Q + + e.g. (y^1 + y^2 rt2) / (a^1 + a^2 rt2) in Q + [(y1 + y2 rt2) * (a1 - a2rt2)] / (a1**2 - 2 a2**2) + [a1 y1 - a2y1rt2 + y2 a1 rt2 - 2 a2y2] / (a1**2 - 2 a2**2) in Q + (a1 y1 - 2 a2 y2)/(a1**2 - 2 a2**2) + rt2 (y2 a1 - a2 y1)/(a1**2 - 2 a2**2) in Q + for each element of the basis that is not rational, get equations + rt2 equation is ==> (a1 y2 - a2 y1)/(a1**2 - 2 a2**2) = 0 + i.e. [a1/(a1**2 - 2 a2**2)] y2 - [a2/(a1**2 - 2 a2**2)] y1 = 0 + + now solve the resulting system of equations + ''' + + + + From aaf1613beccdc654f531a9f7815ac5f34fcfa59f Mon Sep 17 00:00:00 2001 From: Sam Freedman Date: Sun, 29 Jan 2023 20:01:25 -0500 Subject: [PATCH 2/4] start implementing algorithm --- flatsurf/geometry/periodic_points.py | 97 ++++++++++++++++++++++++++-- 1 file changed, 90 insertions(+), 7 deletions(-) diff --git a/flatsurf/geometry/periodic_points.py b/flatsurf/geometry/periodic_points.py index cdfe05991..898771eef 100644 --- a/flatsurf/geometry/periodic_points.py +++ b/flatsurf/geometry/periodic_points.py @@ -23,8 +23,91 @@ class PeriodicPoints(): def __init__(self, T): # we will need hp, vp, hmult, vmult + self._surface = T pass + def list(self): + r'''Return the periodic points''' + + G = self._candidate_graph() + G = self._prune_candidate_graph(G) + return G.vertices() + + def _constraint_segments_in_region(self, R): + r'''Return segments containing all periodic points in region ``R`` ''' + segments = [] + + width, height = R.width(), R.height() + horizontal_cylinder = R.horizontal_cylinder() + vertical_cylinder = R.vertcial_cylinder() + + horizontal_constraint = self._horizontal_constraint(height) + vertical_constraint = self._vertical_constraint(width) + + from sage.all import QQ + if (horizontal_cylinder.circumference() / width) not in QQ: + for R1 in horizontal_cylinder.regions(start=R, multiplicity=horizontal_cylinder.multiplicity()): + connecting_constraint = self._connecting_equation(R, R1) + segments.append(self._solve_constraints(horizontal_constraint, vertical_constraint, connecting_constraint)) + + else: + assert (vertical_cylinder.circumference() / height) not in QQ + for R1 in vertical_cylinder.regions(start=R, multiplicity=vertical_cylinder.multiplicity()): + connecting_constraint = self._connecting_equation(R, R1) + segments.append(self._solve_constraints(horizontal_constraint, vertical_constraint, connecting_constraint)) + + def _constraint_segments(self): + r'''Return set of segments containing all periodic points''' + return {l + for R in self._surface.regions() + for l in self._constraint_segments_in_region(R)} + + def _candidate_points(self): + r'''Return a (potentially large) set of candidate periodic points''' + + candidate_points = set() + S = self._constraint_segments() + + while True: + g = self._surface.veech_group.random_element() + S1 = {l.apply_matrix(g) for l in S} + if {l.slope() for l in S}.isdisjoint({l1.slope() for l1 in S1}): + break + + for l in S: + for l1 in S1: + candidate_points.add(l.intersection(l1)) + + return candidate_points + + def _candidate_graph(self): + r'''Return graph with: + vertices being either a set S of potential periodic points or ``None``, and + p <--> q iff there is a generator of the Veech Group sending p to q, + where p <--> ``None`` when a generator g sends p outside of S + ''' + S = self._candidate_points() + S.add(None) + + G = Graph() + + for p in S: + for gen in self._surface.veech_group.gens(): + q = p.apply_matrix(gen) + if q in S: + G.add_edge(p, q) + else: + G.add_edge(p, None) + + return G + + def _prune_candidate_graph(self, G): + r'''Given candidate graph, return the subgraph of periodic points''' + + # remove the connected component of ``None`` + raise NotImplementedError + + ''' given a rectangle (H cap V), write down the constraints in that rectangle assuming that c(H)/h(V) is irrational @@ -43,16 +126,16 @@ def __init__(self, T): pi_j : (K = Q^n) --> Q pi_j (y / ht(H)) = 0 for all 1 <= j <= n pi_j (x / ht) ... - pi_j (...) + pi_j (...) each constraint --> k1 x + k2 y + k3 in Q and now project k1, k2, k3 to get coefficient for each number field basis elt. let alpha_1, .. alpha_n be a Q-basis for K - then - x = x^i alpha_i - k = a^i alpha_i - (a^i alpha_i) * (x^j alpha_j) + (b^i alpha_i) * (y^j alpha_j) + (c^i alpha_i) in Q + then + x = x^i alpha_i + k = a^i alpha_i + (a^i alpha_i) * (x^j alpha_j) + (b^i alpha_i) * (y^j alpha_j) + (c^i alpha_i) in Q y / ht(H) in Q @@ -64,9 +147,9 @@ def __init__(self, T): rt2 equation is ==> (a1 y2 - a2 y1)/(a1**2 - 2 a2**2) = 0 i.e. [a1/(a1**2 - 2 a2**2)] y2 - [a2/(a1**2 - 2 a2**2)] y1 = 0 - now solve the resulting system of equations + now solve the resulting system of equations ''' - + From f21656aded7df6379b8fcfb747a046313f50e343 Mon Sep 17 00:00:00 2001 From: Sam Freedman Date: Tue, 31 Jan 2023 15:44:13 -0500 Subject: [PATCH 3/4] add constraint equations --- flatsurf/geometry/periodic_points.py | 98 ++++++++++++++-------------- 1 file changed, 49 insertions(+), 49 deletions(-) diff --git a/flatsurf/geometry/periodic_points.py b/flatsurf/geometry/periodic_points.py index 898771eef..2b98f94aa 100644 --- a/flatsurf/geometry/periodic_points.py +++ b/flatsurf/geometry/periodic_points.py @@ -18,7 +18,7 @@ # along with sage-flatsurf. If not, see . ###################################################################### -from thurston_veech import ThurstonVeech +from flatsurf.geometry.thurston_veech import ThurstonVeech class PeriodicPoints(): def __init__(self, T): @@ -33,6 +33,52 @@ def list(self): G = self._prune_candidate_graph(G) return G.vertices() + def _horizontal_constraint(self, height): + r'''Constrain y coordinate of point with rational height lemma + + EXAMPLE:: + + sage: import flatsurf as fs + sage: X = fs.translation_surfaces.mcmullen_genus2_prototype(1, 1, 0, -1) + sage: P = PeriodicPoints(X) + sage: R = X.base_ring() + sage: P._horizontal_constraint(R(R.gen())) + ([0 0 1 0], (0)) + ''' + + A = (~height).matrix().submatrix(1) + A = A.parent().zero().augment(A) + return A, A.column_space().zero() + + def _vertical_constraint(self, width): + r'''Constrain x coordinate of point with rational height lemma + + EXAMPLE:: + + sage: import flatsurf as fs + sage: X = fs.translation_surfaces.mcmullen_genus2_prototype(1, 1, 0, -1) + sage: P = PeriodicPoints(X) + sage: R = X.base_ring() + sage: P._vertical_constraint(R(R.gen())) + ([1 0 0 0], (0)) + + ''' + A = (~width).matrix().submatrix(1) + A = A.augment(A.parent().zero()) + return A, A.column_space().zero() + + def _connecting_constraint_horizontal(self, R, distance): + return self._connecting_constraint_generic(self._surface.horizontal_twist()[0][1], distance, R.width()) + + def _connecting_constraint_vertical(self, R, distance): + return self._connecting_constraint_generic(self._surface.vertical_twist()[0][1], distance, R.height()) + + def _connecting_constraint_generic(self, twist, distance, side_length): + A_t = twist.matrix() + A = (~side_length).matrix() * A_t.parent().one().augment(A_t) + v = (~side_length).matrix() * distance.vector() + return A.submatrix(1), v[1:] + def _constraint_segments_in_region(self, R): r'''Return segments containing all periodic points in region ``R`` ''' segments = [] @@ -46,8 +92,8 @@ def _constraint_segments_in_region(self, R): from sage.all import QQ if (horizontal_cylinder.circumference() / width) not in QQ: - for R1 in horizontal_cylinder.regions(start=R, multiplicity=horizontal_cylinder.multiplicity()): - connecting_constraint = self._connecting_equation(R, R1) + for R1, wrap in horizontal_cylinder.regions(start=R, multiplicity=horizontal_cylinder.multiplicity()): + connecting_constraint = self._connecting_constraint(R1, wrap) segments.append(self._solve_constraints(horizontal_constraint, vertical_constraint, connecting_constraint)) else: @@ -107,49 +153,3 @@ def _prune_candidate_graph(self, G): # remove the connected component of ``None`` raise NotImplementedError - - ''' - given a rectangle (H cap V), write down the constraints in that rectangle - assuming that c(H)/h(V) is irrational - we'd need: - all the rectangles in the horizontal cylinders - global horizontal multitwist - multiplicities of the horizontal cylinder - - y / ht(H) in Q - x / ht(V) in Q - (x + ay - d) / ht(V') in Q, - where d = sum ht(V_k) adding from V --> V_k - and T = [1 a | 0 1] - Note: the left-hand sides are elements of K[x, y] - - pi_j : (K = Q^n) --> Q - pi_j (y / ht(H)) = 0 for all 1 <= j <= n - pi_j (x / ht) ... - pi_j (...) - - each constraint --> k1 x + k2 y + k3 in Q - and now project k1, k2, k3 to get coefficient for each number field basis elt. - - let alpha_1, .. alpha_n be a Q-basis for K - then - x = x^i alpha_i - k = a^i alpha_i - (a^i alpha_i) * (x^j alpha_j) + (b^i alpha_i) * (y^j alpha_j) + (c^i alpha_i) in Q - - y / ht(H) in Q - - e.g. (y^1 + y^2 rt2) / (a^1 + a^2 rt2) in Q - [(y1 + y2 rt2) * (a1 - a2rt2)] / (a1**2 - 2 a2**2) - [a1 y1 - a2y1rt2 + y2 a1 rt2 - 2 a2y2] / (a1**2 - 2 a2**2) in Q - (a1 y1 - 2 a2 y2)/(a1**2 - 2 a2**2) + rt2 (y2 a1 - a2 y1)/(a1**2 - 2 a2**2) in Q - for each element of the basis that is not rational, get equations - rt2 equation is ==> (a1 y2 - a2 y1)/(a1**2 - 2 a2**2) = 0 - i.e. [a1/(a1**2 - 2 a2**2)] y2 - [a2/(a1**2 - 2 a2**2)] y1 = 0 - - now solve the resulting system of equations - ''' - - - - From 647b5f479f68d40fd85eb77c55ae7e5f49069762 Mon Sep 17 00:00:00 2001 From: Sam Freedman Date: Thu, 2 Feb 2023 13:54:38 -0500 Subject: [PATCH 4/4] implemenet solve_constraints --- flatsurf/geometry/periodic_points.py | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/flatsurf/geometry/periodic_points.py b/flatsurf/geometry/periodic_points.py index 2b98f94aa..982aea7ee 100644 --- a/flatsurf/geometry/periodic_points.py +++ b/flatsurf/geometry/periodic_points.py @@ -79,6 +79,18 @@ def _connecting_constraint_generic(self, twist, distance, side_length): v = (~side_length).matrix() * distance.vector() return A.submatrix(1), v[1:] + def _solve_constraints_in_region(self, R, constraints): + M = matrix([c[0] for c in constraints]) + b = vector([c[1] for c in constraints]) + x = M.solve_right(b) + K = M.right_kernel() + + assert K.dimension() <= 1 + + from sage.all import Polyhedron + line = Polyhedron(vertices=[x], lines=K.basis()) + return R.intersection(line) + def _constraint_segments_in_region(self, R): r'''Return segments containing all periodic points in region ``R`` ''' segments = [] @@ -94,7 +106,8 @@ def _constraint_segments_in_region(self, R): if (horizontal_cylinder.circumference() / width) not in QQ: for R1, wrap in horizontal_cylinder.regions(start=R, multiplicity=horizontal_cylinder.multiplicity()): connecting_constraint = self._connecting_constraint(R1, wrap) - segments.append(self._solve_constraints(horizontal_constraint, vertical_constraint, connecting_constraint)) + constraints = [horizontal_constraint, vertical_constraint, connecting_constraint] + segments.append(self._solve_constraints_in_region(R, constraints)) else: assert (vertical_cylinder.circumference() / height) not in QQ @@ -102,6 +115,9 @@ def _constraint_segments_in_region(self, R): connecting_constraint = self._connecting_equation(R, R1) segments.append(self._solve_constraints(horizontal_constraint, vertical_constraint, connecting_constraint)) + return [s for s in segments if not s.is_empty()] + + def _constraint_segments(self): r'''Return set of segments containing all periodic points''' return {l