浮点数的精确值作为有理数
Posted
技术标签:
【中文标题】浮点数的精确值作为有理数【英文标题】：Exact value of a floating-point number as a rational 【发布时间】：2018-07-02 19:00:20 【问题描述】：我正在寻找一种方法来将浮点数的精确值转换为两个整数的有理商，即a / b，其中b 不大于指定的最大分母b_max。如果满足条件b <= b_max 是不可能的，那么结果会退回到仍然满足条件的最佳近似值。
等一下。这里有很多关于截断实数数的最佳有理近似的问题/答案，表示为一个浮点数。但是我对浮点数的 exact 值感兴趣，它本身就是一个具有不同表示的有理数。更具体地说，浮点数的数学集合是有理数的子集。对于 IEEE 754 二进制浮点标准，它是dyadic rationals 的子集。无论如何，任何浮点数都可以转换为两个有限精度整数的有理商a / b。
因此，例如假设 IEEE 754 单精度二进制浮点格式，float f = 1.0f / 3.0f 的有理等价物不是1 / 3，而是11184811 / 33554432。这是f 的精确值，它是来自 IEEE 754 单精度二进制浮点数的数学集合。
根据我的经验，遍历（通过二进制搜索）Stern-Brocot tree 在这里没有用处，因为当浮点数被解释为截断时，它更适合近似浮点数的值真实的 而不是精确的有理。
可能，continued fractions 是要走的路。
这里的另一个问题是整数溢出。想想我们想要将有理数表示为两个int32_t 的商，其中最大分母b_max = INT32_MAX。我们不能依赖像b > b_max 这样的停止标准。所以算法一定不能溢出，否则它必须检测溢出。
我目前发现的是an algorithm from 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];
【问题讨论】：
itu.dk/~sestoft/bachelor/IEEE754_article.pdf @πάνταῥεῖ 我知道这篇文章。它在这里有什么帮助？ 2 基数的表示只能表示所有可能有理分数的子集。 您可能对 Python 的Fraction.limit_denominator 感兴趣，它正是这样做的。当然，这对 C 中的溢出问题没有帮助，但该算法可能有用。具体来说，在 Python 中，对于（Python）floatx，Fraction(x).limit_denominator(b_max) 给出你想要的。
chux 给出了一个可能适用于 IEEE-754 二进制格式的答案（我还没有检查过）。为了完整起见，我将指出 C 提供了FLT_RADIX，它定义了用于浮点格式的基数。乘以或除以FLT_RADIX 应该是精确的（并注意scalbn 函数），因此，给定一些不是整数的浮点值x，可以通过乘以FLT_RADIX 找到所需的分子，直到结果是一个整数，分母是FLT_RADIX乘法次数的幂。
【参考方案1】：
所有有限的double 都是rational numbers，正如 OP 所说..
使用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) 也可以使用一些性能改进。但首先，让我们根据需要获取功能。
【讨论】：
不错！在我接受它之前，我很快就会做一些测试。您知道如何处理“太小”的值吗？ @plasmacelb_max 的局限性甚至类型值得进一步思考。为了更精确地呈现整数 type 解决方案，我会从 intmax_t num, uintmax_t den; 开始，以获得更多范围。
@plasmacel "too small" --> 我认为这个问题可能会变成“integer / power-of-2”，然后范围问题就不再那么重要了（除了如何编码无穷大和 NaN，嗯）。
我认为那些太小的值，更准确地说是1.0 / pow(2, DBL_MAX_EXP) 实际上是不正常的。
@plasmacel“太小的值”对于binary64 来说是不正常的（或低于正常的），但这个答案并不假定 C 格式没有指定。这只是最常见的格式，并且在那里是正确的。以上是关于浮点数的精确值作为有理数的主要内容，如果未能解决你的问题，请参考以下文章