Leo*_*nid 35 c++ java hypotenuse
我用C++ hypot()和做了一些测试Java Math.hypot.它们似乎都明显慢于sqrt(a*a + b*b).那是因为精度更高吗？计算斜边hypot函数的方法有哪些？令人惊讶的是,我在文档中找不到任何性能不佳的迹象.
rkg*_*rkg 27
这不是一个简单的sqrt函数.您应该检查此链接以获取算法的实现:http://www.koders.com/c/fid7D3C8841ADC384A5F8DE0D081C88331E3909BF3A.aspx
它有while循环来检查收敛
/* Slower but safer algorithm due to Moler and Morrison. Never
produces any intermediate result greater than roughly the
larger of X and Y. Should converge to machine-precision
accuracy in 3 iterations. */
double r = ratio*ratio, t, s, p = abig, q = asmall;
do {
t = 4. + r;
if (t == 4.)
break;
s = r / t;
p += 2. * s * p;
q *= s;
r = (q / p) * (q / p);
} while (1);
编辑(JM更新):
这是最初的Moler-Morrison论文,由于Dubrulle ,这是一个很好的后续行动.
小智 5
这是一个更快的实现，其结果也更接近 java.lang.Math.hypot: (注意：对于 Delorie 的实现，需要添加对 NaN 和 +-Infinity 输入的处理)
private static final double TWO_POW_450 = Double.longBitsToDouble(0x5C10000000000000L);
private static final double TWO_POW_N450 = Double.longBitsToDouble(0x23D0000000000000L);
private static final double TWO_POW_750 = Double.longBitsToDouble(0x6ED0000000000000L);
private static final double TWO_POW_N750 = Double.longBitsToDouble(0x1110000000000000L);
public static double hypot(double x, double y) {
x = Math.abs(x);
y = Math.abs(y);
if (y < x) {
double a = x;
x = y;
y = a;
} else if (!(y >= x)) { // Testing if we have some NaN.
if ((x == Double.POSITIVE_INFINITY) || (y == Double.POSITIVE_INFINITY)) {
return Double.POSITIVE_INFINITY;
} else {
return Double.NaN;
}
}
if (y-x == y) { // x too small to substract from y
return y;
} else {
double factor;
if (x > TWO_POW_450) { // 2^450 < x < y
x *= TWO_POW_N750;
y *= TWO_POW_N750;
factor = TWO_POW_750;
} else if (y < TWO_POW_N450) { // x < y < 2^-450
x *= TWO_POW_750;
y *= TWO_POW_750;
factor = TWO_POW_N750;
} else {
factor = 1.0;
}
return factor * Math.sqrt(x*x+y*y);
}
}
小智 5
我发现 Math.hypot() 非常慢。我发现我可以使用相同的算法编写一个快速的 java 版本，并产生相同的结果。可以用这个来代替。
/**
* <b>hypot</b>
* @param x
* @param y
* @return sqrt(x*x +y*y) without intermediate overflow or underflow.
* @Note {@link Math#hypot} is unnecessarily slow. This returns the identical result to
* Math.hypot with reasonable run times (~40 nsec vs. 800 nsec).
* <p>The logic for computing z is copied from "Freely Distributable Math Library"
* fdlibm's e_hypot.c. This minimizes rounding error to provide 1 ulb accuracy.
*/
public static double hypot(double x, double y) {
if (Double.isInfinite(x) || Double.isInfinite(y)) return Double.POSITIVE_INFINITY;
if (Double.isNaN(x) || Double.isNaN(y)) return Double.NaN;
x = Math.abs(x);
y = Math.abs(y);
if (x < y) {
double d = x;
x = y;
y = d;
}
int xi = Math.getExponent(x);
int yi = Math.getExponent(y);
if (xi > yi + 27) return x;
int bias = 0;
if (xi > 510 || xi < -511) {
bias = xi;
x = Math.scalb(x, -bias);
y = Math.scalb(y, -bias);
}
// translated from "Freely Distributable Math Library" e_hypot.c to minimize rounding errors
double z = 0;
if (x > 2*y) {
double x1 = Double.longBitsToDouble(Double.doubleToLongBits(x) & 0xffffffff00000000L);
double x2 = x - x1;
z = Math.sqrt(x1*x1 + (y*y + x2*(x+x1)));
} else {
double t = 2 * x;
double t1 = Double.longBitsToDouble(Double.doubleToLongBits(t) & 0xffffffff00000000L);
double t2 = t - t1;
double y1 = Double.longBitsToDouble(Double.doubleToLongBits(y) & 0xffffffff00000000L);
double y2 = y - y1;
double x_y = x - y;
z = Math.sqrt(t1*y1 + (x_y*x_y + (t1*y2 + t2*y))); // Note: 2*x*y <= x*x + y*y
}
if (bias == 0) {
return z;
} else {
return Math.scalb(z, bias);
}
}
hypot()与幼稚的实现相比，在中间计算中避免上溢和下溢会产生开销sqrt(a*a+b*b)。这涉及缩放操作。在较旧的实现中，缩放可能使用除法，这本身是一个相当慢的操作。即使在通过直接指数操作完成缩放的情况下，在较旧的处理器体系结构上也可能相当慢，因为在整数 ALU 和浮点单元之间传输数据相当慢（例如，涉及通过内存的往返）。数据驱动的分支对于各种扩展方法也很常见；这些很难预测。
数学库设计者的一个常见目标是为简单的数学函数实现忠实的舍入，例如hypot()，最大误差小于 1 ulp。与原始解决方案相比，准确性的这种改进意味着必须以高于本机精度的方式执行中间计算。一种经典方法是将操作数拆分为“高”和“低”部分并模拟扩展精度。这增加了拆分的开销并增加了浮点运算的数量。最后，ISO C99 规范hypot()（后来被 C++ 标准采用）规定了对 NaN 和无穷大的处理，这些处理自然不会从简单的计算中产生。
旧最佳实践的典型例子是__ieee754_hypot()在FDLIBM。虽然它声称最大误差为 1 ulp，但对任意精度库的快速测试表明它实际上并未实现该目标，因为可以观察到高达 1.15 ulp 的误差。
自从提出这个问题以来，处理器体系结构的两项进步已经实现了更高效的hypot()实现。一个是引入了融合乘加 ( FMA ) 操作，它在内部使用全乘积进行加法，并在最后只应用一次舍入。这通常减轻了在中间计算中模拟稍高的精度的需要，并且通过将两个操作合并为一条指令还可以提高性能。另一个架构进步是允许浮点和整数运算在同一组寄存器上运行，从而使浮点操作数的位操作便宜。
因此，在现代高性能数学库中hypot()，吞吐量通常只有 的一半左右sqrt()，例如，从英特尔发布的 MKL性能数据中可以看出。
对于自己的实验，最好从 的单精度实现开始hypot()，因为这允许对所有可能参数值的总组合的较大部分进行测试。这还允许在使用许多现代处理器架构提供的低精度倒数平方根近似值时探索速度和精度之间的权衡。一种这样的探索的示例性代码如下所示。
请注意，ISO C99（以及扩展名 C++）对 NaN 处理的要求fmin() fmax()并不总是与浮点最小/最大机器指令的功能相匹配，这使得这些功能比它们本来的速度要慢。由于下面的代码单独处理 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;
}
一个简单的基于 FMA 的双精度实现hypot()如下所示：
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;
}