I'm looking for a method to convert the exact value of a floating-point number to a rational quotient of two integers, i.e. a / b, where b is not larger than a specified maximum denominator b_max. If satisfying the condition b <= b_max is impossible, then the result falls back to the best approximation which still satisfies the condition.
Hold on. There are a lot of questions/answers here about the best rational approximation of a truncated real number which is represented as a floating-point number. However I'm interested in the exact value of a floating-point number, which is itself a rational number with a different representation. More specifically, the mathematical set of floating-point numbers is a subset of rational numbers. In case of IEEE 754 binary floating-point standard it is a subset of dyadic rationals. Anyway, any floating-point number can be converted to a rational quotient of two finite precision integers as a / b.
So, for example assuming IEEE 754 single-precision binary floating-point format, the rational equivalent of float f = 1.0f / 3.0f is not 1 / 3, but 11184811 / 33554432. This is the exact value of f, which is a number from the mathematical set of IEEE 754 single-precision binary floating-point numbers.
Based on my experience, traversing (by binary search of) the Stern-Brocot tree is not useful here, since that is more suitable for approximating the value of a floating-point number, when it is interpreted as a truncated real instead of an exact rational.
Possibly, continued fractions are the way to go.
The another problem here is integer overflow. Think about that we want to represent the rational as the quotient of two int32_t, where the maximum denominator b_max = INT32_MAX. We cannot rely on a stopping criterion like b > b_max. So the algorithm must never overflow, or it must detect overflow.
What I found so far is an algorithm from Rosetta Code, which is based on continued fractions, but its source mentions it is "still not quite complete". Some basic tests gave good results, but I cannot confirm its overall correctness and I think it can easily overflow.
// https://rosettacode.org/wiki/Convert_decimal_number_to_rational#C
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>
/* f : number to convert.
* num, denom: returned parts of the rational.
* md: max denominator value. Note that machine floating point number
* has a finite resolution (10e-16 ish for 64 bit double), so specifying
* a "best match with minimal error" is often wrong, because one can
* always just retrieve the significand and return that divided by
* 2**52, which is in a sense accurate, but generally not very useful:
* 1.0/7.0 would be "2573485501354569/18014398509481984", for example.
*/
void rat_approx(double f, int64_t md, int64_t *num, int64_t *denom)
{
/* a: continued fraction coefficients. */
int64_t a, h[3] = { 0, 1, 0 }, k[3] = { 1, 0, 0 };
int64_t x, d, n = 1;
int i, neg = 0;
if (md <= 1) { *denom = 1; *num = (int64_t) f; return; }
if (f < 0) { neg = 1; f = -f; }
while (f != floor(f)) { n <<= 1; f *= 2; }
d = f;
/* continued fraction and check denominator each step */
for (i = 0; i < 64; i++) {
a = n ? d / n : 0;
if (i && !a) break;
x = d; d = n; n = x % n;
x = a;
if (k[1] * a + k[0] >= md) {
x = (md - k[0]) / k[1];
if (x * 2 >= a || k[1] >= md)
i = 65;
else
break;
}
h[2] = x * h[1] + h[0]; h[0] = h[1]; h[1] = h[2];
k[2] = x * k[1] + k[0]; k[0] = k[1]; k[1] = k[2];
}
*denom = k[1];
*num = neg ? -h[1] : h[1];
}
f, which is a number from the mathematical set of IEEE 754 single-precision binary floating-point numbers." The mathematical set of floating-point numbers is a subset of rational numbers. All IEEE 754 floating-point numbers can be represented as a quotient of two integers. There is no need for infinite bits. I already wrote a Wolfram Mathematica function to do that. Unfortunately that is not portable to C++. – DevastationFraction.limit_denominator, which does exactly this. Of course, that doesn't help with the overflow problem in C, but the algorithm may be useful. Specifically, in Python, for a (Python)floatx,Fraction(x).limit_denominator(b_max)gives what you want. – Schnabelboost::multiprecision. – Oopskrationals with different denominators, i.e. you have to perform a lot of multiplication and addition (not to mention applying the exponent at the end), which leads to overflow very very soon. Of course you could simplify in each step using the greatest common divisor, but that would make the algorithm very slow. This approach also cannot provide the closest approx if the rational cannot be represented asb <= b_max. – Devastationa / bwill be arbitrary-precision. If they are not, then I am interested in the best approximation which still satisfiesb <= b_max, whereb_maxis the largest value of that integer type. This is why an iterative refinement algorithm is desirable. – DevastationFLT_RADIX, which defines the radix used for the floating-point formats. Multiplying or dividing byFLT_RADIXshould be exact (and note thescalbnfunction), so, given some floating-point valuexthat is not an integer, one can find the desired numerator by multiplying byFLT_RADIXuntil the result is an integer, and the denominator isFLT_RADIXto the power of the number of multiplications. – SetterFLT_RADIX. – DevastationFLT_RADIXuntil a number is an integer will not overflow in any normal format. Theoretically, one could define a floating-point format such that the exponent is limited to a weird range like -20 to -10 instead of a more normal -126 to 127, but, in practice, exponent ranges allow the significand to be scaled to an integer. – Setter