From cedd2de47e2ad3712745de40f0f6b686783433d7 Mon Sep 17 00:00:00 2001 From: John Cremona Date: Fri, 17 Feb 2023 10:02:50 +0000 Subject: [PATCH 01/12] Added Cremona-Sutherland algorithm for recognising Hilbert Class Polynomials, making elliptic curve CM testing much faster --- src/doc/en/reference/references/index.rst | 15 + src/sage/quadratic_forms/binary_qf.py | 6 +- src/sage/quadratic_forms/special_values.py | 2 +- src/sage/rings/integer.pyx | 67 +++++ src/sage/rings/number_field/number_field.py | 11 +- .../number_field/number_field_element.pyx | 8 - .../number_field_element_quadratic.pyx | 17 +- src/sage/schemes/elliptic_curves/cm.py | 281 +++++++++++++++--- .../elliptic_curves/ell_finite_field.py | 159 ++++++++++ src/sage/schemes/elliptic_curves/heegner.py | 21 +- 10 files changed, 505 insertions(+), 82 deletions(-) diff --git a/src/doc/en/reference/references/index.rst b/src/doc/en/reference/references/index.rst index 5a7eb0a163c..c4159f90bb2 100644 --- a/src/doc/en/reference/references/index.rst +++ b/src/doc/en/reference/references/index.rst @@ -1539,6 +1539,10 @@ REFERENCES: and Monographs*. American Mathematical Society, Providence, RI, 1999. +.. [Cox1989] David A. Cox. + Primes of the form `x^2+ny^2`. + Wiley, 1989. + .. [CK2008] Derek G. Corneil and Richard M. Krueger, *A Unified View of Graph Searching*, SIAM Jounal on Discrete Mathematics, 22(4), 1259–-1276, 2008. @@ -1635,6 +1639,11 @@ REFERENCES: .. [CrNa2020] \J.E. Cremona and F. Najman, `\QQ`-curves over odd degree number fields, :arxiv:`2004.10054`. +.. [CreSuth2023] \J.E. Cremona and A.V. Sutherland. + *Computing the endomorphism ring of an elliptic curve + over a number field*. + :arxiv:`2301.11169`. + .. [CoCo1] J.H. Conway, H.S.M. Coxeter *Triangulated polygons and frieze patterns*, The Mathematical Gazette (1973) 57 p.87-94 @@ -5760,6 +5769,12 @@ REFERENCES: Mathematics of Computation **80** (2011), pp. 477-500. :arxiv:`0809.3413v3`. +.. [RouSuthZur2022] Jeremy Rouse, Andrew V. Sutherland, David Zureick-Brown. + *`\ell`-adic images of Galois for elliptic curves over `\Q`* (and an appendix with John Voight). + Forum of Mathematics, Sigma , Volume 10 , 2022. + :doi: `10.1017/fms.2022.38`. + :arxiv:`2106.11141`. + .. [SV1970] \H. Schneider and M. Vidyasagar. Cross-positive matrices. SIAM Journal on Numerical Analysis, 7:508-519, 1970. diff --git a/src/sage/quadratic_forms/binary_qf.py b/src/sage/quadratic_forms/binary_qf.py index e7cd7086f09..31bf51e40b3 100755 --- a/src/sage/quadratic_forms/binary_qf.py +++ b/src/sage/quadratic_forms/binary_qf.py @@ -52,7 +52,7 @@ from functools import total_ordering from sage.libs.pari.all import pari_gen -from sage.rings.all import ZZ, is_fundamental_discriminant +from sage.rings.all import ZZ from sage.arith.all import gcd from sage.structure.sage_object import SageObject from sage.matrix.matrix_space import MatrixSpace @@ -606,7 +606,7 @@ def has_fundamental_discriminant(self): sage: Q.has_fundamental_discriminant() False """ - return is_fundamental_discriminant(self.discriminant()) + return self.discriminant().is_fundamental_discriminant() def is_primitive(self): r""" @@ -1773,7 +1773,7 @@ def BinaryQF_reduced_representatives(D, primitive_only=False, proper=True): # For a fundamental discriminant all forms are primitive so we need not check: if primitive_only: - primitive_only = not is_fundamental_discriminant(D) + primitive_only = not D.is_fundamental_discriminant() form_list = [] diff --git a/src/sage/quadratic_forms/special_values.py b/src/sage/quadratic_forms/special_values.py index af80989032e..4e32d276b51 100644 --- a/src/sage/quadratic_forms/special_values.py +++ b/src/sage/quadratic_forms/special_values.py @@ -158,7 +158,7 @@ def QuadraticBernoulliNumber(k, d): Let us create a list of some odd negative fundamental discriminants:: - sage: test_set = [d for d in range(-163, -3, 4) if is_fundamental_discriminant(d)] + sage: test_set = [d for d in srange(-163, -3, 4) if d.is_fundamental_discriminant()] In general, we have `B_{1, \chi_d} = -2 h/w` for odd negative fundamental discriminants:: diff --git a/src/sage/rings/integer.pyx b/src/sage/rings/integer.pyx index c43319942fb..aa7168df78c 100644 --- a/src/sage/rings/integer.pyx +++ b/src/sage/rings/integer.pyx @@ -6011,6 +6011,73 @@ cdef class Integer(sage.structure.element.EuclideanDomainElement): """ return self.__pari__().issquarefree() + def is_discriminant(self): + """ + Returns True if this integer is a discriminant. + + .. NOTE:: + + A discriminant is an integer congruent to 0 or 1 modulo 4. + + EXAMPLES:: + + sage: (-1).is_discriminant() + False + sage: (-4).is_discriminant() + True + sage: 100.is_discriminant() + True + sage: 101.is_discriminant() + True + + TESTS:: + + sage: 0.is_discriminant() + True + sage: 1.is_discriminant() + True + sage: len([D for D in srange(-100,100) if D.is_discriminant()]) + 100 + """ + return self%4 in [0,1] + + def is_fundamental_discriminant(self): + """ + Returns True if this integer is a fundamental_discriminant. + + .. NOTE:: + + A fundamental discriminant is a discrimimant, not 0 or 1 and not a square multiple of a smaller discriminant. + + EXAMPLES:: + + sage: (-4).is_fundamental_discriminant() + True + sage: (-12).is_fundamental_discriminant() + False + sage: 101.is_fundamental_discriminant() + True + + TESTS:: + + sage: 0.is_fundamental_discriminant() + False + sage: 1.is_fundamental_discriminant() + False + sage: len([D for D in srange(-100,100) if D.is_fundamental_discriminant()]) + 61 + + """ + if self in [0,1]: + return False + Dmod4 = self%4 + if Dmod4 in [2,3]: + return False + if Dmod4 == 1: + return self.is_squarefree() + d = self//4 + return d%4 in [2,3] and d.is_squarefree() + cpdef __pari__(self): """ Returns the PARI version of this integer. diff --git a/src/sage/rings/number_field/number_field.py b/src/sage/rings/number_field/number_field.py index 1dd1bcbc973..e8c256cfd92 100644 --- a/src/sage/rings/number_field/number_field.py +++ b/src/sage/rings/number_field/number_field.py @@ -12188,16 +12188,15 @@ def is_fundamental_discriminant(D): EXAMPLES:: sage: [D for D in range(-15,15) if is_fundamental_discriminant(D)] + ... + DeprecationWarning: is_fundamental_discriminant(D) is deprecated; please use D.is_fundamental_discriminant() + ... [-15, -11, -8, -7, -4, -3, 5, 8, 12, 13] sage: [D for D in range(-15,15) if not is_square(D) and QuadraticField(D,'a').disc() == D] [-15, -11, -8, -7, -4, -3, 5, 8, 12, 13] """ - d = D % 4 - if d not in [0, 1]: - return False - return D != 1 and D != 0 and \ - (arith.is_squarefree(D) or - (d == 0 and (D // 4) % 4 in [2, 3] and arith.is_squarefree(D // 4))) + deprecation(35147, "is_fundamental_discriminant(D) is deprecated; please use D.is_fundamental_discriminant()") + return Integer(D).is_fundamental_discriminant() ################### diff --git a/src/sage/rings/number_field/number_field_element.pyx b/src/sage/rings/number_field/number_field_element.pyx index bbe5b90fb34..b5daa2b2d8d 100644 --- a/src/sage/rings/number_field/number_field_element.pyx +++ b/src/sage/rings/number_field/number_field_element.pyx @@ -4642,10 +4642,6 @@ cdef class NumberFieldElement_absolute(NumberFieldElement): sage: a.absolute_charpoly(algorithm='pari') == a.absolute_charpoly(algorithm='sage') True """ - # this hack is necessary because quadratic fields override - # charpoly(), and they don't take the argument 'algorithm' - if algorithm is None: - return self.charpoly(var) return self.charpoly(var, algorithm) def absolute_minpoly(self, var='x', algorithm=None): @@ -4673,10 +4669,6 @@ cdef class NumberFieldElement_absolute(NumberFieldElement): sage: b.absolute_minpoly(algorithm='pari') == b.absolute_minpoly(algorithm='sage') True """ - # this hack is necessary because quadratic fields override - # minpoly(), and they don't take the argument 'algorithm' - if algorithm is None: - return self.minpoly(var) return self.minpoly(var, algorithm) def charpoly(self, var='x', algorithm=None): diff --git a/src/sage/rings/number_field/number_field_element_quadratic.pyx b/src/sage/rings/number_field/number_field_element_quadratic.pyx index ed662749e35..c3ffbe42e0a 100644 --- a/src/sage/rings/number_field/number_field_element_quadratic.pyx +++ b/src/sage/rings/number_field/number_field_element_quadratic.pyx @@ -2129,10 +2129,16 @@ cdef class NumberFieldElement_quadratic(NumberFieldElement_absolute): # D = 1 mod 4 return mpz_fdiv_ui(self.D.value, 4) == 1 - def charpoly(self, var='x'): + def charpoly(self, var='x', algorithm=None): r""" The characteristic polynomial of this element over `\QQ`. + INPUT: + + - ``var`` -- the minimal polynomial is defined over a polynomial ring + in a variable with this name. If not specified this defaults to ``x`` + - ``algorithm`` -- for compatibility with general number field elements: and ignored + EXAMPLES:: sage: K. = NumberField(x^2-x+13) @@ -2147,15 +2153,16 @@ cdef class NumberFieldElement_quadratic(NumberFieldElement_absolute): R = QQ[var] return R([self.norm(), -self.trace(), 1]) - def minpoly(self, var='x'): + def minpoly(self, var='x', algorithm=None): r""" The minimal polynomial of this element over `\QQ`. INPUT: - - ``var`` -- the minimal polynomial is defined over a polynomial ring - in a variable with this name. If not specified this defaults to - ``x``. + - ``var`` -- the minimal polynomial is defined over a polynomial ring + in a variable with this name. If not specified this defaults to ``x`` + - ``algorithm`` -- for compatibility with general number field elements: and ignored + EXAMPLES:: diff --git a/src/sage/schemes/elliptic_curves/cm.py b/src/sage/schemes/elliptic_curves/cm.py index 6d56371ce93..bbb43bac1d1 100644 --- a/src/sage/schemes/elliptic_curves/cm.py +++ b/src/sage/schemes/elliptic_curves/cm.py @@ -4,11 +4,13 @@ This module implements the functions - ``hilbert_class_polynomial`` +- ``is_HCP`` - ``cm_j_invariants`` - ``cm_orders`` - ``discriminants_with_bounded_class_number`` - ``cm_j_invariants_and_orders`` - ``largest_fundamental_disc_with_class_number`` +- ``is_cm_j_invariant`` AUTHORS: @@ -38,9 +40,10 @@ QQ, ZZ, IntegerRing, - is_fundamental_discriminant, - PolynomialRing) + GF, + RR) +from sage.schemes.elliptic_curves.all import EllipticCurve from sage.misc.cachefunc import cached_function @cached_function @@ -171,6 +174,156 @@ def hilbert_class_polynomial(D, algorithm=None): return IntegerRing()['x'](coeffs) +def is_HCP(f, check_monic_irreducible=True): + r""" + Return ``True``, `D` if ``f`` is the Hilbert Class Polynomial `H_D`, else ``False``, 0. + + INPUT: + + - ``f`` -- a polynomial in `\ZZ[X]`. + - ``check_monic_irreducible`` (boolean, default ``True``) -- if + ``True``, check that ``f`` is a monic, irreducible, integer + polynomial. + + OUTPUT: + + (integer) -- `D` if ``f`` is the Hilbert Class Polynomial (HCP) `H_D`, esle 0. + + ALGORITHM: + + Cremona and Sutherland: Algorithm 2 of [CreSuth2023]_. + + EXAMPLES:: + + Even for large degrees this is fast. We test the largest + discriminant of class number 100, for which the HCP has coefficients + with thousands of digits:: + + sage: from sage.schemes.elliptic_curves.cm import is_HCP + sage: D = -1856563 + sage: D.class_number() + 100 + sage: H = hilbert_class_polynomial(D) + sage: H.degree() + 100 + sage: max(c for c in H).ndigits() + 2774 + sage: is_HCP(H) + (True, -1856563) + + Testing polynomials which are not HCPs is faster:: + + sage: is_HCP(H+1) + (False, 0) + + + TESTS:: + + sage: from sage.schemes.elliptic_curves.cm import is_HCP + sage: all(is_HCP(hilbert_class_polynomial(D))==(True,D) for D in srange(-4,-100,-1) if D.is_discriminant()) + True + sage: all(is_HCP(hilbert_class_polynomial(D)+1)==(False,0) for D in srange(-4,-100,-1) if D.is_discriminant()) + True + + """ + zero = ZZ(0) + # optional check that input is monic and irreducible + if check_monic_irreducible: + try: + if not (all(c in ZZ for c in f) and f.is_monic()): + return (False, zero) + f = f.change_ring(ZZ) + except AttributeError: + return (False, zero) + + h = f.degree() + h2list = [d for d in h.divisors() if (d-h)%2==0 and d.prime_to_m_part(2)==1] + pmin = 33 * (h**2 * (RR(h+2).log().log()+2)**2).ceil() + # Guarantees 4*p > |D| for fundamental D under GRH + p = pmin-1 + n = 0 + from sage.arith.all import next_prime + while True: + p = next_prime(p) + n += 1 + fp = f.change_ring(GF(p)) + # Compute X^p-X mod fp manually, avoiding quotient ring which is slower + r = zpow = z = fp.parent().gen() + m = p>>1 + while m: + zpow = (zpow**2) % fp + if m & 1: + r = (zpow * r) % fp + m >>= 1 + # now r = X^p mod fp + d = (r-z).gcd(fp).degree() # number of roots mod p + if d==0: + continue + if not fp.is_squarefree(): + continue + if d = NumberField(H) + sage: is_cm_j_invariant(j) + (True, (-1856563, 1)) + TESTS:: sage: from sage.schemes.elliptic_curves.cm import is_cm_j_invariant sage: all(is_cm_j_invariant(j) == (True, (d,f)) for d,f,j in cm_j_invariants_and_orders(QQ)) True + """ # First we check that j is an algebraic number: from sage.rings.all import NumberFieldElement, NumberField @@ -641,32 +812,48 @@ def is_cm_j_invariant(j, method='new'): if j in QQ: return False, None - # Now j has degree at least 2. If it is not integral so is not CM: + # Next we find its minimal polynomial of j: - if not j.is_integral(): + if j.parent().absolute_degree()==2: + jpol = j.absolute_minpoly() # no algorithm parameter + else: + jpol = j.absolute_minpoly(algorithm='pari') + + # If it does not have integer coefficients then j is not integral, hence not CM: + + if not all(c in ZZ for c in jpol): return False, None - # Next we find its minimal polynomial and degree h, and if h is - # less than the degree of j.parent() we recreate j as an element - # of Q(j): + # Otherwise test whether it is a Hilbert Class Polynomial + # (using the fact that we know that it is monic and irreducible): - jpol = PolynomialRing(QQ,'x')([-j,1]) if j in QQ else j.absolute_minpoly() - h = jpol.degree() + if method == 'CremonaSutherland': + flag, D = is_HCP(jpol, check_monic_irreducible=False) + if flag: + D0 = D.squarefree_part() + if D0%4 !=1: + D0 *= 4 + f = ZZ(D//D0).isqrt() + return (True, (D0,f)) + else: + return (False, None) - # This will be used as a fall-back if we cannot determine the - # result using local data. For this to be necessary there would - # have to be very few primes of degree 1 and norm under 1000, - # since we only need to find one prime of degree 1, good - # reduction for which a_P is nonzero. - if method=='old': + h = jpol.degree() + if method in ['exhaustive', 'old']: if h>100: raise NotImplementedError("CM data only available for class numbers up to 100") for d,f in cm_orders(h): if jpol == hilbert_class_polynomial(d*f**2): - return True, (d,f) - return False, None + return (True, (d,f)) + return (False, None) + + if method not in ['reduction', 'new']: + raise ValueError("Invalid method {} in is_cm_j_invariant".format(method)) + + # Now we use the reduction method - # replace j by a clone whose parent is Q(j), if necessary: + # If the degree h is less than the degree of j.parent() we recreate j as an element + # of Q(j, and replace j by a clone whose parent is Q(j), if necessary: K = j.parent() if h < K.absolute_degree(): @@ -676,7 +863,6 @@ def is_cm_j_invariant(j, method='new'): # Construct an elliptic curve with j-invariant j, with # integral model: - from sage.schemes.elliptic_curves.all import EllipticCurve E = EllipticCurve(j=j).integral_model() D = E.discriminant() prime_bound = 1000 # test primes of degree 1 up to this norm @@ -718,12 +904,12 @@ def is_cm_j_invariant(j, method='new'): cmd = dP cmf = fP elif cmd != dP: # inconsistent with previous - return False, None + return (False, None) else: # consistent d, so update f cmf = cmf.gcd(fP) - if cmd==0: # no conclusion, we found no degree 1 primes, revert to old method - return is_cm_j_invariant(j, method='old') + if cmd==0: # no conclusion, we found no degree 1 primes, revert to default method + return is_cm_j_invariant(j, method='CremonaSutherland') # it looks like cm by disc cmd * f**2 where f divides cmf @@ -733,11 +919,10 @@ def is_cm_j_invariant(j, method='new'): # Now we must check if h(cmd*f**2)==h for f|cmf; if so we check # whether j is a root of the associated Hilbert class polynomial. + h0 = cmd.class_number() for f in cmf.divisors(): # only positive divisors - d = cmd*f**2 - if h != d.class_number(): + if h != OrderClassNumber(cmd,h0,f): continue - pol = hilbert_class_polynomial(d) - if pol(j) == 0: - return True, (cmd, f) - return False, None + if jpol == hilbert_class_polynomial(cmd*f**2): + return (True, (cmd, f)) + return (False, None) diff --git a/src/sage/schemes/elliptic_curves/ell_finite_field.py b/src/sage/schemes/elliptic_curves/ell_finite_field.py index 419c08cf9ca..dc21b35985a 100644 --- a/src/sage/schemes/elliptic_curves/ell_finite_field.py +++ b/src/sage/schemes/elliptic_curves/ell_finite_field.py @@ -694,6 +694,19 @@ def frobenius_endomorphism(self): """ return self.frobenius_isogeny(self.base_field().degree()) + def frobenius_discriminant(self): + r""" + Return the discriminant of the ring `\ZZ[\pi_E]` where `\pi_E` is the Frobenius endomorphism. + + EXAMPLES:: + + sage: F. = GF(11^4) + sage: E = EllipticCurve([t,t]) + sage: E.frobenius_discriminant() + -57339 + """ + return self.frobenius_polynomial().discriminant() + def cardinality_pari(self): r""" Return the cardinality of ``self`` using PARI. @@ -1353,6 +1366,150 @@ def _fetch_cached_order(self, other): if n is not None: self._order = n + def height_above_floor(self, ell, e): + r""" + Return the height of the `j`-invariant of this ordinary elliptic curve on its `\ell`-volcano. + + INPUT: + + - ``ell`` -- a prime number + - ``e`` -- a non-negative integer + + + .. NOTE:: + + For an ordinary `E/\F_q`, and a prime `\ell`, the height + of the `\ell`-volcano containing `j(E)` is the `\ell`-adic + valuation of the conductor the order generated by the + Frobenius `\pi_E`); the height of `j(E)` on its + ell-volcano is the `\ell`-adic valuation of the conductor + of the order `\End(E)`. + + ALGORITHM: + + See [RouSuthZur2022]_. + + EXAMPLES:: + + sage: F = GF(312401) + sage: E = EllipticCurve(F,(0, 0, 0, 309381, 93465)) + sage: D = E.frobenius_discriminant(); D + -687104 + sage: D.factor() + -1 * 2^10 * 11 * 61 + sage: E.height_above_floor(2,8) + 5 + + """ + if self.is_supersingular(): + raise ValueError("{} is not ordinary".format(self)) + if e == 0: + return 0 + j = self.j_invariant() + if j in [0, 1728]: + return e + F = j.parent() + x = polygen(F) + from sage.rings.polynomial.polynomial_ring import polygens + from sage.libs.pari.convert_sage import gen_to_sage + from sage.libs.pari.all import pari + X,Y = polygens(F,['X', 'Y'],2) + phi = gen_to_sage(pari.polmodular(ell),{'x':X, 'y':Y}) + j1 = phi([x,j]).roots(multiplicities=False) + nj1 = len(j1) + on_floor = self.two_torsion_rank() < 2 if ell==2 else nj1 <= ell + if on_floor: + return 0 + if e == 1 or nj1 != ell+1: # double roots can only happen at the surface + return e + if nj1 < 3: + return 0 + j0 = [j,j,j] + h = 1 + while True: + for i in range(3): + r = (phi([x,j1[i]])//(x-j0[i])).roots(multiplicities=False) + if not r: + return h + j0[i] = j1[i] + j1[i] = r[0] + h += 1 + + def endomorphism_discriminant_from_class_number(self, h): + r""" + Return the endomorphism order discriminant of this ordinary elliptic curve, given its class number ``h``. + + INPUT: + + - ``h`` -- a positive integer + + OUTPUT: + + (integer) The discriminant of the endomorphism ring `\End(E)`, if + this has class number ``h``. If `\End(E)` does not have class + number ``h`` the output is undefined. + + ALGORITHM: + + Compute the trace of Frobenius and hence the discriminant `D_0` + and class number `h_0` of the maximal order containing the + endomorphism order. From the given value of `h`, which must be a + multiple of `h_0`, compute the possible conductors, using + :meth:`height_above_floor` for each prime `\ell` dividing the + quotient `h/h_0`. If exactly one conductor `f` remains, return + `f^2D_0`, otherwise return `0`. + + .. NOTE:: + + Adapted from [RouSuthZur2022]_. The application for which + one knows the class number in advance is in the + recognition of Hilbert Class Polynomials: see + :func:`sage.schemes.elliptic_curves.cm.is_HCP`. + + EXAMPLES:: + + sage: F = GF(312401) + sage: E = EllipticCurve(F,(0, 0, 0, 309381, 93465)) + sage: E.endomorphism_discriminant_from_class_number(30) + -671 + + We check that this is the correct discriminant, and the input value of `h` was correct:: + + sage: H = hilbert_class_polynomial(-671) + sage: H(E.j_invariant()) == 0 and H.degree()==30 + True + + """ + F = self.base_field() + if not F.is_finite(): + raise ValueError("Base field {} must be finite".format(F)) + if self.is_supersingular(): + raise ValueError("Elliptic curve ({}) must be ordinary".format(self)) + D1 = self.frobenius_discriminant() + D0 = D1.squarefree_part() + if D0%4 !=1: + D0 *= 4 + v = ZZ(D1//D0).isqrt() + h0 = D0.class_number() + if h%h0: + return 0 + from sage.schemes.elliptic_curves.cm import OrderClassNumber + cs = [v//f for f in v.divisors() if OrderClassNumber(D0,h0,f) == h] # cofactors c=v/f compatible with h(f**2D0)=h + if not cs: + return 0 + if len(cs) == 1: + return (v//cs[0])**2 * D0 + from sage.sets.set import Set + L = sorted(list(Set(sum([c.prime_factors() for c in cs], [])))) + for ell in L: + e = self.height_above_floor(ell,v.valuation(ell)) + cs = [c for c in cs if c.valuation(ell) == e] + if not cs: + return 0 + if len(cs) == 1: + return (v//cs[0])**2 * D0 + return 0 + # dict to hold precomputed coefficient vectors of supersingular j values (excluding 0, 1728): @@ -1641,3 +1798,5 @@ def is_j_supersingular(j, proof=True): # expensive since it involves counting the number of points on E): return E.trace_of_frobenius() % p == 0 + + diff --git a/src/sage/schemes/elliptic_curves/heegner.py b/src/sage/schemes/elliptic_curves/heegner.py index 298f44fedc6..499a27f9685 100644 --- a/src/sage/schemes/elliptic_curves/heegner.py +++ b/src/sage/schemes/elliptic_curves/heegner.py @@ -5699,9 +5699,9 @@ def red(P, ell): def best_heegner_D(ell_1, ell_2): # return the first Heegner D satisfy all hypothesis such that # both ell_1 and ell_2 are inert - D = -5 + D = ZZ(-5) while True: - if number_field.is_fundamental_discriminant(D) and \ + if D.is_fundamental_discriminant() and \ D % ell_1 and D % ell_2 and \ E.satisfies_heegner_hypothesis(D) and \ is_inert(D, ell_1) and is_inert(D, ell_2) and \ @@ -6070,10 +6070,9 @@ def class_number(D): ... ValueError: D (=-5) must be a fundamental discriminant """ - if not number_field.is_fundamental_discriminant(D): + if not D.is_fundamental_discriminant(): raise ValueError("D (=%s) must be a fundamental discriminant" % D) - return QuadraticField(D, 'a').class_number() - + return D.class_number() def is_inert(D, p): r""" @@ -6221,7 +6220,7 @@ def satisfies_weak_heegner_hypothesis(N, D): sage: EllipticCurve('37a').heegner_discriminants_list(10) [-7, -11, -40, -47, -67, -71, -83, -84, -95, -104] """ - if not number_field.is_fundamental_discriminant(D): + if not D.is_fundamental_discriminant(): return False if D >= 0: return False @@ -6417,7 +6416,7 @@ def ell_heegner_discriminants(self, bound): sage: E.heegner_discriminants(30) # indirect doctest [-7, -8, -19, -24] """ - return [-D for D in range(1, bound) + return [ZZ(-D) for D in range(1, bound) if self.satisfies_heegner_hypothesis(-D)] @@ -6440,7 +6439,7 @@ def ell_heegner_discriminants_list(self, n): [-7, -8, -19, -24] """ v = [] - D = -5 + D = ZZ(-5) while len(v) < n: while not self.satisfies_heegner_hypothesis(D): D -= 1 @@ -7169,7 +7168,7 @@ def _heegner_forms_list(self, D, beta=None, expected_count=None): [389*x^2 + 313*x*y + 63*y^2, 1167*x^2 + 313*x*y + 21*y^2, 3501*x^2 + 313*x*y + 7*y^2] """ if expected_count is None: - expected_count = number_field.QuadraticField(D, 'a').class_number() + expected_count = D.class_number() N = self.conductor() if beta is None: beta = Integers(4*N)(D).sqrt(extend=False) @@ -7229,11 +7228,11 @@ def satisfies_heegner_hypothesis(self, D): sage: E.satisfies_heegner_hypothesis(-11) False """ - if not number_field.is_fundamental_discriminant(D): - return False D = ZZ(D) if D >= 0: return False + if not D.is_fundamental_discriminant(): + return False if D.gcd(self.conductor()) != 1: return False for p, _ in self.conductor().factor(): From 2123c4b7df3f7e495dcaced68b082619956f904d Mon Sep 17 00:00:00 2001 From: John Cremona Date: Fri, 17 Feb 2023 15:40:13 +0000 Subject: [PATCH 02/12] improved efficiency of CM functions through better caching --- src/doc/en/reference/references/index.rst | 10 + src/sage/rings/integer.pyx | 5 +- src/sage/schemes/elliptic_curves/cm.py | 374 ++++++++++++++-------- 3 files changed, 262 insertions(+), 127 deletions(-) diff --git a/src/doc/en/reference/references/index.rst b/src/doc/en/reference/references/index.rst index c4159f90bb2..b829bd36afb 100644 --- a/src/doc/en/reference/references/index.rst +++ b/src/doc/en/reference/references/index.rst @@ -3635,6 +3635,11 @@ REFERENCES: .. [KMOY2007] \M. Kashiwara, K. C. Misra, M. Okado, D. Yamada. *Perfect crystals for* `U_q(D_4^{(3)})`, J. Algebra. **317** (2007). +.. [Klaise2012] Janis Klaise. + *Orders in imaginary quadratic fields of small class number* + University of Warwick Undergraduate Masters thesis, unpublished (2012). + https://warwick.ac.uk/fac/cross_fac/complexity/people/students/dtc/students2013/klaise/janis_klaise_ug_report.pdf + .. [KMR2012] \A. Kleshchev, A. Mathas, and A. Ram, *Universal Specht modules for cyclotomic Hecke algebras*, Proc. London Math. Soc. (2012) 105 (6): 1245-1289. @@ -6085,6 +6090,11 @@ REFERENCES: .. [Watkins] Mark Watkins, *Hypergeometric motives over Q and their L-functions*, http://magma.maths.usyd.edu.au/~watkins/papers/known.pdf +.. [Watkins2004] Mark Watkins. + *Class numbers of imaginary quadratic fields*. + Math. Comp. 73 (2004), 907-938. + https://www.ams.org/journals/mcom/2004-73-246/S0025-5718-03-01517-5/ + .. [Wat2003] Joel Watson. *Strategy: an introduction to game theory*. WW Norton, 2002. diff --git a/src/sage/rings/integer.pyx b/src/sage/rings/integer.pyx index aa7168df78c..e2ac2327e30 100644 --- a/src/sage/rings/integer.pyx +++ b/src/sage/rings/integer.pyx @@ -5567,7 +5567,8 @@ cdef class Integer(sage.structure.element.EuclideanDomainElement): - ``proof`` (boolean, default ``True``) -- if ``False`` then for negative discriminants a faster algorithm is used by the PARI library which is known to give incorrect results - when the class group has many cyclic factors. + when the class group has many cyclic factors. However the + results are correct for discriminants `D` with `|D|\le 2\cdot10^{10}`. OUTPUT: @@ -5576,7 +5577,7 @@ cdef class Integer(sage.structure.element.EuclideanDomainElement): .. NOTE:: - This is not always equal to the number of classes of + For positive `D`, this is not always equal to the number of classes of primitive binary quadratic forms of discriminant `D`, which is equal to the narrow class number. The two notions are the same when `D<0`, or `D>0` and the fundamental unit of diff --git a/src/sage/schemes/elliptic_curves/cm.py b/src/sage/schemes/elliptic_curves/cm.py index bbb43bac1d1..beb15c47d2f 100644 --- a/src/sage/schemes/elliptic_curves/cm.py +++ b/src/sage/schemes/elliptic_curves/cm.py @@ -435,8 +435,10 @@ def cm_j_invariants_and_orders(K, proof=None): # Get the list of CM orders that could possibly have Hilbert class # polynomial F(x) with a root in K. If F(x) has a root alpha in K, # then F is the minimal polynomial of alpha in K, so the degree of - # F(x) is at most [K:QQ]. - dlist = sorted(Df for v in discriminants_with_bounded_class_number(K.degree(), proof=proof).values() for Df in v) + # F(x) divides [K:QQ]. + n = K.absolute_degree() + T = discriminants_with_bounded_class_number(n, proof=proof) + dlist = sorted(sum((Dflist for h,Dflist in T.items() if h.divides(n)), [])) return [(D, f, j) for D, f in dlist for j in hilbert_class_polynomial(D*f*f).roots(K, multiplicities=False)] @@ -456,27 +458,27 @@ def cm_orders(h, proof=None): OUTPUT: - - list of 2-tuples `(D,f)` + - list of 2-tuples `(D,f)` sorted lexicographically by (|D|,f) EXAMPLES:: sage: cm_orders(0) [] sage: v = cm_orders(1); v - [(-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)] + [(-3, 1), (-3, 2), (-3, 3), (-4, 1), (-4, 2), (-7, 1), (-7, 2), (-8, 1), (-11, 1), (-19, 1), (-43, 1), (-67, 1), (-163, 1)] sage: type(v[0][0]), type(v[0][1]) (<... 'sage.rings.integer.Integer'>, <... 'sage.rings.integer.Integer'>) sage: v = cm_orders(2); v - [(-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)] + [(-3, 4), (-3, 5), (-3, 7), (-4, 3), (-4, 4), (-4, 5), (-7, 4), (-8, 2), (-8, 3), (-11, 3), (-15, 1), (-15, 2), (-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)] sage: len(v) 29 sage: set([hilbert_class_polynomial(D*f^2).degree() for D,f in v]) {2} - Any degree up to 100 is implemented, but may be prohibitively slow:: + Any degree up to 100 is implemented, but may be slow:: sage: cm_orders(3) - [(-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)] + [(-3, 6), (-3, 9), (-11, 2), (-19, 2), (-23, 1), (-23, 2), (-31, 1), (-31, 2), (-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)] sage: len(cm_orders(4)) 84 """ @@ -484,13 +486,16 @@ def cm_orders(h, proof=None): if h <= 0: # trivial case return [] - # Get information for all discriminants then throw away everything - # but for h. If this is replaced by a table it will be faster, - # but not now. (David Kohel is rumored to have a large table.) - return discriminants_with_bounded_class_number(h, proof=proof)[h] + + if h in hDf_dict: + return hDf_dict[h] + else: # Get all discriminants for all class numbers up to h (which will + # be stored in hDf_dict), and return just those with class number h. + return discriminants_with_bounded_class_number(h, proof=proof)[h] # Table from Mark Watkins paper "Class numbers of imaginary quadratic fields". -# I extracted this by cutting/pasting from the pdf, and running this program: + +# WAS extracted this by cutting/pasting from the pdf, and running this program: # z = {} # for X in open('/Users/wstein/tmp/a.txt').readlines(): # if len(X.strip()): @@ -498,6 +503,10 @@ def cm_orders(h, proof=None): # for i in range(5): # z[v[3*i]]=(v[3*i+2], v[3*i+1]) +# The keys are integers 1--100 and the value for each h is (|D|,n) +# where |D| is the largest absolute discriminant of an imaginary +# quadratic field with class number h, and n is the number of such +# fields. These are all *unconditional* (not dependent on GRH). watkins_table = {1: (163, 9), 2: (427, 18), 3: (907, 16), 4: (1555, 54), 5: (2683, 25), 6: (3763, 51), 7: (5923, 31), 8: (6307, 131), 9: (10627, 34), 10: @@ -527,13 +536,60 @@ def cm_orders(h, proof=None): 95:(1659067, 241), 96: (1684027, 3283), 97: (1842523, 185), 98: (2383747,580), 99: (1480627, 289), 100: (1856563, 1736)} +# Table from Janis Klaise [Klaise2012]_ + +# Extracted by converting pdf to text via pdf2ps and ps2txt, cutting/pasting +# and running this code: + +# klaise_table = {} +# for X in open('klaise_table.txt').readlines(): +# if len(X.strip()): +# v = [int(a) for a in X.split()] +# for i in range(4): +# klaise_table[v[3*i]]=(v[3*i+2], v[3*i+1]) + +# The keys are integers 1--100 and the value for each h is (|D|,n) +# where |D| is the largest discriminant of an imaginary quadratic +# field with class number h, and n is the number of such fields. +# These are all *unconditional* (not dependent on GRH). + +klaise_table = {1: (163, 13), 2: (427, 29), 3: (907, 25), 4: (1555, 84), 5: (2683, 29), 6: (4075, 101), + 7: (5923, 38), 8: (7987, 208), 9: (10627, 55), 10: (13843, 123), 11: (15667, 46), + 12: (19723, 379), 13: (20563, 43), 14: (30067, 134), 15: (34483, 95), 16: (35275, 531), + 17: (37123, 50), 18: (48427, 291), 19: (38707, 59), 20: (58843, 502), 21: (61483, 118), + 22: (85507, 184), 23: (90787, 78), 24: (111763, 1042), 25: (93307, 101), 26: (103027, 227), + 27: (103387, 136), 28: (126043, 623), 29: (166147, 94), 30: (137083, 473), 31: (133387, 83), + 32: (164803, 1231), 33: (222643, 158), 34: (189883, 262), 35: (210907, 111), 36: (217627, 1306), + 37: (158923, 96), 38: (289963, 284), 39: (253507, 162), 40: (274003, 1418), 41: (296587, 125), + 42: (301387, 596), 43: (300787, 123), 44: (319867, 911), 45: (308323, 231), 46: (462883, 330), + 47: (375523, 117), 48: (335203, 2895), 49: (393187, 146), 50: (389467, 445), 51: (546067, 217), + 52: (457867, 1006), 53: (425107, 130), 54: (532123, 812), 55: (452083, 177), 56: (494323, 1812), + 57: (615883, 237), 58: (586987, 361), 59: (474307, 144), 60: (662803, 2361), 61: (606643, 149), + 62: (647707, 386), 63: (991027, 311), 64: (693067, 2919), 65: (703123, 192), 66: (958483, 861), + 67: (652723, 145), 68: (819163, 1228), 69: (888427, 292), 70: (821683, 704), 71: (909547, 176), + 72: (947923, 4059), 73: (886867, 137), 74: (951043, 474), 75: (916507, 353), 76: (1086187, 1384), + 77: (1242763, 236), 78: (1004347, 925), 79: (1333963, 200), 80: (1165483, 3856), 81: (1030723, 339), + 82: (1446547, 487), 83: (1074907, 174), 84: (1225387, 2998), 85: (1285747, 246), 86: (1534723, 555), + 87: (1261747, 313), 88: (1265587, 2771), 89: (1429387, 206), 90: (1548523, 1516), 91: (1391083, 249), + 92: (1452067, 1591), 93: (1475203, 354), 94: (1587763, 600), 95: (1659067, 273), 96: (1684027, 7276), + 97: (1842523, 208), 98: (2383747, 710), 99: (1480627, 396), 100: (1856563, 2311)} + def largest_fundamental_disc_with_class_number(h): - """ - Return largest absolute value of any fundamental discriminant with - class number `h`, and the number of fundamental discriminants with - that class number. This is known for `h` up to 100, by work of Mark - Watkins. + r""" + Return largest absolute value of any fundamental negative discriminant with + class number `h`, and the number of fundamental negative discriminants with + that class number. This is known (unconditionally) for `h` up to 100, + by work of Mark Watkins ([Watkins2004]_). + + .. NOTE:: + + The class number of a fundamental negative discriminant `D` is + the same as the class number of the imaginary quadratic field + `\QQ(\sqrt{D})`, so this function gives the number of such + fields of each class number `h\le100`. It is easy to extend + this to larger class number conditional on the GRH, but much + harder to obyain unconditional results. INPUT: @@ -541,60 +597,124 @@ class number `h`, and the number of fundamental discriminants with EXAMPLES:: - sage: sage.schemes.elliptic_curves.cm.largest_fundamental_disc_with_class_number(0) + sage: from sage.schemes.elliptic_curves.cm import largest_fundamental_disc_with_class_number + sage: largest_fundamental_disc_with_class_number(0) (0, 0) - sage: sage.schemes.elliptic_curves.cm.largest_fundamental_disc_with_class_number(1) + sage: largest_fundamental_disc_with_class_number(1) (163, 9) - sage: sage.schemes.elliptic_curves.cm.largest_fundamental_disc_with_class_number(2) + sage: largest_fundamental_disc_with_class_number(2) (427, 18) - sage: sage.schemes.elliptic_curves.cm.largest_fundamental_disc_with_class_number(10) + sage: largest_fundamental_disc_with_class_number(10) (13843, 87) - sage: sage.schemes.elliptic_curves.cm.largest_fundamental_disc_with_class_number(100) + sage: largest_fundamental_disc_with_class_number(100) (1856563, 1736) - sage: sage.schemes.elliptic_curves.cm.largest_fundamental_disc_with_class_number(101) + sage: largest_fundamental_disc_with_class_number(101) Traceback (most recent call last): ... - NotImplementedError: largest discriminant not known for class number 101 + NotImplementedError: largest fundamental discriminant not available for class number 101 + """ h = Integer(h) if h <= 0: - # very easy special case return Integer(0), Integer(0) try: - # simply look up the answer in Watkins's table. B, c = watkins_table[h] return (Integer(B), Integer(c)) except KeyError: - # nobody knows, since I guess Watkins's is state of the art. - raise NotImplementedError("largest discriminant not known for class number %s" % h) + raise NotImplementedError("largest fundamental discriminant not available for class number %s" % h) + +def largest_disc_with_class_number(h): + r""" + Return largest absolute value of any negative discriminant with + class number `h`, and the number of fundamental negative + discriminants with that class number. This is known + (unconditionally) for `h` up to 100, by work of Mark Watkins + [Watkins2004]_ for fundamental discriminants, extended to all + discriminants of class number `h\le100` by Klaise [Klaise2012]_. + + .. NOTE:: + + The class number of a negative discriminant `D` is + the same as the class number of the unique imaginary quadratic order + of discriminant `D`, so this function gives the number of such + orders of each class number `h\le100`. It is easy to extend + this to larger class number conditional on the GRH, but much + harder to obyain unconditional results. + + INPUT: + + - `h` -- integer + + EXAMPLES:: + + sage: from sage.schemes.elliptic_curves.cm import largest_disc_with_class_number + sage: largest_disc_with_class_number(0) + (0, 0) + sage: largest_disc_with_class_number(1) + (163, 13) + sage: largest_disc_with_class_number(2) + (427, 29) + sage: largest_disc_with_class_number(10) + (13843, 123) + sage: largest_disc_with_class_number(100) + (1856563, 2311) + sage: largest_disc_with_class_number(101) + Traceback (most recent call last): + ... + NotImplementedError: largest discriminant not available for class number 101 + + For most `h\le100`, the largest fundamental discriminant with class number `h` is also the largest discriminant, but this is not the case for `h=`:: + + sage: from sage.schemes.elliptic_curves.cm import largest_disc_with_class_number, largest_fundamental_disc_with_class_number + sage: [h for h in range(1,101) if largest_disc_with_class_number(h)[0] != largest_fundamental_disc_with_class_number(h)[0]] + [6, 8, 12, 16, 20, 30, 40, 42, 52, 70] + sage: largest_fundamental_disc_with_class_number(6) + (3763, 51) + sage: largest_disc_with_class_number(6) + (4075, 101) + """ + h = Integer(h) + if h <= 0: + return Integer(0), Integer(0) + try: + B, c = klaise_table[h] + return (Integer(B), Integer(c)) + except KeyError: + raise NotImplementedError("largest discriminant not available for class number %s" % h) +# This dict has class numbers h as keys, the value at h is a complete +# list of paits (D0,f) such that D=D0*f**2 has class number h. We +# initialise it with h=1 only; other values will be added by calls to +# discriminants_with_bounded_class_number(). + +hDf_dict = {ZZ(1): [(ZZ(D), ZZ(h)) for D,h in + [(-3, 1), (-3, 2), (-3, 3), (-4, 1), (-4, 2), (-7, 1), (-7, 2), + (-8, 1), (-11, 1), (-19, 1), (-43, 1), (-67, 1), (-163, 1)]]} @cached_function def discriminants_with_bounded_class_number(hmax, B=None, proof=None): r""" - Return dictionary with keys class numbers `h\le hmax` and values the - list of all pairs `(D, f)`, with `D<0` a fundamental discriminant such - that `Df^2` has class number `h`. If the optional bound `B` is given, - return only those pairs with fundamental `|D| \le B`, though `f` can - still be arbitrarily large. + Return a dictionary with keys class numbers `h\le hmax` and values the + list of all pairs `(D_0, f)`, with `D_0<0` a fundamental discriminant such + that `D=D_0f^2` has class number `h`. If the optional bound `B` is given, + return only those pairs with fundamental `|D| \le B`. INPUT: - ``hmax`` -- integer - `B` -- integer or None; if None returns all pairs - ``proof`` -- this code calls the PARI function :pari:`qfbclassno`, so it - could give wrong answers when ``proof``==``False``. The default is - whatever ``proof.number_field()`` is. If ``proof==False`` and `B` is - ``None``, at least the number of discriminants is correct, since it - is double checked with Watkins's table. + could give wrong answers when ``proof``==``False`` (though only for + discriminants greater than `2\cdot10^{10}). The default is + the current value of ``proof.number_field()``. OUTPUT: - dictionary - In case `B` is not given, we use Mark Watkins's: "Class numbers of - imaginary quadratic fields" to compute a `B` that captures all `h` - up to `hmax` (only available for `hmax\le100`). + In case `B` is not given, for ``hmax`` at most 100 we use the + tables from [Watkins2004]_ and [Klaise2012]_ to compute a `B` that + captures all `h` up to `hmax` (only available for `hmax\le100`). EXAMPLES:: @@ -602,11 +722,11 @@ def discriminants_with_bounded_class_number(hmax, B=None, proof=None): sage: sorted(v) [1, 2, 3] sage: v[1] - [(-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)] + [(-3, 1), (-3, 2), (-3, 3), (-4, 1), (-4, 2), (-7, 1), (-7, 2), (-8, 1), (-11, 1), (-19, 1), (-43, 1), (-67, 1), (-163, 1)] sage: v[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), (-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)] + [(-3, 4), (-3, 5), (-3, 7), (-4, 3), (-4, 4), (-4, 5), (-7, 4), (-8, 2), (-8, 3), (-11, 3), (-15, 1), (-15, 2), (-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)] sage: v[3] - [(-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)] + [(-3, 6), (-3, 9), (-11, 2), (-19, 2), (-23, 1), (-23, 2), (-31, 1), (-31, 2), (-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)] sage: v = sage.schemes.elliptic_curves.cm.discriminants_with_bounded_class_number(8, proof=False) sage: sorted(len(v[h]) for h in v) [13, 25, 29, 29, 38, 84, 101, 208] @@ -614,110 +734,114 @@ def discriminants_with_bounded_class_number(hmax, B=None, proof=None): Find all class numbers for discriminant up to 50:: sage: sage.schemes.elliptic_curves.cm.discriminants_with_bounded_class_number(hmax=5, B=50) - {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)]} + {1: [(-3, 1), (-3, 2), (-3, 3), (-4, 1), (-4, 2), (-7, 1), (-7, 2), (-8, 1), (-11, 1), (-19, 1), (-43, 1)], 2: [(-3, 4), (-4, 3), (-8, 2), (-15, 1), (-20, 1), (-24, 1), (-35, 1), (-40, 1)], 3: [(-11, 2), (-23, 1), (-31, 1)], 4: [(-39, 1)], 5: [(-47, 1)]} + """ + hmax = Integer(hmax) + global hDf_dict + + # Easy case where we have already computed and cached the relevant values + if hDf_dict and hmax <= max(hDf_dict.keys()): + T = {h:Dflist for h,Dflist in hDf_dict.items() if h<=hmax} + if B: + for h in T: + T[h] = [Df for Df in T[h] if Df[0].abs()*Df[1]**2<=B] + return T + # imports that are needed only for this function + from sage.arith.srange import xsrange from sage.structure.proof.proof import get_flag - - # deal with input defaults and type checking proof = get_flag(proof, 'number_field') - hmax = Integer(hmax) - - # T stores the output - T = {} - - # Easy case -- instead of giving error, give meaningful output - if hmax < 1: - return T if B is None: - # Determine how far we have to go by applying Watkins's theorem. - v = [largest_fundamental_disc_with_class_number(h) for h in range(1, hmax+1)] + # Determine how far we have to go by applying Watkins + Klaise's results. + v = [largest_disc_with_class_number(h) for h in range(1, hmax+1)] B = max([b for b,_ in v]) - fund_count = [0] + [cnt for _,cnt in v] + #print("Testing all discriminants up to {}".format(B)) + count = [0] + [cnt for _,cnt in v] else: # Nothing to do -- set to None so we can use this later to know not # to do a double check about how many we find. - fund_count = None + count = None B = Integer(B) if B <= 2: # This is an easy special case, since there are no fundamental discriminants # this small. - return T + return {} + + # T stores the values found, to be returned. It will also be used + # to update the global hDf_dict *provided that* the parameter B + # was not provided. Note that we only reach this point if we have + # not already computed data for at least one class number, hmax. + + # To avoid recomputing class numbers, we initialise T with the + # global hDf_dict, only including those (D,f) for which |D|*f**2 + # <= B if B was provided: + + T = hDf_dict.copy() + + if not count: + for h in T: + T[h] = [Df for Df in T[h] if Df[0].abs()*Df[1]**2<=B] + + # h_dict caches the class number h of all discriminants previously + # encountered; we will use the function OrderClassNumber() to + # quicky compute the class number of non-fundamental discriminants + # from the fundamental ones. Note that in the initialisation, the + # keys of h_dict include nonfundamental discriminants, but we only + # update it with fundamental ones. + + h_dict = {} + for h, Dflist in hDf_dict.items(): + for D0,f in Dflist: + h_dict[D0*f**2] = h + + # We do not need to certaify the class number from :pari:`qfbclassno` for discriminants under 2*10^10 + if B > 2*10**10: + proof = False + + for D in xsrange(-3, -B-1, -1): + if not D.is_discriminant(): + continue + D0 = D.squarefree_part() + if D0%4 !=1: + D0 *= 4 + f = (D//D0).isqrt() + + # Now D0 is the fundamental discriminant and f the conductor - # This lower bound gets used in an inner loop below. - from math import log - - def lb(f): - """Lower bound on euler_phi.""" - # 1.79 > e^gamma = 1.7810724... - if f <= 1: - return 0 # don't do log(log(1)) = log(0) - llf = log(log(f)) - return f/(1.79*llf + 3.0/llf) - - for D in range(-B, -2): - D = Integer(D) - if D.is_fundamental_discriminant(): - h_D = D.class_number(proof) - # For each fundamental discriminant D, loop through the f's such - # that h(D*f^2) could possibly be <= hmax. As explained to me by Cremona, - # we have h(D*f^2) >= (1/c)*h(D)*phi_D(f) >= (1/c)*h(D)*euler_phi(f), where - # phi_D(f) is like euler_phi(f) but the factor (1-1/p) is replaced - # by a factor of (1-kr(D,p)*p), where kr(D/p) is the Kronecker symbol. - # The factor c is 1 unless D=-4 and f>1 (when c=2) or D=-3 and f>1 (when c=3). - # Since (1-1/p) <= 1 and (1-1/p) <= (1+1/p), we see that - # euler_phi(f) <= phi_D(f). - # - # We have the following analytic lower bound on euler_phi: - # - # euler_phi(f) >= lb(f) = f / (exp(euler_gamma)*log(log(n)) + 3/log(log(n))). - # - # See Theorem 8 of Peter Clark's - # http://math.uga.edu/~pete/4400arithmeticorders.pdf - # which is a consequence of Theorem 15 of - # [Rosser and Schoenfeld, 1962]. - # - # By Calculus, we see that the lb(f) is an increasing function of f >= 2. - # - # NOTE: You can visibly "see" that it is a lower bound in Sage with - # lb(n) = n/(exp(euler_gamma)*log(log(n)) + 3/log(log(n))) - # plot(lb, (n, 1, 10^4), color='red') + plot(lambda x: euler_phi(int(x)), 1, 10^4).show() - # - # So we consider f=1,2,..., until the first f with lb(f)*h_D > c*h_max. - # (Note that lb(f) is <= 0 for f=1,2, so nothing special is needed there.) - # - # TODO: Maybe we could do better using a bound for phi_D(f). - # - f = Integer(1) - chmax=hmax - if D==-3: - chmax*=3 + if D in h_dict: + h = h_dict[D] + else: + if f == 1: # D itself is fundamental + h = D.class_number(proof) + h_dict[D] = h + else: + h = OrderClassNumber(D0,h_dict[D0],f) + + # If the class number of this order is within the range, then store (D0,f) + if h <= hmax: + z = (D0, f) + if h in T: + if z not in T[h]: + T[h].append(z) else: - if D==-4: - chmax*=2 - while lb(f)*h_D <= chmax: - h = OrderClassNumber(D,h_D,f) - # If the class number of this order is within the range, then use it. - if h <= hmax: - z = (D, f) - if h in T: - T[h].append(z) - else: - T[h] = [z] - f += 1 + T[h] = [z] + + # sort each list of (D,f) pairs by (|D|,f) for h in T: - T[h] = list(reversed(T[h])) + T[h].sort(key=lambda Df: (Df[0].abs(), Df[1])) - if fund_count is not None: - # Double check that we found the right number of fundamental - # discriminants; we might as well, since Watkins provides this - # data. + # count is None unless the user provided a value of B + if count is not None: + # 1. Check that we found the right number of discriminants for h in T: - if len([DD for DD, ff in T[h] if ff == 1]) != fund_count[h]: + if len(T[h]) != count[h]: raise RuntimeError("number of discriminants inconsistent with Watkins's table") + # 2. Update the global dict + hDf_dict.update(T) return T From 988e88b457a956cad5d6ddd2c947638761ca6f2c Mon Sep 17 00:00:00 2001 From: John Cremona Date: Mon, 20 Feb 2023 19:15:55 +0000 Subject: [PATCH 03/12] fix lint and docstring errors --- src/sage/schemes/elliptic_curves/cm.py | 4 ++-- src/sage/schemes/elliptic_curves/ell_finite_field.py | 2 -- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/sage/schemes/elliptic_curves/cm.py b/src/sage/schemes/elliptic_curves/cm.py index beb15c47d2f..f3d01e5ea21 100644 --- a/src/sage/schemes/elliptic_curves/cm.py +++ b/src/sage/schemes/elliptic_curves/cm.py @@ -193,7 +193,7 @@ def is_HCP(f, check_monic_irreducible=True): Cremona and Sutherland: Algorithm 2 of [CreSuth2023]_. - EXAMPLES:: + EXAMPLES: Even for large degrees this is fast. We test the largest discriminant of class number 100, for which the HCP has coefficients @@ -293,7 +293,7 @@ def OrderClassNumber(D0,h0,f): .. MATH:: - h(D) = \frac{h(D_0)f}{[\mathcal{O}_{D_0}^\times:\mathcal{O}_{D}^\times]}\prod_{p\,|\,f}\left(1-\left(\frac{D_0}{p}\right)\frac{1}{p}\right) + h(D) = \frac{h(D_0)f}{[\mathcal{O}_{D_0}^\times:\mathcal{O}_{D}^\times]}\prod_{p\,|\,f}\left(1-\left(\frac{D_0}{p}\right)\frac{1}{p}\right) EXAMPLES:: diff --git a/src/sage/schemes/elliptic_curves/ell_finite_field.py b/src/sage/schemes/elliptic_curves/ell_finite_field.py index dc21b35985a..456a0da6117 100644 --- a/src/sage/schemes/elliptic_curves/ell_finite_field.py +++ b/src/sage/schemes/elliptic_curves/ell_finite_field.py @@ -1798,5 +1798,3 @@ def is_j_supersingular(j, proof=True): # expensive since it involves counting the number of points on E): return E.trace_of_frobenius() % p == 0 - - From e866fe994ba99b61ae5d751b2cc972193c646109 Mon Sep 17 00:00:00 2001 From: John Cremona Date: Mon, 20 Feb 2023 19:33:48 +0000 Subject: [PATCH 04/12] fix another docstring error --- src/sage/schemes/elliptic_curves/cm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sage/schemes/elliptic_curves/cm.py b/src/sage/schemes/elliptic_curves/cm.py index f3d01e5ea21..cf89334b383 100644 --- a/src/sage/schemes/elliptic_curves/cm.py +++ b/src/sage/schemes/elliptic_curves/cm.py @@ -458,7 +458,7 @@ def cm_orders(h, proof=None): OUTPUT: - - list of 2-tuples `(D,f)` sorted lexicographically by (|D|,f) + - list of 2-tuples `(D,f)` sorted lexicographically by `(|D|, f)` EXAMPLES:: From be1b18f689cb862f2bfbd2366ad9151762feca1a Mon Sep 17 00:00:00 2001 From: John Cremona Date: Tue, 21 Feb 2023 08:55:07 +0000 Subject: [PATCH 05/12] fix typos in docstrings --- src/sage/schemes/elliptic_curves/cm.py | 2 +- src/sage/schemes/elliptic_curves/ell_finite_field.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/sage/schemes/elliptic_curves/cm.py b/src/sage/schemes/elliptic_curves/cm.py index cf89334b383..62dfb72906c 100644 --- a/src/sage/schemes/elliptic_curves/cm.py +++ b/src/sage/schemes/elliptic_curves/cm.py @@ -705,7 +705,7 @@ def discriminants_with_bounded_class_number(hmax, B=None, proof=None): - `B` -- integer or None; if None returns all pairs - ``proof`` -- this code calls the PARI function :pari:`qfbclassno`, so it could give wrong answers when ``proof``==``False`` (though only for - discriminants greater than `2\cdot10^{10}). The default is + discriminants greater than `2\cdot10^{10}`). The default is the current value of ``proof.number_field()``. OUTPUT: diff --git a/src/sage/schemes/elliptic_curves/ell_finite_field.py b/src/sage/schemes/elliptic_curves/ell_finite_field.py index 1043eac394b..88117d35be5 100644 --- a/src/sage/schemes/elliptic_curves/ell_finite_field.py +++ b/src/sage/schemes/elliptic_curves/ell_finite_field.py @@ -1136,7 +1136,7 @@ def is_isogenous(self, other, field=None, proof=True): ... ValueError: Curves have different base fields: use the field parameter. - When the field is given: + When the field is given:: sage: E1 = EllipticCurve(GF(13^2,'a'),[2,7]); E1 Elliptic Curve defined by y^2 = x^3 + 2*x + 7 over Finite Field in a of size 13^2 @@ -1439,7 +1439,7 @@ def height_above_floor(self, ell, e): .. NOTE:: - For an ordinary `E/\F_q`, and a prime `\ell`, the height + For an ordinary `E/\GF{q}`, and a prime `\ell`, the height of the `\ell`-volcano containing `j(E)` is the `\ell`-adic valuation of the conductor the order generated by the Frobenius `\pi_E`); the height of `j(E)` on its From 64f5fe43163d0a641770f81bc2962f756fe8a160 Mon Sep 17 00:00:00 2001 From: John Cremona Date: Wed, 22 Feb 2023 11:06:09 +0000 Subject: [PATCH 06/12] docstring typo --- src/sage/schemes/elliptic_curves/cm.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sage/schemes/elliptic_curves/cm.py b/src/sage/schemes/elliptic_curves/cm.py index 62dfb72906c..7fa881d7c1a 100644 --- a/src/sage/schemes/elliptic_curves/cm.py +++ b/src/sage/schemes/elliptic_curves/cm.py @@ -187,7 +187,7 @@ def is_HCP(f, check_monic_irreducible=True): OUTPUT: - (integer) -- `D` if ``f`` is the Hilbert Class Polynomial (HCP) `H_D`, esle 0. + (integer) -- `D` if ``f`` is the Hilbert Class Polynomial (HCP) `H_D`, else 0. ALGORITHM: From b6fa715c4188eb2b27b348e1631ede4c12f06eee Mon Sep 17 00:00:00 2001 From: John Cremona Date: Wed, 22 Feb 2023 11:16:29 +0000 Subject: [PATCH 07/12] fix failing doctests --- .../number_field_element_quadratic.pyx | 20 ++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/src/sage/rings/number_field/number_field_element_quadratic.pyx b/src/sage/rings/number_field/number_field_element_quadratic.pyx index ada2e1a0b65..bcba57457e5 100644 --- a/src/sage/rings/number_field/number_field_element_quadratic.pyx +++ b/src/sage/rings/number_field/number_field_element_quadratic.pyx @@ -2137,7 +2137,7 @@ cdef class NumberFieldElement_quadratic(NumberFieldElement_absolute): - ``var`` -- the minimal polynomial is defined over a polynomial ring in a variable with this name. If not specified this defaults to ``x`` - - ``algorithm`` -- for compatibility with general number field elements: and ignored + - ``algorithm`` -- for compatibility with general number field elements; ignored EXAMPLES:: @@ -2658,11 +2658,18 @@ cdef class OrderElement_quadratic(NumberFieldElement_quadratic): """ return ZZ(NumberFieldElement_quadratic.trace(self)) - def charpoly(self, var='x'): + def charpoly(self, var='x', algorithm=None): r""" The characteristic polynomial of this element, which is over `\ZZ` because this element is an algebraic integer. + INPUT: + + - ``var`` -- the minimal polynomial is defined over a polynomial ring + in a variable with this name. If not specified this defaults to ``x`` + - ``algorithm`` -- for compatibility with general number field elements; ignored + + EXAMPLES:: sage: K. = NumberField(x^2 - 5) @@ -2678,10 +2685,17 @@ cdef class OrderElement_quadratic(NumberFieldElement_quadratic): R = ZZ[var] return R([self.norm(), -self.trace(), 1]) - def minpoly(self, var='x'): + def minpoly(self, var='x', algorithm=None): r""" The minimal polynomial of this element over `\ZZ`. + INPUT: + + - ``var`` -- the minimal polynomial is defined over a polynomial ring + in a variable with this name. If not specified this defaults to ``x`` + - ``algorithm`` -- for compatibility with general number field elements; ignored + + EXAMPLES:: sage: K. = NumberField(x^2 + 163) From 80d8ecf077d306232c336e546721a402fb9c2f7b Mon Sep 17 00:00:00 2001 From: John Cremona Date: Wed, 22 Feb 2023 14:53:32 +0000 Subject: [PATCH 08/12] fix a doctest in a thematic tutorial --- .../explicit_methods_in_number_theory/nf_galois_groups.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/doc/en/thematic_tutorials/explicit_methods_in_number_theory/nf_galois_groups.rst b/src/doc/en/thematic_tutorials/explicit_methods_in_number_theory/nf_galois_groups.rst index fb7697587c5..9e5b3b8fbf0 100644 --- a/src/doc/en/thematic_tutorials/explicit_methods_in_number_theory/nf_galois_groups.rst +++ b/src/doc/en/thematic_tutorials/explicit_methods_in_number_theory/nf_galois_groups.rst @@ -188,7 +188,7 @@ field). :: sage: w = [-1,-2,..,-200] - sage: v = [D for D in w if is_fundamental_discriminant(D)] + sage: v = [D for D in w if D.is_fundamental_discriminant()] sage: v [-3, -4, -7, -8, -11, -15, -19, -20, ..., -195, -199] @@ -222,7 +222,7 @@ class number :math:`1`! :: sage: w = [1..1000] - sage: v = [D for D in w if is_fundamental_discriminant(D)] + sage: v = [D for D in w if D.is_fundamental_discriminant()] sage: len(v) 302 sage: len([D for D in v if QuadraticField(D,'a').class_number() == 1]) From 8622db969769880fd84fec80164d4a67545b0c7b Mon Sep 17 00:00:00 2001 From: roed314 Date: Fri, 3 Mar 2023 10:12:09 -0500 Subject: [PATCH 09/12] Fix documentation typos --- src/sage/schemes/elliptic_curves/ell_finite_field.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/sage/schemes/elliptic_curves/ell_finite_field.py b/src/sage/schemes/elliptic_curves/ell_finite_field.py index 88117d35be5..3ecaab7b068 100644 --- a/src/sage/schemes/elliptic_curves/ell_finite_field.py +++ b/src/sage/schemes/elliptic_curves/ell_finite_field.py @@ -1441,8 +1441,8 @@ def height_above_floor(self, ell, e): For an ordinary `E/\GF{q}`, and a prime `\ell`, the height of the `\ell`-volcano containing `j(E)` is the `\ell`-adic - valuation of the conductor the order generated by the - Frobenius `\pi_E`); the height of `j(E)` on its + valuation of the conductor of the order generated by the + Frobenius `\pi_E`; the height of `j(E)` on its ell-volcano is the `\ell`-adic valuation of the conductor of the order `\End(E)`. From 68847d44e5778d5635a58923f0db2caf31c83872 Mon Sep 17 00:00:00 2001 From: John Cremona Date: Mon, 6 Mar 2023 11:01:23 +0000 Subject: [PATCH 10/12] dealt with most referee suggestions --- src/sage/schemes/elliptic_curves/cm.py | 130 ++++++++++-------- .../elliptic_curves/ell_finite_field.py | 33 ++--- 2 files changed, 86 insertions(+), 77 deletions(-) diff --git a/src/sage/schemes/elliptic_curves/cm.py b/src/sage/schemes/elliptic_curves/cm.py index 7fa881d7c1a..914cf5a159d 100644 --- a/src/sage/schemes/elliptic_curves/cm.py +++ b/src/sage/schemes/elliptic_curves/cm.py @@ -176,7 +176,7 @@ def hilbert_class_polynomial(D, algorithm=None): def is_HCP(f, check_monic_irreducible=True): r""" - Return ``True``, `D` if ``f`` is the Hilbert Class Polynomial `H_D`, else ``False``, 0. + Determine whether a polynomial is a Hilbert Class Polynomial. INPUT: @@ -187,7 +187,8 @@ def is_HCP(f, check_monic_irreducible=True): OUTPUT: - (integer) -- `D` if ``f`` is the Hilbert Class Polynomial (HCP) `H_D`, else 0. + (integer) -- either `D` if ``f`` is the Hilbert Class Polynomial + `H_D` for discriminant `D`, or `0` if not an HCP. ALGORITHM: @@ -209,35 +210,34 @@ def is_HCP(f, check_monic_irreducible=True): sage: max(c for c in H).ndigits() 2774 sage: is_HCP(H) - (True, -1856563) + -1856563 Testing polynomials which are not HCPs is faster:: sage: is_HCP(H+1) - (False, 0) + 0 TESTS:: sage: from sage.schemes.elliptic_curves.cm import is_HCP - sage: all(is_HCP(hilbert_class_polynomial(D))==(True,D) for D in srange(-4,-100,-1) if D.is_discriminant()) + sage: all(is_HCP(hilbert_class_polynomial(D))==D for D in srange(-4,-100,-1) if D.is_discriminant()) True - sage: all(is_HCP(hilbert_class_polynomial(D)+1)==(False,0) for D in srange(-4,-100,-1) if D.is_discriminant()) + sage: all(not is_HCP(hilbert_class_polynomial(D)+1) for D in srange(-4,-100,-1) if D.is_discriminant()) True - """ zero = ZZ(0) # optional check that input is monic and irreducible if check_monic_irreducible: try: if not (all(c in ZZ for c in f) and f.is_monic()): - return (False, zero) + return zero f = f.change_ring(ZZ) except AttributeError: - return (False, zero) + return zero h = f.degree() - h2list = [d for d in h.divisors() if (d-h)%2==0 and d.prime_to_m_part(2)==1] + h2list = [d for d in h.divisors() if (d-h)%2 == 0 and d.prime_to_m_part(2) == 1] pmin = 33 * (h**2 * (RR(h+2).log().log()+2)**2).ceil() # Guarantees 4*p > |D| for fundamental D under GRH p = pmin-1 @@ -247,30 +247,25 @@ def is_HCP(f, check_monic_irreducible=True): p = next_prime(p) n += 1 fp = f.change_ring(GF(p)) - # Compute X^p-X mod fp manually, avoiding quotient ring which is slower - r = zpow = z = fp.parent().gen() - m = p>>1 - while m: - zpow = (zpow**2) % fp - if m & 1: - r = (zpow * r) % fp - m >>= 1 - # now r = X^p mod fp - d = (r-z).gcd(fp).degree() # number of roots mod p - if d==0: + # Compute X^p-X mod fp + z = fp.parent().gen() + r = pow(z, p, fp) - z + d = r.gcd(fp).degree() # number of roots mod p + if d == 0: continue if not fp.is_squarefree(): continue if d100 you must specify a discriminant bound B") else: # Nothing to do -- set to None so we can use this later to know not # to do a double check about how many we find. @@ -797,8 +798,8 @@ def discriminants_with_bounded_class_number(hmax, B=None, proof=None): for D0,f in Dflist: h_dict[D0*f**2] = h - # We do not need to certaify the class number from :pari:`qfbclassno` for discriminants under 2*10^10 - if B > 2*10**10: + # We do not need to certify the class number from :pari:`qfbclassno` for discriminants under 2*10^10 + if B < 2*10**10: proof = False for D in xsrange(-3, -B-1, -1): @@ -834,7 +835,7 @@ def discriminants_with_bounded_class_number(hmax, B=None, proof=None): for h in T: T[h].sort(key=lambda Df: (Df[0].abs(), Df[1])) - # count is None unless the user provided a value of B + # count is None precisely when the user provided a value of B if count is not None: # 1. Check that we found the right number of discriminants for h in T: @@ -847,17 +848,19 @@ def discriminants_with_bounded_class_number(hmax, B=None, proof=None): @cached_function -def is_cm_j_invariant(j, method='CremonaSutherland'): +def is_cm_j_invariant(j, algorithm='CremonaSutherland', method=None): r"""Return whether or not this is a CM `j`-invariant, and the CM discriminant if it is. INPUT: - ``j`` -- an element of a number field `K` - - ``method`` (string, default 'CremonaSutherland') -- the method + - ``algorithm`` (string, default 'CremonaSutherland') -- the algorithm used, either 'CremonaSutherland' (the default, very much faster for all but very small degrees), 'exhaustive' or 'reduction' + - ``method`` (string) -- deprecated name for ``algorithm`` + OUTPUT: A pair (bool, (d,f)) which is either (False, None) if `j` is not a @@ -872,7 +875,7 @@ def is_cm_j_invariant(j, method='CremonaSutherland'): :func:`is_HCP` which implements Algorithm 2 of [CreSuth2023]_ by Cremona and Sutherland. - Two older methods are available, both of which are much slower + Two older algorithms are available, both of which are much slower except for very small degrees. Method 'exhaustive' makes use of the complete and unconditionsl classification of @@ -917,6 +920,11 @@ def is_cm_j_invariant(j, method='CremonaSutherland'): True """ + if method: + if not algorithm: + algorithm = method + raise DeprecationWarning("'method' is deprecated, use 'algorithm instead'") + # First we check that j is an algebraic number: from sage.rings.all import NumberFieldElement, NumberField if not isinstance(j, NumberFieldElement) and j not in QQ: @@ -938,7 +946,7 @@ def is_cm_j_invariant(j, method='CremonaSutherland'): # Next we find its minimal polynomial of j: - if j.parent().absolute_degree()==2: + if j.parent().absolute_degree() == 2: jpol = j.absolute_minpoly() # no algorithm parameter else: jpol = j.absolute_minpoly(algorithm='pari') @@ -951,9 +959,9 @@ def is_cm_j_invariant(j, method='CremonaSutherland'): # Otherwise test whether it is a Hilbert Class Polynomial # (using the fact that we know that it is monic and irreducible): - if method == 'CremonaSutherland': - flag, D = is_HCP(jpol, check_monic_irreducible=False) - if flag: + if algorithm == 'CremonaSutherland': + D = is_HCP(jpol, check_monic_irreducible=False) + if D: D0 = D.squarefree_part() if D0%4 !=1: D0 *= 4 @@ -963,7 +971,7 @@ def is_cm_j_invariant(j, method='CremonaSutherland'): return (False, None) h = jpol.degree() - if method in ['exhaustive', 'old']: + if algorithm in ['exhaustive', 'old']: if h>100: raise NotImplementedError("CM data only available for class numbers up to 100") for d,f in cm_orders(h): @@ -971,10 +979,10 @@ def is_cm_j_invariant(j, method='CremonaSutherland'): return (True, (d,f)) return (False, None) - if method not in ['reduction', 'new']: - raise ValueError("Invalid method {} in is_cm_j_invariant".format(method)) + if algorithm not in ['reduction', 'new']: + raise ValueError("Invalid algorithm {} in is_cm_j_invariant".format(algorithm)) - # Now we use the reduction method + # Now we use the reduction algorithm # If the degree h is less than the degree of j.parent() we recreate j as an element # of Q(j, and replace j by a clone whose parent is Q(j), if necessary: @@ -1024,7 +1032,7 @@ def is_cm_j_invariant(j, method='CremonaSutherland'): DP = aP**2 - 4*P.norm() dP = DP.squarefree_part() fP = ZZ(DP//dP).isqrt() - if cmd==0: # first one, so store d and f + if cmd == 0: # first one, so store d and f cmd = dP cmf = fP elif cmd != dP: # inconsistent with previous @@ -1032,8 +1040,8 @@ def is_cm_j_invariant(j, method='CremonaSutherland'): else: # consistent d, so update f cmf = cmf.gcd(fP) - if cmd==0: # no conclusion, we found no degree 1 primes, revert to default method - return is_cm_j_invariant(j, method='CremonaSutherland') + if cmd == 0: # no conclusion, we found no degree 1 primes, revert to default algorithm + return is_cm_j_invariant(j) # it looks like cm by disc cmd * f**2 where f divides cmf diff --git a/src/sage/schemes/elliptic_curves/ell_finite_field.py b/src/sage/schemes/elliptic_curves/ell_finite_field.py index 3ecaab7b068..4a1718a6d98 100644 --- a/src/sage/schemes/elliptic_curves/ell_finite_field.py +++ b/src/sage/schemes/elliptic_curves/ell_finite_field.py @@ -1434,13 +1434,14 @@ def height_above_floor(self, ell, e): INPUT: - ``ell`` -- a prime number - - ``e`` -- a non-negative integer + - ``e`` -- a non-negative integer, the ell-adic valuation of + the conductor the Frobenius order .. NOTE:: For an ordinary `E/\GF{q}`, and a prime `\ell`, the height - of the `\ell`-volcano containing `j(E)` is the `\ell`-adic + `e` of the `\ell`-volcano containing `j(E)` is the `\ell`-adic valuation of the conductor of the order generated by the Frobenius `\pi_E`; the height of `j(E)` on its ell-volcano is the `\ell`-adic valuation of the conductor @@ -1508,17 +1509,18 @@ def endomorphism_discriminant_from_class_number(self, h): (integer) The discriminant of the endomorphism ring `\End(E)`, if this has class number ``h``. If `\End(E)` does not have class - number ``h`` the output is undefined. + number ``h``, a ``ValueError`` is raised. ALGORITHM: - Compute the trace of Frobenius and hence the discriminant `D_0` - and class number `h_0` of the maximal order containing the - endomorphism order. From the given value of `h`, which must be a - multiple of `h_0`, compute the possible conductors, using - :meth:`height_above_floor` for each prime `\ell` dividing the - quotient `h/h_0`. If exactly one conductor `f` remains, return - `f^2D_0`, otherwise return `0`. + Compute the trace of Frobenius and hence the discriminant + `D_0` and class number `h_0` of the maximal order containing + the endomorphism order. From the given value of `h`, which + must be a multiple of `h_0`, compute the possible conductors, + using :meth:`height_above_floor` for each prime `\ell` + dividing the quotient `h/h_0`. If exactly one conductor `f` + remains, return `f^2D_0`, otherwise raise a ``ValueError``; + this can onlyhappen when the input value of `h` was incorrect. .. NOTE:: @@ -1553,24 +1555,23 @@ def endomorphism_discriminant_from_class_number(self, h): v = ZZ(D1//D0).isqrt() h0 = D0.class_number() if h%h0: - return 0 + raise ValueError("Incorrect class number {}".format(h)) from sage.schemes.elliptic_curves.cm import OrderClassNumber cs = [v//f for f in v.divisors() if OrderClassNumber(D0,h0,f) == h] # cofactors c=v/f compatible with h(f**2D0)=h if not cs: - return 0 + raise ValueError("Incorrect class number {}".format(h)) if len(cs) == 1: return (v//cs[0])**2 * D0 from sage.sets.set import Set - L = sorted(list(Set(sum([c.prime_factors() for c in cs], [])))) + L = sorted(set(sum([c.prime_factors() for c in cs], []))) for ell in L: e = self.height_above_floor(ell,v.valuation(ell)) cs = [c for c in cs if c.valuation(ell) == e] if not cs: - return 0 + raise ValueError("Incorrect class number {}".format(h)) if len(cs) == 1: return (v//cs[0])**2 * D0 - return 0 - + raise ValueError("Incorrect class number {}".format(h)) # dict to hold precomputed coefficient vectors of supersingular j values (excluding 0, 1728): From ccea429cc46311519f3c16ffbc1f7ece014ed006 Mon Sep 17 00:00:00 2001 From: John Cremona Date: Mon, 6 Mar 2023 12:02:51 +0000 Subject: [PATCH 11/12] fix 2 imports --- src/sage/schemes/elliptic_curves/cm.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/sage/schemes/elliptic_curves/cm.py b/src/sage/schemes/elliptic_curves/cm.py index 1d29409f5db..ffb26b13b40 100644 --- a/src/sage/schemes/elliptic_curves/cm.py +++ b/src/sage/schemes/elliptic_curves/cm.py @@ -237,6 +237,9 @@ def is_HCP(f, check_monic_irreducible=True): except AttributeError: return zero + from sage.rings.real_mpfr import RR + from sage.rings.finite_rings.finite_field_constructor import GF + h = f.degree() h2list = [d for d in h.divisors() if (d-h)%2 == 0 and d.prime_to_m_part(2) == 1] pmin = 33 * (h**2 * (RR(h+2).log().log()+2)**2).ceil() From d83aa3a5043c05deb6eeecb4d9d0b601425a05d8 Mon Sep 17 00:00:00 2001 From: John Cremona Date: Mon, 6 Mar 2023 16:12:50 +0000 Subject: [PATCH 12/12] use sets and defaultdict --- src/sage/schemes/elliptic_curves/cm.py | 21 +++++++------------ .../elliptic_curves/ell_finite_field.py | 6 +++--- 2 files changed, 11 insertions(+), 16 deletions(-) diff --git a/src/sage/schemes/elliptic_curves/cm.py b/src/sage/schemes/elliptic_curves/cm.py index ffb26b13b40..cb6a1ef04df 100644 --- a/src/sage/schemes/elliptic_curves/cm.py +++ b/src/sage/schemes/elliptic_curves/cm.py @@ -784,12 +784,6 @@ def discriminants_with_bounded_class_number(hmax, B=None, proof=None): # global hDf_dict, only including those (D,f) for which |D|*f**2 # <= B if B was provided: - T = hDf_dict.copy() - - if not count: - for h in T: - T[h] = [Df for Df in T[h] if Df[0].abs()*Df[1]**2<=B] - # h_dict caches the class number h of all discriminants previously # encountered; we will use the function OrderClassNumber() to # quicky compute the class number of non-fundamental discriminants @@ -797,10 +791,15 @@ def discriminants_with_bounded_class_number(hmax, B=None, proof=None): # keys of h_dict include nonfundamental discriminants, but we only # update it with fundamental ones. + from collections import defaultdict + T = defaultdict(set) h_dict = {} for h, Dflist in hDf_dict.items(): for D0,f in Dflist: h_dict[D0*f**2] = h + if not count: + Dflist = [Df for Df in Dflist if Df[0].abs()*Df[1]**2<=B] + T[h] = set(Dflist) # We do not need to certify the class number from :pari:`qfbclassno` for discriminants under 2*10^10 if B < 2*10**10: @@ -827,16 +826,12 @@ def discriminants_with_bounded_class_number(hmax, B=None, proof=None): # If the class number of this order is within the range, then store (D0,f) if h <= hmax: - z = (D0, f) - if h in T: - if z not in T[h]: - T[h].append(z) - else: - T[h] = [z] + T[h].add((D0,f)) # sort each list of (D,f) pairs by (|D|,f) for h in T: + T[h] = list(T[h]) T[h].sort(key=lambda Df: (Df[0].abs(), Df[1])) # count is None precisely when the user provided a value of B @@ -846,7 +841,7 @@ def discriminants_with_bounded_class_number(hmax, B=None, proof=None): if len(T[h]) != count[h]: raise RuntimeError("number of discriminants inconsistent with Watkins's table") # 2. Update the global dict - hDf_dict.update(T) + hDf_dict.update(dict(T)) return T diff --git a/src/sage/schemes/elliptic_curves/ell_finite_field.py b/src/sage/schemes/elliptic_curves/ell_finite_field.py index f23d6acf217..5b648cf601d 100644 --- a/src/sage/schemes/elliptic_curves/ell_finite_field.py +++ b/src/sage/schemes/elliptic_curves/ell_finite_field.py @@ -1453,7 +1453,7 @@ def height_above_floor(self, ell, e): valuation of the conductor of the order generated by the Frobenius `\pi_E`; the height of `j(E)` on its ell-volcano is the `\ell`-adic valuation of the conductor - of the order `\End(E)`. + of the order `\text{End}(E)`. ALGORITHM: @@ -1515,8 +1515,8 @@ def endomorphism_discriminant_from_class_number(self, h): OUTPUT: - (integer) The discriminant of the endomorphism ring `\End(E)`, if - this has class number ``h``. If `\End(E)` does not have class + (integer) The discriminant of the endomorphism ring `\text{End}(E)`, if + this has class number ``h``. If `\text{End}(E)` does not have class number ``h``, a ``ValueError`` is raised. ALGORITHM: