Convert a float to a rational number that is guaranteed to convert back to the original float

Question

I am looking for an algorithm to convert a float to a rational number, such that the rational number is guaranteed to evaluate back to the original float, and the denominator is minimized.

A naive algorithm can just return the actual value of the float as X / 2^N, but that 2^N tends to be pretty large for anything that is not a finite binary fraction. For example, the number 0.1, when stored in a double-precision float, is actually approximated by ³⁶⁰²⁸⁷⁹⁷⁰¹⁸⁹⁶³⁹⁷⁄₃₆₀₂₈₇₉₇₀₁₈₉₆₃₉₆₈ (the denominator being 2⁵⁵). However, converting 0.1 to ¹⁄₁₀ is obviously better, and ¹⁄₁₀ will evaluate to ³⁶⁰²⁸⁷⁹⁷⁰¹⁸⁹⁶³⁹⁷⁄₃₆₀₂₈₇₉₇₀₁₈₉₆₃₉₆₈ under floating point arithmetic.

A related problem is printing floats in decimal with the least amount of digits (this paper describes some techniques), which can be considered a specialized version of this problem with an additional constraint that the denominator must be a power of 10.

There is an existing questions, and likely more, but they don't have the constraint that the converted rational number must evaluate to the original float.

Answer 1

Let's start with a definition that pins down exactly which fraction we're looking for in any particular case:

Definition. Say that a fraction a/b is simpler than another fraction c/d (with both fractions written in lowest terms, with positive denominators) if b <= d, abs(a) <= abs(c), and at least one of those two inequalities is strict.

So for example 5/7 is simpler than 6/7, and 5/7 is simpler than 5/8, but neither 2/5 nor 3/4 is simpler than the other. (We don't have a total ordering here.)

Then with this definition, there's a not-immediately-obvious theorem, that guarantees that the fraction we're looking for always exists:

Theorem. Given a subinterval J of the real numbers that contains at least one fraction, J contains a unique simplest fraction. In other words, there's a unique fraction f such that

• f is in J and,
• for any other fraction g in J, f is simpler than g.

In particular, the simplest fraction in an interval will always have smallest possible denominator, as required in the question. The "contains at least one fraction" condition in the theorem is there to exclude degenerate cases like the closed interval [√2, √2], which doesn't contain any fractions at all.

Our job is to write a function that takes a finite floating-point input x and returns the simplest fraction n/d for which x is the closest float to n/d, in the target floating-point format. Assuming a reasonably sane floating-point format and rounding mode, the set of real numbers that round to x will form a nonempty subinterval of the real line, with rational endpoints. So our problem breaks naturally into two subproblems:

• Problem 1. Given a float x in the target floating-point format, describe the interval of all values that round to x under the rules for that floating-point format. This involves both identifying the endpoints of that interval and establishing whether the interval is open, closed, or half-open.

• Problem 2. Given a nonempty subinterval J of the real line with rational endpoints, compute the simplest fraction in that subinterval.

The second problem is more interesting and less dependent on platform and language details; let's tackle that one first.

Finding the simplest fraction in an interval

Assuming an IEEE 754 floating-point format and the default round-ties-to-even rounding mode, the interval rounding to a given float will be either open or closed; with other rounding modes or formats, it could potentially be half open (open at one end, closed at the other). So for this section, we only look at open and closed intervals, but adapting to half-open intervals isn't hard.

Suppose that J is a nonempty subinterval of the real line with rational endpoints. For simplicity, let's assume that J is a subinterval of the positive real line. If it's not, then either it contains 0 — in which case 0/1 is the simplest fraction in J — or it's a subinterval of the negative real line and we can negate, find the simplest fraction, and negate back.

Then the following gives a simple recursive algorithm for finding the simplest fraction in J:

• If J contains 1, then 1/1 is the simplest fraction in J
• Otherwise, if J is a subinterval of (0, 1), then the simplest fraction in J is 1/f, where f is the simplest fraction in 1/J. (This is immediate from the definition of 'simplest'.)
• Otherwise, J must be a subinterval of (1, ∞), and the simplest fraction in J is q + f, where q is the largest integer such that J - q is still within the positive reals, and f is the simplest fraction in J - q.

For a sketch of proof of the last statement: if a / b is the simplest fraction in J and c / d is the simplest fraction in J - q, then a / b is simpler than or equal to (c + qd) / d, and c / d is simpler than or equal to (a - qb) / b. So b <= d, a <= c + qd, d <= b and c <= a - qb, and it follows that b = d and a = c + qd, so c / d = a / b - q.

In Python-like pseudocode:

def simplest_in_interval(J):
    # Simplest fraction in a subinterval J of the positive reals
    if J < 1:
        return 1 / simplest_in_interval(1/J)
    else if 1 < J:
        q = largest_integer_to_the_left_of(J)
        return q + simplest_in_interval(J - q)
    else:
        # If we get here then J contains 1.
        return 1

To see that the algorithm must always terminate and can't enter an infinite loop, note that every inversion step must be followed by a J - q step, and every J - q step reduces the numerators of the left and right endpoints of the interval. Concretely, if the endpoints of the interval are a/b and c/d, the sum abs(a) + abs(c) + b + d is a positive integer that steadily decreases as the algorithm progresses.

To translate the above to real Python code, we have to deal with some details. First, let's assume for now that J is a closed interval; we'll adapt to open intervals below.

We'll represent our interval by its endpoints left and right, both of which are positive fraction.Fraction instances. Then the following Python code implements the above algorithm.

from fractions import Fraction
from math import ceil

def simplest_in_closed_interval(left, right):
    """
    Simplest fraction in [left, right], assuming 0 < left <= right < ∞.
    """
    if right < 1:  # J ⊂ (0, 1)
        return 1 / simplest_in_closed_interval(1 / right, 1 / left)
    elif 1 < left:  # J ⊂ (1, ∞):
        q = ceil(left) - 1  # largest q satisfying q < left
        return q + simplest_in_closed_interval(left - q, right - q)
    else:  #  left <= 1 <= right, so 1 ∈ J
        return Fraction(1)

Here's an example run:

>>> simplest_in_closed_interval(Fraction("3.14"), Fraction("3.15"))
Fraction(22, 7)

In principle, the code for open intervals is equally simple, but in practice there's a complication: we may need to deal with infinite intervals. For example, if our original interval is J = (2, 5/2), then the first step shifts that interval by 2 to get (0, 1/2), and then that interval is inverted to give (2, ∞).

So for open intervals, we'll continue to represent our interval by a pair (left, right) of its endpoints, but now right is either a fractions.Fraction instance or a special constant INFINITY. And instead of simply being able to use 1 / left to take the reciprocal of the left endpoint, we'll need an auxiliary function that can compute a reciprocal of something that's either a fraction or INFINITY, and another auxiliary function for subtraction, ensuring that INFINITY - q gives INFINITY. Here are those auxiliary functions:

#: Constant used to represent an unbounded interval
INFINITY = "infinity"

def reciprocal(f):
    """ 1 / f, for f either a nonnegative fraction or ∞ """
    if f == INFINITY:
        return 0
    elif f == 0:
        return INFINITY
    else:
        return 1 / f


def shift(f, q):
    """ f - q, for f either a nonnegative fraction or ∞ """
    if f == INFINITY:
        return INFINITY
    else:
        return f - q

And here's the main function. Note the changes to the inequalities in the if and elif conditions, and the fact that we now want to use floor(left) instead of ceil(left) - 1 to find the largest integer q lying to the left of the interval:

from fractions import Fraction
from math import floor

def simplest_in_open_interval(left, right):
    """
    Simplest fraction in (left, right), assuming 0 <= left < right <= ∞.
    """
    if 1 <= left:  # J ⊆ (1, ∞):
        q = floor(left)
        return q + simplest_in_open_interval(shift(left, q), shift(right, q))
    elif right != INFINITY and right <= 1:  # J ⊆ (0, 1)
        return 1 / simplest_in_open_interval(reciprocal(right), reciprocal(left))
    else:  #  left < 1 < right, so 1 ∈ J
        return Fraction(1)

The above code is optimised for clarity, not efficiency: it's reasonably efficient in terms of big-Oh complexity, but not in terms of the implementation details. I leave it to the reader to convert to something more efficient. The first step would be to work with integer numerators and denominators throughout instead of fractions.Fraction instances. If you're interested in what that looks like, take a look at the implementation in my simplefractions package on PyPI.

Finding the interval that rounds to a given float

Now that we can find the simplest fraction in a given interval, we need to solve the other half of the problem: finding the interval that rounds to a given floating-point number. The details of doing this will depend much more heavily on the language, the floating-point format in use, and even things like the rounding mode being used.

Here we outline one way to do this in Python, assuming IEEE 754 binary64 floating-point format with the default round-ties-to-even rounding mode.

For simplicity, assume that our input float x is positive (and finite).

Python >= 3.9 provides a function math.nextafter that allows us to retrieve the next floats up and down from x. Here's an example of doing that for the float nearest π:

>>> import math
>>> x = 3.141592653589793
>>> x_plus = math.nextafter(x, math.inf)
>>> x_minus = math.nextafter(x, 0.0)
>>> x_plus, x_minus
(3.1415926535897936, 3.1415926535897927)

(Note that to do this in general, we also need to deal with the special case where x is the largest representable float and math.nextafter(x, math.inf) gives infinity.)

The bounds for the interval that rounds to x are midway between x and the neighbouring floats. Python lets us convert floats to the corresponding exact value as a fraction:

>>> from fractions import Fraction
>>> left = (Fraction(x) + Fraction(x_minus)) / 2
>>> right = (Fraction(x) + Fraction(x_plus)) / 2
>>> print(left, right)
14148475504056879/4503599627370496 14148475504056881/4503599627370496

We also need to know whether we have a closed or an open interval. We could look at the bit representation to figure that out (it depends on whether the least significant bit of the float is 0 or 1), or we could just test to see whether our interval endpoints round to x or not:

>>> float(left) == x
True
>>> float(right) == x
True

They do, so we have a closed interval. That's confirmed by looking at the hex representation of the float:

>>> x.hex()
'0x1.921fb54442d18p+1'

So we can find the simplest fraction that rounds to x using simplest_in_closed_interval:

>>> simplest_in_closed_interval(left, right)
Fraction(245850922, 78256779)
>>> 245850922 / 78256779 == x
True

Putting it all together

While the core algorithm is simple, there are enough corner cases to deal with (negative values, open versus closed intervals, sys.float_info.max, etc.) that a complete solution ends up being a bit too messy to post in full in this answer. Some time ago I put together a Python package called simplefractions that deals with all those corner cases; it's available on PyPI. Here it is in action:

>>> from simplefractions import simplest_from_float
>>> simplest_from_float(0.1)
Fraction(1, 10)
>>> simplest_from_float(-3.3333333333333333)
Fraction(-10, 3)
>>> simplest_from_float(22/7)
Fraction(22, 7)
>>> import math
>>> simplest_from_float(math.pi)
Fraction(245850922, 78256779)