hypot() 与幼稚的实现sqrt(a*a+b*b) 相比，在避免中间计算中的上溢和下溢方面会产生开销。这涉及缩放操作。在较旧的实现中，缩放可能使用除法，这本身就是一个相当慢的操作。即使缩放是通过直接指数操作完成的，在旧的处理器架构上它也可能相当慢，因为在整数 ALU 和浮点单元之间传输数据相当慢（例如，涉及到内存的往返）。数据驱动的分支对于各种扩展方法也很常见；这些很难预测。
数学库设计者的一个常见目标是实现对 hypot() 等简单数学函数的忠实舍入，即最大误差小于 1 ulp。与朴素解决方案相比，这种准确性的提高意味着必须以高于原生精度的方式执行中间计算。一种经典的方法是将操作数分成“高”和“低”部分并模拟扩展精度。这增加了拆分的开销并增加了浮点运算的数量。最后，hypot() 的 ISO C99 规范（后来被 C++ 标准采用）规定了对 NaN 和无穷大的处理，这些计算不会自然而然地脱离直接计算。
较早的最佳做法的代表性示例是 __ieee754_hypot() in FDLIBM。虽然它声称最大误差为 1 ulp，但针对任意精度库的快速测试表明它实际上并未实现该目标，因为可以观察到高达 1.15 ulp 的误差。
自从提出这个问题以来，处理器架构的两项进步使hypot() 实现更加高效。一个是引入了融合乘加 (FMA) 操作，它在内部使用完整的产品进行加法，并在最后只应用一次舍入。这通常减轻了在中间计算中模拟稍高精度的需求，并且通过将两个操作组合到一条指令中也可以提高性能。另一个架构上的进步是允许浮点和整数运算在同一组寄存器上进行操作，从而降低了浮点操作数的位操作成本。
因此，在现代高性能数学库中，hypot() 的吞吐量通常约为sqrt() 的一半，例如，英特尔为 MKL 发布的performance data 中可以看到。
对于自己的实验，最好从 hypot() 的单精度实现开始，因为这允许在所有可能的参数值的总组合中的较大部分进行测试。这也允许在使用许多现代处理器架构提供的低精度倒数平方根近似时探索速度和准确性之间的权衡。一种这样的探索的示例性代码如下所示。
请注意，ISO C99（以及扩展的 C++）对 fmin() fmax() 施加的 NaN 处理要求并不总是与浮点最小/最大机器指令的功能相匹配，从而使这些功能变慢比他们可能。由于下面的代码分别处理 NaN，因此应该可以直接使用此类机器指令。代码编译时应进行全面优化，但也要严格遵守 IEEE-754。
基于对 500 亿个随机参数对的测试，以下代码变体观察到的最大错误为：(USE_BEEBE == 1): 1.51 ulp; （USE_BEEBE == 0 && USE_SQRT == 1）：1.02 ulp； （USE_BEEBE == 0 && USE_SQRT == 0 && USE_2ND_ITERATION == 0）：2.77 ulp； (USE_BEEBE == 0 && USE_SQRT == 0 && USE_2ND_ITERATION == 1)：1.02 ulp。
#include <stdint.h>
#include <string.h>
#include <float.h>
#include <math.h>
#include "immintrin.h"
uint32_t __float_as_uint32 (float a);
float __uint32_as_float (uint32_t a);
float sse_rsqrtf (float a);
float my_hypotf (float a, float b)
{
float fa, fb, mn, mx, res, r, t;
fa = fabsf (a);
fb = fabsf (b);
mx = fmaxf (fa, fb);
mn = fminf (fa, fb);
#if USE_BEEBE
/* Nelson H. F. Beebe, "The Mathematical Function Computation Handbook."
Springer 2017
*/
float s, c;
r = mn / mx;
t = fmaf (r, r, 1.0f);
s = sqrtf (t);
c = fmaf (-s, s, t) / (s + s);
res = fmaf (mx, s, mx * c);
if (mx == 0) res = mx; // fixup hypot(0,0)
#else // USE_BEEBE
float scale_in, scale_out, s, v, w;
uint32_t expo;
/* compute scale factors */
expo = __float_as_uint32 (mx) & 0xfc000000;
scale_in = __uint32_as_float (0x7e000000 - expo);
scale_out = __uint32_as_float (expo + 0x01000000);
/* scale mx towards unity */
mn = mn * scale_in;
mx = mx * scale_in;
/* mx in [2**-23, 2**6) */
s = fmaf (mx, mx, mn * mn); // 0.75 ulp
#if USE_SQRT
w = sqrtf (s);
#else // USE_SQRT
r = sse_rsqrtf (s);
/* use A. Schoenhage's coupled iteration for the square root */
v = r * 0.5f;
w = r * s;
w = fmaf (fmaf (w, -w, s), v, w);
#if USE_2ND_ITERATION
v = fmaf (fmaf (r, -w, 1), v, v);
w = fmaf (fmaf (w, -w, s), v, w);
#endif // USE_2ND_ITERATION
if (mx == 0) w = mx; // fixup hypot(0,0)
#endif // USE_SQRT
/* reverse previous scaling */
res = w * scale_out;
#endif // USE_BEEBE
/* check for special cases: infinity and NaN */
t = a + b;
if (!(fabsf (t) <= INFINITY)) res = t; // isnan(t)
if (mx == INFINITY) res = mx; // isinf(mx)
return res;
}
uint32_t __float_as_uint32 (float a)
{
uint32_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;
}
float sse_rsqrtf (float a)
{
__m128 b, t;
float res;
int old_mxcsr;
old_mxcsr = _mm_getcsr();
_MM_SET_FLUSH_ZERO_MODE (_MM_FLUSH_ZERO_OFF);
_MM_SET_ROUNDING_MODE (_MM_ROUND_NEAREST);
b = _mm_set_ss (a);
t = _mm_rsqrt_ss (b);
_mm_store_ss (&res, t);
_mm_setcsr (old_mxcsr);
return res;
}
hypot() 的基于 FMA 的简单双精度实现如下所示：
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;
}
double my_hypot (double a, double b)
{
double fa, fb, mn, mx, scale_in, scale_out, r, s, t;
uint64_t expo;
fa = fabs (a);
fb = fabs (b);
mx = fmax (fa, fb);
mn = fmin (fa, fb);
/* compute scale factors */
expo = __double_as_uint64 (mx) & 0xff80000000000000ULL;
scale_in = __uint64_as_double (0x7fc0000000000000ULL - expo);
scale_out = __uint64_as_double (expo + 0x0020000000000000ULL);
/* scale mx towards unity */
mn = mn * scale_in;
mx = mx * scale_in;
/* mx in [2**-52, 2**6) */
s = fma (mx, mx, mn * mn); // 0.75 ulp
r = sqrt (s);
/* reverse previous scaling */
r = r * scale_out;
/* check for special cases: infinity and NaN */
t = a + b;
if (!(fabs (t) <= INFINITY)) r = t; // isnan(t)
if (mx == INFINITY) r = mx; // isinf(mx)
return r;
}