Project Euler

Lattice Points on a Circle

Let f(n) be the number of lattice points on the circle x^2 + y^2 = n^2. Find the sum of all positive integers n <= 10^11 such that f(n) = 420.

Source sync Apr 19, 2026

Problem #0233

Level Level 09

Solved By 2,654

Languages C++, Python

Answer 271204031455541309

Length 354 words

modular_arithmeticnumber_theorylinear_algebra

Let $f(N)$ be the number of points with integer coordinates that are on a circle passing through $(0,0)$, $(N,0)$,$(0,N)$, and $(N,N)$.

It can be shown that $f(10000) = 36$.

What is the sum of all positive integers $N \le 10^{11}$ such that $f(N) = 420$?

Problem 233: Lattice Points on a Circle

Mathematical Foundation

Theorem (Sum of Two Squares Representation Count). For a positive integer $m$ , the number of representations of $m$ as a sum of two squares (counting order and signs) is

r_2(m) = 4 \sum_{d \mid m} \chi(d)

where $\chi$ is the non-principal Dirichlet character modulo 4: $\chi(d) = 0$ if $2 \mid d$ , $\chi(d) = 1$ if $d \equiv 1 \pmod{4}$ , $\chi(d) = -1$ if $d \equiv 3 \pmod{4}$ .

Proof. This is a classical result following from the factorization of ideals in the Gaussian integers $\mathbb{Z}[i]$ . The norm $N(a + bi) = a^2 + b^2$ is multiplicative, and $r_2(m) = 4 \sum_{d \mid m} \chi(d)$ follows from the identity $\sum_{d \mid m} \chi(d) = \prod_{p \mid m} \sum_{j=0}^{v_p(m)} \chi(p)^j$ and the evaluation of each local factor. Specifically, for $p = 2$ : the factor is 1. For $p \equiv 1 \pmod{4}$ : the factor is $v_p(m) + 1$ . For $p \equiv 3 \pmod{4}$ : the factor is 1 if $v_p(m)$ is even, and 0 if odd. $\square$

Lemma. For $m = n^2$ with $n = 2^{a_0} \prod_{i} p_i^{a_i} \prod_{j} q_j^{b_j}$ where $p_i \equiv 1 \pmod{4}$ and $q_j \equiv 3 \pmod{4}$ :

f(n) = r_2(n^2) = 4 \prod_{i} (2a_i + 1).

Proof. In $n^2$ , each $p_i$ appears with exponent $2a_i$ and each $q_j$ with exponent $2b_j$ (even). By the theorem, the $q_j$ -factors each contribute 1 (since $2b_j$ is even), the factor of 2 contributes 1, and each $p_i$ -factor contributes $2a_i + 1$ . $\square$

Corollary. The condition $f(n) = 420$ is equivalent to $\prod_i (2a_i + 1) = 105$ , where the product is over primes $p_i \equiv 1 \pmod{4}$ dividing $n$ .

Proof. Immediate from $f(n) = 4 \prod_i (2a_i + 1) = 420$ . $\square$

The factorizations of 105 into factors $\geq 3$ (since each $a_i \geq 1$ forces $2a_i + 1 \geq 3$ ) are:

Factorization	Exponent pattern $(a_i)$
$105$	$(52)$
$3 \times 35$	$(1, 17)$
$5 \times 21$	$(2, 10)$
$7 \times 15$	$(3, 7)$
$3 \times 5 \times 7$	$(1, 2, 3)$

For each pattern, $n$ has the form $n = P \cdot k$ where $P = \prod p_i^{a_i}$ uses primes $\equiv 1 \pmod{4}$ and $k$ has no prime factor $\equiv 1 \pmod{4}$ .

Editorial

r_2(n^2) = 4 * product(2a_i + 1) for primes p_i ≡ 1 (mod 4) in factorization of n. We need product(2a_i + 1) = 105 = 357. n = k * P where P = product of p_i^{a_i} with p_i ≡ 1 mod 4, and k has no prime factor ≡ 1 mod 4. We primes ≡ 1 (mod 4) up to N: 5, 13, 17, 29, 37, 41, . We then iterate over each such P. Finally, count sum of all k <= L with no prime factor ≡ 1 (mod 4).

Pseudocode

N = 10^11
Primes ≡ 1 (mod 4) up to N: 5, 13, 17, 29, 37, 41, 
for each such P
Count sum of all k <= L with no prime factor ≡ 1 (mod 4)
Use inclusion-exclusion over primes ≡ 1 (mod 4) up to L

Complexity Analysis

Time: Dominated by the enumeration of products $P$ and the inclusion-exclusion for each. The number of valid $P$ values is moderate (tens of thousands). For each, the inclusion-exclusion over primes $\equiv 1 \pmod{4}$ up to $L$ takes $O(\sqrt{L})$ with careful sieving.
Space: $O(\sqrt{N})$ for sieve arrays.

Answer

\boxed{271204031455541309}

C++ project_euler/problem_233/solution.cpp

#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
typedef unsigned long long ull;
typedef __int128 lll;

const ll LIMIT = 100000000000LL; // 10^11

// We need:
// r_2(n^2) = 420 => product(2*a_i+1) = 105 = 3*5*7
// Exponent patterns (sorted descending for canonical form):
// {52}, {17,1}, {10,2}, {7,3}, {3,2,1}
//
// For {52}: p^52 with p=5 gives 5^52 >> 10^11. No solutions.
// For {17,1}: smallest is 5^17 * 13 = 762939453125*13 >> 10^11.
//   Or 5^1 * 13^17 >> 10^11. No solutions.
// For {10,2}: 5^10 * 13^2 = 9765625*169 = 1,650,390,625. OK.
//   5^2 * 13^10 = 25 * 137,858,491,849 >> 10^11. Only (big exp on small prime) works.
// For {7,3}: 5^7 * 13^3 = 78125*2197 = 171,656,250. OK.
//   5^3 * 13^7 = 125*62748517 = 7,843,564,625. OK.
// For {3,2,1}: many possibilities.
//
// Strategy: for each exponent assignment to primes in increasing order,
// enumerate using DFS. The prime with exp=1 can be huge, so we need a way
// to generate primes ≡ 1 mod 4 up to ~5*10^9.
//
// For the "free multiplier" k (no prime factor ≡ 1 mod 4), we need sum_no1mod4(L).
// When L is small, sieve directly. When L is large, use inclusion-exclusion
// with primes ≡ 1 mod 4 up to L.
//
// Alternative approach: generate primes ≡ 1 mod 4 using segmented sieve,
// or just sieve up to sqrt(LIMIT) ~ 316228 and for the last (largest) prime
// in the product, iterate directly.

// Sieve primes up to some limit
vector<int> small_primes;
vector<int> primes_1mod4_small; // primes ≡ 1 mod 4 up to sieve limit

void sieve(int n) {
    vector<bool> is_p(n + 1, true);
    is_p[0] = is_p[1] = false;
    for (int i = 2; (ll)i * i <= n; i++)
        if (is_p[i])
            for (int j = i * i; j <= n; j += i)
                is_p[j] = false;
    for (int i = 2; i <= n; i++)
        if (is_p[i]) {
            small_primes.push_back(i);
            if (i % 4 == 1) primes_1mod4_small.push_back(i);
        }
}

// Check if a number > sieve_limit is prime using trial division
bool is_prime_check(ll n) {
    if (n < 2) return false;
    if (n < (ll)small_primes.back() * small_primes.back()) {
        for (int p : small_primes) {
            if ((ll)p * p > n) return true;
            if (n % p == 0) return false;
        }
        return true;
    }
    for (int p : small_primes) {
        if ((ll)p * p > n) return true;
        if (n % p == 0) return false;
    }
    return true; // probably prime if we've checked up to sqrt
}

// Sum of all k in [1, L] that have no prime factor ≡ 1 mod 4
// For small L, sieve directly.
// primes_1mod4_for_sieve are the primes ≡ 1 mod 4 up to L.
ll sum_no_1mod4(ll L) {
    if (L <= 0) return 0;
    if (L == 1) return 1;

    // Sieve: mark numbers divisible by any prime ≡ 1 mod 4
    vector<bool> bad(L + 1, false);
    for (int p : primes_1mod4_small) {
        if (p > L) break;
        for (ll j = p; j <= L; j += p)
            bad[j] = true;
    }
    ll s = 0;
    for (ll i = 1; i <= L; i++)
        if (!bad[i]) s += i;
    return s;
}

ll total_sum = 0;

// Enumerate products of primes ≡ 1 mod 4 with given exponent sequence.
// exps[pos..] are exponents to assign to primes in increasing order.
// Primes must be > prev_prime.
// current_P is the product so far.
// For the last exponent, if it's 1, we can iterate over all primes ≡ 1 mod 4
// greater than prev_prime up to LIMIT/current_P using a segmented approach.

void enumerate(const vector<int>& exps, int pos, ll prev_prime, ll current_P) {
    if (pos == (int)exps.size()) {
        ll L = LIMIT / current_P;
        ll s = sum_no_1mod4(L);
        total_sum += s * current_P;
        return;
    }

    int e = exps[pos];
    bool is_last = (pos == (int)exps.size() - 1);

    if (is_last && e == 1) {
        // Last exponent is 1: the prime p contributes p^1 to the product.
        // We need current_P * p <= LIMIT, so p <= LIMIT / current_P.
        // p must be prime, ≡ 1 mod 4, and > prev_prime.
        // For each such p, the multiplier k has L = LIMIT / (current_P * p).
        // We need to sum over all such p: sum_no_1mod4(L) * current_P * p.
        //
        // For primes within our sieve range, iterate directly.
        // For primes beyond sieve range, we'd need segmented sieve.
        // Let's use segmented sieve approach.

        ll max_p = LIMIT / current_P;

        // Use primes from our small list first
        for (int pp : primes_1mod4_small) {
            ll p = pp;
            if (p <= prev_prime) continue;
            if (p > max_p) break;
            ll new_P = current_P * p;
            ll L = LIMIT / new_P;
            ll s = sum_no_1mod4(L);
            total_sum += s * new_P;
        }

        // Now handle primes beyond our sieve limit
        ll sieve_lim = primes_1mod4_small.empty() ? 2 : primes_1mod4_small.back();
        if (max_p > sieve_lim) {
            ll seg_start = max(sieve_lim + 1, prev_prime + 1);
            // Round up to make seg_start odd if needed
            // Use segmented sieve
            const ll SEG_SIZE = 1000000;
            for (ll lo = seg_start; lo <= max_p; lo += SEG_SIZE) {
                ll hi = min(lo + SEG_SIZE - 1, max_p);
                ll sz = hi - lo + 1;
                vector<bool> is_p(sz, true);
                for (int sp : small_primes) {
                    if ((ll)sp * sp > hi) break;
                    ll start = ((lo + sp - 1) / sp) * sp;
                    if (start == sp) start += sp; // don't mark the prime itself...
                    // Actually for segmented sieve, lo > sieve_limit > sp, so start > sp always
                    for (ll j = start; j <= hi; j += sp)
                        is_p[j - lo] = false;
                }
                for (ll i = 0; i < sz; i++) {
                    if (is_p[i]) {
                        ll p = lo + i;
                        if (p <= prev_prime) continue;
                        if (p % 4 != 1) continue;
                        ll new_P = current_P * p;
                        ll L = LIMIT / new_P;
                        ll s = sum_no_1mod4(L);
                        total_sum += s * new_P;
                    }
                }
            }
        }
        return;
    }

    // General case: pick a prime ≡ 1 mod 4, raise to power e
    for (int pp : primes_1mod4_small) {
        ll p = pp;
        if (p <= prev_prime) continue;
        ll pe = 1;
        bool overflow = false;
        for (int j = 0; j < e; j++) {
            if (pe > LIMIT / p) { overflow = true; break; }
            pe *= p;
        }
        if (overflow) break;
        if (current_P > LIMIT / pe) break;
        enumerate(exps, pos + 1, p, current_P * pe);
    }

    // Also check primes beyond sieve for small exponents (e >= 2)
    // p^e <= LIMIT / current_P => p <= (LIMIT/current_P)^(1/e)
    // For e >= 2 and LIMIT = 10^11, p <= 10^(11/2) = ~316227 which is within sieve range.
    // For e = 1 handled above. For e >= 2, sieve up to 316228 is sufficient.
    // Actually for e=2: p <= sqrt(10^11 / current_P). current_P >= 1, so p <= 316227. Covered.
    // For e=3: p <= (10^11)^(1/3) ~ 46415. Covered.
    // So we're fine.
}

int main() {
    // Sieve up to enough to cover:
    // - primes for enumeration: for e=2, up to 316228
    // - primes for sum_no_1mod4: L can be up to ~LIMIT/P_min
    //   For pattern (3,2,1) with smallest P: 5^3*13^2*17 = 358525, L ~ 278854
    //   For (7,3): 5^7*13^3 = 171656250, L ~ 583
    //   For (10,2): 5^10*13^2 = 1650390625, L ~ 60
    //   For last-prime-with-exp-1 cases: L = LIMIT/(current_P * p),
    //     when p is large, L is small.
    //   Maximum L occurs for the (3,2,1) pattern with the last prime = 17 (smallest):
    //     L = 10^11 / (5^3*13^2*17) ≈ 278854. Need primes ≡ 1 mod 4 up to 278854.
    //   But wait, for the last-prime-with-exp-1 optimization, L can be up to
    //     LIMIT / (current_P * (prev_prime+1)). E.g., current_P = 5^3 * 13^2 = 21125,
    //     prev_prime = 13, smallest valid p = 17: L = 10^11 / (21125*17) ≈ 278,416.
    //   So we need sieve up to ~280,000 for sum_no_1mod4 sieving purposes,
    //   PLUS we need segmented sieve for finding large primes.
    //   For trial division in is_prime_check, we need primes up to sqrt(5*10^9) ~ 70711.

    sieve(320000); // covers all needs

    // Exponent patterns: assignments of exponents to primes in increasing order
    // For unordered multiset of exponents, all orderings give different prime assignments.
    // Patterns from 105 = 3*5*7:
    // Unordered exponent sets and their permutations:
    // {52} -> only (52)
    // {17,1} -> (17,1), (1,17)
    // {10,2} -> (10,2), (2,10)
    // {7,3} -> (7,3), (3,7)
    // {3,2,1} -> all 6 permutations

    vector<vector<int>> patterns = {
        {52},
        {17, 1}, {1, 17},
        {10, 2}, {2, 10},
        {7, 3}, {3, 7},
        {3, 2, 1}, {3, 1, 2}, {2, 3, 1}, {2, 1, 3}, {1, 3, 2}, {1, 2, 3}
    };

    total_sum = 0;
    for (auto& exps : patterns) {
        enumerate(exps, 0, 0, 1);
    }

    cout << total_sum << endl;
    return 0;
}

Python project_euler/problem_233/solution.py

"""
Problem 233: Lattice Points on a Circle

Find sum of all n <= 10^11 such that the circle x^2+y^2=n^2 has exactly 420 lattice points.

r_2(n^2) = 4 * product(2*a_i + 1) for primes p_i ≡ 1 (mod 4) in factorization of n.
We need product(2*a_i + 1) = 105 = 3*5*7.

n = k * P where P = product of p_i^{a_i} with p_i ≡ 1 mod 4,
and k has no prime factor ≡ 1 mod 4.
"""

from itertools import permutations
import math

def sieve_primes(n):
    is_p = bytearray(b'\x01') * (n + 1)
    is_p[0] = is_p[1] = 0
    for i in range(2, int(n**0.5) + 1):
        if is_p[i]:
            is_p[i*i::i] = bytearray(len(is_p[i*i::i]))
    return [i for i in range(2, n + 1) if is_p[i]]

def solve():
    LIMIT = 10**11

    primes = sieve_primes(320000)
    primes_1mod4 = [p for p in primes if p % 4 == 1]

    def sum_no_1mod4(L):
        """Sum of k in [1, L] with no prime factor ≡ 1 mod 4."""
        if L <= 0:
            return 0
        bad = bytearray(L + 1)
        for p in primes_1mod4:
            if p > L:
                break
            for j in range(p, L + 1, p):
                bad[j] = 1
        return sum(i for i in range(1, L + 1) if not bad[i])

    def is_prime_check(n):
        if n < 2:
            return False
        for p in primes:
            if p * p > n:
                return True
            if n % p == 0:
                return False
        return True

    # Segmented sieve to find primes in [lo, hi]
    def seg_sieve(lo, hi):
        if lo > hi:
            return []
        sz = hi - lo + 1
        is_p = bytearray(b'\x01') * sz
        for p in primes:
            if p * p > hi:
                break
            start = ((lo + p - 1) // p) * p
            if start < lo:
                start += p
            for j in range(start, hi + 1, p):
                is_p[j - lo] = 0
        result = []
        for i in range(sz):
            if is_p[i] and (lo + i) > 1:
                result.append(lo + i)
        return result

    total = 0

    # Exponent patterns from factorizations of 105:
    # 105 -> (52), 3*35 -> (1,17), 5*21 -> (2,10), 7*15 -> (3,7), 3*5*7 -> (1,2,3)
    base_patterns = [
        (52,),
        (1, 17),
        (2, 10),
        (3, 7),
        (1, 2, 3),
    ]

    for base_pat in base_patterns:
        perms = set(permutations(base_pat))
        for perm in perms:
            def dfs(pos, prev_prime, current_P):
                nonlocal total
                if pos == len(perm):
                    L = LIMIT // current_P
                    s = sum_no_1mod4(L)
                    total += s * current_P
                    return

                e = perm[pos]
                is_last = (pos == len(perm) - 1)

                if is_last and e == 1:
                    # Last exponent is 1: iterate over primes ≡ 1 mod 4
                    max_p = LIMIT // current_P

                    # From small primes list
                    for p in primes_1mod4:
                        if p <= prev_prime:
                            continue
                        if p > max_p:
                            break
                        new_P = current_P * p
                        L = LIMIT // new_P
                        s = sum_no_1mod4(L)
                        total += s * new_P

                    # From segmented sieve for larger primes
                    sieve_lim = primes_1mod4[-1] if primes_1mod4 else 2
                    if max_p > sieve_lim:
                        seg_start = max(sieve_lim + 1, prev_prime + 1)
                        SEG_SIZE = 1000000
                        lo_seg = seg_start
                        while lo_seg <= max_p:
                            hi_seg = min(lo_seg + SEG_SIZE - 1, max_p)
                            seg_primes = seg_sieve(lo_seg, hi_seg)
                            for p in seg_primes:
                                if p <= prev_prime:
                                    continue
                                if p % 4 != 1:
                                    continue
                                new_P = current_P * p
                                L = LIMIT // new_P
                                s = sum_no_1mod4(L)
                                total += s * new_P
                            lo_seg = hi_seg + 1
                    return

                # General case
                for p in primes_1mod4:
                    if p <= prev_prime:
                        continue
                    pe = p ** e
                    if current_P * pe > LIMIT:
                        break
                    dfs(pos + 1, p, current_P * pe)

            dfs(0, 0, 1)

    print(total)

if __name__ == "__main__":
    solve()