Точное значение числа с плавающей точкой как рациональное

Я ищу способ преобразовать точный значение числа с плавающей точкой в рациональный частное от двух целых чисел, т.е. 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, Таким образом, алгоритм никогда не должен переполняться или обнаруживать переполнение.

Что я нашел до сих пор является алгоритм из кода Розетты, который основан на непрерывных дробях, но его источник упоминает, что он «все еще не совсем завершен». Некоторые базовые тесты дали хорошие результаты, но я не могу подтвердить их общую правильность и думаю, что они могут легко переполниться.

// 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];
}