diff --git a/headers/vmath.hpp/vmath_fun.hpp b/headers/vmath.hpp/vmath_fun.hpp
index e2d0bc1..a876bd8 100644
--- a/headers/vmath.hpp/vmath_fun.hpp
+++ b/headers/vmath.hpp/vmath_fun.hpp
@@ -431,25 +431,36 @@ namespace vmath_hpp
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
     constexpr equal_to(T x, T y) noexcept {
-        return x == y;
+        if constexpr ( std::is_floating_point_v<T> ) {
+            // http://www.realtimecollisiondetection.net/pubs/Tolerances
+            const T epsilon = std::numeric_limits<T>::epsilon();
+            return abs(x - y) <= epsilon * max(T(1), abs(x), abs(y));
+        } else {
+            return x == y;
+        }
     }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
     constexpr equal_to(T x, T y, T epsilon) noexcept {
-        return abs(x - y) <= epsilon;
+        if constexpr ( std::is_floating_point_v<T> ) {
+            // http://www.realtimecollisiondetection.net/pubs/Tolerances
+            return abs(x - y) <= epsilon * max(T(1), abs(x), abs(y));
+        } else {
+            return abs(x - y) <= epsilon;
+        }
     }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
     constexpr not_equal_to(T x, T y) noexcept {
-        return x != y;
+        return !equal_to(x, y);
     }
 
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_arithmetic_v<T>, bool>
     constexpr not_equal_to(T x, T y, T epsilon) noexcept {
-        return abs(x - y) > epsilon;
+        return !equal_to(x, y, epsilon);
     }
 
     template < typename T >
diff --git a/headers/vmath.hpp/vmath_fwd.hpp b/headers/vmath.hpp/vmath_fwd.hpp
index 2b330ec..b7a247a 100644
--- a/headers/vmath.hpp/vmath_fwd.hpp
+++ b/headers/vmath.hpp/vmath_fwd.hpp
@@ -11,6 +11,7 @@
 
 #include <initializer_list>
 #include <iterator>
+#include <limits>
 #include <stdexcept>
 #include <type_traits>
 #include <utility>
diff --git a/headers/vmath.hpp/vmath_mat_fun.hpp b/headers/vmath.hpp/vmath_mat_fun.hpp
index 5960e2e..8ffc74e 100644
--- a/headers/vmath.hpp/vmath_mat_fun.hpp
+++ b/headers/vmath.hpp/vmath_mat_fun.hpp
@@ -340,9 +340,7 @@ namespace vmath_hpp
 
     template < typename T, std::size_t Size >
     [[nodiscard]] constexpr bool operator!=(const mat<T, Size>& xs, const mat<T, Size>& ys) {
-        return fold_join([](bool acc, const vec<T, Size>& x, const vec<T, Size>& y){
-            return acc || (x != y);
-        }, false, xs, ys);
+        return !(xs == ys);
     }
 
     // operator<
diff --git a/headers/vmath.hpp/vmath_vec_fun.hpp b/headers/vmath.hpp/vmath_vec_fun.hpp
index ce780dc..d9a0cf0 100644
--- a/headers/vmath.hpp/vmath_vec_fun.hpp
+++ b/headers/vmath.hpp/vmath_vec_fun.hpp
@@ -293,15 +293,13 @@ namespace vmath_hpp
     template < typename T, std::size_t Size >
     [[nodiscard]] constexpr bool operator==(const vec<T, Size>& xs, const vec<T, Size>& ys) {
         return fold_join([](bool acc, T x, T y){
-            return acc && (x == y);
+            return acc && equal_to(x, y);
         }, true, xs, ys);
     }
 
     template < typename T, std::size_t Size >
     [[nodiscard]] constexpr bool operator!=(const vec<T, Size>& xs, const vec<T, Size>& ys) {
-        return fold_join([](bool acc, T x, T y){
-            return acc || (x != y);
-        }, false, xs, ys);
+        return !(xs == ys);
     }
 
     // operator<
diff --git a/untests/vmath_fun_tests.cpp b/untests/vmath_fun_tests.cpp
index 2cfcb89..33787f5 100644
--- a/untests/vmath_fun_tests.cpp
+++ b/untests/vmath_fun_tests.cpp
@@ -165,6 +165,16 @@ TEST_CASE("vmath/fun") {
         STATIC_REQUIRE_FALSE(not_equal_to(1, 1, 0));
         STATIC_REQUIRE_FALSE(not_equal_to(1, 1, 1));
 
+        STATIC_REQUIRE(equal_to(1.f, 1.f + std::numeric_limits<float>::epsilon() * 0.5f));
+        STATIC_REQUIRE_FALSE(equal_to(1.f, 1.f + std::numeric_limits<float>::epsilon() * 1.5f));
+        STATIC_REQUIRE(equal_to(100.f, 100.f + std::numeric_limits<float>::epsilon() * 90.f));
+        STATIC_REQUIRE_FALSE(equal_to(100.f, 100.f + std::numeric_limits<float>::epsilon() * 110.f));
+
+        STATIC_REQUIRE_FALSE(not_equal_to(1.f, 1.f + std::numeric_limits<float>::epsilon() * 0.5f));
+        STATIC_REQUIRE(not_equal_to(1.f, 1.f + std::numeric_limits<float>::epsilon() * 1.5f));
+        STATIC_REQUIRE_FALSE(not_equal_to(100.f, 100.f + std::numeric_limits<float>::epsilon() * 90.f));
+        STATIC_REQUIRE(not_equal_to(100.f, 100.f + std::numeric_limits<float>::epsilon() * 110.f));
+
         STATIC_REQUIRE_FALSE(any(false));
         STATIC_REQUIRE_FALSE(any(0));
         STATIC_REQUIRE(any(true));