vector_x86.h

/*
 *  Copyright (C) 2004-2021 Edward F. Valeev
 *
 *  This file is part of Libint.
 *
 *  Libint is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  Libint is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with Libint.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
 
#ifndef _libint2_src_lib_libint_vectorx86_h_
#define _libint2_src_lib_libint_vectorx86_h_
 
#include <cstring>
#include <cmath>
#include <iostream>
 
#include <libint2/util/cxxstd.h>
#include <libint2/util/type_traits.h>
 
#if defined(__SSE2__) || defined(__SSE__) || defined(__AVX__)
#  include <x86intrin.h>
#endif
 
#ifdef __SSE2__
 
namespace libint2 { namespace simd {
 
  struct VectorSSEDouble {
 
      typedef double T;
      __m128d d;
 
      VectorSSEDouble() {}
 
      VectorSSEDouble(T a) {
        d = _mm_set_pd(a, a);
      }
 
      VectorSSEDouble(T (&a)[2]) {
        d = _mm_loadu_pd(&a[0]);
      }
 
      VectorSSEDouble(T a0, T a1) {
        d = _mm_set_pd(a1, a0);
      }
 
      VectorSSEDouble(__m128d a) {
        d = a;
      }
 
      VectorSSEDouble& operator=(T a) {
        d = _mm_set_pd(a, a);
        return *this;
      }
 
      VectorSSEDouble& operator+=(VectorSSEDouble a) {
        d = _mm_add_pd(d, a.d);
        return *this;
      }
 
      VectorSSEDouble& operator-=(VectorSSEDouble a) {
        d = _mm_sub_pd(d, a.d);
        return *this;
      }
 
      VectorSSEDouble operator-() const {
        // TODO improve using bitflips
        const static __m128d minus_one = _mm_set_pd(-1.0, -1.0);;
        VectorSSEDouble result;
        result.d = _mm_mul_pd(this->d, minus_one);
        return result;
      }
 
#if LIBINT2_CPLUSPLUS_STD >= 2011
      explicit
#endif
      operator double() const {
        double d0[2];
        ::memcpy(&(d0[0]), &d, sizeof(__m128d));
        // this segfaults on OS X
        //_mm_storeu_pd(&(d0[0]), d);
//        // aligned alternative requires C++11's alignas, but efficiency should not matter here
//        alignas(__m128d) double d0[2];
//        _mm_store_pd(&(d0[0]), d);
        return d0[0];
      }
 
      operator __m128d() const {
        return d;
      }
 
      void load(T const* a) {
        d = _mm_loadu_pd(a);
      }
      void load_aligned(T const* a) {
        d = _mm_load_pd(a);
      }
      void convert(T* a) const {
        _mm_storeu_pd(&a[0], d);
      }
      void convert_aligned(T* a) const {
        _mm_store_pd(&a[0], d);
      }
 
  };
 
  inline VectorSSEDouble operator*(double a, VectorSSEDouble b) {
    VectorSSEDouble c;
    VectorSSEDouble _a(a);
    c.d = _mm_mul_pd(_a.d, b.d);
    return c;
  }
 
  inline VectorSSEDouble operator*(VectorSSEDouble a, double b) {
    VectorSSEDouble c;
    VectorSSEDouble _b(b);
    c.d = _mm_mul_pd(a.d, _b.d);
    return c;
  }
 
  inline VectorSSEDouble operator*(int a, VectorSSEDouble b) {
    if (a == 1)
      return b;
    else {
      VectorSSEDouble c;
      VectorSSEDouble _a((double)a);
      c.d = _mm_mul_pd(_a.d, b.d);
      return c;
    }
  }
 
  inline VectorSSEDouble operator*(VectorSSEDouble a, int b) {
    if (b == 1)
      return a;
    else {
      VectorSSEDouble c;
      VectorSSEDouble _b((double)b);
      c.d = _mm_mul_pd(a.d, _b.d);
      return c;
    }
  }
 
  inline VectorSSEDouble operator*(VectorSSEDouble a, VectorSSEDouble b) {
    VectorSSEDouble c;
    c.d = _mm_mul_pd(a.d, b.d);
    return c;
  }
 
  inline VectorSSEDouble operator+(VectorSSEDouble a, VectorSSEDouble b) {
    VectorSSEDouble c;
    c.d = _mm_add_pd(a.d, b.d);
    return c;
  }
 
  inline VectorSSEDouble operator-(VectorSSEDouble a, VectorSSEDouble b) {
    VectorSSEDouble c;
    c.d = _mm_sub_pd(a.d, b.d);
    return c;
  }
 
  inline VectorSSEDouble operator/(VectorSSEDouble a, VectorSSEDouble b) {
    VectorSSEDouble c;
    c.d = _mm_div_pd(a.d, b.d);
    return c;
  }
 
#if defined(__FMA__)
  inline VectorSSEDouble fma_plus(VectorSSEDouble a, VectorSSEDouble b, VectorSSEDouble c) {
    VectorSSEDouble d;
    d.d = _mm_fmadd_pd(a.d, b.d, c.d);
    return d;
  }
  inline VectorSSEDouble fma_minus(VectorSSEDouble a, VectorSSEDouble b, VectorSSEDouble c) {
    VectorSSEDouble d;
    d.d = _mm_fmsub_pd(a.d, b.d, c.d);
    return d;
  }
#elif defined(__FMA4__)
  inline VectorSSEDouble fma_plus(VectorSSEDouble a, VectorSSEDouble b, VectorSSEDouble c) {
    VectorSSEDouble d;
    d.d = _mm_macc_pd(a.d, b.d, c.d);
    return d;
  }
  inline VectorSSEDouble fma_minus(VectorSSEDouble a, VectorSSEDouble b, VectorSSEDouble c) {
    VectorSSEDouble d;
    d.d = _mm_msub_pd(a.d, b.d, c.d);
    return d;
  }
#endif
 
  inline double horizontal_add (VectorSSEDouble const & a) {
//      Agner Fog's version
#if defined(__SSE3__)
    __m128d t1 = _mm_hadd_pd(a,a);
    return _mm_cvtsd_f64(t1);
#else // SSE2 only
    __m128  t0 = _mm_castpd_ps(a);
    __m128d t1 = _mm_castps_pd(_mm_movehl_ps(t0,t0));
    __m128d t2 = _mm_add_sd(a,t1);
    return _mm_cvtsd_f64(t2);
#endif
  }
 
  inline VectorSSEDouble horizontal_add (VectorSSEDouble const & a, VectorSSEDouble const & b) {
#if defined(__SSE3__)
    return _mm_hadd_pd(a, b);
#else // will be very inefficient without SSE3
    return VectorSSEDouble(horizontal_add(a), horizontal_add(b));
#endif
  }
 
 
  inline VectorSSEDouble exp(VectorSSEDouble a) {
#if HAVE_INTEL_SVML
    VectorSSEDouble result;
    result.d = _mm_exp_pd(a.d);
#else
    double a_d[2]; a.convert(a_d);
    for(int i=0; i<2; ++i) a_d[i] = std::exp(a_d[i]);
    VectorSSEDouble result(a_d);
#endif
    return result;
  }
  inline VectorSSEDouble sqrt(VectorSSEDouble a) {
#if HAVE_INTEL_SVML
    VectorSSEDouble result;
    result.d = _mm_sqrt_pd(a.d);
#else
    double a_d[2]; a.convert(a_d);
    for(int i=0; i<2; ++i) a_d[i] = std::sqrt(a_d[i]);
    VectorSSEDouble result(a_d);
#endif
    return result;
  }
  inline VectorSSEDouble erf(VectorSSEDouble a) {
#if HAVE_INTEL_SVML
    VectorSSEDouble result;
    result.d = _mm_erf_pd(a.d);
#else
    double a_d[2]; a.convert(a_d);
    for(int i=0; i<2; ++i) a_d[i] = ::erf(a_d[i]);
    VectorSSEDouble result(a_d);
#endif
    return result;
  }
  inline VectorSSEDouble erfc(VectorSSEDouble a) {
#if HAVE_INTEL_SVML
    VectorSSEDouble result;
    result.d = _mm_erfc_pd(a.d);
#else
    double a_d[2]; a.convert(a_d);
    for(int i=0; i<2; ++i) a_d[i] = ::erfc(a_d[i]);
    VectorSSEDouble result(a_d);
#endif
    return result;
  }
 
};}; // namespace libint2::simd
 
inline std::ostream& operator<<(std::ostream& os, libint2::simd::VectorSSEDouble a) {
  double ad[2];
  a.convert(ad);
  os << "{" << ad[0] << "," << ad[1] << "}";
  return os;
}
 
namespace libint2 {
 
 
  template <>
  struct is_vector<simd::VectorSSEDouble> {
      static const bool value = true;
  };
 
  template <>
  struct vector_traits<simd::VectorSSEDouble> {
      typedef double scalar_type;
      static const size_t extent = 2;
  };
 
 
} // namespace libint2
 
#endif // SSE2-only
 
#ifdef __SSE__
 
namespace libint2 { namespace simd {
 
  struct VectorSSEFloat {
 
      typedef float T;
      __m128 d;
 
      VectorSSEFloat() {}
 
      VectorSSEFloat(T a) {
        d = _mm_set_ps(a, a, a, a);
      }
 
      VectorSSEFloat(T (&a)[4]) {
        d = _mm_loadu_ps(&a[0]);
      }
 
      VectorSSEFloat(T a0, T a1, T a2, T a3) {
        d = _mm_set_ps(a3, a2, a1, a0);
      }
 
      VectorSSEFloat& operator=(T a) {
        d = _mm_set_ps(a, a, a, a);
        return *this;
      }
 
      VectorSSEFloat& operator+=(VectorSSEFloat a) {
        d = _mm_add_ps(d, a.d);
        return *this;
      }
 
      VectorSSEFloat& operator-=(VectorSSEFloat a) {
        d = _mm_sub_ps(d, a.d);
        return *this;
      }
 
      VectorSSEFloat operator-() const {
        // TODO improve using bitflips
        const static __m128 minus_one = _mm_set_ps(-1.0, -1.0, -1.0, -1.0);;
        VectorSSEFloat result;
        result.d = _mm_mul_ps(this->d, minus_one);
        return result;
      }
 
#if LIBINT2_CPLUSPLUS_STD >= 2011
      explicit
#endif
      operator float() const {
        float d0[4];
        ::memcpy(&(d0[0]), &d, sizeof(__m128));
        // this segfaults on OS X
        //_mm_storeu_ps(&(d0[0]), d);
//        // aligned alternative requires C++11's alignas, but efficiency should not matter here
//        alignas(__m128) float d0[4];
//        _mm_store_ps(&(d0[0]), d);
        return d0[0];
      }
 
#if LIBINT2_CPLUSPLUS_STD >= 2011
      explicit
#endif
      operator double() const {
        const float result_flt = this->operator float();
        return result_flt;
      }
 
      operator __m128() const {
        return d;
      }
 
      void load(T const* a) {
        d = _mm_loadu_ps(a);
      }
      void load_aligned(T const* a) {
        d = _mm_load_ps(a);
      }
      void convert(T* a) const {
        _mm_storeu_ps(&a[0], d);
      }
      void convert_aligned(T* a) const {
        _mm_store_ps(&a[0], d);
      }
  };
 
  inline VectorSSEFloat operator*(float a, VectorSSEFloat b) {
    VectorSSEFloat c;
    VectorSSEFloat _a(a);
    c.d = _mm_mul_ps(_a.d, b.d);
    return c;
  }
 
  inline VectorSSEFloat operator*(VectorSSEFloat a, float b) {
    VectorSSEFloat c;
    VectorSSEFloat _b(b);
    c.d = _mm_mul_ps(a.d, _b.d);
    return c;
  }
 
  // narrows a!
  inline VectorSSEFloat operator*(double a, VectorSSEFloat b) {
    VectorSSEFloat c;
    VectorSSEFloat _a((float)a);
    c.d = _mm_mul_ps(_a.d, b.d);
    return c;
  }
 
  // narrows b!
  inline VectorSSEFloat operator*(VectorSSEFloat a, double b) {
    VectorSSEFloat c;
    VectorSSEFloat _b((float)b);
    c.d = _mm_mul_ps(a.d, _b.d);
    return c;
  }
 
 
  inline VectorSSEFloat operator*(int a, VectorSSEFloat b) {
    if (a == 1)
      return b;
    else {
      VectorSSEFloat c;
      VectorSSEFloat _a((float)a);
      c.d = _mm_mul_ps(_a.d, b.d);
      return c;
    }
  }
 
  inline VectorSSEFloat operator*(VectorSSEFloat a, int b) {
    if (b == 1)
      return a;
    else {
      VectorSSEFloat c;
      VectorSSEFloat _b((float)b);
      c.d = _mm_mul_ps(a.d, _b.d);
      return c;
    }
  }
 
  inline VectorSSEFloat operator*(VectorSSEFloat a, VectorSSEFloat b) {
    VectorSSEFloat c;
    c.d = _mm_mul_ps(a.d, b.d);
    return c;
  }
 
  inline VectorSSEFloat operator+(VectorSSEFloat a, VectorSSEFloat b) {
    VectorSSEFloat c;
    c.d = _mm_add_ps(a.d, b.d);
    return c;
  }
 
  inline VectorSSEFloat operator-(VectorSSEFloat a, VectorSSEFloat b) {
    VectorSSEFloat c;
    c.d = _mm_sub_ps(a.d, b.d);
    return c;
  }
 
  inline VectorSSEFloat operator/(VectorSSEFloat a, VectorSSEFloat b) {
    VectorSSEFloat c;
    c.d = _mm_div_ps(a.d, b.d);
    return c;
  }
 
#if defined(__FMA__)
  inline VectorSSEFloat fma_plus(VectorSSEFloat a, VectorSSEFloat b, VectorSSEFloat c) {
    VectorSSEFloat d;
    d.d = _mm_fmadd_ps(a.d, b.d, c.d);
    return d;
  }
  inline VectorSSEFloat fma_minus(VectorSSEFloat a, VectorSSEFloat b, VectorSSEFloat c) {
    VectorSSEFloat d;
    d.d = _mm_fmsub_ps(a.d, b.d, c.d);
    return d;
  }
#elif defined(__FMA4__)
  inline VectorSSEFloat fma_plus(VectorSSEFloat a, VectorSSEFloat b, VectorSSEFloat c) {
    VectorSSEFloat d;
    d.d = _mm_macc_ps(a.d, b.d, c.d);
    return d;
  }
  inline VectorSSEFloat fma_minus(VectorSSEFloat a, VectorSSEFloat b, VectorSSEFloat c) {
    VectorSSEFloat d;
    d.d = _mm_msub_ps(a.d, b.d, c.d);
    return d;
  }
#endif
 
 
  inline VectorSSEFloat exp(VectorSSEFloat a) {
#if HAVE_INTEL_SVML
    VectorSSEFloat result;
    result.d = _mm_exp_ps(a.d);
#else
    float a_d[4]; a.convert(a_d);
    for(int i=0; i<4; ++i) a_d[i] = std::exp(a_d[i]);
    VectorSSEFloat result(a_d);
#endif
    return result;
  }
  inline VectorSSEFloat sqrt(VectorSSEFloat a) {
#if HAVE_INTEL_SVML
    VectorSSEFloat result;
    result.d = _mm_sqrt_ps(a.d);
#else
    float a_d[4]; a.convert(a_d);
    for(int i=0; i<4; ++i) a_d[i] = std::sqrt(a_d[i]);
    VectorSSEFloat result(a_d);
#endif
    return result;
  }
  inline VectorSSEFloat erf(VectorSSEFloat a) {
#if HAVE_INTEL_SVML
    VectorSSEFloat result;
    result.d = _mm_erf_ps(a.d);
#else
    float a_d[4]; a.convert(a_d);
    for(int i=0; i<4; ++i) a_d[i] = ::erf(a_d[i]);
    VectorSSEFloat result(a_d);
#endif
    return result;
  }
  inline VectorSSEFloat erfc(VectorSSEFloat a) {
#if HAVE_INTEL_SVML
    VectorSSEFloat result;
    result.d = _mm_erfc_ps(a.d);
#else
    float a_d[4]; a.convert(a_d);
    for(int i=0; i<4; ++i) a_d[i] = ::erfc(a_d[i]);
    VectorSSEFloat result(a_d);
#endif
    return result;
  }
 
};}; // namespace libint2::simd
 
inline std::ostream& operator<<(std::ostream& os, libint2::simd::VectorSSEFloat a) {
  float ad[4];
  a.convert(ad);
  os << "{" << ad[0] << "," << ad[1] << "," << ad[2] << "," << ad[3] << "}";
  return os;
}
 
namespace libint2 {
 
 
  template <>
  struct is_vector<simd::VectorSSEFloat> {
      static const bool value = true;
  };
 
  template <>
  struct vector_traits<simd::VectorSSEFloat> {
      typedef float scalar_type;
      static const size_t extent = 4;
  };
 
 
} // namespace libint2
 
#endif // SSE-only
 
#ifdef __AVX__
 
namespace libint2 { namespace simd {
 
  struct VectorAVXDouble {
 
      typedef double T;
      __m256d d;
 
      VectorAVXDouble() {}
 
      VectorAVXDouble(T a) {
        d = _mm256_set_pd(a, a, a, a);
      }
 
      VectorAVXDouble(T (&a)[4]) {
        d = _mm256_loadu_pd(&a[0]);
      }
 
      VectorAVXDouble(T a0, T a1, T a2, T a3) {
        d = _mm256_set_pd(a3, a2, a1, a0);
      }
 
      VectorAVXDouble(__m256d a) {
        d = a;
      }
 
      VectorAVXDouble& operator=(T a) {
        d = _mm256_set_pd(a, a, a, a);
        return *this;
      }
 
      VectorAVXDouble& operator+=(VectorAVXDouble a) {
        d = _mm256_add_pd(d, a.d);
        return *this;
      }
 
      VectorAVXDouble& operator-=(VectorAVXDouble a) {
        d = _mm256_sub_pd(d, a.d);
        return *this;
      }
 
      VectorAVXDouble operator-() const {
        // TODO improve using bitflips
        const static __m256d minus_one = _mm256_set_pd(-1.0, -1.0, -1.0, -1.0);;
        VectorAVXDouble result;
        result.d = _mm256_mul_pd(this->d, minus_one);
        return result;
      }
 
#if LIBINT2_CPLUSPLUS_STD >= 2011
      explicit
#endif
      operator double() const {
        double d0[4];
        ::memcpy(&(d0[0]), &d, sizeof(__m256d));
        // this segfaults on OS X
//        _mm256_storeu_pd(&d0[0], d);
//        // aligned alternative requires C++11's alignas, but efficiency should not matter here
//        //        alignas(__m256d) double d0[4];
//        //        _mm256_store_pd(&(d0[0]), d);
        return d0[0];
      }
 
      operator __m256d() const {
        return d;
      }
 
      void load(T const* a) {
        d = _mm256_loadu_pd(a);
      }
      void load_aligned(T const* a) {
        d = _mm256_load_pd(a);
      }
      void convert(T* a) const {
        _mm256_storeu_pd(&a[0], d);
      }
      void convert_aligned(T* a) const {
        _mm256_store_pd(&a[0], d);
      }
  };
 
  inline VectorAVXDouble operator*(double a, VectorAVXDouble b) {
    VectorAVXDouble c;
    VectorAVXDouble _a(a);
    c.d = _mm256_mul_pd(_a.d, b.d);
    return c;
  }
 
  inline VectorAVXDouble operator*(VectorAVXDouble a, double b) {
    VectorAVXDouble c;
    VectorAVXDouble _b(b);
    c.d = _mm256_mul_pd(a.d, _b.d);
    return c;
  }
 
  inline VectorAVXDouble operator*(int a, VectorAVXDouble b) {
    if (a == 1)
      return b;
    else {
      VectorAVXDouble c;
      VectorAVXDouble _a((double)a);
      c.d = _mm256_mul_pd(_a.d, b.d);
      return c;
    }
  }
 
  inline VectorAVXDouble operator*(VectorAVXDouble a, int b) {
    if (b == 1)
      return a;
    else {
      VectorAVXDouble c;
      VectorAVXDouble _b((double)b);
      c.d = _mm256_mul_pd(a.d, _b.d);
      return c;
    }
  }
 
  inline VectorAVXDouble operator*(VectorAVXDouble a, VectorAVXDouble b) {
    VectorAVXDouble c;
    c.d = _mm256_mul_pd(a.d, b.d);
    return c;
  }
 
  inline VectorAVXDouble operator+(VectorAVXDouble a, VectorAVXDouble b) {
    VectorAVXDouble c;
    c.d = _mm256_add_pd(a.d, b.d);
    return c;
  }
 
  inline VectorAVXDouble operator+(int a, VectorAVXDouble b) {
    if(a == 0){
      return b;
    } 
    else {
      VectorAVXDouble c;
      VectorAVXDouble _a = (static_cast<double>(a));
      c.d = _mm256_add_pd(_a.d,b.d);
      return c;
    }
  }
 
 
  inline VectorAVXDouble operator+(VectorAVXDouble a, int b) {
    if(b == 0){
      return a;
    } 
    else {
      VectorAVXDouble c;
      VectorAVXDouble _b = (static_cast<double>(b));
      c.d = _mm256_add_pd(a.d,_b.d);
      return c;
    }
  }
 
 
  inline VectorAVXDouble operator-(VectorAVXDouble a, VectorAVXDouble b) {
    VectorAVXDouble c;
    c.d = _mm256_sub_pd(a.d, b.d);
    return c;
  }
 
  inline VectorAVXDouble operator/(VectorAVXDouble a, VectorAVXDouble b) {
    VectorAVXDouble c;
    c.d = _mm256_div_pd(a.d, b.d);
    return c;
  }
 
#if defined(__FMA__)
  inline VectorAVXDouble fma_plus(VectorAVXDouble a, VectorAVXDouble b, VectorAVXDouble c) {
    VectorAVXDouble d;
    d.d = _mm256_fmadd_pd(a.d, b.d, c.d);
    return d;
  }
  inline VectorAVXDouble fma_minus(VectorAVXDouble a, VectorAVXDouble b, VectorAVXDouble c) {
    VectorAVXDouble d;
    d.d = _mm256_fmsub_pd(a.d, b.d, c.d);
    return d;
  }
#elif defined(__FMA4__)
  inline VectorAVXDouble fma_plus(VectorAVXDouble a, VectorAVXDouble b, VectorAVXDouble c) {
    VectorAVXDouble d;
    d.d = _mm256_macc_pd(a.d, b.d, c.d);
    return d;
  }
  inline VectorAVXDouble fma_minus(VectorAVXDouble a, VectorAVXDouble b, VectorAVXDouble c) {
    VectorAVXDouble d;
    d.d = _mm256_msub_pd(a.d, b.d, c.d);
    return d;
  }
#endif
 
  inline double horizontal_add (VectorAVXDouble const & a) {
      __m256d s = _mm256_hadd_pd(a,a);
      return ((double*)&s)[0] + ((double*)&s)[2];
//      Agner Fog's version
//      __m256d t1 = _mm256_hadd_pd(a,a);
//      __m128d t2 = _mm256_extractf128_pd(t1,1);
//      __m128d t3 = _mm_add_sd(_mm256_castpd256_pd128(t1),t2);
//      return _mm_cvtsd_f64(t3);
  }
 
  inline VectorSSEDouble horizontal_add (VectorAVXDouble const & a, VectorAVXDouble const & b) {
    // solution from http://stackoverflow.com/questions/9775538/fastest-way-to-do-horizontal-vector-sum-with-avx-instructions
    __m256d sum = _mm256_hadd_pd(a, b);
    // extract upper 128 bits of result
    __m128d sum_high = _mm256_extractf128_pd(sum, 1);
    // add upper 128 bits of sum to its lower 128 bits
    return _mm_add_pd(sum_high, _mm256_castpd256_pd128(sum));
  }
 
  inline VectorAVXDouble horizontal_add (VectorAVXDouble const & a,
                                         VectorAVXDouble const & b,
                                         VectorAVXDouble const & c,
                                         VectorAVXDouble const & d) {
    // solution from http://stackoverflow.com/questions/10833234/4-horizontal-double-precision-sums-in-one-go-with-avx?lq=1
    // {a[0]+a[1], b[0]+b[1], a[2]+a[3], b[2]+b[3]}
    __m256d sumab = _mm256_hadd_pd(a, b);
    // {c[0]+c[1], d[0]+d[1], c[2]+c[3], d[2]+d[3]}
    __m256d sumcd = _mm256_hadd_pd(c, d);
 
    // {a[0]+a[1], b[0]+b[1], c[2]+c[3], d[2]+d[3]}
    __m256d blend = _mm256_blend_pd(sumab, sumcd, 0b1100);
    // {a[2]+a[3], b[2]+b[3], c[0]+c[1], d[0]+d[1]}
    __m256d perm = _mm256_permute2f128_pd(sumab, sumcd, 0x21);
 
    return _mm256_add_pd(perm, blend);
  }
 
 
  inline VectorAVXDouble exp(VectorAVXDouble a) {
#if HAVE_INTEL_SVML
    VectorAVXDouble result;
    result.d = _mm256_exp_pd(a.d);
#else
    double a_d[4]; a.convert(a_d);
    for(int i=0; i<4; ++i) a_d[i] = ::exp(a_d[i]);
    VectorAVXDouble result(a_d);
#endif
    return result;
  }
  inline VectorAVXDouble sqrt(VectorAVXDouble a) {
#if HAVE_INTEL_SVML
    VectorAVXDouble result;
    result.d = _mm256_sqrt_pd(a.d);
#else
    double a_d[4]; a.convert(a_d);
    for(int i=0; i<4; ++i) a_d[i] = ::sqrt(a_d[i]);
    VectorAVXDouble result(a_d);
#endif
    return result;
  }
  inline VectorAVXDouble erf(VectorAVXDouble a) {
#if HAVE_INTEL_SVML
    VectorAVXDouble result;
    result.d = _mm256_erf_pd(a.d);
#else
    double a_d[4]; a.convert(a_d);
    for(int i=0; i<4; ++i) a_d[i] = ::erf(a_d[i]);
    VectorAVXDouble result(a_d);
#endif
    return result;
  }
  inline VectorAVXDouble erfc(VectorAVXDouble a) {
#if HAVE_INTEL_SVML
    VectorAVXDouble result;
    result.d = _mm256_erfc_pd(a.d);
#else
    double a_d[4]; a.convert(a_d);
    for(int i=0; i<4; ++i) a_d[i] = ::erfc(a_d[i]);
    VectorAVXDouble result(a_d);
#endif
    return result;
  }
 
 
 
 
  struct VectorAVXFloat {
 
      typedef float T;
      __m256 d;
 
      VectorAVXFloat() {}
 
      VectorAVXFloat(T a) {
        d = _mm256_set_ps(a, a, a, a, a, a, a, a);
      }
 
      VectorAVXFloat(T (&a)[8]) {
        d = _mm256_loadu_ps(&a[0]);
      }
 
      VectorAVXFloat(T a0, T a1, T a2, T a3, T a4, T a5, T a6, T a7) {
        d = _mm256_set_ps(a0, a1, a2, a3, a4, a5, a6, a7);
      }
 
      VectorAVXFloat& operator=(T a) {
        d = _mm256_set_ps(a, a, a, a, a, a, a, a);
        return *this;
      }
 
      VectorAVXFloat& operator+=(VectorAVXFloat a) {
        d = _mm256_add_ps(d, a.d);
        return *this;
      }
 
      VectorAVXFloat& operator-=(VectorAVXFloat a) {
        d = _mm256_sub_ps(d, a.d);
        return *this;
      }
 
      VectorAVXFloat operator-() const {
        // TODO improve using bitflips
        const static __m256 minus_one = _mm256_set_ps(-1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0);;
        VectorAVXFloat result;
        result.d = _mm256_mul_ps(this->d, minus_one);
        return result;
      }
 
      explicit operator float() const {
        float d0[8];
        ::memcpy(&(d0[0]), &d, sizeof(__m256));
        // this segfaults on OS X
//        _mm256_storeu_ps(&d0[0], d);
//        // aligned alternative requires C++11's alignas, but efficiency should not matter here
//        //        alignas(__m256) float d0[8];
//        //        _mm256_store_ps(&(d0[0]), d);
        return d0[0];
      }
 
      explicit operator double() const {
              const float result_flt = this->operator float();
              return result_flt;
      }
 
      void convert(T(&a)[8]) const {
        _mm256_storeu_ps(&a[0], d);
      }
  };
 
  inline VectorAVXFloat operator*(double a, VectorAVXFloat b) {
    VectorAVXFloat c;
    VectorAVXFloat _a(a);
    c.d = _mm256_mul_ps(_a.d, b.d);
    return c;
  }
 
  inline VectorAVXFloat operator*(VectorAVXFloat a, double b) {
    VectorAVXFloat c;
    VectorAVXFloat _b(b);
    c.d = _mm256_mul_ps(a.d, _b.d);
    return c;
  }
 
  inline VectorAVXFloat operator*(int a, VectorAVXFloat b) {
    if (a == 1)
      return b;
    else {
      VectorAVXFloat c;
      VectorAVXFloat _a((float)a);
      c.d = _mm256_mul_ps(_a.d, b.d);
      return c;
    }
  }
 
  inline VectorAVXFloat operator*(VectorAVXFloat a, int b) {
    if (b == 1)
      return a;
    else {
      VectorAVXFloat c;
      VectorAVXFloat _b((float)b);
      c.d = _mm256_mul_ps(a.d, _b.d);
      return c;
    }
  }
 
  inline VectorAVXFloat operator*(VectorAVXFloat a, VectorAVXFloat b) {
    VectorAVXFloat c;
    c.d = _mm256_mul_ps(a.d, b.d);
    return c;
  }
 
  inline VectorAVXFloat operator+(VectorAVXFloat a, VectorAVXFloat b) {
    VectorAVXFloat c;
    c.d = _mm256_add_ps(a.d, b.d);
    return c;
  }
 
  inline VectorAVXFloat operator-(VectorAVXFloat a, VectorAVXFloat b) {
    VectorAVXFloat c;
    c.d = _mm256_sub_ps(a.d, b.d);
    return c;
  }
 
  inline VectorAVXFloat operator/(VectorAVXFloat a, VectorAVXFloat b) {
    VectorAVXFloat c;
    c.d = _mm256_div_ps(a.d, b.d);
    return c;
  }
 
#if defined(__FMA__)
  inline VectorAVXFloat fma_plus(VectorAVXFloat a, VectorAVXFloat b, VectorAVXFloat c) {
    VectorAVXFloat d;
    d.d = _mm256_fmadd_ps(a.d, b.d, c.d);
    return d;
  }
  inline VectorAVXFloat fma_minus(VectorAVXFloat a, VectorAVXFloat b, VectorAVXFloat c) {
    VectorAVXFloat d;
    d.d = _mm256_fmsub_ps(a.d, b.d, c.d);
    return d;
  }
#elif defined(__FMA4__)
  inline VectorAVXFloat fma_plus(VectorAVXFloat a, VectorAVXFloat b, VectorAVXFloat c) {
    VectorAVXFloat d;
    d.d = _mm256_facc_ps(a.d, b.d, c.d);
    return d;
  }
  inline VectorAVXFloat fma_minus(VectorAVXFloat a, VectorAVXFloat b, VectorAVXFloat c) {
    VectorAVXFloat d;
    d.d = _mm256_fsub_ps(a.d, b.d, c.d);
    return d;
  }
#endif
 
 
  inline VectorAVXFloat exp(VectorAVXFloat a) {
#if HAVE_INTEL_SVML
    VectorAVXFloat result;
    result.d = _mm256_exp_ps(a.d);
#else
    float a_d[8]; a.convert(a_d);
    for(int i=0; i<8; ++i) a_d[i] = ::exp(a_d[i]);
    VectorAVXFloat result(a_d);
#endif
    return result;
  }
  inline VectorAVXFloat sqrt(VectorAVXFloat a) {
#if HAVE_INTEL_SVML
    VectorAVXFloat result;
    result.d = _mm256_sqrt_ps(a.d);
#else
    float a_d[8]; a.convert(a_d);
    for(int i=0; i<8; ++i) a_d[i] = ::sqrt(a_d[i]);
    VectorAVXFloat result(a_d);
#endif
    return result;
  }
  inline VectorAVXFloat erf(VectorAVXFloat a) {
#if HAVE_INTEL_SVML
    VectorAVXFloat result;
    result.d = _mm256_erf_ps(a.d);
#else
    float a_d[8]; a.convert(a_d);
    for(int i=0; i<8; ++i) a_d[i] = ::erf(a_d[i]);
    VectorAVXFloat result(a_d);
#endif
    return result;
  }
  inline VectorAVXFloat erfc(VectorAVXFloat a) {
#if HAVE_INTEL_SVML
    VectorAVXFloat result;
    result.d = _mm256_erfc_ps(a.d);
#else
    float a_d[8]; a.convert(a_d);
    for(int i=0; i<8; ++i) a_d[i] = ::erfc(a_d[i]);
    VectorAVXFloat result(a_d);
#endif
    return result;
  }
 
};}; // namespace libint2::simd
 
inline std::ostream& operator<<(std::ostream& os, libint2::simd::VectorAVXDouble a) {
  double ad[4];
  a.convert(ad);
  os << "{" << ad[0] << "," << ad[1] << "," << ad[2] << "," << ad[3] << "}";
  return os;
}
 
inline std::ostream& operator<<(std::ostream& os, libint2::simd::VectorAVXFloat a) {
  float ad[8];
  a.convert(ad);
  os << "{" << ad[0] << "," << ad[1] << "," << ad[2] << "," << ad[3] << "," << ad[4]<< "," << ad[5]<< ","  << ad[6]<< "," << ad[7] << "}";
  return os;
}
 
namespace libint2 {
 
 
  template <>
  struct is_vector<simd::VectorAVXDouble> {
      static const bool value = true;
  };
 
  template <>
  struct vector_traits<simd::VectorAVXDouble> {
      typedef double scalar_type;
      static const size_t extent = 4;
  };
 
 
} // namespace libint2
 
#endif // AVX-only
 
 
#ifdef LIBINT2_HAVE_AGNER_VECTORCLASS
#include <vectorclass.h>
#endif
 
#endif // header guard