21248 Adding reduced binary form algorithm

src/sage/rings/polynomial/multi_polynomial.pyx

@@ -15,9 +15,15 @@ def is_MPolynomial(x):
     return isinstance(x, MPolynomial)
 
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
-
 from sage.categories.map cimport Map
 from sage.categories.morphism cimport Morphism
+#from sage.calculus.functions import jacobian
+#from sage.modules.free_module_element import vector
+#from sage.rings.rational_field import QQ
+# from sage.matrix.constructor import matrix
+#from sage.arith.misc import gcd
+#from sage.rings.complex_interval_field import ComplexIntervalField
+
 
 cdef class MPolynomial(CommutativeRingElement):
 
@@ -1820,12 +1826,174 @@ cdef class MPolynomial(CommutativeRingElement):
 
             return ans
 
+    def reduced_form(self, prec=100, return_conjugation=True):
+        r"""
+        Returns reduced binary form of a polynomial.
+
+        Implimented algorithim from Stoll and Cremona's "On the Reduction Theory of Binary Forms" _[SC].
+        This takes a degree two homogenous polynomial and finds it's roots. Then it finds a unique Q_0 which is
+        in quadratic form, we apply the quadratic folrmula to this in order to find our z. This s is them moved into
+        the fundamental domain. We keep track f the elements of PGL(2) used to move our z. This z gives us a new Q and through
+        the same process a new z. We then apply this element of PGL(2) to our original polynomial, this puts it in reduced form.
+
+        Implimented by Rebecca Lauren Miller as part of GSOC 2016.
+
+        INPUT:
+
+        -``prec`` --  sets the precision (default:100)
+
+        - ``conjugation`` -- A booolean. Returns element of PGL(2) (default:True)
+
+        OUTPUT:
+            - reduced binary form
+
+            - element of PFL(2)
+
+        REFERENCES:
+
+        .. [SC] Michael Stoll and John E. Cremona. On The Reduction Theory of Binary Forms. Journal für die reine und angewandte Mathematik, 565 (2003), 79-99.
+
+        EXAMPLES::
+
+            sage: R.<x,y> = PolynomialRing(QQ)
+            sage: F = - 8*x^4 - 3933*x^3*y - 725085*x^2*y^2 - 59411592*x*y^3 - 1825511633*y^4
+            sage: F.reduced_form(prec=120, return_conjugation=False)
+            x^4 + 9*x^3*y - 3*x*y^3 - 8*y^4
+
+        ::
+
+#sage: R.<x,y,z> = PolynomialRing(QQ)
+# sage: F = x^4 + x^3*y*z + y^2*z
+#sage: F.reduced_form(prec=200, return_conjugation=False)
+#ValueError: need two variable polynomial
+#Failing not sure why?
+        """
+        # need to have error example  with  non-homogenous
+        from sage.calculus.functions import jacobian
+        from sage.modules.free_module_element import vector
+        from sage.rings.rational_field import QQ
+        from sage.matrix.constructor import matrix
+        from sage.arith.misc import gcd
+        from sage.rings.complex_interval_field import ComplexIntervalField
+
+        if self.parent().ngens() != 2:
+            raise ValueError("need two variable polynomial")
+        if self.is_homogeneous() != True:
+            raise ValueError("polynomial must be homogenous")
+        KK = ComplexIntervalField(prec=prec)
+        R = self.parent()
+        S = PolynomialRing(R.base_ring(),'z')
+        phi = R.hom([S.gen(0),1],S)# dehomogenization
+        FF = S(phi(self)/gcd(phi(self),phi(self).derivative())) # removes multiple roots
+        roots = FF.roots(ring=KK, multiplicities=False)
+        dF = FF.derivative()
+        n = FF.degree()
+        R=PolynomialRing(KK,'x,y')
+        x,y = R.gens()
+        Q  =[]
+        for j in range(len(roots)):
+            k = (1/(dF(roots[j]).abs()**(2/(n-2))))*((x-(roots[j]*y))*(x-(roots[j].conjugate()*y)))
+            Q.append(k)
+        q = sum([Q[i] for i in range(len(Q))]) # this is q_o , always positive def as long as F HAS DISTINCT ROOTS
+
+        A = q.monomial_coefficient(x**2)
+        B = q.monomial_coefficient(x*y)
+        C = q.monomial_coefficient(y**2)
+        try:
+            z=(-B + ((B**2)-(4*A*C)).sqrt())/(2*A)# this is z_o
+        except ValueError:
+            raise  ValueError("not enough precision")
+        if z.imag()<0:
+            z=(-B - ((B**2)-(4*A*C)).sqrt())/(2*A)
+        M=matrix(QQ,[[1,0],[0,1]])# this moves Z to fundamental domain
+        while z.real()<-0.5 or z.real()>=0.5 or (z.real()<=0 and z.abs()<1) or (z.real()>0 and z.abs()<=1):
+            if z.real()<-0.5:
+                m = z.real().abs().round() # member abs
+                Qm = QQ(m.center())
+                M = M*matrix(QQ,[[1,-Qm],[0,1]]) # move
+                z += m  # M.inverse()*Z is supposed to move z by m
+            elif z.real()>=(0.5): #else if
+                m = z.real().round()
+                Qm = QQ(m.center())
+                M = M*matrix(QQ,[[1,Qm],[0,1]])  #move z
+                z -= m
+            elif (z.real()<=0 and z.abs()<1) or (z.real()>0 and z.abs()<=1):
+                z = -1/z
+                M = M*matrix(QQ,[[0,-1],[1,0]])
+        q = q(tuple((M*vector([x,y]))))
+        A = q.monomial_coefficient(x**2)
+        B = q.monomial_coefficient(x*y)
+        C = q.monomial_coefficient(y**2)
+        try:
+            z0=(-B + ((B**2)-(4*A*C)).sqrt())/(2*A)# this is z_o
+        except ValueError:
+            raise ValueError("not enough precision")
+        if z.imag()<0:
+            z0=(-B - ((B**2)-(4*A*C)).sqrt())/(2*A)
+        L = [ (-B + ((B**2)-(4*A*C)).sqrt())/(2*A), (-B - ((B**2)-(4*A*C)).sqrt())/(2*A)] # if highr than quadratic need to change?
+        a = 0
+        c = 0
+        RR = PolynomialRing(KK,'u,t')
+        u,t = RR.gens()
+        for j in range(len(L)):
+            b= u**2/((t-(L[j]))*(t-(L[j].conjugate()))+ u**2)
+            a += b
+            d = (t-(L[j].real()))/((t-(L[j]))*(t-(L[j].conjugate()))+ u**2)
+            c += d
+        # return "A"broke here
+        #Newton's Method, error bound is while less than diameter of our z instaed of solve use newtons method
+        err = z0.diameter()
+        zz =  z0.diameter()
+        n = q.degree()
+        g1 = a.numerator() - n/2*a.denominator()
+        g2 = c.numerator()
+        G = vector([g1,g2])
+        J = jacobian(G,[u,t])
+        v0 = vector([z0.imag(),z0.real()]) #z0 as starting point
+        z = z0
+        while err <= zz:
+            NJ = J.subs({u:v0[0],t:v0[1]})
+            NJinv = NJ.inverse()
+            for ii in range(2):
+                for jj in range(2):
+                    tt = NJinv[ii,jj]
+                    try:
+                        NJinv[ii,jj] = (tt.numerator()/tt.denominator()).numerator()
+                    except:
+                        pass
+            w = z
+            v0 = v0 - NJinv*G.subs({u:v0[0],t:v0[1]})
+            z = v0[1].coefficients()[0] + v0[0].coefficients()[0]*KK.gen(0)
+            err = z.diameter() # prec
+            zz = (w - z).abs() #diff in w and z
+        while z.real()<-0.5 or z.real()>=0.5or (z.real()<=0 and z.abs()<1) or (z.real()>0 and z.abs()<=1):
+            if z.real()<-0.5:
+                m = z.real().abs().round() # member abs
+                Qm = QQ(m.center())
+                M = M*matrix(QQ,[[1,-Qm],[0,1]]) # move
+                z += m  # M.inverse()*Z is supposed to move z by m
+            elif z.real()>=(0.5): #else if
+                m = z.real().round()
+                Qm = QQ(m.center())
+                M = M*matrix(QQ,[[1,Qm],[0,1]])  #move z
+                z -= m
+            elif (z.real()<=0 and z.abs()<1) or (z.real()>0 and z.abs()<=1):
+                z = -1/z
+                M = M*matrix(QQ,[[0,-1],[1,0]])
+        err = z.diameter()
+        zc = z.center()
+        zr = zc.real()
+        if (zr.abs()-1/2).abs() <= err or (zc.abs()-1).abs() <= err:
+            raise ValueError("not enough precision")
+        x,y = self.parent().gens()
+        if return_conjugation:
+            return (self(tuple(M*vector([x,y]))), M)
+        return self(tuple(M*vector([x,y])))
 
 cdef remove_from_tuple(e, int ind):
     w = list(e)
     del w[ind]
     if len(w) == 1:
         return w[0]
     else:
-        return tuple(w)
-
+        return tuple(w)

src/sage/schemes/projective/projective_morphism.py

@@ -3414,6 +3414,88 @@ def sigma_invariants(self, n, formal=True, embedding=None):
             polys.append(e([i+1]).expand(N)(multipliers))
         return polys
 
+    def reduced_form(self, prec=100, return_conjugation=True):  # add bool for matrices, set prec? author ref?
+        r"""
+        Returns reduced form of a Projective Morphism.
+
+        Implimeting Reduced fuction from Stoll and Cremona's paper "On the Reduction theory of Binary Forms" _[SC].
+
+        Takes Projective morphism and calles reduced_form fuction on the Dynatomic polynomial in order to return the reduced formof the mapping.
+
+        INPUT::
+
+        - ``prec`` -- desired precision (default: 100)
+
+        -``return_conjuagtion`` -- A Boolean. Returns element of PGL(2). (default: True)
+
+        OUTPUT:
+
+        Reduced mapping
+
+        Conjuagtion element of PGL(2)
+
+        EXAMPLES:
+
+            #sage: PS.<x,y> = ProjectiveSpace(QQ,1)
+            #sage: H = End(PS)
+            #sage: f=H([x^3+x*y^2,y^3])
+            #sage: m=matrix(QQ,2,2,[1,13,0,1])*matrix(QQ,2,2,[0,-1,1,0])*matrix(QQ,2,2,[1,0,234,1])
+            #sage: f=f.conjugate(m)
+            #sage: f.reduced_form(prec=200) #needs 2 periodic
+            #ValueError: not enough precision
+
+        ::
+
+            sage: PS.<x,y> = ProjectiveSpace(QQ,1)
+            sage: H = End(PS)
+            sage: f=H([x^2+ x*y,y^2]) #needs 3 periodic
+            sage: m=matrix(QQ,2,2,[1,13,0,1])*matrix(QQ,2,2,[0,-1,1,0])*matrix(QQ,2,2,[1,0,234,1])
+            sage: f=f.conjugate(m)
+            sage: f.reduced_form(prec=200)
+            (
+            Scheme endomorphism of Projective Space of dimension 1 over Rational
+            Field
+            Defn: Defined on coordinates by sending (x : y) to
+            (-x^2 + x*y - y^2 : -y^2)
+            ,
+            [  0  -1]
+            [  1 220]
+            )
+
+
+        ::
+
+            sage: PS.<X,Y> = ProjectiveSpace(QQ,1)
+            sage: H = End(PS)
+            sage: f = H([7365/2*X^4 + 6282*X^3*Y + 4023*X^2*Y^2 + 1146*X*Y^3 + 245/2*Y^4,-12329/2*X^4 - 10506*X^3*Y - 6723*X^2*Y^2 - 1914*X*Y^3 - 409/2*Y^4])
+            sage: f.reduced_form()
+            (
+            Scheme endomorphism of Projective Space of dimension 1 over Rational
+            Field
+            Defn: Defined on coordinates by sending (X : Y) to
+            (-7/2*X^4 - 6*X^3*Y - 21*X^2*Y^2 - 6*X*Y^3 - 7/2*Y^4 : -1/2*X^4
+            - 2*X^3*Y - 3*X^2*Y^2 - 2*X*Y^3 - 1/2*Y^4),
+            [-1  2]
+            [ 2 -5]
+            )
+
+
+        """
+        R = self.coordinate_ring()
+        F = R(self.dynatomic_polynomial(1))
+        x,y = R.gens()
+        #F.parent() == R #needs ticket 21097,
+        F = R(F/gcd(F,F.derivative(x))) #remove multiple roots
+        n=2
+        while F.degree() <= 2: # checks to have enough distinct roots need at least 3
+            F = self.dynatomic_polynomial(n) # n period fix points
+            F = R(F/gcd(F,F.derivative(x)))#remove multiple roots
+            n+=1
+        G,m = F.reduced_form( prec=prec, return_conjugation=return_conjugation)
+        if return_conjugation:
+            return (self.conjugate(m), m)
+        return self.conjugate(m)
+            # answers to examples have fractions?, 1st example now doenst have enough prec!, problem raising value error
 class SchemeMorphism_polynomial_projective_space_field(SchemeMorphism_polynomial_projective_space):
 
     def lift_to_rational_periodic(self, points_modp, B=None):