Project Euler Problem 880
(x,y) is called a nested radical pair if x and y are non-zero integers such that dfrac{x}{y} is not a cube of a rational
Solution
Answer: 522095328
Let
$$\alpha=\sqrt[3]{x},\qquad \beta=\sqrt[3]{y},$$
and suppose
$$\sqrt{\alpha+\beta} = \sqrt[3]{a}+\sqrt[3]{b}+\sqrt[3]{c}.$$
Squaring gives
$$\alpha+\beta = \sqrt[3]{a^2}+\sqrt[3]{b^2}+\sqrt[3]{c^2} +2\sqrt[3]{ab}+2\sqrt[3]{ac}+2\sqrt[3]{bc}.$$
The right-hand side initially contains six cube-root terms.
For it to collapse to only two cube-root terms, the six radicals must lie in the same cubic field
$$\mathbb Q(\sqrt[3]{k}),$$
so after extracting cubes we may write
$$\sqrt[3]{a}=A,\qquad \sqrt[3]{b}=B,t,\qquad \sqrt[3]{c}=C,t^2,$$
where $t^3=k$ and $A,B,C\in\mathbb Z$.
Then
$$(A+Bt+Ct^2)^2 = (A^2+2BCk) + (2AB+C^2k)t + (B^2+2AC)t^2.$$
To obtain exactly two cube-root terms, exactly one coefficient must vanish.
That produces two primitive families.
1. First family
Set
$$2AB+C^2k=0.$$
After solving the primitiveness constraints and removing cube scalings, every primitive pair is
$$x=4qs^3,\qquad y=p(9p+4s)^3,$$
with
$$q=2p+s,$$
subject to
- $p$ odd,
- $\gcd(p,s)=1$,
- $q>0$,
- $s\neq0$.
2. Second family
Set
$$B^2+2AC=0.$$
This yields
$$x=4p(9p+2s)^3,\qquad y=qs^3,$$
with
$$q=4p+s,$$
subject to
- $q$ odd,
- $\gcd(p,s)=1$,
- $q>0$,
- $s\neq0$.
Every nested radical pair arises uniquely from one of these primitive families, together with the square scaling
$$(x,y)\mapsto (k^2x,k^2y),$$
coming from
$$(a,b,c)\mapsto (ka,kb,kc).$$
Hence if $(x,y)$ is primitive and
$$M=\max(|x|,|y|),$$
then all admissible scaled pairs are obtained for
$$1\le k\le \sqrt{\frac{N}{M}}.$$
Their total contribution is
$$(|x|+|y|) \sum_{k\le \sqrt{N/M}} k^2.$$
Using
$$\sum_{i=1}^n i^2=\frac{n(n+1)(2n+1)}6,$$
we can compute $H(N)$ efficiently.
Python implementation
from __future__ import annotations
import math
N = 10**15
MOD = 1031**3 + 2
# ------------------------------------------------------------
# Integer cube root
# ------------------------------------------------------------
def icbrt(n: int) -> int:
"""Floor cube root."""
if n <= 1:
return n
x = int(round(n ** (1.0 / 3.0)))
while (x + 1) ** 3 <= n:
x += 1
while x ** 3 > n:
x -= 1
return x
# ------------------------------------------------------------
# Perfect cube test
# ------------------------------------------------------------
CUBE_RES_9 = {0, 1, 8}
CUBE_RES_7 = {0, 1, 6}
CUBE_RES_13 = {pow(i, 3, 13) for i in range(13)}
def is_cube_int(z: int) -> bool:
if z == 0:
return True
n = abs(z)
if n % 9 not in CUBE_RES_9:
return False
if n % 7 not in CUBE_RES_7:
return False
if n % 13 not in CUBE_RES_13:
return False
r = icbrt(n)
return r**3 == n
def is_rational_cube_ratio(x: int, y: int) -> bool:
g = math.gcd(abs(x), abs(y))
return is_cube_int(x // g) and is_cube_int(y // g)
# ------------------------------------------------------------
# Sum of squares
# ------------------------------------------------------------
def sumsq(k: int) -> int:
return k * (k + 1) * (2 * k + 1) // 6
# ------------------------------------------------------------
# Canonical ordering
# ------------------------------------------------------------
def canon_pair(x: int, y: int):
ax, ay = abs(x), abs(y)
if ax > ay or (ax == ay and x > y):
return (y, x)
return (x, y)
# ------------------------------------------------------------
# Primitive generators
# ------------------------------------------------------------
def iter_primitive_pairs(limit: int):
gcd = math.gcd
# ---------------- Family A ----------------
p = 1
while True:
if p * (p + 4) ** 3 > limit:
break
t = icbrt(limit // p)
s_hi = (t - 9 * p) // 4
s_lo = 1 - 2 * p
if s_hi >= s_lo:
for s in range(s_lo, s_hi + 1):
if s == 0:
continue
if gcd(p, s) != 1:
continue
q = 2 * p + s
if q <= 0:
continue
x = 4 * q * s**3
y = p * (9 * p + 4 * s) ** 3
if max(abs(x), abs(y)) <= limit:
yield canon_pair(x, y)
p += 2
# ---------------- Family B ----------------
p = 1
while True:
if 4 * p * (p + 2) ** 3 > limit:
break
t = icbrt(limit // (4 * p))
s_hi = (t - 9 * p) // 2
s_lo = 1 - 4 * p
if s_lo % 2 == 0:
s_lo += 1
if s_hi >= s_lo:
for s in range(s_lo, s_hi + 1, 2):
if s == 0:
continue
if gcd(p, s) != 1:
continue
q = 4 * p + s
if q <= 0:
continue
x = 4 * p * (9 * p + 2 * s) ** 3
y = q * s**3
if max(abs(x), abs(y)) <= limit:
yield canon_pair(x, y)
p += 1
# ------------------------------------------------------------
# Main computation
# ------------------------------------------------------------
def H_mod(limit: int) -> int:
seen = set()
total = 0
for x, y in iter_primitive_pairs(limit):
# Exclude rational cubes
if is_rational_cube_ratio(x, y):
continue
if (x, y) in seen:
continue
seen.add((x, y))
ax, ay = abs(x), abs(y)
# Number of square scalings
kmax = math.isqrt(limit // ay)
contribution = (
(ax + ay) % MOD
* (sumsq(kmax) % MOD)
)
total = (total + contribution) % MOD
return total
# ------------------------------------------------------------
# Verification
# ------------------------------------------------------------
assert H_mod(10**3) == 2535
print(H_mod(10**15))
Code walkthrough
icbrtcomputes exact integer cube roots.is_cube_intchecks whether an integer is a perfect cube using modular filters plus exact verification.iter_primitive_pairsgenerates all primitive nested radical pairs from the two parametrized families.- Each primitive pair generates infinitely many scaled pairs by multiplication with $k^2$.
sumsq(k)evaluates
$$1^2+2^2+\cdots+k^2.$$
H_modaccumulates all contributions modulo
$$1031^3+2.$$
The verification
assert H_mod(10**3) == 2535
matches the value given in the problem statement.
The final computation gives
$$H(10^{15}) \equiv 522095328 \pmod{1031^3+2}.$$
Therefore:
Answer: 522095328