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, int v)
{
auto w = lldiv(a, b);

// u * a' - u * b' == a
// v * a' - v * b' == b

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

auto m = u - w.quot * v;
auto n = u - w.quot * v;
u = v;
u = v;

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 = int(m);
v = 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 = int(d);
v = 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;
auto mm = v;

// 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;

approx_brute(dd, &n, &m);
std::cout << n << " / " << m << "  dur: " << (after - before).count() << std::endl;
approx(dd, &n, &m);
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++