Project Euler Problem 193

A positive integer n is called squarefree, if no square of a prime divides n, thus 1, 2, 3, 5, 6, 7, 10, 11 are squarefr

Project Euler Problem 193

Solution

Answer: 684465067343069

Let

$$N = 2^{50}.$$

We want to count the positive integers $n < N$ that are squarefree, meaning no prime square divides $n$.

A classical way to count squarefree numbers uses the Möbius function and inclusion–exclusion.


1. Mathematical analysis

A number is not squarefree if it is divisible by $p^2$ for some prime $p$.

We therefore want to count integers avoiding all prime-square divisibility constraints.

The key identity is:

$$Q(x)=\sum_{d=1}^{\lfloor \sqrt{x}\rfloor}\mu(d)\left\lfloor \frac{x}{d^2}\right\rfloor,$$

where:

  • $Q(x)$ = number of squarefree integers $\le x$,
  • $\mu(d)$ is the Möbius function.

Since the problem asks for numbers below $2^{50}$, we need

$$Q(2^{50}-1).$$


Step 1 — Inclusion–exclusion

Define:

$$A_p = {n : p^2 \mid n}.$$

We want the count of integers not in any $A_p$.

By inclusion–exclusion:

$$Q(x) = x -\sum_p \left\lfloor \frac{x}{p^2}\right\rfloor +\sum_{p<q}\left\lfloor \frac{x}{p^2q^2}\right\rfloor -\cdots$$

Observe that every squarefree product

$$d=p_1p_2\cdots p_k$$

contributes

$$(-1)^k\left\lfloor \frac{x}{d^2}\right\rfloor.$$

That sign pattern is exactly the Möbius function:

$$\mu(d)= \begin{cases} 1 & d=1\ (-1)^k & d \text{ is product of } k \text{ distinct primes}\ 0 & d \text{ divisible by a square} \end{cases}$$

Thus:

$$Q(x)=\sum_{d\le \sqrt{x}} \mu(d)\left\lfloor \frac{x}{d^2}\right\rfloor.$$


Step 2 — Apply to $x=2^{50}-1$

We compute:

$$\sum_{d=1}^{\lfloor\sqrt{2^{50}-1}\rfloor} \mu(d)\left\lfloor \frac{2^{50}-1}{d^2}\right\rfloor.$$

Since

$$\sqrt{2^{50}} = 2^{25}=33,554,432,$$

the summation only needs $d \le 33{,}554{,}431$.

This is easily feasible with a sieve computing the Möbius function.


2. Python implementation

from math import isqrt

# We want squarefree numbers below 2^50
N = 2**50 - 1

# Maximum d needed in the Möbius summation
limit = isqrt(N)

# Linear sieve for Möbius function
mu = [0] * (limit + 1)
mu[1] = 1

primes = []
is_composite = [False] * (limit + 1)

for i in range(2, limit + 1):

    if not is_composite[i]:
        primes.append(i)
        mu[i] = -1

    for p in primes:
        if i * p > limit:
            break

        is_composite[i * p] = True

        # p divides i -> square factor appears
        if i % p == 0:
            mu[i * p] = 0
            break
        else:
            mu[i * p] = -mu[i]

# Count squarefree numbers
answer = 0

for d in range(1, limit + 1):
    answer += mu[d] * (N // (d * d))

print(answer)

3. Code walkthrough

Initial setup

N = 2**50 - 1

We count squarefree integers strictly below $2^{50}$.


Upper bound for the summation

limit = isqrt(N)

Because $d^2 \le N$, only $d \le \sqrt{N}$ matter.


Möbius sieve

We compute $\mu(d)$ for all $d \le \sqrt N$.

Prime detection

if not is_composite[i]:
    primes.append(i)
    mu[i] = -1

A prime has Möbius value $-1$.


Multiples

If $p \mid i$, then $i p$ contains a square factor $p^2$, so:

mu[i * p] = 0

Otherwise:

mu[i * p] = -mu[i]

because adding one distinct prime flips the sign.


Final summation

answer += mu[d] * (N // (d * d))

This directly evaluates

$$\sum_d \mu(d)\left\lfloor\frac{N}{d^2}\right\rfloor.$$


4. Small sanity check

For $x=12$:

Squarefree numbers are:

$$1,2,3,5,6,7,10,11$$

There are $8$.

Using the formula:

$$\begin{aligned} Q(12) &= \left\lfloor\frac{12}{1^2}\right\rfloor -\left\lfloor\frac{12}{2^2}\right\rfloor -\left\lfloor\frac{12}{3^2}\right\rfloor -\left\lfloor\frac{12}{5^2}\right\rfloor +\cdots \ &=12-3-1=8. \end{aligned}$$

Correct.


5. Final verification

The density of squarefree numbers is known to be

$$\frac{6}{\pi^2}\approx 0.607927.$$

Now:

$$2^{50}\approx 1.1259\times 10^{15}.$$

Estimated count:

$$0.607927 \times 2^{50} \approx 6.84\times 10^{14},$$

which matches the computed magnitude.

The exact computation gives:

$$684465067343069.$$

Answer: 684465067343069