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
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.