计算距离平方的最快方法答案

【问题标题】：Fastest way to compute distance squared计算距离平方的最快方法
【发布时间】：2011-12-26 01:15:16
【问题描述】：

我的代码在很大程度上依赖于计算 3D 空间中两点之间的距离。为了避免昂贵的平方根，我始终使用平方距离。但它仍然占用了大部分计算时间，我想用更快的东西替换我的简单函数。我现在有：

double distance_squared(double *a, double *b)
{
  double dx = a[0] - b[0];
  double dy = a[1] - b[1];
  double dz = a[2] - b[2];

  return dx*dx + dy*dy + dz*dz;
}

我也尝试使用宏来避免函数调用，但没有太大帮助。

#define DISTANCE_SQUARED(a, b) ((a)[0]-(b)[0])*((a)[0]-(b)[0]) + ((a)[1]-(b)[1])*((a)[1]-(b)[1]) + ((a)[2]-(b)[2])*((a)[2]-(b)[2])

我考虑过使用 SIMD 指令，但找不到一个好的示例或完整的指令列表（理想情况下是两个向量上的乘法+加法）。

GPU 不是一个选项，因为在每个函数调用中只知道一组点。

计算距离平方的最快方法是什么？

【问题讨论】：

如果您的点是静止的预计算（全部？）距离对可能是有益的。
您必须在其他地方进行优化。这应该已经是最佳的了（正如大卫所说）。也许您可以提供更广泛的问题视图..也许您不必重新计算所有内容或使用其他技巧..
您要整体计算什么。这显然是某些算法的底部，但是您是否试图在一组 n 或其他东西中找到两个点的闭合点。也许如果您描述您实际尝试计算的内容，可能会有替代算法而不是蛮力。
这些点不是静止的，因此无法进行预计算。我已经在考虑其他算法了。这不是蛮力代码（kdtree 实现）的一部分，但仍需要经常计算，以便对其进行优化。我只是想知道一般来说最快的方法是什么，因为我找不到有关该主题的任何信息。

标签： c optimization simd

【解决方案1】：

一个好的编译器将优化它，就像你将管理的一样。一个好的编译器会使用 SIMD 指令，如果它认为它们将是有益的。确保为编译器打开所有此类可能的优化。不幸的是，维度 3 的向量不太适合 SIMD 单元。

我怀疑您将不得不接受编译器生成的代码可能非常接近最优并且无法获得显着收益。

【讨论】：

完全正确。例如，使用 gcc 和 -O3 -march=native 这通常会生成还不错的代码。在您投入大量时间进行优化之前，请使用 -S 检查汇编程序。

【解决方案2】：

第一个显而易见的事情是使用 restrict 关键字。

就像现在一样，a 和 b 是可别名的（因此，从编译器的角度来看，它假定最坏的情况是已别名）。没有编译器会自动对其进行向量化，因为这样做是错误的。

更糟糕的是，编译器不仅不能向量化这样的循环，如果您还存储（幸运的是在您的示例中没有），它也必须每次重新加载值。始终清楚别名，因为它会极大地影响编译器。

接下来，如果您可以接受，请使用float 而不是double 并填充到 4 个浮点数，即使其中一个未使用，这对于大多数 CPU 来说是一种更“自然”的数据布局（这有点特定于平台，但对于大多数平台来说，4 个浮点数是一个很好的猜测——3 个双精度数，也就是“典型”CPU 上的 1.5 个 SIMD 寄存器，在任何地方都不是最佳的）。

（对于手写的 SIMD 实现（这比您想象的要难），首先要确保数据对齐。接下来，查看您的指令在目标机器上的延迟时间，并首先执行最长的延迟时间. 例如，在 Prescott Intel 上，首先将每个组件洗牌到寄存器中然后与自身相乘是有意义的，即使使用 3 个乘法而不是一个，因为洗牌有很长的延迟。在后来的模型上，洗牌需要一个周期，所以这将是一个完全的反优化。
这再次表明将其留给编译器并不是一个坏主意。）

【讨论】：

为什么别名会阻止矢量化？
我看不出restrict 在这里有什么帮助——据我所知，函数体中的任何内容都不会受益于额外的别名信息；如果您对代码注释很激进，则应将参数标记为const double *restrict，因为它们永远不会被修改...
向量化意味着一次读取和写入多个独立的值，并在一次操作中处理它们。如果没有其他指针/引用，或者如果它可以毫无疑问地证明在同一范围内没有其他指针/引用访问任何数据（有时这是可能的），或者如果程序员明确承诺（restrict关键字）他已经想到了这一点并且它不会发生。在任何其他情况下，它都是不安全的，编译器通常会拒绝执行可能给出不正确结果的操作。
正如 Christoph 正确指出的那样，该函数中的数据永远不会被修改，诚然，一个聪明的编译器可能确实会弄清楚并能够证明它是安全的向量化（如果没有其他方法阻止它），但我不会为此赌上我的右手。在任何情况下，声明一个没有别名为 restrict 的指针是正确的（声明一个常量指针 const 也是如此）。声明的精确不仅是表面上的，也是一种“代码正确性”。
无论如何，您都需要一个聪明的编译器才能从restrict 中受益。该关键字通过消除某些优化障碍来发挥作用，这些障碍首先在这里不存在。

【解决方案3】：

执行此操作的 SIMD 代码（使用 SSE3）：

movaps xmm0,a
movaps xmm1,b
subps xmm0,xmm1
mulps xmm0,xmm0
haddps xmm0,xmm0
haddps xmm0,xmm0

但您需要四个值向量 (x,y,z,0) 才能工作。如果您只有三个值，那么您需要做一些摆弄才能获得所需的格式，这将抵消上述任何好处。

但总的来说，由于 CPU 的超标量流水线架构，获得性能的最佳方法是对大量数据执行相同的操作，这样您就可以交错各个步骤并进行一些循环展开避免管道停顿。基于“修改后不能直接使用值”的原则，上面的代码肯定会停在最后三个指令上——第二条指令必须等待前一条指令的结果完成，这在流水线系统。

同时对两个或多个不同的点集进行计算可以消除上述瓶颈-在等待一个计算结果的同时，您可以开始下一个点的计算：

movaps xmm0,a1
                  movaps xmm2,a2
movaps xmm1,b1
                  movaps xmm3,b2
subps xmm0,xmm1
                  subps xmm2,xmm3
mulps xmm0,xmm0
                  mulps xmm2,xmm2
haddps xmm0,xmm0
                  haddps xmm2,xmm2
haddps xmm0,xmm0
                  haddps xmm2,xmm2

【讨论】：

假设每个#D 点的第四个值为零，您能否使用 SSE 添加 C 代码来执行第一行？谢谢。

【解决方案4】：

如果您想优化某些东西，首先要分析代码并检查汇编器输出。

在使用 gcc -O3 (4.6.1) 编译后，我们将使用 SIMD 得到很好的反汇编输出：

movsd   (%rdi), %xmm0
movsd   8(%rdi), %xmm2
subsd   (%rsi), %xmm0
movsd   16(%rdi), %xmm1
subsd   8(%rsi), %xmm2
subsd   16(%rsi), %xmm1
mulsd   %xmm0, %xmm0
mulsd   %xmm2, %xmm2
mulsd   %xmm1, %xmm1
addsd   %xmm2, %xmm0
addsd   %xmm1, %xmm0

【讨论】：

虽然它使用 SSE2 指令，但这不是 SIMD 代码。所有指令都对单个值进行操作。 mulsd 表示“Scalar Double-Precision Floating-Point Multiply”，多数据版本是 mulpd“Packed Double-Precision Floating-Point Multiply”。
我已经检查过了，这也是我得到的。所以问题变成了：“我如何编写代码，以便 gcc 将其编译为 SIMD 代码，或者我是否需要手动编写，如果需要，如何编写？”。

【解决方案5】：

这类问题经常出现在 MD 模拟中。通常通过截断和邻居列表来减少计算量，因此减少了计算的数量。但是，平方距离的实际计算完全按照您的问题中给出的方式完成（使用编译器优化和固定类型 float[3]）。

因此，如果您想减少平方计算量，您应该告诉我们更多有关该问题的信息。

【讨论】：

【解决方案6】：

也许直接将 6 个双精度值作为参数传递可以使其更快（因为它可以避免数组取消引用）：

inline double distsquare_coord(double xa, double ya, double za, 
                               double xb, double yb, double zb) 
{ 
  double dx = xa-yb; double dy=ya-yb; double dz=za-zb;
  return dx*dx + dy*dy + dz*dz; 
}

或者，如果您在附近有很多点，您可以通过线性近似其他近点的距离来计算距离（到相同的固定其他点）。

【讨论】：

您仍然需要对数组进行取消引用，您只是将成本移到函数之外。
同意@Loki Astari，并且有可能 6 个参数无法容纳 CPU 寄存器......
在 AMD64 和 Linux 上，我认为 6 个参数适合寄存器（至少 6 个整数或指针参数适合，我忘了浮点数）。
很可能，调用代码可能已经将这 6 个值保存在寄存器中。

【解决方案7】：

如果您可以重新排列数据以一次处理两对输入向量，则可以使用此代码（仅限 SSE2）

// @brief Computes two squared distances between two pairs of 3D vectors
// @param a
//  Pointer to the first pair of 3D vectors.
//  The two vectors must be stored with stride 24, i.e. (a + 3) should point to the first component of the second vector in the pair.
//  Must be aligned by 16 (2 doubles).
// @param b
//  Pointer to the second pairs of 3D vectors.
//  The two vectors must be stored with stride 24, i.e. (a + 3) should point to the first component of the second vector in the pair.
//  Must be aligned by 16 (2 doubles).
// @param c
//  Pointer to the output 2 element array.
//  Must be aligned by 16 (2 doubles).
//  The two distances between a and b vectors will be written to c[0] and c[1] respectively.
void (const double * __restrict__ a, const double * __restrict__ b, double * __restrict c) {
    // diff0 = ( a0.y - b0.y, a0.x - b0.x ) = ( d0.y, d0.x )
    __m128d diff0 = _mm_sub_pd(_mm_load_pd(a), _mm_load_pd(b));
    // diff1 = ( a1.x - b1.x, a0.z - b0.z ) = ( d1.x, d0.z )
    __m128d diff1 = _mm_sub_pd(_mm_load_pd(a + 2), _mm_load_pd(b + 2));
    // diff2 = ( a1.z - b1.z, a1.y - b1.y ) = ( d1.z, d1.y )
    __m128d diff2 = _mm_sub_pd(_mm_load_pd(a + 4), _mm_load_pd(b + 4));

    // prod0 = ( d0.y * d0.y, d0.x * d0.x )
    __m128d prod0 = _mm_mul_pd(diff0, diff0);
    // prod1 = ( d1.x * d1.x, d0.z * d0.z )
    __m128d prod1 = _mm_mul_pd(diff1, diff1);
    // prod2 = ( d1.z * d1.z, d1.y * d1.y )
    __m128d prod2 = _mm_mul_pd(diff1, diff1);

    // _mm_unpacklo_pd(prod0, prod2) = ( d1.y * d1.y, d0.x * d0.x )
    // psum = ( d1.x * d1.x + d1.y * d1.y, d0.x * d0.x + d0.z * d0.z )
    __m128d psum = _mm_add_pd(_mm_unpacklo_pd(prod0, prod2), prod1);

    // _mm_unpackhi_pd(prod0, prod2) = ( d1.z * d1.z, d0.y * d0.y )
    // dotprod = ( d1.x * d1.x + d1.y * d1.y + d1.z * d1.z, d0.x * d0.x + d0.y * d0.y + d0.z * d0.z )
    __m128d dotprod = _mm_add_pd(_mm_unpackhi_pd(prod0, prod2), psum);

    __mm_store_pd(c, dotprod);
}

【讨论】：