浮点数的精确值作为有理数

9
我正在寻找一种方法,将浮点数的确切值转换为两个整数的有理商,即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 这样的停止准则。因此,算法必须永远不会溢出,或者必须检测溢出。 到目前为止我找到的 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
你可能会对Python的Fraction.limit_denominator感兴趣,它恰好可以做到这一点。当然,这并不能解决C中的溢出问题,但该算法可能会有用。具体而言,在Python中,对于一个(Python)float xFraction(x).limit_denominator(b_max)可以得到你想要的结果。 - Mark Dickinson
你应该专注于算法。通过库函数(如boost::multiprecision),计算的精度可以得到扩展(以内存和运行时间为代价)。 - Red.Wave
最大公约数对于分母是2的幂次方并没有太大帮助,它最多只能移除一些零位。你需要一个大整数表示来获得精确值。 - Dietmar Kühl
1
chux 给出了一个可能适用于 IEEE-754 二进制格式的答案(我还没有检查过)。为了完整起见,我要指出 C 提供了 FLT_RADIX,它定义了浮点格式所使用的基数。乘以或除以 FLT_RADIX 应该是精确的(请注意 scalbn 函数),因此,给定一些不是整数的浮点值 x,可以通过将其乘以 FLT_RADIX 直到结果为整数来找到所需的分子,并且分母是 FLT_RADIX 的幂次方。 - Eric Postpischil
@plasmacel:乘以非规格化数应该没有问题。通过乘以FLT_RADIX直到一个数字成为整数,不会在任何正常格式中溢出。理论上,可以定义一种浮点格式,使指数限制在奇怪的范围内,如-20到-10,而不是更正常的-126到127,但实际上,指数范围允许有效数字被缩放为整数。 - Eric Postpischil
显示剩余14条评论
1个回答

8

所有有限的double都是有理数,正如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的考虑留到以后。


更快捷的代码可以通过用ldexp(1, expo)@gammatesterexp2(expo)@Bob__替换pow(2.0, expo)来实现。 while (*numerator && fmod(*numerator,2) == 0 && fmod(*denominator,2) == 0)也可以进行一些性能优化。但首先,让我们按需求获得所需功能。

https://en.cppreference.com/w/c/numeric/math/exp2 是否可以比通用的 pow(2.0,x) 更高效或更准确地实现? - Bob__
如果我理解正确,小于1.0 / pow(2, MAX_EXP)的数字(即“太小”)会导致分子是2.0的负幂。这种情况可以通过将分子和分母都乘以1.0 / numerator来解决。或者这些情况已经具有最大可表示的分母了吗? - plasmacel
@chux,我的假设是正确的吗?任何浮点数都可以用相同类型的两个浮点数的商来表示吗?还有另一个优化想法:最大公约数(GCD,由欧几里得算法实现)除法是否比末尾的while循环更有效?我认为它会更快地收敛。 - plasmacel
在C++代码中,考虑使用boost::rational<long>,它在其构造函数中对分子和分母进行了归一化(但必须将double值显式转换为long),因此不需要使用while (*numerator && ... - Tomáš Kolárik
@TomášKolárik 尽管 std::numeric_limits<double>::digits 可以作为更好的 C++ 答案,但是 OP 更倾向于使用 <*.h> 文件来获得类似 C 的答案。如果你有更符合惯用法的 C++ 答案,请随意发表。据我所知,boost:: 不是标准 C++ 的一部分,因此它的包含偏离了标准 C++。同样,您可以通过对 boost 的深入了解来回答问题。 - chux - Reinstate Monica
显示剩余13条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接