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 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.

Congruence restrictions from small primes

Several modular arguments remove almost all integers before any primality test is attempted.

If $n$ were odd, then $n^2+1$ would be even and larger than 2, so $n$ must be even.

Modulo 5, the quadratic residues are $0,1,4$. If $n^2\equiv 1\pmod 5$, then $n^2+9\equiv 0\pmod 5$. If $n^2\equiv 4\pmod 5$, then $n^2+1\equiv 0\pmod 5$. Therefore the only surviving case is $n^2\equiv 0\pmod 5$, which forces $5\mid n$.

Modulo 3, a multiple of 3 would make $n^2+3$, $n^2+9$, and $n^2+27$ divisible by 3, so $3\nmid n$.

The strongest restriction comes from modulo 7. The quadratic residues mod 7 are $0,1,2,4$. If $n^2\equiv 0\pmod 7$, then $n^2+7$ is divisible by 7. If $n^2\equiv 1\pmod 7$, then $n^2+13\equiv 0\pmod 7$. If $n^2\equiv 4\pmod 7$, then $n^2+3\equiv 0\pmod 7$. Only $n^2\equiv 2\pmod 7$ survives, which is equivalent to $n\equiv \pm 3\pmod 7$.

Combining $2\mid n$, $5\mid n$, and $n\equiv \pm 3\pmod 7$ leaves exactly two arithmetic progressions: $$n\equiv 10 \pmod{70}\qquad\text{or}\qquad n\equiv 60 \pmod{70}.$$ The same mod-7 argument automatically makes $n^2+5$ and $n^2+19$ divisible by 7, so two forbidden offsets are eliminated for free.

Worked example: $n=10$

The smallest candidate already exhibits the full pattern. For $n=10$, we have $n^2=100$, so the required values are $$101,\ 103,\ 107,\ 109,\ 113,\ 127,$$ and each of them is prime. The forbidden odd offsets give $$105,\ 111,\ 115,\ 117,\ 119,\ 121,\ 123,\ 125,$$ and each of those is composite. So $10$ is a genuine solution and a concrete illustration of the exact condition being tested.

A residue sieve for many small primes

The fixed mod-70 reduction is only the first layer. For each small prime $p$ up to 2000, the implementations precompute the bad residue set $$R_p=\{r\in\{0,\dots,p-1\}: r^2+g\equiv 0\pmod p\text{ for some }g\in G\}.$$ If $n\bmod p\in R_p$, then at least one required value is divisible by $p$, so the candidate can be discarded immediately.

Primes 2 and 5 are omitted from this table because divisibility by them has already been absorbed into the forced form $n\equiv 0\pmod{10}$. For very small $n$, a congruence $n^2+g\equiv 0\pmod p$ could still mean $n^2+g=p$, which is prime rather than composite. The implementations avoid that edge case by activating these precomputed tables only once $n\ge 1000$, a conservative threshold well above $\sqrt{2000}$.

Exact verification after filtering

After the modular screens, the surviving candidates are rare. Each survivor is then checked exactly: the six values $n^2+g$ must be prime and the eight values $n^2+c$ must be composite. Since $$n^2+27 < (150{,}000{,}000)^2+27 < 2.25\times 10^{16},$$ every tested integer fits comfortably in 64-bit arithmetic, so fast Miller-Rabin based certification is practical for the last stage.

How the Code Works

Scanning only the two valid progressions

The C++, Python, and Java implementations iterate only over the sequences $70m+10$ and $70m+60$. This is the direct payoff of the congruence analysis: instead of checking all integers, the search immediately keeps only about one out of every 35 candidates.

Maintaining squares and modular states incrementally

Inside one progression, the step size is always 70, so the square is updated by the recurrence $$ (n+70)^2=n^2+140n+4900.$$ That avoids recomputing $n^2$ from scratch. The same incremental idea is used for every small-prime filter: if the current residue modulo $p$ is known, the next one is obtained by adding $70\bmod p$. The sieve phase therefore becomes a sequence of cheap table lookups instead of repeated modular squaring.

The decision pipeline and parallel search

Each candidate passes through a short pipeline: inexpensive divisibility checks for 3 and 13, lookup in the precomputed bad-residue tables, exact primality tests on the six required offsets, and exact compositeness checks on the eight forbidden offsets. When every condition succeeds, the candidate $n$ is added to the running sum.

The interval splits naturally into independent chunks. The C++ and Java implementations parallelize those chunks directly, and the Python implementation can also distribute chunks across worker processes. The mathematics is identical in all three languages; only the execution strategy changes.

Complexity Analysis

Let $L$ be the search limit and $B$ the upper bound for the small-prime residue tables. Building the tables costs $O(B\log\log B+\sum_{p\le B} p)$ time and $O(\sum_{p\le B} p)$ memory, because each prime $p$ stores a boolean array of length $p$.

The main scan visits only about $L/35$ integers because it uses two residue classes modulo 70. For each candidate it performs constant-time arithmetic plus one residue update per sieve prime, and only a tiny fraction survive to the primality stage. With the fixed choice $B=2000$ used by the implementations, memory usage stays small and the running time is effectively linear in the search limit.

Footnotes and References

Problem page: Project Euler 146
Prime constellations: Wikipedia - Prime k-tuple
Modular arithmetic: Wikipedia - Modular arithmetic
Quadratic residues: Wikipedia - Quadratic residue
Miller-Rabin primality test: Wikipedia - Miller-Rabin primality test
Sieve of Eratosthenes: Wikipedia - Sieve of Eratosthenes

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