Problem 758: Buckets of Water

Project Euler Problem 758 Solution

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

Problem Summary For coprime positive integers $a\le b$, the quantity studied here is built from Bézout relations between $a$ and $b$. In the Project Euler instance, the inputs are Mersenne numbers $M_u=2^u-1$ and $M_v=2^v-1$, where $u=r^5$ and $v=s^5$ for distinct primes $r\lt s\lt 1000$. The goal is to evaluate the contribution for every such prime pair and sum all contributions modulo $10^9+7$. Directly constructing $2^u-1$ and $2^v-1$ is hopeless because the exponents are enormous. The key observation is that the Euclidean algorithm on Mersenne numbers can be simulated on the exponents alone, while every coefficient that matters is tracked only modulo the final modulus. Mathematical Approach Write $M_n=2^n-1$. For two coprime positive integers $a\le b$, choose the least positive inverse $x$ of $a$ modulo $b$: $$a x \equiv 1 \pmod{b},\qquad 1\le x \lt b.$$ Then $$y=\frac{a x - 1}{b}$$ is a positive integer and satisfies $a x - b y = 1$. The complementary positive solution to the opposite sign equation is $(b-x,\ a-y)$, so the implementations evaluate $$P(a,b)=2\min\left(x+y-1,\ (b-x)+(a-y)-1\right).$$ The problem is therefore to compute $P(M_u,M_v)$ quickly for very large exponents $u$ and $v$. Step 1: Replace gigantic integers by Mersenne exponents Let $a=M_u$ and $b=M_v$ with $u\lt v$....

Detailed mathematical approach

Problem Summary

For coprime positive integers $a\le b$, the quantity studied here is built from Bézout relations between $a$ and $b$. In the Project Euler instance, the inputs are Mersenne numbers $M_u=2^u-1$ and $M_v=2^v-1$, where $u=r^5$ and $v=s^5$ for distinct primes $r\lt s\lt 1000$. The goal is to evaluate the contribution for every such prime pair and sum all contributions modulo $10^9+7$.

Directly constructing $2^u-1$ and $2^v-1$ is hopeless because the exponents are enormous. The key observation is that the Euclidean algorithm on Mersenne numbers can be simulated on the exponents alone, while every coefficient that matters is tracked only modulo the final modulus.

Mathematical Approach

Write $M_n=2^n-1$. For two coprime positive integers $a\le b$, choose the least positive inverse $x$ of $a$ modulo $b$:

$$a x \equiv 1 \pmod{b},\qquad 1\le x \lt b.$$

Then

$$y=\frac{a x - 1}{b}$$

is a positive integer and satisfies $a x - b y = 1$. The complementary positive solution to the opposite sign equation is $(b-x,\ a-y)$, so the implementations evaluate

$$P(a,b)=2\min\left(x+y-1,\ (b-x)+(a-y)-1\right).$$

The problem is therefore to compute $P(M_u,M_v)$ quickly for very large exponents $u$ and $v$.

Step 1: Replace gigantic integers by Mersenne exponents

Let $a=M_u$ and $b=M_v$ with $u\lt v$. If $u$ and $v$ are coprime, then the classical identity

$$\gcd(M_u,M_v)=M_{\gcd(u,v)}$$

gives

$$\gcd(M_u,M_v)=M_1=1.$$

This is exactly the situation here, because distinct prime fifth powers remain coprime:

$$\gcd(r^5,s^5)=1\qquad (r\ne s).$$

So every target pair is valid for the Bézout setup, even though the numbers themselves are astronomically large.

Step 2: Run the Euclidean algorithm on exponents

If $A=qB+r$ with $0\le r \lt B$, then

$$2^A-1=2^r\left(2^{qB}-1\right)+\left(2^r-1\right).$$

After factoring $2^{qB}-1=(2^B-1)\sum_{j=0}^{q-1}2^{jB}$, this becomes

$$M_A = Q\,M_B + M_r,$$

with quotient

$$Q = 2^r\sum_{j=0}^{q-1}2^{jB}.$$

Therefore the remainder when dividing $M_A$ by $M_B$ is exactly $M_r$. The full remainder sequence for $(M_A,M_B)$ is obtained by applying ordinary Euclid to the exponent pair $(A,B)$. This is the compression that makes the problem tractable.

Step 3: Compute Euclidean quotients without big integers

Each Euclidean step needs the quotient $Q$ only modulo $10^9+7$. Using the formula above, we only need

$$Q \equiv 2^r\sum_{j=0}^{q-1}\left(2^B\right)^j \pmod{10^9+7}.$$

The implementation evaluates this geometric sum by divide and conquer, simultaneously tracking a power and the corresponding partial sum. As a result, every quotient update is handled with modular arithmetic and logarithmic recursion depth, not with gigantic integers.

Feeding these modular quotients into the extended Euclidean recurrence yields the inverse of a Mersenne remainder modulo another Mersenne number, but represented only modulo $10^9+7$.

Step 4: Recover the inverse of the Mersenne remainder

Let

$$w = v \bmod u,$$

and let $t$ be the inverse of $w$ modulo $u$:

$$w t \equiv 1 \pmod{u},\qquad 1\le t \lt u.$$

Now consider

$$R = 1 + 2^w + 2^{2w} + \cdots + 2^{(t-1)w}.$$

Then

$$\begin{aligned} (2^w-1)R &= 2^{tw}-1 \\ &\equiv 2^1-1 \equiv 1 \pmod{2^u-1}, \end{aligned}$$

so $R$ is an inverse of $M_w$ modulo $M_u$. Because $M_v\equiv M_w \pmod{M_u}$, the same value also acts as an inverse of $M_v$ modulo $M_u$.

The opposite Bézout sign branch is obtained by replacing $R$ with $M_u-R$. The implementation chooses between these two branches with the criterion

$$2t\le u,$$

which is equivalent to choosing between the $t$-term geometric sum and its $(u-t)$-term complement inside $M_u=1+2+\cdots+2^{u-1}$.

Step 5: Convert the inverse into the required Bézout sum

Let $d$ be the branch selected in the previous step. Then one of the two relations

$$M_v d = M_u x \pm 1.$$

holds, depending on the chosen sign. Once $x$ is recovered, the contribution is

$$P(M_u,M_v)=2(x+d-1).$$

This is the same minimum as in the generic formula for $P(a,b)$; the only difference is that the implementation works with the coefficient attached to $M_v$, because that is the quantity naturally produced by the Mersenne inverse computation.

Step 6: Worked example with small exponents

Take $u=3$ and $v=5$. Then

$$M_3=7,\qquad M_5=31,\qquad w=5\bmod 3=2.$$

The inverse of $w=2$ modulo $u=3$ is $t=2$, since

$$2\cdot 2 \equiv 1 \pmod{3}.$$

Hence

$$R = 1 + 2^2 = 5,$$

and indeed

$$M_2 R = 3\cdot 5 = 15 \equiv 1 \pmod{7}.$$

Because $2t=4>3$, the smaller branch uses the complement

$$d = M_3 - R = 7-5=2.$$

Now

$$x=\frac{M_5 d + 1}{M_3}=\frac{31\cdot 2 + 1}{7}=9,$$

$$P(M_3,M_5)=2(x+d-1)=2(9+2-1)=20.$$

This matches the small verification case built into the implementation.

Step 7: Final summation over prime fifth powers

If $p_1,p_2,\dots,p_k$ are the primes below $1000$, the final value is

$$\sum_{1\le i \lt j\le k} P\left(M_{p_i^5},M_{p_j^5}\right)\pmod{10^9+7}.$$

The entire mathematical task is therefore to evaluate each pair contribution via exponent-level Euclid and accumulate the results modulo the final prime modulus.

How the Code Works

The C++, Python, and Java implementations all expose the same computation, but the actual number-theoretic work is performed by the C++ implementation. It begins with several small checkpoints: first on ordinary coprime integer pairs, and then on small Mersenne exponents where a direct comparison against the generic Bézout formula is still possible.

After the checkpoints, the implementation generates all primes below $1000$, raises each of them to the fifth power, and iterates over every unordered pair. For each pair $(u,v)$, it computes $M_u \bmod 10^9+7$ and $M_v \bmod 10^9+7$, runs the Euclidean algorithm on the exponent pair, and reconstructs the required quotient information only modulo $10^9+7$.

The key accelerator is the divide-and-conquer evaluation of geometric sums. Instead of forming the true quotient between two huge Mersenne numbers, the implementation directly produces that quotient modulo $10^9+7$, which is exactly what the extended Euclidean recurrence needs. The branch decision based on the inverse of $w$ modulo $u$ then determines which Bézout sign gives the smaller contribution.

Finally, the contribution of every prime pair is added into a running sum modulo $10^9+7$. The Python and Java implementations are thin launchers around that same compiled computation, so all three language entries stay synchronized on both the mathematics and the final output.

Complexity Analysis

There are $168$ primes below $1000$, so the outer loop evaluates $\binom{168}{2}=14028$ unordered pairs. For one pair, the exponent-level Euclidean algorithm takes $O(\log \max(u,v))$ divisions, just as ordinary Euclid does on integers. Each division performs a constant number of modular exponentiations and one divide-and-conquer geometric-sum evaluation, both logarithmic in the relevant exponent sizes. Therefore the total running time is $O(P^2\operatorname{polylog} U)$ with $P=168$ and $U=\max p^5$, while memory usage is $O(P)$ for the prime list plus $O(\log U)$ recursion depth.

Footnotes and References

Problem page: https://projecteuler.net/problem=758
Bézout's identity: Wikipedia — Bézout's identity
Extended Euclidean algorithm: Wikipedia — Extended Euclidean algorithm
Mersenne number: Wikipedia — Mersenne number
Modular multiplicative inverse: Wikipedia — Modular multiplicative inverse
Geometric series: Wikipedia — Geometric series

Problem 758 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <utility>
#include <vector>

namespace {

using i64 = long long;
using i128 = __int128_t;
using u64 = std::uint64_t;

constexpr i64 MOD = 1'000'000'007LL;

i64 mod_pow(i64 base, u64 exp) {
    i64 result = 1 % MOD;
    i64 cur = base % MOD;
    while (exp > 0) {
        if (exp & 1ULL) {
            result = static_cast<i64>((static_cast<i128>(result) * cur) % MOD);
        }
        cur = static_cast<i64>((static_cast<i128>(cur) * cur) % MOD);
        exp >>= 1ULL;
    }
    return result;
}

i64 mod_inv(i64 x) {
    assert(x % MOD != 0);
    return mod_pow((x % MOD + MOD) % MOD, static_cast<u64>(MOD - 2));
}

std::pair<i64, i64> pow_sum(i64 ratio, u64 n) {
    if (n == 0ULL) {
        return {1LL, 0LL};
    }
    if (n == 1ULL) {
        return {ratio % MOD, 1LL};
    }
    if ((n & 1ULL) == 0ULL) {
        auto [p, s] = pow_sum(ratio, n >> 1ULL);
        const i64 p2 = static_cast<i64>((static_cast<i128>(p) * p) % MOD);
        const i64 s2 = static_cast<i64>((static_cast<i128>(s) * ((1LL + p) % MOD)) % MOD);
        return {p2, s2};
    }
    auto [p, s] = pow_sum(ratio, n - 1ULL);
    const i64 p2 = static_cast<i64>((static_cast<i128>(p) * (ratio % MOD)) % MOD);
    const i64 s2 = (s + p) % MOD;
    return {p2, s2};
}

std::pair<i64, u64> quotient_mod_from_exponents(u64 e_prev, u64 e_cur) {
    const u64 a = e_prev / e_cur;
    const u64 b = e_prev % e_cur;
    const i64 two_b = mod_pow(2LL, b);
    const i64 ratio = mod_pow(2LL, e_cur);
    const i64 geom = pow_sum(ratio, a).second;
    const i64 q_mod = static_cast<i64>((static_cast<i128>(two_b) * geom) % MOD);
    return {q_mod, b};
}

i64 inverse_mersenne_mod_mod(u64 u, u64 w) {
    u64 e_prev = u;
    u64 e_cur = w;
    i64 t_prev = 0;
    i64 t_cur = 1;
    int steps = 0;

    while (e_cur != 0ULL) {
        auto [q_mod, b] = quotient_mod_from_exponents(e_prev, e_cur);
        i64 t_next = static_cast<i64>((t_prev - static_cast<i128>(q_mod) * t_cur) % MOD);
        if (t_next < 0) {
            t_next += MOD;
        }
        e_prev = e_cur;
        e_cur = b;
        t_prev = t_cur;
        t_cur = t_next;
        ++steps;
    }

    const i64 a_mod = (mod_pow(2LL, u) - 1 + MOD) % MOD;
    if (steps % 2 == 1) {
        return t_prev;
    }
    return (a_mod + t_prev) % MOD;
}

i64 ext_gcd(i64 a, i64 b, i64& x, i64& y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    i64 x1 = 0;
    i64 y1 = 0;
    const i64 g = ext_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

u64 inverse_u64_mod(u64 a, u64 m) {
    i64 x = 0;
    i64 y = 0;
    const i64 g = ext_gcd(static_cast<i64>(a), static_cast<i64>(m), x, y);
    assert(g == 1);
    i64 r = x % static_cast<i64>(m);
    if (r < 0) {
        r += static_cast<i64>(m);
    }
    return static_cast<u64>(r);
}

u64 pow_u64(u64 base, int exp) {
    u64 result = 1ULL;
    for (int i = 0; i < exp; ++i) {
        result *= base;
    }
    return result;
}

u64 p_general(u64 a, u64 b) {
    assert(a <= b);
    assert(std::gcd(a, b) == 1ULL);

    const u64 inv = inverse_u64_mod(a % b, b);
    const u64 y = static_cast<u64>((static_cast<i128>(a) * inv - 1) / static_cast<i128>(b));
    const u64 s1 = inv + y - 1;
    const u64 s2 = (b - inv) + (a - y) - 1;
    return 2ULL * std::min(s1, s2);
}

i64 p_mersenne_from_exponents(u64 u, u64 v) {
    assert(u < v);
    assert(std::gcd(u, v) == 1ULL);

    const i64 a_mod = (mod_pow(2LL, u) - 1 + MOD) % MOD;
    const i64 b_mod = (mod_pow(2LL, v) - 1 + MOD) % MOD;
    const i64 inv_a_mod = mod_inv(a_mod);

    const u64 w = v % u;
    const i64 r_mod = inverse_mersenne_mod_mod(u, w);
    const u64 t = inverse_u64_mod(w, u);

    i64 d_mod = 0;
    i64 x_mod = 0;
    if (2ULL * t <= u) {
        d_mod = r_mod;
        i64 rhs = (static_cast<i64>((static_cast<i128>(b_mod) * d_mod) % MOD) - 1 + MOD) % MOD;
        x_mod = static_cast<i64>((static_cast<i128>(rhs) * inv_a_mod) % MOD);
    } else {
        d_mod = (a_mod - r_mod + MOD) % MOD;
        i64 rhs = (static_cast<i64>((static_cast<i128>(b_mod) * d_mod) % MOD) + 1) % MOD;
        x_mod = static_cast<i64>((static_cast<i128>(rhs) * inv_a_mod) % MOD);
    }

    const i64 m_mod = (x_mod + d_mod - 1 + MOD) % MOD;
    return static_cast<i64>((2LL * m_mod) % MOD);
}

std::vector<int> primes_below(int n) {
    std::vector<bool> is_prime(static_cast<std::size_t>(n), true);
    is_prime[0] = false;
    is_prime[1] = false;
    for (int p = 2; static_cast<i64>(p) * p < n; ++p) {
        if (!is_prime[p]) {
            continue;
        }
        for (int x = p * p; x < n; x += p) {
            is_prime[x] = false;
        }
    }
    std::vector<int> primes;
    for (int x = 2; x < n; ++x) {
        if (is_prime[x]) {
            primes.push_back(x);
        }
    }
    return primes;
}

}  // namespace

int main() {
    assert(p_general(3, 5) == 4ULL);
    assert(p_general(7, 31) == 20ULL);
    assert(p_general(1234, 4321) == 2780ULL);

    assert(p_mersenne_from_exponents(2, 5) == 20);
    assert(p_mersenne_from_exponents(3, 5) == 20);
    assert(p_mersenne_from_exponents(5, 8) == 164);
    for (u64 u = 2; u <= 9; ++u) {
        for (u64 v = u + 1; v <= 10; ++v) {
            if (std::gcd(u, v) != 1ULL) {
                continue;
            }
            const u64 a = (1ULL << u) - 1ULL;
            const u64 b = (1ULL << v) - 1ULL;
            assert(p_mersenne_from_exponents(u, v) == static_cast<i64>(p_general(a, b)));
        }
    }

    const std::vector<int> primes = primes_below(1000);
    std::vector<u64> p5;
    p5.reserve(primes.size());
    for (int p : primes) {
        p5.push_back(pow_u64(static_cast<u64>(p), 5));
    }

    i64 answer = 0;
    for (std::size_t i = 0; i < p5.size(); ++i) {
        for (std::size_t j = i + 1; j < p5.size(); ++j) {
            answer += p_mersenne_from_exponents(p5[i], p5[j]);
            if (answer >= MOD) {
                answer -= MOD;
            }
        }
    }

    std::cout << answer << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

import java.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler758 {
    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("Euler758.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(".euler758_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 Euler758 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("Euler758 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("Euler758 C++ bridge produced empty output.");
        }
        return parsed;
    }

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