Merge branch 'u/gh-EnderWannabe/min_res_locus' of git://trac.sagemath…

….org/sage into min_res_locus

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

src/sage/dynamics/arithmetic_dynamics/berkovich_ds.py

@@ -43,6 +43,7 @@
 from sage.rings.integer_ring import ZZ
 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):
@@ -742,7 +743,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()
@@ -874,6 +874,187 @@ 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.
+
+        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`.
+        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``, ``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``] with value ``minimum``.
+
+        - 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``.
+
+        ALGORITHM:
+
+        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()
+        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())
+        #print('all factors:', factorization)
+        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)
+        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]
+        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)
+            if not ret_conjugation:
+                return (min_res_locus, min_ord_res)
+            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()
+                    #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, new_extension, extension, ext_to_new_ext)
+        else:
+            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"""
     A dynamical system of the affine Berkovich line over `\CC_p`.