Самый быстрый способ вычисления расстояния в квадрате

Question

Мой код сильно зависит от вычисления расстояния между двумя точками в трехмерном пространстве. Чтобы избежать дорогого квадратного корня, я использую квадратное расстояние. Но все же это занимает большую часть вычислительного времени, и я хотел бы заменить свою простую функцию чем-то еще более быстрым. Теперь у меня есть:

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, но не смог найти хороший пример или полный список инструкций (в идеале некоторые умножить + добавить на два вектора).

Графические процессоры не являются опцией, поскольку при каждом вызове функции известен только один набор точек.

Какой самый быстрый способ вычислить квадрат расстояния?

Answer 1

Хороший компилятор оптимизирует это так же хорошо, как и когда-либо. Хороший компилятор будет использовать инструкции SIMD, если сочтет их полезными. Убедитесь, что вы включили все возможные оптимизации для вашего компилятора. К сожалению, векторы измерения 3 не склонны хорошо сидеть с SIMD-единицами.

Я подозреваю, что вам просто придется признать, что код, сгенерированный компилятором, вероятно, очень близок к оптимальному, и что никаких существенных преимуществ достичь невозможно.

Answer 2

Первым очевидным было бы использование ключевого слова restrict.

Как и сейчас, a и b являются псевдонимами (и, таким образом, с точки зрения компилятора, который предполагает наихудший возможный случай, имеют псевдонимы). Ни один компилятор не будет автоматически векторизовать это, так как это неправильно.

Хуже того, компилятор не только не может векторизовать такой цикл, если вы также сохраняете (к счастью, не в вашем примере), он также должен каждый раз перезагружать значения. Всегда будьте осторожны с псевдонимами, так как это сильно влияет на компилятор.

Далее, если вы можете жить с этим, используйте float вместо double и pad до 4 с плавающей точкой, даже если один не используется, это более «естественная» компоновка данных для большинства процессоров (это несколько для конкретной платформы, но 4 числа с плавающей запятой - хорошее предположение для большинства платформ - 3 двойные, то есть 1,5 SIMD-регистра на «типичных» процессорах, нигде не оптимальны).

(Для рукописной реализации SIMD (которая сложнее, чем вы думаете), сначала и прежде всего убедитесь, что выровнены данные. Затем посмотрите, какие задержки у ваших инструкторов на целевой машине, и сначала выполните самые длинные. Например, в Intel pre-Prescott имеет смысл сначала перетасовать каждый компонент в регистр, а затем умножить на себя, даже если для этого используется 3 умножения вместо одного, поскольку перемешивание имеет большую задержку. На более поздних моделях перемешивание занимает один цикл, так что это будет полная антиоптимизация.
Что еще раз показывает, что оставить его компилятору - не такая уж плохая идея.)

Answer 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).Если у вас есть только три значения, вам нужно будет немного поиграться, чтобы получить требуемый формат, который бы отменил все вышеперечисленные преимущества.

В общем, из-за конвейерного суперскалярного преобразованияВ архитектуре ЦП наилучшим способом получения производительности является выполнение одной и той же операции с большим количеством данных. Таким образом, вы можете чередовать различные этапы и выполнять небольшую циклическую развертку, чтобы избежать остановок конвейера.Приведенный выше код определенно остановится на последних трех инструкциях, основанных на принципе «не может использовать значение сразу после его изменения» - вторая инструкция должна дождаться завершения результата предыдущей инструкции, что плохо втрубопроводная система.

Выполнение расчета по двум или более различным наборам точек одновременно может устранить вышеуказанное узкое место - во время ожидания результата одного вычисления вы можете начать вычисление следующей точки:

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

Answer 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

Answer 5

Этот тип проблемы часто возникает при моделировании MD.Обычно количество вычислений уменьшается за счет списков и списков соседей, поэтому число вычислений уменьшается.Однако фактический расчет квадратов расстояний выполняется точно (с оптимизацией компилятора и с плавающей точкой фиксированного типа [3]), как указано в вашем вопросе.

Так что, если вы хотите уменьшить количество квадратов вычислений, вы должны сказатьнам больше о проблеме.

Answer 6

Если вы можете переставить ваши данные для обработки двух пар входных векторов одновременно, вы можете использовать этот код (только 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);
}

Answer 7

Возможно, передача 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; 
}

Или, возможно, если у вас много точек поблизости, вы можете вычислить расстояние (до той же фиксированной другой точки) путем линейного приближения расстояний других близких точек.