diff options
| -rw-r--r-- | .gitignore | 1 | ||||
| -rw-r--r-- | mathfuncs_sqrt.h | 8 | ||||
| -rw-r--r-- | test.cc | 6 | ||||
| -rw-r--r-- | vec_double_avx.h | 184 | ||||
| -rw-r--r-- | vec_double_sse2.h | 433 | ||||
| -rw-r--r-- | vec_float_avx.h | 506 | ||||
| -rw-r--r-- | vec_float_sse2.h | 415 | ||||
| -rw-r--r-- | vecmathlib.h | 10 |
8 files changed, 1460 insertions, 103 deletions
@@ -5,3 +5,4 @@ example old test .ninja_log +OUT diff --git a/mathfuncs_sqrt.h b/mathfuncs_sqrt.h index c1ea3f7..050fc91 100644 --- a/mathfuncs_sqrt.h +++ b/mathfuncs_sqrt.h @@ -57,6 +57,14 @@ namespace vecmathlib { template<typename realvec_t> realvec_t mathfuncs<realvec_t>::vml_rsqrt(realvec_t x) { +#if 0 + double const x_2 = 0.5 * x; // x/2 + // Newton method: + // Solve f(r) = 0 for f(r) = x - 1/r^2 + // r <- r - f(r) / f'(r) + // r <- (3 r - r^3 x) / 2 + return r * (1.5 - r*r * x_2); +#endif return rcp(sqrt(x)); } @@ -40,7 +40,7 @@ struct vecmathlib_test { // Test each function with this many random values - static int const imax = 100000; + static int const imax = 10000; // Require that 3/4 of the digits are correct static real_t constexpr accuracy = pow(realvec_t::epsilon(), R(0.75)); @@ -380,7 +380,11 @@ int main(int argc, char** argv) cout << "Testing math functions:\n"; vecmathlib_test<realvec<float,1>>::test(); + vecmathlib_test<realvec<float,4>>::test(); + vecmathlib_test<realvec<float,8>>::test(); + vecmathlib_test<realvec<double,1>>::test(); + vecmathlib_test<realvec<double,2>>::test(); vecmathlib_test<realvec<double,4>>::test(); cout << "\n"; diff --git a/vec_double_avx.h b/vec_double_avx.h index f311095..081045f 100644 --- a/vec_double_avx.h +++ b/vec_double_avx.h @@ -171,34 +171,27 @@ namespace vecmathlib { // Note: not all arithmetic operations are supported! intvec operator+() const { return *this; } - intvec operator-() const - { - __m128i v01 = _mm256_castsi256_si128(v); - __m128i v23 = _mm256_extractf128_si256(v, 1); - v01 = _mm_sub_epi64(_mm_set1_epi64x(0), v01); - v23 = _mm_sub_epi64(_mm_set1_epi64x(0), v23); - return _mm256_insertf128_si256(_mm256_castsi128_si256(v01), v23, 1); - } + intvec operator-() const { return IV(I(0)) - *this; } intvec operator+(intvec x) const { - __m128i v01 = _mm256_castsi256_si128(v); - __m128i v23 = _mm256_extractf128_si256(v, 1); - __m128i xv01 = _mm256_castsi256_si128(x.v); - __m128i xv23 = _mm256_extractf128_si256(x.v, 1); - v01 = _mm_add_epi64(v01, xv01); - v23 = _mm_add_epi64(v23, xv23); - return _mm256_insertf128_si256(_mm256_castsi128_si256(v01), v23, 1); + __m128i vlo = _mm256_castsi256_si128(v); + __m128i vhi = _mm256_extractf128_si256(v, 1); + __m128i xvlo = _mm256_castsi256_si128(x.v); + __m128i xvhi = _mm256_extractf128_si256(x.v, 1); + vlo = _mm_add_epi64(vlo, xvlo); + vhi = _mm_add_epi64(vhi, xvhi); + return _mm256_insertf128_si256(_mm256_castsi128_si256(vlo), vhi, 1); } intvec operator-(intvec x) const { - __m128i v01 = _mm256_castsi256_si128(v); - __m128i v23 = _mm256_extractf128_si256(v, 1); - __m128i xv01 = _mm256_castsi256_si128(x.v); - __m128i xv23 = _mm256_extractf128_si256(x.v, 1); - v01 = _mm_sub_epi64(v01, xv01); - v23 = _mm_sub_epi64(v23, xv23); - return _mm256_insertf128_si256(_mm256_castsi128_si256(v01), v23, 1); + __m128i vlo = _mm256_castsi256_si128(v); + __m128i vhi = _mm256_extractf128_si256(v, 1); + __m128i xvlo = _mm256_castsi256_si128(x.v); + __m128i xvhi = _mm256_extractf128_si256(x.v, 1); + vlo = _mm_sub_epi64(vlo, xvlo); + vhi = _mm_sub_epi64(vhi, xvhi); + return _mm256_insertf128_si256(_mm256_castsi128_si256(vlo), vhi, 1); } intvec& operator+=(intvec const& x) { return *this=*this+x; } @@ -206,11 +199,7 @@ namespace vecmathlib { - intvec operator~() const - { - return _mm256_castpd_si256(_mm256_xor_pd(_mm256_castsi256_pd(IV(~U(0))), - _mm256_castsi256_pd(v))); - } + intvec operator~() const { return IV(~U(0)) ^ *this; } intvec operator&(intvec x) const { @@ -234,118 +223,109 @@ namespace vecmathlib { - // SSE2 shift instructions potentially relevant for 64-bit - // operations: - // _mm_slli_epi64 - // _mm_srli_epi64 - // _mm_sll_epi64 - // _mm_srl_epi64 - // _mm_srai_epi32 - // _mm_sra_epi32 - // _mm_srli_si128 - // _mm_slli_si128 intvec lsr(int_t n) const { - __m128i v01 = _mm256_castsi256_si128(v); - __m128i v23 = _mm256_extractf128_si256(v, 1); - v01 = _mm_srli_epi64(v01, n); - v23 = _mm_srli_epi64(v23, n); - return _mm256_insertf128_si256(_mm256_castsi128_si256(v01), v23, 1); + __m128i vlo = _mm256_castsi256_si128(v); + __m128i vhi = _mm256_extractf128_si256(v, 1); + vlo = _mm_srli_epi64(vlo, n); + vhi = _mm_srli_epi64(vhi, n); + return _mm256_insertf128_si256(_mm256_castsi128_si256(vlo), vhi, 1); } intvec operator>>(int_t n) const { - __m128i v01 = _mm256_castsi256_si128(v); - __m128i v23 = _mm256_extractf128_si256(v, 1); - // There is no _mm_srai_epi64. To emulate it, add 0x8000000 - // before shifting, and subtracted the shifted 0x80000000 after shifting + __m128i vlo = _mm256_castsi256_si128(v); + __m128i vhi = _mm256_extractf128_si256(v, 1); + // There is no _mm_srai_epi64. To emulate it, add 0x80000000 + // before shifting, and subtract the shifted 0x80000000 after + // shifting #if 0 __m128i signmask01 = _mm_sub_epi64(_mm_set1_epi64x(0), - _mm_srli_epi64(v01, 63)); + _mm_srli_epi64(vlo, 63)); __m128i signmask23 = _mm_sub_epi64(_mm_set1_epi64x(0), - _mm_srli_epi64(v23, 63)); - v01 = _mm_xor_si128(signmask01, v01); - v23 = _mm_xor_si128(signmask23, v23); - v01 = _mm_srli_epi64(v01, n); - v23 = _mm_srli_epi64(v23, n); - v01 = _mm_xor_si128(signmask01, v01); - v23 = _mm_xor_si128(signmask23, v23); + _mm_srli_epi64(vhi, 63)); + vlo = _mm_xor_si128(signmask01, vlo); + vhi = _mm_xor_si128(signmask23, vhi); + vlo = _mm_srli_epi64(vlo, n); + vhi = _mm_srli_epi64(vhi, n); + vlo = _mm_xor_si128(signmask01, vlo); + vhi = _mm_xor_si128(signmask23, vhi); #else // Convert signed to unsiged - v01 += _mm_set1_epi64x(U(1) << (bits-1)); - v23 += _mm_set1_epi64x(U(1) << (bits-1)); + vlo += _mm_set1_epi64x(U(1) << (bits-1)); + vhi += _mm_set1_epi64x(U(1) << (bits-1)); // Shift - v01 = _mm_srli_epi64(v01, n); - v23 = _mm_srli_epi64(v23, n); + vlo = _mm_srli_epi64(vlo, n); + vhi = _mm_srli_epi64(vhi, n); // Undo conversion - v01 -= _mm_set1_epi64x(U(1) << (bits-n)); - v23 -= _mm_set1_epi64x(U(1) << (bits-n)); + vlo -= _mm_set1_epi64x(U(1) << (bits-n)); + vhi -= _mm_set1_epi64x(U(1) << (bits-n)); #endif - return _mm256_insertf128_si256(_mm256_castsi128_si256(v01), v23, 1); + return _mm256_insertf128_si256(_mm256_castsi128_si256(vlo), vhi, 1); } intvec operator<<(int_t n) const { - __m128i v01 = _mm256_castsi256_si128(v); - __m128i v23 = _mm256_extractf128_si256(v, 1); - v01 = _mm_slli_epi64(v01, n); - v23 = _mm_slli_epi64(v23, n); - return _mm256_insertf128_si256(_mm256_castsi128_si256(v01), v23, 1); + __m128i vlo = _mm256_castsi256_si128(v); + __m128i vhi = _mm256_extractf128_si256(v, 1); + vlo = _mm_slli_epi64(vlo, n); + vhi = _mm_slli_epi64(vhi, n); + return _mm256_insertf128_si256(_mm256_castsi128_si256(vlo), vhi, 1); } intvec& operator>>=(int_t n) { return *this=*this>>n; } intvec& operator<<=(int_t n) { return *this=*this<<n; } intvec lsr(intvec n) const { - __m128i v01 = _mm256_castsi256_si128(v); - __m128i v23 = _mm256_extractf128_si256(v, 1); - __m128i nv01 = _mm256_castsi256_si128(n.v); - __m128i nv23 = _mm256_extractf128_si256(n.v, 1); - v01 = _mm_srl_epi64(v01, nv01); - v23 = _mm_srl_epi64(v23, nv23); - return _mm256_insertf128_si256(_mm256_castsi128_si256(v01), v23, 1); + __m128i vlo = _mm256_castsi256_si128(v); + __m128i vhi = _mm256_extractf128_si256(v, 1); + __m128i nvlo = _mm256_castsi256_si128(n.v); + __m128i nvhi = _mm256_extractf128_si256(n.v, 1); + vlo = _mm_srl_epi64(vlo, nvlo); + vhi = _mm_srl_epi64(vhi, nvhi); + return _mm256_insertf128_si256(_mm256_castsi128_si256(vlo), vhi, 1); } intvec operator>>(intvec n) const { - __m128i v01 = _mm256_castsi256_si128(v); - __m128i v23 = _mm256_extractf128_si256(v, 1); - __m128i nv01 = _mm256_castsi256_si128(n.v); - __m128i nv23 = _mm256_extractf128_si256(n.v, 1); + __m128i vlo = _mm256_castsi256_si128(v); + __m128i vhi = _mm256_extractf128_si256(v, 1); + __m128i nvlo = _mm256_castsi256_si128(n.v); + __m128i nvhi = _mm256_extractf128_si256(n.v, 1); #if 0 // There is no _mm_srai_epi64. To emulate it, invert all bits // before and after shifting if the sign bit is set. __m128i signmask01 = _mm_sub_epi64(_mm_set1_epi64x(0), - _mm_srli_epi64(v01, 63)); + _mm_srli_epi64(vlo, 63)); __m128i signmask23 = _mm_sub_epi64(_mm_set1_epi64x(0), - _mm_srli_epi64(v23, 63)); - v01 = _mm_xor_si128(signmask01, v01); - v23 = _mm_xor_si128(signmask23, v23); - v01 = _mm_srl_epi64(v01, nv01); - v23 = _mm_srl_epi64(v23, nv23); - v01 = _mm_xor_si128(signmask01, v01); - v23 = _mm_xor_si128(signmask23, v23); + _mm_srli_epi64(vhi, 63)); + vlo = _mm_xor_si128(signmask01, vlo); + vhi = _mm_xor_si128(signmask23, vhi); + vlo = _mm_srl_epi64(vlo, nvlo); + vhi = _mm_srl_epi64(vhi, nvhi); + vlo = _mm_xor_si128(signmask01, vlo); + vhi = _mm_xor_si128(signmask23, vhi); #else // Convert signed to unsiged - v01 += _mm_set1_epi64x(U(1) << (bits-1)); - v23 += _mm_set1_epi64x(U(1) << (bits-1)); + vlo += _mm_set1_epi64x(U(1) << (bits-1)); + vhi += _mm_set1_epi64x(U(1) << (bits-1)); // Shift - v01 = _mm_srl_epi64(v01, nv01); - v23 = _mm_srl_epi64(v23, nv23); + vlo = _mm_srl_epi64(vlo, nvlo); + vhi = _mm_srl_epi64(vhi, nvhi); // Undo conversion - v01 -= _mm_sll_epi64(_mm_set1_epi64x(1), - _mm_sub_epi64(_mm_set1_epi64x(bits), nv01)); - v23 -= _mm_sll_epi64(_mm_set1_epi64x(1), - _mm_sub_epi64(_mm_set1_epi64x(bits), nv23)); + vlo -= _mm_sll_epi64(_mm_set1_epi64x(1), + _mm_sub_epi64(_mm_set1_epi64x(bits), nvlo)); + vhi -= _mm_sll_epi64(_mm_set1_epi64x(1), + _mm_sub_epi64(_mm_set1_epi64x(bits), nvhi)); #endif - return _mm256_insertf128_si256(_mm256_castsi128_si256(v01), v23, 1); + return _mm256_insertf128_si256(_mm256_castsi128_si256(vlo), vhi, 1); } intvec operator<<(intvec n) const { - __m128i v01 = _mm256_castsi256_si128(v); - __m128i v23 = _mm256_extractf128_si256(v, 1); - __m128i nv01 = _mm256_castsi256_si128(n.v); - __m128i nv23 = _mm256_extractf128_si256(n.v, 1); - v01 = _mm_sll_epi64(v01, nv01); - v23 = _mm_sll_epi64(v23, nv23); - return _mm256_insertf128_si256(_mm256_castsi128_si256(v01), v23, 1); + __m128i vlo = _mm256_castsi256_si128(v); + __m128i vhi = _mm256_extractf128_si256(v, 1); + __m128i nvlo = _mm256_castsi256_si128(n.v); + __m128i nvhi = _mm256_extractf128_si256(n.v, 1); + vlo = _mm_sll_epi64(vlo, nvlo); + vhi = _mm_sll_epi64(vhi, nvhi); + return _mm256_insertf128_si256(_mm256_castsi128_si256(vlo), vhi, 1); } intvec& operator>>=(intvec n) { return *this=*this>>n; } intvec& operator<<=(intvec n) { return *this=*this<<n; } @@ -403,7 +383,7 @@ namespace vecmathlib { realvec operator+() const { return *this; } - realvec operator-() const { return _mm256_sub_pd(_mm256_set1_pd(0.0), v); } + realvec operator-() const { return RV(0.0) - *this; } realvec operator+(realvec x) const { return _mm256_add_pd(v, x.v); } realvec operator-(realvec x) const { return _mm256_sub_pd(v, x.v); } diff --git a/vec_double_sse2.h b/vec_double_sse2.h new file mode 100644 index 0000000..1c6648c --- /dev/null +++ b/vec_double_sse2.h @@ -0,0 +1,433 @@ +// -*-C++-*- + +#ifndef VEC_DOUBLE_SSE2_H +#define VEC_DOUBLE_SSE2_H + +#include "floatprops.h" +#include "mathfuncs.h" +#include "vec_base.h" + +#include <cmath> + +// SSE2 intrinsics +#include <emmintrin.h> +#ifdef __SSE4_1__ // Intel's SSE 4.1 +# include <smmintrin.h> +#endif +#ifdef __SSE4A__ // AMD's SSE 4a +# include <ammintrin.h> +#endif + + + +namespace vecmathlib { + + template<> struct boolvec<double,2>; + template<> struct intvec<double,2>; + template<> struct realvec<double,2>; + + + + template<> + struct boolvec<double,2>: floatprops<double> + { + static int const size = 2; + typedef bool scalar_t; + typedef __m128d bvector_t; + + static_assert(size * sizeof(real_t) == sizeof(bvector_t), + "vector size is wrong"); + + private: + // true values have the sign bit set, false values have it unset + static uint_t from_bool(bool a) { return - uint_t(a); } + static bool to_bool(uint_t a) { return int_t(a) < int_t(0); } + public: + + typedef boolvec boolvec_t; + typedef intvec<real_t, size> intvec_t; + typedef realvec<real_t, size> realvec_t; + + // Short names for type casts + typedef real_t R; + typedef int_t I; + typedef uint_t U; + typedef realvec_t RV; + typedef intvec_t IV; + typedef boolvec_t BV; + typedef floatprops<real_t> FP; + typedef mathfuncs<realvec_t> MF; + + + + bvector_t v; + + boolvec() {} + boolvec(boolvec const& x): v(x.v) {} + boolvec& operator=(boolvec const& x) { return v=x.v, *this; } + boolvec(bvector_t x): v(x) {} + boolvec(bool a): + v(_mm_castsi128_pd(_mm_set1_epi64x(from_bool(a)))) {} + boolvec(bool const* as): + v(_mm_castsi128_pd(_mm_set_epi64x(from_bool(as[1]), from_bool(as[0])))) {} + + operator bvector_t() const { return v; } + bool operator[](int n) const { return to_bool(((uint_t const*)&v)[n]); } + boolvec& set_elt(int n, bool a) + { + return ((int_t*)&v)[n] = from_bool(a), *this; + } + + + + auto as_int() const -> intvec_t; // defined after intvec + auto convert_int() const -> intvec_t; // defined after intvec + + + + boolvec operator!() const { return _mm_xor_pd(boolvec(true), v); } + + boolvec operator&&(boolvec x) const { return _mm_and_pd(v, x.v); } + boolvec operator||(boolvec x) const { return _mm_or_pd(v, x.v); } + boolvec operator==(boolvec x) const { return !(*this==x); } + boolvec operator!=(boolvec x) const { return _mm_xor_pd(v, x.v); } + + bool all() const + { + return (*this)[0] && (*this)[1]; + } + bool any() const + { + return (*this)[0] || (*this)[1]; + } + + + + // ifthen(condition, then-value, else-value) + auto ifthen(intvec_t x, + intvec_t y) const -> intvec_t; // defined after intvec + auto ifthen(realvec_t x, + realvec_t y) const -> realvec_t; // defined after realvec + + }; + + + + template<> + struct intvec<double,2>: floatprops<double> + { + static int const size = 2; + typedef int_t scalar_t; + typedef __m128i ivector_t; + + static_assert(size * sizeof(real_t) == sizeof(ivector_t), + "vector size is wrong"); + + typedef boolvec<real_t, size> boolvec_t; + typedef intvec intvec_t; + typedef realvec<real_t, size> realvec_t; + + // Short names for type casts + typedef real_t R; + typedef int_t I; + typedef uint_t U; + typedef realvec_t RV; + typedef intvec_t IV; + typedef boolvec_t BV; + typedef floatprops<real_t> FP; + typedef mathfuncs<realvec_t> MF; + + + + ivector_t v; + + intvec() {} + intvec(intvec const& x): v(x.v) {} + intvec& operator=(intvec const& x) { return v=x.v, *this; } + intvec(ivector_t x): v(x) {} + intvec(int_t a): v(_mm_set1_epi64x(a)) {} + intvec(int_t const* as): v(_mm_set_epi64x(as[1], as[0])) {} + + operator ivector_t() const { return v; } + int_t operator[](int n) const { return ((int_t const*)&v)[n]; } + intvec& set_elt(int n, int_t a) { return ((int_t*)&v)[n]=a, *this; } + + + + auto as_bool() const -> boolvec_t { return _mm_castsi128_pd(v); } + auto convert_bool() const -> boolvec_t + { + // Result: convert_bool(0)=false, convert_bool(else)=true + // There is no intrinsic to compare with zero. Instead, we check + // whether x is positive and x-1 is negative. + intvec x = *this; + intvec xm1 = x - 1; + // We know that boolvec values depend only on the sign bit + // return (~xm1 | x).as_bool(); + return x.as_bool() || !xm1.as_bool(); + } + auto as_float() const -> realvec_t; // defined after realvec + auto convert_float() const -> realvec_t; // defined after realvec + + + + // Note: not all arithmetic operations are supported! + + intvec operator+() const { return *this; } + intvec operator-() const { return IV(I(0)) - *this; } + + intvec operator+(intvec x) const { return _mm_add_epi64(v, x.v); } + intvec operator-(intvec x) const { return _mm_sub_epi64(v, x.v); } + + intvec& operator+=(intvec const& x) { return *this=*this+x; } + intvec& operator-=(intvec const& x) { return *this=*this-x; } + + + + intvec operator~() const { return IV(~U(0)) ^ *this; } + + intvec operator&(intvec x) const + { + return _mm_castpd_si128(_mm_and_pd(_mm_castsi128_pd(v), + _mm_castsi128_pd(x.v))); + } + intvec operator|(intvec x) const + { + return _mm_castpd_si128(_mm_or_pd(_mm_castsi128_pd(v), + _mm_castsi128_pd(x.v))); + } + intvec operator^(intvec x) const + { + return _mm_castpd_si128(_mm_xor_pd(_mm_castsi128_pd(v), + _mm_castsi128_pd(x.v))); + } + + intvec& operator&=(intvec const& x) { return *this=*this&x; } + intvec& operator|=(intvec const& x) { return *this=*this|x; } + intvec& operator^=(intvec const& x) { return *this=*this^x; } + + + + intvec lsr(int_t n) const { return _mm_srli_epi64(v, n); } + intvec operator>>(int_t n) const + { + // There is no _mm_srai_epi64. To emulate it, add 0x80000000 + // before shifting, and subtract the shifted 0x80000000 after + // shifting + intvec x = *this; + // Convert signed to unsiged + x += U(1) << (bits-1); + // Shift + x = x.lsr(x); + // Undo conversion + x -= U(1) << (bits-n); + return x; + } + intvec operator<<(int_t n) const { return _mm_slli_epi64(v, n); } + intvec& operator>>=(int_t n) { return *this=*this>>n; } + intvec& operator<<=(int_t n) { return *this=*this<<n; } + + intvec lsr(intvec n) const { return _mm_srl_epi64(v, n.v); } + intvec operator>>(intvec n) const + { + // There is no _mm_srai_epi64. To emulate it, add 0x80000000 + // before shifting, and subtract the shifted 0x80000000 after + // shifting + intvec x = *this; + // Convert signed to unsiged + x += U(1) << (bits-1); + // Shift + x = x.lsr(n); + // Undo conversion + x -= IV(1) << (IV(bits) - n); + return x; + } + intvec operator<<(intvec n) const { return _mm_sll_epi64(v, n.v); } + intvec& operator>>=(intvec n) { return *this=*this>>n; } + intvec& operator<<=(intvec n) { return *this=*this<<n; } + }; + + + + template<> + struct realvec<double,2>: floatprops<double> + { + static int const size = 2; + typedef real_t scalar_t; + typedef __m128d vector_t; + + static constexpr char const* const name = "<SSE2:2*double>"; + inline void barrier() { asm("": "+x" (v)); } + + static_assert(size * sizeof(real_t) == sizeof(vector_t), + "vector size is wrong"); + + typedef boolvec<real_t, size> boolvec_t; + typedef intvec<real_t, size> intvec_t; + typedef realvec realvec_t; + + // Short names for type casts + typedef real_t R; + typedef int_t I; + typedef uint_t U; + typedef realvec_t RV; + typedef intvec_t IV; + typedef boolvec_t BV; + typedef floatprops<real_t> FP; + typedef mathfuncs<realvec_t> MF; + + + + vector_t v; + + realvec() {} + realvec(realvec const& x): v(x.v) {} + realvec& operator=(realvec const& x) { return v=x.v, *this; } + realvec(vector_t x): v(x) {} + realvec(real_t a): v(_mm_set1_pd(a)) {} + realvec(real_t const* as): v(_mm_set_pd(as[1], as[0])) {} + + operator vector_t() const { return v; } + real_t operator[](int n) const { return ((real_t const*)&v)[n]; } + realvec& set_elt(int n, real_t a) { return ((real_t*)&v)[n]=a, *this; } + + + + intvec_t as_int() const { return _mm_castpd_si128(v); } + intvec_t convert_int() const { return MF::vml_convert_int(*this); } + + + + realvec operator+() const { return *this; } + realvec operator-() const { return RV(0.0) - *this; } + + realvec operator+(realvec x) const { return _mm_add_pd(v, x.v); } + realvec operator-(realvec x) const { return _mm_sub_pd(v, x.v); } + realvec operator*(realvec x) const { return _mm_mul_pd(v, x.v); } + realvec operator/(realvec x) const { return _mm_div_pd(v, x.v); } + + realvec& operator+=(realvec const& x) { return *this=*this+x; } + realvec& operator-=(realvec const& x) { return *this=*this-x; } + realvec& operator*=(realvec const& x) { return *this=*this*x; } + realvec& operator/=(realvec const& x) { return *this=*this/x; } + + real_t prod() const + { + return (*this)[0] * (*this)[1]; + } + real_t sum() const + { + return (*this)[0] + (*this)[1]; + } + + + + boolvec_t operator==(realvec const& x) const + { + return _mm_cmpeq_pd(v, x.v); + } + boolvec_t operator!=(realvec const& x) const + { + return _mm_cmpneq_pd(v, x.v); + } + boolvec_t operator<(realvec const& x) const + { + return _mm_cmplt_pd(v, x.v); + } + boolvec_t operator<=(realvec const& x) const + { + return _mm_cmple_pd(v, x.v); + } + boolvec_t operator>(realvec const& x) const + { + return _mm_cmpgt_pd(v, x.v); + } + boolvec_t operator>=(realvec const& x) const + { + return _mm_cmpge_pd(v, x.v); + } + + + + realvec acos() const { return MF::vml_acos(*this); } + realvec acosh() const { return MF::vml_acosh(*this); } + realvec asin() const { return MF::vml_asin(*this); } + realvec asinh() const { return MF::vml_asinh(*this); } + realvec atan() const { return MF::vml_atan(*this); } + realvec atanh() const { return MF::vml_atanh(*this); } + realvec ceil() const { return _mm_ceil_pd(v); } + realvec copysign(realvec y) const { return MF::vml_copysign(*this, y); } + realvec cos() const { return MF::vml_cos(*this); } + realvec cosh() const { return MF::vml_cosh(*this); } + realvec exp() const { return MF::vml_exp(*this); } + realvec exp10() const { return MF::vml_exp10(*this); } + realvec exp2() const { return MF::vml_exp2(*this); } + realvec expm1() const { return MF::vml_expm1(*this); } + realvec fabs() const { return MF::vml_fabs(*this); } + realvec floor() const { return _mm_floor_pd(v); } + realvec fmod(realvec y) const { return MF::vml_fmod(*this, y); } + intvec_t ilogb() const { return MF::vml_ilogb(*this); } + realvec log() const { return MF::vml_log(*this); } + realvec log10() const { return MF::vml_log10(*this); } + realvec log1p() const { return MF::vml_log1p(*this); } + realvec log2() const { return MF::vml_log2(*this); } + realvec pow(realvec y) const { return MF::vml_pow(*this, y); } + realvec rcp() const { return _mm_div_pd(_mm_set1_pd(1.0), v); } + realvec remainder(realvec y) const { return MF::vml_remainder(*this, y); } + realvec round() const { return _mm_round_pd(v, _MM_FROUND_NINT); } + realvec rsqrt() const { return MF::vml_rsqrt(*this); } + realvec scalbn(intvec_t n) const { return MF::vml_scalbn(*this, n); } + boolvec_t signbit() const { return v; } + realvec sin() const { return MF::vml_sin(*this); } + realvec sinh() const { return MF::vml_sinh(*this); } + realvec sqrt() const { return _mm_sqrt_pd(v); } + realvec tan() const { return MF::vml_tan(*this); } + realvec tanh() const { return MF::vml_tanh(*this); } + }; + + + + // boolvec definitions + + inline + auto boolvec<double,2>::as_int() const -> intvec_t + { + return _mm_castpd_si128(v); + } + + inline + auto boolvec<double,2>::convert_int() const -> intvec_t + { + //return ifthen(v, U(1), U(0)); + return lsr(as_int(), bits-1); + } + + inline + auto boolvec<double,2>::ifthen(intvec_t x, intvec_t y) const -> intvec_t + { + return ifthen(x.as_float(), y.as_float()).as_int(); + } + + inline + auto boolvec<double,2>::ifthen(realvec_t x, realvec_t y) const -> realvec_t + { + return _mm_blendv_pd(y.v, x.v, v); + } + + + + // intvec definitions + + inline auto intvec<double,2>::as_float() const -> realvec_t + { + return _mm_castsi128_pd(v); + } + + inline auto intvec<double,2>::convert_float() const -> realvec_t + { + return MF::vml_convert_float(*this); + } + +} // namespace vecmathlib + +#endif // #ifndef VEC_DOUBLE_SSE2_H diff --git a/vec_float_avx.h b/vec_float_avx.h new file mode 100644 index 0000000..4a19673 --- /dev/null +++ b/vec_float_avx.h @@ -0,0 +1,506 @@ +// -*-C++-*- + +#ifndef VEC_FLOAT_AVX_H +#define VEC_FLOAT_AVX_H + +#include "floatprops.h" +#include "mathfuncs.h" +#include "vec_base.h" + +#include <cmath> + +// AVX intrinsics +#include <immintrin.h> + + + +namespace vecmathlib { + + template<> struct boolvec<float,8>; + template<> struct intvec<float,8>; + template<> struct realvec<float,8>; + + + + template<> + struct boolvec<float,8>: floatprops<float> + { + static int const size = 8; + typedef bool scalar_t; + typedef __m256 bvector_t; + + static_assert(size * sizeof(real_t) == sizeof(bvector_t), + "vector size is wrong"); + + private: + // true values have the sign bit set, false values have it unset + static uint_t from_bool(bool a) { return - uint_t(a); } + static bool to_bool(uint_t a) { return int_t(a) < int_t(0); } + public: + + typedef boolvec boolvec_t; + typedef intvec<real_t, size> intvec_t; + typedef realvec<real_t, size> realvec_t; + + // Short names for type casts + typedef real_t R; + typedef int_t I; + typedef uint_t U; + typedef realvec_t RV; + typedef intvec_t IV; + typedef boolvec_t BV; + typedef floatprops<real_t> FP; + typedef mathfuncs<realvec_t> MF; + + + + bvector_t v; + + boolvec() {} + boolvec(boolvec const& x): v(x.v) {} + boolvec& operator=(boolvec const& x) { return v=x.v, *this; } + boolvec(bvector_t x): v(x) {} + boolvec(bool a): + v(_mm256_castsi256_ps(_mm256_set1_epi32(from_bool(a)))) {} + boolvec(bool const* as): + v(_mm256_castsi256_ps(_mm256_set_epi32(from_bool(as[7]), + from_bool(as[6]), + from_bool(as[5]), + from_bool(as[4]), + from_bool(as[3]), + from_bool(as[2]), + from_bool(as[1]), + from_bool(as[0])))) {} + + operator bvector_t() const { return v; } + bool operator[](int n) const { return to_bool(((uint_t const*)&v)[n]); } + boolvec& set_elt(int n, bool a) + { + return ((int_t*)&v)[n] = from_bool(a), *this; + } + + + + auto as_int() const -> intvec_t; // defined after intvec + auto convert_int() const -> intvec_t; // defined after intvec + + + + boolvec operator!() const { return _mm256_xor_ps(boolvec(true), v); } + + boolvec operator&&(boolvec x) const { return _mm256_and_ps(v, x.v); } + boolvec operator||(boolvec x) const { return _mm256_or_ps(v, x.v); } + boolvec operator==(boolvec x) const { return !(*this==x); } + boolvec operator!=(boolvec x) const { return _mm256_xor_ps(v, x.v); } + + bool all() const + { + return (*this)[0] && (*this)[1] && (*this)[2] && (*this)[3]; + } + bool any() const + { + return (*this)[0] || (*this)[1] || (*this)[2] || (*this)[3]; + } + + + + // ifthen(condition, then-value, else-value) + auto ifthen(intvec_t x, + intvec_t y) const -> intvec_t; // defined after intvec + auto ifthen(realvec_t x, + realvec_t y) const -> realvec_t; // defined after realvec + + }; + + + + template<> + struct intvec<float,8>: floatprops<float> + { + static int const size = 8; + typedef int_t scalar_t; + typedef __m256i ivector_t; + + static_assert(size * sizeof(real_t) == sizeof(ivector_t), + "vector size is wrong"); + + typedef boolvec<real_t, size> boolvec_t; + typedef intvec intvec_t; + typedef realvec<real_t, size> realvec_t; + + // Short names for type casts + typedef real_t R; + typedef int_t I; + typedef uint_t U; + typedef realvec_t RV; + typedef intvec_t IV; + typedef boolvec_t BV; + typedef floatprops<real_t> FP; + typedef mathfuncs<realvec_t> MF; + + + + ivector_t v; + + intvec() {} + intvec(intvec const& x): v(x.v) {} + intvec& operator=(intvec const& x) { return v=x.v, *this; } + intvec(ivector_t x): v(x) {} + intvec(int_t a): v(_mm256_set1_epi32(a)) {} + intvec(int_t const* as): v(_mm256_set_epi32(as[7], as[6], as[5], as[4], + as[3], as[2], as[1], as[0])) {} + + operator ivector_t() const { return v; } + int_t operator[](int n) const { return ((int_t const*)&v)[n]; } + intvec& set_elt(int n, int_t a) { return ((int_t*)&v)[n]=a, *this; } + + + + auto as_bool() const -> boolvec_t { return _mm256_castsi256_ps(v); } + auto convert_bool() const -> boolvec_t + { + // Result: convert_bool(0)=false, convert_bool(else)=true + // There is no intrinsic to compare with zero. Instead, we check + // whether x is positive and x-1 is negative. + intvec x = *this; + intvec xm1 = x - 1; + // We know that boolvec values depend only on the sign bit + // return (~xm1 | x).as_bool(); + return x.as_bool() || !xm1.as_bool(); + } + auto as_float() const -> realvec_t; // defined after realvec + auto convert_float() const -> realvec_t; // defined after realvec + + + + // Note: not all arithmetic operations are supported! + + intvec operator+() const { return *this; } + intvec operator-() const { return IV(0) - *this; } + + intvec operator+(intvec x) const + { + __m128i vlo = _mm256_castsi256_si128(v); + __m128i vhi = _mm256_extractf128_si256(v, 1); + __m128i xvlo = _mm256_castsi256_si128(x.v); + __m128i xvhi = _mm256_extractf128_si256(x.v, 1); + vlo = _mm_add_epi32(vlo, xvlo); + vhi = _mm_add_epi32(vhi, xvhi); + return _mm256_insertf128_si256(_mm256_castsi128_si256(vlo), vhi, 1); + } + intvec operator-(intvec x) const + { + __m128i vlo = _mm256_castsi256_si128(v); + __m128i vhi = _mm256_extractf128_si256(v, 1); + __m128i xvlo = _mm256_castsi256_si128(x.v); + __m128i xvhi = _mm256_extractf128_si256(x.v, 1); + vlo = _mm_sub_epi32(vlo, xvlo); + vhi = _mm_sub_epi32(vhi, xvhi); + return _mm256_insertf128_si256(_mm256_castsi128_si256(vlo), vhi, 1); + } + + intvec& operator+=(intvec const& x) { return *this=*this+x; } + intvec& operator-=(intvec const& x) { return *this=*this-x; } + + + + intvec operator~() const { return IV(~U(0)) ^ *this; } + + intvec operator&(intvec x) const + { + return _mm256_castps_si256(_mm256_and_ps(_mm256_castsi256_ps(v), + _mm256_castsi256_ps(x.v))); + } + intvec operator|(intvec x) const + { + return _mm256_castps_si256(_mm256_or_ps(_mm256_castsi256_ps(v), + _mm256_castsi256_ps(x.v))); + } + intvec operator^(intvec x) const + { + return _mm256_castps_si256(_mm256_xor_ps(_mm256_castsi256_ps(v), + _mm256_castsi256_ps(x.v))); + } + + intvec& operator&=(intvec const& x) { return *this=*this&x; } + intvec& operator|=(intvec const& x) { return *this=*this|x; } + intvec& operator^=(intvec const& x) { return *this=*this^x; } + + + + intvec lsr(int_t n) const + { + __m128i vlo = _mm256_castsi256_si128(v); + __m128i vhi = _mm256_extractf128_si256(v, 1); + vlo = _mm_srli_epi32(vlo, n); + vhi = _mm_srli_epi32(vhi, n); + return _mm256_insertf128_si256(_mm256_castsi128_si256(vlo), vhi, 1); + } + intvec operator>>(int_t n) const + { + __m128i vlo = _mm256_castsi256_si128(v); + __m128i vhi = _mm256_extractf128_si256(v, 1); + // There is no _mm_srai_epi32. To emulate it, add 0x80000000 + // before shifting, and subtract the shifted 0x80000000 after + // shifting + // Convert signed to unsiged + vlo += _mm_set1_epi32(U(1) << (bits-1)); + vhi += _mm_set1_epi32(U(1) << (bits-1)); + // Shift + vlo = _mm_srli_epi32(vlo, n); + vhi = _mm_srli_epi32(vhi, n); + // Undo conversion + vlo -= _mm_set1_epi32(U(1) << (bits-n)); + vhi -= _mm_set1_epi32(U(1) << (bits-n)); + return _mm256_insertf128_si256(_mm256_castsi128_si256(vlo), vhi, 1); + } + intvec operator<<(int_t n) const + { + __m128i vlo = _mm256_castsi256_si128(v); + __m128i vhi = _mm256_extractf128_si256(v, 1); + vlo = _mm_slli_epi32(vlo, n); + vhi = _mm_slli_epi32(vhi, n); + return _mm256_insertf128_si256(_mm256_castsi128_si256(vlo), vhi, 1); + } + intvec& operator>>=(int_t n) { return *this=*this>>n; } + intvec& operator<<=(int_t n) { return *this=*this<<n; } + + intvec lsr(intvec n) const + { + __m128i vlo = _mm256_castsi256_si128(v); + __m128i vhi = _mm256_extractf128_si256(v, 1); + __m128i nvlo = _mm256_castsi256_si128(n.v); + __m128i nvhi = _mm256_extractf128_si256(n.v, 1); + vlo = _mm_srl_epi32(vlo, nvlo); + vhi = _mm_srl_epi32(vhi, nvhi); + return _mm256_insertf128_si256(_mm256_castsi128_si256(vlo), vhi, 1); + } + intvec operator>>(intvec n) const + { + __m128i vlo = _mm256_castsi256_si128(v); + __m128i vhi = _mm256_extractf128_si256(v, 1); + __m128i nvlo = _mm256_castsi256_si128(n.v); + __m128i nvhi = _mm256_extractf128_si256(n.v, 1); + // Convert signed to unsiged + vlo += _mm_set1_epi32(U(1) << (bits-1)); + vhi += _mm_set1_epi32(U(1) << (bits-1)); + // Shift + vlo = _mm_srl_epi32(vlo, nvlo); + vhi = _mm_srl_epi32(vhi, nvhi); + // Undo conversion + vlo -= _mm_sll_epi32(_mm_set1_epi32(1), + _mm_sub_epi32(_mm_set1_epi32(bits), nvlo)); + vhi -= _mm_sll_epi32(_mm_set1_epi32(1), + _mm_sub_epi32(_mm_set1_epi32(bits), nvhi)); + return _mm256_insertf128_si256(_mm256_castsi128_si256(vlo), vhi, 1); + } + intvec operator<<(intvec n) const + { + __m128i vlo = _mm256_castsi256_si128(v); + __m128i vhi = _mm256_extractf128_si256(v, 1); + __m128i nvlo = _mm256_castsi256_si128(n.v); + __m128i nvhi = _mm256_extractf128_si256(n.v, 1); + vlo = _mm_sll_epi32(vlo, nvlo); + vhi = _mm_sll_epi32(vhi, nvhi); + return _mm256_insertf128_si256(_mm256_castsi128_si256(vlo), vhi, 1); + } + intvec& operator>>=(intvec n) { return *this=*this>>n; } + intvec& operator<<=(intvec n) { return *this=*this<<n; } + }; + + + + template<> + struct realvec<float,8>: floatprops<float> + { + static int const size = 8; + typedef real_t scalar_t; + typedef __m256 vector_t; + + static constexpr char const* const name = "<AVX:8*float>"; + inline void barrier() { asm("": "+x" (v)); } + + static_assert(size * sizeof(real_t) == sizeof(vector_t), + "vector size is wrong"); + + typedef boolvec<real_t, size> boolvec_t; + typedef intvec<real_t, size> intvec_t; + typedef realvec realvec_t; + + // Short names for type casts + typedef real_t R; + typedef int_t I; + typedef uint_t U; + typedef realvec_t RV; + typedef intvec_t IV; + typedef boolvec_t BV; + typedef floatprops<real_t> FP; + typedef mathfuncs<realvec_t> MF; + + + + vector_t v; + + realvec() {} + realvec(realvec const& x): v(x.v) {} + realvec& operator=(realvec const& x) { return v=x.v, *this; } + realvec(vector_t x): v(x) {} + realvec(real_t a): v(_mm256_set1_ps(a)) {} + realvec(real_t const* as): v(_mm256_set_ps(as[7], as[6], as[5], as[4], + as[3], as[2], as[1], as[0])) {} + + operator vector_t() const { return v; } + real_t operator[](int n) const { return ((real_t const*)&v)[n]; } + realvec& set_elt(int n, real_t a) { return ((real_t*)&v)[n]=a, *this; } + + + + intvec_t as_int() const { return _mm256_castps_si256(v); } + intvec_t convert_int() const { return _mm256_cvtps_epi32(v); } + + + + realvec operator+() const { return *this; } + realvec operator-() const { return RV(0.0) - *this; } + + realvec operator+(realvec x) const { return _mm256_add_ps(v, x.v); } + realvec operator-(realvec x) const { return _mm256_sub_ps(v, x.v); } + realvec operator*(realvec x) const { return _mm256_mul_ps(v, x.v); } + realvec operator/(realvec x) const { return _mm256_div_ps(v, x.v); } + + realvec& operator+=(realvec const& x) { return *this=*this+x; } + realvec& operator-=(realvec const& x) { return *this=*this-x; } + realvec& operator*=(realvec const& x) { return *this=*this*x; } + realvec& operator/=(realvec const& x) { return *this=*this/x; } + + real_t prod() const + { + return (*this)[0] * (*this)[1] * (*this)[2] * (*this)[3]; + } + real_t sum() const + { + return (*this)[0] + (*this)[1] + (*this)[2] + (*this)[3]; + } + + + + boolvec_t operator==(realvec const& x) const + { + return _mm256_cmp_ps(v, x.v, _CMP_EQ_OQ); + } + boolvec_t operator!=(realvec const& x) const + { + return _mm256_cmp_ps(v, x.v, _CMP_NEQ_OQ); + } + boolvec_t operator<(realvec const& x) const + { + return _mm256_cmp_ps(v, x.v, _CMP_LT_OQ); + } + boolvec_t operator<=(realvec const& x) const + { + return _mm256_cmp_ps(v, x.v, _CMP_LE_OQ); + } + boolvec_t operator>(realvec const& x) const + { + return _mm256_cmp_ps(v, x.v, _CMP_GT_OQ); + } + boolvec_t operator>=(realvec const& x) const + { + return _mm256_cmp_ps(v, x.v, _CMP_GE_OQ); + } + + + + realvec acos() const { return MF::vml_acos(*this); } + realvec acosh() const { return MF::vml_acosh(*this); } + realvec asin() const { return MF::vml_asin(*this); } + realvec asinh() const { return MF::vml_asinh(*this); } + realvec atan() const { return MF::vml_atan(*this); } + realvec atanh() const { return MF::vml_atanh(*this); } + realvec ceil() const { return _mm256_ceil_ps(v); } + realvec copysign(realvec y) const { return MF::vml_copysign(*this, y); } + realvec cos() const { return MF::vml_cos(*this); } + realvec cosh() const { return MF::vml_cosh(*this); } + realvec exp() const { return MF::vml_exp(*this); } + realvec exp10() const { return MF::vml_exp10(*this); } + realvec exp2() const { return MF::vml_exp2(*this); } + realvec expm1() const { return MF::vml_expm1(*this); } + realvec fabs() const { return MF::vml_fabs(*this); } + realvec floor() const { return _mm256_floor_ps(v); } + realvec fmod(realvec y) const { return MF::vml_fmod(*this, y); } + intvec_t ilogb() const { return MF::vml_ilogb(*this); } + realvec log() const { return MF::vml_log(*this); } + realvec log10() const { return MF::vml_log10(*this); } + realvec log1p() const { return MF::vml_log1p(*this); } + realvec log2() const { return MF::vml_log2(*this); } + realvec pow(realvec y) const { return MF::vml_pow(*this, y); } + realvec rcp() const + { + realvec x = *this; + realvec r = _mm256_rcp_ps(x); // this is only an approximation + r *= RV(2.0) - r*x; // one Newton iteration (see vml_rcp) + return r; + } + realvec remainder(realvec y) const { return MF::vml_remainder(*this, y); } + realvec round() const { return _mm256_round_ps(v, _MM_FROUND_NINT); } + realvec rsqrt() const + { + realvec x = *this; + realvec r = _mm256_rsqrt_ps(x); // this is only an approximation + r *= RV(1.5) - RV(0.5)*x * r*r; // one Newton iteration (see vml_rsqrt) + return r; + } + realvec scalbn(intvec_t n) const { return MF::vml_scalbn(*this, n); } + boolvec_t signbit() const { return v; } + realvec sin() const { return MF::vml_sin(*this); } + realvec sinh() const { return MF::vml_sinh(*this); } + realvec sqrt() const { return _mm256_sqrt_ps(v); } + realvec tan() const { return MF::vml_tan(*this); } + realvec tanh() const { return MF::vml_tanh(*this); } + }; + + + + // boolvec definitions + + inline + auto boolvec<float,8>::as_int() const -> intvec_t + { + return _mm256_castps_si256(v); + } + + inline + auto boolvec<float,8>::convert_int() const -> intvec_t + { + return lsr(as_int(), bits-1); + } + + inline + auto boolvec<float,8>::ifthen(intvec_t x, intvec_t y) const -> intvec_t + { + return ifthen(x.as_float(), y.as_float()).as_int(); + } + + inline + auto boolvec<float,8>::ifthen(realvec_t x, realvec_t y) const -> realvec_t + { + return _mm256_blendv_ps(y.v, x.v, v); + } + + + + // intvec definitions + + inline auto intvec<float,8>::as_float() const -> realvec_t + { + return _mm256_castsi256_ps(v); + } + + inline auto intvec<float,8>::convert_float() const -> realvec_t + { + return _mm256_cvtepi32_ps(v); + } + +} // namespace vecmathlib + +#endif // #ifndef VEC_FLOAT_AVX_H diff --git a/vec_float_sse2.h b/vec_float_sse2.h new file mode 100644 index 0000000..f99eae8 --- /dev/null +++ b/vec_float_sse2.h @@ -0,0 +1,415 @@ +// -*-C++-*- + +#ifndef VEC_FLOAT_SSE2_H +#define VEC_FLOAT_SSE2_H + +#include "floatprops.h" +#include "mathfuncs.h" +#include "vec_base.h" + +#include <cmath> + +// SSE2 intrinsics +#include <xmmintrin.h> +#if defined __SSE4_1__ // Intel's SSE 4.1 +# include <smmintrin.h> +#endif +#if defined __SSE4A__ // AMD's SSE 4a +# include <ammintrin.h> +#endif + + + +namespace vecmathlib { + + template<> struct boolvec<float,4>; + template<> struct intvec<float,4>; + template<> struct realvec<float,4>; + + + + template<> + struct boolvec<float,4>: floatprops<float> + { + static int const size = 4; + typedef bool scalar_t; + typedef __m128 bvector_t; + + static_assert(size * sizeof(real_t) == sizeof(bvector_t), + "vector size is wrong"); + + private: + // true values have the sign bit set, false values have it unset + static uint_t from_bool(bool a) { return - uint_t(a); } + static bool to_bool(uint_t a) { return int_t(a) < int_t(0); } + public: + + typedef boolvec boolvec_t; + typedef intvec<real_t, size> intvec_t; + typedef realvec<real_t, size> realvec_t; + + // Short names for type casts + typedef real_t R; + typedef int_t I; + typedef uint_t U; + typedef realvec_t RV; + typedef intvec_t IV; + typedef boolvec_t BV; + typedef floatprops<real_t> FP; + typedef mathfuncs<realvec_t> MF; + + + + bvector_t v; + + boolvec() {} + boolvec(boolvec const& x): v(x.v) {} + boolvec& operator=(boolvec const& x) { return v=x.v, *this; } + boolvec(bvector_t x): v(x) {} + boolvec(bool a): + v(_mm_castsi128_ps(_mm_set1_epi32(from_bool(a)))) {} + boolvec(bool const* as): + v(_mm_castsi128_ps(_mm_set_epi32(from_bool(as[3]), + from_bool(as[2]), + from_bool(as[1]), + from_bool(as[0])))) {} + + operator bvector_t() const { return v; } + bool operator[](int n) const { return to_bool(((uint_t const*)&v)[n]); } + boolvec& set_elt(int n, bool a) + { + return ((int_t*)&v)[n] = from_bool(a), *this; + } + + + + auto as_int() const -> intvec_t; // defined after intvec + auto convert_int() const -> intvec_t; // defined after intvec + + + + boolvec operator!() const { return _mm_xor_ps(boolvec(true), v); } + + boolvec operator&&(boolvec x) const { return _mm_and_ps(v, x.v); } + boolvec operator||(boolvec x) const { return _mm_or_ps(v, x.v); } + boolvec operator==(boolvec x) const { return !(*this==x); } + boolvec operator!=(boolvec x) const { return _mm_xor_ps(v, x.v); } + + bool all() const + { + return (*this)[0] && (*this)[1] && (*this)[2] && (*this)[3]; + } + bool any() const + { + return (*this)[0] || (*this)[1] || (*this)[2] || (*this)[3]; + } + + + + // ifthen(condition, then-value, else-value) + auto ifthen(intvec_t x, + intvec_t y) const -> intvec_t; // defined after intvec + auto ifthen(realvec_t x, + realvec_t y) const -> realvec_t; // defined after realvec + + }; + + + + template<> + struct intvec<float,4>: floatprops<float> + { + static int const size = 4; + typedef int_t scalar_t; + typedef __m128i ivector_t; + + static_assert(size * sizeof(real_t) == sizeof(ivector_t), + "vector size is wrong"); + + typedef boolvec<real_t, size> boolvec_t; + typedef intvec intvec_t; + typedef realvec<real_t, size> realvec_t; + + // Short names for type casts + typedef real_t R; + typedef int_t I; + typedef uint_t U; + typedef realvec_t RV; + typedef intvec_t IV; + typedef boolvec_t BV; + typedef floatprops<real_t> FP; + typedef mathfuncs<realvec_t> MF; + + + + ivector_t v; + + intvec() {} + intvec(intvec const& x): v(x.v) {} + intvec& operator=(intvec const& x) { return v=x.v, *this; } + intvec(ivector_t x): v(x) {} + intvec(int_t a): v(_mm_set1_epi32(a)) {} + intvec(int_t const* as): v(_mm_set_epi32(as[3], as[2], as[1], as[0])) {} + + operator ivector_t() const { return v; } + int_t operator[](int n) const { return ((int_t const*)&v)[n]; } + intvec& set_elt(int n, int_t a) { return ((int_t*)&v)[n]=a, *this; } + + + + auto as_bool() const -> boolvec_t { return _mm_castsi128_ps(v); } + auto convert_bool() const -> boolvec_t + { + // Result: convert_bool(0)=false, convert_bool(else)=true + return ! IV(_mm_cmpeq_epi32(v, IV(0))).as_bool(); + } + auto as_float() const -> realvec_t; // defined after realvec + auto convert_float() const -> realvec_t; // defined after realvec + + + + // Note: not all arithmetic operations are supported! + + intvec operator+() const { return *this; } + intvec operator-() const { return IV(0) - *this; } + + intvec operator+(intvec x) const { return _mm_add_epi32(v, x.v); } + intvec operator-(intvec x) const { return _mm_sub_epi32(v, x.v); } + + intvec& operator+=(intvec const& x) { return *this=*this+x; } + intvec& operator-=(intvec const& x) { return *this=*this-x; } + + + + intvec operator~() const { return IV(~U(0)) ^ *this; } + + intvec operator&(intvec x) const + { + return _mm_castps_si128(_mm_and_ps(_mm_castsi128_ps(v), + _mm_castsi128_ps(x.v))); + } + intvec operator|(intvec x) const + { + return _mm_castps_si128(_mm_or_ps(_mm_castsi128_ps(v), + _mm_castsi128_ps(x.v))); + } + intvec operator^(intvec x) const + { + return _mm_castps_si128(_mm_xor_ps(_mm_castsi128_ps(v), + _mm_castsi128_ps(x.v))); + } + + intvec& operator&=(intvec const& x) { return *this=*this&x; } + intvec& operator|=(intvec const& x) { return *this=*this|x; } + intvec& operator^=(intvec const& x) { return *this=*this^x; } + + + + intvec lsr(int_t n) const { return _mm_srli_epi32(v, n); } + intvec operator>>(int_t n) const { return _mm_srai_epi32(v, n); } + intvec operator<<(int_t n) const { return _mm_slli_epi32(v, n); } + intvec& operator>>=(int_t n) { return *this=*this>>n; } + intvec& operator<<=(int_t n) { return *this=*this<<n; } + + intvec lsr(intvec n) const { return _mm_srl_epi32(v, n.v); } + intvec operator>>(intvec n) const { return _mm_sra_epi32(v, n.v); } + intvec operator<<(intvec n) const { return _mm_sll_epi32(v, n.v); } + intvec& operator>>=(intvec n) { return *this=*this>>n; } + intvec& operator<<=(intvec n) { return *this=*this<<n; } + }; + + + + template<> + struct realvec<float,4>: floatprops<float> + { + static int const size = 4; + typedef real_t scalar_t; + typedef __m128 vector_t; + + static constexpr char const* const name = "<SSE2:4*float>"; + inline void barrier() { asm("": "+x" (v)); } + + static_assert(size * sizeof(real_t) == sizeof(vector_t), + "vector size is wrong"); + + typedef boolvec<real_t, size> boolvec_t; + typedef intvec<real_t, size> intvec_t; + typedef realvec realvec_t; + + // Short names for type casts + typedef real_t R; + typedef int_t I; + typedef uint_t U; + typedef realvec_t RV; + typedef intvec_t IV; + typedef boolvec_t BV; + typedef floatprops<real_t> FP; + typedef mathfuncs<realvec_t> MF; + + + + vector_t v; + + realvec() {} + realvec(realvec const& x): v(x.v) {} + realvec& operator=(realvec const& x) { return v=x.v, *this; } + realvec(vector_t x): v(x) {} + realvec(real_t a): v(_mm_set1_ps(a)) {} + realvec(real_t const* as): v(_mm_set_ps(as[3], as[2], as[1], as[0])) {} + + operator vector_t() const { return v; } + real_t operator[](int n) const { return ((real_t const*)&v)[n]; } + realvec& set_elt(int n, real_t a) { return ((real_t*)&v)[n]=a, *this; } + + + + intvec_t as_int() const { return _mm_castps_si128(v); } + intvec_t convert_int() const { return _mm_cvtps_epi32(v); } + + + + realvec operator+() const { return *this; } + realvec operator-() const { return RV(0.0) - *this; } + + realvec operator+(realvec x) const { return _mm_add_ps(v, x.v); } + realvec operator-(realvec x) const { return _mm_sub_ps(v, x.v); } + realvec operator*(realvec x) const { return _mm_mul_ps(v, x.v); } + realvec operator/(realvec x) const { return _mm_div_ps(v, x.v); } + + realvec& operator+=(realvec const& x) { return *this=*this+x; } + realvec& operator-=(realvec const& x) { return *this=*this-x; } + realvec& operator*=(realvec const& x) { return *this=*this*x; } + realvec& operator/=(realvec const& x) { return *this=*this/x; } + + real_t prod() const + { + return (*this)[0] * (*this)[1] * (*this)[2] * (*this)[3]; + } + real_t sum() const + { + return (*this)[0] + (*this)[1] + (*this)[2] + (*this)[3]; + } + + + + boolvec_t operator==(realvec const& x) const + { + return _mm_cmpeq_ps(v, x.v); + } + boolvec_t operator!=(realvec const& x) const + { + return _mm_cmpneq_ps(v, x.v); + } + boolvec_t operator<(realvec const& x) const + { + return _mm_cmplt_ps(v, x.v); + } + boolvec_t operator<=(realvec const& x) const + { + return _mm_cmple_ps(v, x.v); + } + boolvec_t operator>(realvec const& x) const + { + return _mm_cmpgt_ps(v, x.v); + } + boolvec_t operator>=(realvec const& x) const + { + return _mm_cmpge_ps(v, x.v); + } + + + + realvec acos() const { return MF::vml_acos(*this); } + realvec acosh() const { return MF::vml_acosh(*this); } + realvec asin() const { return MF::vml_asin(*this); } + realvec asinh() const { return MF::vml_asinh(*this); } + realvec atan() const { return MF::vml_atan(*this); } + realvec atanh() const { return MF::vml_atanh(*this); } + realvec ceil() const { return _mm_ceil_ps(v); } + realvec copysign(realvec y) const { return MF::vml_copysign(*this, y); } + realvec cos() const { return MF::vml_cos(*this); } + realvec cosh() const { return MF::vml_cosh(*this); } + realvec exp() const { return MF::vml_exp(*this); } + realvec exp10() const { return MF::vml_exp10(*this); } + realvec exp2() const { return MF::vml_exp2(*this); } + realvec expm1() const { return MF::vml_expm1(*this); } + realvec fabs() const { return MF::vml_fabs(*this); } + realvec floor() const { return _mm_floor_ps(v); } + realvec fmod(realvec y) const { return MF::vml_fmod(*this, y); } + intvec_t ilogb() const { return MF::vml_ilogb(*this); } + realvec log() const { return MF::vml_log(*this); } + realvec log10() const { return MF::vml_log10(*this); } + realvec log1p() const { return MF::vml_log1p(*this); } + realvec log2() const { return MF::vml_log2(*this); } + realvec pow(realvec y) const { return MF::vml_pow(*this, y); } + realvec rcp() const + { + realvec x = *this; + realvec r = _mm_rcp_ps(x); // this is only an approximation + r *= RV(2.0) - r*x; // one Newton iteration (see vml_rcp) + return r; + } + realvec remainder(realvec y) const { return MF::vml_remainder(*this, y); } + realvec round() const { return _mm_round_ps(v, _MM_FROUND_NINT); } + realvec rsqrt() const + { + realvec x = *this; + realvec r = _mm_rsqrt_ps(x); // this is only an approximation + r *= RV(1.5) - RV(0.5)*x * r*r; // one Newton iteration (see vml_rsqrt) + return r; + } + realvec scalbn(intvec_t n) const { return MF::vml_scalbn(*this, n); } + boolvec_t signbit() const { return v; } + realvec sin() const { return MF::vml_sin(*this); } + realvec sinh() const { return MF::vml_sinh(*this); } + realvec sqrt() const { return _mm_sqrt_ps(v); } + realvec tan() const { return MF::vml_tan(*this); } + realvec tanh() const { return MF::vml_tanh(*this); } + }; + + + + // boolvec definitions + + inline + auto boolvec<float,4>::as_int() const -> intvec_t + { + return _mm_castps_si128(v); + } + + inline + auto boolvec<float,4>::convert_int() const -> intvec_t + { + return lsr(as_int(), bits-1); + } + + inline + auto boolvec<float,4>::ifthen(intvec_t x, intvec_t y) const -> intvec_t + { + return ifthen(x.as_float(), y.as_float()).as_int(); + } + + inline + auto boolvec<float,4>::ifthen(realvec_t x, realvec_t y) const -> realvec_t + { + return _mm_blendv_ps(y.v, x.v, v); + } + + + + // intvec definitions + + inline auto intvec<float,4>::as_float() const -> realvec_t + { + return _mm_castsi128_ps(v); + } + + inline auto intvec<float,4>::convert_float() const -> realvec_t + { + return _mm_cvtepi32_ps(v); + } + +} // namespace vecmathlib + +#endif // #ifndef VEC_FLOAT_SSE2_H diff --git a/vecmathlib.h b/vecmathlib.h index b5703bb..83396ad 100644 --- a/vecmathlib.h +++ b/vecmathlib.h @@ -4,7 +4,17 @@ #define VECMATHLIB_H #include "vec_float.h" +#if defined __SSE2__ // Intel SSE 2 +# include "vec_float_sse2.h" +#endif +#if defined __AVX__ // Intel AVX +# include "vec_float_avx.h" +#endif + #include "vec_double.h" +#if defined __SSE2__ // Intel SSE 2 +# include "vec_double_sse2.h" +#endif #if defined __AVX__ // Intel AVX # include "vec_double_avx.h" #endif |
