Problem 355: Maximal Coprime Subset

Project Euler Problem 355 Solution

EulerSolve provides an optimized solution for Project Euler Problem 355, Maximal Coprime Subset, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $\operatorname{Co}(n)$ denote the maximum possible sum of a subset $A\subseteq\{1,2,\dots,n\}$ such that any two distinct chosen numbers are coprime: $$\operatorname{Co}(n)=\max\left\{\sum_{x\in A}x : A\subseteq\{1,\dots,n\},\ \gcd(a,b)=1\ \forall a\ne b\in A\right\}.$$ The local C++, Python, and Java solutions all use the same structure: start from the best one-prime-per-slot baseline, then add only profitable replacements of one small-prime power and one large prime by a larger mixed value. Those replacements form a maximum-weight bipartite matching. Mathematical Approach Step 1: One Prime, One Slot If two selected numbers were both divisible by the same prime $p$, their gcd would be at least $p$, so they could not be pairwise coprime. Therefore each prime may appear in at most one chosen number. For a fixed prime $p\le n$, the best number using only that prime is the largest prime power not exceeding $n$: $$P_p=\max\{p^k : k\ge 1,\ p^k\le n\}.$$ This gives the baseline set $$B_0=\{1\}\cup\{P_p : p\in\mathbb{P},\ p\le n\},$$ Here $\mathbb{P}$ denotes the set of prime numbers. which is pairwise coprime because different members are powers of different primes. Its sum is $$S_0=1+\sum_{p\le n}P_p,$$ where the sum runs over primes $p\le n$....

Detailed mathematical approach

Problem Summary

For a positive integer $n$, $\operatorname{Co}(n)$ is the largest possible sum of a set of mutually coprime elements of $\{1,2,\dots,n\}$: any two different chosen numbers must have greatest common divisor $1$. ^[2] The statement gives $\operatorname{Co}(10)=30$, attained on $\{1,5,7,8,9\}$, together with $\operatorname{Co}(30)=193$ and $\operatorname{Co}(100)=1356$, and asks for $\operatorname{Co}(200000)$. ^[1] Since $1$ is coprime to every integer, it always belongs to an optimal set.

The whole solution in one sentence: give every prime its largest power that does not exceed $n$, then let a few of the smaller primes team up with one large prime each — for $n=100$ a number such as $8\cdot11=88$ replaces the two numbers $64$ and $11$ — and decide who teams up with whom by a maximum-weight matching.

Mathematical Approach

1. One prime, one element

By unique factorisation every integer greater than $1$ has a well defined set of prime divisors. ^[3] If a prime $p$ divided two chosen numbers, their greatest common divisor would be at least $p$. Hence each prime divides at most one element, and the prime sets of the elements greater than $1$ are pairwise disjoint; $1$ has no prime divisors at all, which is why it is free.

2. Every prime is used, so the task is a partition

Write $P_p$ for the largest power of the prime $p$ that fits:

$$P_p=\max\{p^k : k\ge1,\ p^k\le n\}.$$

If a prime $p\le n$ divides no element of a coprime set, then $P_p$ is coprime to every element and not already present, so adding it increases the sum. In an optimal set every prime up to $n$ therefore divides exactly one element, and optimal sets correspond to the ways of splitting the primes $p\le n$ into disjoint groups — we call them parts — one part per element. Write $r_S=\prod_{p\in S}p$ for the product of the primes of a part and $\operatorname{rad}(m)$ for the product of the distinct primes dividing $m$. The best element a part $S$ can produce is then

$$v(S)=\max\{m\le n:\ \operatorname{rad}(m)=r_S\},$$

which exists exactly when $r_S\le n$, because $r_S$ itself is then a candidate. The correspondence works both ways: the prime sets of the elements of a coprime set, plus a singleton for every unused prime, form such a partition, and replacing each element by the largest integer with the same prime divisors keeps the set coprime without lowering the sum; conversely the values $v(S_i)$ of any partition, together with $1$, are pairwise coprime and distinct, since different parts have different prime sets. Consequently

$$\operatorname{Co}(n)=1+\max\sum_i v(S_i),$$

the maximum running over all partitions of the primes $p\le n$ into feasible parts. For $n=10$ the primes are $2,3,5,7$; four singleton parts give $8,9,5,7$, which with the $1$ is precisely the $\{1,5,7,8,9\}$ of the statement.

3. The baseline

The partition into singletons is always feasible and gives the baseline

$$B(n)=1+\sum_{p\le n}P_p.$$

It is already optimal for $n=10$, where $B(10)=30$ equals $\operatorname{Co}(10)$. For $n=30$ we get $B(30)=188$, five short of $193$; for $n=100$, $B(100)=1263$, which is $93$ short of $1356$. The baseline is thus the bulk of the sum, and the correction is the whole difficulty.

One inequality is used repeatedly: if $P_p=p^k$ then $p^{k+1}\gt n$, hence

$$P_p\gt n/p.$$

A small prime therefore already owns a large number. For $n=200000$, $P_2=2^{17}=131072$ and $P_3=3^{11}=177147$; their sum alone exceeds $n$, so no single element can ever repay giving the primes $2$ and $3$ to one part.

4. The gain of a part

Measure each part against the baseline: the gain of a part $S$ is

$$g(S)=v(S)-\sum_{p\in S}P_p,$$

so that $g(S)=0$ for a singleton and

$$\operatorname{Co}(n)=B(n)+\max\sum_i g(S_i),$$

the maximum now running over families of disjoint parts; parts with $g(S)\le0$ can be dropped and replaced by singletons. For $n=30$ exactly one part has a positive gain, namely $S=\{2,7\}$ with $v(S)=4\cdot7=28$ and

$$g(S)=28-16-7=5,$$

which gives $\operatorname{Co}(30)=188+5=193$. The optimal set is $\{1,11,13,17,19,23,25,27,28,29\}$: the merged element $28=2^2\cdot7$ replaces $16$ and $7$.

5. Small primes and large primes

Call a prime small when $p\le\lfloor\sqrt n\rfloor$ and large otherwise. Two large primes $q_1\ne q_2$ satisfy $q_1q_2\gt n$, so a part contains at most one large prime. For a large prime $q$ we have $q^2\gt n$ and therefore $P_q=q$: large primes contribute themselves. For a part $\{p,q\}$ with $p$ small and $q$ large, no second factor $q$ fits either, so

$$v(\{p,q\})=q\cdot\max\{p^a : a\ge1,\ p^aq\le n\}.$$

This is the only shape the three programs consider. For $n=200000$ the split is at $\lfloor\sqrt{200000}\rfloor=447$: there are $86$ small primes, the largest being $443$, and $17898$ large ones, the smallest being $449$, out of $17984$ primes in total.

6. A complete example: $n=100$

Here $\lfloor\sqrt{100}\rfloor=10$, the small primes are $2,3,5,7$ with $P_2=64$, $P_3=81$, $P_5=25$, $P_7=49$, and every large prime $11,\dots,97$ contributes itself. Ten parts have a positive gain, all of the predicted shape except one; written $(p,q)\mapsto v$ with the gain in brackets:

$(5,19)\mapsto95\ [51]$, $(5,17)\mapsto85\ [43]$, $(7,13)\mapsto91\ [29]$, $(5,13)\mapsto65\ [27]$;
$(5,11)\mapsto55\ [19]$, $(7,11)\mapsto77\ [17]$, $(2,11)\mapsto88\ [13]$, $(3,11)\mapsto99\ [7]$, $(2,23)\mapsto92\ [5]$;
and the exception $\{2,5\}\mapsto100\ [11]$, a part made of two small primes.

The gains are easy to check by hand: $88=2^3\cdot11$ and $88-64-11=13$; $100=2^2\cdot5^2$ and $100-64-25=11$. The disjoint parts $\{2,11\}$, $\{5,19\}$ and $\{7,13\}$ have total gain $13+51+29=93$, and a price argument proves that no family does better. Give the primes $2,5,7,11$ the prices $6,51,29,7$ and every other prime the price $0$. Each of the ten gains listed above is at most the sum of the prices of its primes — for instance $13\le6+7$ for $\{2,11\}$ and $11\le6+51$ for $\{2,5\}$ — and disjoint parts use each price at most once, so no family can gain more than $6+51+29+7=93$. Our three pairs attain it, so $\operatorname{Co}(100)=1263+93=1356$, as the statement promises. The set itself is

$$\{1,17,23,29,\dots,97\}\cup\{81,88,91,95\},$$

that is $1$, the merged $88$, $95$, $91$, the untouched $P_3=81$ and every unmerged large prime; its $23$ members are pairwise coprime and sum to $1356$.

7. Why the choice is a matching

Each prime lies in at most one part, which is exactly the matching condition in the bipartite graph with the small primes on one side, the large primes on the other and an edge of weight $g(\{p,q\})\gt0$ for every profitable pair. The total gain of a family of disjoint pairs is the weight of the corresponding matching, so the best correction obtainable from such pairs is a maximum-weight matching; Section 8 shows that for $n=200000$ no other kind of part is needed.

Choosing each small prime's best partner independently is not enough, because two small primes may want the same large prime: at $n=88$ both $2$ and $7$ prefer $11$, which is their only profitable partner. Taking the heaviest admissible edge each time can fail too, first at $n=104$: that rule picks $(5,19)\ [51]$, then $(7,13)\ [29]$, then $(2,11)\ [13]$ and gains $93$, while $(5,19)\mapsto95\ [51]$, $(2,13)\mapsto104\ [27]$ and $(7,11)\mapsto77\ [17]$ gain $95$. The greedy choice blocks $13$, the partner that $2$ needs, so $\operatorname{Co}(104)=1467+95=1562$ is missed by $2$.

The same clash occurs at the top of the range for $n=200000$: since $439\cdot457$ and $443\cdot457$ both exceed $n$, each of $439$ and $443$ has exactly one admissible partner, the large prime $449$. The two gains are

$$439\cdot449-439^2-449=3941,$$

$$443\cdot449-443^2-449=2209,$$

so the two pairs cannot both be used. A matching that takes $\{443,449\}$ leaves $439$ unmatched, since $439$ has no other admissible partner, and exchanging the pair for $\{439,449\}$ then gains $3941-2209=1732$ more. Hence every maximum-weight matching leaves $443$ alone with $P_{443}=196249$, and in the optimal matching the prime $449$ does go to $439$. Over the whole range, taking the heaviest admissible edge each time would lose $416$.

8. The parts the programs leave out

The restriction to one small and one large prime is an assumption, not a theorem. For $n=200000$ a complete enumeration of the $121580$ feasible parts — one per squarefree number in $[2,n]$, the largest having six primes, $2\cdot3\cdot5\cdot7\cdot11\cdot13=30030$ — finds $25381$ with positive gain: $19711$ of the expected shape, $311$ pairs of two small primes and $5359$ parts of three or four primes. Two omitted ones are large:

$\{29,83\}$ with $v=29\cdot83^2=199781$ and gain $199781-29^3-83^2=168503$;
$\{23,89,97\}$ with $v=23\cdot89\cdot97=198559$ and gain $169062$.

They still do not pay. No part is worth more than $n$, while a small prime with a well chosen large partner already reaches almost $n$: the optimal family contains $\{83,2399\}$ and $\{29,6883\}$ with gains $189829$ and $168335$, together more than twice the gain of $\{29,83\}$. Merging two small primes spends two of the $86$ scarce small primes on one value below $n$, so it can pay only when a small prime finds no free large partner — an intuition, not a proof, so the restriction was certified instead.

The exact task is a maximum-weight set packing over all $25381$ profitable parts, NP-hard in general ^[7] but small enough here; its optimum consists of $85$ pairs of one small and one large prime and agrees with the three programs. Optimality is proved in exact integer arithmetic by prices, as for $n=100$ in Section 6: a nonnegative integer price $y_p$ per prime with

$$g(S)\le\sum_{p\in S}y_p$$

for every one of the $25381$ parts, with $\sum_p y_p$ equal to the weight of that matching. Only $4654$ primes get a nonzero price, $128$ inequalities are tight, and all $25381$ were checked with integers. Disjoint parts use each price at most once, so no family of any shape can gain more: the restriction loses nothing for $n=200000$. It is a restriction all the same: at $n=2047$, merging $\{2,3\}$ into $2^3\cdot3^5=1944$ beats $1024+729=1753$ by $191$, and only the better pairs $\{2,127\}$ and $\{3,73\}$, gaining $881$ and $1169$, rescue the restricted method there.

9. From matching to the assignment problem

A maximum-weight matching in a bipartite graph is an assignment problem, and the classical algorithms assign every row to a column of its own. ^[5] To make the matching optional, the programs use $s$ rows for the small primes, $\ell$ real columns for the large primes and $s$ further columns of weight $0$; every entry outside a profitable pair is $0$ too. Keeping only the positive entries of a complete assignment leaves a family of disjoint profitable pairs of the same total weight, and conversely any such family extends to one by sending unused rows to distinct zero columns. The two optima therefore agree, and a row ending on a zero entry is a small prime left alone with $P_p$.

How the Code Works

All three files follow the same five steps, driven by solve_co_optimized in C++, solve_n in Python and solveN in Java; Java spells the helpers below sievePrimes, maxPrimePowerLeqN, highestPowerWithLimit and hungarianMaxWeightBipartite. sieve_primes runs a sieve of Eratosthenes up to $n$, crossing out multiples of $i$ from $i^2$ on. ^[4] C++ uses std::vector<bool>, Java a boolean[] and Python a bytearray.

max_prime_power_leq_n computes $P_p$ of Section 2 by multiplying value by prime while value <= n / prime; the division instead of a multiplication keeps the test free of overflow. These values are accumulated into baseline_sum, which starts at $1$ for the element $1$ and ends as $B(n)$ of Section 3.

The primes are then split at split, which is math.isqrt(n) in Python and a cast of sqrt(n) in C++ and Java, into small_primes and large_primes together with their powers, as in Section 5. The C++ version keeps the four lists plus the baseline in one SolverData struct filled by prepare_solver_data.

Next comes the double loop over candidate pairs. For a fixed small prime p it walks the large primes in increasing order and breaks as soon as p * q > n, since no later q can fit either. For the surviving pairs highest_power_with_limit(p, n / q) returns the largest power of p not exceeding n / q, so pair_value is exactly $v(\{p,q\})$ of Section 5. Pairs with pair_value <= p_power + q_power are skipped, and the rest are stored as $g(\{p,q\})$ in a weight matrix of $s$ rows and $\ell+s$ columns, at index si * right_size + li; the last $s$ columns stay zero and are the padding of Section 9. Python and Java write into that matrix directly, while the C++ build_gain_edges first returns the edges as (small index, large index, gain) tuples that solve_co_optimized copies in.

hungarian_max_weight_bipartite then solves the padded instance. It negates the weights, turning the maximum into a minimum-cost assignment, and runs the standard potential-based Hungarian method: u and v hold the row and column potentials, p[j] the row currently matched to column j, minv[j] the smallest reduced cost reaching column j, way[j] the previous column on that path and used[j] the columns already visited. ^[6] For each row the inner loop extends the cheapest path, shifts the potentials by delta so that an edge into the chosen column has zero reduced cost, and stops at the first unmatched column; the walk along way then flips the path and matches the new row. Summing the matched entries and flipping the sign back gives the total gain. The three versions of this function agree line by line: the weights, the potentials and minv are signed 64-bit integers, p and way are plain int, and used is a bit vector, a boolean[] or a bytearray. No value comes near the signed 64-bit limit, since one gain is below $n$ and the total below $86n$, and Python's integers are unbounded anyway.

Only the C++ file adds options and self-checks: parse_arguments accepts --n= to change $n$ and --skip-checkpoints, and run_checkpoints verifies $\operatorname{Co}(10)=30$, $\operatorname{Co}(30)=193$ and $\operatorname{Co}(100)=1356$ before the real computation, aborting with a message if one fails. ^[9]

Complexity and Verification

The sieve costs $O(n\log\log n)$ time and $O(n)$ memory. ^[4] Writing $\pi(x)$ for the number of primes up to $x$, the number of pairs that pass the product test is

$$\sum_{p\le\sqrt n}\bigl(\pi(n/p)-\pi(\sqrt n)\bigr),$$

which is $41489$ for $n=200000$, each costing $O(\log n)$ multiplications; $19711$ survive as positive gains, involving only $4568$ different large primes, the largest $49999$. The weight matrix has $86\cdot17984=1546624$ entries and dominates the memory of the matching stage: about $12.4$ MB of 64-bit integers in C++ and Java, while Python stores references instead. The Hungarian method inserts $s$ rows; one search can cross up to $s$ occupied columns and each step scans all $\ell+s$ columns, giving $O(s^2(\ell+s))$ worst-case time, ^[6] here some $10^8$ operations. In practice it is far quicker, and that was measured rather than assumed: the $86$ insertions took $109$ column passes in all, at most $6$ for one row and a single pass for $76$ of them, about $2\cdot10^6$ entries read. Compiled with -O2 on one laptop, the C++ needs $0.01$ s including the checkpoints, Java $0.07$ s and Python $0.26$ s.

Beyond those checkpoints, the following was verified for this article. An independent brute force — a dynamic program over subsets of primes that tries every number $2,\dots,n$, with no pair restriction — reproduces $\operatorname{Co}(n)$ for all $n\le60$, agrees with the three programs there and matches the $57$ terms of the sequence published in the OEIS, including $\operatorname{Co}(51)=453$. ^[8] The partition model of Section 2, combined with an exact set packing over all profitable parts, agrees with the three programs for every $n\le2400$ and for $195$ further values sampled up to $200000$, among them $200000$ itself, where the prices of Section 8 close the argument. Every gain quoted above was recomputed from the definition of $v(S)$, not taken from the programs, and the optimal sets for $n=10$, $30$, $100$ and $104$ are pairwise coprime and sum to $30$, $193$, $1356$ and $1562$. Finally, the C++, Java and Python programs print the same value for $n=200000$.

Footnotes and References

The references below supply the background facts used above. The reduction to a partition of the primes, the gain of a part, the matching argument and the prices that justify the restriction to small-and-large pairs are derived or computed explicitly for this article.

Project Euler 355 — Maximal Coprime Subset. The official statement defines Co(n) as the maximal possible sum of a set of mutually coprime elements of {1, 2, …, n}, notes that Co(10) = 30 is attained on {1, 5, 7, 8, 9}, gives Co(30) = 193 and Co(100) = 1356, and asks for Co(200000). The three given values are used here as test cases only.
Wikipedia — Coprime integers. Two integers are coprime when their only common positive divisor is 1, equivalently when their greatest common divisor is 1, and a set is pairwise coprime when every two distinct members are coprime. It also states that 1 is coprime to every integer, as used in the summary.
Wikipedia — Fundamental theorem of arithmetic. Every integer greater than 1 is a product of primes, uniquely up to the order of the factors. Sections 1 and 2 use it to speak of the set of prime divisors of an element.
Wikipedia — Sieve of Eratosthenes. The sieve marks the multiples of each prime starting from its square, needs O(n log log n) operations and O(n) memory in its basic array form. It is the first stage of all three programs.
Wikipedia — Assignment problem. The article describes the linear assignment problem, calls an instance with unequal sides unbalanced and illustrates dummy tasks of cost 0. The zero columns of Section 9 are such dummies, there not to balance the sides — the column side is already larger — but to let a row stay unassigned.
cp-algorithms — Hungarian algorithm for the assignment problem. This page presents the implementation the three files follow: the potentials u and v, the array p of the row matched to each column, minv, way and used, the search that shifts the potentials by delta until it meets a free column, and the complexity O(n³), or O(n²m) for an n × m matrix with n ≤ m.
Wikipedia — Set packing. Set packing asks for pairwise disjoint subsets from a given list and is NP-complete, while the special case of two-element subsets is equivalent to a matching problem, solvable in polynomial time. This is why parts of two primes (Section 7) make the problem easy and why the check of Section 8 needs a packing solver.
OEIS A186736 — maximum sum of relatively prime integers no larger than n. The entry lists 57 terms, beginning 1, 3, 6, 8, 13, 13, 20, 24, 30, 30, and gives the optimal sets for n = 30, namely {1, 11, 13, 17, 19, 23, 25, 27, 28, 29}, and for n = 51. The brute force was compared against all of them.
C++ — Euler355.cpp. The linked immutable C++ revision contains prepare_solver_data, build_gain_edges, hungarian_max_weight_bipartite, solve_co_optimized and the three checkpoints of run_checkpoints. It is the implementation record for this article; the enumeration, the prices and the brute force described above are separate checks, not part of that file.

Problem 355 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <tuple>
#include <vector>

namespace {

using u64 = std::uint64_t;
using i64 = std::int64_t;

constexpr u64 kDefaultN = 200'000ULL;

struct Options {
    u64 n = kDefaultN;
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }

    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        return false;
    }

    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(const 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 (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        u64 parsed_u64 = 0ULL;
        if (parse_u64_after_prefix(arg, "--n=", parsed_u64)) {
            options.n = parsed_u64;
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

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

    return true;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const std::size_t workload_units) {
    if (!allow_multithreading || workload_units < 2ULL) {
        return 1U;
    }

    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) {
            threads = 1U;
        }
    }

    return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload_units)));
}

std::vector<int> sieve_primes(const int n) {
    if (n < 2) {
        return {};
    }

    std::vector<bool> is_prime(static_cast<std::size_t>(n) + 1ULL, true);
    is_prime[0] = false;
    is_prime[1] = false;

    for (int i = 2; static_cast<u64>(i) * static_cast<u64>(i) <= static_cast<u64>(n); ++i) {
        if (!is_prime[static_cast<std::size_t>(i)]) {
            continue;
        }
        for (int j = i * i; j <= n; j += i) {
            is_prime[static_cast<std::size_t>(j)] = false;
        }
    }

    std::vector<int> primes;
    for (int i = 2; i <= n; ++i) {
        if (is_prime[static_cast<std::size_t>(i)]) {
            primes.push_back(i);
        }
    }
    return primes;
}

u64 max_prime_power_leq_n(const u64 n, const int prime) {
    u64 value = static_cast<u64>(prime);
    while (value <= n / static_cast<u64>(prime)) {
        value *= static_cast<u64>(prime);
    }
    return value;
}

u64 highest_power_with_limit(const int prime, const u64 limit) {
    if (limit < static_cast<u64>(prime)) {
        return 0ULL;
    }

    u64 value = static_cast<u64>(prime);
    while (value <= limit / static_cast<u64>(prime)) {
        value *= static_cast<u64>(prime);
    }
    return value;
}

struct SolverData {
    std::vector<int> all_primes;
    std::vector<u64> max_prime_powers;
    std::vector<int> small_primes;
    std::vector<u64> small_prime_powers;
    std::vector<int> large_primes;
    std::vector<u64> large_prime_powers;
    u64 baseline_sum = 1ULL;
};

SolverData prepare_solver_data(const u64 n) {
    SolverData data;
    data.all_primes = sieve_primes(static_cast<int>(n));

    data.max_prime_powers.reserve(data.all_primes.size());
    for (const int p : data.all_primes) {
        const u64 power = max_prime_power_leq_n(n, p);
        data.max_prime_powers.push_back(power);
        data.baseline_sum += power;
    }

    const int split = static_cast<int>(std::sqrt(static_cast<long double>(n)));

    for (std::size_t i = 0; i < data.all_primes.size(); ++i) {
        const int p = data.all_primes[i];
        const u64 power = data.max_prime_powers[i];
        if (p <= split) {
            data.small_primes.push_back(p);
            data.small_prime_powers.push_back(power);
        } else {
            data.large_primes.push_back(p);
            data.large_prime_powers.push_back(power);
        }
    }

    return data;
}

std::vector<std::tuple<int, int, u64>> build_gain_edges_chunk(const SolverData& data,
                                                               const u64 n,
                                                               const std::size_t begin,
                                                               const std::size_t end) {
    std::vector<std::tuple<int, int, u64>> edges;

    for (std::size_t si = begin; si < end; ++si) {
        const int p = data.small_primes[si];
        const u64 p_power = data.small_prime_powers[si];

        for (std::size_t li = 0; li < data.large_primes.size(); ++li) {
            const int q = data.large_primes[li];
            const u64 q_power = data.large_prime_powers[li];

            if (static_cast<u64>(p) * static_cast<u64>(q) > n) {
                break;
            }

            const u64 best_p_factor = highest_power_with_limit(p, n / static_cast<u64>(q));
            if (best_p_factor == 0ULL) {
                continue;
            }

            const u64 pair_value = best_p_factor * static_cast<u64>(q);
            if (pair_value <= p_power + q_power) {
                continue;
            }

            const u64 gain = pair_value - p_power - q_power;
            edges.emplace_back(static_cast<int>(si), static_cast<int>(li), gain);
        }
    }

    return edges;
}

std::vector<std::tuple<int, int, u64>> build_gain_edges(const SolverData& data,
                                                        const u64 n,
                                                        const bool allow_multithreading,
                                                        const unsigned requested_threads) {
    const std::size_t small_count = data.small_primes.size();
    if (small_count == 0ULL || data.large_primes.empty()) {
        return {};
    }

    const unsigned threads = choose_thread_count(allow_multithreading, requested_threads, small_count);
    if (threads == 1U) {
        return build_gain_edges_chunk(data, n, 0ULL, small_count);
    }

    std::vector<std::vector<std::tuple<int, int, u64>>> partial(
        static_cast<std::size_t>(threads));
    std::vector<std::thread> pool;
    pool.reserve(static_cast<std::size_t>(threads));

    const std::size_t chunk = (small_count + static_cast<std::size_t>(threads) - 1ULL) /
                              static_cast<std::size_t>(threads);

    for (unsigned t = 0U; t < threads; ++t) {
        const std::size_t begin = static_cast<std::size_t>(t) * chunk;
        const std::size_t end = std::min<std::size_t>(small_count, begin + chunk);

        if (begin >= end) {
            break;
        }

        pool.emplace_back([&, t, begin, end]() {
            partial[static_cast<std::size_t>(t)] = build_gain_edges_chunk(data, n, begin, end);
        });
    }

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

    std::size_t total = 0ULL;
    for (const auto& part : partial) {
        total += part.size();
    }

    std::vector<std::tuple<int, int, u64>> edges;
    edges.reserve(total);
    for (auto& part : partial) {
        edges.insert(edges.end(), part.begin(), part.end());
    }

    return edges;
}

i64 hungarian_max_weight_bipartite(const int left_size,
                                   const int right_size,
                                   const std::vector<i64>& weights) {
    if (left_size == 0 || right_size == 0) {
        return 0;
    }

    const int n = left_size;
    const int m = right_size;

    auto cost = [&](const int i, const int j) -> i64 {
        const std::size_t idx = static_cast<std::size_t>(i - 1) * static_cast<std::size_t>(m) +
                                static_cast<std::size_t>(j - 1);
        return -weights[idx];
    };

    const i64 inf = std::numeric_limits<i64>::max() / 4;

    std::vector<i64> u(static_cast<std::size_t>(n) + 1ULL, 0);
    std::vector<i64> v(static_cast<std::size_t>(m) + 1ULL, 0);
    std::vector<int> p(static_cast<std::size_t>(m) + 1ULL, 0);
    std::vector<int> way(static_cast<std::size_t>(m) + 1ULL, 0);

    for (int i = 1; i <= n; ++i) {
        p[0] = i;
        int j0 = 0;

        std::vector<i64> minv(static_cast<std::size_t>(m) + 1ULL, inf);
        std::vector<bool> used(static_cast<std::size_t>(m) + 1ULL, false);

        do {
            used[static_cast<std::size_t>(j0)] = true;
            const int i0 = p[static_cast<std::size_t>(j0)];

            i64 delta = inf;
            int j1 = 0;

            for (int j = 1; j <= m; ++j) {
                if (used[static_cast<std::size_t>(j)]) {
                    continue;
                }

                const i64 cur = cost(i0, j) - u[static_cast<std::size_t>(i0)] - v[static_cast<std::size_t>(j)];
                if (cur < minv[static_cast<std::size_t>(j)]) {
                    minv[static_cast<std::size_t>(j)] = cur;
                    way[static_cast<std::size_t>(j)] = j0;
                }

                if (minv[static_cast<std::size_t>(j)] < delta) {
                    delta = minv[static_cast<std::size_t>(j)];
                    j1 = j;
                }
            }

            for (int j = 0; j <= m; ++j) {
                if (used[static_cast<std::size_t>(j)]) {
                    u[static_cast<std::size_t>(p[static_cast<std::size_t>(j)])] += delta;
                    v[static_cast<std::size_t>(j)] -= delta;
                } else {
                    minv[static_cast<std::size_t>(j)] -= delta;
                }
            }

            j0 = j1;
        } while (p[static_cast<std::size_t>(j0)] != 0);

        do {
            const int j1 = way[static_cast<std::size_t>(j0)];
            p[static_cast<std::size_t>(j0)] = p[static_cast<std::size_t>(j1)];
            j0 = j1;
        } while (j0 != 0);
    }

    i64 min_cost = 0;
    for (int j = 1; j <= m; ++j) {
        if (p[static_cast<std::size_t>(j)] != 0) {
            min_cost += cost(p[static_cast<std::size_t>(j)], j);
        }
    }

    return -min_cost;
}

u64 solve_co_optimized(const u64 n,
                       const bool allow_multithreading,
                       const unsigned requested_threads) {
    const SolverData data = prepare_solver_data(n);
    const auto edges = build_gain_edges(data, n, allow_multithreading, requested_threads);

    const int left_size = static_cast<int>(data.small_primes.size());
    const int large_size = static_cast<int>(data.large_primes.size());
    const int right_size = large_size + left_size;  // includes dummy columns

    std::vector<i64> weights(static_cast<std::size_t>(left_size) * static_cast<std::size_t>(right_size), 0);

    for (const auto& edge : edges) {
        const int small_index = std::get<0>(edge);
        const int large_index = std::get<1>(edge);
        const i64 gain = static_cast<i64>(std::get<2>(edge));

        const std::size_t idx = static_cast<std::size_t>(small_index) * static_cast<std::size_t>(right_size) +
                                static_cast<std::size_t>(large_index);
        if (gain > weights[idx]) {
            weights[idx] = gain;
        }
    }

    const i64 max_gain = hungarian_max_weight_bipartite(left_size, right_size, weights);
    return data.baseline_sum + static_cast<u64>(max_gain);
}

bool run_checkpoints() {
    struct Checkpoint {
        u64 n;
        u64 expected;
    };

    const std::vector<Checkpoint> checkpoints = {
        {10ULL, 30ULL},
        {30ULL, 193ULL},
        {100ULL, 1356ULL},
    };

    for (const Checkpoint& cp : checkpoints) {
        const u64 got = solve_co_optimized(cp.n, false, 1U);
        if (got != cp.expected) {
            std::cerr << "Checkpoint failed for n=" << cp.n << ": expected " << cp.expected
                      << ", got " << got << '\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 1;
    }

    const u64 answer = solve_co_optimized(options.n,
                                          options.allow_multithreading,
                                          options.requested_threads);

    std::cout << answer << '\n';
    return 0;
}

Python

import math

def sieve_primes(n):
    if n < 2: return []
    is_prime = bytearray(n + 1)
    for i in range(2, n + 1):
        is_prime[i] = 1
    for i in range(2, math.isqrt(n) + 1):
        if is_prime[i]:
            for j in range(i * i, n + 1, i):
                is_prime[j] = 0
    return [i for i in range(2, n + 1) if is_prime[i]]

def max_prime_power_leq_n(n, prime):
    val = prime
    while val <= n // prime:
        val *= prime
    return val

def highest_power_with_limit(prime, limit):
    if limit < prime: return 0
    val = prime
    while val <= limit // prime:
        val *= prime
    return val

def hungarian_max_weight_bipartite(left_size, right_size, weights):
    if left_size == 0 or right_size == 0: return 0
    
    n = left_size
    m = right_size
    
    INF = 10**18
    u = [0] * (n + 1)
    v = [0] * (m + 1)
    p = [0] * (m + 1)
    way = [0] * (m + 1)
    
    for i in range(1, n + 1):
        p[0] = i
        j0 = 0
        minv = [INF] * (m + 1)
        used = bytearray(m + 1)
        
        while True:
            used[j0] = 1
            i0 = p[j0]
            delta = INF
            j1 = 0
            
            for j in range(1, m + 1):
                if not used[j]:
                    cost = -weights[(i0 - 1) * m + (j - 1)]
                    cur = cost - u[i0] - v[j]
                    if cur < minv[j]:
                        minv[j] = cur
                        way[j] = j0
                    if minv[j] < delta:
                        delta = minv[j]
                        j1 = j
                        
            for j in range(m + 1):
                if used[j]:
                    u[p[j]] += delta
                    v[j] -= delta
                else:
                    minv[j] -= delta
                    
            j0 = j1
            if p[j0] == 0:
                break
                
        while True:
            j1 = way[j0]
            p[j0] = p[j1]
            j0 = j1
            if j0 == 0:
                break
                
    min_cost = 0
    for j in range(1, m + 1):
        if p[j] != 0:
            min_cost += -weights[(p[j] - 1) * m + (j - 1)]
            
    return -min_cost

def solve_n(n):
    all_primes = sieve_primes(n)
    
    max_prime_powers = []
    baseline_sum = 1
    for p in all_primes:
        power = max_prime_power_leq_n(n, p)
        max_prime_powers.append(power)
        baseline_sum += power
        
    split = math.isqrt(n)
    small_primes = []
    small_prime_powers = []
    large_primes = []
    large_prime_powers = []
    
    for i in range(len(all_primes)):
        p = all_primes[i]
        power = max_prime_powers[i]
        if p <= split:
            small_primes.append(p)
            small_prime_powers.append(power)
        else:
            large_primes.append(p)
            large_prime_powers.append(power)
            
    left_size = len(small_primes)
    large_size = len(large_primes)
    right_size = large_size + left_size
    
    weights = [0] * (left_size * right_size)
    
    for si in range(left_size):
        p = small_primes[si]
        p_power = small_prime_powers[si]
        for li in range(large_size):
            q = large_primes[li]
            q_power = large_prime_powers[li]
            if p * q > n:
                break
                
            best_p_factor = highest_power_with_limit(p, n // q)
            if best_p_factor == 0:
                continue
                
            pair_value = best_p_factor * q
            if pair_value <= p_power + q_power:
                continue
                
            gain = pair_value - p_power - q_power
            idx = si * right_size + li
            weights[idx] = max(weights[idx], gain)
            
    max_gain = hungarian_max_weight_bipartite(left_size, right_size, weights)
    return baseline_sum + max_gain

def solve():
    n = 200000
    ans = solve_n(n)
    return str(ans)

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

Java

import java.util.ArrayList;
import java.util.List;

public class Euler355 {

    static List<Integer> sievePrimes(int n) {
        List<Integer> primes = new ArrayList<>();
        if (n < 2)
            return primes;
        boolean[] isPrime = new boolean[n + 1];
        for (int i = 2; i <= n; i++)
            isPrime[i] = true;
        for (int i = 2; i * i <= n; i++) {
            if (isPrime[i]) {
                for (int j = i * i; j <= n; j += i) {
                    isPrime[j] = false;
                }
            }
        }
        for (int i = 2; i <= n; i++) {
            if (isPrime[i])
                primes.add(i);
        }
        return primes;
    }

    static long maxPrimePowerLeqN(long n, int prime) {
        long value = prime;
        while (value <= n / prime) {
            value *= prime;
        }
        return value;
    }

    static long highestPowerWithLimit(int prime, long limit) {
        if (limit < prime)
            return 0;
        long value = prime;
        while (value <= limit / prime) {
            value *= prime;
        }
        return value;
    }

    static long hungarianMaxWeightBipartite(int leftSize, int rightSize, long[] weights) {
        if (leftSize == 0 || rightSize == 0)
            return 0;
        int n = leftSize;
        int m = rightSize;
        long inf = Long.MAX_VALUE / 4;

        long[] u = new long[n + 1];
        long[] v = new long[m + 1];
        int[] p = new int[m + 1];
        int[] way = new int[m + 1];

        for (int i = 1; i <= n; i++) {
            p[0] = i;
            int j0 = 0;
            long[] minv = new long[m + 1];
            for (int k = 0; k <= m; k++)
                minv[k] = inf;
            boolean[] used = new boolean[m + 1];

            do {
                used[j0] = true;
                int i0 = p[j0];
                long delta = inf;
                int j1 = 0;

                for (int j = 1; j <= m; j++) {
                    if (!used[j]) {
                        long cost = -weights[(i0 - 1) * m + (j - 1)];
                        long cur = cost - u[i0] - v[j];
                        if (cur < minv[j]) {
                            minv[j] = cur;
                            way[j] = j0;
                        }
                        if (minv[j] < delta) {
                            delta = minv[j];
                            j1 = j;
                        }
                    }
                }

                for (int j = 0; j <= m; j++) {
                    if (used[j]) {
                        u[p[j]] += delta;
                        v[j] -= delta;
                    } else {
                        minv[j] -= delta;
                    }
                }
                j0 = j1;
            } while (p[j0] != 0);

            do {
                int j1 = way[j0];
                p[j0] = p[j1];
                j0 = j1;
            } while (j0 != 0);
        }

        long minCost = 0;
        for (int j = 1; j <= m; j++) {
            if (p[j] != 0) {
                minCost += -weights[(p[j] - 1) * m + (j - 1)];
            }
        }
        return -minCost;
    }

    static long solveN(long n) {
        List<Integer> allPrimes = sievePrimes((int) n);
        List<Long> maxPrimePowers = new ArrayList<>();
        long baselineSum = 1;

        for (int p : allPrimes) {
            long power = maxPrimePowerLeqN(n, p);
            maxPrimePowers.add(power);
            baselineSum += power;
        }

        int split = (int) Math.sqrt((double) n);
        List<Integer> smallPrimes = new ArrayList<>();
        List<Long> smallPrimePowers = new ArrayList<>();
        List<Integer> largePrimes = new ArrayList<>();
        List<Long> largePrimePowers = new ArrayList<>();

        for (int i = 0; i < allPrimes.size(); i++) {
            int p = allPrimes.get(i);
            long power = maxPrimePowers.get(i);
            if (p <= split) {
                smallPrimes.add(p);
                smallPrimePowers.add(power);
            } else {
                largePrimes.add(p);
                largePrimePowers.add(power);
            }
        }

        int leftSize = smallPrimes.size();
        int largeSize = largePrimes.size();
        int rightSize = largeSize + leftSize;
        long[] weights = new long[leftSize * rightSize];

        for (int si = 0; si < leftSize; si++) {
            int p = smallPrimes.get(si);
            long pPower = smallPrimePowers.get(si);
            for (int li = 0; li < largeSize; li++) {
                int q = largePrimes.get(li);
                long qPower = largePrimePowers.get(li);

                if ((long) p * q > n)
                    break;

                long bestPFactor = highestPowerWithLimit(p, n / q);
                if (bestPFactor == 0)
                    continue;

                long pairValue = bestPFactor * q;
                if (pairValue <= pPower + qPower)
                    continue;

                long gain = pairValue - pPower - qPower;
                weights[si * rightSize + li] = Math.max(weights[si * rightSize + li], gain);
            }
        }

        long maxGain = hungarianMaxWeightBipartite(leftSize, rightSize, weights);
        return baselineSum + maxGain;
    }

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