From 0dae4eeaa26f207c0e36589781548ed05de83b66 Mon Sep 17 00:00:00 2001 From: David Roe Date: Sun, 19 Jul 2026 14:33:20 -0400 Subject: [PATCH 1/3] Identify abelian number fields without polredabs in polynomial lookup (LMFDB#5471) Add a fast path to WebNumberField.from_polynomial implementing jwj61's abelian identification: after polredbest, certify the field abelian with galoisinit over an order maximal at the small primes of the discriminant (never factoring the discriminant), compute the field discriminant from that order, query nf_fields by degree/signature/discriminant, and confirm candidates by exhibiting a verified root of their defining polynomial, so any label returned is provably correct; all failures fall back to the unchanged polredabs path. Inputs whose degree exceeds every field in nf_fields now return "not in database" immediately instead of running polredabs (a degree 94 kernel polynomial went from >300s to 0.15s; the issue's degree 47 field is identified from a non-reduced polynomial in ~2s). Verified: all 36 number_fields tests pass (3 new: identified jumps for the issue's degree 47 field and Q(zeta_32) from non-reduced polynomials, abelian_nf_label unit checks, fast large-degree jump); pyflakes clean. Co-Authored-By: Claude Fable 5 --- lmfdb/number_fields/test_numberfield.py | 44 ++++++++ lmfdb/number_fields/web_number_field.py | 135 +++++++++++++++++++++++- 2 files changed, 177 insertions(+), 2 deletions(-) diff --git a/lmfdb/number_fields/test_numberfield.py b/lmfdb/number_fields/test_numberfield.py index 850f340559..5a72235598 100644 --- a/lmfdb/number_fields/test_numberfield.py +++ b/lmfdb/number_fields/test_numberfield.py @@ -46,6 +46,50 @@ 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 + + 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 e9ab44c86b..a59e9e55bc 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, @@ -510,6 +512,125 @@ 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 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 + # Primes up to 10^5 dividing D + S = D.gcd(_prime_product(10**5)).prime_divisors() + # If the remaining unfactored part of D is a pseudoprime of + # reasonable size, use it too: this catches fields ramified at one + # larger prime. If it is a hard composite we leave it alone. + C = D.abs() + for p in S: + C //= p**C.valuation(p) + if C > 1 and C.ndigits() <= 300 and C.is_pseudoprime(): + S.append(C) + # 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]) + 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: + for rt in pari.nfroots(nf, gpol): + 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'] + except PariError: + continue + return None + + class WebNumberField: """ Class for retrieving number field information from the database @@ -563,8 +684,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 From fd85f3dc831936a468ca00d5a54a50ff3d3e585b Mon Sep 17 00:00:00 2001 From: David Roe Date: Wed, 5 Aug 2026 03:05:54 -0400 Subject: [PATCH 2/3] Harden the abelian fast path: prime-power cofactors and nfroots input (LMFDB#5471) Two review fixes, both of which could make the abelian fast path return None and fall back to the expensive polredabs path (neither could produce a wrong label, which the exact root check still rules out). Recognize a large ramified prime that appears in the discriminant as a prime power: in a field of degree at least 4 a ramified p contributes p^e with e > 1, and the index contributes further powers, so the cofactor left after dividing out the small primes is typically p^e rather than p and is_pseudoprime() rejected it, leaving p out of S and making DK (hence the database query) wrong. Use is_pseudoprime_power(get_data=True) and append the base; a cofactor with several large primes is still left unfactored. The selection of these primes moves to a helper _known_discriminant_primes. Give nfroots the defining polynomial rather than the conditional nf structure: nfinit([T, S]) is only certified maximal at the primes of S, and the PARI manual warns that nfroots may miss a root when handed such a structure, while it recovers in polynomial time from nf.pol. The exact substitution certificate and the polredabs fallback are unchanged. Verified: all 38 number_fields tests pass, including two new regressions (a discriminant whose large ramified prime occurs as p^7, and an abelian field entered with an index divisible by two primes above 10^5, so that nfroots runs against a conditional order); pyflakes and ruff clean. On the degree 47 conductor 283 example nfroots goes from 0.37s to 0.78s and the whole fast path takes 1.5s, against more than 300s for polredabs. Co-Authored-By: Claude Opus 5 --- lmfdb/number_fields/test_numberfield.py | 18 +++++++++++ lmfdb/number_fields/web_number_field.py | 43 ++++++++++++++++++------- 2 files changed, 50 insertions(+), 11 deletions(-) diff --git a/lmfdb/number_fields/test_numberfield.py b/lmfdb/number_fields/test_numberfield.py index 5bb591310e..bea90dd006 100644 --- a/lmfdb/number_fields/test_numberfield.py +++ b/lmfdb/number_fields/test_numberfield.py @@ -76,6 +76,24 @@ def test_abelian_nf_label(self): # 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 diff --git a/lmfdb/number_fields/web_number_field.py b/lmfdb/number_fields/web_number_field.py index cae4b73d5a..de90140e92 100644 --- a/lmfdb/number_fields/web_number_field.py +++ b/lmfdb/number_fields/web_number_field.py @@ -568,6 +568,30 @@ def _probably_galois(T, needed=6, maxp=1000): 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 @@ -601,19 +625,16 @@ def abelian_nf_label(T): D = ZZ(T.poldisc()) if D == 0: return None - # Primes up to 10^5 dividing D - S = D.gcd(_prime_product(10**5)).prime_divisors() - # If the remaining unfactored part of D is a pseudoprime of - # reasonable size, use it too: this catches fields ramified at one - # larger prime. If it is a hard composite we leave it alone. - C = D.abs() - for p in S: - C //= p**C.valuation(p) - if C > 1 and C.ndigits() <= 300 and C.is_pseudoprime(): - S.append(C) + 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, used in place of nf when looking for + # roots below: 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 + field_pol = nf.getattr("pol") gal = pari.galoisinit(nf) if gal == 0 or pari.galoisisabelian(gal) == 0: # not certified Galois, or certified non-abelian @@ -633,7 +654,7 @@ def abelian_nf_label(T): for cand in db.nf_fields.search(query, ['label', 'coeffs']): gpol = pari(coeff_to_poly(cand['coeffs'])) try: - for rt in pari.nfroots(nf, gpol): + for rt in pari.nfroots(field_pol, gpol): 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 From 23f7ab000f2783e8186cd433a57297f1d984fbcf Mon Sep 17 00:00:00 2001 From: David Roe Date: Wed, 5 Aug 2026 03:33:24 -0400 Subject: [PATCH 3/3] Only pay for the conditional-order nfroots retry when it is needed (LMFDB#5471) Passing nf.pol to nfroots for every candidate protects against the root that a conditional nf structure can miss, but it also slows down the candidates that succeed, which is the case the fast path exists for: on the degree 47 conductor 283 example nfroots went from 0.37s to 0.78s. A missed root can only show up as an empty result, so try the nf form first and ask again with the defining polynomial only when nothing came back, which is exactly the situation that would otherwise fall through to polredabs. The degree 47 lookup is back to 1.18s end to end, and a candidate that really is not isomorphic to K costs a second nfroots call that returns immediately (0.00s at degree 8 and 16). Verified: all 38 number_fields tests pass; pyflakes and ruff clean. Co-Authored-By: Claude Opus 5 --- lmfdb/number_fields/web_number_field.py | 34 ++++++++++++++++--------- 1 file changed, 22 insertions(+), 12 deletions(-) diff --git a/lmfdb/number_fields/web_number_field.py b/lmfdb/number_fields/web_number_field.py index de90140e92..0e3ea20423 100644 --- a/lmfdb/number_fields/web_number_field.py +++ b/lmfdb/number_fields/web_number_field.py @@ -629,11 +629,8 @@ def abelian_nf_label(T): # 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, used in place of nf when looking for - # roots below: 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 + # 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: @@ -654,14 +651,27 @@ def abelian_nf_label(T): for cand in db.nf_fields.search(query, ['label', 'coeffs']): gpol = pari(coeff_to_poly(cand['coeffs'])) try: - for rt in pari.nfroots(field_pol, gpol): - 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'] + roots = pari.nfroots(nf, gpol) except PariError: - continue + 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