Problem 146: Investigating a Prime Pattern

Project Euler Problem 146 Solution

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

Problem Summary We seek the sum of all integers $n < 150{,}000{,}000$ for which the six values $$n^2+1,\ n^2+3,\ n^2+7,\ n^2+9,\ n^2+13,\ n^2+27$$ are prime, while the odd values between them, $$n^2+5,\ n^2+11,\ n^2+15,\ n^2+17,\ n^2+19,\ n^2+21,\ n^2+23,\ n^2+25,$$ are composite. Because the primes must be consecutive inside that window, those eight exclusions are part of the definition of a valid $n$, not an optional afterthought. A brute-force search would examine far too many large numbers near $2.25\times 10^{16}$. The successful strategy is to force $n$ into a tiny set of residue classes, precompute modular obstructions for many small primes, and reserve exact primality testing for the few candidates that survive all cheap filters. Mathematical Approach Let $$G=\{1,3,7,9,13,27\},\qquad C=\{5,11,15,17,19,21,23,25\}.$$ Then $n$ is a solution exactly when $n^2+g$ is prime for every $g\in G$ and $n^2+c$ is composite for every $c\in C$. The required and forbidden offsets If $n$ is even, then $n^2$ is even, so every even offset produces an even number greater than 2 and is automatically composite. That is why only the odd offsets matter. The implementation therefore tracks two explicit sets: the six offsets that must produce primes and the eight odd offsets in the same interval that must not....

Detailed mathematical approach

Problem Summary

We look for positive integers $n$ for which the six numbers $n^2+1$, $n^2+3$, $n^2+7$, $n^2+9$, $n^2+13$ and $n^2+27$ are consecutive primes: all six are prime, and no other prime lies between $n^2+1$ and $n^2+27$. The smallest such $n$ is $10$, and the sum of all such $n$ below one million is $1242490$. We must find the sum of all such $n$ below $150$ million. ^[1]

The whole solution in one sentence: divisibility by small primes forces $n$ to be a multiple of $10$ of a special form, makes seven of the eight odd numbers between the six primes composite automatically and cheaply rules out almost all remaining candidates, so that an exact primality test is needed for only a few thousand values of $n$.

Mathematical Approach

1. What "consecutive" demands

For $n=10$ the window runs from $n^2+1=101$ to $n^2+27=127$. The six required numbers $101$, $103$, $107$, $109$, $113$ and $127$ are prime. The even numbers of the window are composite, and so are its other odd numbers: $105=3\cdot5\cdot7$, $111=3\cdot37$, $115=5\cdot23$, $117=3^2\cdot13$, $119=7\cdot17$, $121=11^2$, $123=3\cdot41$ and $125=5^3$. So the six primes are consecutive, and $n=10$ counts.

The word "consecutive" matters: for $n=2$ the six numbers $5$, $7$, $11$, $13$, $17$ and $31$ are all prime, but $19$, $23$ and $29$, which lie between $17$ and $31$, are prime too, so $n=2$ does not count. In general, with the set of required offsets

$$G=\{1,3,7,9,13,27\},$$

a positive integer $n$ is valid exactly when $n^2+k$ is prime for every $k\in G$ and composite for every other $k$ with $1\le k\le27$.

2. Parity

If $n$ is odd, $n^2+1$ is even, so it can only be prime if it equals $2$, that is for $n=1$; but then $n^2+3=4$ is not prime. Hence $n$ is even, and then every $n^2+k$ with an even $k$ is an even number larger than $2$, composite without any test. The demand "no other prime" is left with the eight odd offsets

$$C=\{5,11,15,17,19,21,23,25\}:$$

for a valid $n$, every number $n^2+c$ with $c\in C$ is composite.

3. Squares modulo 5

Whether a prime $p$ divides $n^2+g$ depends only on the remainder of $n^2$ modulo $p$: $p\mid n^2+g$ exactly when $n^2\equiv-g\pmod p$. The remainders that squares leave modulo $p$ are called quadratic residues. ^[2] Modulo $5$ the squares of $0,1,2,3,4$ leave $0,1,4,4,1$, so there are three cases:

if $n^2\equiv1\pmod5$, then $5$ divides $n^2+9$, which exceeds $5$ and is therefore composite;
if $n^2\equiv4\pmod5$, then $5$ divides $n^2+1$, which is prime only if $n^2+1=5$, that is for $n=2$, excluded in Section 1;
otherwise $n^2\equiv0\pmod5$, and since $5$ is prime, $5\mid n$.

Together with Section 2, every valid $n$ is a multiple of $10$.

4. The primes 3 and 7

Modulo $3$ a square leaves $0$ or $1$. If $3\mid n$, then $3$ divides $n^2+3$, which is larger than $3$; so $3\nmid n$ and $n^2\equiv1\pmod3$. Modulo $7$ the squares of $0,1,\dots,6$ leave $0,1,4,2,2,4,1$, and three of the four possible remainders are fatal:

$n^2\equiv0$: $7$ divides $n^2+7$;
$n^2\equiv1$: $7$ divides $n^2+13$ and $n^2+27$;
$n^2\equiv4$: $7$ divides $n^2+3$, which equals $7$ only for $n=2$.

Hence $n^2\equiv2\pmod7$, which holds exactly for $n\equiv3$ and $n\equiv4\pmod7$. In this case the six required numbers leave the remainders $3,5,2,4,1,1$ modulo $7$, so $7$ divides none of them. Similarly, for a multiple of $10$ with $3\nmid n$ none of the six numbers is divisible by $2$, $3$ or $5$.

5. Two arithmetic progressions

The conditions $n\equiv0\pmod{10}$ and $n\equiv3$ or $4\pmod7$ concern the coprime moduli $10$ and $7$, so by the Chinese remainder theorem each of the two combinations has exactly one solution modulo $70$. ^[3] The multiples $0,10,\dots,60$ of $10$ leave the remainders $0,3,6,2,5,1,4$ modulo $7$, so every valid $n$ satisfies $n\equiv10\pmod{70}$ or $n\equiv60\pmod{70}$; the progressions $10,80,150,\dots$ and $60,130,200,\dots$ contain $2$ of every $70$ integers. The condition $3\nmid n$ removes a third of them.

6. The numbers in between

Claim. If $10\mid n$, $3\nmid n$ and $n\equiv\pm3\pmod7$, then $n^2+c$ is composite for every $c\in C$ other than $21$.

Proof. The hypotheses give $n^2\equiv1\pmod3$, $n^2\equiv0\pmod5$ and $n^2\equiv2\pmod7$. Hence $3$ divides $n^2+c$ for $c=5,11,17,23$, because $1+c$ is a multiple of $3$; $5$ divides $n^2+15$ and $n^2+25$; and $7$ divides $n^2+19$, because $2+19=21$. All these numbers exceed $7$, so they are composite. For $c=21$ none of the three primes helps, since $1+21=22$, $0+21=21$ and $2+21=23$ are not divisible by $3$, $5$ and $7$, respectively. $\square$

So for the candidates of Section 5 with $3\nmid n$, "consecutive" adds a single condition to the six primality tests: $n^2+21$ must be composite. It matters: below $150$ million exactly two values of $n$ give six prime numbers that are not consecutive primes, $n=2$ from Section 1 and $n=144774340$, for which $n^2+21$ is also prime.

7. The prime 13

Modulo $13$ the squares leave $0,1,3,4,9,10,12$, and the six offsets leave $1,3,7,9,0,1$. A required number $n^2+g$ is divisible by $13$ when $n^2\equiv-g\pmod{13}$, that is when $n^2$ leaves $12$, $10$, $4$ or $0$; the offset $7$ is harmless, because $-7\equiv6$ is not a square modulo $13$. Only $n^2\equiv1,3,9$ survive, which means $n\equiv\pm1,\pm4,\pm3\pmod{13}$, six of the thirteen remainders. The programs test $13\mid n$ directly and leave the other forbidden remainders to the tables of the next section.

8. Tables of forbidden remainders

The same reasoning works for every prime $p$. Call a remainder $r$ forbidden modulo $p$ if $p$ divides $r^2+g$ for some $g\in G$:

$$R_p=\{\,r\in\{0,1,\dots,p-1\}\ :\ \exists\,g\in G,\ p\mid r^2+g\,\}.$$

For example $R_{11}=\{2,3,8,9\}$, because $2^2+7=11$, $3^2+13=22$, $8^2+13=77$ and $9^2+7=88$. If $n$ leaves a forbidden remainder modulo $p$, then $p$ divides one of the six numbers, which is therefore composite as long as it is larger than $p$. Modulo a prime a number has at most two square roots, ^[2] so $R_p$ contains at most $12$ of the $p$ possible remainders.

The proviso "larger than $p$" matters for small $n$: for $n=10$ the number $n^2+1=101$ is itself a prime below $2000$, and $10\in R_{101}$, so the table for $101$ would throw away the smallest solution. For $n\ge45$ every $n^2+g$ exceeds $2000$. The programs use the tables only from $n=1000$ on, a safe round bound, and send the few smaller candidates straight to the exact tests.

Over a full period of all these moduli the Chinese remainder theorem makes the remainders modulo different primes independent, so the surviving fraction of the progressions is the product of $2/3$ for $p=3$, $1$ for $p=7$ and $1-|R_p|/p$ for $11\le p\le2000$, about $0.00058$. Our range is far shorter than that period, yet this predicts $2471$ survivors between $1000$ and $150$ million, close to the true $2541$. The same factors appear in the constant of the Bateman–Horn conjecture, ^[4] which predicts infinitely many $n$ with six prime values: each $n^2+g$ is irreducible, having no real root, and $|R_p|\lt p$ for every prime $p$, by Sections 2 to 4 and the example $R_{11}$ for $p\le11$ and by $|R_p|\le12$ for larger $p$. The prediction is unproved; even the infinitude of primes of the form $n^2+1$ is an open problem of Landau. ^[5]

9. The checkpoint below one million

Below $10^6$ the two progressions contain $28571$ numbers. The tests $3\nmid n$ and $13\nmid n$ keep $17582$ of them, $18$ of which are below $1000$; the tables remove all but $26$ of the other $17564$, so only $44$ candidates reach the exact tests. Exactly three of them make all six numbers prime, $n=10$, $315410$ and $927070$, and for each of them $n^2+21$ is composite. Their sum $10+315410+927070=1242490$ is the value given in the statement.

10. Exact primality tests

Since $n\le149999990$, every number to be tested is below $2.25\cdot10^{16}\lt2^{55}$; trial division would need divisors up to $1.5\cdot10^8$. The Miller–Rabin test ^[6] is much faster. For an odd number $m$ write $m-1=2^s d$ with $d$ odd and choose a base $a$. The number $m$ passes for the base $a$ if $a^d\equiv1\pmod m$, or if $a^{2^r d}\equiv-1\pmod m$ for some $r$ with $0\le r\lt s$. Every prime $m$ passes for every base $a$ that it does not divide. Indeed, in the sequence

$$a^{d},\ a^{2d},\ a^{4d},\ \dots,\ a^{2^{s}d}=a^{m-1}$$

each number is the square of the one before, and the last one is $\equiv1\pmod m$ by Fermat's little theorem. Modulo a prime the only square roots of $1$ are $1$ and $-1$, so either $a^d\equiv1$, or the number just before the first $1$ is $\equiv-1$. For example $m=101$ has $100=2^2\cdot25$, and for the base $2$ we get $2^{25}\equiv10$ and $2^{50}\equiv10^2\equiv-1\pmod{101}$, so $101$ passes.

A composite number can pass for some bases. The smallest example for the base $2$ is $2047=23\cdot89$: here $2046=2\cdot1023$, and $2^{11}=2048\equiv1$ gives $2^{1023}\equiv1\pmod{2047}$. The base $325$ exposes it, since $325^{1023}\equiv967\pmod{2047}$, which is neither $1$ nor $2046\equiv-1$. For every odd composite number at most a quarter of the bases let it pass, ^[6] and for the numbers below $2^{64}$ Jim Sinclair found seven bases, $2$, $325$, $9375$, $28178$, $450775$, $9780504$ and $1795265022$, such that no composite number below $2^{64}$ passes for all seven. ^[7] With these bases the test is exact for all numbers in this problem.

How the Code Works

The C++ program ^[10] first runs run_checkpoints, unless --skip-checkpoints is given: the primality test must be right for $1$, $2$, $3$ and $21$, and the sum below one million must be $1242490$. Then build_filters(2000) builds the tables and solve searches below the limit, $150$ million unless an option such as --limit=1000000 sets another one. Python has the same functions and options. Java has neither the checkpoint nor the options; its routines are sievePrimes, buildFilters, runSequence and isPrime, and it keeps $n$, its square and the sum in long.

sieve_primes lists the $303$ primes up to $2000$ with the sieve of Eratosthenes. ^[8] build_filters skips $2$ and $5$, whose tables could never reject a multiple of $10$ (Section 4), and stores for each of the other $301$ primes an array bad of length $p$ that marks $R_p$, and step $=10\bmod p$.

run_sequence walks along one progression of Section 5 from the first element that first_with_mod finds. It keeps square $=n^2$ and, for every table, the remainder residues[i] of $n$ modulo $p_i$ and the step step70[i] $=70\bmod p_i$. For each $n$ it requires that $3$ and $13$ do not divide square, that is $3\nmid n$ and $13\nmid n$; from $n=1000$ on (filter_threshold) it looks up the remainders in the tables, in increasing order of $p$, until one is forbidden. A surviving candidate is added to the sum if the six numbers at the offsets in kGood are prime and the eight at the offsets in kBad are not; each check stops at the first failure, and by Section 6, is_prime recognises seven of the eight kBad numbers as composite at once by trial division. Then $n$ grows by $70$ and square by $140(n-70)+4900$, where $n-70$ is the previous value, since $(m+70)^2=m^2+140m+4900$; every remainder grows by its step, minus $p_i$ when it reaches $p_i$, so the remainders are kept up to date without any division. solve runs both progressions and adds the results.

is_prime rejects numbers below $2$, divides by the twelve primes up to $37$ and runs the Miller–Rabin test of Section 10 with the seven bases, skipping a base that the number divides, as recommended in ^[7]. pow_mod computes powers by repeated squaring, and mul_mod multiplies two remainders below $2^{55}$ in unsigned __int128, as such a product has up to $110$ bits; everything else fits in the 64-bit type u64. ^[9] Python runs the same code with unbounded integers. Java keeps the trial division but then calls BigInteger.isProbablePrime(20), a probabilistic test whose positive answers are documented to be correct with probability above $1-2^{-20}$, ^[9] unlike the deterministic test of the other two programs; at this limit it gives the same total.

Complexity and Verification

Below the limit $L=1.5\cdot10^8$ the progressions contain $4285714\approx2L/70$ numbers, and for each of them the loop updates all $301$ remainders: $1289999914$ additions, the dominant cost. The tests for $3$ and $13$ pass $2637362$ candidates, and the $2637344$ of them from $1000$ on need $15235513$ table lookups, about $5.8$ each: the table for $11$ alone rejects $959035$ of them, while those for $3$ and $7$, looked up first, never reject. Only $2559$ candidates reach the primality tests. In C++ and Python they cause $4233$ calls of is_prime, and $4127$ of these reach the Miller–Rabin stage, with at most seven modular powers of about $55$ squarings each. With $F=301$ tables the time is thus $O(LF)$ plus $O(\log L)$ multiplications per base and call. The memory is dominated by the $277043$ table entries, one byte each in C++, filled with at most $6\cdot277043$ offset checks; it does not grow with $L$. On a desktop computer the C++ and Java programs finish in well under a second, the Python program in under a minute.

By default the C++ and Python programs check the value $1242490$ from the statement before the main computation. For this article the result was also confirmed by checks that are not part of the published solutions:

a brute force without modular shortcuts, testing all $27$ numbers of the window with the isprime function of SymPy for every $n$ below $150$ million with $n^2+1$ prime, gives the same sum as the C++ program;
the three programs and a separate reimplementation of the search with NumPy and SymPy print the same total;
the published is_prime agrees with a sieve up to $3\cdot10^6$ and with SymPy, which uses other deterministic base sets at this size, on $340000$ random numbers of the sizes that occur, and it rejects strong pseudoprimes such as $2047$ and $3825123056546413051$.

Footnotes and References

The references below supply the background facts used above. The two progressions, the claim of Section 6, the tables $R_p$ and all counts are derived or computed explicitly in this article.

Project Euler 146 — Investigating a Prime Pattern. The official statement asks for the sum of all n below 150 million for which n² + 1, n² + 3, n² + 7, n² + 9, n² + 13 and n² + 27 are consecutive primes, and gives the smallest such n, 10, and the sum 1242490 below one million.
Wikipedia — Quadratic residue. A quadratic residue modulo n is a number congruent to a perfect square modulo n. Modulo an odd prime p a number a has 1 + (a|p) square roots, never more than two; the article's table lists the residues modulo 5, 7 and 13.
Wikipedia — Chinese remainder theorem. If the moduli are pairwise coprime, the congruences x ≡ aᵢ (mod nᵢ) have exactly one solution x with 0 ≤ x < N, where N is the product of the moduli.
Wikipedia — Bateman–Horn conjecture. For irreducible integer polynomials whose product f has no fixed prime divisor, it predicts how many n below x make all of them prime, with a constant built from the numbers N(p) of solutions of f(n) ≡ 0 (mod p).
Wikipedia — Landau's problems. The fourth problem asks whether there are infinitely many primes of the form n² + 1; all four problems are unresolved.
Wikipedia — Miller–Rabin primality test. With n − 1 = 2ˢ·d and d odd, n is a strong probable prime to the base a if a^d ≡ 1 or a^(2ʳ·d) ≡ −1 (mod n) for some 0 ≤ r < s. Primes always pass, and for an odd composite n at most a quarter of the bases are strong liars.
W. Izykowski — Deterministic variants of the Miller–Rabin primality test. The table of the best known base sets lists 2, 325, 9375, 28178, 450775, 9780504 and 1795265022, found by Jim Sinclair in 2011, with the bound 2⁶⁴; a base that is 0 modulo n is to count as passed.
Wikipedia — Sieve of Eratosthenes. The sieve marks the multiples of each prime as composite, starting at the square of the prime; the numbers left unmarked are the primes.
GCC manual — 128-bit Integers documents unsigned __int128. Java — BigInteger.isProbablePrime returns false only for definitely composite numbers; when it returns true, the number is prime with probability above 1 − 1/2^certainty.
C++ — Euler146.cpp. The linked immutable C++ revision contains is_prime, build_filters, run_sequence, solve and run_checkpoints. It is the implementation record for this article.

Problem 146 source code

C++

#include <array>
#include <algorithm>
#include <cstdint>
#include <iostream>
#include <limits>
#include <pthread.h>
#include <string>
#include <unistd.h>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

struct Options {
    u64 limit = 150000000ULL;
    bool run_checkpoints = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    u64 parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10ULL + static_cast<u64>(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 (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }
        u64 thread_count_value = 0;
        if (parse_u64_after_prefix(arg, "--threads=", thread_count_value)) {
            if (thread_count_value > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
                std::cerr << "Thread count too large: " << thread_count_value << '\n';
                return false;
            }
            options.requested_threads = static_cast<unsigned>(thread_count_value);
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 10;
}

u64 mul_mod(const u64 a, const u64 b, const u64 mod) {
    return static_cast<u64>((static_cast<u128>(a) * static_cast<u128>(b)) % mod);
}

u64 pow_mod(u64 base, u64 exp, const u64 mod) {
    base %= mod;
    u64 result = 1 % mod;
    while (exp > 0) {
        if ((exp & 1ULL) != 0ULL) {
            result = mul_mod(result, base, mod);
        }
        base = mul_mod(base, base, mod);
        exp >>= 1ULL;
    }
    return result;
}

bool is_prime(const u64 n) {
    if (n < 2) {
        return false;
    }
    for (u64 p : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
        if (n == p) {
            return true;
        }
        if (n % p == 0ULL) {
            return false;
        }
    }

    u64 d = n - 1;
    int s = 0;
    while ((d & 1ULL) == 0ULL) {
        d >>= 1ULL;
        ++s;
    }

    static constexpr std::array<u64, 7> kBases{
        2ULL, 325ULL, 9375ULL, 28178ULL, 450775ULL, 9780504ULL, 1795265022ULL};

    for (u64 a : kBases) {
        if (a % n == 0ULL) {
            continue;
        }
        u64 x = pow_mod(a, d, n);
        if (x == 1ULL || x == n - 1) {
            continue;
        }
        bool witness = true;
        for (int r = 1; r < s; ++r) {
            x = mul_mod(x, x, n);
            if (x == n - 1) {
                witness = false;
                break;
            }
        }
        if (witness) {
            return false;
        }
    }
    return true;
}

std::vector<int> sieve_primes(int n) {
    std::vector<bool> is_prime(static_cast<std::size_t>(n + 1), true);
    if (n >= 0) is_prime[0] = false;
    if (n >= 1) is_prime[1] = false;
    for (int p = 2; static_cast<long long>(p) * p <= n; ++p) {
        if (!is_prime[static_cast<std::size_t>(p)]) continue;
        for (int q = p * p; q <= n; q += p) {
            is_prime[static_cast<std::size_t>(q)] = false;
        }
    }
    std::vector<int> primes;
    for (int i = 2; i <= n; ++i) {
        if (is_prime[static_cast<std::size_t>(i)]) primes.push_back(i);
    }
    return primes;
}

struct ModFilter {
    int p = 0;
    int step = 0;
    std::vector<std::uint8_t> bad;
};

std::vector<ModFilter> build_filters(int limit) {
    static constexpr std::array<u64, 6> kGood{1ULL, 3ULL, 7ULL, 9ULL, 13ULL, 27ULL};
    const auto primes = sieve_primes(limit);
    std::vector<ModFilter> filters;
    filters.reserve(primes.size());
    for (int p : primes) {
        if (p == 2 || p == 5) {
            continue;
        }
        ModFilter f;
        f.p = p;
        f.step = 10 % p;
        f.bad.assign(static_cast<std::size_t>(p), 0);
        for (int r = 0; r < p; ++r) {
            const int r2 = static_cast<int>((1LL * r * r) % p);
            bool bad = false;
            for (u64 k : kGood) {
                if ((r2 + static_cast<int>(k % p)) % p == 0) {
                    bad = true;
                    break;
                }
            }
            f.bad[static_cast<std::size_t>(r)] = bad ? 1 : 0;
        }
        filters.push_back(std::move(f));
    }
    return filters;
}

unsigned thread_count(bool allow_multithreading, unsigned requested) {
    if (!allow_multithreading) {
        return 1U;
    }
    if (requested != 0U) {
        return requested;
    }
    long n = sysconf(_SC_NPROCESSORS_ONLN);
    if (n < 1) n = 1;
    return static_cast<unsigned>(n);
}

struct Task {
    u64 start = 0;
    u64 end = 0;
    const std::vector<ModFilter>* filters = nullptr;
    u64 sum = 0;
};

u64 first_with_mod(u64 start, const u64 mod, const u64 residue) {
    const u64 current = start % mod;
    if (current <= residue) {
        return start + (residue - current);
    }
    return start + (mod - (current - residue));
}

void run_sequence(u64 start, const u64 end, const u64 residue, const std::vector<ModFilter>& filters, u64& sum) {
    static constexpr std::array<u64, 6> kGood{1ULL, 3ULL, 7ULL, 9ULL, 13ULL, 27ULL};
    static constexpr std::array<u64, 8> kBad{5ULL, 11ULL, 15ULL, 17ULL, 19ULL, 21ULL, 23ULL, 25ULL};

    constexpr u64 filter_threshold = 1000ULL;
    constexpr u64 step = 70ULL;

    u64 n = first_with_mod(start, step, residue);
    if (n >= end) {
        return;
    }

    u64 square = n * n;

    std::vector<int> residues(filters.size(), 0);
    std::vector<int> step70(filters.size(), 0);
    for (std::size_t i = 0; i < filters.size(); ++i) {
        const int p = filters[i].p;
        residues[i] = static_cast<int>(n % static_cast<u64>(p));
        step70[i] = (filters[i].step * 7) % p;
    }

    while (n < end) {
        if (square % 3ULL != 0ULL && square % 13ULL != 0ULL) {
            bool ok = true;
            if (n >= filter_threshold) {
                for (std::size_t i = 0; i < filters.size(); ++i) {
                    if (filters[i].bad[static_cast<std::size_t>(residues[i])]) {
                        ok = false;
                        break;
                    }
                }
            }
            if (ok) {
                for (u64 k : kGood) {
                    if (!is_prime(square + k)) {
                        ok = false;
                        break;
                    }
                }
                if (ok) {
                    for (u64 k : kBad) {
                        if (is_prime(square + k)) {
                            ok = false;
                            break;
                        }
                    }
                }
                if (ok) {
                    sum += n;
                }
            }
        }

        n += step;
        square += 140ULL * (n - step) + 4900ULL;
        for (std::size_t i = 0; i < filters.size(); ++i) {
            const int p = filters[i].p;
            int next = residues[i] + step70[i];
            if (next >= p) next -= p;
            residues[i] = next;
        }
    }
}

static void* worker_fn(void* arg) {
    auto* task = static_cast<Task*>(arg);
    const auto& filters = *task->filters;

    u64 sum = 0;
    run_sequence(task->start, task->end, 10ULL, filters, sum);
    run_sequence(task->start, task->end, 60ULL, filters, sum);
    task->sum = sum;
    return nullptr;
}

u64 solve(const u64 limit, const std::vector<ModFilter>& filters, bool allow_multithreading, unsigned requested_threads) {
    const unsigned threads = thread_count(allow_multithreading, requested_threads);
    if (threads <= 1) {
        Task single{10, limit, &filters, 0};
        worker_fn(&single);
        return single.sum;
    }

    const u64 start = 10;
    const u64 count = (limit - start + 9) / 10;
    const u64 block = (count + threads - 1) / threads;

    std::vector<Task> tasks(threads);
    std::vector<pthread_t> workers(threads);
    for (unsigned t = 0; t < threads; ++t) {
        const u64 idx_start = static_cast<u64>(t) * block;
        const u64 idx_end = std::min(count, idx_start + block);
        const u64 n_start = start + idx_start * 10ULL;
        const u64 n_end = start + idx_end * 10ULL;
        tasks[t] = {n_start, n_end, &filters, 0};
        pthread_create(&workers[t], nullptr, worker_fn, &tasks[t]);
    }
    u64 sum = 0;
    for (unsigned t = 0; t < threads; ++t) {
        pthread_join(workers[t], nullptr);
        sum += tasks[t].sum;
    }
    return sum;
}

bool run_checkpoints() {
    if (!is_prime(2ULL) || !is_prime(3ULL) || is_prime(1ULL) || is_prime(21ULL)) {
        std::cerr << "Checkpoint failed for primality tester" << '\n';
        return false;
    }
    const auto filters = build_filters(2000);
    if (solve(1000000ULL, filters, false, 1) != 1242490ULL) {
        std::cerr << "Checkpoint failed for limit 1,000,000" << '\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;
    }
    const auto filters = build_filters(2000);
    std::cout << solve(options.limit, filters, options.allow_multithreading, options.requested_threads) << '\n';
    return 0;
}

Python

import sys
import math
from typing import List, Tuple

def mul_mod(a: int, b: int, mod: int) -> int:
    return (a * b) % mod

def pow_mod(base: int, exp: int, mod: int) -> int:
    base %= mod
    result = 1 % mod
    while exp > 0:
        if (exp & 1) != 0:
            result = mul_mod(result, base, mod)
        base = mul_mod(base, base, mod)
        exp >>= 1
    return result

def is_prime(n: int) -> bool:
    if n < 2:
        return False
    small_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
    for p in small_primes:
        if n == p:
            return True
        if n % p == 0:
            return False
    
    d = n - 1
    s = 0
    while (d & 1) == 0:
        d >>= 1
        s += 1
    
    bases = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]
    
    for a in bases:
        if a % n == 0:
            continue
        x = pow_mod(a, d, n)
        if x == 1 or x == n - 1:
            continue
        witness = True
        for r in range(1, s):
            x = mul_mod(x, x, n)
            if x == n - 1:
                witness = False
                break
        if witness:
            return False
    return True

def sieve_primes(n: int) -> List[int]:
    if n < 2:
        return []
    is_prime_arr = [True] * (n + 1)
    is_prime_arr[0] = False
    if n >= 1:
        is_prime_arr[1] = False
    
    for p in range(2, int(math.isqrt(n)) + 1):
        if not is_prime_arr[p]:
            continue
        for q in range(p * p, n + 1, p):
            is_prime_arr[q] = False
    
    return [i for i in range(2, n + 1) if is_prime_arr[i]]

def build_filters(limit: int):
    kGood = [1, 3, 7, 9, 13, 27]
    primes = sieve_primes(limit)
    filters = []
    
    for p in primes:
        if p == 2 or p == 5:
            continue
        
        step = 10 % p
        bad = [False] * p
        
        for r in range(p):
            r2 = (r * r) % p
            is_bad = False
            for k in kGood:
                if (r2 + (k % p)) % p == 0:
                    is_bad = True
                    break
            bad[r] = is_bad
        
        filters.append({'p': p, 'step': step, 'bad': bad})
    
    return filters

def first_with_mod(start: int, mod: int, residue: int) -> int:
    current = start % mod
    if current <= residue:
        return start + (residue - current)
    return start + (mod - (current - residue))

def run_sequence(start: int, end: int, residue: int, filters, sum_ref):
    kGood = [1, 3, 7, 9, 13, 27]
    kBad = [5, 11, 15, 17, 19, 21, 23, 25]
    
    filter_threshold = 1000
    step = 70
    
    n = first_with_mod(start, step, residue)
    if n >= end:
        return
    
    square = n * n
    
    residues = [int(n % f['p']) for f in filters]
    step70 = [(f['step'] * 7) % f['p'] for f in filters]
    
    while n < end:
        if square % 3 != 0 and square % 13 != 0:
            ok = True
            if n >= filter_threshold:
                for i, f in enumerate(filters):
                    if f['bad'][residues[i]]:
                        ok = False
                        break
            
            if ok:
                for k in kGood:
                    if not is_prime(square + k):
                        ok = False
                        break
                
                if ok:
                    for k in kBad:
                        if is_prime(square + k):
                            ok = False
                            break
                
                if ok:
                    sum_ref[0] += n
        
        n += step
        square += 140 * (n - step) + 4900
        
        for i, f in enumerate(filters):
            p = f['p']
            next_val = residues[i] + step70[i]
            if next_val >= p:
                next_val -= p
            residues[i] = next_val

def worker(t, block, limit, filters, tasks_sum_arr):
    start = 10
    count = (limit - start + 9) // 10
    idx_start = t * block
    idx_end = min(count, idx_start + block)
    n_start = start + idx_start * 10
    n_end = start + idx_end * 10
    
    sum_ref = [0]
    run_sequence(n_start, n_end, 10, filters, sum_ref)
    run_sequence(n_start, n_end, 60, filters, sum_ref)
    tasks_sum_arr[t] = sum_ref[0]

def solve(limit: int, filters, allow_multithreading: bool, requested_threads: int) -> int:
    if not allow_multithreading:
        threads = 1
    elif requested_threads > 0:
        threads = requested_threads
    else:
        try:
            import os
            threads = os.cpu_count() or 1
        except:
            threads = 1
    
    if threads <= 1:
        sum_ref = [0]
        run_sequence(10, limit, 10, filters, sum_ref)
        run_sequence(10, limit, 60, filters, sum_ref)
        return sum_ref[0]
    
    start = 10
    count = (limit - start + 9) // 10
    block = (count + threads - 1) // threads
    
    import multiprocessing
    
    tasks_sum = multiprocessing.Array('L', threads)
    processes = []
    
    
    
    for t in range(threads):
        p = multiprocessing.Process(target=worker, args=(t, block, limit, filters, tasks_sum))
        processes.append(p)
        p.start()
    
    for p in processes:
        p.join()
    
    return sum(tasks_sum)

def run_checkpoints() -> bool:
    if not is_prime(2) or not is_prime(3) or is_prime(1) or is_prime(21):
        return False
    
    filters = build_filters(2000)
    result = solve(1000000, filters, False, 1)
    return result == 1242490

def main():
    import argparse
    
    parser = argparse.ArgumentParser()
    parser.add_argument('--skip-checkpoints', action='store_true')
    parser.add_argument('--single-thread', action='store_true')
    parser.add_argument('--limit', type=int, default=150000000)
    parser.add_argument('--threads', type=int, default=0)
    
    args = parser.parse_args()
    
    if not args.skip_checkpoints:
        if not run_checkpoints():
            sys.exit(2)
    
    filters = build_filters(2000)
    result = solve(args.limit, filters, not args.single_thread, args.threads)
    print(result)

if __name__ == '__main__':
    main()

Java

import java.math.BigInteger;
import java.util.ArrayList;
import java.util.List;
import java.util.stream.IntStream;

public class Euler146 {
    static class ModFilter {
        int p;
        int step;
        boolean[] bad;
    }

    static boolean isPrime(long n) {
        if (n < 2) return false;
        long[] small = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
        for (long p : small) {
            if (n == p) return true;
            if (n % p == 0) return false;
        }
        return BigInteger.valueOf(n).isProbablePrime(20);
    }

    static List<Integer> sievePrimes(int limit) {
        boolean[] isPrime = new boolean[limit + 1];
        for (int i = 2; i <= limit; i++) isPrime[i] = true;
        for (int p = 2; p * p <= limit; p++) {
            if (isPrime[p]) {
                for (int q = p * p; q <= limit; q += p) {
                    isPrime[q] = false;
                }
            }
        }
        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= limit; i++) {
            if (isPrime[i]) primes.add(i);
        }
        return primes;
    }

    static List<ModFilter> buildFilters(int limit) {
        long[] kGood = {1, 3, 7, 9, 13, 27};
        List<Integer> primes = sievePrimes(limit);
        List<ModFilter> filters = new ArrayList<>();
        for (int p : primes) {
            if (p == 2 || p == 5) continue;
            ModFilter f = new ModFilter();
            f.p = p;
            f.step = 10 % p;
            f.bad = new boolean[p];
            for (int r = 0; r < p; r++) {
                int r2 = (int) ((1L * r * r) % p);
                boolean bad = false;
                for (long k : kGood) {
                    if ((r2 + (k % p)) % p == 0) {
                        bad = true;
                        break;
                    }
                }
                f.bad[r] = bad;
            }
            filters.add(f);
        }
        return filters;
    }

    static long runSequence(long start, long end, long residue, List<ModFilter> filters) {
        long[] kGood = {1, 3, 7, 9, 13, 27};
        long[] kBad = {5, 11, 15, 17, 19, 21, 23, 25};
        long sum = 0;
        
        long step = 70;
        long currentMod = start % step;
        long n;
        if (currentMod <= residue) {
            n = start + (residue - currentMod);
        } else {
            n = start + (step - (currentMod - residue));
        }

        if (n >= end) return 0;

        long square = n * n;
        
        int[] residues = new int[filters.size()];
        int[] step70 = new int[filters.size()];
        for (int i = 0; i < filters.size(); i++) {
            ModFilter f = filters.get(i);
            residues[i] = (int) (n % f.p);
            step70[i] = (f.step * 7) % f.p;
        }

        while (n < end) {
            if (square % 3 != 0 && square % 13 != 0) {
                boolean ok = true;
                if (n >= 1000) {
                    for (int i = 0; i < filters.size(); i++) {
                        if (filters.get(i).bad[residues[i]]) {
                            ok = false;
                            break;
                        }
                    }
                }
                if (ok) {
                    for (long k : kGood) {
                        if (!isPrime(square + k)) {
                            ok = false;
                            break;
                        }
                    }
                    if (ok) {
                        for (long k : kBad) {
                            if (isPrime(square + k)) {
                                ok = false;
                                break;
                            }
                        }
                    }
                    if (ok) {
                        sum += n;
                    }
                }
            }
            n += step;
            square += 140 * (n - step) + 4900;
            for (int i = 0; i < filters.size(); i++) {
                int p = filters.get(i).p;
                int next = residues[i] + step70[i];
                if (next >= p) next -= p;
                residues[i] = next;
            }
        }
        return sum;
    }

    static long solve(long limit, List<ModFilter> filters) {
        int threads = Runtime.getRuntime().availableProcessors();
        if (threads < 1) threads = 1;
        long start = 10;
        long count = (limit - start + 9) / 10;
        long block = (count + threads - 1) / threads;
        
        return IntStream.range(0, threads).parallel().mapToLong(t -> {
            long idxStart = t * block;
            long idxEnd = Math.min(count, idxStart + block);
            long nStart = start + idxStart * 10;
            long nEnd = start + idxEnd * 10;
            long localSum = 0;
            localSum += runSequence(nStart, nEnd, 10, filters);
            localSum += runSequence(nStart, nEnd, 60, filters);
            return localSum;
        }).sum();
    }

    public static void main(String[] args) {
        long limit = 150000000L;
        List<ModFilter> filters = buildFilters(2000);
        System.out.println(solve(limit, filters));
    }
}