hypot()的伴奏

首先要做的是分解:

b = sqrt(h*h - a*a) = sqrt((h-a)*(h+a))

我们不仅避免了一些溢出,而且获得了准确性。

如果任何因素接近 1E+154 = sqrt(1E+308) (最大使用IEEE 754 64位浮点)然后我们还必须避免溢出:

sqrt((h-a)*(h+a)) = sqrt(h-a) * sqrt(h+a)

这种情况不太可能发生,所以 sqrt 是合理的,即使它比 sqrt公司 。

请注意,如果 h ~ 5E+7 * a 然后 h ~ b 这意味着没有足够的数字来表示 b 不同于 h 。

这个答案假设一个平台使用符合IEEE-754(2008)的浮点算法,并提供融合乘加(FMA)功能。x86-64、ARM64和Power等常见体系结构都满足这两个条件。FMA在ISO C99和更高版本的C标准中作为标准数学函数公开 fma() 。在不提供FMA指令的硬件上,这需要仿真,仿真速度可能较慢 functionally deficient 。

从数学上讲,给定斜边和另一条腿的长度,直角三角形中一条腿(cathetus)的长度可以简单地计算为 √(h²-a²) 哪里 h 是斜边的长度。但是,当使用有限精度浮点算法进行计算时,我们面临两个问题:计算平方时可能会出现上溢或下溢到零的情况,而平方的减法会导致 subtractive cancellation 当正方形的大小相似时。

第一个问题可以通过按2的比例轻松解决 ⁿ 这样,量值较大的一词就更接近于统一。由于可能涉及次正常数,这无法通过操纵指数字段来实现,因为可能需要规范化/非规范化。但我们可以通过指数字段位操作,即乘以因子来计算所需的比例因子。我们知道,对于非特殊情况,斜边必须更长或与给定腿的长度相同,因此可以根据该参数进行缩放。

处理减法对消比较困难,但幸运的是,在其他重要问题中,计算与我们的计算h-a非常相似。例如,浮点计算大师研究了二次公式判别式的精确计算, b²-4ac :

William Kahan,“关于无超精密算法的浮点计算成本”,2004年11月21日( online )

最近,法国研究人员研究了两种产品差异的更一般情况, ad-bc :

Claude Pierre Jeannerod、Nicolas Louvet、Jean-Michel Muller,“进一步分析Kahan算法,以精确计算2 x 2行列式。” 计算数学 ,第82卷,第2842013年10月,第2245-2264页( online )

第二篇文章中基于FMA的算法计算两个乘积的差值,经证明最大误差为1.5 ulp 。有了这个构建块,我们就可以得到下面cathetus计算的简单ISO C99实现。通过与任意精度库的结果进行比较,在10亿个随机试验中观察到最大误差为1.2 ulp:

#include <stdint.h>
#include <string.h>
#include <float.h>
#include <math.h>

uint64_t __double_as_uint64 (double a)
{
    uint64_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

double __uint64_as_double (uint64_t a)
{
    double r;
    memcpy (&r, &a, sizeof r);
    return r;
}

/*
  diff_of_products() computes a*b-c*d with a maximum error < 1.5 ulp

  Claude-Pierre Jeannerod, Nicolas Louvet, and Jean-Michel Muller, 
  "Further Analysis of Kahan's Algorithm for the Accurate Computation 
  of 2x2 Determinants". Mathematics of Computation, Vol. 82, No. 284, 
  Oct. 2013, pp. 2245-2264
*/
double diff_of_products (double a, double b, double c, double d)
{
    double w = d * c;
    double e = fma (-d, c, w);
    double f = fma (a, b, -w);
    return f + e;
}

/* compute sqrt (h*h - a*a) accurately, avoiding spurious overflow */
double my_cathetus (double h, double a)
{
    double fh, fa, res, scale_in, scale_out, d, s;
    uint64_t expo;

    fh = fabs (h);
    fa = fabs (a);

    /* compute scale factors */
    expo = __double_as_uint64 (fh) & 0xff80000000000000ULL;
    scale_in = __uint64_as_double (0x7fc0000000000000ULL - expo);
    scale_out = __uint64_as_double (expo + 0x0020000000000000ULL);

    /* scale fh towards unity */
    fh = fh * scale_in;
    fa = fa * scale_in;

    /* compute sqrt of difference of scaled arguments, avoiding overflow */
    d = diff_of_products (fh, fh, fa, fa);
    s = sqrt (d);

    /* reverse previous scaling */
    res = s * scale_out;

    /* handle special arguments */
    if (isnan (h) || isnan (a)) {
        res = h + a;
    }

    return res;
}

假设IEEE 754基本64位二进制浮点,我会考虑如下算法:

• 设置 s (比例)为2 ^–512 如果2 ¹⁰⁰ -137; 一 2. ^+512个 如果 一 <2. ¹⁰⁰ ,否则为1。
• 允许 一 '是 一 -1美分 s 和 b '是 b -1美分 s 。
• 计算sqrt( 一 '美元 一 ' b ' b ')/ s 。

关于推理的注释:

• 如果 一 大(或小),乘以 s 减少(或增加)值,以便 一 '保持在浮点范围内。
• 比例因子是2的幂,因此在二进制浮点中,乘和除它是精确的。
• b 必须小于(或等于) 一 ,否则返回NaN,这是适当的。如果我们正在增加 一 ,无错误发生; b '和 b '美元 b '保持在范围内。如果我们正在减少 一 ,则, b '可能会失去精度或变为零 b 很小,但是 b 比 一 计算结果不能取决于 b 无论如何。
• 我将浮点范围划分为三个区间,因为两个区间不够。例如,如果设置 s 为2 ⁵¹² 如果1￥ 一 和2 ^+512个 否则,1将缩放为2 ⁵¹² 然后平方为2 ¹⁰²⁴ ,此时a b 略低于1将失去与结果相关的精度。但如果你用较小的量级功率 s ,例如2 ⁵¹¹ ,然后是2 ¹⁰²³ 将缩放到2 ⁵¹² 平方为2 ¹⁰²⁴ ,这是不允许的。因此,我们需要不同的比例因子 一 =1和 一 =2 ¹⁰²³ 。同样地, 一 =2 ^-1049年 需要太大的比例因子 一 =1。所以需要三个。
• 除法是出了名的慢,所以人们可能想用一个准备好的 s ^第1页 而不是除以 s 。

hypot 它的特点是 极少数精选 C标准库函数 不 传播 NaN 哦!(另一个是 pow 对于第一个参数为1的情况。)

放在一边,我倾向于只写

returns sqrt(h * h - a * a); // h is the hypotenuse

作为函数的主体,并让调用者检查输入。如果无法做到这一点,请遵循 海波 准确地