diff --git a/headers/vmath.hpp/vmath_mat.hpp b/headers/vmath.hpp/vmath_mat.hpp
index d9f3972..004cb89 100644
--- a/headers/vmath.hpp/vmath_mat.hpp
+++ b/headers/vmath.hpp/vmath_mat.hpp
@@ -19,108 +19,111 @@ namespace vmath_hpp::detail
     template < typename T >
     class mat_base<T, 2> {
     public:
-        vec<T, 2> rows[2];
+        using row_type = vec<T, 2>;
+        row_type rows[2];
     public:
         constexpr mat_base() : rows{
-            {1, 0},
-            {0, 1},
+            row_type{1, 0},
+            row_type{0, 1},
         } {}
 
         constexpr explicit mat_base(T v)
         : rows{
-            {v, 0},
-            {0, v}} {}
+            row_type{v, 0},
+            row_type{0, v}} {}
 
         constexpr mat_base(
             T m11, T m12,
             T m21, T m22)
         : rows{
-            {m11, m12},
-            {m21, m22}} {}
+            row_type{m11, m12},
+            row_type{m21, m22}} {}
 
         constexpr mat_base(
-            const vec<T, 2>& row0,
-            const vec<T, 2>& row1)
+            const row_type& row0,
+            const row_type& row1)
         : rows{row0, row1} {}
 
         constexpr explicit mat_base(
             const mat_base<T, 3>& other)
         : rows{
-            vec<T, 2>{other.rows[0]},
-            vec<T, 2>{other.rows[1]}} {}
+            row_type{other.rows[0]},
+            row_type{other.rows[1]}} {}
 
         constexpr explicit mat_base(
             const mat_base<T, 4>& other)
         : rows{
-            vec<T, 2>{other.rows[0]},
-            vec<T, 2>{other.rows[1]}} {}
+            row_type{other.rows[0]},
+            row_type{other.rows[1]}} {}
     };
 
     template < typename T >
     class mat_base<T, 3> {
     public:
-        vec<T, 3> rows[3];
+        using row_type = vec<T, 3>;
+        row_type rows[3];
     public:
         constexpr mat_base() : rows{
-            {1, 0, 0},
-            {0, 1, 0},
-            {0, 0, 1},
+            row_type{1, 0, 0},
+            row_type{0, 1, 0},
+            row_type{0, 0, 1},
         } {}
 
         constexpr explicit mat_base(T v)
         : rows{
-            {v, 0, 0},
-            {0, v, 0},
-            {0, 0, v}} {}
+            row_type{v, 0, 0},
+            row_type{0, v, 0},
+            row_type{0, 0, v}} {}
 
         constexpr mat_base(
             T m11, T m12, T m13,
             T m21, T m22, T m23,
             T m31, T m32, T m33)
         : rows{
-            {m11, m12, m13},
-            {m21, m22, m23},
-            {m31, m32, m33}} {}
+            row_type{m11, m12, m13},
+            row_type{m21, m22, m23},
+            row_type{m31, m32, m33}} {}
 
         constexpr mat_base(
-            const vec<T, 3>& row0,
-            const vec<T, 3>& row1,
-            const vec<T, 3>& row2)
+            const row_type& row0,
+            const row_type& row1,
+            const row_type& row2)
         : rows{row0, row1, row2} {}
 
         constexpr explicit mat_base(
             const mat_base<T, 2>& other)
         : rows{
-            {other.rows[0], 0},
-            {other.rows[1], 0},
-            {0, 0, 1}} {}
+            row_type{other.rows[0], 0},
+            row_type{other.rows[1], 0},
+            row_type{0, 0, 1}} {}
 
         constexpr explicit mat_base(
             const mat_base<T, 4>& other)
         : rows{
-            vec<T, 3>{other.rows[0]},
-            vec<T, 3>{other.rows[1]},
-            vec<T, 3>{other.rows[2]}} {}
+            row_type{other.rows[0]},
+            row_type{other.rows[1]},
+            row_type{other.rows[2]}} {}
     };
 
     template < typename T >
     class mat_base<T, 4> {
     public:
-        vec<T, 4> rows[4];
+        using row_type = vec<T, 4>;
+        row_type rows[4];
     public:
         constexpr mat_base() : rows{
-            {1, 0, 0, 0},
-            {0, 1, 0, 0},
-            {0, 0, 1, 0},
-            {0, 0, 0, 1},
+            row_type{1, 0, 0, 0},
+            row_type{0, 1, 0, 0},
+            row_type{0, 0, 1, 0},
+            row_type{0, 0, 0, 1},
         } {}
 
         constexpr explicit mat_base(T v)
         : rows{
-            {v, 0, 0, 0},
-            {0, v, 0, 0},
-            {0, 0, v, 0},
-            {0, 0, 0, v}} {}
+            row_type{v, 0, 0, 0},
+            row_type{0, v, 0, 0},
+            row_type{0, 0, v, 0},
+            row_type{0, 0, 0, v}} {}
 
         constexpr mat_base(
             T m11, T m12, T m13, T m14,
@@ -128,33 +131,33 @@ namespace vmath_hpp::detail
             T m31, T m32, T m33, T m34,
             T m41, T m42, T m43, T m44)
         : rows{
-            {m11, m12, m13, m14},
-            {m21, m22, m23, m24},
-            {m31, m32, m33, m34},
-            {m41, m42, m43, m44}} {}
+            row_type{m11, m12, m13, m14},
+            row_type{m21, m22, m23, m24},
+            row_type{m31, m32, m33, m34},
+            row_type{m41, m42, m43, m44}} {}
 
         constexpr mat_base(
-            const vec<T, 4>& row0,
-            const vec<T, 4>& row1,
-            const vec<T, 4>& row2,
-            const vec<T, 4>& row3)
+            const row_type& row0,
+            const row_type& row1,
+            const row_type& row2,
+            const row_type& row3)
         : rows{row0, row1, row2, row3} {}
 
         constexpr explicit mat_base(
             const mat_base<T, 2>& other)
         : rows{
-            {other.rows[0], 0, 0},
-            {other.rows[1], 0, 0},
-            {0, 0, 1, 0},
-            {0, 0, 0, 1}} {}
+            row_type{other.rows[0], 0, 0},
+            row_type{other.rows[1], 0, 0},
+            row_type{0, 0, 1, 0},
+            row_type{0, 0, 0, 1}} {}
 
         constexpr explicit mat_base(
             const mat_base<T, 3>& other)
         : rows{
-            {other.rows[0], 0},
-            {other.rows[1], 0},
-            {other.rows[2], 0},
-            {0, 0, 0, 1}} {}
+            row_type{other.rows[0], 0},
+            row_type{other.rows[1], 0},
+            row_type{other.rows[2], 0},
+            row_type{0, 0, 0, 1}} {}
     };
 }
 
@@ -166,18 +169,18 @@ namespace vmath_hpp
         using self_type = mat;
         using base_type = detail::mat_base<T, Size>;
     public:
-        using value_type = vec<T, Size>;
+        using row_type = vec<T, Size>;
 
-        using pointer = value_type*;
-        using const_pointer = const value_type*;
+        using pointer = row_type*;
+        using const_pointer = const row_type*;
 
-        using reference = value_type&;
-        using const_reference = const value_type&;
+        using reference = row_type&;
+        using const_reference = const row_type&;
 
         static constexpr std::size_t size = Size;
     public:
-        using base_type::rows;
         using base_type::mat_base;
+        using base_type::rows;
 
         mat() = default;
 
@@ -190,7 +193,7 @@ namespace vmath_hpp
         void swap(mat& other) noexcept(std::is_nothrow_swappable_v<T>) {
             for ( std::size_t i = 0; i < Size; ++i ) {
                 using std::swap;
-                swap(rows[i], other.rows[i]);
+                swap((*this)[i], other[i]);
             }
         }
 
diff --git a/headers/vmath.hpp/vmath_mat_fun.hpp b/headers/vmath.hpp/vmath_mat_fun.hpp
index c194f7a..0f06091 100644
--- a/headers/vmath.hpp/vmath_mat_fun.hpp
+++ b/headers/vmath.hpp/vmath_mat_fun.hpp
@@ -11,34 +11,247 @@
 #include "vmath_fun.hpp"
 #include "vmath_mat.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 mat<A, Size>& a, std::index_sequence<Is...>)
+            -> mat<typename std::invoke_result_t<F, vec<A, Size>>::value_type, Size>
+        {
+            return { f(a[Is])... };
+        }
+
+        template < typename A, typename B, std::size_t Size, typename F, std::size_t... Is >
+        constexpr auto zip_impl(F&& f, const mat<A, Size>& a, const mat<B, Size>& b, std::index_sequence<Is...>)
+            -> mat<typename std::invoke_result_t<F, vec<A, Size>, vec<B, Size>>::value_type, 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 mat<A, Size>& a, const mat<B, Size>& b, const mat<C, Size>& c, std::index_sequence<Is...>)
+            -> mat<typename std::invoke_result_t<F, vec<A, Size>, vec<B, Size>, vec<C, Size>>::value_type, 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 mat<B, Size>& b, std::index_sequence<Is...>)
+            -> A
+        {
+            return ((init = f(std::move(init), b[Is])), ...);
+        }
+
+        template < typename A, typename B, typename C, std::size_t Size, typename F, std::size_t... Is >
+        constexpr auto fold_impl(F&& f, A init, const mat<B, Size>& b, const mat<C, Size>& c, std::index_sequence<Is...>)
+            -> A
+        {
+            return ((init = f(std::move(init), b[Is], c[Is])), ...);
+        }
+    }
+
+    template < typename A, std::size_t Size, typename F >
+    constexpr auto map(F&& f, const mat<A, Size>& a) {
+        return impl::map_impl(std::forward<F>(f), a, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, std::size_t Size, typename F >
+    constexpr auto zip(F&& f, const mat<A, Size>& a, const mat<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 mat<A, Size>& a, const mat<B, Size>& b, const mat<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 mat<B, Size>& b) {
+        return impl::fold_impl(std::forward<F>(f), std::move(init), b, std::make_index_sequence<Size>{});
+    }
+
+    template < typename A, typename B, typename C, std::size_t Size, typename F >
+    constexpr auto fold(F&& f, A init, const mat<B, Size>& b, const mat<C, Size>& c) {
+        return impl::fold_impl(std::forward<F>(f), std::move(init), b, c, std::make_index_sequence<Size>{});
+    }
+}
+
 //
 // Operators
 //
 
 namespace vmath_hpp
 {
+    // -operator
+
     template < typename T, std::size_t Size >
-    constexpr bool operator==(const mat<T, Size>& l, const mat<T, Size>& r) {
-        for ( std::size_t i = 0; i < Size; ++i ) {
-            if ( !(l[i] == r[i]) ) {
-                return false;
-            }
-        }
-        return true;
+    constexpr mat<T, Size> operator-(const mat<T, Size>& xs) {
+        return map(std::negate<>(), xs);
+    }
+
+    // operator+
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size> operator+(const mat<T, Size>& xs, T y) {
+        return map([y](const vec<T, Size>& x){ return x + y; }, xs);
     }
 
     template < typename T, std::size_t Size >
-    constexpr bool operator!=(const mat<T, Size>& l, const mat<T, Size>& r) {
-        return !(l == r);
+    constexpr mat<T, Size> operator+(T x, const mat<T, Size>& ys) {
+        return map([x](const vec<T, Size>& y){ return x + y; }, ys);
     }
 
     template < typename T, std::size_t Size >
-    constexpr bool operator<(const mat<T, Size>& l, const mat<T, Size>& r) {
+    constexpr mat<T, Size> operator+(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return zip(std::plus<>(), xs, ys);
+    }
+
+    // operator-
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size> operator-(const mat<T, Size>& xs, T y) {
+        return map([y](const vec<T, Size>& x){ return x - y; }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size> operator-(T x, const mat<T, Size>& ys) {
+        return map([x](const vec<T, Size>& y){ return x - y; }, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size> operator-(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return zip(std::minus<>(), xs, ys);
+    }
+
+    // operator*
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size> operator*(const mat<T, Size>& xs, T y) {
+        return map([y](const vec<T, Size>& x){ return x * y; }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size> operator*(T x, const mat<T, Size>& ys) {
+        return map([x](const vec<T, Size>& y){ return x * y; }, ys);
+    }
+
+    template < typename T >
+    constexpr vec<T, 2> operator*(const vec<T, 2>& xs, const mat<T, 2>& ys) {
+        return {
+            xs.x * ys[0][0] + xs.y * ys[1][0],
+            xs.x * ys[0][1] + xs.y * ys[1][1]};
+    }
+
+    template < typename T >
+    constexpr mat<T, 2> operator*(const mat<T, 2>& xs, const mat<T, 2>& ys) {
+        return {
+            xs[0][0] * ys[0][0] + xs[0][1] * ys[1][0],
+            xs[0][0] * ys[0][1] + xs[0][1] * ys[1][1],
+
+            xs[1][0] * ys[0][0] + xs[1][1] * ys[1][0],
+            xs[1][0] * ys[0][1] + xs[1][1] * ys[1][1]};
+    }
+
+    template < typename T >
+    constexpr vec<T, 3> operator*(const vec<T, 3>& xs, const mat<T, 3>& ys) {
+        return {
+            xs.x * ys[0][0] + xs.y * ys[1][0] + xs.z * ys[2][0],
+            xs.x * ys[0][1] + xs.y * ys[1][1] + xs.z * ys[2][1],
+            xs.x * ys[0][2] + xs.y * ys[1][2] + xs.z * ys[2][2]};
+    }
+
+    template < typename T >
+    constexpr mat<T, 3> operator*(const mat<T, 3>& xs, const mat<T, 3>& ys) {
+        return {
+            xs[0][0] * ys[0][0] + xs[0][1] * ys[1][0] + xs[0][2] * ys[2][0],
+            xs[0][0] * ys[0][1] + xs[0][1] * ys[1][1] + xs[0][2] * ys[2][1],
+            xs[0][0] * ys[0][2] + xs[0][1] * ys[1][2] + xs[0][2] * ys[2][2],
+
+            xs[1][0] * ys[0][0] + xs[1][1] * ys[1][0] + xs[1][2] * ys[2][0],
+            xs[1][0] * ys[0][1] + xs[1][1] * ys[1][1] + xs[1][2] * ys[2][1],
+            xs[1][0] * ys[0][2] + xs[1][1] * ys[1][2] + xs[1][2] * ys[2][2],
+
+            xs[2][0] * ys[0][0] + xs[2][1] * ys[1][0] + xs[2][2] * ys[2][0],
+            xs[2][0] * ys[0][1] + xs[2][1] * ys[1][1] + xs[2][2] * ys[2][1],
+            xs[2][0] * ys[0][2] + xs[2][1] * ys[1][2] + xs[2][2] * ys[2][2]};
+    }
+
+    template < typename T >
+    constexpr vec<T, 4> operator*(const vec<T, 4>& xs, const mat<T, 4>& ys) {
+        return {
+            xs.x * ys[0][0] + xs.y * ys[1][0] + xs.z * ys[2][0] + xs.w * ys[3][0],
+            xs.x * ys[0][1] + xs.y * ys[1][1] + xs.z * ys[2][1] + xs.w * ys[3][1],
+            xs.x * ys[0][2] + xs.y * ys[1][2] + xs.z * ys[2][2] + xs.w * ys[3][2],
+            xs.x * ys[0][3] + xs.y * ys[1][3] + xs.z * ys[2][3] + xs.w * ys[3][3]};
+    }
+
+    template < typename T >
+    constexpr mat<T, 4> operator*(const mat<T, 4>& xs, const mat<T, 4>& ys) {
+        return {
+            xs[0][0] * ys[0][0] + xs[0][1] * ys[1][0] + xs[0][2] * ys[2][0] + xs[0][3] * ys[3][0],
+            xs[0][0] * ys[0][1] + xs[0][1] * ys[1][1] + xs[0][2] * ys[2][1] + xs[0][3] * ys[3][1],
+            xs[0][0] * ys[0][2] + xs[0][1] * ys[1][2] + xs[0][2] * ys[2][2] + xs[0][3] * ys[3][2],
+            xs[0][0] * ys[0][3] + xs[0][1] * ys[1][3] + xs[0][2] * ys[2][3] + xs[0][3] * ys[3][3],
+
+            xs[1][0] * ys[0][0] + xs[1][1] * ys[1][0] + xs[1][2] * ys[2][0] + xs[1][3] * ys[3][0],
+            xs[1][0] * ys[0][1] + xs[1][1] * ys[1][1] + xs[1][2] * ys[2][1] + xs[1][3] * ys[3][1],
+            xs[1][0] * ys[0][2] + xs[1][1] * ys[1][2] + xs[1][2] * ys[2][2] + xs[1][3] * ys[3][2],
+            xs[1][0] * ys[0][3] + xs[1][1] * ys[1][3] + xs[1][2] * ys[2][3] + xs[1][3] * ys[3][3],
+
+            xs[2][0] * ys[0][0] + xs[2][1] * ys[1][0] + xs[2][2] * ys[2][0] + xs[2][3] * ys[3][0],
+            xs[2][0] * ys[0][1] + xs[2][1] * ys[1][1] + xs[2][2] * ys[2][1] + xs[2][3] * ys[3][1],
+            xs[2][0] * ys[0][2] + xs[2][1] * ys[1][2] + xs[2][2] * ys[2][2] + xs[2][3] * ys[3][2],
+            xs[2][0] * ys[0][3] + xs[2][1] * ys[1][3] + xs[2][2] * ys[2][3] + xs[2][3] * ys[3][3],
+
+            xs[3][0] * ys[0][0] + xs[3][1] * ys[1][0] + xs[3][2] * ys[2][0] + xs[3][3] * ys[3][0],
+            xs[3][0] * ys[0][1] + xs[3][1] * ys[1][1] + xs[3][2] * ys[2][1] + xs[3][3] * ys[3][1],
+            xs[3][0] * ys[0][2] + xs[3][1] * ys[1][2] + xs[3][2] * ys[2][2] + xs[3][3] * ys[3][2],
+            xs[3][0] * ys[0][3] + xs[3][1] * ys[1][3] + xs[3][2] * ys[2][3] + xs[3][3] * ys[3][3]};
+    }
+
+    // operator/
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size> operator/(const mat<T, Size>& xs, T y) {
+        return map([y](const vec<T, Size>& x){ return x / y; }, xs);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size> operator/(T x, const mat<T, Size>& ys) {
+        return map([x](const vec<T, Size>& y){ return x / y; }, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr mat<T, Size> operator/(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return zip(std::divides<>(), xs, ys);
+    }
+
+    // operator==
+
+    template < typename T, std::size_t Size >
+    constexpr bool operator==(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return fold([](bool acc, const vec<T, Size>& x, const vec<T, Size>& y){
+            return acc && (x == y);
+        }, true, xs, ys);
+    }
+
+    template < typename T, std::size_t Size >
+    constexpr bool operator!=(const mat<T, Size>& xs, const mat<T, Size>& ys) {
+        return fold([](bool acc, const vec<T, Size>& x, const vec<T, Size>& y){
+            return acc || (x != y);
+        }, false, xs, ys);
+    }
+
+    // operator<
+
+    template < typename T, std::size_t Size >
+    constexpr bool operator<(const mat<T, Size>& xs, const mat<T, Size>& ys) {
         for ( std::size_t i = 0; i < Size; ++i ) {
-            if ( l[i] < r[i] ) {
+            if ( xs[i] < ys[i] ) {
                 return true;
             }
-            if ( r[i] < l[i] ) {
+            if ( ys[i] < xs[i] ) {
                 return false;
             }
         }
diff --git a/untests/vmath_mat_fun_tests.cpp b/untests/vmath_mat_fun_tests.cpp
index 603955d..ad5ff2c 100644
--- a/untests/vmath_mat_fun_tests.cpp
+++ b/untests/vmath_mat_fun_tests.cpp
@@ -11,8 +11,32 @@
 
 namespace
 {
+    using namespace vmath_hpp;
 }
 
 TEST_CASE("vmath/mat_fun") {
-    using namespace vmath_hpp;
+    SECTION("Operators") {
+        STATIC_REQUIRE(-mat2i(1,2,3,4) == mat2i(-1,-2,-3,-4));
+
+        STATIC_REQUIRE(mat2i(1,2,3,4) + 2 == mat2i(3,4,5,6));
+        STATIC_REQUIRE(mat2i(1,2,3,4) - 2 == mat2i(-1,0,1,2));
+        STATIC_REQUIRE(mat2i(1,2,3,4) * 2 == mat2i(2,4,6,8));
+        STATIC_REQUIRE(mat2i(1,2,3,4) / 2 == mat2i(0,1,1,2));
+
+        STATIC_REQUIRE(4 + mat2i(1,2,3,4) == mat2i(5,6,7,8));
+        STATIC_REQUIRE(4 - mat2i(1,2,3,4) == mat2i(3,2,1,0));
+        STATIC_REQUIRE(4 * mat2i(1,2,3,4) == mat2i(4,8,12,16));
+        STATIC_REQUIRE(4 / mat2i(1,2,3,4) == mat2i(4,2,1,1));
+
+        STATIC_REQUIRE(mat2i(1,2,3,4) + mat2i(5,6,7,8) == mat2i(6,8,10,12));
+        STATIC_REQUIRE(mat2i(1,2,3,4) - mat2i(5,6,7,8) == mat2i(-4,-4,-4,-4));
+        STATIC_REQUIRE(mat2i(5,6,7,8) / mat2i(1,2,3,4) == mat2i(5,3,2,2));
+
+        STATIC_REQUIRE(mat2i() * mat2i() == mat2i());
+        STATIC_REQUIRE(mat3i() * mat3i() == mat3i());
+
+        STATIC_REQUIRE(vec2i(1,2) * mat2i() == vec2i(1,2));
+        STATIC_REQUIRE(vec3i(1,2,3) * mat3i() == vec3i(1,2,3));
+        STATIC_REQUIRE(vec4i(1,2,3,4) * mat4i() == vec4i(1,2,3,4));
+    }
 }