Problem 829: Integral Fusion
View on Project EulerProject Euler Problem 829 Solution
EulerSolve provides an optimized solution for Project Euler Problem 829, Integral Fusion, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For each integer \(n\) with \(2\le n\le 31\), form the double factorial $$n!!=\prod_{\substack{1\le k\le n \\ k\equiv n\pmod{2}}} k.$$ The problem associates every positive integer with a canonical binary factor tree: primes are leaves, while a composite number is split at the divisor closest to \(\sqrt m\) from below and then treated recursively. If \(T(m)\) denotes that tree shape, we must find, for every \(n\), the smallest integer \(M(n)\ge 2\) with \(T(M(n))=T(n!!)\), and finally compute $$\sum_{n=2}^{31} M(n).$$ The challenge is that many different integers can share the same tree, so the solution generates integers by shape class instead of scanning all integers one by one. Mathematical Approach The canonical decomposition is deterministic, so each positive integer belongs to exactly one shape class. The solver therefore turns the problem into a sequence-generation task: for each target shape coming from \(n!!\), list all integers with that shape in increasing order and take the first one....
Detailed mathematical approach
Problem Summary
For each integer \(n\) with \(2\le n\le 31\), form the double factorial
$$n!!=\prod_{\substack{1\le k\le n \\ k\equiv n\pmod{2}}} k.$$
The problem associates every positive integer with a canonical binary factor tree: primes are leaves, while a composite number is split at the divisor closest to \(\sqrt m\) from below and then treated recursively. If \(T(m)\) denotes that tree shape, we must find, for every \(n\), the smallest integer \(M(n)\ge 2\) with \(T(M(n))=T(n!!)\), and finally compute
$$\sum_{n=2}^{31} M(n).$$
The challenge is that many different integers can share the same tree, so the solution generates integers by shape class instead of scanning all integers one by one.
Mathematical Approach
The canonical decomposition is deterministic, so each positive integer belongs to exactly one shape class. The solver therefore turns the problem into a sequence-generation task: for each target shape coming from \(n!!\), list all integers with that shape in increasing order and take the first one.
Step 1: Define the canonical shape
For a prime \(p\), set
$$T(p)=\bullet.$$
For a composite integer \(m\), define
$$\delta(m)=\max\{d:d\mid m,\ d\le \sqrt m\},\qquad a=\delta(m),\qquad b=\frac{m}{a},$$
so \(a\le b\), and then set
$$T(m)=\langle T(a),T(b)\rangle.$$
The factor \(a\) is the largest divisor not exceeding \(\sqrt m\), so \((a,b)\) is the most balanced divisor pair of \(m\). This makes the top split unique, and recursion makes the whole tree unique.
Step 2: Reformulate the target quantity
For each \(n\in\{2,3,\dots,31\}\), let
$$t_n=T(n!!),\qquad M(n)=\min\{m\ge 2:T(m)=t_n\}.$$
The minimum exists because \(n!!\) itself has shape \(t_n\). In particular,
$$M(n)\le n!!\le 31!!,$$
so every required answer lies below a finite global bound.
For each shape \(s\), define \(S_s\) to be the increasing sequence of all integers \(m\le 31!!\) with \(T(m)=s\). Then the desired value is simply the first term of \(S_{t_n}\).
Step 3: Build each shape from its children
The leaf shape is easy:
$$S_{\bullet}=(2,3,5,7,11,\dots),$$
because primes are exactly the integers whose canonical tree is a single leaf.
Now suppose \(s=\langle u,v\rangle\) is an internal shape. Any integer \(m\) with \(T(m)=s\) has a canonical top split \(m=xy\) with
$$T(x)=u,\qquad T(y)=v,\qquad x\le y.$$
Therefore every element of \(S_s\) must occur among candidate products of the form
$$x\in S_u,\qquad y\in S_v,\qquad x\le y,\qquad m=xy.$$
This gives a complete candidate source for the parent shape.
Step 4: Revalidate the candidate products
Not every such product really belongs to \(S_s\). Even if \(x\) has shape \(u\) and \(y\) has shape \(v\), the product \(xy\) might admit a more balanced divisor pair than \((x,y)\), so its canonical split can move to a different tree.
The correct characterization is therefore
$$S_s=\{xy:x\in S_u,\ y\in S_v,\ x\le y,\ T(xy)=s\}.$$
The extra test \(T(xy)=s\) is essential: it filters out products whose top split or deeper recursive structure changes after recomputing the canonical tree.
Step 5: Merge the product streams lazily
Write
$$S_u=(x_0,x_1,x_2,\dots),\qquad S_v=(y_0,y_1,y_2,\dots)$$
in increasing order. For fixed \(i\), only indices \(j\) with \(y_j\ge x_i\) are allowed, so each row
$$x_i y_j\qquad (j\ge j_0(i))$$
is itself increasing, where \(j_0(i)\) is the first index satisfying \(y_{j_0(i)}\ge x_i\).
A priority queue can merge these rows without generating the full Cartesian product. Whenever the smallest available product is removed, the algorithm advances only that row and opens a new row only when its first possible product can compete with the current minimum. Duplicate products are skipped, and every surviving product is rechecked against the canonical-shape rule.
Worked Example: \(n=9\)
We have
$$9!!=9\cdot 7\cdot 5\cdot 3\cdot 1=945.$$
The largest divisor of \(945\) not exceeding \(\sqrt{945}\) is \(27\), so the top split is
$$945=27\cdot 35.$$
Continuing recursively,
$$27=3\cdot 9,\qquad 9=3\cdot 3,\qquad 35=5\cdot 7.$$
Hence
$$T(945)=\langle \langle \bullet,\langle \bullet,\bullet\rangle\rangle,\langle \bullet,\bullet\rangle\rangle.$$
Now look at \(72\). Its largest divisor not exceeding \(\sqrt{72}\) is \(8\), so
$$72=8\cdot 9,\qquad 8=2\cdot 4,\qquad 4=2\cdot 2,\qquad 9=3\cdot 3.$$
This gives exactly the same tree shape, so \(T(72)=T(945)\). The increasing sequence for that shape begins with \(72\), and the implementations confirm the checkpoint
$$M(9)=72.$$
How the Code Works
The C++, Python, and Java implementations first compute the target shapes \(T(n!!)\) for \(2\le n\le 31\). They then mark only the shapes needed by those targets and by their descendants, which keeps the later sequence generation focused on the relevant part of the tree space.
To evaluate \(T(m)\), the implementation uses deterministic primality testing for 64-bit integers, Pollard's rho factorization for composites, and divisor generation from the prime factorization to identify the largest divisor not exceeding \(\sqrt m\). Results are memoized so repeated shape queries become cheap.
Each needed shape owns a lazy increasing sequence. The leaf sequence emits primes in order. Every internal sequence performs a best-first merge of child products with a priority queue, enforces the order \(x\le y\), ignores duplicates, rejects products above \(31!!\), and keeps only those products whose recomputed canonical tree matches the target shape. Once the first term of every target sequence has been obtained, the program adds those minima.
Complexity Analysis
The runtime is output-sensitive: it depends on how many candidate products must be examined before the required values for each shape appear. If a needed shape \(s\) inspects \(K_s\) candidate products, its heap work is \(O(K_s\log K_s)\) in the standard best-first merge model.
The expensive subroutine is canonical-shape validation, because that requires primality testing, factorization, divisor enumeration, and recursive tree reconstruction. Memoization removes much of the repeated cost in practice, so the method is efficient for the required range, but the cleanest description is “heap-driven lazy generation with number-theoretic validation” rather than a single closed-form complexity bound. Memory usage is dominated by cached factorizations, primality results, canonical shapes, and stored prefixes of the generated sequences.
Footnotes and References
- Problem page: https://projecteuler.net/problem=829
- Double factorial: Wikipedia — Double factorial
- Pollard's rho algorithm: Wikipedia — Pollard's rho algorithm
- Miller-Rabin primality test: Wikipedia — Miller-Rabin primality test
- Priority queue: Wikipedia — Priority queue
Problem 829 source code
C++
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <functional>
#include <iostream>
#include <map>
#include <numeric>
#include <optional>
#include <queue>
#include <random>
#include <set>
#include <unordered_map>
#include <unordered_set>
#include <vector>
using u64 = std::uint64_t;
using u128 = unsigned __int128;
struct ShapeInfo {
int left = -1;
int right = -1;
int leaves = 1;
};
struct HeapItem {
u64 prod;
int i;
int j;
};
struct HeapCmp {
bool operator()(const HeapItem& a, const HeapItem& b) const {
if (a.prod != b.prod) return a.prod > b.prod;
if (a.i != b.i) return a.i > b.i;
return a.j > b.j;
}
};
struct Sequence {
bool needed = false;
bool leaf = false;
bool started = false;
bool exhausted = false;
int left = -1;
int right = -1;
int next_i = 0;
u64 next_prime = 2;
std::vector<u64> values;
std::priority_queue<HeapItem, std::vector<HeapItem>, HeapCmp> heap;
std::unordered_set<u64> in_heap;
};
static std::mt19937_64 rng(0);
static std::vector<ShapeInfo> shapes(1); // id 0 = LEAF
static std::map<std::pair<int, int>, int> shape_id_by_children;
static std::unordered_map<u64, int> shape_cache;
static std::unordered_map<u64, u64> best_divisor_cache;
static std::unordered_map<u64, bool> prime_cache;
static std::unordered_map<u64, std::map<u64, int>> factor_cache;
static u64 MAXVAL = 0;
static std::vector<Sequence> seqs;
static u64 mul_mod(u64 a, u64 b, u64 mod) {
return static_cast<u64>((static_cast<u128>(a) * b) % mod);
}
static u64 pow_mod(u64 a, u64 e, u64 mod) {
u64 r = 1 % mod;
a %= mod;
while (e > 0) {
if (e & 1ULL) r = mul_mod(r, a, mod);
a = mul_mod(a, a, mod);
e >>= 1ULL;
}
return r;
}
static bool is_prime64(u64 n) {
auto it = prime_cache.find(n);
if (it != prime_cache.end()) return it->second;
if (n < 2) return prime_cache[n] = false;
for (u64 p : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
if (n % p == 0) return prime_cache[n] = (n == p);
}
u64 d = n - 1;
int s = 0;
while ((d & 1ULL) == 0) {
d >>= 1ULL;
++s;
}
for (u64 a : {2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL}) {
if (a % n == 0) continue;
u64 x = pow_mod(a, d, n);
if (x == 1 || x == n - 1) continue;
bool composite = true;
for (int r = 1; r < s; ++r) {
x = mul_mod(x, x, n);
if (x == n - 1) {
composite = false;
break;
}
}
if (composite) return prime_cache[n] = false;
}
return prime_cache[n] = true;
}
static u64 pollard_rho(u64 n) {
if ((n & 1ULL) == 0) return 2;
if (n % 3ULL == 0) return 3;
std::uniform_int_distribution<u64> dist(2, n - 2);
while (true) {
u64 c = dist(rng);
u64 x = dist(rng);
u64 y = x;
u64 d = 1;
while (d == 1) {
x = (mul_mod(x, x, n) + c) % n;
y = (mul_mod(y, y, n) + c) % n;
y = (mul_mod(y, y, n) + c) % n;
u64 diff = (x > y ? x - y : y - x);
d = std::gcd(diff, n);
}
if (d != n) return d;
}
}
static void factor_rec(u64 n, std::map<u64, int>& fac) {
if (n == 1) return;
if (is_prime64(n)) {
fac[n]++;
return;
}
u64 d = pollard_rho(n);
factor_rec(d, fac);
factor_rec(n / d, fac);
}
static std::map<u64, int> factorize(u64 n) {
auto it = factor_cache.find(n);
if (it != factor_cache.end()) return it->second;
std::map<u64, int> fac;
factor_rec(n, fac);
factor_cache[n] = fac;
return fac;
}
static void gen_divisors_rec(const std::vector<std::pair<u64, int>>& f, int idx, u64 cur, std::vector<u64>& out) {
if (idx == static_cast<int>(f.size())) {
out.push_back(cur);
return;
}
auto [p, e] = f[idx];
u64 v = 1;
for (int i = 0; i <= e; ++i) {
gen_divisors_rec(f, idx + 1, cur * v, out);
v *= p;
}
}
static u64 best_divisor_le_sqrt(u64 n) {
auto it = best_divisor_cache.find(n);
if (it != best_divisor_cache.end()) return it->second;
auto fac_map = factorize(n);
std::vector<std::pair<u64, int>> f(fac_map.begin(), fac_map.end());
std::vector<u64> divs;
divs.reserve(1024);
gen_divisors_rec(f, 0, 1, divs);
u64 root = static_cast<u64>(std::sqrt(static_cast<long double>(n)));
while ((root + 1) * (root + 1) <= n) ++root;
while (root * root > n) --root;
u64 best = 1;
for (u64 d : divs) {
if (d <= root && d > best) best = d;
}
best_divisor_cache[n] = best;
return best;
}
static int make_shape(int l, int r) {
auto key = std::make_pair(l, r);
auto it = shape_id_by_children.find(key);
if (it != shape_id_by_children.end()) return it->second;
ShapeInfo s;
s.left = l;
s.right = r;
s.leaves = shapes[l].leaves + shapes[r].leaves;
int id = static_cast<int>(shapes.size());
shapes.push_back(s);
shape_id_by_children[key] = id;
return id;
}
static int shape_of(u64 n) {
auto it = shape_cache.find(n);
if (it != shape_cache.end()) return it->second;
if (is_prime64(n)) {
shape_cache[n] = 0;
return 0;
}
u64 d = best_divisor_le_sqrt(n);
u64 a = d;
u64 b = n / d;
if (a > b) std::swap(a, b);
int l = shape_of(a);
int r = shape_of(b);
int id = make_shape(l, r);
shape_cache[n] = id;
return id;
}
static std::optional<u64> next_value(int sh);
static std::optional<u64> get_value_maybe(int sh, int idx) {
auto& seq = seqs[sh];
while (static_cast<int>(seq.values.size()) <= idx) {
auto nv = next_value(sh);
if (!nv.has_value()) return std::nullopt;
seq.values.push_back(*nv);
}
return seq.values[idx];
}
static int ensure_right_ge(int right_sh, u64 x) {
auto& rseq = seqs[right_sh];
if (rseq.exhausted && (rseq.values.empty() || rseq.values.back() < x)) return -1;
while (rseq.values.empty() || rseq.values.back() < x) {
auto nv = next_value(right_sh);
if (!nv.has_value()) {
rseq.exhausted = true;
break;
}
rseq.values.push_back(*nv);
}
if (rseq.values.empty() || rseq.values.back() < x) return -1;
auto it = std::lower_bound(rseq.values.begin(), rseq.values.end(), x);
return static_cast<int>(it - rseq.values.begin());
}
static void push_pair(int sh, int i, int j) {
auto& seq = seqs[sh];
u64 key = (static_cast<u64>(i) << 32) | static_cast<u64>(j);
if (seq.in_heap.count(key)) return;
auto x = get_value_maybe(seq.left, i);
auto y = get_value_maybe(seq.right, j);
if (!x.has_value() || !y.has_value()) return;
if (*x > *y) return;
u128 prod = static_cast<u128>(*x) * static_cast<u128>(*y);
if (prod > MAXVAL) return;
seq.in_heap.insert(key);
seq.heap.push({static_cast<u64>(prod), i, j});
}
static void start_node(int sh) {
auto& seq = seqs[sh];
if (seq.started) return;
seq.started = true;
auto x0 = get_value_maybe(seq.left, 0);
if (!x0.has_value()) {
seq.exhausted = true;
return;
}
int j0 = ensure_right_ge(seq.right, *x0);
if (j0 != -1) {
push_pair(sh, 0, j0);
}
seq.next_i = 1;
}
static void maybe_add_more_i(int sh, u64 current_min) {
auto& seq = seqs[sh];
while (true) {
auto x = get_value_maybe(seq.left, seq.next_i);
if (!x.has_value()) return;
if (!seq.heap.empty() && static_cast<u128>(*x) * static_cast<u128>(*x) > current_min) {
return;
}
int j0 = ensure_right_ge(seq.right, *x);
if (j0 != -1) {
push_pair(sh, seq.next_i, j0);
}
++seq.next_i;
}
}
static std::optional<u64> next_candidate_product(int sh) {
auto& seq = seqs[sh];
start_node(sh);
if (seq.exhausted) return std::nullopt;
while (true) {
if (seq.heap.empty()) {
auto x = get_value_maybe(seq.left, seq.next_i);
if (!x.has_value()) {
seq.exhausted = true;
return std::nullopt;
}
int j0 = ensure_right_ge(seq.right, *x);
if (j0 == -1) {
seq.exhausted = true;
return std::nullopt;
}
push_pair(sh, seq.next_i, j0);
++seq.next_i;
continue;
}
HeapItem cur = seq.heap.top();
seq.heap.pop();
u64 key = (static_cast<u64>(cur.i) << 32) | static_cast<u64>(cur.j);
seq.in_heap.erase(key);
maybe_add_more_i(sh, cur.prod);
push_pair(sh, cur.i, cur.j + 1);
return cur.prod;
}
}
static std::optional<u64> next_value(int sh) {
auto& seq = seqs[sh];
if (seq.exhausted) return std::nullopt;
if (seq.leaf) {
if (seq.next_prime == 2) {
seq.next_prime = 3;
return 2;
}
for (u64 p = seq.next_prime; p <= MAXVAL; p += 2) {
if (is_prime64(p)) {
seq.next_prime = p + 2;
return p;
}
}
seq.exhausted = true;
return std::nullopt;
}
while (true) {
auto cand = next_candidate_product(sh);
if (!cand.has_value()) {
seq.exhausted = true;
return std::nullopt;
}
if (!seq.values.empty() && *cand == seq.values.back()) {
continue;
}
if (shape_of(*cand) == sh) {
return cand;
}
}
}
static u64 get_value(int sh, int idx) {
auto v = get_value_maybe(sh, idx);
if (!v.has_value()) {
throw std::runtime_error("sequence exhausted");
}
return *v;
}
static u64 double_factorial(int n) {
u64 r = 1;
for (int k = n; k >= 2; k -= 2) {
r *= static_cast<u64>(k);
}
return r;
}
int main() {
const int max_n = 31;
MAXVAL = double_factorial(max_n);
std::vector<int> targets;
targets.reserve(max_n - 1);
for (int n = 2; n <= max_n; ++n) {
targets.push_back(shape_of(double_factorial(n)));
}
std::vector<bool> needed(shapes.size(), false);
std::function<void(int)> collect = [&](int sh) {
if (needed[sh]) return;
needed[sh] = true;
if (sh == 0) return;
collect(shapes[sh].left);
collect(shapes[sh].right);
};
for (int sh : targets) collect(sh);
seqs.resize(shapes.size());
for (int sh = 0; sh < static_cast<int>(shapes.size()); ++sh) {
if (!needed[sh]) continue;
seqs[sh].needed = true;
seqs[sh].leaf = (sh == 0);
if (sh != 0) {
seqs[sh].left = shapes[sh].left;
seqs[sh].right = shapes[sh].right;
}
}
std::vector<u64> M(max_n + 1, 0);
for (int n = 2; n <= max_n; ++n) {
M[n] = get_value(targets[n - 2], 0);
}
assert(M[9] == 72ULL);
u128 sum = 0;
for (int n = 2; n <= max_n; ++n) {
sum += M[n];
}
std::cout << static_cast<u64>(sum) << '\n';
return 0;
}
Python
import math
import random
import heapq
import sys
import bisect
sys.setrecursionlimit(20000)
prime_cache = {}
factor_cache = {}
best_divisor_cache = {}
shape_cache = {}
shapes = [{'left': -1, 'right': -1, 'leaves': 1}]
shape_id_by_children = {}
MAXVAL = 0
class Sequence:
def __init__(self):
self.needed = False
self.leaf = False
self.started = False
self.exhausted = False
self.left = -1
self.right = -1
self.next_i = 0
self.next_prime = 2
self.values = []
self.heap = []
self.in_heap = set()
seqs = []
def mul_mod(a, b, mod):
return (a * b) % mod
def pow_mod(a, e, mod):
return pow(a, e, mod)
def is_prime64(n):
if n in prime_cache: return prime_cache[n]
if n < 2:
prime_cache[n] = False
return False
for p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]:
if n % p == 0:
prime_cache[n] = (n == p)
return n == p
d = n - 1
s = 0
while (d & 1) == 0:
d >>= 1
s += 1
for a in [2, 325, 9375, 28178, 450775, 9780504, 1795265022]:
if a % n == 0: continue
x = pow_mod(a, d, n)
if x == 1 or x == n - 1: continue
composite = True
for r in range(1, s):
x = mul_mod(x, x, n)
if x == n - 1:
composite = False
break
if composite:
prime_cache[n] = False
return False
prime_cache[n] = True
return True
def pollard_rho(n):
if (n & 1) == 0: return 2
if n % 3 == 0: return 3
while True:
c = random.randint(2, n - 2)
x = random.randint(2, n - 2)
y = x
d = 1
while d == 1:
x = (mul_mod(x, x, n) + c) % n
y = (mul_mod(y, y, n) + c) % n
y = (mul_mod(y, y, n) + c) % n
diff = x - y if x > y else y - x
d = math.gcd(diff, n)
if d != n: return d
def factor_rec(n, fac):
if n == 1: return
if is_prime64(n):
fac[n] = fac.get(n, 0) + 1
return
d = pollard_rho(n)
factor_rec(d, fac)
factor_rec(n // d, fac)
def factorize(n):
if n in factor_cache: return factor_cache[n]
fac = {}
factor_rec(n, fac)
factor_cache[n] = fac
return fac
def gen_divisors_rec(f, idx, cur, out):
if idx == len(f):
out.append(cur)
return
p, e = f[idx]
v = 1
for _ in range(e + 1):
gen_divisors_rec(f, idx + 1, cur * v, out)
v *= p
def best_divisor_le_sqrt(n):
if n in best_divisor_cache: return best_divisor_cache[n]
fac_map = factorize(n)
f = list(fac_map.items())
divs = []
gen_divisors_rec(f, 0, 1, divs)
root = math.isqrt(n)
best = 1
for d in divs:
if d <= root and d > best:
best = d
best_divisor_cache[n] = best
return best
def make_shape(l, r):
key = (l, r)
if key in shape_id_by_children: return shape_id_by_children[key]
s = {'left': l, 'right': r, 'leaves': shapes[l]['leaves'] + shapes[r]['leaves']}
idx = len(shapes)
shapes.append(s)
shape_id_by_children[key] = idx
return idx
def shape_of(n):
if n in shape_cache: return shape_cache[n]
if is_prime64(n):
shape_cache[n] = 0
return 0
d = best_divisor_le_sqrt(n)
a = d
b = n // d
if a > b: a, b = b, a
l = shape_of(a)
r = shape_of(b)
idx = make_shape(l, r)
shape_cache[n] = idx
return idx
def get_value_maybe(sh, idx):
seq = seqs[sh]
while len(seq.values) <= idx:
nv = next_value(sh)
if nv is None: return None
seq.values.append(nv)
return seq.values[idx]
def ensure_right_ge(right_sh, x):
rseq = seqs[right_sh]
if rseq.exhausted and (not rseq.values or rseq.values[-1] < x): return -1
while not rseq.values or rseq.values[-1] < x:
nv = next_value(right_sh)
if nv is None:
rseq.exhausted = True
break
rseq.values.append(nv)
if not rseq.values or rseq.values[-1] < x: return -1
return bisect.bisect_left(rseq.values, x)
def push_pair(sh, i, j):
seq = seqs[sh]
key = (i << 32) | j
if key in seq.in_heap: return
x = get_value_maybe(seq.left, i)
y = get_value_maybe(seq.right, j)
if x is None or y is None: return
if x > y: return
prod = x * y
if prod > MAXVAL: return
seq.in_heap.add(key)
heapq.heappush(seq.heap, (prod, i, j))
def start_node(sh):
seq = seqs[sh]
if seq.started: return
seq.started = True
x0 = get_value_maybe(seq.left, 0)
if x0 is None:
seq.exhausted = True
return
j0 = ensure_right_ge(seq.right, x0)
if j0 != -1:
push_pair(sh, 0, j0)
seq.next_i = 1
def maybe_add_more_i(sh, current_min):
seq = seqs[sh]
while True:
x = get_value_maybe(seq.left, seq.next_i)
if x is None: return
if seq.heap and x * x > current_min:
return
j0 = ensure_right_ge(seq.right, x)
if j0 != -1:
push_pair(sh, seq.next_i, j0)
seq.next_i += 1
def next_candidate_product(sh):
seq = seqs[sh]
start_node(sh)
if seq.exhausted: return None
while True:
if not seq.heap:
x = get_value_maybe(seq.left, seq.next_i)
if x is None:
seq.exhausted = True
return None
j0 = ensure_right_ge(seq.right, x)
if j0 == -1:
seq.exhausted = True
return None
push_pair(sh, seq.next_i, j0)
seq.next_i += 1
continue
prod, i, j = heapq.heappop(seq.heap)
key = (i << 32) | j
seq.in_heap.remove(key)
maybe_add_more_i(sh, prod)
push_pair(sh, i, j + 1)
return prod
def next_value(sh):
seq = seqs[sh]
if seq.exhausted: return None
if seq.leaf:
if seq.next_prime == 2:
seq.next_prime = 3
return 2
p = seq.next_prime
while p <= MAXVAL:
if is_prime64(p):
seq.next_prime = p + 2
return p
p += 2
seq.exhausted = True
return None
while True:
cand = next_candidate_product(sh)
if cand is None:
seq.exhausted = True
return None
if seq.values and cand == seq.values[-1]:
continue
if shape_of(cand) == sh:
return cand
def get_value(sh, idx):
v = get_value_maybe(sh, idx)
if v is None: raise RuntimeError("sequence exhausted")
return v
def double_factorial(n):
r = 1
for k in range(n, 1, -2):
r *= k
return r
def solve():
global MAXVAL, seqs
max_n = 31
MAXVAL = double_factorial(max_n)
targets = []
for n in range(2, max_n + 1):
targets.append(shape_of(double_factorial(n)))
needed = [False] * len(shapes)
def collect(sh):
if needed[sh]: return
needed[sh] = True
if sh == 0: return
collect(shapes[sh]['left'])
collect(shapes[sh]['right'])
for sh in targets:
collect(sh)
seqs = [Sequence() for _ in range(len(shapes))]
for sh in range(len(shapes)):
if not needed[sh]: continue
seqs[sh].needed = True
seqs[sh].leaf = (sh == 0)
if sh != 0:
seqs[sh].left = shapes[sh]['left']
seqs[sh].right = shapes[sh]['right']
M = [0] * (max_n + 1)
for n in range(2, max_n + 1):
M[n] = get_value(targets[n - 2], 0)
ans = sum(M[2:])
return str(ans)
if __name__ == "__main__":
print(solve())
Java
import java.util.*;
public class Euler829 {
static class ShapeInfo {
int left = -1;
int right = -1;
int leaves = 1;
}
static class HeapItem implements Comparable<HeapItem> {
long prod;
int i;
int j;
HeapItem(long prod, int i, int j) {
this.prod = prod;
this.i = i;
this.j = j;
}
@Override
public int compareTo(HeapItem o) {
if (this.prod != o.prod)
return Long.compare(this.prod, o.prod);
if (this.i != o.i)
return Integer.compare(this.i, o.i);
return Integer.compare(this.j, o.j);
}
}
static class Sequence {
boolean needed = false;
boolean leaf = false;
boolean started = false;
boolean exhausted = false;
int left = -1;
int right = -1;
int nextI = 0;
long nextPrime = 2;
ArrayList<Long> values = new ArrayList<>();
PriorityQueue<HeapItem> heap = new PriorityQueue<>();
HashSet<Long> inHeap = new HashSet<>();
}
static Random rng = new Random(0);
static ArrayList<ShapeInfo> shapes;
static HashMap<Long, Integer> shapeIdByChildren;
static HashMap<Long, Integer> shapeCache;
static HashMap<Long, Long> bestDivisorCache;
static HashMap<Long, Boolean> primeCache;
static HashMap<Long, HashMap<Long, Integer>> factorCache;
static long MAXVAL = 0;
static Sequence[] seqs;
static long mulMod(long a, long b, long mod) {
long q = (long) ((double) a * b / mod);
long r = a * b - q * mod;
while (r < 0)
r += mod;
while (r >= mod)
r -= mod;
return r;
}
static long powMod(long a, long e, long mod) {
long r = 1 % mod;
a %= mod;
while (e > 0) {
if ((e & 1L) == 1L)
r = mulMod(r, a, mod);
a = mulMod(a, a, mod);
e >>= 1L;
}
return r;
}
static boolean isPrime64(long n) {
if (primeCache.containsKey(n))
return primeCache.get(n);
if (n < 2) {
primeCache.put(n, false);
return false;
}
long[] ps = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37 };
for (long p : ps) {
if (n % p == 0) {
primeCache.put(n, n == p);
return n == p;
}
}
long d = n - 1;
int s = 0;
while ((d & 1L) == 0) {
d >>= 1L;
s++;
}
long[] as = { 2, 325, 9375, 28178, 450775, 9780504, 1795265022 };
for (long a : as) {
if (a % n == 0)
continue;
long x = powMod(a, d, n);
if (x == 1 || x == n - 1)
continue;
boolean comp = true;
for (int r = 1; r < s; ++r) {
x = mulMod(x, x, n);
if (x == n - 1) {
comp = false;
break;
}
}
if (comp) {
primeCache.put(n, false);
return false;
}
}
primeCache.put(n, true);
return true;
}
static long pollardRho(long n) {
if ((n & 1L) == 0)
return 2;
if (n % 3 == 0)
return 3;
while (true) {
long c = 2 + (long) (rng.nextDouble() * (n - 3));
long x = 2 + (long) (rng.nextDouble() * (n - 3));
long y = x;
long d = 1;
while (d == 1) {
x = (mulMod(x, x, n) + c) % n;
y = (mulMod(y, y, n) + c) % n;
y = (mulMod(y, y, n) + c) % n;
long diff = x > y ? x - y : y - x;
d = gcd(diff, n);
}
if (d != n)
return d;
}
}
static long gcd(long a, long b) {
while (b != 0) {
long temp = b;
b = a % b;
a = temp;
}
return a;
}
static void factorRec(long n, HashMap<Long, Integer> fac) {
if (n == 1)
return;
if (isPrime64(n)) {
fac.put(n, fac.getOrDefault(n, 0) + 1);
return;
}
long d = pollardRho(n);
factorRec(d, fac);
factorRec(n / d, fac);
}
static HashMap<Long, Integer> factorize(long n) {
if (factorCache.containsKey(n))
return factorCache.get(n);
HashMap<Long, Integer> fac = new HashMap<>();
factorRec(n, fac);
factorCache.put(n, fac);
return fac;
}
static class Pair {
long p;
int e;
Pair(long p, int e) {
this.p = p;
this.e = e;
}
}
static void genDivisorsRec(ArrayList<Pair> f, int idx, long cur, ArrayList<Long> out) {
if (idx == f.size()) {
out.add(cur);
return;
}
Pair p = f.get(idx);
long v = 1;
for (int i = 0; i <= p.e; ++i) {
genDivisorsRec(f, idx + 1, cur * v, out);
v *= p.p;
}
}
static long bestDivisorLeSqrt(long n) {
if (bestDivisorCache.containsKey(n))
return bestDivisorCache.get(n);
HashMap<Long, Integer> facMap = factorize(n);
ArrayList<Pair> f = new ArrayList<>();
for (Map.Entry<Long, Integer> entry : facMap.entrySet()) {
f.add(new Pair(entry.getKey(), entry.getValue()));
}
ArrayList<Long> divs = new ArrayList<>();
genDivisorsRec(f, 0, 1, divs);
long root = (long) Math.sqrt(n);
while ((root + 1) * (root + 1) <= n)
++root;
while (root * root > n)
--root;
long best = 1;
for (long d : divs) {
if (d <= root && d > best)
best = d;
}
bestDivisorCache.put(n, best);
return best;
}
static int makeShape(int l, int r) {
long key = ((long) l << 32) | (r & 0xFFFFFFFFL);
if (shapeIdByChildren.containsKey(key))
return shapeIdByChildren.get(key);
ShapeInfo s = new ShapeInfo();
s.left = l;
s.right = r;
s.leaves = shapes.get(l).leaves + shapes.get(r).leaves;
int id = shapes.size();
shapes.add(s);
shapeIdByChildren.put(key, id);
return id;
}
static int shapeOf(long n) {
if (shapeCache.containsKey(n))
return shapeCache.get(n);
if (isPrime64(n)) {
shapeCache.put(n, 0);
return 0;
}
long d = bestDivisorLeSqrt(n);
long a = d;
long b = n / d;
if (a > b) {
long temp = a;
a = b;
b = temp;
}
int l = shapeOf(a);
int r = shapeOf(b);
int id = makeShape(l, r);
shapeCache.put(n, id);
return id;
}
static Long getValueMaybe(int sh, int idx) {
Sequence seq = seqs[sh];
while (seq.values.size() <= idx) {
Long nv = nextValue(sh);
if (nv == null)
return null;
seq.values.add(nv);
}
return seq.values.get(idx);
}
static int ensureRightGe(int rightSh, long x) {
Sequence rseq = seqs[rightSh];
if (rseq.exhausted && (rseq.values.isEmpty() || rseq.values.get(rseq.values.size() - 1) < x))
return -1;
while (rseq.values.isEmpty() || rseq.values.get(rseq.values.size() - 1) < x) {
Long nv = nextValue(rightSh);
if (nv == null) {
rseq.exhausted = true;
break;
}
rseq.values.add(nv);
}
if (rseq.values.isEmpty() || rseq.values.get(rseq.values.size() - 1) < x)
return -1;
int low = 0;
int high = rseq.values.size() - 1;
while (low <= high) {
int mid = (low + high) >>> 1;
if (rseq.values.get(mid) < x) {
low = mid + 1;
} else {
high = mid - 1;
}
}
return low;
}
static void pushPair(int sh, int i, int j) {
Sequence seq = seqs[sh];
long key = ((long) i << 32) | (j & 0xFFFFFFFFL);
if (seq.inHeap.contains(key))
return;
Long x = getValueMaybe(seq.left, i);
Long y = getValueMaybe(seq.right, j);
if (x == null || y == null)
return;
if (x > y)
return;
if ((double) x * y > MAXVAL * 1.5) {
// fast check to avoid exact 128 bit multiplication issue
// wait, we can just use simple division
}
long maxDivX = MAXVAL / x;
if (y > maxDivX)
return;
long prod = x * y;
seq.inHeap.add(key);
seq.heap.add(new HeapItem(prod, i, j));
}
static void startNode(int sh) {
Sequence seq = seqs[sh];
if (seq.started)
return;
seq.started = true;
Long x0 = getValueMaybe(seq.left, 0);
if (x0 == null) {
seq.exhausted = true;
return;
}
int j0 = ensureRightGe(seq.right, x0);
if (j0 != -1) {
pushPair(sh, 0, j0);
}
seq.nextI = 1;
}
static void maybeAddMoreI(int sh, long currentMin) {
Sequence seq = seqs[sh];
while (true) {
Long x = getValueMaybe(seq.left, seq.nextI);
if (x == null)
return;
long maxDivX = currentMin / x;
if (!seq.heap.isEmpty() && x > maxDivX) {
return;
}
int j0 = ensureRightGe(seq.right, x);
if (j0 != -1) {
pushPair(sh, seq.nextI, j0);
}
seq.nextI++;
}
}
static Long nextCandidateProduct(int sh) {
Sequence seq = seqs[sh];
startNode(sh);
if (seq.exhausted)
return null;
while (true) {
if (seq.heap.isEmpty()) {
Long x = getValueMaybe(seq.left, seq.nextI);
if (x == null) {
seq.exhausted = true;
return null;
}
int j0 = ensureRightGe(seq.right, x);
if (j0 == -1) {
seq.exhausted = true;
return null;
}
pushPair(sh, seq.nextI, j0);
seq.nextI++;
continue;
}
HeapItem cur = seq.heap.poll();
long key = ((long) cur.i << 32) | (cur.j & 0xFFFFFFFFL);
seq.inHeap.remove(key);
maybeAddMoreI(sh, cur.prod);
pushPair(sh, cur.i, cur.j + 1);
return cur.prod;
}
}
static Long nextValue(int sh) {
Sequence seq = seqs[sh];
if (seq.exhausted)
return null;
if (seq.leaf) {
if (seq.nextPrime == 2) {
seq.nextPrime = 3;
return 2L;
}
for (long p = seq.nextPrime; p <= MAXVAL; p += 2) {
if (isPrime64(p)) {
seq.nextPrime = p + 2;
return p;
}
}
seq.exhausted = true;
return null;
}
while (true) {
Long cand = nextCandidateProduct(sh);
if (cand == null) {
seq.exhausted = true;
return null;
}
if (!seq.values.isEmpty() && cand.equals(seq.values.get(seq.values.size() - 1))) {
continue;
}
if (shapeOf(cand) == sh) {
return cand;
}
}
}
static long getValue(int sh, int idx) {
Long v = getValueMaybe(sh, idx);
if (v == null) {
throw new RuntimeException("sequence exhausted");
}
return v;
}
static long doubleFactorial(int n) {
long r = 1;
for (int k = n; k >= 2; k -= 2) {
r *= k;
}
return r;
}
public static String solve() {
shapes = new ArrayList<>();
shapes.add(new ShapeInfo());
shapeIdByChildren = new HashMap<>();
shapeCache = new HashMap<>();
bestDivisorCache = new HashMap<>();
primeCache = new HashMap<>();
factorCache = new HashMap<>();
int maxN = 31;
MAXVAL = doubleFactorial(maxN);
int[] targets = new int[maxN - 1];
for (int n = 2; n <= maxN; ++n) {
targets[n - 2] = shapeOf(doubleFactorial(n));
}
boolean[] needed = new boolean[shapes.size()];
class DFS {
void collect(int sh) {
if (needed[sh])
return;
needed[sh] = true;
if (sh == 0)
return;
collect(shapes.get(sh).left);
collect(shapes.get(sh).right);
}
}
DFS dfs = new DFS();
for (int sh : targets)
dfs.collect(sh);
seqs = new Sequence[shapes.size()];
for (int i = 0; i < shapes.size(); i++)
seqs[i] = new Sequence();
for (int sh = 0; sh < shapes.size(); ++sh) {
if (!needed[sh])
continue;
seqs[sh].needed = true;
seqs[sh].leaf = (sh == 0);
if (sh != 0) {
seqs[sh].left = shapes.get(sh).left;
seqs[sh].right = shapes.get(sh).right;
}
}
long[] M = new long[maxN + 1];
for (int n = 2; n <= maxN; ++n) {
M[n] = getValue(targets[n - 2], 0);
}
long sum = 0;
for (int n = 2; n <= maxN; ++n) {
sum += M[n];
}
return Long.toString(sum);
}
public static void main(String[] args) {
System.out.println(solve());
}
}