Problem 198: Ambiguous Numbers

Project Euler Problem 198 Solution

EulerSolve provides an optimized solution for Project Euler Problem 198, Ambiguous Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a denominator bound $B$, a rational $r/s$ with $s \le B$ is called a best approximation to a real number $x$ if it minimizes $|x-r/s|$ among all rationals whose denominators do not exceed $B$. A rational number $x$ is called ambiguous if there exists at least one bound $B$ for which two distinct rationals are tied as best approximations. The problem asks for the number of reduced rationals $x=p/q$ such that $0 \lt x \lt 1/100$, $q \le 10^8$, and $x$ is ambiguous. The successful strategy is not to enumerate all fractions with denominator up to $10^8$, but to characterize exactly which rationals can be ambiguous and count only those. Mathematical Approach The central objects are neighboring Farey fractions. Every ambiguous rational in the target range comes from one Farey gap, and the denominator condition on $x$ becomes a simple product condition on the two endpoint denominators. Ambiguity means a midpoint of Farey neighbors Assume a denominator bound $B$ gives two distinct best approximations $a/b$ and $c/d$ with $a/b \lt x \lt c/d$ and $b,d \le B$. There cannot be any reduced fraction $u/v$ with $v \le B$ strictly between them, because such a fraction would be closer to $x$ than at least one endpoint....

Detailed mathematical approach

Problem Summary

For a real number $x$ and a bound $B$ on denominators (the statement calls it $d$), a best approximation is a fraction $r/s$ in lowest terms with $s\le B$ such that every fraction closer to $x$ than $r/s$ has a denominator larger than $B$. A real number is called ambiguous if at least one bound gives two of them. The statement's example is $x=9/40$, which has the two best approximations $1/4$ and $1/5$ for the bound $6$. We must count the ambiguous numbers $x=p/q$ with $0\lt x\lt 1/100$ whose denominator, in lowest terms, satisfies $q\le Q=10^8$. ^[1]

The whole solution in one sentence: a number is ambiguous exactly when it is the midpoint of two fractions $a/b\lt c/d$ with $bc-ad=1$, its denominator is then $2bd$, and since every such pair appears exactly once in the Stern–Brocot tree, a walk through that tree that turns back as soon as $bd\gt Q/2$ or the interval lies at or above $1/100$ can count, one by one, the pairs whose midpoint lies below $1/100$.

Mathematical Approach

We first find out which numbers are ambiguous, and then how to visit exactly the pairs of fractions that produce them. Throughout, $M=\lfloor Q/2\rfloor$, and every fraction is written in lowest terms.

1. The example of the statement

A direct search over all bounds shows how the tie arises. For $x=9/40=0.225$, the nearest fraction with denominator at most $B$ is $0/1$ for $B\le2$, $1/3$ for $B=3$ and $1/4$ alone for $B=4$. For $B=5,6,7,8$ the fractions $1/5=0.2$ and $1/4=0.25$ are both at the distance $1/40$, and nothing with such a denominator comes closer. From $B=9$ on, the fraction $2/9$, at the distance $1/360$, is closer, and no tie occurs again. Three features of this example will turn out to be general: $9/40$ is the midpoint of $1/5$ and $1/4$; the two fractions satisfy $5\cdot1-1\cdot4=1$; and the tie lasts from the bound $\max(5,4)=5$ up to $5+4-1=8$, ending when $2/9=(1+1)/(5+4)$ becomes available. Note also that $40=2\cdot5\cdot4$.

2. Two best approximations force a midpoint

By the definition, the best approximations for a bound $B$ are exactly the fractions with denominator at most $B$ at the smallest distance $\delta$ from $x$. Two different ones cannot lie on the same side of $x$, so they are $x-\delta$ and $x+\delta$, and $x$ is their midpoint. No fraction with denominator at most $B$ lies strictly between them, because it would be closer than $\delta$. If $0\lt x\lt1$, the fractions $0/1$ and $1/1$ are always available, so $\delta\le\min(x,1-x)$ and both best approximations lie in $[0,1]$. For $B\ge q$, where $q$ is the denominator of $x$, the fraction $x$ itself is the only best approximation. Since the midpoint of two fractions is a fraction, only rational numbers can be ambiguous.

3. Farey sequences and neighbours

The fractions in $[0,1]$ with denominator at most $B$, listed in increasing order, form the Farey sequence $F_B$; for example, $F_4$ is $0/1,1/4,1/3,1/2,2/3,3/4,1/1$. ^[2] By Section 2, two tied best approximations are consecutive in $F_B$. We call two fractions $a/b\lt c/d$ neighbours if $bc-ad=1$, and we call $(a+c)/(b+d)$ their mediant. Two standard facts connect these notions. ^[2] ^[3]

(i) Between neighbours. Let $bc-ad=1$ and $a/b\lt p/q\lt c/d$. Then $bp-aq$ and $cq-dp$ are positive integers, and multiplying $q$ by $bc-ad=1$ gives

$$q=d\,(bp-aq)+b\,(cq-dp)\ge b+d.$$

Equality needs $bp-aq=cq-dp=1$, and the only solution of these two equations is $p=a+c$, $q=b+d$. So every fraction strictly between neighbours has a denominator of at least $b+d$, and the mediant is the only one with exactly $b+d$. The mediant is in lowest terms, because a common factor of $a+c$ and $b+d$ would divide $b(a+c)-a(b+d)=1$.

(ii) Consecutive Farey fractions are neighbours. This holds in $F_1$, which is $0/1,1/1$. Suppose it holds in $F_{B-1}$, and a new fraction $p/B$ falls into the gap between consecutive fractions $a/b\lt c/d$ of $F_{B-1}$. By (i), $B\ge b+d$. The mediant lies in the gap, so it is missing from $F_{B-1}$, which means $b+d\ge B$. Hence $B=b+d$, the new fraction is the mediant, and the two smaller gaps are neighbours again, since $b(a+c)-a(b+d)=1$ and $(b+d)c-(a+c)d=1$. In $F_4$, for instance, $3\cdot1-1\cdot2=1$ for $1/3\lt1/2$.

4. Which numbers are ambiguous

Theorem. A number $x$ with $0\lt x\lt1$ is ambiguous exactly when it is the midpoint of two neighbours $0\le a/b\lt c/d\le1$. The tie then occurs for the bounds $B$ with $\max(b,d)\le B\le b+d-1$, and the two best approximations are $a/b$ and $c/d$.

If $x$ is ambiguous, Section 2 gives two consecutive fractions of some $F_B$ with midpoint $x$, and by (ii) they are neighbours. Conversely, let $x$ be the midpoint of neighbours $a/b\lt c/d$, and let $\max(b,d)\le B\lt b+d$. A fraction closer to $x$ than these two lies strictly between them, so by (i) its denominator is at least $b+d\gt B$. Both endpoints are therefore best approximations for the bound $B$. For $B\ge b+d$ the mediant, which lies strictly inside the interval, is closer than both, and for $B\lt\max(b,d)$ one endpoint is not allowed; Section 5 shows that $x$ belongs to no other pair, so these are the only bounds with a tie. For $9/40$ this is the range $5\le B\le8$ found in Section 1.

5. The denominator of the midpoint

The midpoint of neighbours $a/b\lt c/d$ is

$$x=\frac12\Bigl(\frac ab+\frac cd\Bigr)=\frac{ad+bc}{2bd},$$

and this fraction is already in lowest terms. Since $bc=ad+1$, the numerator $ad+bc=2ad+1$ is odd. A prime dividing both $b$ and $ad+bc$ would divide $ad$, hence also $bc-ad=1$, which is impossible, and the same holds for $d$. So the denominator of $x$ is exactly $q=2bd$, and the condition $q\le Q$ means $bd\le M$. The half-length of the interval is $(bc-ad)/(2bd)=1/q$, so the pair can be read off from $x=p/q$ alone: $a/b=(p-1)/q$ and $c/d=(p+1)/q$. For $9/40$ these are $8/40=1/5$ and $10/40=1/4$. Different pairs therefore have different midpoints, and the ambiguous numbers below $1/100$ with $q\le Q$ correspond one to one to the pairs of neighbours in $[0,1]$ with $bd\le M$ whose midpoint lies below $1/100$.

Over the whole interval $(0,1)$ the pairs can even be counted. For coprime $b$ and $d$, exactly one pair of neighbours in $[0,1]$ has the left denominator $b$ and the right denominator $d$: the condition $bc-ad=1$ means $bc\equiv1\pmod d$, which fixes $c$ among $1,\dots,d$, and then $a=(bc-1)/d$. Hence the ambiguous numbers in $(0,1)$ with denominator $2m$ match the ordered factorisations $m=bd$ with $\gcd(b,d)=1$. Each prime power in $m$ must go entirely into $b$ or entirely into $d$, so there are $2^{\omega(m)}$ such factorisations, where $\omega(m)$ is the number of distinct primes dividing $m$. ^[5] For $m=6$ they are $1/12$, $5/12$, $7/12$ and $11/12$, the midpoints of $(0/1,1/6)$, $(1/3,1/2)$, $(1/2,2/3)$ and $(5/6,1/1)$; summing over $m\le6$ gives the $13$ ambiguous numbers with $q\le12$.

6. The Stern–Brocot tree lists every pair once

All pairs of neighbours in $[0,1]$ grow from the single interval $(0/1,1/1)$ by repeated splitting at the mediant. An interval $(a/b,c/d)$ has the left child $\bigl(a/b,(a+c)/(b+d)\bigr)$ and the right child $\bigl((a+c)/(b+d),c/d\bigr)$. Both children are neighbours again, since $b(a+c)-a(b+d)$ and $(b+d)c-(a+c)d$ both equal $bc-ad=1$. This is the Stern–Brocot tree, or more precisely its part between $0/1$ and $1/1$, which is called the Farey tree. ^[4] For example, the pair of Section 1 is reached along

$$\begin{aligned}\Bigl(\frac01,\frac11\Bigr)&\to\Bigl(\frac01,\frac12\Bigr)\to\Bigl(\frac01,\frac13\Bigr)\\&\to\Bigl(\frac01,\frac14\Bigr)\to\Bigl(\frac15,\frac14\Bigr).\end{aligned}$$

Every pair of neighbours in $[0,1]$ appears in this tree exactly once. Only the root has $b=d$, since then $b(c-a)=1$ forces $b=d=1$. In a child, the new endpoint is the mediant, and it has the larger denominator. So if $b\gt d$, the pair can only be the right child of $\bigl((a-c)/(b-d),c/d\bigr)$, which is again a pair of neighbours in $[0,1]$: its determinant is $(b-d)c-(a-c)d=1$, and its numerator is not negative because $d(a-c)=c(b-d)-1\ge0$. The case $b\lt d$ is symmetric, with the parent $\bigl(a/b,(c-a)/(d-b)\bigr)$. The sum of the two denominators drops at each step, so following parents leads back to the root along a single path, and each pair has exactly one place in the tree.

7. Two rules that cut off whole subtrees

We walk through the tree and count every pair with $bd\le M$ whose midpoint lies below $1/100$. Two rules let the walk skip whole subtrees without missing anything.

Product rule. The children of $(a/b,c/d)$ have the products $b(b+d)$ and $(b+d)d$, both larger than $bd$. Once $bd\gt M$, every pair below fails as well.
Threshold rule. Every interval below $(a/b,c/d)$ lies inside it. If $a/b\ge1/100$, every midpoint below is larger than $1/100$.

A large midpoint alone is no reason to stop: the root $(0/1,1/1)$ has the midpoint $1/2$, yet every answer lies below it. Conversely, nothing is lost. If a pair with $bd\le M$ has its midpoint below $1/100$, every pair on its path from the root has a smaller product and a left endpoint that is at most its own, which is below $1/100$, so the walk reaches it. In integers, the walk keeps an interval when $bd\le M$ and $100a\lt b$, and counts it when $100\,(ad+bc)\lt2bd$.

8. The left edge of the tree

The tree has a very long left edge: the pairs $(0/1,1/n)$ for $n=1,2,3,\dots$, with product $n$ and midpoint $1/(2n)$. The midpoint is below $1/100$ exactly when $n\ge51$, and the product rule ends the edge after $n=M$. So for every $Q\ge100$ the edge contributes the $M-50$ numbers $1/(2n)$ with $51\le n\le M$. For $Q=120$ these are the ten numbers $1/102,1/104,\dots,1/120$.

Every other pair lies below a right child $(1/(n+1),1/n)$ of an edge pair. For $n\le99$ this child starts at $1/(n+1)\ge1/100$, and the threshold rule removes it with everything below it. For $n\ge100$ the whole subtree lies inside $[1/(n+1),1/n]$, so every midpoint in it is below $1/100$, and every pair in it with $bd\le M$ counts. By (i), the fractions in that interval other than its endpoints have denominators of at least $2n+1$, so these pairs have $bd\ge n(n+1)\ge100\cdot101$ and denominators $q\ge20200$. Writing $R(Q)$ for the number of pairs below the edge that count, the answer for the bound $Q$ is

$$A(Q)=(M-50)+R(Q).$$

For $Q\lt20200$ we get $R(Q)=0$, so $A(Q)=\lfloor Q/2\rfloor-50$ for $100\le Q\lt20200$. The first number below the edge is $201/20200$, the midpoint of $1/101$ and $1/100$. For $Q=10^6$ the edge gives $499950$ numbers and $R=9813$, so $A=509763$. A subtree below the edge contains a pair with $bd\le M$ only if its root does, that is, if $n(n+1)\le M$; for $Q=10^8$ this leaves the subtrees with $100\le n\le7070$.

9. Worked example: $Q=40602$

Here $M=20301$, and the left edge gives $20301-50=20251$ numbers. Below the edge, the roots $(1/(n+1),1/n)$ have the product $n(n+1)$, which stays at most $20301$ for $100\le n\le141$, since $141\cdot142=20022$ and $142\cdot143=20306$. They give $42$ numbers $(2n+1)/(2n(n+1))$, from $201/20200$ to $283/40044$. One level deeper, only the children of $(1/101,1/100)$ pass the product rule: $(2/201,1/100)$ with $bd=20100$ and $(1/101,2/201)$ with $bd=20301$, whose midpoints are $401/40200$ and $403/40602$. For $n\ge101$ even the smaller child product $n(2n+1)$ is at least $101\cdot203=20503$, and the children of the two survivors have products of at least $30100$, so the product rule excludes everything deeper. In total $A(40602)=20251+42+2=20295$.

How the Code Works

The C++ program ^[9] stores an interval as a Node of four unsigned 64-bit integers a, b, c and d. The tests of Section 7 are three small functions: exceeds_product_bound checks $b\gt\lfloor M/d\rfloor$, which is equivalent to $bd\gt M$; left_endpoint_not_below_threshold checks $100a\ge b$; and midpoint_below_threshold checks $100\,(ad+bc)\lt2bd$. The last two multiply in unsigned __int128, a compiler extension, ^[7] although 64 bits would do: after the product test $a\le b$, $c\le d$ and $bd\le M$, so $100\,(ad+bc)\le200M=10^{10}$.

count_subtree_iterative is the walk of Section 7, with a std::vector as its stack. It takes an interval from the top, applies the product rule and the threshold rule, counts the interval if the midpoint test succeeds, and pushes the left child and then the right child, so that the right child is examined first. ^[6] This order keeps the stack short: besides the two children just pushed, it holds one waiting left child for each right turn on the path to the current interval. Each right turn on the way to an expanded interval adds the current right denominator, which is at least $100$, to the left one, so an expanded interval reached with $r$ right turns has $b\ge1+100r$ and $d\ge100$, and $bd\le M$ gives $r\le4999$ for $Q=10^8$. The stack therefore never exceeds $5001$ intervals, a size reached at the end of the chain $(k/(100k+1),1/100)$ below $(1/101,1/100)$. In the opposite order every left turn would leave a right child waiting, and the left edge alone would pile up $M+1=50000001$ intervals, about 1.6 GB of C++ nodes.

main first calls run_validation_checkpoints. This function counts by a single walk from the root for $Q=120$, $200000$ and $10^6$ and compares the results with the stored values $10$, $100967$ and $509763$. It also runs count_ambiguous_bruteforce for $Q=120$, which uses only the definition: after storing the fractions of every Farey sequence up to order $120$ in approximants_by_bound, it tries, for each of the $20$ fractions $1/101,\dots,1/120$, the bounds $1,2,\dots,120$ in turn, compares distances by exact cross-multiplication and stops at the first bound where two fractions share the smallest distance. A mismatch ends the program with an error. Below $10^6$, count_ambiguous_fast simply walks from the root; for the real bound it normally first builds a SeedFrontier: a helper processes the top of the tree breadth-first with a std::deque, with the same tests and counting, until the queue is long enough or empty. The intervals left in the queue are the roots of disjoint unexamined subtrees, which together with the processed top and the parts already cut off make up the whole tree, so adding the counts of count_subtree_iterative for these subtrees to the count of the top gives the result.

The Python and Java programs make the same walk from the root, with the same tests and the same push order, and print the count without any checks. Python keeps tuples on a list used as a stack, and Java keeps arrays long[] on an ArrayDeque with push and pop. ^[8] Python's integers are exact, and Java's long is ample by the bound $10^{10}$. Python's product test b > max_product // d if d > 0 else True is the C++ test, since every denominator is positive.

Complexity and Verification

Each interval taken from the stack costs a fixed number of integer operations and amortised constant stack work. Each expanded interval pushes two, and the root starts the stack, so $1+2E$ intervals are taken when $E$ are expanded. By Section 8 the expanded intervals are the $M$ edge pairs and the $R(Q)$ counted pairs below the edge, and only the $50$ edge pairs with $n\le50$ are expanded without being counted. Hence, for every $Q\ge100$, the walk takes exactly

$$1+2E=2A(Q)+101$$

intervals from the stack, as we confirmed for every $Q$ from $100$ to $3000$ and for $Q=10^6$, where it takes $1019627$. The left edge and its right children account for $2M+1$ of these intervals, so the time is at least of order $Q$. An upper bound comes from all pairs of neighbours in $[0,1]$ with $bd\le M$: their number is the sum of $2^{\omega(m)}$ over $m\le M$, which grows like $(6/\pi^2)\,M\ln M$, ^[5] so the time is at most of order $Q\log Q$. For $Q=10^8$ there are $578195943$ such pairs, and with the product rule alone the walk would expand all of them; the threshold rule cuts only $99$ subtrees, but large ones. The stack holds at most $5001$ intervals, and the C++ split also keeps its list of subtree roots. On a desktop computer the C++ program, checks included, takes less than half a second, the Java program about $0.6$ seconds and the Python program about $12$ seconds. Counting the left edge by the formula $M-50$ and walking only the subtrees with $n\ge100$ and $n(n+1)\le M$ would be much faster: for $Q=10^6$ these are the $607$ subtrees with $100\le n\le706$, and the walk takes $20233$ intervals instead of $1019627$.

Besides the packaged C++ checks described above, we checked the following for this article, outside the published solutions:

A search that uses only the definition, over every bound $B\lt q$, finds $577$ ambiguous numbers among all fractions in $(0,1)$ with $q\le300$. Each of them is the midpoint of the pair predicted in Section 5 and ties exactly for the bounds $\max(b,d)\le B\le b+d-1$, and $577$ is the sum of $2^{\omega(m)}$ over $m\le150$.
An enumeration without the tree runs over all coprime $(b,d)$ with $bd\le M$, finds $c$ as the inverse of $b$ modulo $d$ and tests the midpoint. It reproduces $10$, $100967$, $509763$ and $A(40602)=20295$, and for $Q=10^8$ it agrees with the three published programs, which agree with one another.
The formula $A(Q)=\lfloor Q/2\rfloor-50$ holds for every $Q$ from $102$ to $20199$.

Footnotes and References

The references below supply the background facts used above. The characterisation of ambiguous numbers, the denominator $2bd$, the split into the left edge and the subtrees below it and the exact number of examined intervals are derived explicitly in this article.

Project Euler 198 — Ambiguous Numbers. The official statement defines a best approximation to x for a denominator bound d, calls x ambiguous if some bound gives two best approximations, gives 9/40 with the best approximations 1/4 and 1/5 for the bound 6, and asks how many ambiguous x = p/q with 0 < x < 1/100 have q ≤ 10⁸.
Wikipedia — Farey sequence. The Farey sequence of order n lists the completely reduced fractions in [0, 1] with denominators at most n in increasing order. Neighbours a/b < c/d satisfy bc − ad = 1, fractions with bc − ad = 1 are neighbours in the sequence of order max(b, d), and the first term to appear between neighbours is their mediant, in the sequence of order b + d.
Wikipedia — Mediant (mathematics). The mediant of two fractions lies strictly between them, and when the two fractions satisfy the determinant relation, the mediant is the simplest fraction in the interval between them, in the sense of having the smallest denominator.
Wikipedia — Stern–Brocot tree. Introduced independently by Moritz Stern (1858) and the clockmaker Achille Brocot (1861), the tree is built from mediants and contains every positive rational number exactly once; its part between 0/1 and 1/1 is called the Farey tree.
OEIS A034444 — Number of unitary divisors of n. The unitary divisors of n are the divisors d with gcd(d, n/d) = 1; their number is 2^ω(n), where ω(n) counts the distinct primes dividing n, and the sum of these numbers up to N grows like (6/π²) N log N.
Wikipedia — Depth-first search. The article's non-recursive pseudocode keeps the vertices still to be visited on a stack instead of the queue of breadth-first search, and it visits the neighbours of each vertex in the opposite order to its recursive version.
GCC manual — 128-bit Integers. GCC supports the types __int128 and unsigned __int128 on targets that have an integer mode wide enough to hold 128 bits.
Python tutorial — Using lists as stacks explains that append and pop make a list a last-in, first-out stack. Java — ArrayDeque documents push and pop as insertion and removal at the front of the deque, which makes it a stack.
C++ — Euler198.cpp. The linked immutable C++ revision contains count_subtree_iterative, count_ambiguous_fast, count_ambiguous_bruteforce and the checks in run_validation_checkpoints. It is the implementation record for this article.

Problem 198 source code

C++

#include <atomic>
#include <cstdint>
#include <deque>
#include <iostream>
#include <numeric>
#include <thread>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr u64 kThresholdNumerator = 1;     // x < 1/100
constexpr u64 kThresholdDenominator = 100; // x < 1/100
constexpr u64 kTargetQ = 100'000'000ULL;

struct Node {
    u64 a;
    u64 b;
    u64 c;
    u64 d;
};

inline bool exceeds_product_bound(const Node& node, u64 max_product) {
    return node.b > max_product / node.d;
}

inline bool left_endpoint_not_below_threshold(const Node& node) {
    return static_cast<u128>(node.a) * kThresholdDenominator >=
           static_cast<u128>(kThresholdNumerator) * node.b;
}

inline bool midpoint_below_threshold(const Node& node) {
    // (a/b + c/d) / 2 < T  <=>  (a*d + b*c) * T_den < 2*T_num*b*d
    const u128 lhs = (static_cast<u128>(node.a) * node.d +
                      static_cast<u128>(node.b) * node.c) *
                     kThresholdDenominator;
    const u128 rhs = static_cast<u128>(2) * kThresholdNumerator * node.b * node.d;
    return lhs < rhs;
}

u64 count_subtree_iterative(const Node& seed, u64 max_product) {
    std::vector<Node> stack;
    stack.reserve(4096);
    stack.push_back(seed);

    u64 count = 0;

    while (!stack.empty()) {
        const Node node = stack.back();
        stack.pop_back();

        if (exceeds_product_bound(node, max_product)) {
            continue;
        }

        if (left_endpoint_not_below_threshold(node)) {
            continue;
        }

        if (midpoint_below_threshold(node)) {
            ++count;
        }

        const u64 mediant_n = node.a + node.c;
        const u64 mediant_d = node.b + node.d;

        // Push left first, then right; right is processed first with LIFO.
        // This keeps memory low on long left spines like (0/1, 1/n).
        stack.push_back(Node{node.a, node.b, mediant_n, mediant_d});
        stack.push_back(Node{mediant_n, mediant_d, node.c, node.d});
    }

    return count;
}

struct SeedFrontier {
    std::vector<Node> seeds;
    u64 expanded_node_count = 0;
};

SeedFrontier build_seed_frontier_with_prefix_count(u64 max_product,
                                                   std::size_t desired_seeds) {
    std::deque<Node> queue;
    queue.push_back(Node{0, 1, 1, 1});
    u64 expanded_node_count = 0;

    while (!queue.empty() && queue.size() < desired_seeds) {
        const Node node = queue.front();
        queue.pop_front();

        if (exceeds_product_bound(node, max_product) ||
            left_endpoint_not_below_threshold(node)) {
            continue;
        }

        if (midpoint_below_threshold(node)) {
            ++expanded_node_count;
        }

        const u64 mediant_n = node.a + node.c;
        const u64 mediant_d = node.b + node.d;

        queue.push_back(Node{node.a, node.b, mediant_n, mediant_d});
        queue.push_back(Node{mediant_n, mediant_d, node.c, node.d});
    }

    return SeedFrontier{
        std::vector<Node>(queue.begin(), queue.end()),
        expanded_node_count,
    };
}

u64 count_ambiguous_fast(u64 q_limit, bool allow_multithreading, unsigned requested_threads = 0) {
    if (q_limit < 2) {
        return 0;
    }

    const u64 max_product = q_limit / 2;

    unsigned threads = requested_threads;
    if (threads == 0) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0) {
            threads = 1;
        }
    }

    if (!allow_multithreading || threads <= 1 || q_limit < 1'000'000ULL) {
        return count_subtree_iterative(Node{0, 1, 1, 1}, max_product);
    }

    const std::size_t desired_seeds = static_cast<std::size_t>(threads) * 64;
    const SeedFrontier frontier =
        build_seed_frontier_with_prefix_count(max_product, desired_seeds);
    const std::vector<Node>& seeds = frontier.seeds;

    if (seeds.empty()) {
        return frontier.expanded_node_count;
    }

    if (seeds.size() == 1) {
        return frontier.expanded_node_count +
               count_subtree_iterative(seeds[0], max_product);
    }

    threads = static_cast<unsigned>(std::min<std::size_t>(threads, seeds.size()));

    std::atomic<std::size_t> next_index{0};
    std::vector<u64> partial(threads, 0);
    std::vector<std::thread> workers;
    workers.reserve(threads);

    for (unsigned tid = 0; tid < threads; ++tid) {
        workers.emplace_back([&, tid]() {
            u64 local = 0;
            while (true) {
                const std::size_t index = next_index.fetch_add(1, std::memory_order_relaxed);
                if (index >= seeds.size()) {
                    break;
                }
                local += count_subtree_iterative(seeds[index], max_product);
            }
            partial[tid] = local;
        });
    }

    for (std::thread& worker : workers) {
        worker.join();
    }

    u64 total = 0;
    for (u64 value : partial) {
        total += value;
    }
    return frontier.expanded_node_count + total;
}

u64 count_ambiguous_bruteforce(u64 q_limit) {
    struct Frac {
        u64 p;
        u64 q;
    };

    std::vector<Frac> targets;
    for (u64 q = 1; q <= q_limit; ++q) {
        for (u64 p = 1; p < q; ++p) {
            if (std::gcd(p, q) != 1) {
                continue;
            }
            if (p * kThresholdDenominator >= kThresholdNumerator * q) {
                continue;
            }
            targets.push_back(Frac{p, q});
        }
    }

    std::vector<Frac> approximants;
    std::vector<std::vector<Frac>> approximants_by_bound(q_limit + 1);

    for (u64 den = 1; den <= q_limit; ++den) {
        for (u64 num = 0; num <= den; ++num) {
            if (std::gcd(num, den) == 1) {
                approximants.push_back(Frac{num, den});
            }
        }
        approximants_by_bound[den] = approximants;
    }

    u64 ambiguous_count = 0;

    for (const Frac target : targets) {
        bool ambiguous = false;

        for (u64 bound = 1; bound <= q_limit && !ambiguous; ++bound) {
            u64 best_num = 0;
            u64 best_den = 1;
            bool has_best = false;
            int best_count = 0;

            for (const Frac cand : approximants_by_bound[bound]) {
                const std::int64_t diff = static_cast<std::int64_t>(target.p * cand.q) -
                                          static_cast<std::int64_t>(cand.p * target.q);
                const u64 num = static_cast<u64>(diff < 0 ? -diff : diff);
                const u64 den = target.q * cand.q;

                if (!has_best) {
                    has_best = true;
                    best_num = num;
                    best_den = den;
                    best_count = 1;
                    continue;
                }

                const u128 lhs = static_cast<u128>(num) * best_den;
                const u128 rhs = static_cast<u128>(best_num) * den;

                if (lhs < rhs) {
                    best_num = num;
                    best_den = den;
                    best_count = 1;
                } else if (lhs == rhs) {
                    ++best_count;
                }
            }

            if (best_count >= 2) {
                ambiguous = true;
            }
        }

        if (ambiguous) {
            ++ambiguous_count;
        }
    }

    return ambiguous_count;
}

bool run_validation_checkpoints() {
    struct Checkpoint {
        u64 q_limit;
        u64 expected;
    };

    const std::vector<Checkpoint> checkpoints = {
        {120, 10},
        {200'000, 100'967},
        {1'000'000, 509'763},
    };

    for (const Checkpoint cp : checkpoints) {
        const u64 got = count_ambiguous_fast(cp.q_limit, false, 1);
        if (got != cp.expected) {
            std::cerr << "Checkpoint failed for q <= " << cp.q_limit << ": got " << got
                      << ", expected " << cp.expected << "\n";
            return false;
        }
    }

    const u64 brute_q_limit = 120;
    const u64 brute_expected = checkpoints.front().expected;
    const u64 brute_got = count_ambiguous_bruteforce(brute_q_limit);

    if (brute_got != brute_expected) {
        std::cerr << "Brute-force checkpoint failed for q <= " << brute_q_limit << ": got "
                  << brute_got << ", expected " << brute_expected << "\n";
        return false;
    }

    return true;
}

} // namespace

int main() {
    if (!run_validation_checkpoints()) {
        return 1;
    }

    unsigned threads = std::thread::hardware_concurrency();
    if (threads == 0) {
        threads = 4;
    }

    const u64 answer = count_ambiguous_fast(kTargetQ, true, threads);
    std::cout << answer << '\n';
    return 0;
}

Python

def solve():
    T_NUM = 1
    T_DEN = 100
    TARGET_Q = 100_000_000

    def count_subtree(seed_a, seed_b, seed_c, seed_d, max_product):
        stack = [(seed_a, seed_b, seed_c, seed_d)]
        count = 0
        while stack:
            a, b, c, d = stack.pop()
            if b > max_product // d if d > 0 else True:
                continue
            if a * T_DEN >= T_NUM * b:
                continue
            ad_bc = a * d + b * c
            if ad_bc * T_DEN < 2 * T_NUM * b * d:
                count += 1
            mn = a + c
            md = b + d
            stack.append((a, b, mn, md))
            stack.append((mn, md, c, d))
        return count

    max_product = TARGET_Q // 2
    answer = count_subtree(0, 1, 1, 1, max_product)
    return str(answer)

if __name__ == '__main__':
    print(solve())

Java

import java.util.*;

public class Euler198 {
    public static void main(String[] args) {
        long qLimit = 100000000L;
        long maxProduct = qLimit / 2;
        long count = 0;
        Deque<long[]> stack = new ArrayDeque<>();
        stack.push(new long[] { 0, 1, 1, 1 });
        while (!stack.isEmpty()) {
            long[] f = stack.pop();
            long a = f[0], b = f[1], c = f[2], d = f[3];
            if (b > maxProduct / d)
                continue;
            if (a * 100 >= b)
                continue;
            // midpoint check: (a*d + b*c)*100 < 2*b*d
            if ((a * d + b * c) * 100 < 2 * b * d)
                count++;
            long mn = a + c, md = b + d;
            stack.push(new long[] { a, b, mn, md });
            stack.push(new long[] { mn, md, c, d });
        }
        System.out.println(count);
    }
}