hyperelliptic_generic.py

"""
Hyperelliptic curves over a general ring

EXAMPLE::

    sage: P.<x> = GF(5)[]
    sage: f = x^5 - 3*x^4 - 2*x^3 + 6*x^2 + 3*x - 1
    sage: C = HyperellipticCurve(f); C
    Hyperelliptic Curve over Finite Field of size 5 defined by y^2 = x^5 + 2*x^4 + 3*x^3 + x^2 + 3*x + 4
    
EXAMPLE::

    sage: P.<x> = QQ[]
    sage: f = 4*x^5 - 30*x^3 + 45*x - 22
    sage: C = HyperellipticCurve(f); C
    Hyperelliptic Curve over Rational Field defined by y^2 = 4*x^5 - 30*x^3 + 45*x - 22
    sage: C.genus()
    2

    sage: D = C.affine_patch(0)
    sage: D.defining_polynomials()[0].parent()
    Multivariate Polynomial Ring in x0, x1 over Rational Field
"""

#*****************************************************************************
#  Copyright (C) 2006 David Kohel <[email protected]>
#  Distributed under the terms of the GNU General Public License (GPL)
#                  http://www.gnu.org/licenses/
#*****************************************************************************

from sage.rings.all import PolynomialRing, RR, PowerSeriesRing, LaurentSeriesRing, O
from sage.functions.all import log

import sage.schemes.plane_curves.projective_curve as plane_curve

def is_HyperellipticCurve(C):
    """
    EXAMPLES::
    
        sage: R.<x> = QQ[]; C = HyperellipticCurve(x^3 + x - 1); C
        Hyperelliptic Curve over Rational Field defined by y^2 = x^3 + x - 1
        sage: sage.schemes.hyperelliptic_curves.hyperelliptic_generic.is_HyperellipticCurve(C)
        True
    """
    return isinstance(C,HyperellipticCurve_generic)

class HyperellipticCurve_generic(plane_curve.ProjectiveCurve_generic):
    def __init__(self, PP, f, h=None, names=None, genus=None):
        x, y, z = PP.gens()
        df = f.degree()
        F1 = sum([ f[i]*x**i*z**(df-i) for i in range(df+1) ])
        if h is None:
            F = y**2*z**(df-2) - F1
        else:
            dh = h.degree()
            deg = max(df,dh+1)
            F0 = sum([ h[i]*x**i*z**(dh-i) for i in range(dh+1) ])
            F = y**2*z**(deg-2) + F0*y*z**(deg-dh-1) - F1*z**(deg-df)
        plane_curve.ProjectiveCurve_generic.__init__(self,PP,F)
        R = PP.base_ring()
        if names == None:
            names = ["x","y"]
        elif isinstance(names,str):
            names = names.split(",")
        self._names = names
        P1 = PolynomialRing(R,name=names[0])
        P2 = PolynomialRing(P1,name=names[1])
        self._PP = PP
        self._printing_ring = P2
        self._hyperelliptic_polynomials = (f,h)
        self._genus = genus
        
    def change_ring(self, R):
        """
        Returns this HyperEllipticCurve over a new base ring R.
        
        EXAMPLES::

            sage: R.<x> = QQ['x']
            sage: H = HyperellipticCurve(x^3-10*x+9)
            sage: K = Qp(3,5)
            sage: J.<a> = K.extension(x^30-3)
            sage: HK = H.change_ring(K)
            sage: HJ = HK.change_ring(J); HJ
            Hyperelliptic Curve over Eisenstein Extension of 3-adic Field with capped relative precision 5 in a defined by (1 + O(3^5))*x^30 + (O(3^6))*x^29 + (O(3^6))*x^28 + (O(3^6))*x^27 + (O(3^6))*x^26 + (O(3^6))*x^25 + (O(3^6))*x^24 + (O(3^6))*x^23 + (O(3^6))*x^22 + (O(3^6))*x^21 + (O(3^6))*x^20 + (O(3^6))*x^19 + (O(3^6))*x^18 + (O(3^6))*x^17 + (O(3^6))*x^16 + (O(3^6))*x^15 + (O(3^6))*x^14 + (O(3^6))*x^13 + (O(3^6))*x^12 + (O(3^6))*x^11 + (O(3^6))*x^10 + (O(3^6))*x^9 + (O(3^6))*x^8 + (O(3^6))*x^7 + (O(3^6))*x^6 + (O(3^6))*x^5 + (O(3^6))*x^4 + (O(3^6))*x^3 + (O(3^6))*x^2 + (O(3^6))*x + (2*3 + 2*3^2 + 2*3^3 + 2*3^4 + 2*3^5 + O(3^6)) defined by (1 + O(a^150))*y^2 = (1 + O(a^150))*x^3 + (2 + 2*a^30 + a^60 + 2*a^90 + 2*a^120 + O(a^150))*x + a^60 + O(a^210)
        """
        from constructor import HyperellipticCurve
        f, h = self._hyperelliptic_polynomials
        y = self._printing_ring.variable_name()
        x = self._printing_ring.base_ring().variable_name()
        return HyperellipticCurve(f.change_ring(R), h, "%s,%s"%(x,y))

    def _repr_(self):
        """
        String representation of hyperelliptic curves.
        
        EXAMPLE::
        
            sage: P.<x> = QQ[]
            sage: f = 4*x^5 - 30*x^3 + 45*x - 22
            sage: C = HyperellipticCurve(f); C
            Hyperelliptic Curve over Rational Field defined by y^2 = 4*x^5 - 30*x^3 + 45*x - 22
            sage: C = HyperellipticCurve(f,names='u,v'); C
            Hyperelliptic Curve over Rational Field defined by v^2 = 4*u^5 - 30*u^3 + 45*u - 22
        """

        f, h = self._hyperelliptic_polynomials
        R = self.base_ring()
        y = self._printing_ring.gen()
        x = self._printing_ring.base_ring().gen()
        if h == 0: 
            return "Hyperelliptic Curve over %s defined by %s = %s"%(R, y**2, f(x))
        else:
            return "Hyperelliptic Curve over %s defined by %s + %s = %s"%(R, y**2, h(x)*y, f(x))

    def __cmp__(self, other):
        if not isinstance(other, HyperellipticCurve_generic):
            return -1
        return cmp(self._hyperelliptic_polynomials, other._hyperelliptic_polynomials)

    def hyperelliptic_polynomials(self, K=None, var='x'):
        """
        EXAMPLES::
        
            sage: R.<x> = QQ[]; C = HyperellipticCurve(x^3 + x - 1, x^3/5); C
            Hyperelliptic Curve over Rational Field defined by y^2 + 1/5*x^3*y = x^3 + x - 1
            sage: C.hyperelliptic_polynomials()
            (x^3 + x - 1, 1/5*x^3)
        """
        if K is None:
            return self._hyperelliptic_polynomials
        else:
            f, h = self._hyperelliptic_polynomials
            P = PolynomialRing(K, var)
            return (P(f),P(h))

    def is_singular(self):
        r"""
        Returns False, because hyperelliptic curves are smooth projective
        curves, as checked on construction.

        EXAMPLES::

            sage: R.<x> = QQ[]
            sage: H = HyperellipticCurve(x^5+1)
            sage: H.is_singular()
            False

        A hyperelliptic curve with genus at least 2 always has a singularity at
        infinity when viewed as a *plane* projective curve. This can be seen in
        the following example.::

            sage: R.<x> = QQ[]
            sage: H = HyperellipticCurve(x^5+2)
            sage: set_verbose(None)
            sage: H.is_singular()
            False
            sage: from sage.schemes.plane_curves.projective_curve import ProjectiveCurve_generic
            sage: ProjectiveCurve_generic.is_singular(H)
            True
        """
        return False

    def is_smooth(self):
        r"""
        Returns True, because hyperelliptic curves are smooth projective
        curves, as checked on construction.

        EXAMPLES::

            sage: R.<x> = GF(13)[]
            sage: H = HyperellipticCurve(x^8+1)
            sage: H.is_smooth()
            True

        A hyperelliptic curve with genus at least 2 always has a singularity at
        infinity when viewed as a *plane* projective curve. This can be seen in
        the following example.::

            sage: R.<x> = GF(27, 'a')[]
            sage: H = HyperellipticCurve(x^10+2)
            sage: set_verbose(None)
            sage: H.is_smooth()
            True
            sage: from sage.schemes.plane_curves.projective_curve import ProjectiveCurve_generic
            sage: ProjectiveCurve_generic.is_smooth(H)
            False
        """
        return True

    def lift_x(self, x, all=False):
        f, h = self._hyperelliptic_polynomials
        x += self.base_ring()(0)
        one = x.parent()(1)
        if h.is_zero():
            y2 = f(x)
            if y2.is_square():
                if all:
                    return [self.point([x, y, one], check=False) for y in y2.sqrt(all=True)]
                else:
                    return self.point([x, y2.sqrt(), one], check=False)
        else:
            b = h(x)
            D = b*b + 4*f(x)
            if D.is_square():
                if all:
                    return [self.point([x, (-b+d)/2, one], check=False) for d in D.sqrt(all=True)]
                else:
                    return self.point([x, (-b+D.sqrt())/2, one], check=False)
        if all:
            return []
        else:
            raise ValueError, "No point with x-coordinate %s on %s"%(x, self)


    def genus(self):
        return self._genus

    def jacobian(self):
        import jacobian_generic
        return jacobian_generic.HyperellipticJacobian_generic(self)

    def odd_degree_model(self):
        r"""
        Return an odd degree model of self, or raise ValueError if one does not exist over the field of definition.

        EXAMPLES::

            sage: x = QQ['x'].gen()
            sage: H = HyperellipticCurve((x^2 + 2)*(x^2 + 3)*(x^2 + 5)); H
            Hyperelliptic Curve over Rational Field defined by y^2 = x^6 + 10*x^4 + 31*x^2 + 30
            sage: H.odd_degree_model()
            Traceback (most recent call last):
            ...
            ValueError: No odd degree model exists over field of definition

            sage: K2 = QuadraticField(-2, 'a')
            sage: Hp2 = H.change_ring(K2).odd_degree_model(); Hp2
            Hyperelliptic Curve over Number Field in a with defining polynomial x^2 + 2 defined by y^2 = 6*a*x^5 - 29*x^4 - 20*x^2 + 6*a*x + 1

            sage: K3 = QuadraticField(-3, 'b')
            sage: Hp3 = H.change_ring(QuadraticField(-3, 'b')).odd_degree_model(); Hp3
            Hyperelliptic Curve over Number Field in b with defining polynomial x^2 + 3 defined by y^2 = -4*b*x^5 - 14*x^4 - 20*b*x^3 - 35*x^2 + 6*b*x + 1

            Of course, Hp2 and Hp3 are isomorphic over the composite
            extension.  One consequence of this is that odd degree models
            reduced over "different" fields should have the same number of
            points on their reductions.  43 and 67 split completely in the
            compositum, so when we reduce we find:

            sage: P2 = K2.factor(43)[0][0]
            sage: P3 = K3.factor(43)[0][0]
            sage: Hp2.change_ring(K2.residue_field(P2)).frobenius_polynomial()
            x^4 - 16*x^3 + 134*x^2 - 688*x + 1849
            sage: Hp3.change_ring(K3.residue_field(P3)).frobenius_polynomial()
            x^4 - 16*x^3 + 134*x^2 - 688*x + 1849
            sage: H.change_ring(GF(43)).odd_degree_model().frobenius_polynomial()
            x^4 - 16*x^3 + 134*x^2 - 688*x + 1849

            sage: P2 = K2.factor(67)[0][0]
            sage: P3 = K3.factor(67)[0][0]
            sage: Hp2.change_ring(K2.residue_field(P2)).frobenius_polynomial()
            x^4 - 8*x^3 + 150*x^2 - 536*x + 4489
            sage: Hp3.change_ring(K3.residue_field(P3)).frobenius_polynomial()
            x^4 - 8*x^3 + 150*x^2 - 536*x + 4489
            sage: H.change_ring(GF(67)).odd_degree_model().frobenius_polynomial()
            x^4 - 8*x^3 + 150*x^2 - 536*x + 4489
            
        TESTS::
            sage: HyperellipticCurve(x^5 + 1, 1).odd_degree_model()
            Traceback (most recent call last):
            ...
            NotImplementedError: odd_degree_model only implemented for curves in Weierstrass form

            sage: HyperellipticCurve(x^5 + 1, names="U, V").odd_degree_model()
            Hyperelliptic Curve over Rational Field defined by V^2 = U^5 + 1
        """
        f, h = self._hyperelliptic_polynomials
        if h:
            raise NotImplementedError, "odd_degree_model only implemented for curves in Weierstrass form"
        if f.degree() % 2:
            # already odd, so just yield self
            return self

        rts = f.roots(multiplicities=False)
        if not rts:
            raise ValueError, "No odd degree model exists over field of definition"
        rt = rts[0]
        x = f.parent().gen()
        fnew =  f((x*rt + 1)/x).numerator() # move rt to "infinity"

        from constructor import HyperellipticCurve
        return HyperellipticCurve(fnew, 0, names=self._names, PP=self._PP)

    def has_odd_degree_model(self):
        r"""
        Return True if an odd degree model of self exists over the field of definition; False otherwise.

        Use ``odd_degree_model`` to calculate an odd degree model.

        EXAMPLES::
            sage: x = QQ['x'].0
            sage: HyperellipticCurve(x^5 + x).has_odd_degree_model()
            True
            sage: HyperellipticCurve(x^6 + x).has_odd_degree_model()
            True
            sage: HyperellipticCurve(x^6 + x + 1).has_odd_degree_model()
            False
        """
        try:
            return bool(self.odd_degree_model())
        except ValueError:
            return False

    def _magma_init_(self, magma):
        """
        Internal function. Returns a string to initialize this elliptic
        curve in the Magma subsystem.
        
        EXAMPLES::
        
            sage: R.<x> = QQ[]; C = HyperellipticCurve(x^3 + x - 1, x); C
            Hyperelliptic Curve over Rational Field defined by y^2 + x*y = x^3 + x - 1
            sage: magma(C)             # optional - magma
            Hyperelliptic Curve defined by y^2 + x*y = x^3 + x - 1 over Rational Field
            sage: R.<x> = GF(9,'a')[]; C = HyperellipticCurve(x^3 + x - 1, x^10); C
            Hyperelliptic Curve over Finite Field in a of size 3^2 defined by y^2 + x^10*y = x^3 + x + 2
            sage: D = magma(C); D             # optional - magma
            Hyperelliptic Curve defined by y^2 + (x^10)*y = x^3 + x + 2 over GF(3^2)
            sage: D.sage()                    # optional - magma
            Hyperelliptic Curve over Finite Field in a of size 3^2 defined by y^2 + x^10*y = x^3 + x + 2
        """
        f, h = self._hyperelliptic_polynomials
        return 'HyperellipticCurve(%s, %s)'%(f._magma_init_(magma), h._magma_init_(magma))


    def monsky_washnitzer_gens(self):
        import sage.schemes.elliptic_curves.monsky_washnitzer as monsky_washnitzer
        S = monsky_washnitzer.SpecialHyperellipticQuotientRing(self)
        return S.gens()

    def invariant_differential(self):
        """
        Returns $dx/2y$, as an element of the Monsky-Washnitzer cohomology
        of self

        EXAMPLES::

            sage: R.<x> = QQ['x']
            sage: C = HyperellipticCurve(x^5 - 4*x + 4)
            sage: C.invariant_differential()
            1 dx/2y
            
        """
        import sage.schemes.elliptic_curves.monsky_washnitzer as m_w
        S = m_w.SpecialHyperellipticQuotientRing(self)
        MW = m_w.MonskyWashnitzerDifferentialRing(S)
        return MW.invariant_differential()

    def local_coordinates_at_nonweierstrass(self, P, prec = 20, name = 't'):
        """
        For a non-Weierstrass point P = (a,b) on the hyperelliptic                         
        curve y^2 = f(x), returns (x(t), y(t)) such that (y(t))^2 = f(x(t)), 
        where t = x - a is the local parameter.

        INPUT:

        - P = (a,b) a non-Weierstrass point on self
        - prec: desired precision of the local coordinates
        - name: gen of the power series ring (default: 't')

        OUTPUT:
        (x(t),y(t)) such that y(t)^2 = f(x(t)) and t = x - a 
        is the local parameter at P

        EXAMPLES::

            sage: R.<x> = QQ['x']                                                              
            sage: H = HyperellipticCurve(x^5-23*x^3+18*x^2+40*x)                               
            sage: P = H(1,6)                                                                   
            sage: x,y = H.local_coordinates_at_nonweierstrass(P,prec=5)                        
            sage: x                                                                            
            1 + t + O(t^5)                                                                     
            sage: y                                                                            
            6 + t - 7/2*t^2 - 1/2*t^3 - 25/48*t^4 + O(t^5)                                     
            sage: Q = H(-2,12)                                                                 
            sage: x,y = H.local_coordinates_at_nonweierstrass(Q,prec=5)                        
            sage: x                                                                            
            -2 + t + O(t^5)                                                                    
            sage: y                                                                            
            12 - 19/2*t - 19/32*t^2 + 61/256*t^3 - 5965/24576*t^4 + O(t^5)

        AUTHOR:

            - Jennifer Balakrishnan (2007-12)
        """
        d = P[1]
        if d == 0:                                                                             
            raise TypeError, "P = %s is a Weierstrass point. Use local_coordinates_at_weierstrass instead!"%P          
        pol = self.hyperelliptic_polynomials()[0]
        L = PowerSeriesRing(self.base_ring(), name)
        t = L.gen()
        L.set_default_prec(prec)
        K = PowerSeriesRing(L, 'x')
        pol = K(pol)
        x = K.gen()
        b = P[0]
        f = pol(t+b)
        for i in range((RR(log(prec)/log(2))).ceil()):                
            d = (d + f/d)/2
        return t+b+O(t**(prec)), d + O(t**(prec))

    def local_coordinates_at_weierstrass(self, P,  prec = 20, name = 't'):
        """                                                                                                 
        For a finite Weierstrass point on the hyperelliptic                                
        curve y^2 = f(x), returns (x(t), y(t)) such that                                   
        (y(t))^2 = f(x(t)), where t = y is the local parameter.

        INPUT:
            - P a finite Weierstrass point on self
            - prec: desired precision of the local coordinates
            - name: gen of the power series ring (default: 't')

        OUTPUT:

        (x(t),y(t)) such that y(t)^2 = f(x(t)) and t = y
        is the local parameter at P                                             
                                                                                                            
        EXAMPLES:                                                                                           
            sage: R.<x> = QQ['x']                                                              
            sage: H = HyperellipticCurve(x^5-23*x^3+18*x^2+40*x)               
            sage: A = H(4,0)                                                                   
            sage: x,y = H.local_coordinates_at_weierstrass(A,prec =5)                          
            sage: x
            4 + 1/360*t^2 - 191/23328000*t^4 + 7579/188956800000*t^6 + O(t^7)
            sage: y 
            t + O(t^7)
            sage: B = H(-5,0)
            sage: x,y = H.local_coordinates_at_weierstrass(B,prec = 5)
            sage: x
            -5 + 1/1260*t^2 + 887/2000376000*t^4  + 643759/1587898468800000*t^6 + O(t^7)
            sage: y
            t + O(t^7)
                                                                                                            
        AUTHOR:                                                                                             
            - Jennifer Balakrishnan (2007-12)
        """
        if P[1] != 0:                                                                          
            raise TypeError, "P = %s is not a finite Weierstrass point. Use local_coordinates_at_nonweierstrass instead!"%P
        pol = self.hyperelliptic_polynomials()[0]
        L = PowerSeriesRing(self.base_ring(), name)
        t = L.gen()
        L.set_default_prec(prec+2)
        K = PowerSeriesRing(L, 'x')
        pol = K(pol)
        x = K.gen()
        b = P[0]
        g = pol/(x-b)
        c = b+1/g(b)*t**2
        f = pol - t**2
        fprime = f.derivative()
        for i in range((RR(log(prec+2)/log(2))).ceil()):
            c = c-f(c)/fprime(c)
        return c + O(t**(prec+2)),t+O(t**(prec+2))
    
    def local_coordinates_at_infinity(self, prec = 20, name = 't'):
        """
        For the genus g hyperelliptic curve y^2 = f(x), returns (x(t), y(t)) such that
        (y(t))^2 = f(x(t)), where t = x^g/y is the local parameter at infinity

        INPUT:
            - prec: desired precision of the local coordinates
            - name: gen of the power series ring (default: 't')

        OUTPUT:
        (x(t),y(t)) such that y(t)^2 = f(x(t)) and t = x^g/y
        is the local parameter at infinity


        EXAMPLES:
            sage: R.<x> = QQ['x']
            sage: H = HyperellipticCurve(x^5-5*x^2+1)
            sage: x,y = H.local_coordinates_at_infinity(10)
            sage: x
            t^-2 + 5*t^4 - t^8 - 50*t^10 + O(t^12)
            sage: y 
            t^-5 + 10*t - 2*t^5 - 75*t^7 + 50*t^11 + O(t^12)

            sage: R.<x> = QQ['x']
            sage: H = HyperellipticCurve(x^3-x+1)
            sage: x,y = H.local_coordinates_at_infinity(10)
            sage: x
            t^-2 + t^2 - t^4 - t^6 + 3*t^8 + O(t^12)
            sage: y
            t^-3 + t - t^3 - t^5 + 3*t^7 - 10*t^11 + O(t^12)


        AUTHOR:
            - Jennifer Balakrishnan (2007-12)
        """
        g = self.genus()
        pol = self.hyperelliptic_polynomials()[0]
        K = LaurentSeriesRing(self.base_ring(), name)
        t = K.gen()
        K.set_default_prec(prec+2)
        L = PolynomialRing(self.base_ring(),'x')
        x = L.gen()
        i = 0
        w = (x**g/t)**2-pol
        wprime = w.derivative(x)
        x = t**-2
        for i in range((RR(log(prec+2)/log(2))).ceil()):
            x = x-w(x)/wprime(x)
        y = x**g/t
        return x+O(t**(prec+2)) , y+O(t**(prec+2))


    def local_coord(self, P, prec = 20, name = 't'):
        """                                                                                                          
        Calls the appropriate local_coordinates function

        INPUT:
            - P a point on self
            - prec: desired precision of the local coordinates
            - name: gen of the power series ring (default: 't')

        OUTPUT:
        (x(t),y(t)) such that y(t)^2 = f(x(t)), where t
        is the local parameter at P

        EXAMPLES:
            sage: R.<x> = QQ['x']
            sage: H = HyperellipticCurve(x^5-23*x^3+18*x^2+40*x)
            sage: H.local_coord(H(1,6),prec=5)
            (1 + t + O(t^5), 6 + t - 7/2*t^2 - 1/2*t^3 - 25/48*t^4 + O(t^5))
            sage: H.local_coord(H(4,0),prec=5)
            (4 + 1/360*t^2 - 191/23328000*t^4 + 7579/188956800000*t^6 + O(t^7), t + O(t^7))
            sage: H.local_coord(H(0,1,0),prec=5)
            (t^-2 + 23*t^2 - 18*t^4 - 569*t^6 + O(t^7), t^-5 + 46*t^-1 - 36*t - 609*t^3 + 1656*t^5 + O(t^6))
        
        AUTHOR:
            - Jennifer Balakrishnan (2007-12)

                                                                                                    
        """
        if P[1] == 0:
            return self.local_coordinates_at_weierstrass(P, prec, name)
        elif P[2] == 0:
            return self.local_coordinates_at_infinity(prec, name)
        else:
            return self.local_coordinates_at_nonweierstrass(P, prec, name)