Problem 983: Consonant Circle Crossing

Project Euler Problem 983 Solution

EulerSolve provides an optimized solution for Project Euler Problem 983, Consonant Circle Crossing, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Problem 983 asks for the smallest squared radius $R(n)^2$ that allows a sufficiently large consonant family of equal circles. The implementations realize the configuration on the integer lattice: every chosen direction is a lattice point on the circle $x^2+y^2=R(n)^2$, the actual circle centers are built from subset sums of those directions, and every nontrivial crossing of the chosen circles must occur at one of the designated crossing points produced by the same subset-sum construction. The key simplification is that the code does not search directly for $n$ unrelated circles. Instead it builds two parity classes of subset sums, one class for circle centers and one class for allowed crossings. If $d$ directions are chosen, each parity class contains $2^{d-1}$ points, so for $n>2$ the search dimension is $$d=\min\{t\ge 1:2^{t-1}\ge n\}.$$ The answer is then found by increasing $m=1,2,3,\dots$ until the circle $x^2+y^2=m$ supports such a $d$-direction construction. Mathematical Approach The whole method is built around equal-radius lattice vectors, subset sums, and a sharp characterization of which circle crossings are legitimate. From $n$ circles to $d$ lattice directions For $n\le 2$, the problem is trivial and the implementations immediately return $R(n)^2=1$....

Detailed mathematical approach

Problem Summary

Problem 983 asks for the smallest squared radius $R(n)^2$ that allows a sufficiently large consonant family of equal circles. The implementations realize the configuration on the integer lattice: every chosen direction is a lattice point on the circle $x^2+y^2=R(n)^2$, the actual circle centers are built from subset sums of those directions, and every nontrivial crossing of the chosen circles must occur at one of the designated crossing points produced by the same subset-sum construction.

The key simplification is that the code does not search directly for $n$ unrelated circles. Instead it builds two parity classes of subset sums, one class for circle centers and one class for allowed crossings. If $d$ directions are chosen, each parity class contains $2^{d-1}$ points, so for $n>2$ the search dimension is

$$d=\min\{t\ge 1:2^{t-1}\ge n\}.$$

The answer is then found by increasing $m=1,2,3,\dots$ until the circle $x^2+y^2=m$ supports such a $d$-direction construction.

Mathematical Approach

The whole method is built around equal-radius lattice vectors, subset sums, and a sharp characterization of which circle crossings are legitimate.

From $n$ circles to $d$ lattice directions

For $n\le 2$, the problem is trivial and the implementations immediately return $R(n)^2=1$. For larger $n$, the construction uses $d$ selected vectors $v_1,\dots,v_d$, all with the same squared length $m$. From them it forms the two parity classes

$$E=\left\{\sum_{i\in A} v_i:\ |A|\equiv 0 \pmod 2\right\},\qquad O=\left\{\sum_{i\in A} v_i:\ |A|\equiv 1 \pmod 2\right\}.$$

Here $E$ is the family of circle centers actually counted by the search, while $O$ is the family of crossing points that those circles are allowed to share. Both sets have size $2^{d-1}$, which is why the least admissible $d$ is the smallest one with $2^{d-1}\ge n$.

Lattice points on a fixed radius

Fix $m$ and write

$$\Lambda_m=\{(x,y)\in\mathbb{Z}^2:x^2+y^2=m\}.$$

Every chosen direction must come from $\Lambda_m$. Because $v$ and $-v$ generate the same geometric direction up to parity reversal, the search keeps only one representative from each antipodal pair $\{v,-v\}$. Therefore $m$ is immediately impossible if $\Lambda_m$ contains fewer than $d$ antipodal pairs.

The factorization test used in the implementations is the standard two-squares criterion. If

$$m=2^a\prod_i p_i^{\alpha_i}\prod_j q_j^{2\beta_j},\qquad p_i\equiv 1 \pmod 4,\quad q_j\equiv 3 \pmod 4,$$

then

$$|\Lambda_m|=4\prod_i(\alpha_i+1),\qquad \text{antipodal pairs}=\frac{|\Lambda_m|}{2}=2\prod_i(\alpha_i+1).$$

If any prime $q\equiv 3 \pmod 4$ appears to an odd power, then $\Lambda_m$ is empty and the search skips $m$ at once.

Why subset sums create the intended circle incidences

Take any even subset sum $e\in E$. If the subset defining $e$ does not contain $i$, then $e+v_i\in O$; if it does contain $i$, then $e-v_i\in O$. Since every $v_i$ lies on $x^2+y^2=m$, each point of $O$ produced this way lies on the circle of radius $\sqrt m$ centered at $e$. So every circle centered at a point of $E$ comes with $d$ intended crossing points from $O$.

The same relation can be read backwards: every odd subset sum $o\in O$ differs from $d$ even subset sums by vectors in $\Lambda_m$, so each odd point lies on $d$ circles of the family. The subset-sum construction therefore creates a bipartite incidence pattern isomorphic to the $d$-dimensional cube: toggling one chosen vector changes parity and moves by exactly one lattice vector of length $\sqrt m$.

Injectivity inside the parity classes

The implementations require the map $A\mapsto \sum_{i\in A}v_i$ to be injective separately on even subsets and on odd subsets. Equivalently, they reject any nontrivial relation

$$\sum_{i\in A}v_i=\sum_{i\in B}v_i\quad \text{with}\quad |A|\equiv |B|\pmod 2,\ A\ne B.$$

If such a collision existed inside $E$, two circles that are supposed to be different would collapse to the same center. If it existed inside $O$, two designated crossings would merge. This injectivity check is the first exact validation step performed after the directions have been chosen.

Translated circles and forbidden extra crossings

The deeper condition is that no point outside $O$ may lie on two circles centered in $E$. For an even center $e\in E$ and a lattice vector $u\in \Lambda_m$, consider the point $z=e+u$, which certainly lies on the circle centered at $e$. If $z\in O$, then this is one of the intended crossings and it is allowed to belong to several circles. If $z\notin O$, then the translated circle $z+\Lambda_m$ must not hit a second point of $E$. The exact condition checked by the code is

$$z\notin O\ \Longrightarrow\ |(z+\Lambda_m)\cap E|<2.$$

This is the decisive consonance test: every genuine multiple crossing of the chosen circles must come from the designated odd subset sums, and nowhere else.

The displacement set used for pruning

To reject hopeless branches early, the implementations precompute

$$\Delta_m=\{a-b:\ a,b\in \Lambda_m,\ a\ne b\}.$$

If two even centers lie on the same translated radius-$\sqrt m$ circle, then their difference belongs to $\Delta_m$. Differences of even subset sums are signed sums of an even number of chosen directions. Two-term differences correspond to intended odd crossings, so the first genuinely dangerous obstructions arise from four selected directions. During depth-first search, every newly created 4-tuple is tested through the sixteen sums

$$\pm v_i\pm v_j\pm v_k\pm v_\ell,$$

and the branch is pruned as soon as a nonzero value lands in $\Delta_m$. This is not the whole proof of validity, but it is an extremely effective necessary filter before the full subset-sum test.

Worked Example: the checkpoint $R(4)^2=5$

For $n=4$, the implementations need $d=3$ because $2^{3-1}=4$. Take $m=5$, for which

$$\Lambda_5=\{(\pm 2,\pm 1),(\pm 1,\pm 2)\}.$$

Choose three representatives, for example

$$v_1=(2,1),\qquad v_2=(2,-1),\qquad v_3=(1,2).$$

Then the even subset sums are

$$E=\{(0,0),(4,0),(3,3),(3,1)\},$$

and the odd subset sums are

$$O=\{(2,1),(2,-1),(1,2),(5,2)\}.$$

So we obtain four circles of radius $\sqrt 5$, centered at the points of $E$. The circle centered at $(0,0)$ passes through $(2,1)$, $(2,-1)$, and $(1,2)$; the circle centered at $(4,0)$ passes through $(2,1)$, $(2,-1)$, and $(5,2)$; and similarly for the other two even centers. The odd points are exactly the intended shared crossings, and the full validation confirms that no other point lies on two of these four circles. That is why the checkpoint value is $R(4)^2=5$.

How the Code Works

Arithmetic filters before the search

The C++, Python, and Java implementations first convert $n$ into the required dimension $d$, with the trivial shortcut $R(n)^2=1$ for $n\le 2$. They then scan $m=1,2,3,\dots$. For each $m$, they use the prime-factorization formula above to count antipodal pairs quickly, and discard $m$ immediately if there are fewer than $d$ available directions.

Choosing representatives and pruning the DFS

When $m$ survives the arithmetic filter, the implementations enumerate all lattice points on $x^2+y^2=m$, collapse them into antipodal pairs, and keep one canonical representative from each pair. A depth-first search then tries all $d$-element choices of those representatives. The precomputed displacement set $\Delta_m$ is used at every step: whenever a new 4-tuple already creates a forbidden signed sum, the branch is abandoned before any expensive subset-sum work is done.

The three languages use the same mathematical search. The C++ implementation additionally parallelizes the top-level DFS branches over the choice of the first representative, while the Python and Java implementations run the same search serially.

Full candidate verification and termination

Once $d$ directions have been selected, the implementations generate all $2^d$ subset sums, split them into $E$ and $O$, and test injectivity inside each parity class. If that passes, they run the translated-circle condition $z\notin O \Rightarrow |(z+\Lambda_m)\cap E|<2$. The first $m$ that passes every test is returned, so by construction it is the minimal value $R(n)^2$.

The built-in checkpoints confirm the first nontrivial cases $R(2)^2=1$ and $R(4)^2=5$; the C++ version also checks monotonicity for small $n$.

Complexity Analysis

The outer loop over $m$ is open-ended in theory, because the algorithm stops only when it reaches the first feasible squared radius. For one fixed $m$, the trial-division factorization and the enumeration of lattice points both cost $O(\sqrt m)$ in the current implementations.

Let $L=|\Lambda_m|$ and $u=L/2$. Building the displacement set $\Delta_m$ costs $O(L^2)$ time and memory. The search over direction choices is combinatorial in the worst case, up to $O\!\left(\binom{u}{d}\right)$ branches, although the 4-vector pruning removes most of them in practice. For a complete candidate, generating all subset sums costs $O(d\,2^d)$, while the translated-circle verification costs $O(2^{d-1}L^2)$.

For the actual target $n=500$, the dimension is $d=10$, so each parity class has $512$ points. The exponential part is therefore tied to a very small dimension; the practical difficulty is finding a good $m$ and pruning the combinatorial search aggressively enough.

Footnotes and References

Problem page: https://projecteuler.net/problem=983
Fermat's theorem on sums of two squares: Wikipedia - Fermat's theorem on sums of two squares
Lattice point: Wikipedia - Lattice point
Hypercube graph: Wikipedia - Hypercube graph
Backtracking: Wikipedia - Backtracking

Problem 983 source code

C++

#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <functional>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <unordered_set>
#include <vector>

using std::cerr;
using std::cout;
using std::int64_t;
using std::string;
using std::uint64_t;
using std::vector;

struct Point {
    int x;
    int y;
};

static inline int64_t pair_key(int x, int y) {
    return (static_cast<int64_t>(x) << 32) ^ static_cast<uint32_t>(y);
}

static vector<Point> circle_points(int m) {
    const int r = static_cast<int>(std::sqrt(static_cast<double>(m)));
    vector<Point> pts;
    pts.reserve(64);
    for (int x = -r; x <= r; ++x) {
        const int y2 = m - x * x;
        if (y2 < 0) {
            continue;
        }
        const int y = static_cast<int>(std::sqrt(static_cast<double>(y2)));
        if (y * y == y2) {
            pts.push_back({x, y});
            if (y != 0) {
                pts.push_back({x, -y});
            }
        }
    }
    std::sort(pts.begin(), pts.end(), [](const Point& a, const Point& b) {
        return (a.x < b.x) || (a.x == b.x && a.y < b.y);
    });
    pts.erase(std::unique(pts.begin(), pts.end(), [](const Point& a, const Point& b) {
        return a.x == b.x && a.y == b.y;
    }), pts.end());
    return pts;
}

static vector<std::array<Point, 2>> opposite_pairs(const vector<Point>& points) {
    vector<Point> sorted = points;
    std::sort(sorted.begin(), sorted.end(), [](const Point& a, const Point& b) {
        return (a.x < b.x) || (a.x == b.x && a.y < b.y);
    });

    std::unordered_set<int64_t> used;
    used.reserve(points.size() * 2 + 16);

    vector<std::array<Point, 2>> pairs;
    pairs.reserve(points.size() / 2 + 2);

    for (const Point& v : sorted) {
        const int64_t k = pair_key(v.x, v.y);
        if (used.find(k) != used.end()) {
            continue;
        }
        const Point ov{-v.x, -v.y};
        used.insert(k);
        used.insert(pair_key(ov.x, ov.y));
        pairs.push_back({v, ov});
    }

    return pairs;
}

static int opposite_pair_count_by_factorization(uint64_t m) {
    if (m == 0) {
        return 0;
    }

    uint64_t n = m;
    int pairs = 2;

    int e2 = 0;
    while ((n & 1ULL) == 0ULL) {
        n >>= 1;
        ++e2;
    }
    (void)e2;

    for (uint64_t p = 3; p * p <= n; p += 2) {
        if (n % p != 0) {
            continue;
        }
        int exp = 0;
        while (n % p == 0) {
            n /= p;
            ++exp;
        }
        if ((p & 3ULL) == 3ULL) {
            if ((exp & 1) != 0) {
                return 0;
            }
        } else if ((p & 3ULL) == 1ULL) {
            if (pairs > std::numeric_limits<int>::max() / (exp + 1)) {
                return std::numeric_limits<int>::max();
            }
            pairs *= (exp + 1);
        }
    }

    if (n > 1) {
        if ((n & 3ULL) == 3ULL) {
            return 0;
        }
        if ((n & 3ULL) == 1ULL) {
            if (pairs > std::numeric_limits<int>::max() / 2) {
                return std::numeric_limits<int>::max();
            }
            pairs *= 2;
        }
    }

    return pairs;
}

static std::unordered_set<int64_t> build_bad_displacements(const vector<Point>& points) {
    std::unordered_set<int64_t> bad;
    bad.reserve(points.size() * points.size() * 2 + 16);
    for (const Point& a : points) {
        for (const Point& b : points) {
            if (a.x == b.x && a.y == b.y) {
                continue;
            }
            bad.insert(pair_key(a.x - b.x, a.y - b.y));
        }
    }
    return bad;
}

static bool four_tuple_has_forbidden_sum(
    const vector<Point>& selected,
    int used_count,
    const std::unordered_set<int64_t>& bad_disp
) {
    if (used_count < 4) {
        return false;
    }

    const Point& d = selected[used_count - 1];
    for (int i = 0; i < used_count - 1; ++i) {
        const Point& a = selected[i];
        for (int j = i + 1; j < used_count - 1; ++j) {
            const Point& b = selected[j];
            for (int k = j + 1; k < used_count - 1; ++k) {
                const Point& c = selected[k];
                for (int mask = 0; mask < 16; ++mask) {
                    const int sx =
                        ((mask & 1) ? -a.x : a.x) +
                        ((mask & 2) ? -b.x : b.x) +
                        ((mask & 4) ? -c.x : c.x) +
                        ((mask & 8) ? -d.x : d.x);
                    const int sy =
                        ((mask & 1) ? -a.y : a.y) +
                        ((mask & 2) ? -b.y : b.y) +
                        ((mask & 4) ? -c.y : c.y) +
                        ((mask & 8) ? -d.y : d.y);
                    if (sx == 0 && sy == 0) {
                        continue;
                    }
                    if (bad_disp.find(pair_key(sx, sy)) != bad_disp.end()) {
                        return true;
                    }
                }
            }
        }
    }

    return false;
}

static bool test_selected_vectors(
    const vector<Point>& v,
    const vector<Point>& lattice
) {
    const int d = static_cast<int>(v.size());
    const int full = 1 << d;

    vector<Point> sums(full, Point{0, 0});
    for (int mask = 1; mask < full; ++mask) {
        const int bit = __builtin_ctz(static_cast<unsigned>(mask));
        const int pm = mask ^ (1 << bit);
        sums[mask] = {sums[pm].x + v[bit].x, sums[pm].y + v[bit].y};
    }

    std::unordered_set<int64_t> even_keys;
    std::unordered_set<int64_t> odd_keys;
    even_keys.reserve(full + 16);
    odd_keys.reserve(full + 16);
    vector<Point> even_points;
    even_points.reserve(full / 2 + 1);

    for (int mask = 0; mask < full; ++mask) {
        const int64_t key = pair_key(sums[mask].x, sums[mask].y);
        if ((__builtin_popcount(static_cast<unsigned>(mask)) & 1) == 0) {
            if (!even_keys.insert(key).second) {
                return false;
            }
            even_points.push_back(sums[mask]);
        } else {
            if (!odd_keys.insert(key).second) {
                return false;
            }
        }
    }

    for (const Point& c : even_points) {
        for (const Point& dxy : lattice) {
            const int nx = c.x + dxy.x;
            const int ny = c.y + dxy.y;
            if (odd_keys.find(pair_key(nx, ny)) != odd_keys.end()) {
                continue;
            }
            int count = 0;
            for (const Point& p : lattice) {
                if (even_keys.find(pair_key(nx + p.x, ny + p.y)) != even_keys.end()) {
                    ++count;
                    if (count >= 2) {
                        return false;
                    }
                }
            }
        }
    }

    return true;
}

static bool find_for_m_and_dimension(int m, int d, int thread_count) {
    const vector<Point> lattice = circle_points(m);
    if (lattice.empty()) {
        return false;
    }

    const auto pairs = opposite_pairs(lattice);
    if (static_cast<int>(pairs.size()) < d) {
        return false;
    }

    vector<Point> reps;
    reps.reserve(pairs.size());
    for (const auto& pr : pairs) {
        const Point& a = pr[0];
        const Point& b = pr[1];
        const bool take_a = (a.x > b.x) || (a.x == b.x && a.y > b.y);
        reps.push_back(take_a ? a : b);
    }

    const auto bad_disp = build_bad_displacements(lattice);

    const int u = static_cast<int>(reps.size());
    const int first_max = u - d;
    if (first_max < 0) {
        return false;
    }

    if (thread_count <= 0) {
        thread_count = 1;
    }
    thread_count = std::min(thread_count, first_max + 1);

    std::atomic<bool> found{false};

    auto worker = [&](int tid) {
        vector<Point> selected(d);

        std::function<void(int, int)> dfs = [&](int pos, int next_idx) {
            if (found.load(std::memory_order_relaxed)) {
                return;
            }
            if (pos == d) {
                if (test_selected_vectors(selected, lattice)) {
                    found.store(true, std::memory_order_relaxed);
                }
                return;
            }

            const int remaining = d - pos;
            const int limit = u - remaining;
            for (int i = next_idx; i <= limit; ++i) {
                if (found.load(std::memory_order_relaxed)) {
                    return;
                }
                selected[pos] = reps[i];
                if (four_tuple_has_forbidden_sum(selected, pos + 1, bad_disp)) {
                    continue;
                }
                dfs(pos + 1, i + 1);
            }
        };

        for (int first = tid; first <= first_max; first += thread_count) {
            if (found.load(std::memory_order_relaxed)) {
                return;
            }
            selected[0] = reps[first];
            dfs(1, first + 1);
        }
    };

    vector<std::thread> threads;
    threads.reserve(thread_count);
    for (int t = 0; t < thread_count; ++t) {
        threads.emplace_back(worker, t);
    }
    for (auto& th : threads) {
        th.join();
    }

    return found.load(std::memory_order_relaxed);
}

uint64_t solve_r_sq(uint64_t n, int thread_count = 1) {
    if (n <= 2) {
        return 1;
    }

    int d = 1;
    uint64_t capacity = 1;
    while (capacity < n) {
        ++d;
        if (d >= 64) {
            capacity = std::numeric_limits<uint64_t>::max();
            break;
        }
        capacity <<= 1;
    }

    const int dimension = d;

    int m = 1;
    while (true) {
        if (opposite_pair_count_by_factorization(static_cast<uint64_t>(m)) >= dimension &&
            find_for_m_and_dimension(m, dimension, thread_count)) {
            return static_cast<uint64_t>(m);
        }
        ++m;
    }
}

bool run_checkpoints(int thread_count) {
    const uint64_t r2 = solve_r_sq(2, thread_count);
    if (r2 != 1) {
        cerr << "Checkpoint failed: expected R(2)^2 = 1, got " << r2 << "\n";
        return false;
    }

    const uint64_t r4 = solve_r_sq(4, thread_count);
    if (r4 != 5) {
        cerr << "Checkpoint failed: expected R(4)^2 = 5, got " << r4 << "\n";
        return false;
    }

    const uint64_t r4_again = solve_r_sq(4, thread_count);
    if (r4_again != r4) {
        cerr << "Checkpoint failed: recomputation mismatch for n=4 ("
             << r4 << " vs " << r4_again << ")\n";
        return false;
    }

    uint64_t prev = 0;
    for (uint64_t ncheck = 2; ncheck <= 8; ++ncheck) {
        const uint64_t cur = solve_r_sq(ncheck, thread_count);
        if (ncheck > 2 && cur < prev) {
            cerr << "Checkpoint failed: monotonicity violated at n=" << ncheck << "\n";
            return false;
        }
        prev = cur;
    }

    return true;
}

int main(int argc, char** argv) {
    uint64_t n = 500;
    bool skip_checkpoints = false;

    int thread_count = static_cast<int>(std::thread::hardware_concurrency());
    if (thread_count <= 0) {
        thread_count = 1;
    }

    for (int i = 1; i < argc; ++i) {
        const string arg(argv[i]);
        if (arg.rfind("--n=", 0) == 0) {
            n = std::stoull(arg.substr(4));
        } else if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else if (arg.rfind("--threads=", 0) == 0) {
            thread_count = std::stoi(arg.substr(10));
            if (thread_count <= 0) {
                thread_count = 1;
            }
        } else {
            cerr << "Unknown argument: " << arg << "\n";
            return 1;
        }
    }

    if (!skip_checkpoints) {
        if (!run_checkpoints(thread_count)) {
            return 1;
        }
    }

    cout << solve_r_sq(n, thread_count) << "\n";
    return 0;
}

Python

import math
import sys

def circle_points(m):
    r = int(math.isqrt(m))
    pts = []
    for x in range(-r, r + 1):
        y2 = m - x * x
        if y2 < 0:
            continue
        y = math.isqrt(y2)
        if y * y == y2:
            pts.append((x, y))
            if y != 0:
                pts.append((x, -y))
    
    pts.sort()
    unique = []
    for p in pts:
        if not unique or unique[-1] != p:
            unique.append(p)
    return unique

def opposite_pairs(points):
    used = set()
    pairs = []
    for x, y in points:
        k = (x, y)
        if k in used:
            continue
        ov = (-x, -y)
        used.add(k)
        used.add(ov)
        pairs.append(((x, y), ov))
    return pairs

def opposite_pair_count_by_factorization(m):
    if m == 0:
        return 0
    n = m
    pairs = 2
    
    while n & 1 == 0:
        n >>= 1
        
    p = 3
    while p * p <= n:
        if n % p != 0:
            p += 2
            continue
        exp = 0
        while n % p == 0:
            n //= p
            exp += 1
        if (p & 3) == 3:
            if exp & 1:
                return 0
        elif (p & 3) == 1:
            pairs *= (exp + 1)
        p += 2
        
    if n > 1:
        if (n & 3) == 3:
            return 0
        if (n & 3) == 1:
            pairs *= 2
            
    return pairs

def build_bad_displacements(points):
    bad = set()
    for ax, ay in points:
        for bx, by in points:
            if ax == bx and ay == by:
                continue
            bad.add((ax - bx, ay - by))
    return bad

def four_tuple_has_forbidden_sum(selected, used_count, bad_disp):
    if used_count < 4:
        return False
        
    dx, dy = selected[used_count - 1]
    for i in range(used_count - 1):
        ax, ay = selected[i]
        for j in range(i + 1, used_count - 1):
            bx, by = selected[j]
            for k in range(j + 1, used_count - 1):
                cx, cy = selected[k]
                for mask in range(16):
                    sx = (-ax if (mask & 1) else ax) + \
                         (-bx if (mask & 2) else bx) + \
                         (-cx if (mask & 4) else cx) + \
                         (-dx if (mask & 8) else dx)
                    sy = (-ay if (mask & 1) else ay) + \
                         (-by if (mask & 2) else by) + \
                         (-cy if (mask & 4) else cy) + \
                         (-dy if (mask & 8) else dy)
                    if sx == 0 and sy == 0:
                        continue
                    if (sx, sy) in bad_disp:
                        return True
                        
    return False

def test_selected_vectors(v, lattice):
    d = len(v)
    full = 1 << d
    sums = [(0, 0)] * full
    
    for mask in range(1, full):
        bit = (mask & -mask).bit_length() - 1
        pm = mask ^ (1 << bit)
        sums[mask] = (sums[pm][0] + v[bit][0], sums[pm][1] + v[bit][1])
        
    even_keys = set()
    odd_keys = set()
    even_points = []
    
    for mask in range(full):
        key = sums[mask]
        if mask.bit_count() & 1 == 0:
            if key in even_keys:
                return False
            even_keys.add(key)
            even_points.append(sums[mask])
        else:
            if key in odd_keys:
                return False
            odd_keys.add(key)
            
    for cx, cy in even_points:
        for dx, dy in lattice:
            nx = cx + dx
            ny = cy + dy
            if (nx, ny) in odd_keys:
                continue
            count = 0
            for px, py in lattice:
                if (nx + px, ny + py) in even_keys:
                    count += 1
                    if count >= 2:
                        return False
                        
    return True

def dfs(next_idx, pos, d, reps, bad_disp, selected, lattice):
    if pos == d:
        return test_selected_vectors(selected, lattice)
        
    limit = len(reps) - (d - pos)
    for i in range(next_idx, limit + 1):
        selected[pos] = reps[i]
        if four_tuple_has_forbidden_sum(selected, pos + 1, bad_disp):
            continue
        if dfs(i + 1, pos + 1, d, reps, bad_disp, selected, lattice):
            return True
            
    return False

def find_for_m_and_dimension(m, d):
    lattice = circle_points(m)
    if not lattice: return False
    
    pairs = opposite_pairs(lattice)
    if len(pairs) < d: return False
    
    reps = []
    for a, b in pairs:
        take_a = (a[0] > b[0]) or (a[0] == b[0] and a[1] > b[1])
        reps.append(a if take_a else b)
        
    bad_disp = build_bad_displacements(lattice)
    u = len(reps)
    first_max = u - d
    if first_max < 0: return False
    
    selected = [(0, 0)] * d
    for first in range(first_max + 1):
        selected[0] = reps[first]
        if dfs(first + 1, 1, d, reps, bad_disp, selected, lattice):
            return True
            
    return False

def solve_r_sq(n):
    if n <= 2: return 1
    d = 1
    capacity = 1
    while capacity < n:
        d += 1
        if d >= 64: break
        capacity <<= 1
        
    m = 1
    while True:
        if opposite_pair_count_by_factorization(m) >= d:
            if find_for_m_and_dimension(m, d):
                return m
        m += 1

def solve():
    return str(solve_r_sq(500))

def run_checkpoints():
    assert solve_r_sq(2) == 1
    assert solve_r_sq(4) == 5

if __name__ == "__main__":
    run_checkpoints()
    print(solve())

Java

import java.util.ArrayList;
import java.util.Collections;
import java.util.HashSet;
import java.util.List;

public class Euler983 {
    static class Point implements Comparable<Point> {
        int x, y;

        Point(int x, int y) {
            this.x = x;
            this.y = y;
        }

        @Override
        public int compareTo(Point o) {
            if (this.x != o.x)
                return Integer.compare(this.x, o.x);
            return Integer.compare(this.y, o.y);
        }

        @Override
        public boolean equals(Object obj) {
            if (!(obj instanceof Point))
                return false;
            Point p = (Point) obj;
            return x == p.x && y == p.y;
        }

        @Override
        public int hashCode() {
            return x * 31 + y;
        }
    }

    static long pairKey(int x, int y) {
        return (((long) x) << 32) ^ ((long) y & 0xFFFFFFFFL);
    }

    static List<Point> circlePoints(int m) {
        int r = (int) Math.sqrt(m);
        List<Point> pts = new ArrayList<>();
        for (int x = -r; x <= r; ++x) {
            int y2 = m - x * x;
            if (y2 < 0)
                continue;
            int y = (int) Math.sqrt(y2);
            if (y * y == y2) {
                pts.add(new Point(x, y));
                if (y != 0)
                    pts.add(new Point(x, -y));
            }
        }
        Collections.sort(pts);
        List<Point> unique = new ArrayList<>();
        for (Point p : pts) {
            if (unique.isEmpty() || !unique.get(unique.size() - 1).equals(p)) {
                unique.add(p);
            }
        }
        return unique;
    }

    static List<Point[]> oppositePairs(List<Point> points) {
        HashSet<Long> used = new HashSet<>();
        List<Point[]> pairs = new ArrayList<>();
        for (Point v : points) {
            long k = pairKey(v.x, v.y);
            if (used.contains(k))
                continue;
            Point ov = new Point(-v.x, -v.y);
            used.add(k);
            used.add(pairKey(ov.x, ov.y));
            pairs.add(new Point[] { v, ov });
        }
        return pairs;
    }

    static int oppositePairCountByFactorization(long m) {
        if (m == 0)
            return 0;
        long n = m;
        int pairs = 2;

        while ((n & 1) == 0)
            n >>= 1;

        for (long p = 3; p * p <= n; p += 2) {
            if (n % p != 0)
                continue;
            int exp = 0;
            while (n % p == 0) {
                n /= p;
                exp++;
            }
            if ((p & 3) == 3) {
                if ((exp & 1) != 0)
                    return 0;
            } else if ((p & 3) == 1) {
                pairs *= (exp + 1);
            }
        }

        if (n > 1) {
            if ((n & 3) == 3)
                return 0;
            if ((n & 3) == 1)
                pairs *= 2;
        }

        return pairs;
    }

    static HashSet<Long> buildBadDisplacements(List<Point> points) {
        HashSet<Long> bad = new HashSet<>();
        for (Point a : points) {
            for (Point b : points) {
                if (a.equals(b))
                    continue;
                bad.add(pairKey(a.x - b.x, a.y - b.y));
            }
        }
        return bad;
    }

    static boolean fourTupleHasForbiddenSum(Point[] selected, int usedCount, HashSet<Long> badDisp) {
        if (usedCount < 4)
            return false;
        Point d = selected[usedCount - 1];
        for (int i = 0; i < usedCount - 1; ++i) {
            Point a = selected[i];
            for (int j = i + 1; j < usedCount - 1; ++j) {
                Point b = selected[j];
                for (int k = j + 1; k < usedCount - 1; ++k) {
                    Point c = selected[k];
                    for (int mask = 0; mask < 16; ++mask) {
                        int sx = ((mask & 1) != 0 ? -a.x : a.x) +
                                ((mask & 2) != 0 ? -b.x : b.x) +
                                ((mask & 4) != 0 ? -c.x : c.x) +
                                ((mask & 8) != 0 ? -d.x : d.x);
                        int sy = ((mask & 1) != 0 ? -a.y : a.y) +
                                ((mask & 2) != 0 ? -b.y : b.y) +
                                ((mask & 4) != 0 ? -c.y : c.y) +
                                ((mask & 8) != 0 ? -d.y : d.y);
                        if (sx == 0 && sy == 0)
                            continue;
                        if (badDisp.contains(pairKey(sx, sy)))
                            return true;
                    }
                }
            }
        }
        return false;
    }

    static boolean testSelectedVectors(Point[] v, List<Point> lattice) {
        int d = v.length;
        int full = 1 << d;
        Point[] sums = new Point[full];
        sums[0] = new Point(0, 0);

        for (int mask = 1; mask < full; ++mask) {
            int bit = Integer.numberOfTrailingZeros(mask);
            int pm = mask ^ (1 << bit);
            sums[mask] = new Point(sums[pm].x + v[bit].x, sums[pm].y + v[bit].y);
        }

        HashSet<Long> evenKeys = new HashSet<>();
        HashSet<Long> oddKeys = new HashSet<>();
        List<Point> evenPoints = new ArrayList<>();

        for (int mask = 0; mask < full; ++mask) {
            long key = pairKey(sums[mask].x, sums[mask].y);
            if ((Integer.bitCount(mask) & 1) == 0) {
                if (!evenKeys.add(key))
                    return false;
                evenPoints.add(sums[mask]);
            } else {
                if (!oddKeys.add(key))
                    return false;
            }
        }

        for (Point c : evenPoints) {
            for (Point dxy : lattice) {
                int nx = c.x + dxy.x;
                int ny = c.y + dxy.y;
                if (oddKeys.contains(pairKey(nx, ny)))
                    continue;
                int count = 0;
                for (Point p : lattice) {
                    if (evenKeys.contains(pairKey(nx + p.x, ny + p.y))) {
                        count++;
                        if (count >= 2)
                            return false;
                    }
                }
            }
        }
        return true;
    }

    static boolean dfs(int nextIdx, int pos, int d, List<Point> reps, HashSet<Long> badDisp, Point[] selected,
            List<Point> lattice) {
        if (pos == d) {
            return testSelectedVectors(selected, lattice);
        }
        int remaining = d - pos;
        int limit = reps.size() - remaining;
        for (int i = nextIdx; i <= limit; ++i) {
            selected[pos] = reps.get(i);
            if (fourTupleHasForbiddenSum(selected, pos + 1, badDisp))
                continue;
            if (dfs(i + 1, pos + 1, d, reps, badDisp, selected, lattice))
                return true;
        }
        return false;
    }

    static boolean findForMAndDimension(int m, int d) {
        List<Point> lattice = circlePoints(m);
        if (lattice.isEmpty())
            return false;

        List<Point[]> pairs = oppositePairs(lattice);
        if (pairs.size() < d)
            return false;

        List<Point> reps = new ArrayList<>();
        for (Point[] pr : pairs) {
            Point a = pr[0];
            Point b = pr[1];
            boolean takeA = (a.x > b.x) || (a.x == b.x && a.y > b.y);
            reps.add(takeA ? a : b);
        }

        HashSet<Long> badDisp = buildBadDisplacements(lattice);
        int u = reps.size();
        int firstMax = u - d;
        if (firstMax < 0)
            return false;

        Point[] selected = new Point[d];

        for (int first = 0; first <= firstMax; ++first) {
            selected[0] = reps.get(first);
            if (dfs(first + 1, 1, d, reps, badDisp, selected, lattice))
                return true;
        }
        return false;
    }

    static long solveRSq(long n) {
        if (n <= 2)
            return 1;

        int d = 1;
        long capacity = 1;
        while (capacity < n) {
            d++;
            if (d >= 64) {
                capacity = Long.MAX_VALUE;
                break;
            }
            capacity <<= 1;
        }

        int m = 1;
        while (true) {
            if (oppositePairCountByFactorization(m) >= d) {
                if (findForMAndDimension(m, d)) {
                    return m;
                }
            }
            m++;
        }
    }

    public static String solve() {
        return Long.toString(solveRSq(500));
    }

    public static void main(String[] args) {
        if (solveRSq(2) != 1 || solveRSq(4) != 5) {
            System.err.println("Validation failed");
            return;
        }
        System.out.println(solve());
    }
}