Problem 418: Factorisation Triples
View on Project EulerProject 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
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! = 2^{39} 3^{19} 5^9 7^6 11^3 13^3 17^2 19^2 \cdot 23 \cdot 29 \cdot 31 \cdot 37 \cdot 41 \cdot 43.$$
Every admissible triple is therefore equivalent to choosing, for each prime \(p\), three nonnegative exponents \(\alpha_p,\beta_p,\gamma_p\) such that
$$\alpha_p + \beta_p + \gamma_p = e_p,$$
and then setting
$$a = \prod_{p \le 43} p^{\alpha_p},\qquad b = \prod_{p \le 43} p^{\beta_p},\qquad c = \prod_{p \le 43} p^{\gamma_p}.$$
Step 2: The objective in logarithmic form
Because the logarithm is strictly increasing, minimizing \(c/a\) is equivalent to minimizing
$$\Delta = \log c - \log a.$$
This turns the multiplicative balancing problem into an additive one. If the three current factors have sorted logarithms
$$L_1 \le L_2 \le L_3,$$
then the current quality of the triple is exactly the log-range \(L_3 - L_1\). A perfectly balanced triple would make that range as small as possible.
Step 3: A greedy seed gives the first upper bound
Before the expensive search begins, the implementation constructs a valid triple greedily. Prime copies are processed from large to small, and each copy is assigned to the currently smallest logarithmic bucket. This is not guaranteed to be optimal, but it usually produces a fairly balanced triple very quickly.
That first triple provides an initial upper bound on \(\Delta\). Once such a bound exists, later search branches can be discarded as soon as they can no longer beat it.
Step 4: Phase A - branch-and-bound over exponent splits
For one prime power \(p^e\), the number of three-way exponent splits is
$$\binom{e+2}{2},$$
because we are counting nonnegative integer solutions of \(x+y+z=e\). Naively taking the Cartesian product over all primes would be enormous, so the implementation explores these choices with branch-and-bound.
The split candidates for a fixed exponent \(e\) are tried in an order that favors balanced distributions first. The priority is determined by a small spread
$$\max(x,y,z) - \min(x,y,z),$$
and then by closeness to \(e/3\), measured by
$$|3x-e| + |3y-e| + |3z-e|.$$
After each prime is assigned, the three partial factors are re-sorted by size. This is important: only the ordered triple matters in the end, so sorting removes symmetric duplicates and keeps the pruning formulas simple.
Step 5: Lower bound from the remaining logarithmic mass
Suppose the current sorted logarithms are \(L_1 \le L_2 \le L_3\), and the primes that have not been assigned yet contribute total logarithmic mass \(R\). Even if that remaining mass could be split continuously, the smallest possible final range is obtained by filling the smaller bins first. The resulting optimistic lower bound is
$$\operatorname{LB}(L_1,L_2,L_3,R) = \begin{cases} L_3 - (L_1 + R), & R \le L_2 - L_1, \\ L_3 - \left(L_2 + \frac{R-(L_2-L_1)}{2}\right), & L_2 - L_1 \lt R \le L_2 - L_1 + 2(L_3-L_2), \\ 0, & R > L_2 - L_1 + 2(L_3-L_2). \end{cases}$$
The meaning is straightforward. First use the remaining logarithmic mass to raise the smallest factor until it reaches the middle one. Then use any leftover mass to raise the two smaller factors together toward the largest one. If even that is not enough to beat the best known range, the entire branch is pruned.
This bound is safe because the real problem is more restrictive than the continuous relaxation: prime exponents are discrete and tied to specific primes, so the true achievable range can only be larger.
Step 6: Phase B - exact refinement over the smallest factor
Phase A already produces an excellent ratio \(R^\ast = c^\ast / a^\ast\). The second phase turns that ratio into a search window for the smallest factor \(a\).
From \(a \le b \le c\) and \(abc = n\), we always have
$$a \le n^{1/3}.$$
Now assume a candidate triple is at least as good as the current best, so \(c/a \le R^\ast\). Because \(b \le c\),
$$n = abc \le a c^2 \le a(R^\ast a)^2 = (R^\ast)^2 a^3,$$
which yields
$$a \ge \left(\frac{n}{(R^\ast)^2}\right)^{1/3}.$$
Therefore phase B only needs to inspect divisors \(a\) inside the narrow interval
$$\left(\frac{n}{(R^\ast)^2}\right)^{1/3} \le a \le n^{1/3}.$$
The implementation enumerates precisely those divisors by walking through the exponent choices for \(a\), again pruning with logarithmic lower and upper bounds.
Step 7: For fixed \(a\), the optimal \(b\) is the largest divisor below the square root
Once \(a\) is fixed, define
$$q = \frac{n}{a} = bc.$$
For this fixed \(a\), minimizing \(c/a\) is the same as minimizing \(c\), and since \(bc=q\), that is equivalent to maximizing \(b\). The ordering constraint \(b \le c\) becomes
$$b \le \sqrt{q}.$$
So the best possible companion factor is simply the largest divisor of \(q\) that does not exceed \(\sqrt{q}\). Once that \(b\) is found, the third factor is forced:
$$c = \frac{q}{b}.$$
The residual factorization of \(q\) is already known from the chosen exponents of \(a\), so the implementation searches that divisor space directly. It also uses an upper bound from the remaining prime powers to prune branches that cannot improve the current best \(b\).
Step 8: Exact comparison and a small checkpoint
Logs are used only for pruning. The final comparison between candidate triples is exact. To compare the current candidate \((a,b,c)\) with the best-known triple \((a^\ast,b^\ast,c^\ast)\), it is enough to compare
$$c a^\ast \quad \text{and} \quad c^\ast a.$$
If the ratios are equal, the tie is resolved by the smaller sum \(a+b+c\).
A simple checkpoint is \(n=165=3\cdot 5 \cdot 11\). The only ordered factor triple is \((3,5,11)\), so the optimal sum is \(19\). The implementations verify themselves on small checkpoints of this kind before moving to \(43!\).
How the Code Works
The C++, Python, and Java implementations all follow the same structure. First they compute the prime exponents of \(43!\) and precompute prime powers. Then they generate every three-way split for each exponent value and sort those splits so balanced choices are examined first. A greedy construction gives the initial admissible triple.
After that, phase A performs a depth-first branch-and-bound search in log-space, always keeping the three partial factors sorted. Phase B then enumerates only the admissible divisors for the smallest factor, reconstructs the residual factorization, and finds the largest allowed companion divisor below the square-root threshold. The final choice is made with exact integer arithmetic, so floating-point approximation never decides correctness.
Complexity Analysis
If \(n = \prod p^{e_p}\), phase A would naively examine
$$\prod_p \binom{e_p+2}{2}$$
branches, which is exponential in the prime-exponent structure. Phase B is also combinatorial: in the worst case it may inspect many divisors in the admissible interval for \(a\), and for each one it performs a recursive search for the best residual divisor \(b\).
So there is no simple polynomial worst-case bound in \(\log n\). The method is practical here because four ideas remove most of the search space: a strong greedy seed, balanced ordering of exponent splits, the continuous lower bound on the attainable log-range, and pruning inside the residual divisor search. Memory use remains modest and is dominated by the stored prime powers, split tables, and recursion state.
References
- Problem page: https://projecteuler.net/problem=418
- Legendre's formula for prime exponents in \(n!\): Wikipedia - Legendre's formula
- Branch and bound: Wikipedia - Branch and bound
- Prime factorization and divisor structure: Wikipedia - Prime factor
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());
}
}