Я ищу метод для преобразования точного значения числа с плавающей запятой в рациональное отношение двух целых чисел, то есть 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.
Исходя из моего опыта, обход (путем бинарного поиска) дерева Stern-Brocot здесь бесполезен, поскольку это больше подходит для аппроксимации значения числа с плавающей запятой, когда оно интерпретируется как усеченное действительное вместо точного рационального .
Возможно, продолженные дроби - путь.
Другая проблема здесь - целочисленное переполнение. Подумайте о том, что мы хотим представить рациональное как отношение двух int32_t, где максимальный знаменатель b_max = INT32_MAX. Мы не можем полагаться на критерий остановки, такой как b > b_max. Поэтому алгоритм никогда не должен переполняться или обнаруживать переполнение.
То, что я нашел до сих пор , является алгоритмом из кода Розетты , который основан на непрерывных дробях, но его источник упоминает, что он "все еще не совсем завершен". Некоторые базовые тесты дали хорошие результаты, но я не могу подтвердить их общую правильность и думаю, что они могут легко переполниться.
// 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];
}