Мне нужно найти отношение одного числа с плавающей запятой к другому, а отношение должно быть целым числом. Например:
1.5, 3.25"6:13"Кто-нибудь знает об этом? Поиск в Интернете не нашел такого алгоритма, ни алгоритма для наименьшего общего множителя или знаменателя двух чисел с плавающей запятой (просто целые числа).
Это конечная реализация, которую я буду использовать:
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)
}
}
-> указывает следующий этап в вызове цикла или методаХотя я и не закончил использовать это, это более чем заслуживает признания, поэтому я перевел его на 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
}
}
Вышеупомянутый реализованный способ выполнения работы IN THEORY, однако, из-за ошибок округления с плавающей запятой, это приводит к большим количеством неожиданных исключений, ошибок и выходов. Ниже приведена практическая, надежная, но немного грязная реализация алгоритма определения отношения (Javadoc'd для вашего удобства):
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);
}
}
}
}
Предполагая, что у вас есть тип данных, который может обрабатывать сколь угодно большие числовые значения, вы можете сделать что-то вроде этого:
Итак, для примера у вас будет что-то вроде этого:
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
Это довольно нетривиальная задача. Лучший подход, который я знаю, который дает надежные результаты для любых двух плавающих точек, заключается в использовании продолженных дробей.
Сначала разделите два числа, чтобы получить соотношение в плавающей запятой. Затем запустите алгоритм непрерывной дроби, пока он не завершится. Если это не заканчивается, то оно иррационально и не существует решения.
Если он завершается, оцените полученную непрерывную дроби обратно в одну фракцию, и это будет ответ.
Конечно, нет надежного способа определить, есть ли решение или нет, так как это становится проблемой остановки. Но для целей ограниченной точности с плавающей запятой, если последовательность не заканчивается с разумным числом шагов, тогда не принимайте ответа.
РЕДАКТИРОВАТЬ 2:
Здесь обновление моего исходного решения на С++. Эта версия гораздо более надежна и, похоже, работает с любым положительным числом с плавающей запятой, за исключением 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 . . .
Здесь следует отметить, что он дает 245850922 / 78256779 для pi. Очевидно, pi иррационально. Но что касается двойной точности, 245850922 / 78256779 не отличается от pi.
В принципе, любая фракция с 8-9 цифрами в числителе/знаменателе имеет достаточную энтропию для покрытия почти всех значений с плавающей запятой DP (запрет таких угловых случаев, как INF, NAN или очень большие/малые значения).
Если числа с плавающей запятой имеют предел десятичным знакам, то просто умножьте оба числа на 10 ^ n, где n - предел - так что для двух знаков после запятой умножьте на 100, а затем вычислите для целых чисел - отношение будет то же самое для исходных десятичных знаков, потому что это отношение.
В 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))))))))