Problem 691: Long Substring with Many Repetitions

Project Euler Problem 691 Solution

EulerSolve provides an optimized solution for Project Euler Problem 691, Long Substring with Many Repetitions, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For $0 \le i \lt n$, define $$a_i=\operatorname{popcount}(i)\bmod 2,\qquad b_i=\left\lfloor\frac{i+1}{\varphi}\right\rfloor-\left\lfloor\frac{i}{\varphi}\right\rfloor,\qquad c_i=a_i\oplus b_i,$$ where $\varphi=\frac{1+\sqrt{5}}{2}$. The binary word of interest is $$C_n=c_0c_1\dots c_{n-1}.$$ For each $k\ge 1$, let $L_n(k)$ be the maximum length of a substring of $C_n$ that occurs at least $k$ times, allowing overlaps. The task is to evaluate $$\sum_{\substack{k\ge 1\\L_n(k)>0}} L_n(k)$$ for $n=5{,}000{,}000$. Since $C_n$ has millions of positions and $O(n^2)$ distinct substring intervals, the solution must compress the repeated-substring structure instead of enumerating substrings directly. Mathematical Approach Once the word $C_n$ has been generated, the problem becomes a pure substring-frequency question for one fixed binary string. The decisive tool is a suffix automaton, because it stores all distinct substrings in linear space together with enough structure to recover how many times each substring appears. Step 1: Build the Binary Word The first ingredient, $a_i$, is the parity of the number of $1$-bits in the binary expansion of $i$, so it is the Thue-Morse bit at index $i$. The second ingredient, $b_i$, is a Beatty difference....

Detailed mathematical approach

Problem Summary

For $0 \le i \lt n$, define

$$a_i=\operatorname{popcount}(i)\bmod 2,\qquad b_i=\left\lfloor\frac{i+1}{\varphi}\right\rfloor-\left\lfloor\frac{i}{\varphi}\right\rfloor,\qquad c_i=a_i\oplus b_i,$$

where $\varphi=\frac{1+\sqrt{5}}{2}$. The binary word of interest is

$$C_n=c_0c_1\dots c_{n-1}.$$

For each $k\ge 1$, let $L_n(k)$ be the maximum length of a substring of $C_n$ that occurs at least $k$ times, allowing overlaps. The task is to evaluate

$$\sum_{\substack{k\ge 1\\L_n(k)>0}} L_n(k)$$

for $n=5{,}000{,}000$. Since $C_n$ has millions of positions and $O(n^2)$ distinct substring intervals, the solution must compress the repeated-substring structure instead of enumerating substrings directly.

Mathematical Approach

Once the word $C_n$ has been generated, the problem becomes a pure substring-frequency question for one fixed binary string. The decisive tool is a suffix automaton, because it stores all distinct substrings in linear space together with enough structure to recover how many times each substring appears.

Step 1: Build the Binary Word

The first ingredient, $a_i$, is the parity of the number of $1$-bits in the binary expansion of $i$, so it is the Thue-Morse bit at index $i$. The second ingredient, $b_i$, is a Beatty difference. Since $\frac{1}{\varphi}=\varphi-1\in(0,1)$, the sequence

$$\left\lfloor\frac{i}{\varphi}\right\rfloor$$

can increase by at most $1$ when $i$ is incremented, hence

$$b_i\in\{0,1\}.$$

Therefore $c_i=a_i\oplus b_i$ is again binary, and the entire input to the substring problem is a binary word $C_n$. Generating the word costs only $O(n)$ time, one position at a time.

Step 2: Compress All Substrings with a Suffix Automaton

A suffix automaton for $C_n$ groups substrings by their set of end positions. Each state $v$ represents an equivalence class of substrings that all end in exactly the same positions in $C_n$. Let $M(v)$ be the maximum length represented by state $v$. If the suffix link of $v$ points to $u$, then the represented lengths form the interval

$$M(u)+1,\ M(u)+2,\ \dots,\ M(v).$$

So every distinct substring belongs to exactly one state, and within that state the longest candidate has length $M(v)$. For a string of length $n$, a suffix automaton has at most $2n-1$ states, which is why the whole repeated-substring structure can be stored in linear size.

Step 3: Recover the Number of Occurrences of Each State

While the automaton is built online from left to right, each newly created non-clone state corresponds to one new suffix ending at the current position, so it starts with one end position. Clone states start from zero because they only split structure; they do not introduce a new occurrence by themselves.

After the full word has been inserted, the states are processed in decreasing order of $M(v)$. When a state $v$ contributes its total to its suffix-link parent, all longer substrings have already been counted, so the propagated value becomes the size of the end-position set for that entire class. Denoting this final count by $R(v)$, every substring represented by $v$ occurs exactly $R(v)$ times in $C_n$.

This count automatically includes overlapping repetitions, because different end positions are counted separately even if the corresponding substring occurrences overlap in the word.

Step 4: Convert State Data into $L_n(k)$

Fix a state $v$. Every substring represented by $v$ occurs exactly $R(v)$ times, but those substrings have several possible lengths inside the interval described above. For the quantity $L_n(k)$, only the longest one matters, because among all substrings with the same repetition count we want the maximum possible length. Therefore state $v$ contributes the candidate length $M(v)$ to the bucket for the exact repetition count $R(v)$.

If $B(t)$ denotes the largest substring length seen among states with exact occurrence count $t$, we update

$$B(R(v))=\max\bigl(B(R(v)),M(v)\bigr).$$

But the problem asks for “at least $k$ times”, not “exactly $k$ times”. So we transform exact counts into lower bounds by a suffix maximum:

$$L_n(k)=\max_{t\ge k} B(t).$$

A single pass from right to left over the array of counts computes all values $L_n(1),L_n(2),\dots,L_n(n)$.

Step 5: Sum the Positive Values

Once the entire array $L_n(k)$ is known, the required answer is simply

$$S_n=\sum_{\substack{k\ge 1\\L_n(k)>0}} L_n(k).$$

The sequence $L_n(k)$ is non-increasing in $k$: requiring more repetitions can never make the best repeated substring longer. Hence after the first zero appears, all later terms are also zero, so the final accumulation is straightforward.

Worked Example: $n=10$

For $i=0,1,\dots,9$, the two source sequences are

$$a_i=(0,1,1,0,1,0,0,1,1,0),$$

$$b_i=(0,1,0,1,1,0,1,0,1,1).$$

Taking bitwise XOR gives

$$C_{10}=0011001101.$$

Now inspect repeated substrings:

$L_{10}(1)=10$ because the whole word appears once. The substring $00110$ appears twice, at positions $1$ through $5$ and $5$ through $9$ in one-based indexing, so $L_{10}(2)=5$. The substring $01$ appears three times, and no substring of length $3$ appears three times, so $L_{10}(3)=2$. Finally, single symbols appear five times each, giving $L_{10}(4)=L_{10}(5)=1$, and for $k\ge 6$ the value is $0$.

Therefore

$$L_{10}(k)=(10,5,2,1,1,0,\dots),$$

and the sum of positive terms is

$$10+5+2+1+1=19.$$

This agrees with the small checkpoints used by the implementation, in particular $L_{10}(2)=5$ and $L_{10}(3)=2$.

How the Code Works

The C++, Python, and Java implementations generate the bits of $C_n$ on the fly, so no separate quadratic substring table is ever materialized. Each new bit extends a suffix automaton over the binary alphabet $\{0,1\}$; because the alphabet has size two, every state needs only two outgoing transitions, one suffix link, one maximum-length field, and one occurrence counter.

After the automaton is complete, the implementation counting-sorts states by their represented maximum length, then walks that order from longest to shortest so occurrence totals can be pushed through suffix links. Next it records, for each exact occurrence count, the largest maximum length attained by any state with that count. A backward suffix-maximum pass converts those exact-count buckets into the desired values $L_n(k)$. The last step sums the positive entries and stops once they become zero.

Complexity Analysis

Let $n$ be the length of the word. Building the binary sequence and extending the suffix automaton are both linear, and the automaton contains at most $2n-1$ states. Sorting states by length with counting sort is $O(n)$, propagating occurrence counts is $O(n)$, and the final exact-count and suffix-maximum passes are also $O(n)$. Thus the total running time is $O(n)$, and the memory usage is $O(n)$.

Footnotes and References

Problem page: https://projecteuler.net/problem=691
Suffix automaton overview: cp-algorithms — Suffix Automaton
Beatty sequence: Wikipedia — Beatty sequence
Thue-Morse sequence: Wikipedia — Thue-Morse sequence

Problem 691 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <vector>

namespace {

struct State {
    int next[2];
    int link;
    int len;
    int occ;
};

std::vector<int> build_L(int n) {
    std::vector<State> st(static_cast<std::size_t>(2 * n));

    int sz = 1;
    int last = 0;

    st[0].next[0] = st[0].next[1] = -1;
    st[0].link = -1;
    st[0].len = 0;
    st[0].occ = 0;

    const long double invphi = (std::sqrt(5.0L) - 1.0L) / 2.0L;
    std::uint64_t beat_prev = 0ULL;

    for (int i = 0; i < n; ++i) {
        const int a = __builtin_popcount(static_cast<unsigned int>(i)) & 1;

        const std::uint64_t beat_cur =
            static_cast<std::uint64_t>(std::floor((static_cast<long double>(i) + 1.0L) * invphi));
        const int b = static_cast<int>(beat_cur - beat_prev);
        beat_prev = beat_cur;

        const int c = a ^ b;

        const int cur = sz++;
        st[cur].next[0] = st[cur].next[1] = -1;
        st[cur].len = st[last].len + 1;
        st[cur].occ = 1;

        int p = last;
        while (p != -1 && st[p].next[c] == -1) {
            st[p].next[c] = cur;
            p = st[p].link;
        }

        if (p == -1) {
            st[cur].link = 0;
        } else {
            const int q = st[p].next[c];
            if (st[p].len + 1 == st[q].len) {
                st[cur].link = q;
            } else {
                const int clone = sz++;
                st[clone] = st[q];
                st[clone].len = st[p].len + 1;
                st[clone].occ = 0;

                while (p != -1 && st[p].next[c] == q) {
                    st[p].next[c] = clone;
                    p = st[p].link;
                }

                st[q].link = clone;
                st[cur].link = clone;
            }
        }

        last = cur;
    }

    st.resize(static_cast<std::size_t>(sz));

    std::vector<int> cnt_len(static_cast<std::size_t>(n + 1), 0);
    for (int i = 0; i < sz; ++i) {
        ++cnt_len[static_cast<std::size_t>(st[i].len)];
    }
    for (int i = 1; i <= n; ++i) {
        cnt_len[static_cast<std::size_t>(i)] += cnt_len[static_cast<std::size_t>(i - 1)];
    }

    std::vector<int> order(static_cast<std::size_t>(sz), 0);
    for (int i = sz - 1; i >= 0; --i) {
        const int l = st[i].len;
        order[static_cast<std::size_t>(--cnt_len[static_cast<std::size_t>(l)])] = i;
    }

    for (int i = sz - 1; i > 0; --i) {
        const int v = order[static_cast<std::size_t>(i)];
        const int p = st[v].link;
        if (p >= 0) {
            st[p].occ += st[v].occ;
        }
    }

    std::vector<int> best(static_cast<std::size_t>(n + 1), 0);
    for (int i = 1; i < sz; ++i) {
        const int occ = st[i].occ;
        if (occ <= n && st[i].len > best[static_cast<std::size_t>(occ)]) {
            best[static_cast<std::size_t>(occ)] = st[i].len;
        }
    }

    for (int k = n - 1; k >= 1; --k) {
        if (best[static_cast<std::size_t>(k)] < best[static_cast<std::size_t>(k + 1)]) {
            best[static_cast<std::size_t>(k)] = best[static_cast<std::size_t>(k + 1)];
        }
    }

    return best;
}

std::uint64_t sum_non_zero_L(int n) {
    const std::vector<int> L = build_L(n);
    std::uint64_t out = 0ULL;
    for (int k = 1; k <= n; ++k) {
        const int v = L[static_cast<std::size_t>(k)];
        if (v == 0) {
            break;
        }
        out += static_cast<std::uint64_t>(v);
    }
    return out;
}

}  // namespace

int main() {
    {
        const std::vector<int> L10 = build_L(10);
        assert(L10[2] == 5);
        assert(L10[3] == 2);
    }

    {
        const std::vector<int> L100 = build_L(100);
        assert(L100[2] == 14);
        assert(L100[4] == 6);
    }

    {
        const std::vector<int> L1000 = build_L(1000);
        assert(L1000[2] == 86);
        assert(L1000[3] == 45);
        assert(L1000[5] == 31);
    }

    assert(sum_non_zero_L(1000) == 2460ULL);

    std::cout << sum_non_zero_L(5'000'000) << "\n";
    return 0;
}

Python

import math

def solve():
    n = 5000000
    invphi = (math.sqrt(5) - 1) / 2

    # SAM
    class State:
        __slots__ = ['next', 'link', 'length', 'occ']
        def __init__(self):
            self.next = [-1, -1]; self.link = -1; self.length = 0; self.occ = 0

    st = [State() for _ in range(2*n)]
    sz = 1; last = 0
    beat_prev = 0

    for i in range(n):
        a = bin(i).count('1') & 1
        beat_cur = int(math.floor((i + 1) * invphi))
        b = beat_cur - beat_prev; beat_prev = beat_cur
        c = a ^ b

        cur = sz; sz += 1
        st[cur].length = st[last].length + 1; st[cur].occ = 1
        p = last
        while p != -1 and st[p].next[c] == -1:
            st[p].next[c] = cur; p = st[p].link
        if p == -1: st[cur].link = 0
        else:
            q = st[p].next[c]
            if st[p].length + 1 == st[q].length: st[cur].link = q
            else:
                clone = sz; sz += 1
                st[clone].next = st[q].next[:]; st[clone].link = st[q].link
                st[clone].length = st[p].length + 1; st[clone].occ = 0
                while p != -1 and st[p].next[c] == q:
                    st[p].next[c] = clone; p = st[p].link
                st[q].link = clone; st[cur].link = clone
        last = cur

    cnt_len = [0]*(n+1)
    for i in range(sz): cnt_len[st[i].length] += 1
    for i in range(1, n+1): cnt_len[i] += cnt_len[i-1]
    order = [0]*sz
    for i in range(sz-1, -1, -1):
        cnt_len[st[i].length] -= 1; order[cnt_len[st[i].length]] = i

    for i in range(sz-1, 0, -1):
        v = order[i]; p = st[v].link
        if p >= 0: st[p].occ += st[v].occ

    best = [0]*(n+1)
    for i in range(1, sz):
        o = st[i].occ
        if o <= n and st[i].length > best[o]: best[o] = st[i].length

    for k in range(n-1, 0, -1):
        if best[k] < best[k+1]: best[k] = best[k+1]

    return str(sum(v for v in best[1:] if v > 0))

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

Java

public class Euler691 {
    static class State {
        int[] next = { -1, -1 };
        int link;
        int len;
        int occ;
    }

    static int[] buildL(int n) {
        State[] st = new State[2 * n];
        for (int i = 0; i < st.length; i++) {
            st[i] = new State();
        }

        int sz = 1;
        int last = 0;

        st[0].link = -1;
        st[0].len = 0;
        st[0].occ = 0;

        double invphi = (Math.sqrt(5.0) - 1.0) / 2.0;
        long beatPrev = 0L;

        for (int i = 0; i < n; ++i) {
            int a = Integer.bitCount(i) & 1;

            long beatCur = (long) Math.floor((i + 1) * invphi);
            int b = (int) (beatCur - beatPrev);
            beatPrev = beatCur;

            int c = a ^ b;

            int cur = sz++;
            st[cur].len = st[last].len + 1;
            st[cur].occ = 1;

            int p = last;
            while (p != -1 && st[p].next[c] == -1) {
                st[p].next[c] = cur;
                p = st[p].link;
            }

            if (p == -1) {
                st[cur].link = 0;
            } else {
                int q = st[p].next[c];
                if (st[p].len + 1 == st[q].len) {
                    st[cur].link = q;
                } else {
                    int clone = sz++;
                    st[clone].next[0] = st[q].next[0];
                    st[clone].next[1] = st[q].next[1];
                    st[clone].link = st[q].link;
                    st[clone].len = st[p].len + 1;
                    st[clone].occ = 0;

                    while (p != -1 && st[p].next[c] == q) {
                        st[p].next[c] = clone;
                        p = st[p].link;
                    }

                    st[q].link = clone;
                    st[cur].link = clone;
                }
            }

            last = cur;
        }

        int[] cntLen = new int[n + 1];
        for (int i = 0; i < sz; ++i) {
            cntLen[st[i].len]++;
        }
        for (int i = 1; i <= n; ++i) {
            cntLen[i] += cntLen[i - 1];
        }

        int[] order = new int[sz];
        for (int i = sz - 1; i >= 0; --i) {
            int l = st[i].len;
            order[--cntLen[l]] = i;
        }

        for (int i = sz - 1; i > 0; --i) {
            int v = order[i];
            int p = st[v].link;
            if (p >= 0) {
                st[p].occ += st[v].occ;
            }
        }

        int[] best = new int[n + 1];
        for (int i = 1; i < sz; ++i) {
            int occ = st[i].occ;
            if (occ <= n && st[i].len > best[occ]) {
                best[occ] = st[i].len;
            }
        }

        for (int k = n - 1; k >= 1; --k) {
            if (best[k] < best[k + 1]) {
                best[k] = best[k + 1];
            }
        }

        return best;
    }

    public static String solve() {
        int n = 5000000;
        int[] L = buildL(n);
        long out = 0L;
        for (int k = 1; k <= n; ++k) {
            int v = L[k];
            if (v == 0) {
                break;
            }
            out += v;
        }
        return Long.toString(out);
    }

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