Ticket #14755: 14755-reduction-cleaner.patch

File 14755-reduction-cleaner.patch, 25.5 KB (added by florian, 9 years ago)
  • sage/schemes/hyperelliptic_curves/mestre.py

    # HG changeset patch
    # User Florian Bouyer <f.j.s.c.bouyer@gmail.com>
    # Date 1374763285 -3600
    # Node ID 7b9ec6250ec1712ff4036f6d25951f03e694ee3c
    # Parent  c53db919bae2ec678ceb8d356fb222a493f9fe80
    14755 reduction patch
    
    diff --git a/sage/schemes/hyperelliptic_curves/mestre.py b/sage/schemes/hyperelliptic_curves/mestre.py
    a b  
    2828from sage.schemes.plane_conics.constructor import Conic
    2929from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
    3030from sage.schemes.hyperelliptic_curves.constructor import HyperellipticCurve
     31from sage.schemes.hyperelliptic_curves.reduction import reduce_hyperelliptic_curve_polynomial
    3132
    3233
    3334def HyperellipticCurve_from_invariants(i, reduced=True, precision=None,
     
    6667
    6768        sage: HyperellipticCurve_from_invariants([3840,414720,491028480,2437709561856])
    6869        Traceback (most recent call last):
    69         ...
    70         NotImplementedError: Reduction of hyperelliptic curves not yet implemented. See trac #14755 and #14756.
     70        Hyperelliptic Curve over Rational Fields defined by y^2 = -x^6 - 3*x + 1
    7171        sage: HyperellipticCurve_from_invariants([3840,414720,491028480,2437709561856],reduced = False)
    7272        Hyperelliptic Curve over Rational Field defined by y^2 = -x^6 + 4410*x^5 - 540*x^4 + 4320*x^3 - 19440*x^2 + 46656*x - 46656
    7373        sage: HyperellipticCurve_from_invariants([21, 225/64, 22941/512, 1])
    74         Traceback (most recent call last):
    75         ...
    76         NotImplementedError: Reduction of hyperelliptic curves not yet implemented. See trac #14755 and #14756.
     74        Hyperelliptic Curve over Rational Field defined by y^2 = 908*x^6 - 7156*x^5 + 69567*x^4 - 208564*x^3 - 3663441*x^2 - 10179516*x - 59567060
    7775
    7876    An example over a finite field::
    7977
     
    209207        pass
    210208
    211209    if reduced:
    212         raise NotImplementedError("Reduction of hyperelliptic curves not " \
    213                                    "yet implemented. " \
    214                                    "See trac #14755 and #14756.")
     210        f = reduce_hyperelliptic_curve_polynomial(f)
    215211
    216212    return HyperellipticCurve(f)
    217213
  • new file sage/schemes/hyperelliptic_curves/reduction.py

    diff --git a/sage/schemes/hyperelliptic_curves/reduction.py b/sage/schemes/hyperelliptic_curves/reduction.py
    new file mode 100644
    - +  
     1r"""
     2Reduction of binary forms and hyperelliptic curve equations.
     3
     4This file contains functions that convert hyperelliptic curve equations into
     5equivalent ones with hopefully smaller coefficients.
     6
     7The first step is to reduce discriminants of hyperelliptic curve equations at
     8non-archimedean primes of a number field. Doing this globally is generally
     9impossible when the class number is > 1, so this is only approximate.
     10
     11This is followed by reduction as in Stoll-Cremona [SC2003]_, which
     12reduces the archimedean size of the coefficients. See
     13:module:`sage.schemes.hyperelliptic_curves.stoll_cremona`.
     14
     15For details see Bouyer-Streng [BS2013]_
     16
     17REFERENCES:
     18.. [SC2003] M. Stoll, J. Cremona, "On the reduction theory of binary forms"
     19   Journal fur die reine und angewandter Mathematik 565, pages 79 - 99
     20.. [BS2013] F. Bouyer, M. Steng, "Examples of CM curves of genus two defined
     21   over the reflex field", arXiv:1307.0486[math.NT]
     22
     23AUTHORS:
     24
     25- Florian Bouyer
     26- Marco Streng
     27
     28"""
     29#*****************************************************************************
     30#       Copyright (C) 2011,2012,2013
     31#                 Florian Bouyer <f.j.s.c.bouyer@gmail.com>
     32#                 Marco Streng <marco.streng@gmail.com>
     33#
     34#  Distributed under the terms of the GNU General Public License (GPL)
     35#  as published by the Free Software Foundation; either version 2 of
     36#  the License, or (at your option) any later version.
     37#                  http://www.gnu.org/licenses/
     38#*****************************************************************************
     39
     40from sage.rings.arith import gcd
     41from sage.rings.all import ZZ, QQ, RealField, FiniteField
     42from sage.misc.misc import prod
     43from sage.schemes.hyperelliptic_curves.stoll_cremona import stoll_cremona_reduction
     44
     45
     46def reduce_hyperelliptic_curve_polynomial(f, precision=3000, algorithm='default', primes=None):   
     47    r"""
     48
     49    Reduces f by using the three functions reduce_gcd, reduce_discriminant,
     50    stoll_cremona_reduction.
     51    (See their respective documentations for details)
     52   
     53    Currently this is an internal function used by
     54    :func: `HyperellipticCurve_from_invariants`. In the future it will
     55    be a method of the HyperellipticCurve class.
     56   
     57    INPUT:
     58       
     59    - ``f`` - a polynomial
     60   
     61    - ``precision`` - an integer (default: 3000) the bit precision of the
     62    numbers in `\CC` and `\RR` that is going to be used.
     63   
     64    - ``algorithm`` - default/magma. If set to magma in the stoll_cremona
     65    reduction step, the function uses magma to calculate a 'better' covariant
     66    of f
     67
     68    - ``primes`` (default:None) -- a list of primes. If supplied, then reduce
     69    the discriminant only at the given primes
     70
     71    OUTPUT:
     72       
     73    A polynomial
     74   
     75    EXAMPLES:
     76   
     77    Over `\QQ`::
     78
     79        sage: P.<x> = QQ[]
     80        sage: from sage.schemes.hyperelliptic_curves.reduction import reduce_hyperelliptic_curve_polynomial
     81        sage: f = -16*x^6 - 1611*x^5 - 8640*x^4 + 69120*x^3 - 311040*x^2 + 746496*x - 746496
     82        sage: reduce_hyperelliptic_curve_polynomial(f)
     83        -x^6 + 3*x - 2
     84        sage: f =  -5*x^6 + 174*x^5 - 2700*x^4 + 21600*x^3 - 97200*x^2 + 233280*x - 233280
     85        sage: reduce_hyperelliptic_curve_polynomial(f)
     86        -x^6 - x^5 - 5
     87
     88    Over number fields::
     89   
     90        sage: X = var('X')
     91        sage: k = NumberField(X^2 - 41,'a')
     92        sage: a = k.an_element()
     93        sage: P.<x> = k[]
     94        sage: f =  x^6 + 2*x^4 + (64*a + 410)*x^3 + 1311360*a + 8396801
     95        sage: reduce_hyperelliptic_curve_polynomial(f)
     96        (-320*a + 2049)*x^6 + (-640*a + 4098)*x^4 + (-64*a + 410)*x^3 + 320*a + 2049
     97        sage: k = NumberField(X^2 - 73,'a')
     98        sage: a = k.an_element()
     99        sage: P.<x> = k[]
     100        sage: f = (1/778034*a - 1/778034)*x^6 + (3/5329*a - 3/5329)*x^5 + (15/146*a - 9/146)*x^4 + (10*a + 2)*x^3 + (1095/2*a + 1533/2)*x^2 + (15987*a + 47961)*x + 389017/2*a + 2723119/2
     101        sage: reduce_hyperelliptic_curve_polynomial(f)
     102        x^6 + 3*x^2 + 1/2*a - 1/2
     103     
     104    """
     105    if precision is None:
     106        precision = 3000
     107    k = f(0).parent()
     108    f = reduce_gcd(f)
     109    f = reduce_discriminant(f, primes=primes)
     110    f = stoll_cremona_reduction(f, algorithm=algorithm, precision=precision)
     111    try:
     112        f = reduce_unit_in_disc(f)
     113    except NotImplementedError:
     114        pass
     115    return f
     116
     117
     118def reduce_gcd(f):   
     119    r"""
     120    Takes out the common factors of the coefficients of f. In the case of
     121    number fields of class number larger than one, there may still be a small
     122    prime ideal in the denominator after the reduction.
     123
     124    This is an internal function called by :func: `reduce_hyperelliptic_curve_polynomial`
     125    INPUT:
     126   
     127    - ``f`` - a polynomial
     128   
     129    OUTPUT:
     130 
     131    If the base ring is ZZ or QQ, then the output is an integer polynomial
     132    whose coefficients have gcd 1 and which is a rational scaling of f.
     133   
     134    If the base ring is a number field, then it still removes many prime factors
     135    from the gcd of the coefficients (all in the case of class number one).
     136   
     137    EXAMPLES::
     138
     139        sage: P.<x> = QQ[]
     140        sage: from sage.schemes.hyperelliptic_curves.reduction import reduce_gcd
     141        sage: reduce_gcd(3*x^2+6*x+321)
     142        x^2 + 2*x + 107
     143        sage: reduce_gcd(321*x^6 + 1284*x^2 - 2247*x + 5136)
     144        x^6 + 4*x^2 - 7*x + 16
     145       
     146        sage: X = var('X')
     147        sage: k = NumberField(X^2 + 41,'a')
     148        sage: a = k.an_element()
     149        sage: P.<x> = k[]
     150        sage: reduce_gcd((a + 41)*x^6 + (41*a + 41)*x^4 + (42*a + 82)*x^2 + 16*a + 656)
     151        (-a + 1)*x^6 + (-a + 41)*x^4 + (-2*a + 42)*x^2 - 16*a + 16
     152       
     153    """
     154    k = f.base_ring()
     155    f = f.denominator()*f
     156   
     157    if k is QQ or k is ZZ:
     158        gcdofcoeffs = gcd([ZZ(coeff) for coeff in f.coeffs()])
     159        return f/gcdofcoeffs
     160   
     161    A = k.ideal(f.list())
     162    return f / gen_almost(A)
     163
     164
     165def gen_almost(A):
     166    r"""
     167    If A is a principal ideal, returns a generator.
     168   
     169    Otherwise, if A is a number field ideal, returns an element x of the number
     170    field such that x divides A and (A/x).norm() is small.
     171   
     172    This is an internal function, used by :func:`reduce_at_infinity`
     173    and :func:`reduce_gcd`.
     174   
     175    EXAMPLES::
     176   
     177        sage: k = QuadraticField(-5,'a')
     178        sage: p = k.ideal(next_prime(2000)).factor()[0][0]; p
     179        Fractional ideal (2003, a + 1483)
     180        sage: from sage.schemes.hyperelliptic_curves.reduction import gen_almost
     181        sage: x = gen_almost(p); x
     182        -27/2*a + 19/2
     183        sage: p/x
     184        Fractional ideal (2, a + 1)
     185
     186    """
     187    if A.is_principal():
     188        return A.gens_reduced()[0]
     189    # Here is the obvious way:
     190    #B = A.reduce_equiv()
     191    #C = A/B
     192    #if not C.is_principal():
     193    #    raise RuntimeError
     194   
     195    # So we do something that looks stupid and slow,
     196    # but seemed to be faster:
     197    k = A.number_field()
     198    for c in k.class_group():
     199        B = c.ideal()
     200        C = A/B
     201        if C.is_principal():
     202            break
     203
     204    if not B.is_integral():
     205        raise RuntimeError
     206    return C.gens_reduced()[0]
     207
     208
     209def reduce_at_infinity(f, primes=None):
     210    r"""
     211    If there is a 4-fold root at infinity modulo a prime, then try to remove
     212    it.
     213   
     214    This works over any number field. If the class number is one, then the
     215    discriminant of the output divides the discriminant of the input.
     216   
     217    This is an internal function, called by :func:`reduce_discriminant`.
     218   
     219    INPUT:
     220       
     221    - ``f`` -- a polynomial of degree 5 or 6
     222    - ``primes`` (default:None) -- a list of primes,
     223      reduce only at the given primes
     224   
     225    OUTPUT:
     226   
     227    A polynomial ``g`` such that as many prime powers
     228    as possible are removed from the discriminant by
     229    transformations of the form `f(x/u^k)u^l` with
     230    `k,l>=0` and `u` a uniformizer of a prime.
     231
     232    If the class number of the base field is one,
     233    then a small prime may remain or be introduced.
     234   
     235    EXAMPLES::
     236   
     237        sage: P.<x> = QQ[]
     238        sage: h = (prod([x+3^7*5^8*n for n in [1,2,3,4]])*(134*x^2+384*x+378)).reverse(); h
     239        4832126692668678984045982360839843750000*x^6 + 4908827128145711155933141708374023437500*x^5 + 1712976141118022078847492961120605468750*x^4 + 4177353776971384564386468750000*x^3 + 3422870579757550781628*x^2 + 1144757812884*x + 134
     240        sage: log(h.discriminant()).n()
     241        643.624696259543 + 3.14159265358979*I
     242        sage: from sage.schemes.hyperelliptic_curves.reduction import reduce_at_infinity
     243        sage: log(reduce_at_infinity(h).discriminant()).n()
     244        405.008434398004 + 3.14159265358979*I
     245
     246        sage: X = var('X')
     247        sage: k = NumberField(X^2 + 41,'a')
     248        sage: a = k.an_element()
     249        sage: P.<x> = k[]
     250        sage: h = (prod([x+3^7*5^8*n for n in [1,2,3,4]])*(134*x^2+384*x+378)).reverse(); h
     251        4832126692668678984045982360839843750000*x^6 + 4908827128145711155933141708374023437500*x^5 + 1712976141118022078847492961120605468750*x^4 + 4177353776971384564386468750000*x^3 + 3422870579757550781628*x^2 + 1144757812884*x + 134
     252        sage: log(h.discriminant().norm()).n()
     253        1287.24939251909
     254        sage: log(reduce_at_infinity(h).discriminant().norm()).n()
     255        796.153925184809
     256
     257    """
     258    k = f.base_ring()
     259    x = k['x'].gen()
     260   
     261    if k is ZZ or k is QQ:
     262        A = gcd([ZZ(c) for c in f.list()[4:]])
     263        Bz = 1
     264        Cz = 1
     265    else:
     266        A = k.ideal(f.list()[4:])
     267        Bz = k.ideal(1)
     268        Cz = k.ideal(1)
     269    By = 1
     270    Cy = 1
     271    if primes is None:
     272        primes = [p for (p,e) in A.factor()]
     273    for p in primes:
     274        maxpow = ZZ(min([(f[6].valuation(p)-1)/3, \
     275                         (f[5].valuation(p)-1)/2, \
     276                          f[4].valuation(p)-1, \
     277                          f[3].valuation(p)]).floor())
     278        if maxpow > 0:
     279            powers = max([[3*pow + \
     280                         min([f[indx].valuation(p)-indx*pow \
     281                             for indx in range(7)]),pow] \
     282                         for pow in range(maxpow+1)])
     283            if powers[0] > 0:
     284                # The following two operations (taking powers of p) take way
     285                # too long in Sage < 5.7 when p is an ideal, hence we
     286                # use different powering functions, making use of
     287                # the class group.
     288                #B = B * p**powers[1]
     289                #C = C * p**(3*powers[1]-powers[0])
     290                y, z = power_of_ideal(p, powers[1], factors=True)
     291                Bz = Bz * z
     292                By = By * y
     293                y, z = power_of_ideal(p, 3*powers[1]-powers[0], factors=True)
     294                Cz = Cz * z
     295                Cy = Cy * y
     296                # If we knew that B and C were principal, then
     297                # the reduction step would be
     298                # g = f(x/p**powers[1]) * p**(3*powers[1]-powers[0])
     299
     300    if k is ZZ or k is QQ:
     301        return f(x/By/Bz)*Cy*Cz
     302   
     303    f = f(x/By/gen_almost(Bz)) * Cy*gen_almost(Cz)
     304    return reduce_gcd(f)
     305
     306
     307def power_of_ideal(A, e, factors=False):
     308    r"""
     309    Return A**e. Uses the fact that the class group
     310    is small, so goes via principal ideals.
     311   
     312    If factors is True, returns an element
     313    and an ideal, such that A**e is the product.
     314
     315    This is an internal function called by :func: `reduce_at_infinity`
     316    """
     317    if A in QQ:
     318        if factors:
     319            return A**e, 1
     320        return A**e
     321
     322    k = A.number_field()
     323    m = k.class_group()(A).order()
     324    q, r = ZZ(e).quo_rem(m)
     325    assert (A**m).is_principal()
     326    g = (A**m).gens_reduced()[0]
     327    x = g**q
     328    B = A**r
     329    if factors:
     330        return (x, B)
     331    return x*B
     332
     333
     334def reduce_discriminant(f, primes=None):   
     335    r"""
     336    Removes 10th powers from the discriminant by a `GL_2(QQ)` transformation.
     337
     338    This is an internal function called by :func: `reduce_hyperelliptic_curve_polynomial`
     339   
     340    INPUT:
     341       
     342    - ``f`` - a univariate polynomial of degree at most 6, interpreted as the
     343      binary sextic form `F(X,Y) = Y^6 f(X/Y)`. The polynomial ``f`` is not
     344      allowed to have repeated roots.
     345   
     346    OUTPUT:
     347       
     348    a polynomial giving a `GL_2(QQ)` equivalent binary form with 10th powers
     349    removed from the discriminant where possible
     350   
     351    EXAMPLES::
     352
     353        sage: from sage.schemes.hyperelliptic_curves.reduction import reduce_discriminant
     354        sage: P.<x> = QQ[]
     355        sage: f = 32*x^6 + 6*x^2 - 5
     356        sage: factor(f.discriminant())
     357        2^33 * 3^6 * 5 * 13^2
     358        sage: g = reduce_discriminant(f); g
     359        x^6 + 3*x^2 - 10
     360        sage: factor(g.discriminant())
     361        2^13 * 3^6 * 5 * 13^2
     362
     363    An example where the degree of the input is 5::
     364
     365        sage: P.<y> = QQ[]
     366        sage: f = y^5 + 3^20*y + 3^20
     367        sage: g = f(y+1)
     368        sage: f
     369        y^5 + 3486784401*y + 3486784401
     370        sage: reduce_discriminant(f)
     371        x^5 + 81*x + 1
     372        sage: g
     373        y^5 + 5*y^4 + 10*y^3 + 10*y^2 + 3486784406*y + 6973568803
     374        sage: reduce_discriminant(g)
     375        x^5 + 5*x^4 + 10*x^3 + 10*x^2 + 86*x + 83
     376       
     377    An example where a large power is removed::
     378   
     379        sage: x = QQ['x'].gen()
     380        sage: f = ((x-2)*(x-2+3^20)*(x-2+3^10)*(x-2+3^20*7)*(x^2+x+1))
     381        sage: reduce_discriminant(f)
     382        3486784401*x^6 + 1647132543836838*x^5 + 85105305481926930091*x^4 + 85110864542496497113*x^3 + 7206360479055108*x^2 + 170852435649*x
     383        sage: reduce_discriminant(f).discriminant().factor()
     384        -1 * 2^10 * 3^83 * 7^6 * 11^4 * 61^2 * 233^2 * 887^2 * 1549321^2 * 4607887^2 * 11920789^2 * 1964677987^2 * 3486489163^2 * 6188121169^2
     385        sage: f.discriminant().factor()
     386        -1 * 2^10 * 3^183 * 7^6 * 11^4 * 61^2 * 233^2 * 887^2 * 1549321^2 * 4607887^2 * 11920789^2 * 1964677987^2 * 3486489163^2 * 6188121169^2
     387        sage: reduce_discriminant(f.reverse()).discriminant() == reduce_discriminant(f).discriminant()
     388        True
     389
     390    ALGORITHM:
     391       
     392    First we create a list of primes that are factors of the discriminant D of
     393    f. For each prime `p` check that ``valuation(D, p)>=10``
     394    Then we try to maximize `e = min\{v(a_n)+(3-n)k\}` with 0 <= k <=
     395    `min\{(v(a_{n})-1)/3,v(a_{n-1})-1)/2,v(a_{n-2})-1,v(a_{n-3})\}
     396    (where `v(a_n)` is ``valuation(a_n, p)``).
     397    Finally we do the substitution `f = f(x/p^e) and clear denominator
     398    Next we do the same trick but looking at the coefficients in the reverse
     399    order. Note that this comes down to seeing whether f has 0 as a fourfold
     400    root mod `p`. So first we improve a bit by moving any fourfold root to 0
     401    before applying the same trick. (but replacing `a_{n-i}` with `a_i` in the
     402    condition k <= min...).
     403   
     404    """
     405    k = f(0).parent()
     406    OK = k.OK()
     407    P = k['x']
     408    x = P.gen()
     409    f = reduce_gcd(f)
     410    f = reduce_at_infinity(f, primes=primes)
     411    one = ZZ(1)
     412
     413    from sage.schemes.hyperelliptic_curves.invariants import igusa_clebsch_invariants
     414
     415    if k is QQ or k is ZZ:
     416        D = gcd(igusa_clebsch_invariants(f))
     417   
     418    if not (k is QQ or k is ZZ):
     419        D = k.ideal(igusa_clebsch_invariants(f))
     420        c = coprime_representatives(k, D)
     421
     422    if primes is None:
     423        primes = [prime for (prime, power) in D.factor()]
     424    for prime in primes:
     425        if (D.valuation(prime)>=2):
     426            if k is QQ or k is ZZ:
     427                u_unif = prime
     428                l_unif = prime
     429            else:
     430                u_unif = almost_generator(prime, c)
     431                l_unif = one/almost_generator(one/prime, c)
     432            f = reduce_locally(f, prime, u_unif, l_unif)
     433    f = reduce_gcd(f)
     434    f = reduce_at_infinity(f)
     435    return f
     436
     437
     438def reduce_locally(f, prime, u_unif=None, l_unif=None):
     439    r"""
     440    Reduces `C : y^2 = f(x)` locally at ``prime``. Assumes that `C` has no
     441    automorphisms over the algebraic closure other than the identity and the
     442    hyperelliptic involution.
     443   
     444    This is an internal function called by :func: `reduce_discriminant`
     445    INPUT:
     446       
     447    - ``f`` - a univariate polynomial of degree 2g+1 or 2g+2, interpreted as
     448      the binary sextic form `F(X,Y) = Y^(2g+2) f(X/Y)`. The polynomial ``f``
     449      is not allowed to have repeated roots.
     450    - ``prime`` - a prime number
     451    - ``u_unif`` - a uniformizer such that u_unif / prime is integral
     452    - ``l_unif`` - a uniformizer such that prime / l_unif is integral
     453   
     454    OUTPUT:
     455       
     456    a polynomial giving a `GL_2(QQ)` equivalent binary form with this prime
     457    removed from the discriminant where possible
     458    """
     459    g = ((ZZ(f.degree())-1)/2).floor()
     460   
     461    k = f.base_ring()
     462   
     463    if k is QQ:
     464        F = FiniteField(prime)
     465        if u_unif is None:
     466            u_unif = QQ(prime)
     467        if l_unif is None:
     468            l_unif = QQ(prime)
     469    else:
     470        F = prime.residue_field()
     471        if prime.is_principal():
     472            if u_unif is None:
     473                u_unif = prime.gens_reduced()[0]
     474            if l_unif is None:
     475                l_unif = prime.gens_reduced()[0]
     476        elif u_unif is None or l_unif is None:
     477            raise ValueError, "Please supply u_unif and l_unif for non-principal ideal"
     478
     479    x = f.parent().gen()
     480
     481    try_this_prime = True
     482    while try_this_prime:
     483        try_this_prime = False
     484        m = min([c.valuation(prime) for c in f.list()])
     485        f = f / l_unif**m
     486        for k in range(g+2):
     487            if f[2*g+2-k].valuation(prime) < g+2-k:
     488                break
     489        else:
     490            # now val(f[2*g+2-k])>=g+2-k for all k
     491            f = f(x/l_unif)
     492            m = -min([c.valuation(prime) for c in f.list()])
     493            assert m <= g
     494            f = f * u_unif ** m
     495       
     496        rts = f.roots(F)
     497        for ind in range(len(rts)):
     498            if rts[ind][1] >= 4:
     499                root_mod_prime = F.lift(rts[ind][0])
     500                h = f(x + root_mod_prime)
     501                for k in range(g+2):
     502                    if h[k].valuation(prime) < g+2-k:
     503                        break
     504                else:
     505                    # now val(h[k])>=g+2-k for all k
     506                    f = h(x*u_unif)
     507                    try_this_prime = True
     508                   
     509    return f
     510
     511
     512def coprime_generator_representatives(k, A):
     513    r"""
     514    Given a number field k and an ideal A, returns a list of integral ideals,
     515    one for every generator of the class group, all of them coprime to A.
     516   
     517    This is an auxiliary function of :func:`coprime_representatives` and only
     518    used for number fields of class number >1.
     519    """
     520    c = k.class_group()
     521    gens = c.gens()
     522    gens_todo = [(gens[i], i) for i in range(len(gens))]
     523    primes = [None for g in gens]
     524    for p in k.primes_of_degree_one_iter():
     525        if p+A == 1:
     526            for a in gens_todo:
     527                if c(p) == a[0]:
     528                    primes[a[1]] = p
     529                    gens_todo.remove(a)
     530                    if len(gens_todo) == 0:
     531                        return primes
     532                    break
     533
     534
     535def coprime_representatives(k, A):
     536    r"""
     537    Given a number field k and an ideal A, returns a list of integral ideals,
     538    one for every ideal class, all of them coprime to A.
     539   
     540    This provides data used by :func:`almost_generator` and only
     541    used for number fields of class number >1.
     542    """
     543    c = k.class_group()
     544    l = coprime_generator_representatives(k, A)
     545    ret = []
     546    for a in c:
     547        m = a.list()
     548        ret.append(prod([l[i]**m[i] for i in range(len(m))]))
     549    return ret
     550
     551
     552def almost_generator(I, coprime_rep):
     553    r"""
     554    Given an ideal I and the output of :func:`coprime_representatives`, returns
     555    a generator of I times an element of ``coprime_rep``.
     556   
     557    This is only used for number fields of class number >1.
     558    """
     559    for J in coprime_rep:
     560        K = I*J
     561        if K.is_principal():
     562            return K.gens_reduced()[0]
     563
     564
     565def reduce_unit_in_disc(f, precision=100):
     566    """
     567    Reduce the size of the coefficients of f by multiplying f by a unit.
     568   
     569    INPUT:
     570       
     571    -``f``- a polynomial in a number field `K`
     572   
     573    OUTPUT:
     574   
     575    a polynomial of the form u*f, where u is a unit of the maximal order of `K`
     576   
     577    EXAMPLES::
     578   
     579        sage: from sage.schemes.hyperelliptic_curves.reduction import reduce_unit_in_disc
     580        sage: X = var('X')
     581        sage: k = NumberField(X^2-41,'a')
     582        sage: a = k.an_element()
     583        sage: P.<x> = k[]
     584        sage: f = 2*x^4 + (64*a + 410)*x^3 + 1311360*a + 8396801
     585        sage: factor(f.discriminant())
     586        (-369836393348730718080*a - 2368108374136006451201) * (-118*a - 725) * (-118*a + 725) * (-1/2*a + 7/2)^4 * (1/2*a + 7/2)^4
     587        sage: g = reduce_unit_in_disc(f); g
     588        (-640*a + 4098)*x^4 + (-64*a + 410)*x^3 + 320*a + 2049
     589        sage: factor(g.discriminant())
     590        (-1) * (-118*a - 725) * (-118*a + 725) * (-1/2*a + 7/2)^4 * (1/2*a + 7/2)^4
     591        sage: f = (-320*a + 2049)*x^6 + (-10374*a + 66426)*x^5 + (-74076*a + 474318)*x^4 + (228552*a - 1463448)*x^3 + (1510568*a - 9672351)*x^2 + (-4678674*a + 29958126)*x + 3024324*a - 19365120
     592        sage: factor(f.discriminant())
     593        (-6210885195014008862151805440*a + 39769069528107045198367948801) * (-a - 6)^4 * (-2*a + 1)^2 * (-a)^6 * (-a + 2)^2 * (4*a - 27)^2 * (3*a - 20)^2 * (-1/2*a + 7/2)^26 * (a - 6)^6 * (-a + 8)^2 * 3^2 * (1/2*a + 7/2)^22
     594        sage: g = reduce_unit_in_disc(f); g
     595        x^6 + (-6*a - 6)*x^5 + (36*a + 462)*x^4 + (-312*a - 2712)*x^3 + (1512*a + 4961)*x^2 + (-2706*a - 2706)*x + 1476*a
     596        sage: factor(g.discriminant())
     597        (1311360*a + 8396801) * (-a - 6)^4 * (-2*a + 1)^2 * (-a)^6 * (-a + 2)^2 * (4*a - 27)^2 * (3*a - 20)^2 * (-1/2*a + 7/2)^26 * (a - 6)^6 * (-a + 8)^2 * 3^2 * (1/2*a + 7/2)^22
     598
     599    Note that this does nothing over `\QQ`::
     600   
     601        sage: P.<x> = QQ[]
     602        sage: f = x^2 + 3
     603        sage: reduce_unit_in_disc(f)
     604        x^2 + 3
     605
     606    This is not implemented for number fields of degree greater than 2::
     607
     608        sage: P.<y> = NumberField(x^3+2,'a')[]
     609        sage: reduce_unit_in_disc(y^6+y+1)
     610        Traceback (most recent call last):
     611        ...
     612        NotImplementedError
     613    """
     614    k = f.base_ring()
     615    if k is QQ or k.degree() == 1 or (k.degree() == 2 and k.is_totally_imaginary()):
     616        return f
     617    if k.degree() > 2:
     618        raise NotImplementedError
     619    epsilon = k(k.unit_group().gens()[1])
     620    Phi = k.embeddings(RealField(precision))
     621    d = f.discriminant()
     622#    e = abs(Phi[0](epsilon))**30
     623#    while abs(Phi[0](d)/Phi[1](d)) < e:
     624#        f = f*epsilon
     625#        d = f.discriminant()
     626#    e = e**-1
     627#    while abs(Phi[0](d)/Phi[1](d)) > e:
     628#        f = f/epsilon
     629#        d = f.discriminant()
     630#    return f
     631
     632    h = f*epsilon
     633    while max([max([abs(Phi[0](coeff)),abs(Phi[1](coeff))]) for coeff in (h).coefficients()])\
     634          < max([max([abs(Phi[0](coeff)),abs(Phi[1](coeff))]) for coeff in (f).coefficients()]):
     635        f = h
     636        h = f*epsilon
     637    h = f/epsilon
     638    while max([max([abs(Phi[0](coeff)),abs(Phi[1](coeff))]) for coeff in (h).coefficients()])\
     639           < max([max([abs(Phi[0](coeff)),abs(Phi[1](coeff))]) for coeff in (f).coefficients()]):
     640        f = h
     641        h = f/epsilon
     642    return f
     643 No newline at end of file