Problem Analysis
The problem hinges on the notion of an abundant number: a positive integer whose proper divisors (every divisor except the number itself) sum to more than the number. That sum of proper divisors is called the aliquot sum, and it is the central quantity I have to compute. The problem hands me a crucial mathematical fact: every integer greater than 28123 can be expressed as the sum of two abundant numbers. This bounds an otherwise infinite search into the finite range [1, 28123], which is what makes the problem tractable at all. Without that bound there would be no terminating computation.
With the range fixed, the problem decomposes into two clean phases. Phase one classifies every number up to the limit as abundant or not by computing its aliquot sum. Phase two enumerates all sums of two abundant numbers and marks which targets in the range are reachable; the answer is the sum of the unmarked targets. The classification phase can be done either by per-number trial division in O(N√N) total time, or by a sieve in O(N log N) time. The pairwise-sum phase costs O(A²) where A is the number of abundant numbers below the limit - roughly 6965, giving about 24 million pair checks, which a modern machine dispatches in milliseconds. Recognising the two-phase structure is the key to optimising each part independently.
Solution Approach
Algorithm 0: Sieve-Based Aliquot-Sum Construction (s0.py and s0.c)
My first approach computes the aliquot sum for the entire range in a single sieve pass, inspired by the Sieve of Eratosthenes. Instead of asking "what are the divisors of n?" for each n, I invert the question: for each divisor i from 1 to limit/2, I add i to every multiple of i starting at 2*i. Each divisor value is processed exactly once and contributes to all the numbers it divides. After the pass, div_sums[n] holds the sum of all proper divisors of n. The outer loop runs limit/2 times, the inner loop runs about limit/(2i) times per step, and the total work is proportional to the harmonic series, giving O(N log N).
Once the sieve completes, I collect abundant numbers with a single comparison div_sums[i] > i, naturally producing them in ascending order. The pairwise-sum phase iterates over all pairs (i, j) with j >= i, marking each sum that lands within the limit in a flat boolean array, and finally sums the integers whose mark was never set.
The two language versions differ in how the inner loops are expressed. In C I use plain nested loops with direct array indexing - cache-friendly and idiomatic. In Python I lean on NumPy: the entire inner sieve loop becomes a single vectorised slice assignment div_sums[2*i::i] += i, abundant numbers are extracted with np.flatnonzero(div_sums > np.arange(limit + 1)), and even the sum-marking loop uses vectorised slicing and boolean indexing. The Python interpreter only drives the outer loops; the heavy lifting runs in compiled NumPy code.
Algorithm 1: Trial Division per Number (s1.py and s1.c)
My second approach computes each number's aliquot sum independently by trial division. For a given n, I start the accumulator at 1 (every n > 1 has 1 as a proper divisor) and iterate i from 2 up to √n. Whenever i divides n I add both i and its cofactor n/i, taking care not to double-count when i == n/i (the case of a perfect square). This costs O(√n) per number and O(N√N) overall. The divisor-pairing trick - that divisors come in pairs straddling the square root - is what keeps the per-number cost at √n rather than n.
The pairwise-sum phase is identical to algorithm 0: a sorted list of abundant numbers, symmetric pair iteration, a boolean marking array, and a linear final scan. Only the classification phase differs. The C version outperforms the Python version dramatically because the tight square-root loop runs in native code; even so, it is slower than the sieve-based C solution because trial division repeats divisibility work that the sieve shares across multiples.