title: Euclidean distance using intrinsic instruction
body: <p>For a research project, I need to compute a lot of Euclidean distances, where certain dimensions must be selected and others discarded. In the current state of the program, the array of selected dimensions has 100-ish elements, and I compute around 2-3 million distances. My current piece of code is as follow :</p>
  
  <pre><code>float compute_distance(const float* p1, const float* p2) const
  {
      __m256 euclidean = _mm256_setzero_ps();
  
      const uint16_t n = nbr_dimensions;
      const uint16_t aligend_n = n - n % 16;
      const float* local_selected = selected_dimensions;
  
      for (uint16_t i = 0; i &lt; aligend_n; i += 16)
      {
          const __m256 r1 = _mm256_sub_ps(_mm256_load_ps(&amp;p1[i]), _mm256_load_ps(&amp;p2[i]));
          euclidean = _mm256_fmadd_ps(_mm256_mul_ps(r1, r1), _mm256_load_ps(&amp;local_selected[i]), euclidean);
          const __m256 r2 = _mm256_sub_ps(_mm256_load_ps(&amp;p1[i + 8]), _mm256_load_ps(&amp;p2[i + 8]));
          euclidean = _mm256_fmadd_ps(_mm256_mul_ps(r2, r2), _mm256_load_ps(&amp;local_selected[i + 8]), euclidean);
      }
      float distance = hsum256_ps_avx(euclidean);
  
      for (uint16_t i = aligend_n; i &lt; n; ++i)
      {
          const float num = p1[i] - p2[i];
          distance += num * num * local_selected[i];
      }
  
      return distance;
  }
  </code></pre>
  
  <p>The selected dimensions are pre-determined. I could thus pre-compute an array of  <code>__m256</code> to pass to a <code>_mm256_blendv_ps</code> instead of multiplying by 0 or 1 in the line <code>euclidean = _mm256_fmadd_ps(_mm256_mul_ps(r1, r1), _mm256_load_ps(&amp;local_selected[i]), euclidean);</code>. But I'm rather a novice in intrinsic instructions, and I have not yet found a working solution.</p>
  
  <p>I was wondering if you guys could have some advice, or even code suggestions,  to improve the running speed of this function. As a side note, I do not have access to AVX-512 instructions.</p>
  
  <hr>
  
  <p>Update:
  Using the first above-mentioned solution, it comes:</p>
  
  <pre><code>float compute_distance(const float* p1, const float* p2) const
  {
      const size_t n = nbr_dimensions;
      const size_t aligend_n = n - n % 16;
      const unsigned int* local_selected = selected_dimensions;
      const __m256* local_masks = masks;
  
      __m256 euc1 = _mm256_setzero_ps(), euc2 = _mm256_setzero_ps(),
          euc3 = _mm256_setzero_ps(), euc4 = _mm256_setzero_ps();
  
      const size_t n_max = aligend_n/8;
      for (size_t i = 0; i &lt; n_max; i += 4)
      {       
          const __m256 r1 = _mm256_sub_ps(_mm256_load_ps(&amp;p1[i * 8 + 0]), _mm256_load_ps(&amp;p2[i * 8 + 0]));
          const __m256 r1_1 = _mm256_and_ps(r1, local_masks[i + 0]);
          euc1 = _mm256_fmadd_ps(r1_1, r1_1, euc1);
  
          const __m256 r2 = _mm256_sub_ps(_mm256_load_ps(&amp;p1[i * 8 + 8]), _mm256_load_ps(&amp;p2[i * 8 + 8]));
          const __m256 r2_1 = _mm256_and_ps(r2, local_masks[i + 1]);
          euc2 = _mm256_fmadd_ps(r2_1, r2_1, euc2);
  
          const __m256 r3 = _mm256_sub_ps(_mm256_load_ps(&amp;p1[i * 8 + 16]), _mm256_load_ps(&amp;p2[i * 8 + 16]));
          const __m256 r3_1 = _mm256_and_ps(r3, local_masks[i + 2]);
          euc3 = _mm256_fmadd_ps(r3_1, r3_1, euc3);
  
          const __m256 r4 = _mm256_sub_ps(_mm256_load_ps(&amp;p1[i * 8 + 24]), _mm256_load_ps(&amp;p2[i * 8 + 24]));
          const __m256 r4_1 = _mm256_and_ps(r4, local_masks[i + 3]);
          euc4 = _mm256_fmadd_ps(r4_1, r4_1, euc4);
      }
  
      float distance = hsum256_ps_avx(_mm256_add_ps(_mm256_add_ps(euc1, euc2), _mm256_add_ps(euc3, euc4)));
  
      for (size_t i = aligend_n; i &lt; n; ++i)
      {
          const float num = p1[i] - p2[i];
          distance += num * num * local_selected[i];      
      }
  
      return distance;
  }
  </code></pre>