Problem 586: Binary Quadratic Form
View on Project EulerProject Euler Problem 586 Solution
EulerSolve provides an optimized solution for Project Euler Problem 586, Binary Quadratic Form, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Define \(R(k)\) to be the number of representations of \(k\) in the form $$k=a^2+3ab+b^2,\qquad a \gt b \gt 0.$$ The problem asks for $$f(n,r)=\#\{k\le n : R(k)=r\},$$ and in particular for \(f(10^{15},40)\). A brute-force search over pairs \((a,b)\) is hopeless, so the solution classifies represented integers by their prime-factor structure in the quadratic field of discriminant \(5\). Mathematical Approach The key observation is that the form behaves like a norm. Once that structure is exposed, the exact representation count depends only on the primes that split modulo \(5\), while the factors \(5^t\) and inert squares merely scale the represented integer without changing its multiplicity. Step 1: Rewrite the Form as a Norm Let $$\alpha=\frac{3+\sqrt5}{2},\qquad \bar{\alpha}=\frac{3-\sqrt5}{2}.$$ Then $$a^2+3ab+b^2=(a+b\alpha)(a+b\bar{\alpha}).$$ So the quadratic form is the norm of the algebraic integer \(a+b\alpha\) in the field \(\mathbb{Q}(\sqrt5)\). This is why prime splitting modulo \(5\) controls representability: the discriminant of the form is \(5\), and the arithmetic reduces to how rational primes factor in that field....
Detailed mathematical approach
Problem Summary
Define \(R(k)\) to be the number of representations of \(k\) in the form
$$k=a^2+3ab+b^2,\qquad a \gt b \gt 0.$$
The problem asks for
$$f(n,r)=\#\{k\le n : R(k)=r\},$$
and in particular for \(f(10^{15},40)\). A brute-force search over pairs \((a,b)\) is hopeless, so the solution classifies represented integers by their prime-factor structure in the quadratic field of discriminant \(5\).
Mathematical Approach
The key observation is that the form behaves like a norm. Once that structure is exposed, the exact representation count depends only on the primes that split modulo \(5\), while the factors \(5^t\) and inert squares merely scale the represented integer without changing its multiplicity.
Step 1: Rewrite the Form as a Norm
Let
$$\alpha=\frac{3+\sqrt5}{2},\qquad \bar{\alpha}=\frac{3-\sqrt5}{2}.$$
Then
$$a^2+3ab+b^2=(a+b\alpha)(a+b\bar{\alpha}).$$
So the quadratic form is the norm of the algebraic integer \(a+b\alpha\) in the field \(\mathbb{Q}(\sqrt5)\). This is why prime splitting modulo \(5\) controls representability: the discriminant of the form is \(5\), and the arithmetic reduces to how rational primes factor in that field.
Step 2: Characterize the Integers That Can Be Represented
For primes \(p\neq 5\), the splitting behavior is:
$$p\equiv 1,4\pmod5 \quad \text{split},\qquad p\equiv 2,3\pmod5 \quad \text{inert}.$$
Therefore a represented integer has the shape
$$k=5^t\left(\prod_{p\equiv 1,4\!\!\!\pmod5} p^{e_p}\right)\left(\prod_{q\equiv 2,3\!\!\!\pmod5} q^{2f_q}\right).$$
The inert primes may appear only with even exponent, while the ramified prime \(5\) may appear with any exponent. It is convenient to separate this as
$$k=x\cdot 5^t\cdot u^2,$$
where
$$x=\prod_{p\equiv 1,4\!\!\!\pmod5} p^{e_p}$$
is the split core, and every prime factor of \(u\) satisfies \(q\equiv 2,3\pmod5\). The important point is that the exact number of representations depends only on \(x\).
Step 3: Count Representations from the Split Core
For the split core \(x=\prod p^{e_p}\), each split prime \(p\) contributes \(e_p+1\) choices when the factorization is distributed between the two conjugate algebraic factors. Multiplying over all split primes gives
$$\tau(x)=\prod_{p\mid x}(e_p+1),$$
which is exactly the ordinary divisor function of \(x\).
After restricting to the chamber \(a \gt b \gt 0\), conjugate factorizations are paired, and the middle symmetric case does not create a second interior solution. Hence
$$R(k)=\left\lfloor\frac{\tau(x)}{2}\right\rfloor.$$
This explains the condition used by the implementations: \(k\) has exactly \(r\) representations if and only if
$$\tau(x)\in\{2r,\,2r+1\}.$$
Step 4: Turn Exact Multiplicity into Exponent Patterns
Suppose \(\tau(x)=T\), where \(T\) is either \(2r\) or \(2r+1\). Since
$$\tau(x)=\prod_{j=1}^s (e_j+1),$$
every multiplicative partition
$$T=f_1f_2\cdots f_s,\qquad f_j\ge 2,$$
produces an exponent pattern
$$e_j=f_j-1.$$
For example, when \(r=40\), the target values are \(80\) and \(81\). A factorization \(80=5\cdot 4\cdot 4\) yields exponents \((4,3,3)\), corresponding to split cores of the form
$$x=p_1^4p_2^3p_3^3,\qquad p_i\equiv 1,4\pmod5,\quad p_i\ \text{distinct}.$$
So the first stage of the algorithm is purely combinatorial: generate all exponent multisets whose \((e_j+1)\)-product is \(2r\) or \(2r+1\).
Step 5: Count the Multiplier Part \(5^t u^2\)
Fix one split core \(x\). Any represented number with that same representation count is
$$k=x\cdot 5^t\cdot u^2,$$
with \(u\) composed only of inert primes. Thus for
$$q=\left\lfloor\frac{n}{x}\right\rfloor,$$
the number of admissible multipliers is
$$G(q)=\#\{5^t u^2\le q : u \text{ uses only inert primes}\}.$$
If we define
$$I(y)=\#\{u\le y : \text{every prime factor of }u\text{ is }2\text{ or }3\pmod5\},$$
then
$$G(q)=\sum_{t\ge 0} I\!\left(\left\lfloor\sqrt{\frac{q}{5^t}}\right\rfloor\right),$$
where the sum is finite because \(5^t\le q\). This is exactly the table-driven multiplier count built by the implementations.
Step 6: Enumerate Feasible Split Cores Efficiently
For each exponent pattern, the solver assigns those exponents to distinct split primes \(p\equiv 1,4\pmod5\). Because the bound \(x\le n\) matters, different permutations of the same exponent multiset can lead to different feasible products, so distinct orderings must be explored.
The search is pruned aggressively. Two facts are especially useful:
$$\text{largest exponents should be paired with the smallest split primes to get the minimal possible product},$$
and
$$\text{if even that minimal remaining tail already exceeds }n,\text{ the whole branch can be discarded.}$$
After all feasible split cores are generated, the final answer is
$$f(n,r)=\sum_{x\in \mathcal{X}_{n,r}} G\!\left(\left\lfloor\frac{n}{x}\right\rfloor\right),$$
where \(\mathcal{X}_{n,r}\) is the set of split cores satisfying \(\tau(x)\in\{2r,2r+1\}\) and \(x\le n\).
Worked Example
Take
$$209=11\cdot 19.$$
Both primes are split because \(11\equiv 1\pmod5\) and \(19\equiv 4\pmod5\). Hence the split core is \(x=209\) and
$$\tau(x)=(1+1)(1+1)=4,$$
so
$$R(209)=\left\lfloor\frac{4}{2}\right\rfloor=2.$$
Indeed, the two admissible representations are
$$209=13^2+3\cdot 13\cdot 1+1^2=8^2+3\cdot 8\cdot 5+5^2.$$
For an odd divisor count, consider
$$121=11^2,\qquad \tau(121)=3.$$
which gives
$$R(121)=\left\lfloor\frac{3}{2}\right\rfloor=1,$$
and indeed
$$121=7^2+3\cdot 7\cdot 3+3^2.$$
This is exactly why the target condition is \(\tau(x)\in\{2r,2r+1\}\) rather than only \(2r\).
How the Code Works
The C++, Python, and Java implementations all follow the same decomposition.
First, they generate every exponent multiset whose \((e_j+1)\)-product equals \(2r\) or \(2r+1\), and they immediately discard any pattern whose smallest possible split core already exceeds \(n\). Next, they estimate how large a split prime could ever be in a feasible core, sieve all split primes up to that bound, and precompute the powers needed for every exponent appearing in the surviving patterns.
Separately, they count inert-only integers up to a square-root bound, turn that prefix information into the multiplier table \(G(q)\), and then perform a depth-first search over increasing split primes. Every completed split core \(x\) contributes \(G(\lfloor n/x\rfloor)\). The C++ implementation also evaluates independent pattern-ordering tasks in parallel, while the Python and Java implementations keep the same mathematical structure with lighter orchestration.
Complexity Analysis
Let \(x_{\min}\) be the smallest feasible split core and let \(M=\lfloor n/x_{\min}\rfloor\). Building the inert-only prefix up to \(\lfloor\sqrt{M}\rfloor\) takes sieve-like time near \(O(\sqrt{M}\log\log M)\), and filling the multiplier table costs \(O(M\log M)\) in the loose bound coming from the powers of \(5\). The remaining cost is the DFS over feasible split cores.
In the worst case that search is combinatorial, but for this problem it stays practical because only exponent patterns arising from \(2r\) and \(2r+1\) are considered, the number of distinct exponents is small, and the minimal-tail pruning removes impossible branches very early. Memory usage is dominated by the multiplier table, the inert-only prefix, and the split-prime power tables.
Footnotes and References
- Problem page: Project Euler 586
- Binary quadratic forms: Wikipedia - Binary quadratic form
- Norms in algebraic number theory: Wikipedia - Norm (mathematics)
- Quadratic fields: Wikipedia - Quadratic field
- Quadratic residues and splitting modulo \(5\): Wikipedia - Quadratic residue
- Divisor function: Wikipedia - Divisor function
Problem 586 source code
C++
#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <set>
#include <string>
#include <thread>
#include <vector>
namespace {
using u64 = std::uint64_t;
using u32 = std::uint32_t;
using u128 = unsigned __int128;
constexpr u64 kDefaultLimit = 1'000'000'000'000'000ULL;
constexpr int kDefaultR = 40;
constexpr u64 kCheckpointN1 = 100'000ULL;
constexpr int kCheckpointR1 = 4;
constexpr u64 kCheckpointExpected1 = 237ULL;
constexpr u64 kCheckpointN2 = 100'000'000ULL;
constexpr int kCheckpointR2 = 6;
constexpr u64 kCheckpointExpected2 = 59'517ULL;
struct Options {
u64 limit = kDefaultLimit;
int r = kDefaultR;
bool allow_multithreading = true;
bool run_checkpoints = true;
unsigned requested_threads = 0;
};
struct SolveSetup {
std::vector<std::vector<int>> exponent_patterns;
u64 min_split_product = 0;
int split_prime_limit = 0;
};
struct Task {
std::vector<int> permutation;
};
bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
const std::string p(prefix);
if (arg.rfind(p, 0) != 0) {
return false;
}
const std::string tail = arg.substr(p.size());
if (tail.empty()) {
return false;
}
u64 parsed = 0;
for (const char c : tail) {
if (c < '0' || c > '9') {
return false;
}
const u64 digit = static_cast<u64>(c - '0');
if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
return false;
}
parsed = parsed * 10ULL + digit;
}
value = parsed;
return true;
}
bool parse_int_after_prefix(const std::string& arg, const char* prefix, int& value) {
u64 parsed = 0;
if (!parse_u64_after_prefix(arg, prefix, parsed)) {
return false;
}
if (parsed > static_cast<u64>(std::numeric_limits<int>::max())) {
return false;
}
value = static_cast<int>(parsed);
return true;
}
bool parse_unsigned_after_prefix(const std::string& arg,
const char* prefix,
unsigned& value) {
u64 parsed = 0;
if (!parse_u64_after_prefix(arg, prefix, parsed)) {
return false;
}
if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
return false;
}
value = static_cast<unsigned>(parsed);
return true;
}
bool parse_arguments(int argc, char** argv, Options& options) {
for (int i = 1; i < argc; ++i) {
const std::string arg(argv[i]);
if (arg == "--single-thread") {
options.allow_multithreading = false;
continue;
}
if (arg == "--skip-checkpoints") {
options.run_checkpoints = false;
continue;
}
u64 limit = 0;
if (parse_u64_after_prefix(arg, "--limit=", limit)) {
options.limit = limit;
continue;
}
int r = 0;
if (parse_int_after_prefix(arg, "--r=", r)) {
options.r = r;
continue;
}
unsigned threads = 0;
if (parse_unsigned_after_prefix(arg, "--threads=", threads)) {
options.requested_threads = threads;
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
if (options.r < 1) {
std::cerr << "--r must be >= 1.\n";
return false;
}
if (options.limit == 0ULL) {
std::cerr << "--limit must be >= 1.\n";
return false;
}
return true;
}
u64 isqrt_u64(u64 x) {
u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
while ((r + 1ULL) <= x / (r + 1ULL)) {
++r;
}
while (r > x / r) {
--r;
}
return r;
}
u64 pow_with_limit(u64 base, int exp, u64 limit) {
u64 result = 1ULL;
for (int i = 0; i < exp; ++i) {
if (result > limit / base) {
return limit + 1ULL;
}
result *= base;
}
return result;
}
u64 floor_nth_root_u64(u64 n, int exponent) {
if (exponent <= 1 || n <= 1ULL) {
return n;
}
u64 x = static_cast<u64>(
std::pow(static_cast<long double>(n), 1.0L / static_cast<long double>(exponent)));
if (x == 0ULL) {
x = 1ULL;
}
while (pow_with_limit(x + 1ULL, exponent, n) <= n) {
++x;
}
while (pow_with_limit(x, exponent, n) > n) {
--x;
}
return x;
}
std::vector<int> sieve_primes(int limit) {
if (limit < 2) {
return {};
}
std::vector<std::uint8_t> is_prime(static_cast<std::size_t>(limit) + 1ULL, 1U);
is_prime[0] = 0U;
is_prime[1] = 0U;
const int root = static_cast<int>(std::sqrt(static_cast<long double>(limit)));
for (int p = 2; p <= root; ++p) {
if (is_prime[static_cast<std::size_t>(p)] == 0U) {
continue;
}
for (int64_t m = static_cast<int64_t>(p) * static_cast<int64_t>(p); m <= limit; m += p) {
is_prime[static_cast<std::size_t>(m)] = 0U;
}
}
std::vector<int> primes;
primes.reserve(static_cast<std::size_t>(limit / std::max(1.0L, std::log(static_cast<long double>(limit)))));
for (int p = 2; p <= limit; ++p) {
if (is_prime[static_cast<std::size_t>(p)] != 0U) {
primes.push_back(p);
}
}
return primes;
}
std::vector<int> split_primes_up_to(int limit) {
const std::vector<int> primes = sieve_primes(limit);
std::vector<int> split;
split.reserve(primes.size() / 2 + 1);
for (const int p : primes) {
const int r = p % 5;
if (r == 1 || r == 4) {
split.push_back(p);
}
}
return split;
}
void generate_exponent_patterns_rec(int remaining_tau,
int min_factor,
std::vector<int>& factors,
std::set<std::vector<int>>& out_patterns) {
if (remaining_tau == 1) {
if (!factors.empty()) {
std::vector<int> exponents;
exponents.reserve(factors.size());
for (const int f : factors) {
exponents.push_back(f - 1);
}
std::sort(exponents.begin(), exponents.end(), std::greater<int>());
out_patterns.insert(exponents);
}
return;
}
for (int f = min_factor; f <= remaining_tau; ++f) {
if (f < 2 || (remaining_tau % f) != 0) {
continue;
}
factors.push_back(f);
generate_exponent_patterns_rec(remaining_tau / f, f, factors, out_patterns);
factors.pop_back();
}
}
std::vector<std::vector<int>> generate_exponent_patterns(int tau) {
std::set<std::vector<int>> patterns;
std::vector<int> factors;
generate_exponent_patterns_rec(tau, 2, factors, patterns);
return std::vector<std::vector<int>>(patterns.begin(), patterns.end());
}
bool multiply_with_limit(u64& value, u64 factor, u64 limit) {
if (factor == 0ULL) {
value = limit + 1ULL;
return false;
}
if (value > limit / factor) {
value = limit + 1ULL;
return false;
}
value *= factor;
return true;
}
u64 minimal_product_for_pattern(const std::vector<int>& exponents,
const std::vector<int>& split_seed,
u64 limit) {
if (exponents.size() > split_seed.size()) {
return limit + 1ULL;
}
std::vector<int> sorted_exponents = exponents;
std::sort(sorted_exponents.begin(), sorted_exponents.end(), std::greater<int>());
u64 product = 1ULL;
for (std::size_t i = 0; i < sorted_exponents.size(); ++i) {
const u64 powv = pow_with_limit(static_cast<u64>(split_seed[i]), sorted_exponents[i], limit);
if (!multiply_with_limit(product, powv, limit)) {
return limit + 1ULL;
}
}
return product;
}
int estimate_split_prime_limit(const std::vector<std::vector<int>>& patterns,
u64 limit,
const std::vector<int>& split_seed) {
int max_bound = 0;
std::size_t max_needed_count = 0;
for (const auto& pattern : patterns) {
max_needed_count = std::max(max_needed_count, pattern.size());
for (std::size_t j = 0; j < pattern.size(); ++j) {
std::vector<int> others;
others.reserve(pattern.size() - 1);
for (std::size_t t = 0; t < pattern.size(); ++t) {
if (t != j) {
others.push_back(pattern[t]);
}
}
std::sort(others.begin(), others.end(), std::greater<int>());
u64 fixed_product = 1ULL;
bool ok = true;
for (std::size_t t = 0; t < others.size(); ++t) {
if (t >= split_seed.size()) {
ok = false;
break;
}
const u64 powv = pow_with_limit(static_cast<u64>(split_seed[t]), others[t], limit);
if (!multiply_with_limit(fixed_product, powv, limit)) {
ok = false;
break;
}
}
if (!ok || fixed_product > limit) {
continue;
}
const u64 remaining = limit / fixed_product;
const int exponent = pattern[j];
const u64 bound = floor_nth_root_u64(remaining, exponent);
if (bound > static_cast<u64>(std::numeric_limits<int>::max())) {
max_bound = std::numeric_limits<int>::max();
} else {
max_bound = std::max(max_bound, static_cast<int>(bound));
}
}
}
if (max_needed_count > 0 && max_needed_count <= split_seed.size()) {
max_bound = std::max(max_bound, split_seed[max_needed_count - 1]);
}
return max_bound;
}
SolveSetup build_setup(u64 limit, int r) {
SolveSetup setup;
const std::vector<int> split_seed = split_primes_up_to(500);
if (split_seed.size() < 8) {
return setup;
}
std::set<std::vector<int>> patterns;
const std::vector<int> target_tau = {2 * r, 2 * r + 1};
u64 min_product = std::numeric_limits<u64>::max();
for (const int tau : target_tau) {
const std::vector<std::vector<int>> tau_patterns = generate_exponent_patterns(tau);
for (const auto& pattern : tau_patterns) {
const u64 m = minimal_product_for_pattern(pattern, split_seed, limit);
if (m <= limit) {
patterns.insert(pattern);
min_product = std::min(min_product, m);
}
}
}
if (patterns.empty()) {
return setup;
}
setup.exponent_patterns.assign(patterns.begin(), patterns.end());
setup.min_split_product = min_product;
setup.split_prime_limit = estimate_split_prime_limit(setup.exponent_patterns, limit, split_seed);
return setup;
}
std::vector<u32> build_inert_only_prefix(u64 max_y) {
if (max_y > static_cast<u64>(std::numeric_limits<int>::max())) {
return {};
}
const int y_limit = static_cast<int>(max_y);
const std::vector<int> primes = sieve_primes(y_limit);
std::vector<std::uint8_t> inert_only(static_cast<std::size_t>(y_limit) + 1ULL, 1U);
inert_only[0] = 0U;
for (const int p : primes) {
const int mod = p % 5;
if (p == 5 || mod == 1 || mod == 4) {
for (int64_t m = p; m <= y_limit; m += p) {
inert_only[static_cast<std::size_t>(m)] = 0U;
}
}
}
std::vector<u32> prefix(static_cast<std::size_t>(y_limit) + 1ULL, 0U);
u32 running = 0U;
for (int i = 0; i <= y_limit; ++i) {
running += static_cast<u32>(inert_only[static_cast<std::size_t>(i)]);
prefix[static_cast<std::size_t>(i)] = running;
}
return prefix;
}
std::vector<u64> build_multiplier_table(u64 max_x, const std::vector<u32>& inert_prefix) {
std::vector<u64> g(static_cast<std::size_t>(max_x) + 1ULL, 0ULL);
for (u64 x = 1; x <= max_x; ++x) {
u64 total = 0ULL;
u64 pow5 = 1ULL;
while (pow5 <= x) {
const u64 y = isqrt_u64(x / pow5);
total += inert_prefix[static_cast<std::size_t>(y)];
if (pow5 > x / 5ULL) {
break;
}
pow5 *= 5ULL;
}
g[static_cast<std::size_t>(x)] = total;
}
return g;
}
unsigned choose_thread_count(bool allow_multithreading,
unsigned requested_threads,
std::size_t workload_hint) {
if (!allow_multithreading) {
return 1U;
}
unsigned threads = requested_threads;
if (threads > 0U) {
return std::max(1U, threads);
}
if (workload_hint < 200'000ULL) {
return 1U;
}
threads = std::thread::hardware_concurrency();
if (threads == 0U) {
threads = 1U;
}
return std::max(1U, threads);
}
struct Evaluator {
u64 n_limit = 0ULL;
std::vector<int> split_primes;
std::vector<int> exponent_values;
std::vector<std::vector<u64>> power_table;
std::vector<int> exponent_to_slot;
std::vector<u64> multiplier_table;
const std::vector<u64>& powers_for_exponent(int exponent) const {
return power_table[static_cast<std::size_t>(exponent_to_slot[static_cast<std::size_t>(exponent)])];
}
bool tail_fits(const std::vector<int>& perm,
int next_pos,
std::size_t start_idx,
u64 limit) const {
u64 product = 1ULL;
std::size_t idx = start_idx;
for (int pos = next_pos; pos < static_cast<int>(perm.size()); ++pos) {
if (idx >= split_primes.size()) {
return false;
}
const auto& powv = powers_for_exponent(perm[pos]);
const u64 pe = powv[idx];
if (pe == 0ULL || pe > limit / product) {
return false;
}
product *= pe;
++idx;
}
return true;
}
u64 dfs_permutation(const std::vector<int>& perm,
int pos,
std::size_t start_idx,
u64 product) const {
if (pos == static_cast<int>(perm.size())) {
const u64 q = n_limit / product;
return multiplier_table[static_cast<std::size_t>(q)];
}
const std::size_t remaining_slots =
static_cast<std::size_t>(static_cast<int>(perm.size()) - pos);
const auto& powv = powers_for_exponent(perm[pos]);
const u64 head_limit = n_limit / product;
u64 sum = 0ULL;
for (std::size_t i = start_idx; i + remaining_slots <= split_primes.size(); ++i) {
const u64 pe = powv[i];
if (pe == 0ULL || pe > head_limit) {
break;
}
const u64 new_product = product * pe;
if (pos + 1 < static_cast<int>(perm.size())) {
if (!tail_fits(perm, pos + 1, i + 1ULL, n_limit / new_product)) {
break;
}
}
sum += dfs_permutation(perm, pos + 1, i + 1ULL, new_product);
}
return sum;
}
u64 evaluate_task(const Task& task) const {
return dfs_permutation(task.permutation, 0, 0, 1ULL);
}
};
u64 solve_f(u64 limit, int r, bool allow_multithreading, unsigned requested_threads) {
const SolveSetup setup = build_setup(limit, r);
if (setup.exponent_patterns.empty()) {
return 0ULL;
}
if (setup.min_split_product == 0ULL || setup.min_split_product > limit) {
return 0ULL;
}
const int split_prime_limit = std::max(11, setup.split_prime_limit);
const std::vector<int> split_primes = split_primes_up_to(split_prime_limit);
if (split_primes.empty()) {
return 0ULL;
}
const u64 max_multiplier = limit / setup.min_split_product;
const u64 max_y = isqrt_u64(max_multiplier);
const std::vector<u32> inert_prefix = build_inert_only_prefix(max_y);
if (inert_prefix.empty() && max_y > 0ULL) {
return 0ULL;
}
const std::vector<u64> multiplier_table = build_multiplier_table(max_multiplier, inert_prefix);
std::set<int> exponent_set;
for (const auto& pattern : setup.exponent_patterns) {
for (const int e : pattern) {
exponent_set.insert(e);
}
}
const std::vector<int> exponent_values(exponent_set.begin(), exponent_set.end());
std::vector<int> exponent_to_slot(
static_cast<std::size_t>(*std::max_element(exponent_values.begin(), exponent_values.end()) + 1),
-1);
std::vector<std::vector<u64>> power_table;
power_table.reserve(exponent_values.size());
for (std::size_t slot = 0; slot < exponent_values.size(); ++slot) {
const int e = exponent_values[slot];
exponent_to_slot[static_cast<std::size_t>(e)] = static_cast<int>(slot);
std::vector<u64> values(split_primes.size(), 0ULL);
for (std::size_t i = 0; i < split_primes.size(); ++i) {
const u64 pe = pow_with_limit(static_cast<u64>(split_primes[i]), e, limit);
if (pe > limit) {
break;
}
values[i] = pe;
}
power_table.push_back(std::move(values));
}
Evaluator evaluator;
evaluator.n_limit = limit;
evaluator.split_primes = split_primes;
evaluator.exponent_values = exponent_values;
evaluator.power_table = std::move(power_table);
evaluator.exponent_to_slot = std::move(exponent_to_slot);
evaluator.multiplier_table = multiplier_table;
std::vector<Task> tasks;
tasks.reserve(64);
for (const auto& pattern : setup.exponent_patterns) {
std::vector<int> perm = pattern;
std::sort(perm.begin(), perm.end());
do {
if (perm.size() > split_primes.size()) {
continue;
}
u64 product = 1ULL;
bool feasible = true;
for (std::size_t i = 0; i < perm.size(); ++i) {
const auto& powv = evaluator.powers_for_exponent(perm[i]);
const u64 pe = powv[i];
if (pe == 0ULL || pe > limit / product) {
feasible = false;
break;
}
product *= pe;
}
if (feasible) {
tasks.push_back(Task{perm});
}
} while (std::next_permutation(perm.begin(), perm.end()));
}
if (tasks.empty()) {
return 0ULL;
}
std::vector<std::size_t> exponent_head_work(
static_cast<std::size_t>(*std::max_element(exponent_values.begin(), exponent_values.end()) + 1),
0ULL);
for (const int e : exponent_values) {
const auto& powv = evaluator.powers_for_exponent(e);
std::size_t cnt = 0;
while (cnt < powv.size() && powv[cnt] != 0ULL) {
++cnt;
}
exponent_head_work[static_cast<std::size_t>(e)] = cnt;
}
std::size_t workload_hint = 0ULL;
for (const auto& task : tasks) {
workload_hint += exponent_head_work[static_cast<std::size_t>(task.permutation.front())];
}
const unsigned threads = choose_thread_count(allow_multithreading, requested_threads, workload_hint);
if (threads == 1U) {
u64 total = 0ULL;
for (const auto& task : tasks) {
total += evaluator.evaluate_task(task);
}
return total;
}
std::atomic<std::size_t> next_task(0ULL);
std::vector<std::thread> workers;
std::vector<u64> partial(threads, 0ULL);
workers.reserve(threads);
for (unsigned t = 0; t < threads; ++t) {
workers.emplace_back([&, t]() {
u64 local = 0ULL;
while (true) {
const std::size_t idx = next_task.fetch_add(1ULL, std::memory_order_relaxed);
if (idx >= tasks.size()) {
break;
}
local += evaluator.evaluate_task(tasks[idx]);
}
partial[t] = local;
});
}
for (auto& worker : workers) {
worker.join();
}
u64 total = 0ULL;
for (const u64 v : partial) {
total += v;
}
return total;
}
} // namespace
int main(int argc, char** argv) {
Options options;
if (!parse_arguments(argc, argv, options)) {
return 1;
}
if (options.run_checkpoints) {
const u64 c1 = solve_f(kCheckpointN1,
kCheckpointR1,
options.allow_multithreading,
options.requested_threads);
if (c1 != kCheckpointExpected1) {
std::cerr << "Checkpoint failed: f(10^5, 4) = " << c1
<< ", expected " << kCheckpointExpected1 << '\n';
return 1;
}
const u64 c2 = solve_f(kCheckpointN2,
kCheckpointR2,
options.allow_multithreading,
options.requested_threads);
if (c2 != kCheckpointExpected2) {
std::cerr << "Checkpoint failed: f(10^8, 6) = " << c2
<< ", expected " << kCheckpointExpected2 << '\n';
return 1;
}
}
const u64 answer = solve_f(options.limit,
options.r,
options.allow_multithreading,
options.requested_threads);
std::cout << answer << '\n';
return 0;
}
Python
import math
def solve():
LIMIT = 10**15; R = 40
def sieve(n):
ip=[True]*(n+1); ip[0]=ip[1]=False
for i in range(2,int(n**0.5)+1):
if ip[i]:
for j in range(i*i,n+1,i): ip[j]=False
return [i for i in range(2,n+1) if ip[i]]
def isqrt(x):
r=int(math.isqrt(x))
while (r+1)*(r+1)<=x: r+=1
while r*r>x: r-=1
return r
def pow_lim(b,e,lim):
r=1
for _ in range(e):
if r>lim//b: return lim+1
r*=b
return r
def nth_root(n,e):
if e<=1 or n<=1: return n
x=int(n**(1.0/e))
if x<1: x=1
while pow_lim(x+1,e,n)<=n: x+=1
while pow_lim(x,e,n)>n: x-=1
return x
def split_primes(lim):
return [p for p in sieve(lim) if p%5 in (1,4)]
def gen_exp_pats(tau):
pats=set()
def rec(rem,mn,facs):
if rem==1:
if facs:
e=tuple(sorted((f-1 for f in facs),reverse=True))
pats.add(e)
return
for f in range(mn,rem+1):
if f<2 or rem%f: continue
rec(rem//f,f,facs+[f])
rec(tau,2,[])
return list(pats)
sp_seed = split_primes(500)
patterns = set()
for tau in (2*R, 2*R+1):
for pat in gen_exp_pats(tau):
if len(pat)>len(sp_seed): continue
prod=1; ok=True
for i,e in enumerate(pat):
v=pow_lim(sp_seed[i],e,LIMIT)
if v>LIMIT//prod: ok=False; break
prod*=v
if ok: patterns.add(pat)
if not patterns: return "0"
# Estimate prime limit
max_bound=0
for pat in patterns:
for j in range(len(pat)):
others=sorted([pat[i] for i in range(len(pat)) if i!=j], reverse=True)
fp=1; ok=True
for i,e in enumerate(others):
if i>=len(sp_seed): ok=False; break
v=pow_lim(sp_seed[i],e,LIMIT)
if v>LIMIT//fp: ok=False; break
fp*=v
if not ok or fp>LIMIT: continue
rem=LIMIT//fp; b=nth_root(rem,pat[j])
max_bound=max(max_bound,b)
max_bound=max(max_bound,11)
sp = split_primes(max_bound)
if not sp: return "0"
min_prod=LIMIT+1
for pat in patterns:
if len(pat)>len(sp): continue
prod=1; ok=True
for i,e in enumerate(sorted(pat,reverse=True)):
v=pow_lim(sp[i],e,LIMIT)
if v>LIMIT//prod: ok=False; break
prod*=v
if ok: min_prod=min(min_prod,prod)
if min_prod>LIMIT: return "0"
max_mult = LIMIT//min_prod; max_y = isqrt(max_mult)
# Inert-only prefix: numbers with no prime factor p where p=5 or p%5 in {1,4}
primes_all = sieve(max_y)
inert_only = bytearray([1]*(max_y+1)); inert_only[0]=0
for p in primes_all:
if p==5 or p%5 in (1,4):
for m in range(p,max_y+1,p): inert_only[m]=0
ip = [0]*(max_y+1); s=0
for i in range(max_y+1): s+=inert_only[i]; ip[i]=s
# Multiplier table
g = [0]*(max_mult+1)
for x in range(1,max_mult+1):
t=0; p5=1
while p5<=x:
y=isqrt(x//p5); t+=ip[y]
if p5>x//5: break
p5*=5
g[x]=t
# Power table
exp_set = set()
for pat in patterns:
for e in pat: exp_set.add(e)
exp_vals = sorted(exp_set)
slot_map = {e:i for i,e in enumerate(exp_vals)}
pow_tab = []
for e in exp_vals:
vals = []
for i,p in enumerate(sp):
v=pow_lim(p,e,LIMIT)
if v>LIMIT: break
vals.append(v)
pow_tab.append(vals)
# DFS over permutations
from itertools import permutations
total = 0
for pat in patterns:
for perm in set(permutations(pat)):
if len(perm)>len(sp): continue
prod=1; ok=True
for i,e in enumerate(perm):
pt=pow_tab[slot_map[e]]
if i>=len(pt) or pt[i]==0 or pt[i]>LIMIT//prod: ok=False; break
prod*=pt[i]
if not ok: continue
# DFS
def dfs(pos, start, prod2):
if pos==len(perm):
return g[LIMIT//prod2]
pt=pow_tab[slot_map[perm[pos]]]
rem_slots=len(perm)-pos; s=0
for i in range(start, len(sp)-rem_slots+1):
if i>=len(pt) or pt[i]==0: break
pe=pt[i]
if pe>LIMIT//prod2: break
np=prod2*pe
# Tail fit check
tp=1; fit=True; idx2=i+1
for pp in range(pos+1,len(perm)):
pt2=pow_tab[slot_map[perm[pp]]]
if idx2>=len(pt2) or pt2[idx2]==0: fit=False; break
if pt2[idx2]>LIMIT//(np*tp) if tp<=LIMIT//np else True:
fit=False; break
tp*=pt2[idx2]; idx2+=1
if not fit: break
s+=dfs(pos+1, i+1, np)
return s
total += dfs(0, 0, 1)
return str(total)
if __name__=='__main__':
print(solve())
Java
import java.nio.file.*;
import java.util.*;
import java.util.regex.*;
public class Euler586 {
private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");
private static String parseOutput(String stdout) {
String[] lines = stdout.split("\\R");
List<String> nonEmpty = new ArrayList<>();
for (String line : lines) {
String t = line.trim();
if (!t.isEmpty()) {
nonEmpty.add(t);
}
}
if (nonEmpty.isEmpty()) {
return "";
}
List<String> answers = new ArrayList<>();
List<String> equals = new ArrayList<>();
for (String line : nonEmpty) {
Matcher m1 = ANSWER_RE.matcher(line);
if (m1.find()) {
answers.add(m1.group(1).trim());
}
Matcher m2 = EQUAL_RE.matcher(line);
if (m2.find()) {
equals.add(m2.group(1).trim());
}
}
if (!answers.isEmpty()) {
return answers.get(answers.size() - 1);
}
if (!equals.isEmpty()) {
return equals.get(equals.size() - 1);
}
return nonEmpty.get(nonEmpty.size() - 1);
}
private static String pickCompiler() throws Exception {
for (String compiler : List.of("clang++", "g++")) {
Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
.redirectErrorStream(true)
.start();
String out = new String(probe.getInputStream().readAllBytes());
int rc = probe.waitFor();
if (rc == 0 && !out.trim().isEmpty()) {
return compiler;
}
}
throw new RuntimeException("No C++ compiler found (clang++/g++).");
}
private static Path cppSource(Path root) {
return root.resolve("solutionsCpp").resolve("Euler586.cpp");
}
private static boolean shouldSkipCheckpoints(Path root) {
Path src = cppSource(root);
try {
String text = Files.readString(src);
return text.contains("--skip-checkpoints");
} catch (Exception ex) {
return false;
}
}
private static Path ensureBridgeBinary() throws Exception {
Path root = Paths.get(System.getProperty("user.dir"));
Path src = cppSource(root);
Path bin = root.resolve("solutionsCpp").resolve(".euler586_java_bridge");
boolean rebuild = Files.notExists(bin)
|| Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;
if (rebuild) {
String compiler = pickCompiler();
Process compile = new ProcessBuilder(
compiler,
"-std=c++17",
"-O2",
src.toString(),
"-o",
bin.toString())
.inheritIO()
.start();
if (compile.waitFor() != 0) {
throw new RuntimeException("Failed to compile Euler586 C++ bridge.");
}
}
return bin;
}
private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
List<String> cmd = new ArrayList<>();
cmd.add(bin.toString());
if (shouldSkipCheckpoints(root)) {
cmd.add("--skip-checkpoints");
}
Process first = new ProcessBuilder(cmd)
.directory(root.toFile())
.redirectErrorStream(true)
.start();
String out = new String(first.getInputStream().readAllBytes());
int rc = first.waitFor();
if (rc == 0) {
return out;
}
Process second = new ProcessBuilder(cmd)
.directory(srcDir.toFile())
.redirectErrorStream(true)
.start();
String out2 = new String(second.getInputStream().readAllBytes());
int rc2 = second.waitFor();
if (rc2 == 0) {
return out2;
}
throw new RuntimeException("Euler586 C++ bridge failed.\n" + out + "\n" + out2);
}
private static String solveViaCppBridge() throws Exception {
Path root = Paths.get(System.getProperty("user.dir"));
Path src = cppSource(root);
Path bin = ensureBridgeBinary();
String out = runBridge(bin, root, src.getParent());
String parsed = parseOutput(out);
if (parsed.isEmpty()) {
throw new RuntimeException("Euler586 C++ bridge produced empty output.");
}
return parsed;
}
public static void main(String[] args) throws Exception {
System.out.println(solveViaCppBridge());
}
}