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

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. Practical Numbers

A positive integer $m$ is practical if every integer from $1$ to $m$ can be written as a sum of distinct divisors of $m$.

The code only needs the practical-number test for five nearby values around a candidate $n$:

$$n-8,\ n-4,\ n,\ n+4,\ n+8.$$

This is why the candidate must satisfy $n\ge 9$: otherwise $n-8$ would be non-positive.

3. Stewart-Sierpiński Criterion

Let a number have prime factorization

$$N=\prod_{i=1}^{r} q_i^{a_i},\qquad 2=q_1<q_2<\cdots<q_r.$$

Define the prefix product

$$A_i=\prod_{j=1}^{i} q_j^{a_j},$$

and its divisor sum

$$S_i=\sigma(A_i).$$

The Stewart-Sierpiński theorem says that $N$ is practical if and only if $q_1=2$ and

$$q_{i+1}\le S_i+1\qquad(1\le i<r).$$

In words: once the first prime factor is 2, every next prime factor must not jump too far beyond the divisor-sum coverage already built by the prefix.

4. Why the Criterion Is Easy to Check Incrementally

The code does not recompute $\sigma(A_i)$ from scratch at every step. It factors the number in increasing prime order and maintains a running prefix value.

For a prime power,

$$\sigma(p^a)=1+p+p^2+\cdots+p^a=\frac{p^{a+1}-1}{p-1}.$$

So after factoring $N$, the solver only needs to multiply the current prefix by $\sigma(q_i^{a_i})$ and compare the next prime with $S_i+1$.

A simple worked example shows the idea:

$$12=2^2\cdot 3.$$

For the prefix $2^2$,

$$S_1=\sigma(2^2)=1+2+4=7,$$

and the next prime factor is $3$, which satisfies $3\le 7+1$. So 12 is practical.

By contrast,

$$14=2\cdot 7$$

fails because the prefix after $2$ is $S_1=3$, but $7>3+1$. So 14 is not practical.

5. Why a Brute-Force Cross-Check Exists

The implementation also includes a reference test for practical numbers. It enumerates all divisors, runs a subset-sum style reachability DP, and verifies that every value from 1 to $n$ is reachable as a sum of distinct divisors.

This slow method is used only for verification on small inputs. The program checks that the Stewart-Sierpiński criterion and the brute-force definition agree for every $n\le 400$.

This checkpoint is stronger than a casual spot-check. It compares the exact implementation of is_practical against the definition itself for every integer in the whole range $1\le n\le 400$, so mistakes in factor handling, divisor generation, or subset-sum reachability are caught before the large search starts.

6. Segmented Prime Generation

Prime generation is split into two layers.

First, the code computes all base primes up to $\sqrt{\text{limit}+9}$ with an ordinary sieve. Those primes are used both for segmented sieving and for factoring candidate numbers.

Then the search range $[3,\text{limit}+9]$ is processed in segments. Each segment stores only odd numbers, so the sieve memory stays small. Segment batches can be processed in parallel, but the primes inside each batch are still handled in increasing order so the sliding window stays correct.

7. The Paradise Search

The function find_paradises reads the prime stream and maintains a deque of the last four primes. Whenever the last four primes satisfy the gap-6 pattern, the code forms the candidate

$$n=p_0+9$$

and then tests the five practical-number conditions.

If they all pass, the candidate is recorded, added to the running sum, and printed as a found paradise.

The code stops once it has found the requested number of paradises. The first two known values are used as checkpoints:

$$219869980,\qquad 312501820.$$

8. Candidate Validation and Threading

Threading is used only for segmented sieve batches. Candidate validation itself is kept ordered and deterministic because the prime-window scan is sequential. The helper choose_thread_count caps the number of workers by the workload and the hardware.

This separation is important: the sieve can be parallelized safely, while the window scan must preserve prime order.

How the Code Works

The solver parses options such as --limit, --target, --segment, --threads, --single-thread, and --skip-checkpoints. It builds the base-prime list once, runs the practical-number cross-check up to 400 unless disabled, and then calls find_paradises.

The practical test is implemented in is_practical. The brute-force reference is is_practical_bruteforce. The segmented sieve lives in sieve_segment. The final search loop is a four-prime sliding window in find_paradises.

One subtle detail is that segment batches may be sieved in parallel, but the lambda process_prime still consumes the resulting prime lists in increasing segment order. That preserves the meaning of “four consecutive primes”, which would be broken by an unordered merge.

At the end, the program prints the sum of the first requested paradises and also reports elapsed time.

Complexity Analysis

The ordinary sieve up to $\sqrt{\text{limit}+9}$ is tiny compared with the full search. The segmented sieve is close to

$$O(L\log\log L)$$

for limit $L$, while using only segment-sized working memory. The practical-number test for a candidate factors five nearby numbers and applies the multiplicative criterion, so candidate checks are sparse and relatively cheap.

The brute-force practical-number check is much slower, but it runs only on the small checkpoint range $1\le n\le 400$. That makes it a safe validator without affecting the main asymptotics.

Footnotes and References

Problem page: https://projecteuler.net/problem=263
Practical numbers: Wikipedia - Practical number
Stewart-Sierpiński theorem: Wikipedia - Characterization of practical numbers
Segmented sieve: Wikipedia - Segmented sieve
Prime factorization and divisor sums: Wikipedia - Divisor function

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());
    }
}