Ticket #14755: 14755-reduction.patch

File 14755-reduction.patch, 19.6 KB (added by mstreng, 9 years ago)
  • new file sage/schemes/hyperelliptic_curves/reduction.py

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