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
The Thue-Morse sequence \(\{T_n\}\) is the binary sequence defined by \(T_0=0\), \(T_{2n}=T_n\) and \(T_{2n+1}=1-T_n\); it begins \(0110100110010110\ldots\). Let \(\{A_n\}\) be the sorted sequence of all integers whose binary expression occurs inside \(\{T_n\}\). For instance \(18\) is \(10010\) in binary, which occurs as \(T_8\) to \(T_{12}\), so \(18\) is an element, while \(14\) is \(1110\), which never occurs. The statement lists \(A_0\) to \(A_{12}\) as
$$0,1,2,3,4,5,6,9,10,11,12,13,18,$$
and gives the checkpoints \(A_{100}=3251\) and \(A_{1000}=80852364498\). We must find the last nine digits of
$$\sum_{k=1}^{18}A_{10^k},$$
that is, this sum modulo \(10^9\). [1] The official text says "as a subsequence", but its examples fix the meaning as consecutive positions: with gaps allowed, \(1,1,1,0\) would occur and \(14\) would not be excluded. We call a word a block when it occurs in consecutive positions.
The whole solution in one sentence: since the elements are sorted by bit length, we count the blocks of each length to find the length of \(A_N\), pick the required block by a recursion that halves that length at every step, and evaluate it modulo \(10^9\) from its compact description.
Mathematical Approach
1. The word, its morphism and one small lemma
Write \(t=t_0t_1t_2\ldots\) for the infinite word \(T_0T_1T_2\ldots\). The two recurrences make \(t_n\) the parity of the number of ones in the binary expansion of \(n\), [2] and they say something more useful as well. Let \(\mu\) be the substitution
$$\mu(0)=01,\qquad \mu(1)=10,$$
extended to words letter by letter. Then \(\mu(t)=t\): the pair \(\mu\) produces from \(t_n\) is \(t_n\,(1-t_n)\), which is exactly \(t_{2n}t_{2n+1}\). So \(t\) is cut into two-letter frame groups \(\mu(t_0),\mu(t_1),\ldots\) sitting at positions \(0,2,4,\ldots\).
Lemma 1. Neither \(000\) nor \(111\) is a block. Proof. Any three consecutive positions \(i,i+1,i+2\) contain a whole frame group: if \(i\) is even the group is \(\{i,i+1\}\), otherwise \(\{i+1,i+2\}\). The two letters of a group are \(t_k\) and \(1-t_k\), which differ, so three consecutive letters are never all equal. \(\square\) Much stronger facts are known — \(t\) contains no cube \(xxx\) and no overlap \(cxcxc\) — but Lemma 1 is all we need. [2]
2. Counting the blocks of each length
Let \(p(\ell)\) be the number of distinct blocks of length \(\ell\), the complexity of the word. [6] The first values follow by hand. Both letters occur, so \(p(1)=2\); all four pairs occur, so \(p(2)=4\); of the eight words of length \(3\), Lemma 1 removes \(000\) and \(111\) and the other six all occur, so \(p(3)=6\), namely \(001,010,011,100,101,110\). At length \(4\) exactly ten words occur,
$$\begin{aligned}&0010,0011,0100,0101,0110,\\ &1001,1010,1011,1100,1101.\end{aligned}$$
Lemma 2. A word \(w\) is a block if and only if its complement \(\overline{w}\) is. Proof. Since \(\mu(\overline{c})=\overline{\mu(c)}\) for a letter, \(\mu(\overline{v})=\overline{\mu(v)}\) for every word, hence \(\mu^k(1)=\overline{\mu^k(0)}\) for every \(k\). Every block \(w\) occurs inside some \(\mu^k(0)\), because these are the prefixes of \(t\) of length \(2^k\), so \(\overline{w}\) occurs inside \(\mu^k(1)\). And \(\mu^k(1)\) is a block: the letter \(1\) is one, and if \(v\) occurs in \(t\) then \(\mu(v)\) occurs in \(\mu(t)=t\). \(\square\)
No nonempty word equals its own complement, so the blocks of length \(\ell\) split into complementary pairs, one member starting with \(0\) and one with \(1\). Hence \(p(\ell)\) is even and exactly \(p(\ell)/2\) blocks of each length begin with \(1\).
3. From integers to blocks
The binary expression of a positive integer has no leading zero, so it is a word beginning with \(1\). The positive elements of \(\{A_n\}\) with exactly \(\ell\) bits are therefore the blocks of length \(\ell\) beginning with \(1\), and there are \(p(\ell)/2\) of them. Numbers with fewer bits are smaller, and among numbers of equal bit length the numerical order is the lexicographic order of the digit strings, so \(\{A_n\}\) is \(A_0=0\) followed by those blocks, grouped by increasing length and sorted inside each group. Writing
$$s(\ell)=\sum_{j=1}^{\ell}p(j),\qquad c(\ell)=\frac{s(\ell)}{2},$$
there are \(c(\ell)\) positive elements of at most \(\ell\) bits. From \(p(1),\ldots,p(6)=2,4,6,10,12,16\) we get \(s(1),\ldots,s(6)=2,6,12,22,34,50\) and
$$c(1),\ldots,c(6)=1,3,6,11,17,25.$$
This explains the statement's table: with \(c(4)=11\) elements of at most four bits, \(A_{11}=13\) is the largest of them, and \(A_{12}\) is the smallest length-\(5\) block beginning with \(1\), which is \(10010\), so \(A_{12}=18\).
One remark turns a rank in \(\{A_n\}\) into a rank in the full list of blocks: lexicographically the \(p(\ell)/2\) blocks beginning with \(0\) come first, so if \(A_N\) has \(\ell\) bits and is the \(k\)-th \(\ell\)-bit element, the block we want sits at position \(j=p(\ell)/2+k\) of the full length-\(\ell\) list. Everything reduces to one question: what is the \(j\)-th block of length \(\ell\)?
4. The reading frame
A block of length \(\ell\) is a window on \(t\), starting either in step with the frame groups or one letter later. Both cases describe it through a shorter block, its ancestor. For instance the block \(010\) expands to \(\mu(010)=01\,10\,01\); the length-\(4\) window at a group boundary is \(0110=\mu(01)\), and the window one letter later is \(1100\), which keeps half of the first group and half of the last. For even \(\ell=2n\), an even start covers whole groups, so the window is \(\mu(u)\) for a block \(u\) of length \(n\); conversely \(\mu(u)\) is a block whenever \(u\) is one, which gives \(p(n)\) words. An odd start takes the second letter of one group, then \(n-1\) whole groups, then the first letter of one more; if \(a=a_1\cdots a_{n+1}\) is the block covered by those groups, the second letter of the first group is \(1-a_1\) and the first letter of the last is \(a_{n+1}\), so the window is
$$\mathrm{OE}(a)=(1-a_1)\,\mu(a_2\cdots a_n)\,a_{n+1}.$$
For odd \(\ell=2n+1\) the same bookkeeping gives two maps from blocks \(a\) of length \(n+1\),
$$\mathrm{EO}(a)=\mu(a_1\cdots a_n)\,a_{n+1},$$
$$\mathrm{OO}(a)=(1-a_1)\,\mu(a_2\cdots a_{n+1}),$$
for even and odd starts. All four maps are injective: the ancestor is read back from the leading letter, the groups and the trailing letter.
The two families of one length are disjoint once \(\ell\ge4\). Suppose \(w\) of length \(2n\), \(n\ge2\), were in both. From \(w=\mu(u)\) we get \(w_{2i-1}\ne w_{2i}\); from \(w=\mathrm{OE}(a)\) the copy of \(\mu(a_2\cdots a_n)\) occupies positions \(2\) to \(2n-1\), giving \(w_{2i}\ne w_{2i+1}\). All neighbours then differ, so \(w\) alternates, and reading the ancestor back from an alternating \(w\) yields \(a=11\cdots1\) or \(a=00\cdots0\) of length \(n+1\ge3\), which Lemma 1 forbids. In the odd case the two patterns again force \(w\) to alternate, and its even-start ancestor has \(a_j=w_{2j-1}\) for \(j\le n\) and \(a_{n+1}=w_{2n+1}\), all odd positions of an alternating word and hence all equal, so \(a\) is constant again. Counting both families gives, for \(n\ge2\),
$$\begin{aligned}p(2n)&=p(n)+p(n+1),\\ p(2n+1)&=2p(n+1).\end{aligned}$$
The restriction \(n\ge2\) matters: for \(n=1\) the words \(010\) and \(101\) are both even- and odd-aligned, so \(2p(2)=8\) would overcount \(p(3)=6\), and \(p(1),p(2),p(3)\) are therefore the base of the recursion. The resulting \(2,4,6,10,12,16,20,22,\ldots\) is the known complexity of the Thue-Morse word, [3] [5] with the closed form \(p(\ell)=3(\ell-1)+d(\ell-1)\) for \(\ell\ge3\), where \(d(m)\) is the distance from \(m\) to the nearest power of two; [4] the two values \(p(1)=2\) and \(p(2)=4\) are exceptions to it. On \(2^r\le m\le2^{r+1}\) the distance \(d\) rises to \(2^{r-1}\) at \(m=3\cdot2^{r-1}\) and falls back, so \(0\le d(m)\le m/3\) and
$$3(\ell-1)\le p(\ell)\le\tfrac{10}{3}(\ell-1)\quad(\ell\ge3).$$
5. The lexicographic order of the families
Split the odd-start family by the first letter of its ancestor. The three parts then come in a fixed order, and the whole selection rests on that.
Lemma 3. Let \(\ell\ge4\). An odd-start word whose ancestor begins with \(1\) begins with \(001\) or \(0100\); an even-start word begins with \(0101\), \(0110\), \(1001\) or \(1010\); an odd-start word whose ancestor begins with \(0\) begins with \(1011\) or \(110\). Proof. An even-start word is \(\mu(u_1u_2\cdots)\), whose first four letters are \(\mu(u_1)\mu(u_2)\), one of the four patterns listed. For an odd start with \(a_1=1\) the first letter is \(0\) and the next group is \(\mu(a_2)\): if \(a_2=0\) the word begins \(001\); if \(a_2=1\) it begins \(010\) and the next letter is \(a_3=0\), because \(a\) has no \(111\), so it begins \(0100\). The case \(a_1=0\) is the complement. \(\square\)
Because
$$\begin{aligned}001&\lt 0100\lt 0101\lt 0110\\ &\lt 1001\lt 1010\lt 1011\lt 110,\end{aligned}$$
the sorted list of all blocks of length \(\ell\ge4\) is the three parts written one after another in that order. Each outer part has \(p(n+1)/2\) members, and the middle part has \(p(n)\) for even \(\ell\) and \(p(n+1)\) for odd \(\ell\). At \(\ell=4\) the list of Section 2 splits exactly so: \(0010,0011,0100\), then \(0101,0110,1001,1010\), then \(1011,1100,1101\).
Inside a part the order follows the order of the ancestors, because \(\mu\) preserves lexicographic order: two words that first differ at position \(i\) have images that first differ inside the \(i\)-th group, where \(01\lt10\). Within one part the first ancestor letter is fixed, so a first difference between two ancestors falls either inside a substituted group, where that group's first letter decides, or in the letter which \(\mathrm{OE}\) and \(\mathrm{EO}\) append unchanged. Each map is therefore increasing on its part, and the \(i\)-th word of a part comes from the \(i\)-th ancestor of its kind.
6. The selection recursion
To find the \(j\)-th block of length \(\ell\ge4\), put \(n=\lfloor \ell/2\rfloor\) and compare \(j\) with the three sizes. In the first part the ancestor is the \(j\)-th length-\((n+1)\) block beginning with \(1\), at position \(p(n+1)/2+j\) of the full list; in the last part it is at position \(j\) minus the first two sizes, which automatically begins with \(0\); in both we recurse and apply \(\mathrm{OE}\) or \(\mathrm{OO}\). In the middle part we recurse at length \(n\) for even \(\ell\) or \(n+1\) for odd \(\ell\) and apply \(\mu\) or \(\mathrm{EO}\). Lengths \(1,2,3\) come from the explicit lists, and each step roughly halves the length, so there are about \(\log_2\ell\) steps.
For \(N=12\) we saw \(c(4)=11\lt12\le c(5)=17\), so \(\ell=5\), \(k=1\) and \(j=p(5)/2+1=7\). Here \(\ell\) is odd with \(n=2\) and the sizes are \(3,6,3\); since \(3\lt7\le9\) the middle part applies at position \(7-3=4\), so the ancestor is the fourth length-\(3\) block, \(100\), and
$$\mathrm{EO}(100)=\mu(10)\,0=1001\,0=10010,$$
which is \(18=A_{12}\), as the statement says.
The first checkpoint exercises both outer parts. From \(c(11)=88\) and \(c(12)=106\) we get \(\ell=12\), \(k=12\) and \(j=30\); the sizes \(10,16,10\) put \(j\) in the last part at position \(4\), the sizes \(5,10,5\) at \(\ell=7\) put that in the first part at ancestor rank \(5+4=9\), and the sizes \(3,4,3\) at \(\ell=4\) put that in the last part at position \(2\), giving \(010\). Building the images back up,
$$\begin{aligned}\mathrm{OE}(010)&=1\,\mu(1)\,0=1100,\\ \mathrm{OO}(1100)&=0\,\mu(100)=0100101,\\ \mathrm{OE}(0100101)&=1\,\mu(10010)\,1,\end{aligned}$$
and the last line is \(110010110011=3251=A_{100}\). The same descent at \(N=1000\) ends at a block of \(37\) bits whose value is \(80852364498\), the second checkpoint.
7. The cumulative counts
The ranks need \(c(\ell)=s(\ell)/2\) for enormous \(\ell\), so \(s\) needs a halving recursion too. Split off the three base values and apply Section 4 to the rest:
$$\begin{aligned}\sum_{j=4}^{2n}p(j)&=\sum_{m=2}^{n}\bigl(p(m)+p(m+1)\bigr)\\ &\quad+\sum_{m=2}^{n-1}2p(m+1).\end{aligned}$$
The three sums are \(s(n)-2\), \(s(n+1)-6\) and \(2s(n)-12\), because \(p(1)=2\) and \(p(1)+p(2)=6\). Adding \(s(3)=12\) gives the first identity below; adding \(p(2n+1)=2p(n+1)\) and using \(s(n+1)=s(n)+p(n+1)\) gives the second:
$$\begin{aligned}s(2n)&=3s(n)+s(n+1)-8,\\ s(2n+1)&=4s(n)+3p(n+1)-8,\end{aligned}$$
valid for \(n\ge2\) with \(s(1)=2\), \(s(2)=6\), \(s(3)=12\). Like \(p\), each value depends on two values at roughly half the argument. Halving a window of consecutive integers and adding \(1\) leaves a window of at most three, so only three arguments are reachable at each depth and a memoised evaluation of one count touches \(O(\log \ell)\) pairs of function and argument.
8. Finding the length of \(A_N\)
Since \(c\) is strictly increasing, the bit length of \(A_N\) is the least \(\ell\) with \(c(\ell)\ge N\), and then \(k=N-c(\ell-1)\). The program doubles a bound until \(c\) passes \(N\), then bisects. The growth is quadratic. Summing the bounds of Section 4 and adding \(p(1)+p(2)=6\) gives, for \(\ell\ge2\),
$$\begin{aligned}\tfrac32\ell(\ell-1)+3&\le s(\ell),\\ s(\ell)&\le\tfrac53\ell(\ell-1)+\tfrac83,\end{aligned}$$
so \(c(\ell)=\Theta(\ell^2)\) and the bit length grows like \(\sqrt N\); measured, \(c(\ell)/\ell^2\) stays between \(0.785\) and \(0.800\) for \(\ell\le2\cdot10^5\), and \(\ell\) is close to \(1.12\sqrt N\). The search itself uses the exact counts, not this estimate. For \(N=10^{18}\) it returns \(\ell=1\,126\,216\,794\), a window of more than a billion letters: writing it out would cost time and space proportional to that length, which is exactly what the expression avoids.
9. Evaluating a block without expanding it
Let \(M=10^9\) and build a tower of bases
$$\begin{aligned}B_i&=2^{2^i}\bmod M,\\ B_0,B_1,B_2&=2,4,16.\end{aligned}$$
For a binary word \(w=c_1\cdots c_m\) let \(V_i(w)\) be its positional value in base \(B_i\), reduced modulo \(M\); every identity below is an equality of residues. In particular \(V_0(w)\) is the binary value of \(w\) modulo \(M\), which Horner's scheme produces in \(m\) steps, each replacing \(v\) by \(B_iv+c_j\). [7] Two rules then make \(V_0\) computable from the expression instead of the bits. Concatenation obeys
$$V_i(uv)=V_i(u)\,B_i^{|v|}+V_i(v),$$
and the substitution obeys
$$\begin{aligned}V_i(\mu(w))&=G_{i+1}(m)\\ &\quad+(B_i-1)\,V_{i+1}(w),\\ G_j(m)&=\sum_{r=0}^{m-1}B_j^{\,r}.\end{aligned}$$
The second rule carries the whole evaluation and is short to prove. Group the letters of \(\mu(w)\) in pairs: each pair is one group \(\mu(c_j)\), and a pair of bits read in base \(B_i\) behaves as one digit in base \(B_i^2\equiv B_{i+1}\). The group \(\mu(0)=01\) has value \(1\) and \(\mu(1)=10\) has value \(B_i\), that is \(1+(B_i-1)c_j\) in both cases; summing \(\bigl(1+(B_i-1)c_j\bigr)B_{i+1}^{\,m-j}\) over \(j\) separates into \(G_{i+1}(m)\) and \((B_i-1)V_{i+1}(w)\). The index rises by one with every substitution, so a word with \(D\) nested applications of \(\mu\) needs \(B_0,\ldots,B_D\); our deepest nesting is \(29\) and the programs keep \(66\) residues per node.
As a check, follow \(10010=\mathrm{EO}(100)=\mu(10)\,0\). Here \(G_1(2)=1+4=5\) and \(V_1(10)=4\), so \(V_0(\mu(10))=5+1\cdot4=9\), the value of \(1001\), and the concatenation rule gives \(V_0(10010)=9\cdot2+0=18\). Both \(B^m\) and \(G(m)\) come out of one halving recursion from \(B^0=1\), \(G(0)=0\): an even argument uses \(B^{2h}=(B^h)^2\) with \(G(2h)=G(h)\bigl(1+B^h\bigr)\), and an odd one is reduced by a step, \(B^{m}=B\,B^{m-1}\) with \(G(m)=G(m-1)+B^{m-1}\). Each pair therefore costs \(O(\log m)\) modular operations. [8]
10. The expression graph
The recursion of Section 6 never produces bits: it produces a small acyclic expression whose leaves are blocks of length at most \(3\) and whose inner nodes are concatenations and applications of \(\mu\). One detail remains. The maps need \(a_2\cdots a_n\), \(a_1\cdots a_n\) and \(a_2\cdots a_{n+1}\), and the ancestor is itself an expression, so dropping its first or last letter must be symbolic as well. For a concatenation this recurses into one side; for \(\mu(v)\) two rules suffice. Dropping the last letter of \(\mu(v)\) leaves \(\mu(v')\) followed by the last letter of \(v\), where \(v'\) is \(v\) without its last letter, because the final group starts with that letter; dropping the first letter leaves \(1-v_1\) followed by \(\mu(v'')\), with \(v''\) the tail of \(v\). This is why every node records its length and its first and last letter. For all eighteen targets the finished graph has at most \(132\) nodes.
How the Code Works
The three programs share one structure. Each forms the eighteen indices \(10^k\), computes one residue for each and adds them modulo \(10^9\) — in main for C++, in solve for Python and Java — and the C++ program runs run_validation first. Java spells the helpers below in camel case: pCount, sCount, findLengthForIndex, computeAMod. p_count and s_count implement the recurrences of Sections 4 and 7 with the hard-coded bases \(p(1),p(2),p(3)\) and \(s(1),s(2),s(3)\), memoised in hash maps; count_upto returns \(c(\ell)=s(\ell)/2\), an exact halving because every \(p(j)\) is even. find_length_for_index is the search of Section 8: it doubles hi while count_upto(hi) is below \(n\), then bisects. compute_A_mod glues the steps together — it returns \(0\) for \(n=0\), and otherwise computes \(\ell\), the rank \(k\) and the position \(j\) as offset, builds the expression and reports its residue.
The expression lives in a WordBuilder, an arena of EMPTY, LEAF, CONCAT and MU nodes addressed by index, which starts with shared nodes for the empty word and the two letters. make_bit returns one of those, make_leaf reuses them for words of length at most one, and make_concat returns the other operand when one is empty, which makes the shortest base cases harmless; every other node is appended with its len, first and last. prefix, suffix and middle are the letter-dropping operations of Section 10, each with its own cache, and middle of a concatenation is prefix of its suffix. The maps of Section 4 are odd_even_map, even_odd_map and odd_odd_map, and select_node is the recursion of Section 6: it takes size_O1 and size_E from p_count, compares k with them and descends, the outer branches sharing one map with ancestor ranks \(p(n+1)/2+k\) and \(k\) minus the two sizes.
Evaluation is the vector-valued eval, named get_eval in Python and getEval in Java. For a node it fills an array of \(66\) residues, cached on the node: a leaf by Horner in each base, a concatenation by the concatenation rule with pow_base, a MU node by the substitution rule, reading the child's entry at index \(i+1\) and adding sum_geom(i+1, …). A MU node can therefore fill only the indices \(0\) to \(64\) and leaves the last slot at zero; an entry at index \(i\) is correct as long as \(i\) plus the substitution depth below that node fits in the array, which the requested entry \(0\) does with a wide margin. Both helpers come from one memoised pow_sum returning the pair \(\bigl(B_i^{\,m},G_i(m)\bigr)\).
The implementations differ in their numeric types and in their built-in checks. C++ multiplies through mul_mod and forms the cumulative-count intermediates in __int128, keeping the counts themselves in unsigned long long; Java uses plain long and may multiply directly, both factors being below \(10^9\); Python has unbounded integers and needs no care, and its Node uses __slots__. Only C++ carries self-checks: run_validation asserts \(p(1),\ldots,p(4)\) and \(s(3),s(4)\), then rebuilds the blocks for \(N=100\) and \(N=1000\) with materialize, which expands an expression into bits while it stays under a limit, and asserts \(3251\) and \(80852364498\). All three pad the printed result to nine digits, so a leading zero cannot be lost.
Complexity and Verification
Write \(L\) for \(\log_2\ell\) and \(G\) for the number of expression nodes. The length search makes \(O(L)\) calls to count_upto, each visiting \(O(L)\) arguments by Section 7, so \(O(L^2)\) operations on fixed-width integers; for \(N=10^{18}\) it makes exactly \(63\) calls and leaves \(435\) entries for \(p\) and \(901\) for \(s\). Selection descends \(O(L)\) levels, and each node is evaluated once in \(66\) bases with a halving chain of \(O(L)\) exponents per base, which is \(O(66\,G\,L)\) operations; for the eighteen targets \(G\le132\), a few thousand operations per query. Storage is of the same order, the retained memo tables included, so it is \(O(L^2)\) words rather than \(O(L)\) — but nothing depends on the \(10^9\) bits of the block, which is the point of Sections 9 and 10. The largest quantity stored is \(s(2^{31})\approx7.3\cdot10^{18}\), the bound the doubling reaches, and the recurrence forms \(s(2^{31})+8\) before subtracting \(8\); both stay below \(2^{63}-1\approx9.22\cdot10^{18}\), so Java's signed long suffices, with about a fifth of its range to spare. Measured on one machine, all eighteen queries and the checks take about \(9\) ms in C++, \(34\) ms in Python and \(49\) ms in Java.
The claims above were verified independently of the solution files. The first \(2^{21}\) letters of \(t\) were generated from the recurrence, cross-checked against the parity of the binary digit sum, and all blocks of every length \(\ell\le64\) collected by brute force.
- The counts agree with
p_countfor every \(\ell\le64\). Since Section 4 derives \(p(\ell)\) independently of any prefix, a collection of that many blocks is complete, so this prefix holds every block of those lengths; a prefix of \(2^{17}\) letters gives the same sets. - Complement closure and the exact half-count of leading ones hold for every \(\ell\le40\), the prefix classes of Lemma 3 for every word, and the decomposition of Section 5 reproduces the sorted list of blocks part by part and in the ancestor order for every \(4\le\ell\le40\).
- The closed form agrees with the recurrence for \(3\le\ell\le2\cdot10^5\) and \(s(\ell)=\sum_{j\le\ell}p(j)\) with direct summation for \(\ell\le2\cdot10^5\); the substitution rule of Section 9 matches a direct base-\(B_i\) evaluation on \(400\) random words.
- The brute-force list gives \(A_n\) for every \(n\le3196\), including both checkpoints, and
compute_A_modagrees with all of them modulo \(10^9\); the C++ assertions pass when enabled, and the output of the three programs, compared byte for byte, is identical. [9]
Footnotes and References
The references below supply the background facts used above; the reading-frame maps, the three-part order, the counting recurrences and the substitution rule are derived in this article.
- Project Euler 361 — Subsequence of Thue-Morse Sequence. The official statement defines the sequence by T₀ = 0, T₂ₙ = Tₙ and T₂ₙ₊₁ = 1 − Tₙ, defines {Aₙ} as the sorted integers whose binary expression appears in it, with 18 = 10010 appearing and 14 = 1110 not, tabulates A₀ to A₁₂, gives A₁₀₀ = 3251 and A₁₀₀₀ = 80852364498, and asks for the last nine digits of the sum of the terms at the indices 10, 100, …, 10¹⁸.
- Wikipedia — Thue-Morse sequence. It identifies Tₙ as the parity of the number of ones in the binary expansion of n, states that the sequence is the fixed point of 0 → 01, 1 → 10, and records that it contains no cube XXX and no overlap 0X0X0 or 1X1X1.
- OEIS A005942 — subword complexity of the Thue-Morse sequence. It defines a(n) as the number of factors of length n, lists 1, 2, 4, 6, 10, 12, 16, 20, 22, 24, 28, 32, 36, … and states a(2n) = a(n) + a(n+1), a(2n+1) = 2a(n+1) for n ≥ 2, the recurrence of Section 4.
- OEIS A053646 — distance to the nearest power of 2. Its terms 0, 0, 1, 0, 1, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1, 0, … are the d(m) of the closed form in Section 4, which A005942 states for n > 2.
- J.-P. Allouche and J. Shallit, Automatic Sequences. Cambridge University Press, 2003; its Chapter 10, "Subwords" (pages 298 to 344), treats the factors of automatic sequences, and A005942 cites its Problem 10, page 335, for this complexity.
- Wikipedia — Complexity function. It defines the complexity function of an infinite word as the number of distinct factors of each length, the p(ℓ) of Section 2.
- Wikipedia — Horner's method. The nested scheme evaluates a degree-n polynomial with n multiplications and n additions, which is optimal, and it turns a digit string of a positional system into a number.
- Wikipedia — Exponentiation by squaring. The identities x²ʰ = (xʰ)² and x²ʰ⁺¹ = x·(xʰ)² compute a power with about log₂ n squarings, the recursion of Section 9.
- C++ — Euler361.cpp. The linked immutable revision contains
p_count,s_count,find_length_for_index, theWordBuilderarena,select_nodeand the assertions ofrun_validation; it is the implementation record for this article. The brute-force comparisons above are additional validation, not packaged in that file.
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());
}
}