source: sage/schemes/hyperelliptic_curves/jacobian_morphism.py @ 5447:cd60620ddfe5

Revision 5447:cd60620ddfe5, 5.9 KB checked in by David Kohel <kohel@…>, 6 years ago (diff)

Fixes to left multiplication and coercion of 0 in Jacobians and point sets.
Line 
1"""
2Jacobian ``morphism'' as a class in the Picard group
3"""
4
5#*****************************************************************************
6#  Copyright (C) 2005 David Kohel <kohel@maths.usyd.edu.au>
7#  Distributed under the terms of the GNU General Public License (GPL)
8#                  http://www.gnu.org/licenses/
9#*****************************************************************************
10
11from sage.rings.all import PolynomialRing, ZZ
12from sage.structure.element import AdditiveGroupElement
13from sage.schemes.generic.morphism import SchemeMorphism
14
15def cantor_reduction_simple(a1,b1,f,genus):
16    # Divisor reduction.
17    a2 = (f - b1**2).div(a1)
18    a2 *= 1/a2.leading_coefficient()
19    b2 = -b1.mod(a2);
20    if a2.degree() == a1.degree():
21        assert a2.degree() == genus+1
22        print "Returning ambiguous form of degree genus+1."
23        return (a2, b2)
24    elif a2.degree() > genus:
25        return cantor_reduction_simple(a2,b2,f)
26    return (a2, b2)
27
28def cantor_reduction(a,b,f,h,genus):
29    assert a.degree() < 2*genus+1
30    assert b.degree() < a.degree()
31    k = f - h*b - b**2
32    if 2*a.degree() == k.degree():
33        # must adjust b to include the point at infinity
34        g1 = a.degree()
35        x = a.parent().gen()
36        r = (x**2 + h[g1]*x - f[2*g1]).roots()[0][0]
37        b = b + r*(x**g1 - (x**g1).mod(a))
38        k = f - h*b - b**2
39    assert k.mod(a) == 0
40    a = k.div(a)
41    a /= a.leading_coefficient()
42    b = -(b+h).mod(a)
43    if a.degree() > genus:
44        return cantor_reduction(a,b,f,h,genus)
45    return (a, b)
46
47def cantor_composition_simple(D1,D2,f,genus):
48    a1, b1 = D1
49    a2, b2 = D2
50    if a1 == a2 and b1 == b2:
51        # Duplication law:
52        d, h1, h3 = a1.xgcd(2*b1)
53        a = (a1.div(d))**2
54        b = (b1 + h3*((f - b1**2).div(d))).mod(a)
55    else:
56        d0, _, h2 = a1.xgcd(a2)
57        if d0 == 1:
58            a = a1*a2
59            b = (b2 + h2*a2*(b1-b2)).mod(a)
60        else:
61            d, l, h3 = d0.xgcd(b1 + b2)
62            a = (a1*a2).div(d**2)
63            b = (b2 + l*h2*(b1-b2)*a2.div(d)) + h3*((f - b2**2).div(d)).mod(a)
64    if a.degree() > genus:
65        return cantor_reduction_simple(a,b,f,genus)
66    return (a,b)
67
68def cantor_composition(D1,D2,f,h,genus):
69    a1, b1 = D1
70    a2, b2 = D2
71    if a1 == a2 and b1 == b2:
72        # Duplication law:
73        d, h1, h3 = a1.xgcd(2*b1 + h)
74        # NOTE THAT d is not normalised, but this gives a crash:
75        # d *= 1/d.leading_coefficient()
76        # print "d =", d
77        a = (a1.div(d))**2;
78        b = (b1 + h3*((f-h*b1-b1**2).div(d))).mod(a)
79    else:
80        d0, _, h2 = a1.xgcd(a2)
81        if d0 == 1:
82            a = a1*a2;
83            b = (b2 + h2*a2*(b1-b2)).mod(a)
84        else:
85            e0 = b1+b2+h
86            if e0 == 0:
87                a = (a1*a2).div(d0**2);
88                b = (b2 + h2*(b1-b2)*(a2.div(d0))).mod(a)
89            else:
90                d, l, h3 = d0.xgcd(e0)
91                a = (a1*a2).div(d**2);
92                b = (b2 + l*h2*(b1-b2)*(a2.div(d)) + h3*((f-h*b2-b2**2).div(d))).mod(a)
93    a *= 1/a.leading_coefficient()
94    if a.degree() > genus:
95        return cantor_reduction(a,b,f,h,genus)
96    return (a,b)
97
98class JacobianMorphism_divisor_class_field(AdditiveGroupElement, SchemeMorphism):
99    r"""
100    An element of a $J(K) = \Pic^0_K(C)$.
101    """
102    def __init__(self, parent, polys, reduce=True, check=False):
103        SchemeMorphism.__init__(self, parent)
104        if polys == 0:
105            P = PolynomialRing(self.parent(), 'x')
106            self.__polys = (P(1),P(0))
107        if check:
108            C = parent.curve()
109            K = parent.value_ring()
110            f, h = C.hyperelliptic_polynomials(K)
111            a, b = polys
112            if not (b**2 + h*b - f)%a == 0:
113                raise ValueError, \
114                      "Argument polys (= %s) must be reduced divisor on curve %s."%(
115                    polys, C)
116        self.__polys = polys
117       
118    def __repr__(self):
119        a, b = self.__polys
120        if a == 1:
121            return "(1)"
122        P = self.parent()._printing_ring
123        y = P.gen()
124        x = P.base_ring().gen()
125        return "(%s, %s)"%(a(x), y - b(x))
126
127    def scheme(self):
128        return self.codomain()
129
130    def list(self):
131        return self.__polys
132
133    def __add__(self,other):
134        X = self.parent()
135        C = X.curve()
136        K = X.value_ring()
137        f, h = C.hyperelliptic_polynomials(K)
138        if h == 0:
139            D = cantor_composition_simple(self.__polys,other.__polys,f,C.genus())
140        else:
141            D = cantor_composition(self.__polys,other.__polys,f,h,C.genus())
142        return JacobianMorphism_divisor_class_field(X, D, reduce=False, check=False)
143
144    def __cmp__(self, other):
145        if not isinstance(other, JacobianMorphism_divisor_class_field):
146            try:
147                other = self.parent()(other)
148            except TypeError:
149                return -1
150        return cmp(self.__polys, other.__polys)
151
152    def __nonzero__(self):
153        return self.__polys[0] != 1
154
155    def __sub__(self, other):
156        return self + (-other)
157
158    def __neg__(self):
159        if self.is_zero():
160            return self
161        polys = self.__polys
162        X = self.parent()
163        K = X.value_ring()
164        f, h = X.curve().hyperelliptic_polynomials(K)
165        if h == 0:
166            D = (polys[0],-polys[1])
167        else:
168            D = (polys[0],-polys[1]-h.mod(polys[0]))
169        return JacobianMorphism_divisor_class_field(X, D, reduce=False, check=False)
170       
171    def __mul__(self, n):
172        try:
173            n = ZZ(n)
174        except TypeError:
175            raise TypeError, "Argument n (= %s) must be an integer."
176        X = self.parent()
177        if n < 0:
178            return - self * (-n)
179        elif n == 0:
180            return self.parent()(0)
181        elif n == 1:
182            return self
183        elif n < 4:
184            D = self
185        else:       
186            D = self.__mul__(n//2)
187        if n % 2 == 0:
188            return D + D
189        else:     
190            return D + D + self
191
192    def _rmul_(self, n):
193        return self.__mul__(n)
194
195
Note: See TracBrowser for help on using the repository browser.