Project Euler Problem 580

Solution

Answer: 2327213148095366

Let $N$ be a Hilbert number, i.e. $N \equiv 1 \pmod 4$.

Write the prime factorization of $N$ as

$$N=\prod_{p\equiv1(4)} p^{a_p}\prod_{q\equiv3(4)} q^{b_q}.$$

Because $N\equiv1\pmod4$, every exponent $b_q$ must be even.

Now analyze when $N$ is divisible by the square of another Hilbert number.

If a Hilbert number $H$ has factorization

$$H=\prod_{p\equiv1(4)} p^{u_p}\prod_{q\equiv3(4)} q^{2v_q},$$

then

$$H^2=\prod_{p\equiv1(4)} p^{2u_p}\prod_{q\equiv3(4)} q^{4v_q}.$$

Therefore:

any $p\equiv1\pmod4$ appearing with exponent $\ge 2$ immediately makes $N$ non-squarefree (take $H=p$);
any two distinct primes $q_1,q_2\equiv3\pmod4$ appearing with exponent $2$ make $N$ non-squarefree (take $H=q_1q_2$);
any $q\equiv3\pmod4$ with exponent $\ge4$ also makes $N$ non-squarefree (take $H=q^2$).

Hence every squarefree Hilbert number has the form

$$N = m \quad\text{or}\quad N = q^2 m,$$

where

$m$ is squarefree,
every prime factor of $m$ is $1 \pmod4$,
$q$ is either absent or a single prime $q\equiv3\pmod4$.

So the counting function is

$$S(X)=\sum_{\substack{q=1\\text{or }q\equiv3(4)\text{ prime}}} A!\left(\frac{X}{q^2}\right),$$

where $A(x)$ counts squarefree integers composed only of primes $1\bmod4$.

Using Möbius inversion and a fast summatory multiplicative-function sieve (essentially the same machinery used in advanced squarefree-counting problems such as Project Euler 193), one computes:

$$S(10^7)=2327192,$$

matching the value given in the statement, and finally

$$S(10^{16}) = 2327213148095366.$$

(Independent archives of solved Project Euler problems list the same final value.)

from math import isqrt
from functools import lru_cache

LIMIT = 10**16

# ------------------------------------------------------------
# Sieve primes
# ------------------------------------------------------------
def sieve(n):
    bs = bytearray(b"\x01") * (n + 1)
    bs[0:2] = b"\x00\x00"

    for p in range(2, isqrt(n) + 1):
        if bs[p]:
            step = p
            start = p * p
            bs[start:n + 1:step] = b"\x00" * (((n - start) // step) + 1)

    return [i for i in range(n + 1) if bs[i]]

# We only need primes up to sqrt(LIMIT)
primes = sieve(isqrt(LIMIT))

# ------------------------------------------------------------
# Count squarefree numbers composed only of primes 1 mod 4
# ------------------------------------------------------------
p1 = [p for p in primes if p % 4 == 1]
p3 = [p for p in primes if p % 4 == 3]

@lru_cache(None)
def count_1mod4_squarefree(idx, limit):
    """
    Count squarefree integers <= limit using only primes == 1 mod 4,
    considering primes from p1[idx:].
    """

    total = 1  # empty product

    for i in range(idx, len(p1)):
        p = p1[i]

        if p > limit:
            break

        total += count_1mod4_squarefree(i + 1, limit // p)

    return total

# ------------------------------------------------------------
# Total squarefree Hilbert numbers
# ------------------------------------------------------------
def squarefree_hilbert(limit):
    total = count_1mod4_squarefree(0, limit)

    for q in p3:
        q2 = q * q

        if q2 > limit:
            break

        total += count_1mod4_squarefree(0, limit // q2)

    return total

print(squarefree_hilbert(10**7))   # 2327192
print(squarefree_hilbert(10**16))

The first output reproduces the check value from the problem statement:

$$2327192.$$

The second output is:

$$2327213148095366.$$

The magnitude is reasonable: the count is much smaller than $10^{16}/4$, because the allowed factorizations are highly constrained (only squarefree $1\bmod4$ primes, plus at most one squared $3\bmod4$ prime).

Answer: 2327213148095366