Ticket #12356: trac_12356-extended-part2.patch

File trac_12356-extended-part2.patch, 6.4 KB (added by was, 10 years ago)

replacing part 2 with a refreshed patch that has a better reference.

  • sage/schemes/elliptic_curves/cm.py

    # HG changeset patch
    # User William Stein <wstein@gmail.com>
    # Date 1327646756 28800
    # Node ID 010c3b81471e685bc6de1437bc04a56dcc04955f
    # Parent  da68bdfdc1f9b1b6e401b49731d4849a86f4d9bc
    trac 12356 (part 2) -- refactored cm_j_invariants_and_orders; fix significant bug pointed out by Cremona
    
    diff --git a/sage/schemes/elliptic_curves/cm.py b/sage/schemes/elliptic_curves/cm.py
    a b  
    364364    Return dictionary with keys class numbers h<=hmax and values the
    365365    list of all pairs (D, f), with D<0 a fundamental discriminant such
    366366    that D*f^2 has class number h.  If the optional bound B is given,
    367     return only those pairs with |D*f^2| <= B.
     367    return only those pairs with fundamental |D| <= B, though f can
     368    still be arbitrarily large.
    368369
    369370    INPUT:
    370371    - ``hmax`` -- integer
     
    395396        [(-3, 9), (-3, 6), (-11, 2), (-19, 2), (-23, 2), (-23, 1), (-31, 2), (-31, 1), (-43, 2), (-59, 1), (-67, 2), (-83, 1), (-107, 1), (-139, 1), (-163, 2), (-211, 1), (-283, 1), (-307, 1), (-331, 1), (-379, 1), (-499, 1), (-547, 1), (-643, 1), (-883, 1), (-907, 1)]
    396397        sage: v = sage.schemes.elliptic_curves.cm.discriminants_with_bounded_class_number(8, proof=False)
    397398        sage: [len(v[h]) for h in v.keys()]
    398         [13, 29, 25, 84, 29, 101, 38, 207]
     399        [13, 29, 25, 84, 29, 101, 38, 208]
    399400
    400401    Find all class numbers for discriminant up to 50::
    401402   
    402         sage: sage.schemes.elliptic_curves.cm.discriminants_with_bounded_class_number(hmax=10^10, B=50)
    403         {1: [(-3, 3), (-3, 2), (-3, 1), (-4, 2), (-4, 1), (-7, 2), (-7, 1), (-8, 1), (-11, 1), (-19, 1), (-43, 1)], 2: [(-3, 4), (-4, 3), (-8, 2), (-15, 2), (-15, 1), (-20, 1), (-24, 1), (-35, 1), (-40, 1)], 3: [(-11, 2), (-23, 1), (-31, 1)], 4: [(-7, 3), (-39, 1)], 5: [(-47, 1)]}   
     403        sage: sage.schemes.elliptic_curves.cm.discriminants_with_bounded_class_number(hmax=5, B=50)
     404        {1: [(-3, 3), (-3, 2), (-3, 1), (-4, 2), (-4, 1), (-7, 2), (-7, 1), (-8, 1), (-11, 1), (-19, 1), (-43, 1)], 2: [(-3, 7), (-3, 5), (-3, 4), (-4, 5), (-4, 4), (-4, 3), (-7, 4), (-8, 3), (-8, 2), (-11, 3), (-15, 2), (-15, 1), (-20, 1), (-24, 1), (-35, 1), (-40, 1)], 3: [(-3, 9), (-3, 6), (-11, 2), (-19, 2), (-23, 2), (-23, 1), (-31, 2), (-31, 1), (-43, 2)], 4: [(-3, 13), (-3, 11), (-3, 8), (-4, 10), (-4, 8), (-4, 7), (-4, 6), (-7, 8), (-7, 6), (-7, 3), (-8, 6), (-8, 4), (-11, 5), (-15, 4), (-19, 5), (-19, 3), (-20, 3), (-20, 2), (-24, 2), (-35, 3), (-39, 2), (-39, 1), (-40, 2), (-43, 3)], 5: [(-47, 2), (-47, 1)]}
    404405    """
    405406    # imports that are needed only for this function
    406407    from sage.structure.proof.proof import get_flag
     
    435436        # this small.
    436437        return T
    437438
     439    # This lower bound gets used in an inner loop below.
     440    from math import log
     441    def lb(f):
     442        """Lower bound on euler_phi."""
     443        # 1.79 > e^gamma = 1.7810724...
     444        if f <= 1: return 0  # don't do log(log(1)) = log(0)
     445        return f/(1.79*log(log(f)) + 3.0/log(log(f)))
     446
     447    # We define a little function to compute the class number of
     448    # discriminant d quickly, using pari, with proof inherited from
     449    # the containing scope.  Fast classno functionality should
     450    # probably be moved elsewhere in Sage, but I'm not sure where.
     451    # The same line also occurs in the number fields code, but to use
     452    # it requires making a number field, which is very slow, and this
     453    # function must be fast, since it is the main bottleneck.
     454    def classno(d):
     455        """Return the class number of the order of discriminant d."""
     456        # There is currently no qfbclassno method in gen.pyx, hence the string.
     457        return Integer(pari('qfbclassno(%s,%s)'%(d, 1 if proof else 0)))
     458
    438459    for D in range(-B, -2):
    439460        if is_fundamental_discriminant(D):
    440             # For each fundamental discrimant loop through the f's up to the bound.
    441             # Note that class_number(D*f^2)>=class_nubmer(D) by a
    442             # standard exact sequence, so we can use Watkins's bound in this more
    443             # general setting.
    444             for f in range(1, round(math.sqrt(B//(-D)))+1):
    445                 h = Integer(pari('qfbclassno(%s,%s)'%(D*f*f, 1 if proof else 0)))
     461            h_D = classno(D)
     462            # For each fundamental discrimant D, loop through the f's such
     463            # that h(D*f^2) could possibly be <= hmax.  As explained to me by Cremona,
     464            # we have h(D*f^2) >= h(D)*phi_D(f) >= h(D)*euler_phi(f), where
     465            # phi_D(f) is like euler_phi(f) but the factor (1-1/p) is replaced
     466            # by a factor of (1-kr(D,p)*p), where kr(D/p) is the Kronecker symbol.
     467            # Since (1-1/p) <= 1 and (1-1/p) <= (1+1/p), we see that
     468            #     euler_phi(f) <= phi_D(f).
     469            #
     470            # We have the following analytic lower bound on euler_phi:
     471            #
     472            #     euler_phi(f) >= lb(f) = f / (exp(euler_gamma)*log(log(n)) + 3/log(log(n))).
     473            #
     474            # See Theorem 8 of Peter Clark's
     475            #   http://math.uga.edu/~pete/4400arithmeticorders.pdf
     476            # which is a consequence of Theorem 15 of
     477            # [Rosser and Schoenfeld, 1962].
     478            #
     479            # By Calculus, we see that the lb(f) is an increasing function of f >= 2.
     480            #
     481            # NOTE: You can visibly "see" that it is a lower bound in Sage with
     482            #   lb(n) = n/(exp(euler_gamma)*log(log(n)) + 3/log(log(n)))
     483            #   plot(lb, (n, 1, 10^4), color='red') + plot(lambda x: euler_phi(int(x)), 1, 10^4).show()
     484            #
     485            # So we consider f=1,2,..., until the first f with lb(f)*h_D > h_max.
     486            # (Note that lb(f) is <= 0 for f=1,2, so nothing special is needed there.)
     487            #
     488            # TODO: Maybe we could do better using a bound for for phi_D(f).
     489            #
     490            f = 1
     491            while lb(f)*h_D <= hmax:
     492                if f == 1:
     493                    h = h_D
     494                else:
     495                    h = classno(D*f*f)
    446496                # If the class number of this order is within the range, then
    447497                # use it.  (NOTE: In some cases there is a simple relation between
    448498                # the class number for D and D*f^2, and this could be used to
     
    453503                        T[h].append(z)
    454504                    else:
    455505                        T[h] = [z]
     506                f += 1
    456507
    457508    for h in T.keys():
    458509        T[h] = list(reversed(T[h]))