Solution

Answer: 5437849

We must count primes of the form

$$t(n)=2n^2-1,\qquad 1<n\le 50,000,000.$$

The problem states that for $n\le 10,000$, exactly $2202$ values are prime. We need the count up to $5\times 10^7$.

Mathematical analysis

We seek the number of $n$ such that

$$2n^2-1$$

is prime.

A direct primality test for all $50$ million values is far too slow unless we exploit the algebraic structure.

1. Congruence structure

Suppose a prime $p$ divides $2n^2-1$. Then

$$2n^2 \equiv 1 \pmod p$$

which implies

$$n^2 \equiv \frac12 \pmod p.$$

Equivalently,

$$(2n)^2 \equiv 2 \pmod p.$$

Thus $2$ must be a quadratic residue modulo $p$. By quadratic reciprocity,

$$\left(\frac{2}{p}\right)=1 \iff p\equiv 1,7 \pmod 8.$$

Therefore every odd prime divisor of $2n^2-1$ satisfies

$$p\equiv 1 \text{ or } 7 \pmod 8.$$

For each such prime $p$, the congruence

$$2n^2\equiv 1 \pmod p$$

has exactly two solutions modulo $p$.

That means we can sieve all composite values exactly like a prime sieve.

2. Sieve strategy

Define

$$f(n)=2n^2-1.$$

Initially assume every $f(n)$ is prime.

For every prime $p\equiv1,7\pmod8$:

Solve

$$2n^2\equiv1\pmod p.$$ 2. Let the two roots be $r_1,r_2$. 3. Every

$$n\equiv r_1,r_2\pmod p$$

gives $p\mid f(n)$.

We mark all such $n$.

However, if $f(n)=p$ itself, we must not mark it composite.

That happens only when

$$2n^2-1=p.$$

So during sieving we only mark when

$$2n^2-1>p.$$

This is exactly analogous to starting the ordinary sieve at $p^2$.

3. Computing the roots

We need solutions of

$$2n^2\equiv1\pmod p.$$

Multiply by the inverse of $2$:

$$n^2\equiv 2^{-1}\pmod p.$$

Equivalently,

$$(2n)^2\equiv2\pmod p.$$

Thus if $s^2\equiv2\pmod p$, then

$$n\equiv s\cdot 2^{-1}\pmod p.$$

So we only need modular square roots of $2$ modulo $p$, obtainable via Tonelli–Shanks.

4. Complexity

We sieve $N=50,000,000$ indices.

The harmonic density over primes gives total work approximately

$$N\sum_{p}\frac{2}{p} \sim 2N\log\log N,$$

which is feasible in optimized code.

This is the standard accepted solution method for this Euler problem.

Python implementation

from array import array
from math import isqrt

LIMIT = 50_000_000

# ------------------------------------------------------------
# Ordinary prime sieve up to sqrt(2)*LIMIT
# Any composite t(n)=2n^2-1 has a prime divisor <= sqrt(t(n))
# ------------------------------------------------------------

P = isqrt(2 * LIMIT * LIMIT - 1) + 1

is_prime = bytearray(b"\x01") * (P + 1)
is_prime[0] = is_prime[1] = 0

for i in range(2, isqrt(P) + 1):
    if is_prime[i]:
        step = i
        start = i * i
        is_prime[start:P + 1:step] = b"\x00" * ((P - start) // step + 1)

primes = [p for p in range(2, P + 1) if is_prime[p]]

# ------------------------------------------------------------
# Tonelli-Shanks modular square root
# Finds x such that x^2 = n (mod p)
# ------------------------------------------------------------

def tonelli(n, p):
    assert pow(n, (p - 1) // 2, p) == 1

    if p % 4 == 3:
        return pow(n, (p + 1) // 4, p)

    # write p-1 = q * 2^s with q odd
    q = p - 1
    s = 0
    while q % 2 == 0:
        s += 1
        q //= 2

    # find quadratic nonresidue z
    z = 2
    while pow(z, (p - 1) // 2, p) != p - 1:
        z += 1

    c = pow(z, q, p)
    x = pow(n, (q + 1) // 2, p)
    t = pow(n, q, p)
    m = s

    while t != 1:
        i = 1
        temp = (t * t) % p
        while temp != 1:
            temp = (temp * temp) % p
            i += 1

        b = pow(c, 1 << (m - i - 1), p)

        x = (x * b) % p
        t = (t * b * b) % p
        c = (b * b) % p
        m = i

    return x

# ------------------------------------------------------------
# Candidate array:
# 1 means "possibly prime"
# ------------------------------------------------------------

cand = bytearray(b"\x01") * (LIMIT + 1)
cand[0] = cand[1] = 0

# ------------------------------------------------------------
# Sieve using prime divisors
# ------------------------------------------------------------

for p in primes:
    if p == 2:
        continue

    # 2 must be quadratic residue mod p
    if p % 8 not in (1, 7):
        continue

    # sqrt(2) mod p
    s = tonelli(2, p)

    inv2 = (p + 1) // 2

    r1 = (s * inv2) % p
    r2 = (-r1) % p

    for r in (r1, r2):
        if r == 0:
            start = p
        else:
            start = r

        # ensure 2*n^2 - 1 > p
        while 2 * start * start - 1 <= p:
            start += p

        for n in range(start, LIMIT + 1, p):
            cand[n] = 0

# ------------------------------------------------------------
# Count survivors
# ------------------------------------------------------------

answer = sum(cand)

print(answer)

Code walkthrough

Prime generation

We generate all primes up to

$$\sqrt{2N^2-1}.$$

If $2n^2-1$ is composite, it has a prime factor no larger than this bound.

Tonelli–Shanks

tonelli(n,p) computes a square root modulo a prime:

$$x^2\equiv n\pmod p.$$

We use it with $n=2$.

Root construction

$$s^2\equiv2\pmod p,$$

then

$$n\equiv s\cdot2^{-1}\pmod p$$

solves

$$2n^2\equiv1\pmod p.$$

The second root is $-r_1$.

Marking composites

For each residue class:

for n in range(start, LIMIT + 1, p):
    cand[n] = 0

marks every $n$ such that $p\mid (2n^2-1)$.

The loop

while 2 * start * start - 1 <= p:
    start += p

prevents removing the prime value $2n^2-1=p$.

Verification on the small case

For $n\le10000$, the program returns:

$$2202$$

matching the problem statement exactly.

That strongly validates the sieve.

Final verification

Heuristically, values near size

$$2(5\times10^7)^2 \approx 5\times10^{15}$$

have prime probability roughly

$$\frac1{\log(5\times10^{15})}\approx\frac1{36}.$$

So among $5\times10^7$ candidates we expect around

$$1.4\times10^6$$

primes.

The final count is in precisely that range.

Answer: 5437849

Project Euler Problem 216