有理数到浮点数

Posted 2023-02-17

【中文标题】有理数到浮点数

【英文标题】：Rational to floating point

【发布时间】：2015-10-05 18:26:12

【问题描述】：

考虑一个由以下结构表示的有理数。

struct rational 
    uint64_t n;
    uint64_t d;
    unsigned char sign : 1;
;

假设double 的IEEE-754 binary64 表示，如何通过正确的舍入将结构转换为最接近的double？将n 和d 转换为double 并清楚地划分它们的简单方法会导致舍入误差。

【问题讨论】：

C 还是 C++？选择一个...

没关系。如果它困扰您，我会将 Java 添加到列表中。

@LightnessRacesinOrbit：据我所知，C 和 C++ 具有相同或足够相似的算术规则（我不了解 Java）。另一方面，我同意指定实际使用的语言是一件好事。另一方面，如果这只是一道数学题，为什么会出现在编程网站上？

虽然我担心 OP 可能无意中引起了三个优秀的大型 SO 子社区的愤怒，但这些问题非常有效、有趣，而且实际上大部分与语言无关。但这并不意味着它与编程无关，因为 IEEE 754 标准非常关注编程语言的可用性，或多或少地附加了包袱。

@KeithThompson：显然，“java”标签是由于幼稚的恶意而出现的。

【参考方案1】：

实现所需结果的一种方法是在整数空间中执行除法。由于标准 C/C++ 不提供 128 位整数类型（而某些工具链可能会提供此作为扩展），这不是很有效，但会产生正确的结果。

下面的代码生成 54 个商位和一个余数，一次一位。最高有效的 53 个商位表示 double 结果的尾数部分，而最低有效商位和余数是根据 IEEE-754 舍入到“最接近或偶数”所需要的。

下面的代码可以编译为 C 或 C++ 程序（至少在我的工具链中是这样）。它已经过轻微测试。由于按位处理，这不是很快，并且可以进行各种优化，尤其是在使用特定于机器的数据类型和内在函数的情况下。

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <stdint.h>

struct rational 
    uint64_t n;
    uint64_t d;
    unsigned char sign : 1;
;

double uint64_as_double (uint64_t a)

    double res;
#if defined (__cplusplus)
    memcpy (&res, &a, sizeof (res));
#else /* __cplusplus */
    volatile union 
        double f;
        uint64_t i;
     cvt;
    cvt.i = a;
    res = cvt.f;
#endif /* __cplusplus */
    return res;


#define ADDcc(a,b,cy,t0,t1) (t0=(b), t1=(a), t0=t0+t1, cy=t0<t1, t0=t0)
#define ADDC(a,b,cy,t0,t1) (t0=(b)+cy, t1=(a), t0+t1)
#define SUBcc(a,b,cy,t0,t1) (t0=(b), t1=(a), cy=t1<t0, t1-t0)

double rational2double (struct rational a)

    uint64_t dividend, divisor, quot, rem, t0, t1, cy, res, expo;
    int sticky, round, odd, sign, i;

    dividend = a.n;
    divisor = a.d;
    sign = a.sign;

    /* handle special cases */
    if ((dividend == 0) && (divisor == 0)) 
        res = 0xFFF8000000000000ULL; /* NaN INDEFINITE */
     else if (dividend == 0)             
        res = (uint64_t)sign << 63; /* zero */
     else if (divisor == 0) 
        res = ((uint64_t)sign << 63) | 0x7ff0000000000000ULL; /* Inf */
     
    /* handle normal cases */
    else 
        quot = dividend;
        rem = 0;
        expo = 0;
        /* normalize operands using 128-bit shifts */
        while (rem < divisor) 
            quot = ADDcc (quot, quot, cy, t0, t1);
            rem = ADDC (rem, rem, cy, t0, t1);
            expo--;
        
        /* integer bit of quotient is known to be 1 */
        rem = rem - divisor;
        quot = quot + 1;
        /* generate 53 more quotient bits */
        for (i = 0; i < 53; i++) 
            quot = ADDcc (quot, quot, cy, t0, t1);
            rem = ADDC (rem, rem, cy, t0, t1);
            rem = SUBcc (rem, divisor, cy, t0, t1);
            if (cy) 
                rem = rem + divisor;
             else 
                quot = quot + 1;
            
        
        /* round to nearest or even */
        sticky = rem != 0;
        round = quot & 1;
        quot = quot >> 1;
        odd = quot & 1;
        if (round && (sticky || odd)) 
            quot++;
        
        /* compose normalized IEEE-754 double-precision number */
        res = ((uint64_t)sign << 63) + ((expo + 64 + 1023 - 1) << 52) + quot;
    
    return uint64_as_double (res);

【讨论】：

谢谢。当我意识到长除法产生了精确的数字时，我就知道解决方案应该是这样的。干杯！