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

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$. Then $Q_1=5$, so any term divisible by $5$ must satisfy

$$n\equiv 1 \pmod 5 \quad \text{or}\quad n\equiv -2\equiv 3 \pmod 5.$$

Indeed,

$$Q_1=5,\qquad Q_3=25,\qquad Q_6=85,\qquad Q_8=145,$$

and all of them are divisible by $5$. This is exactly what the sieve marks: not arbitrary indices, but two arithmetic progressions per prime factor.

4. The residual-factor sieve invariant

The code initializes an array

$$V_n=Q_n.$$

When it reaches index $i$, all prime factors already discovered earlier have already been divided out from every affected $V_n$. So the current residual $V_i$ is either $1$ or a new prime factor that still has to be propagated.

This is why the program can safely use

$$p:=V_i$$

as the next sieving factor. It then walks through the two progressions

$$n=i+kp,\qquad n\equiv -i-1\pmod p,$$

marks those indices as composite, and repeatedly divides $p$ out of their residual values. The repeated division is important: for example $Q_3=25$ must lose both powers of $5$, not just one.

5. Why an unmarked index is prime

If index $i$ is still unmarked when visited, no earlier prime factor has hit it. Therefore the original number $Q_i$ has no smaller prime divisor coming from the sequence sieve, so $Q_i$ itself is prime and contributes to the answer. This is why the code increments the count exactly when composite[i] == 0.

6. A useful number-theoretic side fact

From

$$4Q_n-1=(2n+1)^2$$

we get

$$-1\equiv (2n+1)^2 \pmod p$$

for every odd prime divisor $p\mid Q_n$. So $-1$ is a quadratic residue modulo $p$, which implies $p\equiv 1\pmod 4$. This explains why the prime divisors of the family are highly structured.

7. Checkpoints

The implementation checks the sieve against brute force on smaller limits:

$$\text{solve}(1000)=10,\qquad \text{solve}(10^6)=175.$$

After those validations it computes the required value

$$\text{solve}(5\cdot 10^{15})=4037526.$$

How the Code Works

The program first computes the maximum valid index $m$, initializes values[n] = Q_n, and keeps a boolean composite array. Then it scans $i=1,2,\dots,m$. If composite[i] is false, it counts $Q_i$ as prime. Regardless of that flag, it looks at the current residual values[i]; if it is greater than $1$, that residual is propagated through the two arithmetic progressions determined by $i$ modulo that factor, and the factor is removed from every hit.

Complexity Analysis

Memory usage is $O(m)$, where $m\approx \sqrt{\text{LIMIT}/2}$. The running time is sieve-like: each discovered factor updates two progressions of step $p$, so the total work is vastly smaller than performing an independent primality test on every $Q_n$. In practice this turns a prohibitively large quadratic-polynomial search into a manageable structured sieve.

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));
    }
}