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

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

Source: Windows Questions C++

LEAVE A COMMENT