Точное значение числа с плавающей точкой как рациональное
Я ищу метод для преобразования точного значения числа с плавающей запятой в рациональное отношение двух целых чисел, т.е. a / b, где b не больше указанного максимального знаменателя b_max, Если удовлетворяет условию b <= b_max невозможно, тогда результат возвращается к наилучшему приближению, которое все еще удовлетворяет условию.
Оставайтесь на линии. Здесь много вопросов / ответов о наилучшем рациональном приближении усеченного действительного числа, представленного в виде числа с плавающей запятой. Однако меня интересует точное значение числа с плавающей запятой, которое само по себе является рациональным числом с другим представлением. Более конкретно, математический набор чисел с плавающей точкой является подмножеством рациональных чисел. В случае двоичного стандарта IEEE 754 это подмножество двоичных рациональных чисел. В любом случае, любое число с плавающей точкой может быть преобразовано в рациональное отношение двух целых чисел конечной точности как a / b,
Так, например, предполагая, что двоичный формат с плавающей запятой IEEE 754 одинарной точности, является рациональным эквивалентом float f = 1.0f / 3.0f не является 1 / 3, но 11184811 / 33554432, Это точное значение f, который представляет собой число из математического набора двоичных чисел с плавающей точкой одинарной точности IEEE 754.
Исходя из моего опыта, обход (путем бинарного поиска) дерева Штерна-Броко здесь бесполезен, поскольку он больше подходит для аппроксимации значения числа с плавающей запятой, когда оно интерпретируется как усеченное действительное число вместо точного рациональный
Возможно, продолжение фракции - это путь.
Другая проблема здесь - целочисленное переполнение. Подумайте о том, что мы хотим представить рациональное как отношение двух int32_tгде максимальный знаменатель b_max = INT32_MAX, Мы не можем полагаться на такой критерий остановки, как b > b_max, Таким образом, алгоритм никогда не должен переполняться или обнаруживать переполнение.
До сих пор я обнаружил алгоритм из Rosetta Code, который основан на непрерывных дробях, но его источник упоминает, что он "все еще не вполне завершен". Некоторые базовые тесты дали хорошие результаты, но я не могу подтвердить их общую правильность и думаю, что они могут легко переполниться.
// https://rosettacode.org/wiki/Convert_decimal_number_to_rational#C
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>
/* f : number to convert.
* num, denom: returned parts of the rational.
* md: max denominator value. Note that machine floating point number
* has a finite resolution (10e-16 ish for 64 bit double), so specifying
* a "best match with minimal error" is often wrong, because one can
* always just retrieve the significand and return that divided by
* 2**52, which is in a sense accurate, but generally not very useful:
* 1.0/7.0 would be "2573485501354569/18014398509481984", for example.
*/
void rat_approx(double f, int64_t md, int64_t *num, int64_t *denom)
{
/* a: continued fraction coefficients. */
int64_t a, h[3] = { 0, 1, 0 }, k[3] = { 1, 0, 0 };
int64_t x, d, n = 1;
int i, neg = 0;
if (md <= 1) { *denom = 1; *num = (int64_t) f; return; }
if (f < 0) { neg = 1; f = -f; }
while (f != floor(f)) { n <<= 1; f *= 2; }
d = f;
/* continued fraction and check denominator each step */
for (i = 0; i < 64; i++) {
a = n ? d / n : 0;
if (i && !a) break;
x = d; d = n; n = x % n;
x = a;
if (k[1] * a + k[0] >= md) {
x = (md - k[0]) / k[1];
if (x * 2 >= a || k[1] >= md)
i = 65;
else
break;
}
h[2] = x * h[1] + h[0]; h[0] = h[1]; h[1] = h[2];
k[2] = x * k[1] + k[0]; k[0] = k[1]; k[1] = k[2];
}
*denom = k[1];
*num = neg ? -h[1] : h[1];
}
1 ответ
Все конечно double рациональные числа как ОП хорошо сформулированы..
использование frexp() разбить число на его дробь и показатель степени. Конечный результат еще нужно использовать double представлять значения целых чисел из-за требований диапазона. Некоторые цифры слишком малы, (x меньше чем 1.0/(2.0,DBL_MAX_EXP)) и бесконечность, а не число являются вопросами.
frexpфункции разбивают число с плавающей запятой на нормализованную дробь и интегральную степень 2. .. интервал [1/2, 1) или ноль...
C11 §7.12.6.4 2/3
#include <math.h>
#include <float.h>
_Static_assert(FLT_RADIX == 2, "TBD code for non-binary FP");
// Return error flag
int split(double x, double *numerator, double *denominator) {
if (!isfinite(x)) {
*numerator = *denominator = 0.0;
if (x > 0.0) *numerator = 1.0;
if (x < 0.0) *numerator = -1.0;
return 1;
}
int bdigits = DBL_MANT_DIG;
int expo;
*denominator = 1.0;
*numerator = frexp(x, &expo) * pow(2.0, bdigits);
expo -= bdigits;
if (expo > 0) {
*numerator *= pow(2.0, expo);
}
else if (expo < 0) {
expo = -expo;
if (expo >= DBL_MAX_EXP-1) {
*numerator /= pow(2.0, expo - (DBL_MAX_EXP-1));
*denominator *= pow(2.0, DBL_MAX_EXP-1);
return fabs(*numerator) < 1.0;
} else {
*denominator *= pow(2.0, expo);
}
}
while (*numerator && fmod(*numerator,2) == 0 && fmod(*denominator,2) == 0) {
*numerator /= 2.0;
*denominator /= 2.0;
}
return 0;
}
void split_test(double x) {
double numerator, denominator;
int err = split(x, &numerator, &denominator);
printf("e:%d x:%24.17g n:%24.17g d:%24.17g q:%24.17g\n",
err, x, numerator, denominator, numerator/ denominator);
}
int main(void) {
volatile float third = 1.0f/3.0f;
split_test(third);
split_test(0.0);
split_test(0.5);
split_test(1.0);
split_test(2.0);
split_test(1.0/7);
split_test(DBL_TRUE_MIN);
split_test(DBL_MIN);
split_test(DBL_MAX);
return 0;
}
Выход
e:0 x: 0.3333333432674408 n: 11184811 d: 33554432 q: 0.3333333432674408
e:0 x: 0 n: 0 d: 9007199254740992 q: 0
e:0 x: 1 n: 1 d: 1 q: 1
e:0 x: 0.5 n: 1 d: 2 q: 0.5
e:0 x: 1 n: 1 d: 1 q: 1
e:0 x: 2 n: 2 d: 1 q: 2
e:0 x: 0.14285714285714285 n: 2573485501354569 d: 18014398509481984 q: 0.14285714285714285
e:1 x: 4.9406564584124654e-324 n: 4.4408920985006262e-16 d: 8.9884656743115795e+307 q: 4.9406564584124654e-324
e:0 x: 2.2250738585072014e-308 n: 2 d: 8.9884656743115795e+307 q: 2.2250738585072014e-308
e:0 x: 1.7976931348623157e+308 n: 1.7976931348623157e+308 d: 1 q: 1.7976931348623157e+308
Оставь b_max рассмотрение на потом.
Более целесообразный код возможен с заменой pow(2.0, expo) с ldexp(1, expo) @gammatester или exp2(expo) @Bob__
while (*numerator && fmod(*numerator,2) == 0 && fmod(*denominator,2) == 0) также можно использовать некоторые улучшения производительности. Но сначала давайте получим функциональность по мере необходимости.