Problem 361: Subsequence of Thue-Morse Sequence
View on Project EulerProject Euler Problem 361 Solution
EulerSolve provides an optimized solution for Project Euler Problem 361, Subsequence of Thue-Morse Sequence, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Let \(t=0110100110010110\ldots\) be the infinite Thue-Morse word, the fixed point of the morphism \(\mu:0\mapsto01,\;1\mapsto10\). For each length \(\ell\), let \(T_\ell\) be the set of distinct factors of \(t\) of length \(\ell\). Sort \(T_\ell\) lexicographically, keep only the words whose first bit is \(1\), and then concatenate these filtered lists for \(\ell=1,2,3,\ldots\). Reading every chosen word as a binary integer produces the sequence \(A(1),A(2),\ldots\). The goal is to compute $$\sum_{k=1}^{18} A(10^k) \pmod{10^9}.$$ Mathematical Approach 1. Count the Number of Factors of Each Length Write $$p(\ell)=|T_\ell|.$$ The solution files use the standard Thue-Morse factor-complexity recurrence with the explicit base values $$p(1)=2,\qquad p(2)=4,\qquad p(3)=6,$$ and for \(n\ge 2\), $$p(2n)=p(n)+p(n+1),\qquad p(2n+1)=2p(n+1).$$ Now define the cumulative count $$s(\ell)=\sum_{j=1}^{\ell} p(j).$$ The code memoizes the derived recurrences $$s(2n)=3s(n)+s(n+1)-8,\qquad s(2n+1)=4s(n)+3p(n+1)-8.$$ Why does this help with the Project Euler sequence? Because the Thue-Morse word is closed under bitwise complement: if \(w\in T_\ell\), then \(\overline{w}\in T_\ell\). No nonempty binary word equals its own complement, so the factors of length \(\ell\) are paired into one word starting with \(0\) and one starting with \(1\)....
Detailed mathematical approach
Problem Summary
Let \(t=0110100110010110\ldots\) be the infinite Thue-Morse word, the fixed point of the morphism \(\mu:0\mapsto01,\;1\mapsto10\). For each length \(\ell\), let \(T_\ell\) be the set of distinct factors of \(t\) of length \(\ell\). Sort \(T_\ell\) lexicographically, keep only the words whose first bit is \(1\), and then concatenate these filtered lists for \(\ell=1,2,3,\ldots\). Reading every chosen word as a binary integer produces the sequence \(A(1),A(2),\ldots\). The goal is to compute
$$\sum_{k=1}^{18} A(10^k) \pmod{10^9}.$$
Mathematical Approach
1. Count the Number of Factors of Each Length
Write
$$p(\ell)=|T_\ell|.$$
The solution files use the standard Thue-Morse factor-complexity recurrence with the explicit base values
$$p(1)=2,\qquad p(2)=4,\qquad p(3)=6,$$
and for \(n\ge 2\),
$$p(2n)=p(n)+p(n+1),\qquad p(2n+1)=2p(n+1).$$
Now define the cumulative count
$$s(\ell)=\sum_{j=1}^{\ell} p(j).$$
The code memoizes the derived recurrences
$$s(2n)=3s(n)+s(n+1)-8,\qquad s(2n+1)=4s(n)+3p(n+1)-8.$$
Why does this help with the Project Euler sequence? Because the Thue-Morse word is closed under bitwise complement: if \(w\in T_\ell\), then \(\overline{w}\in T_\ell\). No nonempty binary word equals its own complement, so the factors of length \(\ell\) are paired into one word starting with \(0\) and one starting with \(1\). Hence exactly half of them start with \(1\), and the number of valid Euler terms up to length \(\ell\) is
$$\operatorname{count\_upto}(\ell)=\frac{s(\ell)}{2}.$$
2. Convert a Global Index into a Length and a Local Rank
Given \(N\), the program first finds the smallest length \(\ell\) such that
$$\frac{s(\ell)}{2}\ge N.$$
This is done by doubling an upper bound and then applying binary search. Once \(\ell\) is known, the rank of the desired word inside the filtered list for this fixed length is
$$k=N-\operatorname{count\_upto}(\ell-1).$$
Inside the full lexicographic list of all words in \(T_\ell\), the leading-\(0\) half comes first and the leading-\(1\) half comes second. Therefore the wanted factor is the
$$\frac{p(\ell)}{2}+k$$
th factor of length \(\ell\) in complete lexicographic order. This is the reason for the line
$$\texttt{offset}=p(\ell)/2+k$$
in all three implementations.
3. Reconstruct the Lexicographic Factor Order Recursively
The difficult part is selecting the \(\texttt{offset}\)-th factor of length \(\ell\) without listing all factors. The implementations stop recursion at lengths \(1,2,3\), where the factor sets are stored explicitly:
$$T_1=\{0,1\},\qquad T_2=\{00,01,10,11\},\qquad T_3=\{001,010,011,100,101,110\}.$$
For \(\ell\ge 4\), the solution uses the unique reading frame induced by the length-2 morphism \(\mu\). Every factor either starts on a block boundary of \(\mu\) or one symbol later, and that leads to three lexicographic classes called \(O1\), \(E\), and \(O0\) in the source code.
4. Even Length Factors
Let \(\ell=2n\). If a factor starts on an even boundary, it is exactly \(\mu(u)\) for a unique \(u\in T_n\). This gives the middle class
$$E(u)=\mu(u),\qquad u\in T_n,$$
with size \(p(n)\).
If the factor starts one position later, it begins inside one \(\mu\)-block and ends inside another. For an ancestor \(a=a_1a_2\cdots a_{n+1}\in T_{n+1}\), the odd-start construction used by the code is
$$\operatorname{OE}(a)=(1-a_1)\,\mu(a_2a_3\cdots a_n)\,a_{n+1}.$$
If \(a_1=1\), the result starts with \(0\); if \(a_1=0\), the result starts with \(1\). Thus the complete lexicographic order for length \(2n\) splits into three consecutive blocks:
$$O1,\qquad E,\qquad O0,$$
with sizes
$$\frac{p(n+1)}{2},\qquad p(n),\qquad \frac{p(n+1)}{2}.$$
A concrete example is \(\ell=4\), where the ten factors are
$$0010,0011,0100 \mid 0101,0110,1001,1010 \mid 1011,1100,1101.$$
The first block is \(O1\), the middle block is \(E\), and the last block is \(O0\).
5. Odd Length Factors
Let \(\ell=2n+1\). There are again two alignment types, but now an odd-start factor crosses only the left block boundary, while an even-start factor crosses only the right one. For \(a=a_1a_2\cdots a_{n+1}\in T_{n+1}\), the code uses
$$\operatorname{OO}(a)=(1-a_1)\,\mu(a_2a_3\cdots a_{n+1}),$$
$$\operatorname{EO}(a)=\mu(a_1a_2\cdots a_n)\,a_{n+1}.$$
Exactly as in the even case, the lexicographic order is
$$O1,\qquad E,\qquad O0,$$
but now the block sizes are
$$\frac{p(n+1)}{2},\qquad p(n+1),\qquad \frac{p(n+1)}{2}.$$
This is why the recursive selector only needs the values \(p(n)\) and \(p(n+1)\) to decide which branch contains the required factor.
6. Evaluate the Selected Word Modulo \(10^9\) Without Expanding It
The chosen factor can be very long, so the program never builds its ordinary binary integer directly. Instead it defines
$$M=10^9,\qquad B_i=2^{2^i}\bmod M,$$
and lets \(V_i(w)\) be the value of a binary word \(w\) when read in base \(B_i\). In particular, \(V_0(w)\) is exactly the ordinary binary value modulo \(M\).
Concatenation satisfies the standard rule
$$V_i(uv)=V_i(u)\,B_i^{|v|}+V_i(v).$$
For the morphism, one block \(\mu(c)\) has value
$$V_i(\mu(c))=1+(B_i-1)c,\qquad c\in\{0,1\}.$$
So for \(w=c_1c_2\cdots c_m\),
$$V_i(\mu(w))=\sum_{j=1}^{m} B_i^{m-j}\bigl(1+(B_i-1)c_j\bigr).$$
Separating the constant part from the digit-dependent part gives
$$V_i(\mu(w))=G_{i+1}(m)+(B_i-1)V_{i+1}(w),$$
where
$$G_j(m)=\sum_{r=0}^{m-1} B_j^r.$$
This identity is exactly what the cached MU node evaluation computes. The helper
pow_sum stores both \(B_i^m\) and \(G_i(m)\), so each evaluation step remains small and reusable.
How the Code Works
The C++, Python, and Java versions all follow the same structure. They memoize p_count and
s_count, use find_length_for_index to locate the correct length, and then call
select_node to build only a symbolic expression DAG for the requested factor. WordBuilder
has four node types: empty, leaf, concatenation, and morphism application. The helper methods
prefix, suffix, and middle implement the ancestor-to-descendant maps behind
the \(O1/E/O0\) split. Finally get_eval computes the cached values \(V_i\) for each node, and
compute_A_mod returns \(V_0\). The C++ file also contains sanity checks proving that the construction
reproduces \(A(100)=3251\) and \(A(1000)=80852364498\).
Complexity Analysis
Because every recurrence halves its argument, computing \(p(\ell)\) or \(s(\ell)\) touches only \(O(\log \ell)\)
memoized states. The length search for \(A(N)\) costs \(O(\log N)\) calls to count_upto. After that,
select_node descends through \(O(\log \ell)\) recursive levels, and the node evaluator works with a
fixed base table of size 64. So each query is handled in polylogarithmic time with \(O(\log \ell)\) symbolic
memory, which is dramatically smaller than generating all factors or building the full binary integer.
References
- Problem page: https://projecteuler.net/problem=361
- Thue-Morse sequence: Wikipedia — Thue-Morse sequence
- Subword complexity: Wikipedia — Subword complexity
- J.-P. Allouche and J. Shallit, Automatic Sequences, Cambridge University Press, 2003.
Problem 361 source code
C++
#include <cassert>
#include <cstdint>
#include <future>
#include <iomanip>
#include <iostream>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>
using namespace std;
namespace {
constexpr uint64_t MOD = 1000000000ULL;
const vector<vector<string>> kBaseFactors = {
{},
{"0", "1"},
{"00", "01", "10", "11"},
{"001", "010", "011", "100", "101", "110"},
};
thread_local unordered_map<long long, unsigned long long> p_cache;
thread_local unordered_map<long long, unsigned long long> s_cache;
// For n>=4: p(2n)=p(n)+p(n+1), p(2n+1)=2*p(n+1).
// For n>=4: S(2n)=3*S(n)+S(n+1)-8, S(2n+1)=4*S(n)+3*p(n+1)-8.
unsigned long long p_count(long long n) {
if (n <= 0) return 0;
if (n == 1) return 2;
if (n == 2) return 4;
if (n == 3) return 6;
auto it = p_cache.find(n);
if (it != p_cache.end()) return it->second;
unsigned long long res = 0;
if ((n & 1LL) == 0) {
long long m = n / 2;
res = p_count(m) + p_count(m + 1);
} else {
long long m = n / 2;
res = 2ULL * p_count(m + 1);
}
p_cache[n] = res;
return res;
}
unsigned long long s_count(long long n) {
if (n <= 0) return 0;
if (n == 1) return 2;
if (n == 2) return 6;
if (n == 3) return 12;
auto it = s_cache.find(n);
if (it != s_cache.end()) return it->second;
unsigned long long res = 0;
if ((n & 1LL) == 0) {
long long m = n / 2;
__int128 tmp = 3;
tmp = tmp * s_count(m) + s_count(m + 1) - 8;
res = static_cast<unsigned long long>(tmp);
} else {
long long m = n / 2;
__int128 tmp = 4;
tmp = tmp * s_count(m) + 3 * static_cast<__int128>(p_count(m + 1)) - 8;
res = static_cast<unsigned long long>(tmp);
}
s_cache[n] = res;
return res;
}
unsigned long long count_upto(long long len) {
if (len <= 0) return 0;
return s_count(len) / 2ULL;
}
long long find_length_for_index(unsigned long long n) {
long long lo = 1;
long long hi = 1;
while (count_upto(hi) < n) {
hi *= 2;
}
while (lo < hi) {
long long mid = lo + (hi - lo) / 2;
if (count_upto(mid) >= n) {
hi = mid;
} else {
lo = mid + 1;
}
}
return lo;
}
struct WordBuilder {
struct Node {
enum Type { EMPTY, LEAF, CONCAT, MU } type = EMPTY;
int left = -1;
int right = -1;
int child = -1;
string bits;
long long len = 0;
int first = -1;
int last = -1;
};
uint64_t mod = MOD;
int max_i = 64;
vector<uint64_t> base;
vector<unordered_map<long long, pair<uint64_t, uint64_t>>> pow_sum_cache;
vector<Node> nodes;
vector<vector<uint64_t>> eval_cache;
unordered_map<uint64_t, int> prefix_cache;
unordered_map<uint64_t, int> suffix_cache;
unordered_map<uint64_t, int> middle_cache;
int empty_id = 0;
int leaf0_id = 0;
int leaf1_id = 0;
explicit WordBuilder(uint64_t mod_, int max_i_ = 64) : mod(mod_), max_i(max_i_) {
base.assign(max_i + 2, 0);
base[0] = 2 % mod;
for (int i = 1; i < static_cast<int>(base.size()); ++i) {
base[i] = mul_mod(base[i - 1], base[i - 1]);
}
pow_sum_cache.resize(base.size());
nodes.reserve(512);
eval_cache.reserve(512);
Node empty;
empty.type = Node::EMPTY;
empty.len = 0;
nodes.push_back(empty);
eval_cache.emplace_back();
empty_id = 0;
Node leaf0;
leaf0.type = Node::LEAF;
leaf0.bits = "0";
leaf0.len = 1;
leaf0.first = 0;
leaf0.last = 0;
nodes.push_back(leaf0);
eval_cache.emplace_back();
leaf0_id = static_cast<int>(nodes.size() - 1);
Node leaf1;
leaf1.type = Node::LEAF;
leaf1.bits = "1";
leaf1.len = 1;
leaf1.first = 1;
leaf1.last = 1;
nodes.push_back(leaf1);
eval_cache.emplace_back();
leaf1_id = static_cast<int>(nodes.size() - 1);
}
uint64_t mul_mod(uint64_t a, uint64_t b) const {
return static_cast<uint64_t>((__int128)a * b % mod);
}
int make_leaf(const string& bits) {
if (bits.empty()) return empty_id;
if (bits.size() == 1) return bits[0] == '0' ? leaf0_id : leaf1_id;
Node node;
node.type = Node::LEAF;
node.bits = bits;
node.len = static_cast<long long>(bits.size());
node.first = bits.front() - '0';
node.last = bits.back() - '0';
nodes.push_back(node);
eval_cache.emplace_back();
return static_cast<int>(nodes.size() - 1);
}
int make_bit(int b) { return b == 0 ? leaf0_id : leaf1_id; }
int make_concat(int left, int right) {
if (nodes[left].len == 0) return right;
if (nodes[right].len == 0) return left;
Node node;
node.type = Node::CONCAT;
node.left = left;
node.right = right;
node.len = nodes[left].len + nodes[right].len;
node.first = nodes[left].len ? nodes[left].first : nodes[right].first;
node.last = nodes[right].len ? nodes[right].last : nodes[left].last;
nodes.push_back(node);
eval_cache.emplace_back();
return static_cast<int>(nodes.size() - 1);
}
int make_mu(int child) {
if (nodes[child].len == 0) return empty_id;
Node node;
node.type = Node::MU;
node.child = child;
node.len = nodes[child].len * 2;
node.first = nodes[child].first;
node.last = 1 - nodes[child].last;
nodes.push_back(node);
eval_cache.emplace_back();
return static_cast<int>(nodes.size() - 1);
}
int prefix(int node_id) {
const Node& node = nodes[node_id];
if (node.len <= 1) return empty_id;
uint64_t key = (static_cast<uint64_t>(node_id) << 2) | 0ULL;
auto it = prefix_cache.find(key);
if (it != prefix_cache.end()) return it->second;
int res = empty_id;
if (node.type == Node::LEAF) {
res = make_leaf(node.bits.substr(0, static_cast<size_t>(node.len - 1)));
} else if (node.type == Node::CONCAT) {
const Node& right = nodes[node.right];
if (right.len > 1) {
res = make_concat(node.left, prefix(node.right));
} else if (right.len == 1) {
res = node.left;
} else {
res = prefix(node.left);
}
} else if (node.type == Node::MU) {
int p = prefix(node.child);
int mu_p = make_mu(p);
int bit = make_bit(nodes[node.child].last);
res = make_concat(mu_p, bit);
}
prefix_cache[key] = res;
return res;
}
int suffix(int node_id) {
const Node& node = nodes[node_id];
if (node.len <= 1) return empty_id;
uint64_t key = (static_cast<uint64_t>(node_id) << 2) | 1ULL;
auto it = suffix_cache.find(key);
if (it != suffix_cache.end()) return it->second;
int res = empty_id;
if (node.type == Node::LEAF) {
res = make_leaf(node.bits.substr(1));
} else if (node.type == Node::CONCAT) {
const Node& left = nodes[node.left];
if (left.len > 1) {
res = make_concat(suffix(node.left), node.right);
} else if (left.len == 1) {
res = node.right;
} else {
res = suffix(node.right);
}
} else if (node.type == Node::MU) {
int s = suffix(node.child);
int mu_s = make_mu(s);
int bit = make_bit(1 - nodes[node.child].first);
res = make_concat(bit, mu_s);
}
suffix_cache[key] = res;
return res;
}
int middle(int node_id) {
const Node& node = nodes[node_id];
if (node.len <= 2) return empty_id;
uint64_t key = (static_cast<uint64_t>(node_id) << 2) | 2ULL;
auto it = middle_cache.find(key);
if (it != middle_cache.end()) return it->second;
int res = empty_id;
if (node.type == Node::LEAF) {
res = make_leaf(node.bits.substr(1, static_cast<size_t>(node.len - 2)));
} else if (node.type == Node::CONCAT) {
int tmp = suffix(node_id);
res = prefix(tmp);
} else if (node.type == Node::MU) {
int m = middle(node.child);
int mu_m = make_mu(m);
int bit_left = make_bit(1 - nodes[node.child].first);
int bit_right = make_bit(nodes[node.child].last);
res = make_concat(bit_left, make_concat(mu_m, bit_right));
}
middle_cache[key] = res;
return res;
}
int odd_even_map(int ancestor) {
int mid = middle(ancestor);
int mu_mid = make_mu(mid);
int bit_left = make_bit(1 - nodes[ancestor].first);
int bit_right = make_bit(nodes[ancestor].last);
return make_concat(bit_left, make_concat(mu_mid, bit_right));
}
int odd_odd_map(int ancestor) {
int suf = suffix(ancestor);
int mu_suf = make_mu(suf);
int bit_left = make_bit(1 - nodes[ancestor].first);
return make_concat(bit_left, mu_suf);
}
int even_odd_map(int ancestor) {
int pref = prefix(ancestor);
int mu_pref = make_mu(pref);
int bit_right = make_bit(nodes[ancestor].last);
return make_concat(mu_pref, bit_right);
}
pair<uint64_t, uint64_t> pow_sum(int i, long long len) {
auto& cache = pow_sum_cache[i];
auto it = cache.find(len);
if (it != cache.end()) return it->second;
pair<uint64_t, uint64_t> res;
if (len == 0) {
res = {1 % mod, 0};
} else if (len == 1) {
res = {base[i] % mod, 1 % mod};
} else if ((len & 1LL) == 0) {
auto half = pow_sum(i, len / 2);
uint64_t p = half.first;
uint64_t s = half.second;
uint64_t p2 = mul_mod(p, p);
uint64_t s2 = mul_mod(s, (p + 1) % mod);
res = {p2, s2};
} else {
auto prev = pow_sum(i, len - 1);
uint64_t p = prev.first;
uint64_t s = prev.second;
uint64_t p2 = mul_mod(p, base[i]);
uint64_t s2 = (s + p) % mod;
res = {p2, s2};
}
cache[len] = res;
return res;
}
uint64_t pow_base(int i, long long len) {
return pow_sum(i, len).first;
}
uint64_t sum_geom(int i, long long len) {
return pow_sum(i, len).second;
}
const vector<uint64_t>& eval(int node_id) {
auto& cached = eval_cache[node_id];
if (!cached.empty()) return cached;
cached.assign(max_i + 2, 0);
const Node& node = nodes[node_id];
if (node.len == 0) {
return cached;
}
if (node.type == Node::LEAF) {
for (int i = 0; i < static_cast<int>(cached.size()); ++i) {
uint64_t v = 0;
for (char c : node.bits) {
v = (mul_mod(v, base[i]) + static_cast<uint64_t>(c - '0')) % mod;
}
cached[i] = v;
}
} else if (node.type == Node::CONCAT) {
const auto& left = eval(node.left);
const auto& right = eval(node.right);
long long len_r = nodes[node.right].len;
for (int i = 0; i < static_cast<int>(cached.size()); ++i) {
uint64_t term = mul_mod(left[i], pow_base(i, len_r));
cached[i] = (term + right[i]) % mod;
}
} else if (node.type == Node::MU) {
const auto& child = eval(node.child);
long long len_c = nodes[node.child].len;
for (int i = 0; i + 1 < static_cast<int>(cached.size()); ++i) {
uint64_t sum = sum_geom(i + 1, len_c);
uint64_t factor = (base[i] + mod - 1) % mod;
uint64_t term = mul_mod(factor, child[i + 1]);
cached[i] = (sum + term) % mod;
}
}
return cached;
}
string materialize(int node_id, size_t limit) {
const Node& node = nodes[node_id];
if (node.len > static_cast<long long>(limit)) return string();
if (node.type == Node::EMPTY) return string();
if (node.type == Node::LEAF) return node.bits;
if (node.type == Node::CONCAT) {
string left = materialize(node.left, limit);
if (left.empty() && nodes[node.left].len > 0) return string();
string right = materialize(node.right, limit);
if (right.empty() && nodes[node.right].len > 0) return string();
return left + right;
}
string child = materialize(node.child, limit);
if (child.empty() && nodes[node.child].len > 0) return string();
string out;
out.reserve(child.size() * 2);
for (char c : child) {
if (c == '0') out += "01";
else out += "10";
}
return out;
}
};
int select_node(long long len, unsigned long long k, WordBuilder& builder) {
if (len <= 3) {
return builder.make_leaf(kBaseFactors[static_cast<size_t>(len)][static_cast<size_t>(k - 1)]);
}
// Lexicographic split: O1 (odd-aligned from leading 1 ancestors), E (even-aligned), O0 (odd-aligned from leading 0 ancestors).
if ((len & 1LL) == 0) {
long long n = len / 2;
unsigned long long size_O1 = p_count(n + 1) / 2ULL;
unsigned long long size_E = p_count(n);
if (k <= size_O1) {
int ancestor = select_node(n + 1, p_count(n + 1) / 2ULL + k, builder);
return builder.odd_even_map(ancestor);
}
if (k <= size_O1 + size_E) {
int ancestor = select_node(n, k - size_O1, builder);
return builder.make_mu(ancestor);
}
int ancestor = select_node(n + 1, k - size_O1 - size_E, builder);
return builder.odd_even_map(ancestor);
}
long long n = len / 2;
unsigned long long size_O1 = p_count(n + 1) / 2ULL;
unsigned long long size_E = p_count(n + 1);
if (k <= size_O1) {
int ancestor = select_node(n + 1, p_count(n + 1) / 2ULL + k, builder);
return builder.odd_odd_map(ancestor);
}
if (k <= size_O1 + size_E) {
int ancestor = select_node(n + 1, k - size_O1, builder);
return builder.even_odd_map(ancestor);
}
int ancestor = select_node(n + 1, k - size_O1 - size_E, builder);
return builder.odd_odd_map(ancestor);
}
uint64_t compute_A_mod(unsigned long long n) {
if (n == 0) return 0;
long long len = find_length_for_index(n);
unsigned long long prev = count_upto(len - 1);
unsigned long long k = n - prev;
unsigned long long offset = p_count(len) / 2ULL + k;
WordBuilder builder(MOD);
int root = select_node(len, offset, builder);
const auto& vals = builder.eval(root);
return vals[0] % MOD;
}
void run_validation() {
assert(p_count(1) == 2);
assert(p_count(2) == 4);
assert(p_count(3) == 6);
assert(p_count(4) == 10);
assert(s_count(3) == 12);
assert(s_count(4) == 22);
{
unsigned long long n = 100;
long long len = find_length_for_index(n);
unsigned long long prev = count_upto(len - 1);
unsigned long long k = n - prev;
unsigned long long offset = p_count(len) / 2ULL + k;
WordBuilder builder(MOD);
int root = select_node(len, offset, builder);
string bits = builder.materialize(root, 128);
assert(!bits.empty());
unsigned long long value = 0;
for (char c : bits) value = (value << 1) + (c - '0');
assert(value == 3251ULL);
}
{
unsigned long long n = 1000;
long long len = find_length_for_index(n);
unsigned long long prev = count_upto(len - 1);
unsigned long long k = n - prev;
unsigned long long offset = p_count(len) / 2ULL + k;
WordBuilder builder(MOD);
int root = select_node(len, offset, builder);
string bits = builder.materialize(root, 256);
assert(!bits.empty());
unsigned long long value = 0;
for (char c : bits) value = (value << 1) + (c - '0');
assert(value == 80852364498ULL);
}
}
} // namespace
int main() {
run_validation();
vector<unsigned long long> targets;
unsigned long long v = 1;
for (int k = 1; k <= 18; ++k) {
v *= 10ULL;
targets.push_back(v);
}
vector<future<uint64_t>> futures;
futures.reserve(targets.size());
for (unsigned long long n : targets) {
futures.emplace_back(async(launch::async, [n]() { return compute_A_mod(n); }));
}
uint64_t sum_mod = 0;
for (auto& fut : futures) {
sum_mod = (sum_mod + fut.get()) % MOD;
}
cout << setw(9) << setfill('0') << sum_mod << "\n";
return 0;
}
Python
MOD = 1000000000
kBaseFactors = [
[],
["0", "1"],
["00", "01", "10", "11"],
["001", "010", "011", "100", "101", "110"]
]
p_cache = {}
s_cache = {}
def p_count(n):
if n <= 0: return 0
if n == 1: return 2
if n == 2: return 4
if n == 3: return 6
if n in p_cache: return p_cache[n]
if n % 2 == 0:
m = n // 2
res = p_count(m) + p_count(m + 1)
else:
m = n // 2
res = 2 * p_count(m + 1)
p_cache[n] = res
return res
def s_count(n):
if n <= 0: return 0
if n == 1: return 2
if n == 2: return 6
if n == 3: return 12
if n in s_cache: return s_cache[n]
if n % 2 == 0:
m = n // 2
res = 3 * s_count(m) + s_count(m + 1) - 8
else:
m = n // 2
res = 4 * s_count(m) + 3 * p_count(m + 1) - 8
s_cache[n] = res
return res
def count_upto(length):
if length <= 0: return 0
return s_count(length) // 2
def find_length_for_index(n):
lo = 1
hi = 1
while count_upto(hi) < n:
hi *= 2
while lo < hi:
mid = (lo + hi) // 2
if count_upto(mid) >= n:
hi = mid
else:
lo = mid + 1
return lo
EMPTY = 0
LEAF = 1
CONCAT = 2
MU = 3
class Node:
__slots__ = ['type', 'left', 'right', 'child', 'bits', 'len', 'first', 'last']
def __init__(self):
self.type = EMPTY
self.left = -1
self.right = -1
self.child = -1
self.bits = ""
self.len = 0
self.first = -1
self.last = -1
class WordBuilder:
def __init__(self, mod=MOD, max_i=64):
self.mod = mod
self.max_i = max_i
self.base = [0] * (max_i + 2)
self.base[0] = 2 % mod
for i in range(1, len(self.base)):
self.base[i] = (self.base[i - 1] * self.base[i - 1]) % mod
self.pow_sum_cache = [{} for _ in range(len(self.base))]
self.nodes = []
self.eval_cache = []
self.prefix_cache = {}
self.suffix_cache = {}
self.middle_cache = {}
empty = Node()
empty.type = EMPTY
empty.len = 0
self.nodes.append(empty)
self.eval_cache.append(None)
self.empty_id = 0
leaf0 = Node()
leaf0.type = LEAF
leaf0.bits = "0"
leaf0.len = 1
leaf0.first = 0
leaf0.last = 0
self.nodes.append(leaf0)
self.eval_cache.append(None)
self.leaf0_id = len(self.nodes) - 1
leaf1 = Node()
leaf1.type = LEAF
leaf1.bits = "1"
leaf1.len = 1
leaf1.first = 1
leaf1.last = 1
self.nodes.append(leaf1)
self.eval_cache.append(None)
self.leaf1_id = len(self.nodes) - 1
def make_leaf(self, bits):
if not bits: return self.empty_id
if len(bits) == 1: return self.leaf0_id if bits[0] == '0' else self.leaf1_id
node = Node()
node.type = LEAF
node.bits = bits
node.len = len(bits)
node.first = int(bits[0])
node.last = int(bits[-1])
self.nodes.append(node)
self.eval_cache.append(None)
return len(self.nodes) - 1
def make_bit(self, b):
return self.leaf0_id if b == 0 else self.leaf1_id
def make_concat(self, left, right):
if self.nodes[left].len == 0: return right
if self.nodes[right].len == 0: return left
node = Node()
node.type = CONCAT
node.left = left
node.right = right
node.len = self.nodes[left].len + self.nodes[right].len
node.first = self.nodes[left].first if self.nodes[left].len else self.nodes[right].first
node.last = self.nodes[right].last if self.nodes[right].len else self.nodes[left].last
self.nodes.append(node)
self.eval_cache.append(None)
return len(self.nodes) - 1
def make_mu(self, child):
if self.nodes[child].len == 0: return self.empty_id
node = Node()
node.type = MU
node.child = child
node.len = self.nodes[child].len * 2
node.first = self.nodes[child].first
node.last = 1 - self.nodes[child].last
self.nodes.append(node)
self.eval_cache.append(None)
return len(self.nodes) - 1
def prefix(self, node_id):
node = self.nodes[node_id]
if node.len <= 1: return self.empty_id
key = (node_id << 2) | 0
if key in self.prefix_cache: return self.prefix_cache[key]
res = self.empty_id
if node.type == LEAF:
res = self.make_leaf(node.bits[:-1])
elif node.type == CONCAT:
right_len = self.nodes[node.right].len
if right_len > 1:
res = self.make_concat(node.left, self.prefix(node.right))
elif right_len == 1:
res = node.left
else:
res = self.prefix(node.left)
elif node.type == MU:
p = self.prefix(node.child)
mu_p = self.make_mu(p)
bit = self.make_bit(self.nodes[node.child].last)
res = self.make_concat(mu_p, bit)
self.prefix_cache[key] = res
return res
def suffix(self, node_id):
node = self.nodes[node_id]
if node.len <= 1: return self.empty_id
key = (node_id << 2) | 1
if key in self.suffix_cache: return self.suffix_cache[key]
res = self.empty_id
if node.type == LEAF:
res = self.make_leaf(node.bits[1:])
elif node.type == CONCAT:
left_len = self.nodes[node.left].len
if left_len > 1:
res = self.make_concat(self.suffix(node.left), node.right)
elif left_len == 1:
res = node.right
else:
res = self.suffix(node.right)
elif node.type == MU:
s = self.suffix(node.child)
mu_s = self.make_mu(s)
bit = self.make_bit(1 - self.nodes[node.child].first)
res = self.make_concat(bit, mu_s)
self.suffix_cache[key] = res
return res
def middle(self, node_id):
node = self.nodes[node_id]
if node.len <= 2: return self.empty_id
key = (node_id << 2) | 2
if key in self.middle_cache: return self.middle_cache[key]
res = self.empty_id
if node.type == LEAF:
res = self.make_leaf(node.bits[1:-1])
elif node.type == CONCAT:
tmp = self.suffix(node_id)
res = self.prefix(tmp)
elif node.type == MU:
m = self.middle(node.child)
mu_m = self.make_mu(m)
bit_left = self.make_bit(1 - self.nodes[node.child].first)
bit_right = self.make_bit(self.nodes[node.child].last)
res = self.make_concat(bit_left, self.make_concat(mu_m, bit_right))
self.middle_cache[key] = res
return res
def odd_even_map(self, ancestor):
mid = self.middle(ancestor)
mu_mid = self.make_mu(mid)
bit_left = self.make_bit(1 - self.nodes[ancestor].first)
bit_right = self.make_bit(self.nodes[ancestor].last)
return self.make_concat(bit_left, self.make_concat(mu_mid, bit_right))
def odd_odd_map(self, ancestor):
suf = self.suffix(ancestor)
mu_suf = self.make_mu(suf)
bit_left = self.make_bit(1 - self.nodes[ancestor].first)
return self.make_concat(bit_left, mu_suf)
def even_odd_map(self, ancestor):
pref = self.prefix(ancestor)
mu_pref = self.make_mu(pref)
bit_right = self.make_bit(self.nodes[ancestor].last)
return self.make_concat(mu_pref, bit_right)
def pow_sum(self, i, length):
cache = self.pow_sum_cache[i]
if length in cache: return cache[length]
if length == 0:
res = (1 % self.mod, 0)
elif length == 1:
res = (self.base[i] % self.mod, 1 % self.mod)
elif length % 2 == 0:
p, s = self.pow_sum(i, length // 2)
p2 = (p * p) % self.mod
s2 = (s * (p + 1)) % self.mod
res = (p2, s2)
else:
p, s = self.pow_sum(i, length - 1)
p2 = (p * self.base[i]) % self.mod
s2 = (s + p) % self.mod
res = (p2, s2)
cache[length] = res
return res
def pow_base(self, i, length):
return self.pow_sum(i, length)[0]
def sum_geom(self, i, length):
return self.pow_sum(i, length)[1]
def get_eval(self, node_id):
cached = self.eval_cache[node_id]
if cached is not None: return cached
cached = [0] * (self.max_i + 2)
self.eval_cache[node_id] = cached
node = self.nodes[node_id]
if node.len == 0:
return cached
if node.type == LEAF:
for i in range(len(cached)):
v = 0
for c in node.bits:
v = ((v * self.base[i]) + int(c)) % self.mod
cached[i] = v
elif node.type == CONCAT:
left = self.get_eval(node.left)
right = self.get_eval(node.right)
len_r = self.nodes[node.right].len
for i in range(len(cached)):
term = (left[i] * self.pow_base(i, len_r)) % self.mod
cached[i] = (term + right[i]) % self.mod
elif node.type == MU:
child = self.get_eval(node.child)
len_c = self.nodes[node.child].len
for i in range(len(cached) - 1):
sum_g = self.sum_geom(i + 1, len_c)
factor = (self.base[i] + self.mod - 1) % self.mod
term = (factor * child[i + 1]) % self.mod
cached[i] = (sum_g + term) % self.mod
return cached
def select_node(length, k, builder):
if length <= 3:
return builder.make_leaf(kBaseFactors[length][k - 1])
if length % 2 == 0:
n = length // 2
size_O1 = p_count(n + 1) // 2
size_E = p_count(n)
if k <= size_O1:
ancestor = select_node(n + 1, p_count(n + 1) // 2 + k, builder)
return builder.odd_even_map(ancestor)
if k <= size_O1 + size_E:
ancestor = select_node(n, k - size_O1, builder)
return builder.make_mu(ancestor)
ancestor = select_node(n + 1, k - size_O1 - size_E, builder)
return builder.odd_even_map(ancestor)
n = length // 2
size_O1 = p_count(n + 1) // 2
size_E = p_count(n + 1)
if k <= size_O1:
ancestor = select_node(n + 1, p_count(n + 1) // 2 + k, builder)
return builder.odd_odd_map(ancestor)
if k <= size_O1 + size_E:
ancestor = select_node(n + 1, k - size_O1, builder)
return builder.even_odd_map(ancestor)
ancestor = select_node(n + 1, k - size_O1 - size_E, builder)
return builder.odd_odd_map(ancestor)
def compute_A_mod(n):
if n == 0: return 0
length = find_length_for_index(n)
prev = count_upto(length - 1)
k = n - prev
offset = p_count(length) // 2 + k
builder = WordBuilder(MOD)
root = select_node(length, offset, builder)
vals = builder.get_eval(root)
return vals[0] % MOD
def solve():
targets = [10**k for k in range(1, 19)]
ans = 0
for target in targets:
ans = (ans + compute_A_mod(target)) % MOD
return "{:09d}".format(ans)
if __name__ == '__main__':
print(solve())
Java
import java.util.*;
public class Euler361 {
static final long MOD = 1000000000L;
static final String[][] kBaseFactors = {
{},
{ "0", "1" },
{ "00", "01", "10", "11" },
{ "001", "010", "011", "100", "101", "110" }
};
static Map<Long, Long> pCache = new HashMap<>();
static Map<Long, Long> sCache = new HashMap<>();
static long pCount(long n) {
if (n <= 0)
return 0;
if (n == 1)
return 2;
if (n == 2)
return 4;
if (n == 3)
return 6;
if (pCache.containsKey(n))
return pCache.get(n);
long res = 0;
if (n % 2 == 0) {
long m = n / 2;
res = pCount(m) + pCount(m + 1);
} else {
long m = n / 2;
res = 2L * pCount(m + 1);
}
pCache.put(n, res);
return res;
}
static long sCount(long n) {
if (n <= 0)
return 0;
if (n == 1)
return 2;
if (n == 2)
return 6;
if (n == 3)
return 12;
if (sCache.containsKey(n))
return sCache.get(n);
long res = 0;
if (n % 2 == 0) {
long m = n / 2;
res = 3L * sCount(m) + sCount(m + 1) - 8L;
} else {
long m = n / 2;
res = 4L * sCount(m) + 3L * pCount(m + 1) - 8L;
}
sCache.put(n, res);
return res;
}
static long countUpto(long len) {
if (len <= 0)
return 0;
return sCount(len) / 2L;
}
static long findLengthForIndex(long n) {
long lo = 1;
long hi = 1;
while (countUpto(hi) < n)
hi *= 2L;
while (lo < hi) {
long mid = lo + (hi - lo) / 2L;
if (countUpto(mid) >= n)
hi = mid;
else
lo = mid + 1L;
}
return lo;
}
static class Node {
static final int EMPTY = 0;
static final int LEAF = 1;
static final int CONCAT = 2;
static final int MU = 3;
int type = EMPTY;
int left = -1;
int right = -1;
int child = -1;
String bits = "";
long len = 0;
int first = -1;
int last = -1;
}
static class WordBuilder {
long mod;
int maxI;
long[] base;
Map<Long, long[]>[] powSumCache;
List<Node> nodes;
List<long[]> evalCache;
Map<Long, Integer> prefixCache;
Map<Long, Integer> suffixCache;
Map<Long, Integer> middleCache;
int emptyId = 0;
int leaf0Id = 0;
int leaf1Id = 0;
@SuppressWarnings("unchecked")
WordBuilder(long mod, int maxI) {
this.mod = mod;
this.maxI = maxI;
base = new long[maxI + 2];
base[0] = 2L % mod;
for (int i = 1; i < base.length; i++) {
base[i] = (base[i - 1] * base[i - 1]) % mod;
}
powSumCache = new Map[base.length];
for (int i = 0; i < base.length; i++)
powSumCache[i] = new HashMap<>();
nodes = new ArrayList<>(512);
evalCache = new ArrayList<>(512);
prefixCache = new HashMap<>();
suffixCache = new HashMap<>();
middleCache = new HashMap<>();
Node empty = new Node();
empty.type = Node.EMPTY;
empty.len = 0;
nodes.add(empty);
evalCache.add(null);
emptyId = 0;
Node leaf0 = new Node();
leaf0.type = Node.LEAF;
leaf0.bits = "0";
leaf0.len = 1;
leaf0.first = 0;
leaf0.last = 0;
nodes.add(leaf0);
evalCache.add(null);
leaf0Id = nodes.size() - 1;
Node leaf1 = new Node();
leaf1.type = Node.LEAF;
leaf1.bits = "1";
leaf1.len = 1;
leaf1.first = 1;
leaf1.last = 1;
nodes.add(leaf1);
evalCache.add(null);
leaf1Id = nodes.size() - 1;
}
long mulMod(long a, long b) {
return (a * b) % mod;
}
int makeLeaf(String bits) {
if (bits.isEmpty())
return emptyId;
if (bits.length() == 1)
return bits.equals("0") ? leaf0Id : leaf1Id;
Node node = new Node();
node.type = Node.LEAF;
node.bits = bits;
node.len = bits.length();
node.first = bits.charAt(0) - '0';
node.last = bits.charAt(bits.length() - 1) - '0';
nodes.add(node);
evalCache.add(null);
return nodes.size() - 1;
}
int makeBit(int b) {
return b == 0 ? leaf0Id : leaf1Id;
}
int makeConcat(int left, int right) {
if (nodes.get(left).len == 0)
return right;
if (nodes.get(right).len == 0)
return left;
Node node = new Node();
node.type = Node.CONCAT;
node.left = left;
node.right = right;
node.len = nodes.get(left).len + nodes.get(right).len;
node.first = nodes.get(left).len > 0 ? nodes.get(left).first : nodes.get(right).first;
node.last = nodes.get(right).len > 0 ? nodes.get(right).last : nodes.get(left).last;
nodes.add(node);
evalCache.add(null);
return nodes.size() - 1;
}
int makeMu(int child) {
if (nodes.get(child).len == 0)
return emptyId;
Node node = new Node();
node.type = Node.MU;
node.child = child;
node.len = nodes.get(child).len * 2L;
node.first = nodes.get(child).first;
node.last = 1 - nodes.get(child).last;
nodes.add(node);
evalCache.add(null);
return nodes.size() - 1;
}
int prefix(int nodeId) {
Node node = nodes.get(nodeId);
if (node.len <= 1)
return emptyId;
long key = ((long) nodeId << 2) | 0L;
if (prefixCache.containsKey(key))
return prefixCache.get(key);
int res = emptyId;
if (node.type == Node.LEAF) {
res = makeLeaf(node.bits.substring(0, (int) (node.len - 1)));
} else if (node.type == Node.CONCAT) {
Node right = nodes.get(node.right);
if (right.len > 1) {
res = makeConcat(node.left, prefix(node.right));
} else if (right.len == 1) {
res = node.left;
} else {
res = prefix(node.left);
}
} else if (node.type == Node.MU) {
int p = prefix(node.child);
int muP = makeMu(p);
int bit = makeBit(nodes.get(node.child).last);
res = makeConcat(muP, bit);
}
prefixCache.put(key, res);
return res;
}
int suffix(int nodeId) {
Node node = nodes.get(nodeId);
if (node.len <= 1)
return emptyId;
long key = ((long) nodeId << 2) | 1L;
if (suffixCache.containsKey(key))
return suffixCache.get(key);
int res = emptyId;
if (node.type == Node.LEAF) {
res = makeLeaf(node.bits.substring(1));
} else if (node.type == Node.CONCAT) {
Node left = nodes.get(node.left);
if (left.len > 1) {
res = makeConcat(suffix(node.left), node.right);
} else if (left.len == 1) {
res = node.right;
} else {
res = suffix(node.right);
}
} else if (node.type == Node.MU) {
int s = suffix(node.child);
int muS = makeMu(s);
int bit = makeBit(1 - nodes.get(node.child).first);
res = makeConcat(bit, muS);
}
suffixCache.put(key, res);
return res;
}
int middle(int nodeId) {
Node node = nodes.get(nodeId);
if (node.len <= 2)
return emptyId;
long key = ((long) nodeId << 2) | 2L;
if (middleCache.containsKey(key))
return middleCache.get(key);
int res = emptyId;
if (node.type == Node.LEAF) {
res = makeLeaf(node.bits.substring(1, (int) (node.len - 1)));
} else if (node.type == Node.CONCAT) {
int tmp = suffix(nodeId);
res = prefix(tmp);
} else if (node.type == Node.MU) {
int m = middle(node.child);
int muM = makeMu(m);
int bitLeft = makeBit(1 - nodes.get(node.child).first);
int bitRight = makeBit(nodes.get(node.child).last);
res = makeConcat(bitLeft, makeConcat(muM, bitRight));
}
middleCache.put(key, res);
return res;
}
int oddEvenMap(int ancestor) {
int mid = middle(ancestor);
int muMid = makeMu(mid);
int bitLeft = makeBit(1 - nodes.get(ancestor).first);
int bitRight = makeBit(nodes.get(ancestor).last);
return makeConcat(bitLeft, makeConcat(muMid, bitRight));
}
int oddOddMap(int ancestor) {
int suf = suffix(ancestor);
int muSuf = makeMu(suf);
int bitLeft = makeBit(1 - nodes.get(ancestor).first);
return makeConcat(bitLeft, muSuf);
}
int evenOddMap(int ancestor) {
int pref = prefix(ancestor);
int muPref = makeMu(pref);
int bitRight = makeBit(nodes.get(ancestor).last);
return makeConcat(muPref, bitRight);
}
long[] powSum(int i, long len) {
Map<Long, long[]> cache = powSumCache[i];
if (cache.containsKey(len))
return cache.get(len);
long[] res = new long[2];
if (len == 0) {
res[0] = 1L % mod;
res[1] = 0;
} else if (len == 1) {
res[0] = base[i] % mod;
res[1] = 1L % mod;
} else if (len % 2 == 0) {
long[] half = powSum(i, len / 2L);
long p = half[0];
long s = half[1];
long p2 = mulMod(p, p);
long s2 = mulMod(s, (p + 1L) % mod);
res[0] = p2;
res[1] = s2;
} else {
long[] prev = powSum(i, len - 1L);
long p = prev[0];
long s = prev[1];
long p2 = mulMod(p, base[i]);
long s2 = (s + p) % mod;
res[0] = p2;
res[1] = s2;
}
cache.put(len, res);
return res;
}
long powBase(int i, long len) {
return powSum(i, len)[0];
}
long sumGeom(int i, long len) {
return powSum(i, len)[1];
}
long[] getEval(int nodeId) {
long[] cached = evalCache.get(nodeId);
if (cached != null)
return cached;
cached = new long[maxI + 2];
evalCache.set(nodeId, cached);
Node node = nodes.get(nodeId);
if (node.len == 0)
return cached;
if (node.type == Node.LEAF) {
for (int i = 0; i < cached.length; i++) {
long v = 0;
for (int c = 0; c < node.bits.length(); c++) {
v = (mulMod(v, base[i]) + (node.bits.charAt(c) - '0')) % mod;
}
cached[i] = v;
}
} else if (node.type == Node.CONCAT) {
long[] leftEval = getEval(node.left);
long[] rightEval = getEval(node.right);
long lenR = nodes.get(node.right).len;
for (int i = 0; i < cached.length; i++) {
long term = mulMod(leftEval[i], powBase(i, lenR));
cached[i] = (term + rightEval[i]) % mod;
}
} else if (node.type == Node.MU) {
long[] childEval = getEval(node.child);
long lenC = nodes.get(node.child).len;
for (int i = 0; i < cached.length - 1; i++) {
long sumG = sumGeom(i + 1, lenC);
long factor = (base[i] + mod - 1L) % mod;
long term = mulMod(factor, childEval[i + 1]);
cached[i] = (sumG + term) % mod;
}
}
return cached;
}
}
static int selectNode(long len, long k, WordBuilder builder) {
if (len <= 3) {
return builder.makeLeaf(kBaseFactors[(int) len][(int) (k - 1L)]);
}
if (len % 2 == 0) {
long n = len / 2L;
long sizeO1 = pCount(n + 1L) / 2L;
long sizeE = pCount(n);
if (k <= sizeO1) {
int ancestor = selectNode(n + 1L, pCount(n + 1L) / 2L + k, builder);
return builder.oddEvenMap(ancestor);
}
if (k <= sizeO1 + sizeE) {
int ancestor = selectNode(n, k - sizeO1, builder);
return builder.makeMu(ancestor);
}
int ancestor = selectNode(n + 1L, k - sizeO1 - sizeE, builder);
return builder.oddEvenMap(ancestor);
} else {
long n = len / 2L;
long sizeO1 = pCount(n + 1L) / 2L;
long sizeE = pCount(n + 1L);
if (k <= sizeO1) {
int ancestor = selectNode(n + 1L, pCount(n + 1L) / 2L + k, builder);
return builder.oddOddMap(ancestor);
}
if (k <= sizeO1 + sizeE) {
int ancestor = selectNode(n + 1L, k - sizeO1, builder);
return builder.evenOddMap(ancestor);
}
int ancestor = selectNode(n + 1L, k - sizeO1 - sizeE, builder);
return builder.oddOddMap(ancestor);
}
}
static long computeAMod(long n) {
if (n == 0)
return 0;
long len = findLengthForIndex(n);
long prev = countUpto(len - 1L);
long k = n - prev;
long offset = pCount(len) / 2L + k;
WordBuilder builder = new WordBuilder(MOD, 64);
int root = selectNode(len, offset, builder);
long[] vals = builder.getEval(root);
return vals[0] % MOD;
}
static String solve() {
long[] targets = new long[18];
long v = 1;
for (int i = 0; i < 18; i++) {
v *= 10;
targets[i] = v;
}
long sumMod = 0;
for (long tgt : targets) {
sumMod = (sumMod + computeAMod(tgt)) % MOD;
}
return String.format("%09d", sumMod);
}
public static void main(String[] args) {
System.out.println(solve());
}
}