Problem 384: Rudin-Shapiro Sequence

View on Project Euler

Project Euler Problem 384 Solution

EulerSolve provides an optimized solution for Project Euler Problem 384, Rudin-Shapiro Sequence, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The Rudin-Shapiro sign sequence used by the solutions is $$b(n)=(-1)^{\operatorname{popcount}(n \mathbin{\&} (n \gg 1))}\in\{-1,+1\}.$$ The bit count here is exactly the number of overlapping occurrences of the binary pattern \(11\) inside the expansion of \(n\). Its prefix walk is $$A(n)=\sum_{k=0}^{n} b(k).$$ For integers \(t\ge 1\) and \(1\le c\le t\), define \(g(t,c)\) as the \(c\)-th nonnegative index \(n\) such that \(A(n)=t\). The program evaluates $$\sum_{t=2}^{45} g(F_t,F_{t-1}),\qquad F_0=F_1=1.$$ A direct simulation would have to generate the walk up to extremely large indices. The local C++, Python, and Java implementations avoid that completely: they count how often level \(t\) has appeared up to a bound \(x\), then invert that monotone counting function by binary search. Mathematical Approach Step 1: Two-bit identities for \(b(n)\) Write $$q(n)=\operatorname{popcount}(n \mathbin{\&} (n \gg 1)) \bmod 2,$$ so that \(b(n)=(-1)^{q(n)}\). Looking at the last two binary digits gives four exact identities: $$q(4m)=q(m),\qquad q(4m+1)=q(m),$$ $$q(4m+2)\equiv q(m)+m \pmod 2,\qquad q(4m+3)\equiv q(m)+m+1 \pmod 2.$$ Therefore $$b(4m)=b(m),\qquad b(4m+1)=b(m),$$ $$b(4m+2)=(-1)^m b(m),\qquad b(4m+3)=-(-1)^m b(m).$$ This 4-block structure is the real mathematical source of the recursion in all three solution files....

Detailed mathematical approach

Problem Summary

For a non-negative integer \(n\), let \(a(n)\) be the number of adjacent pairs of ones in the binary expansion of \(n\), overlaps included: \(a(5)=a(101_2)=0\), \(a(6)=a(110_2)=1\) and \(a(7)=a(111_2)=2\), because \(111_2\) contains \(11\) twice. The Rudin-Shapiro sequence is the sign sequence \(b(n)=(-1)^{a(n)}\), and \(s\) is its summatory sequence

$$s(n)=\sum_{i=0}^{n}b(i).$$

For \(n=0,1,\dots,7\) the statement lists \(a(n)=0,0,0,1,0,0,1,2\), hence \(b(n)=1,1,1,-1,1,1,-1,1\) and \(s(n)=1,2,3,2,3,4,3,4\). Every value of \(s\) is positive, and every positive integer \(k\) occurs in \(s\) exactly \(k\) times, so for \(1\le c\le t\) there is a well-defined index \(g(t,c)\) at which \(t\) occurs in \(s\) for the \(c\)-th time. The statement gives \(g(3,3)=6\), \(g(4,2)=7\) and \(g(54321,12345)=1220847710\). With the Fibonacci numbers \(F(0)=F(1)=1\) and \(F(n)=F(n-1)+F(n-2)\) for \(n\gt1\), and with \(GF(t)=g(F(t),F(t-1))\), we must find the sum of \(GF(t)\) over \(2\le t\le45\). [1]

Since \(F(45)=1836311903\), the largest level involved is about \(1.8\cdot10^{9}\), and the bound \(s(n)\lt\sqrt{6n}\) of Section 11 puts every visit to that level beyond \(F(45)^2/6\), which is more than \(5\cdot10^{17}\). Walking along \(s\) one index at a time is out of the question.

The whole solution in one sentence: for a fixed level \(t\) we count the indices \(n\le x\) with \(s(n)=t\) by a recursion that strips the last two binary digits of \(n\), and then, because that count grows with \(x\), we recover \(g(t,c)\) from it by binary search.

Mathematical Approach

One observation drives everything. Appending two binary digits to \(n\) changes \(a(n)\) by an amount that depends only on the last digit of \(n\), so \(b\) and \(s\) obey exact rules linking index \(4m+r\) to index \(m\). Those rules turn “how often has level \(t\) been reached by index \(x\)?” into the same question for a level near \(t/2\) and a bound near \(x/4\).

1. Counting the pattern 11 with one bit operation

The pairs counted by \(a(n)\) are the positions \(j\) at which the binary digits \(j\) and \(j+1\) of \(n\) are both one, which is exactly where n & (n >> 1) carries a one-bit. So \(a(n)\) is the number of one-bits of that value: for \(n=7=111_2\) the operation turns \(111_2\) and \(11_2\) into \(011_2\), which has two one-bits, matching \(a(7)=2\). The three sequences \(a\), \(b\) and \(s\) are catalogued as A014081, A020985 and A020986. [3]

2. What the last two binary digits do to b

Write \(n=4m+r\) with \(r\in\{0,1,2,3\}\); in binary \(n\) is \(m\) followed by the two digits of \(r\). A new pair can appear only inside those two digits or at the junction with the last digit of \(m\):

  • \(r=0\): the digits appended are \(00\), so no pair appears and \(a(4m)=a(m)\).
  • \(r=1\): the digits are \(01\), again no pair, so \(a(4m+1)=a(m)\).
  • \(r=2\): the digits are \(10\), and the leading one pairs with the last digit of \(m\) when that digit is one, so \(a(4m+2)=a(m)+(m\bmod 2)\).
  • \(r=3\): the digits are \(11\), which is one pair by itself, plus one more at the junction when \(m\) is odd, so \(a(4m+3)=a(m)+(m\bmod 2)+1\).

Translated into signs, these read

$$b(4m)=b(4m+1)=b(m),$$

$$b(4m+2)=(-1)^m b(m),\qquad b(4m+3)=-(-1)^m b(m).$$

These identities have a sequential counterpart. Reading the binary digits of \(n\) from the left, only two bits ever need to be remembered: the previous digit, and the parity of the pairs seen so far. The parity flips exactly when the previous digit and the new one are both one. Those four states, with output \(+1\) for even parity and \(-1\) for odd, are why the Rudin-Shapiro sequence is called 2-automatic. [2] [4]

3. The same two digits turn s into a doubling rule

Adding the four identities over \(r=0,1,2,3\) makes the \(r=2\) and \(r=3\) terms cancel and leaves

$$\sum_{r=0}^{3}b(4m+r)=2b(m).$$

Summing this over \(m=0,1,\dots,M\) covers precisely the indices \(0\) to \(4M+3\), so \(s(4M+3)=2s(M)\). Peeling the block off one term at a time yields all four rules:

$$s(4m)=2s(m)-b(m),\qquad s(4m+1)=2s(m),$$

$$s(4m+2)=2s(m)+(-1)^m b(m),\qquad s(4m+3)=2s(m).$$

For \(m=3\), where \(s(3)=2\) and \(b(3)=-1\), the rules predict \(s(12),\dots,s(15)=5,4,5,4\), matching direct summation at the indices \(12\) to \(15\).

4. s(n) is never zero or negative

We argue by induction on \(n\). First, \(s(0)=1\). For \(n\ge1\) put \(m=\lfloor n/4\rfloor\), so that \(m\lt n\) and \(s(m)\ge1\) holds by induction. Since \(b(m)=\pm1\), each of the four rules gives \(s(n)\ge2s(m)-1\ge1\). The same inequality, as \(s(n)-1\ge2\bigl(s(m)-1\bigr)\), is what Section 11 iterates.

5. A parity rule that fixes which residues can carry a level

\(s(n)\) is a sum of \(n+1\) terms, each \(\pm1\), so

$$s(n)\equiv n+1\pmod 2.$$

Together with Section 3 this splits the four residues by parity: \(s(4m+1)\) and \(s(4m+3)\) equal \(2s(m)\) and are even, while \(s(4m)\) and \(s(4m+2)\) are even plus \(\pm1\) and are odd. So an even level sits only at indices \(\equiv1,3\pmod4\), an odd level only at \(\equiv0,2\pmod4\).

The parity rule also recovers the parity of a hidden index from its level: \(s(m)=u\) gives \(m\equiv u+1\pmod2\), hence \((-1)^m=(-1)^{u+1}\). That is what makes the rules of Section 3 usable when only the level of \(m\) is known.

6. The counting function to be inverted

Fix a level \(t\ge1\) and a bound \(x\). Because the rules for \(s(4m)\) and \(s(4m+2)\) contain \(b(m)\), counting indices by level alone would throw away exactly the information the next step needs, so we count with the sign attached:

$$P_t(x)=\#\{0\le n\le x: s(n)=t,\ b(n)=+1\},$$

$$N_t(x)=\#\{0\le n\le x: s(n)=t,\ b(n)=-1\}.$$

Let \(C_t(x)=P_t(x)+N_t(x)\) count all visits to level \(t\) up to \(x\). As a function of \(x\) it is non-decreasing and grows by \(1\) exactly at the indices where level \(t\) is reached, so inverting it gives what we want:

$$g(t,c)=\min\{x\ge0: C_t(x)\ge c\}.$$

Two boundary values start the recursion. For \(x\lt0\) the range is empty, so \(P_t(x)=N_t(x)=0\). And \(P_1(x)=1\), \(N_1(x)=0\) for \(x\ge0\), because level \(1\) is reached only at \(n=0\), where \(b(0)=+1\): an odd level sits at an index \(4m\) or \(4m+2\) by Section 5, and \(s(n)=1\) with \(s(m)\ge1\) forces \(s(m)=1\) in both rules of Section 3, so a descent from any \(n\ge1\) reaches \(m=0\), and of \(n=0,1,2,3\) only \(n=0\) has \(s(n)=1\).

7. The recursion for an even level

Let \(t=2u\) with \(u\ge1\). By Section 5 only the residues \(1\) and \(3\) can contribute, and both carry \(s(4m+r)=2s(m)\), so the index \(4m+r\le x\) is counted exactly when \(s(m)=u\) and \(m\le x_r\), where

$$x_r=\left\lfloor\frac{x-r}{4}\right\rfloor,\qquad r\in\{0,1,2,3\}.$$

The signs decide where each count lands. Residue \(1\) passes the sign through, since \(b(4m+1)=b(m)\). For residue \(3\) we have \(b(4m+3)=-(-1)^m b(m)\), and \(s(m)=u\) forces \((-1)^m=(-1)^{u+1}\) by Section 5, hence \(b(4m+3)=(-1)^u b(m)\): the sign survives when \(u\) is even and flips when \(u\) is odd. For even \(u\) this gives

$$\begin{aligned}P_{2u}(x)&=P_u(x_1)+P_u(x_3),\\ N_{2u}(x)&=N_u(x_1)+N_u(x_3),\end{aligned}$$

and for odd \(u\),

$$\begin{aligned}P_{2u}(x)&=P_u(x_1)+N_u(x_3),\\ N_{2u}(x)&=N_u(x_1)+P_u(x_3).\end{aligned}$$

8. The recursion for an odd level

Let \(t=2u+1\) with \(u\ge1\). Only the residues \(0\) and \(2\) can contribute, and now the equations involve \(b(m)\) as well as \(s(m)\).

At \(n=4m\) the equation \(2s(m)-b(m)=2u+1\) has exactly two solutions: \(s(m)=u\) with \(b(m)=-1\), or \(s(m)=u+1\) with \(b(m)=+1\). Since \(b(4m)=b(m)\), the first case lands in \(N_{2u+1}\) and is counted by \(N_u(x_0)\), the second lands in \(P_{2u+1}\) and is counted by \(P_{u+1}(x_0)\).

At \(n=4m+2\) the equation \(2s(m)+(-1)^m b(m)=2u+1\) splits in the same way: \(s(m)=u\) with \((-1)^m b(m)=+1\), or \(s(m)=u+1\) with \((-1)^m b(m)=-1\). Here \(b(4m+2)=(-1)^m b(m)\), so the first case lands in \(P_{2u+1}\) and the second in \(N_{2u+1}\). The parity rule turns both conditions into conditions on \(b(m)\) alone: \(s(m)=u\) gives \((-1)^m=(-1)^{u+1}\), so \((-1)^m b(m)=+1\) means \(b(m)=(-1)^{u+1}\), and \(s(m)=u+1\) gives \((-1)^m=(-1)^{u}\), so \((-1)^m b(m)=-1\) means \(b(m)=-(-1)^{u}\). For odd \(u\) both read \(b(m)=+1\); for even \(u\) both read \(b(m)=-1\). Collecting the four pieces, for odd \(u\),

$$\begin{aligned}P_{2u+1}(x)&=P_u(x_2)+P_{u+1}(x_0),\\ N_{2u+1}(x)&=N_u(x_0)+P_{u+1}(x_2),\end{aligned}$$

and for even \(u\),

$$\begin{aligned}P_{2u+1}(x)&=N_u(x_2)+P_{u+1}(x_0),\\ N_{2u+1}(x)&=N_u(x_0)+N_{u+1}(x_2).\end{aligned}$$

9. A worked example: the third visit to level 5

Level \(5\) is odd with \(u=2\), which is even, so \(P_5(x)=N_2(x_2)+P_3(x_0)\) and \(N_5(x)=N_2(x_0)+N_3(x_2)\). Take \(x=14\): then \(x_0=\lfloor14/4\rfloor=3\) and \(x_2=\lfloor12/4\rfloor=3\).

Level \(2\) is even with \(u=1\), which is odd, so \(P_2(x)=P_1(x_1)+N_1(x_3)\) and \(N_2(x)=N_1(x_1)+P_1(x_3)\). At \(x=3\) both \(x_1\) and \(x_3\) are \(0\), and the boundary values \(P_1(0)=1\), \(N_1(0)=0\) give \(P_2(3)=1\) and \(N_2(3)=1\): the visits at \(n=1\) and \(n=3\), where \(b(1)=+1\) and \(b(3)=-1\).

Level \(3\) is odd with \(u=1\), which is odd, so \(P_3(x)=P_1(x_2)+P_2(x_0)\) and \(N_3(x)=N_1(x_0)+P_2(x_2)\). At \(x=3\) we get \(x_0=x_2=0\), and \(P_2(0)=0\) because both of its own arguments are \(-1\), so \(P_3(3)=1\) and \(N_3(3)=0\), the single visit at \(n=2\).

Putting the pieces together, \(P_5(14)=1+1=2\) and \(N_5(14)=1+0=1\), so \(C_5(14)=3\). At \(x=13\) the arguments become \(x_0=3\) and \(x_2=2\), where \(N_2(2)=0\) and \(N_3(2)=0\), so \(C_5(13)=2\). The count first reaches \(3\) at \(x=14\), so \(g(5,3)=14\). A direct walk confirms it: level \(5\) is visited at \(n=8,12,14,16\) and \(26\). Since \(F(4)=5\) and \(F(3)=3\), this is \(GF(4)=14\); the two smaller terms of the required sum are \(GF(2)=g(2,1)=1\) and \(GF(3)=g(3,2)=4\).

10. How often a level is reached, and with which sign

Let \(p(t)\) and \(q(t)\) be the totals of \(P_t(x)\) and \(N_t(x)\) over the whole sequence, that is, for arbitrarily large \(x\). Every argument \(x_r\) then grows beyond all bounds too, so Sections 7 and 8 become recursions in the level alone:

  • \(t=2u\): for even \(u\), \(p(2u)=2p(u)\) and \(q(2u)=2q(u)\); for odd \(u\), \(p(2u)=q(2u)=p(u)+q(u)\).
  • \(t=2u+1\): for odd \(u\), \(p(2u+1)=p(u)+p(u+1)\) and \(q(2u+1)=q(u)+p(u+1)\); for even \(u\), \(p(2u+1)=q(u)+p(u+1)\) and \(q(2u+1)=q(u)+q(u+1)\).

We claim \(p(t)=\lceil t/2\rceil\) and \(q(t)=\lfloor t/2\rfloor\) for every \(t\ge1\), by induction on \(t\). For \(t=1\) the boundary values give \(p(1)=1\) and \(q(1)=0\), as claimed. For even \(t=2u\) the recurrence uses only \(u\lt t\), and for odd \(t=2u+1\ge3\) both \(u\) and \(u+1\) are smaller than \(t\); so the totals on the right are finite by induction, which makes those at \(t\) finite too, and the claimed formulas may be substituted. If \(t=2u\), either \(u\) is even and \(2p(u)=2q(u)=u\), or \(u\) is odd and \(p(u)+q(u)=u\); either way \(p(t)=q(t)=u=\lceil t/2\rceil=\lfloor t/2\rfloor\). If \(t=2u+1\) with \(u\) odd, then \(p(t)=\frac{u+1}2+\frac{u+1}2=u+1\) and \(q(t)=\frac{u-1}2+\frac{u+1}2=u\). If \(t=2u+1\) with \(u\) even, then \(p(t)=\frac u2+\frac{u+2}2=u+1\) and \(q(t)=\frac u2+\frac u2=u\). In both odd cases \(u+1=\lceil t/2\rceil\) and \(u=\lfloor t/2\rfloor\), which closes the induction.

Adding the two formulas gives \(p(t)+q(t)=t\). So every positive integer \(k\) really does occur exactly \(k\) times in \(s\), as the statement asserts, and \(g(t,c)\) is defined precisely for \(1\le c\le t\). The sign split comes for free: of the \(t\) visits to level \(t\), exactly \(\lceil t/2\rceil\) have \(b(n)=+1\).

11. How far the search must go

For \(n\ge1\), Section 4 gives \(s(n)-1\ge2\bigl(s(\lfloor n/4\rfloor)-1\bigr)\). Iterating it \(k\) times, with \(k=\lfloor\log_4 n\rfloor\) so that \(\lfloor n/4^k\rfloor\) lands in \(\{1,2,3\}\) where \(s-1\ge1\), yields

$$s(n)\ge1+2^{\lfloor\log_4 n\rfloor}\gt1+\tfrac12\sqrt n.$$

Hence \(s(n)=t\) implies \(n\lt4t^2\): the indices are quadratic in the level, so the doubling phase needs about \(2\log_2 t\) steps. The quadratic order is right: level \(k\) occurs \(k\) times, and after the \(t\)-th visit to level \(t\) the walk can never return to \(t\) or below, since \(s\) moves by \(\pm1\) and is unbounded above; so the \(t(t+1)/2\) indices with \(s(n)\le t\) are all at most \(g(t,t)\), giving \(g(t,t)\ge t(t+1)/2-1\).

Sharp constants are known: \(\sqrt{3n/5}\lt s(n)\lt\sqrt{6n}\) for \(n\ge1\). [2] [5] The lower bound turns \(s(n)=t\) into \(n\lt5t^2/3\), so \(GF(t)\lt5F(t)^2/3\). With the Fibonacci identity \(\sum_{t=0}^{n}F(t)^2=F(n)F(n+1)\), which is the classical sum-of-squares identity shifted by one index, [8] the whole sum obeys

$$\sum_{t=2}^{45}GF(t)\lt\tfrac53\bigl(F(45)F(46)-2\bigr),$$

and the right-hand side is \(9093462674871523195\), below \(2^{63}-1\). Signed 64-bit arithmetic is therefore provably sufficient, with no appeal to the computed values.

12. Why only a handful of states appear

One evaluation of \(C_t(x)\) calls itself on the levels \(\lfloor t/2\rfloor\) and \(\lfloor t/2\rfloor+1\) and on the bounds \(x_0,\dots,x_3\). Both families collapse at once.

Levels: suppose the levels at some depth all lie in \(\{v,v+1\}\). If \(v=2w\), its children are \(\{w\}\) and those of \(v+1\) are \(\{w,w+1\}\); if \(v=2w+1\), its children are \(\{w,w+1\}\) and those of \(v+1\) are \(\{w+1\}\). Either way the next depth lies in \(\{w,w+1\}\) with \(w=\lfloor v/2\rfloor\): at most two consecutive levels per depth.

Bounds: the four values \(x_0,\dots,x_3\) lie between \(\lfloor(x-3)/4\rfloor\) and \(\lfloor x/4\rfloor\), which differ by at most \(1\). If a depth holds \(\{y,y+1\}\), its children lie between \(\lfloor(y-3)/4\rfloor\) and \(\lfloor(y+1)/4\rfloor\), which differ by exactly \(1\): again at most two consecutive bounds per depth.

Each depth therefore holds at most \(2\cdot2=4\) distinct states \((t,x)\). Every child level is at most \(\lceil t/2\rceil\), so level \(1\) is reached within \(\lceil\log_2 t\rceil\) steps and the depth is governed by the level; the bound gives no comparable guarantee, because \(x_0=0\) when \(x=0\), so a bound of \(0\) can persist through further odd levels. At \(t=F(45)\) and \(x=2^{61}\), the largest bound the search of Section 6 ever tests, the measured depth is \(31=\lceil\log_2 F(45)\rceil\) with \(123\) distinct argument pairs, \(119\) of them cached, inside the bound \(4\cdot32=128\).

How the Code Works

All three programs implement the same two routines. prefix_counts(t, x) in C++ and Python, prefixCounts in Java, returns the pair \(\bigl(P_t(x),N_t(x)\bigr)\). It answers \((0,0)\) when \(t=0\) or \(x\lt0\) and \((1,0)\) when \(t=1\), the boundary values of Section 6; otherwise the parity of \(t\) selects the branch of Section 7 or of Section 8, and inside it the parity of \(u=\lfloor t/2\rfloor\) picks one of the two sign patterns. The children are evaluated, their components added, and the result cached under \((t,x)\). That memoisation is what keeps Section 12's four states per depth from unfolding into a call tree of about \(8\cdot10^{14}\) nodes. [6]

The arguments \(x_r\) must be rounded towards minus infinity, not truncated towards zero: near the leaves \(x-r\) is often negative, and truncation would turn \(-1\) into \(0\) and count an index that does not exist. Python's // floors by definition, Java calls Math.floorDiv, and the C++ file supplies floor_div, because C++ division truncates.

count_total, called countTotal in Java, adds the two components to give \(C_t(x)\), and g(t, c) performs the inversion of Section 6 as an exponential search followed by a binary search. [7] It starts at hi = 1 and doubles while count_total(t, hi) is below \(c\), then halves the interval between lo and hi, keeping hi a bound whose count reaches \(c\), and returns the least \(x\) with \(C_t(x)\ge c\). All three guard the entry with a test for \(c\gt t\) and return \(0\); by Section 10 that never fires here, since \(1\le F(t-1)\le F(t)\). In Python and Java, solve clears the cache, builds the Fibonacci table with \(F(0)=F(1)=1\) up to \(F(45)\) and adds g(fib[t], fib[t-1]) for \(t=2,\dots,45\). C++ has no single function of that name: fibonacci_values returns the table, and the summing routine reached from main accumulates \(g\) over \(t\) from kMinT to kMaxT.

The rest is plumbing. C++ wraps the cache and the three routines in a class RudinShapiroSelector whose memo_ is an unordered_map keyed by a CacheKey and hashed by CacheKeyHash; Java has the matching CacheKey.hashCode and a HashMap; Python uses the tuple (t, x) as a dictionary key. C++ holds counts and the running sum in uint64_t and search bounds in int64_t, Java uses signed long throughout, and Python's integers are unbounded. Only C++ guards its doubling loop, capping the bound at the largest int64_t; Java has no such guard. Both searches are nonetheless safe at the levels required here: by Section 11 an index \(n\ge2^{62}\) forces \(s(n)\ge1+2^{31}\), which exceeds \(F(45)\), so every visit to every required level lies below \(2^{62}\). The recursion depth stays near \(31\), far below Python's default limit.

The C++ file also carries the self-checks. sign_for_b evaluates \(b(n)\) with the bit trick of Section 1 and __builtin_popcountll, a GCC and Clang builtin. run_small_bruteforce_check walks the real sequence index by index, recording the first \(t\) indices at which each level \(t\le18\) appears, and compares all \(171\) of those values against the recursion. run_checkpoints verifies the three checkpoints of the statement, runs that walk, and compares two evaluations of the full sum. The comparison applies the same recursion twice, so it checks consistency; the direct walk is what tests the recursion against the definition. --skip-checkpoints turns all of this off. [9]

Complexity and Verification

By Section 12 one evaluation of \(C_t(x)\) visits \(O(\log t)\) distinct states, with constant work and one cache entry each, so it costs \(O(\log t)\) time. By Section 11 the search tests bounds of order \(t^2\) and therefore makes \(O(\log t)\) evaluations, so one value of \(GF\) costs \(O(\log^2 t)\) time and retains \(O(\log^2 t)\) cache entries, and there are \(44\) values. Instrumenting those \(44\) selections with one cache that starts empty gives \(2801\) calls of \(C_t(x)\) and \(83082\) cached states, each holding \(t\), \(x\), \(P_t(x)\) and \(N_t(x)\) plus hash-table and object overhead, and excluding the self-checks. On the machine used here the C++ process stayed below \(10\) MB and all three programs returned in well under a second, self-checks included.

These facts were checked by running code. The identities of Sections 2 and 3 and the residue split of Section 5 hold for every \(m\lt10^6\), and the parity rule \(s(n)\equiv n+1\pmod2\) for every \(n\lt4\cdot10^6\). The recursion of Sections 7 and 8 reproduces \(P_t(x)\) and \(N_t(x)\) from a direct walk for all levels \(t\le40\) and bounds \(x\le3000\), and \(g(t,c)\) agrees with it for all \(1\le c\le t\le400\), where the largest index needed is \(125695\). Evaluated at the bound \(10^{40}\), the recursion returns \(\lceil t/2\rceil\) and \(\lfloor t/2\rfloor\) for every \(t\lt3000\), confirming Section 10. Over \(1\le n\le3\cdot10^6\) the extremes of \(s(n)/\sqrt n\) are \(0.775353\) and \(2.448892\), tightly inside the Section 11 bounds \(\sqrt{3/5}=0.774597\) and \(\sqrt6=2.449490\).

The three checkpoints of the statement are reproduced by the recursion and by the C++ self-check. The C++, Python and Java programs print the same final value, with and without the C++ checkpoints. One caution: outside \(1\le c\le t\) all three return \(0\), a legitimate index, so a caller straying outside that range gets a silently wrong answer rather than an error.

Footnotes and References

The references supply the background facts used above. The identities for \(b\) and \(s\), the signed counting recursion, the closed forms for \(p(t)\) and \(q(t)\), the elementary bound \(s(n)\gt1+\sqrt n/2\) for \(n\ge1\) and the state count of Section 12 are derived in this article.

  1. Project Euler 384 — Rudin-Shapiro Sequence. The official statement defines a(n) as the number of adjacent, possibly overlapping, pairs of ones in the binary expansion of n, sets b(n) = (−1)^a(n) and s(n) as its partial sums, tabulates all three for n ≤ 7, states that every positive integer k occurs exactly k times in s, gives the three checkpoints for g, and asks for the sum of GF(t) = g(F(t), F(t−1)) over 2 ≤ t ≤ 45.
  2. Wikipedia — Rudin–Shapiro sequence. The article counts the occurrences of the block 11 in the binary expansion of n with overlaps, two of them for 7 = 111₂, defines the sign sequence as (−1) to that power, states that the partial sums taken from index 0 satisfy √(3n/5) < s(n) < √(6n) for n ≥ 1, and records that it is 2-automatic, generated by a four-state automaton.
  3. OEIS A014081 is the number of occurrences of 11 in the binary expansion of n, with a(15) = 3 for 1111₂, which shows that overlapping pairs are counted; OEIS A020985 is the Rudin-Shapiro sequence, equal to (−1) raised to A014081(n), with a(2n) = a(n) and a(2n+1) = (−1)ⁿ a(n); OEIS A020986 is its sequence of partial sums, beginning 1, 2, 3, 2, 3, 4, 3, 4, 5, 6, 7, 6.
  4. Wikipedia — Automatic sequence. A sequence is k-automatic when its n-th term is the output attached to the state a deterministic finite automaton reaches after reading the base-k digits of n; the article lists the Rudin–Shapiro sequence as a 2-automatic example.
  5. N. Rampersad and J. Shallit, Rudin-Shapiro Sums Via Automata Theory and Logic. The paper recovers many classical results of Brillhart and Morton on Rudin-Shapiro sums in one framework of automata theory and logic; the two-sided bound used in Section 11 is one of them.
  6. Wikipedia — Memoization. Storing the results of calls to a pure function and reusing them for repeated inputs is described there as a space–time tradeoff, usually implemented with a hash table, exactly how the cache keyed by (t, x) works.
  7. Wikipedia — Exponential search. The method doubles an index until it passes the target, then binary-searches the bracketed range, at a cost of O(log i) where i is the position of the target: the two-phase search inside g.
  8. Wikipedia — Fibonacci sequence. The article proves that the sum of the squares of the first n Fibonacci numbers equals the product of the n-th and the (n+1)-st; shifted by one index, because the statement starts with F(0) = F(1) = 1, this is the form used in Section 11.
  9. C++ — Euler384.cpp. The linked immutable C++ revision contains RudinShapiroSelector with prefix_counts and g, the helpers floor_div and sign_for_b, and the checks in run_small_bruteforce_check and run_checkpoints. It is the implementation record for this article.

Problem 384 source code

C++

#include <algorithm>
#include <atomic>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <unordered_map>
#include <utility>
#include <vector>

namespace {

using u64 = std::uint64_t;
using i64 = std::int64_t;

constexpr int kMinT = 2;
constexpr int kMaxT = 45;
constexpr u64 kCheckpointT = 54321ULL;
constexpr u64 kCheckpointC = 12345ULL;
constexpr u64 kCheckpointExpected = 1220847710ULL;

struct Options {
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

struct CountPair {
    u64 positive = 0ULL;
    u64 negative = 0ULL;
};

struct CacheKey {
    u64 t = 0ULL;
    i64 x = 0LL;

    bool operator==(const CacheKey& other) const {
        return t == other.t && x == other.x;
    }
};

struct CacheKeyHash {
    std::size_t operator()(const CacheKey& key) const {
        const std::size_t h1 = std::hash<u64>{}(key.t);
        const std::size_t h2 = std::hash<i64>{}(key.x);
        return h1 ^ (h2 + 0x9e3779b97f4a7c15ULL + (h1 << 6) + (h1 >> 2));
    }
};

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }

        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }

        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }

    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        return false;
    }

    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(const int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return true;
}

i64 floor_div(const i64 a, const i64 b) {
    if (a >= 0) {
        return a / b;
    }
    return -(((-a) + b - 1) / b);
}

int sign_for_b(const u64 n) {
    const unsigned parity = static_cast<unsigned>(__builtin_popcountll(n & (n >> 1ULL)) & 1U);
    return parity == 0U ? 1 : -1;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const std::size_t workload_units) {
    if (!allow_multithreading || workload_units < 2ULL) {
        return 1U;
    }

    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) {
            threads = 1U;
        }
    }

    return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload_units)));
}

class RudinShapiroSelector {
public:
    u64 g(const u64 t, const u64 c) {
        if (t == 0ULL || c == 0ULL || c > t) {
            return 0ULL;
        }

        i64 lo = 0;
        i64 hi = 1;

        while (count_total(t, hi) < c) {
            if (hi > std::numeric_limits<i64>::max() / 2LL) {
                hi = std::numeric_limits<i64>::max();
                break;
            }
            hi *= 2LL;
        }

        while (lo < hi) {
            const i64 mid = lo + (hi - lo) / 2LL;
            if (count_total(t, mid) >= c) {
                hi = mid;
            } else {
                lo = mid + 1LL;
            }
        }

        return static_cast<u64>(lo);
    }

private:
    CountPair prefix_counts(const u64 t, const i64 x) {
        if (t == 0ULL || x < 0LL) {
            return {0ULL, 0ULL};
        }

        if (t == 1ULL) {
            return {1ULL, 0ULL};
        }

        const CacheKey key{t, x};
        const auto it = memo_.find(key);
        if (it != memo_.end()) {
            return it->second;
        }

        CountPair result;
        if ((t & 1ULL) == 0ULL) {
            const u64 u = t >> 1U;
            const i64 x1 = floor_div(x - 1LL, 4LL);
            const i64 x3 = floor_div(x - 3LL, 4LL);
            const CountPair a = prefix_counts(u, x1);
            const CountPair b = prefix_counts(u, x3);

            if ((u & 1ULL) == 0ULL) {
                result.positive = a.positive + b.positive;
                result.negative = a.negative + b.negative;
            } else {
                result.positive = a.positive + b.negative;
                result.negative = a.negative + b.positive;
            }
        } else {
            const u64 u = t >> 1U;
            const i64 x0 = floor_div(x, 4LL);
            const i64 x2 = floor_div(x - 2LL, 4LL);

            const CountPair u0 = prefix_counts(u, x0);
            const CountPair u2 = prefix_counts(u, x2);
            const CountPair v0 = prefix_counts(u + 1ULL, x0);
            const CountPair v2 = prefix_counts(u + 1ULL, x2);

            if ((u & 1ULL) == 1ULL) {
                // u odd.
                result.positive = u2.positive + v0.positive;
                result.negative = u0.negative + v2.positive;
            } else {
                // u even.
                result.positive = u2.negative + v0.positive;
                result.negative = u0.negative + v2.negative;
            }
        }

        memo_.emplace(key, result);
        return result;
    }

    u64 count_total(const u64 t, const i64 x) {
        const CountPair pair = prefix_counts(t, x);
        return pair.positive + pair.negative;
    }

    std::unordered_map<CacheKey, CountPair, CacheKeyHash> memo_;
};

std::vector<u64> fibonacci_values() {
    std::vector<u64> fib(static_cast<std::size_t>(kMaxT) + 1ULL, 0ULL);
    fib[0] = 1ULL;
    fib[1] = 1ULL;
    for (int i = 2; i <= kMaxT; ++i) {
        fib[static_cast<std::size_t>(i)] = fib[static_cast<std::size_t>(i - 1)] +
                                           fib[static_cast<std::size_t>(i - 2)];
    }
    return fib;
}

u64 solve_sum_gf_single_thread(const std::vector<u64>& fib) {
    RudinShapiroSelector selector;
    u64 sum = 0ULL;

    for (int t = kMinT; t <= kMaxT; ++t) {
        const u64 value = selector.g(fib[static_cast<std::size_t>(t)],
                                     fib[static_cast<std::size_t>(t - 1)]);
        sum += value;
    }

    return sum;
}

u64 solve_sum_gf_parallel(const std::vector<u64>& fib,
                          const bool allow_multithreading,
                          const unsigned requested_threads) {
    constexpr std::size_t kTaskCount = static_cast<std::size_t>(kMaxT - kMinT + 1);
    const unsigned thread_count = choose_thread_count(allow_multithreading,
                                                      requested_threads,
                                                      kTaskCount);
    if (thread_count == 1U) {
        return solve_sum_gf_single_thread(fib);
    }

    std::vector<u64> values(static_cast<std::size_t>(kMaxT) + 1ULL, 0ULL);
    std::atomic<int> next_t(kMinT);

    std::vector<std::thread> workers;
    workers.reserve(thread_count);
    for (unsigned tid = 0U; tid < thread_count; ++tid) {
        workers.emplace_back([&]() {
            RudinShapiroSelector selector;
            while (true) {
                const int t = next_t.fetch_add(1);
                if (t > kMaxT) {
                    break;
                }

                const u64 value = selector.g(fib[static_cast<std::size_t>(t)],
                                             fib[static_cast<std::size_t>(t - 1)]);
                values[static_cast<std::size_t>(t)] = value;
            }
        });
    }

    for (std::thread& worker : workers) {
        worker.join();
    }

    u64 sum = 0ULL;
    for (int t = kMinT; t <= kMaxT; ++t) {
        sum += values[static_cast<std::size_t>(t)];
    }
    return sum;
}

bool run_small_bruteforce_check() {
    constexpr int kMaxLevel = 18;
    std::vector<std::vector<u64>> positions(static_cast<std::size_t>(kMaxLevel) + 1ULL);
    std::vector<int> done(static_cast<std::size_t>(kMaxLevel) + 1ULL, 0);

    i64 current = 0;
    for (u64 n = 0ULL; n < 1'000'000ULL; ++n) {
        current += sign_for_b(n);
        if (current >= 1 && current <= kMaxLevel) {
            const std::size_t level = static_cast<std::size_t>(current);
            if (done[level] < static_cast<int>(level)) {
                positions[level].push_back(n);
                ++done[level];
            }
        }

        bool all_done = true;
        for (int level = 1; level <= kMaxLevel; ++level) {
            if (done[static_cast<std::size_t>(level)] < level) {
                all_done = false;
                break;
            }
        }
        if (all_done) {
            break;
        }
    }

    RudinShapiroSelector selector;
    for (int t = 1; t <= kMaxLevel; ++t) {
        for (int c = 1; c <= t; ++c) {
            const u64 got = selector.g(static_cast<u64>(t), static_cast<u64>(c));
            const u64 expected = positions[static_cast<std::size_t>(t)][static_cast<std::size_t>(c - 1)];
            if (got != expected) {
                std::cerr << "Small brute-force check failed at g(" << t << "," << c
                          << "): expected " << expected << ", got " << got << '\n';
                return false;
            }
        }
    }

    return true;
}

bool run_checkpoints(const Options& options) {
    RudinShapiroSelector selector;

    if (selector.g(3ULL, 3ULL) != 6ULL) {
        std::cerr << "Checkpoint failed: g(3,3) should be 6.\n";
        return false;
    }

    if (selector.g(4ULL, 2ULL) != 7ULL) {
        std::cerr << "Checkpoint failed: g(4,2) should be 7.\n";
        return false;
    }

    const u64 checkpoint = selector.g(kCheckpointT, kCheckpointC);
    if (checkpoint != kCheckpointExpected) {
        std::cerr << "Checkpoint failed: g(" << kCheckpointT << ',' << kCheckpointC << ") expected "
                  << kCheckpointExpected << ", got " << checkpoint << '\n';
        return false;
    }

    if (!run_small_bruteforce_check()) {
        return false;
    }

    const std::vector<u64> fib = fibonacci_values();
    const u64 single = solve_sum_gf_single_thread(fib);
    const u64 parallel = solve_sum_gf_parallel(fib,
                                               options.allow_multithreading,
                                               options.requested_threads);
    if (single != parallel) {
        std::cerr << "Parallel consistency check failed: single=" << single
                  << ", parallel=" << parallel << '\n';
        return false;
    }

    std::cout << "Checkpoints passed (threads="
              << choose_thread_count(options.allow_multithreading,
                                     options.requested_threads,
                                     static_cast<std::size_t>(kMaxT - kMinT + 1))
              << ").\n";

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    if (options.run_checkpoints && !run_checkpoints(options)) {
        return 1;
    }

    const std::vector<u64> fib = fibonacci_values();
    const u64 answer = solve_sum_gf_parallel(fib,
                                             options.allow_multithreading,
                                             options.requested_threads);

    std::cout << "Answer: " << answer << '\n';
    return 0;
}

Python

memo = {}

def prefix_counts(t, x):
    if t == 0 or x < 0: return (0, 0)
    if t == 1: return (1, 0)
    key = (t, x)
    if key in memo: return memo[key]
    
    if (t % 2) == 0:
        u = t // 2
        x1 = (x - 1) // 4
        x3 = (x - 3) // 4
        a_pos, a_neg = prefix_counts(u, x1)
        b_pos, b_neg = prefix_counts(u, x3)
        if (u % 2) == 0:
            res = (a_pos + b_pos, a_neg + b_neg)
        else:
            res = (a_pos + b_neg, a_neg + b_pos)
    else:
        u = t // 2
        x0 = x // 4
        x2 = (x - 2) // 4
        u0_p, u0_n = prefix_counts(u, x0)
        u2_p, u2_n = prefix_counts(u, x2)
        v0_p, v0_n = prefix_counts(u + 1, x0)
        v2_p, v2_n = prefix_counts(u + 1, x2)
        
        if (u % 2) == 1:
            res = (u2_p + v0_p, u0_n + v2_p)
        else:
            res = (u2_n + v0_p, u0_n + v2_n)
            
    memo[key] = res
    return res

def count_total(t, x):
    pos, neg = prefix_counts(t, x)
    return pos + neg

def g(t, c):
    if t == 0 or c == 0 or c > t: return 0
    lo = 0
    hi = 1
    while count_total(t, hi) < c:
        hi *= 2
    while lo < hi:
        mid = (lo + hi) // 2
        if count_total(t, mid) >= c:
            hi = mid
        else:
            lo = mid + 1
    return lo

def solve():
    memo.clear()
    fib = [0] * 46
    fib[0] = 1
    fib[1] = 1
    for i in range(2, 46):
        fib[i] = fib[i-1] + fib[i-2]
        
    total_val = 0
    for t in range(2, 46):
        total_val += g(fib[t], fib[t-1])
        
    return str(total_val)

if __name__ == '__main__':
    print(solve())

Java

import java.util.HashMap;
import java.util.Map;

public class Euler384 {
    static class CacheKey {
        long t, x;

        CacheKey(long t, long x) {
            this.t = t;
            this.x = x;
        }

        @Override
        public int hashCode() {
            int h1 = Long.hashCode(t);
            int h2 = Long.hashCode(x);
            return h1 ^ (h2 + 0x9e3779b9 + (h1 << 6) + (h1 >> 2));
        }

        @Override
        public boolean equals(Object o) {
            if (!(o instanceof CacheKey))
                return false;
            CacheKey other = (CacheKey) o;
            return this.t == other.t && this.x == other.x;
        }
    }

    static class CountPair {
        long pos, neg;

        CountPair(long p, long n) {
            this.pos = p;
            this.neg = n;
        }
    }

    static Map<CacheKey, CountPair> memo = new HashMap<>();

    static CountPair prefixCounts(long t, long x) {
        if (t == 0 || x < 0)
            return new CountPair(0, 0);
        if (t == 1)
            return new CountPair(1, 0);

        CacheKey key = new CacheKey(t, x);
        CountPair res = memo.get(key);
        if (res != null)
            return res;

        if ((t & 1) == 0) {
            long u = t / 2;
            long x1 = Math.floorDiv(x - 1, 4);
            long x3 = Math.floorDiv(x - 3, 4);
            CountPair a = prefixCounts(u, x1);
            CountPair b = prefixCounts(u, x3);
            if ((u & 1) == 0)
                res = new CountPair(a.pos + b.pos, a.neg + b.neg);
            else
                res = new CountPair(a.pos + b.neg, a.neg + b.pos);
        } else {
            long u = t / 2;
            long x0 = Math.floorDiv(x, 4);
            long x2 = Math.floorDiv(x - 2, 4);
            CountPair u0 = prefixCounts(u, x0);
            CountPair u2 = prefixCounts(u, x2);
            CountPair v0 = prefixCounts(u + 1, x0);
            CountPair v2 = prefixCounts(u + 1, x2);

            if ((u & 1) == 1)
                res = new CountPair(u2.pos + v0.pos, u0.neg + v2.pos);
            else
                res = new CountPair(u2.neg + v0.pos, u0.neg + v2.neg);
        }

        memo.put(key, res);
        return res;
    }

    static long countTotal(long t, long x) {
        CountPair p = prefixCounts(t, x);
        return p.pos + p.neg;
    }

    static long g(long t, long c) {
        if (t == 0 || c == 0 || c > t)
            return 0;
        long lo = 0;
        long hi = 1;
        while (countTotal(t, hi) < c) {
            hi *= 2;
        }
        while (lo < hi) {
            long mid = lo + (hi - lo) / 2;
            if (countTotal(t, mid) >= c)
                hi = mid;
            else
                lo = mid + 1;
        }
        return lo;
    }

    static String solve() {
        memo.clear();
        long[] fib = new long[46];
        fib[0] = 1;
        fib[1] = 1;
        for (int i = 2; i <= 45; i++) {
            fib[i] = fib[i - 1] + fib[i - 2];
        }

        long sum = 0;
        for (int t = 2; t <= 45; t++) {
            sum += g(fib[t], fib[t - 1]);
        }
        return Long.toString(sum);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}