正确舍入两个处理溢出的浮点数之和的 sqrt 计算
Correctly rounded computation of sqrt of sum of two floats handling overflow
有没有好的方法来计算
的正确舍入结果
sqrt(a+b)
对于浮点数 a 和 b（相同精度），其中 0<=a<+inf 和 0<=b<+inf?
特别是对于 a+b 的计算会溢出的输入值？
(此处“正确舍入”的含义与 sqrt 本身的计算相同，即返回最接近以无限精度计算的“真实”结果的可表示值。)
（注意：一种明显的方法是以更大的浮点大小进行计算并避免溢出。不幸的是，这通常不起作用（例如，如果不支持更大的浮点格式))
我试过Herbie，但它完全放弃了。它似乎没有对 a+b 溢出的足够点进行采样以检测问题，而且似乎也没有很好地处理相关采样。不幸的是，它通常是一个很棒的工具。
到目前为止我一直在做的是（伪代码）
if a + b would overflow:
2*sqrt(a/4 + b/4) # Cannot overflow for finite inputs, as f::MAX/4 + f::MAX/4 <= f::MAX
else:
... # handle non-overflow case. Also interesting; not quite the topic of this question.
...这似乎主要在实践中起作用，但是 a) 完全没有原则和 b) 在实践中偶尔 returns 结果在溢出避免部分被 epsilon 关闭（例如在真正的结果是 x + 0.2(x.next_larger()-x) 但这个 returns x.next_larger() 而不是 x)
有关 f32 中“off-by-epsilon”问题的快速示例：
>>> import decimal
>>> decimal.getcontext().prec = 256
>>> from decimal import Decimal as D
>>> from numpy import float32 as f32
>>> a = D(f32("6.0847234e31").astype(float))
>>> b = D(f32("3.4028235e38").astype(float))
>>> res_act = (a+b).sqrt()
>>> res_calc = D(f32("1.8446744e19").astype(float)) # 2*sqrt(a/4 + b/4) in f32 precision
>>> res_best = D(f32("1.8446746e19").astype(float)) # obtained by brute-force
>>> abs(res_calc - res_act) > abs(res_best - res_act)
True # oops
（您必须相信我对 f32 计算结果的承诺，因为 Python 通常以 f64 精度运行。这也是 f32 舞蹈的原因。）
通过适当地按 2 的幂缩放可以很容易地避免溢出，这样量级大的参数被缩放到统一。困难的部分是产生正确的舍入结果。由于双舍入的潜在问题，我什至不完全相信在下一个更大的 IEEE-754 二进制浮点类型中执行中间计算可以保证。
在没有更广泛的浮点类型的情况下，人们将不得不退回到将多个本机精度数字链接在一起以执行具有更高中间精度的操作。 Dekker 提出的一种常见方案称为对精度。它使用成对的浮点数，其中较重要的部分通常称为“头”，不太重要的部分称为“尾”。将这两个部分归一化，使尾部的大小最多为头部大小的半个 ulp。
此方案中的有效有效位数为 2*p+1，其中 p 是基础浮点类型中的有效位数。 “额外”位由尾部的符号位表示。重要的是要注意，与底层基本类型相比，指数范围没有变化，因此我们需要相当积极地向统一性扩展，以避免在中间计算中遇到次正规操作数。对精度计算不能保证正确舍入的结果。使用三胞胎可能会奏效，但需要付出更多的努力，我无法负担得起答案。
但是，对精度可以提供忠实四舍五入且几乎总是正确四舍五入的结果。当 FMA（融合乘加）可用时，可以相当有效地构建基于 Newton-Raphson 的对精度平方根，产生大约 2*p-1 个好位。这就是我在下面的示例性 IS0-C99 代码中使用的，它使用 float 映射到 IEEE-754 binary32 作为本机浮点类型。成对精度代码的编译应最高遵守 IEEE-754 标准，以防止与浮点运算的书面顺序出现意外偏差。就我而言，我使用了 MSVC 2019 的 /fp:strict 命令行开关。
有数百亿个随机测试向量，我的测试程序报告的最大误差为 0.500000179 ulp。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
/* compute square root of sum of two positive floating-point numbers */
float sqrt_sum_pos (float a, float b)
{
float mn, mx, res, scale_in, scale_out;
float r, s, t, u, v, w, x;
/* sort arguments according to magnitude */
mx = a < b ? b : a;
mn = a < b ? a : b;
/* select scale factor: scale argument larger in magnitude towards unity */
scale_in = (mx > 1.0f) ? 0x1.0p-64f : 0x1.0p+64f;
scale_out = (mx > 1.0f) ? 0x1.0p+32f : 0x1.0p-32f;
/* scale input arguments */
mn = mn * scale_in;
mx = mx * scale_in;
/* represent sum as a normalized pair s:t of 'float' */
s = mx + mn; // most significant bits
t = (mx - s) + mn; // least significant bits
/* compute square root of s:t. Based on Alan Karp and Peter Markstein,
"High Precision Division and Square Root", ACM TOMS, vol. 23, no. 4,
December 1997, pp. 561-589
*/
r = sqrtf (1.0f / s);
if (s == 0.0f) r = 0.0f;
x = r * s;
s = fmaf (x, -x, s);
r = 0.5f * r;
u = s + t;
v = (s - u) + t;
s = r * u;
t = fmaf (r, u, -s);
t = fmaf (r, v, t);
r = x + s;
s = (x - r) + s;
s = s + t;
t = r + s;
s = (r - t) + s;
/* Component sum of t:s represents square root with maximum error very close to 0.5 ulp */
w = s + t;
/* compensate scaling of source operands */
res = w * scale_out;
/* handle special cases: NaN, Inf */
t = a + b;
if (isinf (mx)) res = mx;
if (isnan (t)) res = t;
return res;
}
// George Marsaglia's KISS PRNG, period 2**123. Newsgroup sci.math, 21 Jan 1999
// Bug fix: Greg Rose, "KISS: A Bit Too Simple" http://eprint.iacr.org/2011/007
static uint32_t kiss_z=362436069, kiss_w=521288629;
static uint32_t kiss_jsr=123456789, kiss_jcong=380116160;
#define znew (kiss_z=36969*(kiss_z&65535)+(kiss_z>>16))
#define wnew (kiss_w=18000*(kiss_w&65535)+(kiss_w>>16))
#define MWC ((znew<<16)+wnew )
#define SHR3 (kiss_jsr^=(kiss_jsr<<13),kiss_jsr^=(kiss_jsr>>17), \
kiss_jsr^=(kiss_jsr<<5))
#define CONG (kiss_jcong=69069*kiss_jcong+1234567)
#define KISS ((MWC^CONG)+SHR3)
uint32_t float_as_uint32 (float a)
{
uint32_t r;
memcpy (&r, &a, sizeof r);
return r;
}
uint64_t double_as_uint64 (double a)
{
uint64_t r;
memcpy (&r, &a, sizeof r);
return r;
}
float uint32_as_float (uint32_t a)
{
float r;
memcpy (&r, &a, sizeof r);
return r;
}
double floatUlpErr (float res, double ref)
{
uint64_t i, j, err, refi;
int expoRef;
/* ulp error cannot be computed if either operand is NaN, infinity, zero */
if (isnan (res) || isnan (ref) || isinf (res) || isinf (ref) ||
(res == 0.0f) || (ref == 0.0f)) {
return 0.0;
}
/* Convert the float result to an "extended float". This is like a float
with 56 instead of 24 effective mantissa bits
*/
i = ((uint64_t) float_as_uint32 (res)) << 32;
/* Convert the double reference to an "extended float". If the reference is
>= 2^129, we need to clamp to the maximum "extended float". If reference
is < 2^-126, we need to denormalize because of float's limited exponent
range.
*/
refi = double_as_uint64 (ref);
expoRef = (int)(((refi >> 52) & 0x7ff) - 1023);
if (expoRef >= 129) {
j = 0x7fffffffffffffffULL;
} else if (expoRef < -126) {
j = ((refi << 11) | 0x8000000000000000ULL) >> 8;
j = j >> (-(expoRef + 126));
} else {
j = ((refi << 11) & 0x7fffffffffffffffULL) >> 8;
j = j | ((uint64_t)(expoRef + 127) << 55);
}
j = j | (refi & 0x8000000000000000ULL);
err = (i < j) ? (j - i) : (i - j);
return err / 4294967296.0;
}
int main (void)
{
float arga, argb, res, reff;
uint32_t argai, argbi, resi, refi, diff;
double ref, ulp, maxulp = 0;
unsigned long long int count = 0;
do {
/* random positive inputs */
argai = KISS & 0x7fffffff;
argbi = KISS & 0x7fffffff;
/* increase occurence of zero, infinity */
if ((argai & 0xffff) == 0x5555) argai = 0x00000000;
if ((argbi & 0xffff) == 0x3333) argbi = 0x00000000;
if ((argai & 0xffff) == 0xaaaa) argai = 0x7f800000;
if ((argbi & 0xffff) == 0xcccc) argbi = 0x7f800000;
arga = uint32_as_float (argai);
argb = uint32_as_float (argbi);
res = sqrt_sum_pos (arga, argb);
ref = sqrt ((double)arga + (double)argb);
reff = (float)ref;
ulp = floatUlpErr (res, ref);
resi = float_as_uint32 (res);
refi = float_as_uint32 (reff);
diff = (refi > resi) ? (refi - resi) : (resi - refi);
if (diff > 1) {
/* if both source operands were NaNs, result could be either NaN,
quietened if necessary
*/
if (!(isnan (arga) && isnan (argb) &&
((resi == (argai | 0x00400000)) ||
(resi == (argbi | 0x00400000))))) {
printf ("\rerror: refi=%08x resi=%08x a=% 15.8e %08x b=% 15.8e %08x\n",
refi, resi, arga, argai, argb, argbi);
return EXIT_FAILURE;
}
}
if (ulp > maxulp) {
printf ("\rulp = %.9f @ a=%14.8e (%15.6a) b=%14.8e (%15.6a) a+b=%22.13a res=%15.6a ref=%22.13a\n",
ulp, arga, arga, argb, argb, (double)arga + argb, res, ref);
maxulp = ulp;
}
count++;
if (!(count & 0xffffff)) printf ("\r%llu", count);
} while (1);
printf ("\ntest passed\n");
return EXIT_SUCCESS;
}
另一种方法，现在@EricPostpischil 和@njuffa 强调了实际问题（即双舍入）。
(注意：下面是在谈论“行为良好”的数字。它不考虑精度边界或次正规，尽管它可以扩展到这样做。 )
首先，请注意 sqrt(x) 和 a+b 都保证 return 最接近结果的可表示值。问题是双舍入。也就是说，当我们 想要 计算 round(sqrt(a+b)) 时，我们本质上是在计算 round(sqrt(round(a+b)))。注意缺少内圆。
那么，内轮对结果的影响有多大？好吧，内圆加起来就是加法结果的 ±0.5 ULP。所以我们大致有 sqrt((a+b)*(1 ±2**-p))，假设有 p 位尾数。
这减少到 sqrt(a+b)*sqrt(1 ±2**-p)...但是 sqrt(1 ±2**-p) 比 (1 ±2**-p) 更接近 1！ （它很接近，但不完全是 (1 ±2**-(p+1))，因为这是一个有限差分。您可以从 1 (d/dx = 1/2) 附近的泰勒级数中看到这一点。）第二个舍入会影响另一个±0.5ULP的结果。
这意味着我们保证与“真实”结果相差不超过 1 ULP。因此，仅在 {sqrt(a+b)-1ULP, sqrt(a+b), sqrt(a+b)+1ULP} 之间进行选择的修正是一个可行的策略，如果我们可以“只是”弄清楚如何选择...
那么让我们看看我们是否可以想出一种基于比较的方法，该方法可以在有限精度下工作。 （注：以下如无特殊说明，均为无限精度）
resy = float(sqrt(a+b))
resx = resy.prev_nearest()
resz = resy.next_nearest()
注意 resx < resy < resz.
假设我们的浮点数有 p 位精度，那就变成
res = sqrt(a+b) // in infinite precision
resy = float(res)
resx = resy * (1 - 2**(1-p))
resz = resy * (1 + 2**(1-p))
所以让我们比较一下 resx 和 resy：
distx = abs(resx - res)
disty = abs(resy - res)
checkxy: distx < disty
checkxy: abs(resx - res) < abs(resy - res)
checkxy: (resx - res)**2 < (resy - res)**2
checkxy: resx**2 - 2*resx*res - res**2 < resy**2 - 2*resy*res - res**2
checkxy: resx**2 - resy**2 < 2*resx*res - 2*resy*res
checkxy: resx**2 - resy**2 < 2*res*(resx - resy)
// Assuming resx < resy
checkxy: resx+resy > 2*res
checkxy: resx+resy > 2*sqrt(a+b)
// Assuming resx+resy >= 0
checkxy: (resx+resy)**2 > 4*(a+b)
checkxy: (resy*(2 - 2**(1-p)))**2 > 4*(a+b)
checkxy: (resy**2)*((2 - 2**(1-p)))**2 > 4*(a+b)
checkxy: (resy**2)*(4 - 2*2**(1-p) + 2**(2-2p)) > 4*(a+b)
checkxy: (resy**2)*(4 - 4*2**(0-p) + 4*2**(0-2p)) > 4*(a+b)
checkxy: (resy**2)*(1 - 2**-p + 2**-2p) > a+b
...这是我们实际上可以以有限精度进行的检查（尽管它仍然需要 更高 精度，这很烦人）。
同上，对于 checkyz 我们得到
checkxy: disty < distz
checkyz: (resy**2)*(1 + 2**-p + 2**-2p) < a+b
通过这两项检查，您可以select 得到正确的结果。 ...然后“只是”检查/处理我在上面忽略的边缘情况。
现在，在实践中，与一开始只以更高的精度执行 sqrt 相比，我认为这不值得，至少除非有人能想出更好的选择方法。但它仍然是一个有趣的选择。
举个极端的例子。让我们有 u = 2^-p 其中 p 是浮点精度。
我们有 (1+u)^2 = (1+2u) + u^2.
如果我们取 a = 1+2u，我们有 float(a)=a，a 是一个可表示的浮点数（它是 1 之后的下一个浮点数），并且 b= u^2、float(b)=b , b 也可以表示为浮点数（作为 2^(-2p) 的幂）。
精确的sqrt(a+b)是(1+u)，应该四舍五入到float(1+u)=1，由于精确并列，它向下四舍五入到最接近的偶数尾数...
float(a+b)=a 和 float(sqrt(a))=1，所以没关系。
但是让我们改变b的最后一位：b=(1+2*u)*u^2； float(b)=b, b只是按比例缩小了两倍的精度。
我们现在有了确切的 sqrt(a+b) > 1+u，因此它应该四舍五入为 float(sqrt(a+b)) = 1+2u。
我们看到一点到 2^(-3p+1) 位（浮点精度的三倍）可以改变正确的舍入！
这意味着您应该不依赖双精度来执行正确的舍入运算。
有没有好的方法来计算
的正确舍入结果sqrt(a+b)
对于浮点数 a 和 b（相同精度），其中 0<=a<+inf 和 0<=b<+inf?
特别是对于 a+b 的计算会溢出的输入值？
(此处“正确舍入”的含义与 sqrt 本身的计算相同，即返回最接近以无限精度计算的“真实”结果的可表示值。)
（注意：一种明显的方法是以更大的浮点大小进行计算并避免溢出。不幸的是，这通常不起作用（例如，如果不支持更大的浮点格式))
我试过Herbie，但它完全放弃了。它似乎没有对 a+b 溢出的足够点进行采样以检测问题，而且似乎也没有很好地处理相关采样。不幸的是，它通常是一个很棒的工具。
到目前为止我一直在做的是（伪代码）
if a + b would overflow:
2*sqrt(a/4 + b/4) # Cannot overflow for finite inputs, as f::MAX/4 + f::MAX/4 <= f::MAX
else:
... # handle non-overflow case. Also interesting; not quite the topic of this question.
...这似乎主要在实践中起作用，但是 a) 完全没有原则和 b) 在实践中偶尔 returns 结果在溢出避免部分被 epsilon 关闭（例如在真正的结果是 x + 0.2(x.next_larger()-x) 但这个 returns x.next_larger() 而不是 x)
有关 f32 中“off-by-epsilon”问题的快速示例：
>>> import decimal
>>> decimal.getcontext().prec = 256
>>> from decimal import Decimal as D
>>> from numpy import float32 as f32
>>> a = D(f32("6.0847234e31").astype(float))
>>> b = D(f32("3.4028235e38").astype(float))
>>> res_act = (a+b).sqrt()
>>> res_calc = D(f32("1.8446744e19").astype(float)) # 2*sqrt(a/4 + b/4) in f32 precision
>>> res_best = D(f32("1.8446746e19").astype(float)) # obtained by brute-force
>>> abs(res_calc - res_act) > abs(res_best - res_act)
True # oops
（您必须相信我对 f32 计算结果的承诺，因为 Python 通常以 f64 精度运行。这也是 f32 舞蹈的原因。）
通过适当地按 2 的幂缩放可以很容易地避免溢出，这样量级大的参数被缩放到统一。困难的部分是产生正确的舍入结果。由于双舍入的潜在问题，我什至不完全相信在下一个更大的 IEEE-754 二进制浮点类型中执行中间计算可以保证。
在没有更广泛的浮点类型的情况下，人们将不得不退回到将多个本机精度数字链接在一起以执行具有更高中间精度的操作。 Dekker 提出的一种常见方案称为对精度。它使用成对的浮点数，其中较重要的部分通常称为“头”，不太重要的部分称为“尾”。将这两个部分归一化，使尾部的大小最多为头部大小的半个 ulp。
此方案中的有效有效位数为 2*p+1，其中 p 是基础浮点类型中的有效位数。 “额外”位由尾部的符号位表示。重要的是要注意，与底层基本类型相比，指数范围没有变化，因此我们需要相当积极地向统一性扩展，以避免在中间计算中遇到次正规操作数。对精度计算不能保证正确舍入的结果。使用三胞胎可能会奏效，但需要付出更多的努力，我无法负担得起答案。
但是，对精度可以提供忠实四舍五入且几乎总是正确四舍五入的结果。当 FMA（融合乘加）可用时，可以相当有效地构建基于 Newton-Raphson 的对精度平方根，产生大约 2*p-1 个好位。这就是我在下面的示例性 IS0-C99 代码中使用的，它使用 float 映射到 IEEE-754 binary32 作为本机浮点类型。成对精度代码的编译应最高遵守 IEEE-754 标准，以防止与浮点运算的书面顺序出现意外偏差。就我而言，我使用了 MSVC 2019 的 /fp:strict 命令行开关。
有数百亿个随机测试向量，我的测试程序报告的最大误差为 0.500000179 ulp。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
/* compute square root of sum of two positive floating-point numbers */
float sqrt_sum_pos (float a, float b)
{
float mn, mx, res, scale_in, scale_out;
float r, s, t, u, v, w, x;
/* sort arguments according to magnitude */
mx = a < b ? b : a;
mn = a < b ? a : b;
/* select scale factor: scale argument larger in magnitude towards unity */
scale_in = (mx > 1.0f) ? 0x1.0p-64f : 0x1.0p+64f;
scale_out = (mx > 1.0f) ? 0x1.0p+32f : 0x1.0p-32f;
/* scale input arguments */
mn = mn * scale_in;
mx = mx * scale_in;
/* represent sum as a normalized pair s:t of 'float' */
s = mx + mn; // most significant bits
t = (mx - s) + mn; // least significant bits
/* compute square root of s:t. Based on Alan Karp and Peter Markstein,
"High Precision Division and Square Root", ACM TOMS, vol. 23, no. 4,
December 1997, pp. 561-589
*/
r = sqrtf (1.0f / s);
if (s == 0.0f) r = 0.0f;
x = r * s;
s = fmaf (x, -x, s);
r = 0.5f * r;
u = s + t;
v = (s - u) + t;
s = r * u;
t = fmaf (r, u, -s);
t = fmaf (r, v, t);
r = x + s;
s = (x - r) + s;
s = s + t;
t = r + s;
s = (r - t) + s;
/* Component sum of t:s represents square root with maximum error very close to 0.5 ulp */
w = s + t;
/* compensate scaling of source operands */
res = w * scale_out;
/* handle special cases: NaN, Inf */
t = a + b;
if (isinf (mx)) res = mx;
if (isnan (t)) res = t;
return res;
}
// George Marsaglia's KISS PRNG, period 2**123. Newsgroup sci.math, 21 Jan 1999
// Bug fix: Greg Rose, "KISS: A Bit Too Simple" http://eprint.iacr.org/2011/007
static uint32_t kiss_z=362436069, kiss_w=521288629;
static uint32_t kiss_jsr=123456789, kiss_jcong=380116160;
#define znew (kiss_z=36969*(kiss_z&65535)+(kiss_z>>16))
#define wnew (kiss_w=18000*(kiss_w&65535)+(kiss_w>>16))
#define MWC ((znew<<16)+wnew )
#define SHR3 (kiss_jsr^=(kiss_jsr<<13),kiss_jsr^=(kiss_jsr>>17), \
kiss_jsr^=(kiss_jsr<<5))
#define CONG (kiss_jcong=69069*kiss_jcong+1234567)
#define KISS ((MWC^CONG)+SHR3)
uint32_t float_as_uint32 (float a)
{
uint32_t r;
memcpy (&r, &a, sizeof r);
return r;
}
uint64_t double_as_uint64 (double a)
{
uint64_t r;
memcpy (&r, &a, sizeof r);
return r;
}
float uint32_as_float (uint32_t a)
{
float r;
memcpy (&r, &a, sizeof r);
return r;
}
double floatUlpErr (float res, double ref)
{
uint64_t i, j, err, refi;
int expoRef;
/* ulp error cannot be computed if either operand is NaN, infinity, zero */
if (isnan (res) || isnan (ref) || isinf (res) || isinf (ref) ||
(res == 0.0f) || (ref == 0.0f)) {
return 0.0;
}
/* Convert the float result to an "extended float". This is like a float
with 56 instead of 24 effective mantissa bits
*/
i = ((uint64_t) float_as_uint32 (res)) << 32;
/* Convert the double reference to an "extended float". If the reference is
>= 2^129, we need to clamp to the maximum "extended float". If reference
is < 2^-126, we need to denormalize because of float's limited exponent
range.
*/
refi = double_as_uint64 (ref);
expoRef = (int)(((refi >> 52) & 0x7ff) - 1023);
if (expoRef >= 129) {
j = 0x7fffffffffffffffULL;
} else if (expoRef < -126) {
j = ((refi << 11) | 0x8000000000000000ULL) >> 8;
j = j >> (-(expoRef + 126));
} else {
j = ((refi << 11) & 0x7fffffffffffffffULL) >> 8;
j = j | ((uint64_t)(expoRef + 127) << 55);
}
j = j | (refi & 0x8000000000000000ULL);
err = (i < j) ? (j - i) : (i - j);
return err / 4294967296.0;
}
int main (void)
{
float arga, argb, res, reff;
uint32_t argai, argbi, resi, refi, diff;
double ref, ulp, maxulp = 0;
unsigned long long int count = 0;
do {
/* random positive inputs */
argai = KISS & 0x7fffffff;
argbi = KISS & 0x7fffffff;
/* increase occurence of zero, infinity */
if ((argai & 0xffff) == 0x5555) argai = 0x00000000;
if ((argbi & 0xffff) == 0x3333) argbi = 0x00000000;
if ((argai & 0xffff) == 0xaaaa) argai = 0x7f800000;
if ((argbi & 0xffff) == 0xcccc) argbi = 0x7f800000;
arga = uint32_as_float (argai);
argb = uint32_as_float (argbi);
res = sqrt_sum_pos (arga, argb);
ref = sqrt ((double)arga + (double)argb);
reff = (float)ref;
ulp = floatUlpErr (res, ref);
resi = float_as_uint32 (res);
refi = float_as_uint32 (reff);
diff = (refi > resi) ? (refi - resi) : (resi - refi);
if (diff > 1) {
/* if both source operands were NaNs, result could be either NaN,
quietened if necessary
*/
if (!(isnan (arga) && isnan (argb) &&
((resi == (argai | 0x00400000)) ||
(resi == (argbi | 0x00400000))))) {
printf ("\rerror: refi=%08x resi=%08x a=% 15.8e %08x b=% 15.8e %08x\n",
refi, resi, arga, argai, argb, argbi);
return EXIT_FAILURE;
}
}
if (ulp > maxulp) {
printf ("\rulp = %.9f @ a=%14.8e (%15.6a) b=%14.8e (%15.6a) a+b=%22.13a res=%15.6a ref=%22.13a\n",
ulp, arga, arga, argb, argb, (double)arga + argb, res, ref);
maxulp = ulp;
}
count++;
if (!(count & 0xffffff)) printf ("\r%llu", count);
} while (1);
printf ("\ntest passed\n");
return EXIT_SUCCESS;
}
另一种方法，现在@EricPostpischil 和@njuffa 强调了实际问题（即双舍入）。
(注意：下面是在谈论“行为良好”的数字。它不考虑精度边界或次正规，尽管它可以扩展到这样做。 )
首先，请注意 sqrt(x) 和 a+b 都保证 return 最接近结果的可表示值。问题是双舍入。也就是说，当我们 想要 计算 round(sqrt(a+b)) 时，我们本质上是在计算 round(sqrt(round(a+b)))。注意缺少内圆。
那么，内轮对结果的影响有多大？好吧，内圆加起来就是加法结果的 ±0.5 ULP。所以我们大致有 sqrt((a+b)*(1 ±2**-p))，假设有 p 位尾数。
这减少到 sqrt(a+b)*sqrt(1 ±2**-p)...但是 sqrt(1 ±2**-p) 比 (1 ±2**-p) 更接近 1！ （它很接近，但不完全是 (1 ±2**-(p+1))，因为这是一个有限差分。您可以从 1 (d/dx = 1/2) 附近的泰勒级数中看到这一点。）第二个舍入会影响另一个±0.5ULP的结果。
这意味着我们保证与“真实”结果相差不超过 1 ULP。因此，仅在 {sqrt(a+b)-1ULP, sqrt(a+b), sqrt(a+b)+1ULP} 之间进行选择的修正是一个可行的策略，如果我们可以“只是”弄清楚如何选择...
那么让我们看看我们是否可以想出一种基于比较的方法，该方法可以在有限精度下工作。 （注：以下如无特殊说明，均为无限精度）
resy = float(sqrt(a+b))
resx = resy.prev_nearest()
resz = resy.next_nearest()
注意 resx < resy < resz.
假设我们的浮点数有 p 位精度，那就变成
res = sqrt(a+b) // in infinite precision
resy = float(res)
resx = resy * (1 - 2**(1-p))
resz = resy * (1 + 2**(1-p))
所以让我们比较一下 resx 和 resy：
distx = abs(resx - res)
disty = abs(resy - res)
checkxy: distx < disty
checkxy: abs(resx - res) < abs(resy - res)
checkxy: (resx - res)**2 < (resy - res)**2
checkxy: resx**2 - 2*resx*res - res**2 < resy**2 - 2*resy*res - res**2
checkxy: resx**2 - resy**2 < 2*resx*res - 2*resy*res
checkxy: resx**2 - resy**2 < 2*res*(resx - resy)
// Assuming resx < resy
checkxy: resx+resy > 2*res
checkxy: resx+resy > 2*sqrt(a+b)
// Assuming resx+resy >= 0
checkxy: (resx+resy)**2 > 4*(a+b)
checkxy: (resy*(2 - 2**(1-p)))**2 > 4*(a+b)
checkxy: (resy**2)*((2 - 2**(1-p)))**2 > 4*(a+b)
checkxy: (resy**2)*(4 - 2*2**(1-p) + 2**(2-2p)) > 4*(a+b)
checkxy: (resy**2)*(4 - 4*2**(0-p) + 4*2**(0-2p)) > 4*(a+b)
checkxy: (resy**2)*(1 - 2**-p + 2**-2p) > a+b
...这是我们实际上可以以有限精度进行的检查（尽管它仍然需要 更高 精度，这很烦人）。
同上，对于 checkyz 我们得到
checkxy: disty < distz
checkyz: (resy**2)*(1 + 2**-p + 2**-2p) < a+b
通过这两项检查，您可以select 得到正确的结果。 ...然后“只是”检查/处理我在上面忽略的边缘情况。
现在，在实践中，与一开始只以更高的精度执行 sqrt 相比，我认为这不值得，至少除非有人能想出更好的选择方法。但它仍然是一个有趣的选择。
举个极端的例子。让我们有 u = 2^-p 其中 p 是浮点精度。
我们有 (1+u)^2 = (1+2u) + u^2.
如果我们取 a = 1+2u，我们有 float(a)=a，a 是一个可表示的浮点数（它是 1 之后的下一个浮点数），并且 b= u^2、float(b)=b , b 也可以表示为浮点数（作为 2^(-2p) 的幂）。
精确的sqrt(a+b)是(1+u)，应该四舍五入到float(1+u)=1，由于精确并列，它向下四舍五入到最接近的偶数尾数...
float(a+b)=a 和 float(sqrt(a))=1，所以没关系。
但是让我们改变b的最后一位：b=(1+2*u)*u^2； float(b)=b, b只是按比例缩小了两倍的精度。
我们现在有了确切的 sqrt(a+b) > 1+u，因此它应该四舍五入为 float(sqrt(a+b)) = 1+2u。
我们看到一点到 2^(-3p+1) 位（浮点精度的三倍）可以改变正确的舍入！
这意味着您应该不依赖双精度来执行正确的舍入运算。