IEEE 754数的模定点算法的数学基础是什么

问题描述 投票：0回答：1

在问题的最后，给出了来自 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 库 源代码：

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 */
}
c++ c algorithm floating-point ieee-754
1个回答
3
投票
  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
的 
some
倍数。也就是说，即使正常减法的有效数字可能不对齐，由 
hz
表示的数字与由 
y
表示的数字具有相同的剩余模 
hx
，因为所有数字相差的倍数
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
大于或等于零，
hx
没有完全减少，所以我们用
hz
代替它，我们包含通过减去
hx
减少的
hy
。这是最后需要的减法——我们知道没有更多，因为我们在二进制中工作，并且 
hy
的前导位由于归一化而设置，所以 
hx/hy
必须小于二，这意味着 
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;

一些。该代码还从 

y
中减去 
x
的倍数，目的是减少 
x
模 
y
.

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.