Finding the nearest Rational to a double - is there a more efficient mechanism?

I use the following code to find the lowest denominator Rational that is within a certain delta from a double.

The rationale is that the I am pulling float numbers from a database and in many cases summing them. All of the numbers are calculated using simple maths such as +, -, * and /. No transcendental numbers are involved, nor is there any trigonometry. In most cases finding the nearest Rational to the float gets what the original figure is supposed to be rather than the results of adding mashed-up numbers together.

// Create a good rational for the value within the delta supplied.
public static Rational valueOf(double dbl, double delta) {
  // Primary checks.
  if (delta  0.0");
  }
  // Remove the sign and integral part.
  long integral = (long) Math.floor(dbl);
  dbl -= integral;
  // The value we are looking for.
  final Rational d = new Rational((long) ((dbl) / delta), (long) (1 / delta));
  // Min value = d - delta.
  final Rational min = new Rational((long) ((dbl - delta) / delta), (long) (1 / delta));
  // Max value = d + delta.
  final Rational max = new Rational((long) ((dbl + delta) / delta), (long) (1 / delta));
  // Start the fairey sequence.
  Rational l = ZERO;
  Rational h = ONE;
  Rational found = null;
  // Keep slicing until we arrive within the delta range.
  do {
    // Either between min and max -> found it.
    if (found == null && min.compareTo(l) = 0) {
      found = l;
    }
    if (found == null && min.compareTo(h) = 0) {
      found = h;
    }
    if (found == null) {
      // Make the mediant.
      Rational m = mediant(l, h);
      // Replace either l or h with mediant.
      if (m.compareTo(d) < 0) {
        l = m;
      } else {
        h = m;
      }
    }

  } while (found == null);

  // Bring back the sign and the integral.
  if (integral != 0) {
    found = found.plus(new Rational(integral, 1));
  }
  // That's me.
  return found;
}

In a recent test using 0.000001 as my delta this code took 75% of the CPU. Droppi

Naming

You definitely need to improve your variable names. Here are some suggested changes with the reasoning

• dbl -> value (does not change much but repeating the type of a variable as its name does not add value (no pun intended))

• delta -> epsilon (the name epsilon is much more common if you want to define a closeness boundary)

• ABig1 (I have no idea for this but it does not help much. This one is as big as every other double. More importantly, double is not capable of catching all the digits you give there.)

Now to the

Algorithm

Doing a binary search is quite efficient in comparison with other methods but it still has an \$\mathbb{O}(\log n)\$ runtime. Lets think about the representation of doubles in IEEE 754 format.

A double consists of 64 bits of which the first is the sign bit, the next eleven bits store a biased exponent and the remaining 52 bits store the mantissa.

The idea is as follows: We use the mantissa as numerator and the exponent as denominator (if it is negative) or as factor for the numerator when it is positive:

public static long getMantissaBits(double value) {
    // select the 52 lower bits which make up the mantissa
    return Double.doubleToLongBits(value) & 0xFFFFFFFFFFFFFL;
}

public static long getMantissa(double value) {
    // add the "hidden" 1 of normalized doubles
    return (1L > exponentOffset;
    // remove the bias
    return biasedExponent - 1023;
}

public static Rational valueOf(double value) {
    long mantissa = getMantissa(value);
    long exponent = getExponent(value) - 52;
    int numberOfTrailingZeros = Long.numberOfTrailingZeros(mantissa);
    mantissa >>= numberOfTrailingZeros;
    exponent += numberOfTrailingZeros;
    // apply the sign to the numerator
    long numerator = (long) Math.signum(value) * mantissa;
    if(exponent < 0)
        return new Rational(numerator, 1L << -exponent);
    else
        return new Rational(numerator << exponent, 1);
}

As you can see, we don't need the delta anymore as we are as close as possible. Of course there are some caveats. If the double is denormalized getMantissa will give wrong results (but you could detect that and return the correct result). Another problem stems from too big/small exponents where the shift of 1L is greater or equal to 64bits (and thus the result is 0). However, if you think about this problem you will find that these exponents only occur when the number is too big/small to fit into your Rational class anyways.

As noticed in the comments this will return the exact result which is not wanted. I tried to come up with a solution to get the rounding correct but I failed so I shamelessly translated python's solution to Java:

public Rational limitDenominator(long maximumDenominator) {
    if (maximumDenominator  maximumDenominator)
            break;
        long oldP0 = p0;
        p0 = p1;
        q0 = q1;
        p1 = oldP0 + a * p1;
        q1 = q2;
        long oldN = n;
        n = d;
        d = oldN - a * d;
    }
    long k = (maximumDenominator - q0) / q1;
    Rational bound1 = new Rational(p0 + k * p1, q0 + k * q1);
    Rational bound2 = new Rational(p1, q1);
    if(bound2.minus(this).abs().compareTo(bound1.minus(this).abs()) <= 0){
        return bound2;
    } else {
        return bound1;
    }
}

I cannot say much about the exact details in play here because I need to first understand them myself but the idea is that you find the closest fraction with a denominator less or equal than the given maximum. To do so you find an upper and lower closest and choose the one that is closer.

Algorithm explanation (some math?)

I finally found some time to have a closer look at the algorithm. Let us at first note that the problem we are trying to solve is the Diophantine approximation. As noted in the linked article we can use convergents or semiconvergents of the continued fraction representation of the given number to approximate it.


Its best approximations for the second definition are


\$3, \frac{8}{3}, \frac{11}{4}, \frac{19}{7}, \frac{87}{32}, \ldots\,\$ ,


while, for the first definition, they are


\$3, \frac{5}{2}, \frac{8}{3}, \frac{11}{4}, \frac{19}{7}, \frac{30}{11}, \frac{49}{18}, \frac{68}{25}, \frac{87}{32}, \frac{106}{39}, \ldots \$

As you can see the semiconvergents have a much tighter spacing so we should use them (because the target number can lie in a hole between the convergents that is filled by semiconvergents).

If you take a look into the linked python documentation you will find that it also uses semiconvergents.

The Wikipedia article also tells us that


Even-numbered convergents are smaller than the original number, while odd-numbered ones are bigger.

So we corner the number from both sides. The same section also tells us that


If successive convergents are found, with numerators h1, h2, … and denominators k1, k2, … then the relevant recursive relation is:


\$h_n = a_nh_{n − 1} + h_{n − 2}

StackExchange Code Review Q#48577, answer score: 7