diff --git a/lmfdb/number_fields/test_numberfield.py b/lmfdb/number_fields/test_numberfield.py index 1c22946781..bea90dd006 100644 --- a/lmfdb/number_fields/test_numberfield.py +++ b/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 diff --git a/lmfdb/number_fields/web_number_field.py b/lmfdb/number_fields/web_number_field.py index 0751ba8bfa..0e3ea20423 100644 --- a/lmfdb/number_fields/web_number_field.py +++ b/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