merge Libin's upgrade for AVX512

include/sctl/vec.hpp

@@ -25,28 +25,79 @@
 #include <immintrin.h>
 #endif
 
+// TODO: Implement AVX versions of floats, int32_t, int64_t
+
 // TODO: Check alignment when SCTL_MEMDEBUG is defined
 // TODO: Replace pointers with iterators
 
 namespace SCTL_NAMESPACE {
 
+  enum class DataType {
+    Integer,
+    Real,
+    Bool
+  };
+
   template <class ValueType> class TypeTraits {
     public:
-      static constexpr Integer SigBits = 0;
+      static constexpr DataType Type = DataType::Bool;
+      static constexpr Integer Size = sizeof(ValueType);
+      static constexpr Integer SigBits = 1;
+  };
+  template <> class TypeTraits<int32_t> {
+    public:
+      static constexpr DataType Type = DataType::Integer;
+      static constexpr Integer Size = sizeof(int32_t);
+      static constexpr Integer SigBits = Size * 8;
+  };
+  template <> class TypeTraits<int64_t> {
+    public:
+      static constexpr DataType Type = DataType::Integer;
+      static constexpr Integer Size = sizeof(int64_t);
+      static constexpr Integer SigBits = Size * 8;
   };
   template <> class TypeTraits<float> {
     public:
+      static constexpr DataType Type = DataType::Real;
+      static constexpr Integer Size = sizeof(float);
       static constexpr Integer SigBits = 23;
   };
   template <> class TypeTraits<double> {
     public:
+      static constexpr DataType Type = DataType::Real;
+      static constexpr Integer Size = sizeof(double);
       static constexpr Integer SigBits = 52;
   };
 
+  template <DataType type, Integer size> class GetType {
+    public:
+      typedef bool ValueType;
+  };
+  template <> class GetType<DataType::Integer,4> {
+    public:
+      typedef int32_t ValueType;
+  };
+  template <> class GetType<DataType::Integer,8> {
+    public:
+      typedef int64_t ValueType;
+  };
+  template <> class GetType<DataType::Real,4> {
+    public:
+      typedef float ValueType;
+  };
+  template <> class GetType<DataType::Real,8> {
+    public:
+      typedef double ValueType;
+  };
 
   template <class ValueType, Integer N> class alignas(sizeof(ValueType) * N) Vec {
+
     public:
 
+      typedef typename GetType<DataType::Integer,TypeTraits<ValueType>::Size>::ValueType IntegerType;
+      typedef typename GetType<DataType::Real,TypeTraits<ValueType>::Size>::ValueType RealType;
+      typedef Vec<IntegerType,N> IntegerVec;
+      typedef Vec<RealType,N> RealVec;
       typedef ValueType ScalarType;
 
       static constexpr Integer Size() {
@@ -75,7 +126,7 @@ namespace SCTL_NAMESPACE {
         return r;
       }
 
-      Vec() {}
+      Vec() = default;
 
       Vec(const ValueType& a) {
         for (Integer i = 0; i < N; i++) v[i] = a;
@@ -108,11 +159,11 @@ namespace SCTL_NAMESPACE {
       }
 
       // C-style cast
-      template <class RetValueType> explicit operator Vec<RetValueType,N>() const {
-        Vec<RetValueType,N> r;
-        for (Integer i = 0; i < N; i++) r.v[i] = (RetValueType)v[i];
-        return r;
-      }
+      //template <class RetValueType> explicit operator Vec<RetValueType,N>() const {
+      //  Vec<RetValueType,N> r;
+      //  for (Integer i = 0; i < N; i++) r.v[i] = (RetValueType)v[i];
+      //  return r;
+      //}
 
       // Arithmetic operators
       friend Vec operator*(Vec lhs, const Vec& rhs) {
@@ -193,6 +244,13 @@ namespace SCTL_NAMESPACE {
         return lhs & (~rhs);
       }
 
+      // Bitshift
+      friend IntegerVec operator<<(const Vec& lhs, const Integer& rhs) {
+        IntegerVec r = IntegerVec::LoadAligned(&lhs.v[0]);
+        for (Integer i = 0; i < N; i++) r.v[i] = r.v[i] << rhs;
+        return r;
+      }
+
       // Assignment operators
       Vec& operator+=(const Vec& rhs) {
         for (Integer i = 0; i < N; i++) v[i] += rhs.v[i];
@@ -225,6 +283,9 @@ namespace SCTL_NAMESPACE {
         return *this;
       }
 
+      // Conversion operators
+      // /
+
       // Other operators
       friend Vec max(Vec lhs, const Vec& rhs) {
         for (Integer i = 0; i < N; i++) {
@@ -251,6 +312,8 @@ namespace SCTL_NAMESPACE {
         return r;
       }
 
+      template <class Vec1, class Vec2> friend Vec1 reinterpret(const Vec2& x);
+
     private:
 
       static const ValueType const_zero() {
@@ -274,6 +337,11 @@ namespace SCTL_NAMESPACE {
   };
 
   // Other operators
+  template <class RetVec, class Vec> RetVec reinterpret(const Vec& v){
+    static_assert(sizeof(RetVec) == sizeof(Vec));
+    RetVec& r = *(RetVec*)&v;
+    return r;
+  }
   template <class RealVec, class IntVec> RealVec ConvertInt2Real(const IntVec& x) {
     typedef typename RealVec::ScalarType Real;
     typedef typename IntVec::ScalarType Int;
@@ -287,18 +355,20 @@ namespace SCTL_NAMESPACE {
     IntVec l(x + IntVec(Cint));
     return *(RealVec*)&l - RealVec(Creal);
   }
-  template <class IntVec, class RealVec> IntVec RoundReal2Int(const RealVec& x) {
-    typedef typename RealVec::ScalarType Real;
-    typedef typename IntVec::ScalarType Int;
-    assert(sizeof(RealVec) == sizeof(IntVec));
-    assert(sizeof(Real) == sizeof(Int));
-    static constexpr Integer SigBits = TypeTraits<Real>::SigBits;
+  template <class Vec> typename Vec::IntegerVec RoundReal2Int(const Vec& x) {
+    using IntegerType = typename Vec::IntegerType;
+    using RealType = typename Vec::RealType;
+    using IntegerVec = typename Vec::IntegerVec;
+    using RealVec = typename Vec::RealVec;
+    static_assert(std::is_same<RealVec,Vec>::value, "RoundReal2Int: expected real input argument!");
+
+    static constexpr Integer SigBits = TypeTraits<RealType>::SigBits;
     union {
-      Int Cint = (1UL << (SigBits - 1)) + ((SigBits + ((1UL<<(sizeof(Real)*8 - SigBits - 2))-1)) << SigBits);
-      Real Creal;
+      IntegerType Cint = (1UL << (SigBits - 1)) + ((SigBits + ((1UL<<(sizeof(RealType)*8 - SigBits - 2))-1)) << SigBits);
+      RealType Creal;
     };
-    RealVec d(x + RealVec(Creal));
-    return *(IntVec*)&d - IntVec(Cint);
+    RealVec d = x + RealVec(Creal);
+    return reinterpret<IntegerVec>(d) - IntegerVec(Cint);
   }
   template <class Vec> Vec RoundReal2Real(const Vec& x) {
     typedef typename Vec::ScalarType Real;
@@ -408,8 +478,8 @@ namespace SCTL_NAMESPACE {
       Real real_two;
     };
     Vec x_offset(real_zero);
-    Vec xAnd1 = (((x_+x_offset) & Vec(real_one)) == x_offset);
-    Vec xAnd2 = (((x_+x_offset) & Vec(real_two)) == x_offset);
+    auto xAnd1 = (((x_+x_offset) & Vec(real_one)) == x_offset);
+    auto xAnd2 = (((x_+x_offset) & Vec(real_two)) == x_offset);
 
     Vec s2 = AndNot( c1,xAnd1) | (s1 & xAnd1);
     Vec c2 = AndNot(-s1,xAnd1) | (c1 & xAnd1);
@@ -419,6 +489,95 @@ namespace SCTL_NAMESPACE {
     sinx = s3;
     cosx = c3;
   }
+  template <class Vec> void exp_intrin(Vec& expx, const Vec& x) {
+    constexpr Integer ORDER = 10;
+    using IntegerType = typename Vec::IntegerType;
+    using RealType = typename Vec::RealType;
+    using IntegerVec = typename Vec::IntegerVec;
+    using RealVec = typename Vec::RealVec;
+    static_assert(std::is_same<Vec,RealVec>::value, "exp_intrin: expected a real argument");
+
+    using Real = typename RealVec::ScalarType;
+    static constexpr Integer SigBits = TypeTraits<Real>::SigBits;
+    static constexpr Real coeff2  = 1/(((Real)2));
+    static constexpr Real coeff3  = 1/(((Real)2)*3);
+    static constexpr Real coeff4  = 1/(((Real)2)*3*4);
+    static constexpr Real coeff5  = 1/(((Real)2)*3*4*5);
+    static constexpr Real coeff6  = 1/(((Real)2)*3*4*5*6);
+    static constexpr Real coeff7  = 1/(((Real)2)*3*4*5*6*7);
+    static constexpr Real coeff8  = 1/(((Real)2)*3*4*5*6*7*8);
+    static constexpr Real coeff9  = 1/(((Real)2)*3*4*5*6*7*8*9);
+    static constexpr Real coeff10 = 1/(((Real)2)*3*4*5*6*7*8*9*10);
+    static constexpr Real x0 = (Real)0.693147180559945309417232121458l; // ln(2)
+    static constexpr Real invx0 = 1 / x0;
+
+    RealVec x_ = RoundReal2Real(x * invx0); // 4.5 - cycles
+    IntegerVec int_x_ = RoundReal2Int<RealVec>(x_);
+    RealVec x1 = x - x_ * x0; // 2 - cycles
+    RealVec x2, x3, x4, x5, x6, x7, x8, x9, x10;
+
+    RealVec e1 = 1.0 + x1;
+    if (ORDER >= 2) {
+      x2 = x1 * x1;
+      e1 += x2 * coeff2;
+    }
+    if (ORDER >= 3) {
+      x3 = x2 * x1;
+      e1 += x3 * coeff3;
+    }
+    if (ORDER >= 4) {
+      x4 = x2 * x2;
+      e1 += x4 * coeff4;
+    }
+    if (ORDER >= 5) {
+      x5 = x3 * x2;
+      e1 += x5 * coeff5;
+    }
+    if (ORDER >= 6) {
+      x6 = x3 * x3;
+      e1 += x6 * coeff6;
+    }
+    if (ORDER >= 7) {
+      x7 = x4 * x3;
+      e1 += x7 * coeff7;
+    }
+    if (ORDER >= 8) {
+      x8 = x4 * x4;
+      e1 += x8 * coeff8;
+    }
+    if (ORDER >= 9) {
+      x9 = x5 * x4;
+      e1 += x9 * coeff9;
+    }
+    if (ORDER >= 10) {
+      x10 = x5 * x5;
+      e1 += x10 * coeff10;
+    }
+
+    RealVec e2;
+    { // set e2 = 2 ^ x_
+      union {
+        RealType real_one = 1.0;
+        IntegerType int_one;
+      };
+      //__m256i int_e2 = _mm256_add_epi64(
+      //                                  _mm256_set1_epi64x(int_one),
+      //                                  _mm256_slli_epi64(
+      //                                                      _mm256_load_si256((__m256i const*)&int_x_),
+      //                                                      SigBits
+      //                                                    )
+      //                                  ); // int_e2 = int_one + (int_x_ << SigBits);
+      IntegerVec int_e2 = IntegerVec(int_one) + (int_x_ << SigBits);
+
+      // Handle underflow
+      static constexpr IntegerType max_exp = -(IntegerType)(1UL<<((sizeof(Real)*8-SigBits-2)));
+      int_e2 &= (int_x_ > IntegerVec(max_exp));
+
+      e2 = RealVec::LoadAligned((RealType*)&int_e2);
+    }
+
+    expx = e1 * e2;
+  }
 
 #ifdef __AVX__
   template <> class alignas(sizeof(double)*4) Vec<double,4> {
@@ -427,6 +586,10 @@ namespace SCTL_NAMESPACE {
     static constexpr Integer N = 4;
     public:
 
+      typedef typename GetType<DataType::Integer,TypeTraits<ValueType>::Size>::ValueType IntegerType;
+      typedef typename GetType<DataType::Real,TypeTraits<ValueType>::Size>::ValueType RealType;
+      typedef Vec<IntegerType,N> IntegerVec;
+      typedef Vec<RealType,N> RealVec;
       typedef ValueType ScalarType;
 
       static constexpr Integer Size() {
@@ -455,7 +618,7 @@ namespace SCTL_NAMESPACE {
         return r;
       }
 
-      Vec() {}
+      Vec() = default;
 
       Vec(const ValueType& a) {
         v = _mm256_set1_pd(a);
@@ -485,12 +648,8 @@ namespace SCTL_NAMESPACE {
       }
 
       // C-style cast
-      template <class RetValueType> explicit operator Vec<RetValueType,N>() const {
-        Vec<RetValueType,N> r;
-        VecType& ret_v = *(VecType*)&r.v;
-        ret_v = v;
-        return r;
-      }
+      //template <class RetValueType> explicit operator Vec<RetValueType,N>() const {
+      //}
 
       // Arithmetic operators
       friend Vec operator*(Vec lhs, const Vec& rhs) {
@@ -534,7 +693,6 @@ namespace SCTL_NAMESPACE {
       friend Vec operator==(Vec lhs, const Vec& rhs) {
         lhs.v = _mm256_cmp_pd(lhs.v, rhs.v, _CMP_EQ_OS);
         return lhs;
-        return lhs;
       }
       friend Vec operator!=(Vec lhs, const Vec& rhs) {
         lhs.v = _mm256_cmp_pd(lhs.v, rhs.v, _CMP_NEQ_OS);
@@ -610,14 +768,25 @@ namespace SCTL_NAMESPACE {
         return r;
       }
 
+      template <class Vec1, class Vec2> friend Vec1 reinterpret(const Vec2& x);
       template <class Vec> friend Vec RoundReal2Real(const Vec& x);
       template <class Vec> friend void sincos_intrin(Vec& sinx, Vec& cosx, const Vec& x);
+      template <class Vec> friend void exp_intrin(Vec& expx, const Vec& x);
 
     private:
 
       VecType v;
   };
 
+  template <> inline Vec<int64_t,4> reinterpret<Vec<int64_t,4>,Vec<double,4>>(const Vec<double,4>& x){
+    union {
+      Vec<int64_t,4> r;
+      __m256i y;
+    };
+    y = _mm256_castpd_si256(x.v);
+    return r;
+  }
+
   template <> inline Vec<double,4> RoundReal2Real(const Vec<double,4>& x) {
     Vec<double,4> r;
     r.v = _mm256_round_pd(x.v,_MM_FROUND_TO_NEAREST_INT |_MM_FROUND_NO_EXC);
@@ -629,6 +798,270 @@ namespace SCTL_NAMESPACE {
     sinx.v = _mm256_sin_pd(x.v);
     cosx.v = _mm256_cos_pd(x.v);
   }
+
+  template <> inline void exp_intrin(Vec<double,4>& expx, const Vec<double,4>& x) {
+    expx.v = _mm256_exp_pd(x.v);
+  }
+  #endif
+
+#endif
+
+#ifdef __AVX512F__
+  template <> class alignas(sizeof(double)*8) Vec<double,8> {
+    typedef __m512d VecType;
+    typedef double ValueType;
+    static constexpr Integer N = 8;
+    public:
+
+      typedef typename GetType<DataType::Integer,TypeTraits<ValueType>::Size>::ValueType IntegerType;
+      typedef typename GetType<DataType::Real,TypeTraits<ValueType>::Size>::ValueType RealType;
+      typedef Vec<IntegerType,N> IntegerVec;
+      typedef Vec<RealType,N> RealVec;
+      typedef ValueType ScalarType;
+
+      static constexpr Integer Size() {
+        return N;
+      }
+
+      static Vec Zero() {
+        Vec r;
+        r.v = _mm512_setzero_pd();
+        return r;
+      }
+
+      static Vec Load1(ValueType const* p) {
+        Vec r;
+        // TODO: different from _m256d, could make it faster?
+        // r.v = _mm512_broadcast_f64x4(_mm256_broadcast_sd(p));
+        r.v = _mm512_set1_pd(*p);
+        return r;
+      }
+      static Vec Load(ValueType const* p) {
+        Vec r;
+        r.v = _mm512_loadu_pd(p);
+        return r;
+      }
+      static Vec LoadAligned(ValueType const* p) {
+        Vec r;
+        r.v = _mm512_load_pd(p);
+        return r;
+      }
+
+      Vec() = default;
+
+      Vec(const ValueType& a) {
+        v = _mm512_set1_pd(a);
+      }
+
+      Vec(const __mmask8& a) = delete; // disallow implicit conversions
+
+      void Store(ValueType* p) const {
+        _mm512_storeu_pd(p, v);
+      }
+      void StoreAligned(ValueType* p) const {
+        _mm512_store_pd(p, v);
+      }
+
+      // Bitwise NOT
+      Vec operator~() const {
+        Vec r;
+        static constexpr ScalarType Creal = -1.0;
+        r.v = _mm512_xor_pd(v, _mm512_set1_pd(Creal));
+        return r;
+      }
+
+      // Unary plus and minus
+      Vec operator+() const {
+        return *this;
+      }
+      Vec operator-() const {
+        return Zero() - (*this);
+      }
+
+      // C-style cast
+      //template <class RetValueType> explicit operator Vec<RetValueType,N>() const {
+      //}
+
+      // Arithmetic operators
+      friend Vec operator*(Vec lhs, const Vec& rhs) {
+        lhs.v = _mm512_mul_pd(lhs.v, rhs.v);
+        return lhs;
+      }
+      friend Vec operator+(Vec lhs, const Vec& rhs) {
+        lhs.v = _mm512_add_pd(lhs.v, rhs.v);
+        return lhs;
+      }
+      friend Vec operator-(Vec lhs, const Vec& rhs) {
+        lhs.v = _mm512_sub_pd(lhs.v, rhs.v);
+        return lhs;
+      }
+      friend Vec FMA(Vec a, const Vec& b, const Vec& c) {
+        a.v = _mm512_fmadd_pd(a.v, b.v, c.v);
+        //a.v = _mm512_add_pd(_mm512_mul_pd(a.v, b.v), c.v);
+        return a;
+      }
+
+      // Comparison operators
+      //friend Vec operator< (Vec lhs, const Vec& rhs) {
+      //  lhs.v = _mm512_castsi512_pd(_mm512_movm_epi64(_mm512_cmp_pd_mask(lhs.v, rhs.v, _CMP_LT_OS)));
+      //  return lhs;
+      //}
+      //friend Vec operator<=(Vec lhs, const Vec& rhs) {
+      //  lhs.v = _mm512_castsi512_pd(_mm512_movm_epi64(_mm512_cmp_pd_mask(lhs.v, rhs.v, _CMP_LE_OS)));
+      //  return lhs;
+      //}
+      //friend Vec operator>=(Vec lhs, const Vec& rhs) {
+      //  lhs.v = _mm512_castsi512_pd(_mm512_movm_epi64(_mm512_cmp_pd_mask(lhs.v, rhs.v, _CMP_GE_OS)));
+      //  return lhs;
+      //}
+      //friend Vec operator> (Vec lhs, const Vec& rhs) {
+      //  lhs.v = _mm512_castsi512_pd(_mm512_movm_epi64(_mm512_cmp_pd_mask(lhs.v, rhs.v, _CMP_GT_OS)));
+      //  return lhs;
+      //}
+      //friend Vec operator==(Vec lhs, const Vec& rhs) {
+      //  lhs.v = _mm512_castsi512_pd(_mm512_movm_epi64(_mm512_cmp_pd_mask(lhs.v, rhs.v, _CMP_EQ_OS)));
+      //  return lhs;
+      //}
+      //friend Vec operator!=(Vec lhs, const Vec& rhs) {
+      //  lhs.v = _mm512_castsi512_pd(_mm512_movm_epi64(_mm512_cmp_pd_mask(lhs.v, rhs.v, _CMP_NEQ_OS)));
+      //  return lhs;
+      //}
+
+      friend __mmask8 operator< (Vec lhs, const Vec& rhs) {
+        return _mm512_cmp_pd_mask(lhs.v, rhs.v, _CMP_LT_OS);
+      }
+      friend __mmask8 operator<=(Vec lhs, const Vec& rhs) {
+        return _mm512_cmp_pd_mask(lhs.v, rhs.v, _CMP_LE_OS);
+      }
+      friend __mmask8 operator>=(Vec lhs, const Vec& rhs) {
+        return _mm512_cmp_pd_mask(lhs.v, rhs.v, _CMP_GE_OS);
+      }
+      friend __mmask8 operator> (Vec lhs, const Vec& rhs) {
+        return _mm512_cmp_pd_mask(lhs.v, rhs.v, _CMP_GT_OS);
+      }
+      friend __mmask8 operator==(Vec lhs, const Vec& rhs) {
+        return _mm512_cmp_pd_mask(lhs.v, rhs.v, _CMP_EQ_OS);
+      }
+      friend __mmask8 operator!=(Vec lhs, const Vec& rhs) {
+        return _mm512_cmp_pd_mask(lhs.v, rhs.v, _CMP_NEQ_OS);
+      }
+
+      // Bitwise operators
+      friend Vec operator&(Vec lhs, const Vec& rhs) {
+        lhs.v = _mm512_and_pd(lhs.v, rhs.v);
+        return lhs;
+      }
+      friend Vec operator^(Vec lhs, const Vec& rhs) {
+        lhs.v = _mm512_xor_pd(lhs.v, rhs.v);
+        return lhs;
+      }
+      friend Vec operator|(Vec lhs, const Vec& rhs) {
+        lhs.v = _mm512_or_pd(lhs.v, rhs.v);
+        return lhs;
+      }
+      friend Vec AndNot(Vec lhs, const Vec& rhs) {
+        lhs.v = _mm512_andnot_pd(rhs.v, lhs.v);
+        return lhs;
+      }
+      friend Vec operator&(Vec lhs, const __mmask8& rhs) {
+        lhs.v = _mm512_maskz_mov_pd(rhs, lhs.v);
+        return lhs;
+      }
+      friend Vec AndNot(Vec lhs, const __mmask8& rhs) {
+        lhs.v = _mm512_mask_mov_pd(lhs.v, rhs, _mm512_setzero_pd());
+        return lhs;
+      }
+
+      // Assignment operators
+      Vec& operator*=(const Vec& rhs) {
+        v = _mm512_mul_pd(v, rhs.v);
+        return *this;
+      }
+      Vec& operator+=(const Vec& rhs) {
+        v = _mm512_add_pd(v, rhs.v);
+        return *this;
+      }
+      Vec& operator-=(const Vec& rhs) {
+        v = _mm512_sub_pd(v, rhs.v);
+        return *this;
+      }
+      Vec& operator&=(const Vec& rhs) {
+        v = _mm512_and_pd(v, rhs.v);
+        return *this;
+      }
+      Vec& operator^=(const Vec& rhs) {
+        v = _mm512_xor_pd(v, rhs.v);
+        return *this;
+      }
+      Vec& operator|=(const Vec& rhs) {
+        v = _mm512_or_pd(v, rhs.v);
+        return *this;
+      }
+      Vec& operator&=(const __mmask8& rhs) {
+        v = _mm512_maskz_mov_pd(rhs, v);
+        return *this;
+      }
+
+      // Other operators
+      friend Vec max(Vec lhs, const Vec& rhs) {
+        lhs.v = _mm512_max_pd(lhs.v, rhs.v);
+        return lhs;
+      }
+      friend Vec min(Vec lhs, const Vec& rhs) {
+        lhs.v = _mm512_min_pd(lhs.v, rhs.v);
+        return lhs;
+      }
+
+      friend std::ostream& operator<<(std::ostream& os, const Vec& in) {
+        union {
+          VecType vec;
+          ValueType val[N];
+        };
+        vec = in.v;
+        for (Integer i = 0; i < N; i++) os << val[i] << ' ';
+        return os;
+      }
+      friend Vec approx_rsqrt(const Vec& x) {
+        Vec r;
+        r.v = _mm512_cvtps_pd(_mm256_rsqrt_ps(_mm512_cvtpd_ps(x.v)));
+        return r;
+      }
+
+      template <class Vec1, class Vec2> friend Vec1 reinterpret(const Vec2& x);
+      template <class Vec> friend Vec RoundReal2Real(const Vec& x);
+      template <class Vec> friend void sincos_intrin(Vec& sinx, Vec& cosx, const Vec& x);
+      template <class Vec> friend void exp_intrin(Vec& expx, const Vec& x);
+
+    private:
+
+      VecType v;
+  };
+
+  template <> inline Vec<int64_t,8> reinterpret<Vec<int64_t,8>,Vec<double,8>>(const Vec<double,8>& x){
+    union {
+      Vec<int64_t,8> r;
+      __m512i y;
+    };
+    y = _mm512_castpd_si512(x.v);
+    return r;
+  }
+
+  template <> inline Vec<double,8> RoundReal2Real(const Vec<double,8>& x) {
+    Vec<double,8> r;
+    // TODO: need double check
+    r.v = _mm512_roundscale_pd(x.v,_MM_FROUND_TO_NEAREST_INT |_MM_FROUND_NO_EXC);
+    return r;
+  }
+
+  #ifdef SCTL_HAVE_SVML
+  template <> inline void sincos_intrin(Vec<double,8>& sinx, Vec<double,8>& cosx, const Vec<double,8>& x) {
+    sinx.v = _mm512_sin_pd(x.v);
+    cosx.v = _mm512_cos_pd(x.v);
+  }
+
+  template <> inline void exp_intrin(Vec<double,8>& expx, const Vec<double,8>& x) {
+    expx.v = _mm512_exp_pd(x.v);
+  }
   #endif
 
 #endif