Ticket #12356: trac_12356-cm-rebase-sage5.0.patch

File trac_12356-cm-rebase-sage5.0.patch, 25.1 KB (added by was, 10 years ago)

Rebased everything as a single patch against sage-5.0.beta2 -- apply only this.

  • sage/schemes/elliptic_curves/all.py

    # HG changeset patch
    # User William Stein <wstein@gmail.com>
    # Date 1327494911 0
    # Node ID cc7b0b365082532447e026f955040fb3113e525f
    # Parent  0acd49e280b7e8d9ac56065c56b025e69117a6d1
    trac 12356 - Refactor cm_j_invariants_and_orders and extend to more general case
    * * *
    trac 12356 (part 2) -- refactored cm_j_invariants_and_orders; fix significant bug pointed out by Cremona
    
    diff --git a/sage/schemes/elliptic_curves/all.py b/sage/schemes/elliptic_curves/all.py
    a b  
    2727
    2828from ell_rational_field import cremona_curves, cremona_optimal_curves
    2929
    30 from cm import ( cm_j_invariants,
     30from cm import ( cm_orders,
     31                 cm_j_invariants,
    3132                 cm_j_invariants_and_orders,
    3233                 hilbert_class_polynomial )
    3334
  • sage/schemes/elliptic_curves/cm.py

    diff --git a/sage/schemes/elliptic_curves/cm.py b/sage/schemes/elliptic_curves/cm.py
    a b  
    77- ``cm_j_invariants``
    88- ``cm_j_invariants_and_orders``
    99
     10AUTHORS:
     11
     12- Robert Bradshaw
     13- John Cremona
     14- William Stein
    1015"""
    1116
    1217#*****************************************************************************
    13 #       Copyright (C) 2005 William Stein <wstein@gmail.com>
     18#       Copyright (C) 2005-2012 William Stein <wstein@gmail.com>
    1419#
    1520#  Distributed under the terms of the GNU General Public License (GPL)
    1621#
     
    144149    return IntegerRing()['x'](coeffs)
    145150
    146151
    147 def cm_j_invariants(K):
     152def cm_j_invariants(K, proof=None):
    148153    r"""
    149154    Return a list of all CM `j`-invariants in the field `K`.
    150155
    151156    INPUT:
    152157   
    153158    - ``K`` -- a number field (currently only implemented for `K` of degree at most 2)
     159    - `proof` -- (default: proof.number_field())
    154160
    155161    OUTPUT:
    156162   
    157163    (list) -- A list of CM `j`-invariants in the field `K`.
    158164
    159     .. note::
    160 
    161        This is currently only implemented for the rationals and
    162        quadratic fields.  David Kohel has large tables for other
    163        fields, but they are not in Sage yet.
    164 
    165165    EXAMPLE::
    166166   
    167167        sage: cm_j_invariants(QQ)
    168         [0, 54000, -12288000, 1728, 287496, -3375, 16581375, 8000, -32768, -884736, -884736000, -147197952000, -262537412640768000]
     168        [-262537412640768000, -147197952000, -884736000, -12288000, -884736, -32768, -3375, 0, 1728, 8000, 54000, 287496, 16581375]
    169169
    170170    Over imaginary quadratic fields there are no more than over `QQ`::
    171171   
    172         sage: cm_j_invariants(QuadraticField(-1, 'i'))
    173         [0, 54000, -12288000, 1728, 287496, -3375, 16581375, 8000, -32768, -884736, -884736000, -147197952000, -262537412640768000]
     172        sage: cm_j_invariants(QuadraticField(-1, 'i'))   # sorted differently since QQ[i] sorts funny
     173        [-12288000, 54000, 0, 287496, 1728, 16581375, -3375, 8000, -32768, -884736, -884736000, -147197952000, -262537412640768000]
     174       
     175    Over real quadratic fields there may be more, for example::
     176   
     177        sage: len(cm_j_invariants(QuadraticField(5, 'a')))
     178        31
     179       
     180    Over number fields K of many higher degrees this also works::
     181   
     182        sage: K.<a> = NumberField(x^3 - 2)
     183        sage: cm_j_invariants(K)
     184        [-12288000, 54000, 0, 287496, 1728, 16581375, -3375, 8000, -32768, -884736, -884736000, -147197952000, -262537412640768000, 31710790944000*a^2 + 39953093016000*a + 50337742902000]
     185        sage: K.<a> = NumberField(x^4 - 2)
     186        sage: len(cm_j_invariants(K))
     187        23
     188    """
     189    return list(sorted([j for D,f,j in cm_j_invariants_and_orders(K, proof=proof)]))
    174190
    175     Over real quadratic fields there may be more::
    176    
    177         sage: cm_j_invariants(QuadraticField(5, 'a')) 
    178         [0, 54000, -12288000, 1728, 287496, -3375, 16581375, 8000, -32768, -884736, -884736000, -147197952000, -262537412640768000, 282880*a + 632000, -282880*a + 632000, 95178240*a + 212846400, -95178240*a + 212846400, 85995/2*a - 191025/2, -85995/2*a - 191025/2, 16554983445/2*a + 37018076625/2, -16554983445/2*a + 37018076625/2, 26378240*a - 58982400, -26378240*a - 58982400, 95673435586560*a - 213932305612800, -95673435586560*a - 213932305612800, 184068066743177379840*a - 411588709724712960000, -184068066743177379840*a - 411588709724712960000]
    179 
    180     Over fields of higher degree this is not yet implemented::
    181    
    182         sage: K.<a> = NumberField(x^3-2)
    183         sage: cm_j_invariants(K)
    184         Traceback (most recent call last):
    185         ...
    186         NotImplementedError: Enumeration of CM j-invariants over fields of degree 3 not yet implemented
    187 
    188     """
    189     if K == QQ:
    190         return [Integer(x) for x in [0, 54000, -12288000, 1728, \
    191                                287496, -3375, 16581375, 8000, \
    192                                -32768,  -884736, -884736000,\
    193                                -147197952000, -262537412640768000]]
    194 
    195     else:
    196         return [j for D,f,j in cm_j_invariants_and_orders(K)]
    197 
    198 def cm_j_invariants_and_orders(K):
     191def cm_j_invariants_and_orders(K, proof=None):
    199192    r"""
    200193    Return a list of all CM `j`-invariants in the field `K`, together with the associated orders.
    201194
    202195    INPUT:
    203196
    204197    - ``K`` -- a number field (currently only implemented for `K` of degree at most 2)
     198    - `proof` -- (default: proof.number_field())
    205199
    206200    OUTPUT:
    207201
     
    212206    EXAMPLE::
    213207   
    214208        sage: cm_j_invariants_and_orders(QQ)
    215         [(-3, 1, 0), (-3, 2, 54000), (-3, 3, -12288000), (-4, 1, 1728), (-4, 2, 287496), (-7, 1, -3375), (-7, 2, 16581375), (-8, 1, 8000), (-11, 1, -32768), (-19, 1, -884736), (-43, 1, -884736000), (-67, 1, -147197952000), (-163, 1, -262537412640768000)]
     209        [(-3, 3, -12288000), (-3, 2, 54000), (-3, 1, 0), (-4, 2, 287496), (-4, 1, 1728), (-7, 2, 16581375), (-7, 1, -3375), (-8, 1, 8000), (-11, 1, -32768), (-19, 1, -884736), (-43, 1, -884736000), (-67, 1, -147197952000), (-163, 1, -262537412640768000)]
    216210
    217211    Over an imaginary quadratic field there are no more than over `QQ`::
    218212
    219         sage: cm_j_invariants_and_orders(QuadraticField(-1, 'i')) 
    220         [(-3, 1, 0), (-3, 2, 54000), (-3, 3, -12288000), (-4, 1, 1728), (-4, 2, 287496), (-7, 1, -3375), (-7, 2, 16581375), (-8, 1, 8000), (-11, 1, -32768), (-19, 1, -884736), (-43, 1, -884736000), (-67, 1, -147197952000), (-163, 1, -262537412640768000)]
     213        sage: cm_j_invariants_and_orders(QuadraticField(-1, 'i'))
     214        [(-3, 3, -12288000), (-3, 2, 54000), (-3, 1, 0), (-4, 2, 287496), (-4, 1, 1728), (-7, 2, 16581375), (-7, 1, -3375), (-8, 1, 8000), (-11, 1, -32768), (-19, 1, -884736), (-43, 1, -884736000), (-67, 1, -147197952000), (-163, 1, -262537412640768000)]
    221215
    222216    Over real quadratic fields there may be more::
    223217       
    224         sage: cm_j_invariants_and_orders(QuadraticField(5,'a'))
    225         [(-3, 1, 0), (-3, 2, 54000), (-3, 3, -12288000), (-4, 1, 1728), (-4, 2, 287496), (-7, 1, -3375), (-7, 2, 16581375), (-8, 1, 8000), (-11, 1, -32768), (-19, 1, -884736), (-43, 1, -884736000), (-67, 1, -147197952000), (-163, 1, -262537412640768000), (-20, 1, 282880*a + 632000), (-20, 1, -282880*a + 632000), (-40, 1, 95178240*a + 212846400), (-40, 1, -95178240*a + 212846400), (-15, 1, 85995/2*a - 191025/2), (-15, 1, -85995/2*a - 191025/2), (-15, 2, 16554983445/2*a + 37018076625/2), (-15, 2, -16554983445/2*a + 37018076625/2), (-35, 1, 26378240*a - 58982400), (-35, 1, -26378240*a - 58982400), (-115, 1, 95673435586560*a - 213932305612800), (-115, 1, -95673435586560*a - 213932305612800), (-235, 1, 184068066743177379840*a - 411588709724712960000), (-235, 1, -184068066743177379840*a - 411588709724712960000)]
     218        sage: v = cm_j_invariants_and_orders(QuadraticField(5,'a')); len(v)
     219        31
     220        sage: [(D,f) for D,f,j in v if j not in QQ]
     221        [(-3, 5), (-3, 5), (-4, 5), (-4, 5), (-15, 2), (-15, 2), (-15, 1), (-15, 1), (-20, 1), (-20, 1), (-35, 1), (-35, 1), (-40, 1), (-40, 1), (-115, 1), (-115, 1), (-235, 1), (-235, 1)]
    226222
    227         sage: [(D,f) for D,f,j in cm_j_invariants_and_orders(QuadraticField(5,'a')) if j not in QQ]
    228         [(-20, 1), (-20, 1), (-40, 1), (-40, 1), (-15, 1), (-15, 1), (-15, 2), (-15, 2), (-35, 1), (-35, 1), (-115, 1), (-115, 1), (-235, 1), (-235, 1)]
     223    Over number fields K of many higher degrees this also works::
     224   
     225        sage: K.<a> = NumberField(x^3 - 2)
     226        sage: cm_j_invariants_and_orders(K)
     227        [(-3, 3, -12288000), (-3, 2, 54000), (-3, 1, 0), (-4, 2, 287496), (-4, 1, 1728), (-7, 2, 16581375), (-7, 1, -3375), (-8, 1, 8000), (-11, 1, -32768), (-19, 1, -884736), (-43, 1, -884736000), (-67, 1, -147197952000), (-163, 1, -262537412640768000), (-3, 6, 31710790944000*a^2 + 39953093016000*a + 50337742902000)]
     228    """
     229    # Get the list of CM orders that could possibly have Hilbert class
     230    # polynomial F(x) with a root in K.  If F(x) has a root alpha in K,
     231    # then F is the minimal polynomial of alpha in K, so the degree of
     232    # F(x) is at most [K:QQ].
     233    dlist = sum([v for h,v in discriminants_with_bounded_class_number(K.degree(), proof=proof).iteritems()], [])
    229234
    230     Over fields of higher degree this is not yet implemented::
     235    return [(D,f,j) for D, f in dlist
     236             for j in hilbert_class_polynomial(D*f*f).roots(K, multiplicities=False)]
     237       
    231238   
    232         sage: K.<a> = NumberField(x^3-2)
    233         sage: cm_j_invariants_and_orders(K)
     239def cm_orders(h, proof=None):
     240    """
     241    Return a list of all pairs `(D,f)` where there is a CM order of
     242    discriminant `D f^2` with class number h, with `D` a fundamental
     243    discriminant.
     244
     245    INPUT:
     246
     247    - `h` -- positive integer
     248    - `proof` -- (default: proof.number_field())
     249
     250    OUTPUT:
     251
     252    - list of 2-tuples `(D,f)`
     253
     254    EXAMPLES::
     255
     256        sage: cm_orders(0)
     257        []
     258        sage: v = cm_orders(1); v
     259        [(-3, 3), (-3, 2), (-3, 1), (-4, 2), (-4, 1), (-7, 2), (-7, 1), (-8, 1), (-11, 1), (-19, 1), (-43, 1), (-67, 1), (-163, 1)]
     260        sage: type(v[0][0]), type(v[0][1])
     261        (<type 'sage.rings.integer.Integer'>, <type 'sage.rings.integer.Integer'>)
     262        sage: v = cm_orders(2); v
     263         [(-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), (-51, 1), (-52, 1), (-88, 1), (-91, 1), (-115, 1), (-123, 1), (-148, 1), (-187, 1), (-232, 1), (-235, 1), (-267, 1), (-403, 1), (-427, 1)]
     264        sage: len(v)
     265        29
     266        sage: set([hilbert_class_polynomial(D*f^2).degree() for D,f in v])
     267        set([2])
     268
     269    Any degree up to 100 is implemented, but may be prohibitively slow::
     270   
     271        sage: cm_orders(3)
     272        [(-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)]
     273        sage: len(cm_orders(4))
     274        84
     275    """
     276    h = Integer(h)
     277    T = None
     278    if h <= 0:
     279        # trivial case
     280        return []
     281    # Get information for all discriminants then throw away everything
     282    # but for h.  If this is replaced by a table it will be faster,
     283    # but not now.   (David Kohel is rumored to have a large table.)
     284    return discriminants_with_bounded_class_number(h, proof=proof)[h]
     285
     286# Table from Mark Watkins paper "Class numbers of imaginary quadratic fields".
     287# I extracted this by cutting/pasting from the pdf, and running this program:
     288# z = {}
     289# for X in open('/Users/wstein/tmp/a.txt').readlines():
     290#    if len(X.strip()):
     291#        v = [int(a) for a in X.split()]
     292#        for i in range(5):
     293#            z[v[3*i]]=(v[3*i+2], v[3*i+1])
     294watkins_table = {1: (163, 9), 2: (427, 18), 3: (907, 16), 4: (1555, 54), 5: (2683, 25),
     295                 6: (3763, 51), 7: (5923, 31), 8: (6307, 131), 9: (10627, 34), 10:
     296                 (13843, 87), 11: (15667, 41), 12: (17803, 206), 13: (20563, 37), 14:
     297                 (30067, 95), 15: (34483, 68), 16: (31243, 322), 17: (37123, 45), 18:
     298                 (48427, 150), 19: (38707, 47), 20: (58507, 350), 21: (61483, 85), 22:
     299                 (85507, 139), 23: (90787, 68), 24: (111763, 511), 25: (93307, 95), 26:
     300                 (103027, 190), 27: (103387, 93), 28: (126043, 457), 29: (166147, 83),
     301                 30: (134467, 255), 31: (133387, 73), 32: (164803, 708), 33: (222643, 101),
     302                 34: (189883, 219), 35: (210907, 103), 36: (217627, 668), 37:
     303                 (158923, 85), 38: (289963, 237), 39: (253507, 115), 40: (260947, 912),
     304                 41: (296587, 109), 42: (280267, 339), 43: (300787, 106), 44: (319867, 691),
     305                 45: (308323, 154), 46: (462883, 268), 47: (375523, 107), 48:
     306                 (335203, 1365), 49: (393187, 132), 50: (389467, 345), 51: (546067, 159),
     307                 52: (439147, 770), 53: (425107, 114), 54: (532123, 427), 55: (452083,163),
     308                 56: (494323, 1205), 57: (615883, 179), 58: (586987, 291),
     309                 59:(474307, 128), 60: (662803, 1302), 61: (606643, 132), 62: (647707, 323),
     310                 63: (991027, 216), 64: (693067, 1672), 65: (703123, 164), 66: (958483, 530),
     311                 67: (652723, 120), 68: (819163, 976), 69: (888427, 209), 70:(811507, 560),
     312                 71: (909547, 150), 72: (947923, 1930), 73: (886867, 119),
     313                 74: (951043, 407), 75: (916507, 237), 76: (1086187, 1075), 77: (1242763, 216),
     314                 78: (1004347, 561), 79: (1333963, 175), 80: (1165483, 2277), 81: (1030723, 228),
     315                 82: (1446547, 402), 83: (1074907, 150), 84: (1225387,1715),
     316                 85: (1285747, 221), 86: (1534723, 472), 87: (1261747, 222),
     317                 88:(1265587, 1905), 89: (1429387, 192), 90: (1548523, 801),
     318                 91: (1391083,214), 92: (1452067, 1248), 93: (1475203, 262), 94: (1587763, 509),
     319                 95:(1659067, 241), 96: (1684027, 3283), 97: (1842523, 185), 98: (2383747,580),
     320                 99: (1480627, 289), 100: (1856563, 1736)}
     321
     322def largest_fundamental_disc_with_class_number(h):
     323    """
     324    Return largest absolute value of any fundamental discriminant with
     325    class number h, and the number of fundamental discriminants with
     326    that class number.  This is known for `h` up to 100, by work of Mark
     327    Watkins.
     328
     329    INPUT:
     330
     331    - `h` -- integer     
     332   
     333    EXAMPLES::
     334
     335        sage: sage.schemes.elliptic_curves.cm.largest_fundamental_disc_with_class_number(0)
     336        (0, 0)
     337        sage: sage.schemes.elliptic_curves.cm.largest_fundamental_disc_with_class_number(1)
     338        (163, 9)
     339        sage: sage.schemes.elliptic_curves.cm.largest_fundamental_disc_with_class_number(2)
     340        (427, 18)
     341        sage: sage.schemes.elliptic_curves.cm.largest_fundamental_disc_with_class_number(10)
     342        (13843, 87)
     343        sage: sage.schemes.elliptic_curves.cm.largest_fundamental_disc_with_class_number(100)
     344        (1856563, 1736)       
     345        sage: sage.schemes.elliptic_curves.cm.largest_fundamental_disc_with_class_number(101)
    234346        Traceback (most recent call last):
    235347        ...
    236         NotImplementedError: Enumeration of CM j-invariants over fields of degree 3 not yet implemented
     348        NotImplementedError: largest discriminant not known for class number 101
     349    """
     350    h = Integer(h)
     351    if h <= 0:
     352        # very easy special case
     353        return Integer(0), Integer(0)
     354    try:
     355        # simply look up the answer in Watkins's table.
     356        B, c = watkins_table[h]
     357        return (Integer(B), Integer(c))
     358    except KeyError:
     359        # nobody knows, since I guess Watkins's is state of the art.
     360        raise NotImplementedError, "largest discriminant not known for class number %s"%h
    237361
     362def discriminants_with_bounded_class_number(hmax, B=None, proof=None):
     363    """
     364    Return dictionary with keys class numbers h<=hmax and values the
     365    list of all pairs (D, f), with D<0 a fundamental discriminant such
     366    that D*f^2 has class number h.  If the optional bound B is given,
     367    return only those pairs with fundamental |D| <= B, though f can
     368    still be arbitrarily large.
    238369
     370    INPUT:
     371    - ``hmax`` -- integer
     372    - `B` -- integer or None; if None returns all pairs
     373    - ``proof -- this code calls the PARI function qfbclassno, so it
     374      could give wrong answers when proof==False.  The default is
     375      whatever proof.number_field() is.  If proof==False and B is
     376      None, at least the number of discriminants is correct, since it
     377      is double checked with Watkins's table.
     378
     379    OUTPUT:
     380    - dictionary
     381
     382    In case B is not given, we use Mark Watkins's: "Class numbers of
     383    imaginary quadratic fields" to compute a B that captures all h
     384    up to hmax.
     385
     386    EXAMPLES::
     387
     388        sage: v = sage.schemes.elliptic_curves.cm.discriminants_with_bounded_class_number(3)
     389        sage: v.keys()
     390        [1, 2, 3]
     391        sage: v[1]
     392        [(-3, 3), (-3, 2), (-3, 1), (-4, 2), (-4, 1), (-7, 2), (-7, 1), (-8, 1), (-11, 1), (-19, 1), (-43, 1), (-67, 1), (-163, 1)]
     393        sage: v[2]
     394        [(-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), (-51, 1), (-52, 1), (-88, 1), (-91, 1), (-115, 1), (-123, 1), (-148, 1), (-187, 1), (-232, 1), (-235, 1), (-267, 1), (-403, 1), (-427, 1)]
     395        sage: v[3]
     396        [(-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)]
     397        sage: v = sage.schemes.elliptic_curves.cm.discriminants_with_bounded_class_number(8, proof=False)
     398        sage: [len(v[h]) for h in v.keys()]
     399        [13, 29, 25, 84, 29, 101, 38, 208]
     400
     401    Find all class numbers for discriminant up to 50::
     402   
     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)]}
    239405    """
    240     if K == QQ:
    241         T = [ (0,-3, 1), (54000,-3,2), (-12288000, -3,3), (1728,-1, 1), \
    242                (287496,-1, 2), (-3375,-7,1), (16581375, -7, 2), (8000,-2,1), \
    243                (-32768, -11, 1),  (-884736, -19,1), (-884736000,-43,1),\
    244                (-147197952000, -67,1), (-262537412640768000,-163,1)
    245                ]
    246         dis = lambda D:  Integer(D) if D%4 in [0,1] else Integer(4*D)
    247         T = [(dis(D),Integer(f),Integer(j)) for (j,D,f) in T]
     406    # imports that are needed only for this function
     407    from sage.structure.proof.proof import get_flag
     408    from sage.libs.pari.gen import pari
     409    import math
     410    from sage.misc.functional import round
     411
     412    # deal with input defaults and type checking
     413    proof = get_flag(proof, 'number_field')
     414    hmax = Integer(hmax)
     415
     416    # T stores the output
     417    T = {}
     418
     419    # Easy case -- instead of giving error, give meaningful output
     420    if hmax < 1:
    248421        return T
    249422
    250     if K.degree()==2:
    251         # Below we use that the imaginary quadratic orders of class number 2 are the
    252         # maximal orders in Q(sqrt(-d)) for d in
    253         # [-5,-6,-10,-13,-15,-22,-35,-37,-51,-58,-91,-115,-123,-187,-235,-267,-403,-427]
    254         # and the order of index 2 in Q(sqrt(-15)). [Reference: many
    255         # places including J E Cremona, Abelian Varieties with Extra
    256         # Twist, Cusp Forms, and Elliptic Curves Over Imaginary
    257         # Quadratic Fields, Journal of the London Mathematical Society
    258         # 45 (1992) 402-416.]     
     423    if B is None:
     424        # Determine how far we have to go by applying Watkins's theorem.
     425        v = [largest_fundamental_disc_with_class_number(h) for h in range(1, hmax+1)]
     426        B = max([b for b,_ in v])
     427        fund_count = [0] + [cnt for _,cnt in v]
     428    else:
     429        # Nothing to do -- set to None so we can use this later to know not
     430        # to do a double check about how many we find.
     431        fund_count = None
     432        B = Integer(B)
    259433
    260         def data_for_discriminant(d):
    261             D = d if d%4 == 1 else 4*d
    262             for j in hilbert_class_polynomial(D).roots(K,multiplicities=False):
    263                 yield (Integer(D), Integer(1), j)
    264             # the only non-maximal order of class number 2 is for d == -15
    265             if d == -15:
    266                 for j in hilbert_class_polynomial(-60).roots(K,multiplicities=False):
    267                     yield (Integer(-15), Integer(2), j)
     434    if B <= 2:
     435        # This is an easy special case, since there are no fundamental discriminants
     436        # this small.
     437        return T
    268438
    269         # list the discriminants of class number 2:
    270         dlist = [-5,-6,-10,-13,-15,-22,-35,-37,-51,-58,-91,-115,-123,-187,-235,-267,-403,-427]
    271         return sum((list(data_for_discriminant(d)) for d in dlist), cm_j_invariants_and_orders(QQ))
     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)))
    272446
    273     raise NotImplementedError, "Enumeration of CM j-invariants over fields of degree %s not yet implemented"%K.degree()
    274     # If you want to implement this in general, see the remarks at trac 11220.
     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
     459    for D in range(-B, -2):
     460        if is_fundamental_discriminant(D):
     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)
     496                # If the class number of this order is within the range, then
     497                # use it.  (NOTE: In some cases there is a simple relation between
     498                # the class number for D and D*f^2, and this could be used to
     499                # optimize this inner loop a little.)
     500                if h <= hmax:
     501                    z = (Integer(D), Integer(f))
     502                    if T.has_key(h):
     503                        T[h].append(z)
     504                    else:
     505                        T[h] = [z]
     506                f += 1
     507
     508    for h in T.keys():
     509        T[h] = list(reversed(T[h]))
     510
     511    if fund_count is not None:
     512        # Double check that we found the right number of fundamental
     513        # discriminants; we might as well, since Watkins provides this
     514        # data.
     515        for h in T.keys():
     516            if len([D for D,f in T[h] if f==1]) != fund_count[h]:
     517                raise RuntimeError, "number of discriminants inconsistent with Watkins's table"
     518       
     519    return T
     520                   
     521
     522
     523