From b6219b9d8c3df3c24926f6791711d720ccc57d77 Mon Sep 17 00:00:00 2001 From: David Roe Date: Sun, 9 Aug 2026 16:20:58 -0400 Subject: [PATCH] Shape matchers: recognizing families in a parsed equation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The shape matchers: `matchers.run(parsed)` returns `Match(slug, data)` records. This is the file where a human eye is most useful on the mathematics: the degree/genus tests, the reductions to standard forms, and the data each match extracts (`D` for Pell, a-invariants for Weierstrass, the regime of a generalized Fermat equation, ...). Matchers never mutate the parsed equation. **API note:** `Match.transform` is a `CoordinateTransform` (`transforms.py`), not a prose string. It holds both directions of the map between the user's variables and the family's standard coordinates, `push_forward`/`pull_back`, and the structural roles a family assigns — which have moved out of `Match.data["roles"]`, since a solver needs them together with the map that consumes them. A string cannot transport a solution: `5*x^2 - y^2 = 1` is a Pell equation only after reading the user's `y` as the standard `x`, and the answer has to come back. Operations that change the equation but no coordinate (multiplying through by `-1`) are recorded separately, in `.operations`. Two recognition fixes worth a look: a parameter in a coefficient is no longer dropped (a `Term` keeps its parameter factors apart from `coeff`, so `A*2^n` and `2^n` shared `coeff == 1` and Pillai matched both), and orienting `D*x^2 - y^2 = N` as a Pell equation swaps the variables, which negates the right-hand side too — without that, `5*x^2 - y^2 = 1` was matched as the different equation `x^2 - 5*y^2 = 1`. The recognizers cover 44 families, most of which are not registered yet. That is deliberate and safe: a recognizer whose family has not landed is inert, because the classifier (next PR) ranks matches through the registry and drops slugs it does not know. `tests/test_registry.py` checks the flags of the registered ones; `99-polish` tightens that to *every* emitted slug once the registry is complete. Part of the series that splits #1 into reviewable pieces: 1. `01-bibliography` — packaging, docs, annotated bibliography 2. `02-parsing` — equation strings to a term model 3. `03-registry` — the YAML family registry (3 seed families) 4. `04-matchers` — shape recognizers 5. `05-classify` — the classification pipeline 6. `06-solvers` — solver framework, two seed solvers, and the CLI 7. `07..09-backbone` — the 23 parent families of the DAG, by depth 8. one PR per remaining family (38 of them, mutually independent) 9. `99-polish` — restore the full doctests and tighten the invariants --- CLAUDE.md | 4 +- README.md | 14 +- diophantine_classifier/matchers.py | 1384 ++++++++++++++++++++++++++ diophantine_classifier/transforms.py | 566 +++++++++++ tests/test_matchers.py | 116 +++ tests/test_registry.py | 18 + tests/test_transforms.py | 272 +++++ 7 files changed, 2366 insertions(+), 8 deletions(-) create mode 100644 diophantine_classifier/matchers.py create mode 100644 diophantine_classifier/transforms.py create mode 100644 tests/test_matchers.py create mode 100644 tests/test_transforms.py diff --git a/CLAUDE.md b/CLAUDE.md index 26ef3fd..baa5f8c 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -33,8 +33,8 @@ are specified here but land in a later PR of the series. - `data/references.bib` — bibliography; `references.py` parses/formats it; `tools/check_references.py` validates it (and local PDFs in `references/pdf/.pdf`). -- `matchers.py` *(later in the series)* — shape recognizers emitting `Match(slug, data, transform)`. -- `transforms.py` *(later in the series)* — `CoordinateTransform`: the executable, invertible map +- `matchers.py` — shape recognizers emitting `Match(slug, data, transform)`. +- `transforms.py` — `CoordinateTransform`: the executable, invertible map from the user's variables to a family's standard coordinates, plus the structural roles. Solvers work normalized and `pull_back()`. - `classify.py` *(later in the series)* — factor-split, run matchers, rank by DAG depth (most diff --git a/README.md b/README.md index b636129..44343ef 100644 --- a/README.md +++ b/README.md @@ -17,10 +17,10 @@ standalone library. > **Status of this branch.** This is one PR of a stacked series that builds the > classifier layer by layer; the description above is where the series lands. -> What runs *here* is the parser, the family registry and the bibliography — -> the examples below all work on this branch. Matchers, classification, -> solvers and the remaining 61 families arrive in the later PRs of the series; -> the closing PR restores the full README. +> What runs *here* is the parser, the family registry, the bibliography and the +> structural matchers — the examples below all work on this branch. The +> classification pipeline, the solvers and the remaining 61 families arrive in +> the later PRs of the series; the closing PR restores the full README. ## Quick start @@ -91,11 +91,13 @@ reports the most specific match and the full lineage. - **The bibliography and its pipeline**: BibTeX parsing and display formatting, plus `tools/check_references.py` and its monotone verification ledger. +- **Structural matchers**: `matchers.run()` recognizes the shapes of ~45 + families and extracts each one's data. They are inert for families the + registry does not have yet: the classifier ranks matches through the + registry, which is what lets the families land one PR at a time. ## Coming in the rest of the series -- **Structural matchers** for ~45 families, and a genus-based geometry - fallback for irreducible plane curves. - **Classification**: reducible equations split into components, matches are ranked by depth in the family DAG, and `explain()` / `as_dict()` produce the human report and the JSON contract for the website backend. diff --git a/diophantine_classifier/matchers.py b/diophantine_classifier/matchers.py new file mode 100644 index 0000000..5779a68 --- /dev/null +++ b/diophantine_classifier/matchers.py @@ -0,0 +1,1384 @@ +r""" +Structural matchers: recognize which families a parsed equation belongs to. + +Each matcher inspects a :class:`~diophantine_classifier.parsing.ParsedEquation` +and emits :class:`Match` objects (family slug + extracted data). :func:`run` +collects matches from every applicable matcher; ranking by DAG specificity +happens in :mod:`diophantine_classifier.classify`. + +Wave 1 recognizes equations presented in (or very near) a family's standard +form, plus genus-based routing for plane curves. Deeper normalization +(unimodular reduction, completing the square, Nagell's algorithm) is the +wave-2 transformation layer; see ``docs/DESIGN.md``. + +EXAMPLES:: + + sage: from diophantine_classifier.parsing import parse + sage: from diophantine_classifier.matchers import run + sage: [m.slug for m in run(parse("x^2 - 61*y^2 = 1"))][-2:] + ['pell-like', 'pell'] +""" + +from dataclasses import dataclass, field + +from sage.all import QQ, ZZ, Curve, EllipticCurve, prod + +from .transforms import CoordinateTransform, identity, negate, rename + +#: skip genus computations for inputs of absurdly large degree +MAX_GENUS_DEGREE = 20 + + +@dataclass +class Match: + r""" + One structural match: a family together with extracted data. + + ATTRIBUTES: + + - ``slug`` -- string; the matched family's registry slug. + - ``data`` -- dict of family-specific extracted data, with stringified + values (e.g. ``{"D": "61", "N": "1"}`` for a Pell match, a-invariants + for a Weierstrass match). + - ``summary`` -- string; one-line human-readable description with the + parameters filled in. + - ``transform`` -- a + :class:`~diophantine_classifier.transforms.CoordinateTransform`: the + executable, invertible map from the user's variables to the family's + standard coordinates, together with the structural roles the family + assigns them. The identity (with no operations) when the equation was + already in standard form. A solver works in the standard coordinates + and calls ``transform.pull_back`` to state its answer in the user's. + + EXAMPLES:: + + sage: from diophantine_classifier.matchers import run + sage: from diophantine_classifier.parsing import parse + sage: m = run(parse("3*x + 5*y = 1"))[-1] + sage: m.slug, m.data["b"] + ('linear', '1') + sage: m.transform.pull_back({"x": 1766319049, "y": 226153980}) + {'x': 1766319049, 'y': 226153980} + """ + slug: str + data: dict = field(default_factory=dict) + summary: str = "" + transform: CoordinateTransform = field( + default_factory=CoordinateTransform) + + def __repr__(self): + r""" + Terse representation. + + EXAMPLES:: + + sage: from diophantine_classifier.matchers import run + sage: from diophantine_classifier.parsing import parse + sage: run(parse("3*x + 5*y = 1"))[-1] + Match('linear') + """ + return f"Match({self.slug!r})" + + +def _sign_of(c): + r""" + Return -1/0/+1 for constant coefficients, ``None`` if parametric. + + INPUT: + + - ``c`` -- a coefficient: an element of ``QQ`` or of ``QQ[params]`` + + EXAMPLES:: + + sage: from diophantine_classifier.matchers import _sign_of + sage: from sage.all import QQ, PolynomialRing + sage: _sign_of(QQ(-3)), _sign_of(QQ(0)), _sign_of(QQ(2)) + (-1, 0, 1) + sage: R = PolynomialRing(QQ, "k") + sage: _sign_of(R.gen()) is None + True + """ + try: + q = QQ(c) + except (TypeError, ValueError): + return None + return 0 if q == 0 else (1 if q > 0 else -1) + + +def _unit_coefficient(t): + r""" + Whether a term's *full* scalar coefficient is ``+1`` or ``-1``. + + A :class:`~diophantine_classifier.parsing.Term` keeps its parameter + factors apart from ``coeff``, so ``A*2^n`` and ``2^n`` share + ``coeff == 1``. Named classical shapes — Pillai, Catalan, Fermat, + Ramanujan–Nagell — are statements about unit multipliers and must test + the whole coefficient, not just the rational part. + + INPUT: + + - ``t`` -- a :class:`~diophantine_classifier.parsing.Term` + + OUTPUT: boolean + + EXAMPLES:: + + sage: from diophantine_classifier.parsing import parse + sage: from diophantine_classifier.matchers import _unit_coefficient + sage: [_unit_coefficient(t) for t in parse("2^n - 3^m = 1").terms] + [True, True, True] + sage: pe = parse("A*2^n - 3^m = 1", params="A") + sage: [_unit_coefficient(t) for t in pe.exponential_terms()] + [False, True] + """ + return not t.param_powers and abs(t.coeff) == 1 + + +def _concrete_coefficients(terms): + r""" + Whether every term in ``terms`` has a parameter-free coefficient. + + Guards the matchers that read coefficients off a term list as plain + rationals (monicity, a fixed constant multiplier, a known sign); with a + parameter present those readings drop the parameter silently. + + INPUT: + + - ``terms`` -- iterable of :class:`~diophantine_classifier.parsing.Term` + + OUTPUT: boolean + + EXAMPLES:: + + sage: from diophantine_classifier.parsing import parse + sage: from diophantine_classifier.matchers import _concrete_coefficients + sage: _concrete_coefficients(parse("x^2 + 7 = 2^n").terms) + True + sage: _concrete_coefficients(parse("A*x^2 + 7 = 2^n", params="A").terms) + False + """ + return not any(t.param_powers for t in terms) + + +def _as_int(c): + r""" + Coerce a coefficient to a Sage integer, or return ``None``. + + EXAMPLES:: + + sage: from diophantine_classifier.matchers import _as_int + sage: from sage.all import QQ + sage: _as_int(QQ(7)) + 7 + sage: _as_int(QQ(1)/2) is None + True + """ + try: + return ZZ(c) + except (TypeError, ValueError): + return None + + +def _homogeneous_parts(P): + r""" + Decompose a polynomial into its homogeneous components. + + INPUT: + + - ``P`` -- multivariate polynomial + + OUTPUT: dict mapping total degree to the homogeneous component + + EXAMPLES:: + + sage: from diophantine_classifier.matchers import _homogeneous_parts + sage: R. = QQ[] + sage: parts = _homogeneous_parts(y^2 - x^3 + 2) + sage: sorted(parts) + [0, 2, 3] + sage: parts[3] + -x^3 + """ + parts = {} + for c, m in zip(P.coefficients(), P.monomials()): + d = m.degree() + parts[d] = parts.get(d, P.parent().zero()) + c * m + return parts + + +def _diagonal_scan(P, R): + r""" + Detect a diagonal shape: every nonconstant monomial a pure power. + + INPUT: + + - ``P`` -- multivariate polynomial + - ``R`` -- its parent ring + + OUTPUT: ``(entries, const)`` with ``entries`` a list of triples + ``(coefficient, variable name, exponent)``, one per variable used, and + ``const`` the constant term — or ``None`` if some monomial mixes + variables or a variable occurs with two different exponents + + EXAMPLES:: + + sage: from diophantine_classifier.matchers import _diagonal_scan + sage: R. = QQ[] + sage: entries, const = _diagonal_scan(x^3 + y^3 + z^3 - 42, R) + sage: [(str(v), e) for _, v, e in entries] + [('x', 3), ('y', 3), ('z', 3)] + sage: const + -42 + sage: _diagonal_scan(x*y + z, R) is None + True + """ + entries = [] + seen = set() + const = P.parent().base_ring().zero() + gens = R.gens() + for c, m in zip(P.coefficients(), P.monomials()): + degs = m.degrees() + nonzero = [(i, e) for i, e in enumerate(degs) if e] + if not nonzero: + const += c + continue + if len(nonzero) != 1: + return None + i, e = nonzero[0] + if i in seen: + return None + seen.add(i) + entries.append((c, str(gens[i]), ZZ(e))) + return entries, const + + +def _gram(P, R): + r""" + Gram matrix of a quadratic form, as a list of lists of rationals. + + EXAMPLES:: + + sage: from diophantine_classifier.matchers import _gram + sage: R. = QQ[] + sage: _gram(x^2 + 3*x*y - y^2, R) + [[1, 3/2], [3/2, -1]] + + Parametric coefficients stay in their own ring instead of being coerced + to `QQ`, which a parameter cannot survive:: + + sage: S. = QQ[] + sage: T. = S[] + sage: _gram(x^2 + y^2 - D*z^2, T)[2][2] + -D + """ + gens = R.gens() + k = len(gens) + base = R.base_ring() + G = [[base.zero()] * k for _ in range(k)] + for i in range(k): + G[i][i] = base(P.monomial_coefficient(gens[i] ** 2)) + for j in range(i + 1, k): + half = base(P.monomial_coefficient(gens[i] * gens[j])) / 2 + G[i][j] = G[j][i] = half + return G + + +def _plane_curve_genus(P): + r""" + Geometric genus of the plane curve ``P = 0``, or ``None``. + + Returns ``None`` for degrees above :data:`MAX_GENUS_DEGREE` or when the + curve construction fails (reducible input, etc.). + + EXAMPLES:: + + sage: from diophantine_classifier.matchers import _plane_curve_genus + sage: R. = QQ[] + sage: _plane_curve_genus(y^2 - x^5 + x - 1) + 2 + sage: _plane_curve_genus(y^2 - x^3) # cuspidal cubic + 0 + """ + if P.total_degree() > MAX_GENUS_DEGREE: + return None + try: + return ZZ(Curve(P).genus()) + except Exception: + return None + + +def _is_irreducible(P): + r""" + Whether ``P`` is irreducible (up to constants), or ``None`` on failure. + + EXAMPLES:: + + sage: from diophantine_classifier.matchers import _is_irreducible + sage: R. = QQ[] + sage: _is_irreducible(x^2 + y^2 - 1) + True + sage: _is_irreducible(x^2 - y^2) + False + """ + try: + fac = P.factor() + except Exception: + return None + nontrivial = [(g, e) for g, e in fac if g.degree() > 0] + return len(nontrivial) == 1 and nontrivial[0][1] == 1 + + +# -------------------------------------------------------------------------- +# polynomial equations +# -------------------------------------------------------------------------- + +def _match_polynomial(pe): + r""" + Matches for a purely polynomial equation. + + Always emits ``general-polynomial`` with basic invariants, then + dispatches on the number of unknowns. + + INPUT: + + - ``pe`` -- a polynomial :class:`~diophantine_classifier.parsing.ParsedEquation` + + OUTPUT: list of :class:`Match` + + EXAMPLES:: + + sage: from diophantine_classifier.parsing import parse + sage: from diophantine_classifier.matchers import _match_polynomial + sage: [m.slug for m in _match_polynomial(parse("3*x + 5*y = 1"))] + ['general-polynomial', 'linear'] + sage: [m.slug for m in _match_polynomial(parse("x^2 - 5*x + 6 = 0"))] + ['general-polynomial', 'univariate'] + """ + P, R = pe.poly, pe.poly_ring + k = R.ngens() + d = P.total_degree() + out = [Match( + "general-polynomial", + data={"variables": k, "degree": ZZ(d), + "homogeneous": P.is_homogeneous()}, + summary=f"polynomial equation in {k} unknowns of degree {d}", + )] + + if k == 1: + out.append(Match("univariate", data={"f": str(P)}, + summary="single-variable polynomial: solve by factoring")) + return out + + if d == 1: + coeffs = [P.monomial_coefficient(g) for g in R.gens()] + b = -P.constant_coefficient() + out.append(Match("linear", + data={"coeffs": [str(c) for c in coeffs], "b": str(b)}, + summary="linear Diophantine equation")) + return out + + if k == 2: + out.extend(_match_binary(pe, P, R, d)) + else: + out.extend(_match_multivar(pe, P, R, k, d)) + return out + + +def _match_binary(pe, P, R, d): + r""" + Matches for polynomial equations in two unknowns. + + Degree 2 goes through the binary-quadratic battery (Pell, sums of two + squares, representation by a form); degree ≥ 3 through binary forms + (Thue), curve shapes (Weierstrass, quartic, hyperelliptic, + superelliptic), then genus routing. + + INPUT: + + - ``pe`` -- the parsed equation; ``P`` its polynomial in ring ``R``; + ``d`` the total degree + + OUTPUT: list of :class:`Match` + + EXAMPLES:: + + sage: from diophantine_classifier.parsing import parse + sage: from diophantine_classifier.matchers import _match_binary + sage: pe = parse("x^2 - 61*y^2 = 5") + sage: [m.slug for m in _match_binary(pe, pe.poly, pe.poly_ring, 2)] + ['binary-quadratic', 'binary-qf-representation', 'pell-like'] + sage: pe = parse("x^3 + 2*y^3 = 11") + sage: [m.slug for m in _match_binary(pe, pe.poly, pe.poly_ring, 3)] + ['thue'] + """ + x, y = R.gens() + out = [] + if d == 2: + a = P.monomial_coefficient(x ** 2) + b = P.monomial_coefficient(x * y) + c = P.monomial_coefficient(y ** 2) + dd = P.monomial_coefficient(x) + ee = P.monomial_coefficient(y) + f0 = P.constant_coefficient() + disc = b ** 2 - 4 * a * c + out.append(Match( + "binary-quadratic", + data={"a": str(a), "b": str(b), "c": str(c), "d": str(dd), + "e": str(ee), "f": str(f0), "disc": str(disc)}, + summary=f"binary quadratic with discriminant {disc}", + )) + if dd == 0 and ee == 0 and f0 != 0: + n = -f0 + operations = () + swapped = False + # orient: prefer positive coefficient on the first square + if _sign_of(a) == -1: + a, b, c, n = -a, -b, -c, -n + operations = ("multiplied by -1",) + if b == 0 and _sign_of(a) == 1 and _sign_of(c) == -1 \ + and a != 1 and c == -1: + # D x^2 - y^2 = N: read it as y^2 - D x^2 = -N, which is the + # standard form with the variables swapped. Negating the + # right-hand side is part of the swap: dropping it turned + # 5x^2 - y^2 = 1 into the different equation x^2 - 5y^2 = 1 + a, c, n = -c, -a, -n + swapped = True + src = (str(x), str(y)) + transform = rename( + zip(("x", "y"), reversed(src) if swapped else src), + description="swapped variables" if swapped else "", + operations=operations) + out.append(Match( + "binary-qf-representation", + data={"a": str(a), "b": str(b), "c": str(c), "n": str(n), + "disc": str(disc)}, + summary=f"representation of {n} by the form ({a}, {b}, {c})", + transform=transform, + )) + if b == 0 and a == 1: + if c == 1: + out.append(Match("sum-of-two-squares", data={"n": str(n)}, + summary=f"x^2 + y^2 = {n}", + transform=transform)) + else: + D = -c + Dz = _as_int(D) + if Dz is not None and Dz > 0: + if Dz.is_square(): + out.append(Match( + "binary-form-reducible", + data={"factors": + f"(x - {Dz.sqrt()}*y)*(x + {Dz.sqrt()}*y)", + "m": str(n)}, + summary=f"x^2 - {Dz}y^2 factors (D is a " + "square): solve by divisor " + "enumeration", + transform=transform, + )) + else: + Nz = _as_int(n) + out.append(Match( + "pell-like", data={"D": str(Dz), "N": str(n)}, + summary=f"x^2 - {Dz}y^2 = {n}", + transform=transform, + )) + if Nz is not None and Nz in (1, -1): + out.append(Match( + "pell", data={"D": str(Dz), "N": str(Nz)}, + summary=f"Pell equation with D = {Dz}" + + (" (negative Pell)" + if Nz == -1 else ""), + transform=transform, + )) + elif Dz is None: + # parametric D + out.append(Match("pell-like", + data={"D": str(D), "N": str(n)}, + summary=f"x^2 - ({D})y^2 = {n}", + transform=transform)) + if n == 1: + out.append(Match("pell", + data={"D": str(D), "N": "1"}, + summary=f"Pell equation with " + f"D = {D}", + transform=transform)) + elif dd == 0 and ee == 0 and f0 == 0: + out.append(Match( + "binary-form-reducible", + data={"m": "0"}, + summary="homogeneous quadratic = 0: rational lines exist iff " + "the discriminant is a square", + )) + return out + + # degree >= 3 in two variables + parts = _homogeneous_parts(P) + c0 = P.constant_coefficient() + nonconst = sorted(deg for deg in parts if deg > 0) + if len(nonconst) == 1 and nonconst[0] == d: + # F_d(x, y) = m + F = parts[d] + m = -c0 + data = {"form": str(F), "m": str(m), "degree": ZZ(d), + "fx": str(F.subs({y: 1})), "x": str(x), "y": str(y)} + if not pe.is_concrete: + out.append(Match("binary-form", data=data, + summary=f"binary form of degree {d} = {m}")) + return out + fac = F.factor() + nontrivial = [(g, e) for g, e in fac if g.degree() > 0] + irreducible = len(nontrivial) == 1 and nontrivial[0][1] == 1 + if m == 0: + out.append(Match( + "binary-form-reducible", data=data, + summary="homogeneous form = 0: solutions come from rational " + "linear factors" + ("" if not irreducible else + " (none here: only (0,0))"), + )) + elif irreducible: + mz = _as_int(m) + out.append(Match( + "thue", data=dict(data, m=str(mz if mz is not None else m)), + summary=f"Thue equation of degree {d}", + )) + else: + data["factors"] = str(fac) + out.append(Match( + "binary-form-reducible", data=data, + summary="reducible binary form: enumerate factorizations of " + f"{m} across the factors", + )) + return out + + shape = _match_two_var_shapes(pe, P, R, d) + if shape: + out.extend(shape) + return out + + out.extend(_genus_route(pe, P, affine=True)) + return out + + +def _match_two_var_shapes(pe, P, R, d): + r""" + Weierstrass / quartic / hyperelliptic / superelliptic shapes. + + Two passes over both variable orderings: the quadratic-in-``v`` shapes + (Weierstrass, quartic, hyperelliptic) are tried for *both* orientations + before any superelliptic shape, so that e.g. ``x^2 = y^3 - k`` is + recognized as a Mordell equation rather than as ``y^3 = x^2 + k``. + + OUTPUT: list of :class:`Match`, or ``None`` when no shape applies (the + caller then falls through to genus routing) + + EXAMPLES:: + + sage: from diophantine_classifier.parsing import parse + sage: from diophantine_classifier.matchers import _match_two_var_shapes + sage: pe = parse("y^2 = x^3 - 2") + sage: [m.slug for m in + ....: _match_two_var_shapes(pe, pe.poly, pe.poly_ring, 3)] + ['elliptic-weierstrass', 'mordell'] + sage: pe = parse("y^3 = x^4 + 2") + sage: [m.slug for m in + ....: _match_two_var_shapes(pe, pe.poly, pe.poly_ring, 4)] + ['superelliptic'] + sage: pe = parse("y^2 = x^3") # cusp: shapes decline, genus routes + sage: _match_two_var_shapes(pe, pe.poly, pe.poly_ring, 3) is None + True + """ + orientations = (tuple(R.gens()), tuple(reversed(R.gens()))) + for v, u in orientations: + if P.degree(v) != 2: + continue + coeffs = [P.coefficient({v: j}) for j in range(3)] + A = coeffs[2] + if not A.is_constant(): + continue + B1, C = coeffs[1], coeffs[0] + if A == -1: + A, B1, C = 1, -B1, -C + if A != 1: + continue + # v^2 + B1*v + C = 0 <=> v^2 + B1*v = g(u), g = -C + g = -C + n = g.degree(u) + if B1.degree(u) <= 1 and n == 3 \ + and g.monomial_coefficient(u ** 3) == 1: + a1 = B1.monomial_coefficient(u) + a3 = B1.constant_coefficient() + a2 = g.monomial_coefficient(u ** 2) + a4 = g.monomial_coefficient(u) + a6 = g.constant_coefficient() + ainvs = [a1, a2, a3, a4, a6] + if pe.is_concrete: + try: + EllipticCurve(QQ, [QQ(t) for t in ainvs]) + except (ArithmeticError, TypeError, ValueError): + continue # singular: fall through to genus routing + data = {"ainvs": "[%s]" % ", ".join(str(t) for t in ainvs), + "magma_ainvs": "[%s]" % ", ".join(str(t) for t in ainvs), + "x": str(u), "y": str(v)} + curve_map = _curve_transform(u, v) + out = [Match("elliptic-weierstrass", data=data, + transform=curve_map, + summary=f"Weierstrass equation with " + f"a-invariants {data['ainvs']}")] + if a1 == 0 and a2 == 0 and a3 == 0 and a4 == 0: + out.append(Match("mordell", data={"k": str(a6), + "x": str(u), "y": str(v)}, + transform=curve_map, + summary=f"Mordell equation with k = {a6}")) + return out + if B1 == 0 and n == 4: + qu = g.univariate_polynomial() + if pe.is_concrete and qu.discriminant() == 0: + return None + qcoeffs = [qu[i] for i in range(5)] + return [Match( + "elliptic-quartic", + data={"q": str(g), "x": str(u), "y": str(v), + "qcoeffs": str(list(reversed(qcoeffs)))}, + transform=_curve_transform(u, v), + summary=f"genus-one quartic {v}^2 = {g}", + )] + if B1 == 0 and n >= 5: + qu = g.univariate_polynomial() + if pe.is_concrete and not qu.is_squarefree(): + return None + genus = (n - 1) // 2 if n % 2 else (n - 2) // 2 + slug = "genus-two" if genus == 2 else "hyperelliptic" + return [Match( + slug, + data={"f": str(g), "genus": ZZ(genus), "x": str(u), + "y": str(v), + "magma_coeffs": str([qu[i] for i in range(n + 1)])}, + transform=_curve_transform(u, v), + summary=f"hyperelliptic: {v}^2 = degree-{n} polynomial, " + f"genus {genus}", + )] + # second pass: superelliptic shapes v^m = f(u), m >= 3 + for v, u in orientations: + m = P.degree(v) + if m < 3: + continue + coeffs = [P.coefficient({v: j}) for j in range(m + 1)] + A = coeffs[m] + if not A.is_constant(): + continue + if any(c != 0 for c in coeffs[1:m]): + continue + C = coeffs[0] + if A == -1: + A, C = 1, -C + if A != 1: + continue + g = -C + n = g.degree(u) + if n >= 2: + qu = g.univariate_polynomial() + if pe.is_concrete and not qu.is_squarefree(): + continue # e.g. x^3 = y^2: let genus routing handle it + return [Match( + "superelliptic", + data={"m": ZZ(m), "f": str(g), "x": str(u), "y": str(v)}, + transform=_curve_transform(u, v), + summary=f"superelliptic: {v}^{m} = {g}", + )] + return None + + +def _curve_transform(u, v): + r""" + The coordinate map for a two-variable curve shape ``v^m = f(u)``. + + Both orientations of the pair are tried when matching, so the standard + ``x`` (the base) and ``y`` (the power) may be either of the user's + variables. That reversal has to be executable, not merely recorded in + the match data: a solver returns points ``(x, y)`` of the standard model + and the map carries them back. + + INPUT: + + - ``u``, ``v`` -- the ring generators playing the standard ``x`` and ``y`` + + OUTPUT: a + :class:`~diophantine_classifier.transforms.CoordinateTransform` + + EXAMPLES:: + + sage: from diophantine_classifier.parsing import parse + sage: from diophantine_classifier.matchers import _curve_transform + sage: R = parse("y^2 = x^3 - 2").poly_ring + sage: y, x = R.gens() # y appears first here + sage: t = _curve_transform(x, y) + sage: t.pull_back({"x": 3, "y": 5}) + {'x': 3, 'y': 5} + sage: t.roles + {'x': 'x', 'y': 'y'} + """ + reversed_roles = str(u) != "x" or str(v) != "y" + return rename( + [("x", str(u)), ("y", str(v))], + description=(f"read {u} as the base and {v} as the power" + if reversed_roles else "")) + + +def _genus_route(pe, P, affine=True): + r""" + Route an irreducible plane curve by genus. + + Mirrors the Library's classification plan: genus 0 → parametrize, + genus 1 → find a point then reduce to Weierstrass form, genus ≥ 2 → + Faltings finiteness and Chabauty-type methods. + + OUTPUT: list of at most one :class:`Match` (empty for parametric input, + reducible polynomials, or when the genus computation is skipped) + + EXAMPLES:: + + sage: from diophantine_classifier.parsing import parse + sage: from diophantine_classifier.matchers import _genus_route + sage: pe = parse("x^2*y^2 = x^3 + 1") + sage: [m.slug for m in _genus_route(pe, pe.poly)] + ['genus-one-curve'] + sage: pe = parse("y^2 = x^3") + sage: [m.slug for m in _genus_route(pe, pe.poly)] + ['genus-zero-curve'] + """ + if not pe.is_concrete: + return [] + if _is_irreducible(P) is not True: + return [] + g = _plane_curve_genus(P) + if g is None: + return [] + data = {"genus": g, "model": "affine" if affine else "projective"} + if g == 0: + return [Match("genus-zero-curve", data=data, + summary="plane curve of genus 0: parametrize")] + if g == 1: + return [Match("genus-one-curve", data=data, + summary="plane curve of genus 1: find a point, then " + "reduce to Weierstrass form")] + return [Match("general-curve", data=data, + summary=f"plane curve of genus {g}: Faltings finiteness; " + "Chabauty-type methods apply")] + + +def _match_multivar(pe, P, R, k, d): + r""" + Matches for polynomial equations in three or more unknowns. + + Handles quadratic forms (isotropy, representation, affine quadrics), + Markov–Hurwitz shapes, diagonal equations (generalized Fermat, sums of + cubes, Waring, equal sums of like powers), and plane projective curves. + + EXAMPLES:: + + sage: from diophantine_classifier.parsing import parse + sage: from diophantine_classifier.matchers import _match_multivar + sage: pe = parse("x^2 + y^2 = z^2") + sage: [m.slug for m in + ....: _match_multivar(pe, pe.poly, pe.poly_ring, 3, 2)] + ['quadratic-form-zero', 'legendre', 'pythagorean'] + sage: pe = parse("x^3 + y^3 + z^3 = 42") + sage: [m.slug for m in + ....: _match_multivar(pe, pe.poly, pe.poly_ring, 3, 3)] + ['sum-of-three-cubes'] + """ + out = [] + gens = R.gens() + parts = _homogeneous_parts(P) + c0 = P.constant_coefficient() + + if d == 2: + has_linear = 1 in parts + if not has_linear: + Q2 = parts.get(2, P.parent().zero()) + gram = _gram(Q2, R) + if c0 == 0: + out.append(Match("quadratic-form-zero", + data={"gram": str(gram), "k": k}, + summary=f"isotropy of a quadratic form in " + f"{k} variables (Hasse-Minkowski)")) + scan = _diagonal_scan(Q2, R) + if scan and k == 3: + entries, _ = scan + if len(entries) == 3: + coeffs = [c for c, _, _ in entries] + signs = [_sign_of(c) for c in coeffs] + if None not in signs: + if signs.count(1) == 1: + coeffs = [-c for c in coeffs] + signs = [-s for s in signs] + a, b, c = coeffs + out.append(Match( + "legendre", + data={"a": str(a), "b": str(b), "c": str(c)}, + summary=f"Legendre equation " + f"({a})x^2 + ({b})y^2 + ({c})z^2 = 0", + )) + if sorted(coeffs) == [-1, 1, 1]: + legs = [entries[i][1] for i in range(3) + if coeffs[i] == 1] + hyp = [entries[i][1] for i in range(3) + if coeffs[i] == -1][0] + out.append(Match( + "pythagorean", + transform=rename( + [("x", legs[0]), ("y", legs[1]), + ("z", hyp)], + roles={"legs": legs, + "hypotenuse": hyp}), + summary="Pythagorean equation: solutions " + "(m^2-n^2, 2mn, m^2+n^2)", + )) + else: + n = -c0 + out.append(Match("quadratic-form-representation", + data={"gram": str(gram), "n": str(n), + "k": k}, + summary=f"representation of {n} by a " + f"quadratic form in {k} variables")) + scan = _diagonal_scan(P, R) + if scan: + entries, _ = scan + if len(entries) == k and all(c == 1 for c, _, _ in entries): + if k == 3: + out.append(Match("sum-of-three-squares", + data={"n": str(n)}, + summary=f"three squares: " + f"n = {n}")) + elif k == 4: + out.append(Match("sum-of-four-squares", + data={"n": str(n)}, + summary=f"four squares: " + f"n = {n}")) + else: + out.append(Match("quadric", + data={"k": k}, + summary="affine quadric: complete the square, " + "then Hasse-Minkowski / " + "Grunewald-Segal")) + return out + + # --- degree >= 3 --- + + # Markov-Hurwitz: sum of squares = a * product of all variables + prod_mon = prod(gens) + expected = {g ** 2 for g in gens} | {prod_mon} + mons = set(P.monomials()) + if k >= 3 and mons == expected and c0 == 0: + sq = {P.monomial_coefficient(g ** 2) for g in gens} + pcoeff = P.monomial_coefficient(prod_mon) + if len(sq) == 1: + s = sq.pop() + if _sign_of(s) == -1: + s, pcoeff = -s, -pcoeff + a = -pcoeff + az = _as_int(a) + if s == 1 and az is not None and az > 0: + out.append(Match( + "markov-hurwitz", data={"a": str(az), "k": k}, + summary=("Markov equation" if (az == 3 and k == 3) else + f"Hurwitz equation x_1^2+...+x_{k}^2 = " + f"{az} x_1...x_{k}"), + )) + return out + + scan = _diagonal_scan(P, R) + if scan: + entries, _ = scan + if len(entries) == k: + exps = [e for _, _, e in entries] + coeffs = [c for c, _, _ in entries] + signs = [_sign_of(c) for c in coeffs] + if c0 == 0 and k == 3 and all(e >= 2 for e in exps) \ + and None not in signs: + names = pe.unknowns # == the diagonal variables here + transform = identity(names) + if abs(sum(signs)) == 3: + # all terms on one side: for an odd exponent we may flip + # the sign of that variable (x -> -x) + flippable = [i for i, e in enumerate(exps) if e % 2 == 1] + if flippable: + i = flippable[0] + coeffs = list(coeffs) + signs = list(signs) + coeffs[i] = -coeffs[i] + signs[i] = -signs[i] + transform = negate( + names, (entries[i][1],), + description=f"substituted {entries[i][1]} -> " + f"-{entries[i][1]} (odd exponent)") + if abs(sum(signs)) == 1: + # orient to a x^p + b y^q - c z^r = 0 with a, b, c > 0 + operations = () + if sum(signs) == -1: + coeffs = [-c for c in coeffs] + signs = [-s for s in signs] + operations = ("multiplied by -1",) + # sorting by exponent reorders the variables, so carry + # each one's name along: a, p and standard x have to + # keep meaning the same term + quads = list(zip(coeffs, [v for _, v, _ in entries], + exps, signs)) + pos = sorted((e, c, v) for c, v, e, s in quads if s > 0) + neg = [(e, -c, v) for c, v, e, s in quads if s < 0] + (p, ca, xsrc), (q, cb, ysrc) = pos + (r, cc, zsrc) = neg[0] + transform = transform.then(rename( + [("x", xsrc), ("y", ysrc), ("z", zsrc)], + operations=operations)) + chi = QQ(1) / p + QQ(1) / q + QQ(1) / r + regime = ("spherical" if chi > 1 else + "euclidean" if chi == 1 else "hyperbolic") + data = {"signature": f"({p}, {q}, {r})", "chi": str(chi), + "regime": regime, "a": str(ca), "b": str(cb), + "c": str(cc)} + out.append(Match( + "generalized-fermat", data=data, transform=transform, + summary=f"generalized Fermat equation of signature " + f"({p}, {q}, {r}), chi = {chi} ({regime})", + )) + if p == q == r >= 3 and abs(ca) == 1 and abs(cb) == 1 \ + and abs(cc) == 1: + out.append(Match( + "fermat", data={"n": str(p)}, transform=transform, + summary=f"Fermat equation with exponent {p}", + )) + return out + if c0 != 0 and len(set(exps)) == 1 and None not in signs: + e = exps[0] + n = -c0 + if e == 3 and k == 3 and all(abs(c) == 1 for c in coeffs): + out.append(Match( + "sum-of-three-cubes", data={"n": str(n)}, + summary=f"sum of three cubes: n = {n}", + )) + return out + if e >= 3 and all(c == 1 for c in coeffs): + out.append(Match( + "waring", data={"k": ZZ(e), "s": k, "n": str(n)}, + summary=f"Waring-type: sum of {k} {e}-th powers " + f"= {n}", + )) + return out + if e >= 3: + out.append(Match( + "diagonal-form", + data={"k": ZZ(e), "s": k, + "coeffs": [str(c) for c in coeffs], + "n": str(n)}, + summary=f"diagonal equation of degree {e}", + )) + return out + if c0 == 0 and k >= 4 and len(set(exps)) == 1 \ + and None not in signs and exps[0] >= 3: + e = exps[0] + s_count = signs.count(1) + t_count = signs.count(-1) + # "equal sums of like powers" means the powers themselves are + # summed: weighted forms such as 2*x^3 + y^3 = z^3 + 7*w^3 are + # diagonal, not equal sums + unit_coefficients = all(abs(c) == 1 for c in coeffs) + if 0 < s_count and 0 < t_count and unit_coefficients: + out.append(Match( + "equal-sums-like-powers", + data={"k": ZZ(e), "s": min(s_count, t_count), + "t": max(s_count, t_count)}, + summary=f"equal sums of {e}-th powers " + f"({min(s_count, t_count)} vs " + f"{max(s_count, t_count)} terms)", + )) + return out + out.append(Match( + "diagonal-form", + data={"k": ZZ(e), "s": k, + "coeffs": [str(c) for c in coeffs], "n": "0"}, + summary=f"diagonal form of degree {e} = 0" + + (" (only trivial real solutions)" + if e % 2 == 0 else ""), + )) + return out + + if k == 3 and P.is_homogeneous(): + if d == 3 and pe.is_concrete and _is_irreducible(P) is True: + try: + smooth = Curve(P).is_smooth() + except Exception: + smooth = None + if smooth: + out.append(Match( + "plane-cubic", data={"F": str(P)}, + summary="smooth plane cubic: genus 1; find a point, " + "then reduce to Weierstrass form (Nagell)", + )) + return out + routed = _genus_route(pe, P, affine=False) + if routed: + out.extend(routed) + return out + + return out + + +# -------------------------------------------------------------------------- +# exponential equations +# -------------------------------------------------------------------------- + +def _match_exponential(pe): + r""" + Matches for equations with exponential content. + + Dispatches on the composition of the terms: purely exponential (Pillai, + S-unit), variable powers (Catalan, symbolic Fermat, Lebesgue–Nagell, + Schinzel–Tijdeman power values), polynomial + one exponential + (Ramanujan–Nagell, Thue–Mahler), with ``polynomial-exponential`` as the + root fallback. + + EXAMPLES:: + + sage: from diophantine_classifier.parsing import parse + sage: from diophantine_classifier.matchers import _match_exponential + sage: [m.slug for m in _match_exponential(parse("x^2 + 7 = 2^n"))] + ['polynomial-exponential', 'ramanujan-nagell'] + sage: [m.slug for m in _match_exponential(parse("x^p - y^q = 1"))] + ['polynomial-exponential', 'catalan'] + sage: [m.slug for m in _match_exponential(parse("3^m - 2^n = 5"))] + ['polynomial-exponential', 'exponential-diophantine', 'pillai'] + """ + terms = pe.terms + pure_exp = [t for t in terms + if t.exp_terms and not t.var_powers and not t.powers] + var_pow = [t for t in terms if t.var_powers] + poly_terms = [t for t in terms if t.is_polynomial] + const_terms = [t for t in poly_terms if t.is_constant] + nonconst_poly = [t for t in poly_terms if not t.is_constant] + mixed = [t for t in terms + if (t.exp_terms and t.powers) or (t.var_powers and t.powers) + or (t.exp_terms and t.var_powers) or len(t.var_powers) > 1] + + out = [Match( + "polynomial-exponential", + data={"exp_terms": len(pure_exp) + len(var_pow) + len(mixed)}, + summary="mixed polynomial-exponential equation", + )] + + c0 = sum((t.coeff for t in const_terms if not t.param_powers), QQ(0)) + param_const = [t for t in const_terms if t.param_powers] + + # ---- purely exponential: sum of c * prod(b^n) terms and constants + if not var_pow and not nonconst_poly and not mixed and pure_exp: + bases = sorted({b for t in pure_exp for b, _ in t.exp_terms}, + key=str) + out.append(Match( + "exponential-diophantine", + data={"terms": len(pure_exp), "bases": [str(b) for b in bases]}, + summary=f"purely exponential equation in bases {bases}", + )) + single = all(len(t.exp_terms) == 1 for t in pure_exp) + if len(pure_exp) == 2 and single and not param_const and c0 != 0: + t1, t2 = pure_exp + s1, s2 = _sign_of(t1.coeff), _sign_of(t2.coeff) + # a^x - b^y = c is a statement about *unit* multipliers: A*2^n is + # not 2^n, and t.coeff alone does not see the A + if _unit_coefficient(t1) and _unit_coefficient(t2) and s1 == -s2 \ + and t1.exp_terms[0][1] != t2.exp_terms[0][1]: + if s1 < 0: + t1, t2 = t2, t1 + (a, xvar), (b, yvar) = t1.exp_terms[0], t2.exp_terms[0] + out.append(Match( + "pillai", + data={"a": str(a), "b": str(b), "c": str(-c0)}, + transform=rename([("x", xvar), ("y", yvar)]), + summary=f"Pillai equation {a}^{xvar} - {b}^{yvar} " + f"= {-c0}", + )) + else: + out.append(Match( + "s-unit", data={"bases": [str(b) for b in bases]}, + summary="three-term S-unit-type equation", + )) + elif len(pure_exp) + (1 if (c0 != 0 or param_const) else 0) == 3: + out.append(Match( + "s-unit", data={"bases": [str(b) for b in bases]}, + summary="three-term S-unit-type equation " + f"(S generated by {bases})", + )) + return out + + # ---- Thue-Mahler: homogeneous binary form = m * prod(p^N) + if not var_pow and not mixed and len(pure_exp) == 1 and nonconst_poly \ + and not const_terms: + polyvars = sorted({v for t in nonconst_poly for v, _ in t.powers}) + degs = {sum(e for _, e in t.powers) for t in nonconst_poly} + if len(polyvars) == 2 and len(degs) == 1: + d = degs.pop() + # m is a fixed constant multiplier: a parametric one is not it + if d >= 3 and not pure_exp[0].param_powers: + t = pure_exp[0] + m = -t.coeff + primes = [b for b, _ in t.exp_terms] + form = " + ".join( + "%s*%s" % (tt.coefficient_string(), + "*".join(f"{v}^{e}" for v, e in tt.powers)) + for tt in nonconst_poly) + out.append(Match( + "thue-mahler", + data={"form": form, "m": str(m), "degree": ZZ(d), + "primes": [str(p) for p in primes]}, + summary=f"Thue-Mahler equation of degree {d} with " + f"primes {primes}", + )) + return out + + # ---- var^var families + simple_vv = [t for t in var_pow + if len(t.var_powers) == 1 and not t.powers + and not t.exp_terms] + if len(simple_vv) == len(var_pow) and var_pow and not pure_exp \ + and not mixed: + if len(var_pow) == 2 and not nonconst_poly and not param_const: + t1, t2 = var_pow + s1, s2 = _sign_of(t1.coeff), _sign_of(t2.coeff) + if _unit_coefficient(t1) and _unit_coefficient(t2) and s1 == -s2 \ + and c0 != 0: + if s1 < 0: + t1, t2 = t2, t1 + (xb, pv), (yb, qv) = t1.var_powers[0], t2.var_powers[0] + c = -c0 + if xb != yb and pv != qv and abs(c) == 1: + if c == -1: + (xb, pv), (yb, qv) = (yb, qv), (xb, pv) + out.append(Match( + "catalan", + transform=rename([("x", xb), ("p", pv), ("y", yb), + ("q", qv)], + description="oriented as x^p - y^q = 1" + if c == -1 else ""), + summary=f"Catalan equation {xb}^{pv} - {yb}^{qv} " + "= 1: only 3^2 - 2^3 = 1 (Mihailescu)", + )) + return out + if xb != yb: + out.append(Match( + "pillai", + data={"c": str(c), "variable_bases": True}, + transform=rename([("x", xb), ("p", pv), ("y", yb), + ("q", qv)]), + summary=f"perfect-power difference {xb}^{pv} - " + f"{yb}^{qv} = {c} (Pillai's conjecture)", + )) + return out + if len(var_pow) == 3 and not nonconst_poly and c0 == 0 \ + and not param_const: + coeffs = [t.coeff for t in var_pow] + signs = [_sign_of(c) for c in coeffs] + expvars = [t.var_powers[0][1] for t in var_pow] + bases = [t.var_powers[0][0] for t in var_pow] + if all(_unit_coefficient(t) for t in var_pow) and None not in signs \ + and abs(sum(signs)) == 1 and len(set(bases)) == 3: + # the exponents are themselves unknowns, and x^p + y^p = z^q + # reuses one of them. A rename onto independent canonical + # p, q, r would claim a bijection that does not exist, so + # the coordinate map stays the identity and the structure is + # recorded as roles instead + roles = {"bases": list(bases), "exponents": list(expvars)} + if len(set(expvars)) == 1: + out.append(Match( + "fermat", data={"n": expvars[0], "symbolic": True}, + transform=identity(pe.unknowns, roles=roles), + summary=f"Fermat equation with unknown exponent " + f"{expvars[0]} (no solutions for " + f"{expvars[0]} >= 3, Wiles)", + )) + return out + # any other symbolic signature is generalized Fermat: the + # exponents need not be pairwise distinct, so x^p + y^p = z^q + # belongs here just as x^p + y^q = z^r does + out.append(Match( + "generalized-fermat", + data={"signature": f"({', '.join(expvars)})", + "symbolic": True}, + transform=identity(pe.unknowns, roles=roles), + summary="generalized Fermat with unknown exponents " + "(Beal/Fermat-Catalan territory)", + )) + return out + if len(var_pow) == 1 and nonconst_poly and not param_const: + t = var_pow[0] + base, expvar = t.var_powers[0] + polyvars = sorted({v for tt in nonconst_poly + for v, _ in tt.powers}) + # f(u) is read off as plain rational coefficients below (its + # leading one has to be 1 for Lebesgue-Nagell), so a parametric + # coefficient anywhere would be dropped silently + if len(polyvars) == 1 and polyvars[0] != base \ + and _unit_coefficient(t) \ + and _concrete_coefficients(poly_terms): + u = polyvars[0] + # f(u) = c * base^expvar with c = -t.coeff + fcoeffs = {} + for tt in poly_terms: + e = tt.powers[0][1] if tt.powers else 0 + fcoeffs[e] = fcoeffs.get(e, QQ(0)) + tt.coeff + if _sign_of(t.coeff) == 1: + fcoeffs = {e: -c for e, c in fcoeffs.items()} + degf = max(fcoeffs) + fstr = " + ".join(f"({c})*{u}^{e}" if e else f"({c})" + for e, c in sorted(fcoeffs.items(), + reverse=True)) + power_map = rename([("x", u), ("y", base), ("n", expvar)]) + if degf == 2 and fcoeffs.get(2) == 1 and not fcoeffs.get(1): + dval = fcoeffs.get(0, QQ(0)) + out.append(Match( + "lebesgue-nagell", data={"d": str(dval)}, + transform=power_map, + summary=f"Lebesgue-Nagell equation {u}^2 + ({dval}) " + f"= {base}^{expvar}", + )) + return out + if degf >= 2: + out.append(Match( + "power-values", data={"f": fstr}, + transform=power_map, + summary=f"power values of a polynomial: {fstr} = " + f"{base}^{expvar} (Schinzel-Tijdeman)", + )) + return out + + # ---- polynomial part + one true exponential: Ramanujan-Nagell type + if not var_pow and not mixed and len(pure_exp) == 1 and nonconst_poly: + polyvars = sorted({v for t in nonconst_poly for v, _ in t.powers}) + # the multiplier k and the quadratic's coefficients (monicity below) + # are read as plain rationals: a parametric factor would vanish, and + # A*x^2 + 7 = 2^n is not a Ramanujan-Nagell equation + if len(polyvars) == 1 and not pure_exp[0].param_powers \ + and _concrete_coefficients(poly_terms): + u = polyvars[0] + t = pure_exp[0] + k_coeff = -t.coeff + fcoeffs = {} + for tt in poly_terms: + e = tt.powers[0][1] if tt.powers else 0 + fcoeffs[e] = fcoeffs.get(e, QQ(0)) + tt.coeff + degf = max(fcoeffs) + if degf == 2: + if _sign_of(fcoeffs[2]) == -1: + fcoeffs = {e: -c for e, c in fcoeffs.items()} + k_coeff = -k_coeff + base_desc = "*".join(f"{b}^{n}" for b, n in t.exp_terms) + dval = fcoeffs.get(0, QQ(0)) + monic_pure = (fcoeffs.get(2) == 1 and not fcoeffs.get(1)) + data = {"k": str(k_coeff), "d": str(dval), + "base": str(t.exp_terms[0][0]), + "exp": base_desc} + classical = (monic_pure and k_coeff == 1 and dval == 7 + and len(t.exp_terms) == 1 + and t.exp_terms[0][0] == 2) + out.append(Match( + "ramanujan-nagell", data=data, + transform=rename([("x", u), + ("n", t.exp_terms[0][1])]), + summary=("the Ramanujan-Nagell equation x^2 + 7 = 2^n: " + "n in {3, 4, 5, 7, 15}" if classical else + f"generalized Ramanujan-Nagell: quadratic in " + f"{u} = ({k_coeff})*{base_desc}"), + )) + return out + + return out + + +# -------------------------------------------------------------------------- +# unit fractions +# -------------------------------------------------------------------------- + +def _match_fractional(pe): + r""" + Matches from the unit-fraction structure of the *uncleared* input. + + EXAMPLES:: + + sage: from diophantine_classifier.parsing import parse + sage: from diophantine_classifier.matchers import _match_fractional + sage: pe = parse("4/n = 1/x + 1/y + 1/z", params="n") + sage: [m.slug for m in _match_fractional(pe)] + ['egyptian-fractions', 'erdos-straus'] + sage: _match_fractional(parse("x^2 + y^2 = 610")) + [] + """ + fs = pe.fractional + out = [] + if fs is None: + return out + k, a, n = fs["k"], fs["a"], fs["n"] + out.append(Match( + "egyptian-fractions", + data={"k": ZZ(k), "a": str(a), "n": str(n), + "unit_vars": list(fs["unit_vars"])}, + summary=f"sum of {k} unit fractions = {a}/{n}", + )) + if k == 3 and a == 4: + out.append(Match( + "erdos-straus", data={"n": str(n), "a": "4", "k": "3"}, + summary=f"Erdos-Straus equation 4/{n} = 1/x + 1/y + 1/z", + )) + elif k == 3 and a == 5: + out[-1].summary += " (Sierpinski's 5/n problem)" + return out + + +def run(pe): + r""" + Collect matches from all applicable matchers, deduplicated by slug. + + INPUT: + + - ``pe`` -- a :class:`~diophantine_classifier.parsing.ParsedEquation` + + OUTPUT: list of :class:`Match`, in emission order (ranking by + specificity happens in :mod:`~diophantine_classifier.classify`) + + Every emitted match carries a coordinate map; a matcher that found the + equation already in standard form gets the identity on the equation's + own unknowns, so a solver can always call ``pull_back``:: + + sage: from diophantine_classifier.parsing import parse + sage: from diophantine_classifier.matchers import run + sage: m = run(parse("3*x + 5*y = 1"))[-1] + sage: m.transform.is_identity, m.transform.source_variables + (True, ('x', 'y')) + + EXAMPLES:: + + sage: from diophantine_classifier.parsing import parse + sage: from diophantine_classifier.matchers import run + sage: [m.slug for m in run(parse("y^2 = x^3 - 2"))] + ['general-polynomial', 'elliptic-weierstrass', 'mordell'] + sage: [m.slug for m in run(parse("2^a + 3^b = 5^c"))] + ['polynomial-exponential', 'exponential-diophantine', 's-unit'] + """ + matches = [] + matches.extend(_match_fractional(pe)) + if pe.is_polynomial: + matches.extend(_match_polynomial(pe)) + else: + matches.extend(_match_exponential(pe)) + seen = set() + result = [] + for m in matches: + if m.slug in seen: + continue + seen.add(m.slug) + if not m.transform.source_variables: + # a matcher that normalized nothing still owes the solvers a + # map: the identity on this equation's unknowns + m.transform = identity(pe.unknowns, + description=m.transform.description, + operations=m.transform.operations, + roles=m.transform.roles) + result.append(m) + return result diff --git a/diophantine_classifier/transforms.py b/diophantine_classifier/transforms.py new file mode 100644 index 0000000..d866e4f --- /dev/null +++ b/diophantine_classifier/transforms.py @@ -0,0 +1,566 @@ +r""" +Coordinate maps between a user's equation and the form a matcher recognized. + +A matcher rarely finds an equation in exactly the standard form of its +family. ``5*x^2 - y^2 = 1`` is a Pell equation, but only after reading the +user's ``y`` as the Pell ``x``; ``x^3 + y^3 = z^3`` becomes a generalized +Fermat equation of signature `(3, 3, 3)` only after substituting +``z -> -z``. Recording those normalizations as prose is enough to explain +what happened and useless for anything else: a solver that works in the +standard coordinates has to carry its answer *back*, and a string cannot +transport a solution. + +:class:`CoordinateTransform` is that map, in both directions and executable. +It also carries the structural :attr:`~CoordinateTransform.roles` a family +assigns to variables, so solvers can build an assignment in canonical names +and pull it back without reaching into a match's ``data`` dict. + +Operations that change the *equation* but not the coordinates — multiplying +through by `-1`, say — are not coordinate changes and live in +:attr:`~CoordinateTransform.operations` instead. + +EXAMPLES:: + + sage: from diophantine_classifier.transforms import rename, identity + sage: swap = rename([("x", "b"), ("y", "a")], description="swapped variables") + sage: swap.push_forward({"a": 2, "b": 3}) # into standard coordinates + {'x': 3, 'y': 2} + sage: swap.pull_back({"x": 3, "y": 2}) # and back again + {'a': 2, 'b': 3} + sage: identity(("x", "y")).is_identity + True +""" + +from dataclasses import dataclass, field + +from sage.all import QQ, SR, ZZ +from sage.calculus.calculus import symbolic_expression_from_string + + +def _identifiers(text): + r""" + The identifiers occurring in an expression string. + + EXAMPLES:: + + sage: from diophantine_classifier.transforms import _identifiers + sage: sorted(_identifiers("-x + y_1")) + ['x', 'y_1'] + """ + import re + return set(re.findall(r"[A-Za-z_][A-Za-z_0-9]*", text)) + + +def _substitute(text, mapping): + r""" + Rewrite an expression by replacing its identifiers with expressions. + + INPUT: + + - ``text`` -- string; an expression + - ``mapping`` -- dict mapping identifier names to expression strings + + OUTPUT: string; the rewritten expression + + EXAMPLES:: + + sage: from diophantine_classifier.transforms import _substitute + sage: _substitute("x", {"x": "-v"}) + '-v' + sage: _substitute("-z", {"z": "-w"}) + 'w' + sage: _substitute("x + 1", {"y": "u"}) # nothing to do + 'x + 1' + """ + names = _identifiers(text) + replace = {name: str(mapping[name]) for name in names if name in mapping} + if not replace: + return text + table = {name: SR.var(name) for name in names} + expr = symbolic_expression_from_string(text, table) + subs = {} + for name, repl in replace.items(): + rnames = _identifiers(repl) + subs[SR.var(name)] = symbolic_expression_from_string( + repl, {other: SR.var(other) for other in rnames}) + return str(expr.subs(subs)) + + +def _evaluate(text, assignment): + r""" + Evaluate one coordinate expression on an assignment, exactly. + + Identifiers resolve against ``assignment`` only, so Sage's global names + (``e``, ``I``, ``pi``, ...) cannot leak into a coordinate map. Integer + and rational results come back as Sage numbers rather than as symbolic + expressions, so transported solutions compare and print like the ones a + solver produced directly. + + INPUT: + + - ``text`` -- string; an expression in the other coordinate system + - ``assignment`` -- dict mapping variable names to values + + OUTPUT: the value of the expression + + EXAMPLES:: + + sage: from diophantine_classifier.transforms import _evaluate + sage: _evaluate("-x", {"x": 5}) + -5 + sage: _evaluate("x", {"x": QQ(1)/2}) + 1/2 + sage: _evaluate("x + y", {"x": 1}) + Traceback (most recent call last): + ... + KeyError: "coordinate expression 'x + y' needs a value for 'y'" + """ + names = _identifiers(text) + missing = sorted(name for name in names if name not in assignment) + if missing: + raise KeyError(f"coordinate expression {text!r} needs a value for " + f"{missing[0]!r}") + table = {name: SR.var(name) for name in names} + value = symbolic_expression_from_string(text, table) + if names: + value = value.subs({SR.var(name): assignment[name] + for name in names}) + for ring in (ZZ, QQ): + try: + return ring(value) + except (TypeError, ValueError): + continue + return value + + +def _rewrite_role(value, mapping): + r""" + Rewrite a role's variable (or list of variables) through a mapping. + + INPUT: + + - ``value`` -- a variable name, or a list/tuple of them + - ``mapping`` -- dict from those names to expressions + + OUTPUT: the rewritten name, or list of names + + EXAMPLES:: + + sage: from diophantine_classifier.transforms import _rewrite_role + sage: _rewrite_role("v", {"v": "b"}) + 'b' + sage: _rewrite_role(["u", "v"], {"u": "a", "v": "b"}) + ['a', 'b'] + """ + if isinstance(value, (list, tuple)): + return [_substitute(str(item), mapping) for item in value] + return _substitute(str(value), mapping) + + +@dataclass(frozen=True) +class CoordinateTransform: + r""" + An invertible change of variables between source and standard form. + + Both directions are stored explicitly, each keyed by *its own* system's + variables and written in terms of the other's: ``to_normalized["x"]`` is + the formula for the standard coordinate ``x`` in the user's variables, + and ``to_source["a"]`` is the formula for the user's ``a`` in standard + coordinates. Expression maps cover the swaps and sign changes the + matchers perform today and leave room for affine or birational maps + later. + + ATTRIBUTES: + + - ``description`` -- string; what the normalization did, for display. + - ``source_variables`` -- tuple of strings; the unknowns as the user + wrote them. + - ``normalized_variables`` -- tuple of strings; the family's standard + coordinate names. + - ``to_normalized`` -- dict mapping each standard coordinate to an + expression in the source variables. + - ``to_source`` -- dict mapping each source variable to an expression in + the standard coordinates. + - ``conditions`` -- tuple of + :class:`~diophantine_classifier.conditions.NonzeroCondition`; where + the map is defined, when that is not everywhere. + - ``roles`` -- dict mapping a structural role name to the source + variable (or list of source variables) playing it, e.g. + ``{"legs": ["x", "y"], "hypotenuse": "z"}``. Solver-relevant, hence + here rather than buried in a match's ``data``. + - ``operations`` -- tuple of strings; normalizations that changed the + equation but no coordinate, such as multiplying through by `-1`. + + EXAMPLES:: + + sage: from diophantine_classifier.transforms import CoordinateTransform + sage: flip = CoordinateTransform( + ....: description="substituted z -> -z (odd exponent)", + ....: source_variables=("z",), normalized_variables=("z",), + ....: to_normalized={"z": "-z"}, to_source={"z": "-z"}) + sage: flip.push_forward({"z": 4}) + {'z': -4} + sage: flip.pull_back(flip.push_forward({"z": 4})) + {'z': 4} + """ + description: str = "" + source_variables: tuple = () + normalized_variables: tuple = () + to_normalized: dict = field(default_factory=dict) + to_source: dict = field(default_factory=dict) + conditions: tuple = () + roles: dict = field(default_factory=dict) + operations: tuple = () + + @property + def is_identity(self): + r""" + Whether the map leaves every coordinate alone. + + Equation-level :attr:`operations` do not count: they change no + coordinate, so a sign-flipped equation still has an identity map. + + EXAMPLES:: + + sage: from diophantine_classifier.transforms import identity, rename + sage: identity(("x", "y"), operations=("multiplied by -1",)).is_identity + True + sage: rename([("x", "y"), ("y", "x")]).is_identity + False + """ + return all(name == expr for name, expr in self.to_normalized.items()) + + def push_forward(self, assignment): + r""" + Carry an assignment in the user's variables to standard coordinates. + + INPUT: + + - ``assignment`` -- dict mapping source variable names to values + + OUTPUT: dict mapping standard coordinate names to values + + EXAMPLES:: + + sage: from diophantine_classifier.transforms import rename + sage: t = rename([("x", "b"), ("y", "a")]) + sage: t.push_forward({"a": 2, "b": 3}) + {'x': 3, 'y': 2} + """ + return {name: _evaluate(expr, assignment) + for name, expr in self.to_normalized.items()} + + def pull_back(self, assignment): + r""" + Carry an assignment in standard coordinates back to the user's. + + This is what turns a solver's answer into a solution of the problem + that was submitted. + + INPUT: + + - ``assignment`` -- dict mapping standard coordinate names to values + + OUTPUT: dict mapping source variable names to values + + EXAMPLES:: + + sage: from diophantine_classifier.transforms import rename + sage: t = rename([("x", "b"), ("y", "a")]) + sage: t.pull_back({"x": 3, "y": 2}) + {'a': 2, 'b': 3} + """ + return {name: _evaluate(expr, assignment) + for name, expr in self.to_source.items()} + + def conditions_hold(self, assignment): + r""" + Decide this map's own side conditions on a *source* assignment. + + :attr:`conditions` are stored in the transform's source coordinates + (see :meth:`then`), so they are checked against an assignment in the + user's variables — the same one the parser's conditions are checked + against. + + INPUT: + + - ``assignment`` -- dict mapping source variable names to values + + OUTPUT: ``True``, ``False``, or ``None`` when undetermined + + EXAMPLES:: + + sage: from diophantine_classifier.conditions import NonzeroCondition + sage: from diophantine_classifier.transforms import rename + sage: t = rename([("x", "u"), ("y", "v")], + ....: conditions=(NonzeroCondition("u", ("u",)),)) + sage: t.conditions_hold({"u": 1, "v": 0}) + True + sage: t.conditions_hold({"u": 0, "v": 1}) + False + """ + from .conditions import hold + return hold(self.conditions, assignment) + + def then(self, after): + r""" + Compose two maps: ``self`` from `A` to `B`, ``after`` from `B` to + `C`, giving the map from `A` to `C`. + + The composite satisfies, for every assignment ``a`` in `A` and ``c`` + in `C`:: + + self.then(after).push_forward(a) + == after.push_forward(self.push_forward(a)) + self.then(after).pull_back(c) + == self.pull_back(after.pull_back(c)) + + Conditions live in each map's **source** coordinates, so ``after``'s + are pulled back through ``self`` and the composite's conditions are + all checkable on an assignment in `A`. Equation-level + :attr:`operations` are concatenated in the order they happened, and + roles naming `B`-variables are rewritten to their `A`-expressions. + + INPUT: + + - ``after`` -- a :class:`CoordinateTransform` whose source + coordinates are this one's normalized coordinates + + OUTPUT: a :class:`CoordinateTransform` + + EXAMPLES:: + + sage: from diophantine_classifier.transforms import negate, rename + sage: first = negate(("u", "v"), ("v",)) # A -> B + sage: second = rename([("x", "v"), ("y", "u")]) # B -> C + sage: both = first.then(second) + sage: both.to_normalized == {"x": "-v", "y": "u"} + True + sage: both.push_forward({"u": 2, "v": 3}) + {'x': -3, 'y': 2} + sage: both.pull_back(both.push_forward({"u": 2, "v": 3})) + {'u': 2, 'v': 3} + + Composing with an identity changes nothing:: + + sage: from diophantine_classifier.transforms import identity + sage: identity(("u", "v")).then(second).to_normalized + {'x': 'v', 'y': 'u'} + """ + parts = [text for text in (self.description, after.description) + if text] + return CoordinateTransform( + description="; ".join(parts), + source_variables=self.source_variables, + normalized_variables=after.normalized_variables, + # a C-coordinate is written in B; rewrite each B-name in A + to_normalized={name: _substitute(expr, self.to_normalized) + for name, expr in after.to_normalized.items()}, + # an A-variable is written in B; rewrite each B-name in C + to_source={name: _substitute(expr, after.to_source) + for name, expr in self.to_source.items()}, + conditions=(tuple(self.conditions) + + tuple(c.substitute(self.to_normalized) + for c in after.conditions)), + roles={**{key: _rewrite_role(value, self.to_normalized) + for key, value in after.roles.items()}, + **dict(self.roles)}, + operations=tuple(self.operations) + tuple(after.operations), + ) + + def as_dict(self): + r""" + JSON-serializable form (the website-backend contract). + + OUTPUT: dict with plain types only + + EXAMPLES:: + + sage: from diophantine_classifier.transforms import rename + sage: import json + sage: d = rename([("x", "y"), ("y", "x")], + ....: description="swapped variables").as_dict() + sage: d["to_normalized"], d["to_source"] + ({'x': 'y', 'y': 'x'}, {'x': 'y', 'y': 'x'}) + sage: _ = json.dumps(d) + """ + out = { + "description": self.description, + "identity": self.is_identity, + "source_variables": list(self.source_variables), + "normalized_variables": list(self.normalized_variables), + "to_normalized": {str(k): str(v) + for k, v in self.to_normalized.items()}, + "to_source": {str(k): str(v) for k, v in self.to_source.items()}, + } + if self.conditions: + out["conditions"] = [c.as_dict() for c in self.conditions] + if self.roles: + out["roles"] = {str(k): (list(v) if isinstance(v, (list, tuple)) + else str(v)) + for k, v in self.roles.items()} + if self.operations: + out["operations"] = list(self.operations) + return out + + def __str__(self): + r""" + One-line rendering: the coordinate change and the equation-level + operations, or the empty string when neither happened. + + EXAMPLES:: + + sage: from diophantine_classifier.transforms import identity, rename + sage: str(rename([("x", "y"), ("y", "x")], + ....: description="swapped variables", + ....: operations=("multiplied by -1",))) + 'multiplied by -1; swapped variables' + sage: str(identity(("x",))) + '' + """ + parts = list(self.operations) + if self.description: + parts.append(self.description) + return "; ".join(parts) + + def __bool__(self): + r""" + Whether anything happened at all (used for "was there a transform?"). + + EXAMPLES:: + + sage: from diophantine_classifier.transforms import identity, rename + sage: bool(identity(("x", "y"))) + False + sage: bool(rename([("x", "y"), ("y", "x")])) + True + """ + return bool(self.operations) or not self.is_identity + + +def identity(variables, description="", operations=(), roles=None, + normalized_variables=None, conditions=()): + r""" + The transform that renames nothing. + + The default for a match found already in standard form — including one + reached by an equation-level operation, which changes no coordinate. + + INPUT: + + - ``variables`` -- iterable of source variable names + - ``description`` -- (default: ``""``) display text + - ``operations`` -- (default: ``()``) equation-level normalizations + - ``roles`` -- (default: ``None``) structural role map + - ``normalized_variables`` -- (default: ``None``) standard coordinate + names, when they should be recorded as something other than + ``variables`` + - ``conditions`` -- (default: ``()``) where the map is defined + + OUTPUT: a :class:`CoordinateTransform` + + EXAMPLES:: + + sage: from diophantine_classifier.transforms import identity + sage: t = identity(("x", "y"), operations=("multiplied by -1",)) + sage: t.push_forward({"x": 1, "y": 2}) + {'x': 1, 'y': 2} + sage: str(t) + 'multiplied by -1' + """ + variables = tuple(str(v) for v in variables) + same = {v: v for v in variables} + return CoordinateTransform( + description=description, + source_variables=variables, + normalized_variables=tuple(normalized_variables) + if normalized_variables is not None else variables, + to_normalized=dict(same), to_source=dict(same), + conditions=tuple(conditions), + roles=dict(roles or {}), operations=tuple(operations), + ) + + +def rename(pairs, description="", operations=(), roles=None, conditions=()): + r""" + A transform that only renames or permutes coordinates. + + Covers every variable-role reversal the matchers perform: reading the + user's second variable as the standard ``x``, calling the base of a + curve ``x`` and the squared variable ``y``, and so on. + + INPUT: + + - ``pairs`` -- iterable of ``(normalized name, source name)`` + - ``description`` -- (default: ``""``) display text + - ``operations`` -- (default: ``()``) equation-level normalizations + - ``roles`` -- (default: ``None``) structural role map; defaults to the + pairing itself + - ``conditions`` -- (default: ``()``) where the map is defined + + OUTPUT: a :class:`CoordinateTransform` + + EXAMPLES:: + + sage: from diophantine_classifier.transforms import rename + sage: t = rename([("x", "u"), ("y", "v")]) + sage: t.push_forward({"u": 7, "v": 8}) + {'x': 7, 'y': 8} + sage: t.pull_back({"x": 7, "y": 8}) + {'u': 7, 'v': 8} + sage: t.roles + {'x': 'u', 'y': 'v'} + """ + pairs = [(str(n), str(s)) for n, s in pairs] + return CoordinateTransform( + description=description, + source_variables=tuple(s for _, s in pairs), + normalized_variables=tuple(n for n, _ in pairs), + to_normalized={n: s for n, s in pairs}, + to_source={s: n for n, s in pairs}, + conditions=tuple(conditions), + roles=dict(roles) if roles is not None else {n: s for n, s in pairs}, + operations=tuple(operations), + ) + + +def negate(variables, flipped, description="", operations=(), roles=None, + conditions=()): + r""" + A transform substituting ``v -> -v`` for the named variables. + + Used to orient diagonal equations: at an odd exponent, flipping the sign + of a variable moves its term to the other side of the equation, which is + how ``x^3 + y^3 + z^3 = 0`` is read as a generalized Fermat equation. + The map is its own inverse. + + INPUT: + + - ``variables`` -- iterable of all source variable names + - ``flipped`` -- iterable of the names to negate + - ``description`` -- (default: ``""``) display text + - ``operations`` -- (default: ``()``) equation-level normalizations + - ``roles`` -- (default: ``None``) structural role map + - ``conditions`` -- (default: ``()``) where the map is defined + + OUTPUT: a :class:`CoordinateTransform` + + EXAMPLES:: + + sage: from diophantine_classifier.transforms import negate + sage: t = negate(("x", "y", "z"), ("z",)) + sage: t.push_forward({"x": 1, "y": 2, "z": 3}) + {'x': 1, 'y': 2, 'z': -3} + sage: t.pull_back(t.push_forward({"x": 1, "y": 2, "z": 3})) + {'x': 1, 'y': 2, 'z': 3} + """ + variables = tuple(str(v) for v in variables) + flipped = {str(v) for v in flipped} + mapping = {v: (f"-{v}" if v in flipped else v) for v in variables} + return CoordinateTransform( + description=description, + source_variables=variables, normalized_variables=variables, + to_normalized=dict(mapping), to_source=dict(mapping), + conditions=tuple(conditions), + roles=dict(roles or {}), operations=tuple(operations), + ) diff --git a/tests/test_matchers.py b/tests/test_matchers.py new file mode 100644 index 0000000..3e19072 --- /dev/null +++ b/tests/test_matchers.py @@ -0,0 +1,116 @@ +"""Matchers: what a named family requires, and what it must not swallow.""" + +from diophantine_classifier.matchers import run +from diophantine_classifier.parsing import parse + + +def slugs(equation, **kwargs): + return {m.slug for m in run(parse(equation, **kwargs))} + + +def test_parametric_quadratic_form_does_not_crash(): + """Coefficients live in QQ[params]; a Gram matrix must not force QQ.""" + found = slugs("x^2 + y^2 = D*z^2", params="D") + assert "quadratic-form-zero" in found + + +def test_weighted_diagonal_is_not_equal_sums(): + """Equal sums of like powers means the powers are summed, unweighted.""" + found = slugs("2*x^3 + y^3 = z^3 + 7*w^3") + assert "equal-sums-like-powers" not in found + assert "diagonal-form" in found + + +def test_unweighted_equal_sums_still_matches(): + assert "equal-sums-like-powers" in slugs("x^4 + y^4 + z^4 = w^4") + + +def test_sum_of_three_cubes_still_matches(): + assert "sum-of-three-cubes" in slugs("x^3 + y^3 + z^3 = 42") + + +# --- the Pell orientation swap is an equivalence --------------------------- + +def test_swapped_pell_keeps_the_same_equation(): + """5*x^2 - y^2 = 1 reads as y^2 - 5*x^2 = -1, so N flips with the swap.""" + match, = [m for m in run(parse("5*x^2 - y^2 = 1")) if m.slug == "pell"] + assert match.data == {"D": "5", "N": "-1"} + + +def test_unswapped_pell_is_unchanged(): + match, = [m for m in run(parse("x^2 - 61*y^2 = 1")) if m.slug == "pell"] + assert match.data == {"D": "61", "N": "1"} + + +def test_swapped_pell_like_keeps_the_same_equation(): + match, = [m for m in run(parse("3*x^2 - y^2 = 6")) if m.slug == "pell-like"] + assert match.data == {"D": "3", "N": "-6"} + + +# --- parametric coefficients are part of the coefficient (brief 5.1) ------- + +def test_parametric_multiplier_blocks_pillai(): + """A*2^n - 3^m = 1 is not Pillai's equation: the multiplier is unknown.""" + assert "pillai" not in slugs("A*2^n - 3^m = 1", params="A") + + +def test_concrete_multiplier_still_matches_pillai(): + assert "pillai" in slugs("2^n - 3^m = 1") + + +def test_parametric_leading_coefficient_blocks_ramanujan_nagell(): + assert "ramanujan-nagell" not in slugs("A*x^2 + 7 = 2^n", params="A") + + +def test_concrete_leading_coefficient_still_matches_ramanujan_nagell(): + assert "ramanujan-nagell" in slugs("x^2 + 7 = 2^n") + + +def test_parametric_multiplier_blocks_variable_base_pillai(): + """The variable-base branch reads |coeff| = 1 too.""" + assert "pillai" not in slugs("A*x^p - y^q = 2", params="A") + assert "catalan" not in slugs("A*x^p - y^q = 1", params="A") + + +def test_parametric_coefficient_blocks_symbolic_fermat(): + assert "fermat" not in slugs("A*x^n + y^n = z^n", params="A") + + +def test_parametric_coefficient_blocks_lebesgue_nagell(): + assert "lebesgue-nagell" not in slugs("A*x^2 + 3 = y^n", params="A") + assert "lebesgue-nagell" in slugs("x^2 + 3 = y^n") + + +def test_parametric_multiplier_blocks_thue_mahler(): + assert "thue-mahler" not in slugs("x^3 + 2*y^3 = A*2^n", params="A") + assert "thue-mahler" in slugs("x^3 + 2*y^3 = 5*2^n") + + +def test_parametric_coefficient_survives_in_the_broad_fallback(): + """A blocked classical shape still gets its root family.""" + assert "polynomial-exponential" in slugs("A*x^2 + 7 = 2^n", params="A") + + +def test_full_coefficient_is_available_on_terms(): + from diophantine_classifier.parsing import parse + pe = parse("A*x^2 + 7 = 2^n", params="A") + quadratic, = [t for t in pe.terms if t.powers] + assert quadratic.has_parametric_coefficient + assert quadratic.coefficient_string() == "A" + + +# --- symbolic generalized-Fermat signatures (brief 5.4b) ------------------ + +def test_repeated_symbolic_signature_is_generalized_fermat(): + """x^p + y^p = z^q: the exponents need not be pairwise distinct.""" + assert "generalized-fermat" in slugs("x^p + y^p = z^q") + + +def test_distinct_symbolic_signature_is_still_generalized_fermat(): + assert "generalized-fermat" in slugs("x^p + y^q = z^r") + + +def test_all_equal_symbolic_signature_is_fermat(): + found = slugs("x^n + y^n = z^n") + assert "fermat" in found + assert "generalized-fermat" not in found diff --git a/tests/test_registry.py b/tests/test_registry.py index e59e243..c401ba4 100644 --- a/tests/test_registry.py +++ b/tests/test_registry.py @@ -15,6 +15,24 @@ def test_dag_depths(): assert "general-polynomial" in ancestors("linear") +def test_matcher_flags_match_implementation(): + """Every registered slug the matcher code can emit is flagged matcher: true. + + A recognizer whose family has not landed in the registry yet is inert: + :func:`~diophantine_classifier.classify.classify` drops matches naming an + unregistered family. The check tightens to *every* emitted slug once the + registry is complete. + """ + import re + import diophantine_classifier.matchers as m + source = open(m.__file__.replace(".pyc", ".py")).read() + emitted = set(re.findall(r'Match\(\s*\n?\s*"([a-z0-9-]+)"', source)) + emitted |= set(re.findall(r'Match\("([a-z0-9-]+)"', source)) + fams = families() + for slug in sorted(emitted & set(fams)): + assert fams[slug].matcher, f"{slug}: emitted but matcher flag is false" + + # --- ancestry is paths and edges, not a flattened chain (brief 4.1) ------- def test_path_traversal_on_a_diamond_invents_no_edge(): diff --git a/tests/test_transforms.py b/tests/test_transforms.py new file mode 100644 index 0000000..ff212a0 --- /dev/null +++ b/tests/test_transforms.py @@ -0,0 +1,272 @@ +"""Coordinate maps: every normalization a matcher performs must be +executable and invertible, not merely described in prose. + +A solver works in a family's standard coordinates; the transform is what +carries its answer back to the variables the user wrote. +""" + +import json + +import pytest + +from diophantine_classifier.matchers import run +from diophantine_classifier.parsing import parse +from diophantine_classifier.transforms import (CoordinateTransform, identity, + negate, rename) + + +def match_for(equation, slug, **kwargs): + for m in run(parse(equation, **kwargs)): + if m.slug == slug: + return m + raise AssertionError(f"{equation!r} did not match {slug!r}") + + +# --- the map itself ------------------------------------------------------- + +def test_identity_round_trips_and_knows_it_is_one(): + t = identity(("x", "y")) + assert t.is_identity + assert t.push_forward({"x": 1, "y": 2}) == {"x": 1, "y": 2} + assert t.pull_back({"x": 1, "y": 2}) == {"x": 1, "y": 2} + + +def test_equation_level_operations_are_not_coordinate_changes(): + """Multiplying an equation by -1 changes no coordinate.""" + t = identity(("x", "y"), operations=("multiplied by -1",)) + assert t.is_identity + assert bool(t) # something did happen, though + assert str(t) == "multiplied by -1" + assert t.pull_back({"x": 3, "y": 4}) == {"x": 3, "y": 4} + + +def test_rename_is_invertible(): + t = rename([("x", "b"), ("y", "a")]) + assert t.push_forward({"a": 2, "b": 3}) == {"x": 3, "y": 2} + assert t.pull_back(t.push_forward({"a": 2, "b": 3})) == {"a": 2, "b": 3} + + +def test_negate_is_its_own_inverse(): + t = negate(("x", "y", "z"), ("z",)) + source = {"x": 1, "y": 2, "z": 3} + assert t.push_forward(source) == {"x": 1, "y": 2, "z": -3} + assert t.pull_back(t.push_forward(source)) == source + + +def test_missing_coordinate_is_an_error_not_a_guess(): + t = rename([("x", "a"), ("y", "b")]) + with pytest.raises(KeyError): + t.push_forward({"a": 1}) + + +def test_rational_coordinates_survive_transport(): + from sage.all import QQ + t = negate(("x",), ("x",)) + assert t.pull_back({"x": QQ(1) / 2}) == {"x": QQ(-1) / 2} + + +def test_as_dict_is_json_serializable(): + t = rename([("x", "b"), ("y", "a")], description="swapped variables", + operations=("multiplied by -1",)) + blob = json.dumps(t.as_dict()) + assert "swapped variables" in blob + assert json.loads(blob)["to_source"] == {"a": "y", "b": "x"} + + +def test_default_transform_is_the_identity(): + assert CoordinateTransform().is_identity + + +# --- composition (brief 5.1) --------------------------------------------- + +def test_transform_composition_round_trip(): + first = negate(("u", "v"), ("v",)) # source -> intermediate + second = rename([("x", "v"), ("y", "u")]) # intermediate -> standard + combined = first.then(second) + source = {"u": 2, "v": 3} + assert combined.pull_back(combined.push_forward(source)) == source + + +def test_composition_agrees_with_applying_the_steps_in_turn(): + first = negate(("u", "v"), ("v",)) + second = rename([("x", "v"), ("y", "u")]) + combined = first.then(second) + source = {"u": 5, "v": -7} + assert combined.push_forward(source) == second.push_forward( + first.push_forward(source)) + standard = {"x": 11, "y": 13} + assert combined.pull_back(standard) == first.pull_back( + second.pull_back(standard)) + + +def test_composition_tracks_the_variable_lists(): + first = rename([("a", "u"), ("b", "v")]) + second = rename([("x", "a"), ("y", "b")]) + combined = first.then(second) + assert combined.source_variables == ("u", "v") + assert combined.normalized_variables == ("x", "y") + + +def test_composition_with_identity_changes_nothing(): + t = rename([("x", "v"), ("y", "u")]) + assert identity(("u", "v")).then(t).to_normalized == t.to_normalized + assert t.then(identity(("x", "y"))).to_source == t.to_source + + +def test_composition_keeps_operations_in_order(): + first = identity(("u",), operations=("multiplied by -1",)) + second = identity(("u",), operations=("divided by the content",)) + assert first.then(second).operations == ("multiplied by -1", + "divided by the content") + + +def test_composition_transports_conditions(): + """A condition of the second map is stated in *its* source coordinates; + composing must pull it back so the whole map's conditions are checkable + on the original assignment.""" + from diophantine_classifier.conditions import NonzeroCondition + first = negate(("u", "v"), ("v",)) + second = rename([("x", "v"), ("y", "u")], + conditions=(NonzeroCondition("v", ("v",)),)) + combined = first.then(second) + assert [str(c) for c in combined.conditions] == ["-v != 0"] + assert combined.conditions_hold({"u": 1, "v": 3}) is True + assert combined.conditions_hold({"u": 1, "v": 0}) is False + + +def test_composed_conditions_serialize(): + from diophantine_classifier.conditions import NonzeroCondition + first = rename([("a", "u"), ("b", "v")], + conditions=(NonzeroCondition("u", ("u",)),)) + second = rename([("x", "a"), ("y", "b")], + conditions=(NonzeroCondition("b", ("b",)),)) + blob = json.dumps(first.then(second).as_dict()) + assert "conditions" in blob + + +# --- generalized Fermat: data and map agree (brief 5.2) ------------------ + +def test_generalized_fermat_sorting_updates_the_transform(): + m = match_for("u^5 + v^3 = w^7", "generalized-fermat") + assert m.data["signature"] == "(3, 5, 7)" + assert m.transform.to_normalized["x"] == "v" + assert m.transform.to_normalized["y"] == "u" + assert m.transform.to_normalized["z"] == "w" + + +def test_generalized_fermat_coefficients_follow_the_same_variables(): + """a is the coefficient of standard x, which is the user's v.""" + m = match_for("2*u^5 + 3*v^3 = 5*w^7", "generalized-fermat") + assert m.data["signature"] == "(3, 5, 7)" + assert (m.data["a"], m.data["b"], m.data["c"]) == ("3", "2", "5") + assert m.transform.to_normalized["x"] == "v" + + +def test_generalized_fermat_sign_flip_and_permutation_compose(): + m = match_for("u^5 + v^3 + w^7 = 0", "generalized-fermat") + source = {"u": 2, "v": 3, "w": 4} + assert m.transform.pull_back(m.transform.push_forward(source)) == source + negated = [expr for expr in m.transform.to_normalized.values() + if expr.startswith("-")] + assert len(negated) == 1 + + +def test_global_orientation_records_equation_operation(): + m = match_for("-x^3 - y^3 = -z^3", "generalized-fermat") + assert m.transform.operations == ("multiplied by -1",) + + +def test_symbolic_signature_does_not_claim_a_false_bijection(): + """x^p + y^p = z^q reuses an exponent, so there is no invertible rename + onto independent canonical p, q, r; the structure is roles instead.""" + m = match_for("x^p + y^p = z^q", "generalized-fermat") + assert m.transform.is_identity + assert m.transform.roles["exponents"] == ["p", "p", "q"] + assert m.transform.roles["bases"] == ["x", "y", "z"] + + +# --- the maps the matchers actually emit --------------------------------- + +def test_every_match_carries_a_map_over_the_real_unknowns(): + for equation in ["3*x + 5*y = 1", "x^2 - 61*y^2 = 1", "y^2 = x^3 - 2", + "x^2 + y^2 = z^2", "x^2 + 7 = 2^n", "x^p - y^q = 1", + "x^3 + y^3 + z^3 = 0"]: + pe = parse(equation) + for m in run(pe): + assert set(m.transform.source_variables) <= set(pe.unknowns) + assert set(m.transform.to_source) <= set(pe.unknowns) + + +def test_pell_swap_transform_round_trip(): + """5*x^2 - y^2 = 1 is Pell only after reading the user's y as x.""" + m = match_for("5*x^2 - y^2 = 1", "pell") + assert not m.transform.is_identity + assert m.transform.to_normalized == {"x": "y", "y": "x"} + source = {"x": 2, "y": 3} + assert m.transform.pull_back(m.transform.push_forward(source)) == source + + +def test_pell_without_swap_is_the_identity(): + m = match_for("x^2 - 61*y^2 = 1", "pell") + assert m.transform.is_identity + assert m.transform.source_variables == ("x", "y") + + +def test_sign_flip_multiplication_is_an_operation_not_a_swap(): + """-x^2 + 5*y^2 = -1 is oriented by multiplying through by -1.""" + m = match_for("-x^2 + 5*y^2 = -1", "pell") + assert m.transform.operations == ("multiplied by -1",) + assert m.transform.is_identity + + +def test_odd_exponent_sign_flip_round_trip(): + """x^3 + y^3 + z^3 = 0 becomes generalized Fermat by a sign flip, then + a permutation onto the canonical (x, y, z); the two compose into one + invertible map.""" + m = match_for("x^3 + y^3 + z^3 = 0", "generalized-fermat") + assert not m.transform.is_identity + negated = [expr for expr in m.transform.to_normalized.values() + if expr.startswith("-")] + assert len(negated) == 1 + source = {"x": 3, "y": 4, "z": 5} + assert m.transform.pull_back(m.transform.push_forward(source)) == source + + +def test_curve_role_reversal_is_executable(): + """x^3 = y^2 - 2 is Mordell with the user's y as the standard x.""" + m = match_for("x^3 = y^2 - 2", "mordell") + assert m.transform.to_normalized == {"x": "x", "y": "y"} + m2 = match_for("y^2 = x^3 - 2", "mordell") + assert m2.transform.pull_back({"x": 3, "y": 5}) == {"x": 3, "y": 5} + + +def test_curve_roles_follow_the_user_names(): + """u^2 = v^3 - 2: the standard x is the user's v.""" + m = match_for("u^2 = v^3 - 2", "mordell") + assert m.transform.to_normalized == {"x": "v", "y": "u"} + assert m.transform.pull_back({"x": 3, "y": 5}) == {"v": 3, "u": 5} + + +def test_roles_live_on_the_transform_not_in_the_data(): + m = match_for("x^2 + y^2 = z^2", "pythagorean") + assert m.transform.roles["hypotenuse"] == "z" + assert sorted(m.transform.roles["legs"]) == ["x", "y"] + assert "roles" not in m.data + + +def test_catalan_roles_are_a_coordinate_map(): + m = match_for("x^p - y^q = 1", "catalan") + assert m.transform.pull_back({"x": 3, "p": 2, "y": 2, "q": 3}) == { + "x": 3, "p": 2, "y": 2, "q": 3} + + +def test_catalan_orientation_is_invertible(): + """1 = y^q - x^p is the same equation with the roles the other way.""" + m = match_for("y^q - x^p = 1", "catalan") + pulled = m.transform.pull_back({"x": 3, "p": 2, "y": 2, "q": 3}) + assert pulled == {"y": 3, "q": 2, "x": 2, "p": 3} + + +def test_ramanujan_nagell_roles_are_a_coordinate_map(): + m = match_for("u^2 + 7 = 2^k", "ramanujan-nagell") + assert m.transform.pull_back({"x": 181, "n": 15}) == {"u": 181, "k": 15}