Ticket #14756: 14756-sc-reduction.patch

File 14756-sc-reduction.patch, 35.0 KB (added by Marco Streng, 9 years ago)
  • new file sage/schemes/hyperelliptic_curves/stoll_cremona.py

    # HG changeset patch
    # User Marco Streng <marco.streng@gmail.com>
    # Date 1371474259 -7200
    # Node ID 2e8268247f33914a0ce66c06ce6a1bfe2ee2889b
    # Parent  602ccb7e0480afb988289e17adbb2e6695bbf17a
    Reduction of hyperelliptic curves via a covariant
    
    diff --git a/sage/schemes/hyperelliptic_curves/stoll_cremona.py b/sage/schemes/hyperelliptic_curves/stoll_cremona.py
    new file mode 100644
    - +  
     1r"""
     2Reduction of binary forms and hyperelliptic curve equations
     3
     4This file contains functions that reduce the coefficients of a binary form or
     5hyperelliptic curve by integral transformations, leaving the discriminant
     6intact.
     7
     8AUTHORS:
     9
     10- Florian Bouyer (2011 - 09) : covariant and stoll_cremona_reduction
     11- Marco Streng (2012 - 09) : Hilbert fundamental domains in some cases
     12
     13"""
     14#*****************************************************************************
     15#       Copyright (C) 2011,2012,2013 Florian Bouyer <f.j.s.c.bouyer@gmail.com>
     16#                                and Marco Streng <marco.streng@gmail.com>
     17#
     18#  Distributed under the terms of the GNU General Public License (GPL)
     19#  as published by the Free Software Foundation; either version 2 of
     20#  the License, or (at your option) any later version.
     21#                  http://www.gnu.org/licenses/
     22#*****************************************************************************
     23
     24from sage.structure.sage_object import SageObject
     25from sage.matrix.all import Matrix
     26from sage.schemes.plane_conics.constructor import Conic
     27from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
     28from sage.rings.arith import gcd, primes_first_n, factor
     29from sage.rings.all import (ZZ, QQ, CC, ComplexField, RealField,
     30                            is_NumberField, FiniteField)
     31from sage.schemes.hyperelliptic_curves.constructor import HyperellipticCurve
     32from sage.misc.misc import prod
     33from sage.matrix.constructor import identity_matrix
     34
     35def stoll_cremona_reduction(f, degree=None, algorithm='default',
     36                            precision=53, group='gl', transformation=False):
     37    r"""
     38    Reduces the polynomial f using the Stoll-Cremona method
     39   
     40    INPUT:
     41   
     42    - ``f`` - a univariate polynomial
     43   
     44    - ``degree`` - a positive integer (default: None). Interpret ``f`` as a
     45      homogeneous polynomial of degree ``degree`` in two variables, i.e.,
     46      as `F(X,Y) = Y^{degree} f(X/Y)`. If None, use the degree of `f`.
     47   
     48    - ``algorithm`` - ('default' or 'magma'). If set to 'magma', instead of
     49      using :func:`covariant_z0` to find the covariant `z_0`, the algorithm
     50      uses Magma to compute the 'better' covariant `z` (see [SC2002] for more
     51      details) (requires Magma to be installed)
     52     
     53    - ``precision`` - a positive integer (default=53). The precision to use
     54      for the covariants.
     55
     56    - ``group`` (default='sl') -- a string: 'sl' for SL_2(O_K), 'gl' for
     57      GL_2(O_K), 'plus' for the subgroup GL_2(O_K)^+ of elements of totally
     58      positive determinant in GL_2(O_K).
     59   
     60    - ``transformation`` (default: False). Whether to also output a matrix
     61      (a,b;c,d) such that the output g satisfies
     62      g(z) = (cz+d)^degree * f((az+b)/(cz+d)).
     63     
     64    OUTPUT:
     65       
     66    A polynomial whose associated covariant point in the upper half plane
     67    is in a fundamental domain and that is `SL_2(ZZ)`-equivalent to f
     68   
     69    EXAMPLES:
     70   
     71    We check with Stoll and Cremona's examples::
     72       
     73        sage: from sage.schemes.hyperelliptic_curves.stoll_cremona import stoll_cremona_reduction
     74        sage: P.<x> = QQ[]
     75        sage: f = 19*x^8 - 262*x^7 + 1507*x^6 - 4784*x^5 + 9202*x^4 - 10962*x^3 + 7844*x^2 - 3040*x + 475
     76        sage: stoll_cremona_reduction(f)
     77        -x^8 + 6*x^7 - 7*x^6 - 12*x^5 + 27*x^4 - 4*x^3 - 19*x^2 + 10*x - 5
     78       
     79        sage: Q.<y> = QQ[]
     80        sage: f = y^6 + 30*y^5 + 371*y^4 + 2422*y^3 + 8813*y^2 + 16968*y + 13524
     81        sage: stoll_cremona_reduction(f)   
     82        y^6 - 4*y^4 + 2*y^3 + 8*y^2 - 12*y + 9
     83       
     84    An example over a real quadratic number field::
     85   
     86        sage: K.<a> = NumberField(x^2-5)
     87        sage: P.<x> = K[]
     88        sage: f = (x+a)^6+7; f
     89        x^6 + 6*a*x^5 + 75*x^4 + 100*a*x^3 + 375*x^2 + 150*a*x + 132
     90        sage: stoll_cremona_reduction(f)
     91        x^6 + 7
     92
     93    The chosen degree is relevant, in the following example, the Igusa-Clebsch
     94    invariants show that the reduction for degree 5 is not equivalent for
     95    degree 6::
     96   
     97        sage: Q.<x> = QQ[]
     98        sage: f = 4*x^5 + 26*x^4 + 88*x^3 + 144*x^2 + 112*x + 32
     99        sage: C = HyperellipticCurve(f)
     100        sage: i = C.igusa_clebsch_invariants()
     101        sage: g = stoll_cremona_reduction(f); g
     102        2*x^5 + 4*x^4 + 4*x^3 + 24*x^2 - 6*x + 4
     103       
     104        sage: D = HyperellipticCurve(g)
     105        sage: j = D.igusa_clebsch_invariants()
     106        sage: i[0]^5/i[3] == j[0]^5/j[3]
     107        False
     108       
     109        sage: stoll_cremona_reduction(f, degree=6)
     110        Traceback (most recent call last):
     111        ...
     112        NotImplementedError: Please implement covariant_z0 for binary forms with roots at infinity
     113       
     114        sage: h = stoll_cremona_reduction(f, degree=6, algorithm='magma'); h # optional - magma
     115        4*x^5 + 6*x^4 + 24*x^3 - 4*x^2 + 4*x - 2
     116        sage: D = HyperellipticCurve(h) # optional - magma
     117        sage: j = D.igusa_clebsch_invariants() # optional - magma
     118        sage: i[0]^5/i[3] == j[0]^5/j[3]
     119        True
     120
     121
     122    With Magma, we can get a better reduction, as the covariant `z` is
     123    implemented there, rather than just the covariant `z0`::
     124   
     125        sage: P.<x> = QQ[]
     126        sage: f = good example, maybe from [SC2002]
     127        sage: stoll_cremona_reduction(f)
     128       
     129        sage: stoll_cremona_reduction(f, algorithm='magma') # optional - magma
     130       
     131    An example with ``transformation=True``::
     132   
     133        sage: g = 32*x^6 + 688*x^5 + 6144*x^4 + 29176*x^3 + 77714*x^2 + 110104*x + 64830
     134        sage: h, T = stoll_cremona_reduction(g, algorithm='magma', transformation=True)
     135        sage: T
     136        [ 3  4]
     137        [-1 -1]
     138        sage: g((3*x+4)/(-x-1))*(-x-1)^6 == h
     139        True
     140
     141    ALGORITHM:
     142   
     143    This is the algorithm from [SC2002].
     144   
     145    If the base field is `QQ`, then the algorithm is exactly the same. For
     146    real quadratic fields, we implemented a reduction in the Hilbert upper
     147    half space HH^2. If the real quadratic field has class number one, then
     148    this region is the standard fundamental domain. For more general number
     149    fields, a more general reduction algorithm for the covariant point is
     150    needed.
     151   
     152    REFERENCES
     153   
     154    [SC2002] M. Stoll and J. Cremona "On the reduction theory of binary forms"
     155             J. Reine Angew. Math. 565 (2003) pp. 79-99
     156    """
     157#   s = []
     158    K = f.base_ring()
     159    P = f.parent()
     160    x = P.gen()
     161    h = _HilbertFundDomData(K, embs=precision, group=group)
     162    if degree is None:
     163        degree = f.degree()
     164    R = RealField(precision)
     165    C = ComplexField(precision)
     166    T = Matrix([[1,0],[0,1]])
     167    previous_height = R(-2)
     168    height = R(-1)
     169    one_minus_epsilon = R(1) - R(2)**(-5)   
     170    if algorithm == 'default':
     171        while previous_height < height * one_minus_epsilon:
     172            previous_height = height
     173            z = covariant_z0(f, precision=precision, degree=degree)
     174            M, z = h.step(z)
     175            [[a,b],[c,d]] = M
     176            f = P(f((d*x-b)/(-c*x+a)) * (-c*x+a)**degree)
     177#           s.append((f,z,covariant_z0(f)))
     178            T = T * Matrix([[d, -b], [-c, a]])
     179            height = prod([t.imag() for t in covariant_z0(f)]).abs()
     180            print f
     181            print height
     182            print previous_height
     183        if height < previous_height * one_minus_epsilon:
     184            raise RuntimeError, "Possibly too low precision in stoll_cremona_reduction"
     185    elif algorithm == 'magma':
     186        from sage.interfaces.all import magma
     187        z = 0; zp = 1
     188        while z != zp:
     189            Phi = K.embeddings(R)
     190            z = []
     191            for emb in Phi:
     192                gcoeffs = [emb(f[k]) for k in range(degree+1)]
     193                P = R['x']
     194                x = P.gen()
     195                g = sum([gcoeffs[indx]*x**indx for indx in range(degree+1)])
     196                zi = magma(magma.Covariant(degree, magma(g))).sage()
     197                z.append(C(zi))
     198            M,zp = h.step(z)
     199            [a,b],[c,d] = M
     200            P = K['x']
     201            x = P.gen()
     202            f = f.parent()(f((d*x-b)/(-c*x+a)) * (-c*x+a)**degree)
     203            T = T * Matrix([[d, -b], [-c, a]])
     204    else:
     205        raise ValueError("Unkown algorithm '%s'" % algorithm)   
     206    if transformation:
     207        return f, T
     208    return f
     209
     210
     211def covariant_z0(f, precision=53, degree=None):   
     212    r"""
     213    Returns the covariant `z_0` of f in H (upper half plane) from [SC2002].
     214   
     215    INPUT:
     216       
     217     - ``f`` - a polynomial
     218   
     219     - ``precision`` - An integer (Default = 53) What precision for real and
     220       complex numbers should the algorithm use
     221       
     222     - ``degree`` - the degree to homogenize with
     223   
     224    OUTPUT:
     225   
     226    A point z in the upper half plane associated to f. If f is defined over
     227    a totally real number field of degree g, then the output is a sequence
     228    of length g.
     229   
     230    EXAMPLES::
     231   
     232        sage: z0 = sage.schemes.hyperelliptic_curves.reduction.covariant_z0
     233        sage: f = x^6 + 3*x^3 +4*x^2+ 3
     234        sage: z0(f)
     235        -0.0478810014556224 + 1.20469810839717*I
     236        sage: -z0(f)^-1
     237       
     238        sage: z0(f.reverse())
     239       
     240        sage: z0(f)-1
     241       
     242        sage: z0(f(x+1))
     243       
     244    Comparing the different answers using different precision ::
     245
     246        sage: z0 = sage.schemes.hyperelliptic_curves.reduction.covariant_z0
     247        sage: P.<x> = QQ[]
     248        sage: f = x^8 + 24*x^7 + 3*x^5 +4*x^2+ 3
     249        sage: z0(f)
     250        -0.00832825511768667 + 0.835323919995269*I
     251        sage: z0(f,200)
     252        -0.0083282551176867045891624031810609648523736031152551979546157 + 0.83532391999526880130877046174257268976426032236105109987080*I
     253   
     254    If f is not in `\QQ` then ``covariant_z0`` returns a list where every entry
     255    is a z associated to one embeddings of f in `\RR`::
     256   
     257        sage: z0 = sage.schemes.hyperelliptic_curves.reduction.covariant_z0
     258        sage: X = var('X')
     259        sage: k = NumberField(X^2-41,'a')
     260        sage: P.<x> = k[]
     261        sage: a = k.an_element()
     262        sage: f = x^8 + 24*x^7 + 3*x^5 +4*x^2+ 3
     263        sage: z0(f)
     264        [-0.00832825511768667 + 0.835323919995269*I, -0.00832825511768667 + 0.835323919995269*I]
     265        sage: f = x^6 + (a+1)*x^5 + (5*a+2)*x^2 + (a/2+1/2)*x + 5*a
     266        [0.167978514759694 + 1.68636006242584*I, -0.125921323440448 + 1.64323412027265*I]
     267
     268    Even though Stoll and Cremona allow roots of higher multiplicity, and roots
     269    at infinity (i.e. degree > f.degree()), this is currently not fully
     270    supported::
     271       
     272        sage: z0 = sage.schemes.hyperelliptic_curves.reduction.covariant_z0
     273        sage: P.<x> = RR[]
     274        sage: i = (x-3)^2 * (x^2 + 2)^2
     275        sage: z0(i)
     276        Traceback (most recent call last):
     277        ...
     278        ValueError: Cannot convert NaN or infinity to Pari float
     279       
     280        sage: j = x^8 + 3
     281        sage: z0(j, degree=9)
     282        Traceback (most recent call last):
     283        ...
     284        NotImplementedError
     285       
     286    ALGORITHM:
     287       
     288    This algorithm is taken straight from [SC2002] page 6
     289   
     290    Let n be the degree of the polynomial f. For each embeddings of f into
     291    `\RR` let `\alpha_j` be the roots of f. Define a quadratic polynomial
     292    `Q_0` by `Q_0(x) = \sum_{j=1}^{n}
     293                       (x-\alpha_j)(x-\bar{\alpha_j})/|f'(\alpha_j)|^{2/(n-2)}`
     294    Finally take `z_0` to be the root with positive imaginary part of `Q_0`
     295   
     296    REFERENCE:
     297   
     298    [SC2002] M. Stoll and J. Cremona "On the reduction theory of binary forms"
     299             J. Reine Angew. Math. 565 (2003) pp. 79-99
     300    """   
     301    C_out = ComplexField(precision)
     302 
     303    internal_precision = 6*precision # TODO: make an informed choice
     304    C = ComplexField(internal_precision)
     305    R = RealField(internal_precision)
     306   
     307
     308    if degree is None:
     309        degree = f.degree()
     310    elif degree != f.degree():
     311        raise NotImplementedError("Please implement covariant_z0 for "\
     312                                  "binary forms with roots at infinity")
     313   
     314    k = f.base_ring()
     315    x = R['x'].gen()
     316    if k is QQ or is_NumberField(k):
     317        Phi = k.embeddings(C)
     318    else:
     319        Phi = [C.convert_map_from(k)]
     320    I = C.gen()
     321   
     322    # Actually, the roots should be taken with multiplicities,
     323    # but this does not matter, as the relevant coefficients
     324    # are divided by infinity anyway.
     325   
     326    z_list = []
     327    for emb in Phi:
     328        gcoeffs = [emb(coeff) for coeff in f.coeffs()]
     329        g = sum([gcoeffs[indx] * x**indx for indx in range(degree+1)])
     330        roots = g.roots(C, multiplicities=False)
     331        gdiv = g.derivative()
     332        q0 = R(0)
     333        for ai in roots:
     334            q0 += ((x - ai) * (x - ai.conjugate()) /
     335                   (abs(gdiv(ai))**(2 / (degree - 2))))
     336        z1,z2 = q0.roots(C, multiplicities=False)
     337        if z1.imag_part() > 0:
     338            z0 = z1
     339        else:
     340            z0 = z2
     341        z_list.append(z0)
     342
     343    return [C_out(z) for z in z_list]
     344
     345
     346def _HilbertFundDomData(K, embs=None, group='sl'):
     347    """
     348    Create an object containing the data for a fundamental domain
     349    for the action of SL_2(O_K) on (CC\setminus RR)^g for a totally real
     350    field K of degree g.
     351
     352    Warning: assumes deg(K) <= 2 and cl(K) = 1.
     353
     354    Input:
     355
     356     - ``K`` -- A totally real number field.
     357
     358     - ``embs`` (default=None) -- None, or the complete set
     359       of embeddings of K into a real field that should be used,
     360       or the precision to use for creating such a set of embeddings.
     361
     362     - ``group`` (default='sl') -- a string: 'sl' for SL_2(O_K),
     363       'gl' for GL_2(O_K), 'plus' for
     364       the subgroup GL_2(O_K)^+ of elements of totally positive
     365       determinant in GL_2(O_K).
     366       
     367    EXAMPLES::
     368       
     369        sage: from sage.schemes.hyperelliptic_curves.stoll_cremona import _HilbertFundDomData
     370        sage: K = QuadraticField(41,'a')
     371        sage: K(K.unit_group().gen(1)).norm()
     372        -1
     373        sage: z = [0.01*CC.gen()+23.34,0.02+0.03*CC.gen()]
     374        sage: h = _HilbertFundDomData(K, group='sl')
     375        sage: r = h.reduce(z); r
     376        (
     377        [    9/2*a - 11/2 -125/2*a + 801/2]                                     
     378       
     379        [   -3/2*a - 17/2         2*a - 13], [1.03241154280086 +
     380        15.5033759282630*I, -1.33856972357549 + 0.0496716993976175*I]
     381        )
     382        sage: prod([e.imag() for e in r[1]]), max([log(abs(e.imag())) for e in r[1]])
     383        (0.770079028756939, 2.74105780203325)
     384        sage: h = _HilbertFundDomData(K, group='gl')
     385        sage: r = h.reduce(z); r
     386        (
     387        [15/2*a + 99/2  5/2*a - 37/2]                                           
     388       
     389        [-3/2*a - 17/2      2*a - 13], [1.83799715597054 + 0.242181137053011*I,
     390        -2.13249730501399 - 3.17976469236058*I]
     391        )
     392        sage: prod([e.imag() for e in r[1]]), max([log(abs(e.imag())) for e in r[1]])
     393        (-0.770079028756901, 1.15680719794121)
     394    """
     395    if K.degree() == 1:
     396        return _HilbertFundDomData_classical(embs)
     397    if K.degree() == 2:
     398        return _HilbertFundDomData_quadratic(K, embs, group)
     399    return _HilbertFundDomData_generic(K, embs, group)
     400
     401
     402class _HilbertFundDomData_generic(SageObject):
     403    """
     404    An object containing the data for a fundamental domain
     405    for the action of SL_2(O_K) on HH^g for a totally real
     406    field K of degree g.
     407    """
     408
     409    _K = None
     410    _embs = None
     411    _RR = None
     412
     413    def __init__(self, K, embs, group):
     414        """
     415        Create the object.
     416        """
     417        self._K = K
     418        if embs == None:
     419            embs = K.embeddings(RR)
     420            self._RR = RR
     421        elif embs in ZZ:
     422            self._RR = RealField(embs)
     423            embs = K.embeddings(self._RR)
     424        else:
     425            self._RR = embs[0].codomain()
     426        self._embs = embs
     427
     428    def to_additive(self, x):
     429        """
     430        Move x in RR^g to an additive fundamental domain.
     431
     432        Returns an pair (a, x') with a in O_K and x' = x-a in the fundamental domain.
     433        """
     434        raise NotImplementedError, "Sorry reduction to a Hilbert fundamental domain is not implemented for g>2"
     435
     436    def step(self, z, bound_on_c=None):
     437        """
     438        Do one reduction round.
     439
     440        Returns a pair (M, z') with z' = M(z) nearer the fundamental domain if possible.
     441        """
     442        M1 = self.up(z, bound_on_c=bound_on_c)
     443        z = self.act(M1, z)
     444        M2 = self.to_multiplicative([e.imag() for e in z])[0]
     445        z = self.act(M2, z)
     446        t = self.to_additive([e.real() for e in z])[0]
     447        M3 = Matrix([[1, -t], [0, 1]])
     448        z = self.act(M3, z)
     449        return M3*M2*M1, z
     450
     451    def reduce(self, z, max_rounds = 0, bound_on_c=None):
     452        """
     453        Do reduction rounds until reduced.
     454
     455        Returns a pair (M, z') with z' = M(z) in the fundamental domain.
     456
     457        max_rounds = 0: go on until finishes, risking an infinite loop from bugs or rounding errors
     458        max_rounds > 0: do at most that number of rounds
     459        """
     460        M = Matrix([[1,0],[0,1]])
     461        rounds = 0
     462        while true:
     463            rounds = rounds + 1
     464            M1, z = self.step(z, bound_on_c=bound_on_c)
     465            M = M1 * M
     466            if M1 == 1 or rounds == max_rounds:
     467                return M, z
     468
     469    def act(self, M, z):
     470        """
     471        Return M(z)
     472        """
     473        a = M[0,0]
     474        b = M[0,1]
     475        c = M[1,0]
     476        d = M[1,1]
     477        embs = self._embs
     478        return [((embs[k](a)*z[k] + embs[k](b))/(embs[k](c)*z[k] + embs[k](d))) for k in range(len(embs))]
     479
     480    def _up_cd_lll(self, z):
     481        """
     482        As _up_cd, but instead of making prod(Im(z)) maximal
     483        and log(Im(z)) in a fundamental domain, try to make min(Im(z)) large
     484        using LLL.
     485        """
     486        # The new imaginary part Im(z) / |cz+d|^2
     487        # In particular, max(|cz+d|^2/Im(z)) must be as small as possible,
     488        # and less than 1 / mim(Im(z)).
     489        # Let v = ((cz+d)/sqrt(Im(z)))_i.
     490        # Then we want |v|_infty to be as small as possible and less than
     491        # 1 / min(sqrt(Im(z))).
     492        # Note, g*|v|_2 >= g*|v|_infty^2 >= |v|_2, so it is necessary that
     493        # |v|_2 < g/(min(Im(z))) and sufficient that |v|_2 < 1/(min(Im(z)))
     494        g = self._K.degree()
     495        embs = self._embs
     496        m = 2**(z[0].parent().precision())
     497        O = self._K.maximal_order()
     498        vectors = ( [ [(embs[i](b)) for i in range(g)] +
     499        [0 for i in range(g)] for b in O.basis()] +
     500        [ [(embs[i](b)*z[i].real()) for i in range(g)] +
     501        [ (embs[i](b)*z[i].imag()) for i in range(g)]
     502        for b in O.basis()])
     503        id = identity_matrix(2*g)
     504        vectors = [id[j].list() +
     505                   [(m*vectors[j][i]/(z[i%g].imag().abs().sqrt())).round()
     506                      for i in range(len(vectors[j]))] for j in range(2*g)]
     507#        c = []
     508#        for i in range(len(vectors)):
     509#            for j in range(i, len(vectors)):
     510#                c.append(sum([vectors[i][k]*vectors[j][k] for k in range(2*g)]))
     511#        Q = QuadraticForm(ZZ, 2*g, c)
     512#        l = Q.short_vector_list_up_to_length(g/min([z.imag().abs()]))
     513#        for x in l[1:]:
     514        M = Matrix(vectors).LLL()
     515        d = sum([M[0][i]*O.basis()[i] for i in range(g)])
     516        c = sum([M[0][g+i]*O.basis()[i] for i in range(g)])
     517        B = c * O + d * O
     518        if B == 1*O:
     519            return c, d
     520        if B.is_principal():
     521            g = B.gens_reduced()[0]
     522            return c/g, d/g
     523        # TODO: if B is not principal, find another candidate pair c, d
     524
     525    def up(self, z, bound_on_c=None):
     526        """
     527        Move z in HH^g, i.e., increase the norm of its imaginary part,
     528        if possible.
     529
     530        Returns either None or a matrix M in SL_2(O_K) that such that
     531        Mz has a larger norm of its imaginary part.
     532
     533        Assumes K has class number one.
     534
     535        EXAMPLES::
     536
     537            sage: from sage.schemes.hyperelliptic_curves.stoll_cremona import _HilbertFundDomData
     538            sage: h = _HilbertFundDomData(QuadraticField(15))
     539            sage: h.up([123+CC.0*0.05, 456+CC.0*0.0234])
     540            [          0          -1]
     541            [          1 -43*a - 290]
     542            sage: h.act(_, [123+CC.0*0.05, 456+CC.0*0.0234])
     543            [2.14072851429831 + 0.231823024326897*I, 1.85425166669223 + 0.0806070738044104*I]
     544            h.reduce([123+CC.0*0.05, 456+CC.0*0.0234])
     545            (
     546            [      2*a + 3 -708*a - 2158]
     547            [           -2    86*a + 579], [-0.0404815976687969 + 3.15206362011738*I, -0.618912376653362 + 2.90579848705465*I]
     548            )
     549        """
     550        cd = self._up_cd(z, bound_on_c=bound_on_c)
     551        if cd == None:
     552            return Matrix([[1,0],[0,1]])
     553        c, d = cd
     554        a = d.inverse_mod(c)
     555        b = (a*d-1)/c
     556        return Matrix([[a,b],[c,d]])
     557
     558    def _up_cd(self, z):
     559        """
     560        Move z up in HH^g, i.e., increase the norm of its imaginary part,
     561        if possible. Actually, not implemented in this base class: just calls
     562        _up_cd_lll instead.
     563        """
     564        return self._up_cd_lll(z)
     565
     566
     567class _HilbertFundDomData_quadratic(_HilbertFundDomData_generic):
     568    """
     569    As _HilbertFundDomData, but only for quadratic fields.
     570    """
     571
     572    _eps = None
     573    _epsmat = None
     574    _epslog = None
     575    _RR = None
     576
     577    def __init__(self, K, embs, group):
     578        """
     579        Create the object.
     580        """
     581        _HilbertFundDomData_generic.__init__(self, K, embs, group)
     582        eps = K(K.unit_group().gens()[1])
     583        if group == 'sl' or (group == 'plus' and eps.norm() == -1):
     584            self._epsmat = Matrix([[eps, 0], [0, eps**-1]])
     585            eps = eps**2
     586        elif group in ['plus', 'gl']:
     587            if self._embs[0](eps) < 0:
     588                eps = -eps
     589            self._epsmat = Matrix([[eps, 0], [0, 1]])
     590        else:
     591            raise ValueError, "Unknown group in _HilbertFundDomData: %s" %group
     592        self._eps = eps
     593        self._epslog = self._RR(self._embs[1](eps).abs().log())
     594
     595    def to_additive(self, x):
     596        """
     597        Move x in RR^g to an additive fundamental domain.
     598
     599        Returns an pair (a, x') with a in O_K and x' = x-a in the fundamental domain.
     600
     601        Examples::
     602
     603            sage: from sage.schemes.hyperelliptic_curves.stoll_cremona import _HilbertFundDomData
     604            sage: h = _HilbertFundDomData(QuadraticField(5))
     605            sage: h.to_additive([3.2,7.9])
     606            (a + 6, [-0.563932022500210, -0.336067977499789])
     607        """
     608        embs = self._embs
     609        K = self._K
     610        om = K.maximal_order().gens()[1]
     611        y = ((x[1]-x[0]) / (embs[1](om) - embs[0](om))).round()
     612        a = y * om
     613        x = [x[0] - y*embs[0](om), x[1] - y*embs[1](om)]
     614        y = ((x[1]+x[0]) / 2).round()
     615        a = a + y
     616        x = [x[0] - y, x[1] - y]
     617        return a, x
     618
     619    def to_multiplicative(self, y):
     620        """
     621        Move y in RR_+^g to a multiplicative fundamental domain.
     622
     623        Returns a pair (M, y') with M a diagonal matrix and M*y = y' and y' in a fundamental domain.
     624
     625        Examples::
     626
     627            sage: from sage.schemes.hyperelliptic_curves.stoll_cremona import _HilbertFundDomData
     628            sage: h = _HilbertFundDomData(QuadraticField(15))
     629            sage: h.to_multiplicative([1,45678])
     630            (
     631            [-a + 4      0]
     632            [     0  a + 4], [61.9838667696635, 736.933695500892]
     633            )
     634            sage: h.act(_[0], [CC.0, CC.0*45678])
     635            [61.9838667696594*I, 736.933695500891*I]
     636            sage: h = _HilbertFundDomData(QuadraticField(15), sl=False)
     637            sage: h.to_multiplicative([1,45678])
     638            (
     639            [-63*a + 244           0]
     640            [          0           1], [487.997950815860, 93.6028520695257]
     641            )
     642            sage: h.act(_[0], [CC.0, CC.0*45678])
     643            [487.997950811067*I, 93.6028520686063*I]
     644        """
     645        embs = self._embs
     646        K = self._K
     647        eps = self._eps
     648        l = self._epslog
     649        e = ((y[1].abs().log() - y[0].abs().log()) / (2*l)).round()
     650        a = eps**e
     651        y = [y[0]/embs[0](a), y[1]/embs[1](a)]
     652        return self._epsmat**-e, y
     653
     654
     655    def _up_cd(self, z, bound_on_c=None):
     656        """
     657        Move z in HH^g, i.e., increase the norm of its imaginary part,
     658        if possible.
     659
     660        Returns either None or a matrix M in SL_2(O_K) that such that
     661        Mz has a larger norm of its imaginary part.
     662
     663        Assumes K has class number one.
     664
     665        bound_on_c None means try a small one first, double each time until a proven one is reached
     666        """
     667        embs = self._embs
     668        cd = self._up_cd_lll(z)
     669        K = self._K
     670        O = K.maximal_order()
     671        if not cd is None:
     672            (c, d) = cd
     673            cz_plus_d = [embs[i](c)*z[i] + embs[i](d) for i in range(2)]
     674            if abs(prod(cz_plus_d)) < 1:
     675                assert c*O+d*O == 1*O
     676                return c, d
     677        # Assuming class number 1 and the above did not work, we are pretty
     678        # high up in the fundamental domain, so the search ranges below are not
     679        # that large.
     680        # TODO: two different notions of reduced are tested, so something
     681        # may be missed here. In that case, we may be high up after all.
     682
     683        om = O.gens()[1]
     684
     685        # M = (a, b; c, d) we need norm(cz+d)<1
     686        # and we want such c and d only up to units.
     687        # We have 1 > norm(|cz+d|) >= |norm(c)|*norm(im(z)),
     688        # so |norm(c)| < norm(im(z))^-1.
     689        # We can find all such c up to units by enumerating ideals
     690        # and finding generators.
     691        # Given c, we want a coprime d with norm(|cx+d|)^2 minimal.
     692        x = [z[0].real(), z[1].real()]
     693        y = [z[0].imag(), z[1].imag()]
     694        norm_im_z = self._RR(abs(prod(y)))
     695        bound = ZZ((norm_im_z**-1).floor())
     696#        if bound_on_c is None:
     697#            bound_to_try = 1
     698#            while bound_to_try < 2*bound:
     699#                m = self._up_cd(z, bound_on_c = bound_to_try)
     700#                if m != None:
     701#                    return m
     702#                bound_to_try = bound_to_try * 2
     703#            return None
     704        if (not bound_on_c is None) and bound > bound_on_c:
     705            bound = bound_on_c
     706        ids = K.ideals_of_bdd_norm(bound)
     707        for n in ids:
     708            for C in ids[n]:
     709                if C.is_principal():
     710                    eps = K(K.unit_group().gens()[1])
     711                    c = C.gens_reduced()[0]
     712                    c_times_y = [embs[0](c)*y[0], embs[1](c)*y[1]]
     713                    e = ((self._RR(c_times_y[0]/c_times_y[1]).abs().log()) /
     714                         (2*(self._RR(embs[0](eps)).abs().log())))
     715#                    print c, e
     716                    c = c / eps ** (e.round())
     717#                    print c
     718                       
     719                    # Next, we try tofind d with |norm(cz+d)|<1.
     720                    # We start by moving to an additive fundamental domain,
     721                    # and do a and thorough search if that does not work
     722                    # (which it usually does not).
     723                    c_times_x = [embs[0](c)*x[0], embs[1](c)*x[1]]
     724                    d0, xp = self.to_additive(c_times_x)
     725                    d0 = -d0
     726                    # now xp = c*x + d_0 is in the fundamental domain
     727                    cz_plus_d = [embs[0](c)*z[0] + embs[0](d0), embs[1](c)*z[1] + embs[1](d0)]
     728                    # Actually, the following check already uses class number != 1, because
     729                    # otherwise the check is supposed to be < Norm(c*O+d0*O) or something.
     730
     731                    if abs(prod(cz_plus_d)) < 1:
     732                        if c*O+d0*O == 1*O:
     733                           return c, d0
     734                        # In the class number > 1 case, we may want to do something
     735                        # if the test c*O+d0*O == 1*O fails.
     736                               
     737                    # We have 1 > norm(|cz+d|^2) =
     738
     739                    #            |norm(c)|^2*norm(im(z))^2
     740                    #          + |c_1|^2*im(z_1)^2 *|c_0x_0+d_0|^2
     741                    #          + |c_0|^2*im(z_0)^2 *|c_1x_1+d_1|^2
     742                    #          +  norm(|cx+d|)^2,
     743                    # define x_0' = c_0x_0+d_0
     744                    # and    x_1' = c_1x_1+d_1, then they have the following
     745                    # absolute value bounds:
     746                    #
     747                    bound_on_x0p = ((1 - n**2 * norm_im_z**2)
     748                                     / (embs[1](c)**2 * y[1]**2)).sqrt()
     749                    bound_on_x1p = ((1 - n**2 * norm_im_z**2)
     750                                     / (embs[0](c)**2 * y[0]**2)).sqrt()
     751                    d0_interval = [-embs[0](c) * x[0] - bound_on_x0p,
     752                                   -embs[0](c) * x[0] + bound_on_x0p]
     753                    d1_interval = [-embs[1](c) * x[1] - bound_on_x1p,
     754                                   -embs[1](c) * x[1] + bound_on_x1p]
     755                    # Now write d = e + f*om = (g + f*delta)/2,
     756                    # where g = 2*e + f*trace(om),
     757                    # and om = (trace(om) + delta) / 2.
     758                    delta = 2*om - om.trace()
     759                    assert embs[0](delta) < 0
     760                    assert embs[1](delta) > 0
     761                    delta_num = embs[1](delta)
     762                    # Then g = d0 + d1 is in the following interval:
     763                    g_interval = [d0_interval[i] + d1_interval[i]
     764                                  for i in [0,1]]
     765                    g_start = sum([-embs[0](c) * x[i] for i in [0,1]]).round()
     766                    g_bound = (bound_on_x0p + bound_on_x1p + 3).floor()
     767                    t = om.trace()
     768                    for i in xrange(g_bound):
     769                        # TODO: this range can be way too big, something needs
     770                        #       to be done about that.
     771                        if t == 0 or not ((g_start + i) % 2):
     772                            for s in ([-1,1] if i else [1]):
     773                                g = g_start + s * i
     774                                # now f = (-g + 2*d)/delta, and we know
     775                                # bounds on each embedding:
     776                                f0_lower = (-g + 2*d0_interval[1]) / -delta_num
     777                                f1_lower = (-g + 2*d1_interval[0]) /  delta_num
     778                                f0_upper = (-g + 2*d0_interval[0]) / -delta_num
     779                                f1_upper = (-g + 2*d1_interval[1]) /  delta_num
     780                                f_lower = max([f0_lower, f1_lower]).ceil()
     781                                f_upper = min([f0_upper, f1_upper]).floor()
     782                                for f in range(f_lower, f_upper+1):
     783                                    if not ((g - f*t) % 2):
     784                                        d = (g+f*delta)/2
     785                                        assert d in O
     786                                        cz_plus_d = [embs[0](c)*z[0] +
     787                                                     embs[0](d),
     788                                                     embs[1](c)*z[1] +
     789                                                     embs[1](d)]
     790                                        if abs(prod(cz_plus_d)) < 1:
     791                                            if c*O+d*O == 1*O:
     792                                                return c, d
     793                                            # In the class number > 1 case,
     794                                            # we may want to do something
     795                                            # if the test c*O+d0*O == 1*O
     796                                            # fails. For now, the code just
     797                                            # skips c,d in that case.
     798                                        #
     799                    # # if we write d = e + f*om, then d0+d1 = 2*e + f*tr(om)
     800                    # # is in the following interval:
     801                    #
     802                    # d0-d1 = f*(om_0 - om_1) is in the following interval
     803                    # interval = [d0_interval[0]-d1_interval[1], d0_interval[1]-d1_interval[0]]
     804                    # f_interval = [interval[i] / (embs[0](om)-embs[1](om)) for i in [0,1]]
     805                    # f_interval.sort()
     806                    # for f in range(int(self._RR(f_interval[0]).ceil()), int(self._RR(f_interval[1]).floor()+1)):
     807                        # e_lower = self._RR(max([d0_interval[0]-f*embs[0](om), d1_interval[0]-f*embs[1](om)])).ceil()
     808                        # e_upper = self._RR(min([d0_interval[1]-f*embs[0](om), d1_interval[1]-f*embs[1](om)])).floor()
     809                        # for e in range(int(e_lower), int(e_upper+1)):
     810                            # d = e + f*om
     811                            # cz_plus_d = [embs[0](c)*z[0] + embs[0](d), embs[1](c)*z[1] + embs[1](d)]
     812                            # if abs(prod(cz_plus_d)) < 1:
     813                                # if c*O+d*O == 1*O:
     814                                    # return c, d
     815                                # # In the class number > 1 case, we may want to do something
     816                                # # if the test c*O+d0*O == 1*O fails.
     817
     818
     819class _HilbertFundDomData_classical(_HilbertFundDomData_generic):
     820    """
     821    As _HilbertFundDomData, but only for the classical upper half
     822    plane and the action of SL_2(ZZ).
     823    """
     824
     825    _eps = None
     826
     827    def __init__(self, embs=None):
     828        """
     829        Create the object.
     830        """
     831        self._K = QQ
     832
     833    def to_additive(self, x):
     834        """
     835        Move x in RR^1 to an additive fundamental domain.
     836        """
     837        a = x[0].round()
     838        x = [x[0]-a]
     839        return a, x
     840
     841    def to_multiplicative(self, y):
     842        """
     843        Returns (identitymatrix(2), y).
     844        """
     845        return Matrix([[1,0],[0,1]]), y
     846
     847
     848    def up(self, z):
     849        """
     850        Move z up in HH, i.e., increase its imaginary part,
     851        if possible by the transformation z -> -1/z
     852
     853        Returns either 1 or a matrix M in SL_2(ZZ) such that
     854        Mz has a larger imaginary part.
     855
     856        Assumes the real part of x is between -1/2 and 1/2
     857        """
     858        if z[0].abs() < 1:
     859            return Matrix(QQ, [[0,-1],[1,0]])
     860        return Matrix(QQ, [[1,0],[0,1]])
     861
     862    def step(self, z):
     863        """
     864        Do one reduction round.
     865
     866        Returns a pair (M, z') with z' = M(z) nearer the fundamental domain if possible.
     867        """
     868        M1 = self.up(z)
     869        z = self.act(M1, z)
     870        t = self.to_additive([e.real() for e in z])[0]
     871        M3 = Matrix([[1, -t], [0, 1]])
     872        z = self.act(M3, z)
     873        return M3*M1, z
     874
     875
     876    def act(self, M, z):
     877        """
     878        Return M(z)
     879        """
     880        a = M[0,0]
     881        b = M[0,1]
     882        c = M[1,0]
     883        d = M[1,1]
     884        z = z[0]
     885        return [(a*z+b)/(c*z+d)]
     886
     887
     888