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 Apr 19, 2026

Problem #0216

Level Level 07

Solved By 4,760

Languages C++, Python

Answer 5437849

Length 346 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 Foundation

Theorem 1 (Quadratic residue condition). If $p$ is an odd prime dividing $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 $n^2 \equiv 2^{-1} \pmod{p}$ . Since $n^2 \equiv 2^{-1}$ has a solution, $2^{-1}$ is a quadratic residue mod $p$ . Since $\left(\frac{a^{-1}}{p}\right) = \left(\frac{a}{p}\right)$ for $\gcd(a, p) = 1$ , we need $\left(\frac{2}{p}\right) = 1$ . By the second supplement to the law of quadratic reciprocity:

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

This equals $1$ if and only if $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 all integers $k$ . Moreover, $n_0$ and $p - n_0$ are the two distinct roots of $2x^2 \equiv 1 \pmod{p}$ (assuming they exist).

Proof. We have $t(n_0 + kp) = 2(n_0 + kp)^2 - 1 = 2n_0^2 + 4n_0 kp + 2k^2 p^2 - 1 \equiv 2n_0^2 - 1 \equiv 0 \pmod{p}$ . For the roots: if $n_0^2 \equiv 2^{-1} \pmod{p}$ , then $(p - n_0)^2 = p^2 - 2pn_0 + n_0^2 \equiv n_0^2 \equiv 2^{-1} \pmod{p}$ . Since $x^2 \equiv c \pmod{p}$ has at most 2 solutions (the polynomial $x^2 - c$ has degree 2 over $\mathbb{F}_p$ ), $n_0$ and $p - n_0$ are all the roots. $\square$

Lemma 1 (Tonelli-Shanks algorithm). Given $a$ with $\left(\frac{a}{p}\right) = 1$ , the Tonelli-Shanks algorithm computes $r$ with $r^2 \equiv a \pmod{p}$ in $O(\log^2 p)$ time.

Proof. Write $p - 1 = 2^s q$ with $q$ odd. The algorithm maintains invariants $t \equiv a^q \cdot c^{2^{s-M}} \pmod{p}$ (converging to 1) and $R^2 \equiv a \cdot t^{-1} \pmod{p}$ , where $c$ is derived from a quadratic non-residue. At each iteration, $M$ decreases by at least 1, so the algorithm terminates in at most $s \leq \log_2 p$ iterations. Each iteration uses $O(\log p)$ modular multiplications. $\square$

Theorem 3 (Sieve correctness). Maintain an array $T[n] = t(n)$ for $2 \leq n \leq N$ . For each prime $p \leq \sqrt{2N^2}$ with $p \equiv \pm 1 \pmod{8}$ : find the roots $n_0, p - n_0$ of $2x^2 \equiv 1 \pmod{p}$ , and for each root $r \in \{n_0, p - n_0\}$ , divide $T[n]$ by $p$ (repeatedly) for all $n \equiv r \pmod{p}$ in range. After the sieve, $t(n)$ is prime if and only if $T[n] > 1$ .

Proof. The sieve removes all prime factors $\leq \sqrt{2N^2}$ from $T[n]$ . If $t(n)$ is composite, it has a prime factor $\leq \sqrt{t(n)} \leq \sqrt{2N^2}$ , which the sieve removes, leaving $T[n] < t(n)$ . If $t(n)$ is prime, no factor is found, so $T[n] = t(n) > 1$ . The remaining case $T[n] = 1$ means $t(n)$ was fully factored by small primes, hence composite. $\square$

Editorial

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). We initialize t(n) values. We then sieve primes up to sqrt(2*N^2) using standard sieve. Finally, iterate over each prime p in primes.

Pseudocode

Initialize t(n) values
Sieve primes up to sqrt(2*N^2) using standard sieve
for each prime p in primes
Find r such that 2*r^2 ≡ 1 (mod p)
Sieve with root r
Sieve with root p - r
Count primes

Complexity Analysis

Time: $O(N \log \log N)$ for the sieve over polynomial values (analogous to the Sieve of Eratosthenes). The Tonelli-Shanks computation per prime is $O(\log^2 p)$ , negligible in total.
Space: $O(N)$ for the array $T$ , plus $O(\sqrt{N})$ for the prime sieve.

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 = (__int128)result * base % mod;
        base = (__int128)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 = (__int128)t * t % p;
        while (tmp != 1) { tmp = (__int128)tmp * tmp % p; i++; }
        long long b = c;
        for (long long j = 0; j < M - i - 1; j++) b = (__int128)b * b % p;
        M = i;
        c = (__int128)b * b % p;
        t = (__int128)t * c % p;
        R = (__int128)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 Foundation

Editorial

Pseudocode

Complexity Analysis

Answer

Code

Problem Statement

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

Mathematical Foundation

Editorial

Pseudocode

Complexity Analysis

Answer

Code

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