Identify abelian number fields in the jump box without polredabs

lmfdb/number_fields/test_numberfield.py

@@ -46,6 +46,68 @@ def test_search_multiple_fields(self):
         self.check_args('/NumberField/?jump=Qsqrt5%2c+x%5E2-3&search=Go', '2.2.5.1')
         self.check_args('/NumberField/?jump=Qsqrt5%2c+x%5E2-3&search=Go', '2.2.12.1')
 
+    def test_search_abelian_jump(self):
+        # Abelian fields entered by a non-reduced defining polynomial are
+        # identified without running polredabs (issue #5471).
+        from sage.all import pari
+        from urllib.parse import quote
+        # the degree 47 field of conductor 283 from the issue, entered via
+        # the minimal polynomial of z + z^2 for z a root of the stored
+        # polynomial
+        label = "47.47.60558628944427886416035618894711378994697503545758179730765967261479053047453845062877188530544503628007178201769.1"
+        coeffs = self.db.nf_fields.lookup(label, "coeffs")
+        g = pari([int(c) for c in coeffs]).Polrev()
+        T = pari("x + x^2").Mod(g).charpoly()
+        self.check_args('/NumberField/?jump=' + quote(str(T)),
+                        [label, 'uses a different defining polynomial'])
+        # same for a moderate degree: Q(zeta_32), degree 16
+        T = pari("x + x^2").Mod(pari("polcyclo(32)")).charpoly()
+        self.check_args('/NumberField/?jump=' + quote(str(T)),
+                        '16.0.18446744073709551616.1')
+
+    def test_abelian_nf_label(self):
+        # the underlying fast path for issue #5471
+        from sage.all import pari
+        from lmfdb.number_fields.web_number_field import abelian_nf_label
+        # degree 8: Q(zeta_20), entered via the minimal polynomial of z + 3z^3
+        T = (pari("x") + 3 * pari("x^3")).Mod(pari("polcyclo(20)")).charpoly()
+        assert abelian_nf_label(T) == "8.0.4000000.1"
+        # non-Galois and Galois-but-non-abelian inputs are left to the
+        # polredabs path
+        assert abelian_nf_label(pari("x^8 - 2")) is None
+        assert abelian_nf_label(pari("polcompositum(x^4 - 2, x^2 + 1)[1]")) is None
+        # same field, but entered so that the order we can certify maximal is
+        # very far from maximal: the index is divisible by two primes above
+        # 10^5, which stay out of S, so nfroots is called with a conditional
+        # structure (hence gets the defining polynomial, not the nf)
+        m = 100003 * 100019
+        T = (m * pari("x")).Mod(pari("polcyclo(20)")).charpoly()
+        assert abelian_nf_label(T) == "8.0.4000000.1"
+
+    def test_known_discriminant_primes(self):
+        # a large ramified prime shows up in the discriminant as a prime
+        # power, not as a prime (issue #5471)
+        from sage.all import ZZ
+        from lmfdb.number_fields.web_number_field import _known_discriminant_primes
+        q = ZZ(100003)
+        S = _known_discriminant_primes(ZZ(2)**20 * ZZ(5)**10 * q**7)
+        assert S == [ZZ(2), ZZ(5), q]
+        # a cofactor with two large prime factors is left unfactored
+        assert _known_discriminant_primes(ZZ(2)**20 * q * ZZ(100019)) == [ZZ(2)]
+
+    def test_jump_degree_too_large(self):
+        # for degrees beyond anything in the database the jump returns
+        # quickly instead of attempting polredabs (issue #5471): here a
+        # degree 94 subfield of Q(zeta_283), for which polredabs takes
+        # more than five minutes
+        from sage.all import pari
+        from urllib.parse import quote
+        T = pari.polsubcyclo(283, 94)
+        if T.type() == 't_VEC':
+            T = T[0]
+        self.check_args('/NumberField/?jump=' + quote(str(T)),
+                        'does not define a number field in the database')
+
     def test_search_disc(self):
         self.check_args('/NumberField/?discriminant=1988-2014', '401') # factor of one of the discriminants
 

lmfdb/number_fields/web_number_field.py

@@ -3,9 +3,11 @@
 import yaml
 
 from flask import url_for
+from cypari2 import PariError
 from sage.all import (
     Set, ZZ, RR, pi, gcd, euler_phi, CyclotomicField, gap, RealField, sqrt, prod,
-    QQ, NumberField, QuadraticField, PolynomialRing, latex, pari, cached_function, Permutation)
+    QQ, NumberField, QuadraticField, PolynomialRing, latex, pari, cached_function,
+    Permutation, prime_range)
 
 from lmfdb import db
 from lmfdb.utils import (web_latex, coeff_to_poly,
@@ -523,6 +525,156 @@ def get_local_field(lab):
     return LF
 
 
+@cached_function
+def max_nf_degree():
+    return db.nf_fields.max('degree')
+
+
+# Minimal degree for which we attempt the abelian lookup below; for
+# smaller degrees polredabs is cheap anyway
+ABELIAN_LOOKUP_MIN_DEGREE = 8
+
+
+@cached_function
+def _prime_product(bound):
+    # The product of all primes up to bound, used to extract the small
+    # prime factors of a discriminant with a single gcd (much faster
+    # than trial division when the discriminant is huge)
+    return prod(prime_range(bound))
+
+
+def _probably_galois(T, needed=6, maxp=1000):
+    """
+    Cheap necessary condition for the irreducible pari polynomial ``T`` to
+    define a Galois number field: modulo any prime not dividing its
+    discriminant, all irreducible factors have the same degree.  Returns
+    False only if ``T`` is provably not Galois (hence not abelian); True
+    means "maybe Galois".
+    """
+    count = 0
+    lead = ZZ(T.pollead())
+    for p in prime_range(maxp):
+        if lead % p == 0:
+            continue
+        fm = T.factormod(p)
+        if any(int(e) > 1 for e in fm[1]):
+            # p divides the discriminant of T
+            continue
+        if len({int(f.poldegree()) for f in fm[0]}) > 1:
+            return False
+        count += 1
+        if count >= needed:
+            break
+    return True
+
+
+def _known_discriminant_primes(D):
+    """
+    The prime divisors of the nonzero integer ``D`` that can be found
+    cheaply: those below 10^5, together with the base of the remaining
+    cofactor when that cofactor is a power of a single pseudoprime of
+    reasonable size.  A hard composite cofactor is left alone; ``D`` is
+    never fully factored.
+    """
+    S = D.gcd(_prime_product(10**5)).prime_divisors()
+    C = D.abs()
+    for p in S:
+        C //= p**C.valuation(p)
+    if C > 1 and C.ndigits() <= 300:
+        # This catches a field ramified at one larger prime p.  Such a p
+        # occurs in the discriminant with exponent greater than one as soon
+        # as the degree is at least 4, and the index of Z[x]/(T) contributes
+        # further powers of p, so the cofactor is a prime power rather than
+        # a prime.
+        p, e = C.is_pseudoprime_power(get_data=True)
+        if e:
+            S.append(p)
+    return S
+
+
+def abelian_nf_label(T):
+    """
+    Attempt to find the label of the number field K defined by ``T`` (an
+    integral irreducible pari polynomial in x, typically the output of
+    polredbest) without running polredabs, using the strategy suggested by
+    jwj61 for abelian fields (see issue #5471):
+
+    - certify that K is abelian, using galoisinit on an order that is
+      maximal at the "known" primes of disc(T); no attempt is ever made to
+      fully factor the discriminant, which can be infeasible;
+    - read off the field discriminant from the valuations of the
+      discriminant of that order at the known primes (correct as soon as
+      they include all ramified primes; if not, the database query below
+      simply comes back empty);
+    - query nf_fields by degree, signature and discriminant, and confirm
+      each candidate with nfroots, by exhibiting a root of its defining
+      polynomial in K.
+
+    Returns the label, or None (not applicable / not certified abelian /
+    no confirmed match), in which case the caller should fall back to the
+    polredabs path.  A non-None answer is always correct: the root of the
+    candidate's defining polynomial giving the isomorphism is verified by
+    an exact polynomial computation.
+    """
+    n = int(T.poldegree())
+    if n < ABELIAN_LOOKUP_MIN_DEGREE:
+        return None
+    try:
+        if not _probably_galois(T):
+            return None
+        D = ZZ(T.poldisc())
+        if D == 0:
+            return None
+        S = _known_discriminant_primes(D)
+        # Order maximal (at least) at the primes in S; written in the
+        # variable y so that we can factor polynomials in x over K below
+        nf = pari.nfinit([T.subst("x", "y"), S])
+        # The defining polynomial of K, kept for the second attempt at root
+        # finding below
+        field_pol = nf.getattr("pol")
+        gal = pari.galoisinit(nf)
+        if gal == 0 or pari.galoisisabelian(gal) == 0:
+            # not certified Galois, or certified non-abelian
+            return None
+        # The order is p-maximal for every p in S, so the valuations of its
+        # discriminant at those primes are those of disc(K); primes outside
+        # S are ignored (they are index primes, unless one of them ramifies,
+        # in which case DK is wrong at it and the query just finds no match)
+        d_ord = ZZ(nf.disc())
+        DK = prod(p**d_ord.valuation(p) for p in S)
+        r1 = int(T.polsturm())
+    except PariError:
+        return None
+    r2 = (n - r1) // 2
+    query = {'degree': n, 'r2': r2, 'disc_abs': int(DK),
+             'disc_sign': 1 if r2 % 2 == 0 else -1}
+    for cand in db.nf_fields.search(query, ['label', 'coeffs']):
+        gpol = pari(coeff_to_poly(cand['coeffs']))
+        try:
+            roots = pari.nfroots(nf, gpol)
+        except PariError:
+            roots = []
+        if len(roots) == 0:
+            # nf is only conditional (its order is certified maximal at the
+            # primes of S and nowhere else), and PARI warns that nfroots can
+            # miss a root when handed such a structure, while it recovers in
+            # polynomial time when handed nf.pol.  Ask again that way before
+            # discarding the candidate: this costs nothing when the first
+            # call already found a root, and the alternative to a second try
+            # is the polredabs path, which is orders of magnitude slower
+            try:
+                roots = pari.nfroots(field_pol, gpol)
+            except PariError:
+                continue
+        for rt in roots:
+            if gpol.subst("x", rt) == 0:
+                # certified: K contains a root of gpol, which is irreducible
+                # of the same degree n, so K is isomorphic to the field of
+                # this candidate
+                return cand['label']
+    return None
+
+
 class WebNumberField:
     """
      Class for retrieving number field information from the database
@@ -576,8 +728,18 @@ def from_polynomial(cls, pol):
         # For some reason the error raised by Pari on a constant polynomial is not being caught
         if pol.degree() < 1:
             raise ValueError("Polynomial cannot be constant")
+        if pol.degree() > max_nf_degree():
+            # there is no field of this degree in the database, so we can
+            # skip the (potentially very expensive) canonicalization below
+            return cls('a')  # will initialize data to None
         R = pol.parent()
-        pol = R(pari(pol).polredbest().polredabs())
+        pol = pari(pol).polredbest()
+        # For abelian fields we may be able to identify the field without
+        # running polredabs, which can be very expensive (issue #5471)
+        label = abelian_nf_label(pol)
+        if label is not None:
+            return cls(label)
+        pol = R(pol.polredabs())
         return cls.from_coeffs([int(c) for c in pol.coefficients(sparse=False)])
 
     # If we already have the database entry