floating-point - 正确舍入两个处理溢出的 float 之和的 sqrt 计算

Question

是否有一个好的方法来计算正确舍入的结果

sqrt(a+b)

对于 float 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) 在实践中偶尔会返回一个结果，该结果在溢出避免部分中被 epsilon 关闭(例如，真实结果是 x + 0.2(x.next_larger()-x)但这返回 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 计算结果的 promise ，因为 Python 通常以 f64 精度运行。这也是 f32 舞蹈的原因。)

Answer 1

通过适当缩放 2 的幂可以很容易地避免溢出，这样量级大的参数就可以统一缩放。困难的部分是产生正确的舍入结果。由于双舍入的潜在问题，我什至不完全相信在下一个更大的 IEEE-754 二进制浮点类型中执行中间计算可以保证这一点。

在没有更广泛的浮点类型的情况下，人们将不得不回退到将多个 native 精度数字链接在一起以执行具有更高中间精度的操作。 Dekker 提出的一种常见方案称为对精度。它使用成对的 float ，其中较重要的部分通常称为“头”，不太重要的部分称为“尾”。对这两个部分进行归一化处理，使尾部的大小最多为头部大小的一半 ulp。

此方案中的有效有效位数为 2*p+1，其中 p 是基础浮点类型中的有效位数。 “额外”位由尾部的符号位表示。重要的是要注意，与底层基本类型相比，指数范围没有变化，因此我们需要相当积极地向统一性扩展，以避免在中间计算中遇到次正规操作数。对精度计算不能保证正确舍入的结果。使用三胞胎可能会奏效，但需要付出更多的努力，我无法负担得起答案。

但是，对精度可以提供忠实四舍五入且几乎总是正确四舍五入的结果。当 FMA(融合乘加)可用时，可以相当有效地构建基于 Newton-Raphson 的对精度平方根，产生大约 2*p-1 个好位。这就是我在下面的示例性 IS0-C99 代码中使用的，它使用映射到 IEEE-754 binary32 的 float 作为 native 浮点类型。成对精度代码的编译应最高遵守 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;
}