有理数到浮点数

Rational to floating point

本文关键字：浮点数 有理数      更新时间：2023-10-16

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

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

假设double的IEEE-754二进制64表示，如何将结构转换为最接近的double并正确舍入?将n和d转换为double并将它们分开的简单方法会导致舍入误差。

实现预期结果的一种方法是在整数空间中执行除法。由于标准C/c++不提供128位整数类型(虽然一些工具链可能作为扩展提供此类型)，因此这不是很有效，但它将产生正确的结果。

下面的代码产生54个商位和一个余数位，每次1位。最有效的53个商位代表double结果的尾数部分，而最不有效的商位和余数需要根据IEEE-754进行舍入到"最接近或偶数"。

下面的代码可以被编译为C或c++程序(至少它可以用我的工具链)。它只经受了轻微的考验。由于按位处理，这不是非常快，并且可以进行各种优化，特别是在使用特定于机器的数据类型和intrinsic时。

#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);
}
相关文章：