Problem 375: Minimum of Subsequences

Project Euler Problem 375 Solution

EulerSolve provides an optimized solution for Project Euler Problem 375, Minimum of Subsequences, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The sequence in the repository is generated by $$S_{n+1}\equiv S_n^2 \pmod{50515093},\qquad S_0=290797,$$ and the working array starts at $S_1=629527$. For the first $N=2\cdot 10^9$ terms we must evaluate $$M(N)=\sum_{1\le l\le r\le N}\min(S_l,S_{l+1},\ldots,S_r).$$ A direct scan over all subarrays is hopeless, and even the standard linear-time subarray-minimum technique cannot be run on $2\cdot 10^9$ explicit values. The solution used in the C++, Python, and Java files therefore compresses the sequence to one period and then counts how often each position of that period owns a minimum. Mathematical Approach Step 1: Reduce the recurrence to one period The map $x\mapsto x^2 \bmod 50515093$ is deterministic on a finite state space, so once a value repeats, the future repeats as well. The implementation starts from $S_1$, keeps generating values until $S_1$ appears again, and obtains one full cycle $$a_1,a_2,\ldots,a_p,\qquad p=6308948.$$ Hence the infinite sequence relevant to the problem is the periodic extension $$b_x=a_{((x-1)\bmod p)+1}\qquad (x\ge 1).$$ Write the target length as $$N=qp+r,\qquad 0\le r<p.$$ Then period index $i$ appears $$\operatorname{occ}_i=q+\mathbf{1}_{i\le r}$$ times inside the first $N$ terms....

Detailed mathematical approach

Problem Summary

A sequence is produced by a small pseudo-random generator,

$$\begin{aligned}S_0&=290797,\\ S_{n+1}&=S_n^2\bmod 50515093 .\end{aligned}$$

For $i\le j$ let $A(i,j)$ be the smallest of the numbers $S_i,S_{i+1},\dots,S_j$, and let

$$M(N)=\sum_{1\le i\le j\le N}A(i,j).$$

The statement gives $M(10)=432256955$ and $M(10\,000)=3264567774119$ as checkpoints and asks for $M(N)$ at $N=2\,000\,000\,000$. ^[1] The seed $S_0$ never takes part in a comparison; the terms that occur are $S_1,\dots,S_N$, which begin

$$629527,\ 13339144,\ 15552512,\ 17939732,$$

$$3034546,\ 22608053,\ 23794117,\ \dots$$

The number of windows,

$$\frac{N(N+1)}2=2\,000\,000\,001\,000\,000\,000,$$

is far too large to enumerate. A linear scan is thinkable, since the terms can be generated as they are consumed, but it processes $2\cdot10^9$ of them, about $317$ times the work of the method below.

The whole solution in one sentence: the sequence repeats from $S_1$ with period $p=6308948$, so we study one period only, work out for each of its $p$ entries how far a window may reach left and right while that entry stays its smallest, and add the contributions of its $317$ or $318$ copies with a closed formula.

Mathematical Approach

1. Every window has exactly one owner

Write $b_x=S_x$ for $x\ge1$. Call the owner of a window $[l,r]$ the largest position $x\in[l,r]$ at which the minimum $A(l,r)$ is attained, that is its rightmost smallest entry. Every window has exactly one owner, so the minima can be summed position by position.

Position $x$ is the rightmost minimum of $[l,r]$ precisely when no entry of the window strictly before $x$ is smaller than $b_x$ and no entry after $x$ is smaller than or equal to $b_x$. With

$$P(x)=\max\{y\lt x:b_y\lt b_x\},$$

$$Q(x)=\min\{y\gt x:b_y\le b_x\},$$

where $P(x)=0$ and $Q(x)=N+1$ when no such $y$ exists in $1,\dots,N$, that condition reads $P(x)\lt l\le x\le r\lt Q(x)$. Hence $x$ owns $L_xR_x$ windows, where

$$L_x=x-P(x),\qquad R_x=Q(x)-x,$$

and, counting each window once at its owner,

$$M(N)=\sum_{x=1}^{N}b_xL_xR_x .$$

The two comparisons must differ in strictness, or a window whose minimum occurs twice would be counted twice or not at all. For $(2,2)$, position $2$ owns $[2,2]$ and $[1,2]$, position $1$ only $[1,1]$, and the three minima sum to $2+2\cdot2=6$.

2. A first check: the ten pairs behind $M(10)$

The last three of the ten terms are $S_8=25904962$, $S_9=10474246$ and $S_{10}=38959070$. The pairs $(L_x,R_x)$ for $x=1,\dots,10$ are, in order,

$$\begin{aligned}&(1,10),(1,3),(1,2),(1,1),(4,6),\\ &(1,3),(1,2),(1,1),(4,2),(1,1).\end{aligned}$$

The first term $629527$ is the smallest of the ten, so $P(1)=0$, $Q(1)=11$ and it owns all ten windows starting at position $1$. At $x=5$ the value $3034546$ is smaller than the three entries just before it but larger than $S_1$, so $P(5)=1$; nothing to its right is smaller, so $Q(5)=11$ and it owns $4\cdot6=24$ windows. The ten products $L_xR_x$ add up to $55=10\cdot11/2$, and

$$\sum_{x=1}^{10}b_xL_xR_x=432256955,$$

the first checkpoint of the statement.

3. The sequence repeats from $S_1$

Iterating a fixed map on a finite set must repeat, so the orbit is a tail running into a cycle. ^[5] Here the tail is one step long, provably so. The modulus factors as

$$50515093=5807\cdot 8699$$

with both factors prime and congruent to $3\pmod 4$, so the recurrence is the state iteration of the Blum–Blum–Shub generator. ^[2] Now $S_n\equiv S_1^{2^{n-1}}$ modulo $50515093$, so everything depends on the multiplicative order of $S_1=629527$ — the least $k\ge1$ with $S_1^k\equiv1$, which exists because $S_1$ is coprime to the modulus. ^[4] The group orders are $5806=2\cdot2903$ and $8698=2\cdot4349$ with $2903$, $4349$ prime, so the order of the square $S_1=S_0^2$ is odd and divides $2903$, respectively $4349$; neither is $1$, since $S_1\equiv2371\pmod{5807}$ and $S_1\equiv3199\pmod{8699}$. A power is $1$ modulo the product exactly when it is $1$ modulo both factors, so

$$d=\operatorname{ord}_{50515093}(S_1)=2903\cdot4349$$

is odd. Consequently $2$ is invertible modulo $d$, and

$$\begin{aligned}S_{1+T}=S_1&\iff S_1^{2^T-1}\equiv1\\ &\iff 2^T\equiv1\pmod d.\end{aligned}$$

The smallest such $T$ is $\operatorname{ord}_d(2)$, so the sequence is purely periodic from $n=1$ with period

$$\begin{aligned}p=\operatorname{ord}_d(2)&=\mathrm{lcm}(1451,4348)\\&=6308948,\end{aligned}$$

because $\operatorname{ord}_{2903}(2)=1451$ and $\operatorname{ord}_{4349}(2)=4348$; both follow from a few modular exponentiations, the only candidates being the divisors of $2902=2\cdot1451$ and $4348=4\cdot1087$. Every value on the cycle likewise has odd order, dividing the Carmichael value $\lambda(50515093)=2\cdot2903\cdot4349$ ^[3] and hence dividing $d$, so its period divides $\mathrm{lcm}(2902,4348)=p$ ^[4]: our start attains the longest period available here. The seed itself has the even order $8698$ modulo $8699$, while a value $x$ on the cycle satisfies $x^{2^p-1}\equiv1$ with $2^p-1$ odd and so has odd order; hence $S_0$ is off the cycle and the tail is exactly one step, as Floyd's tortoise and hare confirms. ^[5]

Write $a_i=S_i$ for $1\le i\le p$ and extend periodically, $b_x=a_{((x-1)\bmod p)+1}$. The entries $a_1,\dots,a_p$ are pairwise distinct: if $a_i=a_j$ with $i\lt j$, then applying the map $p+1-j$ times to both gives $S_{i+p+1-j}=S_{p+1}=S_1$ with $2\le i+p+1-j\le p$, so $S_1$ would return in fewer than $p$ steps, contradicting the minimality of $p$. The smallest entry is $a_{2633997}=3$, the largest $50515083$.

4. The two spans of a period index

For a period index $i$ define, in the two-sided periodic extension,

$$\rho_i=\min\{k\ge1:b_{i+k}\le a_i\},$$

$$\ell_i=\min\{k\ge1:b_{i-k}\lt a_i\},$$

the distance forward to the next entry at most $a_i$ and the distance back to the previous strictly smaller entry; the asymmetry is the tie-break of Section 1. Two bounds keep the work finite.

$\rho_i\le p$, because $b_{i+p}=a_i$. Equality holds exactly when $a_i$ is the smallest value of the period, and as the values are distinct that is a single index: $\rho_{2633997}=p$, while the next largest span $p-1$ belongs to the second smallest value $9$.
If $a_i$ is not the smallest value, a strictly smaller entry lies among the $p-1$ positions $i-1,\dots,i-p+1$, so $\ell_i\le p-1$; the largest that really occurs is $\ell_{2049738}=5724688$, at the index whose value is $22$. For the smallest value no strictly smaller entry exists at all, and the programs record that case as $\ell_i=0$: a window may then reach back to position $1$.

5. Both spans in linear time

Finding for every position the nearest earlier smaller entry is the classical all-nearest-smaller-values problem, solved by one stack sweep: pop while the stacked value is at least the current one, read the surviving top as the answer, then push the current position. ^[6] Popping is safe: a new position whose value is no larger makes an older stacked entry useless, being nearer to every later position and not larger, so the survivors increase in value downwards and the top is the answer. Each position is pushed once and popped at most once, so a sweep costs $O(p)$. The programs run two such sweeps over two consecutive copies of the period, positions $1,\dots,2p$: the backward one, popping while the stacked value is larger than the current, records $\rho_i$ for $i\le p$; the forward one, popping while it is larger than or equal, records $\ell_i$ from the second copy, where a full period of predecessors is available. It writes $\ell_i=0$ when the stack is empty or its top lies further back than one period; by the bounds above the second case cannot occur here, and the first occurs exactly once, at the smallest entry of the period. Two copies suffice precisely because both spans are at most $p$.

The backward sweep never finds its stack empty while $i\le p$: after the pops the top is the nearest position to the right with value at most $a_i$, and the twin position $i+p\le2p$ is such a position — which is why the C++ may read the top without a test.

6. From one period to $N$ terms

Write $N=qp+r$ with $0\le r\lt p$. For the target, $q=317$ and $r=63484$, since $317\cdot6308948+63484=2\cdot10^9$. Period index $i$ occurs

$$m_i=q+[\,i\le r\,]$$

times among the first $N$ terms, the bracket being $1$ for $i\le r$ and $0$ otherwise: $63484$ indices occur $318$ times and the other $6245464$ occur $317$ times. An index with $m_i=0$, possible only for $N\lt p$, contributes nothing. The $t$-th copy sits at $x_{i,t}=(t-1)p+i$, and every copy has the spans of Section 4, except where the ends of the prefix cut them off.

On the right only the last copy can be cut. Its tail is $\tau_i=N-x_{i,m_i}+1$, which is $r-i+1$ for $i\le r$ and $p+r-i+1$ otherwise, so its right count is $c_i=\min(\rho_i,\tau_i)$. An earlier copy lies at least $p$ positions further left, so it has at least $p+\tau_i\ge p+1$ admissible right endpoints, more than $\rho_i$, and is not cut.
On the left only the first copy can be cut, and only if $\ell_i\gt0$: its left count is $f_i=\min(\ell_i,i)$, while a later copy starts at $x_{i,t}\ge p+i\gt p-1\ge\ell_i$. If $\ell_i=0$, every copy may reach back to position $1$, so its left count is the position itself, $x_{i,t}$.

7. The contribution of one period index

Multiplying the counts of Section 6 by the value gives the contribution $C_i$ of all copies of index $i$. For $\ell_i\gt0$ and $m_i\ge2$,

$$C_i=a_i\bigl(f_i\rho_i+(m_i-2)\ell_i\rho_i+\ell_ic_i\bigr),$$

the three summands being the first copy, the $m_i-2$ copies in the middle that are cut on neither side, and the last copy; for $m_i=1$ the single copy is cut on both sides and $C_i=a_if_ic_i$. For the one index with $\ell_i=0$ the left count grows with the position, so the copies before the last form an arithmetic progression,

$$\begin{aligned}\Sigma_i&=\sum_{t=1}^{m_i-1}x_{i,t}\\ &=(m_i-1)i+\frac{p(m_i-2)(m_i-1)}{2},\end{aligned}$$

$$C_i=a_i\bigl(\rho_i\Sigma_i+((m_i-1)p+i)\,c_i\bigr),$$

and $C_i=a_i\,i\,c_i$ when $m_i=1$. Finally

$$M(N)=\sum_{i=1}^{p}C_i .$$

The huge $N$ survives only inside $q$ and $r$. For $N=10$ and $N=10\,000$, $q=0$, so every index above $r$ is skipped and the same formula reproduces both checkpoints.

8. A toy period worked out completely

Take $p=3$ with the period $a=(3,1,2)$ and $N=7$, so the sequence is $3,1,2,3,1,2,3$. The spans are $\rho=(1,3,2)$ and $\ell=(1,0,1)$: one step after index $1$ comes an entry at most $3$; index $2$ carries the smallest value, so $\rho_2=p=3$ and $\ell_2=0$; after index $3$ the next entry at most $2$ is two steps away. With $N=2\cdot3+1$ we get $q=2$, $r=1$, $m=(3,2,2)$, $\tau=(1,3,2)$, $c=(1,3,2)$ and $f_1=f_3=1$. The contributions are

$$C_1=3\,(1\cdot1+1\cdot1\cdot1+1\cdot1)=9,$$

$$C_2=1\,(3\cdot2+(1\cdot3+2)\cdot3)=21,$$

$$C_3=2\,(1\cdot2+1\cdot2)=8,$$

with $\Sigma_2=1\cdot2+3\cdot0\cdot1/2=2$ in the middle line, and with no middle summand in $C_3$ because $m_3=2$. So $M(7)=9+21+8=38$. Adding the minima directly, the windows of length $1,2,\dots,7$ contribute $15,8,5,4,3,2,1$, and $15+8+5+4+3+2+1=38$ as well.

No span is cut there. The period $a=(2,1,3)$ with $N=3$ shows both cuts: $\ell=(2,0,1)$, $\rho=(1,3,1)$, every $m_i=1$, $\tau=(3,2,1)$ and $c=(1,2,1)$. The start of the prefix cuts the left count of index $1$ from $2$ to $f_1=1$, its end cuts the right count of index $2$ from $3$ to $2$, and the contributions are $2\cdot1\cdot1=2$, $1\cdot2\cdot2=4$ and $3\cdot1\cdot1=3$. The sum $9$ matches the direct sums $6$, $2$, $1$ over the windows of length $1$, $2$, $3$.

How the Code Works

build_period_sequence() implements Section 3: it squares the seed once to get $S_1=629527$, keeps it as first, then appends values until the next equals first again. The vector reserves $6\,500\,000$ entries, more than the true length $p=6308948$. It terminates only because $S_1$ lies on the cycle; from a value off the cycle it would run forever.

build_spans() implements Sections 4 and 5, with value_at(pos) returning a[(pos - 1) % p], the periodic extension. Its backward loop stores d_right[pos - 1] = stack.back() - pos for $\mathrm{pos}\le p$; its forward loop writes d_left[idx - 1] with idx = pos - p for $\mathrm{pos}\gt p$, using $0$ when the stack is empty — the case of the smallest entry — or when its top lies more than one period back, which cannot happen here. The two strictness conventions are exactly the $\lt$ and $\le$ of Section 1, so spans and ownership rule agree.

sum_subarray_min_prefix(a, d_left, d_right, n) is Sections 6 and 7: it computes q = n / p and r = n % p, loops over idx from $1$ to $p$, skips an index whose occ is zero (possible only for $n\lt p$), forms rem_last and r_last = min(dr, rem_last), and takes one of four branches: dl == 0 or not, each with occ == 1 or occ >= 2. Each factor is cast to unsigned __int128, a GCC extension ^[7], before multiplication, and to_string_u128 turns the total into a decimal string that main prints. main runs run_checkpoints unless --skip-checkpoints is given, rejects any other argument and returns a non-zero status if a checkpoint fails. ^[9] The period and the spans are kept in $32$-bit words, but the generator state is $64$-bit, because a residue near $5\cdot10^7$ squares to about $2.6\cdot10^{15}$.

Python puts everything in solve(), with a nested val_at and plain lists; its integers have unlimited precision ^[8]. Java's solve() calls buildPeriodSequence(), buildSpans() and sumSubarrayMinPrefix(); it holds the period in a boxed List<Integer>, hence its larger memory use, and uses a preallocated int[2 * p + 8] with a top counter instead of a growable stack. Its $0$ fallback for top == 0 is unreachable in the backward loop by Section 5, while in the forward loop it is exactly the $\ell_i=0$ case. Neither port runs the checkpoints. Java's accumulator is a single long, and as Java fixes long at $64$ signed bits and never signals integer overflow ^[8], both its products and its total must be shown to fit.

Complexity and Verification

Building the period costs one initial modular squaring and $p$ more in the loop; each sweep visits $2p$ positions, pushing and popping each at most once; the final loop performs $p$ iterations of constant work. The method uses $O(p)$ arithmetic operations and $O(p)$ memory with $p=6308948$, and $N$ enters only through $q$ and $r$; with fixed-width arithmetic that is $O(p)$ time. The three arrays of $p$ $32$-bit words hold $76$ MB, the measured peak of the C++ run, since reserved but untouched capacity never becomes resident and the stack stays under $40$ entries. Java's boxed period and Python's integer objects push their peaks to about $370$ MB. In single runs on one machine the C++ took about $0.6$ s, the Java $0.5$ s and the Python $7$ s, all printing the same value.

run_checkpoints verifies the period length $6308948$ and the statement values $M(10)=432256955$ and $M(10\,000)=3264567774119$, both through the full periodic machinery. For this article the following was verified independently.

A separate program that generates $S_1,\dots,S_N$ straight from the recurrence, with no periodicity, and sums the minima by the sweeps of Section 1 agrees with the periodic formula at $N=10$, $10\,000$, $6308947$, $6308948$, $6308949$, $12617896$, $12617900$, $18926844$, $18928078$ and $20000000$ — on both sides of the first three period boundaries. The definition itself, a double loop over all pairs, agrees with both for every $N\le400$.
Keeping the computed spans and replacing only the value factor $a_i$ by $1$ turns the same bookkeeping into a count of windows, and it returns exactly $N(N+1)/2$ at $N=10$, $10\,000$, $p$, $2p+7$ and $2\cdot10^9$. This tests the occurrence numbers $m_i$, both cuts and the partition of the windows, with no value of the sequence involved.
Floyd's algorithm on the raw recurrence returns a tail of $1$ and a cycle of $6308948$, and a hash set confirms that the $p$ values are pairwise distinct. The spans obey Section 4: $1\le\rho_i\le p$ with equality once, every non-zero $\ell_i$ at most $p-1$, and exactly one $\ell_i=0$.
Every term is at least $3$, so $M(N)\ge 3N(N+1)/2$, two thirds of the signed $64$-bit range. Recomputing the sum with unlimited precision confirms that it, and every single product the code forms, stays below $2^{63}-1$, so the Java accumulator is sound here; but the margin is thin, and as $M$ grows like $N^2$ the same program would wrap silently at a moderately larger target. The C++ needs no such computation: every window minimum is at most $50515083$ and there are $N(N+1)/2$ windows, so $M(N)\lt2^{87}$, far inside the $128$-bit range.

Footnotes and References

The sources below supply the background facts; the ownership rule, the span bounds, the occurrence counts and the closed form for $C_i$ are derived in this article.

Project Euler 375 — Minimum of Subsequences. The official statement defines S₀ = 290797 and S_{n+1} = S_n² mod 50515093, sets A(i, j) to the minimum of S_i, …, S_j for i ≤ j and M(N) to the sum of A(i, j) over 1 ≤ i ≤ j ≤ N, and gives M(10) = 432256955 and M(10 000) = 3264567774119 before asking for M(2·10⁹).
Wikipedia — Blum Blum Shub. Defines the generator x_{n+1} = x_n² mod M with M = pq, both primes congruent to 3 (mod 4), and computes a term directly as x_i = x₀^(2^i mod λ(M)) mod M with λ(M) = lcm(p − 1, q − 1). It is the background for this recurrence, whose period Section 3 derives from multiplicative orders.
Wikipedia — Carmichael function. λ(n) is the smallest positive m with a^m ≡ 1 (mod n) for all a coprime to n, the exponent of the multiplicative group modulo n, and λ of coprime factors is the lcm of their values, so λ(50515093) = lcm(5806, 8698).
Wikipedia — Multiplicative order. For a coprime to n, ord_n(a) is the smallest positive k with a^k ≡ 1 (mod n), the order of a in the group of units modulo n, and it divides λ(n) and φ(n).
Wikipedia — Cycle detection. The orbit of a repeated map on a finite set looks like the letter rho, a tail leading into a cycle, and Floyd's tortoise-and-hare algorithm finds both lengths with two pointers of different speed, in constant memory.
Wikipedia — All nearest smaller values. The problem of finding, for every position, the nearest earlier position with a smaller value; its stack sweep pops entries greater than or equal to the current value, reads the remaining top as the answer, and is linear because each element is popped at most once.
GCC manual — 128-bit integers. Documents __int128 and unsigned __int128 as a GCC extension on targets with an integer mode wide enough for 128 bits. The C++ accumulator is the unsigned variant, whose range ends at 2¹²⁸ − 1, far above the bound in the verification section.
Java Language Specification — primitive types. Section 4.2.1 fixes long at −9223372036854775808 to 9223372036854775807, and 4.2.2 states that the integer operators do not indicate overflow or underflow in any way, so a Java sum that leaves that range is wrong silently. Python — numeric types documents the unlimited precision of Python integers.
C++ — Euler375.cpp. The immutable revision holding build_period_sequence, build_spans, sum_subarray_min_prefix and run_checkpoints; the implementation record for this article. The comparisons above are extra checks, not tests packaged in that file.

Problem 375 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr u64 kMod = 50515093ULL;
constexpr u64 kTargetN = 2000000000ULL;

std::vector<u32> build_period_sequence() {
    u64 s = 290797ULL;
    s = (s * s) % kMod;  // S_1
    const u32 first = static_cast<u32>(s);

    std::vector<u32> seq;
    seq.reserve(6500000);
    while (true) {
        seq.push_back(static_cast<u32>(s));
        s = (s * s) % kMod;
        if (static_cast<u32>(s) == first) {
            break;
        }
    }
    return seq;
}

void build_spans(const std::vector<u32>& a, std::vector<u32>& d_left, std::vector<u32>& d_right) {
    const u32 p = static_cast<u32>(a.size());
    d_left.assign(static_cast<std::size_t>(p), 0U);   // 0 means no strictly smaller in previous period
    d_right.assign(static_cast<std::size_t>(p), 0U);  // always positive

    auto value_at = [&](const u32 pos) -> u32 {
        return a[static_cast<std::size_t>((pos - 1U) % p)];
    };

    std::vector<u32> stack;
    stack.reserve(static_cast<std::size_t>(2U * p + 8U));

    // Next position with value <= current.
    for (u32 pos = 2U * p; pos >= 1U; --pos) {
        const u32 v = value_at(pos);
        while (!stack.empty() && value_at(stack.back()) > v) {
            stack.pop_back();
        }
        if (pos <= p) {
            d_right[static_cast<std::size_t>(pos - 1U)] = stack.back() - pos;
        }
        stack.push_back(pos);
        if (pos == 1U) {
            break;
        }
    }

    // Previous position with value < current, limited to at most one period back.
    stack.clear();
    for (u32 pos = 1U; pos <= 2U * p; ++pos) {
        const u32 v = value_at(pos);
        while (!stack.empty() && value_at(stack.back()) >= v) {
            stack.pop_back();
        }
        if (pos > p) {
            const u32 idx = pos - p;  // 1-based index in one period
            if (!stack.empty() && stack.back() >= idx) {
                d_left[static_cast<std::size_t>(idx - 1U)] = pos - stack.back();
            } else {
                d_left[static_cast<std::size_t>(idx - 1U)] = 0U;
            }
        }
        stack.push_back(pos);
    }
}

u128 sum_subarray_min_prefix(const std::vector<u32>& a,
                             const std::vector<u32>& d_left,
                             const std::vector<u32>& d_right,
                             const u64 n) {
    const u64 p = static_cast<u64>(a.size());
    const u64 q = n / p;
    const u64 r = n % p;

    u128 total = 0U;
    for (u64 idx = 1ULL; idx <= p; ++idx) {
        const u64 occ = q + ((idx <= r) ? 1ULL : 0ULL);
        if (occ == 0ULL) {
            continue;
        }

        const u64 v = static_cast<u64>(a[static_cast<std::size_t>(idx - 1ULL)]);
        const u64 dr = static_cast<u64>(d_right[static_cast<std::size_t>(idx - 1ULL)]);
        const u64 dl = static_cast<u64>(d_left[static_cast<std::size_t>(idx - 1ULL)]);

        const u64 rem_last = (idx <= r) ? (r - idx + 1ULL) : (p + r - idx + 1ULL);
        const u64 r_last = std::min(dr, rem_last);

        if (dl == 0ULL) {
            // No strictly smaller element exists to the left in any previous period.
            if (occ == 1ULL) {
                total += static_cast<u128>(v) * static_cast<u128>(idx) * static_cast<u128>(r_last);
            } else {
                const u64 n1 = occ - 1ULL;
                const u128 sum_l = static_cast<u128>(n1) * static_cast<u128>(idx) +
                                   static_cast<u128>(p) * static_cast<u128>(n1 - 1ULL) * static_cast<u128>(n1) / 2U;
                total += static_cast<u128>(v) * static_cast<u128>(dr) * sum_l;

                const u64 l_last = (occ - 1ULL) * p + idx;
                total += static_cast<u128>(v) * static_cast<u128>(l_last) * static_cast<u128>(r_last);
            }
        } else {
            if (occ == 1ULL) {
                const u64 l0 = std::min(dl, idx);
                total += static_cast<u128>(v) * static_cast<u128>(l0) * static_cast<u128>(r_last);
            } else {
                const u64 l0 = std::min(dl, idx);
                total += static_cast<u128>(v) * static_cast<u128>(l0) * static_cast<u128>(dr);
                if (occ > 2ULL) {
                    total += static_cast<u128>(v) * static_cast<u128>(dl) * static_cast<u128>(dr) *
                             static_cast<u128>(occ - 2ULL);
                }
                total += static_cast<u128>(v) * static_cast<u128>(dl) * static_cast<u128>(r_last);
            }
        }
    }

    return total;
}

std::string to_string_u128(u128 value) {
    if (value == 0U) {
        return "0";
    }
    std::string out;
    while (value > 0U) {
        const unsigned digit = static_cast<unsigned>(value % 10U);
        out.push_back(static_cast<char>('0' + digit));
        value /= 10U;
    }
    std::reverse(out.begin(), out.end());
    return out;
}

bool run_checkpoints(const std::vector<u32>& a, const std::vector<u32>& d_left, const std::vector<u32>& d_right) {
    if (a.size() != 6308948U) {
        std::cerr << "Checkpoint failed: period length\n";
        return false;
    }

    if (sum_subarray_min_prefix(a, d_left, d_right, 10ULL) != 432256955ULL) {
        std::cerr << "Checkpoint failed: M(10)\n";
        return false;
    }

    if (sum_subarray_min_prefix(a, d_left, d_right, 10000ULL) != 3264567774119ULL) {
        std::cerr << "Checkpoint failed: M(10000)\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    const std::vector<u32> period = build_period_sequence();
    std::vector<u32> d_left;
    std::vector<u32> d_right;
    build_spans(period, d_left, d_right);

    if (!skip_checkpoints && !run_checkpoints(period, d_left, d_right)) {
        return 2;
    }

    const u128 answer = sum_subarray_min_prefix(period, d_left, d_right, kTargetN);
    std::cout << to_string_u128(answer) << '\n';
    return 0;
}

Python

def solve():
    MOD = 50515093
    TARGET_N = 2_000_000_000

    # Build period sequence
    s = 290797
    s = s * s % MOD  # S_1
    first = s
    seq = []
    while True:
        seq.append(s)
        s = s * s % MOD
        if s == first:
            break

    p = len(seq)

    # Build spans (d_left, d_right) for the periodic sequence
    d_left = [0] * p
    d_right = [0] * p

    def val_at(pos):
        return seq[(pos - 1) % p]

    # d_right: next position with value <= current (within 2 periods)
    stack = []
    for pos in range(2 * p, 0, -1):
        v = val_at(pos)
        while stack and val_at(stack[-1]) > v:
            stack.pop()
        if pos <= p:
            d_right[pos - 1] = stack[-1] - pos
        stack.append(pos)

    # d_left: prev position with value < current (within 2 periods)
    stack = []
    for pos in range(1, 2 * p + 1):
        v = val_at(pos)
        while stack and val_at(stack[-1]) >= v:
            stack.pop()
        if pos > p:
            idx = pos - p
            if stack and stack[-1] >= idx:
                d_left[idx - 1] = pos - stack[-1]
        stack.append(pos)

    # Compute sum of subarray minimums for prefix of length n
    n = TARGET_N
    q_full = n // p
    r = n % p

    total = 0
    for idx in range(1, p + 1):
        occ = q_full + (1 if idx <= r else 0)
        if occ == 0:
            continue
        v = seq[idx - 1]
        dr = d_right[idx - 1]
        dl = d_left[idx - 1]

        rem_last = (r - idx + 1) if idx <= r else (p + r - idx + 1)
        r_last = min(dr, rem_last)

        if dl == 0:
            if occ == 1:
                total += v * idx * r_last
            else:
                n1 = occ - 1
                sum_l = n1 * idx + p * (n1 - 1) * n1 // 2
                total += v * dr * sum_l
                l_last = (occ - 1) * p + idx
                total += v * l_last * r_last
        else:
            if occ == 1:
                l0 = min(dl, idx)
                total += v * l0 * r_last
            else:
                l0 = min(dl, idx)
                total += v * l0 * dr
                if occ > 2:
                    total += v * dl * dr * (occ - 2)
                total += v * dl * r_last

    return str(total)

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

Java

import java.util.ArrayList;
import java.util.List;

public class Euler375 {
    static final long kMod = 50515093L;
    static final long kTargetN = 2000000000L;

    static List<Integer> buildPeriodSequence() {
        long s = 290797L;
        s = (s * s) % kMod;
        int first = (int) s;

        List<Integer> seq = new ArrayList<>(6500000);
        while (true) {
            seq.add((int) s);
            s = (s * s) % kMod;
            if ((int) s == first) {
                break;
            }
        }
        return seq;
    }

    static int valueAt(List<Integer> a, int pos) {
        int p = a.size();
        return a.get((pos - 1) % p);
    }

    static void buildSpans(List<Integer> a, int[] dLeft, int[] dRight) {
        int p = a.size();
        int[] stack = new int[2 * p + 8];
        int top = 0;

        for (int pos = 2 * p; pos >= 1; pos--) {
            int v = valueAt(a, pos);
            while (top > 0 && valueAt(a, stack[top - 1]) > v) {
                top--;
            }
            if (pos <= p) {
                dRight[pos - 1] = top > 0 ? (stack[top - 1] - pos) : 0;
            }
            stack[top++] = pos;
            if (pos == 1)
                break;
        }

        top = 0;
        for (int pos = 1; pos <= 2 * p; pos++) {
            int v = valueAt(a, pos);
            while (top > 0 && valueAt(a, stack[top - 1]) >= v) {
                top--;
            }
            if (pos > p) {
                int idx = pos - p;
                if (top > 0 && stack[top - 1] >= idx) {
                    dLeft[idx - 1] = pos - stack[top - 1];
                } else {
                    dLeft[idx - 1] = 0;
                }
            }
            stack[top++] = pos;
        }
    }

    static long sumSubarrayMinPrefix(List<Integer> a, int[] dLeft, int[] dRight, long n) {
        long p = a.size();
        long q = n / p;
        long r = n % p;

        long total = 0;
        for (long idx = 1; idx <= p; idx++) {
            long occ = q + ((idx <= r) ? 1 : 0);
            if (occ == 0)
                continue;

            long v = a.get((int) (idx - 1));
            long dr = dRight[(int) (idx - 1)];
            long dl = dLeft[(int) (idx - 1)];

            long remLast = (idx <= r) ? (r - idx + 1) : (p + r - idx + 1);
            long rLast = Math.min(dr, remLast);

            if (dl == 0) {
                if (occ == 1) {
                    total += v * idx * rLast;
                } else {
                    long n1 = occ - 1;
                    long sumL = n1 * idx + p * (n1 - 1) * n1 / 2;
                    total += v * dr * sumL;

                    long lLast = (occ - 1) * p + idx;
                    total += v * lLast * rLast;
                }
            } else {
                if (occ == 1) {
                    long l0 = Math.min(dl, idx);
                    total += v * l0 * rLast;
                } else {
                    long l0 = Math.min(dl, idx);
                    total += v * l0 * dr;
                    if (occ > 2) {
                        total += v * dl * dr * (occ - 2);
                    }
                    total += v * dl * rLast;
                }
            }
        }
        return total;
    }

    static String solve() {
        List<Integer> period = buildPeriodSequence();
        int p = period.size();
        int[] dLeft = new int[p];
        int[] dRight = new int[p];
        buildSpans(period, dLeft, dRight);

        long answer = sumSubarrayMinPrefix(period, dLeft, dRight, kTargetN);
        return Long.toString(answer); // Works for BigInteger since it fits in Long
    }

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