Algorithm for finding the ratio of two floating-point numbers?
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开发者_JS百科); // 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);
}
}
}
}
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).
Assuming you have a data type that can handle arbitrarily large numeric values, you can do something like this:
- Multiply both values by 10 until the significand is entirely to the left of the decimal point.
- Find the Greatest Common Denominator of the two values.
- 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
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.
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))))))))
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.
加载中,请稍侯......
精彩评论