All Euler problems
Project Euler

Ambiguous Numbers

A best rational approximation to a real number x for denominator bound d is a fraction r/s (fully reduced, s <= d) such that |x - r/s| < |x - p/q| for every p/q with q <= d and q ne s. A real numbe...

Source sync May 21, 2026
Problem #0198
Level Level 13
Solved By 1,378
Languages C++, Python
Answer 52374425
Length 490 words
modular_arithmeticnumber_theoryanalytic_math

Problem Statement

This archive keeps the full statement, math, and original media on the page.

A best approximation to a real number \(x\) for the denominator bound \(d\) is a rational number \(\frac r s\) (in reduced form) with \(s \le d\), so that any rational number \(\frac p q\) which is closer to \(x\) than \(\frac r s\) has \(q > d\).

Usually the best approximation to a real number is uniquely determined for all denominator bounds. However, there are some exceptions, e.g. \(\frac 9 {40}\) has the two best approximations \(\frac 1 4\) and \(\frac 1 5\) for the denominator bound \(6\). We shall call a real number \(x\) ambiguous, if there is at least one denominator bound for which \(x\) possesses two best approximations. Clearly, an ambiguous number is necessarily rational.

How many ambiguous numbers \(x=\frac p q, 0 < x < \frac 1 {100}\), are there whose denominator \(q\) does not exceed \(10^8\)?

Problem 198: Ambiguous Numbers

Mathematical Development

Theorem 1. (Midpoint Characterization.) A rational number xx is ambiguous if and only if xx is the midpoint of two Farey neighbours. That is, x=a/b+c/d2x = \frac{a/b + c/d}{2} where a/b<c/da/b < c/d are adjacent fractions in some Farey sequence FDF_D (meaning bc−ad=1bc - ad = 1).

Proof. (⇒\Rightarrow) If xx is ambiguous, there exists a bound DD such that two distinct fractions r1/s1r_1/s_1 and r2/s2r_2/s_2 (both with denominator ≤D\leq D) are equidistant from xx and closer than all others. By the properties of best approximations, these two fractions must be adjacent in FDF_D (no fraction with denominator ≤D\leq D lies between them), and xx is their midpoint.

(⇐\Leftarrow) If x=(a/b+c/d)/2x = (a/b + c/d)/2 with bc−ad=1bc - ad = 1, then for D=max⁡(b,d)−1D = \max(b, d) - 1 (or more precisely, for the largest DD such that no Farey fraction lies strictly between a/ba/b and c/dc/d), both a/ba/b and c/dc/d are equidistant from xx, making xx ambiguous. □\square

Theorem 2. (Lowest-Terms Form of Midpoints.) For Farey neighbours a/b<c/da/b < c/d with bc−ad=1bc - ad = 1, their midpoint is

x=ad+bc2bd=2ad+12bd,x = \frac{ad + bc}{2bd} = \frac{2ad + 1}{2bd},

and this fraction is already in lowest terms, i.e., gcd⁡(2ad+1,2bd)=1\gcd(2ad + 1, 2bd) = 1, so the denominator in lowest terms is q=2bdq = 2bd.

Proof. Since bc−ad=1bc - ad = 1, we have ad+bc=ad+(ad+1)=2ad+1ad + bc = ad + (ad + 1) = 2ad + 1, which is odd. Thus gcd⁡(2ad+1,2)=1\gcd(2ad + 1, 2) = 1.

For gcd⁡(2ad+1,b)\gcd(2ad + 1, b): since ad≡−1(modb)ad \equiv -1 \pmod{b} (from bc−ad=1bc - ad = 1 and bc≡0(modb)bc \equiv 0 \pmod{b}), we get 2ad+1≡−2+1=−1(modb)2ad + 1 \equiv -2 + 1 = -1 \pmod{b}, so gcd⁡(2ad+1,b)=gcd⁡(−1,b)=1\gcd(2ad + 1, b) = \gcd(-1, b) = 1.

For gcd⁡(2ad+1,d)\gcd(2ad + 1, d): since bc≡1(modd)bc \equiv 1 \pmod{d}, we get ad+bc≡0+1=1(modd)ad + bc \equiv 0 + 1 = 1 \pmod{d}, so gcd⁡(2ad+1,d)=gcd⁡(1,d)=1\gcd(2ad + 1, d) = \gcd(1, d) = 1.

Therefore gcd⁡(2ad+1,bd)=1\gcd(2ad + 1, bd) = 1, and since gcd⁡(2ad+1,2)=1\gcd(2ad + 1, 2) = 1, we have gcd⁡(2ad+1,2bd)=1\gcd(2ad + 1, 2bd) = 1. □\square

Lemma 1. (Counting Reduction.) The condition q≤N=108q \leq N = 10^8 becomes bd≤N/2bd \leq N/2. For each coprime pair (b,d)(b, d) with b≥1b \geq 1, d≥1d \geq 1, bd≤N/2bd \leq N/2, there is a unique a∈{0,1,…,b−1}a \in \{0, 1, \ldots, b-1\} satisfying ad≡−1(modb)ad \equiv -1 \pmod{b} (namely a=b−d−1 mod ba = b - d^{-1} \bmod b). This gives a midpoint x∈[0,1)x \in [0, 1).

Proof. Since gcd⁡(b,d)=1\gcd(b, d) = 1, the inverse d−1 mod bd^{-1} \bmod b exists and is unique. The value a=(−d−1) mod ba = (-d^{-1}) \bmod b is the unique solution in {0,…,b−1}\{0, \ldots, b-1\}. □\square

Lemma 2. (Range Constraint.) The midpoint x=(2ad+1)/(2bd)x = (2ad+1)/(2bd) lies in (0,1/100)(0, 1/100) if and only if

0<a<b100−12d.0 < a < \frac{b}{100} - \frac{1}{2d}.

For b<100b < 100, the only possibility is a=0a = 0, which gives x=1/(2bd)<1/100x = 1/(2bd) < 1/100 iff bd>50bd > 50.

Proof. x>0x > 0 requires a≥0a \geq 0 (with a=0a = 0 giving x=1/(2bd)x = 1/(2bd)). x<1/100x < 1/100 requires 2ad+1<2bd/1002ad + 1 < 2bd/100, i.e., a<b/100−1/(2d)a < b/100 - 1/(2d). For b<100b < 100, the bound b/100<1b/100 < 1 forces a=0a = 0. □\square

Editorial

An ambiguous rational appears exactly at the midpoint of two Farey neighbours, so the search space is much smaller than all fractions with denominator up to 10^8. The midpoint formula shows that the reduced denominator is always 2bd, which immediately turns the denominator cap into the product bound bd <= 5 * 10^7. After that, the remaining test is whether the unique residue class for a places the midpoint below 1/100.

The implementation splits the count at B = floor(sqrt(N/2)). For small b, there are many valid d values, so it is efficient to fix b, enumerate the feasible a values, and count matching d values in one arithmetic progression determined by a modular inverse. For large b, the product bound forces d to be small, so the loops are reversed: the code fixes d, works with the near-terminal values of the auxiliary parameter in dr-kb=1, and counts compatible b values in residue classes modulo d. Those two passes cover the whole search space without ever scanning every denominator up to 5 * 10^7.

Pseudocode

Set HALF = 10^8 / 2 and B = floor(sqrt(HALF)).
Initialize the answer with the special case b = 1, namely all d > 50.

For each small denominator b from 100 up to B:
    compute the largest possible d from HALF / b;
    iterate over the feasible values of a below b / 100;
    skip a when gcd(a, b) is not 1;
    use the modular inverse of a modulo b to determine the residue class for d;
    count how many d in that arithmetic progression satisfy both the range bound
    and the midpoint inequality.

For each small d from 1 up to B:
    compute the largest compatible b from HALF / d;
    iterate over the admissible near-terminal values of the auxiliary parameter m;
    skip m when gcd(m, d) is not 1;
    use the modular inverse of m modulo d to recover the residue class for b;
    count the valid b values in that progression.

Return the combined total.

Complexity Analysis

  • Time: Cases 2 and 3 each involve O(N/2)O(\sqrt{N/2}) outer iterations, with inner work proportional to O(N)O(\sqrt{N}) total. Overall: O(N)O(N) arithmetic operations.
  • Space: O(1)O(1) beyond loop variables and small lookup tables for modular inverses.

Answer

52374425\boxed{52374425}

Code

Each problem page includes the exact C++ and Python source files from the local archive.

C++ project_euler/problem_198/solution.cpp
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

// Problem 198: Ambiguous Numbers
//
// x = p/q is ambiguous iff it equals (ad+bc)/(2bd) for Farey neighbors a/b, c/d
// with bc-ad=1. The denominator is always 2bd in lowest terms.
// So we need bd <= N/2 and the midpoint in (0, 1/100).
//
// Split into: b=1 case, small b (100..B), and large b (B+1..N/2) via d-iteration.

ll mod_inverse(ll a, ll m) {
    // Extended Euclidean
    ll g = __gcd(a, m);
    if (g != 1) return -1;
    // Using pow with modular inverse for prime... but m may not be prime.
    // Use extended gcd.
    ll x = 1, y = 0, x1 = 0, y1 = 1, a0 = a, m0 = m;
    while (m0 != 0) {
        ll q = a0 / m0;
        ll tmp;
        tmp = x - q * x1; x = x1; x1 = tmp;
        tmp = y - q * y1; y = y1; y1 = tmp;
        tmp = a0 - q * m0; a0 = m0; m0 = tmp;
    }
    return ((x % m) + m) % m;
}

int main() {
    ios_base::sync_with_stdio(false);

    const ll N = 100000000LL;
    const ll HALF = N / 2;
    const ll B = (ll)sqrt((double)HALF);

    ll count = 0;

    // b = 1: d > 50, d <= HALF
    count += HALF - 50;

    // Small b: 100 to B
    for (ll b = 100; b <= B; b++) {
        ll d_max = HALF / b;
        if (d_max < 1) break;

        ll a_upper;
        if (b % 100 == 0)
            a_upper = b / 100 - 1;
        else
            a_upper = b / 100;

        for (ll a = 1; a <= a_upper; a++) {
            if (__gcd(a, b) != 1) continue;

            ll a_inv = mod_inverse(a, b);
            ll d_res = (b - a_inv) % b;  // (-a_inv) mod b
            if (d_res <= 0) d_res += b;

            ll gap = b - 100 * a;
            if (gap <= 0) continue;

            ll d_min = 50 / gap + 1; // strict >

            ll first_d;
            if (d_res >= d_min)
                first_d = d_res;
            else {
                ll k = (d_min - d_res + b - 1) / b;
                first_d = d_res + k * b;
            }

            if (first_d > d_max) continue;
            count += (d_max - first_d) / b + 1;
        }
    }

    // Large b: iterate d from 1 to B, then m values
    for (ll d = 1; d <= B; d++) {
        ll b_max = HALF / d;
        if (b_max <= B) continue;

        ll m_min;
        if ((99 * d) % 100 == 0)
            m_min = 99 * d / 100;
        else
            m_min = 99 * d / 100 + 1;

        ll m_max = d - 1;

        for (ll m = m_min; m <= m_max; m++) {
            if (__gcd(m, d) != 1) continue;

            ll m_inv = mod_inverse(m, d);
            ll b_res = (d - m_inv) % d; // (-m^{-1}) mod d
            if (b_res <= 0) b_res += d;

            ll first_b;
            if (b_res >= B + 1)
                first_b = b_res;
            else {
                ll k = (B + 1 - b_res + d - 1) / d;
                first_b = b_res + k * d;
            }

            if (first_b > b_max) continue;
            count += (b_max - first_b) / d + 1;
        }
    }

    cout << count << endl;
    return 0;
}