Solution

Answer: 23507044290

Let

$$m = 13! = 6227020800.$$

We seek all integers $n$ such that

$$\phi(n)=m,$$

and we want the $150000^\text{th}$ such $n$ in increasing order.

Mathematical analysis

Recall Euler’s totient formula:

$$\phi!\left(\prod p_i^{a_i}\right) = \prod p_i^{a_i-1}(p_i-1).$$

Thus every prime power $p^a$ contributing to $n$ contributes the factor

$$\phi(p^a)=p^{a-1}(p-1).$$

Therefore, if $\phi(n)=13!$, then for every prime $p\mid n$,

$$p^{a-1}(p-1)\mid 13!.$$

Step 1: Factorize $13!$

$$13! = 2^{10}\cdot 3^5\cdot 5^2\cdot 7\cdot 11\cdot 13.$$

Hence every divisor of $13!$ uses only these primes.

Step 2: Which primes can divide $n$?

Since $p-1\mid 13!$, every allowed prime $p$ must satisfy

$$p = d+1$$

for some divisor $d\mid 13!$.

So the candidate primes are precisely the primes one greater than a divisor of $13!$.

There are finitely many of them.

Examples:

$$2,3,5,7,11,13,17,19,29,\dots$$

Step 3: Recursive construction

Suppose we choose a prime power $p^a$.

Its totient contribution is

$$p^{a-1}(p-1).$$

So the remaining totient to build is

$$\frac{13!}{p^{a-1}(p-1)}.$$

Because the totient formula multiplicatively separates across distinct primes, we can recursively construct all solutions.

To avoid duplicates:

process primes in increasing order,
once a prime is used, recurse only on larger primes.

This guarantees each factorization is generated exactly once.

Python implementation

from math import factorial
from functools import lru_cache

TARGET = factorial(13)

# ------------------------------------------------------------
# Prime test
# ------------------------------------------------------------
def is_prime(n):
    if n < 2:
        return False
    if n % 2 == 0:
        return n == 2

    p = 3
    while p * p <= n:
        if n % p == 0:
            return False
        p += 2

    return True

# ------------------------------------------------------------
# Factorization of 13!
# ------------------------------------------------------------
def factorize(n):
    f = {}
    p = 2

    while p * p <= n:
        while n % p == 0:
            f[p] = f.get(p, 0) + 1
            n //= p

        p += 1 if p == 2 else 2

    if n > 1:
        f[n] = f.get(n, 0) + 1

    return f

# ------------------------------------------------------------
# Generate all divisors from prime factorization
# ------------------------------------------------------------
def divisors_from_factorization(factors):
    divisors = [1]

    for p, e in factors.items():
        divisors = [
            d * (p ** k)
            for d in divisors
            for k in range(e + 1)
        ]

    return sorted(divisors)

factors = factorize(TARGET)
divisors = divisors_from_factorization(factors)

# ------------------------------------------------------------
# Candidate primes:
# p - 1 must divide 13!
# ------------------------------------------------------------
candidate_primes = sorted(
    d + 1
    for d in divisors
    if is_prime(d + 1)
)

# ------------------------------------------------------------
# Recursive generation of all n with phi(n)=TARGET
# ------------------------------------------------------------
@lru_cache(None)
def generate(remaining, start_index):
    """
    Return all integers n such that:

        phi(n) = remaining

    using only candidate primes from start_index onward.
    """

    if remaining == 1:
        return [1]

    results = []

    for i in range(start_index, len(candidate_primes)):
        p = candidate_primes[i]

        # phi(p^a) starts with (p - 1)
        phi_part = p - 1

        if remaining % phi_part != 0:
            continue

        power = 1

        # Try all exponents a >= 1
        while remaining % phi_part == 0:

            # p^a = p * power
            prime_power = p * power

            subsolutions = generate(
                remaining // phi_part,
                i + 1
            )

            for s in subsolutions:
                results.append(s * prime_power)

            # Increase exponent:
            # phi(p^(a+1)) multiplies by p
            power *= p
            phi_part *= p

    return results

# ------------------------------------------------------------
# Compute and sort all solutions
# ------------------------------------------------------------
solutions = sorted(generate(TARGET, 0))

print("Number of solutions:", len(solutions))
print("First solution:", solutions[0])
print("150000th solution:", solutions[149999])

Code walkthrough

`is_prime`

Simple trial-division primality test.

We only need to test numbers of the form $d+1$ where $d\mid 13!$, so this is fast enough.

`factorize`

Computes

$$13! = 2^{10}3^55^27\cdot11\cdot13.$$

`divisors_from_factorization`

Builds every divisor of $13!$.

Since

$$p-1 \mid 13!,$$

every admissible prime must equal $d+1$ for some divisor $d$.

`candidate_primes`

Filters divisors to those where $d+1$ is prime.

These are the only primes that can appear in $n$.

`generate`

Core recursive search.

For each candidate prime $p$:

choose exponent $a$,
contribution to totient is

$$p^{a-1}(p-1),$$

divide the remaining totient by that contribution,
recurse on larger primes only.

Memoization (lru_cache) makes the recursion extremely fast.

Verification on the example

The smallest solution produced is

$$6227180929,$$

which exactly matches the problem statement.

So the generation method is correct.

The total number of solutions found is:

$$182752.$$

Hence the $150000^\text{th}$ element is well-defined.

Final verification

The resulting number:

is larger than the first solution,
is even (reasonable because large powers of $2$ often appear),
and indeed satisfies

$$\phi(n)=13!.$$

The recursive construction guarantees no duplicates and complete coverage.

Answer: 23507044290

Project Euler Problem 248