Problem 238: Infinite String Tour

Project Euler Problem 238 Solution

EulerSolve provides an optimized solution for Project Euler Problem 238, Infinite String Tour, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Start with $s_0=14025256$ and iterate $$s_{n+1}=s_n^2 \pmod{20300713}.$$ Write the decimal expansions of $s_0,s_1,s_2,\dots$ one after another to obtain an infinite digit string $w=w_1w_2w_3\cdots$. For each positive integer $k$, define $p(k)$ to be the smallest starting position $z\ge 1$ for which there exists an ending position $y\ge z$ such that the digit sum of the block $w_z w_{z+1}\dots w_y$ is exactly $k$. The problem asks for $$\sum_{k=1}^{2\cdot 10^{15}} p(k).$$ A direct attack is impossible: the outer summation is enormous, and the string itself is infinite. The solution used by the implementations turns the problem into finite data extracted from one period of the quadratic generator. Mathematical Approach The key objects are the state cycle of the modular recurrence, the resulting periodic digit block, and the prefix sums of that block taken modulo one full-period digit sum. The modular squaring sequence is purely periodic The recurrence evolves in a finite state space, so some state must eventually repeat. Once a state repeats, the future trajectory repeats as well. In this problem the repetition is especially clean: the first repeated state is the initial seed itself, so there is no transient prefix before the cycle starts. Therefore the sequence of generated integers is a pure cycle from $s_0$ onward....

Detailed mathematical approach

Problem Summary

Start with $s_0=14025256$ and square repeatedly modulo $20300713$, so that $s_{n+1}=s_n^2\bmod 20300713$; this is the Blum Blum Shub generator. ^[2] Writing $s_0,s_1,s_2,\dots$ one after another in ordinary decimal notation, without leading zeros, gives an infinite string of digits $w=w_1w_2w_3\cdots$ that begins $14025256\,7410149\,5847003\,8053646\cdots$. For a positive integer $k$, $p(k)$ is the position at which the earliest substring, a block of consecutive digits, with digit sum $k$ starts, and $p(k)=0$ if no substring has digit sum $k$. The statement gives $p(1)=p(5)=p(7)=1$ from the substrings $1$, $14$ and $1402$, then $p(4)=p(6)=p(11)=2$ and $p(2)=p(9)=3$, and the checkpoint $\sum_{0\lt k\le10^3}p(k)=4742$. We must find the sum of $p(k)$ over $0\lt k\le2\times10^{15}$. ^[1]

The whole solution in one sentence: $w$ repeats one long block of digits, so whether a substring starting at a given position can have digit sum $k$ depends only on the remainder of $k$ modulo the digit sum $S$ of the block, which makes $p$ repeat with period $S$, lets a few table lookups find each $p(k)$, and turns the sum over $2\times10^{15}$ targets into whole copies of one period plus a remainder.

Mathematical Approach

We first turn substring sums into differences of running totals, then show that $w$ is periodic, and finally reduce every question about the infinite string to a question about one period.

1. Substring sums as differences of running totals

Let $T_0=0$ and $T_n=w_1+\dots+w_n$; these running totals are the prefix sums of the digits. ^[7] The substring from position $z$ to position $y\ge z$ has digit sum

$$w_z+w_{z+1}+\dots+w_y=T_y-T_{z-1}.$$

Call $z$ a good start for $k$ if some substring starting at $z$ has digit sum $k$, so that $p(k)$ is the smallest good start, or $0$ if there is none. Digits are never negative, so the totals never decrease, and for $k\ge1$ an equation $T_y=T_{z-1}+k$ forces $y\ge z$. Hence $z$ is a good start for $k$ exactly when $T_{z-1}+k$ is itself one of the totals $T_0,T_1,T_2,\dots$.

The first sixteen digits $1,4,0,2,5,2,5,6,7,4,1,0,1,4,9,5$ give the totals $T_0,\dots,T_{16}$:

$0,1,5,5,7,12,14,19,25,$
$32,36,37,37,38,42,51,56.$

For $k=7$ the start $1$ is good, because $T_4=7$: the substring is $1402$. For $k=6$ it is not, because the totals jump from $5$ to $7$, but $T_1+6=7=T_4$ makes $z=2$ good, so $p(6)=2$ with the substring $402$. For $k=2$ the starts $1$ and $2$ would need the totals $2$ and $3$, which never occur, while $T_2+2=7=T_4$ gives $p(2)=3$, the substring $02$ of the statement.

2. A small target that starts late

A substring with digit sum $3$ can contain only the digits $0$ to $3$. Among the first $19$ digits $1402525674101495847$ the maximal runs of such digits are $1$, $02$, $2$ and $101$, and any substring with digit sum $3$ would have to lie inside one of them; but their digit sums are only $1,2,2,2$. The next three digits are $w_{20}w_{21}w_{22}=003$, so $p(3)=20$: the starts $20$, $21$ and $22$ all reach the sum $3$, and the earliest of them counts. This is a general rule: if $w_z=0$ then $T_z=T_{z-1}$, so the starts $z$ and $z+1$ are good for the same targets, and a start that comes right after a zero digit is never the smallest good start. Small targets can therefore have late starts, but Section 10 shows that no search is long.

3. The generator returns to its seed

Squaring modulo $M=20300713$ maps the finite set $\{0,1,\dots,M-1\}$ into itself, so the sequence must repeat a value sooner or later and from then on runs around a cycle; in general a tail may lead into the cycle, a picture shaped like the Greek letter rho. ^[3] Here there is no tail. The modulus is $M=4127\cdot4919$, and both primes leave the remainder $3$ on division by $4$. Such a product is a Blum integer, and squaring permutes its quadratic residues, the numbers coprime to $M$ that are congruent to a square modulo $M$. ^[4] The seed is one of them: $\gcd(s_0,M)=1$, and by Euler's criterion the congruences $s_0^{2063}\equiv1\pmod{4127}$ and $s_0^{2459}\equiv1\pmod{4919}$ say that $s_0$ is a square modulo each prime; ^[5] the Chinese remainder theorem combines the two square roots into one modulo $M$. ^[6] So every $s_n$ is such a residue, and if $s_j=s_i$ with $0\lt i\lt j$, then $s_{j-1}^2\equiv s_{i-1}^2\pmod M$ and the permutation property give $s_{j-1}=s_{i-1}$. The first value to repeat is therefore $s_0$ itself, and the sequence is a pure cycle. The simulation agrees: $s_0,\dots,s_{2534197}$ are distinct and $s_{2534198}=s_0$.

4. Why the cycle has $2\,534\,198$ numbers

By induction $s_n\equiv s_0^{2^n}\pmod M$. ^[2] For $a$ coprime to $m$, the order $\operatorname{ord}_m(a)$ of $a$ modulo $m$ is the smallest $e\ge1$ with $a^e\equiv1\pmod m$; for example, the powers of $4$ modulo $7$ are $4,2,1,4,2,1,\dots$, so $\operatorname{ord}_7(4)=3$. ^[6] If $a^n\equiv1$, dividing $n$ by the order leaves a remainder $t$ with $a^t\equiv1$, so $t=0$: thus $a^n\equiv1$ exactly when the order divides $n$. Both primes are safe primes, $4127=2\cdot2063+1$ and $4919=2\cdot2459+1$ with $2063$ and $2459$ prime. By Euler's criterion the order of $s_0$ modulo $4127$ divides the prime $2063$, and it is not $1$ because $s_0\equiv1710\pmod{4127}$, so it equals $2063$; likewise the order modulo $4919$ is $2459$, since $s_0\equiv1187\pmod{4919}$. Hence $s_n=s_0$, that is $s_0^{2^n-1}\equiv1\pmod M$, holds exactly when $2063$ and $2459$ both divide $2^n-1$, which means that $n$ is a multiple of $\operatorname{ord}_{2063}(2)$ and of $\operatorname{ord}_{2459}(2)$.

Because $2063\equiv7\pmod8$, the second supplement of quadratic reciprocity makes $2$ a square modulo $2063$, so $2^{1031}\equiv1$ by Euler's criterion; the order divides the prime $1031$ and is not $1$, so it equals $1031$. Because $2459\equiv3\pmod8$, $2$ is not a square modulo $2459$, so $2^{1229}\equiv-1$. As $1229$ is prime, the divisors of $2458=2\cdot1229$ are $1$, $2$, $1229$ and $2458$; this congruence rules out $1$ and $1229$, and $2^2=4$ rules out $2$, so the order is $2458$. ^[5] The length of the cycle is therefore

$$\operatorname{lcm}(1031,2458)=1031\cdot2458=2\,534\,198,$$

the number found by the simulation. Safe primes with a small $\gcd\bigl((p-3)/2,(q-3)/2\bigr)$, here $\gcd(2062,2458)=2$, are what the description of the generator recommends for a long cycle. ^[2]

5. One period of digits

Let $d_1d_2\cdots d_L$ be the digits of $s_0,s_1,\dots,s_{2534197}$ written one after another. Since the sequence is a pure cycle, $w$ is this block repeated forever: $w_{mL+i}=d_i$ for $m\ge0$ and $1\le i\le L$. The numbers of the cycle have between $1$ and $8$ digits, and the block has $L=18\,886\,117$ digits with digit sum $S=80\,846\,691$, about $4.28$ per digit. With $P_0=0$ and $P_i=d_1+\dots+d_i$ we have $P_L=S$ and

$$T_{mL+i}=mS+P_i\qquad(m\ge0,\ 0\le i\le L).$$

In particular any $L$ consecutive digits of $w$ have digit sum $S$, because they meet every position of the block exactly once. The block ends with the digit $6$, the last digit of $s_{2534197}=782226$, so $P_{L-1}=S-6$, and $0\le P_i\lt S$ for every $i\lt L$.

6. The residue test

Let $R=\{P_0,P_1,\dots,P_{L-1}\}$, viewed as a set of residues modulo $S$. By Section 5 the totals of $w$ are exactly the numbers $mS+r$ with $m\ge0$ and $r\in R$. Every integer $x\ge0$ can be written in exactly one way as $x=mS+r$ with $0\le r\lt S$, namely with $r=x\bmod S$, and all elements of $R$ lie in that range, so

$$x\in\{T_0,T_1,T_2,\dots\}\iff(x\bmod S)\in R\qquad(x\ge0).$$

Now take a start $z$ and write $z-1=mL+b$ with $0\le b\lt L$. Then $T_{z-1}=mS+P_b$, and adding $mS$ does not change a residue, so by Section 1 the start $z$ is good for $k$ exactly when

$$(P_b+k)\bmod S\in R,\qquad b=(z-1)\bmod L.$$

All three programs evaluate exactly this test, one addition, one reduction and one table lookup, even when the substring it certifies is far longer than one period.

7. Three consequences

Periodicity in $k$. The test depends on $k$ only through $k\bmod S$, so $k$ and $k+S$ have the same good starts, and $p(k+S)=p(k)$ for every $k\ge1$, zeros included. Concretely, a substring with digit sum $k$ followed by the next $L$ digits has digit sum $k+S$: the substring $003$ at positions $20$ to $22$, extended to position $22+L$, has digit sum $S+3$, and indeed $p(S+3)=p(3)=20$.
One period of starts suffices. The test depends on $z$ only through $(z-1)\bmod L$, so $z$ and $z+L$ are good for the same targets. If none of the starts $1,\dots,L$ is good, then none at all is and $p(k)=0$; otherwise $1\le p(k)\le L$. Trying $z=1,2,3,\dots$ until a start is good therefore finds $p(k)$, provided that some start is good.
The first start. Since $P_0=0$, we have $p(k)=1$ exactly when $k\bmod S\in R$; in particular $p(mS)=1$ for every $m\ge1$.

8. From $2\times10^{15}$ targets to one period

Divide $N=2\times10^{15}$ by $S$: $N=qS+r$ with $q=24\,738\,180$ and $r=5\,637\,620$. The targets $1\le k\le qS$ form $q$ runs of $S$ consecutive values, and by Section 7 every run has the same sum as the first one, so

$$\sum_{k=1}^{N}p(k)=q\sum_{k=1}^{S}p(k)+\sum_{k=1}^{r}p(k).$$

Only the $S$ values $p(1),\dots,p(S)$ are needed, about $8\cdot10^7$ of them instead of $2\times10^{15}$.

9. Worked example with a tiny generator

The same steps work for any purely periodic digit string whose period has a positive digit sum, and a tiny generator lets us check them completely by hand. Take $M=77=7\cdot11$, again a Blum integer, and the seed $4$. The cycle is $4\to16\to25\to9\to4$, since $9^2=81\equiv4\pmod{77}$, so $w=416259\,416259\cdots$ with $L=6$, $S=27$, $P_0,\dots,P_5=0,4,5,11,13,18$ and $R=\{0,4,5,11,13,18\}$. A start $z\le6$ has $b=z-1$, and by Section 6 it is good for the targets $1\le k\le27$ with $(P_b+k)\bmod27\in R$:

$z=1$, $P_0=0$: $k=4,5,11,13,18,27$;
$z=2$, $P_1=4$: $k=1,7,9,14,23,27$;
$z=3$, $P_2=5$: $k=6,8,13,22,26,27$;
$z=4$, $P_3=11$: $k=2,7,16,20,21,27$;
$z=5$, $P_4=13$: $k=5,14,18,19,25,27$;
$z=6$, $P_5=18$: $k=9,13,14,20,22,27$.

The first start in whose list each $k$ appears gives $p(1),\dots,p(27)$:

$2,4,0,1,1,3,2,3,2,0,1,0,1,2,$
$0,4,0,1,5,4,4,3,2,0,5,3,1.$

For instance $p(19)=5$: the substring $5941$ at positions $5$ to $8$ has digit sum $19$, while the starts $1$ to $4$ need the residues $19$, $23$, $24$ and $3$, none of which is in $R$. The six targets $3,10,12,15,17,24$ appear in no list, so $p(k)=0$ for them and for every $k$ that differs from one of them by a multiple of $27$. For $k=3$ this is clear directly: a substring with digit sum $3$ could use only the digits $1$ and $2$, and these are never next to each other. One period contributes $p(1)+\dots+p(27)=54$, and since $100=3\cdot27+19$, Section 8 gives $\sum_{k\le100}p(k)=3\cdot54+32=194$, where $32=p(1)+\dots+p(19)$; a direct search over the string gives the same value.

10. How far the searches go

The set $R$ has $17\,207\,333$ elements: consecutive $P_i$ are equal exactly when the digit between them is $0$, and $1\,678\,784$ of the $L$ digits are zeros, none of them the last one. So $R$ contains about $21.3\%$ of the residues modulo $S$, and a start is good for roughly one target in five. The toy example shows that a target can be unattainable, but here a finite computation rules this out: the good residues of the start $b+1$ are the elements of $R$ minus $P_b$, reduced modulo $S$, and for $0\le b\lt89$ these sets together cover all $S$ residues, so by Section 7 every positive $k$ is attainable and $p(k)\le89$. The searches are short: in one period $p(k)=1$ for $17\,207\,333$ values of $k$, as Section 7 predicts, $p(k)\le10$ for $72\,793\,414$ of them, about $90\%$, $p(k)\gt50$ for only $1298$, and the value $89$ is reached first at $k=68\,290\,826$.

How the Code Works

The C++ program ^[9] follows Sections 3 to 8 in order. BBSSequence::generate iterates $s\mapsto s^2\bmod M$ in unsigned 64-bit arithmetic, which is exact because $s^2\lt M^2\approx4.1\cdot10^{14}$. It records the index of every value in the std::unordered_map seen and stops at the first value that is already present, the hash-table method of cycle detection. ^[3] Each value is converted with std::to_string, so no leading zeros appear, and its digits are appended to w_digits (and to a string w_str that is not used later). The program prints the index of the repeat, $2534198$, and the index at which the repeated value was first seen, $0$, which confirms Section 3; for a nonzero index it would print a warning, because everything that follows assumes a pure cycle. The vector prefix_sum then holds $P_0,\dots,P_L$, and its last entry is $S$.

solve_full builds the residue test of Section 6. R is a std::vector<bool> with $S+1$ entries, which common implementations pack into single bits, about 10 MB here; ^[8] the loop over $0\le i\le L$ sets R[prefix_sum[i] % S], where $i=L$ marks $0$ a second time. P_mod_S[i] stores $P_i\bmod S$ for $0\le i\lt L$, which by Section 5 is $P_i$ itself. For every $k$ from $1$ to $S$ the program then tries $z=1,2,\dots$: it reads P_mod_S[(z - 1) % L], adds $k$, subtracts $S$ once if the sum is at least $S$, which suffices because $0\le P_b\lt S$ and $1\le k\le S$, so the sum is below $2S$, and stops at the first $z$ whose target is set in R, storing $z$ in p_vals[k]. The loop has no upper limit on $z$: it relies on every $k$ being attainable, which Section 10 confirms for this string. For the toy string of Section 9 it would never stop at $k=3$.

Finally the program sets LIMIT to $2\cdot10^{15}$, computes num_periods $=q$ and remainder $=r$, adds p_vals over $1\le k\le S$ and over $1\le k\le r$, and combines the two sums as in Section 8. It prints $S$, $L$, both sums and the result. Every $p(k)$ is at most $89$, so the result is at most $89N\lt1.8\cdot10^{17}$, far below $2^{63}$, and the unsigned 64-bit sums cannot overflow. The constants TARGET_K_SMALL and EXPECTED_SMALL_SUM, $1000$ and $4742$, are declared but never used: the program does not test the checkpoint of the statement.

The Python program has the same structure in a single function solve: a dictionary seen, lists w_digits, prefix_sum and P_mod_S, a bytearray with one byte per entry of R, and a list p_vals; its integers have unlimited precision. The Java program uses a HashMap<Long, Integer>, lists of boxed values for the digits and prefix sums, a byte[] for R and long[] arrays for the residues and the values $p(k)$; its signed 64-bit long is safe by the same bound. Both print only the result. They also stop at the first repeated value and use all digits generated before it, so they too rely on Section 3, but without the warning.

Complexity and Verification

Generating the cycle takes $2\,534\,198$ squarings and hash-table operations, the digits and prefix sums take $O(L)$ further steps, and the flag table R takes $O(S)$. The search for $p(k)$ makes $p(k)$ tests, so one period costs $p(1)+\dots+p(S)$ tests of constant cost. The crude bound from Section 7 is $SL\approx1.5\cdot10^{15}$ and Section 10 gives at most $89S\approx7.2\cdot10^9$, but the actual number is about $4\cdot10^8$, just under five tests per $k$. A rough model explains this: if every new start were good independently with probability $\rho=|R|/S\approx0.213$, about $1/\rho\approx4.7$ different starts would be needed on average, the largest of $S$ such waiting times would be near $\ln S/(-\ln(1-\rho))\approx76$, and zero digits, which repeat a start, add about $10\%$ to the positions, in line with the measured values. The two final sums take $O(S)$ steps. Memory is $O(S+L)$. The largest structure is p_vals, $S+1$ entries of 8 bytes, about 647 MB in C++ and Java; prefix sums, residues, digits and the hash map come on top, Java's boxed values and Python's integer objects add overhead, and Python's final sum over a slice copies the list again. On a desktop computer with an Apple M4 Max (Clang -O2, OpenJDK 17, CPython 3.11) the C++ and Java programs took about one second and the Python program about $22$ seconds, with peak memory of about 1.2 GB, 2.3 GB and 3.4 GB. Since every $p(k)$ is at most $89$, one byte per entry would do, and since $r\lt S$, a running sum recorded when $k$ reaches $r$ would avoid the table altogether.

Apart from the warning about a tail, the published programs contain no checks. For this article we verified the following separately:

A direct search from the definition over the explicit string reproduces the examples of the statement and the checkpoint $4742$, gives $p(3)=20$, and agrees with the residue test for every $k\le200\,000$.
An independent sweep over the starts $z=1,2,\dots$, which marks with word-wide bit operations all residues $k\bmod S$ for which $z$ is good, has assigned every residue after $z=89$ and agrees entry by entry with a re-implementation of the published search over all $S$ values of $k$.
A two-pointer search from the definition over an explicitly generated prefix of $w$, four full cycles produced by continued squaring, uses neither the residue test nor any periodicity. It agrees with the residue test for $57$ targets up to about $2.95S$, among them $k=1,\dots,30$, $S-2,\dots,S+3$, $S+458$, $2S+3$ and $68\,290\,826+S$, where the value is $89$.
The simulated cycle length agrees with Section 4, and for the toy generator of Section 9 the residue test agrees with a direct search for all $k\le108$, including the sum $194$.
The C++, Python and Java programs print the same result as the sweep.

Footnotes and References

The references below supply the background facts used above. The residue test, the periodicity of $p$, the length of the cycle and the statistics of the searches are derived or computed explicitly in this article.

Project Euler 238 — Infinite String Tour. The official statement defines the sequence s₀ = 14025256, sₙ₊₁ = sₙ² mod 20300713, the string w of the concatenated numbers and p(k), gives the examples p(1) = p(5) = p(7) = 1, p(4) = p(6) = p(11) = 2 and p(2) = p(9) = 3 and the checkpoint 4742 for 0 < k ≤ 10³, and asks for the sum of p(k) over 0 < k ≤ 2 × 10¹⁵.
Wikipedia — Blum Blum Shub. The generator is xₙ₊₁ = xₙ² mod M with M = pq for two primes congruent to 3 modulo 4, which guarantees that each quadratic residue has one square root that is also a quadratic residue. The article gives xᵢ = x₀^(2ⁱ mod λ(M)) mod M with the Carmichael function λ, and it recommends safe primes with a small gcd((p − 3)/2, (q − 3)/2) because this makes the cycle long. The generator goes back to L. Blum, M. Blum and M. Shub, SIAM Journal on Computing 15 (1986).
Wikipedia — Cycle detection. Iterating a function from a finite set to itself must eventually produce the same value twice, after which the values cycle along a path into a loop shaped like the Greek letter ρ. Storing the computed values in a hash table finds the first repeat, at the cost of memory that grows with the number of values.
Wikipedia — Blum integer. A Blum integer is a product of two distinct primes that are both congruent to 3 modulo 4. Every quadratic residue coprime to it has four square roots, exactly one of which is again such a residue, and squaring is a permutation of these residues.
Wikipedia — Euler's criterion: for an odd prime p and a coprime to p, a^((p − 1)/2) ≡ 1 (mod p) if a is a quadratic residue and ≡ −1 (mod p) if it is not. Wikipedia — Quadratic reciprocity: by the second supplement, x² ≡ 2 (mod p) is solvable exactly when p ≡ ±1 (mod 8).
Wikipedia — Multiplicative order: for a coprime to n, the order of a modulo n is the smallest positive k with aᵏ ≡ 1 (mod n); the powers of 4 modulo 7 give the order 3, and the order always divides φ(n). Wikipedia — Chinese remainder theorem: the remainders of an integer modulo pairwise coprime integers determine its remainder modulo their product; Sections 3 and 4 use this to combine statements modulo 4127 and 4919.
Wikipedia — Prefix sum. The prefix sums of a sequence are its running totals; Section 1 writes every substring sum as the difference of two of them.
cppreference — std::vector<bool>. The specialisation for bool is possibly space-efficient; how this is done is implementation-defined, and one possible optimisation stores each element in a single bit.
C++ — Euler238.cpp. The linked immutable C++ revision contains BBSSequence::generate, which builds one period and its prefix sums, and solve_full, which builds the residue table, finds p(k) for one period and combines the sums. It is the implementation record for this article.

Problem 238 source code

C++

#include <atomic>
#include <iostream>
#include <string>
#include <thread>
#include <unordered_map>
#include <vector>

// Type aliases for easier reading
using u64 = unsigned long long;

// BBS Parameters
const u64 MODULUS = 20300713;
const u64 SEED = 14025256;
const u64 TARGET_K_SMALL = 1000;
const u64 EXPECTED_SMALL_SUM = 4742;

// Function to generate the sequence of digits and find the period
struct BBSSequence {
  std::string w_str;
  std::vector<int> w_digits;
  std::vector<u64> prefix_sum;
  u64 total_sum_per_period;

  void generate() {
    u64 s = SEED;
    std::unordered_map<u64, int> seen;
    int index = 0;

    // We know it will cycle eventually. Since modulus is ~2e7, cycle is shorter
    // than that. Actually, we need to generate until we hit a repeated state
    // s_n. The problem statement says s_0 = 14025256, s_{n+1} = s_n^2 mod
    // 20300713. Digits are concatenated.

    while (seen.find(s) == seen.end()) {
      seen[s] = index;

      // Append digits of s to w
      std::string s_str = std::to_string(s);
      w_str += s_str;
      for (char c : s_str) {
        w_digits.push_back(c - '0');
      }

      s = (s * s) % MODULUS;
      index++;
    }

    // Note: The problem implies w starts from s0.
    // "Concatenate these numbers s0s1s2... to create a string w"
    // Cycle detected at state 's'.
    // However, for the purpose of the sum of digits, we treat w as the sequence
    // of digits. The "infinite length" string is formed by repeating the cycle.
    // Wait, does the string w start repeating exactly?
    // s0 -> s1 -> ... -> sk -> ... -> sk+L -> ... where sk = sk+L
    // The digits from s0 to sk-1 are the pre-period.
    // But usually BBS cycles are pure if s0 is a quadratic residue?
    // Or if we enter a cycle, the digits corresponding to the cycle repeat.
    // Let's verify if s0 is part of the cycle or pre-period.

    std::cout << "Cycle detected at index " << index << " (value: " << s << ")"
              << std::endl;
    std::cout << "First seen at index " << seen[s] << std::endl;

    // For this specific problem, let's assume valid range logic later.
    // Just construct the full digit sequence for one period if it's pure,
    // or enough to cover the non-periodic part + one period.

    // Prefix sums
    prefix_sum.push_back(0);
    u64 current_sum = 0;
    for (int d : w_digits) {
      current_sum += d;
      prefix_sum.push_back(current_sum);
    }
    total_sum_per_period =
        current_sum; // This might be slightly wrong if pre-period exists.
    // If pre-period exists, the "infinite string" isn't just repeating
    // 'total_sum_per_period' from the start. But let's check if pre-period is
    // empty (0).
    if (seen[s] != 0) {
      std::cerr
          << "Warning: Pre-period detected. Simplification assumption failed."
          << std::endl;
    }
  }
};

// Full solver
void solve_full() {
  BBSSequence bbs;
  std::cout << "Generating sequence..." << std::endl;
  bbs.generate();

  u64 S = bbs.total_sum_per_period;
  size_t L = bbs.w_digits.size();

  std::cout << "S (Sum): " << S << std::endl;
  std::cout << "L (Length): " << L << std::endl;

  // Build set R of residues present in prefix sums
  // std::vector<bool> is space efficient (1 bit per bool)
  std::cout << "Building residue table..." << std::endl;
  std::vector<bool> R(S + 1, false);
  std::vector<u64> P_mod_S(L); // Cache P[i] % S for fast access

  // P[0] = 0
  R[0] = true;
  // Fill R and P_mod_S
  // prefix_sum has L+1 entries: 0, w1, w1+w2...
  // prefix_sum[i] corresponds to sum of first i digits.
  // P_mod_S[i] should store prefix_sum[i] % S.
  // NOTE: In the search logic, we use P[z-1].
  // z=1 -> P[0]. z=2 -> P[1].
  // So we need P_mod_S to store the sequence of prefix sums modulo S.

  // Reconstruct just the mods to save memory if needed (though bbs has them)
  // We already have bbs.prefix_sum.
  for (size_t i = 0; i < L; ++i) {
    u64 val = bbs.prefix_sum[i] % S;
    P_mod_S[i] = val;
    // The residue is present.
    // Wait, R represents the set of prefix sums that exist.
    // We iterate through period.
    // prefix_sum[0]...prefix_sum[L-1].
    // prefix_sum[L] = S == 0 mod S.
    // So we just mark all prefix_sum[i] % S.
  }

  for (size_t i = 0; i <= L; ++i) { // Include L to ensure 0/S is marked
    R[bbs.prefix_sum[i] % S] = true;
  }

  std::cout << "Computing p(k) for all k..." << std::endl;

  // We need to sum p(k) for k = 1 to S.
  // Store them to handle the full range multiplication.

  std::vector<u64> p_vals(S + 1);
  std::atomic<u64> processed_count(0);

  unsigned int num_threads = std::thread::hardware_concurrency();
  std::vector<std::thread> threads;
  u64 chunk_size = S / num_threads;

  auto worker = [&](u64 start_k, u64 end_k) {
    for (u64 k = start_k; k <= end_k; ++k) {
      // Find minimal z such that (P[z-1] + k) % S is in R.
      u64 z = 1;
      while (true) {
        // accessing P_mod_S[(z-1) % L]
        // and comparing with k
        u64 start_residue = P_mod_S[(z - 1) % L];
        u64 target = (start_residue + k);
        if (target >= S)
          target -= S; // manually mod S

        if (R[target]) {
          p_vals[k] = z;
          break;
        }
        z++;
      }
    }
    processed_count += (end_k - start_k + 1);
  };

  for (unsigned int i = 0; i < num_threads; ++i) {
    u64 start = 1 + i * chunk_size;
    u64 end = (i == num_threads - 1) ? S : (start + chunk_size - 1);
    if (start <= end) {
      threads.emplace_back(worker, start, end);
    }
  }

  // Monitor progress
  // Move to join
  for (auto &t : threads) {
    if (t.joinable())
      t.join();
  }

  std::cout << "All p(k) computed." << std::endl;

  // Calculate total answer
  // Limit: 2 * 10^15
  u64 LIMIT = 2000000000000000ULL;
  u64 num_periods = LIMIT / S;
  u64 remainder = LIMIT % S;

  u64 sum_p_k_one_period = 0;
  u64 sum_p_k_remainder = 0;

  // Parallel reduction for sum? Or just linear is fast enough (10^8 adds ~
  // 0.1s)
  for (u64 k = 1; k <= S; ++k) {
    sum_p_k_one_period += p_vals[k];
  }

  for (u64 k = 1; k <= remainder; ++k) {
    sum_p_k_remainder += p_vals[k];
  }

  u64 total_ans = num_periods * sum_p_k_one_period + sum_p_k_remainder;

  std::cout << "Sum for one period: " << sum_p_k_one_period << std::endl;
  std::cout << "Remainder sum (" << remainder << "): " << sum_p_k_remainder
            << std::endl;
  std::cout << "Total Answer: " << total_ans << std::endl;
}

int main() {
  // Run validation if argument provided or just default
  // For now, let's run full solver which includes generation
  // Maybe verify small case again inside logic?
  // The previous main verified small case.
  // Let's keep a quick check.

  // Re-instantiate for small check? No, waste of time.
  // Just run the fast solver.
  solve_full();
  return 0;
}

Python

def solve():
    MODULUS = 20300713
    SEED = 14025256

    s = SEED
    seen = {}
    index = 0
    w_digits = []

    while s not in seen:
        seen[s] = index
        s_str = str(s)
        for c in s_str:
            w_digits.append(int(c))
        s = (s * s) % MODULUS
        index += 1

    L = len(w_digits)
    prefix_sum = [0] * (L + 1)
    for i in range(L):
        prefix_sum[i + 1] = prefix_sum[i] + w_digits[i]
    S = prefix_sum[L]

    R = bytearray(S + 1)
    for i in range(L + 1):
        R[prefix_sum[i] % S] = 1

    P_mod_S = [prefix_sum[i] % S for i in range(L)]

    p_vals = [0] * (S + 1)
    for k in range(1, S + 1):
        z = 1
        while True:
            start_residue = P_mod_S[(z - 1) % L]
            target = start_residue + k
            if target >= S:
                target -= S
            if R[target]:
                p_vals[k] = z
                break
            z += 1

    LIMIT = 2_000_000_000_000_000
    num_periods = LIMIT // S
    remainder = LIMIT % S

    sum_one = sum(p_vals[1:S+1])
    sum_rem = sum(p_vals[1:remainder+1])

    answer = num_periods * sum_one + sum_rem
    return str(answer)

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

Java

import java.util.*;
import java.util.concurrent.*;

public class Euler238 {
    static final long MODULUS = 20300713;
    static final long SEED = 14025256;
    static final long LIMIT = 2000000000000000L;

    public static String solve() {
        List<Integer> wDigits = new ArrayList<>();
        List<Long> prefixSum = new ArrayList<>();
        Map<Long, Integer> seen = new HashMap<>();

        long s = SEED;
        int index = 0;

        while (!seen.containsKey(s)) {
            seen.put(s, index);
            String sStr = String.valueOf(s);
            for (int i = 0; i < sStr.length(); i++) {
                wDigits.add(sStr.charAt(i) - '0');
            }
            s = (s * s) % MODULUS;
            index++;
        }

        prefixSum.add(0L);
        long currentSum = 0;
        for (int d : wDigits) {
            currentSum += d;
            prefixSum.add(currentSum);
        }

        long S = currentSum;
        int L = wDigits.size();

        byte[] R = new byte[(int) S + 1];
        long[] PModS = new long[L];

        for (int i = 0; i < L; ++i) {
            PModS[i] = prefixSum.get(i) % S;
        }

        for (int i = 0; i <= L; ++i) {
            R[(int) (prefixSum.get(i) % S)] = 1;
        }

        long[] pVals = new long[(int) S + 1];
        int numThreads = Math.max(1, Runtime.getRuntime().availableProcessors());
        long chunkSize = S / numThreads;

        ExecutorService executor = Executors.newFixedThreadPool(numThreads);
        List<Future<?>> futures = new ArrayList<>();

        for (int i = 0; i < numThreads; ++i) {
            final long startK = 1 + i * chunkSize;
            final long endK = (i == numThreads - 1) ? S : (startK + chunkSize - 1);

            if (startK <= endK) {
                futures.add(executor.submit(() -> {
                    for (long k = startK; k <= endK; ++k) {
                        long z = 1;
                        while (true) {
                            long startResidue = PModS[(int) ((z - 1) % L)];
                            long target = startResidue + k;
                            if (target >= S) {
                                target -= S;
                            }

                            if (R[(int) target] == 1) {
                                pVals[(int) k] = z;
                                break;
                            }
                            z++;
                        }
                    }
                }));
            }
        }

        for (Future<?> f : futures) {
            try {
                f.get();
            } catch (Exception e) {
            }
        }
        executor.shutdown();

        long numPeriods = LIMIT / S;
        long remainder = LIMIT % S;

        long sumPKOnePeriod = 0;
        long sumPKRemainder = 0;

        for (int k = 1; k <= S; ++k) {
            sumPKOnePeriod += pVals[k];
        }

        for (int k = 1; k <= remainder; ++k) {
            sumPKRemainder += pVals[k];
        }

        long totalAns = numPeriods * sumPKOnePeriod + sumPKRemainder;
        return String.valueOf(totalAns);
    }

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