Project Euler Problem 196

Build a triangle from all positive integers in the following way: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

Project Euler Problem 196

Solution

Answer: 322303240771079935

Let the entries of row $n$ be

$$T(n,k)=\frac{n(n-1)}2+k,\qquad 1\le k\le n.$$

So row $n$ contains the consecutive integers

$$\frac{n(n-1)}2+1,\dots,\frac{n(n+1)}2.$$

We must compute

$$S(5678027)+S(7208785),$$

where $S(n)$ is the sum of the primes in row $n$ that belong to at least one prime triplet.


Mathematical analysis

1. Neighbour structure

For a number at position $(n,k)$, the neighbours are all valid positions among

$$(n+\delta_n,;k+\delta_k),$$

with

$$\delta_n,\delta_k\in{-1,0,1},$$

excluding $(0,0)$.

Thus an interior element has up to $8$ neighbours.


2. What does it mean to belong to a prime triplet?

A prime triplet is a set of three primes such that one of them is adjacent to the other two.

Define the prime degree of a prime $p$ to be the number of neighbouring primes.

Then:

  • if a prime has degree $\ge 2$, it is itself the centre of a prime triplet;
  • if a prime has degree $1$, it can still belong to a triplet if its unique prime neighbour has degree $\ge 2$.

Therefore:

A prime $p$ belongs to some prime triplet iff

  • $p$ has at least two prime neighbours, or
  • one of its neighbouring primes has at least two prime neighbours.

This reduces the problem to a purely local computation.


3. Locality

To determine whether a prime in row $n$ belongs to a triplet, we only need information from rows

$$n-2,\ n-1,\ n,\ n+1,\ n+2.$$

Indeed:

  • neighbours of row $n$ lie in rows $n-1,n,n+1$;
  • to know the neighbour-degrees of those primes, we need one more layer.

Thus the computation is feasible even for $n\approx 7\times10^6$.


4. Prime generation

The numbers involved are about

$$\frac{n^2}{2}\approx 2.6\times10^{13},$$

far too large for ordinary sieving over the whole range.

But only about $5n$ consecutive integers are needed, so a segmented sieve is ideal.


Python implementation

from math import isqrt

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

def value(n, k):
    return row_start(n) + k - 1

def neighbours(n, k):
    """
    Generate all neighbouring coordinates of (n,k).
    """
    out = []

    for dn in (-1, 0, 1):
        rn = n + dn
        if rn < 1:
            continue

        for dk in (-1, 0, 1):
            if dn == 0 and dk == 0:
                continue

            kk = k + dk

            if 1 <= kk <= rn:
                out.append((rn, kk))

    return out

def segmented_sieve(lo, hi):
    """
    Return a boolean array is_prime[x - lo]
    for all x in [lo, hi].
    """

    size = hi - lo + 1
    is_prime = [True] * size

    if lo <= 1 <= hi:
        is_prime[1 - lo] = False

    limit = isqrt(hi)

    # ordinary sieve up to sqrt(hi)
    base = [True] * (limit + 1)
    base[0] = base[1] = False

    for p in range(2, isqrt(limit) + 1):
        if base[p]:
            for m in range(p * p, limit + 1, p):
                base[m] = False

    primes = [i for i, v in enumerate(base) if v]

    # segmented marking
    for p in primes:
        start = max(p * p, ((lo + p - 1) // p) * p)

        for m in range(start, hi + 1, p):
            is_prime[m - lo] = False

    return is_prime

def S(n):
    """
    Sum of primes in row n that belong to some prime triplet.
    """

    # enough rows to determine all local structures
    lo = row_start(n - 2)
    hi = row_start(n + 3) - 1

    is_prime_segment = segmented_sieve(lo, hi)

    def is_prime(x):
        return is_prime_segment[x - lo]

    # prime neighbour counts
    degree = {}

    for r in range(n - 1, n + 2):
        for k in range(1, r + 1):

            x = value(r, k)

            if not is_prime(x):
                continue

            cnt = 0

            for rr, kk in neighbours(r, k):
                y = value(rr, kk)

                if is_prime(y):
                    cnt += 1

            degree[(r, k)] = cnt

    total = 0

    for k in range(1, n + 1):

        x = value(n, k)

        if not is_prime(x):
            continue

        ok = False

        # case 1: x itself is centre of a triplet
        if degree.get((n, k), 0) >= 2:
            ok = True

        # case 2: some neighbour is centre
        else:
            for rr, kk in neighbours(n, k):

                y = value(rr, kk)

                if is_prime(y) and degree.get((rr, kk), 0) >= 2:
                    ok = True
                    break

        if ok:
            total += x

    return total

answer = S(5678027) + S(7208785)

print(answer)

Code walkthrough

Coordinate system

row_start(n)

computes the first number in row $n$:

$$1+\frac{n(n-1)}2.$$

Then

value(n, k)

gives the actual integer at position $(n,k)$.


Neighbour generation

neighbours(n, k)

enumerates all adjacent cells using offsets $-1,0,+1$ in both row and column directions.

Invalid coordinates are discarded.


Segmented sieve

Instead of sieving all integers up to $2.6\times10^{13}$, we sieve only the short interval actually needed.

This makes the computation practical.


Degree computation

For every prime in rows $n-1,n,n+1$, we count how many neighbouring entries are prime.

That number is stored in:

degree[(r, k)]

Membership test

A prime belongs to a triplet if either:

degree[(n,k)] >= 2

or one of its neighbouring primes has degree at least $2$.


Small checks

The program reproduces the examples:

$$S(8)=60$$

and

$$S(9)=37.$$

It also reproduces the given verification value:

$$S(10000)=950007619.$$

So the implementation matches the problem statement exactly.


Final computation

The program gives

$$S(5678027)+S(7208785) = 322303240771079935.$$

The magnitude is reasonable:

  • each row contains millions of numbers near $10^{13}$,
  • about $1/\log n$ are prime,
  • a substantial fraction participate in triplets.

So the final value being around $3\times10^{17}$ is consistent.


Answer: 322303240771079935