text: /// \copyright Unlicense
  
  #pragma once
  #ifndef RRASSEMBLY_EUCLIDEAN_DISTANCE_INL_
  #define RRASSEMBLY_EUCLIDEAN_DISTANCE_INL_
  
  #include "euclidean_distance.hxx"
  
  #include <cassert>
  
  #include <immintrin.h>
  
  RRASSEMBLY_BEGIN_NAMESPACE
  
  namespace details {
  
  inline float squared_euclidean_distance(const float *pt1, const float *pt2, std::ptrdiff_t pt_dim)
  {
      assert(pt1);
      assert(pt2);
      assert(0 < pt_dim);
  
      float sum_of_sqr {};  // 平方和（sum of square）
      for (std::ptrdiff_t i = 0; i < pt_dim; ++i)
          sum_of_sqr += (pt1[i] - pt2[i]) * (pt1[i] - pt2[i]);
      return sum_of_sqr;
  }
  
  #ifdef __AVX__
  inline float squared_euclidean_distance_128(const float *pt1, const float *pt2, std::ptrdiff_t pt_dim)
  {
      constexpr std::ptrdiff_t step_len {sizeof(__m128) / sizeof(float)};
      assert(pt1);
      assert(pt2);
      assert(0 < pt_dim);
      assert(pt_dim % step_len == 0);
  
      float sum_of_sqr {};
      for (std::ptrdiff_t offset = 0; offset < pt_dim; offset += step_len) {
          const __m128 xmm0 = _mm_loadu_ps(pt1 + offset);
          const __m128 xmm1 = _mm_loadu_ps(pt2 + offset);
          const __m128 xmm2 = _mm_sub_ps(xmm0, xmm1);
          const __m128 xmm3 = _mm_dp_ps(xmm2, xmm2, 0b1111'0001);
          sum_of_sqr += xmm3.m128_f32[0];
      }
      return sum_of_sqr;
  }
  #endif
  
  #ifdef __AVX2__
  inline float squared_euclidean_distance_256(const float *pt1, const float *pt2, std::ptrdiff_t pt_dim)
  {
      constexpr std::ptrdiff_t step_len {sizeof(__m256) / sizeof(float)};
      assert(pt1);
      assert(pt2);
      assert(0 < pt_dim);
      assert(pt_dim % step_len == 0);
  
      float sum_of_sqr {};
      for (std::ptrdiff_t offset = 0; offset < pt_dim; offset += step_len) {
          const __m256 ymm0 = _mm256_loadu_ps(pt1 + offset);
          const __m256 ymm1 = _mm256_loadu_ps(pt2 + offset);
          const __m256 ymm2 = _mm256_sub_ps(ymm0, ymm1);
          const __m256 ymm3 = _mm256_dp_ps(ymm2, ymm2, 0b1111'0001);
          sum_of_sqr += ymm3.m256_f32[0] + ymm3.m256_f32[4];
      }
      return sum_of_sqr;
  }
  #endif
  
  }
  
  RRASSEMBLY_API float squared_euclidean_distance(const float *pt1, const float *pt2, std::ptrdiff_t pt_dim)
  {
      assert(pt1);
      assert(pt2);
      assert(0 < pt_dim);
  
      auto &it1 = pt1;
      auto &it2 = pt2;
      auto &left = pt_dim;
  
      float sum_of_sqr {};
  #ifdef __AVX2__
      do {
          constexpr std::ptrdiff_t step_len {sizeof(__m256) / sizeof(float)};
          const std::ptrdiff_t delta {left - left % step_len};
          if (delta == 0)
              break;
          sum_of_sqr += details::squared_euclidean_distance_256(it1, it2, delta);
          it1 += delta;
          it2 += delta;
          left -= delta;
          if (left == 0)
              return sum_of_sqr;
      } while (false);
  #endif
  #ifdef __AVX__
      do {
          constexpr std::ptrdiff_t step_len {sizeof(__m128) / sizeof(float)};
          const std::ptrdiff_t delta {left - left % step_len};
          if (delta == 0)
              break;
          sum_of_sqr += details::squared_euclidean_distance_128(it1, it2, delta);
          it1 += delta;
          it2 += delta;
          left -= delta;
          if (left == 0)
              return sum_of_sqr;
      } while (false);
  #endif
      sum_of_sqr += details::squared_euclidean_distance(pt1, pt2, pt_dim);
      return sum_of_sqr;
  }
  
  RRASSEMBLY_END_NAMESPACE
  
  #endif
  