diff --git a/include/vectorwrapper/vectorwrapper.hpp b/include/vectorwrapper/vectorwrapper.hpp
index 4892f96..6e157a8 100644
--- a/include/vectorwrapper/vectorwrapper.hpp
+++ b/include/vectorwrapper/vectorwrapper.hpp
@@ -291,6 +291,11 @@ namespace vwr {
 			scalar_type& x ( void );
 			const scalar_type& x ( void ) const;
 		};
+
+		template <bool LastVal, typename V1, typename V2, typename ComposeOp, typename Op>
+		constexpr bool compare ( const Vec<V1>& parLeft, const Vec<V2>& parRight, Op parComposeOp, Op parOp, bt::number_seq<size_type> );
+		template <bool LastVal, typename V1, typename V2, typename ComposeOp, typename Op, size_type I1, size_type... I>
+		bool compare ( const Vec<V1>& parLeft, const Vec<V2>& parRight, ComposeOp parComposeOp, Op parOp, bt::number_seq<size_type, I1, I...> );
 	} //namespace implem
 
 	template <typename V, size_type S>
@@ -392,7 +397,15 @@ namespace vwr {
 	template <typename V1, typename V2>
 	bool operator== ( const Vec<V1>& parLeft, const Vec<V2>& parRight );
 	template <typename V1, typename V2>
+	bool operator!= ( const Vec<V1>& parLeft, const Vec<V2>& parRight );
+	template <typename V1, typename V2>
 	bool operator< ( const Vec<V1>& parLeft, const Vec<V2>& parRight );
+	template <typename V1, typename V2>
+	bool operator> ( const Vec<V1>& parLeft, const Vec<V2>& parRight );
+	template <typename V1, typename V2>
+	bool operator<= ( const Vec<V1>& parLeft, const Vec<V2>& parRight );
+	template <typename V1, typename V2>
+	bool operator>= ( const Vec<V1>& parLeft, const Vec<V2>& parRight );
 
 	template <typename V>
 	bool operator== ( const Vec<V>& parLeft, const typename VectorWrapperInfo<V>::scalar_type& parRight );
diff --git a/include/vectorwrapper/vectorwrapper.inl b/include/vectorwrapper/vectorwrapper.inl
index e68883f..b3f0f66 100644
--- a/include/vectorwrapper/vectorwrapper.inl
+++ b/include/vectorwrapper/vectorwrapper.inl
@@ -279,6 +279,18 @@ namespace vwr {
 		{
 			static_assert(sizeof...(I) == S, "Bug?");
 		}
+
+		template <bool LastVal, typename V1, typename V2, typename ComposeOp, typename Op>
+		inline constexpr
+		bool compare (const Vec<V1>&, const Vec<V2>&, ComposeOp, Op, bt::number_seq<size_type>) {
+			return LastVal;
+		}
+		template <bool LastVal, typename V1, typename V2, typename ComposeOp, typename Op, size_type I1, size_type... I>
+		inline
+		bool compare (const Vec<V1>& parLeft, const Vec<V2>& parRight, ComposeOp parComposeOp, Op parOp, bt::number_seq<size_type, I1, I...>) {
+			static_assert(I1 < VectorWrapperInfo<V1>::dimensions, "Index out of range");
+			return parComposeOp(parOp(parLeft[I1], parRight[I1]), compare<LastVal>(parLeft, parRight, parComposeOp, parOp, bt::number_seq<size_type, I...>()));
+		}
 	} //namespace implem
 
 	template <typename V> const Vec<V, 1> Vec<V, 1>::unit_x(scalar_type(1));
@@ -296,20 +308,74 @@ namespace vwr {
 	template <typename V1, typename V2>
 	inline bool operator== (const Vec<V1>& parLeft, const Vec<V2>& parRight) {
 		static_assert(static_cast<int>(VectorWrapperInfo<V1>::dimensions) == static_cast<int>(VectorWrapperInfo<V2>::dimensions), "Dimensions mismatch");
-		bool retval = true;
-		for (size_type z = 0; z < VectorWrapperInfo<V1>::dimensions; ++z) {
-			retval &= (parLeft[z] == parRight[z]);
-		}
-		return retval;
+		typedef typename std::common_type<typename VectorWrapperInfo<V1>::scalar_type, typename VectorWrapperInfo<V2>::scalar_type>::type scalar_type;
+		return implem::compare<true>(
+			parLeft,
+			parRight,
+			std::logical_and<bool>(),
+			std::equal_to<scalar_type>(),
+			bt::number_range<size_type, 0, VectorWrapperInfo<V1>::dimensions>()
+		);
+	}
+	template <typename V1, typename V2>
+	inline bool operator!= (const Vec<V1>& parLeft, const Vec<V2>& parRight) {
+		static_assert(static_cast<int>(VectorWrapperInfo<V1>::dimensions) == static_cast<int>(VectorWrapperInfo<V2>::dimensions), "Dimensions mismatch");
+		typedef typename std::common_type<typename VectorWrapperInfo<V1>::scalar_type, typename VectorWrapperInfo<V2>::scalar_type>::type scalar_type;
+		return implem::compare<false>(
+			parLeft,
+			parRight,
+			std::logical_or<bool>(),
+			std::not_equal_to<scalar_type>(),
+			bt::number_range<size_type, 0, VectorWrapperInfo<V1>::dimensions>()
+		);
 	}
 	template <typename V1, typename V2>
 	inline bool operator< (const Vec<V1>& parLeft, const Vec<V2>& parRight) {
 		static_assert(static_cast<int>(VectorWrapperInfo<V1>::dimensions) == static_cast<int>(VectorWrapperInfo<V2>::dimensions), "Dimensions mismatch");
-		bool retval = true;
-		for (size_type z = 0; z < VectorWrapperInfo<V1>::dimensions; ++z) {
-			retval &= (parLeft[z] < parRight[z]);
-		}
-		return retval;
+		typedef typename std::common_type<typename VectorWrapperInfo<V1>::scalar_type, typename VectorWrapperInfo<V2>::scalar_type>::type scalar_type;
+		return implem::compare<true>(
+			parLeft,
+			parRight,
+			std::logical_and<bool>(),
+			std::less<scalar_type>(),
+			bt::number_range<size_type, 0, VectorWrapperInfo<V1>::dimensions>()
+		);
+	}
+	template <typename V1, typename V2>
+	inline bool operator> (const Vec<V1>& parLeft, const Vec<V2>& parRight) {
+		static_assert(static_cast<int>(VectorWrapperInfo<V1>::dimensions) == static_cast<int>(VectorWrapperInfo<V2>::dimensions), "Dimensions mismatch");
+		typedef typename std::common_type<typename VectorWrapperInfo<V1>::scalar_type, typename VectorWrapperInfo<V2>::scalar_type>::type scalar_type;
+		return implem::compare<true>(
+			parLeft,
+			parRight,
+			std::logical_and<bool>(),
+			std::greater<scalar_type>(),
+			bt::number_range<size_type, 0, VectorWrapperInfo<V1>::dimensions>()
+		);
+	}
+	template <typename V1, typename V2>
+	inline bool operator<= (const Vec<V1>& parLeft, const Vec<V2>& parRight) {
+		static_assert(static_cast<int>(VectorWrapperInfo<V1>::dimensions) == static_cast<int>(VectorWrapperInfo<V2>::dimensions), "Dimensions mismatch");
+		typedef typename std::common_type<typename VectorWrapperInfo<V1>::scalar_type, typename VectorWrapperInfo<V2>::scalar_type>::type scalar_type;
+		return implem::compare<true>(
+			parLeft,
+			parRight,
+			std::logical_and<bool>(),
+			std::less_equal<scalar_type>(),
+			bt::number_range<size_type, 0, VectorWrapperInfo<V1>::dimensions>()
+		);
+	}
+	template <typename V1, typename V2>
+	inline bool operator>= (const Vec<V1>& parLeft, const Vec<V2>& parRight) {
+		static_assert(static_cast<int>(VectorWrapperInfo<V1>::dimensions) == static_cast<int>(VectorWrapperInfo<V2>::dimensions), "Dimensions mismatch");
+		typedef typename std::common_type<typename VectorWrapperInfo<V1>::scalar_type, typename VectorWrapperInfo<V2>::scalar_type>::type scalar_type;
+		return implem::compare<true>(
+			parLeft,
+			parRight,
+			std::logical_and<bool>(),
+			std::greater_equal<scalar_type>(),
+			bt::number_range<size_type, 0, VectorWrapperInfo<V1>::dimensions>()
+		);
 	}
 
 	template <typename V>
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index cbbf5c0..0b65f23 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -6,6 +6,7 @@ add_executable(${PROJECT_NAME}
 	test_ops.cpp
 	example.cpp
 	test_get_at.cpp
+	test_operators.cpp
 )
 
 target_link_libraries(${PROJECT_NAME}
diff --git a/test/unit/test_operators.cpp b/test/unit/test_operators.cpp
new file mode 100644
index 0000000..f773616
--- /dev/null
+++ b/test/unit/test_operators.cpp
@@ -0,0 +1,44 @@
+#include "sample_vectors.hpp"
+#include <gtest/gtest.h>
+
+TEST(vwr, operators) {
+	using namespace vwr;
+
+	{
+		ivec3 a(5, 6, 7);
+		ivec3 b(6, 7, 8);
+
+		EXPECT_LT(a, b);
+		EXPECT_LE(a, b);
+		EXPECT_NE(a, b);
+		EXPECT_FALSE(a == b);
+		EXPECT_FALSE(a > b);
+		EXPECT_FALSE(a >= b);
+		EXPECT_GT(b, a);
+		EXPECT_GE(b, a);
+	}
+	{
+		ivec3 a(6, 6, 7);
+		ivec3 b(6, 7, 8);
+
+		EXPECT_FALSE(a < b);
+		EXPECT_LE(a, b);
+		EXPECT_NE(a, b);
+		EXPECT_FALSE(a == b);
+		EXPECT_FALSE(a > b);
+		EXPECT_FALSE(a >= b);
+		EXPECT_GE(b, a);
+	}
+	{
+		ivec3 a(0xAABB, 0xAABB, 0xAABB);
+		ivec3 b(0xAABB, 0xAABB, 0xAABB);
+
+		EXPECT_FALSE(a < b);
+		EXPECT_LE(a, b);
+		EXPECT_FALSE(a != b);
+		EXPECT_EQ(a, b);
+		EXPECT_FALSE(a > b);
+		EXPECT_GE(a, b);
+		EXPECT_GE(b, a);
+	}
+}