notes.html
problem 73 · 32 lines
<article>
<h3>Problem Analysis</h3>
<p>
I am asked to count the reduced proper fractions n/d - fractions with n < d and gcd(n, d) = 1 - that lie strictly between 1/3 and 1/2 when the denominator is bounded by some maximum D (here D = 12000). The natural setting is the <a class="reference-source" href="https://en.wikipedia.org/wiki/Farey_sequence" rel="noopener noreferrer" target="_blank">Farey sequence</a> F_D: the ordered list of every reduced fraction in [0, 1] whose denominator is at most D. Two structural facts about <a class="reference-source" href="https://en.wikipedia.org/wiki/Farey_sequence" rel="noopener noreferrer" target="_blank">Farey sequences</a> make this tractable without ever factoring anything. First, the <a class="reference-source" href="https://en.wikipedia.org/wiki/Mediant_(mathematics)" rel="noopener noreferrer" target="_blank">mediant</a> property: if a/b and c/d are adjacent in some F_n, the first fraction to appear between them as n grows is their mediant (a + c)/(b + d). Second, the neighbour identity: adjacent fractions a/b and c/d satisfy bc - ad = 1. These two facts let me walk the interval (1/3, 1/2) one fraction at a time.
</p>
<p>
There is a faster, completely different angle: counting by rank. The rank of a fraction p/q in F_D is its position - the number of reduced fractions in [0, 1] with denominator at most D that do not exceed p/q. If I can compute rank(1/2) and rank(1/3) cheaply, the answer is simply rank(1/2) - rank(1/3) - 1 (the -1 removes 1/2 itself, which is a boundary and not counted). A rank reduces to counting coprime pairs, which is exactly the job of <a class="reference-source" href="https://en.wikipedia.org/wiki/M%C3%B6bius_function" rel="noopener noreferrer" target="_blank">Mobius inversion</a> implemented as a divisor sieve. That yields an O(D log D) algorithm. The two enumeration approaches, by contrast, do O(1) work per counted fraction and therefore run in O(answer) time, which grows like O(D^2). At D = 12000 the answer is about 7.3 million; at D = 100000 it is about 500 million - which is precisely why the sieve wins so decisively.
</p>
<h3>Solution Approach</h3>
<h4>Algorithm s0: Iterative Farey Neighbour Traversal</h4>
<p>
This approach (<code>p0073_s0.py</code> and <code>p0073_s0.c</code>) walks the <a class="reference-source" href="https://en.wikipedia.org/wiki/Farey_sequence" rel="noopener noreferrer" target="_blank">Farey sequence</a> one fraction at a time using the neighbour recurrence. Given two consecutive fractions whose denominators are <code>prev_d</code> and <code>d</code>, the denominator of the next one is <code>d_next = D - ((D + prev_d) mod d)</code>. This formula follows directly from the neighbour identity bc - ad = 1 combined with the requirement that the next denominator be as large as possible while staying at most D. The crucial simplification is that I track <em>only denominators</em>. Both boundary fractions 1/3 and 1/2 have numerator 1, and once the boundaries are fixed, the <a class="reference-source" href="https://en.wikipedia.org/wiki/Mediant_(mathematics)" rel="noopener noreferrer" target="_blank">mediant</a> structure forces every intermediate numerator uniquely - so the numerators carry no information and I drop them entirely.
</p>
<p>
The traversal starts at the fraction just above 1/3 with the largest denominator at most D. That denominator is <code>2 + 3 * floor((D - 2) / 3)</code> - the largest value of the form 2 + 3k not exceeding D - and its predecessor in the walk is 1/3, denominator 3. The loop counts one fraction per iteration and stops when <code>d</code> reaches 2, the denominator of the upper boundary 1/2, which is not counted. Each step is O(1), so the total cost is O(answer). The Python and C versions are line-for-line equivalent; at D = 12000 the C version runs in 0.0247 s against 0.3569 s for Python, a roughly 14x speedup purely from compilation.
</p>
<h4>Algorithm s1: Stern-Brocot Recursive Mediant Search</h4>
<p>
This approach (<code>p0073_s1.py</code> and <code>p0073_s1.c</code>) enumerates the same fractions but via the recursive structure of the <a class="reference-source" href="https://en.wikipedia.org/wiki/Stern%E2%80%93Brocot_tree" rel="noopener noreferrer" target="_blank">Stern-Brocot tree</a>. Given lower and upper boundaries with denominators <code>b_lo</code> and <code>b_hi</code>, their <a class="reference-source" href="https://en.wikipedia.org/wiki/Mediant_(mathematics)" rel="noopener noreferrer" target="_blank">mediant</a> has denominator <code>b_lo + b_hi</code>. If that sum exceeds D, no fraction with a small enough denominator lies between the boundaries and the function returns 0. Otherwise the mediant is one counted fraction, and I recurse into the two sub-intervals (lower, mediant) and (mediant, upper). The initial call is <code>recursion(3, 2)</code> because the boundaries are 1/3 and 1/2 and, as in s0, only denominators matter.
</p>
<p>
This is correct and elegant, but it generates one recursive call per counted fraction plus the failing leaves, so it too runs in O(answer) time with O(depth) stack usage. The C version stores D in a file-scope variable <code>g_max_d</code> to avoid threading it through every call. The Python version raises the recursion limit to 10^6, since the tree can become very deep at large D. Interestingly, the Python recursion is the slowest of all approaches at scale (48.98 s at D = 100000) because of Python's per-call overhead and lack of tail-call optimisation; in C, however, the recursion (0.8119 s) beats the iterator (1.7096 s) at the same scale.
</p>
<h4>Algorithm s2: Mobius Sieve Rank Computation</h4>
<p>
This is the fast algorithm (<code>p0073_s2.py</code> and <code>p0073_s2.c</code>). I define <code>rank(n, d)</code> as the number of reduced fractions with denominator at most D that are at most n/d, and return <code>rank(1, 2) - rank(1, 3) - 1</code>. To compute a rank, I first build an array where <code>data[i] = floor(i * n / d)</code>. At this point <code>data[i]</code> counts <em>all</em> fractions k/i (for k from 1 to floor(i*n/d)), including the ones not in lowest terms. To isolate the reduced fractions of each exact denominator, I run an additive sieve: for each i, I subtract <code>data[i]</code> from every multiple <code>data[j]</code> with j = 2i, 3i, .... This is <a class="reference-source" href="https://en.wikipedia.org/wiki/M%C3%B6bius_function" rel="noopener noreferrer" target="_blank">Mobius inversion</a> expressed as repeated subtraction: it strips out the contribution of fractions that share a common factor with their denominator, leaving <code>data[i]</code> equal to the count of fractions in lowest terms with denominator exactly i that are at most n/d. Summing the array gives the rank.
</p>
<p>
The inner loop visits D/i elements for each i, and the sum of D/i over all i is the <a class="reference-source" href="https://en.wikipedia.org/wiki/Harmonic_series_(mathematics)" rel="noopener noreferrer" target="_blank">harmonic series</a>, which is O(log D). The total work is therefore O(D log D) with O(D) memory. The C version allocates the array with <code>malloc</code> and uses <code>long long</code> to keep the running sum safely in range; the Python version builds the array with a list comprehension. This algorithm completes in 0.0028 s (C) and 0.1095 s (Python) at D = 100000 - the only approaches comfortable at that scale.
</p>
</article>