notes.html

problem 10 · 99 lines

<article>
<h3>Problem Analysis</h3>
<p>
    The task is to add up every prime below two million. The arithmetic is trivial — a single
    running sum — so the entire difficulty lies in <em>enumerating the primes</em> quickly. By the
    <a class="reference-source" href="https://en.wikipedia.org/wiki/Prime_number_theorem" rel="noopener noreferrer" target="_blank">prime number theorem</a>
    there are about <code>n / ln n</code> primes below <code>n</code>, so for <code>n = 2,000,000</code>
    we are summing roughly 148,933 numbers; the answer, 142,913,828,922, overflows 32 bits and must
    be accumulated in a 64-bit integer (Python integers are unbounded, so this only bites the C side).
  </p>
<p>
    The decisive question is how cheaply each prime can be found. Testing each candidate independently
    by
    <a class="reference-source" href="https://en.wikipedia.org/wiki/Trial_division" rel="noopener noreferrer" target="_blank">trial division</a>
    costs <code>O(√n)</code> per number and <code>O(n√n / ln n)</code> overall — correct
    but slow. A sieve flips the economics: instead of interrogating each number, it crosses off the
    multiples of primes it has already found, sharing work across the whole range. The classical
    <a class="reference-source" href="https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes" rel="noopener noreferrer" target="_blank">Sieve of Eratosthenes</a>
    runs in <code>O(n log log n)</code> time — effectively linear — which is the complexity class
    that makes this problem easy at two million and still comfortable at ten million. I implemented
    five solutions to contrast the approaches: a
    <a class="reference-source" href="https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes" rel="noopener noreferrer" target="_blank">Sieve of Eratosthenes</a>,
    a
    <a class="reference-source" href="https://en.wikipedia.org/wiki/Sieve_of_Sundaram" rel="noopener noreferrer" target="_blank">Sieve of Sundaram</a>,
    an unbounded incremental sieve, plain
    <a class="reference-source" href="https://en.wikipedia.org/wiki/Trial_division" rel="noopener noreferrer" target="_blank">trial division</a>,
    and a call into a specialised C library.
  </p>
<h3>Solution Approach</h3>
<h4>Solution 0 — Sieve of Eratosthenes</h4>
<p>
    My baseline is the
    <a class="reference-source" href="https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes" rel="noopener noreferrer" target="_blank">Sieve of Eratosthenes</a>.
    I allocate a <code>bytearray</code> of <code>max_num + 1</code> flags initialised to 1, clear
    indices 0 and 1, then for each <code>i</code> up to <code>√max_num</code> that is still
    marked prime I clear every multiple starting at <code>i*i</code> (smaller multiples were already
    cleared by smaller primes). The remaining set indices are the primes, and I sum them. The Python
    version leans on a key trick: the inner loop is replaced by the slice assignment
    <code>sieve[i*i::i] = bytearray(len(...))</code>, which clears a whole arithmetic progression in
    one C-level operation. The C sibling writes the inner loop explicitly over an
    <code>unsigned char</code> array. Time is <code>O(n log log n)</code>, space <code>O(n)</code>.
  </p>
<h4>Solution 1 — Sieve of Sundaram</h4>
<p>
    The
    <a class="reference-source" href="https://en.wikipedia.org/wiki/Sieve_of_Sundaram" rel="noopener noreferrer" target="_blank">Sieve of Sundaram</a>
    sieves only the odd numbers. It works on indices <code>i ≤ n = (max_num - 1) / 2</code>, where
    index <code>k</code> represents the odd number <code>2k + 1</code>. The insight is that an odd
    number <code>2k + 1</code> is composite exactly when <code>k</code> can be written as
    <code>i + j + 2ij</code> for some <code>1 ≤ i ≤ j</code>, because
    <code>2(i + j + 2ij) + 1 = (2i + 1)(2j + 1)</code>. So I mark every <code>i + j + 2ij</code>, then
    every unmarked index <code>i</code> yields the prime <code>2i + 1</code> (with 2 prepended). It is
    a tidy reformulation of the same sieving idea on a half-sized array; in practice its bound is a
    touch looser than Eratosthenes, around <code>O(n log n)</code>. Python and C implement the marking
    loop identically.
  </p>
<h4>Solution 2 — Unbounded incremental sieve</h4>
<p>
    Solution 2 is an unbounded
    <a class="reference-source" href="https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes" rel="noopener noreferrer" target="_blank">Sieve of Eratosthenes</a>
    expressed as a
    <a class="reference-source" href="https://docs.python.org/3/glossary.html#term-generator" rel="noopener noreferrer" target="_blank">generator function</a>.
    Rather than allocating the full table up front, it keeps a dictionary mapping the next composite
    each known prime will strike to that prime's step. When the running counter <code>n</code> is
    absent from the map it is prime, and I schedule its first multiple <code>n*n</code>; when it is
    present I advance that prime's entry to its next free odd multiple. This yields primes lazily in
    order, so <code>solve()</code> can stop the instant a prime reaches <code>max_num</code> without
    knowing the bound in advance. The C translation reproduces the dictionary with a hand-written
    hash table using
    <a class="reference-source" href="https://en.wikipedia.org/wiki/Open_addressing" rel="noopener noreferrer" target="_blank">open addressing</a>,
    including the cluster-rehash needed after a deletion. The time complexity matches the array sieve,
    <code>O(n log log n)</code>, but the constant factor is far worse because of per-number map churn.
  </p>
<h4>Solution 3 — Trial division</h4>
<p>
    Solution 3 is the deliberately naive baseline:
    <a class="reference-source" href="https://en.wikipedia.org/wiki/Trial_division" rel="noopener noreferrer" target="_blank">trial division</a>.
    For each candidate below <code>max_num</code> I test divisibility by 2 and then by odd numbers up
    to <code>√num</code>; the first divisor found proves compositeness. It needs no extra memory
    beyond the loop variables, which is its only virtue here. The cost is <code>O(n√n / ln n)</code>,
    and the timings below show exactly how badly that scales against the sieves. Python and C share
    the same <code>is_prime</code> structure.
  </p>
<h4>Solution 4 — primesieve library</h4>
<p>
    Solution 4 delegates to
    <a class="reference-source" href="https://github.com/kimwalisch/primesieve" rel="noopener noreferrer" target="_blank">primesieve</a>,
    a highly optimised segmented sieve that uses
    <a class="reference-source" href="https://en.wikipedia.org/wiki/Wheel_factorization" rel="noopener noreferrer" target="_blank">wheel factorization</a>
    and cache-friendly segmentation to enumerate primes at close to memory bandwidth. The Python code
    is a one-liner over <code>primesieve.primes(max_num - 1)</code>; the C code calls
    <code>primesieve_generate_primes</code> over <code>[2, max_num - 1]</code> and sums the array. Both
    pass <code>max_num - 1</code> so the upper endpoint stays strictly below <code>max_num</code>, as
    the problem demands. This is the fastest of the five and shows what a purpose-built library buys
    over a hand-rolled
    <a class="reference-source" href="https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes" rel="noopener noreferrer" target="_blank">Sieve of Eratosthenes</a>.
  </p>
</article>