Family: sum-of-four-squares — Sum of four squares (Lagrange)

diophantine_classifier/data/families/sum-of-four-squares.yaml

@@ -0,0 +1,22 @@
+# Family: Sum of four squares (Lagrange)
+# 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: sum-of-four-squares
+name: Sum of four squares (Lagrange)
+priority: 1
+status: solved
+class: quadratic
+form: x^2 + y^2 + z^2 + w^2 = n
+parents:
+- quadratic-form-representation
+matcher: true
+methods:
+- always solvable (Lagrange 1770)
+- Rabin-Shallit
+software:
+  sage: four_squares(n)
+code:
+  sage: four_squares({n})
+references:
+- key: RabinShallit1986
+  why: randomized polynomial-time four-squares algorithm

diophantine_classifier/solvers.py

@@ -943,6 +943,48 @@ def stream():
                        complete=True, stream=stream)
 
 
+def _solve_four_squares(cls, match):
+    r"""
+    All representations ``n = x^2 + y^2 + z^2 + w^2``, ``x <= y <= z <= w``.
+
+    EXAMPLES::
+
+        sage: from diophantine_classifier import solve
+        sage: solve("x^2 + y^2 + z^2 + w^2 = 7").solutions
+        [(1, 1, 1, 2)]
+    """
+    n = _zz(match.data, "n")
+    if n is None:
+        raise SolverUnavailable("needs concrete n")
+    if n < 0:
+        return SolutionSet(_normalized(match), [], "empty", "n < 0",
+                           complete=True)
+    if n > MAX_FOUR_SQUARES:
+        sol = four_squares(n)
+        return SolutionSet(_normalized(match), [tuple(sol)], "witness",
+                           "one representation (Lagrange: always solvable; "
+                           "n too large for full enumeration)",
+                           complete=False)
+    sols = []
+    x = 0
+    while 4 * x * x <= n:
+        y = x
+        while x * x + 3 * y * y <= n:
+            z = y
+            while x * x + y * y + 2 * z * z <= n:
+                w2 = n - x * x - y * y - z * z
+                w = isqrt(w2)
+                if w * w == w2 and w >= z:
+                    sols.append((ZZ(x), ZZ(y), ZZ(z), ZZ(w)))
+                z += 1
+            y += 1
+        x += 1
+    return SolutionSet(_normalized(match), sols, "finite-complete",
+                       "all representations with 0 <= x <= y <= z <= w; the "
+                       "rest differ by signs and order (Lagrange: always "
+                       "solvable)", complete=True)
+
+
 def _solve_bqf(cls, match):
     r"""
     Representations by a binary quadratic form.
@@ -1158,6 +1200,7 @@ def rec(k_left, target, minimum, acc):
     "univariate": _solve_univariate,
     "linear": _solve_linear,
     "pell-like": _solve_pell_like,
+    "sum-of-four-squares": _solve_four_squares,
     "binary-qf-representation": _solve_bqf,
     "quadratic-form-zero": _solve_qf_zero,
     "legendre": _solve_legendre,

docs/FAMILIES.md

@@ -118,6 +118,13 @@ uniform results: 15-theorem (Conway–Schneeberger–Bhargava), 290-theorem
 `qfsolve`; Magma: `RepresentationNumber`, ternary form machinery.
 **References.** Cassels; Bhargava 2000; Bhargava–Hanke 2005.
 
+### `sum-of-four-squares` — Four squares (Lagrange) — P1, solved
+**Form.** x² + y² + z² + w² = n.
+**Status.** Always solvable (Lagrange 1770); Jacobi's formula counts representations;
+randomized polynomial-time algorithms.
+**Software.** Sage: `four_squares(n)`.
+**References.** Lagrange 1770; Jacobi 1834; Rabin–Shallit 1986.
+
 ### `quadric` — General quadratic Diophantine equation — P2, algorithmic
 **Form.** Q(x₁,…,x_k) + L(x₁,…,x_k) + c = 0 (arbitrary quadratic, k ≥ 3).
 **Status.** Decidable in general — the deepest case of the quadratic theory

tests/test_classify.py

@@ -19,6 +19,7 @@
     ("3*x^2 + 7*y^2 = 19", "", "binary-qf-representation"),
     # quadratic, more variables
     ("x^2 + 3*y^2 = 7*z^2", "", "legendre"),
+    ("x^2 + y^2 + z^2 + w^2 = n", "n", "sum-of-four-squares"),
     ("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"),
     ("x^2 + y^2 - z^2 + 3*x - 7 = 0", "", "quadric"),

tests/test_solvers.py

@@ -33,6 +33,12 @@ def test_pell_like():
     assert s.complete    # orbit representatives + automorph action
 
 
+def test_four_squares():
+    s = solve("x^2 + y^2 + z^2 + w^2 = 7")
+    assert s.solutions == [(1, 1, 1, 2)]
+    assert s.complete
+
+
 def test_bqf():
     s = solve("3*x^2 + 7*y^2 = 19")
     assert s.solutions == [(-2, -1), (-2, 1), (2, -1), (2, 1)]