伴侣hypot()

Question

hypot 函数在 1999 年的 C 语言修订版中引入，计算以其他边为参数的直角三角形的斜边，但要注意避免因幼稚而导致的上溢/下溢实现为

double hypot(double a, double b)
{
  return sqrt(a*a + b*b);
}

我发现自己需要配套功能：给定三角形的边和斜边，找到第三边（避免下溢/溢出）。我可以想到几种方法来做到这一点，但想知道是否存在现有的“最佳实践”？

我的目标是 Python，但实际上我正在寻找算法指针。

---

感谢您的回复。如果有人对结果感兴趣，可以找到我的 C99 实现 here 和 Python 版本 here，这是 Hypothesis 项目的一部分。

Answer 1

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

在数学上，直角三角形中一条腿（导管）的长度，给定斜边和另一条腿的长度，简单地计算为√(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，284，2013年10月，第2245-2264页（online）

第二篇论文中基于 FMA 的算法计算两个产品的差异，经证明的最大误差为 1.5 ulp。有了这个构建块，我们就可以直接实现下面的导管计算的 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;
}

Answer 2

首先要做的是分解：

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 不同。

Answer 3

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

• 如果 2¹⁰⁰ ≤ a，则将 s（用于比例）设置为 2⁻⁵¹²，2+512 如果 a -100，否则为 1。
• 设 a' 为 a•s 和 b' 为 b •s。
• 计算 sqrt(a'•a' - b'•b') / s.

推理说明：

• 如果 a 大（或小），乘以 s 会减小（或增大）值，以便 a' 的平方保持在浮点范围内。
• 比例因子是 2 的幂，因此乘以它和除以它是精确的二进制浮点数。
• b 必须小于（或等于）a，否则我们返回 NaN，这是合适的。在我们增加a的情况下，不会发生错误； b' 和 b'•b' 仍在范围内。在我们减少 a 的情况下，如果 b 很小，b' 可能会丢失精度或变为零，但随后 b 比 a 小得多，以至于计算结果在任何情况下都不能依赖于 b 的精确值。
• 我将浮点范围划分为三个区间，因为两个是不够的。例如，如果您将 s 设置为 2⁻⁵¹² 如果 1 ≤ a 则设置为 2⁺⁵¹² 否则，然后 1 将缩放到 2⁻⁵¹² 然后平方到 2⁻¹⁰²⁴，此时 a b 略低于 1 将失去相关的精度到结果。但是如果你对 s 使用较小的幂，例如 2⁻⁵¹¹，那么 2¹⁰²³ 将缩放到 2⁵¹² 并平方到 2¹⁰²⁴，超出范围。因此，对于 a = 1 和 a = 2¹⁰²³，我们需要不同的比例因子。同样，a = 2⁻¹⁰⁴⁹ 需要一个对于 a = 1 来说太大的比例因子。所以需要三个。
• 除法速度非常慢，因此可能需要乘以准备好的 s^-1 而不是除以 s。

Answer 4

hypot 的独特之处在于它是极少数不 传播NaN 的C 标准库函数之一！ （另一个是pow，用于第一个参数为 1 的情况。）

抛开这些，我倾向于只写

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

作为函数的主体，并让调用者检查输入。如果你不能这样做，那么请忠实地遵循hypot的规范。