寻找两个浮点数比率的算法？

Question

寻找两个浮点数比率的算法？

algorithmfloating-pointinteger

17

我需要找到一个浮点数相对于另一个浮点数的比率，并且这个比率需要是两个整数。例如:

输入: 1.5, 3.25
输出: "6:13"

有人知道如何实现吗？我在搜索引擎上没有找到这样的算法，也没有找到用于两个浮点数（而不是整数）的最小公倍数或最小公约数的算法。

Java实现代码:

这是最终的实现代码：

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)
  }
}

-> 表示循环或方法调用中的下一个阶段

Mystical的Java实现：

尽管我最终没有使用它，但它值得被认可，因此我将其翻译成了Java，以便我能够理解：

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
    }
}

还有一件事：

上述实现方法在理论上是可行的，但由于浮点数舍入误差，会导致很多意外的异常、错误和输出。下面是一个实用、健壮但有点“脏”的比率查找算法的实现方式（为方便起见，已进行了Javadoc注释）：

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);
      }
    }
  }
}

- Ky -

1

将它们都转换为分数，然后输出该比例（在简化之后）？ - BlackJack

@Supuhstar：我刚刚用C++实现了我的伪代码，应该不难转换成Java。 - Mysticial

我刚刚更全面地测试了我的版本。它能够完美地检测最多有4-5位数的分子和分母之间的分数。如果超过这个长度，结果就会不尽如人意。（对于Java版本，请将“4294967296”更改为“2147483648”，否则在转换过程中可能会发生整数溢出。） - Mysticial

@Supuhstar：对于你更新的代码，将比例因子从10改为2可以消除那些舍入误差。但这会得到不同（但仍然有效）的答案。我的实现也严重受到舍入误差的影响。（因此，如果分数大于5位数字，则无法正常工作。）有一种方法可以解决这个问题，但它更加混乱... - Mysticial

5个回答

9

假设您有一个可以处理任意大的数字值的数据类型，您可以按照以下方式执行操作：

将两个数字值乘以10，直到小数点完全位于左侧。
找出这两个值的最大公约数。
除以GCD

因此，对于您的示例，您将得到以下结果：

a = 1.5
b = 3.25
乘以10:  15, 32.5
乘以10:  150, 325
找到GCD:  25
除以GCD:  6, 13

- Austin Salonen

7

看我的另一条评论。具体取决于数字和你想要的分数是多少。这些缩放方法并不总能给你像“1/3”这样的分数，如果这是所需的答案。 - Mysticial

@Mysticial：说得好。没想到那种情况。 - Austin Salonen

@Mysticial，那很好；这对我的需求来说足够了。不过感谢您提醒任何看到这个问题的人作为参考！ - Ky -

为什么要乘以10？为什么不是2、16或37？这种方法只有在两个数字都在十进制下终止时才有效；即它们自己的分母除了2和5之外没有其他因子。10是你可以在这里使用的最任意的值...至少使用2会基于内部表示有一些合理性。（但@Mystical的解决方案是正确的。） - Nemo

@Nemo 好在10是2的倍数，但仍然不是一个好的数字。 - Pascal Cuoq

1

如果浮点数有小数位的限制 - 那么只需将两个数字乘以10^n，其中n是限制 - 因此对于2个小数位，乘以100，然后计算整数 - 比率将与原始小数相同，因为它是比率。

- Billy Moon

好主意 - 如果我假设它们都达到了语言的最大限度，那么我可以从那里开始工作。 - Ky -

4

这种方法可能有效，也可能无效，具体取决于你想要的分数是什么。如果你的两个数字分别是1.0和0.333333333333333，那么这种放大10倍的方法将得到一个分数为“333333333333333 / 1000000000000000”。但如果你想得到的是1/3这个分数，那么这种方法就行不通了。 - Mysticial

你应该乘以二的幂，而不是十的幂。由于这些被假定为二进制浮点数，因此很明显（忽略溢出），乘以二是无损的，而乘以十则不是。只要没有溢出，那么这种技术总是有效的。1.0和0.333...的假设问题将完美地解决。它不会给你1/3，但那是因为比率实际上不是1/3。 - Bruce Dawson

0

在Maxima CAS中，只需简单地：

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

来自 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))))))))

- Adam

这个回答有什么其他答案没有的优点？ - Ky -

0

我正在使用以下算法。它快速简单，利用了这个事实：10^N = 2^N * 5^N，并且还可以处理数字的重复模式！希望它能对你有所帮助。

分数转比例转换器

该网站还提供了一些演示。

- Matej Gomboc

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Mysticial · Accepted Answer

这是一个相当复杂的任务。我知道的最好方法是使用连分数，可以为任何两个浮点数提供可靠的结果。

首先，将两个数字除以得到浮点比率。然后运行连分数算法，直到它终止。如果它不终止，则它是无理数，没有解决方案。

如果它终止，将结果连分数计算回单个分数，这将是答案。

当然，没有可靠的方法确定是否有解决方案，因为这变成了停机问题。但对于有限精度浮点数而言，如果序列不能在合理的步骤内终止，则假设没有答案。

编辑2：这是我在C++中原始解决方案的更新版本。这个版本更加健壮，并且似乎适用于除“INF”、“NAN”或极大或极小的值（会导致整数溢出）之外的任何正浮点数。

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");
}

输出：

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 . . .

这里需要注意的是，它给出了pi的值为245850922 / 78256779。显然，pi是无理数。但就双精度而言，245850922 / 78256779与pi没有任何区别。基本上，任何分子/分母中有8-9位数字的分数都具有足够的熵来覆盖几乎所有DP浮点值（除了像INF、NAN或极大/极小值这样的特殊情况）。