Solution

Answer: 172023848

Let

$$\pi(n)=#{p\le n : p \text{ prime}}.$$

A $\pi$-sequence is a finite sequence

$$u=(u_0,u_1,\dots,u_m)$$

with $m\ge1$, $u_{i+1}=\pi(u_i)$, and all terms positive.

For each sequence $u$, let $c(u)$ be the number of non-prime elements in $u$.

We define:

$$p(n,k)=#{u : u_0\le n,\ c(u)=k},$$

and

$$P(n)=\prod_{k:,p(n,k)>0} p(n,k).$$

We are asked for:

$$P(10^8)\pmod{1,000,000,007}.$$

1. Mathematical analysis

1.1 Structure of a $\pi$-sequence

Starting from $x=u_0$, the sequence is forced:

$$x,\ \pi(x),\ \pi(\pi(x)),\ \dots$$

until it reaches $1$, because $\pi(1)=0$ and sequences must stay positive.

Thus every starting value $x\ge2$ generates exactly one maximal chain.

Example:

$$10\to4\to2\to1.$$

The valid $\pi$-sequences are all prefixes of length at least $2$:

$$(10,4),\quad (10,4,2),\quad (10,4,2,1).$$

1.2 Recursive description

Define $V(x)$ to be the multiset of all values $c(u)$ for sequences starting at $x$.

If $x>1$, let

$$y=\pi(x).$$

Let

$$b(x)= \begin{cases} 0,&x\text{ prime}\ 1,&x\text{ non-prime} \end{cases}$$

Then:

the shortest sequence $(x,y)$ contributes

$$b(x)+b(y),$$

every longer sequence is obtained by prepending $x$ to a sequence from $y$, adding $b(x)$ to the non-prime count.

Therefore:

$$V(x)={b(x)+b(y)}\cup{b(x)+t:t\in V(y)}.$$

This gives a very small recursion depth, because repeated application of $\pi$ shrinks numbers extremely quickly.

For $10^8$:

$$10^8\to5761455\to412849\to\cdots\to1$$

so the maximum depth is tiny (well under 16).

1.3 Grouping by equal $\pi(x)$

A crucial observation:

$$\pi(x)=y$$

holds for all

$$x\in[p_y,\ p_{y+1}-1],$$

where $p_y$ is the $y$-th prime.

Inside this interval:

exactly one number is prime ($p_y$),
the rest are composite.

Thus all numbers with the same $\pi(x)$ fall into only two types:

the prime $p_y$,
the composites between $p_y$ and $p_{y+1}$.

So we only need to compute $V(y)$ once for each $y\le \pi(10^8)$.

Since

$$\pi(10^8)=5,761,455,$$

this is completely feasible.

2. Efficient algorithm

We sieve all primes up to $10^8$.

Then for every $y\le\pi(10^8)$:

compute $V(y)$ from the recursion,
add contributions from:

the prime $p_y$,
the composite numbers with $\pi(x)=y$.

Because the recursion depth is tiny, we store only a short vector of counts.

Complexity:

sieve: $O(n\log\log n)$,
DP: essentially linear in $\pi(n)$.

3. Python implementation

MOD = 1_000_000_007
N = 10**8
K = 16                      # safely larger than maximum depth

# ------------------------------------------------------------
# Odd-only sieve
# ------------------------------------------------------------

size = N // 2
sieve = bytearray(b'\x01') * size
sieve[0] = 0

limit = int(N**0.5)

for i in range(3, limit + 1, 2):
    if sieve[i // 2]:
        start = i * i // 2
        sieve[start::i] = b'\x00' * (((size - start - 1) // i) + 1)

primes = [2] + [2 * i + 1 for i in range(1, size) if sieve[i]]

m = len(primes)             # pi(10^8)

# ------------------------------------------------------------
# primality up to m
# ------------------------------------------------------------

is_prime_small = bytearray(m + 1)

for p in primes:
    if p <= m:
        is_prime_small[p] = 1

# ------------------------------------------------------------
# pi(x) for x <= m
# ------------------------------------------------------------

pi_small = [0] * (m + 1)

c = 0
for i in range(1, m + 1):
    c += is_prime_small[i]
    pi_small[i] = c

# ------------------------------------------------------------
# DP storage
#
# C[y][k] = number of sequences starting at y
#            having exactly k non-primes
# ------------------------------------------------------------

C = [[0] * K for _ in range(m + 1)]

# final p(n,k)
P = [0] * K

for y in range(1, m + 1):

    # build V(y)
    if y > 1:
        z = pi_small[y]

        by = 0 if is_prime_small[y] else 1
        bz = 0 if is_prime_small[z] else 1

        # shortest sequence (y, z)
        C[y][by + bz] += 1

        # longer sequences
        for k in range(K - by):
            C[y][by + k] += C[z][k]

    # --------------------------------------------------------
    # Contributions from x with pi(x)=y
    # --------------------------------------------------------

    bz = 0 if is_prime_small[y] else 1

    # prime x = p_y
    P[bz] += 1
    for k in range(K):
        P[k] += C[y][k]

    # composite x with pi(x)=y
    if y < m:
        gap = primes[y] - primes[y - 1] - 1
    else:
        gap = N - primes[-1]

    P[1 + bz] += gap

    for k in range(K - 1):
        P[k + 1] += gap * C[y][k]

# ------------------------------------------------------------
# Final product
# ------------------------------------------------------------

ans = 1

for v in P:
    if v:
        ans = (ans * (v % MOD)) % MOD

print(ans)

4. Code walkthrough

Sieve

We use an odd-only sieve:

index $i$ represents number $2i+1$,
memory stays manageable even for $10^8$.

This produces all primes up to $10^8$.

Computing $V(y)$

For each $y$:

z = pi_small[y]

Then:

C[y][by + bz] += 1

adds the shortest sequence $(y,z)$.

C[y][by + k] += C[z][k]

extends every sequence from $z$.

Grouped counting

All numbers with $\pi(x)=y$ are handled together.

There is:

one prime $x=p_y$,
many composites between consecutive primes.

This reduces the work enormously.

5. Verification on small cases

The program reproduces:

$$P(10)=648,$$

and

$$P(100)=31038676032,$$

matching the statement.

6. Final computation

Running the algorithm for $N=10^8$ gives:

$$P(10^8)\equiv 172023848 \pmod{1,000,000,007}.$$

Answer: 172023848

Project Euler Problem 609