Project Euler

Investigating the Primality of $2n^2 - 1$

Consider the sequence t(n) = 2n^2 - 1 for n >= 1: 1, 7, 17, 31, 49, 71, 97, 127, 161, 199,... How many values of t(n) are prime for 2 <= n <= 50,000,000?

Source sync May 21, 2026

Problem #0216

Level Level 07

Solved By 4,760

Languages C++, Python

Answer 5437849

Length 340 words

modular_arithmeticnumber_theoryalgebra

Consider numbers $t(n)$ of the form $t(n) = 2n^2 - 1$ with $n > 1$.

The first such numbers are $7, 17, 31, 49, 71, 97, 127$ and $161$.

It turns out that only $49 = 7 \cdot 7$ and $161 = 7 \cdot 23$ are not prime.

For $n \le 10000$ there are $2202$ numbers $t(n)$ that are prime.

How many numbers $t(n)$ are prime for $n \le 50\,000\,000$?

Problem 216: Investigating the Primality of $2n^2 - 1$

Mathematical Development

Theorem 1 (Quadratic residue condition). If an odd prime $p$ divides $t(n) = 2n^2 - 1$ , then

p \equiv \pm 1 \pmod{8}.

Proof. From $p \mid 2n^2 - 1$ we get

2n^2 \equiv 1 \pmod{p},

so $2^{-1}$ is a quadratic residue modulo $p$ . Equivalently,

\left(\frac{2}{p}\right) = 1.

By the second supplement to quadratic reciprocity,

\left(\frac{2}{p}\right) = (-1)^{(p^2 - 1)/8},

which equals 1 exactly for $p \equiv \pm 1 \pmod{8}$ . $\square$

Theorem 2 (Root periodicity). If $p \mid t(n_0)$ , then

p \mid t(n_0 + kp)

for every integer $k$ . Moreover, when the congruence

2x^2 \equiv 1 \pmod{p}

has a solution, its two roots are $r$ and $p-r$ .

Proof. Expanding gives

t(n_0 + kp) = 2(n_0 + kp)^2 - 1 \equiv 2n_0^2 - 1 \equiv 0 \pmod{p}.

If $r^2 \equiv 2^{-1} \pmod{p}$ , then $(p-r)^2 \equiv r^2 \pmod{p}$ as well. Since a quadratic congruence over $\mathbb{F}_p$ has at most two roots, these are the only solutions. $\square$

Lemma 1 (Tonelli-Shanks). For an odd prime $p$ and a quadratic residue $a \pmod p$ , the Tonelli-Shanks algorithm finds a square root of $a$ modulo $p$ in polylogarithmic time.

Theorem 3 (Boolean residue-class sieve). Maintain an array is_composite[n], initially false for $2 \leq n \leq N$ . For each prime

p \leq \sqrt{2N^2}

with $p \equiv \pm 1 \pmod{8}$ : solve

2x^2 \equiv 1 \pmod{p},

obtaining roots $r$ and $p-r$ , and mark every $n \equiv r \pmod p$ or $n \equiv -r \pmod p$ as composite, except for the unique case $t(n) = p$ . After all such primes are processed, is_composite[n] is false exactly when $t(n)$ is prime.

Proof. If $n$ is marked for some prime $p$ , then $p \mid t(n)$ . The skipped exceptional case is precisely $t(n)=p$ , which is prime; every other marked value is composite.

Conversely, if $t(n)$ is composite, it has some prime divisor

q \leq \sqrt{t(n)} \leq \sqrt{2N^2}.

By Theorem 1, $q \equiv \pm 1 \pmod{8}$ , so the congruence $2x^2 \equiv 1 \pmod q$ has the relevant roots, and by Theorem 2 the value of $n$ lies in one of the marked residue classes modulo $q$ . Therefore every composite $t(n)$ is marked. $\square$

Editorial

The sieve works on the parameter $n$ , not on the values $2n^2-1$ themselves. For a prime $p$ to divide $2n^2-1$ , the congruence

2n^2 \equiv 1 \pmod p

must be solvable, so only primes with $p \equiv \pm 1 \pmod 8$ matter. Once such a prime is fixed, its solutions form two residue classes modulo $p$ , and every $n$ in those classes produces a composite value unless $2n^2-1$ happens to equal $p$ itself.

That turns the problem into a quadratic-residue sieve. Generate the relevant primes, use Tonelli-Shanks to find the two roots of the congruence, mark the corresponding arithmetic progressions in $n$ , and count the values left unmarked.

Pseudocode

Set LIMIT = 50,000,000.
Sieve all primes up to floor(sqrt(2 * LIMIT^2)).
Create a boolean array is_composite[0..LIMIT], initially false.

For each odd prime p from the sieve:
    If p is not congruent to 1 or 7 modulo 8:
        continue

    Solve 2 * r^2 ≡ 1 (mod p) using Tonelli-Shanks.
    Let the two roots be r and p - r.

    Check whether 2 * n^2 - 1 = p has an integer solution n.
    If it does, remember that exceptional n so it is not marked.

    For each root root in {r, p - r}:
        Mark every n ≡ root (mod p) with 2 <= n <= LIMIT
        as composite, except the exceptional n if it exists.

Count the integers n from 2 to LIMIT with is_composite[n] = false.
Return that count.

Complexity Analysis

Time: The prime sieve costs $O(M \log \log M)$ with $M = \lfloor \sqrt{2}N \rfloor$ . The marking phase costs $O\!\left(\sum_{p \equiv \pm 1 \bmod 8} \frac{N}{p}\right) = O(N \log \log N).$
Space: $O(M)$ for the prime sieve and $O(N)$ for the composite-marker array.

Answer

\boxed{5437849}

C++ project_euler/problem_216/solution.cpp

#include <bits/stdc++.h>
using namespace std;

// Problem 216: Count primes in t(n) = 2n^2 - 1 for 2 <= n <= 50,000,000
// Uses a sieve approach: for each prime p, find roots of 2x^2 = 1 (mod p)
// and sieve out multiples.

const int N = 50000000;

// We store t[n] and divide out small prime factors via sieve.
// After sieving, if remaining value > 1, t(n) is prime.
// We use long long since t(N) ~ 5e15.

// For memory, we store the "remaining" factor of t(n).
// t(n) = 2n^2 - 1. We sieve primes p up to sqrt(2*N^2) ~ 1e8.

// Actually, a simpler approach: store t[n] as long long, sieve by dividing.
// But 50M long longs = 400MB, too much.

// Better approach: segmented sieve on n. For each segment of n-values,
// compute t(n), sieve small primes, check if remainder is 1 or t(n).

// Even simpler: use a boolean array. Mark n as "composite t(n)" when we find
// a prime dividing t(n). But t(n) could have all prime factors > sieve limit
// only if t(n) is prime or a product of two large primes.
// If t(n) has a factor <= sqrt(t(n)) ~ 7e7, we'll find it.
// Sieve limit: sqrt(2 * 50000000^2 - 1) ~ 70710678

// We sieve primes up to ~70710678 using standard sieve, then for each such prime p,
// find roots of 2x^2 = 1 mod p, and mark those n as having a small factor.
// But this doesn't fully work because t(n) might be a prime power or product of
// primes all > sqrt(t(n)). Actually if t(n) = a*b with a,b > sqrt(t(n)), impossible.
// So if no prime <= sqrt(t(n)) divides t(n), then t(n) is prime.
// sqrt(t(n)) varies with n. sqrt(t(N)) ~ 7.07e7.
// For smaller n, sqrt(t(n)) is smaller, so sieving up to 7.07e7 covers everything.

// Memory for boolean sieve of primes up to 7.07e7: ~70MB (bitset or vector<bool>).
// Boolean array for n: 50MB.

// Approach:
// 1) Sieve primes up to ~70710678.
// 2) For each prime p, find roots r1, r2 of 2x^2 = 1 mod p.
//    This requires p | (2x^2 - 1), so 2 must be a QR mod p, i.e., p = +/-1 mod 8.
//    (But we also need p odd and p >= 3.)
//    Actually, p=2: t(n) = 2n^2-1 is always odd, so p=2 never divides t(n).
// 3) For each root r, mark r, r+p, r+2p, ... as composite.
//    Also mark p-r, p-r+p, ... as composite.
// 4) Count n in [2, N] not marked.

// But wait: marking n as composite just because p | t(n) is wrong if p^2 | t(n)
// or if t(n) = p. We need: t(n) is NOT prime, so we should mark n only if t(n) != p.
// Actually, if p | t(n) and t(n) > p, then t(n) is composite.
// If t(n) = p, then t(n) is prime and we should NOT mark it.
// t(n) = p means 2n^2 - 1 = p, so n = sqrt((p+1)/2). This only happens for
// specific small n/p values. We can handle this edge case.

// For p | t(n) with p < t(n): t(n) is composite.
// For p | t(n) with p = t(n): t(n) is prime. Don't mark.
// t(n) = p iff 2n^2 - 1 = p, i.e., n^2 = (p+1)/2.

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

// Tonelli-Shanks: find x such that x^2 = a mod p
long long sqrt_mod(long long a, long long p) {
    if (a == 0) return 0;
    if (p == 2) return a & 1;

    // Check if a is QR
    if (modpow(a, (p - 1) / 2, p) != 1) return -1;

    if (p % 4 == 3) {
        return modpow(a, (p + 1) / 4, p);
    }

    // Factor p-1 = 2^s * q
    long long s = 0, q = p - 1;
    while (q % 2 == 0) { s++; q /= 2; }

    // Find non-residue
    long long z = 2;
    while (modpow(z, (p - 1) / 2, p) != p - 1) z++;

    long long M = s;
    long long c = modpow(z, q, p);
    long long t = modpow(a, q, p);
    long long R = modpow(a, (q + 1) / 2, p);

    while (true) {
        if (t == 1) return R;
        long long i = 1;
        long long tmp = t * t % p;
        while (tmp != 1) { tmp = tmp * tmp % p; i++; }
        long long b = c;
        for (long long j = 0; j < M - i - 1; j++) b = b * b % p;
        M = i;
        c = b * b % p;
        t = t * c % p;
        R = R * b % p;
    }
}

int main() {
    const int LIMIT = N;
    const int SIEVE_LIMIT = 70710679; // > sqrt(2 * N^2)

    // Step 1: Sieve primes up to SIEVE_LIMIT
    vector<bool> is_prime_sieve(SIEVE_LIMIT + 1, true);
    is_prime_sieve[0] = is_prime_sieve[1] = false;
    for (long long i = 2; i * i <= SIEVE_LIMIT; i++) {
        if (is_prime_sieve[i]) {
            for (long long j = i * i; j <= SIEVE_LIMIT; j += i)
                is_prime_sieve[j] = false;
        }
    }

    // Step 2: For each prime, sieve t(n)
    // is_composite[n] = true means t(n) is composite
    vector<bool> is_composite(LIMIT + 1, false);

    for (long long p = 3; p <= SIEVE_LIMIT; p++) {
        if (!is_prime_sieve[p]) continue;
        // Check if 2 is QR mod p: p = +/-1 mod 8
        if (p % 8 != 1 && p % 8 != 7) continue;

        // Find n such that 2n^2 = 1 mod p, i.e., n^2 = (p+1)/2 mod p
        long long half_inv = (p + 1) / 2; // This is 2^{-1} mod p
        long long r = sqrt_mod(half_inv, p);
        if (r < 0) continue;

        // Two roots: r and p - r
        // For each root, sieve n = root, root+p, root+2p, ...
        // But skip if t(n) = p (meaning n^2 = (p+1)/2 exactly)

        long long p_check = -1;
        // t(n) = p iff 2n^2-1 = p iff n^2 = (p+1)/2
        {
            long long val = (p + 1) / 2;
            long long sq = (long long)sqrt((double)val);
            while (sq * sq < val) sq++;
            while (sq * sq > val) sq--;
            if (sq * sq == val && sq >= 2 && sq <= LIMIT) p_check = sq;
        }

        for (int sign = 0; sign < 2; sign++) {
            long long root = (sign == 0) ? r : p - r;
            if (root == 0) continue;
            // start from root, step by p
            for (long long n = root; n <= LIMIT; n += p) {
                if (n >= 2 && n != p_check) {
                    is_composite[n] = true;
                }
            }
        }
    }

    // Step 3: Count primes
    int count = 0;
    for (int n = 2; n <= LIMIT; n++) {
        if (!is_composite[n]) count++;
    }

    cout << count << endl;
    return 0;
}

Python project_euler/problem_216/solution.py

"""
Problem 216: Investigating the Primality of Values of 2n^2 - 1

Count how many t(n) = 2n^2 - 1 are prime for 2 <= n <= 50,000,000.

Approach: Sieve - for each prime p (where 2 is a QR mod p), find roots of
2x^2 = 1 (mod p) and mark those n-values as having composite t(n).

Answer: 5437849
"""

import math

def solve():
    LIMIT = 50_000_000
    SIEVE_LIMIT = int(math.isqrt(2 * LIMIT * LIMIT)) + 1

    # Sieve of Eratosthenes for primes up to SIEVE_LIMIT
    is_prime = bytearray(b'\x01') * (SIEVE_LIMIT + 1)
    is_prime[0] = is_prime[1] = 0
    for i in range(2, int(math.isqrt(SIEVE_LIMIT)) + 1):
        if is_prime[i]:
            is_prime[i*i::i] = bytearray(len(is_prime[i*i::i]))

    # is_composite[n] = 1 means t(n) is composite
    is_composite = bytearray(LIMIT + 1)

    def modpow(base, exp, mod):
        result = 1
        base %= mod
        while exp > 0:
            if exp & 1:
                result = result * base % mod
            base = base * base % mod
            exp >>= 1
        return result

    def tonelli_shanks(a, p):
        """Find x such that x^2 = a mod p."""
        if a == 0:
            return 0
        if p == 2:
            return a % 2

        # Check QR
        if modpow(a, (p - 1) // 2, p) != 1:
            return -1

        if p % 4 == 3:
            return modpow(a, (p + 1) // 4, p)

        # Factor p-1 = 2^s * q
        s, q = 0, p - 1
        while q % 2 == 0:
            s += 1
            q //= 2

        # Find non-residue
        z = 2
        while modpow(z, (p - 1) // 2, p) != p - 1:
            z += 1

        M = s
        c = modpow(z, q, p)
        t = modpow(a, q, p)
        R = modpow(a, (q + 1) // 2, p)

        while True:
            if t == 1:
                return R
            i = 1
            tmp = t * t % p
            while tmp != 1:
                tmp = tmp * tmp % p
                i += 1
            b = c
            for _ in range(M - i - 1):
                b = b * b % p
            M = i
            c = b * b % p
            t = t * c % p
            R = R * b % p

    for p in range(3, SIEVE_LIMIT + 1):
        if not is_prime[p]:
            continue
        # 2 is QR mod p iff p = +/-1 mod 8
        if p % 8 not in (1, 7):
            continue

        # Find root of n^2 = (p+1)/2 mod p
        half_inv = (p + 1) // 2
        r = tonelli_shanks(half_inv, p)
        if r < 0:
            continue

        # Check if t(n) = p for some n (don't mark those)
        p_check = -1
        val = (p + 1) // 2
        sq = int(math.isqrt(val))
        if sq * sq == val and 2 <= sq <= LIMIT:
            p_check = sq

        for root in (r, p - r):
            if root == 0:
                continue
            # Mark all n = root, root+p, root+2p, ... as composite
            start = root if root >= 2 else root + p * ((2 - root + p - 1) // p)
            for n in range(start, LIMIT + 1, p):
                if n != p_check:
                    is_composite[n] = 1

    count = 0
    for n in range(2, LIMIT + 1):
        if not is_composite[n]:
            count += 1

    print(count)

if __name__ == "__main__":
    solve()

Investigating the Primality of $2n^2 - 1$

Problem Statement

Problem 216: Investigating the Primality of $2n^2 - 1$

Mathematical Development

Editorial

Pseudocode

Complexity Analysis

Answer

Code

Problem Statement

Problem 216: Investigating the Primality of 2n2−12n^2 - 12n2−1

Mathematical Development

Editorial

Pseudocode

Complexity Analysis

Answer

Code

Problem 216: Investigating the Primality of $2n^2 - 1$