Detailed mathematical approach

Problem Summary

Starting from \(S_0=0\), \(S_1=01\), and \(S_n=S_{n-1}S_{n-2}\), every finite factor of the resulting infinite Fibonacci word is called a Fibonacci subword. For each positive \(k\), exactly \(k+1\) distinct factors have length \(k\). Reading each factor as a base-10 integer, leading zeros included harmlessly, we must compute the sum of their squares

$$\Psi(k)=\sum_{u\in\mathcal F_k}\operatorname{val}_{10}(u)^2$$

for \(k=10^{18}\), modulo \(M=101001001\). Neither a word of that length nor its \(k+1\) factors can be constructed explicitly.

Mathematical Approach

Choose the first Fibonacci word longer than the window

Let \(f_n=|S_n|\). Then \(f_0=1\), \(f_1=2\), and

$$f_n=f_{n-1}+f_{n-2}.$$

Choose the smallest \(n\) for which \(L=f_n\gt k\), and put \(W=S_n=w_0w_1\cdots w_{L-1}\). A standard factor property of the Fibonacci word says that the \(L\) cyclic windows of length \(k\) in \(W\) contain all \(k+1\) distinct length-\(k\) factors. Since \(L\) windows represent only \(k+1\) values,

$$\delta=L-k-1$$

of those occurrences are redundant. With the standard indexing used here, one extra copy of each of the first \(\delta\) consecutive windows must be removed. Thus the problem becomes

$$\Psi(k)=\text{sum over all cyclic windows of }W -\text{sum of squares of the first }\delta\text{ windows}.$$

For \(k=10^{18}\), the index \(n\) is only \(O(\log k)\), although \(L\) itself is enormous.

Store words as a concatenation DAG

The implementation never expands \(S_n\). A leaf stores one bit, and an internal node stores two child identifiers and their total length. The recurrence \(S_n=S_{n-1}S_{n-2}\) therefore adds just one concatenation node per Fibonacci level. A prefix of any required length is obtained recursively: it lies wholly in the left child, or it is the entire left child followed by a prefix of the right child.

Concatenations and prefixes are memoized. Consequently, a prefix whose numerical length may be near \(10^{18}\) is still represented by only \(O(\log k)\) nodes.

Four summaries of a binary word

Work modulo \(M\), with \(b=10\). Because \(\gcd(10,M)=\gcd(99,M)=1\), both \(b^{-1}\) and \((b^2-1)^{-1}\) exist modulo \(M\). For a word \(X=x_0x_1\cdots x_{m-1}\), define

$$H_b(X)=\sum_{i=0}^{m-1}x_i b^i,\qquad R_b(X)=\sum_{i=0}^{m-1}x_i b^{m-1-i},$$

$$P_b(X)=\sum_{0\le i\lt j\lt m}x_i x_j b^{j-i},\qquad O(X)=\sum_{i=0}^{m-1}x_i.$$

\(R_b(X)\) is exactly the usual decimal value of \(X\) modulo \(M\). The pair summary \(P_b\) groups every pair of 1-bits by their distance.

If \(X=UV\), \(a=|U|\), and \(c=|V|\), all four values combine without inspecting a digit:

$$H_b(UV)=H_b(U)+b^aH_b(V),$$

$$R_b(UV)=b^cR_b(U)+R_b(V),$$

$$P_b(UV)=P_b(U)+P_b(V)+bR_b(U)H_b(V),$$

$$O(UV)=O(U)+O(V).$$

The cross term in \(P_b\) is correct because a bit \(i\) in \(U\) and a bit \(j\) in \(V\) are separated by \(a+j-i\) positions. The code caches these summaries for both bases \(b\) and \(b^{-1}\).

Range correlation queries on compressed words

Two additional recursive queries filter pairs without opening the word. The difference query computes

$$D_b(X,Y;\ell,h)= \sum_{\substack{i,j\\ \ell\le j-i\le h}}x_i y_j b^{j-i},$$

and the sum-index query computes

$$A_b(X,Y;\ell,h)= \sum_{\substack{i,j\\ \ell\le i+j\le h}}x_i y_j b^{i+j}.$$

When the requested range covers an entire node pair, the answer factors into two cached polynomial summaries. Otherwise the longer node is split, the interval is shifted by the child offset, and the two answers are combined. Memoization prevents the same node-pair/range state from being solved twice.

Cyclic correlations and the useful distances

For \(1\le d\lt L\), define the cyclic correlation

$$C_d=\sum_{i=0}^{L-1}w_iw_{(i+d)\bmod L}.$$

All nonzero cyclic distances can be packed into one polynomial:

$$Q_b=\sum_{d=1}^{L-1}C_db^d =P_b(W)+b^LP_{b^{-1}}(W).$$

The first term counts pairs that do not cross the end of \(W\); the second turns a linear distance \(j-i\) into its wrapped distance \(L-(j-i)\).

A length-\(k\) window can contain two positions only when their forward distance is \(1,\dots,k-1\). Put \(g=L-k\). The unwanted distances \(k,\dots,L-1\) are the reversals \(L-r\) of \(r=1,\dots,g\). Difference queries compute

$$E_b=\sum_{r=1}^{g}C_rb^r,$$

including both the ordinary and wrapped pieces. Hence the correlations actually used by windows are

$$U_b=Q_b-b^LE_{b^{-1}},\qquad U_{b^{-1}}=Q_{b^{-1}}-b^{-L}E_b.$$

Sum the squares of every cyclic window

Let \(V_s\) be the decimal value of the cyclic length-\(k\) window starting at \(s\). On expanding \(V_s^2\), diagonal digit terms and pairs of distinct positions separate cleanly.

Every 1-bit of \(W\) occupies every decimal place once across all cyclic windows. Its diagonal contribution is therefore

$$O(W)\sum_{t=0}^{k-1}b^{2t} =O(W)\frac{b^{2k}-1}{b^2-1}.$$

For a fixed cyclic distance \(d\), \(1\le d\lt k\), a pair appears in \(k-d\) relative placements. The sum of its place-value products is

$$\sum_{t=0}^{k-d-1}b^{2t+d} =\frac{b^{2k-d}-b^d}{b^2-1}.$$

After summing over all correlations, the square sum of all \(L\) cyclic windows is

$$T_{\mathrm{cyc}}=O(W)\frac{b^{2k}-1}{b^2-1} +\frac{2}{b^2-1}\left(b^{2k}U_{b^{-1}}-U_b\right)\pmod M.$$

This identity is the central compression step: exponentially many digit products collapse into two correlation evaluations and modular exponentiation.

Subtract the redundant windows without sliding one by one

If \(\delta=0\), the cyclic total already contains each factor once. Otherwise let \(V_0,\dots,V_{\delta-1}\) be the redundant prefix windows. The first value is \(R_b\) of a compressed prefix of length \(k\).

Put \(m=\delta-1\), and call the outgoing prefix \(x_0,\dots,x_{m-1}\). On this redundant run the entering block is the reversal of that prefix, so the bit entering at transition \(t\) is \(x_{m-1-t}\). The rolling decimal recurrence is therefore

$$V_{t+1}=bV_t-b^kx_t+x_{m-1-t},\qquad 0\le t\lt m.$$

Expanding this recurrence expresses every \(V_t\) in terms of \(V_0\), the outgoing prefix, and its reversal. Squaring and summing needs only counts of 1-bits, weighted bit pairs, the reversal diagonal \(\sum_t x_tx_{m-1-t}\), and two triangular index ranges. Those are precisely the cached summaries and \(A_b\) range queries. The method duplicated_window_sum evaluates the resulting telescoping/geometric expression modulo \(M\), then the final answer is

$$\Psi(k)\equiv T_{\mathrm{cyc}}-\sum_{t=0}^{\delta-1}V_t^2\pmod M.$$

Worked example: \(k=3\)

The first standard word longer than 3 is \(W=S_3=01001\), with \(L=5\) and \(\delta=5-3-1=1\). Its five cyclic windows are

$$010, 100, 001, 010, 101.$$

The first window \(010\) is the one redundant occurrence. Subtracting one copy leaves \(001,010,100,101\), so

$$\Psi(3)=1^2+10^2+100^2+101^2=20302.$$

How the Code Works

FibonacciSubwords builds the logarithmic standard-word DAG and caches powers of \(10\) and \(10^{-1}\). summarize implements the four concatenation formulas; difference_query filters correlations by \(j-i\), while sum_query handles the triangular products needed by the duplicate correction.

solve selects \(L\), forms \(Q_b\), removes separations too large for a \(k\)-window, evaluates \(T_{\mathrm{cyc}}\), and subtracts duplicated_window_sum. All arithmetic is reduced modulo \(101001001\).

The checkpoint suite constructs actual Fibonacci strings only for small inputs. It confirms \(\Psi(3)=20302\), compares the compressed solver with direct set enumeration for every \(1\le k\le50\), and verifies the supplied value \(\Psi(10)\equiv10699667\pmod M\).

Complexity Analysis

There are \(O(\log k)\) Fibonacci levels and prefix nodes. The memoized two-word range queries visit at most a quadratic number of relevant node-pair states, giving \(O(\log^2 k)\) structural work and \(O(\log^2 k)\) cached memory; modular exponentiation contributes \(O(\log k)\) multiplications per previously unseen exponent.

No operation is linear in \(k\), and no length-\(k\) string is stored. This is what makes \(k=10^{18}\) practical.

Footnotes and References

1. Problem page: Project Euler 1006 - Fibonacci Subwords
2. Fibonacci word: Wikipedia - Fibonacci word
3. Sturmian word and factor complexity: Wikipedia - Sturmian word
4. Modular multiplicative inverse: Wikipedia - Modular multiplicative inverse