Problem 124: Ordered Radicals

Project Euler Problem 124 Solution

EulerSolve provides an optimized solution for Project Euler Problem 124, Ordered Radicals, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each positive integer $n$, define $$\operatorname{rad}(n)=\prod_{p\mid n} p,$$ the product of the distinct prime divisors of $n$, with $\operatorname{rad}(1)=1$. The problem orders all integers $1 \le n \le 100{,}000$ by the pair $(\operatorname{rad}(n),n)$: smaller radical first, and if two numbers have the same radical, the smaller integer comes first. If $E(k)$ denotes the $k$-th term of that ordered list, the goal is to determine $E(10000)$. The only real computational issue is that the order depends on all $100{,}000$ radical values. The implementations avoid factoring each number separately. Instead, they compute every radical in one sieve-like pass and then perform a single sort using the exact order required by the statement. Mathematical Approach The structure of the problem is simple but very specific: compute $\operatorname{rad}(n)$ for every $n$, sort the resulting pairs, and read one position. The mathematical content lies in finding the right way to build the whole radical table efficiently. The ordering that defines $E(k)$ Write $$a \prec b \quad\Longleftrightarrow\quad \bigl(\operatorname{rad}(a),a\bigr) \lt \bigl(\operatorname{rad}(b),b\bigr)$$ in lexicographic order. Then $E(1),E(2),\dots,E(100000)$ is just the list of integers $1,2,\dots,100000$ sorted by $\prec$....

Detailed mathematical approach

Problem Summary

The radical $\operatorname{rad}(n)$ of a positive integer $n$ is the product of the distinct prime factors of $n$; for example $504=2^3\cdot3^2\cdot7$, so $\operatorname{rad}(504)=2\cdot3\cdot7=42$. The numbers $1\le n\le N$ are sorted by their radicals, and numbers with the same radical by their own size; $E(k)$ is the $k$-th number of the sorted list. For $N=10$ the statement shows the sorted list $1,2,4,8,3,9,5,6,7,10$, so $E(4)=8$ and $E(6)=9$. We must find $E(10000)$ for $N=100000$. ^[1]

The whole solution in one sentence: compute the radicals of all numbers up to $100000$ in one sweep that multiplies every prime $p$ into all multiples of $p$, sort the numbers by the pair $(\operatorname{rad}(n),n)$, and read the $10000$-th entry of the sorted list.

Mathematical Approach

We need three things: what the radical does to a number, the exact order with its tie-break, and all $100000$ radicals without factoring each number on its own.

1. The radical forgets exponents

By the fundamental theorem of arithmetic, every $n\gt1$ can be written in exactly one way, up to the order of the factors, as a product $n=p_1^{e_1}\cdots p_s^{e_s}$ of powers of distinct primes with exponents $e_i\ge1$; for $n=1$ the product is empty and equals $1$. ^[3] The radical keeps the primes and drops the exponents:

$$\operatorname{rad}(n)=p_1p_2\cdots p_s,\qquad\operatorname{rad}(1)=1.$$

So $2$, $4$, $8$ and $16$ all have the radical $2$, and $12=2^2\cdot3$, $18=2\cdot3^2$ and $72=2^3\cdot3^2$ all have the radical $6$. Three facts follow at once. The radical divides $n$, so $\operatorname{rad}(n)\le n$. Equality holds exactly when every exponent is $1$, that is, when $n$ is square-free: divisible by no square other than $1$. And $\operatorname{rad}(n)$ itself is square-free; it is even the largest square-free divisor of $n$, because a square-free divisor of $n$ is a product of some of the primes $p_i$, each taken once, and so divides $\operatorname{rad}(n)$. ^[2] ^[4]

2. The order and its blocks

The statement compares two numbers by their radicals and, only when the radicals are equal, by the numbers themselves. This is the lexicographic order of the pairs $(\operatorname{rad}(n),n)$: ^[5]

$$(r,a)\lt(r',b)\iff r\lt r'\ \lor\ \bigl(r=r'\land a\lt b\bigr).$$

Two different numbers never have the same pair, so the sorted list is unique and every correct sorting method produces it. It consists of blocks of numbers with equal radical, in increasing order of the radical, and each block is in increasing order. For $N=10$ the radicals of $1,\dots,10$ are $1,2,3,2,5,6,7,2,3,10$, and the blocks are $\{1\}$, $\{2,4,8\}$, $\{3,9\}$, $\{5\}$, $\{6\}$, $\{7\}$ and $\{10\}$. So $E(4)=8$ is the last number with radical $2$, and $E(6)=9$ the last with radical $3$.

By Section 1, $\operatorname{rad}(n)=r$ holds exactly when $r$ is square-free and $n=r\,t$ with a number $t$ whose prime factors all divide $r$, because then $n$ and $r$ have the same prime factors. For $N=100$ the block of $6$ is $6,12,18,24,36,48,54,72,96$, from $t=1,2,3,4,6,8,9,12,16$. The radicals that occur are exactly the square-free $r\le N$: every radical is square-free and at most $N$, and each such $r$ lies in its own block. So the $100000$ numbers form $60794$ blocks. ^[4]

The blocks, and with them $E(k)$, depend on $N$. For $N=100000$ the block of $2$ holds the $16$ powers $2,4,\dots,65536$, so the list begins $1,2,4,8,16,32,\dots$ and $E(6)=32$, not $9$: the order has to be built for $N=100000$ itself.

3. A sieve for all radicals

One sweep gives all radicals together. Start with $R[n]=1$ for every $n\le N$ and walk through $p=2,3,\dots,N$. If $p$ is unmarked, it is prime: multiply $R[m]$ by $p$ for every multiple $m=p,2p,3p,\dots\le N$, and mark the multiples larger than $p$. This is the sieve of Eratosthenes with one change. ^[6] The usual sieve may start marking at $p^2$, since the smaller multiples are already marked; here every multiple of $p$, starting with $p$ itself, needs its factor $p$.

For $N=12$: the pass of $2$ sets $R[m]=2$ for the six even $m$ and marks $4,6,8,10,12$; the pass of $3$ turns $R[3]$, $R[6]$, $R[9]$, $R[12]$ into $3,6,3,6$ and marks $6,9,12$. The marked numbers $4,6,8,9,10,12$ are skipped, and the passes of $5$, $7$ and $11$ give $R[5]=5$, $R[10]=10$, $R[7]=7$ and $R[11]=11$. The final table

$$R[1],\dots,R[12]=1,2,3,2,5,6,7,2,3,10,11,6$$

holds exactly the radicals. The number $8=2^3$ receives the factor $2$ only once, because the pass of $2$ visits $8$ once: this is how the exponents are forgotten.

4. Why the sieve is correct

Claim. When the walk reaches $p$, the number $p$ is unmarked exactly when it is prime, and at the end $R[m]=\operatorname{rad}(m)$ for every $m\le N$.

Proof. We use induction along the walk. A pass is made only for an unmarked number, which by induction is a prime $q$, and it marks only multiples of $q$ larger than $q$. A prime $p$ has no prime divisor $q\lt p$, so it is never marked. A composite $c$ has a prime divisor $q\lt c$, whose pass was made earlier and marked $c$. So exactly the primes are processed.

Next, once the walk has passed $x$, $R[m]$ is the product of the prime divisors of $m$ that are at most $x$. At the start, $x=1$, this is the empty product $1$. The pass of a prime $p$ visits each multiple of $p$ exactly once, stepping through $p,2p,3p,\dots$, and no other entry, so each multiple gains its missing factor $p$, whatever the exponent of $p$ in $m$ is; a composite $x$ changes nothing. At $x=N$ all prime divisors are included, so $R[m]=\operatorname{rad}(m)$. $\square$

The proof also shows that $R[m]$ always divides $\operatorname{rad}(m)$, so $R[m]\le m\le N$ at every moment: no product exceeds $N$, and for the required $N=100000$ the radicals fit easily into 32-bit integers.

5. Sorting and reading the entry

With all radicals known, sort the pairs and read position $k$; arrays count from $0$, so $E(k)$ has index $k-1$. A sort is called stable if items with equal keys keep their original order. ^[9] Since the list $1,2,\dots,N$ starts in increasing order, a stable sort on the radical alone gives the same result as a sort on the pair; for $N=100000$ we checked this. A sort that is not stable, such as std::sort in C++, ^[8] needs the full pair.

6. An alternative without a full sort

Only one position is needed, so counting can replace the sort. Put the numbers $1,\dots,N$ into one bucket per radical; filled in increasing order, every bucket is a sorted block. Then walk through the blocks for $r=1,2,3,\dots$: while $k$ exceeds the size of the current block, subtract that size from $k$ and go on; otherwise the answer is the $k$-th member of the current block. For $N=100$ the first blocks are

$r=1$: $1$; $r=2$: $2,4,8,16,32,64$; $r=3$: $3,9,27,81$; $r=5$: $5,25$;
$r=6$: $6,12,18,24,36,48,54,72,96$; $r=7$: $7,49$; $r=10$: $10,20,40,50,80,100$,

with running totals $1,7,11,13,22,24,30$. The first six blocks hold $24$ numbers, so $E(25)=10$, the first member of the block of $10$; likewise $E(30)=100$ and $E(31)=11$. Every number enters one bucket once and at most $N$ buckets are visited, so after the sieve this costs $O(N)$ time and memory; for $N=100000$ a sort is fast enough, so the published programs sort.

How the Code Works

The C++ program ^[11] starts in main with limit $=100000$, index $=10000$ and the checkpoint on. parse_arguments accepts --limit=, --index= and --skip-checkpoints; an unknown argument, or values outside $1\le k\le N$, give exit code $1$. run_checkpoints checks that solve(10, 4) returns the statement's $E(4)=8$, else the exit code is $2$; then main prints solve(limit, index). Two overflows of the signed 32-bit int, undefined in C++, affect only huge arguments: parse_int_after_prefix reads digits without a range check, so --limit=3000000000 overflows, and the accepted --limit=2147483647 overflows in limit + 1.

solve implements Sections 3 to 5. rad is a std::vector<int> of ones and is_prime a std::vector<char> of true values, each with limit + 1 entries. The outer loop over p skips marked numbers; the inner loop over m = p, 2p, … multiplies rad[m] by p and marks m when m > p. The vector values receives the pairs (rad[n], n), std::sort orders it with a lambda that compares the radicals first and the numbers second, as in Section 2, and values[index - 1].second is returned.

The Python program runs the same sieve on lists and sorts the numbers with sorted(range(1, limit + 1), key=lambda n: (rad[n], n)); tuples are compared lexicographically, ^[9] which is the order of Section 2, and values[9999] is printed. The Java program keeps int[] rad and boolean[] isPrime; since Arrays.sort takes a comparator only for arrays of objects, idx holds the numbers as Integer objects, sorted with rad[a] != rad[b] ? rad[a] - rad[b] : a - b. A comparator returns a negative number, zero or a positive number, ^[10] and these differences cannot overflow, because all values lie between $1$ and $100000$. Both sorts are stable, ^[9] ^[10] so there the tie-break is redundant but harmless. Python and Java fix the limit and the index and run no check.

Complexity and Verification

The pass of a prime $p$ visits its $\lfloor N/p\rfloor$ multiples up to $N$, where $\lfloor x\rfloor$ is $x$ rounded down, so besides the $N-1$ steps of the outer loop the sieve makes

$$\begin{aligned}\sum_{p\le N}\Bigl\lfloor\frac Np\Bigr\rfloor&\le N\sum_{p\le N}\frac1p\\&=N\bigl(\ln\ln N+M+o(1)\bigr)\end{aligned}$$

steps, where $p$ runs over primes, $o(1)$ stands for a term that tends to $0$ as $N$ grows, and $M\approx0.2615$ is the Meissel–Mertens constant of Mertens' second theorem; ^[7] the time is $O(N\log\log N)$. For $N=100000$ there are $9592$ primes up to $N$, and the sieve makes $266400$ multiplications, close to $N(\ln\ln N+M)\approx270500$. Sorting takes $O(N\log N)$ comparisons, ^[8] ^[10] with $N\log_2N\approx1.7\cdot10^6$ here, so it dominates. The memory is $O(N)$: the radicals, the marks, the list of $N$ pairs or numbers and the workspace of the sort. In tests on a desktop computer, each program finished in less than a tenth of a second.

The C++ program checks $E(4)=8$ before the main run, but not $E(6)=9$. For this article all three programs were run and print the same number, and the C++ program with --limit=10 and the indices $1$ to $10$ reproduces the sorted column of the statement. A separate Python check confirms that the radicals of the sieve agree with trial division for every $n\le100000$ and its prime marks for $2\le n\le100000$, and that the counting method of Section 6 agrees with the full sort for every $k$ when $N\le300$ and when $N=1000$. With a detector for undefined behaviour, the default C++ run is clean and both overflows above are reported.

Footnotes and References

The references below supply the background facts; the description of the blocks, the correctness proof of the sieve and the operation counts are derived here.

Project Euler 124 — Ordered Radicals. The official statement defines rad(n), sorts 1 ≤ n ≤ 10 by rad(n) and then by n, gives E(4) = 8 and E(6) = 9, and asks for E(10000) for 1 ≤ n ≤ 100000.
Wikipedia — Radical of an integer. The radical of n is the product of the distinct primes dividing n, for example rad(504) = 42, and it is the largest square-free divisor of n.
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; with the empty product equal to 1, this covers 1 as well.
Wikipedia — Square-free integer. An integer is square-free when no prime exponent in its factorisation exceeds one, exactly when it equals its radical; there are Q(10⁵) = 60794 square-free numbers up to 10⁵.
Wikipedia — Lexicographic order. On pairs, (a, b) ≤ (a′, b′) holds exactly when a < a′, or a = a′ and b ≤ b′; for two totally ordered sets the result is a total order.
Wikipedia — Sieve of Eratosthenes. Marking the multiples of each prime as composite leaves the primes unmarked; marking may start at p², and the total work, a sum of n/p over the primes, is O(n log log n).
Wikipedia — Meissel–Mertens constant. By Mertens' second theorem, the sum of 1/p over the primes p ≤ n, minus ln ln n, tends to the constant M ≈ 0.2614972128.
cppreference — std::sort. The order of equivalent elements is not guaranteed to be preserved, unlike in std::stable_sort, and the comparator is applied O(N·log N) times.
Python — Sorting Techniques. Sorts are guaranteed to be stable, and tuples are compared lexicographically, which allows sorting by several keys.
Java — Arrays.sort with a comparator. The sort is guaranteed to be stable and is implemented as a merge sort. By Comparator.compare, a comparator returns a negative integer, zero or a positive integer as the first argument is less than, equal to or greater than the second.
C++ — Euler124.cpp. The linked immutable revision, with solve and run_checkpoints, is the implementation record for this article.

Problem 124 source code

C++

#include <algorithm>
#include <iostream>
#include <numeric>
#include <string>
#include <utility>
#include <vector>

namespace {

struct Options {
    int limit = 100000;
    int index = 10000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--limit=", options.limit) ||
            parse_int_after_prefix(arg, "--index=", options.index)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.limit >= 1 && options.index >= 1 && options.index <= options.limit;
}

int solve(const int limit, const int index) {
    std::vector<int> rad(static_cast<std::size_t>(limit + 1), 1);
    std::vector<char> is_prime(static_cast<std::size_t>(limit + 1), true);
    if (limit >= 0) {
        is_prime[0] = false;
    }
    if (limit >= 1) {
        is_prime[1] = false;
    }

    for (int p = 2; p <= limit; ++p) {
        if (!is_prime[static_cast<std::size_t>(p)]) {
            continue;
        }
        for (int m = p; m <= limit; m += p) {
            rad[static_cast<std::size_t>(m)] *= p;
            if (m > p) {
                is_prime[static_cast<std::size_t>(m)] = false;
            }
        }
    }

    std::vector<std::pair<int, int>> values;
    values.reserve(static_cast<std::size_t>(limit));
    for (int n = 1; n <= limit; ++n) {
        values.push_back({rad[static_cast<std::size_t>(n)], n});
    }

    std::sort(values.begin(), values.end(), [](const auto& a, const auto& b) {
        if (a.first != b.first) {
            return a.first < b.first;
        }
        return a.second < b.second;
    });

    return values[static_cast<std::size_t>(index - 1)].second;
}

bool run_checkpoints() {
    if (solve(10, 4) != 8) {
        std::cerr << "Checkpoint failed for limit=10,index=4" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << solve(options.limit, options.index) << '\n';
    return 0;
}

Python

# Problem 124: Ordered radicals
# Sort 1..100000 by rad(n), then by n. Find E(10000).

def solve():
    limit = 100000
    rad = [1] * (limit + 1)
    is_prime = [True] * (limit + 1)
    for p in range(2, limit + 1):
        if not is_prime[p]: continue
        for m in range(p, limit + 1, p):
            rad[m] *= p
            if m > p: is_prime[m] = False
    
    values = sorted(range(1, limit + 1), key=lambda n: (rad[n], n))
    print(values[9999])  # E(10000), 0-indexed

solve()

Java

import java.util.*;

public class Euler124 {
    public static void main(String[] args) {
        int limit = 100000;
        int[] rad = new int[limit + 1];
        Arrays.fill(rad, 1);
        boolean[] isPrime = new boolean[limit + 1];
        Arrays.fill(isPrime, true);
        for (int p = 2; p <= limit; p++) {
            if (!isPrime[p])
                continue;
            for (int m = p; m <= limit; m += p) {
                rad[m] *= p;
                if (m > p)
                    isPrime[m] = false;
            }
        }
        Integer[] idx = new Integer[limit];
        for (int i = 0; i < limit; i++)
            idx[i] = i + 1;
        Arrays.sort(idx, (a, b) -> rad[a] != rad[b] ? rad[a] - rad[b] : a - b);
        System.out.println(idx[9999]);
    }
}