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

Question

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) {
    throw new IllegalArgumentException("Delta must be > 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 && max.compareTo(l) >= 0) {
      found = l;
    }
    if (found == null && min.compareTo(h) <= 0 && max.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. Dropping it to 0.0001 reduced that dramatically but it is still a significant bottleneck.

Is there a quicker way of doing this?

My implementation of Rational forces the numerator and denominator to be fully reduced at all times. I accept that that is likely the biggest overhead but as mentioned in Wikipedia the rationals must be fully reduced for the mediant function to work correctly.

Here is the full class - borrowed from the mentioned site and enhanced:

/**
 * ***********************************************************************
 * Immutable ADT for Rational numbers.
 *
 * Invariants
 * -----------
 * - gcd(num, den) = 1, i.e, the rational number is in reduced form
 * - den >= 1, the denominator is always a positive integer
 * - 0/1 is the unique representation of 0
 *
 * We employ some tricks to stave of overflow, but if you
 * need arbitrary precision rationals, use BigRational.java.
 *
 * Borrowed from http://introcs.cs.princeton.edu/java/92symbolic/Rational.java.html
 * because it has a mediant method.
 *
 ************************************************************************
 */
public class Rational extends Number implements Comparable<Rational> {
  public static final Rational ZERO = new Rational(0, 1);
  public static final Rational ONE = new Rational(1, 1);

  private long num;   // the numerator
  private long den;   // the denominator

  // create and initialize a new Rational object
  public Rational(long numerator, long denominator) {

    // deal with x/0
    if (denominator == 0) {
      throw new IllegalArgumentException("Denominator cannot be 0.");
    }
    // reduce fraction
    long g = gcd(numerator, denominator);
    num = numerator / g;
    den = denominator / g;

    // only needed for negative numbers
    if (den < 0) {
      den = -den;
      num = -num;
    }
  }

  public Rational(Rational from) {
    num = from.num;
    den = from.den;
  }

  // return the numerator and denominator of (this)
  public long numerator() {
    return num;
  }

  public long denominator() {
    return den;
  }

  // return double precision representation of (this)
  public double toDouble() {
    return (double) num / den;
  }

  public BigDecimal toBigDecimal() {
    // Do it to just 4 decimal places.
    return toBigDecimal(4);
  }

  public BigDecimal toBigDecimal(int digits) {
    // Do it to n decimal places.
    return new BigDecimal(num).divide(new BigDecimal(den), digits, RoundingMode.DOWN).stripTrailingZeros();
  }

  // return string representation of (this)
  @Override
  public String toString() {
    if (den == 1) {
      return num + "";
    } else {
      return num + "/" + den;
    }
  }

  public int compareTo(Rational b) {
    // return { -1, 0, +1 } if a < b, a = b, or a > b
    Rational a = this;
    long lhs = a.num * b.den;
    long rhs = a.den * b.num;
    if (lhs < rhs) {
      return -1;
    }
    if (lhs > rhs) {
      return +1;
    }
    return 0;
  }

  @Override
  public boolean equals(Object y) {
    // is this Rational object equal to y?
    if (y == null) {
      return false;
    }
    if (y.getClass() != this.getClass()) {
      return false;
    }
    Rational b = (Rational) y;
    return compareTo(b) == 0;
  }

  @Override
  public int hashCode() {
    int hash = 5;
    hash = 97 * hash + (int) (this.num ^ (this.num >>> 32));
    hash = 97 * hash + (int) (this.den ^ (this.den >>> 32));
    return hash;
  }

  // create and return a new rational (r.num + s.num) / (r.den + s.den)
  public static Rational mediant(Rational r, Rational s) {
    return new Rational(r.num + s.num, r.den + s.den);
  }

  // return gcd(|m|, |n|)
  private static long gcd(long m, long n) {
    if (m < 0) {
      m = -m;
    }
    if (n < 0) {
      n = -n;
    }
    if (0 == n) {
      return m;
    } else {
      return gcd(n, m % n);
    }
  }

  // return lcm(|m|, |n|)
  private static long lcm(long m, long n) {
    if (m < 0) {
      m = -m;
    }
    if (n < 0) {
      n = -n;
    }
    return m * (n / gcd(m, n));    // parentheses important to avoid overflow
  }

  // return a * b, staving off overflow as much as possible by cross-cancellation
  public Rational times(Rational b) {
    Rational a = this;

    // reduce p1/q2 and p2/q1, then multiply, where a = p1/q1 and b = p2/q2
    Rational c = new Rational(a.num, b.den);
    Rational d = new Rational(b.num, a.den);
    return new Rational(c.num * d.num, c.den * d.den);
  }

  // return a + b, staving off overflow
  public Rational plus(Rational b) {
    Rational a = this;

    // special cases
    if (a.compareTo(ZERO) == 0) {
      return b;
    }
    if (b.compareTo(ZERO) == 0) {
      return a;
    }

    // Find gcd of numerators and denominators
    long f = gcd(a.num, b.num);
    long g = gcd(a.den, b.den);

    // add cross-product terms for numerator
    Rational s = new Rational((a.num / f) * (b.den / g) + (b.num / f) * (a.den / g),
                              lcm(a.den, b.den));

    // multiply back in
    s.num *= f;
    return s;
  }

  // return -a
  public Rational negate() {
    return new Rational(-num, den);
  }

  // return a - b
  public Rational minus(Rational b) {
    return plus(b.negate());
  }

  public Rational reciprocal() {
    return new Rational(den, num);
  }

  // return a / b
  public Rational divides(Rational b) {
    Rational a = this;
    return a.times(b.reciprocal());
  }

  // Default delta to apply.
  public static final double DELTA = 0.0001;

  public static Rational valueOf(double dbl) {
    return valueOf(dbl, DELTA);
  }

  public static Rational valueOf(BigDecimal dbl) {
    return valueOf(dbl.doubleValue(), DELTA);
  }

  public static Rational valueOf(double dbl, int digits) {
    return valueOf(dbl, Math.pow(10, -digits));
  }

  public static Rational valueOf(BigDecimal dbl, int digits) {
    return valueOf(dbl.doubleValue(), Math.pow(10, -digits));
  }

  // 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) {
      throw new IllegalArgumentException("Delta must be > 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 && max.compareTo(l) >= 0) {
        found = l;
      }
      if (found == null && min.compareTo(h) <= 0 && max.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;
  }

  private static void print(String name, Rational r) {
    System.out.println(name + "=" + r + "(" + r.toDouble() + ")");
  }

  private enum TestNumber {
    OneTenth(0.100000001490116119384765625),
    Pi(Math.PI),
    E(Math.E),
    OneThird(0.3333333333333),
    MinusOneThird(-0.3333333333333),
    ABig1(1.87344227533222758141533568138280569154340745619495504034120344898213260187710089517712780269958755185722145694193999220);

    final double v;

    TestNumber(double v) {
      this.v = v;
    }

  }

  private static void test2() {
    for (TestNumber n : TestNumber.values()) {
      print(n.name(), Rational.valueOf(n.v));
    }
  }

  private static void test1() {
    Rational x, y, z;

    // 1/2 + 1/3 = 5/6
    x = new Rational(1, 2);
    y = new Rational(1, 3);
    z = x.plus(y);
    System.out.println(z);

    // 8/9 + 1/9 = 1
    x = new Rational(8, 9);
    y = new Rational(1, 9);
    z = x.plus(y);
    System.out.println(z);

    // 1/200000000 + 1/300000000 = 1/120000000
    x = new Rational(1, 200000000);
    y = new Rational(1, 300000000);
    z = x.plus(y);
    System.out.println(z);

    // 1073741789/20 + 1073741789/30 = 1073741789/12
    x = new Rational(1073741789, 20);
    y = new Rational(1073741789, 30);
    z = x.plus(y);
    System.out.println(z);

    //  4/17 * 17/4 = 1
    x = new Rational(4, 17);
    y = new Rational(17, 4);
    z = x.times(y);
    System.out.println(z);

    // 3037141/3247033 * 3037547/3246599 = 841/961 
    x = new Rational(3037141, 3247033);
    y = new Rational(3037547, 3246599);
    z = x.times(y);
    System.out.println(z);

    // 1/6 - -4/-8 = -1/3
    x = new Rational(1, 6);
    y = new Rational(-4, -8);
    z = x.minus(y);
    System.out.println(z);
  }

  // test client
  public static void main(String[] args) {
    //test1();
    test2();
  }

  // Implement Number.
  @Override
  public int intValue() {
    return (int) doubleValue();
  }

  @Override
  public long longValue() {
    return (long) doubleValue();
  }

  @Override
  public float floatValue() {
    return (float) doubleValue();
  }

  @Override
  public double doubleValue() {
    return toDouble();
  }

}

Answer 1

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 << 52) + getMantissaBits(value);
}

public static long getExponent(double value) {
    int exponentOffset = 52;
    long lowest11Bits = 0x7FFL;
    long shiftedBiasedExponent = Double.doubleToLongBits(value) & (lowest11Bits << exponentOffset);
    long biasedExponent = shiftedBiasedExponent >> 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 < 1) {
        throw new IllegalArgumentException("Denominator cannot be less than 1.");
    }
    if(this.den <= maximumDenominator)
        // we can't get closer than the current value
        return this;
    long p0 = 0;
    long q0 = 1;
    long p1 = 1;
    long q1 = 0;
    long n = this.num;
    long d = this.den;
    while(true) {
        long a = n / d;
        long q2 = q0 + a * q1;
        if(q2 > 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}\$,
\$k_n = a_nk_{n − 1} + k_{n − 2}\$.

The successive convergents are given by the formula

\$\frac{h_n}{k_n} = \frac{a_nh_{n − 1} + h_{n − 2}}{a_nk_{n − 1} + k_{n − 2}}\$

Looking closer in the code you will note that \$h_{n-1}\$ corresponds to p1 and \$h_{n-2}\$ to p0 (similarly for the \$k\$s and qs), \$a_n\$ is a in the code.

So the loop calculates convergents that approximate better and better and always increase in their denominator. This gives a convenient break condition when the next convergent would have a denominator greater than the allowed maximum. So we hit the end of the road and stop there.

Now we found the nearest convergents but there might be a closer semiconvergent. We know that the convergent with the denominator q0 + a * q1 is too big and the one with the denominator q1 is smaller (or equal) so we only have to look if there are semiconvergents between the two that are closer. Any semiconvergent with bigger denominator will be a better approximation ^{[Citation Needed]} (say closer). So we need to find the semiconvergent with the biggest denominator that is still smaller than the maximum. The following lines do this exactly:

long k = (maximumDenominator - q0) / q1;
denominator = q0 + k * q1;
numerator = p0 + k * p1;

The python code does compare the closeness of this with \$\frac{p1}{q1}\$ but I am not sure if this is really necessary. This concludes my discussion of the algorithm. The next step would be to explain why the semi-/convergents work that way but I will leave that out just believe Wikipedia/some houndred years of mathematics.

Answer 2

Why not simply have java.util.BigDecimal do the heavy lifting :

public Rational valueOf(double d) {
    BigDecimal bigDecimal = BigDecimal.valueOf(d);
    long numerator = bigDecimal.unscaledValue().longValue();
    long denominator = pow(10L, bigDecimal.scale());
    return new Rational(numerator, denominator);
}

private static long pow(long base, int exponent) {
    long result = 1;
    for (int i = 0; i < exponent; i++) {
        result *= base;
    }
    return result;
}

This simply needs some extra checks for values that are out of range.