Problem 196: Prime Triplets

View on Project Euler

Project Euler Problem 196 Solution

EulerSolve provides an optimized solution for Project Euler Problem 196, Prime Triplets, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The positive integers are written in triangular rows: row \(r\) contains \(r\) consecutive values, from \(T_{r-1}+1\) to \(T_r\), where \(T_r=\frac{r(r+1)}{2}\). For a fixed row \(n\), the quantity \(S(n)\) is the sum of those primes in row \(n\) that belong to at least one prime triplet under the problem's neighborhood rule. The full task asks for \(S(5678027)+S(7208785)\). The key fact extracted from the implementations is that triplet membership is a purely local property: even for million-scale rows, the answer can be determined from a five-row strip centered on the target row. Mathematical Approach It is convenient to index the triangle by row and column. For \(0 \le c < r\), let $$a(r,c)=\frac{r(r-1)}{2}+1+c,$$ so \(a(r,c)\) is the value in row \(r\), column \(c\). Triangular Coordinates and the Neighborhood For each cell \((r,c)\), the implementations examine the clipped \(3\times3\) window $$N(r,c)=\{(r+\Delta r,c+\Delta c): \Delta r,\Delta c\in\{-1,0,1\},\ 0\le c+\Delta c<r+\Delta r\}.$$ This set contains the center cell itself, the two horizontal neighbors in the same row when they exist, and up to three cells in the row above and three in the row below. Interior cells therefore have eight neighbors; boundary cells have fewer because invalid coordinates are discarded. Let \(P(r,c)\) be 1 when \(a(r,c)\) is prime and 0 otherwise....

Detailed mathematical approach

Problem Summary

All positive integers are written in rows: row \(1\) holds \(1\), row \(2\) holds \(2,3\), row \(3\) holds \(4,5,6\), and in general row \(r\) holds the next \(r\) consecutive integers. With the rows aligned on the left, as in the picture of the statement, every number has up to eight neighbours: the numbers next to it in its row and the three nearest numbers in the rows above and below. A set of three primes is a prime triplet if one of the three primes has the other two as neighbours. \(S(n)\) is the sum of the primes in row \(n\) that belong to at least one prime triplet. The statement notes that \(2\) and \(3\) belong to triplets, gives \(S(8)=60\), from \(29\) and \(31\), \(S(9)=37\) and \(S(10000)=950007619\), and asks for \(S(5678027)+S(7208785)\). [1] The name has nothing to do with the classical prime triplets, the sets of primes \((p,p+2,p+6)\) or \((p,p+4,p+6)\). [2]

The whole solution in one sentence: whether a prime of row \(n\) belongs to a triplet depends only on which numbers of the five rows \(n-2\) to \(n+2\) are prime, so we find the primes among these \(5n\) consecutive integers with a segmented sieve, mark every prime that is, or lies next to, a prime with at least two prime neighbours, and add up the marked primes of row \(n\).

Mathematical Approach

We describe the neighbours by rows and columns, find a local test for membership in a triplet, and show that five rows of primes decide \(S(n)\).

1. Rows, columns and neighbours

Row \(r\) ends with the triangular number \(T_r=1+2+\dots+r=r(r+1)/2\), [3] so it holds the numbers from \(T_{r-1}+1\) to \(T_r\). Numbering its columns from the left by \(c=0,1,\dots,r-1\), the number in row \(r\) and column \(c\) is

$$a(r,c)=\frac{r(r-1)}{2}+1+c.$$

For example \(a(8,0)=29\) and \(a(9,0)=37\): row \(8\) is \(29,\dots,36\) and row \(9\) is \(37,\dots,45\). The neighbours of the cell \((r,c)\) are the cells \((r,c\pm1)\), \((r-1,c+\delta)\) and \((r+1,c+\delta)\) with \(\delta\in\{-1,0,1\}\) that exist, that is, whose column lies between \(0\) and the length of their row minus \(1\). Since \(a(r+1,c)-a(r,c)=r\), the neighbours of \(x=a(r,c)\) are

  • in its own row: \(x-1\) and \(x+1\);
  • in the row above: \(x-r\), \(x-r+1\) and \(x-r+2\);
  • in the row below: \(x+r-1\), \(x+r\) and \(x+r+1\),

as far as these cells exist. The differences must not be used blindly at the ends of a row: for \(29=a(8,0)\), the number \(x-r=21\) is the last number of row \(6\), not a neighbour, and for \(45=a(9,8)\), the number \(x-r+1=37\) is the first number of row \(9\) itself. A program should therefore work with rows and columns, as the published code does.

2. Odd numbers have at most three odd neighbours

Apart from \(2\), every prime is odd. For an odd \(x\), the neighbours \(x\pm1\) in its own row are even. The three numbers above it are consecutive, and so are the three below, so in each of these rows either the middle number or the two outer ones are odd, where they exist:

  • if \(r\) is even, the odd neighbours are \(x-r\) and \(x-r+2\) above and \(x+r\) below;
  • if \(r\) is odd, they are \(x-r+1\) above and \(x+r-1\) and \(x+r+1\) below.

So every odd number has at most three odd neighbours. In row \(7\) the prime \(23\) has the odd neighbours \(17\), \(29\) and \(31\), all three prime; in row \(9\) the prime \(41\) has the odd neighbours \(33\), \(49\) and \(51\), none of them prime. The prime \(2\) has only the neighbours \(1\), \(3\), \(4\) and \(5\), so from row \(4\) on a prime has at most three prime neighbours.

3. Which primes belong to a triplet

Call a prime \(q\) a centre if at least two of its neighbours are prime, and call the cell of \(q\) together with the cells of its neighbours the window of \(q\), a block of at most \(3\times3\) cells.

Claim. A prime \(p\) belongs to a prime triplet if and only if \(p\) lies in the window of some centre. Since a window contains its own centre, \(q\) is a centre exactly when its window contains at least three primes.

Proof. If \(\{x,q,y\}\) is a prime triplet in which \(q\) has \(x\) and \(y\) as neighbours, then \(q\) is a centre and all three primes lie in its window. Conversely, let \(q\) be a centre and \(p\) a prime in its window. If \(p=q\), then \(q\) and two of its prime neighbours form a triplet. Otherwise \(p\) is a neighbour of \(q\); since \(q\) has at least two prime neighbours, it has a prime neighbour \(p'\ne p\), and \(\{p,q,p'\}\) is a triplet with \(q\) in the middle. \(\square\)

For example, \(23\) is a centre, so the primes \(17\), \(23\), \(29\) and \(31\) of its window all belong to triplets, such as \(\{17,23,29\}\) and \(\{23,29,31\}\).

4. Clusters of primes

Link two primes when they are neighbours. The primes then fall into clusters, the connected components of this graph: two primes share a cluster when a chain of links joins them. [4]

Claim. A prime belongs to a triplet if and only if its cluster contains at least three primes.

Proof. A triplet is held together by the links of its centre, so it lies inside one cluster. Conversely, if the cluster of \(p\) has at least three primes, \(p\) has a prime neighbour \(q\). If \(p\) has a second one, \(p\) is a centre. Otherwise \(q\) has a prime neighbour other than \(p\), since otherwise \(p\) and \(q\) would form a cluster of two; then \(q\) is a centre and \(p\) lies in its window. \(\square\)

So a prime neighbour alone is not enough: \(43\) and \(53\), in rows \(9\) and \(10\), are each other's only prime neighbour, a cluster of two, and neither belongs to a triplet.

5. Five rows are enough

Let \(p\) be a prime of row \(n\) that belongs to a triplet. By Section 3 it lies in the window of a centre \(q\), so \(q\) is \(p\) or a neighbour of \(p\) and lies in row \(n-1\), \(n\) or \(n+1\), and the window of \(q\) lies in rows \(n-2\) to \(n+2\). Hence the primes of these five rows determine \(S(n)\), and only the windows of the primes of rows \(n-1\) to \(n+1\) need to be examined. Consecutive rows hold consecutive integers, so the five rows form the single interval

$$T_{n-3}+1\le x\le T_{n+2},$$

with \(T_{n+2}-T_{n-3}=5n\) numbers, as the five rows hold \(n-2\), \(n-1\), \(n\), \(n+1\) and \(n+2\) of them. For \(n=7\,208\,785\) it runs from \(25\,983\,272\,566\,154\) to \(25\,983\,308\,610\,078\), \(36\,043\,925\) numbers near \(2.6\cdot10^{13}\); for \(n=5\,678\,027\) it holds \(28\,390\,135\) numbers up to \(16\,120\,009\,501\,435\).

6. Worked example: \(S(8)\) and \(S(9)\)

For \(n=8\) we need rows \(6\) to \(10\), the numbers \(16\) to \(55\), whose primes are \(17\), \(19\), \(23\), \(29\), \(31\), \(37\), \(41\), \(43\), \(47\) and \(53\). The candidates for centres are the primes of rows \(7\) to \(9\):

  • \(23\) has the prime neighbours \(17\), \(29\) and \(31\): a centre;
  • \(29\) has \(23\) and \(37\): a centre;
  • \(37\) has \(29\) and \(47\): a centre;
  • \(31\) has only \(23\), \(43\) only \(53\), and \(41\) none: no centres.

The windows of the three centres contain the primes \(17,23,29,31\), then \(23,29,37\), then \(29,37,47\). Both primes of row \(8\) occur, so \(S(8)=29+31=60\); \(29\) lies in all three windows but is added once, which is why the programs mark members with a flag. For \(n=9\) we need rows \(7\) to \(11\), the numbers \(22\) to \(66\), and the windows of the primes of rows \(8\) to \(10\). Row \(10\) adds \(47\), whose only prime neighbour is \(37\), and \(53\), whose only one is \(43\), so the centres are \(29\) and \(37\). Of the three primes \(37\), \(41\) and \(43\) of row \(9\), only \(37\) lies in their windows, and \(S(9)=37\).

7. The primes of a far-away interval

To find the primes of an interval \(L\le x\le U\), we take every prime \(p\le\sqrt U\), cross out its multiples from \(\max(p^2,\lceil L/p\rceil p)\) up to \(U\), and also cross out \(1\) if \(L=1\). This is the segmented form of the sieve of Eratosthenes. [5] [6]

Claim. The numbers that are not crossed out are exactly the primes of the interval.

Proof. Apart from \(1\), which is not prime, a crossed-out number is \(jp\) with \(j\ge p\ge2\), hence composite. A composite \(m\) of the interval has a smallest prime factor \(p\) with \(p^2\le m\le U\), so \(p\le\sqrt U\), and \(m\) is a multiple of \(p\) with \(m\ge p^2\) and \(m\ge L\), so it is crossed out. \(\square\)

Starting at \(p^2\) matters only for small rows, whose interval may contain sieving primes: for \(n=2\), with an empty row \(0\), it is \(1\) to \(10\), and \(2\) and \(3\) must survive. For the worked example, the primes \(2\), \(3\), \(5\) and \(7\) up to \(\sqrt{55}\approx7.4\) leave exactly the ten primes listed in Section 6.

8. The cost of the sieve

The prime \(p\) crosses out at most \((U-L)/p+1\) numbers, so the sieve makes at most \(5n\sum_{p\le\sqrt U}1/p+\pi(\sqrt U)\) steps, where \(\pi(y)\) counts the primes up to \(y\). By Mertens' second theorem, [7]

$$\sum_{p\le y}\frac1p=\ln\ln y+M+O\Bigl(\frac{1}{\ln y}\Bigr),$$

with the Meissel–Mertens constant \(M=0.2614972\ldots\), and \(\sqrt U\approx n/\sqrt2\) because \(U=T_{n+2}\approx n^2/2\), so the sieve takes \(O(n\log\log n)\) steps. For \(n=7\,208\,785\) the \(354\,802\) primes up to \(\lfloor\sqrt U\rfloor=5\,097\,382\) make \(108\,087\,448\) crossings, \(2.9988\) per number of the interval, against \(\ln\ln5\,097\,382+M\approx2.9987\). Few cells need a window: row \(10000\) has \(565\) primes, near the \(10000/\ln(5\cdot10^{7})\approx564\) suggested by the prime number theorem, [8] a heuristic for so short an interval, and only \(19\) of them belong to triplets.

How the Code Works

The C++ program [9] keeps its settings in Options. parse_arguments reads --row-a=, --row-b=, --row= (\(0\) for no single row) and --skip-checkpoints, with parse_int_after_prefix turning digits into an int without an overflow check; other arguments, or a row below \(2\) among those to be computed, give exit code \(1\). Unless switched off, run_checkpoints requires \(60\), \(37\) and \(950007619\) for the rows \(8\), \(9\) and \(10000\), else the exit code is \(2\), and main prints the value of the single row or the sum for the two rows.

solve_row(n) implements Sections 3, 5 and 7. row_start(r) computes \(r(r-1)/2+1\) in unsigned 64-bit arithmetic, and the interval runs from row_start(n - 2) to row_start(n + 3) - 1, that is from \(T_{n-3}+1\) to \(T_{n+2}\). A std::vector<std::uint8_t> named is_prime holds one flag per number, with the number \(1\) cleared if the interval starts there. Every prime from small_primes_up_to(max_p), where max_p is \(\lfloor\sqrt U\rfloor+1\), clears its multiples from the larger of \(p^2\) and its first multiple \(\ge L\), and the loop stops at the first prime with \(p^2\gt U\). The flags are copied into five row-shaped arrays row_prime for rows \(n-2\) to \(n+2\), beside five zeroed arrays row_good. For every prime cell of rows \(n-1\) to \(n+1\), the program counts the primes of its window in prime_neighbors_including_self, skipping columns outside the triangle by the test against rr_len, as Section 1 demands; if the count is at least \(3\), it sets row_good for every prime of the window. Finally the values row_start(n) + col of the marked cells of row \(n\) are added in an unsigned 64-bit sum.

small_primes_up_to is a linear sieve, also known as Euler's sieve. [5] [10] For each \(i\) from \(2\) to the bound, an \(i\) that is still unmarked is appended to the list of primes, and then \(i\cdot p\) is marked for the listed primes \(p\) in increasing order, stopping after the first \(p\) that divides \(i\) or once \(i\cdot p\) exceeds the bound. The loop for \(i\) runs through the primes up to the smallest prime factor of \(i\) while \(i\cdot p\) stays within the bound, so every marked \(i\cdot p\) has the smallest prime factor \(p\), and a composite \(m\) with smallest prime factor \(p\) is marked from \(i=m/p\), whose smallest prime factor is at least \(p\), and from no other \(i\). The work is proportional to the bound, about \(5.1\) million for the larger row.

The Python program follows the C++ function by function, with bytearray flags, math.isqrt and unbounded integers, but without options or checkpoints. The Java program has no checkpoints either; its sieveRange runs the classical sieve of Eratosthenes up to \(\lfloor\sqrt U\rfloor+1\) in an array sc, and each prime \(i\) it finds at once crosses out its multiples in the interval, from \(\max(i^2,\lceil L/i\rceil i)\), in an array comp that is then copied, negated, into the prime flags. Java keeps values and sums in long and array indices in int: the interval has fewer than \(2^{31}\) numbers, but its values, about \(2^{44.6}\), need 64 bits. The two rows contain \(187\,270\) and \(233\,126\) primes, each at most \(T_n\), so all of them together sum to less than

$$\begin{aligned}&187\,270\,T_{5678027}+233\,126\,T_{7208785}\\&\qquad\approx9.08\cdot10^{18},\end{aligned}$$

which is below \(2^{63}-1\approx9.22\cdot10^{18}\), so neither 64-bit sum can overflow.

Complexity and Verification

For one row, the segmented sieve takes \(O(n\log\log n)\) time by Section 8, the base primes \(O(\sqrt U)=O(n)\) with the linear sieve, or \(O(n\log\log n)\) with Java's classical sieve, and the copy into rows \(O(n)\). The window loop visits the \(3n\) cells of rows \(n-1\) to \(n+1\) and reads at most nine flags for each prime cell and nine more for a centre, again \(O(n)\). The memory is one byte per number of the interval in each of is_prime, row_prime and row_good, that is \(15n\) bytes of flags, about \(108\) MB for \(n=7\,208\,785\), plus the base primes; Java has ip, rp and rg and, inside sieveRange, a temporary comp. These arrays serve one row only. Measured with the system's time command on an Apple M4 Max laptop, with Apple clang -O2, OpenJDK 17 and CPython 3.11, a full run takes about \(0.2\) s in C++, checkpoints included, \(0.3\) s in Java and \(7.4\) s in Python.

For this article, a brute force that tests the definition literally, trying every set of three primes within two steps of each prime of the row, agreed with the cluster description of Section 4 for every row \(2\le n\le300\). The cluster description, computed with union–find over rows \(1\) to \(1503\), gave the same \(S(n)\) as the C++ program, run with --row=n, as the Java method solveRow and as the Python solve_row, for every \(2\le n\le1500\). An independent program that uses only the three odd neighbours of Section 2, the centre test of Section 3 and its own segmented sieve agreed for every \(5\le n\le1500\), reproduced \(S(10000)=950007619\) and gave the same final result as the three published programs. These checks are not part of the published solutions.

Footnotes and References

The references below supply the background facts; the criteria for triplets, the five-row interval and the operation counts are derived in this article.

  1. Project Euler 196 — Prime Triplets. The official statement builds the triangle of all positive integers, calls a set of three primes a prime triplet if one of them has the other two as neighbours, gives S(8) = 60, S(9) = 37 and S(10000) = 950007619, and asks for S(5678027) + S(7208785).
  2. Wikipedia — Prime triplet. The classical prime triplets are sets of three primes of the form (p, p + 2, p + 6) or (p, p + 4, p + 6); this problem uses the name differently.
  3. Wikipedia — Triangular number. Tₙ = 1 + 2 + ⋯ + n = n(n + 1)/2 counts the points of a triangle with rows of 1, 2, …, n points; here it is the last number of row n.
  4. Wikipedia — Component (graph theory). A component of an undirected graph is a maximal connected subgraph, every vertex lies in exactly one, and a disjoint-set (union–find) structure can maintain the components as edges are added.
  5. Wikipedia — Sieve of Eratosthenes. Marking the multiples of p may start at p², and the primes up to the square root of the limit suffice; the segmented sieve treats a range piece by piece with the primes up to the square root of its top, in O(n log log n) operations. Euler's sieve removes each composite once, in linear time.
  6. cp-algorithms — Sieve of Eratosthenes. For a short range [L, R] with large R, the primes up to √R are generated first, and each marks its multiples from the larger of p² and its first multiple in the range.
  7. Wikipedia — Mertens' theorems. The second theorem: the sum of 1/p over the primes p ≤ n is log log n + M + O(1/log n), where M ≈ 0.2614972128 is the Meissel–Mertens constant.
  8. Wikipedia — Prime number theorem. A random integer not greater than N is prime with probability very close to 1/log N.
  9. C++ — Euler196.cpp. The linked immutable revision contains parse_arguments, row_start, small_primes_up_to, solve_row and run_checkpoints. It is the implementation record for this article.
  10. cp-algorithms — Linear Sieve. Each i marks i·p for the primes p up to the least prime factor of i, so every composite is reached exactly once, through its least prime factor, in O(n) total time.

Problem 196 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

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

struct Options {
    int row_a = 5678027;
    int row_b = 7208785;
    int single_row = 0;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = 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 == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--row-a=", options.row_a) ||
            parse_int_after_prefix(arg, "--row-b=", options.row_b) ||
            parse_int_after_prefix(arg, "--row=", options.single_row)) {
            continue;
        }

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

    if (options.single_row != 0) {
        return options.single_row >= 2;
    }
    return options.row_a >= 2 && options.row_b >= 2;
}

u64 row_start(const int r) {
    return static_cast<u64>(r) * static_cast<u64>(r - 1) / 2ULL + 1ULL;
}

std::vector<int> small_primes_up_to(const int n) {
    std::vector<std::uint8_t> is_composite(static_cast<std::size_t>(n + 1), 0);
    std::vector<int> primes;

    for (int i = 2; i <= n; ++i) {
        if (!is_composite[static_cast<std::size_t>(i)]) {
            primes.push_back(i);
        }
        for (int p : primes) {
            const long long v = 1LL * i * p;
            if (v > n) {
                break;
            }
            is_composite[static_cast<std::size_t>(v)] = 1;
            if (i % p == 0) {
                break;
            }
        }
    }

    return primes;
}

u64 solve_row(const int n) {
    const int r_min = n - 2;
    const int r_max = n + 2;

    const u64 start_value = row_start(r_min);
    const u64 end_value = row_start(r_max + 1) - 1;

    const u64 range_len = end_value - start_value + 1;
    std::vector<std::uint8_t> is_prime(static_cast<std::size_t>(range_len), 1);

    if (start_value == 1ULL) {
        is_prime[0] = 0;
    }

    const int max_p = static_cast<int>(std::sqrt(static_cast<long double>(end_value))) + 1;
    const std::vector<int> primes = small_primes_up_to(max_p);

    for (int p : primes) {
        const u64 pp = static_cast<u64>(p) * static_cast<u64>(p);
        if (pp > end_value) {
            break;
        }

        u64 first = (start_value + static_cast<u64>(p) - 1ULL) / static_cast<u64>(p);
        first *= static_cast<u64>(p);
        if (first < pp) {
            first = pp;
        }

        for (u64 x = first; x <= end_value; x += static_cast<u64>(p)) {
            is_prime[static_cast<std::size_t>(x - start_value)] = 0;
        }
    }

    std::vector<std::vector<std::uint8_t>> row_prime(5);
    std::vector<std::vector<std::uint8_t>> row_good(5);

    for (int idx = 0; idx < 5; ++idx) {
        const int row = r_min + idx;
        row_prime[static_cast<std::size_t>(idx)].assign(static_cast<std::size_t>(row), 0);
        row_good[static_cast<std::size_t>(idx)].assign(static_cast<std::size_t>(row), 0);

        const u64 rs = row_start(row);
        for (int col = 0; col < row; ++col) {
            const u64 value = rs + static_cast<u64>(col);
            row_prime[static_cast<std::size_t>(idx)][static_cast<std::size_t>(col)] =
                is_prime[static_cast<std::size_t>(value - start_value)];
        }
    }

    for (int rid = 1; rid <= 3; ++rid) {
        const int row_len = r_min + rid;
        for (int col = 0; col < row_len; ++col) {
            if (!row_prime[static_cast<std::size_t>(rid)][static_cast<std::size_t>(col)]) {
                continue;
            }

            int prime_neighbors_including_self = 0;
            for (int dr = -1; dr <= 1; ++dr) {
                const int rr = rid + dr;
                const int rr_len = r_min + rr;
                for (int dc = -1; dc <= 1; ++dc) {
                    const int cc = col + dc;
                    if (cc < 0 || cc >= rr_len) {
                        continue;
                    }
                    prime_neighbors_including_self += row_prime[static_cast<std::size_t>(rr)][static_cast<std::size_t>(cc)] != 0;
                }
            }

            if (prime_neighbors_including_self >= 3) {
                for (int dr = -1; dr <= 1; ++dr) {
                    const int rr = rid + dr;
                    const int rr_len = r_min + rr;
                    for (int dc = -1; dc <= 1; ++dc) {
                        const int cc = col + dc;
                        if (cc < 0 || cc >= rr_len) {
                            continue;
                        }
                        if (row_prime[static_cast<std::size_t>(rr)][static_cast<std::size_t>(cc)]) {
                            row_good[static_cast<std::size_t>(rr)][static_cast<std::size_t>(cc)] = 1;
                        }
                    }
                }
            }
        }
    }

    u64 sum = 0;
    const int mid_idx = 2;
    const u64 mid_start = row_start(n);
    for (int col = 0; col < n; ++col) {
        if (row_good[static_cast<std::size_t>(mid_idx)][static_cast<std::size_t>(col)]) {
            sum += mid_start + static_cast<u64>(col);
        }
    }
    return sum;
}

bool run_checkpoints() {
    if (solve_row(8) != 60ULL) {
        std::cerr << "Checkpoint failed for S(8)" << '\n';
        return false;
    }
    if (solve_row(9) != 37ULL) {
        std::cerr << "Checkpoint failed for S(9)" << '\n';
        return false;
    }
    if (solve_row(10000) != 950007619ULL) {
        std::cerr << "Checkpoint failed for S(10000)" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    if (options.single_row != 0) {
        std::cout << solve_row(options.single_row) << '\n';
    } else {
        std::cout << (solve_row(options.row_a) + solve_row(options.row_b)) << '\n';
    }
    return 0;
}

Python

import math

def solve():
    row_a = 5678027
    row_b = 7208785

    def row_start(r):
        return r * (r - 1) // 2 + 1

    def small_primes_up_to(n):
        is_composite = bytearray(n + 1)
        primes = []
        for i in range(2, n + 1):
            if not is_composite[i]:
                primes.append(i)
            for p in primes:
                v = i * p
                if v > n:
                    break
                is_composite[v] = 1
                if i % p == 0:
                    break
        return primes

    def solve_row(n):
        r_min = n - 2
        r_max = n + 2
        start_value = row_start(r_min)
        end_value = row_start(r_max + 1) - 1
        range_len = end_value - start_value + 1
        is_prime = bytearray(b'\x01' * range_len)
        if start_value == 1:
            is_prime[0] = 0

        max_p = math.isqrt(end_value) + 1
        primes = small_primes_up_to(max_p)

        for p in primes:
            pp = p * p
            if pp > end_value:
                break
            first = ((start_value + p - 1) // p) * p
            if first < pp:
                first = pp
            for x in range(first - start_value, range_len, p):
                is_prime[x] = 0

        row_prime = []
        row_good = []
        for idx in range(5):
            row = r_min + idx
            rp = bytearray(row)
            rg = bytearray(row)
            rs = row_start(row)
            for col in range(row):
                value = rs + col
                rp[col] = is_prime[value - start_value]
            row_prime.append(rp)
            row_good.append(rg)

        for rid in range(1, 4):
            row_len = r_min + rid
            for col in range(row_len):
                if not row_prime[rid][col]:
                    continue
                cnt = 0
                for dr in range(-1, 2):
                    rr = rid + dr
                    rr_len = r_min + rr
                    for dc in range(-1, 2):
                        cc = col + dc
                        if 0 <= cc < rr_len:
                            cnt += (row_prime[rr][cc] != 0)
                if cnt >= 3:
                    for dr in range(-1, 2):
                        rr = rid + dr
                        rr_len = r_min + rr
                        for dc in range(-1, 2):
                            cc = col + dc
                            if 0 <= cc < rr_len and row_prime[rr][cc]:
                                row_good[rr][cc] = 1

        mid_idx = 2
        mid_start = row_start(n)
        s = 0
        for col in range(n):
            if row_good[mid_idx][col]:
                s += mid_start + col
        return s

    return str(solve_row(row_a) + solve_row(row_b))

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

Java

public class Euler196 {
    static long rowStart(int r) {
        return (long) r * (r - 1) / 2 + 1;
    }

    static boolean[] sieveRange(long start, long end) {
        int rng = (int) (end - start + 1);
        boolean[] comp = new boolean[rng];
        if (start <= 1)
            comp[(int) (1 - start)] = true;
        int maxP = (int) Math.sqrt((double) end) + 1;
        boolean[] sc = new boolean[maxP + 1];
        for (int i = 2; i <= maxP; i++) {
            if (sc[i])
                continue;
            for (int j = i * 2; j <= maxP; j += i)
                sc[j] = true;
            long first = ((start + i - 1) / i) * i;
            if (first < (long) i * i)
                first = (long) i * i;
            for (long x = first; x <= end; x += i)
                comp[(int) (x - start)] = true;
        }
        boolean[] ip = new boolean[rng];
        for (int i = 0; i < rng; i++)
            ip[i] = !comp[i];
        if (start == 0)
            ip[0] = false;
        return ip;
    }

    static long solveRow(int n) {
        int rMin = n - 2, rMax = n + 2;
        long sv = rowStart(rMin), ev = rowStart(rMax + 1) - 1;
        boolean[] ip = sieveRange(sv, ev);
        boolean[][] rp = new boolean[5][], rg = new boolean[5][];
        for (int idx = 0; idx < 5; idx++) {
            int row = rMin + idx;
            rp[idx] = new boolean[row];
            rg[idx] = new boolean[row];
            long rs = rowStart(row);
            for (int col = 0; col < row; col++)
                rp[idx][col] = ip[(int) (rs + col - sv)];
        }
        for (int rid = 1; rid <= 3; rid++) {
            int rl = rMin + rid;
            for (int col = 0; col < rl; col++) {
                if (!rp[rid][col])
                    continue;
                int cnt = 0;
                for (int dr = -1; dr <= 1; dr++) {
                    int rr = rid + dr, rrl = rMin + rr;
                    for (int dc = -1; dc <= 1; dc++) {
                        int cc = col + dc;
                        if (cc >= 0 && cc < rrl && rp[rr][cc])
                            cnt++;
                    }
                }
                if (cnt >= 3) {
                    for (int dr = -1; dr <= 1; dr++) {
                        int rr = rid + dr, rrl = rMin + rr;
                        for (int dc = -1; dc <= 1; dc++) {
                            int cc = col + dc;
                            if (cc >= 0 && cc < rrl && rp[rr][cc])
                                rg[rr][cc] = true;
                        }
                    }
                }
            }
        }
        long sum = 0;
        long ms = rowStart(n);
        for (int col = 0; col < n; col++)
            if (rg[2][col])
                sum += ms + col;
        return sum;
    }

    public static void main(String[] args) {
        System.out.println(solveRow(5678027) + solveRow(7208785));
    }
}