Problem 167: Investigating Ulam Sequences
View on Project EulerProject Euler Problem 167 Solution
EulerSolve provides an optimized solution for Project Euler Problem 167, Investigating Ulam Sequences, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary An Ulam sequence \(U(a,b)\) starts with \(a\) and \(b\). Every later term is the smallest integer larger than the previous term that can be written in exactly one way as a sum of two distinct earlier terms. In this problem we must evaluate $$S(k)=\sum_{n=2}^{10} u_k\bigl(U(2,2n+1)\bigr),$$ so the nine relevant sequences are \(U(2,5),U(2,7),\dots,U(2,21)\). The target index is enormous in the original problem, so a term-by-term generator is useless. The solution works only because this particular family of Ulam sequences has a highly constrained structure. Mathematical Approach Fix one sequence \(U(2,v)\) with \(v=2n+1\) odd. The core idea is that the even terms become trivial, and once that happens the odd terms are governed by a linear recurrence over \(\mathbb{F}_2\). The special even term \(E=2(v+1)\) For the nine values \(v=5,7,\dots,21\), the implementation exploits a specific fact about \(U(2,v)\): after the initial term \(2\), there is exactly one further even term, namely $$E=2(v+1).$$ Beyond that point the new members are odd. This matters because every odd candidate \(x>E\) can only be formed as $$x=2+(x-2),\qquad x=E+(x-E).$$ Odd numbers cannot come from odd plus odd, and there are no other even Ulam terms available. Therefore an odd number \(x>E\) enters the sequence exactly when one of \(x-2\) and \(x-E\) is already present and the other is not....
Detailed mathematical approach
Problem Summary
An Ulam sequence \(U(a,b)\) starts with \(a\) and \(b\). Every later term is the smallest integer larger than the previous term that can be written in exactly one way as a sum of two distinct earlier terms. In this problem we must evaluate
$$S(k)=\sum_{n=2}^{10} u_k\bigl(U(2,2n+1)\bigr),$$
so the nine relevant sequences are \(U(2,5),U(2,7),\dots,U(2,21)\). The target index is enormous in the original problem, so a term-by-term generator is useless. The solution works only because this particular family of Ulam sequences has a highly constrained structure.
Mathematical Approach
Fix one sequence \(U(2,v)\) with \(v=2n+1\) odd. The core idea is that the even terms become trivial, and once that happens the odd terms are governed by a linear recurrence over \(\mathbb{F}_2\).
The special even term \(E=2(v+1)\)
For the nine values \(v=5,7,\dots,21\), the implementation exploits a specific fact about \(U(2,v)\): after the initial term \(2\), there is exactly one further even term, namely
$$E=2(v+1).$$
Beyond that point the new members are odd. This matters because every odd candidate \(x>E\) can only be formed as
$$x=2+(x-2),\qquad x=E+(x-E).$$
Odd numbers cannot come from odd plus odd, and there are no other even Ulam terms available. Therefore an odd number \(x>E\) enters the sequence exactly when one of \(x-2\) and \(x-E\) is already present and the other is not.
Odd membership becomes an XOR recurrence
Define an indicator for odd membership by
$$b_t=\chi(2t+1),$$
where \(b_t=1\) if \(2t+1\) belongs to \(U(2,v)\), and \(b_t=0\) otherwise. Also write
$$g=v+1,\qquad E=2g.$$
For odd \(x=2t+1>E\), the previous observation becomes
$$\chi(x)=\chi(x-2)\oplus\chi(x-E),$$
which is the same as
$$b_t=b_{t-1}\oplus b_{t-g}\qquad (t\ge g).$$
So one window of \(g\) consecutive bits determines every later odd term. Instead of thinking about huge integers directly, we can think about a \(g\)-bit state evolving by a deterministic update rule.
Worked example: \(U(2,5)\)
For \(v=5\) we have \(g=6\) and \(E=12\). The beginning of the sequence is
$$2,5,7,9,11,12,13,15,19,23,\dots$$
Look at the odd numbers in the initial window \(\{1,3,5,7,9,11\}\). Their membership bits are
$$b_0,\dots,b_5=(0,0,1,1,1,1).$$
Now apply the recurrence \(b_t=b_{t-1}\oplus b_{t-6}\):
$$b_6=b_5\oplus b_0=1\oplus 0=1,$$
so \(13=2\cdot 6+1\) is in the sequence. Next,
$$b_7=b_6\oplus b_1=1\oplus 0=1,$$
so \(15\) is in. Then
$$b_8=b_7\oplus b_2=1\oplus 1=0,$$
so \(17\) is skipped, and
$$b_9=b_8\oplus b_3=0\oplus 1=1,$$
so \(19\) appears. This is exactly the pattern produced by the brute-force prefix and then continued by the fast method.
Why the bit sequence is eventually periodic
At step \(t\), the recurrence only needs the last \(g\) bits, for example the state
$$\bigl(b_{t-g+1},b_{t-g+2},\dots,b_t\bigr).$$
The next bit is determined by the oldest and newest entries, so the transition is deterministic. There are only \(2^g\) possible states, hence some state must repeat. Once the same \(g\)-bit state appears again, every later update repeats as well, so the odd-membership sequence has a finite preperiod followed by a pure cycle.
Converting the huge index into an odd rank
The full Ulam sequence is almost all odd, but the exceptional even term \(E\) must still be inserted at the correct position. Its index is
$$r_E=2+\#\{x<E:x\equiv 1 \pmod 2,\ x\in U(2,v)\}.$$
The first position is occupied by \(2\), then come all odd members below \(E\), and then \(E\) itself. Therefore
$$u_1(v)=2,\qquad u_{r_E}(v)=E,$$
and every other query is converted to an odd rank
$$r_{\mathrm{odd}}= \begin{cases} k-1,&k<r_E,\\ k-2,&k>r_E. \end{cases}$$
Prefix sums of the bit sequence tell how many odd Ulam numbers have appeared up to a given odd value. Once the preperiod length and the number of 1-bits per cycle are known, the algorithm can skip whole cycles arithmetically and then locate the exact odd value with one binary search.
How the Code Works
Bootstrapping the recurrence
The C++, Python, and Java implementations first generate only a short brute-force prefix, just far enough to know which odd numbers up to \(2g-1\) are members. That gives one complete seed window \(b_0,\dots,b_{g-1}\) and also determines how many odd members lie below the special even term \(E\).
From there the implementation builds a \(g\)-bit state, computes prefix counts of 1-bits, and stops using ordinary Ulam generation. All later odd membership decisions come from the XOR recurrence.
Detecting the cycle and answering \(u_k\)
Each new bit updates the sliding \(g\)-bit state. The first repeated state gives the beginning and end of the eventual cycle. The implementation stores the number of odd members before the cycle, the number contributed by one cycle, and prefix counts inside the cycle itself.
To answer the query, the implementation handles the special cases \(2\) and \(E\), converts the requested position into an odd rank, skips as many full cycles as possible, and then binary-searches the relevant prefix table to recover the matching odd number \(2t+1\). The C++ and Java implementations do this independently for the nine values of \(v\) in parallel; the Python implementation applies the same logic serially.
Complexity Analysis
For one sequence \(U(2,v)\), the main cost is cycle detection on the \(g=v+1\) bit states. In the worst case this takes \(O(2^g)\) time and \(O(2^g)\) memory for the table that records the first visit to each state. Here \(v\le 21\), so \(g\le 22\) and the state space is at most \(2^{22}=4{,}194{,}304\), which is entirely manageable.
After that preprocessing, finding \(u_k(v)\) is \(O(\log(P+\lambda))\), where \(P\) is the preperiod length and \(\lambda\) is the cycle length, because the remaining work is arithmetic plus one binary search on prefix sums. The full problem repeats this for only nine independent sequences and adds the results.
Footnotes and References
- Project Euler problem page: Problem 167 - Investigating Ulam Sequences
- Ulam numbers and Ulam sequences: Wikipedia - Ulam number
- Recurrence relations: Wikipedia - Recurrence relation
- Linear recurrences over \(\mathbb{F}_2\): Wikipedia - Linear feedback shift register
Problem 167 source code
C++
#include <algorithm>
#include <atomic>
#include <cstdint>
#include <exception>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <thread>
#include <unordered_map>
#include <unordered_set>
#include <vector>
namespace {
using u8 = std::uint8_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;
struct Options {
u64 target_k = 100000000000ULL;
bool run_checkpoints = true;
unsigned requested_threads = 0U;
};
bool parse_unsigned_after_prefix(const std::string& arg,
const std::string& prefix,
unsigned& value) {
if (arg.rfind(prefix, 0U) != 0U) {
return false;
}
const std::string tail = arg.substr(prefix.size());
if (tail.empty()) {
return false;
}
u64 parsed = 0ULL;
for (const char c : tail) {
if (c < '0' || c > '9') {
return false;
}
parsed = parsed * 10ULL + static_cast<u64>(c - '0');
if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
return false;
}
}
value = static_cast<unsigned>(parsed);
return true;
}
bool parse_u64_after_prefix(const std::string& arg,
const std::string& prefix,
u64& value) {
if (arg.rfind(prefix, 0U) != 0U) {
return false;
}
const std::string tail = arg.substr(prefix.size());
if (tail.empty()) {
return false;
}
u64 parsed = 0ULL;
for (const char c : tail) {
if (c < '0' || c > '9') {
return false;
}
const u64 digit = static_cast<u64>(c - '0');
if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
throw std::overflow_error("--k overflow");
}
parsed = parsed * 10ULL + digit;
}
value = parsed;
return true;
}
bool parse_arguments(const int argc, char** argv, Options& options) {
for (int i = 1; i < argc; ++i) {
const std::string arg(argv[i]);
if (arg == "--skip-checkpoints") {
options.run_checkpoints = false;
continue;
}
if (parse_unsigned_after_prefix(arg, "--threads=", options.requested_threads)) {
continue;
}
if (parse_u64_after_prefix(arg, "--k=", options.target_k)) {
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
return true;
}
unsigned pick_thread_count(const unsigned requested) {
if (requested > 0U) {
return requested;
}
unsigned hw = std::thread::hardware_concurrency();
if (hw == 0U) {
hw = 4U;
}
return hw;
}
std::vector<u64> brute_ulam_terms(const u64 a, const u64 b, const std::size_t count) {
if (count == 0U) {
return {};
}
if (count == 1U) {
return {a};
}
std::vector<u64> seq;
seq.reserve(count);
seq.push_back(a);
seq.push_back(b);
std::unordered_map<u64, u8> representation_count;
representation_count.reserve(count * count / 2 + 64U);
representation_count[a + b] = 1U;
u64 candidate = b + 1U;
while (seq.size() < count) {
while (true) {
const auto it = representation_count.find(candidate);
const u8 ways = (it == representation_count.end()) ? 0U : it->second;
if (ways == 1U) {
break;
}
++candidate;
}
const u64 next = candidate;
for (const u64 x : seq) {
if (x == next) {
continue;
}
const u64 sum = x + next;
u8& ways = representation_count[sum];
if (ways < 2U) {
++ways;
}
}
seq.push_back(next);
candidate = next + 1U;
}
return seq;
}
std::vector<u64> brute_ulam_until_value(const u64 a, const u64 b, const u64 min_last_value) {
std::vector<u64> seq;
seq.reserve(256U);
seq.push_back(a);
seq.push_back(b);
std::unordered_map<u64, u8> representation_count;
representation_count.reserve(32768U);
representation_count[a + b] = 1U;
u64 candidate = b + 1U;
while (seq.back() < min_last_value) {
while (true) {
const auto it = representation_count.find(candidate);
const u8 ways = (it == representation_count.end()) ? 0U : it->second;
if (ways == 1U) {
break;
}
++candidate;
}
const u64 next = candidate;
for (const u64 x : seq) {
if (x == next) {
continue;
}
const u64 sum = x + next;
u8& ways = representation_count[sum];
if (ways < 2U) {
++ways;
}
}
seq.push_back(next);
candidate = next + 1U;
}
return seq;
}
struct UlamOddCycle {
int v = 0;
int gap = 0;
u64 second_even = 0ULL;
std::vector<u8> bits;
std::vector<u64> prefix_ones;
u64 cycle_t_start = 0ULL;
u64 cycle_t_end = 0ULL;
u64 cycle_len = 0ULL;
u64 ones_before_cycle = 0ULL;
u64 ones_per_cycle = 0ULL;
std::vector<u64> cycle_prefix_ones;
u64 odd_less_than_second_even = 0ULL;
u64 second_even_index = 0ULL;
explicit UlamOddCycle(const int vv) : v(vv), gap(vv + 1), second_even(2ULL * static_cast<u64>(vv + 1)) {
if (v < 1) {
throw std::invalid_argument("v must be positive");
}
if (gap <= 0 || gap > 30) {
throw std::runtime_error("Unexpected recurrence gap");
}
const u64 seed_limit = 2ULL * static_cast<u64>(gap) + 1ULL;
const std::vector<u64> seed = brute_ulam_until_value(2ULL, static_cast<u64>(v), seed_limit);
std::unordered_set<u64> members;
members.reserve(seed.size() * 2U + 16U);
for (const u64 x : seed) {
members.insert(x);
}
bits.resize(static_cast<std::size_t>(gap), 0U);
for (int t = 0; t < gap; ++t) {
const u64 odd_value = 2ULL * static_cast<u64>(t) + 1ULL;
bits[static_cast<std::size_t>(t)] = (members.find(odd_value) != members.end()) ? 1U : 0U;
}
prefix_ones.reserve(1U << 14U);
prefix_ones.push_back(0ULL);
for (const u8 bit : bits) {
prefix_ones.push_back(prefix_ones.back() + static_cast<u64>(bit));
}
odd_less_than_second_even = prefix_ones[static_cast<std::size_t>(gap)];
second_even_index = 2ULL + odd_less_than_second_even;
const u32 state_count = 1U << static_cast<u32>(gap);
std::vector<int> seen(state_count, -1);
u32 state = 0U;
for (int i = 0; i < gap; ++i) {
state = (state << 1U) | static_cast<u32>(bits[static_cast<std::size_t>(i)]);
}
seen[state] = 0;
const u32 lower_mask = (gap == 1) ? 0U : ((1U << static_cast<u32>(gap - 1)) - 1U);
int step = 0;
int cycle_step_start = -1;
int cycle_step_end = -1;
while (true) {
const u32 oldest = (state >> static_cast<u32>(gap - 1)) & 1U;
const u32 newest = state & 1U;
const u32 next_bit = oldest ^ newest;
bits.push_back(static_cast<u8>(next_bit));
prefix_ones.push_back(prefix_ones.back() + static_cast<u64>(next_bit));
state = ((state & lower_mask) << 1U) | next_bit;
++step;
const int prev = seen[state];
if (prev >= 0) {
cycle_step_start = prev;
cycle_step_end = step;
break;
}
seen[state] = step;
}
cycle_t_start = static_cast<u64>(gap - 1 + cycle_step_start);
cycle_t_end = static_cast<u64>(gap - 1 + cycle_step_end);
cycle_len = cycle_t_end - cycle_t_start;
ones_before_cycle = prefix_ones[static_cast<std::size_t>(cycle_t_start + 1ULL)];
ones_per_cycle = prefix_ones[static_cast<std::size_t>(cycle_t_end + 1ULL)] - ones_before_cycle;
if (cycle_len == 0ULL || ones_per_cycle == 0ULL) {
throw std::runtime_error("Degenerate cycle detected");
}
cycle_prefix_ones.assign(static_cast<std::size_t>(cycle_len + 1ULL), 0ULL);
for (u64 i = 0ULL; i < cycle_len; ++i) {
const u64 t = cycle_t_start + 1ULL + i;
cycle_prefix_ones[static_cast<std::size_t>(i + 1ULL)] =
cycle_prefix_ones[static_cast<std::size_t>(i)] + static_cast<u64>(bits[static_cast<std::size_t>(t)]);
}
}
u64 odd_term_by_rank(const u64 odd_rank) const {
if (odd_rank == 0ULL) {
throw std::invalid_argument("odd rank must be positive");
}
u64 t = 0ULL;
if (odd_rank <= ones_before_cycle) {
const auto begin_it = prefix_ones.begin() + 1;
const auto end_it = prefix_ones.begin() + static_cast<std::ptrdiff_t>(cycle_t_start + 2ULL);
const auto it = std::lower_bound(begin_it, end_it, odd_rank);
t = static_cast<u64>(std::distance(prefix_ones.begin(), it)) - 1ULL;
} else {
const u64 rem = odd_rank - ones_before_cycle;
const u64 full_cycles = (rem - 1ULL) / ones_per_cycle;
const u64 rem_inside_cycle = rem - full_cycles * ones_per_cycle;
const auto begin_it = cycle_prefix_ones.begin() + 1;
const auto end_it = cycle_prefix_ones.end();
const auto it = std::lower_bound(begin_it, end_it, rem_inside_cycle);
const u64 d = static_cast<u64>(std::distance(cycle_prefix_ones.begin(), it));
t = cycle_t_start + full_cycles * cycle_len + d;
}
return 2ULL * t + 1ULL;
}
u64 kth_term(const u64 k) const {
if (k == 0ULL) {
throw std::invalid_argument("k must be positive");
}
if (k == 1ULL) {
return 2ULL;
}
if (k == second_even_index) {
return second_even;
}
if (k < second_even_index) {
return odd_term_by_rank(k - 1ULL);
}
return odd_term_by_rank(k - 2ULL);
}
};
bool run_checkpoints() {
{
const std::vector<u64> sample = brute_ulam_terms(1ULL, 2ULL, 7U);
const std::vector<u64> expected{1ULL, 2ULL, 3ULL, 4ULL, 6ULL, 8ULL, 11ULL};
if (sample != expected) {
std::cerr << "Checkpoint failed: U(1,2) sample mismatch\n";
return false;
}
}
const std::vector<u64> ks{1ULL, 2ULL, 3ULL, 4ULL, 5ULL, 6ULL, 10ULL, 20ULL,
50ULL, 100ULL, 250ULL, 500ULL, 1000ULL, 1500ULL, 2000ULL};
for (int n = 2; n <= 10; ++n) {
const int v = 2 * n + 1;
const UlamOddCycle fast(v);
const std::vector<u64> brute = brute_ulam_terms(2ULL, static_cast<u64>(v), 2000U);
for (const u64 k : ks) {
const u64 a = fast.kth_term(k);
const u64 b = brute[static_cast<std::size_t>(k - 1ULL)];
if (a != b) {
std::cerr << "Checkpoint failed: mismatch at v=" << v << ", k=" << k
<< ", fast=" << a << ", brute=" << b << '\n';
return false;
}
}
std::unordered_set<u64> members;
members.reserve(brute.size() * 2U + 16U);
for (const u64 x : brute) {
members.insert(x);
}
const u64 second_even = 2ULL * static_cast<u64>(v + 1);
for (const u64 x : brute) {
if ((x & 1ULL) == 0ULL && x != 2ULL && x != second_even) {
std::cerr << "Checkpoint failed: unexpected extra even term for v=" << v
<< ", value=" << x << '\n';
return false;
}
}
for (u64 odd = second_even + 1ULL; odd <= brute.back(); odd += 2ULL) {
const bool lhs = (members.find(odd) != members.end());
const bool rhs = (members.find(odd - 2ULL) != members.end()) ^
(members.find(odd - second_even) != members.end());
if (lhs != rhs) {
std::cerr << "Checkpoint failed: odd recurrence mismatch at v=" << v
<< ", odd=" << odd << '\n';
return false;
}
}
}
return true;
}
u64 solve_target_sum(const u64 target_k, const unsigned requested_threads) {
constexpr int first_n = 2;
constexpr int last_n = 10;
constexpr int task_count = last_n - first_n + 1;
std::vector<u64> values(task_count, 0ULL);
const unsigned thread_count =
std::max(1U, std::min<unsigned>(pick_thread_count(requested_threads), task_count));
if (thread_count == 1U) {
for (int i = 0; i < task_count; ++i) {
const int n = first_n + i;
const int v = 2 * n + 1;
const UlamOddCycle solver(v);
values[static_cast<std::size_t>(i)] = solver.kth_term(target_k);
}
} else {
std::atomic<int> next_task(0);
std::vector<std::thread> workers;
workers.reserve(thread_count);
for (unsigned t = 0; t < thread_count; ++t) {
workers.emplace_back([&]() {
while (true) {
const int task = next_task.fetch_add(1, std::memory_order_relaxed);
if (task >= task_count) {
return;
}
const int n = first_n + task;
const int v = 2 * n + 1;
const UlamOddCycle solver(v);
values[static_cast<std::size_t>(task)] = solver.kth_term(target_k);
}
});
}
for (std::thread& worker : workers) {
worker.join();
}
}
u64 sum = 0ULL;
for (const u64 value : values) {
sum += value;
}
return sum;
}
} // namespace
int main(int argc, char** argv) {
try {
Options options;
if (!parse_arguments(argc, argv, options)) {
return 1;
}
if (options.run_checkpoints && !run_checkpoints()) {
return 1;
}
const u64 answer = solve_target_sum(options.target_k, options.requested_threads);
std::cout << answer << '\n';
} catch (const std::exception& ex) {
std::cerr << "Error: " << ex.what() << '\n';
return 1;
}
return 0;
}
Python
def solve():
target_k = 100000000000
def brute_ulam(a, b, min_last):
seq = [a, b]
reps = {a + b: 1}
cand = b + 1
while seq[-1] < min_last:
while True:
w = reps.get(cand, 0)
if w == 1: break
cand += 1
nxt = cand
for x in seq:
if x == nxt: continue
s = x + nxt
reps[s] = min(reps.get(s, 0) + 1, 2)
seq.append(nxt)
cand = nxt + 1
return seq
def kth_term(v, k):
gap = v + 1
second_even = 2 * (v + 1)
seed_limit = 2 * gap + 1
seed = brute_ulam(2, v, seed_limit)
members = set(seed)
bits = [1 if (2*t+1) in members else 0 for t in range(gap)]
prefix = [0] * (gap + 1)
for i in range(gap): prefix[i+1] = prefix[i] + bits[i]
odd_lt_se = prefix[gap]
se_idx = 2 + odd_lt_se
state_count = 1 << gap
seen = [-1] * state_count
state = 0
for i in range(gap): state = (state << 1) | bits[i]
seen[state] = 0
lower_mask = (1 << (gap-1)) - 1 if gap > 1 else 0
bit_list = list(bits)
po_list = list(prefix)
step = 0
while True:
oldest = (state >> (gap-1)) & 1
newest = state & 1
nb = oldest ^ newest
bit_list.append(nb)
po_list.append(po_list[-1] + nb)
state = ((state & lower_mask) << 1) | nb
step += 1
prev = seen[state]
if prev >= 0:
cycle_start = prev; cycle_end = step; break
seen[state] = step
cT_start = gap - 1 + cycle_start
cT_end = gap - 1 + cycle_end
cycle_len = cT_end - cT_start
ones_before = po_list[cT_start + 1]
ones_per_cycle = po_list[cT_end + 1] - ones_before
cpf = [0] * (cycle_len + 1)
for i in range(cycle_len):
cpf[i+1] = cpf[i] + bit_list[cT_start + 1 + i]
if k == 1: return 2
if k == se_idx: return second_even
odd_rank = k - 1 if k < se_idx else k - 2
if odd_rank <= ones_before:
lo, hi = 0, cT_start + 1
while lo < hi:
mid = (lo + hi) // 2
if po_list[mid + 1] >= odd_rank: hi = mid
else: lo = mid + 1
return 2 * lo + 1
else:
rem = odd_rank - ones_before
full_cycles = (rem - 1) // ones_per_cycle
rem_inside = rem - full_cycles * ones_per_cycle
lo, hi = 0, cycle_len
while lo < hi:
mid = (lo + hi) // 2
if cpf[mid + 1] >= rem_inside: hi = mid
else: lo = mid + 1
t = cT_start + full_cycles * cycle_len + lo + 1
return 2 * t + 1
total = 0
for n in range(2, 11):
v = 2 * n + 1
total += kth_term(v, target_k)
return str(total)
if __name__ == '__main__':
print(solve())
Java
import java.util.*;
public class Euler167 {
static List<Long> bruteUlamUntilValue(long a, long b, long minLast) {
List<Long> seq = new ArrayList<>();
seq.add(a);
seq.add(b);
Map<Long, Integer> reps = new HashMap<>();
reps.put(a + b, 1);
long cand = b + 1;
while (seq.get(seq.size() - 1) < minLast) {
while (true) {
int w = reps.getOrDefault(cand, 0);
if (w == 1)
break;
cand++;
}
long next = cand;
for (long x : seq) {
if (x == next)
continue;
reps.merge(x + next, 1, (o, n) -> Math.min(o + n, 2));
}
seq.add(next);
cand = next + 1;
}
return seq;
}
static long kthTerm(int v, long targetK) {
int gap = v + 1;
long secondEven = 2L * (v + 1);
long seedLimit = 2L * gap + 1;
List<Long> seed = bruteUlamUntilValue(2, v, seedLimit);
Set<Long> members = new HashSet<>(seed);
byte[] bits = new byte[gap];
long[] prefixOnes = new long[gap + 1];
for (int t = 0; t < gap; t++) {
bits[t] = members.contains(2L * t + 1) ? (byte) 1 : 0;
prefixOnes[t + 1] = prefixOnes[t] + bits[t];
}
long oddLessThanSecondEven = prefixOnes[gap];
long secondEvenIndex = 2 + oddLessThanSecondEven;
// Extend bits via XOR recurrence
int stateCount = 1 << gap;
int[] seen = new int[stateCount];
Arrays.fill(seen, -1);
int state = 0;
for (int i = 0; i < gap; i++)
state = (state << 1) | bits[i];
seen[state] = 0;
int lowerMask = gap == 1 ? 0 : (1 << (gap - 1)) - 1;
List<Byte> bitList = new ArrayList<>();
for (byte b : bits)
bitList.add(b);
List<Long> poList = new ArrayList<>();
for (long p : prefixOnes)
poList.add(p);
int step = 0, cycleStart = -1, cycleEnd = -1;
while (true) {
int oldest = (state >> (gap - 1)) & 1, newest = state & 1, nextBit = oldest ^ newest;
bitList.add((byte) nextBit);
poList.add(poList.get(poList.size() - 1) + nextBit);
state = ((state & lowerMask) << 1) | nextBit;
step++;
int prev = seen[state];
if (prev >= 0) {
cycleStart = prev;
cycleEnd = step;
break;
}
seen[state] = step;
}
long cycleTStart = gap - 1 + cycleStart, cycleTEnd = gap - 1 + cycleEnd;
long cycleLen = cycleTEnd - cycleTStart;
long onesBeforeCycle = poList.get((int) (cycleTStart + 1));
long onesPerCycle = poList.get((int) (cycleTEnd + 1)) - onesBeforeCycle;
long[] cyclePrefixOnes = new long[(int) (cycleLen + 1)];
for (long i = 0; i < cycleLen; i++) {
long t = cycleTStart + 1 + i;
cyclePrefixOnes[(int) (i + 1)] = cyclePrefixOnes[(int) i] + bitList.get((int) t);
}
// kth term
if (targetK == 1)
return 2;
if (targetK == secondEvenIndex)
return secondEven;
long oddRank = targetK < secondEvenIndex ? targetK - 1 : targetK - 2;
// Find t such that prefix_ones[t+1] = oddRank
long t;
if (oddRank <= onesBeforeCycle) {
int lo = 0, hi = (int) (cycleTStart + 1);
while (lo < hi) {
int mid = (lo + hi) / 2;
if (poList.get(mid + 1) >= oddRank)
hi = mid;
else
lo = mid + 1;
}
t = lo;
} else {
long rem = oddRank - onesBeforeCycle;
long fullCycles = (rem - 1) / onesPerCycle;
long remInside = rem - fullCycles * onesPerCycle;
int lo = 0, hi = (int) cycleLen;
while (lo < hi) {
int mid = (lo + hi) / 2;
if (cyclePrefixOnes[mid + 1] >= remInside)
hi = mid;
else
lo = mid + 1;
}
t = cycleTStart + fullCycles * cycleLen + lo + 1;
}
return 2 * t + 1;
}
public static void main(String[] args) {
long targetK = 100000000000L;
long sum = 0;
for (int n = 2; n <= 10; n++) {
int v = 2 * n + 1;
sum += kthTerm(v, targetK);
}
System.out.println(sum);
}
}