Algorithm for finding the ratio of two floating-point numbers?

Question

I need to find the ratio of one floating-point number to another, and the ratio needs to be two integers. For example:

• input: 1.5, 3.25
• output: "6:13"

Does anyone know of one? Searching the internet, I found no such algorithm, nor an algorithm for the least common multiple or denominator of two floating-point numbers (just integers).

Java implementation:

---

This is the final implementation that I will use:

public class RatioTest
{
  public static String getRatio(double d1, double d2)//1.5, 3.25
  {
    while(Math.max(d1,d2) < Long.MAX_VALUE && d1 != (long)d1 && d2 != (long)d2)
    {
      d1 *= 10;//15 -> 150
      d2 *= 10;//32.5 -> 325
    }
    //d1 == 150.0
    //d2 == 325.0
    try
    {
      double gcd = getGCD(d1,d2);//gcd == 25
      return ((long)(d1 / gcd)) + ":" + ((long)(d2 / gcd));//"6:13"
    }
    catch (StackOverflowError er)//in case getGDC (a recursively looping method) repeats too many times
    {
      throw new ArithmeticException("Irrational ratio: " + d1 + " to " + d2);
    }
  }

  public static double getGCD(double i1, double i2)//(150,325) -> (150,175) -> (150,25) -> (125,25) -> (100,25) -> (75,25) -> (50,25) -> (25,25)
  {
    if (i1 == i2)
      return i1;//25
    if (i1 > i2)
      return getGCD(i1 - i2, i2);//(125,25) -> (100,25) -> (75,25) -> (50,25) -> (25,25)
    return getGCD(i1, i2 - i1);//(150,175) -> (150,25)
  }
}

• -> indicates the next stage in the loop or method call

Mystical's implementation as Java:

---

Though I did not end up using this, it more than deserves to be recognized, so I translated it to Java, so I could understand it:

import java.util.Stack;

public class RatioTest
{
    class Fraction{
        long num;
        long den;
        double val;
    };

    Fraction build_fraction(Stack<long> cf){
        long term = cf.size();
        long num = cf[term - 1];
        long den = 1;
        while (term-- > 0){
            long tmp = cf[term];

            long new_num = tmp * num + den;
            long new_den = num;

            num = new_num;
            den = new_den;
        }

        Fraction f;
        f.num = num;
        f.den = den;
        f.val = (double)num / (double)den;

        return f;
    }

    void get_fraction(double x){
        System.out.println("x = " + x);

        //  Generate Continued Fraction
        System.out.print("Continued Fraction: ");
        double t = Math.abs(x);
        double old_error = x;
        Stack<long> cf;
        Fraction f;
        do{
            //  Get next term.
            long tmp = (long)t;
            cf.push(tmp);

            //  Build the current convergent
            f = build_fraction(cf);

            //  Check error
            double new_error = Math.abs(f.val - x);
            if (tmp != 0 && new_error >= old_error){
                //  New error is bigger than old error.
                //  This means that the precision limit has been reached.
                //  Pop this (useless) term and break out.
                cf.pop();
                f = build_fraction(cf);
                break;
            }
            old_error = new_error;
            System.out.print(tmp + ", ");

            //  Error is zero. Break out.
            if (new_error == 0)
                break;

            t -= tmp;
            t = 1/t;
        }while (cf.size() < 39); //  At most 39 terms are needed for double-precision.
        System.out.println();System.out.println();

        //  Print Results
        System.out.println("The fraction is:   " + f.num + " / " + f.den);
        System.out.println("Target x = " + x);
        System.out.println("Fraction = " + f.val);
        System.out.println("Relative error is: " + (Math.abs(f.val - x) / x));System.out.println();
        System.out.println();
    }
    public static void main(String[] args){
        get_fraction(15.38 / 12.3);
        get_fraction(0.3333333333333333333);    //  1 / 3
        get_fraction(0.4184397163120567376);    //  59 / 141
        get_fraction(0.8323518818409020299);    //  1513686 / 1818565
        get_fraction(3.1415926535897932385);    //  pi
    }
}

One MORE thing:

---

The above mentioned implemented way of doing this works IN THEORY, however, due to floating-point rounding errors, this results in alot of unexpected exceptions, errors, and outputs. Below is a practical, robust, but a bit dirty implementation of a ratio-finding algorithm (Javadoc'd for your convenience):

public class RatioTest
{
  /** Represents the radix point */
  public static final char RAD_POI = '.';

  /**
   * Finds the ratio of the two inputs and returns that as a <tt>String</tt>
   * <h4>Examples:</h4>
   * <ul>
   * <li><tt>getRatio(0.5, 12)</tt><ul>
     *   <li>returns "<tt>24:1</tt>"</li></ul></li>
   * <li><tt>getRatio(3, 82.0625)</tt><ul>
   *   <li>returns "<tt>1313:48</tt>"</li></ul></li>
   * </ul>
   * @param d1 the first number of the ratio
   * @param d2 the second number of the ratio
   * @return the resulting ratio, in the format "<tt>X:Y</tt>"
   */
  public static strictfp String getRatio(double d1, double d2)
  {
    while(Math.max(d1,d2) < Long.MAX_VALUE && (!Numbers.isCloseTo(d1,(long)d1) || !Numbers.isCloseTo(d2,(long)d2)))
    {
      d1 *= 10;
      d2 *= 10;
    }
    long l1=(long)d1,l2=(long)d2;
    try
    {
      l1 = (long)teaseUp(d1); l2 = (long)teaseUp(d2);
      double gcd = getGCDRec(l1,l2);
      return ((long)(d1 / gcd)) + ":" + ((long)(d2 / gcd));
    }
    catch(StackOverflowError er)
    {
      try
      {
        double gcd = getGCDItr(l1,l2);
        return ((long)(d1 / gcd)) + ":" + ((long)(d2 / gcd));
      }
      catch (Throwable t)
      {
        return "Irrational ratio: " + l1 + " to " + l2;
      }
    }
  }


  /**
   * <b>Recursively</b> finds the Greatest Common Denominator (GCD)
   * @param i1 the first number to be compared to find the GCD
   * @param i2 the second number to be compared to find the GCD
   * @return the greatest common denominator of these two numbers
   * @throws StackOverflowError if the method recurses to much
   */
  public static long getGCDRec(long i1, long i2)
  {
    if (i1 == i2)
      return i1;
    if (i1 > i2)
      return getGCDRec(i1 - i2, i2);
    return getGCDRec(i1, i2 - i1);
  }

  /**
   * <b>Iteratively</b> finds the Greatest Common Denominator (GCD)
   * @param i1 the first number to be compared to find the GCD
   * @param i2 the second number to be compared to find the GCD
   * @return the greatest common denominator of these two numbers
   */
  public static long getGCDItr(long i1, long i2)
  {
    for (short i=0; i < Short.MAX_VALUE &&  i1 != i2; i++)
    {
      while (i1 > i2)
        i1 = i1 - i2;
      while (i2 > i1)
        i2 = i2 - i1;
    }
      return i1;
  }

  /**
   * Calculates and returns whether <tt>d1</tt> is close to <tt>d2</tt>
   * <h4>Examples:</h4>
   * <ul>
   * <li><tt>d1 == 5</tt>, <tt>d2 == 5</tt>
   *   <ul><li>returns <tt>true</tt></li></ul></li>
   * <li><tt>d1 == 5.0001</tt>, <tt>d2 == 5</tt>
   *   <ul><li>returns <tt>true</tt></li></ul></li>
   * <li><tt>d1 == 5</tt>, <tt>d2 == 5.0001</tt>
   *   <ul><li>returns <tt>true</tt></li></ul></li>
   * <li><tt>d1 == 5.24999</tt>, <tt>d2 == 5.25</tt>
   *   <ul><li>returns <tt>true</tt></li></ul></li>
   * <li><tt>d1 == 5.25</tt>, <tt>d2 == 5.24999</tt>
   *   <ul><li>returns <tt>true</tt></li></ul></li>
   * <li><tt>d1 == 5</tt>, <tt>d2 == 5.1</tt>
   *   <ul><li>returns <tt>false</tt></li></ul></li>
   * </ul>
   * @param d1 the first number to compare for closeness
   * @param d2 the second number to compare for closeness
   * @return <tt>true</tt> if the two numbers are close, as judged by this method
   */
  public static boolean isCloseTo(double d1, double d2)
  {
    if (d1 == d2)
      return true;
    double t;
    String ds = Double.toString(d1);
    if ((t = teaseUp(d1-1)) == d2 || (t = teaseUp(d2-1)) == d1)
      return true;
    return false;
  }

  /**
   * continually increases the value of the last digit in <tt>d1</tt> until the length of the double changes
   * @param d1
   * @return
   */
  public static double teaseUp(double d1)
  {
    String s = Double.toString(d1), o = s;
    byte b;
    for (byte c=0; Double.toString(extractDouble(s)).length() >= o.length() && c < 100; c++)
      s = s.substring(0, s.length() - 1) + ((b = Byte.parseByte(Character.toString(s.charAt(s.length() - 1)))) == 9 ? 0 : b+1);
    return extractDouble(s);
  }

  /**
   * Works like Double.parseDouble, but ignores any extraneous characters. The first radix point (<tt>.</tt>) is the only one treated as such.<br/>
   * <h4>Examples:</h4>
   * <li><tt>extractDouble("123456.789")</tt> returns the double value of <tt>123456.789</tt></li>
   * <li><tt>extractDouble("1qw2e3rty4uiop[5a'6.p7u8&9")</tt> returns the double value of <tt>123456.789</tt></li>
   * <li><tt>extractDouble("123,456.7.8.9")</tt> returns the double value of <tt>123456.789</tt></li>
   * <li><tt>extractDouble("I have $9,862.39 in the bank.")</tt> returns the double value of <tt>9862.39</tt></li>
   * @param str The <tt>String</tt> from which to extract a <tt>double</tt>.
   * @return the <tt>double</tt> that has been found within the string, if any.
   * @throws NumberFormatException if <tt>str</tt> does not contain a digit between 0 and 9, inclusive.
   */
  public static double extractDouble(String str) throws NumberFormatException
  {
    try
    {
      return Double.parseDouble(str);
    }
    finally
    {
      boolean r = true;
      String d = "";
      for (int i=0; i < str.length(); i++)
        if (Character.isDigit(str.charAt(i)) || (str.charAt(i) == RAD_POI && r))
        {
          if (str.charAt(i) == RAD_POI && r)
            r = false;
          d += str.charAt(i);
        }
      try
      {
        return Double.parseDouble(d);
      }
      catch (NumberFormatException ex)
      {
        throw new NumberFormatException("The input string could not be parsed to a double: " + str);
      }
    }
  }
}

Forevermore

Convert them both to fractions, and then output that ratio (after reducing)? – Sordid

@Supuhstar: I just replaced my pseudo-code with a fully working implementation in C++. It shouldn't be too hard to convert to Java. – Weaks

EDIT: I just tested my version more thoroughly. It's able to perfectly detect fractions that are at most 4 - 5 digits long in both numerator and denominator. Anything longer and it's hit and miss. (For the Java version, change the 4294967296 to 2147483648 or there could be an integer overflow during the cast.) – Weaks

@Supuhstar: For your updated code, scaling by 2 instead of 10 will get rid of those round-off errors. But it will give different (but still valid) answers. My implementation also suffers heavily from round-off. (hence why it doesn't work if the fraction is larger than 5 digits.) There's a way to fix this, but it's even messier... – Weaks

Answer 1

Assuming you have a data type that can handle arbitrarily large numeric values, you can do something like this:

1. Multiply both values by 10 until the significand is entirely to the left of the decimal point.
2. Find the Greatest Common Denominator of the two values.
3. Divide by GCD

So for you example you would have something like this:

a = 1.5
b = 3.25

multiply by 10:  15, 32.5
multiply by 10:  150, 325

find GCD:  25

divide by GCD:  6, 13

Shelleyshellfire

See my other comment. It depends on the numbers and what fraction you want. These scaling methods will not always give you fractions such as 1/3 if that's what the desired answer is. – Weaks

@Mysticial: Well said. Didn't think of that case. – Shelleyshellfire

@Weaks that's just fine; it will do for my needs. Than you for reminding anyone who will see this question for reference, though! – Forevermore

Why multiply by 10? Why not 2, or 16, or 37? This method only works if both numbers actually terminate in base 10; i.e. if their own denominators have no factors other than 2 and 5. 10 is about the most arbitrary value you can use here... At least using 2 would have some rationale based on the internal representation. (But @Mystical's solution is the right one.) – Roundel

@Roundel Well, at least 10 is a multiple of 2. Still a terrible number to use. – Asur

Answer 2

This is a fairly non-trivial task. The best approach I know that gives reliable results for any two floating-point is to use continued fractions.

First, divide the two numbers to get the ratio in floating-point. Then run the continued fraction algorithm until it terminates. If it doesn't terminate, then it's irrational and there is no solution.

If it terminates, evaluate the resulting continued fraction back into a single fraction and that will be the answer.

Of course, there is no reliable way to determine if there is a solution or not since this becomes the halting problem. But for the purposes of limited precision floating-point, if the sequence doesn't terminate with a reasonable number of steps, then assume there's no answer.

EDIT 2: Here's an update to my original solution in C++. This version is much more robust and seems to work with any positive floating-point number except for INF, NAN, or extremely large or small values that would overflow an integer.

typedef unsigned long long  uint64;
struct Fraction{
    uint64 num;
    uint64 den;
    double val;
};
Fraction build_fraction(vector<uint64> &cf){
    uint64 term = cf.size();
    uint64 num = cf[--term];
    uint64 den = 1;
    while (term-- > 0){
        uint64 tmp = cf[term];

        uint64 new_num = tmp * num + den;
        uint64 new_den = num;

        num = new_num;
        den = new_den;
    }

    Fraction f;
    f.num = num;
    f.den = den;
    f.val = (double)num / den;

    return f;
}
void get_fraction(double x){
    printf("x = %0.16f\n",x);

    //  Generate Continued Fraction
    cout << "Continued Fraction: ";
    double t = abs(x);
    double old_error = x;
    vector<uint64> cf;
    Fraction f;
    do{
        //  Get next term.
        uint64 tmp = (uint64)t;
        cf.push_back(tmp);

        //  Build the current convergent
        f = build_fraction(cf);

        //  Check error
        double new_error = abs(f.val - x);
        if (tmp != 0 && new_error >= old_error){
            //  New error is bigger than old error.
            //  This means that the precision limit has been reached.
            //  Pop this (useless) term and break out.
            cf.pop_back();
            f = build_fraction(cf);
            break;
        }
        old_error = new_error;
        cout << tmp << ", ";

        //  Error is zero. Break out.
        if (new_error == 0)
            break;

        t -= tmp;
        t = 1/t;
    }while (cf.size() < 39); //  At most 39 terms are needed for double-precision.
    cout << endl << endl;

    //  Print Results
    cout << "The fraction is:   " << f.num << " / " << f.den << endl;
    printf("Target x = %0.16f\n",x);
    printf("Fraction = %0.16f\n",f.val);
    cout << "Relative error is: " << abs(f.val - x) / x << endl << endl;
    cout << endl;
}
int main(){
    get_fraction(15.38 / 12.3);
    get_fraction(0.3333333333333333333);    //  1 / 3
    get_fraction(0.4184397163120567376);    //  59 / 141
    get_fraction(0.8323518818409020299);    //  1513686 / 1818565
    get_fraction(3.1415926535897932385);    //  pi
    system("pause");
}

Output:

x = 1.2504065040650407
Continued Fraction: 1, 3, 1, 152, 1,

The fraction is:   769 / 615
Target x = 1.2504065040650407
Fraction = 1.2504065040650407
Relative error is: 0


x = 0.3333333333333333
Continued Fraction: 0, 3,

The fraction is:   1 / 3
Target x = 0.3333333333333333
Fraction = 0.3333333333333333
Relative error is: 0


x = 0.4184397163120567
Continued Fraction: 0, 2, 2, 1, 1, 3, 3,

The fraction is:   59 / 141
Target x = 0.4184397163120567
Fraction = 0.4184397163120567
Relative error is: 0


x = 0.8323518818409020
Continued Fraction: 0, 1, 4, 1, 27, 2, 7, 1, 2, 13, 3, 5,

The fraction is:   1513686 / 1818565
Target x = 0.8323518818409020
Fraction = 0.8323518818409020
Relative error is: 0


x = 3.1415926535897931
Continued Fraction: 3, 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, 3,

The fraction is:   245850922 / 78256779
Target x = 3.1415926535897931
Fraction = 3.1415926535897931
Relative error is: 0


Press any key to continue . . .

The thing to notice here is that it gives 245850922 / 78256779 for pi. Obviously, pi is irrational. But as far as double-precision allows, 245850922 / 78256779 isn't any different from pi.

Basically, any fraction with 8 - 9 digits in the numerator/denominator has enough entropy to cover nearly all DP floating-point values (barring corner cases like INF, NAN, or extremely large/small values).

Weaks

Interesting... can I get a pseudo-code snippet? I can't quite visualize it in my head ^^; – Forevermore

Generally you would impose an upper limit on the size of the denominator you're willing to use..perhaps denominators less than 100, or able to fit in a 16-bit or 32-bit integer.. The continued fraction method generates successively better and better approximations with monotonically increasing denominators, so you just stop when the approximation is "good enough". – Hydra

so, where did the "large" variable come from? – Forevermore

It's just arbitrary. For double precision, I would imagine it to be about 80 - 90 terms. But in any case, read Jim Lewis' comment above this. The more terms you need, the longer the fraction gets, but the better the approximation becomes. If the answer is something simple like 6/13 in your example, the continued fraction sequence will be very short: for 6/13, it is [0, 2, 6] - which is only 3 terms. – Weaks

I just replaced the pseudo-code with a fully working implementation in C++. It didn't turn out as bad as I thought it would be. :) – Weaks

Thanks for that - I translated it to Java as best I could, so I could understand and compare it. – Forevermore

Irrational? You've dividing two rational numbers, how is the result ever going to be irrational? – Ario

@harold since the user is inputing two floating-point numbers, one or both could be INFINITY, NEGATIVE_INFINITY, or NaN. Also, there's not always a rational fraction for two decimal numbers. What's pi over 7, for instance? – Forevermore

@Supuhstar ok there is infinity and NaN, but apart from those (if they show up here, there is no sane answer anyway) floats are always rationals, always (though not all rationals fit in a float, obviously). There is no pi in floats, just an approximation. This is even true for decimal floats (which shouldn't exist but whatever) – Ario

point is it's always good to have a fallback, just in case your user finds that one situation where you need it. – Forevermore

Answer 3

If the floating point numbers have a limit to decimal places - then just multiply both numbers by 10^n where n is limit - so for 2 decimal places, multiply by 100, then calculate for whole numbers - the ratio will be the same for the original decimals because it is a ratio.

Horsehair

good idea - If I assume they are both carried out to the language's maximum limit, then I could work from there. – Forevermore

This approach may or may not work depending on what you want your fraction to be. If your two numbers are 1.0 and 0.333333333333333. This method of scaling by 10 will give a fraction of '333333333333333 / 1000000000000000'. But this method won't give you 1/3, if that's what you desire. – Weaks

You should multiply by a power of two, not a power of ten. Since these are presumed to be binary-floating point numbers it is clear that (ignoring overflow) multiplying by two is lossless, whereas multiplying by ten is not. As long as there is no overflow then this technique will always work. The hypothetical problem of 1.0 and 0.333... will work perfectly. It won't give you 1/3 but that's because the ratio is not actually 1/3. – Honeymoon

Answer 4

In Maxima CAS simply :

(%i1) rationalize(1.5/3.5);
(%o1) 7720456504063707/18014398509481984

Code from numeric.lisp :

;;; This routine taken from CMUCL, which, in turn is a routine from
;;; CLISP, which is GPL.
;;;
;;; I (rtoy) have modified it from CMUCL so that it only handles bigfloats.
;;;
;;; RATIONALIZE  --  Public
;;;
;;; The algorithm here is the method described in CLISP.  Bruno Haible has
;;; graciously given permission to use this algorithm.  He says, "You can use
;;; it, if you present the following explanation of the algorithm."
;;;
;;; Algorithm (recursively presented):
;;;   If x is a rational number, return x.
;;;   If x = 0.0, return 0.
;;;   If x < 0.0, return (- (rationalize (- x))).
;;;   If x > 0.0:
;;;     Call (integer-decode-float x). It returns a m,e,s=1 (mantissa,
;;;     exponent, sign).
;;;     If m = 0 or e >= 0: return x = m*2^e.
;;;     Search a rational number between a = (m-1/2)*2^e and b = (m+1/2)*2^e
;;;     with smallest possible numerator and denominator.
;;;     Note 1: If m is a power of 2, we ought to take a = (m-1/4)*2^e.
;;;       But in this case the result will be x itself anyway, regardless of
;;;       the choice of a. Therefore we can simply ignore this case.
;;;     Note 2: At first, we need to consider the closed interval [a,b].
;;;       but since a and b have the denominator 2^(|e|+1) whereas x itself
;;;       has a denominator <= 2^|e|, we can restrict the seach to the open
;;;       interval (a,b).
;;;     So, for given a and b (0 < a < b) we are searching a rational number
;;;     y with a <= y <= b.
;;;     Recursive algorithm fraction_between(a,b):
;;;       c := (ceiling a)
;;;       if c < b
;;;         then return c       ; because a <= c < b, c integer
;;;         else
;;;           ; a is not integer (otherwise we would have had c = a < b)
;;;           k := c-1          ; k = floor(a), k < a < b <= k+1
;;;           return y = k + 1/fraction_between(1/(b-k), 1/(a-k))
;;;                             ; note 1 <= 1/(b-k) < 1/(a-k)
;;;
;;; You can see that we are actually computing a continued fraction expansion.
;;;
;;; Algorithm (iterative):
;;;   If x is rational, return x.
;;;   Call (integer-decode-float x). It returns a m,e,s (mantissa,
;;;     exponent, sign).
;;;   If m = 0 or e >= 0, return m*2^e*s. (This includes the case x = 0.0.)
;;;   Create rational numbers a := (2*m-1)*2^(e-1) and b := (2*m+1)*2^(e-1)
;;;   (positive and already in lowest terms because the denominator is a
;;;   power of two and the numerator is odd).
;;;   Start a continued fraction expansion
;;;     p[-1] := 0, p[0] := 1, q[-1] := 1, q[0] := 0, i := 0.
;;;   Loop
;;;     c := (ceiling a)
;;;     if c >= b
;;;       then k := c-1, partial_quotient(k), (a,b) := (1/(b-k),1/(a-k)),
;;;            goto Loop
;;;   finally partial_quotient(c).
;;;   Here partial_quotient(c) denotes the iteration
;;;     i := i+1, p[i] := c*p[i-1]+p[i-2], q[i] := c*q[i-1]+q[i-2].
;;;   At the end, return s * (p[i]/q[i]).
;;;   This rational number is already in lowest terms because
;;;   p[i]*q[i-1]-p[i-1]*q[i] = (-1)^i.
;;;
(defmethod rationalize ((x bigfloat))
  (multiple-value-bind (frac expo sign)
      (integer-decode-float x)
    (cond ((or (zerop frac) (>= expo 0))
       (if (minusp sign)
           (- (ash frac expo))
           (ash frac expo)))
      (t
       ;; expo < 0 and (2*m-1) and (2*m+1) are coprime to 2^(1-e),
       ;; so build the fraction up immediately, without having to do
       ;; a gcd.
       (let ((a (/ (- (* 2 frac) 1) (ash 1 (- 1 expo))))
         (b (/ (+ (* 2 frac) 1) (ash 1 (- 1 expo))))
         (p0 0)
         (q0 1)
         (p1 1)
         (q1 0))
         (do ((c (ceiling a) (ceiling a)))
         ((< c b)
          (let ((top (+ (* c p1) p0))
            (bot (+ (* c q1) q0)))
            (/ (if (minusp sign)
               (- top)
               top)
               bot)))
           (let* ((k (- c 1))
              (p2 (+ (* k p1) p0))
              (q2 (+ (* k q1) q0)))
         (psetf a (/ (- b k))
            b (/ (- a k)))
         (setf p0 p1
               q0 q1
               p1 p2
               q1 q2))))))))

Bettyebettzel

What does this add that the other answers don't? – Forevermore

Answer 5

I'm using the folowing algorithm. It's fast and simple. It uses the fact that 10^N = 2^N * 5^N and it also handles recurring patterns of digits ! I hope it will help you.

Fraction-to-ratio-converter

Some demos are also provided at that side.

Homovec