Problem 216: Investigating the Primality of 2n²-1

Project Euler Problem 216 Solution

EulerSolve provides an optimized solution for Project Euler Problem 216, Investigating the Primality of 2n²-1, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For integers $n \ge 2$, define $$t_n=2n^2-1.$$ Problem 216 asks for the number of indices $2 \le n \le 50{,}000{,}000$ for which $t_n$ is prime. A naive solution would evaluate fifty million numbers of size about $5\times10^{15}$ and run a primality test on each of them. The implementations take a more number-theoretic route: instead of proving primality one value at a time, they sieve the indices $n$ for which some smaller prime is known to divide $2n^2-1$. The survivors are exactly the prime values. Mathematical Approach The key observation is that divisibility of $2n^2-1$ by a fixed prime turns into a modular square-root problem, so the sequence can be attacked with a residue-class sieve on $n$. Why sieving only up to $\sqrt{2N^2-1}$ is enough Let $N=50{,}000{,}000$. If $t_n$ is composite, then it has a prime divisor $p$ with $$p \le \sqrt{t_n} \le \sqrt{2N^2-1}.$$ Therefore every composite value $2n^2-1$ is detected by some prime at most $\sqrt{2N^2-1}$. Conversely, if no prime in that range divides $t_n$, then $t_n$ has no proper prime factor and must itself be prime. This is the invariant that makes the sieve correct. Since $2n^2-1$ is always odd, the prime $2$ can never divide it, so only odd primes need to be examined. Which primes can divide $2n^2-1$?...

Detailed mathematical approach

Problem Summary

For every integer $n\gt1$ let $t(n)=2n^2-1$. The first values, for $n=2,3,\dots,9$, are $7$, $17$, $31$, $49$, $71$, $97$, $127$ and $161$; all of them are prime except $49=7\cdot7$ and $161=7\cdot23$, and for $n\le10000$ exactly $2202$ of the numbers $t(n)$ are prime. We must count the integers $n$ with $2\le n\le N$ for which $t(n)$ is prime, where $N=50\,000\,000$. ^[1] The largest of these numbers, $t(N)=5\cdot10^{15}-1$, has sixteen digits.

The whole solution in one sentence: for each prime $p\le\sqrt{2N^2-1}$, the indices $n$ for which $p$ divides $2n^2-1$ form at most two arithmetic progressions with difference $p$, found with a square root modulo $p$; we cross them all out, except an index whose value is $p$ itself, and count the indices that remain.

Mathematical Approach

Testing each of the fifty million numbers separately is slow: trial division alone needs up to about four million divisions per prime value. Instead we ask, for a fixed prime $p$, which indices $n$ make $p$ a divisor of $t(n)$. The answer repeats with period $p$, and crossing out these patterns is a sieve in the spirit of the sieve of Eratosthenes. ^[2]

1. The first values

For $n=2,\dots,20$, twelve of the numbers $t(n)$ are prime, namely $7, 17, 31, 71, 97, 127,$ $199, 241, 337, 449, 577, 647$, and seven are composite:

$t(5)=49=7\cdot7$, $t(9)=161=7\cdot23$, $t(12)=287=7\cdot41$;
$t(14)=391=17\cdot23$, $t(16)=511=7\cdot73$;
$t(19)=721=7\cdot103$, $t(20)=799=17\cdot47$.

Two patterns stand out. First, the prime factors $7, 17, 23, 41, 47, 73, 103$ all leave the remainder $1$ or $7$ on division by $8$, while $3$, $5$, $11$ and $13$ never appear; the same holds for every prime factor of every $t(n)$ with $n\le1000$. Second, $7$ divides $t(n)$ for $n=2, 5, 9, 12, 16, 19$, that is, exactly when $n$ leaves the remainder $2$ or $5$ on division by $7$. Sections 2 to 4 explain both patterns.

2. Divisibility by $p$ repeats with period $p$

For every prime $p$ and every $n$,

$$t(n+p)=t(n)+p\,(4n+2p),$$

so $t(n+p)$ and $t(n)$ leave the same remainder on division by $p$; in the notation of congruences, $t(n+p)\equiv t(n)\pmod p$. Whether $p$ divides $t(n)$ therefore depends only on the remainder $r$ of $n$ modulo $p$, and we need the remainders $r\in\{0,1,\dots,p-1\}$ with $2r^2\equiv1\pmod p$. For $p=7$ the values $2r^2-1$ for $r=0,1,\dots,6$ are $-1, 1, 7, 17, 31, 49, 71$, and only $7$ and $49$, at $r=2$ and $r=5$, are multiples of $7$: this is the second pattern of Section 1. The prime $2$ never divides $t(n)$, which is odd.

3. Only primes $p\equiv\pm1\pmod 8$ divide $t(n)$

Let $p$ be an odd prime that divides $t(n)$. Then

$$(2n)^2=2\,t(n)+2\equiv2\pmod p,$$

so $2$ is congruent to a square modulo $p$: it is a quadratic residue modulo $p$. ^[3] The squares modulo $5$ are $0, 1, 4$, so $2$ is not a residue modulo $5$; the squares modulo $7$ are $0, 1, 2, 4$, and indeed $3^2=9\equiv2\pmod 7$. Which primes have $2$ as a residue is the content of the second supplement to the law of quadratic reciprocity: ^[3]

Claim. For an odd prime $p$, the number $2$ is a quadratic residue modulo $p$ if and only if $p\equiv\pm1\pmod 8$, that is, $p\equiv1$ or $p\equiv7\pmod 8$.

Proof. By Gauss's lemma, ^[4] $2$ is a residue modulo $p$ exactly when an even number of the multiples $2\cdot1, 2\cdot2, \dots, 2\cdot\frac{p-1}{2}$, reduced modulo $p$, exceed $p/2$. These multiples are the even numbers $2, 4, \dots, p-1$, which are already reduced, and $2k\gt p/2$ means $k\gt p/4$, so their number is $\mu=\frac{p-1}{2}-\lfloor p/4\rfloor$. Writing $p=8j+c$:

$c=1$: $\mu=4j-2j=2j$ is even;
$c=3$: $\mu=(4j+1)-2j=2j+1$ is odd;
$c=5$: $\mu=(4j+2)-(2j+1)=2j+1$ is odd;
$c=7$: $\mu=(4j+3)-(2j+1)=2j+2$ is even. $\square$

For $p=7$ the multiples $4$ and $6$ exceed $3.5$, so $\mu=2$. This is the first pattern of Section 1: no prime $p\equiv\pm3\pmod 8$ ever divides $t(n)$.

4. Two progressions for every prime $p\equiv\pm1\pmod 8$

Conversely, let $p\equiv\pm1\pmod 8$, choose $s$ with $s^2\equiv2\pmod p$, and put $h=(p+1)/2$. Since $2h=p+1\equiv1\pmod p$, the number $h$ is the inverse of $2$ modulo $p$. The remainder $r$ of $sh$ on division by $p$ satisfies

$$2r^2\equiv2\,s^2h^2\equiv(2h)^2\equiv1\pmod p,$$

so $p$ divides $t(r)$; equivalently, $r^2\equiv h\pmod p$. Apart from $-r$ there is no other solution: if also $2x^2\equiv1$, then $p$ divides $2\,(x-r)(x+r)$, and a prime that divides a product divides one of its factors, so $x\equiv\pm r$. The two solutions differ, because $p$ is odd and $r\not\equiv0$. Hence, for every prime $p\equiv\pm1\pmod 8$,

$$p\mid t(n)\iff n\equiv\pm r\pmod p,$$

and for every other prime there is no such $n$ at all. Neither $r$ nor $p-r$ is $0$, $1$ or $p-1$, because $2\cdot0^2-1=-1$ and $2\cdot1^2-1=1$ are not multiples of $p$, so both lie between $2$ and $p-2$. For $p=7$ we may take $s=3$ and $h=4$: the remainder of $12$ is $r=5$, and $p-r=2$, the two remainders found in Section 2.

5. Square roots modulo $p$ when $p\equiv7\pmod 8$

To find $s$ we use Euler's criterion: ^[5] for an odd prime $p$ and an integer $a$ not divisible by $p$,

$$a^{(p-1)/2}\equiv\pm1\pmod p,$$

with $+1$ exactly when $a$ is a quadratic residue. If $p\equiv7\pmod 8$, then $(p+1)/4$ is an integer, and $s=2^{(p+1)/4}$ satisfies

$$s^2=2^{(p+1)/2}=2\cdot2^{(p-1)/2}\equiv2\pmod p,$$

because $2$ is a residue by Section 3. So a single power gives $s$, and then $r\equiv s\,(p+1)/2$. For $p=7$: $s=2^2=4$, $4^2=16\equiv2$, and $r\equiv4\cdot4=16\equiv2$. For $p=23$: $s=2^6=64\equiv18$, $18^2=324=14\cdot23+2$, and $r\equiv18\cdot12=216\equiv9$. So $23$ divides $t(n)$ exactly when $n\equiv9$ or $n\equiv14\pmod{23}$, in agreement with $t(9)=161=7\cdot23$ and $t(14)=391=17\cdot23$.

6. Square roots modulo $p$ when $p\equiv1\pmod 8$

For $p\equiv1\pmod 8$ the exponent $(p+1)/4$ is not an integer, and in general no power of $a$ is a square root of $a$: modulo $17$ the powers of $9$ are $9, 13, 15, 16, 8, 4, 2, 1$, while the square roots of $9$ are $3$ and $14$. The Tonelli–Shanks algorithm ^[6] computes a square root $x$ of any quadratic residue $a$, and the code applies it directly to $a=h=(p+1)/2$, which is a residue because $h\equiv2h^2$ is $2$ times a square. Write $p-1=2^SQ$ with $Q$ odd, here with $S\ge3$, and find a non-residue $z$ by testing $z=2,3,4,\dots$ with Euler's criterion. The algorithm starts from $e=S$, $c=z^Q$, $w=a^Q$ and $x=a^{(Q+1)/2}$, and keeps three congruences modulo $p$:

$$\begin{aligned}x^2&\equiv a\,w,\\ c^{2^{e-1}}&\equiv-1,\\ w^{2^{e-1}}&\equiv1\pmod p.\end{aligned}$$

They hold at the start, because $x^2=a^{Q+1}=a\,w$, $c^{2^{S-1}}=z^{(p-1)/2}\equiv-1$ and $w^{2^{S-1}}=a^{(p-1)/2}\equiv1$ by Euler's criterion. If $w\equiv1$, then $x^2\equiv a$ and we are done. Otherwise let $i$ be the least exponent with $w^{2^i}\equiv1$; then $1\le i\lt e$, and $w^{2^{i-1}}\equiv-1$, because its square is $1$, and, as in Section 4, $y^2\equiv1$ forces $y\equiv\pm1$. Put $b=c^{2^{e-i-1}}$, so that $b^{2^i}=c^{2^{e-1}}\equiv-1$, and replace

$$x\leftarrow xb,\ \ w\leftarrow wb^2,\ \ c\leftarrow b^2,\ \ e\leftarrow i.$$

The first congruence survives because both sides are multiplied by $b^2$, the second because $(b^2)^{2^{i-1}}=b^{2^i}\equiv-1$, and the third because $(wb^2)^{2^{i-1}}\equiv(-1)\cdot(-1)=1$. Every round lowers $e$, and for $e=1$ the third congruence says $w\equiv1$, so after at most $S-1$ rounds $x$ is a square root of $a$.

For $p=17$ we need a square root of $a=9$. Here $p-1=16=2^4\cdot1$, so $S=4$ and $Q=1$. The candidate $z=2$ is a residue, since $2^8\equiv1$, as Section 3 predicts for $p\equiv1\pmod 8$; $z=3$ is not, since $3^8\equiv16\equiv-1$. The start is $e=4$, $c=3$, $w=9$, $x=9$, and three rounds follow:

$9^2\equiv13$, $13^2\equiv16$, $16^2\equiv1$, so $i=3$ and $b=3$: $x=27\equiv10$, $w=81\equiv13$, $c=9$, $e=3$;
$13^2\equiv16$, $16^2\equiv1$, so $i=2$ and $b=9$: $x=90\equiv5$, $w=13\cdot81\equiv16$, $c=81\equiv13$, $e=2$;
$16^2\equiv1$, so $i=1$ and $b=13$: $x=65\equiv14$ and $w=16\cdot169\equiv1$.

Indeed $14^2=196=11\cdot17+9$. So $17$ divides $t(n)$ exactly when $n\equiv14$ or $n\equiv3\pmod{17}$, as $t(3)=17$, $t(14)=391$ and $t(20)=799=17\cdot47$ confirm.

7. The one index that must be spared

Crossing out all $n\equiv\pm r\pmod p$ removes every $n$ with $p\mid t(n)$, and one of them may have $t(n)=p$, which is prime. Since $t$ is increasing, there is at most one such $n$: it exists exactly when $p=2m^2-1$ for an integer $m$, that is, when $(p+1)/2=m^2$ is a perfect square, and then $n=m$. For $p=7=t(2)$ and $p=17=t(3)$ the indices $2$ and $3$ must be spared; for $p=23$ the number $12$ is not a square and nothing is spared.

8. Why the sieve finds exactly the primes

Let $B=\lfloor\sqrt{2N^2-1}\rfloor$. The sieve takes every prime $p\le B$ with $p\equiv\pm1\pmod 8$, crosses out all $n$ with $2\le n\le N$ and $n\equiv\pm r\pmod p$, except the spared index of Section 7, and counts the indices that remain.

Claim. An index $n$ with $2\le n\le N$ is never crossed out if and only if $t(n)$ is prime.

Proof. If $t(n)$ is prime and a prime $p$ crosses $n$ out, then $p$ divides $t(n)$, so $p=t(n)$, and $n$ is the spared index of $p$, a contradiction. If $t(n)$ is composite, write $t(n)=uv$ with $1\lt u\le v$. Then $u^2\le t(n)\le2N^2-1$, so $u\le B$, and every prime factor $p$ of $u$ satisfies $p\le B$ and divides $t(n)$. By Sections 2 and 3, $p$ is odd with $p\equiv\pm1\pmod 8$; by Section 4, $n\equiv\pm r\pmod p$; and $n$ is not spared, because $t(n)\ne p$. So $n$ is crossed out. $\square$

The proof never needs $p\le\sqrt{t(n)}$, so using the same primes for every $n$ is harmless. For $N=50\,000\,000$ the bound is $B=70\,710\,678$, as $\sqrt{5\cdot10^{15}-1}=70\,710\,678.1186\ldots$

9. Worked example: $N=20$

Here $2N^2-1=799$ and $B=28$. Of the primes $2, 3, 5, 7, 11, 13, 17, 19, 23$ up to $28$, only $7$, $17$ and $23$ are $\equiv\pm1\pmod 8$, with the remainders found in Sections 4 to 6:

$p=7$, remainders $2$ and $5$: it crosses out $5, 9, 12, 16, 19$ and spares $n=2$;
$p=17$, remainders $3$ and $14$: it crosses out $14, 20$ and spares $n=3$;
$p=23$, remainders $9$ and $14$: it crosses out $9$ and $14$ a second time.

The crossed-out indices $5, 9, 12, 14, 16, 19, 20$ are exactly the indices of the seven composite values of Section 1, and the $12$ survivors $2, 3, 4, 6, 7, 8, 10, 11, 13, 15, 17, 18$ are the $n\le20$ with $t(n)$ prime. Stopping at $N=8$ leaves $2, 3, 4, 6, 7, 8$, the six primes $7, 17, 31, 71, 97, 127$ of the statement.

10. What is known about these primes

The survivors $2, 3, 4, 6, 7, 8, 10, 11, \dots$ form the sequence A066049 of the OEIS, which is conjectured to be infinite. ^[7] Nobody can prove this. Bunyakovsky's conjecture says that a polynomial with a positive leading coefficient, irreducible over the rationals and with values that have no common prime factor, takes infinitely many prime values; $2x^2-1$ qualifies, since its leading coefficient is positive, it has no rational root and $t(1)=1$, but no case of degree above $1$ has been proved. ^[8] The Bateman–Horn conjecture predicts the count $\frac C2\int_2^N\frac{dx}{\ln x}$ asymptotically, ^[8] where

$$C=\prod_{p}\frac{1-\omega(p)/p}{1-1/p}\approx3.700$$

and $\omega(p)$ is the number of remainders $n$ modulo $p$ with $p\mid t(n)$: by Sections 2 to 4, $\omega(2)=0$, $\omega(p)=2$ for $p\equiv\pm1\pmod 8$ and $\omega(p)=0$ for the other odd primes. Each factor is the ratio of the chance $1-\omega(p)/p$ that $p$ does not divide $t(n)$ to the chance $1-1/p$ for a random integer, for example $2$ for $p=2$ and $5/6$ for $p=7$. The asymptotically equal sum $C\sum_{n=2}^{N}1/\ln t(n)$ gives about $2206$ for $N=10\,000$, close to the true $2202$; but this is an estimate, and the exact count needs a computation such as the sieve.

How the Code Works

The C++ program ^[9] starts in main: parse_arguments accepts --limit-n=N, at least $2$ with the default $50\,000\,000$, and --skip-checkpoints, although its digit parser does not check for overflow of the int range; run_checkpoints requires solve(8) to be $6$ and solve(10000) to be $2202$, the values of the statement, and otherwise ends the program with exit code $2$; finally the result of solve is printed.

solve computes $2N^2-1$ as a 64-bit unsigned integer and $B$ by converting a long double square root to int, without the correction steps of integer_sqrt_int; for the default $N$ this is safe, since the fractional part $0.1186\ldots$ of the root is far from $0$ and $1$. sieve_primes(B) is the sieve of Eratosthenes on a byte array and returns the $4\,157\,407$ primes up to $B$; a second byte array composite records the crossings-out for every $n\le N$. The main loop skips $p=2$ and every prime whose remainder p & 7 is not $1$ or $7$, which leaves $2\,078\,410$ primes. With inv2 $=(p+1)/2$, a prime $p\equiv7$ gets its root as mod_pow(2, (p + 1) / 4, p) times inv2 (Section 5), and a prime $p\equiv1$ from tonelli_shanks(inv2, p) (Section 6), whose variables m, c, t, r are the $e$, $c$, $w$, $x$ of that section. mod_pow uses square-and-multiply: it reads the binary digits of the exponent from the lowest, multiplies the current power of the base, at first the base itself, into the result at every digit $1$, and squares that power after every digit. The branches of tonelli_shanks for $0$, $p=2$, non-residues and $p\equiv3\pmod 4$ are never reached here.

The two remainders are r1, the root, and r2 $=p-{}$r1. integer_sqrt_int corrects a floating-point square root of $(p+1)/2$ until it is exact, and a perfect square gives skip_n, the spared index of Section 7. The lambda mark_progression sets composite[n] for $n=r, r+p, r+2p, \dots$ up to $N$, except skip_n. By Section 4 its adjustment for a start below $2$ is never needed, and the test r2 != r1 always passes. Finally solve counts the $n$ from $2$ to $N$ whose entry is still $0$. The products in mod_pow and tonelli_shanks are formed in unsigned __int128, although $p\lt2^{27}$ keeps them below $2^{54}$.

The Python program follows the same steps inside one function, with the three-argument pow, a bytearray sieve, math.isqrt for $B$ and the square test, and a while loop per progression. The Java program does the same in main with boolean arrays, long arithmetic, (int) Math.sqrt for $B$ and a method mark; its prime array of length maxPrime / 5 is large enough for this $N$, though too short for limits up to $254$. Both fix $N=50\,000\,000$ and have neither options nor checkpoints.

Complexity and Verification

With $B\approx\sqrt2\,N$, the prime sieve takes $O(B\log\log B)$ time and $B+1$ bytes. ^[2] A progression of $p$ has at most $N/p+1$ terms, so the marking takes at most

$$\sum_{p\le B}2\Bigl(\frac Np+1\Bigr)=O(N\log\log N)$$

steps, because the sum of $1/p$ over the primes up to $x$ grows like $\log\log x$; ^[2] for $N=50\,000\,000$ there are $104\,672\,841$ marking steps, about $2.09\,N$. A root costs $O(\log p)$ multiplications for $p\equiv7\pmod 8$; for $p\equiv1\pmod 8$ it costs a few powers of that size, one more for every trial of $z$, and $O(S^2)$ for the rounds. Over all primes up to $B$ the rounds add up to $O(N)$, since $S^2=1+3+\dots+(2S-1)$ and at most $B/2^k+1$ of these primes have $2^k\mid p-1$, and so do the trials, since the least non-residue is bounded on average, a theorem of Erdős. ^[6] In the run for $N=50\,000\,000$, $z$ takes $4.33$ trials on average and never goes beyond $67$, and there are $1.5$ rounds on average and at most $15$. The total time is therefore $O(N\log\log N)$, and the memory $O(N)$.

On the test machine (clang -O2, OpenJDK 17, Python 3.11) the C++, Java and Python programs take about $2.5$, $1.0$ and $12$ seconds, with peak resident memory of 105, 232 and 381 MB. The C++ program trails the Java one because in this build each 128-bit remainder calls the library routine __umodti3; with 64-bit products it also takes about $1.0$ second.

Only the C++ program checks itself, with the two checkpoints. For this article, the published C++ program, run with --limit-n, was compared with an independent count that tests every $t(n)$ with the Miller–Rabin test using the first twelve prime bases, which is deterministic below $2^{64}$. ^[10] The counts agree for every $N$ from $2$ to $3000$, hence index by index for $n\le3000$, and for $N=10^4$, $10^5$, $10^6$, $10^7$ and $5\cdot10^7$; for $N=5\cdot10^7$ the three published programs also print the same number. A copy of the C++ program with an added check confirmed for all $2\,078\,410$ sieving primes that both remainders satisfy $2r^2\equiv1\pmod p$ and lie between $2$ and $p-2$. None of these checks is part of the published solutions.

Footnotes and References

The references below supply the background facts used above. The sieve over residue classes, its correctness proof, the worked examples and all counts are derived or computed in this article.

Project Euler 216 — The Primality of 2n² − 1. The official statement defines t(n) = 2n² − 1 for n > 1, gives its first values and the count 2202 for n ≤ 10000, and asks how many t(n) are prime for n ≤ 50 000 000.
Wikipedia — Sieve of Eratosthenes. Multiples are crossed out from the square of each prime on, primes up to the square root of the limit suffice, and the cost O(n log log n) comes from the sum of the reciprocals of the primes growing like log log n.
Wikipedia — Quadratic reciprocity. A quadratic residue modulo p is a number congruent to a square; by the second supplement, x² ≡ 2 (mod p) is solvable exactly when p ≡ ±1 (mod 8).
Wikipedia — Gauss's lemma (number theory). For an odd prime p and a coprime to p, (a/p) = (−1)ⁿ, where n counts the least positive residues of a, 2a, …, ((p − 1)/2)a above p/2; the page notes that it gives probably the simplest proofs of the second supplementary law.
Wikipedia — Euler's criterion. For an odd prime p and a coprime to p, a^((p − 1)/2) ≡ 1 (mod p) if a is a quadratic residue and ≡ −1 (mod p) otherwise.
Wikipedia — Tonelli–Shanks algorithm. The page gives the algorithm, with M, c, t, R for the e, c, w, x of Section 6, its proof by the same three invariants, and the roots ±n^((p + 1)/4) for p ≡ 3 (mod 4). OEIS A053760 — Smallest positive quadratic nonresidue modulo p gives the mean of the least non-residue over the primes as Σ prime(k)/2^k ≈ 3.67 (Erdős, 1961).
OEIS A066049 — Numbers k such that 2k² − 1 is a prime. The sequence begins 2, 3, 4, 6, 7, 8, 10, 11, 13, 15, 17, 18, the survivors of Section 9, and is conjectured to be infinite.
Wikipedia — Bunyakovsky conjecture and Wikipedia — Bateman–Horn conjecture. The first states the three conditions and notes that no case of degree greater than 1 is proved; the second predicts the count (C/D) ∫₂ˣ dt / log t for a polynomial of degree D, with C the product over p of (1 − N(p)/p)/(1 − 1/p), N(p) being the ω(p) of Section 10.
C++ — Euler216.cpp. The linked immutable revision contains the functions and the two checkpoints described above. It is the implementation record for this article.
Wikipedia — Miller–Rabin primality test. The twelve prime bases from 2 to 37 decide primality correctly for every n < 2⁶⁴, which covers all t(n) < 5 · 10¹⁵.

Problem 216 source code

C++

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

namespace {

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

struct Options {
    int limit_n = 50000000;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--limit-n=", options.limit_n)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit_n >= 2;
}

u64 mod_pow(u64 base, u64 exp, u64 mod) {
    u64 result = 1 % mod;
    u64 cur = base % mod;
    while (exp > 0) {
        if ((exp & 1ULL) != 0ULL) {
            result = static_cast<u64>((static_cast<u128>(result) * cur) % mod);
        }
        cur = static_cast<u64>((static_cast<u128>(cur) * cur) % mod);
        exp >>= 1U;
    }
    return result;
}

u64 tonelli_shanks(u64 n, u64 p) {
    if (n == 0) {
        return 0;
    }
    if (p == 2) {
        return n;
    }
    if (mod_pow(n, (p - 1) / 2, p) != 1) {
        return 0;
    }
    if (p % 4 == 3) {
        return mod_pow(n, (p + 1) / 4, p);
    }

    u64 q = p - 1;
    int s = 0;
    while ((q & 1ULL) == 0ULL) {
        q >>= 1U;
        ++s;
    }

    u64 z = 2;
    while (mod_pow(z, (p - 1) / 2, p) != p - 1) {
        ++z;
    }

    u64 m = static_cast<u64>(s);
    u64 c = mod_pow(z, q, p);
    u64 t = mod_pow(n, q, p);
    u64 r = mod_pow(n, (q + 1) / 2, p);

    while (t != 1) {
        u64 i = 1;
        u64 t2 = static_cast<u64>((static_cast<u128>(t) * t) % p);
        while (t2 != 1) {
            t2 = static_cast<u64>((static_cast<u128>(t2) * t2) % p);
            ++i;
        }
        const u64 b = mod_pow(c, 1ULL << (m - i - 1), p);
        r = static_cast<u64>((static_cast<u128>(r) * b) % p);
        t = static_cast<u64>((static_cast<u128>(t) * b % p * b) % p);
        c = static_cast<u64>((static_cast<u128>(b) * b) % p);
        m = i;
    }
    return r;
}

std::vector<int> sieve_primes(const int limit) {
    std::vector<std::uint8_t> is_prime(static_cast<std::size_t>(limit + 1), 1U);
    is_prime[0] = 0U;
    is_prime[1] = 0U;
    for (int p = 2; static_cast<u64>(p) * p <= static_cast<u64>(limit); ++p) {
        if (is_prime[static_cast<std::size_t>(p)] == 0U) {
            continue;
        }
        for (int q = p * p; q <= limit; q += p) {
            is_prime[static_cast<std::size_t>(q)] = 0U;
        }
    }
    std::vector<int> primes;
    primes.reserve(static_cast<std::size_t>(limit / std::log(static_cast<double>(limit))));
    for (int p = 2; p <= limit; ++p) {
        if (is_prime[static_cast<std::size_t>(p)] != 0U) {
            primes.push_back(p);
        }
    }
    return primes;
}

int integer_sqrt_int(const int x) {
    int r = static_cast<int>(std::sqrt(static_cast<double>(x)));
    while (static_cast<long long>(r + 1) * (r + 1) <= x) {
        ++r;
    }
    while (static_cast<long long>(r) * r > x) {
        --r;
    }
    return r;
}

int solve(const int limit_n) {
    const u64 max_t = 2ULL * static_cast<u64>(limit_n) * static_cast<u64>(limit_n) - 1ULL;
    const int max_prime = static_cast<int>(std::sqrt(static_cast<long double>(max_t)));
    const std::vector<int> primes = sieve_primes(max_prime);

    std::vector<std::uint8_t> composite(static_cast<std::size_t>(limit_n + 1), 0U);

    for (int p : primes) {
        if (p == 2) {
            continue;
        }
        const int mod8 = p & 7;
        if (mod8 != 1 && mod8 != 7) {
            continue;
        }

        const u64 inv2 = static_cast<u64>((p + 1) / 2);
        u64 root = 0;
        if (mod8 == 7) {
            const u64 sqrt2 = mod_pow(2ULL, static_cast<u64>((p + 1) / 4), static_cast<u64>(p));
            root = static_cast<u64>((static_cast<u128>(sqrt2) * inv2) % static_cast<u64>(p));
        } else {
            root = tonelli_shanks(inv2, static_cast<u64>(p));
            if (root == 0) {
                continue;
            }
        }
        const int r1 = static_cast<int>(root);
        const int r2 = (r1 == 0) ? 0 : (p - r1);

        int skip_n = -1;
        const int half = (p + 1) / 2;
        const int sr = integer_sqrt_int(half);
        if (sr * sr == half) {
            skip_n = sr;
        }

        const auto mark_progression = [&](const int residue) {
            int n = residue;
            if (n < 2) {
                const int delta = 2 - n;
                n += ((delta + p - 1) / p) * p;
            }
            for (; n <= limit_n; n += p) {
                if (n == skip_n) {
                    continue;
                }
                composite[static_cast<std::size_t>(n)] = 1U;
            }
        };

        mark_progression(r1);
        if (r2 != r1) {
            mark_progression(r2);
        }
    }

    int count = 0;
    for (int n = 2; n <= limit_n; ++n) {
        if (composite[static_cast<std::size_t>(n)] == 0U) {
            ++count;
        }
    }
    return count;
}

bool run_checkpoints() {
    if (solve(8) != 6) {
        std::cerr << "Checkpoint failed for n <= 8 prefix sample" << '\n';
        return false;
    }
    if (solve(10000) != 2202) {
        std::cerr << "Checkpoint failed for n <= 10000 sample" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << solve(options.limit_n) << '\n';
    return 0;
}

Python

import math

def solve():
    limit_n = 50000000

    def mod_pow(base, exp, mod):
        return pow(base, exp, mod)

    def tonelli_shanks(n, p):
        if n == 0:
            return 0
        if p == 2:
            return n
        if mod_pow(n, (p - 1) // 2, p) != 1:
            return 0
        if p % 4 == 3:
            return mod_pow(n, (p + 1) // 4, p)
        q = p - 1
        s = 0
        while q % 2 == 0:
            q //= 2
            s += 1
        z = 2
        while mod_pow(z, (p - 1) // 2, p) != p - 1:
            z += 1
        m = s
        c = mod_pow(z, q, p)
        t = mod_pow(n, q, p)
        r = mod_pow(n, (q + 1) // 2, p)
        while t != 1:
            i = 1
            t2 = (t * t) % p
            while t2 != 1:
                t2 = (t2 * t2) % p
                i += 1
            b = mod_pow(c, 1 << (m - i - 1), p)
            r = (r * b) % p
            t = (t * b % p * b) % p
            c = (b * b) % p
            m = i
        return r

    def sieve_primes(lim):
        is_p = bytearray(b'\x01' * (lim + 1))
        is_p[0] = 0
        is_p[1] = 0
        for p in range(2, math.isqrt(lim) + 1):
            if is_p[p]:
                is_p[p*p::p] = bytearray(len(is_p[p*p::p]))
        return [p for p in range(2, lim + 1) if is_p[p]]

    max_t = 2 * limit_n * limit_n - 1
    max_prime = math.isqrt(max_t)
    primes = sieve_primes(max_prime)

    composite = bytearray(limit_n + 1)

    for p in primes:
        if p == 2:
            continue
        mod8 = p & 7
        if mod8 != 1 and mod8 != 7:
            continue

        inv2 = (p + 1) // 2
        if mod8 == 7:
            sqrt2 = mod_pow(2, (p + 1) // 4, p)
            root = (sqrt2 * inv2) % p
        else:
            root = tonelli_shanks(inv2, p)
            if root == 0:
                continue

        r1 = root
        r2 = (p - r1) if r1 != 0 else 0

        half = (p + 1) // 2
        sr = math.isqrt(half)
        skip_n = sr if sr * sr == half else -1

        for residue in (r1, r2):
            if residue == r2 and r2 == r1:
                continue
            n = residue
            if n < 2:
                delta = 2 - n
                n += ((delta + p - 1) // p) * p
            while n <= limit_n:
                if n != skip_n:
                    composite[n] = 1
                n += p

    count = 0
    for n in range(2, limit_n + 1):
        if not composite[n]:
            count += 1
    return str(count)

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

Java

public class Euler216 {
    public static void main(String[] args) {
        int limitN = 50000000;
        long maxT = 2L * limitN * limitN - 1;
        int maxPrime = (int) Math.sqrt(maxT);
        boolean[] sieve = new boolean[maxPrime + 1];
        for (int i = 2; (long) i * i <= maxPrime; i++)
            if (!sieve[i])
                for (int j = i * i; j <= maxPrime; j += i)
                    sieve[j] = true;
        int[] primes = new int[maxPrime / 5];
        int pc = 0;
        for (int p = 2; p <= maxPrime; p++)
            if (!sieve[p])
                primes[pc++] = p;

        boolean[] composite = new boolean[limitN + 1];
        for (int pi = 0; pi < pc; pi++) {
            int p = primes[pi];
            if (p == 2)
                continue;
            int mod8 = p & 7;
            if (mod8 != 1 && mod8 != 7)
                continue;
            long inv2 = (p + 1L) / 2;
            long root;
            if (mod8 == 7) {
                long sqrt2 = modPow(2, (p + 1L) / 4, p);
                root = sqrt2 * inv2 % p;
            } else {
                root = tonelliShanks(inv2, p);
                if (root == 0)
                    continue;
            }
            int r1 = (int) root, r2 = r1 == 0 ? 0 : p - r1;
            long half = (p + 1L) / 2;
            int sr = (int) Math.sqrt(half);
            while ((long) (sr + 1) * (sr + 1) <= half)
                sr++;
            while ((long) sr * sr > half)
                sr--;
            int skipN = ((long) sr * sr == half) ? sr : -1;
            mark(composite, r1, p, limitN, skipN);
            if (r2 != r1)
                mark(composite, r2, p, limitN, skipN);
        }
        int count = 0;
        for (int n = 2; n <= limitN; n++)
            if (!composite[n])
                count++;
        System.out.println(count);
    }

    static void mark(boolean[] composite, int residue, int p, int limitN, int skipN) {
        int n = residue;
        if (n < 2)
            n += ((2 - n + p - 1) / p) * p;
        for (; n <= limitN; n += p)
            if (n != skipN)
                composite[n] = true;
    }

    static long modPow(long base, long exp, long mod) {
        long result = 1;
        base %= mod;
        while (exp > 0) {
            if ((exp & 1) == 1)
                result = result * base % mod;
            base = base * base % mod;
            exp >>= 1;
        }
        return result;
    }

    static long tonelliShanks(long n, int p) {
        if (n == 0)
            return 0;
        if (p == 2)
            return n;
        if (modPow(n, (p - 1) / 2, p) != 1)
            return 0;
        if (p % 4 == 3)
            return modPow(n, (p + 1) / 4, p);
        long q = p - 1;
        int s = 0;
        while ((q & 1) == 0) {
            q >>= 1;
            s++;
        }
        long z = 2;
        while (modPow(z, (p - 1) / 2, p) != p - 1)
            z++;
        long m = s, c = modPow(z, q, p), t = modPow(n, q, p), r = modPow(n, (q + 1) / 2, p);
        while (t != 1) {
            long i = 1, t2 = t * t % p;
            while (t2 != 1) {
                t2 = t2 * t2 % p;
                i++;
            }
            long b = modPow(c, 1L << (m - i - 1), p);
            r = r * b % p;
            t = t * b % p * b % p;
            c = b * b % p;
            m = i;
        }
        return r;
    }
}