Problem 87 — Prime Power Triples

Problem Analysis

We are asked to count the distinct integers below N = 50,000,000 that can be written as p² + q³ + r⁴, where p, q, and r are each prime numbers. The word "distinct" is the crux of the problem: many different triples (p, q, r) can produce the same sum, so we cannot simply multiply together the counts of valid bases - we must deduplicate. The fact that makes this tractable is that the upper bound N tightly constrains how large each prime can be in each role. Since p² alone must be below N, the largest prime base we could ever need is at most the square root of N, roughly 7071. The cube and fourth-power roles are far more restrictive: r⁴ < N means r is at most the fourth root of N (about 84), and q³ < N means q is at most the cube root (about 368). The single most demanding constraint is the square term, so one prime list up to sqrt(N) is a superset of every prime we need for all three roles.

The algorithm is a triple nested loop over (fourth-power prime, cube prime, square prime), computing each sum and inserting it into a deduplication structure. Naively this is O(P₄ · P₃ · P₂), where P_k is the number of primes whose k-th power fits under N; by the prime number theorem each P_k is approximately N^1/k / ln N. Because the prime lists are sorted in ascending order, we can break out of each inner loop the moment a partial sum exceeds the budget, which prunes the vast majority of iterations. Empirically the number of valid triples is only a few million, so the dominant costs end up being sieve construction and the deduplication structure, both linear in their respective sizes.

Solution Approach

Triple-loop enumeration with sieve-based prime generation (Python and C)

Both files implement the identical algorithm in two languages, so I describe them together and note the language-specific differences where they matter.

Step 1 - Generate primes with the Sieve of Sundaram. I generate all primes up to sqrt(N) using the Sieve of Sundaram. This sieve marks integers of the form i + j + 2ij (for positive integers i, j ≥ 1) in a half-sized array; any position k that remains unmarked corresponds to the odd prime 2k+1, and I prepend 2 by hand. I chose it over the more familiar Sieve of Eratosthenes largely as a variation - both run in similar time and produce the same primes. The single sqrt(N) bound serves all three roles, as argued above.

Step 2 - Pre-compute early-exit thresholds. Before looping I compute max_quadruple_cube = max_num - 4 (reserving room for the smallest square 2² = 4) and max_quadruple = max_quadruple_cube - 8 (reserving room additionally for the smallest cube 2³ = 8). These constants become the break conditions at the outer two loop levels, so I exit as soon as further iterations could not possibly yield a valid sum.

Step 3 - Triple nested loop. The outermost loop iterates over fourth powers, the middle over cubes, and the innermost over squares. At each level I test the running partial sum against its threshold and break (not continue) when it is exceeded. This is correct precisely because the prime lists are ascending: once a partial sum overshoots, every later term in that loop overshoots too.

Step 4 - Deduplication. In Python I add each valid sum to a set and return its length. In C I allocate a calloc-zeroed unsigned char array of length max_num used as a presence table, incrementing a counter only the first time a given index is set. The C array gives guaranteed O(1) lookup with no hashing overhead at the cost of about 50 MB of memory; the Python set uses memory proportional only to the number of distinct sums (~1 million) but pays hash-table overhead per insertion.