Problem 375: Minimum of Subsequences
View on Project EulerProject 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
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.
Step 2: Give every subarray minimum a unique owner
For a fixed absolute position \(x\), the classical contribution method counts how many subarrays have \(b_x\) as their representative minimum. Equal values require a tie-breaking rule. The code uses a strict comparison on the left and a non-strict comparison on the right, so ties are assigned to the rightmost occurrence of the minimum.
For a period index \(i\), with cyclic indexing modulo \(p\), define
$$d_R(i)=\min\{k\ge 1: a_{i+k}\le a_i\},$$
and define \(d_L(i)\) as the smallest \(k\ge 1\) with \(a_{i-k}<a_i\) if such a value exists within one full period to the left; otherwise set \(d_L(i)=0\).
The quantity \(d_R(i)\) is always at most \(p\), because the same value reappears one full period later. For \(d_L(i)\), looking back only one period is enough: if no strictly smaller value occurs there, periodicity implies that no strictly smaller value exists anywhere earlier.
Step 3: Build the span arrays in linear time
The implementations scan two consecutive copies of the period. A right-to-left monotone stack produces the next position with value \(\le a_i\), which gives \(d_R(i)\). A left-to-right monotone stack with the opposite strictness produces the previous position with value \(<a_i\), which gives \(d_L(i)\).
This is the usual subarray-minimum ownership trick, but adapted to a cyclic array. Two copies are sufficient because no relevant comparison can jump farther than one period.
Step 4: Contribution of one periodic index
The \(t\)-th occurrence of period index \(i\) in the prefix of length \(N\) sits at
$$x_{i,t}=(t-1)p+i,\qquad 1\le t\le \operatorname{occ}_i.$$
If \(L_{i,t}\) is the number of admissible left endpoints and \(R_{i,t}\) is the number of admissible right endpoints for subarrays owned by that occurrence, then its contribution is
$$a_i\,L_{i,t}\,R_{i,t}.$$
Every occurrence except the last one has the full right span \(d_R(i)\). Only the last occurrence may be cut by the end of the prefix. Its remaining tail length is
$$\operatorname{tail}_i=\begin{cases} r-i+1, & i\le r,\\ p+r-i+1, & i>r, \end{cases}$$
so the actual right factor on the last copy is
$$R_i^{\text{last}}=\min(d_R(i),\operatorname{tail}_i).$$
Step 5: Closed forms matching the code
If \(d_L(i)=0\), then no strictly smaller value exists anywhere to the left, so the left span is simply the absolute position itself:
$$L_{i,t}=x_{i,t}=(t-1)p+i.$$
When \(\operatorname{occ}_i=1\), this gives
$$C_i=a_i\cdot i\cdot R_i^{\text{last}}.$$
When \(\operatorname{occ}_i\ge 2\), the first \(\operatorname{occ}_i-1\) copies use the full right span \(d_R(i)\), so the code sums an arithmetic progression:
$$\sum_{t=1}^{\operatorname{occ}_i-1}\big((t-1)p+i\big) =(\operatorname{occ}_i-1)i+p\frac{(\operatorname{occ}_i-2)(\operatorname{occ}_i-1)}{2}.$$
Therefore
$$C_i=a_i\left[ d_R(i)\left((\operatorname{occ}_i-1)i+p\frac{(\operatorname{occ}_i-2)(\operatorname{occ}_i-1)}{2}\right) +\big((\operatorname{occ}_i-1)p+i\big)R_i^{\text{last}} \right].$$
If \(d_L(i)>0\), the first copy may be truncated on the left, but every later copy has the fixed left span \(d_L(i)\). Let
$$L_i^{\text{first}}=\min(d_L(i),i).$$
Then, for \(\operatorname{occ}_i=1\),
$$C_i=a_i\,L_i^{\text{first}}\,R_i^{\text{last}},$$
and for \(\operatorname{occ}_i\ge 2\),
$$C_i=a_i\left[L_i^{\text{first}}d_R(i)+(\operatorname{occ}_i-2)d_L(i)d_R(i)+d_L(i)R_i^{\text{last}}\right].$$
The required value is the sum over one period:
$$\boxed{M(N)=\sum_{i=1}^{p} C_i.}$$
How the Code Works
The C++ solution separates the work into build_period_sequence(), build_spans(), and sum_subarray_min_prefix(). The Python and Java versions mirror the same formulas and case split. The C++ file also checks the key intermediate facts used by the derivation:
$$p=6308948,\qquad M(10)=432256955,\qquad M(10000)=3264567774119.$$
Those checkpoints are important because they validate both the period construction and the strict/non-strict ownership convention before the final \(N=2\cdot 10^9\) computation is performed.
Complexity Analysis
Building one period takes \(O(p)\) time. Each monotone-stack pass is \(O(p)\), because every position is pushed and popped at most once. The final aggregation over all \(p\) period indices is again \(O(p)\). Thus the full algorithm runs in \(O(p)\) time and \(O(p)\) memory, with \(p=6308948\), while the huge input \(N\) only appears through the two integers \(q\) and \(r\).
Footnotes and References
- Problem page: https://projecteuler.net/problem=375
- Monotone stacks and subarray minimum ownership: cp-algorithms
- Cycle detection in deterministic sequences: Wikipedia — Cycle detection
- Modular arithmetic background: Wikipedia — Modular arithmetic
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());
}
}