Rational approximation of double using int numerator and denominator in C++

Question

A real world third party API takes a parameter of type fraction which is a struct of an int numerator and denominator. The value that I need to pass is known to me as a decimal string that is converted to a double.

The range of possible values are, let's say 10K to 300M but if there is a fraction part after the decimal point, it's significant.

I have here code for two approximation approaches, one uses the extended euclidean algorithm while the other is brute-force. Both methods find a rational approximation using int types for a given double.

The brute-force is of course the more accurate of the two and is actually faster when the converted numbers are large. My questions is, can I say anything clever about the quality of the approximation using the euclidean algorithm. More formally, can I put a bound on the approximation using the euclidean algorithm vs. the approximation of the brute-force algorithm (which I believe to be optimal).

An example for a bound:
If the error of the optimal approximation is r, then the euclidean algorithm approximation would produce an error that is less than 2*r.
(I'm not claiming this is the bound and I certainly can't prove it, it's just an example for what a good bound may look like).

Here's the code an a test program:

#include <iostream>
#include <iomanip>
#include <cmath>
#include <limits>
#include <chrono>
#include <random>

// extended euclidian algorithm
// finds the coefficients that produce the gcd
// in u, we store m,n the coefficients that produce m*a - n*b == gcd.
// in v, we store m,n the coefficients that produce m*a - n*b == 0.
// breaks early if the coefficients become larger than INT_MAX
int gcd_e(uint64_t a, int b, int u[2], int v[2])
{
    auto w = lldiv(a, b);

    // u[0] * a' - u[1] * b' == a
    // v[0] * a' - v[1] * b' == b

    // a - w.quot * b == w.rem
    // (u[0] * a' - u[1] * b') - w.quot * (v[0] * a' - v[1] * b') == w.rem
    // (u[0] - w.quot * v[0]) * a' - u[1] * b' + w.quot * v[1] * b' == w.rem
    // (u[0] - w.quot * v[0]) * a' + (w.quot * v[1] - u[1]) * b' == w.rem
    // (u[0] - w.quot * v[0]) * a' - (u[1] - w.quot * v[1]) * b' == w.rem


    auto m = u[0] - w.quot * v[0];
    auto n = u[1] - w.quot * v[1];
    u[0] = v[0];
    u[1] = v[1];

    constexpr auto L = std::numeric_limits<int>::max();
    if (m > L || n > L)
        throw 0;  // break early
    if (m < -L || n < -L)
        throw 0;  // break early

    v[0] = int(m);
    v[1] = int(n);

    if (w.rem == 0)
        return b;

    return gcd_e(b, int(w.rem), u, v);
}


inline double helper_pre(double d, bool* negative, bool* inverse)
{
    bool v = (d < 0);
    *negative = v;
    if (v)
        d = -d;

    v = (d < 1);
    *inverse = v;
    if (v)
        d = 1 / d;

    return d;
}

inline void helper_post(int* m, int* n, bool negative, bool inverse)
{
    if (inverse)
        std::swap(*n, *m);

    if (negative)
        *n = -(*n);
}

// gets a rational approximation for double d
// numerator is stored in n
// denominator is stored in m
void approx(double d, int* n, int *m)
{
    int u[] = { 1, 0 };  // 1*a - 0*b == a
    int v[] = { 0, -1 }; // 0*a - (-1)*b == b

    bool negative, inverse;
    d = helper_pre(d, &negative, &inverse);

    constexpr int q = 1 << 30;
    auto round_d = std::round(d);
    if (d == round_d)
    {
        // nothing to do, it's an integer.
        v[1] = int(d);
        v[0] = 1;
    }
    else try
    {
        uint64_t k = uint64_t(std::round(d*q));
        gcd_e(k, q, u, v);
    }
    catch (...)
    {
        // OK if we got here.
        // int limits
    }

    // get the approximate numerator and denominator
    auto nn = v[1];
    auto mm = v[0];

    // make them positive
    if (mm < 0)
    {
        mm = -mm;
        nn = -nn;
    }

    helper_post(&mm, &nn, negative, inverse);

    *m = mm;
    *n = nn;
}


// helper to test a denominator
// returns the magnitude of the error
double helper_rattest(double x, int tryDenom, int* numerator)
{
    double r = x * tryDenom;
    double rr = std::round(r);
    auto num = int(rr);
    auto err = std::abs(r - rr) / tryDenom;
    *numerator = num;
    return err;
}

// helper to reduce the rational number
int gcd(int a, int b)
{
    auto c = a % b;
    if (c == 0)
        return b;
    return gcd(b, int(c));
}

// gets a rational approximation for double d
// numerator is stored in n
// denominator is stored in m
// uses brute force by scanning denominator range
void approx_brute(double d, int* n, int* m)
{
    bool negative, inverse;
    d = helper_pre(d, &negative, &inverse);

    int upto = int(std::numeric_limits<int>::max() / d);
    int bestNumerator;
    int bestDenominator = 1;
    auto bestErr = helper_rattest(d, 1, &bestNumerator);
    for (int kk = 2; kk < upto; ++kk)
    {
        int n;
        auto e = helper_rattest(d, kk, &n);
        if (e < bestErr)
        {
            bestErr = e;
            bestNumerator = n;
            bestDenominator = kk;
        }
        if (bestErr == 0)
            break;
    }

    // reduce, just in case
    auto g = gcd(bestNumerator, bestDenominator);
    bestNumerator /= g;
    bestDenominator /= g;

    helper_post(&bestDenominator, &bestNumerator, negative, inverse);

    *n = bestNumerator;
    *m = bestDenominator;
}

int main()
{
    int n, m;

    auto re = std::default_random_engine();
    std::random_device rd;
    re.seed(rd());

    for (auto& u : {
        std::uniform_real_distribution<double>(10000,    15000),
        std::uniform_real_distribution<double>(100000,   150000),
        std::uniform_real_distribution<double>(200000,   250000),
        std::uniform_real_distribution<double>(400000,   450000),
        std::uniform_real_distribution<double>(800000,   850000),
        std::uniform_real_distribution<double>(1000000,  1500000),
        std::uniform_real_distribution<double>(2000000,  2500000),
        std::uniform_real_distribution<double>(4000000,  4500000),
        std::uniform_real_distribution<double>(8000000,  8500000),
        std::uniform_real_distribution<double>(10000000, 15000000)
        })
    {
        auto dd = u(re);
        std::cout << "approx: " << std::setprecision(14) << dd << std::endl;

        auto before = std::chrono::steady_clock::now();
        approx_brute(dd, &n, &m);
        auto after = std::chrono::steady_clock::now();
        std::cout << n << " / " << m << "  dur: " << (after - before).count() << std::endl;
        before = std::chrono::steady_clock::now();
        approx(dd, &n, &m);
        after = std::chrono::steady_clock::now();
        std::cout << n << " / " << m << "  dur: " << (after - before).count()
            << std::endl
            << std::endl;
    }
}

Here's some sample output:

approx: 13581.807792679
374722077 / 27590  dur: 3131300
374722077 / 27590  dur: 15000

approx: 103190.31976517
263651267 / 2555  dur: 418700
263651267 / 2555  dur: 6300

approx: 223753.78683426
1726707973 / 7717  dur: 190100
1726707973 / 7717  dur: 5800

approx: 416934.79214075
1941665327 / 4657  dur: 102100
403175944 / 967  dur: 5700

approx: 824300.61241502
1088901109 / 1321  dur: 51900
1088901109 / 1321  dur: 5900

approx: 1077460.29557
1483662827 / 1377  dur: 39600
1483662827 / 1377  dur: 5600

approx: 2414781.364653
1079407270 / 447  dur: 17900
1079407270 / 447  dur: 7300

approx: 4189869.294816
1776504581 / 424  dur: 10600
1051657193 / 251  dur: 9900

approx: 8330270.2432111
308219999 / 37  dur: 5400
308219999 / 37  dur: 10300

approx: 11809264.006453
1830435921 / 155  dur: 4000
1830435921 / 155  dur: 10500

Answer 1

Thanks to all who commented and drew my attention to the concept of continued fractions. According to this paper by (William F. Hammond) There is equivalence between the euclidean algorithm and the continued fractions method.

The sub-optimal results are due to the fact that the numerator is constrained as well as the denominator so if the non brute force algorithm only produces "convergents" it means that it neglects the range of denominators between the first convergent to violate the constraints and the one just before it.

The denominators after the returned convergent and the one that follows may approximate close to the latter convergent and the difference between subsequent convergents can be shown to be:

So I suppose this would be the bound on the difference between the brute-force and the euclidean algorithm. The ratio of the error between them can be practically anything.
(can find examples of error ratios of more than 100 easily)

I hope I read everything correctly. I'm no authority on this.