计算距离平方的最快方法
我的代码在很大程度上依赖于计算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不是一个选项，因为每个函数调用只知道一组点。
计算距离平方的最快方法是什么？
一个优秀的编译器将优化您所能管理的内容。 如果一个好的编译器认为它们将是有益的，它将使用SIMD指令。 确保为编译器打开所有这些可能的优化。 不幸的是，尺寸为3的矢量不适合与SIMD单元配合使用。
我怀疑你只需要接受编译器生成的代码可能非常接近最优并且不能获得显着的增益。
第一个显而易见的事情是使用restrict关键字。
就像现在一样， a和b是可混淆的（因此，从编译器的角度来看，假设最坏的情况是别名）。 没有编译器会自动向量化这个，因为这样做是错误的。
更糟糕的是，编译器不仅无法对这样的循环进行矢量化，以防您也存储（幸运的是不在您的示例中），它还必须每次都重新加载值。 始终要清楚别名，因为它会极大地影响编译器。
接下来，如果你可以使用它，使用float而不是double和pad到4个浮点数，即使一个未使用，这对于大多数CPU来说是一个更“自然”的数据布局（这有点特定于平台，但是4个浮点数是对于大多数平台来说是一个很好的猜测 – 在“典型”CPU上，3个双打，也就是1.5个SIMD寄存器，在任何地方都不是最佳的。
（对于手写的SIMD实现（比你想象的更难），首先要确保对齐数据。然后，查看你的指令在目标机器上的延迟，并先做最长的延迟。例如在预Prescott Intel上，首先将每个组件拖放到一个寄存器然后与它自身相乘是有意义的，即使它使用3个乘法而不是1个，因为shuffle具有长延迟。在后面的模型中，shuffle需要一个周期，这将是一个完全的反优化。
这又表明将它留给编译器并不是一个坏主意。）
执行此操作的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
如果您想优化某些内容，请首先查看配置文件代码并检查汇编程序输出。
用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
这种类型的问题经常发生在MD模拟中。 通常，通过截止和邻居列表减少计算量，因此减少了计算的数量。 然而，平方距离的实际计算是完全完成的（使用编译器优化和固定类型float [3]），如问题中给出的那样。
因此，如果您想减少平方计算量，您应该告诉我们更多有关该问题的信息。
也许直接将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; }
或者，如果你附近有很多点，你可以通过线性近似其他近点的距离来计算距离（到同一个固定的其他点）。
如果您可以重新排列数据以一次处理两对输入向量，则可以使用此代码（仅限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, ie (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, ie (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); }