Problem 418: Factorisation Triples

Project Euler Problem 418 Solution

EulerSolve provides an optimized solution for Project Euler Problem 418, Factorisation Triples, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For $n = 43!$, we seek ordered factor triples $(a,b,c)$ satisfying $$a \le b \le c,\qquad abc = n.$$ The primary objective is to make the triple as balanced as possible, which means minimizing the ratio $c/a$. If several triples achieve the same optimal ratio, the required output is the smallest possible value of $a+b+c$. The implementations do not enumerate every factor triple directly. Instead they work with the prime factorization of $43!$, search balanced exponent distributions first, and then refine the best result with an exact divisor search. Mathematical Approach Step 1: Prime exponents of $43!$ For each prime $p \le 43$, the exponent of $p$ in $43!$ is given by Legendre's formula $$e_p = \sum_{k \ge 1} \left\lfloor \frac{43}{p^k} \right\rfloor.$$ Hence $$43!...

Detailed mathematical approach

Problem Summary

A factorisation triple of a positive integer $n$ is a triple of integers $(a,b,c)$ with $1\le a\le b\le c$ and $abc=n$. Exactly one factorisation triple of $n$ minimises the ratio $c/a$, and $f(n)$ is its sum $a+b+c$. The statement gives the checkpoints $f(165)=19$, $f(100100)=142$ and $f(20!)=4034872$, and asks for $f(43!)$. ^[1]

The whole solution in one sentence: the smallest member $a$ is a divisor of $n$ no larger than $\sqrt[3]{n}$, fixing $a$ determines the best $b$ and $c$, and any ratio already achieved squeezes $a$ into a very short interval below $\sqrt[3]{n}$, so of the $516\,096\,000$ divisors of $43!$ only a handful ever needs its companions computed.

Mathematical Approach

1. The objective is the ratio, not the sum

The quantity to minimise is $c/a$; the sum $a+b+c$ is only reported afterwards. The two goals really differ, and the smallest $n$ at which they disagree is $n=40768=2^6\cdot7^2\cdot13$: the ratio is smallest for $(28,28,52)$, with $c/a=13/7$ and sum $108$, while the smallest sum is $107$, attained by $(26,32,49)$ with the larger ratio $49/26$. The catalogued sequence of smallest sums therefore agrees with $f(n)$ only up to $n=40767$. ^[3] Every comparison below tests $c/a$ first; only a genuine tie is settled by the sum.

The statement asserts that the minimising triple is unique. Exhaustive enumeration confirms this for every $n\le 200\,000$, a finite check rather than a proof; the programs keep a tie-break rule anyway, so they do not rely on the assertion.

2. The smallest member lives just below the cube root

From $a\le b\le c$ we get $a^3\le abc=n$, so every candidate for $a$ is a divisor of $n$ with

$$a\le\big\lfloor\sqrt[3]{n}\big\rfloor.$$

That bound is $10$ for $n=1008$, $46$ for $n=100100$, and $392\,387\,855\,098\,338\,513$ for $43!$, a number of $53$ digits. Exactly $36\,214\,968$ divisors of $43!$ do not exceed it, so the bound alone leaves far too much work.

3. With $a$ fixed, the best companions are determined

Fix a divisor $a$ of $n$ and put $q=n/a$, so $bc=q$. Minimising $c/a$ now means minimising $c=q/b$, that is making $b$ as large as possible, and $b\le c$ is exactly $b^2\le q$. The best companion is therefore

$$\beta(q)=\max\{\,d:\ d\mid q,\ d\le\sqrt{q}\,\},$$

the largest divisor of $q$ at most its square root, ^[4] with $c=q/\beta(q)$. Each $a$ thus contributes at most one candidate, and finding it is a second, smaller divisor problem rather than a search over pairs.

4. When a divisor is the smallest member of no triple

That triple is admissible only if $\beta(q)\ge a$, and the converse is clean: if $\beta(q)\lt a$ then $n$ has no factorisation triple whose smallest member is $a$, since any such triple has $b\mid q$ and $b\le\sqrt q$, hence $b\le\beta(q)\lt a$, contradicting $a\le b$. Such divisors can be discarded.

For $n=100100$ and $a=44$, for instance, $q=2275$ and $\beta(2275)=35\lt44$: no divisor of $2275$ lies between $44$ and $\sqrt{2275}=47.6\ldots$, so no triple of $100100$ begins with $44$.

5. A complete worked example: $n=1008$

Here $1008=2^4\cdot3^2\cdot7$ has $30$ divisors, of which the eight at most $\lfloor\sqrt[3]{1008}\rfloor=10$ are $1,2,3,4,6,7,8,9$. Section 3 gives $b=\beta(1008/a)$ and $c$ for each:

$a=1$: $b=28$, $c=36$, ratio $36$;
$a=2$: $b=21$, $c=24$, ratio $12$;
$a=3$: $b=16$, $c=21$, ratio $7$;
$a=4$: $b=14$, $c=18$, ratio $9/2$;
$a=6$: $b=12$, $c=14$, ratio $7/3$;
$a=7$: $b=12$, $c=12$, ratio $12/7$;
$a=8$: $b=9$, $c=14$, ratio $7/4$;
$a=9$: no triple, by Section 4.

The winner is $(7,12,12)$ with $12/7$, so $f(1008)=31$. The last two lines matter: the largest admissible divisor $a=8$ is not optimal, because $7/4\gt12/7$. The ratio is not a decreasing function of $a$, so the search cannot stop at the first admissible divisor it meets from above; it must try every candidate that a valid bound does not exclude.

The checkpoints of the statement come out the same way. For $n=165=3\cdot5\cdot11$ only $a\in\{1,3,5\}$ is possible: $a=5$ fails by Section 4, $a=1$ gives ratio $15$, and $a=3$ gives $(3,5,11)$ with ratio $11/3$, so $f(165)=19$. For $n=100100=2^2\cdot5^2\cdot7\cdot11\cdot13$, of the $16$ divisors up to $46$ the best is $a=35$, with $q=2860$, $\beta(2860)=52$, $c=55$, ratio $11/7$ and $f(100100)=142$.

6. A ratio already achieved confines $a$

This is the step that makes $43!$ tractable. Let a triple of ratio $R=c^{*}/a^{*}$ be known and let $(a,b,c)$ be at least as good, so $c/a\le R$. From $b\le c$ and $c\le Ra$,

$$n=abc\le ac^2\le a(Ra)^2=R^2a^3,$$

so every triple at least as good as the known one has

$$a\ \ge\ \Big(\frac{n}{R^{2}}\Big)^{1/3}.$$

With Section 2 the range for $a$ runs from $(n/R^2)^{1/3}$ to $\sqrt[3]{n}$, of relative width about $\tfrac23\log R$ when $R$ is near $1$. Two remarks. The inequality must stay non-strict, because both of its steps become equalities when $b=c$ and $c/a=R$, and that happens: for $n=1008$ and $R=12/7$ the bound is $(1008\cdot49/144)^{1/3}=343^{1/3}=7$ exactly, which is the optimal $a$, so a strict test would discard it. And the bound may be recomputed whenever $R$ improves, since a smaller $R$ only raises the lower end.

The gain shows on the small cases: for $1008$ the interval keeps $\{7,8,9\}$ of eight divisors, and for $100100$ with $R=11/7$ the bound is $34.3517\ldots$ and keeps $\{35,44\}$ of sixteen. For $20!$, with $41\,040$ divisors of which $5327$ lie below $\lfloor\sqrt[3]{20!}\rfloor=1\,344\,956$, the optimal ratio $168399/167960\approx1.0026137$ leaves three: $1\,343\,034$, $1\,343\,680$ and $1\,344\,000$. The optimum $(1\,343\,680,\ 1\,344\,000,\ 1\,347\,192)$ sums to the $4034872$ of the statement.

7. $43!$ in exponents

The exponent of a prime $p$ in $n!$ is $\sum_{i\ge1}\lfloor n/p^i\rfloor$, Legendre's formula, since each multiple of $p^i$ up to $n$ adds one more factor $p$. ^[2] This gives

$$\begin{aligned}43!={}&2^{39}3^{19}5^{9}7^{6}11^{3}13^{3}\\&{}\cdot17^{2}19^{2}\cdot23\cdot29\\&{}\cdot31\cdot37\cdot41\cdot43 .\end{aligned}$$

The number of divisors is the product of the exponents increased by one, ^[7] that is $40\cdot20\cdot10\cdot7\cdot4\cdot4\cdot3\cdot3\cdot2^6$, which equals $516\,096\,000$. A triple splits each $e_p$ into three parts $\alpha_p+\beta_p+\gamma_p=e_p$: choosing $\alpha_p$ leaves $e_p-\alpha_p+1$ choices for $\beta_p$, $\gamma_p$ follows, and these counts add to $\binom{e_p+2}{2}$. Over fourteen primes that is $695\,959\,387\,200\,000$ ordered splits, about $7\cdot10^{14}$ — far too many to enumerate, so a cheap usable $R$ is needed.

8. Logarithms, and why perfect balance is hard

Logarithms turn the three products into three sums: with $L_i$ the logarithm of the $i$-th member, minimising $c/a$ is minimising $L_3-L_1$ for the sorted values. The task becomes splitting a multiset — $39$ copies of $\log 2$, $19$ of $\log 3$, and so on — into three groups whose largest and smallest totals are as close as possible: the gap objective of multiway number partitioning, NP-hard in general. ^[5] That motivates the heuristic and the search below, but it is no claim that these instances, whose items are logarithms of primes, inherit the hardness; proving that would need a reduction the other way.

The standard heuristic takes the items from large to small and puts each into the group with the smallest total so far. It is cheap and often excellent: for $1008$ it returns $\{7,12,12\}$, for $100100$ it returns $\{35,52,55\}$ and for $165$ $\{3,5,11\}$, optimal in all three. For $43!$ it does much worse, a log-range of $0.669426$ and ratio $1.953115$; for $20!$ it reaches $1.491228$. Its weakness is that it never revisits an assignment: the thirty-nine copies of $\log2=0.693147\ldots$ come last, each moving one group by that much, only partly repairing the imbalance built up by the large primes. That is no lower bound on what is attainable — for $n=8$ the rule gives $(2,2,2)$ and gap $0$, and the $1008$ seed has log-range $0.538997$, below $\log2$. But by Section 6 a ratio near $2$ leaves an interval reaching a third of the way down from the cube root.

9. Branch and bound over exponent splits

The first phase improves the heuristic triple by a depth-first search over the splits of Section 7, one prime at a time, keeping the three partial logarithms sorted so the bound below can read off the smallest, middle and largest directly. Sorting removes no duplicates: every ordered split is entered, and at an empty node $(1,0,0)$, $(0,1,0)$ and $(0,0,1)$ give the same sorted state. The splits of a prime are tried by increasing spread $\max-\min$ of the exponents, ties broken by closeness to $e_p/3$, so balanced choices come first.

Pruning needs a quantity that can only underestimate what a branch may still achieve. ^[6] Let the sorted partial logarithms be $L_1\le L_2\le L_3$, let $M$ be the total logarithm of the unassigned primes, and put $D_1=L_2-L_1$, $D_2=D_1+2(L_3-L_2)$. Allowing $M$ to be poured out continuously rather than in fixed lumps, define

$\Lambda=L_1+M$ when $M\le D_1$;
$\Lambda=L_2+\tfrac12(M-D_1)$ when $D_1\lt M\le D_2$;
$\Lambda=L_3$ when $M\gt D_2$.

In the first two cases the pouring lifts $L_1$ to the level of $L_2$ and then both towards $L_3$, so $\Lambda$ is the highest level the minimum can reach; as the maximum stays at least $L_3$, the range is at least $L_3-\Lambda$. For $(0,1,3)$ with $M=3$ that gives $\Lambda=2$ and the bound $1$. In the third case the minimum can pass $L_3$ — from $(0,0,0)$ with $M=3$ all three reach $1$ — so $\Lambda$ no longer bounds it, yet the returned $L_3-\Lambda=0$ is still valid, a range never being negative. A relaxation can only be optimistic, so discarding a branch whose bound matches the best known range never discards the optimum. On $1500$ random partial states of five small numbers the bound never exceeded the true best range reachable from the state.

The search is capped at three million admitted nodes, its purpose being only a good $R$; whatever it returns is feasible and hence a valid bound. Which triple it finds is implementation-dependent, since the two sort keys tie — three splits of $e=2$ share both — and a non-stable sort may order those either way. The C++ build measured here reaches a log-range of $0.00138609$, ratio $1.0013871$, after $3\,001\,128$ nodes, counting calls that return at the cap; Python, whose sort is stable and which counts only admitted calls, stops at $3\,000\,000$ nodes with $0.0023320015$, ratio $1.0023347$.

10. What the two phases are worth

With the first of those ratios, Section 6 confines $a$ to the interval from $392\,025\,432\,652\,787\,841$ to $392\,387\,855\,098\,338\,513$, of relative width $9.24\cdot10^{-4}$, holding $4920$ divisors of $43!$ instead of $36\,214\,968$; the weaker Python seed gives $8272$, the heuristic alone $2\,319\,856$. The second phase enumerates those divisors, tightening the lower end at every improvement. The final ratio exceeds $1$ by less than $4.5\cdot10^{-6}$ and its interval holds only $14$ divisors; the measured run evaluates $28$ candidates, of which $18$ fail the test of Section 4.

11. What is exact and what is not

In the second phase C++ and Java compare finished triples exactly, testing $c/a\lt c^{*}/a^{*}$ as $c\,a^{*}\lt c^{*}a$ in integers; the first phase compares log-ranges, as Python does in both phases. An inexact first-phase comparison cannot change the result: whatever triple it keeps is feasible and only serves as the bound $R$, and the exact second-phase comparisons decide among the survivors.

The pruning always reads logarithms. Its margin is a candidate's distance above the current lower end, at least $\tfrac13\log(c/b)$ by Section 6 and equal to $2.2\cdot10^{-7}$ for the optimum of $43!$, while the spacing of floating-point numbers near $\log(43!)=121.53$ is $1.42\cdot10^{-14}$ — on this machine also for long double, whose significand is the $53$ bits of a double. That leaves about seven orders of magnitude of room. But the margin falls to zero when $b=c$, and the tolerances are fixed rather than certified, so this is evidence and not a proof; of the $2999$ instances swept below, $125$ have $b=c$ at the optimum and all were found.

How the Code Works

The C++ entry point ^[9] first runs run_checkpoints unless --skip-checkpoints is given, then calls factorize_factorial(43), which sieves the primes and applies Legendre's formula, and hands the primes and exponents to Solver. Its constructor caches rem_log, the suffix sums of $e_p\log p$ serving as the $M$ of Section 9; the prime powers in three widths, pow_u64, pow_u128 and pow_cpp; desc_idx, the primes in decreasing order; and split_by_exp, listing for every $e$ up to $39$ its $\binom{e+2}{2}$ splits sorted by spread and then by $|3x-e|+|3y-e|+|3z-e|$, the $11\,480$ triples of Section 9.

greedy_seed implements Section 8, adding $\log p$ to the smallest of the three totals for each copy of each prime, large primes first; sort_bins keeps the partial triples ordered and update_best_from_bins records a new best. phase_a_dfs is the search of Section 9, counting nodes against phase_a_cap and calling lower_bound_range for $L_3-\Lambda$.

phase_b_search is Section 10. cube_root_floor gives the upper end from $\exp(\log n/3)$ corrected by exact cubing, and lower_a_bound finds the lower end by binary search, as the smallest $a$ with $(c^{*})^2a^3\ge n(a^{*})^2$ — Section 6 cleared of roots and fractions, hence exact. The recursion walks the exponents of $a$ prime by prime, abandoning a branch when $a$ exceeds the upper end or when all the remaining primes cannot lift it to the lower end; that integer bound is computed once, while its logarithmic form, used at every node, is recomputed from the current best, so an improvement tightens the sweep at once. evaluate_candidate then does Section 3: sqrt_floor_cpp gives $\lfloor\sqrt q\rfloor$, and max_divisor_leq finds $\beta(q)$ by a second branch and bound over the residual exponents, taking the primes in decreasing order and cutting a branch once the partial product times the largest possible remainder cannot beat the best divisor so far.

The three versions share these reductions but differ in their arithmetic and in how the first triple is recorded. C++ carries $43!$ and the triple in Boost's cpp_int, ^[8] uses unsigned __int128 for the overflow-free products inside max_divisor_leq, and is the only version that verifies itself, through run_checkpoints. Java mirrors it with BigInteger, replaces the $128$-bit products by multiplyClamp, which saturates at Long.MAX_VALUE and so errs towards not pruning, and computes $\lfloor\sqrt q\rfloor$ with the Newton iteration sqrtFloorLong before narrowing to long — safe only because the root stays below $4\cdot10^{17}$ here, as longValue silently keeps the low bits of a larger value. ^[8] Python needs no special integer type and takes roots from math.isqrt, but replaces the exact binary search by the logarithmic formula and compares log-ranges with a tolerance of $10^{-18}$. That tolerance is its one real fragility: update recomputes the log-range from the integers while the heuristic stored it from accumulated logarithms, and when the two differ by more than $10^{-18}$ the first triple is never recorded. This does not affect $43!$, but with the factorisation step parameterised the same search returns $0$ for $19$ inputs in $2\le n\le3000$, the smallest $250$; C++ and Java are immune, recording unconditionally while no triple is stored.

Complexity and Verification

No bound polynomial in $\log n$ is available, but the shape of the work can be bounded. With $r$ distinct primes and $\tau$ divisors, the walk over $a$ traverses the divisor tree of $n$, whose node count $\sum_{k=0}^{r}\prod_{i\lt k}(e_i+1)$ is at most $(r+1)\tau$ and equals $1\,028\,000\,841$ for $43!$; every surviving $a$ starts a second such walk over the divisors of $n/a$, so a crude ceiling is $O(r\tau^2)$ large-integer operations. Pruning does far better: the outer walk shrinks to $54\,455\,301$ nodes and $9\,951\,504$ completed divisors, of which $28$ pass the interval test, and those $28$ inner walks explore $2\,748\,227$ nodes in all. A short interval does not, though, make the outer walk proportional to the survivors. Memory stays small: $11\,480$ split triples, $\sum_i(e_i+1)$ cached powers per representation and a recursion $r=14$ deep, a few hundred kilobytes in C++ and more in Java and Python, where every entry is an object.

Measured on one machine, the heuristic takes $0.013$ ms, the first phase $29.8$ ms and the second $170.8$ ms, so the second dominates although it inspects the fewest candidates. Lowering the node cap to $100\,000$, $1000$ or $0$ leaves the result and the time unchanged to within a few milliseconds and only raises the evaluated candidates from $28$ to $33$: the recomputed lower end does nearly all the pruning itself. The whole C++ run, checkpoints included, takes about $0.23$ s, Java $0.6$ s and Python $20.5$ s.

Several independent checks were run for this article. An exhaustive search over all factorisation triples agrees with Sections 3 and 4 on $165$, $100100$ and every $n$ from $2$ to $3000$, that last comparison made against the published solver class with only its entry point replaced by a loop. All three programs reproduce the three checkpoints and agree on $43!$, and $f(1),\dots,f(75)$ match the catalogued sums term by term. ^[3] The bound of Section 9 survived $1500$ partial states, the minimum of $c/a$ proved unique for every $n\le200\,000$, and Java's Newton root matched BigInteger.sqrt on $982\,637$ values, perfect squares and their neighbours included.

Footnotes and References

The references supply the background facts used above. The reduction of a fixed $a$ to the largest divisor below a square root, the criterion of Section 4, the interval bound of Section 6 and the relaxation bound of Section 9 are derived in this article.

Project Euler 418 — Factorisation Triples. The official statement calls (a, b, c) a factorisation triple of n when 1 ≤ a ≤ b ≤ c and abc = n, defines f(n) as a + b + c for the triple minimising c/a, states that it is unique, gives f(165) = 19, f(100100) = 142 and f(20!) = 4034872, and asks for f(43!).
Wikipedia — Legendre's formula. The exponent of a prime p in n! equals the sum of ⌊n/p^i⌋ over i ≥ 1, whose terms vanish once p^i exceeds n, as used in Section 7; the equivalent form (n − s_p(n))/(p − 1) with s_p(n) the base-p digit sum is given too.
OEIS A227215 — smallest sum of the three sides of a box of volume n. Its terms 3, 4, 5, 5, 7, 6, 9, 6, 7, 8, 13, … minimise the sum a + b + c rather than the ratio c/a. The first 75 agree with f(n) as computed here; the first disagreement is the n = 40768 of Section 1.
OEIS A033676 — largest divisor of n at most √n. This is the β of Section 3, the smaller central divisor of n; the entry notes that it is √n exactly when n is a square and 1 when n is 1 or prime. Its terms 1, 1, 1, 2, 1, 2, 1, 2, 3, 2, … agree with the values computed here.
Wikipedia — Multiway number partitioning. Splitting a multiset of numbers into k subsets with sums as equal as possible is NP-hard under each of its usual objectives, among them minimising the gap between the largest and smallest subset sum, the objective of Section 8. The article also describes the greedy rule used as a seed here: each item, taken from large to small, goes into the currently smallest subset.
Wikipedia — Branch and bound. A subproblem may be discarded when its bound is no better than the best solution already found, and the bounding function must be optimistic, never worse than the true value of any solution in that subproblem. Section 9 verifies this for the water-level bound.
Wikipedia — Divisor function. For n with prime exponents a_1, …, a_r, the number of divisors σ₀(n) is the product of the a_i + 1, applied to 43! in Section 7.
Boost.Multiprecision — cpp_int documents cpp_int as a signed arbitrary-precision integer that allocates more storage as values grow; C++ uses it for 43! and the exact ratio comparison. Java — BigInteger documents the immutable arbitrary-precision integers of the Java version and warns that longValue keeps only the low-order 64 bits of a value too large for a long, the narrowing discussed above.
C++ — Euler418.cpp. The linked immutable revision contains greedy_seed, lower_bound_range, phase_a_dfs, lower_a_bound, max_divisor_leq, evaluate_candidate and run_checkpoints. It is the implementation record for this article; the sweep over n ≤ 3000 and the measurements above are additional validation, not packaged in that file.

Problem 418 source code

C++

#include <algorithm>
#include <array>
#include <boost/multiprecision/cpp_int.hpp>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <vector>

namespace {

using boost::multiprecision::cpp_int;
using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr int MAXP = 20;

struct Triple {
    cpp_int a = 0;
    cpp_int b = 0;
    cpp_int c = 0;
    cpp_int sum = 0;
};

cpp_int cube_cpp(const u64 x) {
    cpp_int t = x;
    return t * t * t;
}

long double lower_bound_range(const long double l1, const long double l2, const long double l3, long double rem) {
    const long double t = l2 - l1;
    long double min_value;
    if (rem <= t) {
        min_value = l1 + rem;
    } else {
        rem -= t;
        const long double need = 2.0L * (l3 - l2);
        if (rem <= need) {
            min_value = l2 + rem / 2.0L;
        } else {
            min_value = l3;
        }
    }
    return l3 - min_value;
}

std::string to_string_cpp(const cpp_int& v) {
    return v.convert_to<std::string>();
}

struct Solver {
    std::vector<u64> primes;
    std::vector<int> exps;
    int pcount = 0;

    cpp_int n_value = 1;
    long double log_n = 0.0L;
    std::vector<long double> log_p;
    std::vector<long double> rem_log;

    std::vector<std::vector<u64>> pow_u64;
    std::vector<std::vector<u128>> pow_u128;
    std::vector<std::vector<cpp_int>> pow_cpp;

    std::vector<std::vector<std::array<int, 3>>> split_by_exp;
    std::vector<int> desc_idx;

    long long phase_a_nodes = 0;
    long long phase_a_cap = 3000000;
    long double best_log_range = 1e100L;
    Triple best;

    explicit Solver(std::vector<u64> primes_in, std::vector<int> exps_in)
        : primes(std::move(primes_in)), exps(std::move(exps_in)), pcount(static_cast<int>(primes.size())) {
        log_p.resize(pcount);
        rem_log.assign(pcount + 1, 0.0L);
        pow_u64.resize(pcount);
        pow_u128.resize(pcount);
        pow_cpp.resize(pcount);
        desc_idx.resize(pcount);
        std::iota(desc_idx.begin(), desc_idx.end(), 0);
        std::sort(desc_idx.begin(), desc_idx.end(), [&](const int i, const int j) { return primes[i] > primes[j]; });

        int max_e = 0;
        for (int i = 0; i < pcount; ++i) {
            log_p[i] = std::log(static_cast<long double>(primes[i]));
            log_n += static_cast<long double>(exps[i]) * log_p[i];
            n_value *= cpp_int(1);
            max_e = std::max(max_e, exps[i]);

            pow_u64[i].assign(static_cast<std::size_t>(exps[i] + 1), 1ULL);
            pow_u128[i].assign(static_cast<std::size_t>(exps[i] + 1), 1U);
            pow_cpp[i].assign(static_cast<std::size_t>(exps[i] + 1), cpp_int(1));
            for (int e = 1; e <= exps[i]; ++e) {
                pow_u64[i][static_cast<std::size_t>(e)] =
                    pow_u64[i][static_cast<std::size_t>(e - 1)] * primes[i];
                pow_u128[i][static_cast<std::size_t>(e)] =
                    pow_u128[i][static_cast<std::size_t>(e - 1)] * static_cast<u128>(primes[i]);
                pow_cpp[i][static_cast<std::size_t>(e)] =
                    pow_cpp[i][static_cast<std::size_t>(e - 1)] * cpp_int(primes[i]);
            }
            n_value *= pow_cpp[i][static_cast<std::size_t>(exps[i])];
        }

        for (int i = pcount - 1; i >= 0; --i) {
            rem_log[static_cast<std::size_t>(i)] =
                rem_log[static_cast<std::size_t>(i + 1)] + static_cast<long double>(exps[i]) * log_p[i];
        }

        split_by_exp.assign(static_cast<std::size_t>(max_e + 1), {});
        for (int e = 0; e <= max_e; ++e) {
            auto& v = split_by_exp[static_cast<std::size_t>(e)];
            for (int x = 0; x <= e; ++x) {
                for (int y = 0; y <= e - x; ++y) {
                    int z = e - x - y;
                    v.push_back({x, y, z});
                }
            }
            std::sort(v.begin(), v.end(), [&](const auto& a, const auto& b) {
                const int sa = std::max({a[0], a[1], a[2]}) - std::min({a[0], a[1], a[2]});
                const int sb = std::max({b[0], b[1], b[2]}) - std::min({b[0], b[1], b[2]});
                if (sa != sb) return sa < sb;
                const int ma = std::abs(3 * a[0] - e) + std::abs(3 * a[1] - e) + std::abs(3 * a[2] - e);
                const int mb = std::abs(3 * b[0] - e) + std::abs(3 * b[1] - e) + std::abs(3 * b[2] - e);
                return ma < mb;
            });
        }
    }

    static void sort_bins(std::array<long double, 3>& logs,
                          std::array<std::array<std::uint8_t, MAXP>, 3>& bins) {
        std::array<int, 3> ord = {0, 1, 2};
        std::sort(ord.begin(), ord.end(), [&](const int i, const int j) { return logs[i] < logs[j]; });
        std::array<long double, 3> nl;
        std::array<std::array<std::uint8_t, MAXP>, 3> nb;
        for (int i = 0; i < 3; ++i) {
            nl[i] = logs[ord[i]];
            nb[i] = bins[ord[i]];
        }
        logs = nl;
        bins = nb;
    }

    cpp_int value_from_bin(const std::array<std::uint8_t, MAXP>& bin_exp) const {
        cpp_int v = 1;
        for (int i = 0; i < pcount; ++i) {
            const int e = static_cast<int>(bin_exp[i]);
            if (e > 0) v *= pow_cpp[i][static_cast<std::size_t>(e)];
        }
        return v;
    }

    void update_best_from_bins(const std::array<long double, 3>& logs,
                               const std::array<std::array<std::uint8_t, MAXP>, 3>& bins) {
        const long double range = logs[2] - logs[0];
        if (range > best_log_range + 1e-18L) return;

        std::array<cpp_int, 3> v = {value_from_bin(bins[0]), value_from_bin(bins[1]), value_from_bin(bins[2])};
        std::sort(v.begin(), v.end());
        const cpp_int sum = v[0] + v[1] + v[2];

        if (range < best_log_range - 1e-18L || best.sum == 0 || sum < best.sum) {
            best_log_range = range;
            best.a = v[0];
            best.b = v[1];
            best.c = v[2];
            best.sum = sum;
        }
    }

    void greedy_seed() {
        std::array<long double, 3> logs = {0.0L, 0.0L, 0.0L};
        std::array<std::array<std::uint8_t, MAXP>, 3> bins{};
        for (auto& row : bins) row.fill(0U);

        std::vector<int> idx(pcount);
        std::iota(idx.begin(), idx.end(), 0);
        std::sort(idx.begin(), idx.end(), [&](const int i, const int j) { return primes[i] > primes[j]; });

        for (const int pi : idx) {
            for (int c = 0; c < exps[pi]; ++c) {
                int target = 0;
                if (logs[1] < logs[target]) target = 1;
                if (logs[2] < logs[target]) target = 2;
                logs[target] += log_p[pi];
                ++bins[target][pi];
            }
        }
        sort_bins(logs, bins);
        best_log_range = logs[2] - logs[0];
        update_best_from_bins(logs, bins);
    }

    void phase_a_dfs(int idx, const std::array<long double, 3>& logs,
                     const std::array<std::array<std::uint8_t, MAXP>, 3>& bins) {
        if (phase_a_nodes++ >= phase_a_cap) return;
        if (idx == pcount) {
            update_best_from_bins(logs, bins);
            return;
        }
        if (lower_bound_range(logs[0], logs[1], logs[2], rem_log[static_cast<std::size_t>(idx)]) >=
            best_log_range - 1e-18L) {
            return;
        }

        const int e = exps[idx];
        const auto& splits = split_by_exp[static_cast<std::size_t>(e)];
        for (const auto& s : splits) {
            std::array<long double, 3> nl = logs;
            std::array<std::array<std::uint8_t, MAXP>, 3> nb = bins;

            nl[0] += static_cast<long double>(s[0]) * log_p[idx];
            nl[1] += static_cast<long double>(s[1]) * log_p[idx];
            nl[2] += static_cast<long double>(s[2]) * log_p[idx];

            nb[0][idx] = static_cast<std::uint8_t>(nb[0][idx] + s[0]);
            nb[1][idx] = static_cast<std::uint8_t>(nb[1][idx] + s[1]);
            nb[2][idx] = static_cast<std::uint8_t>(nb[2][idx] + s[2]);

            sort_bins(nl, nb);
            phase_a_dfs(idx + 1, nl, nb);
        }
    }

    u64 cube_root_floor() const {
        long double approx = std::exp(log_n / 3.0L);
        if (approx < 1.0L) approx = 1.0L;
        u64 x = static_cast<u64>(approx);
        while (cube_cpp(x + 1) <= n_value) ++x;
        while (cube_cpp(x) > n_value) --x;
        return x;
    }

    u64 sqrt_floor_cpp(const cpp_int& x) const {
        long double approx = std::sqrt(x.convert_to<long double>());
        if (approx < 0.0L) approx = 0.0L;
        u64 r = static_cast<u64>(approx);
        while (cpp_int(r + 1U) * cpp_int(r + 1U) <= x) ++r;
        while (cpp_int(r) * cpp_int(r) > x) --r;
        return r;
    }

    u64 lower_a_bound(const u64 high) const {
        const cpp_int rhs = n_value * best.a * best.a;
        const cpp_int c2 = best.c * best.c;
        u64 lo = 1;
        u64 hi = high;
        while (lo < hi) {
            const u64 mid = lo + (hi - lo) / 2;
            const cpp_int lhs = c2 * cube_cpp(mid);
            if (lhs >= rhs) {
                hi = mid;
            } else {
                lo = mid + 1;
            }
        }
        return lo;
    }

    u64 max_divisor_leq(const std::array<int, MAXP>& residual, const u64 limit) const {
        std::array<int, MAXP> e_desc{};
        for (int i = 0; i < pcount; ++i) {
            e_desc[i] = residual[desc_idx[i]];
        }

        std::array<u128, MAXP + 1> rem_max{};
        rem_max[pcount] = 1;
        for (int i = pcount - 1; i >= 0; --i) {
            const int pi = desc_idx[i];
            rem_max[i] = rem_max[i + 1] * pow_u128[pi][static_cast<std::size_t>(e_desc[i])];
        }

        u64 best_b = 0;
        auto dfs = [&](auto&& self, const int idx, const u64 cur) -> void {
            if (cur > limit) return;
            if (idx == pcount) {
                if (cur > best_b) best_b = cur;
                return;
            }
            if (static_cast<u128>(cur) * rem_max[idx] <= static_cast<u128>(best_b)) return;

            const int pi = desc_idx[idx];
            const auto& pw = pow_u64[pi];
            const int e = e_desc[idx];
            const u64 lim = limit / cur;
            int kmax = e;
            while (kmax > 0 && pw[static_cast<std::size_t>(kmax)] > lim) --kmax;

            for (int k = kmax; k >= 0; --k) {
                const u64 nxt = cur * pw[static_cast<std::size_t>(k)];
                if (static_cast<u128>(nxt) * rem_max[idx + 1] <= static_cast<u128>(best_b)) break;
                self(self, idx + 1, nxt);
            }
        };
        dfs(dfs, 0, 1ULL);
        return best_b;
    }

    void evaluate_candidate(const u64 a, const std::array<std::uint8_t, MAXP>& a_exp) {
        cpp_int qa = n_value / a;
        const u64 u = sqrt_floor_cpp(qa);

        std::array<int, MAXP> residual{};
        for (int i = 0; i < pcount; ++i) residual[i] = exps[i] - static_cast<int>(a_exp[i]);
        const u64 b = max_divisor_leq(residual, u);
        if (b < a) return;

        cpp_int c = qa / b;
        if (c < b) return;

        const cpp_int lhs = c * best.a;
        const cpp_int rhs = best.c * a;
        if (lhs < rhs) {
            best.a = a;
            best.b = b;
            best.c = c;
            best.sum = cpp_int(a) + cpp_int(b) + c;
            best_log_range = std::log(best.c.convert_to<long double>()) - std::log(best.a.convert_to<long double>());
        } else if (lhs == rhs) {
            const cpp_int s = cpp_int(a) + cpp_int(b) + c;
            if (s < best.sum) {
                best.a = a;
                best.b = b;
                best.c = c;
                best.sum = s;
            }
        }
    }

    void phase_b_search() {
        const u64 a_high = cube_root_floor();
        const u64 base_a_low = lower_a_bound(a_high);
        const long double base_log_low = std::log(static_cast<long double>(base_a_low));
        const long double log_high = std::log(static_cast<long double>(a_high));

        std::array<std::uint8_t, MAXP> a_exp{};
        a_exp.fill(0U);

        auto dfs = [&](auto&& self, const int idx, const u64 cur, const long double lv) -> void {
            const long double dynamic_log_low =
                std::max(base_log_low, (log_n - 2.0L * best_log_range) / 3.0L);
            if (idx == pcount) {
                if (lv + 1e-18L >= dynamic_log_low && cur <= a_high) {
                    evaluate_candidate(cur, a_exp);
                }
                return;
            }
            if (lv > log_high + 1e-18L) return;
            if (lv + rem_log[static_cast<std::size_t>(idx)] < dynamic_log_low - 1e-18L) return;

            const u64 p = primes[idx];
            const int e = exps[idx];
            u64 mul = 1;
            for (int k = 0; k <= e; ++k) {
                if (cur > a_high / mul) break;
                const u64 nxt = cur * mul;
                a_exp[idx] = static_cast<std::uint8_t>(k);
                self(self, idx + 1, nxt, lv + static_cast<long double>(k) * log_p[idx]);
                if (k < e) mul *= p;
            }
            a_exp[idx] = 0U;
        };
        dfs(dfs, 0, 1ULL, 0.0L);
    }

    cpp_int solve() {
        greedy_seed();

        std::array<long double, 3> logs = {0.0L, 0.0L, 0.0L};
        std::array<std::array<std::uint8_t, MAXP>, 3> bins{};
        for (auto& row : bins) row.fill(0U);
        phase_a_nodes = 0;
        phase_a_dfs(0, logs, bins);

        phase_b_search();
        return best.sum;
    }
};

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

std::pair<std::vector<u64>, std::vector<int>> factorize_factorial(const int n) {
    std::vector<u64> primes;
    std::vector<int> exps;
    std::vector<bool> sieve(static_cast<std::size_t>(n + 1), true);
    sieve[0] = false;
    if (n >= 1) sieve[1] = false;
    for (int i = 2; i <= n; ++i) {
        if (!sieve[static_cast<std::size_t>(i)]) continue;
        for (int j = i + i; j <= n; j += i) sieve[static_cast<std::size_t>(j)] = false;
        primes.push_back(static_cast<u64>(i));
        int e = 0;
        int t = n;
        while (t > 0) {
            t /= i;
            e += t;
        }
        exps.push_back(e);
    }
    return {primes, exps};
}

bool run_checkpoints() {
    {
        auto [p, e] = factorize_number(165ULL);
        Solver solver(p, e);
        if (solver.solve() != cpp_int(19)) return false;
    }
    {
        auto [p, e] = factorize_number(100100ULL);
        Solver solver(p, e);
        if (solver.solve() != cpp_int(142)) return false;
    }
    {
        auto [p, e] = factorize_factorial(20);
        Solver solver(p, e);
        if (solver.solve() != cpp_int(4034872)) return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        std::cerr << "Checkpoint failed\n";
        return 2;
    }

    auto [primes, exps] = factorize_factorial(43);
    Solver solver(primes, exps);
    std::cout << to_string_cpp(solver.solve()) << '\n';
    return 0;
}

Python

import math

def solve():
    # Factorize 43!
    n_fact = 43
    sieve = [True]*(n_fact+1); sieve[0] = sieve[1] = False
    for i in range(2, n_fact+1):
        if not sieve[i]: continue
        for j in range(i+i, n_fact+1, i): sieve[j] = False
    primes = []; exps = []
    for i in range(2, n_fact+1):
        if not sieve[i]: continue
        primes.append(i); e = 0; t = n_fact
        while t > 0: t //= i; e += t
        exps.append(e)

    pc = len(primes)
    log_p = [math.log(p) for p in primes]
    log_n = sum(e*lp for e, lp in zip(exps, log_p))
    rem_log = [0.0]*(pc+1)
    for i in range(pc-1, -1, -1): rem_log[i] = rem_log[i+1] + exps[i]*log_p[i]

    # Precompute powers
    pow_val = [[1]*(exps[i]+1) for i in range(pc)]
    for i in range(pc):
        for e in range(1, exps[i]+1): pow_val[i][e] = pow_val[i][e-1] * primes[i]

    n_value = 1
    for i in range(pc): n_value *= pow_val[i][exps[i]]

    # Splits sorted by balance
    splits_by = {}
    for e in range(max(exps)+1):
        sp = []
        for x in range(e+1):
            for y in range(e+1-x):
                z = e-x-y; sp.append((x,y,z))
        sp.sort(key=lambda s: (max(s)-min(s), abs(3*s[0]-e)+abs(3*s[1]-e)+abs(3*s[2]-e)))
        splits_by[e] = sp

    def lb_range(l1, l2, l3, rem):
        t = l2 - l1
        if rem <= t: mv = l1+rem
        else:
            rem -= t; need = 2*(l3-l2)
            mv = l2+rem/2 if rem <= need else l3
        return l3 - mv

    best_lr = [1e100]; best_sum = [0]; best_abc = [None]

    def update(bins):
        vals = sorted(1 if all(b[i]==0 for i in range(pc)) else
                      math.prod(pow_val[i][b[i]] for i in range(pc)) for b in bins)
        s = sum(vals)
        lr = math.log(vals[2]) - math.log(vals[0]) if vals[0] > 0 else 1e100
        if lr < best_lr[0] - 1e-18 or (abs(lr-best_lr[0]) < 1e-18 and (best_sum[0] == 0 or s < best_sum[0])):
            best_lr[0] = lr; best_sum[0] = s; best_abc[0] = vals[:]

    # Greedy seed
    logs = [0.0, 0.0, 0.0]; bins = [[0]*pc, [0]*pc, [0]*pc]
    desc = sorted(range(pc), key=lambda i: -primes[i])
    for pi in desc:
        for _ in range(exps[pi]):
            t = min(range(3), key=lambda j: logs[j])
            logs[t] += log_p[pi]; bins[t][pi] += 1
    # Sort
    order = sorted(range(3), key=lambda i: logs[i])
    logs = [logs[i] for i in order]; bins = [bins[i] for i in order]
    best_lr[0] = logs[2] - logs[0]; update(bins)

    # Phase A: DFS
    nodes = [0]; cap = 3000000
    def phase_a(idx, logs, bins):
        if nodes[0] >= cap: return
        nodes[0] += 1
        if idx == pc: update(bins); return
        if lb_range(logs[0], logs[1], logs[2], rem_log[idx]) >= best_lr[0] - 1e-18: return
        for s in splits_by[exps[idx]]:
            nl = [logs[j] + s[j]*log_p[idx] for j in range(3)]
            nb = [bins[j][:] for j in range(3)]
            for j in range(3): nb[j][idx] += s[j]
            order = sorted(range(3), key=lambda j: nl[j])
            nl = [nl[j] for j in order]; nb = [nb[j] for j in order]
            phase_a(idx+1, nl, nb)

    phase_a(0, [0.0,0.0,0.0], [[0]*pc,[0]*pc,[0]*pc])

    # Phase B: enumerate a
    cr = int(round(n_value ** (1/3)))
    while (cr+1)**3 <= n_value: cr += 1
    while cr**3 > n_value: cr -= 1
    a_high = cr

    def max_div_leq(residual, limit):
        best_b = [0]
        desc2 = sorted(range(pc), key=lambda i: -primes[i])
        rem_max = [1]*(pc+1)
        for i in range(pc-1, -1, -1):
            pi = desc2[i]; rem_max[i] = rem_max[i+1] * pow_val[pi][residual[pi]]
        def dfs(idx, cur):
            if cur > limit: return
            if idx == pc:
                if cur > best_b[0]: best_b[0] = cur; return
            else:
                if cur * rem_max[idx] <= best_b[0]: return
                pi = desc2[idx]; pw = pow_val[pi]
                lim = limit // cur; e = residual[pi]
                kmax = e
                while kmax > 0 and pw[kmax] > lim: kmax -= 1
                for k in range(kmax, -1, -1):
                    nxt = cur * pw[k]
                    if nxt * rem_max[idx+1] <= best_b[0]: break
                    dfs(idx+1, nxt)
        dfs(0, 1)
        return best_b[0]

    def isqrt_big(x):
        if x <= 0: return 0
        r = int(math.isqrt(x))
        while (r+1)*(r+1) <= x: r += 1
        while r*r > x: r -= 1
        return r

    def eval_cand(a, a_exp):
        qa = n_value // a
        u = isqrt_big(qa)
        res = [exps[i] - a_exp[i] for i in range(pc)]
        b = max_div_leq(res, u)
        if b < a: return
        c = qa // b
        if c < b: return
        la = math.log(a); lc = math.log(c)
        lr = lc - la
        if lr < best_lr[0] - 1e-18 or (abs(lr - best_lr[0]) < 1e-18 and a + b + c < best_sum[0]):
            best_lr[0] = lr; best_sum[0] = a + b + c; best_abc[0] = [a, b, c]

    # Lower bound for a
    base_log_low = (log_n - 2*best_lr[0]) / 3.0
    log_high = math.log(a_high)

    def phase_b(idx, cur, lv, a_exp):
        dyn_low = max(base_log_low, (log_n - 2*best_lr[0])/3.0)
        if idx == pc:
            if lv + 1e-18 >= dyn_low and cur <= a_high:
                eval_cand(cur, a_exp)
            return
        if lv > log_high + 1e-18: return
        if lv + rem_log[idx] < dyn_low - 1e-18: return
        mul = 1
        for k in range(exps[idx]+1):
            if cur > a_high // mul: break
            a_exp[idx] = k
            phase_b(idx+1, cur*mul, lv + k*log_p[idx], a_exp)
            if k < exps[idx]: mul *= primes[idx]
        a_exp[idx] = 0

    phase_b(0, 1, 0.0, [0]*pc)
    return str(best_sum[0])

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

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.Comparator;
import java.util.List;

public class Euler418 {

    static class Triple {
        BigInteger a = BigInteger.ZERO;
        BigInteger b = BigInteger.ZERO;
        BigInteger c = BigInteger.ZERO;
        BigInteger sum = BigInteger.ZERO;
    }

    static BigInteger cube(long x) {
        BigInteger t = BigInteger.valueOf(x);
        return t.multiply(t).multiply(t);
    }

    static long sqrtFloorLong(BigInteger x) {
        if (x.compareTo(BigInteger.ZERO) <= 0)
            return 0;
        BigInteger a = BigInteger.ONE.shiftLeft(x.bitLength() / 2);
        boolean decrease = false;
        for (;;) {
            BigInteger b = x.divide(a).add(a).shiftRight(1);
            if (a.compareTo(b) == 0 || (a.compareTo(b) < 0 && decrease)) {
                return a.longValue();
            }
            decrease = a.compareTo(b) > 0;
            a = b;
        }
    }

    static long multiplyClamp(long a, long b) {
        if (a == 0 || b == 0)
            return 0;
        long max = Long.MAX_VALUE / a;
        if (b > max)
            return Long.MAX_VALUE;
        return a * b;
    }

    static double lowerBoundRange(double l1, double l2, double l3, double rem) {
        double t = l2 - l1;
        double minValue;
        if (rem <= t) {
            minValue = l1 + rem;
        } else {
            rem -= t;
            double need = 2.0 * (l3 - l2);
            if (rem <= need) {
                minValue = l2 + rem / 2.0;
            } else {
                minValue = l3;
            }
        }
        return l3 - minValue;
    }

    static class Solver {
        List<Long> primes;
        List<Integer> exps;
        int pcount;

        BigInteger nValue = BigInteger.ONE;
        double logN = 0.0;
        double[] logP;
        double[] remLog;

        long[][] powLong;
        List<List<int[]>> splitByExp;
        int[] descIdx;

        long phaseANodes = 0;
        long phaseACap = 3000000;
        double bestLogRange = 1e100;
        Triple best = new Triple();

        Solver(List<Long> primes, List<Integer> exps) {
            this.primes = primes;
            this.exps = exps;
            this.pcount = primes.size();

            logP = new double[pcount];
            remLog = new double[pcount + 1];
            powLong = new long[pcount][];
            descIdx = new int[pcount];

            Integer[] idx = new Integer[pcount];
            for (int i = 0; i < pcount; i++)
                idx[i] = i;
            Arrays.sort(idx, (i, j) -> Long.compare(primes.get(j), primes.get(i)));
            for (int i = 0; i < pcount; i++)
                descIdx[i] = idx[i];

            int maxE = 0;
            for (int i = 0; i < pcount; i++) {
                logP[i] = Math.log(primes.get(i));
                logN += exps.get(i) * logP[i];
                maxE = Math.max(maxE, exps.get(i));

                powLong[i] = new long[exps.get(i) + 1];
                powLong[i][0] = 1;
                for (int e = 1; e <= exps.get(i); e++) {
                    powLong[i][e] = powLong[i][e - 1] * primes.get(i);
                }

                BigInteger piPow = BigInteger.valueOf(primes.get(i)).pow(exps.get(i));
                nValue = nValue.multiply(piPow);
            }

            for (int i = pcount - 1; i >= 0; i--) {
                remLog[i] = remLog[i + 1] + exps.get(i) * logP[i];
            }

            splitByExp = new ArrayList<>(maxE + 1);
            for (int e = 0; e <= maxE; e++) {
                List<int[]> v = new ArrayList<>();
                for (int x = 0; x <= e; x++) {
                    for (int y = 0; y <= e - x; y++) {
                        int z = e - x - y;
                        v.add(new int[] { x, y, z });
                    }
                }

                final int currentE = e;
                v.sort((a, b) -> {
                    int minA = Math.min(a[0], Math.min(a[1], a[2]));
                    int maxA = Math.max(a[0], Math.max(a[1], a[2]));
                    int minB = Math.min(b[0], Math.min(b[1], b[2]));
                    int maxB = Math.max(b[0], Math.max(b[1], b[2]));
                    int sa = maxA - minA;
                    int sb = maxB - minB;
                    if (sa != sb)
                        return Integer.compare(sa, sb);

                    int ma = Math.abs(3 * a[0] - currentE) + Math.abs(3 * a[1] - currentE)
                            + Math.abs(3 * a[2] - currentE);
                    int mb = Math.abs(3 * b[0] - currentE) + Math.abs(3 * b[1] - currentE)
                            + Math.abs(3 * b[2] - currentE);
                    return Integer.compare(ma, mb);
                });
                splitByExp.add(v);
            }
        }

        void sortBins(double[] logs, int[][] bins) {
            Integer[] ord = { 0, 1, 2 };
            Arrays.sort(ord, (i, j) -> Double.compare(logs[i], logs[j]));
            double[] nl = { logs[ord[0]], logs[ord[1]], logs[ord[2]] };
            int[][] nb = { bins[ord[0]], bins[ord[1]], bins[ord[2]] };
            for (int i = 0; i < 3; i++) {
                logs[i] = nl[i];
                bins[i] = nb[i];
            }
        }

        BigInteger valueFromBin(int[] binExp) {
            BigInteger v = BigInteger.ONE;
            for (int i = 0; i < pcount; i++) {
                int e = binExp[i];
                if (e > 0)
                    v = v.multiply(BigInteger.valueOf(powLong[i][e]));
            }
            return v;
        }

        void updateBestFromBins(double[] logs, int[][] bins) {
            double range = logs[2] - logs[0];
            if (range > bestLogRange + 1e-15)
                return;

            BigInteger[] v = { valueFromBin(bins[0]), valueFromBin(bins[1]), valueFromBin(bins[2]) };
            Arrays.sort(v);
            BigInteger sum = v[0].add(v[1]).add(v[2]);

            if (range < bestLogRange - 1e-15 || best.sum.equals(BigInteger.ZERO) || sum.compareTo(best.sum) < 0) {
                bestLogRange = range;
                best.a = v[0];
                best.b = v[1];
                best.c = v[2];
                best.sum = sum;
            }
        }

        void greedySeed() {
            double[] logs = new double[3];
            int[][] bins = new int[3][pcount];

            Integer[] idx = new Integer[pcount];
            for (int i = 0; i < pcount; i++)
                idx[i] = i;
            Arrays.sort(idx, (i, j) -> Long.compare(primes.get(j), primes.get(i)));

            for (int pi : idx) {
                for (int c = 0; c < exps.get(pi); c++) {
                    int target = 0;
                    if (logs[1] < logs[target])
                        target = 1;
                    if (logs[2] < logs[target])
                        target = 2;
                    logs[target] += logP[pi];
                    bins[target][pi]++;
                }
            }
            sortBins(logs, bins);
            bestLogRange = logs[2] - logs[0];
            updateBestFromBins(logs, bins);
        }

        void phaseADfs(int idx, double[] logs, int[][] bins) {
            if (phaseANodes++ >= phaseACap)
                return;
            if (idx == pcount) {
                updateBestFromBins(logs, bins);
                return;
            }
            if (lowerBoundRange(logs[0], logs[1], logs[2], remLog[idx]) >= bestLogRange - 1e-15) {
                return;
            }

            int e = exps.get(idx);
            List<int[]> splits = splitByExp.get(e);
            for (int[] s : splits) {
                double[] nl = logs.clone();
                int[][] nb = new int[][] { bins[0].clone(), bins[1].clone(), bins[2].clone() };

                nl[0] += s[0] * logP[idx];
                nl[1] += s[1] * logP[idx];
                nl[2] += s[2] * logP[idx];

                nb[0][idx] += s[0];
                nb[1][idx] += s[1];
                nb[2][idx] += s[2];

                sortBins(nl, nb);
                phaseADfs(idx + 1, nl, nb);
            }
        }

        long cubeRootFloor() {
            double approx = Math.exp(logN / 3.0);
            if (approx < 1.0)
                approx = 1.0;
            long x = (long) approx;
            while (cube(x + 1).compareTo(nValue) <= 0)
                x++;
            while (cube(x).compareTo(nValue) > 0)
                x--;
            return x;
        }

        long lowerABound(long high) {
            BigInteger rhs = nValue.multiply(best.a).multiply(best.a);
            BigInteger c2 = best.c.multiply(best.c);
            long lo = 1;
            long hi = high;
            while (lo < hi) {
                long mid = lo + (hi - lo) / 2;
                BigInteger lhs = c2.multiply(cube(mid));
                if (lhs.compareTo(rhs) >= 0) {
                    hi = mid;
                } else {
                    lo = mid + 1;
                }
            }
            return lo;
        }

        long maxDivisorLeq(int[] residual, long limit) {
            int[] eDesc = new int[pcount];
            for (int i = 0; i < pcount; i++)
                eDesc[i] = residual[descIdx[i]];

            long[] remMax = new long[pcount + 1];
            remMax[pcount] = 1;
            for (int i = pcount - 1; i >= 0; i--) {
                int pi = descIdx[i];
                remMax[i] = multiplyClamp(remMax[i + 1], powLong[pi][eDesc[i]]);
            }

            long[] bestB = new long[] { 0 };
            dfsDiv(0, 1L, limit, eDesc, remMax, bestB);
            return bestB[0];
        }

        void dfsDiv(int idx, long cur, long limit, int[] eDesc, long[] remMax, long[] bestB) {
            if (cur > limit)
                return;
            if (idx == pcount) {
                if (cur > bestB[0])
                    bestB[0] = cur;
                return;
            }
            if (multiplyClamp(cur, remMax[idx]) <= bestB[0])
                return;

            int pi = descIdx[idx];
            long[] pw = powLong[pi];
            int e = eDesc[idx];
            long lim = limit / cur;
            int kmax = e;
            while (kmax > 0 && pw[kmax] > lim)
                kmax--;

            for (int k = kmax; k >= 0; k--) {
                long nxt = cur * pw[k];
                if (multiplyClamp(nxt, remMax[idx + 1]) <= bestB[0])
                    break;
                dfsDiv(idx + 1, nxt, limit, eDesc, remMax, bestB);
            }
        }

        void evaluateCandidate(long a, int[] aExp) {
            BigInteger A = BigInteger.valueOf(a);
            BigInteger qa = nValue.divide(A);
            long u = sqrtFloorLong(qa);

            int[] residual = new int[pcount];
            for (int i = 0; i < pcount; i++) {
                residual[i] = exps.get(i) - aExp[i];
            }

            long b = maxDivisorLeq(residual, u);
            if (b < a)
                return;

            BigInteger B = BigInteger.valueOf(b);
            BigInteger c = qa.divide(B);
            if (c.compareTo(B) < 0)
                return;

            BigInteger lhs = c.multiply(best.a);
            BigInteger rhs = best.c.multiply(A);

            if (lhs.compareTo(rhs) < 0) {
                best.a = A;
                best.b = B;
                best.c = c;
                best.sum = A.add(B).add(c);
                // Math.log for biginteger approx
                bestLogRange = (c.bitLength() - A.bitLength()) * Math.log(2.0); // Safe enough approx for next dynamic
                                                                                // pruning!
                // Actually safer to compute precise log:
                double lC = c.doubleValue();
                double lA = A.doubleValue();
                // But doubleValue() may map to Infinity if > 10^308. But here limit is ~10^20
                // max.
                bestLogRange = Math.log(lC) - Math.log(lA);
            } else if (lhs.compareTo(rhs) == 0) {
                BigInteger s = A.add(B).add(c);
                if (s.compareTo(best.sum) < 0) {
                    best.a = A;
                    best.b = B;
                    best.c = c;
                    best.sum = s;
                }
            }
        }

        void phaseBSearch() {
            long aHigh = cubeRootFloor();
            long baseALow = lowerABound(aHigh);
            double baseLogLow = Math.log(baseALow);
            double logHigh = Math.log(aHigh);

            int[] aExp = new int[pcount];

            dfsBSearch(0, 1L, 0.0, aExp, baseLogLow, aHigh, logHigh);
        }

        void dfsBSearch(int idx, long cur, double lv, int[] aExp, double baseLogLow, long aHigh, double logHigh) {
            double dynamicLogLow = Math.max(baseLogLow, (logN - 2.0 * bestLogRange) / 3.0);
            if (idx == pcount) {
                if (lv + 1e-15 >= dynamicLogLow && cur <= aHigh) {
                    evaluateCandidate(cur, aExp);
                }
                return;
            }
            if (lv > logHigh + 1e-15)
                return;
            if (lv + remLog[idx] < dynamicLogLow - 1e-15)
                return;

            long p = primes.get(idx);
            int e = exps.get(idx);
            long mul = 1;
            for (int k = 0; k <= e; k++) {
                if (cur > aHigh / mul)
                    break;
                long nxt = cur * mul;
                aExp[idx] = k;
                dfsBSearch(idx + 1, nxt, lv + k * logP[idx], aExp, baseLogLow, aHigh, logHigh);
                if (k < e)
                    mul *= p;
            }
            aExp[idx] = 0;
        }

        BigInteger solve() {
            greedySeed();

            double[] logs = new double[3];
            int[][] bins = new int[3][pcount];
            phaseANodes = 0;
            phaseADfs(0, logs, bins);

            phaseBSearch();
            return best.sum;
        }
    }

    static Object[] factorizeFactorial(int n) {
        List<Long> primes = new ArrayList<>();
        List<Integer> exps = new ArrayList<>();
        boolean[] sieve = new boolean[n + 1];
        Arrays.fill(sieve, true);
        if (n >= 0)
            sieve[0] = false;
        if (n >= 1)
            sieve[1] = false;

        for (int i = 2; i <= n; i++) {
            if (!sieve[i])
                continue;
            for (int j = i * 2; j <= n; j += i)
                sieve[j] = false;
            primes.add((long) i);
            int e = 0;
            int t = n;
            while (t > 0) {
                t /= i;
                e += t;
            }
            exps.add(e);
        }
        return new Object[] { primes, exps };
    }

    static String solve() {
        Object[] factors = factorizeFactorial(43);
        @SuppressWarnings("unchecked")
        List<Long> primes = (List<Long>) factors[0];
        @SuppressWarnings("unchecked")
        List<Integer> exps = (List<Integer>) factors[1];

        Solver solver = new Solver(primes, exps);
        return solver.solve().toString();
    }

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