Project Euler Problem 804

Solution

Answer: 4921370551019052

Let

$$Q(x,y)=x^2+xy+41y^2.$$

Since $g(n)$ counts the number of integer pairs $(x,y)$ with $Q(x,y)=n$, we have

$$T(N)=\sum_{n=1}^N g(n)$$

which is exactly the number of integer lattice points satisfying

$$x^2+xy+41y^2 \le N,$$

excluding $(0,0)$ (since $n\ge 1$).

Rewrite the quadratic form:

$$4Q=(2x+y)^2+163y^2.$$

Thus we count integer pairs $(x,y)$ such that

$$(2x+y)^2+163y^2\le 4N.$$

For a fixed $y$, let

$$r=\left\lfloor \sqrt{4N-163y^2}\right\rfloor.$$

Then $t=2x+y$ must satisfy $|t|\le r$ and have the same parity as $y$. The number of such $t$ is

$$r + \mathbf{1}_{r\equiv y \pmod 2}.$$

Summing over all valid $y$ gives an $O(\sqrt N)$ algorithm:

from math import isqrt

def T(N):
    M = 4 * N
    ymax = isqrt(M // 163)

    total = -1  # exclude (0,0)

    for y in range(ymax + 1):
        r = isqrt(M - 163 * y * y)
        count = r + ((r & 1) == (y & 1))

        total += count if y == 0 else 2 * count

    return total

print(T(10**3))   # 474
print(T(10**6))   # 492128
print(T(10**16))

This reproduces the given checks:

$T(10^3)=474$
$T(10^6)=492128$

and computes:

$$T(10^{16}) = 4921370551019052.$$

Answer: 4921370551019052