Problem 248: Euler's Totient Function Equals 13!

Project Euler Problem 248 Solution

EulerSolve provides an optimized solution for Project Euler Problem 248, Euler's Totient Function Equals 13!, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The problem asks for all positive integers $n$ such that $$\varphi(n)=13!,$$ then orders those integers increasingly and asks for the $150000$-th entry of that list. This is an inverse-totient problem: instead of evaluating $\varphi(n)$ for a known $n$, we must recover every $n$ whose totient equals the fixed target $$13!=6227020800=2^{10}\cdot 3^5\cdot 5^2\cdot 7\cdot 11\cdot 13.$$ A brute-force scan over $n$ is not realistic. The workable route is to use the multiplicative structure of $\varphi$, restrict the possible prime divisors of $n$, and enumerate only exact factorizations of the target totient. Mathematical Approach If $$n=\prod_{k=1}^{m} p_k^{e_k},$$ then Euler's product formula gives $$\varphi(n)=\prod_{k=1}^{m}\varphi\!\left(p_k^{e_k}\right)=\prod_{k=1}^{m}(p_k-1)p_k^{e_k-1}.$$ The entire search can therefore be phrased as: find every collection of prime powers whose individual totient contributions multiply exactly to $13!$. Factor the target once and exploit it everywhere The fixed target $$T=13!=2^{10}\cdot 3^5\cdot 5^2\cdot 7\cdot 11\cdot 13$$ already tells us that every admissible contribution to $\varphi(n)$ must divide $T$. Because the prime factorization of $T$ is tiny, all divisor-based pruning is cheap....

Detailed mathematical approach

Problem Summary

Euler's totient function $\varphi(n)$ counts the integers $k$ with $1\le k\le n$ that are coprime to $n$, that is, with $\gcd(k,n)=1$. We consider all positive integers $n$ with $\varphi(n)=13!=6\,227\,020\,800$, listed in increasing order. The statement gives the first of them, $6\,227\,180\,929$, and asks for the $150\,000$th. ^[1]

The whole solution in one sentence: only primes $p$ with $p-1$ dividing $13!$ can divide such an $n$, so we combine powers of these primes in a depth-first search until their totients multiply to exactly $13!$, sort the numbers found and pick the $150\,000$th.

Mathematical Approach

We write $T$ for the prescribed value of $\varphi$; the problem has $T=13!$, and the small examples use $T=12$.

1. The totient of a factorisation

The values $\varphi(1)$ to $\varphi(12)$ are $1,1,2,2,4,2,6,4,6,4,10,4$; for instance $\varphi(12)=4$ counts $1,5,7,11$. For a prime power $p^e$ with $e\ge1$, the numbers up to $p^e$ that are not coprime to it are the $p^{e-1}$ multiples of $p$, so

$$\varphi(p^e)=p^e-p^{e-1}=(p-1)\,p^{e-1};$$

for example $\varphi(8)=4$ and $\varphi(9)=6$. For coprime $a$ and $b$, the Chinese remainder theorem matches the residues modulo $ab$ one to one with the pairs of residues modulo $a$ and modulo $b$, and a number is coprime to $ab$ exactly when it is coprime to both. Hence $\varphi(ab)=\varphi(a)\,\varphi(b)$, as in $\varphi(12)=\varphi(4)\,\varphi(3)=2\cdot2$. ^[2] By the fundamental theorem of arithmetic, every $n\gt1$ is, in exactly one way, a product $p_1^{e_1}\cdots p_k^{e_k}$ with primes $p_1\lt\dots\lt p_k$ and exponents $e_j\ge1$, ^[3] and then

$$\varphi(n)=\prod_{j=1}^{k}(p_j-1)\,p_j^{\,e_j-1}.$$

We call $\varphi(p^e)$ the piece of $p^e$. Inverting $\varphi$ means choosing powers of distinct primes whose pieces multiply to exactly $T$. Two consequences are used below. First, $\varphi(n)$ is even for $n\ge3$: an odd prime factor contributes the even factor $p-1$, and otherwise $n=2^e$ with $e\ge2$; so $\varphi(n)=1$ only for $n=1$ and $n=2$. Second, if $n$ has $s$ distinct odd prime factors, then $2^s$ divides $\varphi(n)$. ^[2]

2. Which primes can occur

Let $\varphi(n)=T$, and let $p^e$ be the exact power of a prime $p$ in $n$. Its piece $(p-1)p^{e-1}$ is a factor of the product above, so it divides $T$. Hence $p-1$ divides $T$, and if $e\ge2$, then $p$ divides $T$ too; more precisely $e\le\nu_p(T)+1$, where $\nu_p(T)$ is the exponent of $p$ in $T$. ^[4] So the only primes that can occur are the candidates: the numbers $d+1$ that are prime, where $d$ runs through the divisors of $T$. For $T=12$, the divisors $1,2,3,4,6,12$ give $2,3,4,5,7,13$, and the candidates are $2,3,5,7,13$. For the problem,

$$13!=2^{10}\cdot3^{5}\cdot5^{2}\cdot7\cdot11\cdot13$$

has $(10+1)(5+1)(2+1)(1+1)^3=1584$ divisors, and $459$ of the numbers $d+1$ are prime, from $2$ up to $3\,113\,510\,401=13!/2+1$. Since $13!+1=83\cdot75\,024\,347$ is not prime, no solution is a prime. Among the primes below $100$, only $47$, $59$ and $83$ are not candidates, because $46=2\cdot23$, $58=2\cdot29$ and $82=2\cdot41$ contain a prime larger than $13$. Only the six primes up to $13$ divide $13!$, so all other candidates occur at most to the first power, and the rule allows at most $2^{11}$, $3^6$, $5^3$, $7^2$, $11^2$ and $13^2$.

3. A search over prime powers

List the candidates as $q_0\lt q_1\lt q_2\lt\dots$. A state of the search is a triple $(r,i,m)$: $m$ is the product of the prime powers chosen so far, $r=T/\varphi(m)$ is the part of $T$ still to be produced, and $i$ is the index of the smallest candidate that may still be used. The search starts at $(T,0,1)$. A state with $r=1$ is a leaf, and its $m$ is recorded as a solution. Every other state has one child for each index $j\ge i$ and each exponent $e\ge1$ whose piece divides $r$:

$$\Bigl(\frac{r}{(q_j-1)\,q_j^{\,e-1}},\ j+1,\ m\,q_j^{\,e}\Bigr).$$

A state with $r\gt1$ and no children is a dead end. The search walks through this tree depth first and abandons a branch as soon as no piece fits, a method called backtracking. ^[7] Each piece of $q_j$ is $q_j$ times the previous one, so the exponents that fit form a run $e=1,2,\dots$, and the loop over $e$ can stop at the first failure.

4. Worked example: $\varphi(n)=12$

The candidates $2,3,5,7,13$ have the pieces $\varphi(2)=1$, $\varphi(4)=2$, $\varphi(8)=4$, $\varphi(3)=2$, $\varphi(9)=6$, $\varphi(5)=4$, $\varphi(7)=6$ and $\varphi(13)=12$, while $\varphi(16)=8$, $\varphi(27)=18$ and $\varphi(25)=20$ do not divide $12$. The root $(12,0,1)$ has eight children, treated in this order:

$m=2$, $r=12$: below it, $2\cdot3$ leaves $r=6$, and then $2\cdot3\cdot7=42$ is a leaf; $2\cdot9$ and $2\cdot7$ leave $r=2$, and $2\cdot5$ leaves $r=3$: three dead ends; finally $2\cdot13=26$ is a leaf.
$m=4$, $r=6$: $4\cdot3$ leaves $r=3$, a dead end, while $4\cdot9=36$ and $4\cdot7=28$ are leaves.
$m=8$, $r=3$: a dead end.
$m=3$, $r=6$: $3\cdot7=21$ is a leaf. $m=9$, $r=2$: a dead end.
$m=5$, $r=3$ and $m=7$, $r=2$: dead ends. $m=13$, $r=1$: a leaf.

So the search visits $19$ states, among them $6$ leaves and $8$ dead ends, and finds $42,26,36,28,21,13$, in increasing order $13,21,26,28,36,42$. ^[4] The dead ends with the odd value $r=3$ could never be completed, because the piece of every odd prime is even. Each solution is reached along its factorisation with increasing primes: $42$ only as $2\cdot3\cdot7$, never as $3\cdot2\cdot7$.

5. Every solution exactly once

Claim. For $T\ge2$, the leaves of the search are exactly the integers $n$ with $\varphi(n)=T$, and each of them is reached along exactly one path.

Proof. In every state, $\varphi(m)\,r=T$, and every prime factor of $m$ is a candidate $q_j$ with $j\lt i$. Both hold at the root. A child multiplies $m$ by a power of a candidate $q_j$ with $j\ge i$, which is larger than all prime factors of $m$, hence coprime to $m$, so by Section 1 $\varphi$ of the new product is $\varphi(m)$ times the piece, while $r$ is divided by the piece. At a leaf $r=1$, hence $\varphi(m)=T$. Conversely, let $\varphi(n)=T$ with $n=p_1^{e_1}\cdots p_k^{e_k}$ and $p_1\lt\dots\lt p_k$. By Section 2 all $p_j$ are candidates. Follow the path that chooses $p_1^{e_1}$, then $p_2^{e_2}$, and so on. After $j$ steps, $r$ is the product of the pieces of $p_{j+1}^{e_{j+1}},\dots,p_k^{e_k}$, so the next piece divides it, and the step exists. For $j\lt k$ this $r$ is at least $2$: for $j=0$ it is $T$, and for $j\ge1$ it is a product of pieces of odd primes, each of them even. So the path does not end in a leaf too early, and it ends in the leaf with $m=n$. Two different paths choose different sets of prime powers, whose products differ by unique factorisation, so no solution is found twice. $\square$

For $T=1$ the root itself is a leaf: the search records $1$ and never reaches the second solution $2$, whose piece is $\varphi(2)=1$. The target $13!$ is not affected.

6. How large the search tree is

The same argument describes all states: for $T\ge2$ they correspond one to one to the integers $m\ge1$ with $\varphi(m)\mid T$. A state has $\varphi(m)=T/r$, a divisor of $T$. Conversely, each initial part $p_1^{e_1}\cdots p_j^{e_j}$ of the factorisation of such an $m$ has a totient dividing $\varphi(m)$, so every step along the path of $m$ is allowed, and as in Section 5 the path does not stop early. For $T=12$, the values $\varphi(m)=1,2,3,4,6,12$ are taken by $2,3,0,4,4,6$ integers $m$, which are the $19$ states of Section 4.

For $13!$ there are $2\,795\,333$ states: $182\,752$ leaves, $2\,239\,009$ dead ends and $373\,572$ states with children. Of the dead ends, $923\,338$ have an odd $r\gt1$; choosing $2^{11}$ at the root, for example, uses up the factor $2^{10}$ of $13!$ and leaves the odd number $6\,081\,075$. The tree is shallow: $13!$ contains the factor $2$ exactly ten times, and by Section 1 each distinct odd prime factor of $m$ contributes a factor $2$ to $\varphi(m)$, which divides $13!$ in every state; so $m$ has at most ten distinct odd prime factors. In fact no state has more than nine distinct prime factors, so the recursion is at most ten calls deep.

7. Where the solutions lie

Dividing the formula of Section 1 by $n$ gives

$$\frac{n}{\varphi(n)}=\prod_{p\mid n}\frac{p}{p-1}.$$

Every factor exceeds $1$, so every solution is larger than $T$. All prime factors are candidates, so the product over all $459$ candidates, about $9.55$, bounds the ratio, and every solution lies below $9.55\cdot13!\lt6\cdot10^{10}\lt2^{36}$. The largest solution is $37\,020\,293\,310$${}=2\cdot3\cdot5\cdot7^2\cdot11$${}\cdot13^2\cdot19\cdot23\cdot31$, about $5.9$ times $13!$.

The smallest solution is close to $T$ for a reason. If $n$ has at least three distinct prime factors, the smallest of them, $p$, satisfies $(p-1)^3\le T$, so $p-1\le T^{1/3}\approx1840$ and

$$\begin{aligned}n&\ge T\cdot\frac{p}{p-1}=T+\frac{T}{p-1}\\&\ge T+T^{2/3},\end{aligned}$$

which exceeds $T+3\cdot10^6$. If a prime occurs squared, it is at most $13$, and $n\ge\frac{13}{12}T$. A prime power alone is impossible: $13!+1$ is not prime, and a piece $(p-1)p^{e-1}$ with $p\le13$ has at most three distinct prime factors, while $13!$ has six. What remains are products $n=pq$ of two primes with $(p-1)(q-1)=T$, for which $n=T+p+q-1$. There are $71$ such pairs. The one with the smallest sum, $p=66\,529$ and $q=93\,601$, gives $n=T+160\,129$, far below the bounds of the other cases, so it is the smallest solution, the value $6\,227\,180\,929$ of the statement.

A factor $2$ contributes $2$ to the ratio, so every even solution is at least $2T$: the first $42\,970$ solutions are odd, and the smallest even one is $4\cdot3\,113\,510\,401$, built from the largest candidate. Since $\varphi(2n)=\varphi(n)$ for odd $n$, ^[2] the $52\,827$ odd solutions pair off with the $52\,827$ solutions that are twice an odd number, and the other $77\,098$ solutions are multiples of $4$.

8. Sorting, and counting without listing

The search does not produce the solutions in increasing order, as $T=12$ showed, so all $182\,752$ solutions are collected and sorted, and the $150\,000$th entry is read off. A second method counts the solutions without listing them, by dynamic programming over the divisors of $T$. ^[4] It processes the candidates one at a time and keeps, for every divisor $d$ of $T$, the number $c(d)$ of integers $m$ with $\varphi(m)=d$ whose prime factors are among the candidates processed so far; at the start, $c(1)=1$ and all other counts are $0$. Adding the candidate $p$ replaces every count by

$$c'(d)=c(d)+\sum_{e\ge1}c\Bigl(\frac{d}{(p-1)\,p^{e-1}}\Bigr),$$

where a term counts only when the fraction is an integer. The term for $e$ appends $p^e$ to a product formed before $p$ was available, so at most one power of $p$ is chosen, with any allowed exponent. For $T=12$, processing $2$ first gives $c'(1)=2$ for the products $1$ and $2$, and $c'(2)=c'(4)=1$ for $4$ and $8$; after all five candidates the counts are those of Section 6. For $13!$, the count at $13!$ is $182\,752$, the number listed in OEIS A055506, ^[5] and the counts of all divisors add up to $2\,795\,333$, the number of states of Section 6. Storing the smallest product instead of a count, and taking the minimum of the old entry and the values $p^e$ times the existing old entries at $d/((p-1)p^{e-1})$ instead of the sum, gives the least solution at $T$, listed in OEIS A055487; ^[6] storing lists of products gives every solution, which we use below as an independent check.

How the Code Works

The C++ program ^[9] first reads an optional first command-line argument as the index; only a positive value replaces the default $150\,000$. Then validate_small_case repeats the whole method for $T=120$ and compares the result with a brute force in which phi_u64 computes $\varphi(n)$ by trial division for every $n\le200\,000$, far more than needed: the $9$ candidates $2,3,5,7,11,13,31,41,61$ give the bound $n\le673$ of Section 7, and the $17$ solutions run from $143$ to $462$.

Next, the program builds the candidates of Section 2: factorize_u64 factors $13!$ by trial division, build_divisors_dfs and divisors_of produce the $1584$ divisors in increasing order, and is_prime_u64 tests every $d+1$; the closing sort and std::unique change nothing. The primality test divides by the twelve primes up to $37$ and then runs the strong probable prime test of Miller and Rabin: with $n-1=2^t u$ and $u$ odd, a base $a$ passes if $a^u\equiv1\pmod n$ or $a^{2^v u}\equiv-1\pmod n$ for some $0\le v\lt t$. ^[8] It uses the seven bases 2, 325, 9375, 28178, 450775, 9780504 and 1795265022, which make the test exact for every number below $2^{64}$ when a base that is a multiple of $n$ is skipped, as the code does. The numbers tested here are at most $13!+1\lt2^{33}$ but can exceed $2^{32}$, so the product of two residues can exceed $2^{64}$, and mul_mod_u64 multiplies in unsigned __int128. ^[10]

enumerate_inverse_totient runs the search of Section 3 in dfs_inverse_totient, whose arguments rem_phi, start_index and current_n are the state $(r,i,m)$. At $r=1$ it appends $m$ to the output vector. Otherwise it loops over the candidates from index $i$ on, skips those with $p-1\nmid r$, and runs the exponent loop of Section 3 with phi_piece for $(p-1)p^{e-1}$ and n_piece for $p^e$. The test phi_piece > rem_phi / p ends that loop as soon as the next piece would exceed $r$, so phi_piece *= p cannot overflow. In the default run, the search starts not at the root but at its $479$ children, one for each prime power whose piece divides $13!$: the $459$ candidates and $20$ higher powers of the six primes up to $13$. build_seed_tasks lists them as SeedTask records, the recursion runs below each, and the partial lists are concatenated. With $1$ as its third argument, it runs one undivided search from the root instead; there, the special line that adds the solution $1$ for the target $1$ only repeats the $1$ recorded at the root, and $2$ is missed (Section 5).

normalize_solutions sorts the vector and removes duplicates, of which there are none by Section 5. validate_main_problem checks that the list is not empty and starts with the value $6\,227\,180\,929$ of the statement. It recomputes $\varphi(n)$ for every solution with phi_from_candidate_factorization, which divides $n$ by the candidates, multiplies the pieces and returns $0$ if a factor outside the list remains. After a search that started from the children of the root, it runs the undivided search once more and requires the same sorted list. Finally, main rejects an index beyond the end of the list and prints the entry with to_string_u128.

The solutions are stored as unsigned __int128, although Section 7 shows that $36$ bits suffice; $r$ and the pieces never exceed $13!\lt2^{33}$ and are kept in std::uint64_t. The Python function solve contains the same steps as nested functions. Its is_prime handles numbers below $4$ and multiples of $2$ and $3$ first and then uses the built-in pow with the same seven bases, candidate_primes collects the candidates through a set, and dfs is the same recursion, with an unused variable out_len_before. The list is sorted, made duplicate-free with sorted(set(...)), and the entry is returned as a string; Python integers have unlimited precision. ^[10] The Java version has the same structure: modMul multiplies through BigInteger, dfsInverseTotient keeps $m$ as a BigInteger, and solve sorts with Collections.sort, removes duplicates and returns "0" for an index beyond the end of the list. Neither port runs the checks of validate_small_case or validate_main_problem.

Complexity and Verification

The preparation of Section 2 takes less than a millisecond in C++. With $C$ candidates, $S$ states and $M$ solutions, the loops over the candidates make at most $SC$ divisibility tests, the exponent loops fewer than $2S$, and the sort takes $O(M\log M)$ time. Here $S=2\,795\,333$ (Section 6) and $SC\approx1.28\cdot10^9$: each state with $r\gt1$ scans all remaining candidates, including those with $p-1\gt r$, which can never fit, and these $983\,885\,142$ tests dominate the running time. Since the candidates are sorted, leaving the loop at the first candidate with $p-1\gt r$ would cut them to $41\,926\,639$, about $23$ times fewer; the published code does not do this. The memory is $O(M)$ for the list of $M=182\,752$ solutions, plus small tables and a recursion of at most ten calls. On an Apple M4 Max, a full C++ run with all its checks, compiled with -O2, takes about $0.9$ s; Java takes about $1.3$ s and Python about $27$ s, with a peak memory of about $16$ MB in C++ and $38$ MB in Python.

Besides the checks of the C++ program described above, we verified for this article that:

the C++, Python and Java programs print the same number;
the functions of the C++ file, compiled unchanged into a separate test program, agree for every target $T$ from $2$ to $20\,000$ with a brute force that computes $\varphi(n)$ by a sieve for all $n$ up to the bound of Section 7, which stays below $140\,000$; for $T=1$ the undivided search indeed returns only $1$;
the dynamic programming of Section 8, written separately in Python, reproduces exactly the same $182\,752$ numbers, and its counts and least solutions for $1!,2!,\dots,13!$ agree with OEIS A055506 and A055487; ^[5] ^[6]
the counts in Sections 4, 6 and 7 come from an instrumented copy of the search.

These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts used above; the description of the search tree and its counts, the bounds of Section 7 and the argument for the smallest solution are derived in this article.

Project Euler 248 — Euler's Totient Function Equals 13!. The official statement says that the first number n with φ(n) = 13! is 6227180929 and asks for the 150 000th such number.
Wikipedia — Euler's totient function. The article gives φ(pᵏ) = pᵏ⁻¹(p − 1), proves φ(mn) = φ(m)φ(n) for coprime m and n with the Chinese remainder theorem, and states Euler's product formula. It notes that φ(n) is even for n ≥ 3, that 2ʳ divides φ(n) when n has r distinct odd prime factors, and that φ(2m) = φ(m) for odd m.
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 is a product of primes in exactly one way up to the order of the factors, written canonically with increasing primes and positive exponents.
Max A. Alekseyev — Computing the inverses, their power sums, and extrema for Euler's totient and other multiplicative functions (arXiv:1401.6054). Journal of Integer Sequences 19 (2016), Article 16.5.2. Section 5.1 notes that φ(pᵉ) | n implies that p − 1 divides n and e ≤ νₚ(n) + 1, and finds the possible primes by testing d + 1 for every divisor d of n; for n = 12 it obtains the primes 2, 3, 5, 7, 13 and the inverses 13, 21, 26, 28, 36, 42. Section 4 multiplies one series per prime over the divisors of n, which gives the number or the extrema of the inverses without listing them.
OEIS A055506 — Number of solutions to the equation phi(x) = n!. The sequence begins 2, 3, 4, 10, 17, 49, 93, 359, 1138, 3802, 12124, 52844, 182752; its 13th term counts the solutions for 13!.
OEIS A055487 — Least m such that phi(m) = n!. Its 13th term is 6227180929, the first solution given in the statement.
Wikipedia — Backtracking. Backtracking builds solutions step by step, abandons a partial solution as soon as it cannot be completed, and traverses the resulting search tree recursively, from the root down, in depth-first order.
Wikipedia — Miller–Rabin primality test. With n − 1 = 2ˢ·d and d odd, n is a strong probable prime to the base a if a^d ≡ 1 (mod n) or a^(2ʳ·d) ≡ −1 (mod n) for some 0 ≤ r < s, and a suitable fixed set of bases makes the test deterministic below 2⁶⁴. The table of best known SPRP base sets lists the seven bases 2, 325, 9375, 28178, 450775, 9780504 and 1795265022, found by Jim Sinclair in 2011, for all n below 2⁶⁴, and remarks that a base that reduces to 0 modulo n must count as a pass.
C++ — Euler248.cpp. The linked immutable revision contains dfs_inverse_totient, build_seed_tasks, enumerate_inverse_totient, validate_small_case and validate_main_problem; it is the implementation record for this article.
GCC manual — 128-bit Integers: unsigned __int128 is an extension on targets with an integer mode wide enough for 128 bits. Python — numeric types: integers have unlimited precision. Java — BigInteger: immutable arbitrary-precision integers.

Problem 248 source code

C++

#include <algorithm>
#include <atomic>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <limits>
#include <thread>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr u64 kTargetPhi = 6'227'020'800ULL;  // 13!
constexpr u64 kKnownFirstSolution = 6'227'180'929ULL;
constexpr std::size_t kDefaultIndex = 150'000;

struct SeedTask {
    u64 rem_phi = 1;
    std::size_t next_index = 0;
    u128 current_n = 1;
};

u64 pow_u64(u64 base, int exp) {
    u64 out = 1;
    for (int i = 0; i < exp; ++i) out *= base;
    return out;
}

u64 mul_mod_u64(u64 a, u64 b, u64 mod) {
    return static_cast<u64>((static_cast<u128>(a) * static_cast<u128>(b)) % mod);
}

u64 pow_mod_u64(u64 base, u64 exp, u64 mod) {
    u64 out = 1 % mod;
    u64 cur = base % mod;
    while (exp > 0) {
        if (exp & 1ULL) out = mul_mod_u64(out, cur, mod);
        cur = mul_mod_u64(cur, cur, mod);
        exp >>= 1ULL;
    }
    return out;
}

bool is_prime_u64(u64 n) {
    if (n < 2) return false;
    for (u64 p : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
        if (n == p) return true;
        if (n % p == 0) return false;
    }

    u64 d = n - 1;
    int s = 0;
    while ((d & 1ULL) == 0) {
        d >>= 1ULL;
        ++s;
    }

    // Deterministic bases for 64-bit integers.
    constexpr u64 kBases[] = {
        2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL,
    };

    for (u64 a : kBases) {
        if (a % n == 0) continue;
        u64 x = pow_mod_u64(a, d, n);
        if (x == 1 || x == n - 1) continue;

        bool witness = true;
        for (int r = 1; r < s; ++r) {
            x = mul_mod_u64(x, x, n);
            if (x == n - 1) {
                witness = false;
                break;
            }
        }
        if (witness) return false;
    }

    return true;
}

std::vector<std::pair<u64, int>> factorize_u64(u64 n) {
    std::vector<std::pair<u64, int>> fac;
    for (u64 p = 2; p * p <= n; ++p) {
        if (n % p != 0) continue;
        int e = 0;
        while (n % p == 0) {
            n /= p;
            ++e;
        }
        fac.push_back({p, e});
    }
    if (n > 1) fac.push_back({n, 1});
    return fac;
}

void build_divisors_dfs(const std::vector<std::pair<u64, int>>& fac,
                        std::size_t idx,
                        u64 cur,
                        std::vector<u64>& out) {
    if (idx == fac.size()) {
        out.push_back(cur);
        return;
    }

    const u64 p = fac[idx].first;
    const int e = fac[idx].second;
    u64 value = 1;
    for (int i = 0; i <= e; ++i) {
        build_divisors_dfs(fac, idx + 1, cur * value, out);
        value *= p;
    }
}

std::vector<u64> divisors_of(u64 n) {
    const auto fac = factorize_u64(n);
    std::vector<u64> out;
    out.reserve(2048);
    build_divisors_dfs(fac, 0, 1, out);
    std::sort(out.begin(), out.end());
    return out;
}

std::vector<u64> candidate_primes_for_inverse_totient(u64 phi_target) {
    std::vector<u64> candidates;
    const std::vector<u64> divs = divisors_of(phi_target);
    candidates.reserve(divs.size());

    for (u64 d : divs) {
        const u64 p = d + 1;
        if (is_prime_u64(p)) candidates.push_back(p);
    }

    std::sort(candidates.begin(), candidates.end());
    candidates.erase(std::unique(candidates.begin(), candidates.end()), candidates.end());
    return candidates;
}

void dfs_inverse_totient(u64 rem_phi,
                         std::size_t start_index,
                         u128 current_n,
                         const std::vector<u64>& candidates,
                         std::vector<u128>& out) {
    if (rem_phi == 1) {
        out.push_back(current_n);
        return;
    }

    for (std::size_t i = start_index; i < candidates.size(); ++i) {
        const u64 p = candidates[i];
        u64 phi_piece = p - 1;
        if (rem_phi % phi_piece != 0) continue;

        u128 n_piece = static_cast<u128>(p);
        while (rem_phi % phi_piece == 0) {
            dfs_inverse_totient(rem_phi / phi_piece, i + 1, current_n * n_piece, candidates, out);

            if (phi_piece > rem_phi / p) break;
            phi_piece *= p;
            n_piece *= static_cast<u128>(p);
        }
    }
}

std::vector<SeedTask> build_seed_tasks(u64 phi_target, const std::vector<u64>& candidates) {
    std::vector<SeedTask> tasks;
    tasks.reserve(candidates.size() * 2);

    for (std::size_t i = 0; i < candidates.size(); ++i) {
        const u64 p = candidates[i];
        u64 phi_piece = p - 1;
        u128 n_piece = static_cast<u128>(p);

        while (phi_target % phi_piece == 0) {
            tasks.push_back(SeedTask{phi_target / phi_piece, i + 1, n_piece});
            if (phi_piece > phi_target / p) break;
            phi_piece *= p;
            n_piece *= static_cast<u128>(p);
        }
    }

    return tasks;
}

std::vector<u128> enumerate_inverse_totient(u64 phi_target,
                                            const std::vector<u64>& candidates,
                                            int threads) {
    if (threads < 1) threads = 1;

    std::vector<u128> out;
    if (threads == 1) {
        dfs_inverse_totient(phi_target, 0, 1, candidates, out);
        if (phi_target == 1) out.push_back(1);
        return out;
    }

    const std::vector<SeedTask> tasks = build_seed_tasks(phi_target, candidates);
    std::atomic<std::size_t> next_task{0};
    std::vector<std::vector<u128>> partial(static_cast<std::size_t>(threads));
    std::vector<std::thread> pool;
    pool.reserve(static_cast<std::size_t>(threads));

    for (int t = 0; t < threads; ++t) {
        pool.emplace_back([&, t]() {
            std::vector<u128>& local = partial[static_cast<std::size_t>(t)];
            while (true) {
                const std::size_t id = next_task.fetch_add(1, std::memory_order_relaxed);
                if (id >= tasks.size()) break;
                const SeedTask& task = tasks[id];
                dfs_inverse_totient(task.rem_phi, task.next_index, task.current_n, candidates, local);
            }
        });
    }

    for (std::thread& th : pool) th.join();

    std::size_t total = (phi_target == 1) ? 1 : 0;
    for (const auto& vec : partial) total += vec.size();
    out.reserve(total);

    if (phi_target == 1) out.push_back(1);
    for (auto& vec : partial) {
        out.insert(out.end(), vec.begin(), vec.end());
    }
    return out;
}

void normalize_solutions(std::vector<u128>& solutions) {
    std::sort(solutions.begin(), solutions.end());
    solutions.erase(std::unique(solutions.begin(), solutions.end()), solutions.end());
}

u64 phi_u64(u64 n) {
    if (n == 0) return 0;
    u64 result = n;

    for (u64 p = 2; p * p <= n; ++p) {
        if (n % p != 0) continue;
        while (n % p == 0) n /= p;
        result = result / p * (p - 1);
    }
    if (n > 1) result = result / n * (n - 1);
    return result;
}

u64 phi_from_candidate_factorization(u128 n, const std::vector<u64>& candidates) {
    u128 result = 1;
    u128 rem = n;

    for (u64 p : candidates) {
        if (rem % p != 0) continue;

        int exp = 0;
        while (rem % p == 0) {
            rem /= p;
            ++exp;
        }

        u128 term = static_cast<u128>(p - 1);
        for (int i = 1; i < exp; ++i) term *= static_cast<u128>(p);
        result *= term;
    }

    if (rem != 1) return 0;
    if (result > std::numeric_limits<u64>::max()) return 0;
    return static_cast<u64>(result);
}

std::string to_string_u128(u128 value) {
    if (value == 0) return "0";
    std::string out;
    while (value > 0) {
        const int digit = static_cast<int>(value % 10);
        out.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    std::reverse(out.begin(), out.end());
    return out;
}

bool validate_small_case() {
    constexpr u64 small_target = 120;
    constexpr u64 brute_limit = 200'000;

    const std::vector<u64> candidates = candidate_primes_for_inverse_totient(small_target);
    std::vector<u128> by_solver = enumerate_inverse_totient(small_target, candidates, 1);
    normalize_solutions(by_solver);

    std::vector<u128> by_bruteforce;
    for (u64 n = 1; n <= brute_limit; ++n) {
        if (phi_u64(n) == small_target) by_bruteforce.push_back(static_cast<u128>(n));
    }

    if (by_solver != by_bruteforce) {
        std::cerr << "Validation failed for phi=120 (solver vs brute force mismatch).\n";
        std::cerr << "solver size=" << by_solver.size() << ", brute size=" << by_bruteforce.size() << "\n";
        return false;
    }

    return true;
}

bool validate_main_problem(const std::vector<u64>& candidates,
                           const std::vector<u128>& threaded_solutions,
                           int threads) {
    if (threaded_solutions.empty()) {
        std::cerr << "Validation failed: no solutions found for phi=13!.\n";
        return false;
    }

    if (threaded_solutions.front() != static_cast<u128>(kKnownFirstSolution)) {
        std::cerr << "Validation failed: first solution mismatch. got="
                  << to_string_u128(threaded_solutions.front())
                  << " expected=" << kKnownFirstSolution << "\n";
        return false;
    }

    for (u128 n : threaded_solutions) {
        const u64 phi_value = phi_from_candidate_factorization(n, candidates);
        if (phi_value != kTargetPhi) {
            std::cerr << "Validation failed: phi(n) mismatch for n=" << to_string_u128(n)
                      << ", got phi=" << phi_value << "\n";
            return false;
        }
    }

    if (threads > 1) {
        std::vector<u128> single = enumerate_inverse_totient(kTargetPhi, candidates, 1);
        normalize_solutions(single);
        if (single != threaded_solutions) {
            std::cerr << "Validation failed: single-thread and multi-thread outputs differ.\n";
            std::cerr << "single size=" << single.size() << ", multi size=" << threaded_solutions.size() << "\n";
            return false;
        }
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    std::size_t index = kDefaultIndex;
    int threads = static_cast<int>(std::thread::hardware_concurrency());
    if (threads <= 0) threads = 4;

    if (argc > 1) {
        const long long parsed = std::atoll(argv[1]);
        if (parsed > 0) index = static_cast<std::size_t>(parsed);
    }
    if (argc > 2) {
        const int parsed_threads = std::atoi(argv[2]);
        if (parsed_threads > 0) threads = parsed_threads;
    }

    if (!validate_small_case()) return 1;

    const std::vector<u64> candidates = candidate_primes_for_inverse_totient(kTargetPhi);
    std::vector<u128> solutions = enumerate_inverse_totient(kTargetPhi, candidates, threads);
    normalize_solutions(solutions);

    if (!validate_main_problem(candidates, solutions, threads)) return 1;

    if (index == 0 || index > solutions.size()) {
        std::cerr << "Requested index out of range: " << index
                  << " (number of solutions=" << solutions.size() << ")\n";
        return 1;
    }

    const u128 answer = solutions[index - 1];
    std::cout << to_string_u128(answer) << '\n';
    return 0;
}

Python

def solve():
    TARGET_PHI = 6_227_020_800  # 13!
    INDEX = 150_000

    def is_prime(n):
        if n < 2:
            return False
        if n < 4:
            return True
        if n % 2 == 0 or n % 3 == 0:
            return False
        d = n - 1
        s = 0
        while d % 2 == 0:
            d //= 2
            s += 1
        for a in [2, 325, 9375, 28178, 450775, 9780504, 1795265022]:
            if a % n == 0:
                continue
            x = pow(a, d, n)
            if x == 1 or x == n - 1:
                continue
            witness = True
            for _ in range(s - 1):
                x = (x * x) % n
                if x == n - 1:
                    witness = False
                    break
            if witness:
                return False
        return True

    def factorize(n):
        fac = []
        p = 2
        while p * p <= n:
            if n % p == 0:
                e = 0
                while n % p == 0:
                    n //= p
                    e += 1
                fac.append((p, e))
            p += 1
        if n > 1:
            fac.append((n, 1))
        return fac

    def divisors_of(n):
        fac = factorize(n)
        divs = [1]
        for p, e in fac:
            new_divs = []
            pe = 1
            for _ in range(e + 1):
                for d in divs:
                    new_divs.append(d * pe)
                pe *= p
            divs = new_divs
        divs.sort()
        return divs

    def candidate_primes(phi_target):
        divs = divisors_of(phi_target)
        cands = sorted(set(d + 1 for d in divs if is_prime(d + 1)))
        return cands

    candidates = candidate_primes(TARGET_PHI)

    def dfs(rem_phi, start_index, current_n, out):
        if rem_phi == 1:
            out.append(current_n)
            return
        for i in range(start_index, len(candidates)):
            p = candidates[i]
            phi_piece = p - 1
            if rem_phi % phi_piece != 0:
                continue
            n_piece = p
            while rem_phi % phi_piece == 0:
                out_len_before = len(out)
                dfs(rem_phi // phi_piece, i + 1, current_n * n_piece, out)
                if phi_piece > rem_phi // p:
                    break
                phi_piece *= p
                n_piece *= p

    solutions = []
    dfs(TARGET_PHI, 0, 1, solutions)
    solutions.sort()
    # Remove duplicates
    solutions = sorted(set(solutions))

    return str(solutions[INDEX - 1])

if __name__ == '__main__':
    print(solve())

Java

import java.math.BigInteger;
import java.util.*;

public class Euler248 {
    static final long TARGET_PHI = 6227020800L;
    static final int DEFAULT_INDEX = 150000;

    static long modMul(long a, long b, long mod) {
        return BigInteger.valueOf(a).multiply(BigInteger.valueOf(b)).mod(BigInteger.valueOf(mod)).longValue();
    }

    static long modPow(long base, long exp, long mod) {
        long out = 1 % mod;
        long cur = base % mod;
        while (exp > 0) {
            if ((exp & 1) != 0)
                out = modMul(out, cur, mod);
            cur = modMul(cur, cur, mod);
            exp >>= 1;
        }
        return out;
    }

    static boolean isPrime(long n) {
        if (n < 2)
            return false;
        long[] smallPrimes = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 };
        for (long p : smallPrimes) {
            if (n == p)
                return true;
            if (n % p == 0)
                return false;
        }

        long d = n - 1;
        int s = 0;
        while ((d & 1) == 0) {
            d >>= 1;
            ++s;
        }

        long[] bases = { 2L, 325L, 9375L, 28178L, 450775L, 9780504L, 1795265022L };

        for (long a : bases) {
            if (a % n == 0)
                continue;
            long x = modPow(a, d, n);
            if (x == 1 || x == n - 1)
                continue;

            boolean witness = true;
            for (int r = 1; r < s; ++r) {
                x = modMul(x, x, n);
                if (x == n - 1) {
                    witness = false;
                    break;
                }
            }
            if (witness)
                return false;
        }
        return true;
    }

    static class Factor {
        long p;
        int e;

        Factor(long p, int e) {
            this.p = p;
            this.e = e;
        }
    }

    static List<Factor> factorize(long n) {
        List<Factor> fac = new ArrayList<>();
        for (long p = 2; p * p <= n; ++p) {
            if (n % p != 0)
                continue;
            int e = 0;
            while (n % p == 0) {
                n /= p;
                ++e;
            }
            fac.add(new Factor(p, e));
        }
        if (n > 1)
            fac.add(new Factor(n, 1));
        return fac;
    }

    static void buildDivisorsDFS(List<Factor> fac, int idx, long cur, List<Long> out) {
        if (idx == fac.size()) {
            out.add(cur);
            return;
        }
        long p = fac.get(idx).p;
        int e = fac.get(idx).e;
        long value = 1;
        for (int i = 0; i <= e; ++i) {
            buildDivisorsDFS(fac, idx + 1, cur * value, out);
            value *= p;
        }
    }

    static List<Long> divisorsOf(long n) {
        List<Factor> fac = factorize(n);
        List<Long> out = new ArrayList<>();
        buildDivisorsDFS(fac, 0, 1, out);
        Collections.sort(out);
        return out;
    }

    static List<Long> candidatePrimes(long phiTarget) {
        List<Long> divs = divisorsOf(phiTarget);
        List<Long> candidates = new ArrayList<>();
        for (long d : divs) {
            long p = d + 1;
            if (isPrime(p)) {
                candidates.add(p);
            }
        }
        Collections.sort(candidates);
        List<Long> uniqueCandidates = new ArrayList<>();
        for (Long c : candidates) {
            if (uniqueCandidates.isEmpty() || !uniqueCandidates.get(uniqueCandidates.size() - 1).equals(c)) {
                uniqueCandidates.add(c);
            }
        }
        return uniqueCandidates;
    }

    static void dfsInverseTotient(long remPhi, int startIdx, BigInteger currentN, List<Long> candidates,
            List<BigInteger> out) {
        if (remPhi == 1) {
            out.add(currentN);
            return;
        }

        for (int i = startIdx; i < candidates.size(); ++i) {
            long p = candidates.get(i);
            long phiPiece = p - 1;
            if (remPhi % phiPiece != 0)
                continue;

            BigInteger nPiece = BigInteger.valueOf(p);
            while (remPhi % phiPiece == 0) {
                dfsInverseTotient(remPhi / phiPiece, i + 1, currentN.multiply(nPiece), candidates, out);

                if (phiPiece > remPhi / p)
                    break;
                phiPiece *= p;
                nPiece = nPiece.multiply(BigInteger.valueOf(p));
            }
        }
    }

    public static String solve() {
        List<Long> candidates = candidatePrimes(TARGET_PHI);
        List<BigInteger> out = new ArrayList<>();

        dfsInverseTotient(TARGET_PHI, 0, BigInteger.ONE, candidates, out);

        Collections.sort(out);
        List<BigInteger> uniqueOut = new ArrayList<>();
        for (BigInteger c : out) {
            if (uniqueOut.isEmpty() || !uniqueOut.get(uniqueOut.size() - 1).equals(c)) {
                uniqueOut.add(c);
            }
        }

        if (DEFAULT_INDEX == 0 || DEFAULT_INDEX > uniqueOut.size())
            return "0";
        return uniqueOut.get(DEFAULT_INDEX - 1).toString();
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}