From 4f11a22311b8db3130eace2478c553050a2b7ec0 Mon Sep 17 00:00:00 2001
From: Ahmet Inan <inan@aicodix.de>
Date: Tue, 11 Dec 2018 09:05:12 +0100
Subject: [PATCH] added SIMD wrappers for ARM NEON, Intel SSE4.1 and AVX2

---
 README.md |    7 +
 avx2.hh   |  990 ++++++++++++++++++++++++++++++++++++++++++++
 neon.hh   |  843 +++++++++++++++++++++++++++++++++++++
 simd.hh   | 1189 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 sse4_1.hh |  975 +++++++++++++++++++++++++++++++++++++++++++
 5 files changed, 4004 insertions(+)
 create mode 100644 avx2.hh
 create mode 100644 neon.hh
 create mode 100644 simd.hh
 create mode 100644 sse4_1.hh

diff --git a/README.md b/README.md
index 2f292bf..967da60 100644
--- a/README.md
+++ b/README.md
@@ -65,3 +65,10 @@ It computes the following, but having only O(N) complexity and using O(1) extra
 				output[i] = op(output[i], input[j]);
 ```
 
+### [simd.hh](simd.hh)
+
+Single instruction, multiple data ([SIMD](https://en.wikipedia.org/wiki/SIMD)) wrappers for:
+* [ARM NEON](https://en.wikipedia.org/wiki/ARM_architecture#Advanced_SIMD_(NEON)) ([neon.hh](neon.hh))
+* [Intel SSE4.1](https://en.wikipedia.org/wiki/SSE4) ([sse4_1.hh](sse4_1.hh))
+* [Intel AVX2](https://en.wikipedia.org/wiki/Advanced_Vector_Extensions) ([avx2.hh](avx2.hh))
+
diff --git a/avx2.hh b/avx2.hh
new file mode 100644
index 0000000..281db56
--- /dev/null
+++ b/avx2.hh
@@ -0,0 +1,990 @@
+/*
+Intel AVX2 acceleration
+
+Copyright 2018 Ahmet Inan <inan@aicodix.de>
+*/
+
+#ifndef AVX2_HH
+#define AVX2_HH
+
+#include <immintrin.h>
+
+template <>
+union SIMD<float, 8>
+{
+	static const int SIZE = 8;
+	typedef float value_type;
+	typedef uint32_t uint_type;
+	__m256 m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<double, 4>
+{
+	static const int SIZE = 4;
+	typedef double value_type;
+	typedef uint64_t uint_type;
+	__m256d m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<int8_t, 32>
+{
+	static const int SIZE = 32;
+	typedef int8_t value_type;
+	typedef uint8_t uint_type;
+	__m256i m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<int16_t, 16>
+{
+	static const int SIZE = 16;
+	typedef int16_t value_type;
+	typedef uint16_t uint_type;
+	__m256i m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<int32_t, 8>
+{
+	static const int SIZE = 8;
+	typedef int32_t value_type;
+	typedef uint32_t uint_type;
+	__m256i m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<int64_t, 4>
+{
+	static const int SIZE = 4;
+	typedef int64_t value_type;
+	typedef uint64_t uint_type;
+	__m256i m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<uint8_t, 32>
+{
+	static const int SIZE = 32;
+	typedef uint8_t value_type;
+	typedef uint8_t uint_type;
+	__m256i m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<uint16_t, 16>
+{
+	static const int SIZE = 16;
+	typedef uint16_t value_type;
+	typedef uint16_t uint_type;
+	__m256i m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<uint32_t, 8>
+{
+	static const int SIZE = 8;
+	typedef uint32_t value_type;
+	typedef uint32_t uint_type;
+	__m256i m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<uint64_t, 4>
+{
+	static const int SIZE = 4;
+	typedef uint64_t value_type;
+	typedef uint64_t uint_type;
+	__m256i m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+inline SIMD<float, 8> vreinterpret(SIMD<uint32_t, 8> a)
+{
+	SIMD<float, 8> tmp;
+	tmp.m = (__m256)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 8> vreinterpret(SIMD<float, 8> a)
+{
+	SIMD<uint32_t, 8> tmp;
+	tmp.m = (__m256i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 4> vreinterpret(SIMD<uint64_t, 4> a)
+{
+	SIMD<double, 4> tmp;
+	tmp.m = (__m256d)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 4> vreinterpret(SIMD<double, 4> a)
+{
+	SIMD<uint64_t, 4> tmp;
+	tmp.m = (__m256i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 32> vreinterpret(SIMD<int8_t, 32> a)
+{
+	SIMD<uint8_t, 32> tmp;
+	tmp.m = (__m256i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 32> vreinterpret(SIMD<uint8_t, 32> a)
+{
+	SIMD<int8_t, 32> tmp;
+	tmp.m = (__m256i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 16> vreinterpret(SIMD<int16_t, 16> a)
+{
+	SIMD<uint16_t, 16> tmp;
+	tmp.m = (__m256i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 16> vreinterpret(SIMD<uint16_t, 16> a)
+{
+	SIMD<int16_t, 16> tmp;
+	tmp.m = (__m256i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 8> vreinterpret(SIMD<int32_t, 8> a)
+{
+	SIMD<uint32_t, 8> tmp;
+	tmp.m = (__m256i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 8> vreinterpret(SIMD<uint32_t, 8> a)
+{
+	SIMD<int32_t, 8> tmp;
+	tmp.m = (__m256i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 4> vreinterpret(SIMD<int64_t, 4> a)
+{
+	SIMD<uint64_t, 4> tmp;
+	tmp.m = (__m256i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<int64_t, 4> vreinterpret(SIMD<uint64_t, 4> a)
+{
+	SIMD<int64_t, 4> tmp;
+	tmp.m = (__m256i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 8> vdup<SIMD<float, 8>>(float a)
+{
+	SIMD<float, 8> tmp;
+	tmp.m = _mm256_set1_ps(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 4> vdup<SIMD<double, 4>>(double a)
+{
+	SIMD<double, 4> tmp;
+	tmp.m = _mm256_set1_pd(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 32> vdup<SIMD<int8_t, 32>>(int8_t a)
+{
+	SIMD<int8_t, 32> tmp;
+	tmp.m = _mm256_set1_epi8(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 16> vdup<SIMD<int16_t, 16>>(int16_t a)
+{
+	SIMD<int16_t, 16> tmp;
+	tmp.m = _mm256_set1_epi16(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 8> vdup<SIMD<int32_t, 8>>(int32_t a)
+{
+	SIMD<int32_t, 8> tmp;
+	tmp.m = _mm256_set1_epi32(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<int64_t, 4> vdup<SIMD<int64_t, 4>>(int64_t a)
+{
+	SIMD<int64_t, 4> tmp;
+	tmp.m = _mm256_set1_epi64x(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 8> vzero()
+{
+	SIMD<float, 8> tmp;
+	tmp.m = _mm256_setzero_ps();
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 4> vzero()
+{
+	SIMD<double, 4> tmp;
+	tmp.m = _mm256_setzero_pd();
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 32> vzero()
+{
+	SIMD<int8_t, 32> tmp;
+	tmp.m = _mm256_setzero_si256();
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 16> vzero()
+{
+	SIMD<int16_t, 16> tmp;
+	tmp.m = _mm256_setzero_si256();
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 8> vzero()
+{
+	SIMD<int32_t, 8> tmp;
+	tmp.m = _mm256_setzero_si256();
+	return tmp;
+}
+
+template <>
+inline SIMD<int64_t, 4> vzero()
+{
+	SIMD<int64_t, 4> tmp;
+	tmp.m = _mm256_setzero_si256();
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 8> vadd(SIMD<float, 8> a, SIMD<float, 8> b)
+{
+	SIMD<float, 8> tmp;
+	tmp.m = _mm256_add_ps(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 4> vadd(SIMD<double, 4> a, SIMD<double, 4> b)
+{
+	SIMD<double, 4> tmp;
+	tmp.m = _mm256_add_pd(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 32> vadd(SIMD<int8_t, 32> a, SIMD<int8_t, 32> b)
+{
+	SIMD<int8_t, 32> tmp;
+	tmp.m = _mm256_add_epi8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 16> vadd(SIMD<int16_t, 16> a, SIMD<int16_t, 16> b)
+{
+	SIMD<int16_t, 16> tmp;
+	tmp.m = _mm256_add_epi16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 8> vadd(SIMD<int32_t, 8> a, SIMD<int32_t, 8> b)
+{
+	SIMD<int32_t, 8> tmp;
+	tmp.m = _mm256_add_epi32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int64_t, 4> vadd(SIMD<int64_t, 4> a, SIMD<int64_t, 4> b)
+{
+	SIMD<int64_t, 4> tmp;
+	tmp.m = _mm256_add_epi64(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 32> vqadd(SIMD<int8_t, 32> a, SIMD<int8_t, 32> b)
+{
+	SIMD<int8_t, 32> tmp;
+	tmp.m = _mm256_adds_epi8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 16> vqadd(SIMD<int16_t, 16> a, SIMD<int16_t, 16> b)
+{
+	SIMD<int16_t, 16> tmp;
+	tmp.m = _mm256_adds_epi16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 8> vsub(SIMD<float, 8> a, SIMD<float, 8> b)
+{
+	SIMD<float, 8> tmp;
+	tmp.m = _mm256_sub_ps(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 4> vsub(SIMD<double, 4> a, SIMD<double, 4> b)
+{
+	SIMD<double, 4> tmp;
+	tmp.m = _mm256_sub_pd(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 32> vsub(SIMD<int8_t, 32> a, SIMD<int8_t, 32> b)
+{
+	SIMD<int8_t, 32> tmp;
+	tmp.m = _mm256_sub_epi8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 16> vsub(SIMD<int16_t, 16> a, SIMD<int16_t, 16> b)
+{
+	SIMD<int16_t, 16> tmp;
+	tmp.m = _mm256_sub_epi16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 8> vsub(SIMD<int32_t, 8> a, SIMD<int32_t, 8> b)
+{
+	SIMD<int32_t, 8> tmp;
+	tmp.m = _mm256_sub_epi32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int64_t, 4> vsub(SIMD<int64_t, 4> a, SIMD<int64_t, 4> b)
+{
+	SIMD<int64_t, 4> tmp;
+	tmp.m = _mm256_sub_epi64(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 32> vqsub(SIMD<int8_t, 32> a, SIMD<int8_t, 32> b)
+{
+	SIMD<int8_t, 32> tmp;
+	tmp.m = _mm256_subs_epi8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 16> vqsub(SIMD<int16_t, 16> a, SIMD<int16_t, 16> b)
+{
+	SIMD<int16_t, 16> tmp;
+	tmp.m = _mm256_subs_epi16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 32> vqsub(SIMD<uint8_t, 32> a, SIMD<uint8_t, 32> b)
+{
+	SIMD<uint8_t, 32> tmp;
+	tmp.m = _mm256_subs_epu8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 16> vqsub(SIMD<uint16_t, 16> a, SIMD<uint16_t, 16> b)
+{
+	SIMD<uint16_t, 16> tmp;
+	tmp.m = _mm256_subs_epu16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 8> vabs(SIMD<float, 8> a)
+{
+	SIMD<float, 8> tmp;
+	tmp.m = _mm256_andnot_ps(_mm256_set1_ps(-0.f), a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 4> vabs(SIMD<double, 4> a)
+{
+	SIMD<double, 4> tmp;
+	tmp.m = _mm256_andnot_pd(_mm256_set1_pd(-0.), a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 32> vqabs(SIMD<int8_t, 32> a)
+{
+	SIMD<int8_t, 32> tmp;
+	tmp.m = _mm256_abs_epi8(_mm256_max_epi8(a.m, _mm256_set1_epi8(-INT8_MAX)));
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 16> vqabs(SIMD<int16_t, 16> a)
+{
+	SIMD<int16_t, 16> tmp;
+	tmp.m = _mm256_abs_epi16(_mm256_max_epi16(a.m, _mm256_set1_epi16(-INT16_MAX)));
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 8> vqabs(SIMD<int32_t, 8> a)
+{
+	SIMD<int32_t, 8> tmp;
+	tmp.m = _mm256_abs_epi32(_mm256_max_epi32(a.m, _mm256_set1_epi32(-INT32_MAX)));
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 8> vsign(SIMD<float, 8> a, SIMD<float, 8> b)
+{
+	SIMD<float, 8> tmp;
+	tmp.m = _mm256_andnot_ps(
+		_mm256_cmp_ps(b.m, _mm256_setzero_ps(), _CMP_EQ_OQ),
+		_mm256_xor_ps(a.m, _mm256_and_ps(_mm256_set1_ps(-0.f), b.m)));
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 4> vsign(SIMD<double, 4> a, SIMD<double, 4> b)
+{
+	SIMD<double, 4> tmp;
+	tmp.m = _mm256_andnot_pd(
+		_mm256_cmp_pd(b.m, _mm256_setzero_pd(), _CMP_EQ_OQ),
+		_mm256_xor_pd(a.m, _mm256_and_pd(_mm256_set1_pd(-0.), b.m)));
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 32> vsign(SIMD<int8_t, 32> a, SIMD<int8_t, 32> b)
+{
+	SIMD<int8_t, 32> tmp;
+	tmp.m = _mm256_sign_epi8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 16> vsign(SIMD<int16_t, 16> a, SIMD<int16_t, 16> b)
+{
+	SIMD<int16_t, 16> tmp;
+	tmp.m = _mm256_sign_epi16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 8> vsign(SIMD<int32_t, 8> a, SIMD<int32_t, 8> b)
+{
+	SIMD<int32_t, 8> tmp;
+	tmp.m = _mm256_sign_epi32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 32> vorr(SIMD<uint8_t, 32> a, SIMD<uint8_t, 32> b)
+{
+	SIMD<uint8_t, 32> tmp;
+	tmp.m = _mm256_or_si256(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 16> vorr(SIMD<uint16_t, 16> a, SIMD<uint16_t, 16> b)
+{
+	SIMD<uint16_t, 16> tmp;
+	tmp.m = _mm256_or_si256(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 8> vorr(SIMD<uint32_t, 8> a, SIMD<uint32_t, 8> b)
+{
+	SIMD<uint32_t, 8> tmp;
+	tmp.m = _mm256_or_si256(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 4> vorr(SIMD<uint64_t, 4> a, SIMD<uint64_t, 4> b)
+{
+	SIMD<uint64_t, 4> tmp;
+	tmp.m = _mm256_or_si256(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 32> vand(SIMD<uint8_t, 32> a, SIMD<uint8_t, 32> b)
+{
+	SIMD<uint8_t, 32> tmp;
+	tmp.m = _mm256_and_si256(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 16> vand(SIMD<uint16_t, 16> a, SIMD<uint16_t, 16> b)
+{
+	SIMD<uint16_t, 16> tmp;
+	tmp.m = _mm256_and_si256(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 8> vand(SIMD<uint32_t, 8> a, SIMD<uint32_t, 8> b)
+{
+	SIMD<uint32_t, 8> tmp;
+	tmp.m = _mm256_and_si256(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 4> vand(SIMD<uint64_t, 4> a, SIMD<uint64_t, 4> b)
+{
+	SIMD<uint64_t, 4> tmp;
+	tmp.m = _mm256_and_si256(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 32> veor(SIMD<uint8_t, 32> a, SIMD<uint8_t, 32> b)
+{
+	SIMD<uint8_t, 32> tmp;
+	tmp.m = _mm256_xor_si256(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 16> veor(SIMD<uint16_t, 16> a, SIMD<uint16_t, 16> b)
+{
+	SIMD<uint16_t, 16> tmp;
+	tmp.m = _mm256_xor_si256(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 8> veor(SIMD<uint32_t, 8> a, SIMD<uint32_t, 8> b)
+{
+	SIMD<uint32_t, 8> tmp;
+	tmp.m = _mm256_xor_si256(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 4> veor(SIMD<uint64_t, 4> a, SIMD<uint64_t, 4> b)
+{
+	SIMD<uint64_t, 4> tmp;
+	tmp.m = _mm256_xor_si256(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 32> vbic(SIMD<uint8_t, 32> a, SIMD<uint8_t, 32> b)
+{
+	SIMD<uint8_t, 32> tmp;
+	tmp.m = _mm256_andnot_si256(b.m, a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 16> vbic(SIMD<uint16_t, 16> a, SIMD<uint16_t, 16> b)
+{
+	SIMD<uint16_t, 16> tmp;
+	tmp.m = _mm256_andnot_si256(b.m, a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 8> vbic(SIMD<uint32_t, 8> a, SIMD<uint32_t, 8> b)
+{
+	SIMD<uint32_t, 8> tmp;
+	tmp.m = _mm256_andnot_si256(b.m, a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 4> vbic(SIMD<uint64_t, 4> a, SIMD<uint64_t, 4> b)
+{
+	SIMD<uint64_t, 4> tmp;
+	tmp.m = _mm256_andnot_si256(b.m, a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 32> vbsl(SIMD<uint8_t, 32> a, SIMD<uint8_t, 32> b, SIMD<uint8_t, 32> c)
+{
+	SIMD<uint8_t, 32> tmp;
+	tmp.m = _mm256_or_si256(_mm256_and_si256(a.m, b.m), _mm256_andnot_si256(a.m, c.m));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 16> vbsl(SIMD<uint16_t, 16> a, SIMD<uint16_t, 16> b, SIMD<uint16_t, 16> c)
+{
+	SIMD<uint16_t, 16> tmp;
+	tmp.m = _mm256_or_si256(_mm256_and_si256(a.m, b.m), _mm256_andnot_si256(a.m, c.m));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 8> vbsl(SIMD<uint32_t, 8> a, SIMD<uint32_t, 8> b, SIMD<uint32_t, 8> c)
+{
+	SIMD<uint32_t, 8> tmp;
+	tmp.m = _mm256_or_si256(_mm256_and_si256(a.m, b.m), _mm256_andnot_si256(a.m, c.m));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 4> vbsl(SIMD<uint64_t, 4> a, SIMD<uint64_t, 4> b, SIMD<uint64_t, 4> c)
+{
+	SIMD<uint64_t, 4> tmp;
+	tmp.m = _mm256_or_si256(_mm256_and_si256(a.m, b.m), _mm256_andnot_si256(a.m, c.m));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 8> vceqz(SIMD<float, 8> a)
+{
+	SIMD<uint32_t, 8> tmp;
+	tmp.m = (__m256i)_mm256_cmp_ps(a.m, _mm256_setzero_ps(), _CMP_EQ_OQ);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 4> vceqz(SIMD<double, 4> a)
+{
+	SIMD<uint64_t, 4> tmp;
+	tmp.m = (__m256i)_mm256_cmp_pd(a.m, _mm256_setzero_pd(), _CMP_EQ_OQ);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 32> vceqz(SIMD<int8_t, 32> a)
+{
+	SIMD<uint8_t, 32> tmp;
+	tmp.m = _mm256_cmpeq_epi8(a.m, _mm256_setzero_si256());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 16> vceqz(SIMD<int16_t, 16> a)
+{
+	SIMD<uint16_t, 16> tmp;
+	tmp.m = _mm256_cmpeq_epi16(a.m, _mm256_setzero_si256());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 8> vceqz(SIMD<int32_t, 8> a)
+{
+	SIMD<uint32_t, 8> tmp;
+	tmp.m = _mm256_cmpeq_epi32(a.m, _mm256_setzero_si256());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 4> vceqz(SIMD<int64_t, 4> a)
+{
+	SIMD<uint64_t, 4> tmp;
+	tmp.m = _mm256_cmpeq_epi64(a.m, _mm256_setzero_si256());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 8> vceq(SIMD<float, 8> a, SIMD<float, 8> b)
+{
+	SIMD<uint32_t, 8> tmp;
+	tmp.m = (__m256i)_mm256_cmp_ps(a.m, b.m, _CMP_EQ_OQ);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 4> vceq(SIMD<double, 4> a, SIMD<double, 4> b)
+{
+	SIMD<uint64_t, 4> tmp;
+	tmp.m = (__m256i)_mm256_cmp_pd(a.m, b.m, _CMP_EQ_OQ);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 32> vceq(SIMD<int8_t, 32> a, SIMD<int8_t, 32> b)
+{
+	SIMD<uint8_t, 32> tmp;
+	tmp.m = _mm256_cmpeq_epi8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 16> vceq(SIMD<int16_t, 16> a, SIMD<int16_t, 16> b)
+{
+	SIMD<uint16_t, 16> tmp;
+	tmp.m = _mm256_cmpeq_epi16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 8> vceq(SIMD<int32_t, 8> a, SIMD<int32_t, 8> b)
+{
+	SIMD<uint32_t, 8> tmp;
+	tmp.m = _mm256_cmpeq_epi32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 4> vceq(SIMD<int64_t, 4> a, SIMD<int64_t, 4> b)
+{
+	SIMD<uint64_t, 4> tmp;
+	tmp.m = _mm256_cmpeq_epi64(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 8> vcgtz(SIMD<float, 8> a)
+{
+	SIMD<uint32_t, 8> tmp;
+	tmp.m = (__m256i)_mm256_cmp_ps(a.m, _mm256_setzero_ps(), _CMP_GT_OQ);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 4> vcgtz(SIMD<double, 4> a)
+{
+	SIMD<uint64_t, 4> tmp;
+	tmp.m = (__m256i)_mm256_cmp_pd(a.m, _mm256_setzero_pd(), _CMP_GT_OQ);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 32> vcgtz(SIMD<int8_t, 32> a)
+{
+	SIMD<uint8_t, 32> tmp;
+	tmp.m = _mm256_cmpgt_epi8(a.m, _mm256_setzero_si256());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 16> vcgtz(SIMD<int16_t, 16> a)
+{
+	SIMD<uint16_t, 16> tmp;
+	tmp.m = _mm256_cmpgt_epi16(a.m, _mm256_setzero_si256());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 8> vcgtz(SIMD<int32_t, 8> a)
+{
+	SIMD<uint32_t, 8> tmp;
+	tmp.m = _mm256_cmpgt_epi32(a.m, _mm256_setzero_si256());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 4> vcgtz(SIMD<int64_t, 4> a)
+{
+	SIMD<uint64_t, 4> tmp;
+	tmp.m = _mm256_cmpgt_epi64(a.m, _mm256_setzero_si256());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 8> vcltz(SIMD<float, 8> a)
+{
+	SIMD<uint32_t, 8> tmp;
+	tmp.m = (__m256i)_mm256_cmp_ps(a.m, _mm256_setzero_ps(), _CMP_LT_OQ);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 4> vcltz(SIMD<double, 4> a)
+{
+	SIMD<uint64_t, 4> tmp;
+	tmp.m = (__m256i)_mm256_cmp_pd(a.m, _mm256_setzero_pd(), _CMP_LT_OQ);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 32> vcltz(SIMD<int8_t, 32> a)
+{
+	SIMD<uint8_t, 32> tmp;
+	tmp.m = _mm256_cmpgt_epi8(_mm256_setzero_si256(), a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 16> vcltz(SIMD<int16_t, 16> a)
+{
+	SIMD<uint16_t, 16> tmp;
+	tmp.m = _mm256_cmpgt_epi16(_mm256_setzero_si256(), a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 8> vcltz(SIMD<int32_t, 8> a)
+{
+	SIMD<uint32_t, 8> tmp;
+	tmp.m = _mm256_cmpgt_epi32(_mm256_setzero_si256(), a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 4> vcltz(SIMD<int64_t, 4> a)
+{
+	SIMD<uint64_t, 4> tmp;
+	tmp.m = _mm256_cmpgt_epi64(_mm256_setzero_si256(), a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 8> vmin(SIMD<float, 8> a, SIMD<float, 8> b)
+{
+	SIMD<float, 8> tmp;
+	tmp.m = _mm256_min_ps(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 4> vmin(SIMD<double, 4> a, SIMD<double, 4> b)
+{
+	SIMD<double, 4> tmp;
+	tmp.m = _mm256_min_pd(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 32> vmin(SIMD<int8_t, 32> a, SIMD<int8_t, 32> b)
+{
+	SIMD<int8_t, 32> tmp;
+	tmp.m = _mm256_min_epi8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 16> vmin(SIMD<int16_t, 16> a, SIMD<int16_t, 16> b)
+{
+	SIMD<int16_t, 16> tmp;
+	tmp.m = _mm256_min_epi16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 8> vmin(SIMD<int32_t, 8> a, SIMD<int32_t, 8> b)
+{
+	SIMD<int32_t, 8> tmp;
+	tmp.m = _mm256_min_epi32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int64_t, 4> vmin(SIMD<int64_t, 4> a, SIMD<int64_t, 4> b)
+{
+	SIMD<int64_t, 4> tmp;
+	tmp.m = _mm256_min_epi64(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 8> vmax(SIMD<float, 8> a, SIMD<float, 8> b)
+{
+	SIMD<float, 8> tmp;
+	tmp.m = _mm256_max_ps(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 4> vmax(SIMD<double, 4> a, SIMD<double, 4> b)
+{
+	SIMD<double, 4> tmp;
+	tmp.m = _mm256_max_pd(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 32> vmax(SIMD<int8_t, 32> a, SIMD<int8_t, 32> b)
+{
+	SIMD<int8_t, 32> tmp;
+	tmp.m = _mm256_max_epi8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 16> vmax(SIMD<int16_t, 16> a, SIMD<int16_t, 16> b)
+{
+	SIMD<int16_t, 16> tmp;
+	tmp.m = _mm256_max_epi16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 8> vmax(SIMD<int32_t, 8> a, SIMD<int32_t, 8> b)
+{
+	SIMD<int32_t, 8> tmp;
+	tmp.m = _mm256_max_epi32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int64_t, 4> vmax(SIMD<int64_t, 4> a, SIMD<int64_t, 4> b)
+{
+	SIMD<int64_t, 4> tmp;
+	tmp.m = _mm256_max_epi64(a.m, b.m);
+	return tmp;
+}
+
+#endif
diff --git a/neon.hh b/neon.hh
new file mode 100644
index 0000000..ac5aa92
--- /dev/null
+++ b/neon.hh
@@ -0,0 +1,843 @@
+/*
+ARM NEON acceleration
+
+Copyright 2018 Ahmet Inan <inan@aicodix.de>
+*/
+
+#ifndef NEON_HH
+#define NEON_HH
+
+#include <arm_neon.h>
+
+template <>
+union SIMD<float, 4>
+{
+	static const int SIZE = 4;
+	typedef float value_type;
+	typedef uint32_t uint_type;
+	float32x4_t m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<int8_t, 16>
+{
+	static const int SIZE = 16;
+	typedef int8_t value_type;
+	typedef uint8_t uint_type;
+	int8x16_t m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<int16_t, 8>
+{
+	static const int SIZE = 8;
+	typedef int16_t value_type;
+	typedef uint16_t uint_type;
+	int16x8_t m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<int32_t, 4>
+{
+	static const int SIZE = 4;
+	typedef int32_t value_type;
+	typedef uint32_t uint_type;
+	int32x4_t m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<int64_t, 2>
+{
+	static const int SIZE = 2;
+	typedef int64_t value_type;
+	typedef uint64_t uint_type;
+	int64x2_t m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<uint8_t, 16>
+{
+	static const int SIZE = 16;
+	typedef uint8_t value_type;
+	typedef uint8_t uint_type;
+	uint8x16_t m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<uint16_t, 8>
+{
+	static const int SIZE = 8;
+	typedef uint16_t value_type;
+	typedef uint16_t uint_type;
+	uint16x8_t m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<uint32_t, 4>
+{
+	static const int SIZE = 4;
+	typedef uint32_t value_type;
+	typedef uint32_t uint_type;
+	uint32x4_t m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<uint64_t, 2>
+{
+	static const int SIZE = 2;
+	typedef uint64_t value_type;
+	typedef uint64_t uint_type;
+	uint64x2_t m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+inline SIMD<float, 4> vreinterpret(SIMD<uint32_t, 4> a)
+{
+	SIMD<float, 4> tmp;
+	tmp.m = (float32x4_t)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vreinterpret(SIMD<float, 4> a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = (uint32x4_t)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vreinterpret(SIMD<uint8_t, 16> a)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = (int8x16_t)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vreinterpret(SIMD<int8_t, 16> a)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = (uint8x16_t)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vreinterpret(SIMD<uint16_t, 8> a)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = (int16x8_t)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vreinterpret(SIMD<int16_t, 8> a)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = (uint16x8_t)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 4> vreinterpret(SIMD<uint32_t, 4> a)
+{
+	SIMD<int32_t, 4> tmp;
+	tmp.m = (int32x4_t)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vreinterpret(SIMD<int32_t, 4> a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = (uint32x4_t)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 4> vdup(float a)
+{
+	SIMD<float, 4> tmp;
+	tmp.m = vdupq_n_f32(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vdup(int8_t a)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = vdupq_n_s8(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vdup(int16_t a)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = vdupq_n_s16(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 4> vdup(int32_t a)
+{
+	SIMD<int32_t, 4> tmp;
+	tmp.m = vdupq_n_s32(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<int64_t, 2> vdup(int64_t a)
+{
+	SIMD<int64_t, 2> tmp;
+	tmp.m = vdupq_n_s64(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vdup(uint8_t a)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = vdupq_n_u8(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vdup(uint16_t a)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = vdupq_n_u16(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vdup(uint32_t a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = vdupq_n_u32(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vdup(uint64_t a)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = vdupq_n_u64(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 4> vzero()
+{
+	SIMD<float, 4> tmp;
+	tmp.m = (float32x4_t)veorq_u32((uint32x4_t)tmp.m, (uint32x4_t)tmp.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vzero()
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = veorq_s8(tmp.m, tmp.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vzero()
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = veorq_s16(tmp.m, tmp.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 4> vzero()
+{
+	SIMD<int32_t, 4> tmp;
+	tmp.m = veorq_s32(tmp.m, tmp.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int64_t, 2> vzero()
+{
+	SIMD<int64_t, 2> tmp;
+	tmp.m = veorq_s64(tmp.m, tmp.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vzero()
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = veorq_u8(tmp.m, tmp.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vzero()
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = veorq_u16(tmp.m, tmp.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vzero()
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = veorq_u32(tmp.m, tmp.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vzero()
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = veorq_u64(tmp.m, tmp.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 4> vadd(SIMD<float, 4> a, SIMD<float, 4> b)
+{
+	SIMD<float, 4> tmp;
+	tmp.m = vaddq_f32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vadd(SIMD<int8_t, 16> a, SIMD<int8_t, 16> b)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = vaddq_s8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vadd(SIMD<int16_t, 8> a, SIMD<int16_t, 8> b)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = vaddq_s16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 4> vadd(SIMD<int32_t, 4> a, SIMD<int32_t, 4> b)
+{
+	SIMD<int32_t, 4> tmp;
+	tmp.m = vaddq_s32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int64_t, 2> vadd(SIMD<int64_t, 2> a, SIMD<int64_t, 2> b)
+{
+	SIMD<int64_t, 2> tmp;
+	tmp.m = vaddq_s64(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vqadd(SIMD<int8_t, 16> a, SIMD<int8_t, 16> b)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = vqaddq_s8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vqadd(SIMD<int16_t, 8> a, SIMD<int16_t, 8> b)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = vqaddq_s16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 4> vsub(SIMD<float, 4> a, SIMD<float, 4> b)
+{
+	SIMD<float, 4> tmp;
+	tmp.m = vsubq_f32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vsub(SIMD<int8_t, 16> a, SIMD<int8_t, 16> b)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = vsubq_s8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vsub(SIMD<int16_t, 8> a, SIMD<int16_t, 8> b)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = vsubq_s16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 4> vsub(SIMD<int32_t, 4> a, SIMD<int32_t, 4> b)
+{
+	SIMD<int32_t, 4> tmp;
+	tmp.m = vsubq_s32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int64_t, 2> vsub(SIMD<int64_t, 2> a, SIMD<int64_t, 2> b)
+{
+	SIMD<int64_t, 2> tmp;
+	tmp.m = vsubq_s64(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vqsub(SIMD<int8_t, 16> a, SIMD<int8_t, 16> b)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = vqsubq_s8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vqsub(SIMD<int16_t, 8> a, SIMD<int16_t, 8> b)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = vqsubq_s16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vqsub(SIMD<uint8_t, 16> a, SIMD<uint8_t, 16> b)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = vqsubq_u8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vqsub(SIMD<uint16_t, 8> a, SIMD<uint16_t, 8> b)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = vqsubq_u16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 4> vabs(SIMD<float, 4> a)
+{
+	SIMD<float, 4> tmp;
+	tmp.m = vabsq_f32(a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vqabs(SIMD<int8_t, 16> a)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = vqabsq_s8(a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vqabs(SIMD<int16_t, 8> a)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = vqabsq_s16(a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 4> vsign(SIMD<float, 4> a, SIMD<float, 4> b)
+{
+	SIMD<float, 4> tmp;
+	tmp.m = (float32x4_t)vbicq_u32(
+		veorq_u32((uint32x4_t)a.m, vandq_u32((uint32x4_t)vdupq_n_f32(-0.f), (uint32x4_t)b.m)),
+		vceqq_f32(b.m, vdupq_n_f32(0.f)));
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vsign(SIMD<int8_t, 16> a, SIMD<int8_t, 16> b)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = (int8x16_t)vorrq_u8(
+		vandq_u8(vcgtq_s8(vdupq_n_s8(0), b.m), (uint8x16_t)vnegq_s8(a.m)),
+		vandq_u8(vcgtq_s8(b.m, vdupq_n_s8(0)), (uint8x16_t)a.m));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vorr(SIMD<uint8_t, 16> a, SIMD<uint8_t, 16> b)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = vorrq_u8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vorr(SIMD<uint16_t, 8> a, SIMD<uint16_t, 8> b)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = vorrq_u16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vorr(SIMD<uint32_t, 4> a, SIMD<uint32_t, 4> b)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = vorrq_u32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vorr(SIMD<uint64_t, 2> a, SIMD<uint64_t, 2> b)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = vorrq_u64(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vand(SIMD<uint8_t, 16> a, SIMD<uint8_t, 16> b)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = vandq_u8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vand(SIMD<uint16_t, 8> a, SIMD<uint16_t, 8> b)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = vandq_u16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vand(SIMD<uint32_t, 4> a, SIMD<uint32_t, 4> b)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = vandq_u32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vand(SIMD<uint64_t, 2> a, SIMD<uint64_t, 2> b)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = vandq_u64(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> veor(SIMD<uint8_t, 16> a, SIMD<uint8_t, 16> b)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = veorq_u8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> veor(SIMD<uint16_t, 8> a, SIMD<uint16_t, 8> b)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = veorq_u16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> veor(SIMD<uint32_t, 4> a, SIMD<uint32_t, 4> b)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = veorq_u32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> veor(SIMD<uint64_t, 2> a, SIMD<uint64_t, 2> b)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = veorq_u64(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vbic(SIMD<uint8_t, 16> a, SIMD<uint8_t, 16> b)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = vbicq_u8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vbic(SIMD<uint16_t, 8> a, SIMD<uint16_t, 8> b)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = vbicq_u16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vbic(SIMD<uint32_t, 4> a, SIMD<uint32_t, 4> b)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = vbicq_u32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vbic(SIMD<uint64_t, 2> a, SIMD<uint64_t, 2> b)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = vbicq_u64(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vbsl(SIMD<uint8_t, 16> a, SIMD<uint8_t, 16> b, SIMD<uint8_t, 16> c)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = vbslq_u8(a.m, b.m, c.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vbsl(SIMD<uint16_t, 8> a, SIMD<uint16_t, 8> b, SIMD<uint16_t, 8> c)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = vbslq_u16(a.m, b.m, c.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vbsl(SIMD<uint32_t, 4> a, SIMD<uint32_t, 4> b, SIMD<uint32_t, 4> c)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = vbslq_u32(a.m, b.m, c.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vbsl(SIMD<uint64_t, 2> a, SIMD<uint64_t, 2> b, SIMD<uint64_t, 2> c)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = vbslq_u64(a.m, b.m, c.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vceqz(SIMD<float, 4> a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = vceqq_f32(a.m, vdupq_n_f32(0.f));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vceqz(SIMD<int8_t, 16> a)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = vceqq_s8(a.m, vdupq_n_s8(0));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vceqz(SIMD<int16_t, 8> a)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = vceqq_s16(a.m, vdupq_n_s16(0));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vceqz(SIMD<int32_t, 4> a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = vceqq_s32(a.m, vdupq_n_s32(0));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vceq(SIMD<float, 4> a, SIMD<float, 4> b)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = vceqq_f32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vceq(SIMD<int8_t, 16> a, SIMD<int8_t, 16> b)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = vceqq_s8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vceq(SIMD<int16_t, 8> a, SIMD<int16_t, 8> b)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = vceqq_s16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vceq(SIMD<int32_t, 4> a, SIMD<int32_t, 4> b)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = vceqq_s32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vcgtz(SIMD<float, 4> a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = vcgtq_f32(a.m, vdupq_n_f32(0.f));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vcgtz(SIMD<int8_t, 16> a)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = vcgtq_s8(a.m, vdupq_n_s8(0));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vcgtz(SIMD<int16_t, 8> a)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = vcgtq_s16(a.m, vdupq_n_s16(0));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vcgtz(SIMD<int32_t, 4> a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = vcgtq_s32(a.m, vdupq_n_s32(0));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vcltz(SIMD<float, 4> a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = vcltq_f32(a.m, vdupq_n_f32(0.f));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vcltz(SIMD<int8_t, 16> a)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = vcltq_s8(a.m, vdupq_n_s8(0));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vcltz(SIMD<int16_t, 8> a)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = vcltq_s16(a.m, vdupq_n_s16(0));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vcltz(SIMD<int32_t, 4> a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = vcltq_s32(a.m, vdupq_n_s32(0));
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 4> vmin(SIMD<float, 4> a, SIMD<float, 4> b)
+{
+	SIMD<float, 4> tmp;
+	tmp.m = vminq_f32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vmin(SIMD<int8_t, 16> a, SIMD<int8_t, 16> b)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = vminq_s8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vmin(SIMD<int16_t, 8> a, SIMD<int16_t, 8> b)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = vminq_s16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 4> vmin(SIMD<int32_t, 4> a, SIMD<int32_t, 4> b)
+{
+	SIMD<int32_t, 4> tmp;
+	tmp.m = vminq_s32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 4> vmax(SIMD<float, 4> a, SIMD<float, 4> b)
+{
+	SIMD<float, 4> tmp;
+	tmp.m = vmaxq_f32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vmax(SIMD<int8_t, 16> a, SIMD<int8_t, 16> b)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = vmaxq_s8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vmax(SIMD<int16_t, 8> a, SIMD<int16_t, 8> b)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = vmaxq_s16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 4> vmax(SIMD<int32_t, 4> a, SIMD<int32_t, 4> b)
+{
+	SIMD<int32_t, 4> tmp;
+	tmp.m = vmaxq_s32(a.m, b.m);
+	return tmp;
+}
+
+#endif
diff --git a/simd.hh b/simd.hh
new file mode 100644
index 0000000..197ca75
--- /dev/null
+++ b/simd.hh
@@ -0,0 +1,1189 @@
+/*
+Single instruction, multiple data
+
+Copyright 2018 Ahmet Inan <inan@aicodix.de>
+*/
+
+#ifndef SIMD_HH
+#define SIMD_HH
+
+#include <cstdint>
+#include <cstdlib>
+#include <cmath>
+
+template <typename TYPE, int WIDTH>
+union SIMD;
+
+template <int WIDTH>
+union SIMD<float, WIDTH>
+{
+	static const int SIZE = WIDTH;
+	typedef float value_type;
+	typedef uint32_t uint_type;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <int WIDTH>
+union SIMD<double, WIDTH>
+{
+	static const int SIZE = WIDTH;
+	typedef double value_type;
+	typedef uint64_t uint_type;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <int WIDTH>
+union SIMD<int8_t, WIDTH>
+{
+	static const int SIZE = WIDTH;
+	typedef int8_t value_type;
+	typedef uint8_t uint_type;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <int WIDTH>
+union SIMD<int16_t, WIDTH>
+{
+	static const int SIZE = WIDTH;
+	typedef int16_t value_type;
+	typedef uint16_t uint_type;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <int WIDTH>
+union SIMD<int32_t, WIDTH>
+{
+	static const int SIZE = WIDTH;
+	typedef int32_t value_type;
+	typedef uint32_t uint_type;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <int WIDTH>
+union SIMD<int64_t, WIDTH>
+{
+	static const int SIZE = WIDTH;
+	typedef int64_t value_type;
+	typedef uint64_t uint_type;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <int WIDTH>
+union SIMD<uint8_t, WIDTH>
+{
+	static const int SIZE = WIDTH;
+	typedef uint8_t value_type;
+	typedef uint8_t uint_type;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <int WIDTH>
+union SIMD<uint16_t, WIDTH>
+{
+	static const int SIZE = WIDTH;
+	typedef uint16_t value_type;
+	typedef uint16_t uint_type;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <int WIDTH>
+union SIMD<uint32_t, WIDTH>
+{
+	static const int SIZE = WIDTH;
+	typedef uint32_t value_type;
+	typedef uint32_t uint_type;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <int WIDTH>
+union SIMD<uint64_t, WIDTH>
+{
+	static const int SIZE = WIDTH;
+	typedef uint64_t value_type;
+	typedef uint64_t uint_type;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <typename TYPE>
+static inline TYPE vdup(typename TYPE::value_type a)
+{
+	TYPE tmp;
+	for (int i = 0; i < TYPE::SIZE; ++i)
+		tmp.v[i] = a;
+	return tmp;
+}
+
+template <typename TYPE>
+static inline TYPE vzero()
+{
+	TYPE tmp;
+	for (int i = 0; i < TYPE::SIZE; ++i)
+		tmp.u[i] ^= tmp.u[i];
+	return tmp;
+}
+
+template <typename DST, typename SRC>
+static inline DST vreinterpret(SRC a)
+{
+	static_assert(SRC::SIZE == DST::SIZE, "source and destination width must be same");
+	static_assert(sizeof(typename SRC::value_type) == sizeof(typename DST::value_type), "source and destination value type sizes must be same");
+	DST tmp;
+	for (int i = 0; i < DST::SIZE; ++i)
+		tmp.u[i] = a.u[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vmask(SIMD<float, WIDTH> a)
+{
+	return vreinterpret<SIMD<uint32_t, WIDTH>>(a);
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vmask(SIMD<double, WIDTH> a)
+{
+	return vreinterpret<SIMD<uint64_t, WIDTH>>(a);
+}
+
+template <int WIDTH>
+static inline SIMD<uint8_t, WIDTH> vmask(SIMD<int8_t, WIDTH> a)
+{
+	return vreinterpret<SIMD<uint8_t, WIDTH>>(a);
+}
+
+template <int WIDTH>
+static inline SIMD<uint16_t, WIDTH> vmask(SIMD<int16_t, WIDTH> a)
+{
+	return vreinterpret<SIMD<uint16_t, WIDTH>>(a);
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vmask(SIMD<int32_t, WIDTH> a)
+{
+	return vreinterpret<SIMD<uint32_t, WIDTH>>(a);
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vmask(SIMD<int64_t, WIDTH> a)
+{
+	return vreinterpret<SIMD<uint64_t, WIDTH>>(a);
+}
+
+template <int WIDTH>
+static inline SIMD<uint8_t, WIDTH> vunsigned(SIMD<int8_t, WIDTH> a)
+{
+	return vreinterpret<SIMD<uint8_t, WIDTH>>(a);
+}
+
+template <int WIDTH>
+static inline SIMD<uint16_t, WIDTH> vunsigned(SIMD<int16_t, WIDTH> a)
+{
+	return vreinterpret<SIMD<uint16_t, WIDTH>>(a);
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vunsigned(SIMD<int32_t, WIDTH> a)
+{
+	return vreinterpret<SIMD<uint32_t, WIDTH>>(a);
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vunsigned(SIMD<int64_t, WIDTH> a)
+{
+	return vreinterpret<SIMD<uint64_t, WIDTH>>(a);
+}
+
+template <int WIDTH>
+static inline SIMD<int8_t, WIDTH> vsigned(SIMD<uint8_t, WIDTH> a)
+{
+	return vreinterpret<SIMD<int8_t, WIDTH>>(a);
+}
+
+template <int WIDTH>
+static inline SIMD<int16_t, WIDTH> vsigned(SIMD<uint16_t, WIDTH> a)
+{
+	return vreinterpret<SIMD<int16_t, WIDTH>>(a);
+}
+
+template <int WIDTH>
+static inline SIMD<int32_t, WIDTH> vsigned(SIMD<uint32_t, WIDTH> a)
+{
+	return vreinterpret<SIMD<int32_t, WIDTH>>(a);
+}
+
+template <int WIDTH>
+static inline SIMD<int64_t, WIDTH> vsigned(SIMD<uint64_t, WIDTH> a)
+{
+	return vreinterpret<SIMD<int64_t, WIDTH>>(a);
+}
+
+template <int WIDTH>
+static inline SIMD<float, WIDTH> vneg(SIMD<float, WIDTH> a)
+{
+	SIMD<float, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<double, WIDTH> vneg(SIMD<double, WIDTH> a)
+{
+	SIMD<double, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int8_t, WIDTH> vneg(SIMD<int8_t, WIDTH> a)
+{
+	SIMD<int8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int16_t, WIDTH> vneg(SIMD<int16_t, WIDTH> a)
+{
+	SIMD<int16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int32_t, WIDTH> vneg(SIMD<int32_t, WIDTH> a)
+{
+	SIMD<int32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int64_t, WIDTH> vneg(SIMD<int64_t, WIDTH> a)
+{
+	SIMD<int64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -a.v[i];
+	return tmp;
+}
+
+template <typename TYPE, int WIDTH>
+static inline SIMD<TYPE, WIDTH> vabs(SIMD<TYPE, WIDTH> a)
+{
+	SIMD<TYPE, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::abs(a.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int8_t, WIDTH> vqabs(SIMD<int8_t, WIDTH> a)
+{
+	SIMD<int8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::abs(std::max<int8_t>(a.v[i], -INT8_MAX));
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int16_t, WIDTH> vqabs(SIMD<int16_t, WIDTH> a)
+{
+	SIMD<int16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::abs(std::max<int16_t>(a.v[i], -INT16_MAX));
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int32_t, WIDTH> vqabs(SIMD<int32_t, WIDTH> a)
+{
+	SIMD<int32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::abs(std::max<int32_t>(a.v[i], -INT32_MAX));
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int64_t, WIDTH> vqabs(SIMD<int64_t, WIDTH> a)
+{
+	SIMD<int64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::abs(std::max<int64_t>(a.v[i], -INT64_MAX));
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint8_t, WIDTH> vnot(SIMD<uint8_t, WIDTH> a)
+{
+	SIMD<uint8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = ~a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint16_t, WIDTH> vnot(SIMD<uint16_t, WIDTH> a)
+{
+	SIMD<uint16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = ~a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vnot(SIMD<uint32_t, WIDTH> a)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = ~a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vnot(SIMD<uint64_t, WIDTH> a)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = ~a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint8_t, WIDTH> vorr(SIMD<uint8_t, WIDTH> a, SIMD<uint8_t, WIDTH> b)
+{
+	SIMD<uint8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] | b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint16_t, WIDTH> vorr(SIMD<uint16_t, WIDTH> a, SIMD<uint16_t, WIDTH> b)
+{
+	SIMD<uint16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] | b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vorr(SIMD<uint32_t, WIDTH> a, SIMD<uint32_t, WIDTH> b)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] | b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vorr(SIMD<uint64_t, WIDTH> a, SIMD<uint64_t, WIDTH> b)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] | b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint8_t, WIDTH> vand(SIMD<uint8_t, WIDTH> a, SIMD<uint8_t, WIDTH> b)
+{
+	SIMD<uint8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] & b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint16_t, WIDTH> vand(SIMD<uint16_t, WIDTH> a, SIMD<uint16_t, WIDTH> b)
+{
+	SIMD<uint16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] & b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vand(SIMD<uint32_t, WIDTH> a, SIMD<uint32_t, WIDTH> b)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] & b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vand(SIMD<uint64_t, WIDTH> a, SIMD<uint64_t, WIDTH> b)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] & b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint8_t, WIDTH> veor(SIMD<uint8_t, WIDTH> a, SIMD<uint8_t, WIDTH> b)
+{
+	SIMD<uint8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] ^ b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint16_t, WIDTH> veor(SIMD<uint16_t, WIDTH> a, SIMD<uint16_t, WIDTH> b)
+{
+	SIMD<uint16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] ^ b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> veor(SIMD<uint32_t, WIDTH> a, SIMD<uint32_t, WIDTH> b)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] ^ b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> veor(SIMD<uint64_t, WIDTH> a, SIMD<uint64_t, WIDTH> b)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] ^ b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint8_t, WIDTH> vbic(SIMD<uint8_t, WIDTH> a, SIMD<uint8_t, WIDTH> b)
+{
+	SIMD<uint8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] & ~b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint16_t, WIDTH> vbic(SIMD<uint16_t, WIDTH> a, SIMD<uint16_t, WIDTH> b)
+{
+	SIMD<uint16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] & ~b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vbic(SIMD<uint32_t, WIDTH> a, SIMD<uint32_t, WIDTH> b)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] & ~b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vbic(SIMD<uint64_t, WIDTH> a, SIMD<uint64_t, WIDTH> b)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] & ~b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint8_t, WIDTH> vbsl(SIMD<uint8_t, WIDTH> a, SIMD<uint8_t, WIDTH> b, SIMD<uint8_t, WIDTH> c)
+{
+	SIMD<uint8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = (a.v[i] & b.v[i]) | (~a.v[i] & c.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint16_t, WIDTH> vbsl(SIMD<uint16_t, WIDTH> a, SIMD<uint16_t, WIDTH> b, SIMD<uint16_t, WIDTH> c)
+{
+	SIMD<uint16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = (a.v[i] & b.v[i]) | (~a.v[i] & c.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vbsl(SIMD<uint32_t, WIDTH> a, SIMD<uint32_t, WIDTH> b, SIMD<uint32_t, WIDTH> c)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = (a.v[i] & b.v[i]) | (~a.v[i] & c.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vbsl(SIMD<uint64_t, WIDTH> a, SIMD<uint64_t, WIDTH> b, SIMD<uint64_t, WIDTH> c)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = (a.v[i] & b.v[i]) | (~a.v[i] & c.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vcgtz(SIMD<float, WIDTH> a)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] > 0.f);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vcgtz(SIMD<double, WIDTH> a)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] > 0.);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint8_t, WIDTH> vcgtz(SIMD<int8_t, WIDTH> a)
+{
+	SIMD<uint8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] > 0);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint16_t, WIDTH> vcgtz(SIMD<int16_t, WIDTH> a)
+{
+	SIMD<uint16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] > 0);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vcgtz(SIMD<int32_t, WIDTH> a)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] > 0);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vcgtz(SIMD<int64_t, WIDTH> a)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] > 0);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vceqz(SIMD<float, WIDTH> a)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] == 0.f);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vceqz(SIMD<double, WIDTH> a)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] == 0.);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint8_t, WIDTH> vceqz(SIMD<int8_t, WIDTH> a)
+{
+	SIMD<uint8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -!a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint16_t, WIDTH> vceqz(SIMD<int16_t, WIDTH> a)
+{
+	SIMD<uint16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -!a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vceqz(SIMD<int32_t, WIDTH> a)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -!a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vceqz(SIMD<int64_t, WIDTH> a)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -!a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vcltz(SIMD<float, WIDTH> a)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] < 0.f);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vcltz(SIMD<double, WIDTH> a)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] < 0.);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint8_t, WIDTH> vcltz(SIMD<int8_t, WIDTH> a)
+{
+	SIMD<uint8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] < 0);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint16_t, WIDTH> vcltz(SIMD<int16_t, WIDTH> a)
+{
+	SIMD<uint16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] < 0);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vcltz(SIMD<int32_t, WIDTH> a)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] < 0);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vcltz(SIMD<int64_t, WIDTH> a)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] < 0);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vcgt(SIMD<float, WIDTH> a, SIMD<float, WIDTH> b)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] > b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vcgt(SIMD<double, WIDTH> a, SIMD<double, WIDTH> b)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] > b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint8_t, WIDTH> vcgt(SIMD<int8_t, WIDTH> a, SIMD<int8_t, WIDTH> b)
+{
+	SIMD<uint8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] > b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint16_t, WIDTH> vcgt(SIMD<int16_t, WIDTH> a, SIMD<int16_t, WIDTH> b)
+{
+	SIMD<uint16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] > b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vcgt(SIMD<int32_t, WIDTH> a, SIMD<int32_t, WIDTH> b)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] > b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vcgt(SIMD<int64_t, WIDTH> a, SIMD<int64_t, WIDTH> b)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] > b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint8_t, WIDTH> vcgt(SIMD<uint8_t, WIDTH> a, SIMD<uint8_t, WIDTH> b)
+{
+	SIMD<uint8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] > b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint16_t, WIDTH> vcgt(SIMD<uint16_t, WIDTH> a, SIMD<uint16_t, WIDTH> b)
+{
+	SIMD<uint16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] > b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vcgt(SIMD<uint32_t, WIDTH> a, SIMD<uint32_t, WIDTH> b)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] > b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vcgt(SIMD<uint64_t, WIDTH> a, SIMD<uint64_t, WIDTH> b)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] > b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vceq(SIMD<float, WIDTH> a, SIMD<float, WIDTH> b)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] == b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vceq(SIMD<double, WIDTH> a, SIMD<double, WIDTH> b)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] == b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint8_t, WIDTH> vceq(SIMD<int8_t, WIDTH> a, SIMD<int8_t, WIDTH> b)
+{
+	SIMD<uint8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] == b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint16_t, WIDTH> vceq(SIMD<int16_t, WIDTH> a, SIMD<int16_t, WIDTH> b)
+{
+	SIMD<uint16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] == b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint32_t, WIDTH> vceq(SIMD<int32_t, WIDTH> a, SIMD<int32_t, WIDTH> b)
+{
+	SIMD<uint32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] == b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint64_t, WIDTH> vceq(SIMD<int64_t, WIDTH> a, SIMD<int64_t, WIDTH> b)
+{
+	SIMD<uint64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = -(a.v[i] == b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<float, WIDTH> vmin(SIMD<float, WIDTH> a, SIMD<float, WIDTH> b)
+{
+	SIMD<float, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::min(a.v[i], b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<double, WIDTH> vmin(SIMD<double, WIDTH> a, SIMD<double, WIDTH> b)
+{
+	SIMD<double, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::min(a.v[i], b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int8_t, WIDTH> vmin(SIMD<int8_t, WIDTH> a, SIMD<int8_t, WIDTH> b)
+{
+	SIMD<int8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::min(a.v[i], b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int16_t, WIDTH> vmin(SIMD<int16_t, WIDTH> a, SIMD<int16_t, WIDTH> b)
+{
+	SIMD<int16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::min(a.v[i], b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int32_t, WIDTH> vmin(SIMD<int32_t, WIDTH> a, SIMD<int32_t, WIDTH> b)
+{
+	SIMD<int32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::min(a.v[i], b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int64_t, WIDTH> vmin(SIMD<int64_t, WIDTH> a, SIMD<int64_t, WIDTH> b)
+{
+	SIMD<int64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::min(a.v[i], b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<float, WIDTH> vmax(SIMD<float, WIDTH> a, SIMD<float, WIDTH> b)
+{
+	SIMD<float, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::max(a.v[i], b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<double, WIDTH> vmax(SIMD<double, WIDTH> a, SIMD<double, WIDTH> b)
+{
+	SIMD<double, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::max(a.v[i], b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int8_t, WIDTH> vmax(SIMD<int8_t, WIDTH> a, SIMD<int8_t, WIDTH> b)
+{
+	SIMD<int8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::max(a.v[i], b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int16_t, WIDTH> vmax(SIMD<int16_t, WIDTH> a, SIMD<int16_t, WIDTH> b)
+{
+	SIMD<int16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::max(a.v[i], b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int32_t, WIDTH> vmax(SIMD<int32_t, WIDTH> a, SIMD<int32_t, WIDTH> b)
+{
+	SIMD<int32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::max(a.v[i], b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int64_t, WIDTH> vmax(SIMD<int64_t, WIDTH> a, SIMD<int64_t, WIDTH> b)
+{
+	SIMD<int64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::max(a.v[i], b.v[i]);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<float, WIDTH> vadd(SIMD<float, WIDTH> a, SIMD<float, WIDTH> b)
+{
+	SIMD<float, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] + b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<double, WIDTH> vadd(SIMD<double, WIDTH> a, SIMD<double, WIDTH> b)
+{
+	SIMD<double, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] + b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int8_t, WIDTH> vadd(SIMD<int8_t, WIDTH> a, SIMD<int8_t, WIDTH> b)
+{
+	SIMD<int8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] + b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int16_t, WIDTH> vadd(SIMD<int16_t, WIDTH> a, SIMD<int16_t, WIDTH> b)
+{
+	SIMD<int16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] + b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int32_t, WIDTH> vadd(SIMD<int32_t, WIDTH> a, SIMD<int32_t, WIDTH> b)
+{
+	SIMD<int32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] + b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int64_t, WIDTH> vadd(SIMD<int64_t, WIDTH> a, SIMD<int64_t, WIDTH> b)
+{
+	SIMD<int64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] + b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int8_t, WIDTH> vqadd(SIMD<int8_t, WIDTH> a, SIMD<int8_t, WIDTH> b)
+{
+	SIMD<int8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::min<int16_t>(std::max<int16_t>(int16_t(a.v[i]) + int16_t(b.v[i]), INT8_MIN), INT8_MAX);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int16_t, WIDTH> vqadd(SIMD<int16_t, WIDTH> a, SIMD<int16_t, WIDTH> b)
+{
+	SIMD<int16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::min<int32_t>(std::max<int32_t>(int32_t(a.v[i]) + int32_t(b.v[i]), INT16_MIN), INT16_MAX);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<float, WIDTH> vsub(SIMD<float, WIDTH> a, SIMD<float, WIDTH> b)
+{
+	SIMD<float, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] - b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<double, WIDTH> vsub(SIMD<double, WIDTH> a, SIMD<double, WIDTH> b)
+{
+	SIMD<double, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] - b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int8_t, WIDTH> vsub(SIMD<int8_t, WIDTH> a, SIMD<int8_t, WIDTH> b)
+{
+	SIMD<int8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] - b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int16_t, WIDTH> vsub(SIMD<int16_t, WIDTH> a, SIMD<int16_t, WIDTH> b)
+{
+	SIMD<int16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] - b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int32_t, WIDTH> vsub(SIMD<int32_t, WIDTH> a, SIMD<int32_t, WIDTH> b)
+{
+	SIMD<int32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] - b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int64_t, WIDTH> vsub(SIMD<int64_t, WIDTH> a, SIMD<int64_t, WIDTH> b)
+{
+	SIMD<int64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = a.v[i] - b.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int8_t, WIDTH> vqsub(SIMD<int8_t, WIDTH> a, SIMD<int8_t, WIDTH> b)
+{
+	SIMD<int8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::min<int16_t>(std::max<int16_t>(int16_t(a.v[i]) - int16_t(b.v[i]), INT8_MIN), INT8_MAX);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int16_t, WIDTH> vqsub(SIMD<int16_t, WIDTH> a, SIMD<int16_t, WIDTH> b)
+{
+	SIMD<int16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::min<int32_t>(std::max<int32_t>(int32_t(a.v[i]) - int32_t(b.v[i]), INT16_MIN), INT16_MAX);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint8_t, WIDTH> vqsub(SIMD<uint8_t, WIDTH> a, SIMD<uint8_t, WIDTH> b)
+{
+	SIMD<uint8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::max<int16_t>(int16_t(a.v[i]) - int16_t(b.v[i]), 0);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<uint16_t, WIDTH> vqsub(SIMD<uint16_t, WIDTH> a, SIMD<uint16_t, WIDTH> b)
+{
+	SIMD<uint16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = std::max<int32_t>(int32_t(a.v[i]) - int32_t(b.v[i]), 0);
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<float, WIDTH> vsign(SIMD<float, WIDTH> a, SIMD<float, WIDTH> b)
+{
+	SIMD<float, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = ((b.v[i] > 0.f) - (b.v[i] < 0.f)) * a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<double, WIDTH> vsign(SIMD<double, WIDTH> a, SIMD<double, WIDTH> b)
+{
+	SIMD<double, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = ((b.v[i] > 0.) - (b.v[i] < 0.)) * a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int8_t, WIDTH> vsign(SIMD<int8_t, WIDTH> a, SIMD<int8_t, WIDTH> b)
+{
+	SIMD<int8_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = ((b.v[i] > 0) - (b.v[i] < 0)) * a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int16_t, WIDTH> vsign(SIMD<int16_t, WIDTH> a, SIMD<int16_t, WIDTH> b)
+{
+	SIMD<int16_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = ((b.v[i] > 0) - (b.v[i] < 0)) * a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int32_t, WIDTH> vsign(SIMD<int32_t, WIDTH> a, SIMD<int32_t, WIDTH> b)
+{
+	SIMD<int32_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = ((b.v[i] > 0) - (b.v[i] < 0)) * a.v[i];
+	return tmp;
+}
+
+template <int WIDTH>
+static inline SIMD<int64_t, WIDTH> vsign(SIMD<int64_t, WIDTH> a, SIMD<int64_t, WIDTH> b)
+{
+	SIMD<int64_t, WIDTH> tmp;
+	for (int i = 0; i < WIDTH; ++i)
+		tmp.v[i] = ((b.v[i] > 0) - (b.v[i] < 0)) * a.v[i];
+	return tmp;
+}
+
+#if 1
+#ifdef __AVX2__
+#include "avx2.hh"
+#else
+#ifdef __SSE4_1__
+#include "sse4_1.hh"
+#endif
+#endif
+
+#ifdef __ARM_NEON__
+#include "neon.hh"
+#endif
+#endif
+
+#endif
diff --git a/sse4_1.hh b/sse4_1.hh
new file mode 100644
index 0000000..42ec37d
--- /dev/null
+++ b/sse4_1.hh
@@ -0,0 +1,975 @@
+/*
+Intel SSE4.1 acceleration
+
+Copyright 2018 Ahmet Inan <inan@aicodix.de>
+*/
+
+#ifndef SSE4_1_HH
+#define SSE4_1_HH
+
+#include <smmintrin.h>
+
+template <>
+union SIMD<float, 4>
+{
+	static const int SIZE = 8;
+	typedef float value_type;
+	typedef uint32_t uint_type;
+	__m128 m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<double, 2>
+{
+	static const int SIZE = 4;
+	typedef double value_type;
+	typedef uint64_t uint_type;
+	__m128d m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<int8_t, 16>
+{
+	static const int SIZE = 32;
+	typedef int8_t value_type;
+	typedef uint8_t uint_type;
+	__m128i m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<int16_t, 8>
+{
+	static const int SIZE = 16;
+	typedef int16_t value_type;
+	typedef uint16_t uint_type;
+	__m128i m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<int32_t, 4>
+{
+	static const int SIZE = 8;
+	typedef int32_t value_type;
+	typedef uint32_t uint_type;
+	__m128i m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<int64_t, 2>
+{
+	static const int SIZE = 4;
+	typedef int64_t value_type;
+	typedef uint64_t uint_type;
+	__m128i m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<uint8_t, 16>
+{
+	static const int SIZE = 32;
+	typedef uint8_t value_type;
+	typedef uint8_t uint_type;
+	__m128i m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<uint16_t, 8>
+{
+	static const int SIZE = 16;
+	typedef uint16_t value_type;
+	typedef uint16_t uint_type;
+	__m128i m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<uint32_t, 4>
+{
+	static const int SIZE = 8;
+	typedef uint32_t value_type;
+	typedef uint32_t uint_type;
+	__m128i m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+union SIMD<uint64_t, 2>
+{
+	static const int SIZE = 4;
+	typedef uint64_t value_type;
+	typedef uint64_t uint_type;
+	__m128i m;
+	value_type v[SIZE];
+	uint_type u[SIZE];
+};
+
+template <>
+inline SIMD<float, 4> vreinterpret(SIMD<uint32_t, 4> a)
+{
+	SIMD<float, 4> tmp;
+	tmp.m = (__m128)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vreinterpret(SIMD<float, 4> a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = (__m128i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 2> vreinterpret(SIMD<uint64_t, 2> a)
+{
+	SIMD<double, 2> tmp;
+	tmp.m = (__m128d)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vreinterpret(SIMD<double, 2> a)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = (__m128i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vreinterpret(SIMD<int8_t, 16> a)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = (__m128i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vreinterpret(SIMD<uint8_t, 16> a)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = (__m128i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vreinterpret(SIMD<int16_t, 8> a)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = (__m128i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vreinterpret(SIMD<uint16_t, 8> a)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = (__m128i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vreinterpret(SIMD<int32_t, 4> a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = (__m128i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 4> vreinterpret(SIMD<uint32_t, 4> a)
+{
+	SIMD<int32_t, 4> tmp;
+	tmp.m = (__m128i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vreinterpret(SIMD<int64_t, 2> a)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = (__m128i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<int64_t, 2> vreinterpret(SIMD<uint64_t, 2> a)
+{
+	SIMD<int64_t, 2> tmp;
+	tmp.m = (__m128i)a.m;
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 4> vdup<SIMD<float, 4>>(float a)
+{
+	SIMD<float, 4> tmp;
+	tmp.m = _mm_set1_ps(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 2> vdup<SIMD<double, 2>>(double a)
+{
+	SIMD<double, 2> tmp;
+	tmp.m = _mm_set1_pd(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vdup<SIMD<int8_t, 16>>(int8_t a)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = _mm_set1_epi8(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vdup<SIMD<int16_t, 8>>(int16_t a)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = _mm_set1_epi16(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 4> vdup<SIMD<int32_t, 4>>(int32_t a)
+{
+	SIMD<int32_t, 4> tmp;
+	tmp.m = _mm_set1_epi32(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<int64_t, 2> vdup<SIMD<int64_t, 2>>(int64_t a)
+{
+	SIMD<int64_t, 2> tmp;
+	tmp.m = _mm_set1_epi64x(a);
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 4> vzero()
+{
+	SIMD<float, 4> tmp;
+	tmp.m = _mm_setzero_ps();
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 2> vzero()
+{
+	SIMD<double, 2> tmp;
+	tmp.m = _mm_setzero_pd();
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vzero()
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = _mm_setzero_si128();
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vzero()
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = _mm_setzero_si128();
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 4> vzero()
+{
+	SIMD<int32_t, 4> tmp;
+	tmp.m = _mm_setzero_si128();
+	return tmp;
+}
+
+template <>
+inline SIMD<int64_t, 2> vzero()
+{
+	SIMD<int64_t, 2> tmp;
+	tmp.m = _mm_setzero_si128();
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 4> vadd(SIMD<float, 4> a, SIMD<float, 4> b)
+{
+	SIMD<float, 4> tmp;
+	tmp.m = _mm_add_ps(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 2> vadd(SIMD<double, 2> a, SIMD<double, 2> b)
+{
+	SIMD<double, 2> tmp;
+	tmp.m = _mm_add_pd(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vadd(SIMD<int8_t, 16> a, SIMD<int8_t, 16> b)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = _mm_add_epi8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vadd(SIMD<int16_t, 8> a, SIMD<int16_t, 8> b)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = _mm_add_epi16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 4> vadd(SIMD<int32_t, 4> a, SIMD<int32_t, 4> b)
+{
+	SIMD<int32_t, 4> tmp;
+	tmp.m = _mm_add_epi32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int64_t, 2> vadd(SIMD<int64_t, 2> a, SIMD<int64_t, 2> b)
+{
+	SIMD<int64_t, 2> tmp;
+	tmp.m = _mm_add_epi64(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vqadd(SIMD<int8_t, 16> a, SIMD<int8_t, 16> b)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = _mm_adds_epi8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vqadd(SIMD<int16_t, 8> a, SIMD<int16_t, 8> b)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = _mm_adds_epi16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 4> vsub(SIMD<float, 4> a, SIMD<float, 4> b)
+{
+	SIMD<float, 4> tmp;
+	tmp.m = _mm_sub_ps(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 2> vsub(SIMD<double, 2> a, SIMD<double, 2> b)
+{
+	SIMD<double, 2> tmp;
+	tmp.m = _mm_sub_pd(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vsub(SIMD<int8_t, 16> a, SIMD<int8_t, 16> b)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = _mm_sub_epi8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vsub(SIMD<int16_t, 8> a, SIMD<int16_t, 8> b)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = _mm_sub_epi16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 4> vsub(SIMD<int32_t, 4> a, SIMD<int32_t, 4> b)
+{
+	SIMD<int32_t, 4> tmp;
+	tmp.m = _mm_sub_epi32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int64_t, 2> vsub(SIMD<int64_t, 2> a, SIMD<int64_t, 2> b)
+{
+	SIMD<int64_t, 2> tmp;
+	tmp.m = _mm_sub_epi64(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vqsub(SIMD<int8_t, 16> a, SIMD<int8_t, 16> b)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = _mm_subs_epi8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vqsub(SIMD<int16_t, 8> a, SIMD<int16_t, 8> b)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = _mm_subs_epi16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vqsub(SIMD<uint8_t, 16> a, SIMD<uint8_t, 16> b)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = _mm_subs_epu8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vqsub(SIMD<uint16_t, 8> a, SIMD<uint16_t, 8> b)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = _mm_subs_epu16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 4> vabs(SIMD<float, 4> a)
+{
+	SIMD<float, 4> tmp;
+	tmp.m = _mm_andnot_ps(_mm_set1_ps(-0.f), a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 2> vabs(SIMD<double, 2> a)
+{
+	SIMD<double, 2> tmp;
+	tmp.m = _mm_andnot_pd(_mm_set1_pd(-0.), a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vqabs(SIMD<int8_t, 16> a)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = _mm_abs_epi8(_mm_max_epi8(a.m, _mm_set1_epi8(-INT8_MAX)));
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vqabs(SIMD<int16_t, 8> a)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = _mm_abs_epi16(_mm_max_epi16(a.m, _mm_set1_epi16(-INT16_MAX)));
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 4> vqabs(SIMD<int32_t, 4> a)
+{
+	SIMD<int32_t, 4> tmp;
+	tmp.m = _mm_abs_epi32(_mm_max_epi32(a.m, _mm_set1_epi32(-INT32_MAX)));
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 4> vsign(SIMD<float, 4> a, SIMD<float, 4> b)
+{
+	SIMD<float, 4> tmp;
+	tmp.m = _mm_andnot_ps(
+		_mm_cmpeq_ps(b.m, _mm_setzero_ps()),
+		_mm_xor_ps(a.m, _mm_and_ps(_mm_set1_ps(-0.f), b.m)));
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 2> vsign(SIMD<double, 2> a, SIMD<double, 2> b)
+{
+	SIMD<double, 2> tmp;
+	tmp.m = _mm_andnot_pd(
+		_mm_cmpeq_pd(b.m, _mm_setzero_pd()),
+		_mm_xor_pd(a.m, _mm_and_pd(_mm_set1_pd(-0.), b.m)));
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vsign(SIMD<int8_t, 16> a, SIMD<int8_t, 16> b)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = _mm_sign_epi8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vsign(SIMD<int16_t, 8> a, SIMD<int16_t, 8> b)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = _mm_sign_epi16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 4> vsign(SIMD<int32_t, 4> a, SIMD<int32_t, 4> b)
+{
+	SIMD<int32_t, 4> tmp;
+	tmp.m = _mm_sign_epi32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vorr(SIMD<uint8_t, 16> a, SIMD<uint8_t, 16> b)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = _mm_or_si128(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vorr(SIMD<uint16_t, 8> a, SIMD<uint16_t, 8> b)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = _mm_or_si128(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vorr(SIMD<uint32_t, 4> a, SIMD<uint32_t, 4> b)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = _mm_or_si128(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vorr(SIMD<uint64_t, 2> a, SIMD<uint64_t, 2> b)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = _mm_or_si128(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vand(SIMD<uint8_t, 16> a, SIMD<uint8_t, 16> b)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = _mm_and_si128(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vand(SIMD<uint16_t, 8> a, SIMD<uint16_t, 8> b)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = _mm_and_si128(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vand(SIMD<uint32_t, 4> a, SIMD<uint32_t, 4> b)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = _mm_and_si128(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vand(SIMD<uint64_t, 2> a, SIMD<uint64_t, 2> b)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = _mm_and_si128(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> veor(SIMD<uint8_t, 16> a, SIMD<uint8_t, 16> b)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = _mm_xor_si128(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> veor(SIMD<uint16_t, 8> a, SIMD<uint16_t, 8> b)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = _mm_xor_si128(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> veor(SIMD<uint32_t, 4> a, SIMD<uint32_t, 4> b)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = _mm_xor_si128(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> veor(SIMD<uint64_t, 2> a, SIMD<uint64_t, 2> b)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = _mm_xor_si128(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vbic(SIMD<uint8_t, 16> a, SIMD<uint8_t, 16> b)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = _mm_andnot_si128(b.m, a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vbic(SIMD<uint16_t, 8> a, SIMD<uint16_t, 8> b)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = _mm_andnot_si128(b.m, a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vbic(SIMD<uint32_t, 4> a, SIMD<uint32_t, 4> b)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = _mm_andnot_si128(b.m, a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vbic(SIMD<uint64_t, 2> a, SIMD<uint64_t, 2> b)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = _mm_andnot_si128(b.m, a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vbsl(SIMD<uint8_t, 16> a, SIMD<uint8_t, 16> b, SIMD<uint8_t, 16> c)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = _mm_or_si128(_mm_and_si128(a.m, b.m), _mm_andnot_si128(a.m, c.m));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vbsl(SIMD<uint16_t, 8> a, SIMD<uint16_t, 8> b, SIMD<uint16_t, 8> c)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = _mm_or_si128(_mm_and_si128(a.m, b.m), _mm_andnot_si128(a.m, c.m));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vbsl(SIMD<uint32_t, 4> a, SIMD<uint32_t, 4> b, SIMD<uint32_t, 4> c)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = _mm_or_si128(_mm_and_si128(a.m, b.m), _mm_andnot_si128(a.m, c.m));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vbsl(SIMD<uint64_t, 2> a, SIMD<uint64_t, 2> b, SIMD<uint64_t, 2> c)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = _mm_or_si128(_mm_and_si128(a.m, b.m), _mm_andnot_si128(a.m, c.m));
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vceqz(SIMD<float, 4> a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = (__m128i)_mm_cmpeq_ps(a.m, _mm_setzero_ps());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vceqz(SIMD<double, 2> a)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = (__m128i)_mm_cmpeq_pd(a.m, _mm_setzero_pd());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vceqz(SIMD<int8_t, 16> a)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = _mm_cmpeq_epi8(a.m, _mm_setzero_si128());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vceqz(SIMD<int16_t, 8> a)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = _mm_cmpeq_epi16(a.m, _mm_setzero_si128());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vceqz(SIMD<int32_t, 4> a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = _mm_cmpeq_epi32(a.m, _mm_setzero_si128());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vceqz(SIMD<int64_t, 2> a)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = _mm_cmpeq_epi64(a.m, _mm_setzero_si128());
+	return tmp;
+}
+
+
+template <>
+inline SIMD<uint32_t, 4> vceq(SIMD<float, 4> a, SIMD<float, 4> b)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = (__m128i)_mm_cmpeq_ps(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vceq(SIMD<double, 2> a, SIMD<double, 2> b)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = (__m128i)_mm_cmpeq_pd(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vceq(SIMD<int8_t, 16> a, SIMD<int8_t, 16> b)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = _mm_cmpeq_epi8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vceq(SIMD<int16_t, 8> a, SIMD<int16_t, 8> b)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = _mm_cmpeq_epi16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vceq(SIMD<int32_t, 4> a, SIMD<int32_t, 4> b)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = _mm_cmpeq_epi32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vceq(SIMD<int64_t, 2> a, SIMD<int64_t, 2> b)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = _mm_cmpeq_epi64(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vcgtz(SIMD<float, 4> a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = (__m128i)_mm_cmpgt_ps(a.m, _mm_setzero_ps());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vcgtz(SIMD<double, 2> a)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = (__m128i)_mm_cmpgt_pd(a.m, _mm_setzero_pd());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vcgtz(SIMD<int8_t, 16> a)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = _mm_cmpgt_epi8(a.m, _mm_setzero_si128());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vcgtz(SIMD<int16_t, 8> a)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = _mm_cmpgt_epi16(a.m, _mm_setzero_si128());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vcgtz(SIMD<int32_t, 4> a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = _mm_cmpgt_epi32(a.m, _mm_setzero_si128());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vcgtz(SIMD<int64_t, 2> a)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = _mm_cmpgt_epi64(a.m, _mm_setzero_si128());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vcltz(SIMD<float, 4> a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = (__m128i)_mm_cmplt_ps(a.m, _mm_setzero_ps());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vcltz(SIMD<double, 2> a)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = (__m128i)_mm_cmplt_pd(a.m, _mm_setzero_pd());
+	return tmp;
+}
+
+template <>
+inline SIMD<uint8_t, 16> vcltz(SIMD<int8_t, 16> a)
+{
+	SIMD<uint8_t, 16> tmp;
+	tmp.m = _mm_cmpgt_epi8(_mm_setzero_si128(), a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint16_t, 8> vcltz(SIMD<int16_t, 8> a)
+{
+	SIMD<uint16_t, 8> tmp;
+	tmp.m = _mm_cmpgt_epi16(_mm_setzero_si128(), a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint32_t, 4> vcltz(SIMD<int32_t, 4> a)
+{
+	SIMD<uint32_t, 4> tmp;
+	tmp.m = _mm_cmpgt_epi32(_mm_setzero_si128(), a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<uint64_t, 2> vcltz(SIMD<int64_t, 2> a)
+{
+	SIMD<uint64_t, 2> tmp;
+	tmp.m = _mm_cmpgt_epi64(_mm_setzero_si128(), a.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 4> vmin(SIMD<float, 4> a, SIMD<float, 4> b)
+{
+	SIMD<float, 4> tmp;
+	tmp.m = _mm_min_ps(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 2> vmin(SIMD<double, 2> a, SIMD<double, 2> b)
+{
+	SIMD<double, 2> tmp;
+	tmp.m = _mm_min_pd(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vmin(SIMD<int8_t, 16> a, SIMD<int8_t, 16> b)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = _mm_min_epi8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vmin(SIMD<int16_t, 8> a, SIMD<int16_t, 8> b)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = _mm_min_epi16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 4> vmin(SIMD<int32_t, 4> a, SIMD<int32_t, 4> b)
+{
+	SIMD<int32_t, 4> tmp;
+	tmp.m = _mm_min_epi32(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<float, 4> vmax(SIMD<float, 4> a, SIMD<float, 4> b)
+{
+	SIMD<float, 4> tmp;
+	tmp.m = _mm_max_ps(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<double, 2> vmax(SIMD<double, 2> a, SIMD<double, 2> b)
+{
+	SIMD<double, 2> tmp;
+	tmp.m = _mm_max_pd(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int8_t, 16> vmax(SIMD<int8_t, 16> a, SIMD<int8_t, 16> b)
+{
+	SIMD<int8_t, 16> tmp;
+	tmp.m = _mm_max_epi8(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int16_t, 8> vmax(SIMD<int16_t, 8> a, SIMD<int16_t, 8> b)
+{
+	SIMD<int16_t, 8> tmp;
+	tmp.m = _mm_max_epi16(a.m, b.m);
+	return tmp;
+}
+
+template <>
+inline SIMD<int32_t, 4> vmax(SIMD<int32_t, 4> a, SIMD<int32_t, 4> b)
+{
+	SIMD<int32_t, 4> tmp;
+	tmp.m = _mm_max_epi32(a.m, b.m);
+	return tmp;
+}
+
+#endif