21248 Reduced Binary form, added examples and documentation

src/sage/rings/polynomial/multi_polynomial.pyx

@@ -31,6 +31,7 @@ from sage.modules.free_module_element import vector
 from sage.rings.rational_field import QQ
 from sage.arith.misc import gcd
 from sage.rings.complex_interval_field import ComplexIntervalField
+from sage.rings.real_mpfr import RealField_class,RealField
 
 cdef class MPolynomial(CommutativeRingElement):
 
@@ -1874,6 +1875,16 @@ cdef class MPolynomial(CommutativeRingElement):
             Traceback (most recent call last):
             ...
             ValueError: (=x^3*y*z + x^4 + y^2*z) must have two variables
+
+        ::
+
+            sage: R.<x,y> = PolynomialRing(ZZ)
+            sage: F = - 8*x^6 - 3933*x^3*y - 725085*x^2*y^2 - 59411592*x*y^3 - 99*y^6
+            sage: F.reduced_form(return_conjugation=False)
+            Traceback (most recent call last):
+            ...
+            ValueError: (=-8*x^6 - 99*y^6 - 3933*x^3*y - 725085*x^2*y^2 -
+            59411592*x*y^3) must be homogenous
         """
         from sage.matrix.constructor import matrix
         from sage.calculus.functions import jacobian
@@ -1882,7 +1893,8 @@ cdef class MPolynomial(CommutativeRingElement):
             raise ValueError("(=%s) must have two variables"%self)
         if self.is_homogeneous() != True:
             raise ValueError("(=%s) must be homogenous"%self)
-        KK = ComplexIntervalField(prec=prec)
+        KK = ComplexIntervalField(prec=prec) # keeps trac of our precision error
+        JJ=RealField()
         R = self.parent()
         S = PolynomialRing(R.base_ring(),'z')
         phi = R.hom([S.gen(0),1],S)# dehomogenization
@@ -1893,6 +1905,7 @@ cdef class MPolynomial(CommutativeRingElement):
         R=PolynomialRing(KK,'x,y')
         x,y = R.gens()
         Q  =[]
+        # finds Stoll and Cremona's Q_0'
         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)
@@ -1901,13 +1914,15 @@ cdef class MPolynomial(CommutativeRingElement):
         A = q.monomial_coefficient(x**2)
         B = q.monomial_coefficient(x*y)
         C = q.monomial_coefficient(y**2)
+        # need positive root
         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
+        M = matrix(QQ, [[1,0], [0,1]]) # used to keep track of how our z is moved.
+        # this moves Z to our fundamental domain using the three steps laid out in the algorithim
         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:
@@ -1923,7 +1938,7 @@ cdef class MPolynomial(CommutativeRingElement):
             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]))))
+        q = q(tuple((M * vector([x, y])))) #our new Q
         A = q.monomial_coefficient(x**2)
         B = q.monomial_coefficient(x*y)
         C = q.monomial_coefficient(y**2)
@@ -1944,7 +1959,7 @@ cdef class MPolynomial(CommutativeRingElement):
             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
+        #Newton's Method, error bound is while less than diameter of our z
         err = z0.diameter()
         zz = z0.diameter()
         n = q.degree()
@@ -1954,6 +1969,7 @@ cdef class MPolynomial(CommutativeRingElement):
         J = jacobian(G, [u,t])
         v0 = vector([z0.imag(), z0.real()]) #z0 as starting point
         z = z0
+        #finds our correct z
         while err <= zz:
             NJ = J.subs({u:v0[0], t:v0[1]})
             NJinv = NJ.inverse()
@@ -1967,8 +1983,9 @@ cdef class MPolynomial(CommutativeRingElement):
             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
+            err = z.diameter() # precision
+            zz = (w - z).abs() #difference in w and z
+        # moves our z ro fundamental domain as before
         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:
@@ -1987,7 +2004,7 @@ cdef class MPolynomial(CommutativeRingElement):
         err = z.diameter()
         zc = z.center()
         zr = zc.real()
-        if (zr.abs() - KK(0.5)).abs() <= err or (zc.abs() - 1).abs() <= err:
+        if (zr.abs() - JJ(0.5)).abs() <= err or (zc.abs() - 1).abs() <= err:
             raise ValueError("not enough precision")
         x,y = self.parent().gens()
         if return_conjugation:

src/sage/schemes/projective/projective_morphism.py

@@ -3525,6 +3525,8 @@ def reduced_form(self, prec=100, return_conjugation=True):  # add bool for matri
         the binary form reduction algorithm from Stoll and Cremona _[SC]
         to the periodic points (the dynatomic polynomial).
 
+        Implimented by Rebecca Lauren Miller as part of GSOC 2016.
+
         INPUT::
 
         - ``prec`` -- desired precision (default: 100)
@@ -3559,16 +3561,15 @@ def reduced_form(self, prec=100, return_conjugation=True):  # add bool for matri
             )
 
         ::
-
-            sage: PS.<x,y> = ProjectiveSpace(QQ,1)
+            # failing why?
+            sage: PS.<x,y> = ProjectiveSpace(ZZ,1)
             sage: H = End(PS)
             sage: f = H([x^2+ x*y,y^2]) #needs 3 periodic
             sage: m = matrix(QQ, 2, 2,[-221, -1, 1, 0])
             sage: f = f.conjugate(m)
             sage: f.reduced_form(prec=200)
             (
-            Scheme endomorphism of Projective Space of dimension 1 over Rational
-            Field
+            Scheme endomorphism of Projective Space of dimension 1 over Integer Ring
             Defn: Defined on coordinates by sending (x : y) to
             (-x^2 + x*y - y^2 : -y^2)
             ,
@@ -3592,22 +3593,34 @@ def reduced_form(self, prec=100, return_conjugation=True):  # add bool for matri
             [-1  2]
             [ 2 -5]
             )
+
+         ::
+            # failing why?
+            sage: PS.<x,y> = ProjectiveSpace(ZZ,1)
+            sage: H = End(PS)
+            sage: f = H([7*x^3 + 5*x*y^2 + y, 6*y^3])
+            sage: m = matrix(QQ, 2, 2,[1 ,35, 0, 1])*matrix(QQ, 2, 2,[0, 1, -1, 0])*matrix(QQ,2 ,2, [1,0,131,1])
+            sage: f = f.conjugate(m)
+            sage: f.reduced_form(prec=200)
+            Traceback (most recent call last):
+            ...
+            ValueError: polys (=[7*x^3 + 5*x*y^2 + y, 6*y^3]) must be homogeneous
         """
         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
+        F = R(F/gcd(F,F.derivative(x))) #removes 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
+        # Checks to make sure there are enough distinct, roots we need 3
+        # if there are not if finds the nth period fix points until there are enough
+        while F.degree() <= 2:
+            F = self.dynatomic_polynomial(n) # finds n period fix points
+            F = R(F/gcd(F,F.derivative(x)))#removes 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):