Detailed mathematical approach

Problem Summary

A Cardano triplet is a triple of positive integers \((a,b,c)\) satisfying

$$\sqrt[3]{a+b\sqrt{c}}+\sqrt[3]{a-b\sqrt{c}}=1.$$

The task is to count all such triplets with \(a+b+c\le L\). The implementation does not search directly in \((a,b,c)\)-space; it first converts the radical identity into a Diophantine equation and then enumerates the resulting factorization patterns.

Mathematical Approach

Let

$$x=\sqrt[3]{a+b\sqrt{c}},\qquad y=\sqrt[3]{a-b\sqrt{c}}.$$

Then \(x+y=1\). The whole solution comes from extracting arithmetic consequences of this identity.

From the Cube-Root Identity to an Integer Equation

Using

$$x^3+y^3=(x+y)^3-3xy(x+y),$$

we get

$$2a=x^3+y^3=1-3xy,$$

so

$$xy=\frac{1-2a}{3}.$$

Multiplying the two cube-root arguments also gives

$$x^3y^3=(a+b\sqrt{c})(a-b\sqrt{c})=a^2-b^2c,$$

hence

$$a^2-b^2c=\left(\frac{1-2a}{3}\right)^3.$$

After rearranging,

$$27b^2c=27a^2+(2a-1)^3=(a+1)^2(8a-1).$$

This factorization is the key algebraic simplification behind the code.

Why \(a=3u-1\)

The left-hand side is divisible by \(27\), so \((a+1)^2(8a-1)\) must also be divisible by \(27\). Modulo \(3\), this forces

$$a\equiv 2\pmod 3.$$

Therefore we write

$$a=3u-1,\qquad u\ge 1.$$

Then \(a+1=3u\) and \(8a-1=3(8u-3)\), so the previous equation becomes

$$b^2c=u^2(8u-3).$$

Thus every Cardano triplet corresponds to a positive integer \(u\) together with a decomposition of \(u^2(8u-3)\) into a square part and a residual part.

Complete Parameterization Counted by the Implementation

Set

$$m=8u-3.$$

We need positive integers \(b,c\) such that

$$b^2c=u^2m.$$

Let

$$g=\gcd(b,u),\qquad u=gh,\qquad b=gb_1,$$

so that \(\gcd(h,b_1)=1\). Substituting into \(b^2c=u^2m\) gives

$$g^2b_1^2c=g^2h^2m\qquad\Longrightarrow\qquad b_1^2c=h^2m.$$

Because \(\gcd(h,b_1)=1\), every prime dividing \(b_1\) must come from \(m\), so \(b_1^2\mid m\). Write

$$m=b_1^2r.$$

Then automatically

$$c=h^2r,\qquad b=\frac{u}{h}b_1.$$

Hence every solution enumerated by the code has the form

$$a=3u-1,\qquad m=8u-3=b_1^2r,\qquad h\mid u,\qquad \gcd(h,b_1)=1,$$

$$b=\frac{u}{h}b_1,\qquad c=h^2r.$$

Conversely, any such choice satisfies \(b^2c=u^2m\), so this parameterization is complete.

Local Bounds from the Budget Constraint

For fixed \(u\), the value of \(a=3u-1\) is already known. Define the remaining budget

$$R=L-a.$$

Any candidate must satisfy

$$b+c=\frac{u}{h}b_1+h^2r\le R.$$

The code derives two strong bounds for \(h\). Since \(u/h\ge 1\), we have \(b\ge b_1\), hence

$$c=h^2r\le R-b_1,$$

which yields

$$h\le\left\lfloor\sqrt{\frac{R-b_1}{r}}\right\rfloor.$$

Also \(h\ge 1\) implies \(c=h^2r\ge r\), so \(b\le R-r\). Therefore

$$\frac{u}{h}b_1\le R-r\qquad\Longrightarrow\qquad h\ge\left\lceil\frac{ub_1}{R-r}\right\rceil$$

whenever \(R-r\gt 0\). The implementation checks only divisors \(h\mid u\) inside this interval.

A Global Upper Bound for \(u\)

Before the main enumeration starts, the code computes a safe maximal \(u\) by lower-bounding \(b+c\). Let

$$x=\frac{u}{h}b_1,\qquad y=h^2r.$$

Then \(x+y=b+c\) and

$$x^2y=u^2b_1^2r=u^2(8u-3).$$

For fixed \(x^2y=K\), the minimum of \(x+y\) is \(3\sqrt[3]{K/4}\), by calculus or the weighted AM-GM inequality. Therefore every valid triplet must satisfy

$$a+b+c\ge 3u-1+3\sqrt[3]{\frac{u^2(8u-3)}{4}}.$$

If this lower bound already exceeds \(L\), then that \(u\) cannot contribute. Since the bound grows with \(u\), binary search yields the cutoff used by max_u_bound.

Worked Example: \((2,1,5)\)

Take \(u=1\). Then

$$a=3\cdot 1-1=2,\qquad m=8\cdot 1-3=5.$$

The only square divisor of \(m=5\) is \(b_1=1\), and the only divisor of \(u=1\) is \(h=1\). Thus

$$b=\frac{1}{1}\cdot 1=1,\qquad c=1^2\cdot 5=5.$$

This reproduces the example \((2,1,5)\), and the parameterization counts it exactly once.

How the Code Works

The implementation first computes max_u_bound(limit), then builds an SPF table for fast factorization of every \(u\), and a prime table up to \(\sqrt{8u_{\max}-3}\) for factoring \(m=8u-3\). For each \(u\), it enumerates all divisors \(h\mid u\), all square-divisor choices \(b_1\) of \(m\), applies the interval bounds for \(h\), checks \(\gcd(h,b_1)=1\), and finally verifies \(b+c\le R\). The outer loop over \(u\) is split across threads because different \(u\)-values are independent.

Complexity Analysis

Precomputation costs \(O(u_{\max})\) time and memory for the SPF table, plus a smaller prime sieve up to \(\sqrt{8u_{\max}-3}\). The counting phase is dominated by divisor generation: for each \(u\), the algorithm enumerates divisors of \(u\) and square divisors of \(8u-3\). A simple closed form is awkward because divisor counts fluctuate, but the effective runtime is far below naive search thanks to the global \(u\)-bound, the \(h\)-interval pruning, and the coprimality filter.

Further Reading

1. Problem page: https://projecteuler.net/problem=251
2. Smallest prime factor sieve: https://cp-algorithms.com/algebra/prime-sieve-linear.html
3. Divisor generation from prime exponents: https://en.wikipedia.org/wiki/Divisor_function
4. Arithmetic-geometric mean inequality: https://en.wikipedia.org/wiki/Inequality_of_arithmetic_and_geometric_means