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库源代码：

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;做了两件事：首先，它用hz替换hx，这是一个有效的替换，因为hx和hz在模y上有相同的残基，但hz较小，所以我们对其进行了一些缩减。其次，它将hx左移一位，为下一次循环迭代做准备。

}

这将继续对n进行迭代。

hz=hx-hy;

这是又一次尝试性减法；在测试了符号之后，我们可能会使用结果，也可能不会使用结果。在这一点上，hx已经被前面的循环以y为模进行了减法运算，我们正在测试是否需要最后一次减法运算，或者我们已经完成了减法运算。

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

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

2. 如何将其更改为应用于十进制表示法中的浮点数

必须对算法进行大量重写才能应用于浮点数。特别地，上面的代码依赖于二进制，因为每一步只需要一个减法；商CCD_ 60/CCD_。对于十进制表示，必须考虑到商的整数部分可能在0到9之间，因此必须减去hy的各种倍数。

3. 它与这个算法在寻找余数方面有一些共同点吗

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

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