Problem 291: Panaitopol Primes

Project Euler Problem 291 Solution

EulerSolve provides an optimized solution for Project Euler Problem 291, Panaitopol Primes, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary We must count the prime values of $$Q_n=2n^2+2n+1=n^2+(n+1)^2$$ under the bound $Q_n \lt \text{LIMIT}=5\cdot 10^{15}$. Testing each value independently would be far too slow, so the solution uses a sieve specialized to this quadratic family. Mathematical Approach 1. Bounding the index range The sequence is strictly increasing for $n\ge 1$, so it is enough to search indices $$1\le n\le m,\qquad m=\max\{n: 2n^2+2n+1\lt \text{LIMIT}\}.$$ Solving the quadratic inequality gives the approximation $$m\approx \left\lfloor \frac{\sqrt{2\cdot \text{LIMIT}-1}-1}{2}\right\rfloor,$$ and the code then adjusts by a couple of integer checks to respect the strict inequality. 2. Why every sequence prime generates two residue classes Assume $p=Q_i$ is prime. Then, modulo $p$, $$Q_n-Q_i=2(n-i)(n+i+1).$$ Since $Q_i\equiv 0 \pmod p$, we obtain $$Q_n\equiv 0\pmod p \iff 2(n-i)(n+i+1)\equiv 0\pmod p.$$ Because $p$ is odd, $2$ is invertible modulo $p$, so the only roots are $$\boxed{n\equiv i\pmod p \quad \text{or}\quad n\equiv -i-1\pmod p.}$$ The second class is the same symmetry as $$Q_{-n-1}=Q_n.$$ 3. Worked example of the residue classes Take $i=1$....

Detailed mathematical approach

Problem Summary

A prime $p$ is called a Panaitopol prime if

$$p=\frac{x^4-y^4}{x^3+y^3}$$

for some positive integers $x$ and $y$. We must count the Panaitopol primes below $5\cdot10^{15}$. The statement gives no examples; for instance, $x=6$ and $y=3$ give $1215/243=5$, so $5$ is a Panaitopol prime. ^[1]

The whole solution in one sentence: cancelling the common factor of $x$ and $y$ shows that a prime can be written in this way exactly when it equals $n^2+(n+1)^2$ for some $n\ge1$, so we count the primes among the $49\,999\,999$ numbers $2n^2+2n+1$ below the bound, with a sieve over the index $n$ that takes its sieving primes from these numbers themselves.

Mathematical Approach

Sections 1–4 find out which primes the fraction can produce, and Sections 5–9 count them with a sieve.

1. Small cases

The fraction is positive only when $x\gt y$. A search over all pairs with $y\lt x\le400$ finds $109$ different integer values, such as $5$, $10$, $13$, $15$ and $20$, but only four of them are prime:

$5$ from $(x,y)=(6,3)=3\cdot(2,1)$;
$13$ from $(21,14)=7\cdot(3,2)$;
$41$ from $(105,84)=21\cdot(5,4)$;
$61$ from $(186,155)=31\cdot(6,5)$.

Two patterns stand out. Each prime is the sum of two consecutive squares: $5=1^2+2^2$, $13=2^2+3^2$, $41=4^2+5^2$ and $61=5^2+6^2$. Each pair is a common factor times two consecutive numbers $(n+1,n)$, and the factors $3$, $7$, $21$ and $31$ equal $n^2+n+1$ for $n=1,2,4,5$. These pairs grow like $n^3$: the next prime, $113=7^2+8^2$, needs $(456,399)$, and primes near $5\cdot10^{15}$ need pairs of about $24$ digits, so instead of searching we prove both patterns.

2. Cancelling the common factor

Let $d=\gcd(x,y)$ and write $x=da$, $y=db$ with coprime integers $a\gt b\ge1$. The difference of two squares, used twice, and the sum of two cubes give the factorisations ^[3]

$$\begin{aligned}x^4-y^4&=(x-y)(x+y)(x^2+y^2),\\x^3+y^3&=(x+y)(x^2-xy+y^2).\end{aligned}$$

The factor $x+y$ cancels. The remaining factors are $d$ or $d^2$ times the same expressions in $a$ and $b$, so

$$\frac{x^4-y^4}{x^3+y^3}=\frac{d\,(a-b)(a^2+b^2)}{a^2-ab+b^2}.$$

For $(6,3)$ we have $d=3$, $a=2$, $b=1$, and the right side is $3\cdot1\cdot5/3=5$.

3. The denominator divides the common factor

Put $D=a^2-ab+b^2$. By Euclid's lemma, a prime that divides a product divides one of the factors. ^[4] A prime $q$ dividing $D$ and $a$ divides $D-a(a-b)=b^2$, hence $b$, which is impossible for coprime $a$ and $b$; likewise no prime divides both $D$ and $b$. Since $a^2+b^2=D+ab$ and $D=(a-b)^2+ab$, a prime dividing $D$ and $a^2+b^2$, or $D$ and $a-b$, divides $ab$, hence $a$ or $b$, which was just excluded. So $D$ is coprime to $(a-b)(a^2+b^2)$.

If the fraction is an integer, $D$ divides $d\,(a-b)(a^2+b^2)$, and the general form of Euclid's lemma (a divisor of a product that is coprime to one factor divides the other) gives $D\mid d$. ^[4] With $d=kD$ for an integer $k\ge1$,

$$\frac{x^4-y^4}{x^3+y^3}=k\,(a-b)(a^2+b^2).$$

Conversely, every such product with coprime $a\gt b\ge1$ and $k\ge1$ occurs, for $x=kDa$ and $y=kDb$. For example, $10=2\cdot1\cdot5$ comes from $(12,6)$, and $20=1\cdot2\cdot10$ from $(21,7)$.

4. The Panaitopol primes

A product $k\,(a-b)(a^2+b^2)$ of positive integers is a prime only if all factors but one equal $1$. Since $a\ge2$, the factor $a^2+b^2$ is at least $5$, so $k=1$ and $a=b+1$. With $b=n$ the prime is

$$p=n^2+(n+1)^2=2n^2+2n+1,$$

and $d=D=n^2+n+1$. Conversely, for every $n\ge1$ the pair $x=(n+1)(n^2+n+1)$, $y=n(n^2+n+1)$ has $k=1$ and $a-b=1$, so its fraction is $n^2+(n+1)^2$. For $n=2$ this is the pair $(21,14)$:

$$\frac{21^4-14^4}{21^3+14^3}=\frac{156065}{12005}=13.$$

Hence the Panaitopol primes are exactly the primes among the numbers $Q_n=2n^2+2n+1$ with $n\ge1$, and each of them comes from exactly one pair $(x,y)$. They form the sequence A027862 of the OEIS, which begins $5,13,41,61,113,181,313,421,613,761$. ^[2] Whether there are infinitely many is unknown: the Bunyakovsky conjecture predicts it, but that conjecture has not been proved for any polynomial of degree $2$ or more. ^[9]

5. Which indices to test

Doubling $Q_n$ and subtracting $1$ gives a square, $2Q_n-1=(2n+1)^2$, so $Q_n\lt L$ exactly when $(2n+1)^2\lt2L-1$. ^[2] For $L=5\cdot10^{15}$ we get $2L-1=10^{16}-1$, and the largest odd number with a smaller square is $10^8-1$, so the indices are $1\le n\le m$ with $m=49\,999\,999$. For $L=1000$ the same rule gives $m=21$, and for $L=10^6$ it gives $m=706$.

Trial division of these $5\cdot10^7$ numbers, with divisors up to $\sqrt{5\cdot10^{15}}\approx7\cdot10^7$, or a sieve of Eratosthenes over all numbers up to $5\cdot10^{15}$ would be far too slow. ^[6] Instead, we sieve over the indices $n$.

6. Where a prime factor comes back

For any two indices $n$ and $i$,

$$Q_n-Q_i=2(n-i)(n+i+1).$$

Every $Q_n=2n(n+1)+1$ is odd, so a prime $p$ dividing $Q_i$ is odd. By Euclid's lemma, $p$ then divides $Q_n$ exactly when it divides $n-i$ or $n+i+1$, that is, when $n\equiv i\pmod p$ or $n\equiv-i-1\pmod p$; here $u\equiv v\pmod p$ means that $p$ divides $u-v$. The multiples of $p$ among the $Q_n$ therefore sit at the indices of two arithmetic progressions with step $p$. For $p=5=Q_1$ these are $n\equiv1$ and $n\equiv3\pmod5$: $Q_3=25$, $Q_6=85$, $Q_8=145$, $Q_{11}=265$ and $Q_{13}=365$ are multiples of $5$, while $Q_2=13$, $Q_4=41$, $Q_5=61$, $Q_7=113$, $Q_9=181$ and $Q_{10}=221$ are not.

The two progressions never coincide: they would only if $p$ divided $2i+1$, but $2Q_i-(2i+1)^2=1$, so $Q_i$ and $2i+1$ have no common prime factor. The same identity gives $(2i+1)^2\equiv-1\pmod p$, and by the first supplement to quadratic reciprocity $-1$ is a square modulo an odd prime $p$ exactly when $p\equiv1\pmod4$. ^[5] Conversely, if $p\equiv1\pmod4$, some $t$ has $t^2\equiv-1\pmod p$, and an index $n$ with $2n+1\equiv t\pmod p$ gives $p\mid2Q_n$, hence $p\mid Q_n$. So the primes that divide some $Q_n$ are exactly the primes $5,13,17,29,37,41,53,\dots$ of the form $4k+1$. ^[2]

7. The first index of a prime

For a prime $p$ that divides some $Q_n$, let $s(p)$ be the smallest index $n\ge1$ with $p\mid Q_n$; for example $s(5)=1$, $s(13)=2$, $s(17)=6$, $s(29)=8$ and $s(37)=15$. The key estimate is

$$p\ge2\,s(p)+3.$$

By Section 6, applied with $i=s(p)$, the indices $n\ge0$ with $p\mid Q_n$ are those whose remainder on division by $p$ is $r$ or $p-1-r$, where $r$ is the remainder of $s(p)$. Neither remainder is $0$, because $Q_0=1$. So both lie between $1$ and $p-2$, the smallest positive index of the two progressions is the smaller remainder, and $s(p)=\min(r,p-1-r)\le(p-1)/2$. Equality $p=2s(p)+1$ is impossible by Section 6, and $p$ is odd, so $p\ge2s(p)+3$. The examples agree: $5\ge5$, $13\ge7$, $17\ge15$, $29\ge19$ and $37\ge33$.

Call a prime factor $p$ of $Q_i$ new if $s(p)=i$, and old if $s(p)\lt i$. The estimate has a strong consequence: $Q_i$ has at most one new prime factor, and it divides $Q_i$ only once. Indeed, two new prime factors $p$ and $q$, or one new prime with $p^2\mid Q_i$, would give a divisor $pq$ of $Q_i$, with $q=p$ in the second case, and

$$pq\ge(2i+3)^2\gt2i^2+2i+1=Q_i,$$

which is impossible. For example, $Q_6=85=5\cdot17$ has the old factor $5$ and the new factor $17$, and $Q_{10}=221=13\cdot17$ has two old factors and no new one.

8. The sieve and why it works

Keep a number $V_n$ and a mark for every index $1\le n\le m$, starting with $V_n=Q_n$ and no marks, and go through $i=1,2,\dots,m$ in increasing order:

if $i$ is not marked, count $Q_i$ as a prime;
if $V_i\gt1$, let $p=V_i$; for every index $j$ with $i\lt j\le m$ in the two progressions $j\equiv i$ and $j\equiv-i-1\pmod p$, mark $j$ and divide $V_j$ by $p$ as often as possible.

Unlike in the sieve of Eratosthenes, ^[6] each prime crosses out two progressions, and the sieving primes are not known in advance: they are the leftover values $V_i$, and even a marked index can deliver one, as index $6$ delivers $17$.

Claim. When the loop reaches $i$, $V_i$ is $Q_i$ with every old prime factor divided out completely, so by Section 7 it is $1$ or the new prime factor of $Q_i$; and $i$ is marked exactly when $Q_i$ has an old prime factor.

The proof is by induction on $i$. Let $q$ be an old prime factor of $Q_i$, and $s=s(q)\lt i$. By the claim for $s$, when the loop was at $s$, the value $V_s$ was the new prime factor of $Q_s$, which is $q$. The loop then visited every index $j\le m$ above $s$ in both progressions modulo $q$: the first starts at $s+q$, the second at its smallest positive member $q-s-1$, which is at least $s+2$ because $q\ge2s+3$. Index $i$ lies in one of them by Section 6, so it was marked and lost every factor $q$. Conversely, the loop at an index $s$ only marks indices $j$ for which the prime $V_s$ divides $Q_j$, and divides only by that prime. So exactly the old prime factors are gone from $V_i$, and $i$ is marked exactly when there is one.

It follows that the sieve counts exactly the primes. If $Q_i=p$ is a prime, then $s(p)=i$, because $0\lt Q_s\lt p$ for every $s\lt i$; so $Q_i$ has no old factor, and $i$ is counted. If $i$ is not marked, every prime factor of $Q_i$ is new, so by Section 7 the number $Q_i\ge5$ is a single prime. Dividing as often as possible is needed for the claim: $Q_{119}=28561=13^4$, and with only one division per visit the loop would reach index $119$ with $V_{119}=169=13^2$, which is not a prime.

9. Worked example: the bound $1000$

For $L=1000$ we have $m=21$, and the loop does the following.

$i=1$: $V_1=5$ is counted. Its progressions are $6,11,16,21$ and $3,8,13,18$; all are marked and divided by $5$, so $V_3=25$ becomes $1$ and $V_{21}=925=5^2\cdot37$ becomes $37$.
$i=2$: $V_2=13$ is counted and marks $15$ and $10\equiv-3\pmod{13}$, so $V_{15}=481$ becomes $37$ and $V_{10}=221$ becomes $17$.
$i=6$: marked, with the new prime $V_6=85/5=17$; its progressions start at $23$ and at $17-7=10$, so $V_{10}$ becomes $1$.
$i=8$: marked, with $V_8=145/5=29$, which reaches $29-9=20$: $V_{20}=841=29^2$ is divided twice and becomes $1$.
$i=15$: marked, with $V_{15}=37$, which reaches $37-16=21$, so $V_{21}$ becomes $1$.
At $i=3,10,20,21$ the value is $1$. The indices not yet listed hold new primes, from $41$ at $i=4$ to $761$ at $i=19$, and none of them reaches an index up to $21$.

The unmarked indices are $1,2,4,5,7,9,12,14,17,19$, so the count is $10$: the primes $5,13,41,61,113,181,313,421,613,761$, the first ten terms of A027862. ^[2]

How the Code Works

The C++ program ^[10] follows Sections 5–8. parse_arguments accepts --limit=L, which replaces the default bound $5\cdot10^{15}$ and must be at least $5$, and --skip-checkpoints; the program reads the bound and later doubles it without overflow checks, which matters only above $2^{63}$. q_value(n) returns $2n^2+2n+1$. isqrt_u64 takes a floating-point square root in long double and corrects it with exact 128-bit products until it equals $\lfloor\sqrt x\rfloor$. max_index_for_limit starts from the estimate $m=(\lfloor\sqrt{2L-1}\rfloor-1)/2$ of Section 5, with integer division, and corrects it with two loops until $Q_m\lt L\le Q_{m+1}$. The estimate is exact unless $L$ is itself some $Q_n$, and then it is one too large: for $L=13$ it gives $2$, and the first loop steps down to $1$.

solve keeps the numbers $V_n$ in the vector values of 64-bit integers and the marks in composite, one byte per index. Its loop over i is the sieve of Section 8: it increments count when composite[i] is $0$, skips the index when values[i] is $1$, and otherwise sets prime = values[i], the $p$ of Section 8. The first inner loop runs over $j=i+p,i+2p,\dots$ up to $m$. The second starts at prime - mod with mod = (i + 1) % prime, the smallest positive $j\equiv-i-1\pmod p$. The code also handles mod == 0 and a start equal to i, but these guards never fire: $p\ge2i+3$ makes mod equal to $i+1$ and the start equal to $p-i-1\ge i+2$. Both inner loops set composite[j] and divide values[j] by prime while it is divisible. All numbers fit in 64 bits: $Q_m\lt5\cdot10^{15}\lt2^{63}$, and an index never exceeds $m+p$.

brute_small(L) tests every $Q_n\lt L$ by trial division. run_checkpoints requires solve and brute_small to agree for $L=1000$ and $L=10^6$, where both give $10$ and $175$; then main prints solve for the chosen bound. The Python version nests the same helpers inside solve and uses math.isqrt, a list of Python integers for values and a bytearray for the marks; its bound is fixed and it runs no checkpoints. The Java method solveLimit keeps $m$ and the outer index in int variables, enough for $49\,999\,999$, uses long[] and byte[] arrays, and casts its long inner index to int; it also skips the checkpoints.

Complexity and Verification

The outer loop makes $m$ steps. A prime $p$ found at index $s$ visits at most $(m-s)/p$ indices in its first progression and $(m+s+1)/p$ in its second, together at most $(2m+1)/p$. Every prime is found only once, at its first index, and a prime $p\gt2m$ visits nothing, because its progressions start at $s+p\gt m$ and at $p-s-1\ge(p+1)/2\gt m$ by Section 7. The sum of $1/p$ over the primes up to $x$ grows like $\ln\ln x$, ^[7] so the sieve makes $O(m\log\log m)$ visits, as the sieve of Eratosthenes does. ^[6]

A visit to $j$ divides by $p$ once for each factor $p$ of $Q_j$, so its $k$-th division needs $p^k\mid Q_j$. If $p^k\mid Q_r$, then $p^k\mid Q_j$ exactly when $p^k$ divides $j-r$ or $j+r+1$, since by Section 6 $p$ cannot divide both; so these indices form two progressions with step $p^k$. Summed over the prime powers $p^k\le Q_m$ with $k\ge2$, this allows at most $\sum(2m/p^k+2)=O(m)$ extra divisions, as $\sum_p\sum_{k\ge2}p^{-k}\lt1$ and there are $O(\sqrt{Q_m})=O(m)$ such prime powers. With $m\approx\sqrt{L/2}$, the time is $O(\sqrt L\log\log L)$ arithmetic operations and the memory $O(\sqrt L)$.

For the target bound only the primes $p\equiv1\pmod4$ below $2m$ can visit anything, and the sum of their reciprocals is about $1.17$, so there are at most about $1.17\cdot10^8$ visits; an instrumented copy of the program counted $114\,119\,423$ visits and $120\,581\,657$ successful divisions. On a desktop computer the C++ program, compiled with -O2, takes about $0.9$ s with a peak memory of about $450$ MB, the $9$ bytes per index of its two arrays; the Java program takes about $1.3$ s and $500$ MB, and the Python program about $42$ s and $4.2$ GB, because its list holds references to integer objects.

Besides the brute-force comparison in run_checkpoints, we checked, with code that is not part of the published solutions:

the search of Section 1, the coprimality of Section 3 for all coprime $a\gt b$ with $a\lt300$, and the pair of Section 4 for $n\lt30$;
for every prime below $20000$, that it divides some $Q_n$ exactly when it is of the form $4k+1$, and that $p\ge2s(p)+3$;
at the full target bound, with an instrumented copy of the sieve, that every leftover $V_i$ is $1$ or a prime, that every sieving prime satisfies $p\ge2i+3$, that the two guards never fire, and that the mark of each of the $49\,999\,999$ indices agrees with a deterministic Miller–Rabin test of $Q_i$; ^[8]
the counts $10$, $175$, $3349$ and $76334$ below $10^3$, $10^6$, $10^9$ and $10^{12}$, each confirmed by a direct primality test of every $Q_n$;
that the C++, Java and Python programs print the same count.

Footnotes and References

The references below supply the background facts used above. The characterisation of the Panaitopol primes, the first-index estimate, the correctness proof of the sieve and its operation counts are derived explicitly in this article.

Project Euler 291 — Panaitopol Primes. The official statement calls a prime p a Panaitopol prime if p = (x⁴ − y⁴)/(x³ + y³) for some positive integers x and y, and asks how many Panaitopol primes are less than 5 × 10¹⁵.
OEIS A027862 — Primes of the form j² + (j+1)². The entry begins 5, 13, 41, 61, 113, 181, 313, 421, 613, 761, also describes these primes as the primes p for which 2p − 1 is a square, and lists them as a subsequence of the primes of the form 4k + 1.
Wikipedia — Factorization. The section on recognisable patterns gives the difference of two squares E² − F² = (E + F)(E − F) and the sum of two cubes E³ + F³ = (E + F)(E² − EF + F²).
Wikipedia — Euclid's lemma. If a prime divides the product of two integers, it divides at least one of them; more generally, if an integer n divides ab and is coprime to a, then n divides b.
Wikipedia — Quadratic reciprocity. The first supplement states that x² ≡ −1 (mod p) is solvable for an odd prime p if and only if p ≡ 1 (mod 4).
Wikipedia — Sieve of Eratosthenes. The sieve generates the multiples of each prime as an arithmetic progression with that prime as its step; it needs O(n log log n) operations for the primes below n, because the sum of the reciprocals of the primes grows like log log n.
Wikipedia — Mertens' theorems. By Mertens' second theorem, the sum of 1/p over the primes p ≤ n minus ln ln n tends to a constant, the Meissel–Mertens constant.
Wikipedia — Miller–Rabin primality test. For n < 2⁶⁴ it is enough to test the bases 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31 and 37.
Wikipedia — Bunyakovsky conjecture. A polynomial with positive leading coefficient, irreducible over the rationals and with no common factor of all its values, should take infinitely many prime values. Only degree 1, Dirichlet's theorem, is proved; 2n² + 2n + 1 meets the three conditions.
C++ — Euler291.cpp. The linked immutable C++ revision contains max_index_for_limit, solve, brute_small and run_checkpoints. It is the implementation record for this article.

Problem 291 source code

C++

#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <vector>

namespace {

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

struct Options {
    u64 limit = 5000000000000000ULL;
    bool run_checkpoints = 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.rfind("--threads=", 0U) == 0U) {
            const std::string tail = arg.substr(10);
            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');
                if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) {
                    return false;
                }
            }
            options.requested_threads = static_cast<unsigned>(parsed);
            continue;
        }
        if (parse_u64_after_prefix(arg, "--limit=", options.limit)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 5ULL;
}

u64 isqrt_u64(const u64 x) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(x)));
    while ((r + 1ULL) > r && static_cast<u128>(r + 1ULL) * static_cast<u128>(r + 1ULL) <= x) {
        ++r;
    }
    while (static_cast<u128>(r) * static_cast<u128>(r) > x) {
        --r;
    }
    return r;
}

inline u64 q_value(const u64 n) {
    return 2ULL * n * n + 2ULL * n + 1ULL;
}

u64 max_index_for_limit(const u64 limit) {
    if (limit <= 5ULL) {
        return 0ULL;
    }
    u64 m = (isqrt_u64(2ULL * limit - 1ULL) - 1ULL) / 2ULL;
    while (m > 0ULL && q_value(m) >= limit) {
        --m;
    }
    while (q_value(m + 1ULL) < limit) {
        ++m;
    }
    return m;
}

u64 solve(const u64 limit) {
    const u64 m = max_index_for_limit(limit);
    if (m == 0ULL) {
        return 0ULL;
    }

    std::vector<u64> values(static_cast<std::size_t>(m + 1ULL), 0ULL);
    std::vector<std::uint8_t> composite(static_cast<std::size_t>(m + 1ULL), 0U);
    for (u64 i = 1ULL; i <= m; ++i) {
        values[static_cast<std::size_t>(i)] = q_value(i);
    }

    u64 count = 0ULL;
    for (u64 i = 1ULL; i <= m; ++i) {
        if (composite[static_cast<std::size_t>(i)] == 0U) {
            ++count;
        }
        const u64 prime = values[static_cast<std::size_t>(i)];
        if (prime == 1ULL) {
            continue;
        }

        for (u64 j = i + prime; j <= m; j += prime) {
            composite[static_cast<std::size_t>(j)] = 1U;
            u64 v = values[static_cast<std::size_t>(j)];
            while (v % prime == 0ULL) {
                v /= prime;
            }
            values[static_cast<std::size_t>(j)] = v;
        }

        const u64 mod = (i + 1ULL) % prime;
        u64 start = (mod == 0ULL) ? prime : (prime - mod);
        if (start == i) {
            start += prime;
        }
        for (u64 j = start; j <= m; j += prime) {
            composite[static_cast<std::size_t>(j)] = 1U;
            u64 v = values[static_cast<std::size_t>(j)];
            while (v % prime == 0ULL) {
                v /= prime;
            }
            values[static_cast<std::size_t>(j)] = v;
        }
    }

    return count;
}

u64 brute_small(const u64 limit) {
    auto is_prime = [](u64 x) {
        if (x < 2ULL) {
            return false;
        }
        if ((x & 1ULL) == 0ULL) {
            return x == 2ULL;
        }
        for (u64 p = 3ULL; p * p <= x; p += 2ULL) {
            if (x % p == 0ULL) {
                return false;
            }
        }
        return true;
    };

    u64 count = 0ULL;
    for (u64 n = 1ULL;; ++n) {
        const u64 q = q_value(n);
        if (q >= limit) {
            break;
        }
        if (is_prime(q)) {
            ++count;
        }
    }
    return count;
}

bool run_checkpoints() {
    if (solve(1000ULL) != brute_small(1000ULL)) {
        std::cerr << "Checkpoint failed for limit=1000" << '\n';
        return false;
    }
    if (solve(1000000ULL) != brute_small(1000000ULL)) {
        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;
    }
    (void)options.requested_threads;
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << solve(options.limit) << '\n';
    return 0;
}

Python

import math

def solve():
    LIMIT = 5_000_000_000_000_000

    def isqrt(x):
        r = math.isqrt(x)
        return r

    def q_value(n):
        return 2 * n * n + 2 * n + 1

    def max_index_for_limit(limit):
        if limit <= 5:
            return 0
        m = (isqrt(2 * limit - 1) - 1) // 2
        while m > 0 and q_value(m) >= limit:
            m -= 1
        while q_value(m + 1) < limit:
            m += 1
        return m

    m = max_index_for_limit(LIMIT)
    if m == 0:
        return '0'

    values = [0] * (m + 1)
    composite = bytearray(m + 1)
    for i in range(1, m + 1):
        values[i] = q_value(i)

    count = 0
    for i in range(1, m + 1):
        if composite[i] == 0:
            count += 1
        prime = values[i]
        if prime == 1:
            continue

        for j in range(i + prime, m + 1, prime):
            composite[j] = 1
            v = values[j]
            while v % prime == 0:
                v //= prime
            values[j] = v

        mod = (i + 1) % prime
        start = prime if mod == 0 else prime - mod
        if start == i:
            start += prime
        for j in range(start, m + 1, prime):
            composite[j] = 1
            v = values[j]
            while v % prime == 0:
                v //= prime
            values[j] = v

    return str(count)

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

Java

public class Euler291 {

    static long isqrt(long x) {
        if (x == 0)
            return 0;
        long r = (long) Math.sqrt(x);
        while (r * r > x)
            r--;
        while ((r + 1) * (r + 1) <= x)
            r++;
        return r;
    }

    static long qValue(long n) {
        return 2 * n * n + 2 * n + 1;
    }

    static long maxIndexForLimit(long limit) {
        if (limit <= 5)
            return 0;
        long m = (isqrt(2 * limit - 1) - 1) / 2;
        while (m > 0 && qValue(m) >= limit)
            m--;
        while (qValue(m + 1) < limit)
            m++;
        return m;
    }

    static long solveLimit(long limit) {
        int m = (int) maxIndexForLimit(limit);
        if (m == 0)
            return 0;

        long[] values = new long[m + 1];
        for (int i = 1; i <= m; i++) {
            values[i] = qValue(i);
        }

        byte[] composite = new byte[m + 1];
        long count = 0;

        for (int i = 1; i <= m; i++) {
            if (composite[i] == 0) {
                count++;
            }
            long prime = values[i];
            if (prime == 1)
                continue;

            for (long j = i + prime; j <= m; j += prime) {
                int ji = (int) j;
                composite[ji] = 1;
                long v = values[ji];
                while (v % prime == 0)
                    v /= prime;
                values[ji] = v;
            }

            long mod = (i + 1) % prime;
            long start = (mod == 0) ? prime : (prime - mod);
            if (start == i) {
                start += prime;
            }

            for (long j = start; j <= m; j += prime) {
                int ji = (int) j;
                composite[ji] = 1;
                long v = values[ji];
                while (v % prime == 0)
                    v /= prime;
                values[ji] = v;
            }
        }

        return count;
    }

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