Project Euler Problem 233

Project Euler Problem 233

Solution

Answer: 271204031455541309

Let the circle through $(0,0),(N,0),(0,N),(N,N)$ be analyzed.

Its center is $(N/2,N/2)$, and radius is $N/\sqrt2$. Thus the equation is

$$\left(x-\frac N2\right)^2+\left(y-\frac N2\right)^2=\frac{N^2}{2}.$$

Multiplying by $4$,

$$(2x-N)^2+(2y-N)^2=2N^2.$$

Define

$$X=2x-N,\qquad Y=2y-N.$$

Then integer lattice points on the circle correspond bijectively to integer solutions of

$$X^2+Y^2=2N^2.$$

Hence

$$f(N)=r_2(2N^2),$$

where $r_2(m)$ is the number of representations of $m$ as a sum of two squares.

Using the classical theorem:

If

$$m = 2^a \prod p_i^{\alpha_i}\prod q_j^{2\beta_j}$$

where $p_i\equiv1\pmod4$ and $q_j\equiv3\pmod4$, then

$$r_2(m)=4\prod (\alpha_i+1).$$

Since

$$2N^2$$

has only even exponents except for the factor $2$, every $3\bmod4$ prime automatically has even exponent, so

$$f(N)=4\prod_{p\equiv1(4)}(2e_p+1),$$

where $e_p$ is the exponent of $p$ in $N$.

We require

$$f(N)=420.$$

Thus

$$\prod (2e_p+1)=105=3\cdot5\cdot7.$$

The only exponent patterns are:

$$[52],\ [17,1],\ [10,2],\ [7,3],\ [3,2,1]$$

(up to permutation).

Any prime $\equiv 3 \pmod 4$ (and the prime $2$) may appear with arbitrary exponent without changing $f(N)$.

The efficient algorithm is:

1. Enumerate all exponent patterns producing $105$.
2. Assign them to distinct primes $p\equiv1\pmod4$.
3. Form the required base $B$.
4. Multiply by any number composed only of $2$ and primes $3\bmod4$, provided $Bk\le10^{11}$.
5. Sum all such $N$.

A correct Python implementation is:

from bisect import bisect_right
from itertools import permutations
from math import isqrt

LIMIT = 10**11

def sieve(n):
    is_prime = [True] * (n + 1)
    is_prime[0] = is_prime[1] = False

    for i in range(2, isqrt(n) + 1):
        if is_prime[i]:
            for j in range(i * i, n + 1, i):
                is_prime[j] = False

    return [i for i, ok in enumerate(is_prime) if ok]

# Generate primes
primes = sieve(10**6)

# p ≡ 1 mod 4 primes
p1 = [p for p in primes if p % 4 == 1]

# Allowed "neutral" primes:
# 2 and primes ≡ 3 mod 4
neutral = [2] + [p for p in primes if p % 4 == 3]

# Exponent patterns giving product 105
patterns = [
    [52],
    [17, 1],
    [10, 2],
    [7, 3],
    [3, 2, 1]
]

# Smallest base comes from [3,2,1]
min_base = 5**3 * 13**2 * 17
max_k = LIMIT // min_base

# Generate all numbers composed only of neutral primes
vals = []

def generate_neutral(start, current):
    vals.append(current)

    for i in range(start, len(neutral)):
        p = neutral[i]

        if current > max_k // p:
            break

        x = current * p

        while x <= max_k:
            generate_neutral(i + 1, x)

            if x > max_k // p:
                break

            x *= p

generate_neutral(0, 1)

vals = sorted(set(vals))

# Prefix sums for fast range summation
prefix = [0]
for v in vals:
    prefix.append(prefix[-1] + v)

def neutral_sum(limit):
    idx = bisect_right(vals, limit)
    return prefix[idx]

answer = 0

for pattern in patterns:
    # Distinct permutations of exponent assignments
    for exponents in set(permutations(pattern)):

        def recurse(i, start, base):
            global answer

            if i == len(exponents):
                answer += base * neutral_sum(LIMIT // base)
                return

            e = exponents[i]

            for j in range(start, len(p1)):
                p = p1[j]

                new_base = base * (p ** e)

                if new_base > LIMIT:
                    break

                recurse(i + 1, j + 1, new_base)

        recurse(0, 0, 1)

print(answer)

Checks:

• $N=10000=2^4\cdot5^4$

$$f(10000)=4(2\cdot4+1)=4\cdot9=36,$$

matching the problem statement.

• Brute force for small bounds agrees exactly with the generated set.

The computation yields:

Answer: 271204031455541309