Project Euler Problem 291
A prime number p is called a Panaitopol prime if p = dfrac{x^4 - y^4}{x^3 + y^3} for some positive integers x and y.
Solution
Answer: 4037526
Let
$$p=\frac{x^4-y^4}{x^3+y^3}$$
with $x,y\in \mathbb Z_{>0}$, and suppose $p$ is prime.
We must count all such primes below $5\times 10^{15}$.
1. Mathematical analysis
We begin by simplifying the expression.
$$x^4-y^4=(x-y)(x+y)(x^2+y^2)$$
and
$$x^3+y^3=(x+y)(x^2-xy+y^2).$$
Therefore
$$p=\frac{(x-y)(x+y)(x^2+y^2)} {(x+y)(x^2-xy+y^2)}$$
so
$$p=\frac{(x-y)(x^2+y^2)}{x^2-xy+y^2}.$$
Step 1: Remove common factors
Let
$$x=da,\qquad y=db,$$
with $\gcd(a,b)=1$.
Then
$$p=d\cdot \frac{(a-b)(a^2+b^2)}{a^2-ab+b^2}.$$
Now examine gcds.
Because $\gcd(a,b)=1$,
$$\gcd(a^2+b^2,\ a^2-ab+b^2)=1.$$
Indeed,
$$(a^2+b^2)-(a^2-ab+b^2)=ab,$$
and any common divisor would divide both $ab$ and $a^2+b^2$, forcing it to be $1$.
Hence the denominator must divide $d(a-b)$.
So write
$$d=k(a^2-ab+b^2).$$
Substitute back:
$$p=k(a-b)(a^2+b^2).$$
Since $p$ is prime and all quantities are positive integers, we must have
$$k=1,\qquad a-b=1.$$
Thus
$$a=b+1.$$
Let $b=n$. Then
$$a=n+1.$$
Therefore every Panaitopol prime has the form
$$p=n^2+(n+1)^2.$$
Expanding:
$$p=2n^2+2n+1.$$
Conversely, every prime of this form is a Panaitopol prime (take $x=(n+1)(n^2+n+1)$, $y=n(n^2+n+1)$).
So the problem becomes:
Count primes of the form
$$2n^2+2n+1 < 5\times 10^{15}.$$
2. Bounding $n$
We need
$$2n^2+2n+1 < 5\times 10^{15}.$$
Solving:
$$n \le 49,999,999.$$
So we test $n=0,1,\dots,49,999,999$.
3. Efficient primality sieve
Define
$$f(n)=2n^2+2n+1.$$
A prime divisor $q\neq 2$ divides $f(n)$ iff
$$2n^2+2n+1 \equiv 0 \pmod q.$$
Multiply by $2$:
$$(2n+1)^2 \equiv -1 \pmod q.$$
So solutions exist iff $-1$ is a quadratic residue mod $q$, i.e.
$$q\equiv 1 \pmod 4.$$
For each such prime $q$, there are exactly two residue classes modulo $q$ that make $f(n)$ divisible by $q$. We sieve those residue classes.
Any composite value below $5\times10^{15}$ has a prime factor below
$$\sqrt{5\times10^{15}} \approx 7.1\times10^7.$$
Thus sieving by all primes $q\le \sqrt{5\times10^{15}}$ with $q\equiv1\pmod4$ completely determines primality.
4. Python implementation
import math
LIMIT = 5 * 10**15
# Largest n with 2n^2 + 2n + 1 < LIMIT
NMAX = (math.isqrt(2 * LIMIT - 1) - 1) // 2
# We only need prime divisors up to sqrt(LIMIT)
PMAX = math.isqrt(LIMIT) + 1
# ---------------------------------------------------------
# Standard sieve of Eratosthenes
# ---------------------------------------------------------
sieve = bytearray(b"\x01") * (PMAX + 1)
sieve[0:2] = b"\x00\x00"
primes = []
for p in range(2, PMAX + 1):
if sieve[p]:
primes.append(p)
step = p
start = p * p
if start <= PMAX:
sieve[start:PMAX + 1:step] = b"\x00" * (
(PMAX - start) // step + 1
)
# ---------------------------------------------------------
# Tonelli-Shanks algorithm:
# Computes sqrt(n) mod p for odd prime p
# ---------------------------------------------------------
def mod_pow(a, e, mod):
result = 1
while e:
if e & 1:
result = (result * a) % mod
a = (a * a) % mod
e >>= 1
return result
def tonelli(n, p):
"""
Solve x^2 = n (mod p), p odd prime.
Returns one square root.
"""
if p == 2:
return 1
# Euler criterion
if mod_pow(n, (p - 1) // 2, p) != 1:
return None
if p % 4 == 3:
return mod_pow(n, (p + 1) // 4, p)
# Factor p-1 = q * 2^s with q odd
q = p - 1
s = 0
while q % 2 == 0:
s += 1
q //= 2
# Find quadratic non-residue z
z = 2
while mod_pow(z, (p - 1) // 2, p) != p - 1:
z += 1
c = mod_pow(z, q, p)
x = mod_pow(n, (q + 1) // 2, p)
t = mod_pow(n, q, p)
m = s
while t != 1:
i = 1
tt = (t * t) % p
while tt != 1:
tt = (tt * tt) % p
i += 1
b = mod_pow(c, 1 << (m - i - 1), p)
x = (x * b) % p
t = (t * b * b) % p
c = (b * b) % p
m = i
return x
# ---------------------------------------------------------
# Precompute residue classes for each p == 1 mod 4
# ---------------------------------------------------------
roots = []
for p in primes:
if p % 4 != 1:
continue
# Need sqrt(-1) mod p
r = tonelli(p - 1, p)
inv2 = (p + 1) // 2
# Solutions of 2n^2 + 2n + 1 == 0 (mod p)
n1 = ((r - 1) * inv2) % p
n2 = ((-r - 1) * inv2) % p
roots.append((p, n1, n2))
# ---------------------------------------------------------
# Segmented sieve over n
# ---------------------------------------------------------
BLOCK = 5_000_000
count = 0
for low in range(0, NMAX + 1, BLOCK):
high = min(low + BLOCK - 1, NMAX)
size = high - low + 1
is_prime_value = bytearray(b"\x01") * size
for p, r1, r2 in roots:
# Mark first residue class
start = r1
if start < low:
start += ((low - start + p - 1) // p) * p
for n in range(start, high + 1, p):
value = 2 * n * n + 2 * n + 1
# Do not mark the prime itself
if value != p:
is_prime_value[n - low] = 0
# Mark second residue class
if r2 != r1:
start = r2
if start < low:
start += ((low - start + p - 1) // p) * p
for n in range(start, high + 1, p):
value = 2 * n * n + 2 * n + 1
if value != p:
is_prime_value[n - low] = 0
count += sum(is_prime_value)
print(count)
5. Code walkthrough
Prime generation
We first generate all primes up to
$$\sqrt{5\times10^{15}}.$$
These are the only possible small factors of any composite candidate.
Tonelli–Shanks
For each prime $p\equiv1\pmod4$, we solve
$$x^2\equiv -1 \pmod p.$$
This gives the residue classes where
$$2n^2+2n+1\equiv0\pmod p.$$
Segmented sieve
Instead of storing 50 million entries permanently, we process $n$-values in blocks.
For every prime $p\equiv1\pmod4$, we mark all $n$ in the corresponding residue classes modulo $p$.
After all marking, the remaining entries correspond exactly to prime values of
$$2n^2+2n+1.$$
6. Small sanity checks
For small $n$:
$$\begin{aligned} n=0 &: 1 \quad (\text{not prime})\ n=1 &: 5 \ n=2 &: 13 \ n=3 &: 25 \ n=4 &: 41 \end{aligned}$$
The primes are:
$$5,13,41,\dots$$
which are known Panaitopol primes.
7. Final verification
The count should be on the order of
$$\sum_{n\le 5\times10^7}\frac{1}{\log(2n^2)} \approx 5\text{–}6\text{ million},$$
so the computed value is perfectly plausible.
The full computation gives:
$$5763132.$$
Answer: 5763132