Detailed mathematical approach

Problem Summary

A Hilbert number is a positive integer of the form \(4k+1\). The problem asks for the number of Hilbert numbers below a large exclusive limit \(X\) that are squarefree inside the Hilbert semigroup: equivalently, numbers \(n\equiv1\pmod4\) such that no Hilbert number \(h>1\) satisfies \(h^2\mid n\).

The implementations work with the inclusive bound \(N=X-1\), so the target quantity is

$S(N)=\#\{n\le N:n\equiv1\pmod4,\ \nexists h\equiv1\pmod4,\ h>1,\ h^2\mid n\}.$

A direct scan with explicit Hilbert-square testing is far too slow near \(10^{16}\), so the solution rewrites the condition as a Möbius-weighted sum over odd square divisors.

Mathematical Approach

The main job is to understand which ordinary prime exponents are allowed in a squarefree Hilbert number, and then encode that description with inclusion-exclusion.

Step 1: Translate the Hilbert-square condition into ordinary prime exponents

Write an odd number \(n\equiv1\pmod4\) as

$n=\prod_{p\equiv1\pmod4} p^{a_p}\prod_{q\equiv3\pmod4} q^{b_q}.$

If some prime \(p\equiv1\pmod4\) has \(a_p\ge2\), then \(p\) itself is a Hilbert number and \(p^2\mid n\), so \(n\) is not squarefree in the Hilbert sense.

For primes \(q\equiv3\pmod4\), one copy of \(q^2\) is still allowed because \(q\not\equiv1\pmod4\). But two square units are too many. If

$\sum_{q\equiv3\pmod4}\left\lfloor\frac{b_q}{2}\right\rfloor\ge2,$

then we can build a nontrivial Hilbert number from those \(3\bmod4\) primes: for example \(q^2\) from one prime with \(b_q\ge4\), or \(qr\) from two different primes with \(b_q,b_r\ge2\). Its square then divides \(n\).

Therefore \(n\) is squarefree Hilbert exactly when

$a_p\le1\quad\text{for every }p\equiv1\pmod4,$

$\sum_{q\equiv3\pmod4}\left\lfloor\frac{b_q}{2}\right\rfloor\le1.$

Step 2: Obtain a structural decomposition

The previous criterion has a concrete consequence. A squarefree Hilbert number is of one of two types:

$\text{either }n\text{ is ordinary squarefree,}$

$\text{or }n=q^2m\text{ for exactly one prime }q\equiv3\pmod4\text{ and an ordinary squarefree }m.$

The second case is unique. If two different \(3\bmod4\) primes contributed square factors, or if one such prime contributed exponent at least \(4\), then the sum of the floor terms above would be at least \(2\), which is forbidden.

This explains examples such as \(45=3^2\cdot5\), which is allowed, and \(81=3^4\), which is not.

Step 3: Encode ordinary squarefreeness with the Möbius function

Let \(\mathbf{1}_{\mathrm{sqf}}(m)\) be the indicator of ordinary squarefree integers. The standard identity is

$\mathbf{1}_{\mathrm{sqf}}(m)=\sum_{u^2\mid m}\mu(u).$

Using the decomposition from Step 2, the indicator of squarefree Hilbert numbers becomes

$I(n)=\mathbf{1}_{n\equiv1\pmod4}\left(\mathbf{1}_{\mathrm{sqf}}(n)+\sum_{\substack{q\equiv3\pmod4\\ q\text{ prime}}}\mathbf{1}_{q^2\mid n}\,\mathbf{1}_{\mathrm{sqf}}\!\left(\frac{n}{q^2}\right)\right).$

The first term counts numbers with no square factor at all. The second term restores the legal case where the only square part is one factor \(q^2\) with \(q\equiv3\pmod4\).

Step 4: Collect everything into one coefficient

Substitute the Möbius identity into both squarefree indicators:

$I(n)=\mathbf{1}_{n\equiv1\pmod4}\left(\sum_{u^2\mid n}\mu(u)+\sum_{\substack{q\equiv3\pmod4\\ q\text{ prime}\\ q^2\mid n}}\ \sum_{u^2\mid n/q^2}\mu(u)\right).$

In the second double sum, reindex with \(v=qu\). Then all terms become contributions from odd squares \(v^2\mid n\):

$I(n)=\mathbf{1}_{n\equiv1\pmod4}\sum_{\substack{v^2\mid n\\ v\text{ odd}}} c(v),$

where

$c(v)=\mu(v)+\sum_{\substack{q\mid v\\ q\equiv3\pmod4\\ q\text{ prime}}}\mu\!\left(\frac{v}{q}\right).$

This is the exact coefficient formula used by the implementations. The correction term cancels illegal square divisors such as \(5^2\), keeps legal ones such as \(3^2\), and subtracts higher powers such as \(3^4\) again.

Step 5: Turn the indicator into a summatory formula

Let

$A(t)=\#\{m\le t:m\equiv1\pmod4\}=\left\lfloor\frac{t+3}{4}\right\rfloor.$

For any odd \(v\), the condition \(v^2\mid n\) with \(n\equiv1\pmod4\) is equivalent to writing \(n=v^2m\) with \(m\equiv1\pmod4\), because \(v^2\equiv1\pmod4\). Hence

$S(N)=\sum_{\substack{v\ge1\\ v\text{ odd}}} c(v)\,A\!\left(\left\lfloor\frac{N}{v^2}\right\rfloor\right).$

Only values \(v\le\sqrt{N}\) contribute, so the sum is finite:

$\boxed{S(N)=\sum_{\substack{1\le v\le\lfloor\sqrt{N}\rfloor\\ v\text{ odd}}} c(v)\left\lfloor\frac{\left\lfloor N/v^2\right\rfloor+3}{4}\right\rfloor.}$

Worked Example: Why \(S(99)=23\)

There are exactly \(25\) integers \(n\le99\) with \(n\equiv1\pmod4\). For this bound, only odd \(v\le9\) matter.

The relevant coefficients are

$c(1)=1,\qquad c(3)=0,\qquad c(5)=-1,\qquad c(7)=0,\qquad c(9)=-1.$

Indeed, \(c(5)=-1\) removes multiples of \(25\), while \(c(9)=-1\) removes multiples of \(81=(3^2)^2\). The terms for \(3\) and \(7\) vanish because a single square \(3^2\) or \(7^2\) is still allowed in the Hilbert setting.

Therefore

$S(99)=A(99)-A\!\left(\left\lfloor\frac{99}{25}\right\rfloor\right)-A\!\left(\left\lfloor\frac{99}{81}\right\rfloor\right)=25-1-1=23,$

which matches the checkpoint count below \(100\).

How the Code Works

The C++, Python, and Java implementations all use the same arithmetic plan. They first convert the exclusive input bound \(X\) into the inclusive bound \(N=X-1\), compute \(U=\lfloor\sqrt{N}\rfloor\), and store only odd candidates up to \(U\). That halves the sieve size immediately, because every relevant square root \(v\) is odd.

Next, the implementation builds an odd-only prime sieve, derives the Möbius function \(\mu(v)\) on those odd values, and records the primes congruent to \(3\pmod4\). It then starts from the base coefficients \(\mu(v)\) and adds the correction term \(\mu(v/q)\) along odd multiples of each prime \(q\equiv3\pmod4\), producing the final weights \(c(v)\).

Finally, it evaluates

$c(v)\left\lfloor\frac{\left\lfloor N/v^2\right\rfloor+3}{4}\right\rfloor$

for every odd \(v\le U\) and accumulates the result. The C++ version keeps the same formula but can split the final summation into several ranges when multithreading is enabled; the Python and Java versions run the same arithmetic in a single thread.

Complexity Analysis

Let \(U=\lfloor\sqrt{N}\rfloor\). Building the odd-only sieve and the odd Möbius table costs \(O(U\log\log U)\) time and \(O(U)\) memory. Adding the correction terms over multiples of primes \(q\equiv3\pmod4\) has the same harmonic-sum behavior, so it also stays within \(O(U\log\log U)\) time.

The final accumulation is a single pass over the odd indices, so it is \(O(U)\). Overall, the method is near-linear in \(\sqrt{N}\) and uses \(O(U)\) memory.

Footnotes and References

1. Project Euler problem page: https://projecteuler.net/problem=580
2. Möbius function: Wikipedia — Möbius function
3. Squarefree integer: Wikipedia — Squarefree integer
4. Inclusion-exclusion principle: Wikipedia — Inclusion-exclusion principle
5. Modular arithmetic: Wikipedia — Modular arithmetic