Skip to content

Complete seed-orbit enumeration for generalized Pell equations #80

Open
Open
Complete seed-orbit enumeration for generalized Pell equations#80
Labels
enhancementNew feature or request

Description

@roed-math
roed-math
opened on Aug 9, 2026

Goal

For a single generalized Pell equation x^2 - D*y^2 = N (D a concrete positive nonsquare integer, N a concrete nonzero integer), produce the finite set of orbit seeds together with a certified completeness statement, as a reusable primitive: a new module diophantine_classifier/pell_orbits.py with

seed_orbits(D, N) -> SeedOrbits

"Complete" here carries its supporting contract: the full integer solution set of the equation is exactly { ±eps^k · s : k in ZZ, s in seeds } with eps the fundamental automorph, and completeness of the seed list is certified on every call by the agreement of two complete enumeration methods that are independent of each other: PARI's class enumeration, and an exhaustive scan over the strip 0 <= y <= Y whose per-class completeness is proved inline in this issue (Lemmas A and B below; the inline proofs are normative, the book citation is attribution). This is the enumeration primitive behind the whole of #58; it knows nothing about simultaneous systems.

Target base

Dependencies

Supported inputs

  • D: an element of ZZ, D > 0, not a perfect square. Squarefree is not required (parent contract): everything below lives in the order ZZ[sqrt(D)] — the fundamental automorph is the fundamental solution of u^2 - D*v^2 = 1, which exists for every positive nonsquare D.
  • N: an element of ZZ, N != 0, either sign. N = 0 is rejected as degenerate (parent contract: irrationality of sqrt(D) forces (x, y) = (0, 0); there are no seed orbits).
  • Invalid inputs — D <= 0, D a perfect square, N = 0, non-integers — raise ValueError with a message naming the violated condition. This is a library primitive, not a solver: the SolverUnavailable outcome vocabulary belongs to the solver layer, with the single exception of external-software failure noted under "Failure and fallback semantics".

Output and API contract

The SeedOrbits record

A frozen dataclass with fields, all exact (Sage ZZ; never coerced through int(...) or floats):

  • D, N — the inputs;
  • fundamental — the pair (u0, v0) with u0, v0 > 0, the fundamental solution of u^2 - D*v^2 = 1 (u0 >= 2 always, since u0^2 = 1 + D*v0^2 >= 3);
  • seeds — a tuple of pairs (x, y) in ZZ^2, possibly empty, sorted by the key (|y|, y, x);
  • evidence — one typed evidence.Evidence record (below).

Completeness statement (the contract of the record). Let eps act by (x, y) -> (u0*x + D*v0*y, v0*x + u0*y). Then the full integer solution set of x^2 - D*y^2 = N is exactly

{ ±eps^k · s  :  k in ZZ,  s in seeds }

and the seeds are pairwise inequivalent under the automorph group <±eps> — which is exactly the proper automorphism group SO(Q, ZZ) of the form Q = x^2 - D*y^2, i.e. the class notion qfbsolve enumerates (proved in the lemma section below). An empty seeds tuple is the (certified) statement that the equation has no integer solutions. Verified instance (Sage 10.7 / PARI 2.17.2): for D = 13, N = 27 the four class representatives (-40, -11), (-40, 11), (12, -3), (12, 3) are pairwise inequivalent under <eps> and under negation composed with <eps>, the normalized seed tuple is exactly ((-12, 3), (12, 3), (-40, 11), (40, 11)), and the list includes the imprimitive representative (12, 3) (gcd 3) — imprimitive classes are part of the enumeration.

Normalized seed (deterministic). Each class is represented by the class element with minimal |y|, tie-broken by y >= 0 then x > 0. This is well defined: two solutions with the same y satisfy x^2 = N + D*y^2, so they agree up to signs, and the tie-breaks select a unique element of the class. ("Canonical" may be used for this seed only with this rule stated alongside.) Note that a normalized seed may have x < 0 — for D = 13, N = 27 two of the four seeds do ((-12, 3) and (-40, 11); verified) — which is why the scan below collects both signs of x.

as_dict()

JSON-serializable, exact tagged values via serialize.encode_exact, keys exactly:

{"D": ..., "N": ..., "fundamental": [..., ...], "seeds": [[..., ...], ...],
 "evidence": <the evidence record's to_dict()>}

The evidence record

Evidence(
    <one-line statement: what was enumerated, and that the PARI class
     enumeration and the Lemma-A/B strip scan agreed>,
    kind="external-computation",
    software="pari",
    version=software_version("pari"),
    command=<the executed PARI call, e.g. "qfbsolve(Qfb(1,0,-13),27,3)">,
    certificate={"seeds": [...], "fundamental": [...],
                 "crosscheck": "nagell-bound-exhaustive"},
)

certificate["seeds"] and certificate["fundamental"] carry the same exact values as the record fields; Evidence.to_dict() already encodes the integer content recursively through encode_exact. The statement may describe the scan as an independent complete method only in the exact terms of Lemmas A/B (per-class strip coverage + exhaustive strip enumeration); e.g. "all solution classes of x^2 - 13*y^2 = 27 modulo <±eps>: PARI class enumeration, in agreement with the exhaustive strip scan (complete per class by the proved strip lemmas)".

Algorithm and conventions

  1. Fundamental automorph. (u0, v0) = the fundamental solution of u^2 - D*v^2 = 1 from the existing continued-fraction machinery. Move solvers._pell_unit and solvers._plus_one_unit unchanged (docstrings and doctests included) into pell_orbits.py, and have solvers.py import them back (from .pell_orbits import _pell_unit, _plus_one_unit), so every existing call site and doctest — including from diophantine_classifier.solvers import _pell_unit — keeps working, and the layering is acyclic: solvers imports pell_orbits, never the reverse at module level.
  2. Class representatives. pari(f"qfbsolve(Qfb(1,0,{-D}),{N},3)"). Flag 3 returns representatives of all solution classes modulo the automorph group <±eps> = SO(Q, ZZ) (identification proved below), including imprimitive classes. Each representative is re-verified exactly (x^2 - D*y^2 == N in ZZ); a failing representative raises InternalConsistencyError (a PARI contract violation is never silently dropped).
  3. Certified cross-check — required, runs on every call. This is the certification, not merely a test: the scan of step 4 is a second, complete enumeration of the classes, by Lemmas A and B below, and step 5 checks exact agreement with PARI's enumeration.

The strip lemmas (proved here; these inline proofs are normative)

Throughout: alpha = x + y*sqrt(D) for a solution (x, y), conj(alpha) = x - y*sqrt(D), alpha * conj(alpha) = N; eps = u0 + v0*sqrt(D) > 1; multiplication by eps realizes the action in the completeness statement; the <±eps>-class of alpha is {±eps^k * alpha : k in ZZ}. u0 >= 2 since u0^2 = 1 + D*v0^2 >= 1 + 2 = 3.

(C1) <±eps> is the proper automorphism group of the form, so these classes are the classes qfbsolve enumerates. A matrix M = [[a, b], [c, d]] in SL2(ZZ) preserves Q(x, y) = x^2 - D*y^2 iff, comparing coefficients of Q(a*x+b*y, c*x+d*y) = Q(x, y):

a^2 - D*c^2 = 1,    a*b - D*c*d = 0,    b^2 - D*d^2 = -D.

If c = 0: a = ±1, then b = 0, d^2 = 1, det = a*d = 1 forces M = ±I. If c != 0: a != 0 (as a^2 = 1 + D*c^2 > 1), so b = D*c*d/a; substituting into the third equation gives d^2*(D*c^2 - a^2)/a^2 = -1, i.e. d^2 = a^2, so d = ±a, b = ±D*c, and det = ±(a^2 - D*c^2) = ±1 forces the + sign: M = [[a, D*c], [c, a]] with a^2 - D*c^2 = 1, which is exactly multiplication by the norm-one element a + c*sqrt(D). Finally every norm-one w = a + c*sqrt(D) in ZZ[sqrt(D)] is ±eps^k: negating, assume w > 0; choose k with eps^k <= w < eps^(k+1); then w' = w*eps^(-k) = s + t*sqrt(D) is a norm-one element with 1 <= w' < eps, so conj(w') = 1/w' in (0, 1], giving s = (w' + 1/w')/2 > 0 and t = (w' - 1/w')/(2*sqrt(D)) >= 0; if t > 0 then (s, t) is a positive solution of u^2 - D*v^2 = 1 with s + t*sqrt(D) < eps, contradicting the minimality of the fundamental solution (among positive solutions, v = sqrt((u^2-1)/D) grows with u, so least u0 = least u + v*sqrt(D)); hence t = 0, w' = 1, w = eps^k.

(C2) Exact class test (also Nagell's book relation, so the attribution matches our classes). Solutions alpha_1, alpha_2 are <±eps>-equivalent iff N divides both x1*x2 - D*y1*y2 and x2*y1 - x1*y2. Proof: alpha_1 * conj(alpha_2) = (x1*x2 - D*y1*y2) + (x2*y1 - x1*y2)*sqrt(D). If alpha_2 = w*alpha_1 with w = ±eps^k, then alpha_1*conj(alpha_2) = alpha_1*conj(alpha_1)*conj(w) = N*conj(w), both coordinates divisible by N. Conversely if alpha_1*conj(alpha_2) = N*theta with theta in ZZ[sqrt(D)], then taking norms gives Norm(theta) = 1, and alpha_2 = conj(theta)*alpha_1 with conj(theta) = ±eps^k by (C1). (Usable as an exact pairwise-inequivalence check in the fixtures; verified to agree with explicit orbit walks on 8 pairs for D = 13, N = 27.)

Lemma A (N > 0). Every <±eps>-class of solutions of x^2 - D*y^2 = N contains a representative (x, y), attaining the minimum of |y'| over the class, with

0 <= y <= v0 * sqrt( N / (2*(u0 + 1)) )       and
sqrt(N) <= |x| <= sqrt( N * (u0 + 1) / 2 ).

All four inequalities are weak, and each can hold with equality (verified attainment examples below).

Proof. Pick gamma = x + y*sqrt(D) in the class minimizing |y| (the values |y'| form a nonempty set of nonnegative integers); replacing gamma by -gamma (same class), assume y >= 0. The translates eps*gamma and eps^(-1)*gamma lie in the same class by definition, with y-parts v0*x + u0*y and u0*y - v0*x; minimality gives |v0*x + u0*y| >= y and |u0*y - v0*x| >= y. Replacing x by -x swaps this pair of quantities, so as a set {|u0*y + v0*|x||, |u0*y - v0*|x||} — and since u0*y + v0*|x| >= y always holds, the binding condition is

|u0*y - v0*|x|| >= y,   i.e.   v0*|x| <= (u0 - 1)*y   or   v0*|x| >= (u0 + 1)*y.

From the equation and D*v0^2 = u0^2 - 1: v0^2*x^2 = v0^2*N + (u0^2 - 1)*y^2. In the first case, squaring (both sides nonnegative) gives v0^2*N + (u0^2 - 1)*y^2 <= (u0 - 1)^2*y^2, i.e. v0^2*N <= -2*(u0 - 1)*y^2 <= 0 — impossible for N > 0 (also when y = 0). So the second case holds; squaring gives v0^2*N + (u0^2 - 1)*y^2 >= (u0 + 1)^2*y^2, i.e. 2*(u0 + 1)*y^2 <= v0^2*N, which is the y-bound. For x: x^2 = N + D*y^2 >= N (equality iff y = 0), and x^2 <= N + D*v0^2*N/(2*(u0 + 1)) = N*(1 + (u0 - 1)/2) = N*(u0 + 1)/2. ∎

Lemma B (N < 0). Every <±eps>-class of solutions of x^2 - D*y^2 = N contains a representative (x, y), attaining the minimum of |y'| over the class, with

sqrt(|N|/D) <= y <= v0 * sqrt( |N| / (2*(u0 - 1)) )       and
0 <= |x| <= sqrt( |N| * (u0 - 1) / 2 ).

All inequalities weak; y >= 1 automatically (y = 0 would force x^2 = N < 0), and u0 - 1 >= 1 > 0.

Proof. As above pick the minimal-|y| representative with y > 0 (sign flip via -gamma; y != 0). The same dichotomy holds: v0*|x| <= (u0 - 1)*y or v0*|x| >= (u0 + 1)*y. Now v0^2*x^2 = (u0^2 - 1)*y^2 - v0^2*|N|. The second case would give (u0^2 - 1)*y^2 - v0^2*|N| >= (u0 + 1)^2*y^2, i.e. -v0^2*|N| >= 2*(u0 + 1)*y^2 > 0 — impossible. So the first case holds: (u0^2 - 1)*y^2 - v0^2*|N| <= (u0 - 1)^2*y^2, i.e. 2*(u0 - 1)*y^2 <= v0^2*|N|, the upper y-bound. Lower: D*y^2 = x^2 + |N| >= |N| (equality iff x = 0). For x: x^2 = D*y^2 - |N| <= |N|*((u0^2 - 1)/(2*(u0 - 1)) - 1) = |N|*(u0 - 1)/2. ∎

|x| cannot be strengthened to x. The x-window bounds the absolute value only: for D = 13, N = 27 (with eps = 649 + 180*sqrt(13)), the class of (12, -3) contains exactly one element with 0 <= y <= 25, namely (-12, 3) with x < 0 — its positive-x elements have y-parts ..., -4107, -3, 213, ..., skipping the strip entirely (verified by brute-force class partition to |y| <= 6000; two of the four classes behave this way). Any claim shaped sqrt(N) <= x per class is therefore false, and the scan must collect both signs of x.

Endpoint behavior (verified; pins the weak inequalities and the inclusive loop). D = 3, N = 6: the y-bound is exactly 1 (v0^2*N = 6 = 2*(u0+1)), attained by the single class's strip elements (±3, 1), whose |x| = 3 also attains sqrt(N*(u0+1)/2) = 3. D = 5, N = 4: (±2, 0) attain y = 0 and |x| = sqrt(N) = 2. D = 3, N = -2: y-bound exactly 1, attained by (±1, 1), whose |x| = 1 attains sqrt(|N|*(u0-1)/2) = 1. D = 13, N = -13: (0, 1) attains x = 0 and the lower endpoint y = sqrt(|N|/D) = 1. Because equality occurs, every bound must be implemented weakly and the scan loop must include y = Y.

Attribution (not load-bearing). The classical origin of Lemmas A/B is T. Nagell, Introduction to Number Theory, John Wiley & Sons, New York / Almqvist & Wiksell, Stockholm, 1951; 2nd ed., Chelsea, New York, 1964 — the bounds on fundamental solutions of the classes of x^2 - D*y^2 = ±N. Exact theorem numbers could not be verified against an open copy while writing this revision and are deliberately omitted. The statements are independently cross-checked against the literature: R. Boumahdi, O. Kihel and S. Mavecha, Proof of the conjecture of Keskin, Siar and Karaatli (arXiv:1601.04045), Lemmas 2.1–2.2, quoted verbatim with the source open: for u + v*sqrt(d) the fundamental solution of a class K of u^2 - d*v^2 = N (N a positive integer) and x1 + y1*sqrt(d) the fundamental solution of x^2 - d*y^2 = 1, "0 <= v <= (y1/sqrt(2(x1+1)))*sqrt(N), and 0 < |u| <= sqrt((1/2)(x1+1)N)"; and for u^2 - d*v^2 = -N, "0 < v <= (y1/sqrt(2(x1-1)))*sqrt(N), and 0 <= |u| <= sqrt((1/2)(x1-1)N)" — for the proofs the paper says "see [5]", where [5] is the Nagell 1951 book. Note the literature statement bounds |u|, exactly as Lemmas A/B bound |x|. The inline proofs above are what the implementation and tests rely on; the citations are attribution.

  1. The scan, in exact arithmetic only. y runs over 0 <= y <= Y with

    Y = isqrt( (v0^2 * |N|) // (2*(u0 + 1)) )   for N > 0,
Y = isqrt( (v0^2 * |N|) // (2*(u0 - 1)) )   for N < 0.

    This integer loop bound is the exact floor of the real Lemma-A/B bound — no floating point enters the certification, and no slack is needed: for integers a >= 0, b > 0 and an integer k >= 0, k <= sqrt(a/b) iff k^2 <= a/b iff k^2*b <= a iff k^2 <= a // b; hence isqrt(a // b) = floor(sqrt(a/b)), and since the scanned y are integers, y <= floor(bound) is equivalent to y <= bound — the certified rounding of the real bound to an integer loop bound loses nothing. (For N < 0 the strip [0, Y] is a superset of the proved window [sqrt(|N|/D), Y]; scanning from 0 is harmless and kept for simplicity.) For each y in the strip, test whether N + D*y^2 is a perfect square x^2 (exact isqrt); collect all (±x, y) — both signs, per the |x| note above.

  2. Certification. Normalize both the qfbsolve representatives and the scan hits to their classes' normalized seeds and compare the two seed sets; inequality raises InternalConsistencyError. Normalization walks the orbit with eps^{±1} while |y| strictly decreases, then applies the tie-breaks: along an orbit, y_k = A*eps^k - B*eps^{-k} (real A, B determined by the seed), so |y_k| is unimodal in k with a unique minimum region, and the greedy walk reaches the class minimum. Every class has a scan hit (Lemmas A/B: the guaranteed representative has 0 <= y <= Y, and the scan collects it whatever the sign of its x) and every scan hit lies in a class, so set equality is exactly the statement that the two independent complete enumerations agree. On success, seeds is the agreed set sorted by (|y|, y, x). (The exact divisibility test (C2) may additionally be used in tests as a second equivalence check.)

  3. Docstrings. Every function, private helpers included, gets a Sage-convention docstring with INPUT/OUTPUT where nontrivial and EXAMPLES that pass sage -t; doctests import what they need explicitly; docstring coverage stays 100%.

Failure and fallback semantics

Acceptance criteria

All numeric expectations below are verified against Sage 10.7 / PARI 2.17.2 (fundamental solutions, qfbsolve representative lists, class partitions by brute force to |y| <= 6000, strip bounds, and normalized seeds).

  1. seed_orbits(13, 27): fundamental == (649, 180); the qfbsolve representatives are [(-40, -11), (-40, 11), (12, -3), (12, 3)] (including the imprimitive (12, 3)); len(seeds) == 4; seeds == ((-12, 3), (12, 3), (-40, 11), (40, 11)) (sorted by (|y|, y, x); the multiset {|y|} is {3, 3, 11, 11}); the scan strip is 0 <= y <= 25 (assert Y == 25; isqrt(874800 // 1300) = isqrt(672) = 25); each of the four classes has exactly one strip element, and for the classes of (-12, 3) and (-40, 11) that element has x < 0 — asserting the scan collects both signs of x.
  2. seed_orbits(13, -27) (the N < 0 mirror fixture): fundamental == (649, 180); qfbsolve representatives [(5, -2), (5, 2), (21, -6), (21, 6)]; len(seeds) == 4; seeds == ((-5, 2), (5, 2), (-21, 6), (21, 6)); Y == 25 (isqrt(874800 // 1296) = isqrt(675) = 25, and 874800 = 1296 * 675 exactly); every strip hit has y in {2, 6}, consistent with the proved lower endpoint y >= sqrt(27/13) > 1.
  3. seed_orbits(3, 1): fundamental == (2, 1); seeds == ((1, 0),) (strip Y == 0; the two y = 0 hits (±1, 0) are one class).
  4. seed_orbits(2, 7): fundamental == (3, 2); qfbsolve representatives [(-3, -1), (-3, 1)]; seeds == ((-3, 1), (3, 1)) (both |y| == 1).
  5. seed_orbits(5, 4): fundamental == (9, 4); qfbsolve representatives [(-3, -1), (-3, 1), (2, 0)]; seeds == ((2, 0), (-3, 1), (3, 1)); (2, 0) is its own normalized seed (|y| = 0 minimal, x > 0) and attains both lower endpoints y = 0, |x| = sqrt(N).
  6. Endpoint fixtures (equality cases; pin weak bounds and the inclusive y = Y loop iteration): seed_orbits(3, 6): fundamental == (2, 1), one class, seeds == ((3, 1),), Y == 1 with the real Lemma-A bound exactly 1 (v0^2*N = 6 = 2*(u0+1)), strip hits (±3, 1) attaining y = Y and |x| = 3 = sqrt(N*(u0+1)/2). seed_orbits(3, -2): one class, seeds == ((1, 1),), Y == 1 with the real Lemma-B bound exactly 1, hits (±1, 1) attaining y = Y and |x| = 1 = sqrt(|N|*(u0-1)/2). seed_orbits(13, -13): seeds == ((0, 1),), Y == 18, the hit (0, 1) attaining x = 0 and the lower endpoint y = sqrt(|N|/D) = 1.
  7. Empty case: seed_orbits(3, -1) has seeds == () (x^2 - 3*y^2 = -1 has no integer solutions; qfbsolve returns the empty list and the strip scan finds no hit — both verified; the existing doctest _pell_solutions_stream(3, -1) is None pins the same fact); the evidence record is still present with certificate["crosscheck"] == "nagell-bound-exhaustive".
  8. ValueError: seed_orbits(4, 3) (square D), seed_orbits(-3, 5) (nonpositive D), seed_orbits(13, 0) (zero N) each raise ValueError.
  9. Evidence: kind == "external-computation", software == "pari", version == software_version("pari"), command is the executed qfbsolve call; certificate has exactly the keys {"seeds", "fundamental", "crosscheck"}; as_dict() round-trips through json.dumps/json.loads unchanged, with every integer appearing only as a tagged exact value.
  10. The cross-check runs on every call: a unit test that monkeypatches the qfbsolve wrapper to return a corrupted representative list asserts InternalConsistencyError.
  11. The helper move leaves solvers.py fully working: all existing doctests (including the moved ones) pass under sage -t, and no module-level import cycle exists.
  12. make test, make doctest, make coverage (docstring coverage stays 100%), make registry-docs clean.

Out of scope

References

The completeness of the strip scan rests on Lemmas A and B proved inline above; T. Nagell, Introduction to Number Theory (Wiley/Almqvist & Wiksell 1951; 2nd ed., Chelsea 1964) is cited as the classical origin of those bounds, and Boumahdi–Kihel–Mavecha (arXiv:1601.04045, Lemmas 2.1–2.2, quoting Nagell's book for the proofs) as the open cross-check of the transcription — attribution in both cases, not load-bearing pointers. No data/references.bib change is required by this child: the enumeration evidence carries its cross-check inside the certificate payload rather than as a bibliography reference, and the epic's bibliography additions (BakerDavenport1969, DujellaPetho1998) are cited by #81/#83. Anglin1996 and Bennett1998 remain the family's references.

Metadata

Metadata

Assignees

No one assigned

    Labels

    enhancementNew feature or request

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions