Problem 278: Linear Combinations of Semiprimes

View on Project Euler

Project Euler Problem 278 Solution

EulerSolve provides an optimized solution for Project Euler Problem 278, Linear Combinations of Semiprimes, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For distinct primes \(p \lt q \lt r\), define the three semiprimes $$a=pq,\qquad b=pr,\qquad c=qr.$$ For each prime triple, we want the largest integer that cannot be written as a nonnegative linear combination of \(a,b,c\). Then we sum that Frobenius value over all triples below the prime limit. As usual, the final Project Euler number is omitted; the purpose of this page is to explain the mathematics behind the C++ solution. Mathematical Approach 1) The semigroup we need to understand Let $$S=\langle pq,pr,qr\rangle=\{x\cdot pq+y\cdot pr+z\cdot qr:\ x,y,z\ge0\}.$$ The problem asks for the Frobenius number of this semigroup, namely the largest integer not belonging to \(S\). For three arbitrary generators, no simple closed formula usually exists. The reason this problem is solvable is that the generators share the prime factors \(p,q,r\) in a very rigid way. 2) Rewrite every representable number by residue modulo \(p\) Because $$pq=p\cdot q,\qquad pr=p\cdot r,$$ every representable number has the form $$N=z\cdot qr+p\cdot m,$$ where $$m\in \langle q,r\rangle=\{u\cdot q+v\cdot r:\ u,v\ge0\}.$$ So the whole three-generator problem reduces to this question: For each residue class modulo \(p\), how large can \(N\) be before the quotient part \(m\) is forced into the two-generator semigroup \(\langle q,r\rangle\)?...

Detailed mathematical approach

Problem Summary

Take whole numbers \(1\lt a_1\lt a_2\lt\dots\lt a_n\) and call a number \(b\) representable if \(b=q_1a_1+q_2a_2+\dots+q_na_n\) for some whole numbers \(q_k\ge0\); the number \(0\) is representable, with all \(q_k=0\). The largest number that is not representable, if there is one, is called \(f(a_1,\dots,a_n)\); for the numbers of this problem it exists, as Section 5 shows. For \(5\) and \(7\) the numbers that are not representable are \(1,2,3,4,6,8,9,11,13,16,18\) and \(23\), so \(f(5,7)=23\); the statement also gives \(f(6,10,15)=29\) and \(f(14,22,77)=195\). We must find the sum of \(f(pq,pr,qr)\) over all primes \(p\lt q\lt r\lt5000\); the three numbers are then in increasing order, \(pq\lt pr\lt qr\). [1] With coins of the values \(a_1,\dots,a_n\), \(f\) is the largest amount that cannot be paid exactly, the Frobenius number of the coin problem. [2]

The whole solution in one sentence: when the numbers are sorted by their remainder after division by the largest coin \(qr\), the smallest representable number with each remainder can be written down directly, which gives the formula \(f(pq,pr,qr)=2pqr-pq-pr-qr\), and the program adds this formula up over all triples of primes below \(5000\).

Mathematical Approach

We start with two coins, where everything can be checked by hand, and then apply the same tool, sorting by remainder, to the three semiprimes.

1. Coins of 5 and 7

With coins of \(5\) and \(7\) the representable numbers up to \(30\) are \(0,5,7,10,12,14,15,17,19,20,21,22\) and every number from \(24\) to \(30\). The missing positive numbers, the gaps, are exactly the twelve numbers of the statement. Once the five consecutive numbers \(24,25,26,27,28\) are representable, every larger number is representable too: \(29=24+5\), \(30=25+5\), and so on. In general, as soon as \(a_1\) consecutive numbers are representable, so are all larger numbers, because adding the coin \(a_1\) to each of them gives the next \(a_1\) numbers.

2. Sorting by remainder

Sort the numbers by their remainder after division by \(5\). The numbers with remainder \(3\) are \(3,8,13,18,23,28,33,\dots\); the first five are gaps, \(28=4\cdot7\) is the first representable one, and adding coins of \(5\) reaches all later ones. For the remainders \(0,1,2,3,4\) the first representable numbers are \(0,21,7,28,14\), that is \(0\), \(3\cdot7\), \(1\cdot7\), \(4\cdot7\) and \(2\cdot7\), and the largest gap is \(28-5=23\).

The same works for any coins. Let \(S\) be the set of representable numbers and let \(m\) be one of the coins. Suppose that every remainder \(i=0,1,\dots,m-1\) occurs in \(S\), and let \(w_i\) be the smallest number in \(S\) with remainder \(i\) after division by \(m\); if some remainder never occurs, all numbers with that remainder are gaps and there is no largest gap. The numbers of \(S\) with remainder \(i\) are then exactly \(w_i,w_i+m,w_i+2m,\dots\): each of these is \(w_i\) plus some coins \(m\), and every number of \(S\) with remainder \(i\) is at least \(w_i\) and differs from it by a multiple of \(m\). So the gaps with remainder \(i\) are the positive numbers among \(w_i-m,w_i-2m,\dots\), and

$$f=\max_{0\le i\lt m}w_i-m.$$

The numbers \(w_i\) form the Apéry set of \(m\), and this is Selmer's formula for the largest gap. [3] [5]

3. Two coins: Sylvester's formula

Let the two coins \(a\) and \(b\) be coprime, \(\gcd(a,b)=1\), and take \(m=a\). The numbers \(0,b,2b,\dots,(a-1)b\) have different remainders after division by \(a\): if \(ib\) and \(jb\) had the same remainder, \(a\) would divide \((i-j)b\), hence \(i-j\) by Euclid's lemma, [4] and \(|i-j|\lt a\) leaves only \(i=j\). So these \(a\) numbers meet every remainder exactly once. Each \(kb\) is also the smallest representable number with its remainder, because \(a\) coins \(b\) are worth as much as \(b\) coins \(a\): any \(xa+yb\) with \(x,y\ge0\) can be written with \(y=k+ta\), \(0\le k\lt a\), as \(kb+(x+tb)\,a\), so it has the same remainder as \(kb\), the only one of the \(a\) numbers with this remainder, and it is at least \(kb\). Section 2 now gives

$$f(a,b)=(a-1)b-a=ab-a-b,$$

Sylvester's formula of 1882; [2] for \(5\) and \(7\) it gives \(35-5-7=23\). Coprimality is essential: a common divisor \(d\gt1\) of the coins divides every representable number, and then there is no largest gap.

4. The coins 6, 10 and 15

For three or more coins no closed formula is known in general, [2] but products of three primes have a special structure. For \(p,q,r=2,3,5\) the coins are \(pq=6\), \(pr=10\) and \(qr=15\). Any two of them share a prime factor, while the three together have no common divisor greater than \(1\). A direct search finds fifteen gaps, \(1,2,3,4,5,7,8,9,11,13\), \(14\), \(17\), \(19\), \(23\) and \(29\), so \(f(6,10,15)=29\), as in the statement. Now sort by the remainder after division by the largest coin, \(m=qr=15\), and look at the fifteen sums \(\lambda\cdot10+\mu\cdot6\) with \(0\le\lambda\le2\) and \(0\le\mu\le4\):

  • \(\lambda=0\): \(0,6,12,18,24\), with remainders \(0,6,12,3,9\);
  • \(\lambda=1\): \(10,16,22,28,34\), with remainders \(10,1,7,13,4\);
  • \(\lambda=2\): \(20,26,32,38,44\), with remainders \(5,11,2,8,14\).

Every remainder from \(0\) to \(14\) appears exactly once, and a search confirms that each sum is the smallest representable number with its remainder. The largest sum is \(44\), so Selmer's formula gives \(f=44-15=29\).

5. The formula for every triple

The example is not a coincidence. Let \(p\), \(q\) and \(r\) be pairwise coprime numbers greater than \(1\), for instance three distinct primes. As in Section 3, coins can be exchanged: \(q\) coins \(pr\), \(r\) coins \(pq\) and \(p\) coins \(qr\) are all worth \(pqr\). This suggests looking at the amounts paid with fewer than \(q\) coins \(pr\), fewer than \(r\) coins \(pq\) and no coin \(qr\), that is, at the \(qr\) numbers

$$w(\lambda,\mu)=\lambda\,pr+\mu\,pq,$$

with \(0\le\lambda\lt q\) and \(0\le\mu\lt r\).

They have different remainders modulo \(qr\). Modulo \(q\) the term \(\mu\,pq\) vanishes, so \(w(\lambda,\mu)\equiv\lambda\,pr\pmod q\), and since \(\gcd(q,pr)=1\), the argument of Section 3 with Euclid's lemma shows that this remainder determines \(\lambda\). Modulo \(r\), where \(w(\lambda,\mu)\equiv\mu\,pq\) and \(\gcd(r,pq)=1\), the remainder determines \(\mu\). Two of the numbers with the same remainder modulo \(qr\) agree modulo \(q\) and modulo \(r\), so they are equal, and the \(qr\) numbers meet each of the \(qr\) remainders exactly once.

Each is the smallest representable number with its remainder. Let \(n=u\,pq+v\,pr+z\,qr\) with \(u,v,z\ge0\), and write \(v=\lambda+tq\) and \(u=\mu+sr\) with \(0\le\lambda\lt q\) and \(0\le\mu\lt r\). The exchanges give \(n=w(\lambda,\mu)+(z+pt+ps)\,qr\), so \(n\) has the same remainder modulo \(qr\) as \(w(\lambda,\mu)\), the only one of the \(qr\) numbers with this remainder, and \(n\ge w(\lambda,\mu)\).

So every remainder occurs, the largest gap exists, and the \(w(\lambda,\mu)\) are the numbers \(w_i\) of Section 2 for \(m=qr\). The largest of them is \(w(q-1,r-1)\), and Selmer's formula gives

$$\begin{aligned}f(pq,pr,qr)&=(q-1)pr+(r-1)pq-qr\\&=2pqr-pq-pr-qr.\end{aligned}$$

The value is symmetric in \(p\), \(q\) and \(r\): the proof with the letters renamed, sorting by \(pq\) or by \(pr\), gives the same result. Only pairwise coprimality was used, not primality. In the theory of numerical semigroups, the set \(S\) generated by \(pq\), \(pr\) and \(qr\) is a free semigroup, and the Apéry sets of free semigroups have exactly this shape. [5]

6. The checkpoints of the statement

For \(p,q,r=2,3,5\) the formula gives \(60-6-10-15=29\). For \(p,q,r=2,7,11\) the coins are \(14\), \(22\) and \(77\); the largest of the numbers \(w(\lambda,\mu)\) is \(w(6,10)=6\cdot22+10\cdot14=272\), so \(f=272-77=195\), in agreement with \(2\cdot154-14-22-77=195\). Both values match the statement.

7. Adding up over the prime triples

Below \(10\) there are four primes, \(2,3,5,7\), and four triples, with

  • \(f(6,10,15)\)\({}=60-6-10-15=29\),
  • \(f(6,14,21)\)\({}=84-6-14-21=43\),
  • \(f(10,14,35)\)\({}=140-10-14-35=81\),
  • \(f(15,21,35)\)\({}=210-15-21-35=139\),

so the sum for the limit \(10\) is \(292\). For the limit \(20\) the sum is \(68928\), over \(56\) triples, and for the limit \(100\) it is \(318227794\), over \(2300\) triples. Below \(5000\) there are \(669\) primes, so the program adds the formula over \(\binom{669}3=49\,679\,494\) triples.

8. A shortcut with symmetric sums

The sum can also be found without visiting the triples. Let \(n\) be the number of primes, \(e_2\) the sum of all products of two different primes and \(e_3\) the sum of all products of three different primes; these are elementary symmetric polynomials of the primes. [6] The terms \(2pqr\) add up to \(2e_3\). A pair of two primes lies in exactly \(n-2\) triples, one for each other prime, so the terms \(pq+pr+qr\) add up to \((n-2)\,e_2\), and

$$\sum_{p\lt q\lt r}f(pq,pr,qr)=2e_3-(n-2)\,e_2.$$

One pass over the primes computes these sums. Start from \(e_1=e_2=e_3=0\), where \(e_1\) becomes the sum of the primes, and let each new prime \(x\) replace \(e_3\) by \(e_3+x\,e_2\), then \(e_2\) by \(e_2+x\,e_1\), then \(e_1\) by \(e_1+x\): a new triple is \(x\) with an old pair and a new pair is \(x\) with an old prime, so each update uses the old values. This is the recursive definition of these polynomials. [6] For the limit \(10\), the primes \(2,3,5\) give \(e_1=10\), \(e_2=31\) and \(e_3=30\); then \(7\) gives \(e_3=30+7\cdot31=247\) and \(e_2=31+7\cdot10=101\), and with \(n=4\) indeed \(2\cdot247-2\cdot101=292\). This takes \(O(n)\) operations instead of about \(n^3/6\); the published programs do not use it.

How the Code Works

The C++ program [9] rests on two small functions. primes_below is the sieve of Eratosthenes: [7] in a std::vector<bool> it crosses out the multiples of every prime \(p\) with \(p^2\) below the limit, starting at \(p^2\), and then collects the numbers that remain. frobenius_semiprime_triple evaluates \(2pqr-pq-pr-qr\), the formula of Section 5, in unsigned 64-bit arithmetic. The formula grows with each of \(p\), \(q\) and \(r\), so its largest value here belongs to the three largest primes \(4987,4993,4999\); it is \(248\,876\,319\,707\), far below \(2^{64}\), and every intermediate result of the subtractions is positive.

main first calls parse_arguments, which accepts --skip-checkpoints and --prime-limit= followed by digits, read into an int without an overflow check; an unknown argument or a limit below \(3\) ends the program with exit code \(1\). Unless the checkpoints are skipped, run_checkpoints then compares the formula with the values \(29\) for \(2,3,5\) and \(195\) for \(2,7,11\) from the statement, and with brute_frobenius_three for all \(56\) triples of the primes \(2,3,5,\dots,19\). That function fills a table can with \(4abc+1\) entries for the coins \(a=pq\), \(b=pr\) and \(c=qr\): can[0] is true, and scanning \(x\) upwards, every marked \(x\) marks \(x+a\), \(x+b\) and \(x+c\). A representable \(x\gt0\) is a smaller representable number plus one coin, so it is marked before the scan reaches it. The function returns the largest unmarked entry. It does not check the length of the table, but in all \(56\) test triples at least \(3571\) marked entries follow that entry, far more than \(a\), so by Section 1 it is the exact largest gap. If a comparison fails, the program reports it and exits with code \(2\).

Then solve_sum sieves the primes below the limit and runs three nested loops over indices \(i\lt j\lt k\) of the prime list, so that every triple \(p\lt q\lt r\) of Section 7 is visited exactly once. It adds the values of the formula in a 128-bit unsigned integer, the type __uint128_t of GCC and Clang, [8] and to_string_u128 turns the total into decimal digits by repeated division by \(10\) before main prints it.

The Python and Java versions contain only the computation. Python builds the same sieve in a bytearray, keeps the product \(pq\) of the middle loop for the inner loop, and adds with integers of unlimited precision, so overflow cannot occur. Java uses a boolean[] sieve, an ArrayList<Integer> of primes and long arithmetic for the terms and for the total. A signed long holds numbers up to \(2^{63}-1\approx9.2\cdot10^{18}\). [8] The total is less than \(2e_3\le(\sum p)^3/3\), because the cube \((\sum p)^3\) contains \(6e_3\) and otherwise only nonnegative terms, and the \(669\) primes below \(5000\) add up to \(1\,548\,136\lt2\cdot10^6\); so the total is below \(2.7\cdot10^{18}\) and fits, and the 128-bit total of the C++ version is only a safety margin.

Complexity and Verification

For the limit \(L\) the sieve costs \(O(L\log\log L)\) operations. [7] The triple loop evaluates the formula \(\binom n3\approx n^3/6\) times for the \(n\) primes below \(L\), with a few multiplications each time, so it needs \(O(n^3)\) time; for \(L=5000\) these are the \(49\,679\,494\) evaluations of Section 7. The memory is \(O(L)\) for the sieve and the prime list; only the brute-force check needs more: its largest table, for \(13,17,19\), has \(70\,526\,405\) entries, which std::vector<bool> stores as bits, about \(8.8\) MB, although the largest gap there is only \(7607\). On a desktop computer (Clang with -O2, Java 17, Python 3.11) the full C++ run takes about one second, almost all of it in the checkpoints, whose \(56\) tables have \(291\,642\,676\) entries in total; the sum alone takes a few hundredths of a second in C++ with --skip-checkpoints, less than a tenth in Java and a few seconds in Python.

Independent checks, not part of the published programs, support every step. The formula \(f=2pqr-pq-pr-qr\) agrees with a direct search for the largest gap for all \(220\) triples of primes below \(40\), and with a shortest-path computation of the smallest representable number with each remainder for all \(15\,180\) triples of primes below \(200\) and all \(2518\) pairwise coprime triples \(2\le a\lt b\lt c\le40\), with the coins \(ab\), \(ac\) and \(bc\). For all \(680\) triples of primes below \(60\) the smallest representable numbers modulo \(qr\) are exactly the \(w(\lambda,\mu)\) of Section 5. The shortcut of Section 8 agrees with the triple loop for every limit from \(3\) to \(399\) and for the limit \(5000\), and all three published programs print the same number; with --prime-limit=10, 20 and 100 the C++ program prints \(292\), \(68928\) and \(318227794\).

Footnotes and References

The references below supply the background facts used above. The smallest representable numbers modulo \(qr\), the formula for \(f(pq,pr,qr)\) and the integer bounds are derived explicitly in this article.

  1. Project Euler 278 — Linear Combinations of Semiprimes. The official statement defines f(a₁, …, aₙ) as the largest b that cannot be written as q₁a₁ + … + qₙaₙ with whole numbers qₖ ≥ 0, gives f(5, 7) = 23, f(6, 10, 15) = 29 and f(14, 22, 77) = 195, and asks for the sum of f(pq, pr, qr) over all primes p < q < r < 5000.
  2. Wikipedia — Coin problem. The coin problem asks for the largest amount that cannot be paid with coins of given values, the Frobenius number, whose existence depends on the values having greatest common divisor 1. The article gives Sylvester's formula a₁a₂ − a₁ − a₂ of 1882 for two values and states that no closed-form solution is known for more than two.
  3. Wikipedia — Numerical semigroup. The Apéry set Ap(S, m) of a numerical semigroup S contains the smallest elements of S in each residue class modulo m, and the Frobenius number of S, its largest gap, is max Ap(S, m) − m.
  4. Wikipedia — Euclid's lemma. In its general form: if an integer n divides a product ab and is coprime with a, then n divides b. Sections 3 and 5 use it to cancel a factor that is coprime to the modulus.
  5. A. Assi and P. A. García-Sánchez — Numerical semigroups and applications (arXiv:1411.6093). Proposition 12 gives Selmer's formula F(S) = max Ap(S, n) − n. Section 2.2 defines free semigroups through an arrangement (a₀, …, aₕ) of the generators, and Corollary 60 shows that the Apéry set of a₀ is then the set of sums λ₁a₁ + … + λₕaₕ with 0 ≤ λₖ < eₖ; for the free arrangement (qr, pr, pq), with e₁ = q and e₂ = r, these are the w(λ, μ) of Section 5.
  6. Wikipedia — Elementary symmetric polynomial. The polynomial eₖ is the sum of all products of k distinct variables, and its recursive definition eₖ(X₁, …, Xₙ) = eₖ(X₁, …, Xₙ₋₁) + Xₙ·eₖ₋₁(X₁, …, Xₙ₋₁) is the one-pass update of Section 8.
  7. Wikipedia — Sieve of Eratosthenes. The sieve crosses out the multiples of each prime; it may start at p², because all smaller multiples of p have already been marked, it may stop once p² exceeds the limit, and it finds all primes below n in O(n log log n) operations.
  8. GCC manual — 128-bit Integers. GCC provides the types __int128 and unsigned __int128 on targets that support them; __uint128_t in the C++ program is the unsigned one. The Java Tutorials — Primitive Data Types: long is a 64-bit integer with maximum 2⁶³ − 1; Python — Numeric Types: integers have unlimited precision.
  9. C++ — Euler278.cpp. The linked immutable C++ revision contains primes_below, frobenius_semiprime_triple, solve_sum, brute_frobenius_three and run_checkpoints. It is the implementation record for this article.

Problem 278 source code

C++

#include <algorithm>
#include <array>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

struct Options {
    int prime_limit = 5000;
    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, "--prime-limit=", options.prime_limit)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.prime_limit >= 3;
}

std::vector<int> primes_below(const int limit) {
    std::vector<bool> is_prime(static_cast<std::size_t>(limit), true);
    if (limit > 0) {
        is_prime[0] = false;
    }
    if (limit > 1) {
        is_prime[1] = false;
    }

    for (int p = 2; p * p < limit; ++p) {
        if (!is_prime[static_cast<std::size_t>(p)]) {
            continue;
        }
        for (int x = p * p; x < limit; x += p) {
            is_prime[static_cast<std::size_t>(x)] = false;
        }
    }

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

u64 frobenius_semiprime_triple(const u64 p, const u64 q, const u64 r) {
    return 2ULL * p * q * r - p * q - p * r - q * r;
}

std::string to_string_u128(u128 x) {
    if (x == 0) {
        return "0";
    }
    std::string s;
    while (x > 0) {
        const int digit = static_cast<int>(x % 10);
        s.push_back(static_cast<char>('0' + digit));
        x /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

u128 solve_sum(const int prime_limit) {
    const std::vector<int> primes = primes_below(prime_limit);

    u128 answer = 0;
    const int n = static_cast<int>(primes.size());
    for (int i = 0; i < n; ++i) {
        const u64 p = static_cast<u64>(primes[static_cast<std::size_t>(i)]);
        for (int j = i + 1; j < n; ++j) {
            const u64 q = static_cast<u64>(primes[static_cast<std::size_t>(j)]);
            for (int k = j + 1; k < n; ++k) {
                const u64 r = static_cast<u64>(primes[static_cast<std::size_t>(k)]);
                answer += static_cast<u128>(frobenius_semiprime_triple(p, q, r));
            }
        }
    }
    return answer;
}

u64 brute_frobenius_three(const u64 a, const u64 b, const u64 c) {
    const u64 cap = 4 * a * b * c;
    std::vector<bool> can(static_cast<std::size_t>(cap + 1), false);
    can[0] = true;

    for (u64 x = 0; x <= cap; ++x) {
        if (!can[static_cast<std::size_t>(x)]) {
            continue;
        }
        if (x + a <= cap) {
            can[static_cast<std::size_t>(x + a)] = true;
        }
        if (x + b <= cap) {
            can[static_cast<std::size_t>(x + b)] = true;
        }
        if (x + c <= cap) {
            can[static_cast<std::size_t>(x + c)] = true;
        }
    }

    u64 best = 0;
    for (u64 x = 1; x <= cap; ++x) {
        if (!can[static_cast<std::size_t>(x)]) {
            best = x;
        }
    }
    return best;
}

bool run_checkpoints() {
    if (frobenius_semiprime_triple(2, 3, 5) != 29ULL) {
        std::cerr << "Checkpoint failed for (2,3,5)" << '\n';
        return false;
    }
    if (frobenius_semiprime_triple(2, 7, 11) != 195ULL) {
        std::cerr << "Checkpoint failed for (2,7,11)" << '\n';
        return false;
    }

    // Verify formula against brute force on a small prime pool.
    const std::vector<int> small_primes = primes_below(20);
    for (std::size_t i = 0; i < small_primes.size(); ++i) {
        for (std::size_t j = i + 1; j < small_primes.size(); ++j) {
            for (std::size_t k = j + 1; k < small_primes.size(); ++k) {
                const u64 p = static_cast<u64>(small_primes[i]);
                const u64 q = static_cast<u64>(small_primes[j]);
                const u64 r = static_cast<u64>(small_primes[k]);
                const u64 a = p * q;
                const u64 b = p * r;
                const u64 c = q * r;
                const u64 brute = brute_frobenius_three(a, b, c);
                const u64 formula = frobenius_semiprime_triple(p, q, r);
                if (brute != formula) {
                    std::cerr << "Formula mismatch for (" << p << ',' << q << ',' << r
                              << "): brute=" << brute << ", formula=" << formula << '\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;
    }

    std::cout << to_string_u128(solve_sum(options.prime_limit)) << '\n';
    return 0;
}

Python

def solve():
    prime_limit = 5000

    def primes_below(limit):
        is_p = bytearray(b'\x01' * limit)
        is_p[0] = 0
        is_p[1] = 0
        for p in range(2, limit):
            if p * p >= limit:
                break
            if is_p[p]:
                for q in range(p*p, limit, p):
                    is_p[q] = 0
        return [p for p in range(2, limit) if is_p[p]]

    primes = primes_below(prime_limit)
    n = len(primes)
    answer = 0
    for i in range(n):
        p = primes[i]
        for j in range(i + 1, n):
            q = primes[j]
            pq = p * q
            for k in range(j + 1, n):
                r = primes[k]
                answer += 2 * pq * r - pq - p * r - q * r
    return str(answer)

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

Java

import java.util.*;

public class Euler278 {
    static List<Integer> primesBelow(int limit) {
        boolean[] isPrime = new boolean[limit];
        Arrays.fill(isPrime, true);
        if (limit > 0)
            isPrime[0] = false;
        if (limit > 1)
            isPrime[1] = false;

        for (int p = 2; p * p < limit; ++p) {
            if (!isPrime[p])
                continue;
            for (int x = p * p; x < limit; x += p) {
                isPrime[x] = false;
            }
        }

        List<Integer> primes = new ArrayList<>();
        for (int p = 2; p < limit; ++p) {
            if (isPrime[p]) {
                primes.add(p);
            }
        }
        return primes;
    }

    static long frobeniusSemiprimeTriple(long p, long q, long r) {
        return 2L * p * q * r - p * q - p * r - q * r;
    }

    public static String solve() {
        int primeLimit = 5000;
        List<Integer> primes = primesBelow(primeLimit);
        int n = primes.size();

        // 128-bit emulation for summing large numbers is necessary if the sum exceeds
        // Long.MAX_VALUE.
        // However, max prime is 5000.
        // There are ~669 primes. Number of combinations is (669 choose 3) ~ 50,000,000
        // Average value of frobenius is 2 * p * q * r where p,q,r <= 5000.
        // Max value is 2 * 5000^3 = 2.5e11
        // Expected sum: 2.5e11 * 5e7 = 1.25e19, which fits in an unsigned 64-bit int,
        // but Long.MAX_VALUE is 9e18. So it will overflow a signed long!
        // Wait, the output from C++ is 1228215747273908452.
        // Does 1228215747273908452 fit in a signed long?
        // 1.22e18 fits safely in a signed 64-bit int (Long.MAX_VALUE is 9.22e18).

        long answer = 0;
        for (int i = 0; i < n; ++i) {
            long p = primes.get(i);
            for (int j = i + 1; j < n; ++j) {
                long q = primes.get(j);
                for (int k = j + 1; k < n; ++k) {
                    long r = primes.get(k);
                    answer += frobeniusSemiprimeTriple(p, q, r);
                }
            }
        }
        return String.valueOf(answer);
    }

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