Problem 693: Finite Sequence Generator
View on Project EulerProject Euler Problem 693 Solution
EulerSolve provides an optimized solution for Project Euler Problem 693, Finite Sequence Generator, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Start from integers \(x \ge 2\) and \(y\), and define $$a_0=y,\qquad z_0=x,\qquad a_{t+1}=a_t^2 \bmod z_t,\qquad z_{t+1}=z_t+1.$$ The sequence is counted until the first time \(a_t \le 1\). Since both \(0\) and \(1\) stay fixed under squaring, they are terminal values. Let \(\ell(x,y)\) be the number of terms from \(a_0\) through that first terminal term. Then $$g(x)=\max_y \ell(x,y),\qquad f(n)=\max_{2 \le x \le n} g(x).$$ The direct search over all starting values is far too expensive, so the implementations reorganize the computation around the distinct residues that can still be alive after each step. Mathematical Approach The key observation is that the process quickly merges many different starts into the same residue. Once we track those shared residues instead of individual starts, the problem becomes manageable. Step 1: Compress the Starting Values into First-Step Residues For fixed \(x\), the first update depends only on \(y \bmod x\), so it is enough to consider residues \(0 \le y < x\). Moreover, $$ (x-y)^2 \equiv y^2 \pmod{x}, $$ so the starts \(y\) and \(x-y\) merge immediately. For \(x>2\), the only starts worth keeping are \(2 \le y \le \lfloor x/2 \rfloor\), and only if the first residue is greater than \(1\)....
Detailed mathematical approach
Problem Summary
Start from integers \(x \ge 2\) and \(y\), and define
$$a_0=y,\qquad z_0=x,\qquad a_{t+1}=a_t^2 \bmod z_t,\qquad z_{t+1}=z_t+1.$$
The sequence is counted until the first time \(a_t \le 1\). Since both \(0\) and \(1\) stay fixed under squaring, they are terminal values. Let \(\ell(x,y)\) be the number of terms from \(a_0\) through that first terminal term. Then
$$g(x)=\max_y \ell(x,y),\qquad f(n)=\max_{2 \le x \le n} g(x).$$
The direct search over all starting values is far too expensive, so the implementations reorganize the computation around the distinct residues that can still be alive after each step.
Mathematical Approach
The key observation is that the process quickly merges many different starts into the same residue. Once we track those shared residues instead of individual starts, the problem becomes manageable.
Step 1: Compress the Starting Values into First-Step Residues
For fixed \(x\), the first update depends only on \(y \bmod x\), so it is enough to consider residues \(0 \le y < x\). Moreover,
$$ (x-y)^2 \equiv y^2 \pmod{x}, $$
so the starts \(y\) and \(x-y\) merge immediately. For \(x>2\), the only starts worth keeping are \(2 \le y \le \lfloor x/2 \rfloor\), and only if the first residue is greater than \(1\). Define the first active frontier by
$$A_1(x)=\left\{y^2 \bmod x : 2 \le y \le \left\lfloor \frac{x}{2} \right\rfloor,\ y^2 \bmod x > 1\right\}.$$
If \(A_1(x)=\varnothing\), then every nontrivial start reaches \(0\) or \(1\) after one update, so \(g(x)=2\). The only special case is \(x=2\), where even the initial value is already forced to be terminal and \(g(2)=1\).
Step 2: Propagate the Entire Frontier at Once
For \(t \ge 1\), let \(A_t(x)\) be the set of residues still alive after \(t\) updates. The recurrence is
$$A_{t+1}(x)=\left\{a^2 \bmod (x+t) : a \in A_t(x),\ a^2 \bmod (x+t) > 1\right\}.$$
Each set update simultaneously represents every starting value whose trajectory has merged into one of those residues. When \(A_t(x)\) becomes empty, every trajectory has terminated, so the total length is \(t+1\). This replaces many separate simulations by one deduplicated state update per time step.
Step 3: Exploit Singleton Tails
If the frontier shrinks to a single residue, say
$$A_t(x)=\{v\},$$
then all surviving starting values now share exactly the same future. From that moment onward, there is no benefit in keeping a set representation: the rest is just the single trajectory starting from value \(v\) with current modulus \(x+t\). In other words, after coalescence the problem stops being a many-start search and becomes one ordinary modular-squaring chain.
Step 4: Derive an Interval Upper Bound for the Outer Maximum
To compute \(f(n)\), the implementations do not scan every \(x\) from \(2\) to \(n\). Instead they use a rigorous upper bound on intervals. Fix \(L \le u \le R\). Take any start \(y\) for the modulus \(u\), and advance the sequence for exactly \(R-u\) updates. At that point the current modulus is \(R\), and the current value is some residue \(b<R\). Therefore the remaining lifetime is at most \(g(R)\), which gives
$$\ell(u,y) \le (R-u) + g(R).$$
Taking the maximum over all \(y\) yields
$$g(u) \le (R-u) + g(R) \le (R-L) + g(R).$$
So every interval \([L,R]\) has the valid bound
$$UB([L,R]) = g(R) + (R-L).$$
If this bound is no better than the best value already known, the whole interval can be discarded safely.
Worked Example: \(x=10\)
The first distinct nonterminal residues are
$$A_1(10)=\{2^2,3^2,4^2,5^2\}\bmod 10 \setminus \{0,1\}=\{4,5,6,9\}.$$
Now propagate the frontier:
$$A_2(10)=\{a^2 \bmod 11 : a \in A_1(10)\}\setminus\{0,1\}=\{3,4,5\},$$
$$A_3(10)=\{a^2 \bmod 12 : a \in A_2(10)\}\setminus\{0,1\}=\{4,9\},$$
$$A_4(10)=\{a^2 \bmod 13 : a \in A_3(10)\}\setminus\{0,1\}=\{3\}.$$
At this point the frontier is a singleton, so only one tail remains. Continuing from value \(3\) with moduli \(14,15,16,\dots\) gives
$$3 \to 9 \to 6 \to 4 \to 16 \to 4 \to 16 \to 16 \to \cdots \to 8 \to 0.$$
The terminal \(0\) appears right after the step modulo \(32\), so the sequence has \(24\) terms in total and therefore
$$g(10)=24.$$
This example shows all the main ideas: first-step symmetry, deduplication of equal residues, and the switch to a single tail once every surviving start has merged.
How the Code Works
The C++, Python, and Java implementations follow the same logic. First they build the initial frontier by scanning only half of the possible starting residues, using the symmetry \(y \leftrightarrow x-y\), and discarding first-step residues \(0\) and \(1\). To avoid recomputing each square from scratch during that initial scan, the compiled implementations update consecutive squares incrementally.
While the frontier contains more than one residue, the implementation advances the whole set one modulus at a time. Large frontiers use a dense seen-array for fast deduplication; small frontiers use a sparse hash-based container to avoid touching every residue class of the modulus. Once the frontier becomes empty, the current length is the answer for that \(x\). Once the frontier becomes a singleton, the remaining steps are computed by direct modular squaring of that one value.
For the outer function \(f(n)\), the implementation first samples a coarse grid of \(x\)-values, caches the computed \(g(x)\), and places the intervals between sampled points into a max-priority queue keyed by the bound \(UB([L,R])\). It then repeatedly bisects only those intervals whose upper bound can still beat the current best answer. The C++ and Java implementations evaluate the initial grid points in parallel, while the Python implementation uses the same search strategy sequentially.
Complexity Analysis
For a fixed \(x\), building the first frontier costs \(O(x)\) time, because the scan goes through \(2,3,\dots,\lfloor x/2 \rfloor\). After that, one sparse frontier update costs expected \(O(|A_t(x)|)\), while one dense update costs \(O(x+t)\) because the deduplication array spans the whole current modulus. Therefore the running time for \(g(x)\) is output-sensitive: it is governed by the cumulative frontier sizes until extinction or singleton collapse, rather than by the number of possible starts alone.
The memory usage for one evaluation of \(g(x)\) is \(O(x+t_{\max})\) in dense phases and \(O(|A_t(x)|)\) in sparse phases. For \(f(n)\), if the interval search evaluates \(m\) different values of \(x\), the total cost is
$$O\left(m + \sum_{i=1}^{m} T(x_i)\right),$$
where \(T(x_i)\) is the cost of computing \(g(x_i)\). In the worst case \(m\) can still be \(O(n)\), but the interval bound usually prunes most subintervals long before that.
Footnotes and References
- Problem page: https://projecteuler.net/problem=693
- Modular arithmetic: Wikipedia — Modular arithmetic
- Quadratic residues: Wikipedia — Quadratic residue
- Branch and bound: Wikipedia — Branch and bound
- Priority queues: Wikipedia — Priority queue
Problem 693 source code
C++
#include <algorithm>
#include <atomic>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <queue>
#include <thread>
#include <unordered_map>
#include <unordered_set>
#include <vector>
namespace {
using u32 = std::uint32_t;
using u64 = std::uint64_t;
u32 l_xy(const u32 x, const u32 y) {
u32 z = x;
u32 a = y;
u32 length = 1;
while (a > 1U) {
a = static_cast<u32>((static_cast<u64>(a) * static_cast<u64>(a)) % z);
++z;
++length;
}
return length;
}
std::vector<u32> initial_active_after_first_step(const u32 x) {
if (x <= 2U) {
return {};
}
std::vector<std::uint8_t> seen(x, 0U);
std::vector<u32> active;
active.reserve(x / 4U);
const u64 mod = static_cast<u64>(x);
const u64 limit = mod / 2ULL;
u64 y = 2ULL;
u64 sq = (y * y) % mod;
u64 delta = 2ULL * y + 1ULL;
while (y <= limit) {
if (sq > 1ULL && !seen[static_cast<std::size_t>(sq)]) {
seen[static_cast<std::size_t>(sq)] = 1U;
active.push_back(static_cast<u32>(sq));
}
sq += delta;
if (sq >= mod) {
sq -= mod;
if (sq >= mod) {
sq -= mod;
}
}
delta += 2ULL;
++y;
}
return active;
}
std::vector<u32> step_active_dense(const std::vector<u32>& active, const u32 mod) {
std::vector<std::uint8_t> seen(mod, 0U);
std::vector<u32> next;
next.reserve(active.size());
for (const u32 a : active) {
const u32 v = static_cast<u32>((static_cast<u64>(a) * static_cast<u64>(a)) % mod);
if (v > 1U && !seen[static_cast<std::size_t>(v)]) {
seen[static_cast<std::size_t>(v)] = 1U;
next.push_back(v);
}
}
return next;
}
std::vector<u32> step_active_sparse(const std::vector<u32>& active, const u32 mod) {
std::unordered_set<u32> next_set;
next_set.reserve(active.size() * 2ULL);
for (const u32 a : active) {
const u32 v = static_cast<u32>((static_cast<u64>(a) * static_cast<u64>(a)) % mod);
if (v > 1U) {
next_set.insert(v);
}
}
std::vector<u32> next;
next.reserve(next_set.size());
for (const u32 v : next_set) {
next.push_back(v);
}
return next;
}
u32 g_value(const u32 x, const std::size_t dense_threshold = 120000U) {
if (x < 2U) {
return 0U;
}
if (x == 2U) {
return 1U;
}
std::vector<u32> active = initial_active_after_first_step(x);
u32 length = 2U;
u32 mod = x + 1U;
while (active.size() > 1U) {
if (active.size() >= dense_threshold) {
active = step_active_dense(active, mod);
} else {
active = step_active_sparse(active, mod);
}
++length;
++mod;
if (active.empty()) {
return length;
}
}
if (active.empty()) {
return length;
}
u32 value = active[0];
while (value > 1U) {
value = static_cast<u32>((static_cast<u64>(value) * static_cast<u64>(value)) % mod);
++length;
++mod;
}
return length;
}
u32 snapped_grid_step(const u32 raw_step) {
if (raw_step <= 1U) {
return 1U;
}
const int power = static_cast<int>(std::floor(std::log10(static_cast<double>(raw_step))));
u32 base = 1U;
for (int i = 0; i < power; ++i) {
base *= 10U;
}
for (const u32 m : {1U, 2U, 5U, 10U}) {
if (m * base >= raw_step) {
return m * base;
}
}
return raw_step;
}
struct Interval {
u32 upper_bound;
u32 left;
u32 right;
u32 g_right;
};
struct IntervalOrder {
bool operator()(const Interval& a, const Interval& b) const {
return a.upper_bound < b.upper_bound;
}
};
u32 f_value(const u32 n, const u32 target_points = 16U) {
if (n < 2U) {
return 0U;
}
const u32 points_count = std::max(2U, target_points);
const u32 raw_step = std::max(1U, (n - 2U) / (points_count - 1U));
const u32 grid_step = snapped_grid_step(raw_step);
std::vector<u32> points;
points.reserve(points_count + 2U);
for (u32 x = 2U; x <= n; x += grid_step) {
points.push_back(x);
if (x > n - grid_step) {
break;
}
}
if (points.back() != n) {
points.push_back(n);
}
std::vector<u32> g_points(points.size(), 0U);
{
const unsigned hw = std::thread::hardware_concurrency();
const unsigned thread_count =
std::max(1U, std::min<unsigned>(hw == 0U ? 1U : hw, points.size()));
std::atomic<std::size_t> next_index(0U);
std::vector<std::thread> pool;
pool.reserve(thread_count);
for (unsigned t = 0U; t < thread_count; ++t) {
pool.emplace_back([&]() {
while (true) {
const std::size_t i = next_index.fetch_add(1U, std::memory_order_relaxed);
if (i >= points.size()) {
break;
}
g_points[i] = g_value(points[i]);
}
});
}
for (auto& th : pool) {
th.join();
}
}
std::unordered_map<u32, u32> cache;
cache.reserve(points.size() * 4ULL);
u32 best = 0U;
for (std::size_t i = 0; i < points.size(); ++i) {
cache[points[i]] = g_points[i];
best = std::max(best, g_points[i]);
}
auto G = [&](const u32 x) -> u32 {
const auto it = cache.find(x);
if (it != cache.end()) {
return it->second;
}
const u32 gx = g_value(x);
cache.emplace(x, gx);
return gx;
};
std::priority_queue<Interval, std::vector<Interval>, IntervalOrder> pq;
for (std::size_t i = 1; i < points.size(); ++i) {
const u32 left = points[i - 1U];
const u32 right = points[i];
const u32 g_right = g_points[i];
const u32 upper_bound = g_right + (right - left);
pq.push(Interval{upper_bound, left, right, g_right});
}
while (!pq.empty()) {
const Interval top = pq.top();
if (top.upper_bound <= best) {
break;
}
pq.pop();
if (top.right - top.left <= 1U) {
continue;
}
const u32 mid = (top.left + top.right) / 2U;
const u32 g_mid = G(mid);
best = std::max(best, g_mid);
if (mid - top.left > 1U) {
const u32 ub_left = g_mid + (mid - top.left);
if (ub_left > best) {
pq.push(Interval{ub_left, top.left, mid, g_mid});
}
}
if (top.right - mid > 1U) {
const u32 ub_right = top.g_right + (top.right - mid);
if (ub_right > best) {
pq.push(Interval{ub_right, mid, top.right, top.g_right});
}
}
}
return best;
}
} // namespace
int main() {
assert(l_xy(5U, 3U) == 29U);
assert(g_value(5U) == 29U);
assert(f_value(100U) == 145U);
assert(f_value(10'000U) == 8'824U);
std::cout << f_value(3'000'000U) << '\n';
return 0;
}
Python
def solve():
N = 3000000
def g_value(x):
if x < 2: return 0
if x == 2: return 1
# Initial active set after first step
mod = x; seen = bytearray(mod); active = []
limit = mod // 2
for y in range(2, limit+1):
sq = y*y % mod
if sq > 1 and not seen[sq]:
seen[sq] = 1; active.append(sq)
length = 2; mod = x + 1
while len(active) > 1:
if len(active) >= 120000:
nseen = bytearray(mod); nxt = []
for a in active:
v = a*a % mod
if v > 1 and not nseen[v]: nseen[v] = 1; nxt.append(v)
else:
ns = set()
for a in active:
v = a*a % mod
if v > 1: ns.add(v)
nxt = list(ns)
active = nxt; length += 1; mod += 1
if not active: return length
if not active: return length
v = active[0]
while v > 1: v = v*v % mod; length += 1; mod += 1
return length
import math
def f_value(n):
if n < 2: return 0
pts = list(range(2, n+1, max(1, (n-2)//15)))
if pts[-1] != n: pts.append(n)
gpts = [g_value(x) for x in pts]
cache = dict(zip(pts, gpts))
best = max(gpts)
import heapq
pq = []
for i in range(1, len(pts)):
ub = gpts[i] + (pts[i]-pts[i-1])
heapq.heappush(pq, (-ub, pts[i-1], pts[i], gpts[i]))
while pq:
neg_ub, left, right, gr = heapq.heappop(pq)
if -neg_ub <= best: break
if right - left <= 1: continue
mid = (left+right)//2
gm = cache.get(mid)
if gm is None: gm = g_value(mid); cache[mid] = gm
if gm > best: best = gm
if mid-left > 1:
ub_l = gm + (mid-left)
if ub_l > best: heapq.heappush(pq, (-ub_l, left, mid, gm))
if right-mid > 1:
ub_r = gr + (right-mid)
if ub_r > best: heapq.heappush(pq, (-ub_r, mid, right, gr))
return best
return str(f_value(N))
if __name__ == '__main__':
print(solve())
Java
import java.util.*;
import java.util.concurrent.*;
public class Euler693 {
static List<Integer> initialActiveAfterFirstStep(int x) {
if (x <= 2)
return new ArrayList<>();
byte[] seen = new byte[x];
List<Integer> active = new ArrayList<>();
long mod = x;
long limit = mod / 2;
long y = 2;
long sq = (y * y) % mod;
long delta = 2 * y + 1;
while (y <= limit) {
if (sq > 1 && seen[(int) sq] == 0) {
seen[(int) sq] = 1;
active.add((int) sq);
}
sq += delta;
if (sq >= mod) {
sq -= mod;
if (sq >= mod)
sq -= mod;
}
delta += 2;
y++;
}
return active;
}
static List<Integer> stepActiveDense(List<Integer> active, int mod) {
byte[] seen = new byte[mod];
List<Integer> next = new ArrayList<>();
for (int a : active) {
long v = ((long) a * a) % mod;
if (v > 1 && seen[(int) v] == 0) {
seen[(int) v] = 1;
next.add((int) v);
}
}
return next;
}
static List<Integer> stepActiveSparse(List<Integer> active, int mod) {
HashSet<Integer> nextSet = new HashSet<>();
for (int a : active) {
long v = ((long) a * a) % mod;
if (v > 1) {
nextSet.add((int) v);
}
}
return new ArrayList<>(nextSet);
}
static int gValue(int x) {
if (x < 2)
return 0;
if (x == 2)
return 1;
List<Integer> active = initialActiveAfterFirstStep(x);
int length = 2;
int mod = x + 1;
int denseThreshold = 120000;
while (active.size() > 1) {
if (active.size() >= denseThreshold) {
active = stepActiveDense(active, mod);
} else {
active = stepActiveSparse(active, mod);
}
length++;
mod++;
if (active.isEmpty())
return length;
}
if (active.isEmpty())
return length;
long value = active.get(0);
while (value > 1) {
value = (value * value) % mod;
length++;
mod++;
}
return length;
}
static int snappedGridStep(int rawStep) {
if (rawStep <= 1)
return 1;
int power = (int) Math.floor(Math.log10(rawStep));
int base = 1;
for (int i = 0; i < power; i++)
base *= 10;
int[] mults = { 1, 2, 5, 10 };
for (int m : mults) {
if ((long) m * base >= rawStep) {
return m * base;
}
}
return rawStep;
}
static class Interval implements Comparable<Interval> {
int upperBound, left, right, gRight;
public Interval(int ub, int l, int r, int gr) {
this.upperBound = ub;
this.left = l;
this.right = r;
this.gRight = gr;
}
@Override
public int compareTo(Interval o) {
return Integer.compare(o.upperBound, this.upperBound); // descending
}
}
public static String solve() {
int n = 3000000;
int targetPoints = 16;
int pointsCount = Math.max(2, targetPoints);
int rawStep = Math.max(1, (n - 2) / (pointsCount - 1));
int gridStep = snappedGridStep(rawStep);
List<Integer> points = new ArrayList<>();
for (long x = 2; x <= n; x += gridStep) {
points.add((int) x);
if (x > (long) n - gridStep)
break;
}
if (points.get(points.size() - 1) != n) {
points.add(n);
}
int threads = Runtime.getRuntime().availableProcessors();
if (threads < 1)
threads = 1;
int[] gPoints = new int[points.size()];
ExecutorService pool = Executors.newFixedThreadPool(threads);
List<Future<?>> futures = new ArrayList<>();
for (int i = 0; i < points.size(); i++) {
final int idx = i;
futures.add(pool.submit(() -> {
gPoints[idx] = gValue(points.get(idx));
}));
}
for (Future<?> f : futures) {
try {
f.get();
} catch (Exception e) {
}
}
pool.shutdown();
ConcurrentHashMap<Integer, Integer> cache = new ConcurrentHashMap<>();
int best = 0;
for (int i = 0; i < points.size(); i++) {
cache.put(points.get(i), gPoints[i]);
if (gPoints[i] > best)
best = gPoints[i];
}
PriorityQueue<Interval> pq = new PriorityQueue<>();
for (int i = 1; i < points.size(); i++) {
int left = points.get(i - 1);
int right = points.get(i);
int gRight = gPoints[i];
int upperBound = gRight + (right - left);
pq.add(new Interval(upperBound, left, right, gRight));
}
while (!pq.isEmpty()) {
Interval top = pq.poll();
if (top.upperBound <= best)
break;
if (top.right - top.left <= 1)
continue;
int mid = top.left + (top.right - top.left) / 2;
int gMid = cache.computeIfAbsent(mid, m -> gValue(m));
if (gMid > best)
best = gMid;
if (mid - top.left > 1) {
int ubLeft = gMid + (mid - top.left);
if (ubLeft > best)
pq.add(new Interval(ubLeft, top.left, mid, gMid));
}
if (top.right - mid > 1) {
int ubRight = top.gRight + (top.right - mid);
if (ubRight > best)
pq.add(new Interval(ubRight, mid, top.right, top.gRight));
}
}
return Integer.toString(best);
}
public static void main(String[] args) {
System.out.println(solve());
}
}