最快的算法,用于识别制作双精度方程x + a == b true的最小和最大x

在静态分析的上下文中,我感兴趣的是在以下条件的then-branch中确定x的值:

 double x; x = …; if (x + a == b) { … 

ab可以假设为双精度常量(推广到任意表达式是问题中最容易的部分),并且可以假设编译器严格遵循IEEE 754( FLT_EVAL_METHOD为0)。 运行时的舍入模式可以假设为最接近均匀。

如果用有理数计算是便宜的,那就很简单: x的值将是理性区间中包含的双精度数(b – a – 0.5 * ulp1(b)… b – a + 0.5 * ulp2(b) )。 如果b是偶数,则应该包括边界,如果b是奇数则排除,并且ulp1和ulp2是两个稍微不同的“ULP”定义,如果不介意在2的幂上失去一点精度,则可以采用相同的定义。

不幸的是,使用有理数的计算可能很昂贵。 考虑另一种可能性是通过二分法获得每个边界,在64个双精度加法中(每个操作决定结果的一位)。 获得下限和上限的128个浮点加法可能比基于数学的任何解更快。

我想知道是否有办法改进“128浮点添加”的想法。 实际上我有自己的解决方案,包括更改舍入模式和nextafter调用,但我不想吝啬任何人的风格,并使他们错过比我现有的更优雅的解决方案。 另外我不确定两次更改舍入模式实际上比64个浮点加法更便宜。

您已经在问题中提供了一个漂亮而优雅的解决方案:

如果用有理数计算是便宜的,那就很简单:x的值将是理性区间中包含的双精度数(b – a – 0.5 * ulp1(b)… b – a + 0.5 * ulp2(b) )。 如果b是偶数,则应该包括边界,如果b是奇数则排除,并且ulp1和ulp2是两个稍微不同的“ULP”定义,如果不介意在2的幂上失去一点精度,则可以采用相同的定义。

以下是基于该段的部分解决问题的半合理草图。 希望我很快就有机会充实它。 要获得真正的解决方案,您必须处理次正规,零,NaN以及所有其他有趣的东西。 我会假设ab比如说1e-300 < |a| < 1e300 1e-300 < |a| < 1e3001e-300 < |b| < 1e300 1e-300 < |b| < 1e300以便在任何时候都不会发生疯狂。

如果没有上溢和下溢,你可以从b - nextafter(b, -1.0/0.0)得到ulp1(b) b - nextafter(b, -1.0/0.0) 。 你可以从nextafter(b, 1.0/0.0) - b得到ulp2(b) nextafter(b, 1.0/0.0) - b

如果b/2 <= a <= 2b ,那么Sterbenz定理告诉你b - a是精确的。 所以(b - a) - ulp1 / 2将是最接近下限的double ,而(b - a) + ulp2 / 2将是最接近上限的double 。 尝试这些值,以及之前和之后的值,并选择最有效的最大间隔。

如果b > 2a ,则b - a > b/2b - a的计算值最多偏离半个ulp。 一个ulp1最多是两个ulp,也就是一个ulp2 ,所以你给出的理性间隔最多是两个ulp宽。 找出ba五个最接近的值中的哪一个。

如果a > 2b ,则ba的ul至少与b的ulp一样大; 如果有什么工作,我打赌它必须是与ba最接近的三个值之一。 我想象ab有不同符号的情况类似。

我写了一小堆实现这个想法的C ++代码。 在我厌倦等待之前,它没有失败随机模糊测试(在几个不同的范围内)。 这里是:

 void addeq_range(double a, double b, double &xlo, double &xhi) { if (a != a) return; // empty interval if (b != b) { if (aa != 0) { xlo = xhi = -a; return; } else return; // empty interval } if (bb != 0) { // TODO: handle me. } // b is now guaranteed to be finite. if (aa != 0) return; // empty interval if (b < 0) { addeq_range(-a, -b, xlo, xhi); xlo = -xlo; xhi = -xhi; return; } // b is now guaranteed to be zero or positive finite and a is finite. if (a >= b/2 && a <= 2*b) { double upulp = nextafter(b, 1.0/0.0) - b; double downulp = b - nextafter(b, -1.0/0.0); xlo = (ba) - downulp/2; xhi = (ba) + upulp/2; if (xlo + a == b) { xlo = nextafter(xlo, -1.0/0.0); if (xlo + a != b) xlo = nextafter(xlo, 1.0/0.0); } else xlo = nextafter(xlo, 1.0/0.0); if (xhi + a == b) { xhi = nextafter(xhi, 1.0/0.0); if (xhi + a != b) xhi = nextafter(xhi, -1.0/0.0); } else xhi = nextafter(xhi, -1.0/0.0); } else { double xmid = ba; if (xmid + a < b) { xhi = xlo = nextafter(xmid, 1.0/0.0); if (xhi + a != b) xhi = xmid; } else if (xmid + a == b) { xlo = nextafter(xmid, -1.0/0.0); xhi = nextafter(xmid, 1.0/0.0); if (xlo + a != b) xlo = xmid; if (xhi + a != b) xhi = xmid; } else { xlo = xhi = nextafter(xmid, -1.0/0.0); if (xlo + a != b) xlo = xmid; } } }