Project Euler Problem 668

A positive integer is called square root smooth if all of its prime factors are strictly less than its square root.

Project Euler Problem 668

Solution

Answer: 2811077773

Let $P(n)$ denote the largest prime factor of $n$ (with $P(1)$ undefined, but $1$ is declared square-root smooth by definition).

A number $n$ is square-root smooth iff every prime factor is strictly smaller than $\sqrt n$.

Equivalently,

$$P(n)<\sqrt n$$

or

$$P(n)^2<n.$$

We want the number of integers $n\le 10^{10}$ satisfying this.


1. Mathematical analysis

Step 1: Characterize the complement

It is easier to count the numbers that are not square-root smooth.

A number $n>1$ is not square-root smooth iff

$$P(n)\ge \sqrt n.$$

Write

$$n = pk,$$

where $p=P(n)$ is the largest prime factor.

Then

$$p \ge \sqrt n \quad\Longrightarrow\quad p^2 \ge n = pk \quad\Longrightarrow\quad p \ge k.$$

So every non-smooth number has a unique representation

$$n = pk$$

with:

  • $p$ prime,
  • $k \le p$,
  • $pk \le x$.

Conversely, every such pair produces a non-smooth number, because $p$ is a prime factor and

$$p \ge k ;\Rightarrow; p^2 \ge pk=n ;\Rightarrow; p\ge \sqrt n.$$

Thus the non-smooth numbers are in bijection with pairs

$$(p,k)$$

such that

$$p \text{ prime},\qquad k\le p,\qquad pk\le x.$$


Step 2: Count the complement

Fix $k$.

Then $p$ must satisfy

$$k \le p \le \frac{x}{k}.$$

So the number of valid primes is

$$\pi!\left(\frac{x}{k}\right)-\pi(k-1),$$

where $\pi(n)$ is the prime-counting function.

Also, necessarily $k\le \sqrt x$, otherwise $k^2>x$.

Hence the number of non-smooth integers is

$$N(x) = \sum_{k\le \sqrt x} \left( \pi!\left(\frac{x}{k}\right)-\pi(k-1) \right).$$

Therefore the desired count is

$$S(x)=x-N(x).$$

So:

$$\boxed{ S(x) = x- \sum_{k\le \sqrt x} \left( \pi!\left(\frac{x}{k}\right)-\pi(k-1) \right) }$$

with $x=10^{10}$.


Step 3: Check the example $x=100$

Using the formula:

$$N(100) = \sum_{k=1}^{10} \left( \pi(100/k)-\pi(k-1) \right) =71.$$

Hence

$$S(100)=100-71=29,$$

matching the problem statement.


2. Python implementation

from functools import lru_cache
from math import isqrt

# ------------------------------------------------------------
# Sieve for small primes
# ------------------------------------------------------------

LIMIT = 5_000_000

sieve = bytearray(b"\x01") * (LIMIT + 1)
sieve[0:2] = b"\x00\x00"

for i in range(2, isqrt(LIMIT) + 1):
    if sieve[i]:
        sieve[i * i : LIMIT + 1 : i] = b"\x00" * (
            ((LIMIT - i * i) // i) + 1
        )

primes = [i for i, v in enumerate(sieve) if v]

# pi_small[n] = number of primes <= n
pi_small = [0] * (LIMIT + 1)

count = 0
for i in range(LIMIT + 1):
    if sieve[i]:
        count += 1
    pi_small[i] = count

# ------------------------------------------------------------
# Lehmer prime counting function
# Computes pi(n) quickly for large n
# ------------------------------------------------------------

@lru_cache(None)
def phi(x, s):
    """Count integers <= x not divisible by first s primes."""
    if s == 0:
        return x
    return phi(x, s - 1) - phi(x // primes[s - 1], s - 1)

@lru_cache(None)
def pi(n):
    """Lehmer prime counting function."""
    if n <= LIMIT:
        return pi_small[n]

    a = pi(int(n ** 0.25))
    b = pi(isqrt(n))
    c = pi(int(n ** (1 / 3)))

    # Main Lehmer formula
    result = phi(n, a) + ((b + a - 2) * (b - a + 1)) // 2

    for i in range(a + 1, b + 1):
        p = primes[i - 1]
        w = n // p

        result -= pi(w)

        if i <= c:
            bi = pi(isqrt(w))

            for j in range(i, bi + 1):
                result -= pi(w // primes[j - 1]) - (j - 1)

    return result

# ------------------------------------------------------------
# Count square-root smooth numbers
# ------------------------------------------------------------

def square_root_smooth_count(x):
    root = isqrt(x)

    non_smooth = 0

    for k in range(1, root + 1):
        non_smooth += pi(x // k) - pi(k - 1)

    return x - non_smooth

# ------------------------------------------------------------
# Compute the Project Euler answer
# ------------------------------------------------------------

x = 10_000_000_000

answer = square_root_smooth_count(x)

print(answer)

3. Code walkthrough

Sieve section

We first generate all primes up to $5\times10^6$ using the sieve of Eratosthenes.

sieve = bytearray(...)

stores primality.

We also build:

pi_small[n]

which gives $\pi(n)$ instantly for small $n$.


phi(x, s)

This function counts numbers $\le x$ not divisible by the first $s$ primes.

It satisfies the recursion:

$$\phi(x,s)=\phi(x,s-1)-\phi!\left(\left\lfloor \frac{x}{p_s}\right\rfloor,s-1\right).$$

This is a standard ingredient in Lehmer’s prime-counting algorithm.


pi(n)

This implements the Lehmer prime-counting function.

It computes $\pi(n)$ in roughly $O(n^{2/3})$ time, fast enough for $10^{10}$.


Final counting loop

We use the formula

$$N(x) = \sum_{k\le \sqrt x} \left( \pi(x/k)-\pi(k-1) \right).$$

implemented as:

for k in range(1, root + 1):
    non_smooth += pi(x // k) - pi(k - 1)

Then:

return x - non_smooth

gives the number of square-root smooth integers.


4. Final verification

For the sample:

$$S(100)=29,$$

matching the statement exactly.

For $10^{10}$, the computed value is:

$$2,811,077,773.$$

This is plausible:

  • it is significantly smaller than $10^{10}$,
  • but still billions, since many integers have all prime factors below $\sqrt n$.

The computation is internally consistent and matches the sample test.


Answer: 2811077773