Problem 196: Prime Triplets

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

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. Then the local prime count around $(r,c)$ is

$$M(r,c)=\sum_{(u,v)\in N(r,c)} P(u,v).$$

An Equivalent Criterion for Prime Triplets

The implementations rely on the following equivalent characterization. A prime at $(r,c)$ can act as the center of a prime triplet exactly when $P(r,c)=1$ and $M(r,c)\ge 3$. The reason is that the count $M(r,c)$ includes the center itself, so a value of at least 3 means there are at least two additional primes in the same neighborhood.

This immediately gives the marking rule used in code: if a prime center has $M(r,c)\ge 3$, then every prime in $N(r,c)$ belongs to at least one valid triplet. The center is already adjacent to all of them, and there is always at least one more prime in the same neighborhood to complete a three-prime configuration. So the problem is not just to find centers; it is to mark every prime lying inside a qualifying neighborhood.

Why Only Five Rows Matter

Take any prime in row $n$. If it belongs to a triplet, then either it is the center or it is adjacent to the center. In either case that center must lie in row $n-1$, $n$, or $n+1$. Once the center is fixed, every cell in its neighborhood lies within one extra row, so all relevant values are contained in rows $n-2$ through $n+2$.

This is the central invariant of the solution. Instead of processing all values up to roughly $\frac{n^2}{2}$, it is enough to inspect the exact interval covered by those five rows:

$$\left[\frac{(n-2)(n-3)}{2}+1,\ \frac{(n+2)(n+3)}{2}\right].$$

The total number of cells in that strip is

$$ (n-2)+(n-1)+n+(n+1)+(n+2)=5n,$$

so the relevant search space is linear in $n$.

Segmented Sieving on the Exact Interval

Because the five-row strip is a contiguous block of integers, primality is determined with a segmented sieve on that interval. Only primes up to

$$\sqrt{\frac{(n+2)(n+3)}{2}}$$

are needed to cross out all composites. After sieving, the flat boolean array is projected back into five jagged rows, one for each triangular row from $n-2$ to $n+2$.

With that notation, the answer can be expressed as

$$S(n)=\sum_{c=0}^{n-1} a(n,c)\,G(n,c),$$

where $G(n,c)=1$ exactly when $a(n,c)$ is prime and lies inside $N(r,d)$ for some prime center $(r,d)$ with $r\in\{n-1,n,n+1\}$ and $M(r,d)\ge 3$.

Worked Example: $S(8)=60$ and $S(9)=37$

The small checkpoint rows show the method clearly. Row 8 is $29,30,31,32,33,34,35,36$, while row 9 is $37,38,39,40,41,42,43,44,45$. Consider the prime 23 in row 7. Its clipped neighborhood contains the numbers $16,17,18,22,23,24,29,30,31$, and the primes among them are $17,23,29,31$. Since that neighborhood already contains at least three primes, both 29 and 31 are certified as members of prime triplets, so $S(8)=29+31=60$.

For row 9, the prime 37 lies at the left edge, so its neighborhood is smaller: $29,30,37,38,46,47$. The prime subset is $29,37,47$, which is still enough to form a triplet. By contrast, the primes 41 and 43 never appear inside any qualifying neighborhood of size at least three primes, so they do not contribute. Hence $S(9)=37$. The large target rows are solved by exactly the same local logic.

How the Code Works

Sieving the Five-Row Strip

The C++, Python, and Java implementations first compute the start and end of the interval covering rows $n-2$ through $n+2$. They generate all base primes up to the square root of the interval endpoint, then mark composites inside the interval itself. The result is a primality table for every value that can possibly affect $S(n)$.

Rebuilding the Triangle Shape

Next, the flat segmented-sieve array is mapped back into five row-shaped boolean arrays. One array records which cells are prime, and a second array records which prime cells have already been proved to belong to at least one triplet. This layout mirrors the triangular geometry, so neighborhood scans become simple bounded loops instead of repeated coordinate arithmetic on absolute values.

Scanning Candidate Centers and Marking Members

The implementations inspect only the middle three rows, because only those rows can contain a triplet center relevant to row $n$. For each prime cell in those rows, the code counts primes in the clipped $3\times3$ neighborhood. If the count is at least three, every prime in that neighborhood is marked. A final pass over row $n$ sums exactly the marked primes.

The three languages differ only in low-level representation and in how the small supporting primes are generated. The mathematical workflow is the same in every version.

Complexity Analysis

The active interval has length exactly $5n$, and every later step works only inside that strip. The neighborhood phase is linear as well: there are $3n$ candidate-center positions in rows $n-1$, $n$, and $n+1$, and each one inspects at most nine cells.

Using a standard segmented sieve, the running time is $O(n \log \log n)$ in the usual sieve model, followed by an $O(n)$ marking and summation pass. Extra space is $O(n)$ for the interval primality array and the two five-row masks. The decisive optimization is locality: the algorithm never touches the whole triangle up to row $n$, only the narrow strip that can influence the answer.

Footnotes and References

Project Euler problem page: https://projecteuler.net/problem=196
Prime number: Wikipedia - Prime number
Triangular number: Wikipedia - Triangular number
Sieve of Eratosthenes: Wikipedia - Sieve of Eratosthenes
Segmented sieve overview: cp-algorithms - Sieve of Eratosthenes and segmented sieving

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