Family: pythagorean — Pythagorean triples

diophantine_classifier/data/families/pythagorean.yaml

@@ -0,0 +1,20 @@
+# Family: Pythagorean triples
+# Schema: see data/families/README.md; prose documentation in docs/FAMILIES.md.
+# References are keys into data/references.bib, each with a 'why' annotation.
+slug: pythagorean
+name: Pythagorean triples
+priority: 1
+status: solved
+class: quadratic
+form: x^2 + y^2 = z^2
+parents:
+- legendre
+matcher: true
+methods:
+- Euclid parametrization
+notes: Primitive solutions (m^2-n^2, 2mn, m^2+n^2), gcd(m,n)=1, m != n mod 2.
+examples:
+- x^2 + y^2 = z^2
+references:
+- key: Dickson1920
+  why: history and the classical parametrization of rational right triangles

diophantine_classifier/solvers.py

@@ -1052,6 +1052,38 @@ def _solve_legendre(cls, match):
                           gram=[[a, 0, 0], [0, b, 0], [0, 0, c]])
 
 
+def _solve_pythagorean(cls, match):
+    r"""
+    Pythagorean triples: Euclid's parametrization as a stream.
+
+    Iteration yields the primitive triples ``(m^2 - n^2, 2mn, m^2 + n^2)``
+    for ``m > n >= 1`` coprime of opposite parity, ordered by ``m``; the full
+    solution set consists of their multiples, sign changes, and leg swaps.
+
+    EXAMPLES::
+
+        sage: from diophantine_classifier import solve
+        sage: S = solve("x^2 + y^2 = z^2")
+        sage: S.first(4)
+        [(3, 4, 5), (5, 12, 13), (15, 8, 17), (7, 24, 25)]
+    """
+    def stream():
+        # (leg, leg, hypotenuse) in the family's standard coordinates; which
+        # of the user's variables those are is the transform's business
+        for m in itertools.count(2):
+            for n in range(1, m):
+                if gcd(m, n) == 1 and (m - n) % 2 == 1:
+                    yield (m * m - n * n, 2 * m * n, m * m + n * n)
+
+    first = list(itertools.islice(stream(), 4))
+    return SolutionSet(
+        _normalized(match), first, "parametrized",
+        "primitive solutions (m^2 - n^2, 2mn, m^2 + n^2), gcd(m, n) = 1, "
+        "m ≢ n (mod 2); all solutions are multiples, sign changes and swaps "
+        "of these; iteration enumerates the primitive triples",
+        complete=True, stream=stream)
+
+
 def _solve_weierstrass(cls, match):
     r"""
     Integral points on a Weierstrass model via ``E.integral_points``.
@@ -1161,6 +1193,7 @@ def rec(k_left, target, minimum, acc):
     "binary-qf-representation": _solve_bqf,
     "quadratic-form-zero": _solve_qf_zero,
     "legendre": _solve_legendre,
+    "pythagorean": _solve_pythagorean,
     "elliptic-weierstrass": _solve_weierstrass,
     "thue": _solve_thue,
     "egyptian-fractions": _solve_egyptian,

docs/FAMILIES.md

@@ -99,6 +99,13 @@ PARI: `qfsolve`; Magma: `IsLocallySolvable`, `HasRationalPoint`.
 **References.** Legendre 1785; Holzer 1950; Cremona–Rusin, "Efficient solution of
 rational conics" (Math. Comp. 72, 2003); Simon 2005 (PARI `qfsolve`).
 
+### `pythagorean` — Pythagorean triples — P1, solved
+**Form.** x² + y² = z².
+**Status.** Completely parametrized: primitive solutions (m²−n², 2mn, m²+n²),
+gcd(m,n)=1, m ≢ n (mod 2). The template example of "solved by parametrization"
+(genus 0 with a rational point).
+**References.** Euclid, *Elements* X.29; Dickson, *History*, vol. II ch. IV.
+
 ### `quadratic-form-zero` — Isotropy of a quadratic form — P1, algorithmic
 **Form.** Q(x₁,…,x_k) = 0, Q a nondegenerate integral quadratic form, k ≥ 3.
 **Status.** Hasse–Minkowski: solvable iff solvable over ℝ and all ℚ_p (finite

tests/test_classify.py

@@ -18,6 +18,7 @@
     ("2*x^2 + 3*x*y - 5*y^2 + x - 7 = 0", "", "binary-quadratic"),
     ("3*x^2 + 7*y^2 = 19", "", "binary-qf-representation"),
     # quadratic, more variables
+    ("x^2 + y^2 = z^2", "", "pythagorean"),
     ("x^2 + 3*y^2 = 7*z^2", "", "legendre"),
     ("x^2 - 3*y^2 + 5*z^2 - 7*w^2 = 0", "", "quadratic-form-zero"),
     ("x^2 + x*y + y^2 + z^2 = 14", "", "quadratic-form-representation"),

tests/test_iterators.py

@@ -19,6 +19,15 @@ def test_linear_stream():
     assert all(3 * x + 5 * y == 1 for x, y in sols)
 
 
+def test_pythagorean_stream():
+    s = solve("x^2 + y^2 = z^2")
+    sols = s.first(6)
+    assert sols[0] == (3, 4, 5)
+    for x, y, z in sols:
+        assert x ** 2 + y ** 2 == z ** 2
+        assert gcd(gcd(x, y), z) == 1     # primitive
+
+
 def test_egyptian_complete():
     s = solve("1/x + 1/y + 1/z = 1")
     assert s.solutions == [(2, 3, 6), (2, 4, 4), (3, 3, 3)]

tests/test_solvers.py

@@ -207,6 +207,14 @@ def capture(cls, match):
     assert "coeffs" not in seen["data"]    # ... not the linear match's
 
 
+def test_pythagorean_roles_are_transported_back():
+    S = solve("a^2 + b^2 = c^2")
+    assert S.variables == ("a", "b", "c")
+    sols = S.first(4)
+    assert all(a ** 2 + b ** 2 == c ** 2 for a, b, c in sols)
+    assert_valid_solutions("a^2 + b^2 = c^2", sols)
+
+
 def test_component_of_a_rational_equation_keeps_the_parent_conditions():
     """(x - 2)(x - y)/y = 0 factors; y != 0 still binds each component."""
     cls = classify("(x - 2)*(x - y)/y = 0")