Project Euler Problem 578

Solution

Answer: 9219696799346

Let

$$n = p_1^{a_1}p_2^{a_2}\cdots p_k^{a_k}, \qquad p_1<p_2<\cdots<p_k,$$

with all $a_i>0$.

We want to count integers for which

$$a_1 \ge a_2 \ge \cdots \ge a_k.$$

Define $C(N)$ as the number of such integers $\le N$.

The problem gives:

$$C(100)=94,\qquad C(10^6)=922052.$$

We must compute $C(10^{13})$.

---

Key mathematical insight

Split the factorization into:

• a powerful part (all exponents $\ge 2$),
• and a squarefree tail (all remaining exponents equal to $1$).

Example:

$$360 = 2^3\cdot 3^2\cdot 5$$

has powerful part $2^3 3^2$, and squarefree tail $5$.

Because exponents must be nonincreasing, once exponent $1$ appears, every later exponent must also be $1$.

So every valid integer has the form

$$p_1^{e_1}p_2^{e_2}\cdots p_t^{e_t}\cdot q_1q_2\cdots q_r,$$

where

$$e_1\ge e_2\ge \cdots \ge e_t \ge 2,$$

and all $q_i$ are distinct primes larger than $p_t$.

Thus the problem naturally separates into:

1. recursively enumerating the powerful part;
2. counting squarefree tails efficiently.

---

Counting squarefree numbers

Let

$$Q(x)$$

be the number of squarefree integers $\le x$.

Using Möbius inversion,

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

We precompute the Möbius function with a sieve.

---

Squarefree tails with restricted primes

Suppose the powerful part ended at prime index start_idx.

Then the squarefree tail may only use primes

$$\ge p_{\text{start_idx}}.$$

Define

$$f(x,a)$$

= number of squarefree integers $\le x$ whose prime factors are all at least the $a$-th prime.

We use a smallest-prime-factor decomposition:

$$f(x,a) = Q(x) - \sum_{i=0}^{a-1} f!\left(\left\lfloor \frac{x}{p_i}\right\rfloor, i+1\right).$$

This recursively removes squarefree numbers whose smallest prime factor is forbidden.

---

Recursive enumeration

Now recursively build the powerful part.

State:

$$\text{count}(L,\text{start},\text{maxExp})$$

means:

• remaining limit $L$,
• allowed primes begin at index start,
• next exponent may not exceed maxExp.

If maxExp = 1, the remainder must be squarefree, so the answer is simply

$$f(L,\text{start}).$$

Otherwise:

• count all squarefree tails immediately;
• then try every prime $p$ and exponent $e\in[2,\text{maxExp}]$;
• recurse on

$$\left\lfloor \frac{L}{p^e}\right\rfloor$$

with exponent bound $e$.

This visits only feasible exponent patterns, and is fast enough for $10^{13}$.

---

Python implementation

import math
from functools import lru_cache
from array import array

N_TARGET = 10**13
SIEVE_LIMIT = int(math.isqrt(N_TARGET)) + 10

# ------------------------------------------------------------
# Prime sieve + Möbius function
# ------------------------------------------------------------

def build_primes_and_mobius(n):
    lp = array("I", [0]) * (n + 1)
    mu = array("b", [0]) * (n + 1)
    mu[1] = 1

    primes = []

    for i in range(2, n + 1):
        if lp[i] == 0:
            lp[i] = i
            primes.append(i)
            mu[i] = -1

        for p in primes:
            ip = i * p
            if ip > n:
                break

            lp[ip] = p

            if p == lp[i]:
                mu[ip] = 0
                break

            mu[ip] = -mu[i]

    prefix_mu = array("i", [0]) * (n + 1)

    s = 0
    for i in range(1, n + 1):
        s += mu[i]
        prefix_mu[i] = s

    return primes, prefix_mu

PRIMES, PREFIX_MU = build_primes_and_mobius(SIEVE_LIMIT)

# ------------------------------------------------------------
# Count squarefree numbers <= x
# ------------------------------------------------------------

@lru_cache(maxsize=None)
def squarefree_upto(x):
    if x <= 0:
        return 0

    r = int(math.isqrt(x))

    res = 0
    i = 1

    while i <= r:
        t = x // (i * i)
        j = int(math.isqrt(x // t))

        res += t * (PREFIX_MU[j] - PREFIX_MU[i - 1])

        i = j + 1

    return res

# ------------------------------------------------------------
# Count squarefree numbers whose prime factors are all
# >= PRIMES[start_idx]
# ------------------------------------------------------------

@lru_cache(maxsize=None)
def squarefree_min_prime_index(x, start_idx):

    if x <= 0:
        return 0

    if x == 1:
        return 1

    if start_idx == 0:
        return squarefree_upto(x)

    if start_idx < len(PRIMES) and PRIMES[start_idx] > x:
        return 1

    total = squarefree_upto(x)

    for i in range(start_idx):
        p = PRIMES[i]

        if p > x:
            break

        total -= squarefree_min_prime_index(x // p, i + 1)

    return total

# ------------------------------------------------------------
# Main recursion
# ------------------------------------------------------------

@lru_cache(maxsize=None)
def count_dpowers(limit, start_idx, max_exp):

    if limit <= 0:
        return 0

    if limit == 1:
        return 1

    # only exponent 1 allowed now
    if max_exp <= 1:
        return squarefree_min_prime_index(limit, start_idx)

    # all-squarefree continuation
    res = squarefree_min_prime_index(limit, start_idx)

    for i in range(start_idx, len(PRIMES)):

        p = PRIMES[i]

        if p * p > limit:
            break

        pe = p * p
        e = 2

        while e <= max_exp and pe <= limit:

            res += count_dpowers(limit // pe, i + 1, e)

            e += 1
            pe *= p

    return res

def C(n):
    return count_dpowers(n, 0, int(math.log2(n)))

# ------------------------------------------------------------
# Verification
# ------------------------------------------------------------

print(C(100))      # 94
print(C(10**6))    # 922052

# Final answer
print(C(10**13))

---

Code walkthrough

1. Möbius sieve

We compute:

• all primes up to $\sqrt{10^{13}}$,
• the Möbius function $\mu(n)$,
• prefix sums of $\mu$.

These are needed for fast squarefree counting.

---

2. squarefree_upto(x)

Implements

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

Grouping equal quotient values makes it fast.

---

3. squarefree_min_prime_index

Counts squarefree numbers using only sufficiently large primes.

This avoids explicitly enumerating enormous numbers of squarefree products.

---

4. count_dpowers

Recursively builds valid exponent patterns:

• exponents never increase,
• once we reach exponent $1$, everything else is squarefree.

This exactly matches the definition of decreasing prime powers.

---

Verification

The program reproduces the given values:

$$C(100)=94$$

and

$$C(10^6)=922052.$$

Then it computes:

$$C(10^{13}) = 9219696799346.$$

This is plausible:

• total numbers: $10^{13}$,
• almost all numbers satisfy the condition,
• but a nontrivial fraction fail due to exponent increases such as

$2^1 3^2$, $2^1 5^2$, etc.

The result is slightly below $10^{13}$, as expected.

Answer: 9219696799346