溢出可以通过适当的二次幂缩放来轻松避免，使具有大量幅度的参数缩放到接近1。难点在于产生正确舍入的结果。我甚至不能完全确信在下一个更大的IEEE-754二进制浮点类型中执行中间计算会保证这一点，由于可能存在双重舍入问题。
在没有更广泛的浮点类型的情况下，必须退回到将多个本地精度数字链接在一起以执行具有更高中间精度的操作的方法。由Dekker提出的常见方案称为配对精度。它使用浮点数对，其中更重要的部分通常称为“头部”，不太重要的部分称为“尾部”。这两部分被归一化，使得尾部的大小最多为头部大小的半个ulp。
在此方案中，有效有效位数为2 * p + 1，其中p是基础浮点类型中的有效位数。 "额外"位由尾部的符号位表示。重要的是要注意，指数范围与底层基本类型相比不变，因此我们需要积极地朝向1进行缩放，以避免在中间计算中遇到次正常操作数。 配对精度计算无法保证正确舍入的结果。使用三元组可能有效，但需要更多的工作量，我无法投入答案。
然而，配对精度可以提供忠实舍入并几乎总是正确舍入的结果。当FMA（合并乘加）可用时，可以相当有效地构建基于牛顿-拉弗森的配对精度平方根，产生约2 * p-1个好位。这是我在下面的示例ISO-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>
float sqrt_sum_pos (float a, float b)
{
float mn, mx, res, scale_in, scale_out;
float r, s, t, u, v, w, x;
mx = a < b ? b : a;
mn = a < b ? a : b;
scale_in = (mx > 1.0f) ? 0x1.0p-64f : 0x1.0p+64f;
scale_out = (mx > 1.0f) ? 0x1.0p+32f : 0x1.0p-32f;
mn = mn * scale_in;
mx = mx * scale_in;
s = mx + mn;
t = (mx - s) + mn;
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;
w = s + t;
res = w * scale_out;
t = a + b;
if (isinf (mx)) res = mx;
if (isnan (t)) res = t;
return res;
}
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;
if (isnan (res) || isnan (ref) || isinf (res) || isinf (ref) ||
(res == 0.0f) || (ref == 0.0f)) {
return 0.0;
}
i = ((uint64_t) float_as_uint32 (res)) << 32;
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 {
argai = KISS & 0x7fffffff;
argbi = KISS & 0x7fffffff;
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 (!(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;
}
sqrt(a)接近舍入上限时的情况，而添加一个非常小的b会使“精确”的sqrt(a+b)舍入到更大的下一个值？你是否有一个可用的精确融合乘加操作？ - undefineda+b会溢出，那么将每个数除以浮点格式的基数的平方（对于二进制来说是4），然后取平方根，并乘以基数。这样做不会引入新的舍入误差，尽管仍然存在一个问题，即由于存在两个舍入误差，sqrt(a+b)可能与理想计算结果不同。 - undefined2*sqrt(a/4 + b/4)这部分，而是在于sqrt(a+b)本身。 - undefined