From 4a4469fa6908419cd806f0a3689c2442ebde196d Mon Sep 17 00:00:00 2001 From: Alexander Galarraga Date: Tue, 4 Aug 2020 10:59:08 -0400 Subject: [PATCH 1/3] initial commit --- src/doc/en/reference/references/index.rst | 3 + .../arithmetic_dynamics/berkovich_ds.py | 175 +++++++++++++++++- 2 files changed, 177 insertions(+), 1 deletion(-) diff --git a/src/doc/en/reference/references/index.rst b/src/doc/en/reference/references/index.rst index 32a0d4053ca..8b04dfdfd78 100644 --- a/src/doc/en/reference/references/index.rst +++ b/src/doc/en/reference/references/index.rst @@ -4725,6 +4725,9 @@ REFERENCES: .. [Rud1958] \M. E. Rudin. *An unshellable triangulation of a tetrahedron*. Bull. Amer. Math. Soc. 64 (1958), 90-91. + +.. [Rum2013] Robert Rumely. *The Minimal Resultant Locus*. Acta + Arithmetica, Vol. 169 (2015). :arxiv:`1304.1201` .. [Rus2003] Frank Ruskey. *Combinatorial Generation*. (2003). http://www.1stworks.com/ref/ruskeycombgen.pdf diff --git a/src/sage/dynamics/arithmetic_dynamics/berkovich_ds.py b/src/sage/dynamics/arithmetic_dynamics/berkovich_ds.py index 62cf5e28ff8..ba4122f7de0 100644 --- a/src/sage/dynamics/arithmetic_dynamics/berkovich_ds.py +++ b/src/sage/dynamics/arithmetic_dynamics/berkovich_ds.py @@ -42,6 +42,7 @@ from sage.categories.number_fields import NumberFields from sage.rings.rational_field import QQ from sage.rings.infinity import Infinity +from sage.matrix.constructor import Matrix @add_metaclass(InheritComparisonClasscallMetaclass) class DynamicalSystem_Berkovich(Element): @@ -704,7 +705,6 @@ def __call__(self, x, type_3_pole_check=True): raise NotImplementedError('action on Type IV points not implemented') f = self._system if x.type_of_point() == 2: - from sage.matrix.constructor import Matrix from sage.modules.free_module_element import vector if self.domain().is_number_field_base(): ideal = self.domain().ideal() @@ -836,6 +836,179 @@ def __call__(self, x, type_3_pole_check=True): new_radius = max(new_radius, p**(-valuation/prime.absolute_ramification_index())*r**i) return self.domain()(new_center, new_radius) + def min_res_locus(self, ret_conjugation=False): + r""" + Returns the minimal resultant locus and the minimal order of the resultant. + + Let the action of this map `\phi` be given by + `[X : Y] \to [F(X,Y) : G(X,Y)]`, where `F` and `G` are normalized. + The resultant of `\phi` is the resultant of `F` and `G`. + Define the function `\text{ordRes}(\phi) = \text{ord}(\text{Res}(F, G))`. + `\text{ordRes}` is not constant under conjugacy, hence there is + some pair `(F, G)` for which `\text{ordRes}` achieves a minimum. + + The minimal resultant locus is the segment of Berkovich + space where the minimal value of the resultant is achieved. + See [Rum2013]_ for details. + + INPUT: + + - ``ret_conjugation`` -- (default: ``False``) If the conjugation to achieve + the minimal resultant should be returned. + + OUTPUT: + + - If ``ret_conjugation`` is ``False``, then a tuple (``start``, ``end``) is returned, + where ``start`` and ``end`` are points of Berkovich space. The minimal resultant + is then achieved on the interval defined by [``start``, ``end``]. + + - If ``ret_conjugation`` is ``True``, then a tuple ((``start``, ``end``), ``conj``) + is returned, where [``start``, ``end``] defines the minimal resultant locus + and ``conj`` is a matrix such that this dynamical system achieves the minimal + order of the resultant after conjugation by ``conj``. + + ALGORITHM: + + See [Rum2013]_. + """ + system = self._system + QQ_prime = self.domain().prime() + base_prime_ideal = self.domain().ideal() + affine_system = system.dehomogenize(1) + base = affine_system.base_ring() + T = base['w'] + w = T.gens()[0] + var = affine_system[0].parent().gens()[0] + phi = affine_system[0].numerator().parent().hom([w]) + num, dem = affine_system[0].numerator(), affine_system[0].denominator() + dem = phi(dem) + num = phi(num) + #print('num:', num) + #print('dem:', dem) + fixed_point_equation = (num-w*dem) + factorization = list(fixed_point_equation.factor()) + if(system([1,0]) == system.domain()([1,0])): + Q = dem + else: + value = system([1,0]).dehomogenize(1)[0] + Q = num - value*dem + factorization += list(Q.factor()) + min_list = [] + print('all factors:', factorization) + for i in range(len(factorization)): + irreducible = factorization[i][0] + extension = base.extension(irreducible, 's%s'%i) + new_primes = extension.primes_above(base_prime_ideal) + for new_prime in new_primes: + ramification_index = new_prime.absolute_ramification_index() + valuation = lambda x: x.valuation(new_prime)/ramification_index + s = extension.gens()[0] + transformation_matrix = Matrix([[1, s],[0, 1]]) + new_system = system.change_ring(extension) + new_system = new_system.conjugate(transformation_matrix) + new_system.normalize_coordinates() #very long call + #print('res = ', new_system.resultant()) + res = valuation(new_system.resultant()) + #print('ordres = ', res) + d = new_system.degree() + #print('system = ', new_system) + C_list = [] + D_list = [] + for i in range(d+1): + C_i = new_system[0][i,d-i] + D_i = new_system[1][i,d-i] + #print('C_%s = ' %i, C_i) + #print('val = ', valuation(C_i)) + if C_i != 0: + C_list.append(QQ(res-2*d*valuation(C_i))) + else: + C_list.append(None) + if D_i != 0: + D_list.append(QQ(res-2*d*valuation(D_i))) + else: + D_list.append(None) + #print('C_list =', C_list) + D = QQ['t'] + t = D.gens()[0] + C_lines = [(C_list[i] + (QQ(d^2+d-2*d*i))*t) for i in range(len(C_list)) if C_list[i] != None] + D_lines = [(D_list[i] + (QQ(d^2+d-2*d*(i+1)))*t) for i in range(len(D_list)) if D_list[i] != None] + all_lines = C_lines + D_lines #crude minimization + minimum = None + #print(C_lines,D_lines) + all_intersections = [] + for line_1 in all_lines: + for line_2 in all_lines: + if line_1 != line_2 and line_2[1] != line_1[1]: + intersection = QQ((line_1[0]-line_2[0])/(line_2[1]-line_1[1])) + if intersection not in all_intersections: + all_intersections.append(intersection) + #print(all_intersections) + for intersection in all_intersections: + C_max = max([(i(intersection)) for i in C_lines]) + D_max = max([(i(intersection)) for i in D_lines]) + #print('intersection = ', intersection) + #print('C_max = ', C_max) + #print('D_max = ', D_max) + #print('value of X(i) = ', max(C_max, D_max)) + X_i = (max(C_max,D_max)) + if X_i < minimum or minimum == None: + minimum = (X_i) + intersection_final = QQ(intersection) + #print('minimum: ', minimum) + interval_check = False + for line in all_lines: #check if minimum is achieved on an interval + if line.is_constant(): + if line[0] == minimum: + interval_check = True + break + min_list.append((minimum, intersection_final, interval_check, extension, new_prime)) + minimizing_list = [i[0] for i in min_list] + minimum = min(minimizing_list) + min_index = minimizing_list.index(minimum) + min_res = min_list[min_index][0] + intersection_final = min_list[min_index][1] + interval_check = min_list[min_index][2] + extension = min_list[min_index][3] + prime_above = min_list[min_index][4] + if not interval_check: #minimum is not achieved on an interval + B = Berkovich_Cp_Projective(extension, prime_above) + min_res_locus = B(extension.gens()[0], power=-1*intersection_final) + min_res_locus = min_list[min_index][0] + if not ret_conjugation: + return min_res_locus + else: + #print('intersection_final:', intersection_final) + root = intersection_final.denominator() + if root == 1: + #print('in no denominator case') + #print('extension:', extension) + conj_matrix = Matrix([[QQ_prime**(intersection_final), extension.gens()[0]], [0,1]]) + new_system = system.change_ring(extension) + else: + C = extension['c'] + c = C.gens()[0] + defining_polynomial = c**(root) - QQ_prime**(intersection_final.numerator()) + factorization = defining_polynomial.factor() + degree = None + for factor in factorization: + if degree == None or factor[0].degree() < degree: #slow symbolic calcuation + final_factor = factor[0] + degree = factor[0].degree() + if degree != 1: + new_extension = extension.extension(final_factor, 'w0') + A = new_extension.gens()[0] + B = new_extension.gens()[1] + else: + A = final_factor[0] + B = extension.gens()[0] + conj_matrix = Matrix([[A, B],[0,1]]) + new_system = system.change_ring(new_extension) + new_system = new_system.conjugate(conj_matrix) + return (min_res_locus, conj_matrix, new_system) + else: + print('not implemented') + return [min_list[min_index][0], min_list[min_index][1]] + class DynamicalSystem_Berkovich_affine(DynamicalSystem_Berkovich): r""" A dynamical system of the affine Berkovich line over `\CC_p`. From 519b0cca8563a01cc70448cd94f70ea70d3d51a0 Mon Sep 17 00:00:00 2001 From: Alexander Galarraga Date: Tue, 4 Aug 2020 11:45:56 -0400 Subject: [PATCH 2/3] fixed power error --- .../arithmetic_dynamics/berkovich_ds.py | 29 +++++++++++-------- 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/src/sage/dynamics/arithmetic_dynamics/berkovich_ds.py b/src/sage/dynamics/arithmetic_dynamics/berkovich_ds.py index 722bb7cebbd..cd82dd9e8a4 100644 --- a/src/sage/dynamics/arithmetic_dynamics/berkovich_ds.py +++ b/src/sage/dynamics/arithmetic_dynamics/berkovich_ds.py @@ -859,12 +859,13 @@ def min_res_locus(self, ret_conjugation=False): OUTPUT: - - If ``ret_conjugation`` is ``False``, then a tuple (``start``, ``end``) is returned, + - If ``ret_conjugation`` is ``False``, then a tuple (``start``, ``end``, ``minimum``) is returned, where ``start`` and ``end`` are points of Berkovich space. The minimal resultant - is then achieved on the interval defined by [``start``, ``end``]. + is then achieved on the interval defined by [``start``, ``end``] with value ``minimum``. - - If ``ret_conjugation`` is ``True``, then a tuple ((``start``, ``end``), ``conj``) - is returned, where [``start``, ``end``] defines the minimal resultant locus + - If ``ret_conjugation`` is ``True``, then a tuple (``start``, ``end``, ``minimum``, ``conj``) + is returned, where [``start``, ``end``] defines the minimal resultant locus, + ``minimum`` is the minimum value of order of the resultant, and ``conj`` is a matrix such that this dynamical system achieves the minimal order of the resultant after conjugation by ``conj``. @@ -895,12 +896,13 @@ def min_res_locus(self, ret_conjugation=False): Q = num - value*dem factorization += list(Q.factor()) min_list = [] - print('all factors:', factorization) + #print('all factors:', factorization) for i in range(len(factorization)): irreducible = factorization[i][0] extension = base.extension(irreducible, 's%s'%i) new_primes = extension.primes_above(base_prime_ideal) for new_prime in new_primes: + #print('current prime = ', new_prime) ramification_index = new_prime.absolute_ramification_index() valuation = lambda x: x.valuation(new_prime)/ramification_index s = extension.gens()[0] @@ -912,6 +914,7 @@ def min_res_locus(self, ret_conjugation=False): res = valuation(new_system.resultant()) #print('ordres = ', res) d = new_system.degree() + #print('d = ', d) #print('system = ', new_system) C_list = [] D_list = [] @@ -929,13 +932,16 @@ def min_res_locus(self, ret_conjugation=False): else: D_list.append(None) #print('C_list =', C_list) + #print('D_list =', D_list) D = QQ['t'] t = D.gens()[0] - C_lines = [(C_list[i] + (QQ(d^2+d-2*d*i))*t) for i in range(len(C_list)) if C_list[i] != None] - D_lines = [(D_list[i] + (QQ(d^2+d-2*d*(i+1)))*t) for i in range(len(D_list)) if D_list[i] != None] + C_lines = [(C_list[i] + QQ(d**2 + d - 2*d*i)*t) for i in range(len(C_list)) if C_list[i] != None] + D_lines = [(D_list[i] + QQ(d**2 + d - 2*d*(i + 1))*t) for i in range(len(D_list)) if D_list[i] != None] + #print('C_lines = ', C_lines) + #print('D_lines = ', D_lines) all_lines = C_lines + D_lines #crude minimization minimum = None - #print(C_lines,D_lines) + #print('all lines = ', C_lines, D_lines) all_intersections = [] for line_1 in all_lines: for line_2 in all_lines: @@ -966,7 +972,7 @@ def min_res_locus(self, ret_conjugation=False): minimizing_list = [i[0] for i in min_list] minimum = min(minimizing_list) min_index = minimizing_list.index(minimum) - min_res = min_list[min_index][0] + min_ord_res = min_list[min_index][0] intersection_final = min_list[min_index][1] interval_check = min_list[min_index][2] extension = min_list[min_index][3] @@ -974,9 +980,8 @@ def min_res_locus(self, ret_conjugation=False): if not interval_check: #minimum is not achieved on an interval B = Berkovich_Cp_Projective(extension, prime_above) min_res_locus = B(extension.gens()[0], power=-1*intersection_final) - min_res_locus = min_list[min_index][0] if not ret_conjugation: - return min_res_locus + return (min_res_locus, min_ord_res) else: #print('intersection_final:', intersection_final) root = intersection_final.denominator() @@ -1005,7 +1010,7 @@ def min_res_locus(self, ret_conjugation=False): conj_matrix = Matrix([[A, B],[0,1]]) new_system = system.change_ring(new_extension) new_system = new_system.conjugate(conj_matrix) - return (min_res_locus, conj_matrix, new_system) + return (min_res_locus, min_ord_res, conj_matrix) else: print('not implemented') return [min_list[min_index][0], min_list[min_index][1]] From c21fb94a1633c2ba79eda2d4f1f52404761bcd9b Mon Sep 17 00:00:00 2001 From: Alexander Galarraga Date: Wed, 5 Aug 2020 10:54:21 -0400 Subject: [PATCH 3/3] 30294: beginning work on step 4.b. --- .../arithmetic_dynamics/berkovich_ds.py | 191 +++++++++--------- 1 file changed, 97 insertions(+), 94 deletions(-) diff --git a/src/sage/dynamics/arithmetic_dynamics/berkovich_ds.py b/src/sage/dynamics/arithmetic_dynamics/berkovich_ds.py index b0f3110d5e9..88a8841ec36 100644 --- a/src/sage/dynamics/arithmetic_dynamics/berkovich_ds.py +++ b/src/sage/dynamics/arithmetic_dynamics/berkovich_ds.py @@ -878,6 +878,8 @@ def min_res_locus(self, ret_conjugation=False): r""" Returns the minimal resultant locus and the minimal order of the resultant. + Only implemented for Berkovich spaces backed by number fields. + Let the action of this map `\phi` be given by `[X : Y] \to [F(X,Y) : G(X,Y)]`, where `F` and `G` are normalized. The resultant of `\phi` is the resultant of `F` and `G`. @@ -908,8 +910,87 @@ def min_res_locus(self, ret_conjugation=False): ALGORITHM: - See [Rum2013]_. + Algorithm A of [Rum2013]_ (p.28). """ + + # as step 2 of Rumely's algorithm may be run twice, we define it as a seperate function + def step_2(factorization, base): + min_list = [] + for i in range(len(factorization)): + irreducible = factorization[i][0] + extension = base.extension(irreducible, 's%s'%i) + new_primes = extension.primes_above(base_prime_ideal) + for new_prime in new_primes: + #print('current prime = ', new_prime) + ramification_index = new_prime.absolute_ramification_index() + valuation = lambda x: x.valuation(new_prime)/ramification_index + s = extension.gens()[0] + transformation_matrix = Matrix([[1, s],[0, 1]]) + new_system = system.change_ring(extension) + new_system = new_system.conjugate(transformation_matrix) + print('before normalizing coordinates = ', new_system) + new_system.normalize_coordinates() #very long call + #print('res = ', new_system.resultant()) + res = valuation(new_system.resultant()) + #print('ordres = ', res) + d = new_system.degree() + #print('d = ', d) + #print('system = ', new_system) + C_list = [] + D_list = [] + for i in range(d + 1): + C_i = new_system[0][i, d - i] + D_i = new_system[1][i, d - i] + #print('C_%s = ' %i, C_i) + #print('val = ', valuation(C_i)) + if C_i != 0: + C_list.append(QQ(res - 2*d*valuation(C_i))) + else: + C_list.append(None) + if D_i != 0: + D_list.append(QQ(res - 2*d*valuation(D_i))) + else: + D_list.append(None) + #print('C_list =', C_list) + #print('D_list =', D_list) + D = QQ['t'] + t = D.gens()[0] + C_lines = [(C_list[i] + QQ(d**2 + d - 2*d*i)*t) for i in range(len(C_list)) if C_list[i] != None] + D_lines = [(D_list[i] + QQ(d**2 + d - 2*d*(i + 1))*t) for i in range(len(D_list)) if D_list[i] != None] + #print('C_lines = ', C_lines) + #print('D_lines = ', D_lines) + all_lines = C_lines + D_lines #crude minimization + minimum = None + print('all lines = ', C_lines, D_lines) + all_intersections = [] + for line_1 in all_lines: + for line_2 in all_lines: + if line_1 != line_2 and line_2[1] != line_1[1]: + intersection = QQ((line_1[0]-line_2[0])/(line_2[1]-line_1[1])) + if intersection not in all_intersections: + all_intersections.append(intersection) + #print(all_intersections) + for intersection in all_intersections: + C_max = max([(i(intersection)) for i in C_lines]) + D_max = max([(i(intersection)) for i in D_lines]) + print('intersection = ', intersection) + #print('C_max = ', C_max) + #print('D_max = ', D_max) + print('value of X(i) = ', max(C_max, D_max)) + X_i = (max(C_max,D_max)) + if X_i < minimum or minimum == None: + minimum = X_i + intersection_final = QQ(intersection) + print('minimum: ', minimum) + interval_check = False + for line in all_lines: #check if minimum is achieved on an interval + if line.is_constant(): + if line[0] == minimum: + interval_check = True + break + min_list.append((minimum, intersection_final, interval_check, extension, new_prime)) + return min_list + system = self._system QQ_prime = self.domain().prime() base_prime_ideal = self.domain().ideal() @@ -924,7 +1005,7 @@ def min_res_locus(self, ret_conjugation=False): num = phi(num) #print('num:', num) #print('dem:', dem) - fixed_point_equation = (num-w*dem) + fixed_point_equation = num - w*dem factorization = list(fixed_point_equation.factor()) if(system([1,0]) == system.domain()([1,0])): Q = dem @@ -932,80 +1013,8 @@ def min_res_locus(self, ret_conjugation=False): value = system([1,0]).dehomogenize(1)[0] Q = num - value*dem factorization += list(Q.factor()) - min_list = [] #print('all factors:', factorization) - for i in range(len(factorization)): - irreducible = factorization[i][0] - extension = base.extension(irreducible, 's%s'%i) - new_primes = extension.primes_above(base_prime_ideal) - for new_prime in new_primes: - #print('current prime = ', new_prime) - ramification_index = new_prime.absolute_ramification_index() - valuation = lambda x: x.valuation(new_prime)/ramification_index - s = extension.gens()[0] - transformation_matrix = Matrix([[1, s],[0, 1]]) - new_system = system.change_ring(extension) - new_system = new_system.conjugate(transformation_matrix) - new_system.normalize_coordinates() #very long call - #print('res = ', new_system.resultant()) - res = valuation(new_system.resultant()) - #print('ordres = ', res) - d = new_system.degree() - #print('d = ', d) - #print('system = ', new_system) - C_list = [] - D_list = [] - for i in range(d+1): - C_i = new_system[0][i,d-i] - D_i = new_system[1][i,d-i] - #print('C_%s = ' %i, C_i) - #print('val = ', valuation(C_i)) - if C_i != 0: - C_list.append(QQ(res-2*d*valuation(C_i))) - else: - C_list.append(None) - if D_i != 0: - D_list.append(QQ(res-2*d*valuation(D_i))) - else: - D_list.append(None) - #print('C_list =', C_list) - #print('D_list =', D_list) - D = QQ['t'] - t = D.gens()[0] - C_lines = [(C_list[i] + QQ(d**2 + d - 2*d*i)*t) for i in range(len(C_list)) if C_list[i] != None] - D_lines = [(D_list[i] + QQ(d**2 + d - 2*d*(i + 1))*t) for i in range(len(D_list)) if D_list[i] != None] - #print('C_lines = ', C_lines) - #print('D_lines = ', D_lines) - all_lines = C_lines + D_lines #crude minimization - minimum = None - #print('all lines = ', C_lines, D_lines) - all_intersections = [] - for line_1 in all_lines: - for line_2 in all_lines: - if line_1 != line_2 and line_2[1] != line_1[1]: - intersection = QQ((line_1[0]-line_2[0])/(line_2[1]-line_1[1])) - if intersection not in all_intersections: - all_intersections.append(intersection) - #print(all_intersections) - for intersection in all_intersections: - C_max = max([(i(intersection)) for i in C_lines]) - D_max = max([(i(intersection)) for i in D_lines]) - #print('intersection = ', intersection) - #print('C_max = ', C_max) - #print('D_max = ', D_max) - #print('value of X(i) = ', max(C_max, D_max)) - X_i = (max(C_max,D_max)) - if X_i < minimum or minimum == None: - minimum = (X_i) - intersection_final = QQ(intersection) - #print('minimum: ', minimum) - interval_check = False - for line in all_lines: #check if minimum is achieved on an interval - if line.is_constant(): - if line[0] == minimum: - interval_check = True - break - min_list.append((minimum, intersection_final, interval_check, extension, new_prime)) + min_list = step_2(factorization, base) minimizing_list = [i[0] for i in min_list] minimum = min(minimizing_list) min_index = minimizing_list.index(minimum) @@ -1025,32 +1034,26 @@ def min_res_locus(self, ret_conjugation=False): if root == 1: #print('in no denominator case') #print('extension:', extension) - conj_matrix = Matrix([[QQ_prime**(intersection_final), extension.gens()[0]], [0,1]]) + conj_matrix = Matrix([[QQ_prime**(intersection_final), extension.gens()[0]], [0, 1]]) new_system = system.change_ring(extension) else: C = extension['c'] c = C.gens()[0] - defining_polynomial = c**(root) - QQ_prime**(intersection_final.numerator()) - factorization = defining_polynomial.factor() - degree = None - for factor in factorization: - if degree == None or factor[0].degree() < degree: #slow symbolic calcuation - final_factor = factor[0] - degree = factor[0].degree() - if degree != 1: - new_extension = extension.extension(final_factor, 'w0') - A = new_extension.gens()[0] - B = new_extension.gens()[1] - else: - A = final_factor[0] - B = extension.gens()[0] - conj_matrix = Matrix([[A, B],[0,1]]) + defining_polynomial = c**root - QQ_prime**intersection_final.numerator() + #print('defining polynomial = ', defining_polynomial) + new_extension, ext_to_new_ext = defining_polynomial.splitting_field('w0', map=True) + A = defining_polynomial.roots(new_extension)[0][0] + B = ext_to_new_ext(extension.gens()[0]) + conj_matrix = Matrix([[A, B], [0, 1]]) new_system = system.change_ring(new_extension) new_system = new_system.conjugate(conj_matrix) - return (min_res_locus, min_ord_res, conj_matrix) + return (min_res_locus, min_ord_res, conj_matrix, new_extension, extension, ext_to_new_ext) else: - print('not implemented') - return [min_list[min_index][0], min_list[min_index][1]] + list_of_minimums = [tup for tup in min_list if tup[0] == min_ord_res] + max_radius = max([QQ_prime**(-1*tup[1]) for tup in list_of_minimums]) + maximum_radii = [tup for tup in list_of_minimums if QQ_prime**(-1*tup[1]) == max_radius] + return maximum_radii + class DynamicalSystem_Berkovich_affine(DynamicalSystem_Berkovich): r"""