Problem 950: Pirate Treasure
View on Project EulerProject Euler Problem 950 Solution
EulerSolve provides an optimized solution for Project Euler Problem 950, Pirate Treasure, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Problem 950 asks for the quantity \(T(N,C,D)\) at an extreme scale. The implementations evaluate \[ \sum_{k=1}^{6} T(10^{16},10^k+1,10^k+1) \pmod{10^9}, \] so even one parameter set is far beyond direct simulation. What the code actually sums is a sequence \((c_n)_{n\ge 1}\) with \[ T(N,C,D)=\sum_{n=1}^{N} c_n, \qquad c_1=C. \] The hard part is that the next exceptional value does not depend only on the current term; it depends on a whole multiset of carried values. The successful strategy is to keep that multiset in compressed form, jump directly to the next index where a nontrivial update occurs, and sum the long linear stretches in between with a closed formula. Mathematical Approach The Event-State Invariant At a landing position \(s\), the algorithm maintains two pieces of information: the current sequence value \(c_s\), and a multiset \(\mathcal{M}_s\) of size \(s\) that controls every future transition. Only positive entries are stored explicitly. Zeros are tracked by a counter \(z_s\). Write the multiset as \[ \mathcal{M}_s=\{0^{z_s}\}\cup \{v_1^{m_1},v_2^{m_2},\dots,v_r^{m_r}\}, \qquad 0 \lt v_1 \lt \cdots \lt v_r. \] Here the exponent means multiplicity: \(v_i^{m_i}\) means that the value \(v_i\) occurs \(m_i\) times. The invariant is \[ z_s+\sum_{i=1}^{r} m_i=s....
Detailed mathematical approach
Problem Summary
Problem 950 asks for the quantity \(T(N,C,D)\) at an extreme scale. The implementations evaluate
\[ \sum_{k=1}^{6} T(10^{16},10^k+1,10^k+1) \pmod{10^9}, \]
so even one parameter set is far beyond direct simulation. What the code actually sums is a sequence \((c_n)_{n\ge 1}\) with
\[ T(N,C,D)=\sum_{n=1}^{N} c_n, \qquad c_1=C. \]
The hard part is that the next exceptional value does not depend only on the current term; it depends on a whole multiset of carried values. The successful strategy is to keep that multiset in compressed form, jump directly to the next index where a nontrivial update occurs, and sum the long linear stretches in between with a closed formula.
Mathematical Approach
The Event-State Invariant
At a landing position \(s\), the algorithm maintains two pieces of information: the current sequence value \(c_s\), and a multiset \(\mathcal{M}_s\) of size \(s\) that controls every future transition. Only positive entries are stored explicitly. Zeros are tracked by a counter \(z_s\).
Write the multiset as
\[ \mathcal{M}_s=\{0^{z_s}\}\cup \{v_1^{m_1},v_2^{m_2},\dots,v_r^{m_r}\}, \qquad 0 \lt v_1 \lt \cdots \lt v_r. \]
Here the exponent means multiplicity: \(v_i^{m_i}\) means that the value \(v_i\) occurs \(m_i\) times. The invariant is
\[ z_s+\sum_{i=1}^{r} m_i=s. \]
The key summary statistic derived from this state is
\[ P_s(b)=\text{the sum of the } b \text{ smallest elements of }\mathcal{M}_s. \]
Because zeros come first and the positive values are stored in increasing groups, \(P_s(b)\) can be read off by a short scan of the grouped state instead of by expanding all \(s\) entries.
Feasible Jumps and the Budget Test
From a landing at index \(s\), the next nontrivial landing is searched in the form \(s+d\) with \(1\le d\le s\). For each candidate distance \(d\), the implementations define
\[ b=\left\lfloor\frac{s-d+1}{2}\right\rfloor, \qquad g=\left\lceil\frac{d}{\sqrt D}\right\rceil. \]
The quantity \(b\) tells us how many of the smallest carried values must be consumed, while \(g\) is the common surcharge attached to a jump of length \(d\). The total cost of using that jump is
\[ \operatorname{cost}(d)=P_s(b)+b\,g. \]
The jump is feasible exactly when
\[ \operatorname{cost}(d)\le C. \]
The first feasible \(d\) is chosen. This is the decisive observation in the code: once the earliest feasible landing is known, every index before it follows the trivial linear rule and does not need to be processed one by one.
Why the Search Starts at \(s-2C\) and Always Terminates
The scan does not start from \(d=1\) blindly. Since \(g\ge 1\) and \(P_s(b)\ge 0\), any feasible jump must satisfy \(b\le C\). Using
\[ b=\left\lfloor\frac{s-d+1}{2}\right\rfloor, \]
this implies
\[ d\ge s-2C. \]
So the candidate range can be shortened to
\[ d\in[\max(1,s-2C),\,s]. \]
There is also a guaranteed fallback. When \(d=s\), we get
\[ b=\left\lfloor\frac{s-s+1}{2}\right\rfloor=0, \]
hence \(\operatorname{cost}(s)=0\), which is always feasible. Therefore every state has a next landing. This special case produces a full reset: the next landing is at \(2s\), the new sequence value becomes \(C\), and the positive part of the state collapses to a single copy of \(C\).
Exact Landing-State Update
Assume the chosen jump has \(b>0\). Let \(z_{\mathrm{sel}}=\min(b,z_s)\); these are the selected zeros among the \(b\) smallest elements. The remaining \(b-z_{\mathrm{sel}}\) selected elements are the smallest positive values of \(\mathcal{M}_s\), say
\[ x_1,\dots,x_{b-z_{\mathrm{sel}}}. \]
The landing value is then
\[ c_{s+d}=C-\operatorname{cost}(d). \]
The new positive entries are built as follows:
\[ \underbrace{\{g,\dots,g\}}_{z_{\mathrm{sel}}\text{ copies}} \cup \{x_1+g,\dots,x_{b-z_{\mathrm{sel}}}+g\}. \]
If \(c_{s+d}>0\), one further copy of \(c_{s+d}\) is inserted. Everything else is zero, so the new zero count is simply whatever is needed to make the total size equal to \(s+d\). Finally, equal positive values are merged back into grouped multiplicities. That normalization step is what keeps the state compact.
Arithmetic Blocks Between Landings
If the first feasible jump length is \(d\), then there is no special landing at any of the intermediate indices \(s+1,s+2,\dots,s+d-1\). Their values follow the linear pattern
\[ c_{s+t}=c_s+t \qquad (1\le t \lt d). \]
So the skipped contribution is a simple arithmetic progression:
\[ \sum_{t=1}^{d-1}(c_s+t) = (d-1)c_s+\frac{(d-1)d}{2}. \]
This is the reason the algorithm can jump over huge stretches of indices. Once the next landing is known, the entire block before it is summed in \(O(1)\) time.
Worked Example: \(C=D=3\)
The checkpoint case \(T(30,3,3)=190\) already shows the two essential behaviors. After a few updates the state at \(s=4\) is
\[ c_4=2, \qquad \mathcal{M}_4=\{0,0,1,2\}. \]
Try \(d=1\). Then
\[ b=\left\lfloor\frac{4-1+1}{2}\right\rfloor=2, \qquad g=\left\lceil\frac{1}{\sqrt 3}\right\rceil=1, \qquad P_4(2)=0, \]
because the two smallest multiset elements are both zero. Hence
\[ \operatorname{cost}(1)=2, \qquad c_5=3-2=1. \]
The two selected zeros become two copies of \(1\), and the landing value contributes one more \(1\). So the next grouped state is
\[ \mathcal{M}_5=\{0,0,1,1,1\}. \]
Later, at \(s=8\), every \(d<8\) fails the budget test, so the fallback \(d=8\) is used. That means indices \(9\) through \(15\) are just the arithmetic run \(1,2,3,4,5,6,7\) above \(c_8=0\), and index \(16\) resets to \(c_{16}=3\). This single trace explains why the code alternates between exact landing updates and long linear blocks.
How the Code Works
The C++ and Python implementations maintain the landing index \(s\), the current value \(c_s\), the zero count, and the sorted positive groups. For each landing they scan candidate distances \(d\) from \(\max(1,s-2C)\) upward. To compute \(g=\lceil d/\sqrt D\rceil\) safely, they start from a numerical estimate and then correct it with the exact integer inequality \(g^2D\ge d^2\). Since \(g\) is monotone in \(d\), it can then be updated incrementally.
For each candidate, the implementation evaluates \(P_s(b)\) by consuming zeros first and then walking through the positive groups from smallest to largest. The scan stops early as soon as the partial sum exceeds \(C\), because such a candidate can never be feasible. Once the first feasible \(d\) is found, the code either performs the reset case \(d=s\) or constructs the next grouped multiset from the selected smallest elements, then sorts and merges equal positive values again.
The outer accumulation loop never expands all \(N\) indices. It adds the skipped arithmetic block in closed form, adds the landing value \(c_{s+d}\), and repeats from the new state. The Python version is a direct compact transcription of the same logic. The Java version serves as a thin launcher around the same computation, compiling and invoking the C++ implementation so that all three languages produce the same numerical result.
Complexity Analysis
Let \(J\) be the number of landing states actually visited. At landing \(j\), let \(R_j\) be the number of candidate distances examined before the first feasible one is found, and let \(G_j\) be the number of distinct positive groups in the compressed state. Both the prefix-sum test and the landing update inspect at most \(G_j\) groups, so a safe worst-case bound is
\[ O\!\left(\sum_{j=1}^{J} R_j G_j\right). \]
The memory usage is \(O(\max_j G_j)\), because zeros are implicit and equal positive values are merged. That is the decisive compression: a state with enormous index \(s\) may still be represented by only a few grouped values. The final program repeats this computation for six parameter pairs, which is only a constant-factor multiplier.
The crucial point is qualitative rather than asymptotic. A naive method would need \(10^{16}\) sequence updates per parameter set, while this method processes only the landing states and treats the vast linear stretches analytically.
Footnotes and References
- Problem page: Project Euler 950 - Pirate Treasure
- Multiset: Wikipedia - Multiset
- Arithmetic progression: Wikipedia - Arithmetic progression
- Floor and ceiling functions: Wikipedia - Floor and ceiling functions
Problem 950 source code
C++
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <utility>
#include <vector>
namespace {
using u32 = std::uint32_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;
u64 ceil_div_sqrt(u64 d, u64 D) {
long double approx = static_cast<long double>(d) / std::sqrt(static_cast<long double>(D));
u64 g = static_cast<u64>(approx);
if (g == 0) {
g = 1;
}
const u128 dd = static_cast<u128>(d) * static_cast<u128>(d);
while (static_cast<u128>(g) * static_cast<u128>(g) * static_cast<u128>(D) < dd) {
++g;
}
while (g > 1 && static_cast<u128>(g - 1) * static_cast<u128>(g - 1) * static_cast<u128>(D) >= dd) {
--g;
}
return g;
}
struct Entry {
u32 value;
u64 count;
};
struct State {
u64 s = 0;
u64 C = 0;
u64 D = 0;
u64 z = 0;
std::vector<Entry> groups;
State() = default;
State(u64 s_, u64 C_, u64 D_, const std::vector<Entry>& input_groups)
: s(s_), C(C_), D(D_) {
groups = input_groups;
normalize();
}
void normalize() {
std::vector<Entry> cleaned;
cleaned.reserve(groups.size());
for (const auto& e : groups) {
if (e.value > 0 && e.count > 0) {
cleaned.push_back(e);
}
}
std::sort(cleaned.begin(), cleaned.end(),
[](const Entry& a, const Entry& b) { return a.value < b.value; });
groups.clear();
groups.reserve(cleaned.size());
for (const auto& e : cleaned) {
if (!groups.empty() && groups.back().value == e.value) {
groups.back().count += e.count;
} else {
groups.push_back(e);
}
}
u64 pos_count = 0;
for (const auto& e : groups) {
pos_count += e.count;
}
assert(pos_count <= s);
z = s - pos_count;
}
u64 prefix_smallest(u64 b) const {
if (b <= z) {
return 0ULL;
}
u64 remaining = b - z;
u64 sum = 0;
for (const auto& e : groups) {
if (remaining == 0) {
break;
}
const u64 take = std::min(remaining, e.count);
sum += take * static_cast<u64>(e.value);
if (sum > C) {
return C + 1ULL;
}
remaining -= take;
}
if (remaining > 0) {
return C + 1ULL;
}
return sum;
}
std::pair<State, u64> advance(u64 d, u64 b, u64 g, u64 cost) const {
if (b == 0) {
return {State(s + d, C, D, {{static_cast<u32>(C), 1ULL}}), C};
}
const u64 z_sel = std::min(b, z);
u64 rem = b - z_sel;
std::vector<Entry> selected;
selected.reserve(groups.size());
for (const auto& e : groups) {
if (rem == 0) {
break;
}
const u64 take = std::min(rem, e.count);
if (take > 0) {
selected.push_back({e.value, take});
rem -= take;
}
}
assert(rem == 0);
const u64 c = C - cost;
std::vector<Entry> next_entries;
next_entries.reserve(selected.size() + 2);
if (c > 0) {
next_entries.push_back({static_cast<u32>(c), 1ULL});
}
if (z_sel > 0) {
next_entries.push_back({static_cast<u32>(g), z_sel});
}
for (const auto& e : selected) {
next_entries.push_back({static_cast<u32>(static_cast<u64>(e.value) + g), e.count});
}
return {State(s + d, C, D, next_entries), c};
}
};
u128 compute_T(u64 N, u64 C, u64 D) {
State st(1ULL, C, D, {{static_cast<u32>(C), 1ULL}});
u64 current_c = C;
u128 total = C;
while (st.s < N) {
const u64 s = st.s;
const u64 lo = (s > 2ULL * C) ? (s - 2ULL * C) : 1ULL;
u64 found_d = s;
u64 found_b = 0;
u64 found_g = 1;
u64 found_cost = 0;
u64 g = ceil_div_sqrt(lo, D);
for (u64 d = lo; d <= s; ++d) {
const u128 dd = static_cast<u128>(d) * static_cast<u128>(d);
while (static_cast<u128>(g) * static_cast<u128>(g) * static_cast<u128>(D) < dd) {
++g;
}
const u64 b = (s - d + 1ULL) / 2ULL;
if (b == 0) {
found_d = d;
found_b = 0;
found_g = g;
found_cost = 0;
break;
}
const u64 ps = st.prefix_smallest(b);
if (ps > C) {
continue;
}
const u128 cost128 = static_cast<u128>(b) * static_cast<u128>(g) + static_cast<u128>(ps);
if (cost128 <= C) {
found_d = d;
found_b = b;
found_g = g;
found_cost = static_cast<u64>(cost128);
break;
}
}
const u64 next_s = s + found_d;
if (next_s > N) {
const u64 cnt = N - s;
total += static_cast<u128>(cnt) * static_cast<u128>(current_c);
total += static_cast<u128>(cnt) * static_cast<u128>(cnt + 1ULL) / 2ULL;
break;
}
const u64 cnt = next_s - s - 1ULL;
if (cnt > 0) {
total += static_cast<u128>(cnt) * static_cast<u128>(current_c);
total += static_cast<u128>(cnt) * static_cast<u128>(cnt + 1ULL) / 2ULL;
}
auto [next_state, next_c] = st.advance(found_d, found_b, found_g, found_cost);
total += next_c;
st = std::move(next_state);
current_c = next_c;
}
return total;
}
void run_validations() {
assert(compute_T(30ULL, 3ULL, 3ULL) == 190ULL);
assert(compute_T(50ULL, 3ULL, 31ULL) == 385ULL);
assert(compute_T(1'000ULL, 101ULL, 101ULL) == 142'427ULL);
}
} // namespace
int main() {
run_validations();
constexpr u64 kN = 10'000'000'000'000'000ULL;
u128 sum = 0;
u64 pow10 = 1;
for (int k = 1; k <= 6; ++k) {
pow10 *= 10ULL;
const u64 C = pow10 + 1ULL;
sum += compute_T(kN, C, C);
}
const u64 answer = static_cast<u64>(sum % 1'000'000'000ULL);
std::cout << std::setw(9) << std::setfill('0') << answer << '\n';
return 0;
}
Python
import math
def solve():
def ceil_div_sqrt(d, D):
approx = d / math.sqrt(D); g = max(1, int(approx))
dd = d*d
while g*g*D < dd: g += 1
while g > 1 and (g-1)*(g-1)*D >= dd: g -= 1
return g
def compute_T(N, C, D):
groups = [(C, 1)]; s = 1; z = 0; current_c = C; total = C
while s < N:
lo = max(1, s - 2*C)
g = ceil_div_sqrt(lo, D)
found_d = s; found_b = 0; found_g = 1; found_cost = 0
for d in range(lo, s+1):
dd = d*d
while g*g*D < dd: g += 1
b = (s-d+1)//2
if b == 0:
found_d = d; found_b = 0; found_g = g; found_cost = 0; break
# prefix_smallest(b)
ps = 0; rem = b
if rem <= z: ps = 0; rem = 0
else: rem -= z
for val, cnt in sorted(groups):
if rem == 0: break
take = min(rem, cnt); ps += take*val; rem -= take
if ps > C: ps = C+1; break
if rem > 0: ps = C+1
if ps > C: continue
cost128 = b*g + ps
if cost128 <= C:
found_d = d; found_b = b; found_g = g; found_cost = cost128; break
next_s = s + found_d
if next_s > N:
cnt = N-s; total += cnt*current_c + cnt*(cnt+1)//2; break
cnt = next_s-s-1
if cnt > 0: total += cnt*current_c + cnt*(cnt+1)//2
# advance
if found_b == 0:
groups = [(C, 1)]; z = next_s-1; current_c = C
else:
z_sel = min(found_b, z); rem = found_b - z_sel
selected = []
for val, cnt2 in sorted(groups):
if rem == 0: break
take = min(rem, cnt2); selected.append((val, take)); rem -= take
c_new = C - found_cost
new_entries = []
if c_new > 0: new_entries.append((c_new, 1))
if z_sel > 0: new_entries.append((found_g, z_sel))
for val, cnt2 in selected: new_entries.append((val+found_g, cnt2))
# normalize
from collections import defaultdict
d2 = defaultdict(int)
for v, c in new_entries:
if v > 0 and c > 0: d2[v] += c
groups = sorted(d2.items())
pos_count = sum(c for _, c in groups)
z = next_s - pos_count; current_c = c_new
total += current_c; s = next_s
return total
N = 10**16; ans = 0
pow10 = 1
for k in range(1, 7):
pow10 *= 10; C = pow10+1
ans += compute_T(N, C, C)
return str(ans % 1000000000).zfill(9)
if __name__ == '__main__':
print(solve())
Java
import java.nio.file.*;
import java.util.*;
import java.util.regex.*;
public class Euler950 {
private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");
private static String parseOutput(String stdout) {
String[] lines = stdout.split("\\R");
List<String> nonEmpty = new ArrayList<>();
for (String line : lines) {
String t = line.trim();
if (!t.isEmpty()) {
nonEmpty.add(t);
}
}
if (nonEmpty.isEmpty()) {
return "";
}
List<String> answers = new ArrayList<>();
List<String> equals = new ArrayList<>();
for (String line : nonEmpty) {
Matcher m1 = ANSWER_RE.matcher(line);
if (m1.find()) {
answers.add(m1.group(1).trim());
}
Matcher m2 = EQUAL_RE.matcher(line);
if (m2.find()) {
equals.add(m2.group(1).trim());
}
}
if (!answers.isEmpty()) {
return answers.get(answers.size() - 1);
}
if (!equals.isEmpty()) {
return equals.get(equals.size() - 1);
}
return nonEmpty.get(nonEmpty.size() - 1);
}
private static String pickCompiler() throws Exception {
for (String compiler : List.of("clang++", "g++")) {
Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
.redirectErrorStream(true)
.start();
String out = new String(probe.getInputStream().readAllBytes());
int rc = probe.waitFor();
if (rc == 0 && !out.trim().isEmpty()) {
return compiler;
}
}
throw new RuntimeException("No C++ compiler found (clang++/g++).");
}
private static Path cppSource(Path root) {
return root.resolve("solutionsCpp").resolve("Euler950.cpp");
}
private static boolean shouldSkipCheckpoints(Path root) {
Path src = cppSource(root);
try {
String text = Files.readString(src);
return text.contains("--skip-checkpoints");
} catch (Exception ex) {
return false;
}
}
private static Path ensureBridgeBinary() throws Exception {
Path root = Paths.get(System.getProperty("user.dir"));
Path src = cppSource(root);
Path bin = root.resolve("solutionsCpp").resolve(".euler950_java_bridge");
boolean rebuild = Files.notExists(bin)
|| Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;
if (rebuild) {
String compiler = pickCompiler();
Process compile = new ProcessBuilder(
compiler,
"-std=c++17",
"-O2",
src.toString(),
"-o",
bin.toString())
.inheritIO()
.start();
if (compile.waitFor() != 0) {
throw new RuntimeException("Failed to compile Euler950 C++ bridge.");
}
}
return bin;
}
private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
List<String> cmd = new ArrayList<>();
cmd.add(bin.toString());
if (shouldSkipCheckpoints(root)) {
cmd.add("--skip-checkpoints");
}
Process first = new ProcessBuilder(cmd)
.directory(root.toFile())
.redirectErrorStream(true)
.start();
String out = new String(first.getInputStream().readAllBytes());
int rc = first.waitFor();
if (rc == 0) {
return out;
}
Process second = new ProcessBuilder(cmd)
.directory(srcDir.toFile())
.redirectErrorStream(true)
.start();
String out2 = new String(second.getInputStream().readAllBytes());
int rc2 = second.waitFor();
if (rc2 == 0) {
return out2;
}
throw new RuntimeException("Euler950 C++ bridge failed.\n" + out + "\n" + out2);
}
private static String solveViaCppBridge() throws Exception {
Path root = Paths.get(System.getProperty("user.dir"));
Path src = cppSource(root);
Path bin = ensureBridgeBinary();
String out = runBridge(bin, root, src.getParent());
String parsed = parseOutput(out);
if (parsed.isEmpty()) {
throw new RuntimeException("Euler950 C++ bridge produced empty output.");
}
return parsed;
}
public static void main(String[] args) throws Exception {
System.out.println(solveViaCppBridge());
}
}