Problem 263: An Engineers' Dream Come True

Project Euler Problem 263 Solution

EulerSolve provides an optimized solution for Project Euler Problem 263, An Engineers' Dream Come True, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The code searches for special integers called engineers' paradises . A candidate $n$ must satisfy two independent conditions: 1. four consecutive primes centered around $n$ must form a prime chain with gap 6, and 2. the five numbers $n-8, n-4, n, n+4, n+8$ must all be practical numbers. The solver finds the first $k$ such paradises and sums them. Mathematical Approach 1. The Prime-Window Filter The first filter is purely prime-based. The program keeps a sliding window of four consecutive primes $p_0,p_1,p_2,p_3$ and checks $$p_1-p_0=p_2-p_1=p_3-p_2=6.$$ If this holds, then the four primes are exactly $$p_0,\ p_0+6,\ p_0+12,\ p_0+18,$$ so the candidate center is $$n=p_0+9.$$ Equivalently, the prime pattern is centered at $$n-9,\ n-3,\ n+3,\ n+9.$$ This is the expensive combinatorial pattern the code looks for first, because only such windows can possibly produce a paradise. A concrete miniature example is the first prime window with this shape: $$251,\ 257,\ 263,\ 269.$$ It yields the center $$n=251+9=260.$$ The practical-number stage then checks $$252,\ 256,\ 260,\ 264,\ 268.$$ The first four values are practical, but 268 is not, so 260 is rejected. This shows the full logic of the search: the prime window is only a necessary condition, and the five local practical tests decide whether the candidate is really a paradise. 2....

Detailed mathematical approach

Problem Summary

A positive integer $n$ is practical if every integer from $1$ to $n$ is a sum of distinct divisors of $n$. The statement shows this for $6$: $1=1$, $2=2$, $3=1+2$, $4=1+3$, $5=2+3$ and $6=6$. A sexy pair is a pair of consecutive primes that differ by $6$; the first one is $(23,29)$. A triple-pair consists of three consecutive sexy pairs in which the second member of each pair is the first member of the next. A number $n$ is called an engineers' paradise if $(n-9,n-3)$, $(n-3,n+3)$ and $(n+3,n+9)$ form a triple-pair and the five numbers $n-8$, $n-4$, $n$, $n+4$ and $n+8$ are all practical. We must find the sum of the first four engineers' paradises. ^[1]

The whole solution in one sentence: list the primes in increasing order with a sieve, stop at every four consecutive primes that are $6$ apart, and test the five even numbers around their centre with a rule that decides from the prime factorisation alone whether a number is practical.

Mathematical Approach

1. Practical numbers by hand

The divisors of $12$ are $1,2,3,4,6,12$. Every number up to $12$ is a sum of distinct ones, for example $5=1+4$, $8=2+6$ and $11=1+4+6$, and this even works up to $1+2+3+4+6+12=28$, for instance $23=1+4+6+12$. The number $10$ is not practical: its smallest divisors $1$ and $2$ reach only $3$, and the next divisor $5$ is too large to help with $4$. There are $30$ practical numbers up to $100$, beginning 1, 2, 4, 6, 8, 12, 16, 18, 20, 24, 28, 30; they form OEIS A005153. ^[3]

2. The Stewart–Sierpiński criterion

Write $\sigma(m)$ for the sum of the divisors of $m$. It is multiplicative, $\sigma(ab)=\sigma(a)\sigma(b)$ for coprime $a$ and $b$, and $\sigma(p^a)=1+p+\dots+p^a$ for a prime power. ^[4] Let

$$n=p_1^{a_1}\cdots p_k^{a_k},\qquad p_1\lt\dots\lt p_k,$$

and let $A_i=p_1^{a_1}\cdots p_i^{a_i}$ be the product of the first $i$ prime powers, with $A_0=1$. Stewart (1954) and Sierpiński (1955) proved that $n\ge2$ is practical if and only if ^[2]

$$p_i\le\sigma(A_{i-1})+1\qquad(1\le i\le k).$$

For $i=1$ the condition reads $p_1\le\sigma(1)+1=2$, so $n$ must be even. For example $12=2^2\cdot3$ passes, with $2\le2$ and $3\le\sigma(4)+1=8$, while $10=2\cdot5$ fails, because $5\gt\sigma(2)+1=4$; and $4$ was exactly the number that $10$ could not form.

3. Why the condition is necessary

Suppose $p_i\gt\sigma(A_{i-1})+1$ for some $i$, and put $x=\sigma(A_{i-1})+1$, so that $x\lt p_i\le n$. A divisor of $n$ smaller than $p_i$ contains none of the primes $p_i,\dots,p_k$, so it divides $A_{i-1}$. A sum of distinct divisors equal to $x$ can only use divisors not larger than $x$, hence only divisors of $A_{i-1}$, and all of these together add up to $\sigma(A_{i-1})=x-1$. So $x$ cannot be formed and $n$ is not practical. For $268=2^2\cdot67$ the missing number is $x=\sigma(4)+1=8$: the divisors below $67$ are $1,2,4$, and they reach only $7$.

4. Why the condition is sufficient

We prove a little more. Call $m$ complete if every integer $k$ with $0\le k\le\sigma(m)$ is a sum of distinct divisors of $m$, the empty sum giving $0$; by Section 1, $12$ is complete. A complete number is practical, since $\sigma(m)\ge m$.

For example, the divisors $1,2,4$ of $4$ form every sum from $0$ to $7$. The divisors of $20$ are these three and their multiples $5,10,20$, so the sums of distinct divisors of $20$ are the numbers $u_0+5u_1$ with $0\le u_0,u_1\le7$. For fixed $u_1$ they fill the interval from $5u_1$ to $5u_1+7$, and since $5\le7+1$ these intervals overlap and cover everything from $0$ to $42=\sigma(20)$, so $20$ is complete.

Lemma. Let $m$ be complete, put $s=\sigma(m)$, and let $p$ be a prime with $p\nmid m$ and $p\le s+1$. Then $mp^a$ is complete for every $a\ge1$.

Proof. Every divisor of $mp^a$ is $dp^j$ for exactly one pair with $d\mid m$ and $0\le j\le a$. Grouping a sum of distinct divisors by the exponent $j$ gives $u_0+u_1p+\dots+u_ap^a$, where each $u_j$ is a sum of distinct divisors of $m$; since $m$ is complete, every choice of integers $u_j$ from $0$ to $s$ arises in this way. As $\sigma(mp^a)=s(1+p+\dots+p^a)$ by multiplicativity, it suffices to show, by induction on $a$, that every integer $k$ with $0\le k\le s(1+p+\dots+p^a)$ can be written as

$$k=u_0+u_1p+\dots+u_ap^a,\qquad 0\le u_j\le s.$$

For $a=0$ take $u_0=k$. For $a\ge1$ take $u_a=\min\bigl(s,\lfloor k/p^a\rfloor\bigr)$ and $r=k-u_ap^a\ge0$. If $u_a=\lfloor k/p^a\rfloor$, then $r\le p^a-1$, and $p^a-1=(p-1)(1+p+\dots+p^{a-1})$ is at most $s(1+p+\dots+p^{a-1})$ because $p-1\le s$. Otherwise $u_a=s$, and $r\le s(1+p+\dots+p^{a-1})$ follows from the bound on $k$. In both cases the induction hypothesis writes $r$ with digits $u_0,\dots,u_{a-1}$. $\square$

Now $A_0=1$ is complete, and if $A_{i-1}$ is complete and $p_i\le\sigma(A_{i-1})+1$, the lemma makes $A_i=A_{i-1}p_i^{a_i}$ complete. After $k$ steps $n=A_k$ is complete, hence practical. With Section 3 this proves the criterion, and it shows that a practical number forms every integer up to $\sigma(n)$, which is how OEIS A005153 defines them. ^[3]

5. The criterion on five neighbours

The first triple-pair, found in Section 6, has the centre $n=260$. The criterion decides its five neighbours as follows.

$252=2^2\cdot3^2\cdot7$: $3\le\sigma(4)+1=8$ and $7\le\sigma(36)+1=92$, practical.
$256=2^8$: no prime after $2$, practical.
$260=2^2\cdot5\cdot13$: $5\le8$ and $13\le\sigma(20)+1=43$, practical.
$264=2^3\cdot3\cdot11$: $3\le\sigma(8)+1=16$ and $11\le\sigma(24)+1=61$, practical.
$268=2^2\cdot67$: $67\gt8$, not practical.

So $260$ is not a paradise; it fails only at its last neighbour.

6. Triple-pairs are four consecutive primes

The three sexy pairs share their members, so a triple-pair is the same as four consecutive primes $n-9$, $n-3$, $n+3$, $n+9$: these four numbers are prime, and the odd numbers $n-7$, $n-5$, $n-1$, $n+1$, $n+5$, $n+7$ between them are composite. The word consecutive matters: $5,11,17,23$ are all prime, a sexy prime quadruplet in the terminology of MathWorld, ^[5] but $7$, $13$ and $19$ lie between them. The first four consecutive primes with the gaps $6,6,6$ are $251,257,263,269$, and the next runs start at $1741$, $3301$, $5101$, $5381$, $6311$, …, the sequence A033451 of the OEIS. ^[6] There are $7$ triple-pairs with $n\le10^4$, $189$ with $n\le10^6$, $6698$ with $n\le10^8$ and $75\,291$ with $n\le2\cdot10^9$.

7. What the congruences force

All four primes are odd and at least $251$, so $n$ is even. They are all congruent to $n$ modulo $3$, so $3\nmid n$, and modulo $5$ they are $n+1$, $n+2$, $n+3$, $n+4$, none of them $0$, so $5\mid n$. Hence $n\equiv10$ or $n\equiv20\pmod{30}$. The practical condition adds a factor $4$: if $n\equiv2\pmod4$, then $n/2$ is odd, so by the criterion the second prime factor of $n$ is at most $\sigma(2)+1=4$, that is $3\mid n$, which is excluded. Therefore every paradise satisfies $n\equiv20\pmod{60}$ or $n\equiv40\pmod{60}$. The five neighbours are congruent to $n$ modulo $4$, so for $n\equiv2\pmod4$ the argument applies to each of them, and as only two of them are divisible by $3$, at most two are practical; this concerns $37\,497$ of the $75\,291$ triple-pairs with $n\le2\cdot10^9$. The programs do not use these congruences.

8. Listing the primes with a segmented sieve

The programs search the centres $n\le L=2\cdot10^9$ and so need the primes up to $L+9$. A plain sieve of Eratosthenes first lists the $4648$ base primes up to $\lfloor\sqrt{L+9}\rfloor+1=44\,722$. The numbers from $3$ to $L+9$ are then cut into segments of $8\,000\,000$ integers with one byte per odd number. For every odd base prime $p$ with $p^2$ not above the end of a segment, its odd multiples are crossed off in steps of $2p$, from the first odd multiple in the segment or from $p^2$, whichever is larger. ^[7] Only multiples that are at least $p^2$ are marked, so no prime is crossed off, while an odd composite $x$ is crossed off as a multiple of its smallest prime factor $q$, for which $q^2\le x$; exactly the odd primes stay unmarked.

9. Factoring the five neighbours

Each neighbour $m\le L+8$ is factored by trial division with the base primes in increasing order, dividing out each prime completely and stopping as soon as $p^2$ exceeds the remaining cofactor. ^[8] A remaining cofactor larger than $1$ is prime: otherwise it would have a prime factor not above its square root, which would have been divided out already, because the base primes reach $44\,722$ and $44\,722^2\gt L+8$. The primes thus come out in increasing order, as the criterion needs. The test keeps the divisor sum $\sigma(A_{i-1})$, compares the next prime with this sum plus $1$, and then multiplies the sum by $\sigma(p_i^{a_i})$. These sums are at most $\sigma(m)$, and

$$\begin{aligned}\sigma(m)&=\sum_{d\mid m}\frac md\le m\sum_{j=1}^{m}\frac1j\\&\le m(1+\ln m),\end{aligned}$$

which is below $5\cdot10^{10}$ for $m\le L+8$, far inside the range of 64-bit integers.

10. The first paradise

The first paradise, which the C++ program uses as a checkpoint, is $n=219\,869\,980$, the centre of the $12\,629$th triple-pair. Its primes are $219\,869\,971$, $219\,869\,977$, $219\,869\,983$ and $219\,869\,989$, and every prime factor of a neighbour after $2$ is at most its bound $\sigma(A_{i-1})+1$:

$219\,869\,972$${}=2^2\cdot7\cdot23\cdot223\cdot1531$: $7\le8$, $23\le57$, $223\le1345$, $1531\le301057$;
$219\,869\,976$${}=2^3\cdot3\cdot17\cdot19\cdot113\cdot251$: $3\le16$, $17\le61$, $19\le1081$, $113\le21601$, $251\le2462401$;
$219\,869\,980$${}=2^2\cdot5\cdot11\cdot31\cdot103\cdot313$: $5\le8$, $11\le43$, $31\le505$, $103\le16129$, $313\le1677313$;
$219\,869\,984$${}=2^5\cdot37\cdot233\cdot797$: $37\le64$, $233\le2395$, $797\le560197$;
$219\,869\,988$${}=2^2\cdot3\cdot13\cdot83\cdot16981$: $3\le8$, $13\le29$, $83\le393$, $16981\le32929$.

Such neighbours are rare: of the $75\,291$ triple-pairs with $n\le2\cdot10^9$, exactly $36\,622$, $26\,169$, $10\,422$, $1950$, $124$ and $4$ have $0$, $1$, $2$, $3$, $4$ and $5$ practical neighbours, and the four with five are the first four paradises, so the limit $L$ suffices.

How the Code Works

The C++ program ^[10] starts in main. parse_arguments fills an Options record with the defaults $L=2\cdot10^9$, four paradises and segments of $8\,000\,000$ numbers; its options include --limit=, --target=, --segment= and --skip-checkpoints. isqrt_u64 computes $\lfloor\sqrt{L+9}\rfloor$ exactly, and sieve_primes lists the base primes on a byte array (Section 8).

Unless the checkpoints are skipped, check_practical_criterion compares is_practical with is_practical_bruteforce for $1\le n\le400$. The brute force lists the divisors with generate_divisors_rec and marks the reachable sums in a byte table reachable of length $n+1$, starting from reachable[0]; for each divisor $d$ the index $s$ runs down from $n$ to $d$, so that $d$ is used at most once. This is the standard dynamic program for subset sums, ^[9] and every value from $1$ to $n$ must be reachable.

find_paradises keeps the last four primes in a std::deque called window. Its helper process_prime appends a prime, drops the oldest one beyond four and, when the three gaps are all $6$, takes $n$ as the oldest prime plus $9$. If $n\le L$, is_engineers_paradise_candidate runs is_practical on the five neighbours and stops at the first failure; a paradise is stored, added to the sum and printed. After the prime $2$, sieve_segment(low, high) supplies the primes of consecutive segments from $3$ upward in increasing order, so the window holds consecutive primes also across segment boundaries, and the scan ends when the requested number of paradises is found.

is_practical implements Sections 2 and 9: $1$ is practical and odd numbers greater than $1$ are not; otherwise factorize returns the factorisation in increasing order, its first prime must be $2$, and sigma_prefix starts at $\sigma(2^a)$, computed by sigma_prime_power in unsigned 128-bit arithmetic. A further prime $p$ fails the test if it exceeds sigma_prefix + 1 and otherwise multiplies sigma_prefix by $\sigma(p^a)$. Finally check_known_first_two requires the first two paradises to be $219\,869\,980$ and $312\,501\,820$ when the limit and the target allow it. If fewer paradises than requested were found, the program asks for a larger limit; otherwise it prints the sum and the elapsed time.

Python and Java follow the same steps with the limit, the target $4$ and the segment length fixed, without options or checkpoints, and print the sum without checking that four paradises were found. Python uses bytearray sieves, steps through the odd-number index by $p$, a step of $2p$ in value, and has unlimited integers. Java returns the primes of a segment as a List<Long> of boxed integers and computes $\sigma$ in a 64-bit long, which is safe by the bound of Section 9; it ignores any exception raised while the primes of a batch of segments are fetched or scanned, which would silently skip primes.

Complexity and Verification

Let $X\le L+9$ be the end of the last segment that is sieved, which may lie a little beyond the point where the scan stops. The sieve crosses off every odd multiple $jp\ge p^2$ up to $X$ of every odd base prime once, about $\tfrac X2\sum_{p\le\sqrt X}\tfrac1p=O(X\log\log X)$ operations, ^[7] plus a loop over the $4648$ base primes per segment; the full range up to $L+9$ would take $2\,083\,999\,874$ crossings. The window costs $O(1)$ per prime, and a practical test factors its number completely with at most $4648$ trial divisors; over all $75\,291$ triple-pairs with $n\le L$, the early exit leaves only $87\,362$ calls of is_practical with $36\,205\,233$ trial divisors in total. A segment needs $4\,000\,000$ bytes for its odd numbers and a list of at most $539\,776$ primes; C++ and Java hold several segments at a time. On an Apple M4 Max (clang++ -O2, OpenJDK 17, Python 3.11) the C++, Java and Python programs take about $0.7$, $1.1$ and $39$ seconds of wall-clock time. The peak resident memory is about $112$ MB for C++, $73$ MB for Python and $3.6$ GB for Java under the default heap settings, although Java also finishes with a $256$ MB heap.

Beyond the two checkpoints of the C++ program, this article checked the following. The criterion agreed with the subset-sum definition for every $n\le10^5$, which includes $11\,751$ practical numbers, and every practical $n\le10^4$ is complete. An independent program with a bit-packed sieve and its own practical test produced the counts of Sections 6, 7 and 10 and the same first four paradises, and its sum equals the sums printed by the C++, Java and Python programs; a separate Python scan reproduced the first twenty starting primes of A033451. ^[6] For each paradise, trial division confirmed the primes and the composites between them, and a subset-sum bitset confirmed from the definition alone that all five neighbours are practical. These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts used above. The proof of the criterion, the congruences of Section 7 and all counts and timings are derived or computed in this article.

Project Euler 263 — An Engineers' Dream Come True. The official statement defines practical numbers, sexy pairs, triple-pairs and engineers' paradises as in the summary above, with the examples 6 and (23, 29), and asks for the sum of the first four engineers' paradises.
Wikipedia — Practical number. A practical number is a positive integer n such that all smaller positive integers are sums of distinct divisors of n. Stewart (1954) and Sierpiński (1955) showed that an integer n > 1 with prime factorisation n = p₁^α₁ ⋯ p_k^α_k, p₁ < ⋯ < p_k, is practical if and only if p₁ = 2 and p_i ≤ 1 + σ(p₁^α₁ ⋯ p_{i−1}^α_{i−1}) for 2 ≤ i ≤ k.
OEIS A005153 — Practical numbers. The entry defines the practical numbers as the positive integers m such that every k ≤ σ(m) is a sum of distinct divisors of m, and lists 1, 2, 4, 6, 8, 12, 16, 18, 20, 24, ….
Wikipedia — Divisor function. The sum-of-divisors function σ is multiplicative, and for a prime power σ(pᵃ) = 1 + p + p² + ⋯ + pᵃ.
MathWorld — Sexy Primes. Sexy quadruplets are primes p, p + 6, p + 12, p + 18, not necessarily consecutive, such as (5, 11, 17, 23), (11, 17, 23, 29) and (251, 257, 263, 269).
OEIS A033451 — Initial prime in set of 4 consecutive primes with common difference 6. The entry lists 251, 1741, 3301, 5101, 5381, 6311, …, the smallest primes of the triple-pairs.
Wikipedia — Sieve of Eratosthenes. Crossing off may start at p², an odd-only variant steps by 2p, a segmented sieve treats one segment at a time with the primes up to the square root, and the sieve takes O(n log log n) operations.
Wikipedia — Trial division. Trial factors need go no further than the square root of n, and it suffices to try primes, each repeatedly until it no longer divides.
Wikipedia — Subset sum problem. Dynamic programming over the reachable sums solves the problem in pseudo-polynomial time.
C++ — Euler263.cpp. The linked immutable revision contains is_practical, sieve_segment, find_paradises and both checkpoints; it is the implementation record for this article, and the file shown on this page adds one output line at the end.

Problem 263 source code

C++

#include <algorithm>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <deque>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <utility>
#include <vector>

namespace {

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

constexpr u64 kDefaultLimit = 2'000'000'000ULL;
constexpr unsigned kDefaultTargetCount = 4U;
constexpr u64 kDefaultSegmentSpan = 8'000'000ULL;

constexpr int kPracticalCrossCheckMax = 400;
constexpr u64 kKnownParadise1 = 219'869'980ULL;
constexpr u64 kKnownParadise2 = 312'501'820ULL;

struct Options {
    u64 limit = kDefaultLimit;
    unsigned target_count = kDefaultTargetCount;
    u64 segment_span = kDefaultSegmentSpan;
    bool allow_multithreading = true;
    bool run_checkpoints = true;
    unsigned requested_threads = 0;
};

struct SearchResult {
    std::vector<u64> paradises;
    u64 sum = 0ULL;
};

bool parse_u64_after_prefix(const std::string& arg, const char* prefix, u64& value) {
    const std::string p(prefix);
    if (arg.rfind(p, 0) != 0) {
        return false;
    }

    const std::string tail = arg.substr(p.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            return false;
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg, const char* prefix, unsigned& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) {
        return false;
    }
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
        return false;
    }
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }

        u64 parsed_u64 = 0ULL;
        if (parse_u64_after_prefix(arg, "--limit=", parsed_u64)) {
            options.limit = parsed_u64;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--segment=", parsed_u64)) {
            options.segment_span = parsed_u64;
            continue;
        }

        unsigned parsed_unsigned = 0U;
        if (parse_unsigned_after_prefix(arg, "--target=", parsed_unsigned)) {
            options.target_count = parsed_unsigned;
            continue;
        }
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_unsigned)) {
            options.requested_threads = parsed_unsigned;
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    if (options.limit < 30ULL) {
        std::cerr << "--limit must be at least 30.\n";
        return false;
    }
    if (options.segment_span < 10ULL) {
        std::cerr << "--segment must be at least 10.\n";
        return false;
    }
    if (options.target_count == 0U) {
        std::cerr << "--target must be at least 1.\n";
        return false;
    }

    return true;
}

u64 isqrt_u64(u64 n) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(n)));
    while ((r + 1ULL) <= n / (r + 1ULL)) {
        ++r;
    }
    while (r > n / r) {
        --r;
    }
    return r;
}

std::vector<int> sieve_primes(int limit) {
    if (limit < 2) {
        return {};
    }

    std::vector<std::uint8_t> is_prime(static_cast<std::size_t>(limit) + 1ULL, 1U);
    is_prime[0] = 0U;
    is_prime[1] = 0U;

    const int root = static_cast<int>(std::sqrt(static_cast<long double>(limit)));
    for (int p = 2; p <= root; ++p) {
        if (is_prime[static_cast<std::size_t>(p)] == 0U) {
            continue;
        }
        for (std::int64_t m = static_cast<std::int64_t>(p) * p; m <= limit; m += p) {
            is_prime[static_cast<std::size_t>(m)] = 0U;
        }
    }

    std::vector<int> primes;
    for (int p = 2; p <= limit; ++p) {
        if (is_prime[static_cast<std::size_t>(p)] != 0U) {
            primes.push_back(p);
        }
    }
    return primes;
}

std::vector<std::pair<u64, int>> factorize(u64 n, const std::vector<int>& primes) {
    std::vector<std::pair<u64, int>> factors;
    if (n <= 1ULL) {
        return factors;
    }

    u64 m = n;
    for (const int p_int : primes) {
        const u64 p = static_cast<u64>(p_int);
        if (p > m / p) {
            break;
        }
        if ((m % p) != 0ULL) {
            continue;
        }

        int exponent = 0;
        while ((m % p) == 0ULL) {
            m /= p;
            ++exponent;
        }
        factors.emplace_back(p, exponent);
    }

    if (m > 1ULL) {
        factors.emplace_back(m, 1);
    }

    return factors;
}

u128 sigma_prime_power(u64 prime, int exponent) {
    u128 sum = 1;
    u128 power = 1;
    for (int i = 0; i < exponent; ++i) {
        power *= static_cast<u128>(prime);
        sum += power;
    }
    return sum;
}

bool is_practical(u64 n, const std::vector<int>& factor_primes) {
    if (n == 1ULL) {
        return true;
    }
    if ((n & 1ULL) == 1ULL) {
        return false;
    }

    const std::vector<std::pair<u64, int>> factors = factorize(n, factor_primes);
    if (factors.empty() || factors[0].first != 2ULL) {
        return false;
    }

    u128 sigma_prefix = sigma_prime_power(factors[0].first, factors[0].second);
    for (std::size_t i = 1; i < factors.size(); ++i) {
        const u64 p = factors[i].first;
        if (static_cast<u128>(p) > sigma_prefix + 1U) {
            return false;
        }
        sigma_prefix *= sigma_prime_power(factors[i].first, factors[i].second);
    }

    return true;
}

void generate_divisors_rec(const std::vector<std::pair<u64, int>>& factors,
                           std::size_t index,
                           u64 current,
                           std::vector<u64>& divisors) {
    if (index == factors.size()) {
        divisors.push_back(current);
        return;
    }

    const u64 p = factors[index].first;
    const int exponent = factors[index].second;

    u64 value = 1ULL;
    for (int e = 0; e <= exponent; ++e) {
        generate_divisors_rec(factors, index + 1ULL, current * value, divisors);
        if (e != exponent) {
            value *= p;
        }
    }
}

bool is_practical_bruteforce(u64 n, const std::vector<int>& factor_primes) {
    if (n == 0ULL) {
        return false;
    }
    if (n == 1ULL) {
        return true;
    }

    const std::vector<std::pair<u64, int>> factors = factorize(n, factor_primes);
    std::vector<u64> divisors;
    divisors.reserve(128);
    generate_divisors_rec(factors, 0ULL, 1ULL, divisors);
    std::sort(divisors.begin(), divisors.end());

    std::vector<std::uint8_t> reachable(static_cast<std::size_t>(n) + 1ULL, 0U);
    reachable[0] = 1U;
    for (const u64 d : divisors) {
        for (u64 s = n; s >= d; --s) {
            if (reachable[static_cast<std::size_t>(s - d)] != 0U) {
                reachable[static_cast<std::size_t>(s)] = 1U;
            }
            if (s == d) {
                break;
            }
        }
    }

    for (u64 k = 1ULL; k <= n; ++k) {
        if (reachable[static_cast<std::size_t>(k)] == 0U) {
            return false;
        }
    }
    return true;
}

bool check_practical_criterion(const std::vector<int>& factor_primes) {
    for (int n = 1; n <= kPracticalCrossCheckMax; ++n) {
        const bool criterion = is_practical(static_cast<u64>(n), factor_primes);
        const bool brute = is_practical_bruteforce(static_cast<u64>(n), factor_primes);
        if (criterion != brute) {
            std::cerr << "Practical criterion mismatch at n=" << n
                      << " criterion=" << criterion
                      << " brute=" << brute << '\n';
            return false;
        }
    }
    return true;
}

unsigned choose_thread_count(bool allow_multithreading,
                             unsigned requested_threads,
                             std::size_t workload) {
    if (!allow_multithreading || workload < 2ULL) {
        return 1U;
    }

    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) {
            threads = 1U;
        }
    }

    return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload)));
}

std::vector<u64> sieve_segment(u64 low, u64 high, const std::vector<int>& base_primes) {
    if (high < 2ULL || low > high) {
        return {};
    }

    const u64 odd_low = std::max<u64>(3ULL, low | 1ULL);
    if (odd_low > high) {
        return {};
    }

    const std::size_t odd_count =
        static_cast<std::size_t>((high - odd_low) / 2ULL + 1ULL);
    std::vector<std::uint8_t> is_prime(odd_count, 1U);

    for (const int p_int : base_primes) {
        const u64 p = static_cast<u64>(p_int);
        if (p == 2ULL) {
            continue;
        }
        if (p > high / p) {
            break;
        }

        u64 start = odd_low;
        const u64 rem = start % p;
        if (rem != 0ULL) {
            start += (p - rem);
        }
        if (start < p * p) {
            start = p * p;
        }
        if ((start & 1ULL) == 0ULL) {
            start += p;
        }

        for (u64 m = start; m <= high; m += 2ULL * p) {
            const std::size_t idx = static_cast<std::size_t>((m - odd_low) / 2ULL);
            is_prime[idx] = 0U;
        }
    }

    std::vector<u64> primes;
    primes.reserve(odd_count / 8ULL + 16ULL);
    for (std::size_t i = 0; i < odd_count; ++i) {
        if (is_prime[i] != 0U) {
            primes.push_back(odd_low + 2ULL * static_cast<u64>(i));
        }
    }
    return primes;
}

bool is_engineers_paradise_candidate(u64 n, const std::vector<int>& factor_primes) {
    if (n < 9ULL) {
        return false;
    }

    static constexpr int kOffsets[5] = {-8, -4, 0, 4, 8};
    for (const int offset : kOffsets) {
        const u64 value = static_cast<u64>(static_cast<std::int64_t>(n) + offset);
        if (!is_practical(value, factor_primes)) {
            return false;
        }
    }
    return true;
}

SearchResult find_paradises(u64 limit,
                            unsigned target_count,
                            u64 segment_span,
                            bool allow_multithreading,
                            unsigned requested_threads,
                            const std::vector<int>& base_primes,
                            const std::vector<int>& factor_primes) {
    SearchResult result;
    result.paradises.reserve(target_count);

    const u64 prime_limit = limit + 9ULL;
    u64 segment_low = 3ULL;
    const std::size_t total_segments =
        static_cast<std::size_t>((prime_limit >= segment_low)
                                     ? ((prime_limit - segment_low) / segment_span + 1ULL)
                                     : 0ULL);

    const unsigned threads =
        choose_thread_count(allow_multithreading, requested_threads, total_segments);

    std::deque<u64> window;
    auto process_prime = [&](u64 p) {
        window.push_back(p);
        if (window.size() > 4ULL) {
            window.pop_front();
        }
        if (window.size() != 4ULL) {
            return;
        }
        if ((window[1] - window[0]) != 6ULL ||
            (window[2] - window[1]) != 6ULL ||
            (window[3] - window[2]) != 6ULL) {
            return;
        }

        const u64 n = window[0] + 9ULL;
        if (n > limit) {
            return;
        }
        if (!is_engineers_paradise_candidate(n, factor_primes)) {
            return;
        }

        result.paradises.push_back(n);
        result.sum += n;
        std::cout << "Found paradise #" << result.paradises.size() << ": " << n << '\n';
    };

    if (prime_limit >= 2ULL) {
        process_prime(2ULL);
    }

    while (segment_low <= prime_limit && result.paradises.size() < target_count) {
        std::vector<std::pair<u64, u64>> batch_ranges;
        batch_ranges.reserve(threads);
        for (unsigned t = 0; t < threads && segment_low <= prime_limit; ++t) {
            const u64 low = segment_low;
            const u64 high = std::min<u64>(prime_limit, low + segment_span - 1ULL);
            batch_ranges.emplace_back(low, high);
            segment_low = high + 1ULL;
        }

        std::vector<std::vector<u64>> batch_primes(batch_ranges.size());
        if (batch_ranges.size() == 1ULL) {
            batch_primes[0] = sieve_segment(batch_ranges[0].first, batch_ranges[0].second, base_primes);
        } else {
            std::vector<std::thread> workers;
            workers.reserve(batch_ranges.size());
            for (std::size_t i = 0; i < batch_ranges.size(); ++i) {
                workers.emplace_back([&, i]() {
                    batch_primes[i] =
                        sieve_segment(batch_ranges[i].first, batch_ranges[i].second, base_primes);
                });
            }
            for (std::thread& worker : workers) {
                worker.join();
            }
        }

        for (const std::vector<u64>& primes : batch_primes) {
            for (const u64 p : primes) {
                process_prime(p);
                if (result.paradises.size() >= target_count) {
                    break;
                }
            }
            if (result.paradises.size() >= target_count) {
                break;
            }
        }
    }

    return result;
}

bool check_known_first_two(const SearchResult& result, const Options& options) {
    if (options.limit < kKnownParadise2 || options.target_count < 2U) {
        return true;
    }
    if (result.paradises.size() < 2ULL) {
        std::cerr << "Known-value checkpoint failed: fewer than two paradises found.\n";
        return false;
    }
    if (result.paradises[0] != kKnownParadise1 || result.paradises[1] != kKnownParadise2) {
        std::cerr << "Known-value checkpoint failed.\n";
        std::cerr << "Expected first two paradises: "
                  << kKnownParadise1 << ", " << kKnownParadise2 << '\n';
        std::cerr << "Observed: "
                  << result.paradises[0] << ", " << result.paradises[1] << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }

    const auto start_time = std::chrono::steady_clock::now();

    const u64 sieve_max = options.limit + 9ULL;
    const int root = static_cast<int>(isqrt_u64(sieve_max)) + 1;
    const std::vector<int> base_primes = sieve_primes(root);

    if (options.run_checkpoints) {
        if (!check_practical_criterion(base_primes)) {
            return 1;
        }
        std::cout << "Checkpoint passed: practical criterion matches brute force up to "
                  << kPracticalCrossCheckMax << ".\n";
    }

    SearchResult result = find_paradises(options.limit,
                                         options.target_count,
                                         options.segment_span,
                                         options.allow_multithreading,
                                         options.requested_threads,
                                         base_primes,
                                         base_primes);

    if (options.run_checkpoints) {
        if (!check_known_first_two(result, options)) {
            return 1;
        }
        if (options.limit >= kKnownParadise2 && options.target_count >= 2U) {
            std::cout << "Checkpoint passed: first two paradises match known values.\n";
        }
    }

    if (result.paradises.size() < options.target_count) {
        std::cerr << "Only found " << result.paradises.size()
                  << " paradises up to limit " << options.limit
                  << ". Increase --limit.\n";
        return 1;
    }

    const auto end_time = std::chrono::steady_clock::now();
    const std::chrono::duration<double> elapsed = end_time - start_time;

    std::cout << "Sum of first " << options.target_count
              << " engineers' paradises: " << result.sum << '\n';
    std::cout << "Elapsed: " << elapsed.count() << " seconds\n";
    std::cout << "Answer: " << result.sum << '\n';

    return 0;
}

Python

import math

def solve():
    LIMIT = 2_000_000_000
    TARGET_COUNT = 4

    def isqrt(n):
        return math.isqrt(n)

    def sieve_primes(limit):
        if limit < 2:
            return []
        is_p = bytearray(b'\x01' * (limit + 1))
        is_p[0] = 0
        is_p[1] = 0
        for p in range(2, isqrt(limit) + 1):
            if is_p[p]:
                is_p[p*p::p] = bytearray(len(is_p[p*p::p]))
        return [p for p in range(2, limit + 1) if is_p[p]]

    def factorize(n, primes):
        factors = []
        if n <= 1:
            return factors
        m = n
        for p in primes:
            if p * p > m:
                break
            if m % p != 0:
                continue
            e = 0
            while m % p == 0:
                m //= p
                e += 1
            factors.append((p, e))
        if m > 1:
            factors.append((m, 1))
        return factors

    def sigma_pp(prime, exp):
        s = 1
        pw = 1
        for _ in range(exp):
            pw *= prime
            s += pw
        return s

    def is_practical(n, primes):
        if n == 1:
            return True
        if n & 1:
            return False
        factors = factorize(n, primes)
        if not factors or factors[0][0] != 2:
            return False
        sigma_prefix = sigma_pp(factors[0][0], factors[0][1])
        for i in range(1, len(factors)):
            p = factors[i][0]
            if p > sigma_prefix + 1:
                return False
            sigma_prefix *= sigma_pp(factors[i][0], factors[i][1])
        return True

    def is_paradise_candidate(n, primes):
        if n < 9:
            return False
        for off in [-8, -4, 0, 4, 8]:
            if not is_practical(n + off, primes):
                return False
        return True

    sieve_max = LIMIT + 9
    root = isqrt(sieve_max) + 1
    base_primes = sieve_primes(root)

    def sieve_segment(low, high):
        if high < 2 or low > high:
            return []
        odd_low = max(3, low | 1)
        if odd_low > high:
            return []
        odd_count = (high - odd_low) // 2 + 1
        is_p = bytearray(b'\x01' * odd_count)
        for p in base_primes:
            if p == 2:
                continue
            if p * p > high:
                break
            start = ((odd_low + p - 1) // p) * p
            if start < p * p:
                start = p * p
            if start % 2 == 0:
                start += p
            idx = (start - odd_low) // 2
            while idx < odd_count:
                is_p[idx] = 0
                idx += p
        return [odd_low + 2 * i for i in range(odd_count) if is_p[i]]

    prime_limit = LIMIT + 9
    window = []
    paradises = []
    total_sum = 0

    def process_prime(p):
        nonlocal total_sum
        window.append(p)
        if len(window) > 4:
            window.pop(0)
        if len(window) != 4:
            return
        if (window[1]-window[0] != 6 or window[2]-window[1] != 6 or window[3]-window[2] != 6):
            return
        n = window[0] + 9
        if n > LIMIT:
            return
        if not is_paradise_candidate(n, base_primes):
            return
        paradises.append(n)
        total_sum += n

    process_prime(2)
    seg_size = 8_000_000
    seg_low = 3
    while seg_low <= prime_limit and len(paradises) < TARGET_COUNT:
        seg_high = min(prime_limit, seg_low + seg_size - 1)
        primes = sieve_segment(seg_low, seg_high)
        for p in primes:
            process_prime(p)
            if len(paradises) >= TARGET_COUNT:
                break
        seg_low = seg_high + 1

    return str(total_sum)

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

Java

import java.util.*;
import java.util.concurrent.*;

public class Euler263 {
    static int isqrt(long n) {
        if (n < 0)
            return 0;
        return (int) Math.sqrt(n);
    }

    static List<Integer> sievePrimes(int limit) {
        if (limit < 2)
            return new ArrayList<>();
        boolean[] isPrime = new boolean[limit + 1];
        Arrays.fill(isPrime, true);
        isPrime[0] = isPrime[1] = false;
        for (int p = 2; p <= isqrt(limit); ++p) {
            if (isPrime[p]) {
                for (int i = p * p; i <= limit; i += p) {
                    isPrime[i] = false;
                }
            }
        }
        List<Integer> primes = new ArrayList<>();
        for (int p = 2; p <= limit; ++p) {
            if (isPrime[p])
                primes.add(p);
        }
        return primes;
    }

    static class Factor {
        long prime;
        int exp;

        Factor(long p, int e) {
            prime = p;
            exp = e;
        }
    }

    static List<Factor> factorize(long n, List<Integer> primes) {
        List<Factor> factors = new ArrayList<>();
        long m = n;
        for (int p : primes) {
            if ((long) p * p > m)
                break;
            if (m % p == 0) {
                int exp = 0;
                while (m % p == 0) {
                    m /= p;
                    exp++;
                }
                factors.add(new Factor(p, exp));
            }
        }
        if (m > 1) {
            factors.add(new Factor(m, 1));
        }
        return factors;
    }

    static long sigmaPrimePower(long p, int exp) {
        long s = 1;
        long power = 1;
        for (int i = 0; i < exp; ++i) {
            power *= p;
            s += power;
        }
        return s;
    }

    static boolean isPractical(long n, List<Integer> basePrimes) {
        if (n == 1)
            return true;
        if (n % 2 != 0)
            return false;
        List<Factor> factors = factorize(n, basePrimes);
        if (factors.isEmpty() || factors.get(0).prime != 2)
            return false;
        long sigmaPrefix = sigmaPrimePower(factors.get(0).prime, factors.get(0).exp);
        for (int i = 1; i < factors.size(); ++i) {
            long p = factors.get(i).prime;
            int exp = factors.get(i).exp;
            if (p > sigmaPrefix + 1)
                return false;
            sigmaPrefix *= sigmaPrimePower(p, exp);
        }
        return true;
    }

    static boolean isEngineersParadiseCandidate(long n, List<Integer> basePrimes) {
        if (n < 9)
            return false;
        int[] offsets = { -8, -4, 0, 4, 8 };
        for (int offset : offsets) {
            if (!isPractical(n + offset, basePrimes))
                return false;
        }
        return true;
    }

    static List<Long> sieveSegment(long low, long high, List<Integer> basePrimes) {
        if (high < 2 || low > high)
            return new ArrayList<>();
        long oddLow = Math.max(3, low | 1);
        if (oddLow > high)
            return new ArrayList<>();
        int oddCount = (int) ((high - oddLow) / 2 + 1);
        byte[] isPrime = new byte[oddCount];
        Arrays.fill(isPrime, (byte) 1);

        for (int p : basePrimes) {
            if (p == 2)
                continue;
            if ((long) p * p > high)
                break;
            long start = oddLow;
            long rem = start % p;
            if (rem != 0)
                start += (p - rem);
            if (start < (long) p * p)
                start = (long) p * p;
            if (start % 2 == 0)
                start += p;

            for (long m = start; m <= high; m += 2L * p) {
                int idx = (int) ((m - oddLow) / 2);
                isPrime[idx] = 0;
            }
        }

        List<Long> primes = new ArrayList<>();
        for (int i = 0; i < oddCount; ++i) {
            if (isPrime[i] == 1) {
                primes.add(oddLow + 2L * i);
            }
        }
        return primes;
    }

    public static String solve() {
        long limit = 2000000000L;
        int targetCount = 4;
        long segmentSpan = 8000000L;

        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());
        long primeLimit = limit + 9;
        int root = isqrt(primeLimit) + 1;
        List<Integer> basePrimes = sievePrimes(root);

        List<Long> paradises = new CopyOnWriteArrayList<>();
        List<Long> window = new ArrayList<>();

        java.util.function.Consumer<Long> processPrime = (p) -> {
            window.add(p);
            if (window.size() > 4) {
                window.remove(0);
            }
            if (window.size() != 4)
                return;
            if (window.get(1) - window.get(0) != 6 ||
                    window.get(2) - window.get(1) != 6 ||
                    window.get(3) - window.get(2) != 6) {
                return;
            }

            long n = window.get(0) + 9;
            if (n > limit)
                return;
            if (isEngineersParadiseCandidate(n, basePrimes)) {
                paradises.add(n);
            }
        };

        processPrime.accept(2L);
        long segmentLow = 3;

        ExecutorService executor = Executors.newFixedThreadPool(threads);

        while (segmentLow <= primeLimit && paradises.size() < targetCount) {
            List<long[]> batchRanges = new ArrayList<>();
            for (int t = 0; t < threads; ++t) {
                if (segmentLow > primeLimit)
                    break;
                long high = Math.min(primeLimit, segmentLow + segmentSpan - 1);
                batchRanges.add(new long[] { segmentLow, high });
                segmentLow = high + 1;
            }

            if (batchRanges.size() == 1) {
                List<Long> primes = sieveSegment(batchRanges.get(0)[0], batchRanges.get(0)[1], basePrimes);
                for (long p : primes) {
                    processPrime.accept(p);
                    if (paradises.size() >= targetCount)
                        break;
                }
            } else {
                List<Future<List<Long>>> futures = new ArrayList<>();
                for (long[] range : batchRanges) {
                    futures.add(executor.submit(() -> sieveSegment(range[0], range[1], basePrimes)));
                }

                for (Future<List<Long>> future : futures) {
                    try {
                        List<Long> primes = future.get();
                        for (long p : primes) {
                            processPrime.accept(p);
                            if (paradises.size() >= targetCount)
                                break;
                        }
                    } catch (Exception e) {
                    }
                    if (paradises.size() >= targetCount)
                        break;
                }
            }
        }
        executor.shutdown();

        long sum = 0;
        for (long p : paradises)
            sum += p;
        return String.valueOf(sum);
    }

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