notes.html

problem 23 · 28 lines

<article>
<h3>Problem Analysis</h3>
<p>
The problem hinges on the notion of an <a class="reference-source" href="https://en.wikipedia.org/wiki/Abundant_number" rel="noopener noreferrer" target="_blank">abundant number</a>: 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 <a class="reference-source" href="https://en.wikipedia.org/wiki/Aliquot_sum" rel="noopener noreferrer" target="_blank">aliquot sum</a>, 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 <a class="reference-source" href="https://en.wikipedia.org/wiki/Abundant_number" rel="noopener noreferrer" target="_blank">abundant number</a>s. 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.
</p>
<p>
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 <a class="reference-source" href="https://en.wikipedia.org/wiki/Aliquot_sum" rel="noopener noreferrer" target="_blank">aliquot sum</a>. Phase two enumerates all sums of two <a class="reference-source" href="https://en.wikipedia.org/wiki/Abundant_number" rel="noopener noreferrer" target="_blank">abundant number</a>s 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 <a class="reference-source" href="https://en.wikipedia.org/wiki/Trial_division" rel="noopener noreferrer" target="_blank">trial division</a> 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 <a class="reference-source" href="https://en.wikipedia.org/wiki/Abundant_number" rel="noopener noreferrer" target="_blank">abundant number</a>s 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.
</p>
<h3>Solution Approach</h3>
<h4>Algorithm 0: Sieve-Based Aliquot-Sum Construction (s0.py and s0.c)</h4>
<p>
My first approach computes the <a class="reference-source" href="https://en.wikipedia.org/wiki/Aliquot_sum" rel="noopener noreferrer" target="_blank">aliquot sum</a> for the entire range in a single sieve pass, inspired by the <a class="reference-source" href="https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes" rel="noopener noreferrer" target="_blank">Sieve of Eratosthenes</a>. Instead of asking "what are the divisors of n?" for each n, I invert the question: for each divisor <code>i</code> from 1 to limit/2, I add <code>i</code> to every multiple of <code>i</code> starting at <code>2*i</code>. Each divisor value is processed exactly once and contributes to all the numbers it divides. After the pass, <code>div_sums[n]</code> holds the sum of all proper divisors of <code>n</code>. 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 <a class="reference-source" href="https://en.wikipedia.org/wiki/Harmonic_series_(mathematics)" rel="noopener noreferrer" target="_blank">harmonic series</a>, giving O(N log N).
</p>
<p>
Once the sieve completes, I collect <a class="reference-source" href="https://en.wikipedia.org/wiki/Abundant_number" rel="noopener noreferrer" target="_blank">abundant number</a>s with a single comparison <code>div_sums[i] &gt; i</code>, naturally producing them in ascending order. The pairwise-sum phase iterates over all pairs <code>(i, j)</code> with <code>j &gt;= i</code>, marking each sum that lands within the limit in a flat boolean array, and finally sums the integers whose mark was never set.
</p>
<p>
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 <a class="reference-source" href="https://numpy.org/" rel="noopener noreferrer" target="_blank">NumPy</a>: the entire inner sieve loop becomes a single vectorised slice assignment <code>div_sums[2*i::i] += i</code>, <a class="reference-source" href="https://en.wikipedia.org/wiki/Abundant_number" rel="noopener noreferrer" target="_blank">abundant number</a>s are extracted with <code>np.flatnonzero(div_sums &gt; np.arange(limit + 1))</code>, 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 <a class="reference-source" href="https://numpy.org/" rel="noopener noreferrer" target="_blank">NumPy</a> code.
</p>
<h4>Algorithm 1: Trial Division per Number (s1.py and s1.c)</h4>
<p>
My second approach computes each number's <a class="reference-source" href="https://en.wikipedia.org/wiki/Aliquot_sum" rel="noopener noreferrer" target="_blank">aliquot sum</a> independently by <a class="reference-source" href="https://en.wikipedia.org/wiki/Trial_division" rel="noopener noreferrer" target="_blank">trial division</a>. For a given <code>n</code>, I start the accumulator at 1 (every n &gt; 1 has 1 as a proper divisor) and iterate <code>i</code> from 2 up to √n. Whenever <code>i</code> divides <code>n</code> I add both <code>i</code> and its cofactor <code>n/i</code>, taking care not to double-count when <code>i == n/i</code> (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.
</p>
<p>
The pairwise-sum phase is identical to algorithm 0: a sorted list of <a class="reference-source" href="https://en.wikipedia.org/wiki/Abundant_number" rel="noopener noreferrer" target="_blank">abundant number</a>s, 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 <a class="reference-source" href="https://en.wikipedia.org/wiki/Trial_division" rel="noopener noreferrer" target="_blank">trial division</a> repeats divisibility work that the sieve shares across multiples.
</p>
</article>