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

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$.

Step 2: Why Only Two-Prime Replacements Can Beat the Baseline

Suppose a candidate number $x\le n$ has distinct prime divisors in a set $T$. Choosing $x$ forbids us from keeping the baseline representatives $P_r$ for $r\in T$.

If $P_r=r^k$, then $r^{k+1}>n$, hence

$$P_r>\frac{n}{r}.$$

Therefore

$$\sum_{r\in T}P_r>n\sum_{r\in T}\frac{1}{r}.$$

If $x$ has at least three distinct prime factors, the smallest possible reciprocal sum is

$$\frac{1}{2}+\frac{1}{3}+\frac{1}{5}>1,$$

so $\sum_{r\in T}P_r>n\ge x$. Such a number can never improve the baseline. This is why the implementations only search upgrades involving exactly two primes.

Step 3: Small and Large Primes

Split the primes at $\lfloor\sqrt{n}\rfloor$: small primes satisfy $p\le \sqrt{n}$, while large primes satisfy $q>\sqrt{n}$.

Two large primes cannot appear in the same chosen number because then their product would exceed $n$. Also, for a large prime $q$, we have $P_q=q$ since $q^2>n$.

The local solver therefore encodes profitable mixed replacements as numbers of the form

$$M(p,q)=q\cdot \max\{p^a : a\ge 1,\ p^a q\le n\},\qquad p\le \sqrt{n}\lt q.$$

The helper highest_power_with_limit(p, n / q) computes the inner maximum. Replacing $P_p$ and $P_q$ by $M(p,q)$ changes the sum by

$$w(p,q)=M(p,q)-P_p-P_q.$$

Only positive values of $w(p,q)$ are useful, so only those become edges in the optimization graph.

Step 4: Reduction to Maximum-Weight Matching

Each small prime can be used in at most one mixed number, and each large prime can also be used at most once. That is exactly a matching constraint.

Create a bipartite graph with small primes on the left and large primes on the right. Add an edge $(p,q)$ with weight $w(p,q)$ whenever $w(p,q)>0$. The implementations add extra zero-weight dummy columns so that a small prime may remain unmatched and simply keep $P_p$.

If $\mathcal{M}$ is a matching, its total gain is

$$G(\mathcal{M})=\sum_{(p,q)\in \mathcal{M}}w(p,q).$$

The answer computed by the code is

$$\boxed{\operatorname{Co}(n)=S_0+\max_{\mathcal{M}}G(\mathcal{M}).}$$

The maximum is obtained with the Hungarian algorithm on a rectangular assignment matrix.

Worked Example: $n=30$

The baseline representatives are

$$1,\ 16,\ 27,\ 25,\ 7,\ 11,\ 13,\ 17,\ 19,\ 23,\ 29,$$

$$S_0=188.$$

Here $\lfloor\sqrt{30}\rfloor=5$, so the small primes are $2,3,5$. The only positive gain edge is $(2,7)$:

$$M(2,7)=7\cdot 4=28,\qquad w(2,7)=28-16-7=5.$$

Hence

$$\operatorname{Co}(30)=188+5=193,$$

which matches the checkpoint embedded in the C++ program.

Why Matching Is Needed: $n=100$

For $n=100$, the baseline sum is $S_0=1263$. Several small primes have multiple profitable partners, so a greedy choice is not enough. The maximum matching selected by the implementations uses

$$64+11\to 88,\qquad 25+19\to 95,\qquad 49+13\to 91,$$

with gains $13$, $51$, and $29$. The total gain is $93$, giving

$$\operatorname{Co}(100)=1263+93=1356.$$

How the Code Works

The three solution files implement the same mathematics. They sieve primes up to $n$, compute every $P_p$ with max_prime_power_leq_n, sum the baseline, and split primes at $\lfloor\sqrt{n}\rfloor$.

For each feasible pair $(p,q)$ with $p$ small, $q$ large, and $pq\le n$, they evaluate the best mixed value $M(p,q)$. Positive gains are written into a weight matrix of size

$$s\times(\ell+s),$$

where $s$ is the number of small primes and $\ell$ is the number of large primes. The extra $s$ columns are the dummy unmatched choices. The function hungarian_max_weight_bipartite returns the maximum additional gain.

The C++ version exposes this pipeline through prepare_solver_data, build_gain_edges, and solve_co_optimized. It also parallelizes edge generation over the small-prime side and verifies the checkpoints $\operatorname{Co}(10)=30$, $\operatorname{Co}(30)=193$, and $\operatorname{Co}(100)=1356$. The Python and Java versions follow the same matrix construction in a single thread.

Complexity Analysis

Let $s=\pi(\lfloor\sqrt{n}\rfloor)$ and $\ell=\pi(n)-s$. The sieve costs $O(n\log\log n)$ time and $O(n)$ memory. Building the gain matrix costs $O(E\log n)$, where $E$ is the number of feasible pairs $(p,q)$ with $p\le\sqrt{n}\lt q$ and $pq\le n$; the logarithmic factor comes from multiplying $p$ repeatedly up to the limit $n/q$.

The Hungarian algorithm on an $s\times(\ell+s)$ matrix costs $O(s^2(\ell+s))$ time and $O(s(\ell+s))$ memory. For the actual input $n=200000$, we have $s=\pi(447)=86$, so the matching stage is small compared with the full prime range.

References

Problem page: https://projecteuler.net/problem=355
Coprime integers: Wikipedia - Coprime integers
Prime power: Wikipedia - Prime power
Sieve of Eratosthenes: Wikipedia - Sieve of Eratosthenes
Hungarian algorithm: Wikipedia - Hungarian algorithm

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));
    }
}