From 0544a244e4199c6d89e1d76ceae593c71bc57c71 Mon Sep 17 00:00:00 2001 From: rlmiller Date: Tue, 16 Aug 2016 13:58:22 -0500 Subject: [PATCH] 21248 Adding reduced binary form algorithm --- .../rings/polynomial/multi_polynomial.pyx | 174 +++++++++++++++++- .../schemes/projective/projective_morphism.py | 82 +++++++++ 2 files changed, 253 insertions(+), 3 deletions(-) diff --git a/src/sage/rings/polynomial/multi_polynomial.pyx b/src/sage/rings/polynomial/multi_polynomial.pyx index a53d4a93166..0ebeb7be837 100644 --- a/src/sage/rings/polynomial/multi_polynomial.pyx +++ b/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,6 +1826,169 @@ 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. = 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. = 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) @@ -1827,5 +1996,4 @@ cdef remove_from_tuple(e, int ind): if len(w) == 1: return w[0] else: - return tuple(w) - + return tuple(w) \ No newline at end of file diff --git a/src/sage/schemes/projective/projective_morphism.py b/src/sage/schemes/projective/projective_morphism.py index 0df098923a5..8b768925507 100644 --- a/src/sage/schemes/projective/projective_morphism.py +++ b/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. = 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. = 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. = 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):