diff --git a/libs/math/include/psemek/math/detail/array.hpp b/libs/math/include/psemek/math/detail/array.hpp
index a86f9f13..fd835628 100644
--- a/libs/math/include/psemek/math/detail/array.hpp
+++ b/libs/math/include/psemek/math/detail/array.hpp
@@ -2,7 +2,6 @@
 
 #include <psemek/math/dynamic.hpp>
 #include <psemek/util/exception.hpp>
-#include <psemek/util/array.hpp>
 
 namespace psemek::math::detail
 {
@@ -15,10 +14,42 @@ namespace psemek::math::detail
 		{}
 	};
 
+	template <std::size_t N>
+	struct dynamic_dimensions
+	{};
+
+	template <>
+	struct dynamic_dimensions<dynamic>
+	{
+		std::size_t size = 0;
+	};
+
+	template <std::size_t N>
+	void check_dynamic_size(dynamic_dimensions<N> d1, dynamic_dimensions<N> d2)
+	{
+		if constexpr (N == dynamic)
+		{
+			if (d1.size != d2.size)
+				throw dynamic_size_mismatch(d1.size, d2.size);
+		}
+	}
+
 	template <typename T, std::size_t N>
 	struct array
 	{
-		using type = T[N];
+		struct type
+		{
+			static constexpr std::size_t size = N;
+
+			T data[N];
+
+			type(dynamic_dimensions<N>){}
+
+			dynamic_dimensions<N> dimensions() const { return {}; }
+
+			T const & operator[](std::size_t i) const { return data[i]; }
+			T & operator[](std::size_t i) { return data[i]; }
+		};
 	};
 
 	template <typename T>
@@ -26,17 +57,35 @@ namespace psemek::math::detail
 	{
 		struct type
 		{
+			static constexpr std::size_t size = 0;
+
+			type(dynamic_dimensions<0>){}
+
+			dynamic_dimensions<0> dimensions() const { return {}; }
+
 			T const & operator[](std::size_t) const { throw empty_array_exception{}; }
 			T & operator[](std::size_t) { throw empty_array_exception{}; }
-
-			// type operator + (std::size_t) const { return *this; }
 		};
 	};
 
 	template <typename T>
 	struct array<T, dynamic>
 	{
-		using type = util::array<T>;
+		struct type
+		{
+			std::size_t size;
+			std::unique_ptr<T[]> data;
+
+			type(dynamic_dimensions<dynamic> dimensions)
+				: size(dimensions.size)
+				, data(std::make_unique_for_overwrite<T[]>(dimensions.size))
+			{}
+
+			dynamic_dimensions<dynamic> dimensions() const { return {.size = size}; }
+
+			T const & operator[](std::size_t i) const { return data[i]; }
+			T & operator[](std::size_t i) { return data[i]; }
+		};
 	};
 
 	template <typename T, typename ... Args>
diff --git a/libs/math/include/psemek/math/detail/array_2d.hpp b/libs/math/include/psemek/math/detail/array_2d.hpp
index bd257e33..98f0075f 100644
--- a/libs/math/include/psemek/math/detail/array_2d.hpp
+++ b/libs/math/include/psemek/math/detail/array_2d.hpp
@@ -5,10 +5,37 @@
 namespace psemek::math::detail
 {
 
+	template <std::size_t R, std::size_t C>
+	struct dynamic_dimensions_2d
+	{
+		dynamic_dimensions<R> rows = {};
+		dynamic_dimensions<C> columns = {};
+	};
+
+	template <std::size_t R, std::size_t C>
+	void check_dynamic_size(dynamic_dimensions_2d<R, C> d1, dynamic_dimensions_2d<R, C> d2)
+	{
+		check_dynamic_size(d1.rows, d2.rows);
+		check_dynamic_size(d1.columns, d2.columns);
+	}
+
 	template <typename T, std::size_t R, std::size_t C>
 	struct array_2d
 	{
-		using type = T[R][C];
+		struct type
+		{
+			T data[R * C];
+
+			static constexpr std::size_t rows = R;
+			static constexpr std::size_t columns = C;
+
+			type(dynamic_dimensions_2d<R, C>) {}
+
+			dynamic_dimensions_2d<R, C> dimensions() const { return {}; }
+
+			T const * operator[](std::size_t row) const { return data + C * row; }
+			T * operator[](std::size_t row) { return data + C * row; }
+		};
 	};
 
 	template <typename T, std::size_t R>
@@ -19,6 +46,10 @@ namespace psemek::math::detail
 			static constexpr std::size_t rows = R;
 			static constexpr std::size_t columns = 0;
 
+			type(dynamic_dimensions_2d<R, 0>) {}
+
+			dynamic_dimensions_2d<R, 0> dimensions() const { return {}; }
+
 			T const * operator[](std::size_t) const { throw empty_array_exception{}; }
 			T * operator[](std::size_t) { throw empty_array_exception{}; }
 		};
@@ -32,10 +63,12 @@ namespace psemek::math::detail
 			std::size_t rows;
 			static constexpr std::size_t columns = 0;
 
-			type(std::size_t rows)
-				: rows(rows)
+			type(dynamic_dimensions_2d<dynamic, 0> dimensions)
+				: rows(dimensions.rows.size)
 			{}
 
+			dynamic_dimensions_2d<dynamic, 0> dimensions() const { return {.rows = rows}; }
+
 			T const * operator[](std::size_t) const { throw empty_array_exception{}; }
 			T * operator[](std::size_t) { throw empty_array_exception{}; }
 		};
@@ -49,6 +82,10 @@ namespace psemek::math::detail
 			static constexpr std::size_t rows = 0;
 			static constexpr std::size_t columns = C;
 
+			type(dynamic_dimensions_2d<0, C>) {}
+
+			dynamic_dimensions_2d<0, C> dimensions() const { return {}; }
+
 			T const * operator[](std::size_t) const { throw empty_array_exception{}; }
 			T * operator[](std::size_t) { throw empty_array_exception{}; }
 		};
@@ -62,10 +99,12 @@ namespace psemek::math::detail
 			static constexpr std::size_t rows = 0;
 			std::size_t columns;
 
-			type(std::size_t columns)
-				: columns(columns)
+			type(dynamic_dimensions_2d<0, dynamic> dimensions)
+				: columns(dimensions.columns.size)
 			{}
 
+			dynamic_dimensions_2d<0, dynamic> dimensions() const { return {.columns = columns}; }
+
 			T const * operator[](std::size_t) const { throw empty_array_exception{}; }
 			T * operator[](std::size_t) { throw empty_array_exception{}; }
 		};
@@ -79,6 +118,10 @@ namespace psemek::math::detail
 			static constexpr std::size_t rows = 0;
 			static constexpr std::size_t columns = 0;
 
+			type(dynamic_dimensions_2d<0, 0>) {}
+
+			dynamic_dimensions_2d<0, 0> dimensions() const { return {}; }
+
 			T const * operator[](std::size_t) const { throw empty_array_exception{}; }
 			T * operator[](std::size_t) { throw empty_array_exception{}; }
 		};
@@ -93,20 +136,12 @@ namespace psemek::math::detail
 			std::size_t columns;
 			std::unique_ptr<T[]> data;
 
-			type()
-				: columns(0)
-			{}
-
-			type(std::size_t columns)
-				: columns(columns)
+			type(dynamic_dimensions_2d<R, dynamic> dimensions)
+				: columns(dimensions.columns.size)
 				, data(std::make_unique_for_overwrite<T[]>(rows * columns))
 			{}
 
-			type([[maybe_unused]] std::size_t rows, std::size_t columns)
-				: type(columns)
-			{
-				assert(rows == R);
-			}
+			dynamic_dimensions_2d<R, dynamic> dimensions() const { return {.columns = columns}; }
 
 			T * operator[] (std::size_t row)
 			{
@@ -138,20 +173,12 @@ namespace psemek::math::detail
 			static constexpr std::size_t columns = C;
 			std::unique_ptr<T[]> data;
 
-			type()
-				: rows(0)
-			{}
-
-			type(std::size_t rows)
-				: rows(rows)
+			type(dynamic_dimensions_2d<dynamic, C> dimensions)
+				: rows(dimensions.rows.size)
 				, data(std::make_unique_for_overwrite<T[]>(rows * columns))
 			{}
 
-			type(std::size_t rows, [[maybe_unused]] std::size_t columns)
-				: type(rows)
-			{
-				assert(columns == C);
-			}
+			dynamic_dimensions_2d<dynamic, C> dimensions() const { return {.rows = rows}; }
 
 			T * operator[] (std::size_t row)
 			{
@@ -183,17 +210,14 @@ namespace psemek::math::detail
 			std::size_t columns;
 			std::unique_ptr<T[]> data;
 
-			type()
-				: rows(0)
-				, columns(0)
-			{}
-
-			type(std::size_t rows, std::size_t columns)
-				: rows(rows)
-				, columns(columns)
+			type(dynamic_dimensions_2d<dynamic, dynamic> dimensions)
+				: rows(dimensions.rows.size)
+				, columns(dimensions.columns.size)
 				, data(std::make_unique_for_overwrite<T[]>(rows * columns))
 			{}
 
+			dynamic_dimensions_2d<dynamic, dynamic> dimensions() const { return {.rows = rows, .columns = columns}; }
+
 			T * operator[] (std::size_t row)
 			{
 				return data.get() + row * columns;
diff --git a/libs/math/include/psemek/math/dynamic.hpp b/libs/math/include/psemek/math/dynamic.hpp
index 1a18b928..37f17942 100644
--- a/libs/math/include/psemek/math/dynamic.hpp
+++ b/libs/math/include/psemek/math/dynamic.hpp
@@ -18,10 +18,4 @@ namespace psemek::math
 		{}
 	};
 
-	inline void check_dynamic_size(std::size_t size1, std::size_t size2)
-	{
-		if (size1 != size2)
-			throw dynamic_size_mismatch(size1, size2);
-	}
-
 }
diff --git a/libs/math/include/psemek/math/matrix.hpp b/libs/math/include/psemek/math/matrix.hpp
index 15d9d032..21662f82 100644
--- a/libs/math/include/psemek/math/matrix.hpp
+++ b/libs/math/include/psemek/math/matrix.hpp
@@ -19,30 +19,30 @@ namespace psemek::math
 		static constexpr std::size_t static_rows = R;
 		static constexpr std::size_t static_columns = C;
 
+		using dynamic_dimensions_type = detail::dynamic_dimensions_2d<R, C>;
+
 		typename detail::array_2d<T, R, C>::type coords;
 
+		// Internal
+		matrix(dynamic_dimensions_type dimensions);
+
+		matrix();
+
+		template <typename ... Args>
+		matrix(Args && ... args) requires (R != dynamic && C != dynamic && sizeof...(Args) == R * C && detail::all_convertible_to<T, Args...>::value);
+
+		explicit matrix(std::size_t rows) requires(R == dynamic && C != dynamic);
+		explicit matrix(std::size_t columns) requires(R != dynamic && C == dynamic);
+		matrix(std::size_t rows, std::size_t columns) requires(R == dynamic && C == dynamic);
+
 		std::size_t rows() const
 		{
-			if constexpr (R == dynamic)
-			{
-				return coords.rows;
-			}
-			else
-			{
-				return R;
-			}
+			return coords.rows;
 		}
 
 		std::size_t columns() const
 		{
-			if constexpr (C == dynamic)
-			{
-				return coords.columns;
-			}
-			else
-			{
-				return C;
-			}
+			return coords.columns;
 		}
 
 		auto operator[](std::size_t i)
@@ -85,32 +85,72 @@ namespace psemek::math
 		matrix & operator += (matrix const & v);
 		matrix & operator -= (matrix const & v);
 
-		static matrix zero();
-		static matrix identity();
-		static matrix scalar(T const & s);
+		static matrix zero(dynamic_dimensions_type dimensions = {});
+		static matrix identity(dynamic_dimensions_type dimensions = {});
+		static matrix scalar(T const & s, dynamic_dimensions_type dimensions = {});
 		static matrix diagonal(vector<T, std::min(R, C)> const & d);
 	};
 
 	template <typename T, std::size_t R, std::size_t C>
-	matrix<T, R, C> matrix<T, R, C>::zero()
+	void check_dynamic_size(matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
 	{
-		matrix<T, R, C> m;
+		check_dynamic_size(m1.coords.dimensions(), m2.coords.dimensions());
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C>::matrix(dynamic_dimensions_type dimensions)
+		: coords(dimensions)
+	{}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C>::matrix()
+		: matrix(dynamic_dimensions_type{})
+	{}
+
+	template <typename T, std::size_t R, std::size_t C>
+	template <typename ... Args>
+	matrix<T, R, C>::matrix(Args && ... args) requires (R != dynamic && C != dynamic && sizeof...(Args) == R * C && detail::all_convertible_to<T, Args...>::value)
+		: coords({})
+	{
+		auto out = values().begin();
+		((*out++ = args), ...);
+	}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C>::matrix(std::size_t rows) requires(R == dynamic && C != dynamic)
+		: coords({.rows = rows})
+	{}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C>::matrix(std::size_t columns) requires(R != dynamic && C == dynamic)
+		: coords({.columns = columns})
+	{}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C>::matrix(std::size_t rows, std::size_t columns) requires(R == dynamic && C == dynamic)
+		: coords({.rows = rows, .columns = columns})
+	{}
+
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C> matrix<T, R, C>::zero(dynamic_dimensions_type dimensions)
+	{
+		matrix<T, R, C> m(dimensions);
 		for (auto & v : m.values())
 			v = T(0);
 		return m;
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
-	matrix<T, R, C> matrix<T, R, C>::identity()
+	matrix<T, R, C> matrix<T, R, C>::identity(dynamic_dimensions_type dimensions)
 	{
-		return scalar(T(1));
+		return scalar(T(1), dimensions);
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
-	matrix<T, R, C> matrix<T, R, C>::scalar(T const & s)
+	matrix<T, R, C> matrix<T, R, C>::scalar(T const & s, dynamic_dimensions_type dimensions)
 	{
-		matrix<T, R, C> m = zero();
-		for (std::size_t i = 0; i < std::min(R, C); ++i)
+		matrix<T, R, C> m = zero(dimensions);
+		for (std::size_t i = 0; i < std::min(m.rows(), m.columns()); ++i)
 			m[i][i] = s;
 		return m;
 	}
@@ -118,7 +158,17 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	matrix<T, R, C> matrix<T, R, C>::diagonal(vector<T, std::min(R, C)> const & d)
 	{
-		matrix<T, R, C> m = zero();
+		dynamic_dimensions_type dimensions;
+		if constexpr (R == dynamic)
+		{
+			dimensions.rows = d.dimension();
+		}
+		if constexpr (C == dynamic)
+		{
+			dimensions.columns = d.dimension();
+		}
+
+		matrix<T, R, C> m = zero(dimensions);
 		for (std::size_t i = 0; i < std::min(R, C); ++i)
 			m[i][i] = d[i];
 		return m;
@@ -127,25 +177,33 @@ namespace psemek::math
 	template <typename T1, typename T, std::size_t R, std::size_t C>
 	matrix<T1, R, C> cast(matrix<T, R, C> const & m)
 	{
-		matrix<T1, R, C> r;
-		for (std::size_t i = 0; i < R; ++i)
-			for (std::size_t j = 0; j < C; ++j)
-				r[i][j] = static_cast<T1>(m[i][j]);
+		matrix<T1, R, C> r(m.coords.dimensions());
+
+		auto out = r.values().begin();
+		for (auto const & value : m.values())
+			*out++ = static_cast<T1>(value);
+
 		return r;
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
 	std::strong_ordering operator <=> (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
 	{
-		for (std::size_t i = 0; i < R; ++i)
+		check_dynamic_size(m1, m2);
+		auto begin1 = m1.values().begin();
+		auto begin2 = m2.values().begin();
+		auto end1 = m1.values().end();
+
+		// Intentionally don't use std::lexicographical_compare_three_way to fake strong ordering
+		// (actual ordering might be partial & weak due to NaN's and signed zeros
+		while (begin1 != end1)
 		{
-			for (std::size_t j = 0; j < C; ++j)
-			{
-				if (m1[i][j] < m2[i][j])
-					return std::strong_ordering::less;
-				else if (m1[i][j] > m2[i][j])
-					return std::strong_ordering::greater;
-			}
+			if (*begin1 < *begin2)
+				return std::strong_ordering::less;
+			else if (*begin1 > *begin2)
+				return std::strong_ordering::greater;
+			++begin1;
+			++begin2;
 		}
 		return std::strong_ordering::equal;
 	}
@@ -189,30 +247,30 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	matrix<T, R, C> operator * (matrix<T, R, C> const & m, T const & s)
 	{
-		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R; ++i)
-			for (std::size_t j = 0; j < C; ++j)
-				r[i][j] = m[i][j] * s;
+		matrix<T, R, C> r(m.coords.dimensions());
+		auto out = r.values().begin();
+		for (auto const & value : m.values())
+			*out++ = value * s;
 		return r;
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
 	matrix<T, R, C> operator * (T const & s, matrix<T, R, C> const & m)
 	{
-		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R; ++i)
-			for (std::size_t j = 0; j < C; ++j)
-				r[i][j] = m[i][j] * s;
+		matrix<T, R, C> r(m.coords.dimensions());
+		auto out = r.values().begin();
+		for (auto const & value : m.values())
+			*out++ = s * value;
 		return r;
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
 	matrix<T, R, C> operator / (matrix<T, R, C> const & m, T const & s)
 	{
-		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R; ++i)
-			for (std::size_t j = 0; j < C; ++j)
-				r[i][j] = m[i][j] / s;
+		matrix<T, R, C> r(m.coords.dimensions());
+		auto out = r.values().begin();
+		for (auto const & value : m.values())
+			*out++ = value / s;
 		return r;
 	}
 
@@ -233,30 +291,34 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	matrix<T, R, C> operator - (matrix<T, R, C> const & m)
 	{
-		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R; ++i)
-			for (std::size_t j = 0; j < C; ++j)
-				r[i][j] = -m[i][j];
+		matrix<T, R, C> r(m.coords.dimensions());
+		auto out = r.values().begin();
+		for (auto const & value : m.values())
+			*out++ = -value;
 		return r;
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
 	matrix<T, R, C> operator + (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
 	{
-		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R; ++i)
-			for (std::size_t j = 0; j < C; ++j)
-				r[i][j] = m1[i][j] + m2[i][j];
+		check_dynamic_size(m1, m2);
+		matrix<T, R, C> r(m1.coords.dimensions());
+		auto in1 = m1.values().begin();
+		auto in2 = m2.values().begin();
+		for (auto & value : r.values())
+			value = (*in1++) + (*in2++);
 		return r;
 	}
 
 	template <typename T, std::size_t R, std::size_t C>
 	matrix<T, R, C> operator - (matrix<T, R, C> const & m1, matrix<T, R, C> const & m2)
 	{
-		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R; ++i)
-			for (std::size_t j = 0; j < C; ++j)
-				r[i][j] = m1[i][j] - m2[i][j];
+		check_dynamic_size(m1, m2);
+		matrix<T, R, C> r(m1.coords.dimensions());
+		auto in1 = m1.values().begin();
+		auto in2 = m2.values().begin();
+		for (auto & value : r.values())
+			value = (*in1++) - (*in2++);
 		return r;
 	}
 
@@ -277,11 +339,11 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	vector<T, R> operator * (matrix<T, R, C> const & m, vector<T, C> const & v)
 	{
-		vector<T, R> r;
-		for (std::size_t i = 0; i < R; ++i)
+		detail::check_dynamic_size(m.coords.dimensions().columns, v.coords.dimensions());
+		auto r = vector<T, R>::zero(m.coords.dimensions().rows);
+		for (std::size_t i = 0; i < m.rows(); ++i)
 		{
-			r[i] = T{};
-			for (std::size_t j = 0; j < C; ++j)
+			for (std::size_t j = 0; j < m.columns(); ++j)
 				r[i] += m[i][j] * v[j];
 		}
 		return r;
@@ -290,12 +352,12 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	vector<T, C> operator * (vector<T, R> const & v, matrix<T, R, C> const & m)
 	{
-		vector<T, C> r;
-		for (std::size_t j = 0; j < C; ++j)
+		detail::check_dynamic_size(v.coords.dimensions(), m.coords.dimensions().rows);
+		auto r = vector<T, C>::zero(m.coords.dimensions().columns);
+		for (std::size_t i = 0; i < m.rows(); ++i)
 		{
-			r[j] = T{};
-			for (std::size_t i = 0; i < R; ++i)
-				r[i] += v[j] * m[i][j];
+			for (std::size_t j = 0; j < m.columns(); ++j)
+				r[j] += v[i] * m[i][j];
 		}
 		return r;
 	}
@@ -303,14 +365,13 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t K, std::size_t C>
 	matrix<T, R, C> operator * (matrix<T, R, K> const & m1, matrix<T, K, C> const & m2)
 	{
-		matrix<T, R, C> r;
-		for (std::size_t i = 0; i < R; ++i)
+		detail::check_dynamic_size(m1.coords.dimensions().columns, m2.coords.dimensions().rows);
+		auto r = matrix<T, R, C>::zero();
+		for (std::size_t i = 0; i < m1.rows(); ++i)
 		{
-			for (std::size_t j = 0; j < C; ++j)
+			for (std::size_t j = 0; j < m2.columns(); ++j)
 			{
-				r[i][j] = T{};
-
-				for (std::size_t k = 0; k < K; ++k)
+				for (std::size_t k = 0; k < m1.columns(); ++k)
 					r[i][j] += m1[i][k] * m2[k][j];
 			}
 		}
@@ -320,9 +381,9 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	matrix<T, C, R> transpose(matrix<T, R, C> const & m)
 	{
-		matrix<T, C, R> r;
-		for (std::size_t i = 0; i < R; ++i)
-			for (std::size_t j = 0; j < C; ++j)
+		matrix<T, C, R> r({.rows = m.coords.dimensions().columns, .columns = m.coords.dimensions().rows});
+		for (std::size_t i = 0; i < m.rows(); ++i)
+			for (std::size_t j = 0; j < m.columns(); ++j)
 				r[j][i] = m[i][j];
 		return r;
 	}
@@ -330,9 +391,11 @@ namespace psemek::math
 	template <typename T, std::size_t N, typename ... Rows>
 	auto by_rows(vector<T, N> const & v, Rows const & ... rows)
 	{
+		(check_dynamic_size(v, rows), ...);
+
 		vector<T, N> m[] = {v, rows...};
 
-		matrix<T, sizeof...(Rows) + 1, N> r;
+		matrix<T, sizeof...(Rows) + 1, N> r({.columns = v.coords.dimensions()});
 		for (std::size_t i = 0; i < r.rows(); ++i)
 			for (std::size_t j = 0; j < r.columns(); ++j)
 				r[i][j] = m[i][j];
@@ -342,9 +405,10 @@ namespace psemek::math
 	template <typename T, std::size_t N, typename ... Columns>
 	auto by_columns(vector<T, N> const & v, Columns const & ... columns)
 	{
+		(check_dynamic_size(v, columns), ...);
 		vector<T, N> m[] = {v, columns...};
 
-		matrix<T, N, sizeof...(Columns) + 1> r;
+		matrix<T, N, sizeof...(Columns) + 1> r({.rows = v.coords.dimensions()});
 		for (std::size_t i = 0; i < r.rows(); ++i)
 			for (std::size_t j = 0; j < r.columns(); ++j)
 				r[i][j] = m[j][i];
@@ -354,8 +418,8 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	vector<T, C> row(matrix<T, R, C> const & m, std::size_t i)
 	{
-		vector<T, C> r;
-		for (std::size_t j = 0; j < C; ++j)
+		vector<T, C> r(m.coords.dimensions().columns);
+		for (std::size_t j = 0; j < m.columns(); ++j)
 			r[j] = m[i][j];
 		return r;
 	}
@@ -363,8 +427,8 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	vector<T, R> column(matrix<T, R, C> const & m, std::size_t j)
 	{
-		vector<T, R> r;
-		for (std::size_t i = 0; i < R; ++i)
+		vector<T, R> r(m.coords.dimensions().rows);
+		for (std::size_t i = 0; i < m.rows(); ++i)
 			r[i] = m[i][j];
 		return r;
 	}
@@ -373,8 +437,8 @@ namespace psemek::math
 	T frobenius_norm_sqr(matrix<T, R, C> const & m)
 	{
 		T r{0};
-		for (std::size_t i = 0; i < R; ++i)
-			for (std::size_t j = 0; j < C; ++j)
+		for (std::size_t i = 0; i < m.rows(); ++i)
+			for (std::size_t j = 0; j < m.columns(); ++j)
 				r += sqr(m[i][j]);
 		return r;
 	}
@@ -389,12 +453,13 @@ namespace psemek::math
 	T trace(matrix<T, R, C> const & m)
 	{
 		T r{0};
-		for (std::size_t i = 0; i < std::min(R, C); ++i)
+		for (std::size_t i = 0; i < std::min(m.rows(), m.columns()); ++i)
 			r += m[i][i];
 		return r;
 	}
 
 	template <std::size_t R1, std::size_t C1, std::size_t R2, std::size_t C2, typename T, std::size_t R, std::size_t C>
+	requires (R1 != dynamic && C1 != dynamic && R2 != dynamic && C2 != dynamic)
 	matrix<T, R2 - R1, C2 - C1> submatrix(matrix<T, R, C> const & m)
 	{
 		static_assert(R1 < R2);
@@ -403,13 +468,9 @@ namespace psemek::math
 		static_assert(C2 <= C);
 
 		matrix<T, R2 - R1, C2 - C1> result;
-		for (std::size_t i = 0; i < R2 - R1; ++i)
-		{
-			for (std::size_t j = 0; j < C2 - C1; ++j)
-			{
+		for (std::size_t i = 0; i < result.rows(); ++i)
+			for (std::size_t j = 0; j < result.columns(); ++j)
 				result[i][j] = m[i + R1][j + C1];
-			}
-		}
 		return result;
 	}
 
@@ -426,12 +487,12 @@ namespace psemek::math
 		return r;
 	}
 
-	template <typename T, std::size_t N>
-	matrix<T, N, N> outer_product(vector<T, N> const & v1, vector<T, N> const & v2)
+	template <typename T, std::size_t R, std::size_t C>
+	matrix<T, R, C> outer_product(vector<T, R> const & v1, vector<T, C> const & v2)
 	{
-		matrix<T, N, N> r;
-		for (std::size_t i = 0; i < N; ++i)
-			for (std::size_t j = 0; j < N; ++j)
+		matrix<T, R, C> r({.rows = v1.coords.dimensions(), .columns = v2.coords.dimensions()});
+		for (std::size_t i = 0; i < r.rows(); ++i)
+			for (std::size_t j = 0; j < r.columns(); ++j)
 				r[i][j] = v1[i] * v2[j];
 		return r;
 	}
@@ -439,9 +500,9 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	std::ostream & operator << (std::ostream & os, matrix<T, R, C> const & m)
 	{
-		for (std::size_t i = 0; i < R; ++i)
+		for (std::size_t i = 0; i < m.rows(); ++i)
 		{
-			for (std::size_t j = 0; j < C; ++j)
+			for (std::size_t j = 0; j < m.columns(); ++j)
 				os << m[i][j] << ' ';
 			os << '\n';
 		}
@@ -458,9 +519,9 @@ namespace psemek::math
 	template <typename T, std::size_t R, std::size_t C>
 	std::ostream & operator << (std::ostream & os, setw<T, R, C> const & w)
 	{
-		for (std::size_t i = 0; i < R; ++i)
+		for (std::size_t i = 0; i < w.m.rows(); ++i)
 		{
-			for (std::size_t j = 0; j < C; ++j)
+			for (std::size_t j = 0; j < w.m.columns(); ++j)
 				os << std::fixed << std::right << std::setw(w.w) << w.m[i][j] << ' ';
 			os << '\n';
 		}
diff --git a/libs/math/include/psemek/math/point.hpp b/libs/math/include/psemek/math/point.hpp
index afb4b004..7aa26924 100644
--- a/libs/math/include/psemek/math/point.hpp
+++ b/libs/math/include/psemek/math/point.hpp
@@ -19,17 +19,42 @@ namespace psemek::math
 
 		typename detail::array<T, N>::type coords;
 
-		point() = default;
+		using dynamic_dimensions_type = detail::dynamic_dimensions<N>;
+
+		// Internal
+		point(dynamic_dimensions_type dimensions)
+			: coords(dimensions)
+		{}
+
+		point()
+			: coords({})
+		{}
+
+		explicit point(std::size_t size) requires (N == dynamic)
+			: coords({.size = size})
+		{}
 
 		template <typename ... Args>
-		requires((sizeof...(Args) == N) && detail::all_convertible_to<T, Args...>::value)
+		requires(N != dynamic && sizeof...(Args) == N && detail::all_convertible_to<T, Args...>::value)
 		point(Args && ... args)
-			: coords{ static_cast<T>(std::forward<Args>(args))... }
-		{}
+			: coords({})
+		{
+			auto out = values().begin();
+			((*out++ = args), ...);
+		}
+
+		template <typename ... Args>
+		requires(N == dynamic && detail::all_convertible_to<T, Args...>::value)
+		point(Args && ... args)
+			: coords({.size = sizeof...(Args)})
+		{
+			auto out = values().begin();
+			((*out++ = args), ...);
+		}
 
 		std::size_t dimension() const
 		{
-			return N;
+			return coords.size;
 		}
 
 		T & operator[](std::size_t i)
@@ -42,20 +67,56 @@ namespace psemek::math
 			return coords[i];
 		}
 
+		util::span<T> values()
+		{
+			return {&coords[0], dimension()};
+		}
+
+		util::span<T const> values() const
+		{
+			return {&coords[0], dimension()};
+		}
+
+		point copy() const
+		{
+			if constexpr (N == dynamic)
+			{
+				point result;
+				result.coords = coords.copy();
+				return result;
+			}
+			else
+			{
+				return *this;
+			}
+		}
+
 		point & operator += (vector<T, N> const & v);
 		point & operator -= (vector<T, N> const & v);
 
-		static point zero();
+		static point zero(dynamic_dimensions_type dimensions = {});
 	};
 
 	template <typename ... Args>
 	point(Args && ...) -> point<std::common_type_t<Args...>, sizeof...(Args)>;
 
 	template <typename T, std::size_t N>
-	point<T, N> point<T, N>::zero()
+	void check_dynamic_size(point<T, N> const & p1, point<T, N> const & p2)
 	{
-		point<T, N> p;
-		for (std::size_t i = 0; i < N; ++i)
+		check_dynamic_size(p1.coords.dimensions(), p2.coords.dimensions());
+	}
+
+	template <typename T, std::size_t N>
+	void check_dynamic_size(point<T, N> const & p, vector<T, N> const & v)
+	{
+		check_dynamic_size(p.coords.dimensions(), v.coords.dimensions());
+	}
+
+	template <typename T, std::size_t N>
+	point<T, N> point<T, N>::zero(dynamic_dimensions_type dimensions)
+	{
+		point<T, N> p(dimensions);
+		for (std::size_t i = 0; i < p.dimension(); ++i)
 			p[i] = T(0);
 		return p;
 	}
@@ -63,8 +124,8 @@ namespace psemek::math
 	template <typename T1, typename T, std::size_t N>
 	point<T1, N> cast(point<T, N> const & p)
 	{
-		point<T1, N> r;
-		for (std::size_t i = 0; i < N; ++i)
+		point<T1, N> r(p.coords.dimensions());
+		for (std::size_t i = 0; i < r.dimension(); ++i)
 			r[i] = T1(p[i]);
 		return r;
 	}
@@ -108,7 +169,9 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	std::strong_ordering operator <=> (point<T, N> const & p1, point<T, N> const & p2)
 	{
-		for (std::size_t i = 0; i < N; ++i)
+		check_dynamic_size(p1, p2);
+
+		for (std::size_t i = 0; i < p1.dimension(); ++i)
 		{
 			if (p1[i] < p2[i])
 				return std::strong_ordering::less;
@@ -121,7 +184,8 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	point<T, N> operator + (point<T, N> const & p, vector<T, N> const & v)
 	{
-		point<T, N> r;
+		check_dynamic_size(p, v);
+		point<T, N> r(p.coords.dimensions());
 		for (std::size_t i = 0; i < N; ++i)
 			r[i] = p[i] + v[i];
 		return r;
@@ -130,7 +194,8 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	point<T, N> operator + (vector<T, N> const & v, point<T, N> const & p)
 	{
-		point<T, N> r;
+		check_dynamic_size(p, v);
+		point<T, N> r(p.coords.dimensions());
 		for (std::size_t i = 0; i < N; ++i)
 			r[i] = v[i] + p[i];
 		return r;
@@ -139,7 +204,8 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	point<T, N> operator - (point<T, N> const & p, vector<T, N> const & v)
 	{
-		point<T, N> r;
+		check_dynamic_size(p, v);
+		point<T, N> r(p.coords.dimensions());
 		for (std::size_t i = 0; i < N; ++i)
 			r[i] = p[i] - v[i];
 		return r;
@@ -148,7 +214,8 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	vector<T, N> operator - (point<T, N> const & p1, point<T, N> const & p2)
 	{
-		vector<T, N> r;
+		check_dynamic_size(p1, p2);
+		vector<T, N> r(p1.coords.dimensions());
 		for (std::size_t i = 0; i < N; ++i)
 			r[i] = p1[i] - p2[i];
 		return r;
@@ -199,8 +266,14 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	std::ostream & operator << (std::ostream & os, point<T, N> const & p)
 	{
+		if constexpr (N == 0)
+		{
+			os << "()";
+			return os;
+		}
+
 		os << '(' << p[0];
-		for (std::size_t i = 1; i < N; ++i)
+		for (std::size_t i = 1; i < p.dimension(); ++i)
 			os << ", " << p[i];
 		os << ')';
 		return os;
@@ -209,7 +282,7 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	bool isfinite(point<T, N> const & p)
 	{
-		for (std::size_t i = 0; i < N; ++i)
+		for (std::size_t i = 0; i < p.dimension(); ++i)
 			if (!std::isfinite(p[i]))
 				return false;
 		return true;
diff --git a/libs/math/include/psemek/math/vector.hpp b/libs/math/include/psemek/math/vector.hpp
index 0be0e1e7..dc067e0f 100644
--- a/libs/math/include/psemek/math/vector.hpp
+++ b/libs/math/include/psemek/math/vector.hpp
@@ -2,11 +2,10 @@
 
 #include <psemek/math/detail/array.hpp>
 #include <psemek/util/hash.hpp>
-#include <psemek/util/assert.hpp>
+#include <psemek/util/span.hpp>
 
 #include <iostream>
 #include <cstddef>
-#include <utility>
 #include <type_traits>
 #include <cmath>
 #include <compare>
@@ -22,33 +21,62 @@ namespace psemek::math
 
 		typename detail::array<T, N>::type coords;
 
-		vector() = default;
+		using dynamic_dimensions_type = detail::dynamic_dimensions<N>;
 
-		explicit vector(std::size_t size) requires (N != dynamic)
-		{
-			assert(size == N);
-		}
+		// Internal
+		vector(dynamic_dimensions_type dimensions)
+			: coords(dimensions)
+		{}
+
+		vector()
+			: coords({})
+		{}
 
 		explicit vector(std::size_t size) requires (N == dynamic)
-			: coords(size, T{})
+			: coords({.size = size})
 		{}
 
 		template <typename ... Args>
-		requires((sizeof...(Args) == N || N == dynamic) && detail::all_convertible_to<T, Args...>::value)
+		requires(N != dynamic && sizeof...(Args) == N && detail::all_convertible_to<T, Args...>::value)
 		vector(Args && ... args)
-			: coords{ static_cast<T>(std::forward<Args>(args))... }
-		{}
+			: coords({})
+		{
+			auto out = values().begin();
+			((*out++ = args), ...);
+		}
+
+		template <typename ... Args>
+		requires(N == dynamic && detail::all_convertible_to<T, Args...>::value)
+		vector(Args && ... args)
+			: coords({.size = sizeof...(Args)})
+		{
+			auto out = values().begin();
+			((*out++ = args), ...);
+		}
 
 		std::size_t dimension() const
 		{
-			if constexpr (N == dynamic)
-			{
-				return coords.size();
-			}
-			else
-			{
-				return N;
-			}
+			return coords.size;
+		}
+
+		T & operator[](std::size_t i)
+		{
+			return coords[i];
+		}
+
+		T const & operator[](std::size_t i) const
+		{
+			return coords[i];
+		}
+
+		util::span<T> values()
+		{
+			return {&coords[0], dimension()};
+		}
+
+		util::span<T const> values() const
+		{
+			return {&coords[0], dimension()};
 		}
 
 		vector copy() const
@@ -65,33 +93,28 @@ namespace psemek::math
 			}
 		}
 
-		T & operator[](std::size_t i)
-		{
-			return coords[i];
-		}
-
-		T const & operator[](std::size_t i) const
-		{
-			return coords[i];
-		}
-
 		vector & operator *= (T const & s);
 		vector & operator /= (T const & s);
 
 		vector & operator += (vector const & v);
 		vector & operator -= (vector const & v);
 
-		static vector zero() requires (N != dynamic);
-		static vector zero(std::size_t size);
+		static vector zero(dynamic_dimensions_type dimensions = {});
 	};
 
 	template <typename ... Args>
 	vector(Args && ...) -> vector<std::common_type_t<Args...>, sizeof...(Args)>;
 
+	template <typename T, std::size_t N>
+	void check_dynamic_size(vector<T, N> const & v1, vector<T, N> const & v2)
+	{
+		check_dynamic_size(v1.coords.dimensions(), v2.coords.dimensions());
+	}
+
 	template <typename T1, typename T, std::size_t N>
 	vector<T1, N> cast(vector<T, N> const & v)
 	{
-		vector<T1, N> r(v.dimension());
+		vector<T1, N> r(v.coords.dimensions());
 		for (std::size_t i = 0; i < v.dimension(); ++i)
 			r[i] = static_cast<T1>(v[i]);
 		return r;
@@ -100,7 +123,7 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	std::strong_ordering operator <=> (vector<T, N> const & v1, vector<T, N> const & v2)
 	{
-		check_dynamic_size(v1.dimension(), v2.dimension());
+		check_dynamic_size(v1, v2);
 
 		for (std::size_t i = 0; i < v1.dimension(); ++i)
 		{
@@ -151,7 +174,7 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	vector<T, N> operator * (vector<T, N> const & v, T const & s)
 	{
-		vector<T, N> r(v.dimension());
+		vector<T, N> r(v.coords.dimensions());
 		for (std::size_t i = 0; i < v.dimension(); ++i)
 			r[i] = v[i] * s;
 		return r;
@@ -160,7 +183,7 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	vector<T, N> operator * (T const & s, vector<T, N> const & v)
 	{
-		vector<T, N> r(v.dimension());
+		vector<T, N> r(v.coords.dimensions());
 		for (std::size_t i = 0; i < v.dimension(); ++i)
 			r[i] = s * v[i];
 		return r;
@@ -169,7 +192,7 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	vector<T, N> operator / (vector<T, N> const & v, T const & s)
 	{
-		vector<T, N> r(v.dimension());
+		vector<T, N> r(v.coords.dimensions());
 		for (std::size_t i = 0; i < v.dimension(); ++i)
 			r[i] = v[i] / s;
 		return r;
@@ -178,7 +201,7 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	vector<T, N> operator - (vector<T, N> const & v)
 	{
-		vector<T, N> r(v.dimension());
+		vector<T, N> r(v.coords.dimensions());
 		for (std::size_t i = 0; i < v.dimension(); ++i)
 			r[i] = -v[i];
 		return r;
@@ -187,9 +210,9 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	vector<T, N> operator + (vector<T, N> const & v1, vector<T, N> const & v2)
 	{
-		check_dynamic_size(v1.dimension(), v2.dimension());
+		check_dynamic_size(v1, v2);
 
-		vector<T, N> r(v1.dimension());
+		vector<T, N> r(v1.coords.dimensions());
 		for (std::size_t i = 0; i < v1.dimension(); ++i)
 			r[i] = v1[i] + v2[i];
 		return r;
@@ -198,9 +221,9 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	vector<T, N> operator - (vector<T, N> const & v1, vector<T, N> const & v2)
 	{
-		check_dynamic_size(v1.dimension(), v2.dimension());
+		check_dynamic_size(v1, v2);
 
-		vector<T, N> r(v1.dimension());
+		vector<T, N> r(v1.coords.dimensions());
 		for (std::size_t i = 0; i < v1.dimension(); ++i)
 			r[i] = v1[i] - v2[i];
 		return r;
@@ -231,24 +254,10 @@ namespace psemek::math
 	}
 
 	template <typename T, std::size_t N>
-	vector<T, N> vector<T, N>::zero() requires (N != dynamic)
+	vector<T, N> vector<T, N>::zero(dynamic_dimensions_type dimensions)
 	{
-		vector<T, N> r;
-		for (std::size_t i = 0; i < N; ++i)
-			r[i] = T(0);
-		return r;
-	}
-
-	template <typename T, std::size_t N>
-	vector<T, N> vector<T, N>::zero(std::size_t size)
-	{
-		if constexpr (N != dynamic)
-		{
-			assert(N == size);
-		}
-
-		vector<T, N> r(size);
-		for (std::size_t i = 0; i < size; ++i)
+		vector<T, N> r(dimensions);
+		for (std::size_t i = 0; i < r.dimension(); ++i)
 			r[i] = T(0);
 		return r;
 	}
@@ -256,7 +265,7 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	T dot(vector<T, N> const & v1, vector<T, N> const & v2)
 	{
-		check_dynamic_size(v1.dimension(), v2.dimension());
+		check_dynamic_size(v1, v2);
 
 		T r{};
 		for (std::size_t i = 0; i < v1.dimension(); ++i)
@@ -332,8 +341,8 @@ namespace psemek::math
 	template <typename T, std::size_t N>
 	vector<T, N> ort(vector<T, N> const & v)
 	{
-		vector<T, N> result = vector<T, N>::zero(v.dimension());
-		if constexpr ((N % 2) == 0)
+		vector<T, N> result = vector<T, N>::zero(v.coords.dimensions());
+		if ((N != dynamic && (N % 2) == 0) || (N == dynamic && (v.dimension() % 2) == 0))
 		{
 			for (std::size_t i = 0; i < v.dimension(); i += 2)
 			{
@@ -414,9 +423,9 @@ namespace psemek::math
 	{
 		using R = decltype(f(v0[0], vs[0]...));
 
-		(check_dynamic_size(v0.dimension(), vs.dimension()), ...);
+		(check_dynamic_size(v0, vs), ...);
 
-		vector<R, N> result(v0.dimension());
+		vector<R, N> result(v0.coords.dimensions());
 		for (std::size_t i = 0; i < v0.dimension(); ++i)
 			result[i] = f(v0[i], vs[i]...);
 		return result;