解释 IEEE 754 数字的模定点算法如何工作

Question

在问题的最后，给出了一个从 GNU C 库中计算浮点模运算 (__ieee754_fmod ) 的代码。我对这个算法背后的基本思想很感兴趣。

特别感兴趣的是标记为/ * fix point fmod * /的代码。

    /* fix point fmod */
    n = ix - iy;
    while(n--) {
        hz=hx-hy;
        if(hz<0){hx = hx+hx;}
        else {
            if(hz==0)                /* return sign(x)*0 */
                return Zero[(uint64_t)sx>>63];
            hx = hz+hz;
        }
    }
    hz=hx-hy;
    if(hz>=0) {hx=hz;}

1. 定点运算如何应用于浮点数？

2. 如何更改它以应用于十进制的浮点数？

3. 这种求余数的算法有什么共同点吗？

while(x >= y){
    auto scaled_y = y;
    while(scaled_y*2 < x)
        scaled_y *= 2;
    x -= scaled_y;
}

鼓励使用简化示例或详细说明的链接。

来自 GNU C 库source code：

typedef union
{
    double value;
    struct
    {
        uint32_t lsw;
        uint32_t msw;
    } parts;
    uint64_t word;
} ieee_double_shape_type;


/* Get all in one, efficient on 64-bit machines.  */
#ifndef EXTRACT_WORDS64
# define EXTRACT_WORDS64(i,d)                                        \
do {                                                                \
  ieee_double_shape_type gh_u;                                        \
  gh_u.value = (d);                                                \
  (i) = gh_u.word;                                                \
} while (0)
#endif

/* Get all in one, efficient on 64-bit machines.  */
#ifndef INSERT_WORDS64
# define INSERT_WORDS64(d,i)                                        \
do {                                                                \
  ieee_double_shape_type iw_u;                                        \
  iw_u.word = (i);                                                \
  (d) = iw_u.value;                                                \
} while (0)
#endif


static const double one = 1.0, Zero[] = {0.0, -0.0,};
double __ieee754_fmod (double x, double y)
{
    int32_t n,ix,iy;
    int64_t hx,hy,hz,sx,i;
    EXTRACT_WORDS64(hx,x);
    EXTRACT_WORDS64(hy,y);
    sx = hx&UINT64_C(0x8000000000000000);        /* sign of x */
    hx ^=sx;                                /* |x| */
    hy &= UINT64_C(0x7fffffffffffffff);        /* |y| */
    /* purge off exception values */
    if(__builtin_expect(hy==0
                        || hx >= UINT64_C(0x7ff0000000000000)
                        || hy > UINT64_C(0x7ff0000000000000), 0))
        /* y=0,or x not finite or y is NaN */
        return (x*y)/(x*y);
    if(__builtin_expect(hx<=hy, 0)) {
        if(hx<hy) return x;        /* |x|<|y| return x */
        return Zero[(uint64_t)sx>>63];        /* |x|=|y| return x*0*/
    }
    /* determine ix = ilogb(x) */
    if(__builtin_expect(hx<UINT64_C(0x0010000000000000), 0)) {
        /* subnormal x */
        for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
    } else ix = (hx>>52)-1023;
    /* determine iy = ilogb(y) */
    if(__builtin_expect(hy<UINT64_C(0x0010000000000000), 0)) {        /* subnormal y */
        for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
    } else iy = (hy>>52)-1023;
    /* set up hx, hy and align y to x */
    if(__builtin_expect(ix >= -1022, 1))
        hx = UINT64_C(0x0010000000000000)|(UINT64_C(0x000fffffffffffff)&hx);
    else {                /* subnormal x, shift x to normal */
        n = -1022-ix;
        hx<<=n;
    }
    if(__builtin_expect(iy >= -1022, 1))
        hy = UINT64_C(0x0010000000000000)|(UINT64_C(0x000fffffffffffff)&hy);
    else {                /* subnormal y, shift y to normal */
        n = -1022-iy;
        hy<<=n;
    }
    /* fix point fmod */
    n = ix - iy;
    while(n--) {
        hz=hx-hy;
        if(hz<0){hx = hx+hx;}
        else {
            if(hz==0)                /* return sign(x)*0 */
                return Zero[(uint64_t)sx>>63];
            hx = hz+hz;
        }
    }
    hz=hx-hy;
    if(hz>=0) {hx=hz;}
    /* convert back to floating value and restore the sign */
    if(hx==0)                        /* return sign(x)*0 */
        return Zero[(uint64_t)sx>>63];
    while(hx<UINT64_C(0x0010000000000000)) {        /* normalize x */
        hx = hx+hx;
        iy -= 1;
    }
    if(__builtin_expect(iy>= -1022, 1)) {        /* normalize output */
        hx = ((hx-UINT64_C(0x0010000000000000))|((uint64_t)(iy+1023)<<52));
        INSERT_WORDS64(x,hx|sx);
    } else {                /* subnormal output */
        n = -1022 - iy;
        hx>>=n;
        INSERT_WORDS64(x,hx|sx);
        x *= one;                /* create necessary signal */
    }
    return x;                /* exact output */
}

Answer 1

1. 定点运算如何应用于浮点数？

见下文：

    /* fix point fmod */

此时，我假设（尚未彻底检查先前的代码）ix 和 iy 是输入 x 和 y、hx 和 hy 的浮点指数他们的有效数字和|y| ≤ |x|。注释“fix point fmod”指的是这段代码正在处理浮点数的有效数字。然而，这是一个不恰当的描述；它不是真正的定点fmod，因为它确实使用指数进行一些缩放，正如我们将看到的。

    n = ix - iy;

这会将n 设置为指数差。

    while(n--) {

这将启动一个循环，处理因指数而导致的数字缩放比例不同。

        hz=hx-hy;

这会执行一个试减法；在测试标志之后，我们可能会也可能不会使用这个结果。请注意，无论x 和y 之间的指数有何差异，hy 是y 相对于hx 的一些 倍数。也就是说，即使有效数字可能不会为正常减法对齐，hz 表示的数字与hx 表示的数字具有相同的余数模 y，因为所有数字相差y 有相同的残基。

        if(hz<0){hx = hx+hx;}

如果hz 是负数，我们还不想处理hx-hy。相反，hx 向左移动一位，使其缩放更接近hy。 （回想一下，我们在 n 上循环，因此，最终将 hx 向左移动 n 位将使 hx 和 hy 处于相同的比例。）

        else {
            if(hz==0)                /* return sign(x)*0 */
                return Zero[(uint64_t)sx>>63];
            hx = hz+hz;
        }

如果hz 为零，我们已经证明x 是y 的倍数，所以余数为零，并返回带有适当符号位的余数。否则，hz 为正。在这种情况下，hx = hz+hz; 做了两件事：首先，它将hx 替换为hz，这是一个有效的替换，因为hx 和hz 具有相同的余数模y，但hz 是更小，所以我们稍微减少了它。其次，它将hx 左移一位，为下一次循环迭代做准备。

    }

这会在n 上继续迭代。

    hz=hx-hy;

这是另一个试减法；在测试标志之后，我们可能会或可能不会使用结果。此时，hx 已通过前面的循环以 y 为模减少，我们正在测试是否需要最后一个减法，或者我们已经完成了。

    if(hz>=0) {hx=hz;}

如果hz小于零，循环完全减少hx；它产生的残留物少于hy，我们完成了：hx 包含要使用的残留物。如果hz大于或等于0，则hx没有完全减少，所以我们用hz替换它，我们包含hx减去hy。这是最后需要的减法——我们知道没有更多的了，因为我们在二进制中工作并且hy 的前导位由于规范化而被设置，所以hx/hy 必须小于2，这意味着hx-hy < hy， fmod 的目标是产生小于 |y| 的残基。

1. 如何更改它以应用于十进制表示法的浮点数？

必须对算法进行大量重写才能应用于浮点数。特别是，上面的代码依赖于二进制，因为每一步只需要一个减法；商hx/hy 的整数部分始终为 0 或 1。对于十进制表示，必须对商的整数部分范围从 0 到 9 的事实进行调整，因此可以使用各种倍数必须减去 hy。

1. 这种求余数的算法有什么共同点吗？

while(x >= y){
    auto scaled_y = y;
    while(scaled_y*2 < x)
        scaled_y *= 2;
    x -= scaled_y;

一些。该代码还从x 中减去y 的倍数，以减少x 模y。