Parker Square: searching for a magic square of squares

Status: Numerically checked. No magic square of nine squares exists with a centre of m² for any m up to 1,000,000. We ran the code ourselves. It found nothing new: the problem is still open.

The problem

Can you fill a 3×3 grid with nine different square numbers so that every row, column and both diagonals add up to the same number? Nobody knows. The question has been open since the 1980s. Matt Parker's famous attempt, the "Parker Square", uses nine squares but fails on one diagonal.

The idea

Fable 5.1 (High) did not try all combinations of nine numbers. It used the structure every 3×3 magic square has (the Lucas form):

c+a     c-a-b   c+b
c-a+b   c       c+a-b
c-b     c+a+b   c-a

The nine entries are c and c ± d for d in {a, b, a+b, a−b}. If all nine are squares:

1. The centre must be a square: c = m².
2. For each d, c + d = p² and c − d = q². With u = (p−q)/2 and v = (p+q)/2 this becomes u² + v² = m² and d = 2uv.
3. So each of a, b, a+b and a−b must come from a Pythagorean triple with hypotenuse m.

The search therefore only looks at values of m with at least four Pythagorean triples and checks whether two values x > y from S = {2uv} also have x+y and x−y in S. We checked the algebra: it is correct.

Result

The one primitive 7-square magic square it found is a known one, with centre 425²:

565²     289²     373²
 23²     425²   360721
222121   527²     205²

Review by a second model

Claude Opus 5.5 checked the method and ran the code.

• The maths is correct. The reduction from magic squares to Pythagorean triples holds, and the search for full solutions is complete inside its range.
• Nothing new was found. No solution is expected in this range: earlier published searches have gone far beyond 10¹². The 7-square square with centre 425² is a known result from the late 1990s.
• The 2,352 near misses are mostly copies. 850², 1275², 1700² … are the 425² square multiplied by 2, 3, 4. Only one is primitive.
• The near-miss search is incomplete. It only looks at m with at least four triples and only at squares whose centre is a square, so it can miss other 7-square squares. This does not affect the search for full solutions.

Verdict: a clean, correct and fast search that confirms what is known. Status 2, numerically checked, for the statement "no solution with centre m², m ≤ 1,000,000".

How it was made

• Model: Claude Fable 5.1 (High), Anthropic.
• Check: run again and reviewed by Claude Opus 5.5.
• Tools: Python 3, standard library only.

The code

Anyone can run it: python parker.py 1000000.

"""
Suche nach einem 3x3 magischen Quadrat aus Quadratzahlen ("Parker Square"-Problem).

Idee (keine Brute-Force über 9 Zahlen, sondern Struktur ausnutzen):

Jedes 3x3 magische Quadrat hat die Lucas-Form

        c+a     c-a-b   c+b
        c-a+b   c       c+a-b
        c-b     c+a+b   c-a

mit magischer Summe 3c. Die 9 Einträge sind c und c±d für d in {a, b, a+b, a-b}.

Sollen alle Einträge Quadrate sein, muss c = m² sein. Für jede Differenz d gilt
    c+d = p²,  c-d = q²   =>   p² + q² = 2m²
und mit u = (p-q)/2, v = (p+q)/2 wird daraus
    u² + v² = m²   und   d = 2uv.

Jede der vier Differenzen a, b, a+b, a-b muss also 2·u·v für ein pythagoreisches
Tripel (u, v, m) mit Hypotenuse m sein. Suche daher:

    m mit >= 4 Tripeln,  S = {2uv},  x > y in S mit  x+y in S  und  x-y in S.

Falls nur x+y in S (aber x-y nicht): magisches Quadrat mit 7 Quadrat-Einträgen
("Near Miss").
"""

import sys
from math import gcd, isqrt
from collections import defaultdict

N = int(sys.argv[1]) if len(sys.argv) > 1 else 1_000_000


def triples_by_hypotenuse(limit):
    """Alle pythagoreischen Tripel (u, v, m) mit m <= limit, gebucketet nach m."""
    buckets = defaultdict(list)
    pmax = isqrt(limit) + 1
    for p in range(2, pmax):
        for q in range(1, p):
            if (p - q) % 2 == 0 or gcd(p, q) != 1:
                continue
            h = p * p + q * q
            if h > limit:
                break
            u, v = p * p - q * q, 2 * p * q
            for k in range(1, limit // h + 1):
                buckets[k * h].append((k * u, k * v))
    return buckets


def build_square(c, a, b):
    return [
        [c + a, c - a - b, c + b],
        [c - a + b, c, c + a - b],
        [c - b, c + a + b, c - a],
    ]


def is_square(n):
    return n >= 0 and isqrt(n) ** 2 == n


def check(sq):
    """Gibt (ist_magisch, anzahl_quadrate, alle_verschieden) zurück."""
    s = sum(sq[0])
    lines = [sq[i] for i in range(3)] + [[sq[i][j] for i in range(3)] for j in range(3)]
    lines.append([sq[i][i] for i in range(3)])
    lines.append([sq[i][2 - i] for i in range(3)])
    magic = all(sum(l) == s for l in lines)
    flat = [x for row in sq for x in row]
    return magic, sum(is_square(x) for x in flat), len(set(flat)) == 9


def fmt(sq):
    out = []
    for row in sq:
        cells = []
        for x in row:
            cells.append(f"{isqrt(x)}²" if is_square(x) else f"{x}")
        out.append("  ".join(f"{c:>12}" for c in cells))
    return "\n".join(out)


def main():
    print(f"Enumeriere pythagoreische Tripel bis Hypotenuse m = {N:,} ...", flush=True)
    buckets = triples_by_hypotenuse(N)
    print(f"  {sum(len(v) for v in buckets.values()):,} Tripel, "
          f"{len(buckets):,} verschiedene Hypotenusen", flush=True)

    candidates = {m: t for m, t in buckets.items() if len(t) >= 4}
    print(f"  {len(candidates):,} Werte von m haben >= 4 Tripel (nur die können funktionieren)\n", flush=True)

    # Volle Lösung: alle vier Differenzen a, b, a+b, a-b liegen in S.
    # Near Miss (7 Quadrate): zwei Differenzen "voll" (c+d und c-d Quadrat),
    # bei den anderen beiden ist nur eine Seite ein Quadrat. Dafür probieren wir
    # pro Paar (x, y) aus S alle Rollen durch und zählen einfach die Quadrate.
    full = []
    seven = []
    for m in sorted(candidates):
        c = m * m
        S = sorted({2 * u * v for u, v in candidates[m]})
        Sset = set(S)
        for i, y in enumerate(S):
            for x in S[i + 1:]:
                if (x + y) in Sset and (x - y) in Sset and x != 2 * y:
                    full.append((m, x, y))
                    continue
                # mögliche Rollen von (x, y) unter {a, b, a+b, a-b}
                roles = [(x, y), (y, x - y), (x, y - x)]
                if (x + y) % 2 == 0:
                    roles.append(((x + y) // 2, (x - y) // 2))
                for a, b in roles:
                    sq = build_square(c, a, b)
                    if min(min(r) for r in sq) <= 0:
                        continue
                    magic, nsq, distinct = check(sq)
                    if magic and distinct and nsq >= 7:
                        seven.append((m, a, b, nsq))

    print(f"Volle Lösungen (9 Quadrate):           {len(full)}")
    print(f"Near Misses (>=7 Quadrate, magisch):   {len(seven)}\n")

    for m, a, b in full:
        sq = build_square(m * m, a, b)
        print(f"!!! LÖSUNG GEFUNDEN  m={m}  a={a}  b={b}  check={check(sq)}")
        print(fmt(sq), "\n")

    print("Kleinste Near Misses (7 von 9 Einträgen sind Quadrate):")
    for m, a, b, nsq in seven[:5]:
        sq = build_square(m * m, a, b)
        print(f"\n m = {m}  (Zentrum {m}² = {m*m}, magische Summe {3*m*m}, {nsq} Quadrate)")
        print(fmt(sq))
    eight = [t for t in seven if t[3] == 8]
    print(f"\nQuadrate mit genau 8 Quadrat-Einträgen: {len(eight)}")


if __name__ == "__main__":
    main()

Next steps

1. Rewrite the search in C or Rust and push the bound far beyond 10¹².
2. Extend the near-miss search to all values of m, to find every 7-square square in the range.
3. Look for an 8-square magic square, which nobody has found yet.