diff --git a/headers/vmath.hpp/vmath.hpp b/headers/vmath.hpp/vmath.hpp
index 6c98170..25faaf7 100644
--- a/headers/vmath.hpp/vmath.hpp
+++ b/headers/vmath.hpp/vmath.hpp
@@ -6,6 +6,7 @@
 
 #pragma once
 
+#include "vmath_fun.hpp"
 #include "vmath_mat.hpp"
-#include "vmath_vec_ops.hpp"
 #include "vmath_vec.hpp"
+#include "vmath_vec_fun.hpp"
diff --git a/headers/vmath.hpp/vmath_vec_fun.hpp b/headers/vmath.hpp/vmath_vec_fun.hpp
new file mode 100644
index 0000000..c605411
--- /dev/null
+++ b/headers/vmath.hpp/vmath_vec_fun.hpp
@@ -0,0 +1,544 @@
+/*******************************************************************************
+ * This file is part of the "https://github.com/blackmatov/vmath.hpp"
+ * For conditions of distribution and use, see copyright notice in LICENSE.md
+ * Copyright (C) 2020, by Matvey Cherevko (blackmatov@gmail.com)
+ ******************************************************************************/
+
+#pragma once
+
+#include "vmath_fwd.hpp"
+
+#include "vmath_fun.hpp"
+#include "vmath_vec.hpp"
+
+namespace vmath_hpp::detail
+{
+    namespace impl
+    {
+        template < typename A, std::size_t Size, typename F, std::size_t... Is >
+        constexpr auto map_impl(F&& f, const vec<A, Size>& v, std::index_sequence<Is...>)
+            -> vec<std::invoke_result_t<F, A>, Size>
+        {
+            return { f(v[Is])... };
+        }
+
+        template < typename A, typename B, std::size_t Size, typename F, std::size_t... Is >
+        constexpr auto zip_impl(F&& f, const vec<A, Size>& a, const vec<B, Size>& b, std::index_sequence<Is...>)
+            -> vec<std::invoke_result_t<F, A, B>, Size>
+        {
+            return { f(a[Is], b[Is])... };
+        }
+
+        template < typename A, typename B, typename C, std::size_t Size, typename F, std::size_t... Is >
+        constexpr auto zip_impl(F&& f, const vec<A, Size>& a, const vec<B, Size>& b, const vec<C, Size>& c, std::index_sequence<Is...>)
+            -> vec<std::invoke_result_t<F, A, B, C>, Size>
+        {
+            return { f(a[Is], b[Is], c[Is])... };
+        }
+
+        template < typename A, typename B, std::size_t Size, typename F, std::size_t... Is >
+        constexpr auto fold_impl(F&& f, A init, const vec<B, Size>& v, std::index_sequence<Is...>)
+            -> A
+        {
+            return ((init = f(std::move(init), v[Is])), ...);
+        }
+    }
+
+    template < typename T, std::size_t Size, typename F >
+    constexpr auto map(F&& f, const vec<T, Size>& v) {
+        return impl::map_impl(std::forward<F>(f), v, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F >
+    constexpr auto zip(F&& f, const vec<A, Size>& a, const vec<B, Size>& b) {
+        return impl::zip_impl(std::forward<F>(f), a, b, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, typename C, std::size_t Size, typename F >
+    constexpr auto zip(F&& f, const vec<A, Size>& a, const vec<B, Size>& b, const vec<C, Size>& c) {
+        return impl::zip_impl(std::forward<F>(f), a, b, c, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F >
+    constexpr auto fold(F&& f, A init, const vec<B, Size>& v) {
+        return impl::fold_impl(std::forward<F>(f), std::move(init), v, std::make_index_sequence<Size>{});
+    }
+}
+
+//
+// Operators
+//
+
+namespace vmath_hpp
+{
+    // -operator
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> operator-(const vec<T, Size>& xs) {
+        return map(std::negate<>(), xs);
+    }
+
+    // operator+
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> operator+(const vec<T, Size>& xs, T y) {
+        return map([y](T x){ return x + y; }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> operator+(T x, const vec<T, Size>& ys) {
+        return map([x](T y){ return x + y; }, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> operator+(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return zip(std::plus<>(), xs, ys);
+    }
+
+    // operator-
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> operator-(const vec<T, Size>& xs, T y) {
+        return map([y](T x){ return x - y; }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> operator-(T x, const vec<T, Size>& ys) {
+        return map([x](T y){ return x - y; }, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> operator-(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return zip(std::minus<>(), xs, ys);
+    }
+
+    // operator*
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> operator*(const vec<T, Size>& xs, T y) {
+        return map([y](T x){ return x * y; }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> operator*(T x, const vec<T, Size>& ys) {
+        return map([x](T y){ return x * y; }, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> operator*(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return zip(std::multiplies<>(), xs, ys);
+    }
+
+    // operator/
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> operator/(const vec<T, Size>& xs, T y) {
+        return map([y](T x){ return x / y; }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> operator/(T x, const vec<T, Size>& ys) {
+        return map([x](T y){ return x / y; }, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> operator/(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return zip(std::divides<>(), xs, ys);
+    }
+}
+
+//
+// Angle and Trigonometry Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> radians(const vec<T, Size>& degrees) {
+        return map([](T x) { return radians(x); }, degrees);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> degrees(const vec<T, Size>& radians) {
+        return map([](T x) { return degrees(x); }, radians);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> sin(const vec<T, Size>& xs) {
+        return map([](T x) { return sin(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> cos(const vec<T, Size>& xs) {
+        return map([](T x) { return cos(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> tan(const vec<T, Size>& xs) {
+        return map([](T x) { return tan(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> asin(const vec<T, Size>& xs) {
+        return map([](T x) { return asin(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> acos(const vec<T, Size>& xs) {
+        return map([](T x) { return acos(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> atan(const vec<T, Size>& xs) {
+        return map([](T x) { return atan(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> atan2(const vec<T, Size>& ys, const vec<T, Size>& xs) {
+        return zip([](T y, T x) { return atan2(y, x); }, ys, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> sinh(const vec<T, Size>& xs) {
+        return map([](T x) { return sinh(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> cosh(const vec<T, Size>& xs) {
+        return map([](T x) { return cosh(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> tanh(const vec<T, Size>& xs) {
+        return map([](T x) { return tanh(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> asinh(const vec<T, Size>& xs) {
+        return map([](T x) { return asinh(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> acosh(const vec<T, Size>& xs) {
+        return map([](T x) { return acosh(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> atanh(const vec<T, Size>& xs) {
+        return map([](T x) { return atanh(x); }, xs);
+    }
+}
+
+//
+// Exponential Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T, std::size_t Size >
+    vec<T, Size> pow(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return zip([](T x, T y) { return pow(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> exp(const vec<T, Size>& xs) {
+        return map([](T x) { return exp(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> log(const vec<T, Size>& xs) {
+        return map([](T x) { return log(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> exp2(const vec<T, Size>& xs) {
+        return map([](T x) { return exp2(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> log2(const vec<T, Size>& xs) {
+        return map([](T x) { return log2(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> sqrt(const vec<T, Size>& xs) {
+        return map([](T x) { return sqrt(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> invsqrt(const vec<T, Size>& xs) {
+        return map([](T x) { return invsqrt(x); }, xs);
+    }
+}
+
+//
+// Common Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T, std::size_t Size >
+    vec<T, Size> abs(const vec<T, Size>& xs) {
+        return map([](T x) { return abs(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> sign(const vec<T, Size>& xs) {
+        return map([](T x) { return sign(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> floor(const vec<T, Size>& xs) {
+        return map([](T x) { return floor(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> trunc(const vec<T, Size>& xs) {
+        return map([](T x) { return trunc(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> round(const vec<T, Size>& xs) {
+        return map([](T x) { return round(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> ceil(const vec<T, Size>& xs) {
+        return map([](T x) { return ceil(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> fract(const vec<T, Size>& xs) {
+        return map([](T x) { return fract(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> fmod(const vec<T, Size>& xs, T y) {
+        return map([y](T x) { return fmod(x, y); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> fmod(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return zip([](T x, T y) { return fmod(x, y); }, xs, ys);
+    }
+
+    namespace impl
+    {
+        template < typename T, std::size_t Size, std::size_t... Is >
+        vec<T, Size> modf_impl(const vec<T, Size>& xs, vec<T, Size>* is, std::index_sequence<Is...>) {
+            return { modf(xs[Is], &(*is)[Is])... };
+        }
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> modf(const vec<T, Size>& xs, vec<T, Size>* is) {
+        return impl::modf_impl(xs, is, std::make_index_sequence<Size>{});
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> min(const vec<T, Size>& xs, T y) {
+        return map([y](T x) { return min(x, y); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> min(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return zip([](T x, T y) { return min(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> max(const vec<T, Size>& xs, T y) {
+        return map([y](T x) { return max(x, y); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> max(const vec<T, Size>& xs, const vec<T, Size>& ys) {
+        return zip([](T x, T y) { return max(x, y); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> clamp(const vec<T, Size>& xs, T min_x, T max_x) {
+        return map([min_x, max_x](T x) { return clamp(x, min_x, max_x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> clamp(const vec<T, Size>& xs, const vec<T, Size>& min_xs, const vec<T, Size>& max_xs) {
+        return zip([](T x, T min_x, T max_x) { return clamp(x, min_x, max_x); }, xs, min_xs, max_xs);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> mix(const vec<T, Size>& xs, const vec<T, Size>& ys, T a) {
+        return zip([a](T x, T y) { return mix(x, y, a); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> mix(const vec<T, Size>& xs, const vec<T, Size>& ys, bool a) {
+        return zip([a](T x, T y) { return mix(x, y, a); }, xs, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> mix(const vec<T, Size>& xs, const vec<T, Size>& ys, const vec<T, Size>& as) {
+        return zip([](T x, T y, T a) { return mix(x, y, a); }, xs, ys, as);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> mix(const vec<T, Size>& xs, const vec<T, Size>& ys, const vec<bool, Size>& as) {
+        return zip([](T x, T y, bool a) { return mix(x, y, a); }, xs, ys, as);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> step(T edge, const vec<T, Size>& xs) {
+        return map([edge](T x) { return step(edge, x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> step(const vec<T, Size>& edges, const vec<T, Size>& xs) {
+        return zip([](T edge, T x) { return step(edge, x); }, edges, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> smoothstep(T edge0, T edge1, const vec<T, Size>& xs) {
+        return map([edge0, edge1](T x) { return smoothstep(edge0, edge1, x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<T, Size> smoothstep(const vec<T, Size>& edges0, const vec<T, Size>& edges1, const vec<T, Size>& xs) {
+        return zip([](T edge0, T edge1, T x) { return smoothstep(edge0, edge1, x); }, edges0, edges1, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<bool, Size> isnan(const vec<T, Size>& xs) {
+        return map([](T x) { return isnan(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<bool, Size> isinf(const vec<T, Size>& xs) {
+        return map([](T x) { return isinf(x); }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> fma(const vec<T, Size>& as, const vec<T, Size>& bs, const vec<T, Size>& cs) {
+        return zip([](T a, T b, T c) { return fma(a, b, c); }, as, bs, cs);
+    }
+
+    namespace impl
+    {
+        template < typename T, std::size_t Size, std::size_t... Is >
+        vec<T, Size> frexp_impl(const vec<T, Size>& xs, vec<int, Size>* exps, std::index_sequence<Is...>) {
+            return { frexp(xs.data[Is], &exps->data[Is])... };
+        }
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> frexp(const vec<T, Size>& xs, vec<int, Size>* exps) {
+        return impl::frexp_impl(xs, exps, std::make_index_sequence<Size>{});
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> ldexp(const vec<T, Size>& xs, const vec<int, Size>& exps) {
+        return zip([](T x, int exp) { return ldexp(x, exp); }, xs, exps);
+    }
+}
+
+//
+// Geometric Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T, std::size_t Size >
+    T length(const vec<T, Size>& x) noexcept {
+        return sqrt(dot(x, x));
+    }
+
+    template < typename T, std::size_t Size >
+    T distance(const vec<T, Size>& p0, const vec<T, Size>& p1) noexcept {
+        return length(p0 - p1);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr T dot(const vec<T, Size>& x, const vec<T, Size>& y) noexcept {
+        return fold(std::plus<>(), T(0), zip(std::multiplies<>(), x, y));
+    }
+
+    template < typename T >
+    constexpr vec<T, 3> cross(const vec<T, 3>& x, const vec<T, 3>& y) noexcept {
+        return {
+            x.y * y.z - x.z * y.y,
+            x.z * y.x - x.x * y.z,
+            x.x * y.y - x.y * y.x};
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> normalize(const vec<T, Size>& x) noexcept {
+        return x * invsqrt(dot(x, x));
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> faceforward(const vec<T, Size>& n, const vec<T, Size>& i, const vec<T, Size>& nref) noexcept {
+        return dot(nref, i) < T(0) ? n : -n;
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> reflect(const vec<T, Size>& i, const vec<T, Size>& n) noexcept {
+        return i - n * dot(n, i) * T(2);
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> refract(const vec<T, Size>& i, const vec<T, Size>& n, T eta) noexcept {
+        const T d = dot(n, i);
+        const T k = T(1) - eta * eta * (T(1) - d * d);
+        return T(k >= T(0)) * (eta * i - (eta * d + sqrt(k)) * n);
+    }
+}
+
+//
+// Vector Relational Functions
+//
+
+namespace vmath_hpp
+{
+    template < typename T, std::size_t Size >
+    constexpr vec<bool, Size> less(const vec<T, Size>& x, const vec<T, Size>& y) {
+        return detail::zip(std::less<>(), x, y);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<bool, Size> less_equal(const vec<T, Size>& x, const vec<T, Size>& y) {
+        return detail::zip(std::less_equal<>(), x, y);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<bool, Size> greater(const vec<T, Size>& x, const vec<T, Size>& y) {
+        return detail::zip(std::greater<>(), x, y);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<bool, Size> greater_equal(const vec<T, Size>& x, const vec<T, Size>& y) {
+        return detail::zip(std::greater_equal<>(), x, y);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<bool, Size> equal_to(const vec<T, Size>& x, const vec<T, Size>& y) {
+        return detail::zip(std::equal_to<>(), x, y);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr vec<bool, Size> not_equal_to(const vec<T, Size>& x, const vec<T, Size>& y) {
+        return detail::zip(std::not_equal_to<>(), x, y);
+    }
+
+    template < std::size_t Size >
+    constexpr bool any(const vec<bool, Size>& x) {
+        return detail::fold(std::logical_or<>(), false, x);
+    }
+
+    template < std::size_t Size >
+    constexpr bool all(const vec<bool, Size>& x) {
+        return detail::fold(std::logical_and<>(), true, x);
+    }
+
+    template < std::size_t Size >
+    constexpr vec<bool, Size> not_(const vec<bool, Size>& x) {
+        return detail::map(std::logical_not<>(), x);
+    }
+}
diff --git a/headers/vmath.hpp/vmath_vec_ops.hpp b/headers/vmath.hpp/vmath_vec_ops.hpp
deleted file mode 100644
index a20bfd5..0000000
--- a/headers/vmath.hpp/vmath_vec_ops.hpp
+++ /dev/null
@@ -1,91 +0,0 @@
-/*******************************************************************************
- * This file is part of the "https://github.com/blackmatov/vmath.hpp"
- * For conditions of distribution and use, see copyright notice in LICENSE.md
- * Copyright (C) 2020, by Matvey Cherevko (blackmatov@gmail.com)
- ******************************************************************************/
-
-#pragma once
-
-#include "vmath_fwd.hpp"
-#include "vmath_vec.hpp"
-
-namespace vmath_hpp::detail
-{
-    template < typename T, std::size_t Size, typename F >
-    constexpr auto map(const vec<T, Size>& v, F&& f) {
-        vec<std::invoke_result_t<F, T>, Size> result;
-        for ( std::size_t i = 0; i < Size; ++i ) {
-            result[i] = f(v[i]);
-        }
-        return result;
-    }
-
-    template < typename T, typename U, std::size_t Size, typename F >
-    constexpr auto zip(const vec<T, Size>& l, const vec<U, Size>& r, F&& f) {
-        vec<std::invoke_result_t<F, T, U>, Size> result;
-        for ( std::size_t i = 0; i < Size; ++i ) {
-            result[i] = f(l[i], r[i]);
-        }
-        return result;
-    }
-
-    template < typename T, typename U, std::size_t Size, typename F >
-    constexpr auto fold(const vec<T, Size>& v, U init, F&& f) {
-        for ( std::size_t i = 0; i < Size; ++i ) {
-            init = f(init, v[i]);
-        }
-        return init;
-    }
-}
-
-//
-// Vector Relational Functions
-//
-
-namespace vmath_hpp
-{
-    template < typename T, std::size_t Size >
-    constexpr vec<bool, Size> less(const vec<T, Size>& x, const vec<T, Size>& y) {
-        return detail::zip(x, y, std::less<>());
-    }
-
-    template < typename T, std::size_t Size >
-    constexpr vec<bool, Size> less_equal(const vec<T, Size>& x, const vec<T, Size>& y) {
-        return detail::zip(x, y, std::less_equal<>());
-    }
-
-    template < typename T, std::size_t Size >
-    constexpr vec<bool, Size> greater(const vec<T, Size>& x, const vec<T, Size>& y) {
-        return detail::zip(x, y, std::greater<>());
-    }
-
-    template < typename T, std::size_t Size >
-    constexpr vec<bool, Size> greater_equal(const vec<T, Size>& x, const vec<T, Size>& y) {
-        return detail::zip(x, y, std::greater_equal<>());
-    }
-
-    template < typename T, std::size_t Size >
-    constexpr vec<bool, Size> equal(const vec<T, Size>& x, const vec<T, Size>& y) {
-        return detail::zip(x, y, std::equal_to<>());
-    }
-
-    template < typename T, std::size_t Size >
-    constexpr vec<bool, Size> not_equal(const vec<T, Size>& x, const vec<T, Size>& y) {
-        return detail::zip(x, y, std::not_equal_to<>());
-    }
-
-    template < std::size_t Size >
-    constexpr bool any(const vec<bool, Size>& x) {
-        return detail::fold(x, false, std::logical_or<>());
-    }
-
-    template < std::size_t Size >
-    constexpr bool all(const vec<bool, Size>& x) {
-        return detail::fold(x, true, std::logical_and<>());
-    }
-
-    template < std::size_t Size >
-    constexpr vec<bool, Size> not_(const vec<bool, Size>& x) {
-        return detail::map(x, std::logical_not<>());
-    }
-}
diff --git a/untests/vmath_vec_fun_tests.cpp b/untests/vmath_vec_fun_tests.cpp
new file mode 100644
index 0000000..b1b5730
--- /dev/null
+++ b/untests/vmath_vec_fun_tests.cpp
@@ -0,0 +1,205 @@
+/*******************************************************************************
+ * This file is part of the "https://github.com/blackmatov/vmath.hpp"
+ * For conditions of distribution and use, see copyright notice in LICENSE.md
+ * Copyright (C) 2020, by Matvey Cherevko (blackmatov@gmail.com)
+ ******************************************************************************/
+
+#include <vmath.hpp/vmath_vec_fun.hpp>
+
+#define CATCH_CONFIG_FAST_COMPILE
+#include <catch2/catch.hpp>
+
+namespace
+{
+    using namespace vmath_hpp;
+
+    template < typename T >
+    class approx2 {
+    public:
+        constexpr explicit approx2(T v) : value_(v) {}
+        constexpr explicit approx2(T x, T y) : value_(x, y) {}
+
+        friend constexpr bool operator==(const vec<T, 2>& l, const approx2& r) {
+            return (r.value_.x < l.x + epsilon)
+                && (l.x < r.value_.x + epsilon)
+                && (r.value_.y < l.y + epsilon)
+                && (l.y < r.value_.y + epsilon);
+        }
+    private:
+        vec<T, 2> value_;
+        static constexpr T epsilon = std::numeric_limits<T>::epsilon() * 100;
+    };
+
+    template < typename T >
+    class approx3 {
+    public:
+        constexpr explicit approx3(T v) : value_(v) {}
+        constexpr explicit approx3(T x, T y, T z) : value_(x, y, z) {}
+
+        friend constexpr bool operator==(const vec<T, 3>& l, const approx3& r) {
+            return (r.value_.x < l.x + epsilon)
+                && (l.x < r.value_.x + epsilon)
+                && (r.value_.y < l.y + epsilon)
+                && (l.y < r.value_.y + epsilon)
+                && (r.value_.z < l.z + epsilon)
+                && (l.z < r.value_.z + epsilon);
+        }
+    private:
+        vec<T, 3> value_;
+        static constexpr T epsilon = std::numeric_limits<T>::epsilon() * 100;
+    };
+}
+
+TEST_CASE("vmath/vec_fun") {
+    SECTION("Operators") {
+        STATIC_REQUIRE(-vec2i(1,-2) == vec2i(-1,2));
+
+        STATIC_REQUIRE(vec2i(1,2) + 3 == vec2i(4,5));
+        STATIC_REQUIRE(vec2i(1,2) - 3 == vec2i(-2,-1));
+        STATIC_REQUIRE(vec2i(1,2) * 3 == vec2i(3,6));
+        STATIC_REQUIRE(vec2i(2,4) / 2 == vec2i(1,2));
+
+        STATIC_REQUIRE(3 + vec2i(1,2) == vec2i(4,5));
+        STATIC_REQUIRE(3 - vec2i(1,2) == vec2i(2,1));
+        STATIC_REQUIRE(3 * vec2i(1,2) == vec2i(3,6));
+        STATIC_REQUIRE(4 / vec2i(2,4) == vec2i(2,1));
+
+        STATIC_REQUIRE(vec2i(1,2) + vec2i(3,4) == vec2i(4,6));
+        STATIC_REQUIRE(vec2i(1,2) - vec2i(3,4) == vec2i(-2,-2));
+        STATIC_REQUIRE(vec2i(1,2) * vec2i(3,4) == vec2i(3,8));
+        STATIC_REQUIRE(vec2i(3,4) / vec2i(1,2) == vec2i(3,2));
+    }
+
+    SECTION("Angle and Trigonometry Functions") {
+        STATIC_REQUIRE(radians(degrees(vec2f(12.13f))) == approx2(12.13f));
+        STATIC_REQUIRE(degrees(radians(vec2f(12.13f))) == approx2(12.13f));
+
+        sin(vec2f(1.f));
+        cos(vec2f(1.f));
+        tan(vec2f(1.f));
+
+        asin(vec2f(1.f));
+        acos(vec2f(1.f));
+        atan(vec2f(1.f));
+        atan2(vec2f(1.f), vec2f(1.f));
+
+        sinh(vec2f(1.f));
+        cosh(vec2f(1.f));
+        tanh(vec2f(1.f));
+
+        asinh(vec2f(1.f));
+        acosh(vec2f(1.f));
+        atanh(vec2f(1.f));
+    }
+
+    SECTION("Exponential Functions") {
+        pow(vec2f(1.f), vec2f(2.f));
+        exp(vec2f(1.f));
+        log(vec2f(1.f));
+        exp2(vec2f(1.f));
+        log2(vec2f(1.f));
+        sqrt(vec2f(1.f));
+        invsqrt(vec2f(1.f));
+    }
+
+    SECTION("Common Functions") {
+        REQUIRE(abs(vec2f(1.f, -1.f)) == approx2(1.f,1.f));
+        REQUIRE(sign(vec3f(1.f, -1.f, 0.f)) == approx3(1.f,-1.f,0.f));
+
+        floor(vec2f(1.f, -1.f));
+        trunc(vec2f(1.f, -1.f));
+        round(vec2f(1.f, -1.f));
+        ceil(vec2f(1.f, -1.f));
+        fract(vec2f(1.f, -1.f));
+
+        REQUIRE(fmod(vec2f(1.7f), 1.2f) == approx2(0.5f));
+        REQUIRE(fmod(vec2f(1.7f), vec2f(1.2f)) == approx2(0.5f));
+
+        {
+            vec2f out_i{};
+            REQUIRE(modf(vec2f(1.7f), &out_i) == approx2(0.7f));
+            REQUIRE(out_i.x == Approx(1.f));
+        }
+
+        STATIC_REQUIRE(min(vec2i(1,2), 1) == vec2i(1,1));
+        STATIC_REQUIRE(min(vec2i(1,1), vec2i(0,2)) == vec2i(0,1));
+
+        STATIC_REQUIRE(max(vec2i(1,2), 1) == vec2i(1,2));
+        STATIC_REQUIRE(max(vec2i(1,1), vec2i(0,2)) == vec2i(1,2));
+
+        STATIC_REQUIRE(clamp(vec2i(1,2), 0, 1) == vec2i(1,1));
+        STATIC_REQUIRE(clamp(vec2i(1,2), vec2i(0), vec2i(1)) == vec2i(1,1));
+
+        STATIC_REQUIRE(mix(vec2f(0.f), vec2f(10.f), 0.5f) == approx2(5.f));
+        STATIC_REQUIRE(mix(vec2f(0.f), vec2f(10.f), vec2f(0.5f)) == approx2(5.f));
+
+        STATIC_REQUIRE(mix(vec2f(0.f), vec2f(10.f), true) == approx2(10.f));
+        STATIC_REQUIRE(mix(vec2f(0.f), vec2f(10.f), false) == approx2(0.f));
+
+        STATIC_REQUIRE(mix(vec2f(0.f,5.f), vec2f(10.f,20.f), vec2f(true,false)) == approx2(10.f, 5.f));
+        STATIC_REQUIRE(mix(vec2f(0.f,5.f), vec2f(10.f,20.f), vec2f(true,false)) == approx2(10.f, 5.f));
+
+        STATIC_REQUIRE(step(0.5f, vec2f(0.4f)) == approx2(0.f));
+        STATIC_REQUIRE(step(0.5f, vec2f(0.6f)) == approx2(1.f));
+        STATIC_REQUIRE(step(vec2f(0.5f), vec2f(0.4f)) == approx2(0.f));
+        STATIC_REQUIRE(step(vec2f(0.5f), vec2f(0.6f)) == approx2(1.f));
+
+        STATIC_REQUIRE(smoothstep(0.f, 1.f, vec2f(0.1f)) == approx2(0.028f));
+        STATIC_REQUIRE(smoothstep(vec2f(0.f), vec2f(1.f), vec2f(0.1f)) == approx2(0.028f));
+
+        REQUIRE_FALSE(isnan(vec2f(1.f)).x);
+        REQUIRE_FALSE(isinf(vec2f(1.f)).x);
+
+        REQUIRE_FALSE(fma(vec2f(2.f), vec2f(3.f), vec2f(4.f)).x == Approx(12.f));
+
+        {
+            vec2i out_exp{};
+            REQUIRE(frexp(vec2f(1.7f), &out_exp).x == Approx(0.85f));
+            REQUIRE(out_exp == vec2i(1));
+        }
+
+        REQUIRE(ldexp(vec2f(0.85f), vec2i(1)).x == Approx(1.7f));
+    }
+
+    SECTION("Geometric Functions") {
+        REQUIRE(length(vec2f(10.f,0.f)) == Approx(10.f));
+        REQUIRE(length(vec2f(-10.f,0.f)) == Approx(10.f));
+
+        REQUIRE(distance(vec2f(5.f,0.f), vec2f(10.f,0.f)) == Approx(5.f));
+        REQUIRE(distance(vec2f(-5.f,0.f), vec2f(-10.f,0.f)) == Approx(5.f));
+
+        STATIC_REQUIRE(dot(vec2i(1,2),vec2i(3,4)) == 11);
+        STATIC_REQUIRE(cross(vec3i(1,0,0),vec3i(0,1,0)) == vec3i(0,0,1));
+        REQUIRE(normalize(vec2f(0.5f,0.f)).x == Approx(1.f));
+
+        REQUIRE(faceforward(vec2f(1.f), vec2f(2.f), vec2f(3.f)).x == Approx(-1.f));
+        REQUIRE(reflect(vec2f(1.f), vec2f(2.f)).x == Approx(-15.f));
+        REQUIRE(refract(vec2f(1.f), vec2f(2.f), 1.f).x == Approx(-15.f));
+    }
+
+    SECTION("Vector Relational Functions") {
+        STATIC_REQUIRE(less(vec3i(1,1,1), vec3i(0,1,2)) == vec3b(false, false, true));
+        STATIC_REQUIRE(less_equal(vec3i(1,1,1), vec3i(0,1,2)) == vec3b(false, true, true));
+
+        STATIC_REQUIRE(greater(vec3i(1,1,1), vec3i(0,1,2)) == vec3b(true, false, false));
+        STATIC_REQUIRE(greater_equal(vec3i(1,1,1), vec3i(0,1,2)) == vec3b(true, true, false));
+
+        STATIC_REQUIRE(equal_to(vec3i(1,1,1), vec3i(0,1,2)) == vec3b(false, true, false));
+        STATIC_REQUIRE(not_equal_to(vec3i(1,1,1), vec3i(0,1,2)) == vec3b(true, false, true));
+
+        STATIC_REQUIRE_FALSE(any(vec2b(false, false)));
+        STATIC_REQUIRE(any(vec2b(true, false)));
+        STATIC_REQUIRE(any(vec2b(false, true)));
+        STATIC_REQUIRE(any(vec2b(true, true)));
+
+        STATIC_REQUIRE_FALSE(all(vec2b(false, false)));
+        STATIC_REQUIRE_FALSE(all(vec2b(true, false)));
+        STATIC_REQUIRE_FALSE(all(vec2b(false, true)));
+        STATIC_REQUIRE(all(vec2b(true, true)));
+
+        STATIC_REQUIRE(not_(vec2b(false, false)) == vec2b(true, true));
+        STATIC_REQUIRE(not_(vec2b(true, false)) == vec2b(false, true));
+        STATIC_REQUIRE(not_(vec2b(false, true)) == vec2b(true, false));
+        STATIC_REQUIRE(not_(vec2b(true, true)) == vec2b(false, false));
+    }
+}
diff --git a/untests/vmath_vec_ops_tests.cpp b/untests/vmath_vec_ops_tests.cpp
deleted file mode 100644
index 92bb3c3..0000000
--- a/untests/vmath_vec_ops_tests.cpp
+++ /dev/null
@@ -1,44 +0,0 @@
-/*******************************************************************************
- * This file is part of the "https://github.com/blackmatov/vmath.hpp"
- * For conditions of distribution and use, see copyright notice in LICENSE.md
- * Copyright (C) 2020, by Matvey Cherevko (blackmatov@gmail.com)
- ******************************************************************************/
-
-#include <vmath.hpp/vmath_vec_ops.hpp>
-
-#define CATCH_CONFIG_FAST_COMPILE
-#include <catch2/catch.hpp>
-
-namespace
-{
-}
-
-TEST_CASE("vmath/vec_ops") {
-    using namespace vmath_hpp;
-
-    SECTION("Vector Relational Functions") {
-        STATIC_REQUIRE(less(vec3i(1,1,1), vec3i(0,1,2)) == vec3b(false, false, true));
-        STATIC_REQUIRE(less_equal(vec3i(1,1,1), vec3i(0,1,2)) == vec3b(false, true, true));
-
-        STATIC_REQUIRE(greater(vec3i(1,1,1), vec3i(0,1,2)) == vec3b(true, false, false));
-        STATIC_REQUIRE(greater_equal(vec3i(1,1,1), vec3i(0,1,2)) == vec3b(true, true, false));
-
-        STATIC_REQUIRE(equal(vec3i(1,1,1), vec3i(0,1,2)) == vec3b(false, true, false));
-        STATIC_REQUIRE(not_equal(vec3i(1,1,1), vec3i(0,1,2)) == vec3b(true, false, true));
-
-        STATIC_REQUIRE_FALSE(any(vec2b(false, false)));
-        STATIC_REQUIRE(any(vec2b(true, false)));
-        STATIC_REQUIRE(any(vec2b(false, true)));
-        STATIC_REQUIRE(any(vec2b(true, true)));
-
-        STATIC_REQUIRE_FALSE(all(vec2b(false, false)));
-        STATIC_REQUIRE_FALSE(all(vec2b(true, false)));
-        STATIC_REQUIRE_FALSE(all(vec2b(false, true)));
-        STATIC_REQUIRE(all(vec2b(true, true)));
-
-        STATIC_REQUIRE(not_(vec2b(false, false)) == vec2b(true, true));
-        STATIC_REQUIRE(not_(vec2b(true, false)) == vec2b(false, true));
-        STATIC_REQUIRE(not_(vec2b(false, true)) == vec2b(true, false));
-        STATIC_REQUIRE(not_(vec2b(true, true)) == vec2b(false, false));
-    }
-}