diff --git a/src/entity/commands.cpp b/src/entity/commands.cpp
index 4ad4faa..2c141b1 100644
--- a/src/entity/commands.cpp
+++ b/src/entity/commands.cpp
@@ -157,7 +157,7 @@ void place(entity::registry& registry, entity::id eid, entity::id celestial_body
 		}
 		
 		game::component::transform& transform = registry.get<game::component::transform>(eid);
-		transform.local.translation = math::type_cast<float>(double3{x, y, z});
+		transform.local.translation = math::vector<float, 3>(double3{x, y, z});
 		transform.warp = true;
 	}
 	*/
diff --git a/src/game/state/splash.cpp b/src/game/state/splash.cpp
index 616c41d..e69d9fc 100644
--- a/src/game/state/splash.cpp
+++ b/src/game/state/splash.cpp
@@ -29,6 +29,8 @@
 #include "debug/logger.hpp"
 #include "resources/resource-manager.hpp"
 #include "render/material-flags.hpp"
+#include "math/matrix.hpp"
+#include <iostream>
 
 namespace game {
 namespace state {
diff --git a/src/game/system/astronomy.cpp b/src/game/system/astronomy.cpp
index dbcdbd0..098a3dd 100644
--- a/src/game/system/astronomy.cpp
+++ b/src/game/system/astronomy.cpp
@@ -157,7 +157,7 @@ void astronomy::update(double t, double dt)
 		// Update local transform
 		if (orbit.parent != entt::null)
 		{
-			transform.local.translation = math::normalize(math::type_cast<float>(r_eus));
+			transform.local.translation = math::normalize(float3(r_eus));
 			transform.local.rotation = math::type_cast<float>(rotation_eus);
 			transform.local.scale = {1.0f, 1.0f, 1.0f};
 		}
@@ -206,12 +206,12 @@ void astronomy::update(double t, double dt)
 			(
 				math::look_rotation
 				(
-					math::type_cast<float>(-observer_blackbody_direction_eus),
-					math::type_cast<float>(blackbody_up_eus)
+					float3(-observer_blackbody_direction_eus),
+					float3(blackbody_up_eus)
 				)
 			);
 			
-			sun_light->set_color(math::type_cast<float>(observer_blackbody_transmitted_illuminance));
+			sun_light->set_color(float3(observer_blackbody_transmitted_illuminance));
 			
 			// Bounce sun light
 			bounce_illuminance += std::max(0.0, math::dot(bounce_normal, -observer_blackbody_direction_eus)) * observer_blackbody_transmitted_illuminance * bounce_albedo;
@@ -231,7 +231,7 @@ void astronomy::update(double t, double dt)
 			sky_light_illuminance += starlight_illuminance;
 			
 			// Update sky light
-			sky_light->set_color(math::type_cast<float>(sky_light_illuminance));
+			sky_light->set_color(float3(sky_light_illuminance));
 			
 			// Bounce sky light
 			bounce_illuminance += sky_light_illuminance * bounce_albedo;
@@ -240,9 +240,9 @@ void astronomy::update(double t, double dt)
 		// Upload blackbody params to sky pass
 		if (this->sky_pass)
 		{
-			this->sky_pass->set_sun_position(math::type_cast<float>(blackbody_position_eus));
-			this->sky_pass->set_sun_luminance(math::type_cast<float>(blackbody.luminance));
-			this->sky_pass->set_sun_illuminance(math::type_cast<float>(observer_blackbody_illuminance), math::type_cast<float>(observer_blackbody_transmitted_illuminance));
+			this->sky_pass->set_sun_position(float3(blackbody_position_eus));
+			this->sky_pass->set_sun_luminance(float3(blackbody.luminance));
+			this->sky_pass->set_sun_illuminance(float3(observer_blackbody_illuminance), float3(observer_blackbody_transmitted_illuminance));
 			this->sky_pass->set_sun_angular_radius(static_cast<float>(observer_blackbody_angular_radius));
 		}
 		
@@ -308,23 +308,23 @@ void astronomy::update(double t, double dt)
 				this->sky_pass->set_moon_position(transform.local.translation);
 				this->sky_pass->set_moon_rotation(transform.local.rotation);
 				this->sky_pass->set_moon_angular_radius(static_cast<float>(observer_reflector_angular_radius));
-				this->sky_pass->set_moon_sunlight_direction(math::type_cast<float>(-reflector_blackbody_direction_eus));
-				this->sky_pass->set_moon_sunlight_illuminance(math::type_cast<float>(reflector_blackbody_illuminance * observer_reflector_transmittance));
-				this->sky_pass->set_moon_planetlight_direction(math::type_cast<float>(observer_reflector_direction_eus));
-				this->sky_pass->set_moon_planetlight_illuminance(math::type_cast<float>(reflector_observer_illuminance * observer_reflector_transmittance));
-				this->sky_pass->set_moon_illuminance(math::type_cast<float>(observer_reflector_illuminance / observer_reflector_transmittance), math::type_cast<float>(observer_reflector_illuminance));
+				this->sky_pass->set_moon_sunlight_direction(float3(-reflector_blackbody_direction_eus));
+				this->sky_pass->set_moon_sunlight_illuminance(float3(reflector_blackbody_illuminance * observer_reflector_transmittance));
+				this->sky_pass->set_moon_planetlight_direction(float3(observer_reflector_direction_eus));
+				this->sky_pass->set_moon_planetlight_illuminance(float3(reflector_observer_illuminance * observer_reflector_transmittance));
+				this->sky_pass->set_moon_illuminance(float3(observer_reflector_illuminance / observer_reflector_transmittance), float3(observer_reflector_illuminance));
 			}
 			
 			if (this->moon_light)
 			{
-				const float3 reflector_up_eus = math::type_cast<float>(icrf_to_eus.r * double3{0, 0, 1});
+				const float3 reflector_up_eus = float3(icrf_to_eus.r * double3{0, 0, 1});
 				
-				this->moon_light->set_color(math::type_cast<float>(observer_reflector_illuminance));
+				this->moon_light->set_color(float3(observer_reflector_illuminance));
 				this->moon_light->set_rotation
 				(
 					math::look_rotation
 					(
-						math::type_cast<float>(-observer_reflector_direction_eus),
+						float3(-observer_reflector_direction_eus),
 						reflector_up_eus
 					)
 				);
@@ -337,7 +337,7 @@ void astronomy::update(double t, double dt)
 	
 	if (bounce_light)
 	{
-		bounce_light->set_color(math::type_cast<float>(bounce_illuminance));
+		bounce_light->set_color(float3(bounce_illuminance));
 	}
 }
 
@@ -580,7 +580,7 @@ void astronomy::update_icrf_to_eus(const game::component::celestial_body& body,
 		(
 			math::transformation::se3<float>
 			{
-				math::type_cast<float>(icrf_to_eus.t),
+				float3(icrf_to_eus.t),
 				math::type_cast<float>(icrf_to_eus.r)
 			}
 		);
diff --git a/src/game/system/atmosphere.cpp b/src/game/system/atmosphere.cpp
index b03761d..22874e3 100644
--- a/src/game/system/atmosphere.cpp
+++ b/src/game/system/atmosphere.cpp
@@ -143,10 +143,10 @@ void atmosphere::update_sky_pass()
 		return;
 	
 	sky_pass->set_atmosphere_upper_limit(static_cast<float>(component->upper_limit));
-	sky_pass->set_rayleigh_parameters(static_cast<float>(component->rayleigh_scale_height), math::type_cast<float>(component->rayleigh_scattering));
+	sky_pass->set_rayleigh_parameters(static_cast<float>(component->rayleigh_scale_height), math::vector<float, 3>(component->rayleigh_scattering));
 	sky_pass->set_mie_parameters(static_cast<float>(component->mie_scale_height), static_cast<float>(component->mie_scattering), static_cast<float>(component->mie_extinction), static_cast<float>(component->mie_anisotropy));
-	sky_pass->set_ozone_parameters(static_cast<float>(component->ozone_lower_limit), static_cast<float>(component->ozone_upper_limit), static_cast<float>(component->ozone_mode), math::type_cast<float>(component->ozone_absorption));
-	sky_pass->set_airglow_illuminance(math::type_cast<float>(component->airglow_illuminance));
+	sky_pass->set_ozone_parameters(static_cast<float>(component->ozone_lower_limit), static_cast<float>(component->ozone_upper_limit), static_cast<float>(component->ozone_mode), math::vector<float, 3>(component->ozone_absorption));
+	sky_pass->set_airglow_illuminance(math::vector<float, 3>(component->airglow_illuminance));
 }
 
 void atmosphere::on_atmosphere_construct(entity::registry& registry, entity::id entity_id)
diff --git a/src/geom/plane.hpp b/src/geom/plane.hpp
index fb8b236..4fa408e 100644
--- a/src/geom/plane.hpp
+++ b/src/geom/plane.hpp
@@ -99,7 +99,7 @@ plane<T>::plane(const vector_type& a, const vector_type& b, const vector_type& c
 template <class T>
 plane<T>::plane(const math::vector<T, 4>& coefficients)
 {
-	const vector_type abc = math::resize<3>(coefficients);
+	const vector_type abc = math::vector<T, 3>(coefficients);
 	const float inverse_length = T(1) / math::length(abc);
 	
 	normal = abc * inverse_length;
diff --git a/src/math/matrix.hpp b/src/math/matrix.hpp
index 2c83a4f..4d46974 100644
--- a/src/math/matrix.hpp
+++ b/src/math/matrix.hpp
@@ -134,7 +134,7 @@ struct matrix
 	/**
 	 * Returns a pointer to the first element.
 	 *
-	 * @warning If column_vector_type is not a POD type, elements may not be stored contiguously.
+	 * @warning If matrix::element_type is not a POD type, elements may not be stored contiguously.
 	 */
 	/// @{
 	constexpr inline element_type* data() noexcept
@@ -242,18 +242,18 @@ struct matrix
 	constexpr inline matrix<T, P, O> size_cast(std::index_sequence<I...>) const noexcept
 	{
 		if constexpr (O == M)
-			return {((I < N) ? columns[I] : vector<T, O>::zero()) ...};
+			return {((I < N) ? columns[I] : matrix<T, P, O>::identity()[I]) ...};
 		else
-			return {((I < N) ? vector<T, O>(columns[I]) : vector<T, O>::zero()) ...};
+			return {((I < N) ? vector<T, O>(columns[I]) : matrix<T, P, O>::identity()[I]) ...};
 	}
 	
 	/**
-	 * Size-casts this matrix to a matrix with different dimensions. Casting to greater dimensions causes new elements to be set to zero.
+	 * Size-casts this matrix to a matrix with different dimensions. Casting to greater dimensions causes new elements to be set to identity matrix elements.
 	 *
 	 * @tparam P Target number of columns.
 	 * @tparam O Target number of rows.
 	 *
-	 * @return *p by *o* matrix.
+	 * @return *p* by *o* matrix.
 	 */
 	template <std::size_t P, std::size_t O>
 	constexpr inline explicit operator matrix<T, P, O>() const noexcept
@@ -331,63 +331,48 @@ using matrix4x4 = matrix<T, 4, 4>;
 /**
  * Adds two matrices.
  *
- * @param x First matrix.
- * @param y Second matrix.
+ * @param a First matrix.
+ * @param b Second matrix.
+ *
  * @return Sum of the two matrices.
  */
-template <class T>
-constexpr matrix<T, 2, 2> add(const matrix<T, 2, 2>& x, const matrix<T, 2, 2>& y);
-
-/// @copydoc add(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&)
-template <class T>
-constexpr matrix<T, 3, 3> add(const matrix<T, 3, 3>& x, const matrix<T, 3, 3>& y);
-
-/// @copydoc add(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&)
-template <class T>
-constexpr matrix<T, 4, 4> add(const matrix<T, 4, 4>& x, const matrix<T, 4, 4>& y);
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> add(const matrix<T, N, M>& a, const matrix<T, N, M>& b) noexcept;
 
 /**
- * Reinterprets data as an `NxM` matrix of type `T`.
+ * Adds a matrix and a scalar.
  *
- * @tparam N Number of columns.
- * @tparam M Number of rows.
- * @tparam T Element type.
- * @param data Data to reinterpret.
+ * @param a Matrix.
+ * @param b scalar.
+ *
+ * @return Sum of the matrix and scalar.
  */
-template <std::size_t N, std::size_t M, typename T>
-constexpr matrix<T, N, M>& as_matrix(T& data);
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> add(const matrix<T, N, M>& a, T b) noexcept;
 
 /**
- * Calculates the determinant of a matrix.
+ * Calculates the determinant of a square matrix.
  *
- * @param m Matrix of which to take the determinant. 
+ * @param m Matrix of which to take the determinant.
+ *
+ * @return Determinant of @p m.
+ *
+ * @warning Currently only implemented for 2x2, 3x3, and 4x4 matrices.
  */
-template <class T>
-constexpr T determinant(const matrix<T, 2, 2>& m);
-
-/// @copydoc determinant(const matrix<T, 2, 2>&)
-template <class T>
-constexpr T determinant(const matrix<T, 3, 3>& m);
-
-/// @copydoc determinant(const matrix<T, 2, 2>&)
-template <class T>
-constexpr T determinant(const matrix<T, 4, 4>& m);
+template <class T, std::size_t N>
+constexpr T determinant(const matrix<T, N, N>& m) noexcept;
 
 /**
- * Calculates the inverse of a matrix.
+ * Calculates the inverse of a square matrix.
  *
  * @param m Matrix of which to take the inverse.
+ *
+ * @return Inverse of matrix @p m.
+ *
+ * @warning Currently only implemented for 2x2, 3x3, and 4x4 matrices.
  */
-template <class T>
-constexpr matrix<T, 2, 2> inverse(const matrix<T, 2, 2>& m);
-
-/// @copydoc inverse(const matrix<T, 2, 2>&)
-template <class T>
-constexpr matrix<T, 3, 3> inverse(const matrix<T, 3, 3>& m);
-
-/// @copydoc inverse(const matrix<T, 2, 2>&)
-template <class T>
-constexpr matrix<T, 4, 4> inverse(const matrix<T, 4, 4>& m);
+template <class T, std::size_t N>
+constexpr matrix<T, N, N> inverse(const matrix<T, N, N>& m) noexcept;
 
 /**
  * Performs a component-wise multiplication of two matrices.
@@ -395,151 +380,97 @@ constexpr matrix<T, 4, 4> inverse(const matrix<T, 4, 4>& m);
  * @param x First matrix multiplicand.
  * @param y Second matrix multiplicand.
  */
-template <class T>
-constexpr matrix<T, 2, 2> componentwise_mul(const matrix<T, 2, 2>& x, const matrix<T, 2, 2>& y);
-
-/// @copydoc componentwise_mul(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&)
-template <class T>
-constexpr matrix<T, 3, 3> componentwise_mul(const matrix<T, 3, 3>& x, const matrix<T, 3, 3>& y);
-
-/// @copydoc componentwise_mul(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&)
-template <class T>
-constexpr matrix<T, 4, 4> componentwise_mul(const matrix<T, 4, 4>& x, const matrix<T, 4, 4>& y);
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> componentwise_mul(const matrix<T, N, M>& a, const matrix<T, N, M>& b) noexcept;
 
 /**
- * Creates a viewing transformation matrix.
+ * Divides a matrix by a matrix.
  *
- * @param position Position of the view point.
- * @param target Position of the target.
- * @param up Normalized direction of the up vector.
- * @return Viewing transformation matrix.
+ * @param a First matrix.
+ * @param b Second matrix.
+ * @return Result of the division.
  */
-template <class T>
-constexpr matrix<T, 4, 4> look_at(const vector<T, 3>& position, const vector<T, 3>& target, vector<T, 3> up);
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> div(const matrix<T, N, M>& a, const matrix<T, N, M>& b) noexcept;
 
 /**
- * Multiplies two matrices.
+ * Divides a matrix by a scalar.
  *
- * @param x First matrix.
- * @param y Second matrix.
- * @return Product of the two matrices.
+ * @param a Matrix.
+ * @param b Scalar.
+ * @return Result of the division.
  */
-template <class T>
-constexpr matrix<T, 2, 2> mul(const matrix<T, 2, 2>& x, const matrix<T, 2, 2>& y);
-
-/// @copydoc mul(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&);
-template <class T>
-constexpr matrix<T, 3, 3> mul(const matrix<T, 3, 3>& x, const matrix<T, 3, 3>& y);
-
-/// @copydoc mul(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&);
-template <class T>
-constexpr matrix<T, 4, 4> mul(const matrix<T, 4, 4>& x, const matrix<T, 4, 4>& y);
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> div(const matrix<T, N, M>& a, T b) noexcept;
 
 /**
- * Multiplies a matrix by a scalar.
+ * Divides a scalar by a matrix.
  *
- * @param m Matrix.
- * @param s Scalar.
- * @return Product of the matrix and the scalar..
+ * @param a Scalar.
+ * @param b Matrix.
+ * @return Result of the division.
  */
 template <class T, std::size_t N, std::size_t M>
-constexpr matrix<T, N, M> mul(const matrix<T, N, M>& m, T s);
+constexpr matrix<T, N, M> div(T a, const matrix<T, N, M>& b) noexcept;
 
 /**
- * Transforms a vector by a matrix.
+ * Creates a viewing transformation matrix.
  *
- * @param m Matrix.
- * @param v Vector.
- * @return Transformed vector.
+ * @param position Position of the view point.
+ * @param target Position of the target.
+ * @param up Normalized direction of the up vector.
+ * @return Viewing transformation matrix.
  */
 template <class T>
-constexpr vector<T, 2> mul(const matrix<T, 2, 2>& m, const vector<T, 2>& v);
-
-/// @copydoc mul(const matrix<T, 2, 2>&, const vector<T, 2>&)
-template <class T>
-constexpr vector<T, 3> mul(const matrix<T, 3, 3>& m, const vector<T, 3>& v);
-
-/// @copydoc mul(const matrix<T, 2, 2>&, const vector<T, 2>&)
-template <class T>
-constexpr vector<T, 4> mul(const matrix<T, 4, 4>& m, const vector<T, 4>& v);
+constexpr matrix<T, 4, 4> look_at(const vector<T, 3>& position, const vector<T, 3>& target, vector<T, 3> up);
 
 /**
- * Creates an orthographic projection matrix which will transform the near and far clipping planes to `[-1, 1]`, respectively.
+ * Multiplies two matrices
  *
- * @param left Signed distance to the left clipping plane.
- * @param right Signed distance to the right clipping plane.
- * @param bottom Signed distance to the bottom clipping plane.
- * @param top Signed distance to the top clipping plane.
- * @param z_near Signed distance to the near clipping plane.
- * @param z_far Signed distance to the far clipping plane.
- * @return Orthographic projection matrix.
- */
-template <class T>
-constexpr matrix<T, 4, 4> ortho(T left, T right, T bottom, T top, T z_near, T z_far);
-
-/**
- * Creates an orthographic projection matrix which will transform the near and far clipping planes to `[0, 1]`, respectively.
+ * @tparam T Matrix element type.
+ * @tparam N Number of columns in matrix @p a, and rows in matrix @p b.
+ * @tparam M Number of rows in matrix @p a.
+ * @tparam P Number of columns in matrix @p b.
  *
- * @param left Signed distance to the left clipping plane.
- * @param right Signed distance to the right clipping plane.
- * @param bottom Signed distance to the bottom clipping plane.
- * @param top Signed distance to the top clipping plane.
- * @param z_near Signed distance to the near clipping plane.
- * @param z_far Signed distance to the far clipping plane.
- * @return Orthographic projection matrix.
- */
-template <class T>
-constexpr matrix<T, 4, 4> ortho_half_z(T left, T right, T bottom, T top, T z_near, T z_far);
-
-/**
- * Calculates the outer product of a pair of vectors.
+ * @param a First matrix.
+ * @param b Second matrix.
  *
- * @param c Parameter to be treated as a column vector.
- * @param r Parameter to be treated as a row vector.
+ * @return Product of `a * b`.
  */
-template <class T>
-constexpr matrix<T, 2, 2> outer_product(const vector<T, 2>& c, const vector<T, 2>& r);
-
-/// @copydoc outer_product(const vector<T, 2>&, const vector<T, 2>&)
-template <class T>
-constexpr matrix<T, 3, 3> outer_product(const vector<T, 3>& c, const vector<T, 3>& r);
-
-/// @copydoc outer_product(const vector<T, 2>&, const vector<T, 2>&)
-template <class T>
-constexpr matrix<T, 4, 4> outer_product(const vector<T, 4>& c, const vector<T, 4> r);
+template <typename T, std::size_t N, std::size_t M, std::size_t P>
+constexpr matrix<T, P, M> mul(const matrix<T, N, M>& a, const matrix<T, P, N>& b) noexcept;
 
 /**
- * Creates a perspective projection matrix which will transform the near and far clipping planes to `[-1, 1]`, respectively.
+ * Multiplies a matrix by a scalar.
  *
- * @param vertical_fov Vertical field of view angle, in radians.
- * @param aspect_ratio Aspect ratio which determines the horizontal field of view.
- * @param z_near Distance to the near clipping plane.
- * @param z_far Distance to the far clipping plane.
- * @return Perspective projection matrix.
+ * @param a Matrix.
+ * @param b Scalar.
+ * @return Product of the matrix and the scalar.
  */
-template <class T>
-matrix<T, 4, 4> perspective(T vertical_fov, T aspect_ratio, T z_near, T z_far);
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> mul(const matrix<T, N, M>& a, T b) noexcept;
 
 /**
- * Creates a perspective projection matrix which will transform the near and far clipping planes to `[0, 1]`, respectively.
+ * Calculates the product of a matrix and a row vector.
  *
- * @param vertical_fov Vertical field of view angle, in radians.
- * @param aspect_ratio Aspect ratio which determines the horizontal field of view.
- * @param z_near Distance to the near clipping plane.
- * @param z_far Distance to the far clipping plane.
- * @return Perspective projection matrix.
+ * @param a Matrix.
+ * @param b Row vector
+ *
+ * @return Product of the matrix and the row vector.
  */
-template <class T>
-matrix<T, 4, 4> perspective_half_z(T vertical_fov, T aspect_ratio, T z_near, T z_far);
+template <typename T, std::size_t N, std::size_t M>
+constexpr typename matrix<T, N, M>::column_vector_type mul(const matrix<T, N, M>& a, const typename matrix<T, N, M>::row_vector_type& b) noexcept;
 
 /**
- * Resizes a matrix. Any new elements will be set to `1` if in the diagonal, and `0` otherwise.
+ * Calculates the product of a column vector and a matrix.
  *
- * @param m Matrix to resize.
- * @return Resized matrix.
+ * @param a Column vector.
+ * @param b Matrix.
+ *
+ * @return Product of the column vector and the matrix.
  */
-template <std::size_t N1, std::size_t M1, class T, std::size_t N0, std::size_t M0>
-constexpr matrix<T, N1, M1> resize(const matrix<T, N0, M0>& m);
+template <typename T, std::size_t N, std::size_t M>
+constexpr typename matrix<T, N, M>::row_vector_type mul(const typename matrix<T, N, M>::column_vector_type& a, const matrix<T, N, M>& b) noexcept;
 
 /**
  * Constructs a rotation matrix.
@@ -591,20 +522,35 @@ constexpr matrix<T, 4, 4> scale(const matrix<T, 4, 4>& m, const vector<T, 3>& v)
 /**
  * Subtracts a matrix from another matrix.
  *
- * @param x First matrix.
- * @param y Second matrix.
+ * @param a First matrix.
+ * @param b Second matrix.
+ *
  * @return Difference between the two matrices.
  */
-template <class T>
-constexpr matrix<T, 2, 2> sub(const matrix<T, 2, 2>& x, const matrix<T, 2, 2>& y);
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> sub(const matrix<T, N, M>& a, const matrix<T, N, M>& b) noexcept;
 
-/// @copydoc sub(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&)
-template <class T>
-constexpr matrix<T, 3, 3> sub(const matrix<T, 3, 3>& x, const matrix<T, 3, 3>& y);
+/**
+ * Subtracts a scalar from matrix.
+ *
+ * @param a Matrix.
+ * @param b Scalar.
+ *
+ * @return Difference between the matrix and scalar.
+ */
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> sub(const matrix<T, N, M>& a, T b) noexcept;
 
-/// @copydoc sub(const matrix<T, 2, 2>&, const matrix<T, 2, 2>&)
-template <class T>
-constexpr matrix<T, 4, 4> sub(const matrix<T, 4, 4>& x, const matrix<T, 4, 4>& y);
+/**
+ * Subtracts a matrix from a scalar.
+ *
+ * @param a Scalar.
+ * @param b Matrix.
+ *
+ * @return Difference between the scalar and matrix.
+ */
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> sub(T a, const matrix<T, N, M>& b) noexcept;
 
 /**
  * Translates a matrix.
@@ -619,82 +565,52 @@ constexpr matrix<T, 4, 4> translate(const matrix<T, 4, 4>& m, const vector<T, 3>
 /**
  * Calculates the transpose of a matrix.
  *
- * @param m Matrix of which to take the transpose.
- */
-template <class T>
-constexpr matrix<T, 2, 2> transpose(const matrix<T, 2, 2>& m);
-
-/// @copydoc transpose(const matrix<T, 2, 2>&)
-template <class T>
-constexpr matrix<T, 3, 3> transpose(const matrix<T, 3, 3>& m);
-
-/// @copydoc transpose(const matrix<T, 2, 2>&)
-template <class T>
-constexpr matrix<T, 4, 4> transpose(const matrix<T, 4, 4>& m);
-
-/**
- * Types casts each matrix element and returns a matrix of the casted type.
+ * @param m Matrix to transpose.
  *
- * @tparam T2 Target matrix element type.
- * @tparam T1 Source matrix element type.
- * @tparam N Number of columns.
- * @tparam M Number of rows.
- * @param m Matrix to type cast.
- * @return Type-casted matrix.
+ * @return Transposed matrix.
  */
-template <class T2, class T1, std::size_t N, std::size_t M>
-constexpr matrix<T2, N, M> type_cast(const matrix<T1, N, M>& m);
+template <typename T, std::size_t N, std::size_t M>
+constexpr matrix<T, M, N> transpose(const matrix<T, N, M>& m) noexcept;
 
-template <class T>
-constexpr matrix<T, 2, 2> add(const matrix<T, 2, 2>& x, const matrix<T, 2, 2>& y)
+/// @private
+template <class T, std::size_t N, std::size_t M, std::size_t... I>
+constexpr inline matrix<T, N, M> add(const matrix<T, N, M>& a, const matrix<T, N, M>& b, std::index_sequence<I...>) noexcept
 {
-	return
-		{{
-			x[0] + y[0],
-			x[1] + y[1]
-		}};
+	return {(a[I] + b[I]) ...};
 }
 
-template <class T>
-constexpr matrix<T, 3, 3> add(const matrix<T, 3, 3>& x, const matrix<T, 3, 3>& y)
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> add(const matrix<T, N, M>& a, const matrix<T, N, M>& b) noexcept
 {
-	return
-		{{
-			x[0] + y[0],
-			x[1] + y[1],
-			x[2] + y[2]
-		}};
+	return add(a, b, std::make_index_sequence<N>{});
 }
 
-template <class T>
-constexpr matrix<T, 4, 4> add(const matrix<T, 4, 4>& x, const matrix<T, 4, 4>& y)
+/// @private
+template <class T, std::size_t N, std::size_t M, std::size_t... I>
+constexpr inline matrix<T, N, M> add(const matrix<T, N, M>& a, T b, std::index_sequence<I...>) noexcept
 {
-	return
-		{{
-			x[0] + y[0],
-			x[1] + y[1],
-			x[2] + y[2],
-			x[3] + y[3]
-		}};
+	return {(a[I] + b) ...};
 }
 
-template <std::size_t N, std::size_t M, typename T>
-constexpr inline matrix<T, N, M>& as_matrix(T& data)
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> add(const matrix<T, N, M>& a, T b) noexcept
 {
-	static_assert(std::is_pod<matrix<T, N, M>>::value);
-	return reinterpret_cast<matrix<T, N, M>&>(data);
+	return add(a, b, std::make_index_sequence<N>{});
 }
 
 template <class T>
-constexpr T determinant(const matrix<T, 2, 2>& m)
+constexpr T determinant(const matrix<T, 2, 2>& m) noexcept
 {
-	return m[0][0] * m[1][1] - m[0][1] * m[1][0];
+	return
+		m[0][0] * m[1][1] -
+		m[0][1] * m[1][0];
 }
 
 template <class T>
-constexpr T determinant(const matrix<T, 3, 3>& m)
+constexpr T determinant(const matrix<T, 3, 3>& m) noexcept
 {
-	return m[0][0] * m [1][1] * m[2][2] +
+	return
+		m[0][0] * m[1][1] * m[2][2] +
 		m[0][1] * m[1][2] * m[2][0] +
 		m[0][2] * m[1][0] * m[2][1] -
 		m[0][0] * m[1][2] * m[2][1] -
@@ -703,9 +619,10 @@ constexpr T determinant(const matrix<T, 3, 3>& m)
 }
 
 template <class T>
-constexpr T determinant(const matrix<T, 4, 4>& m)
+constexpr T determinant(const matrix<T, 4, 4>& m) noexcept
 {
-	return m[0][3] * m[1][2] * m[2][1] * m[3][0] - m[0][2] * m[1][3] * m[2][1] * m[3][0] -
+	return
+		m[0][3] * m[1][2] * m[2][1] * m[3][0] - m[0][2] * m[1][3] * m[2][1] * m[3][0] -
 		m[0][3] * m[1][1] * m[2][2] * m[3][0] + m[0][1] * m[1][3] * m[2][2] * m[3][0] +
 		m[0][2] * m[1][1] * m[2][3] * m[3][0] - m[0][1] * m[1][2] * m[2][3] * m[3][0] -
 		m[0][3] * m[1][2] * m[2][0] * m[3][1] + m[0][2] * m[1][3] * m[2][0] * m[3][1] +
@@ -720,317 +637,195 @@ constexpr T determinant(const matrix<T, 4, 4>& m)
 }
 
 template <class T>
-constexpr matrix<T, 2, 2> inverse(const matrix<T, 2, 2>& m)
-{
-	static_assert(std::is_floating_point<T>::value);
-
-	const T rd(T(1) / determinant(m));
-	return
-		{{
-			{ m[1][1] * rd, -m[0][1] * rd},
-			{-m[1][0] * rd,  m[0][0] * rd}
-		}};
-}
-
-template <class T>
-constexpr matrix<T, 3, 3> inverse(const matrix<T, 3, 3>& m)
+constexpr matrix<T, 2, 2> inverse(const matrix<T, 2, 2>& m) noexcept
 {
-	static_assert(std::is_floating_point<T>::value);
-
-	const T rd(T(1) / determinant(m));
-
+	const T inv_det = T{1} / determinant(m);
+	
 	return
-		{{
-			{
-				(m[1][1] * m[2][2] - m[1][2] * m[2][1]) * rd,
-				(m[0][2] * m[2][1] - m[0][1] * m[2][2]) * rd,
-				(m[0][1] * m[1][2] - m[0][2] * m[1][1]) * rd
-			},
-
-			{
-				(m[1][2] * m[2][0] - m[1][0] * m[2][2]) * rd,
-				(m[0][0] * m[2][2] - m[0][2] * m[2][0]) * rd,
-				(m[0][2] * m[1][0] - m[0][0] * m[1][2]) * rd
-			},
-
-			{
-				(m[1][0] * m[2][1] - m[1][1] * m[2][0]) * rd,
-				(m[0][1] * m[2][0] - m[0][0] * m[2][1]) * rd,
-				(m[0][0] * m[1][1] - m[0][1] * m[1][0]) * rd
-			}
-		}};
+	{
+		 m[1][1] * inv_det,
+		-m[0][1] * inv_det,
+		-m[1][0] * inv_det,
+		 m[0][0] * inv_det
+	};
 }
 
 template <class T>
-constexpr matrix<T, 4, 4> inverse(const matrix<T, 4, 4>& m)
+constexpr matrix<T, 3, 3> inverse(const matrix<T, 3, 3>& m) noexcept
 {
-	static_assert(std::is_floating_point<T>::value);
-
-	const T rd(T(1) / determinant(m));
-
-	return
-		{{
-			{
-				(m[1][2] * m[2][3] * m[3][1] - m[1][3] * m[2][2] * m[3][1] + m[1][3] * m[2][1] * m[3][2] - m[1][1] * m[2][3] * m[3][2] - m[1][2] * m[2][1] * m[3][3] + m[1][1] * m[2][2] * m[3][3]) * rd,
-				(m[0][3] * m[2][2] * m[3][1] - m[0][2] * m[2][3] * m[3][1] - m[0][3] * m[2][1] * m[3][2] + m[0][1] * m[2][3] * m[3][2] + m[0][2] * m[2][1] * m[3][3] - m[0][1] * m[2][2] * m[3][3]) * rd,
-				(m[0][2] * m[1][3] * m[3][1] - m[0][3] * m[1][2] * m[3][1] + m[0][3] * m[1][1] * m[3][2] - m[0][1] * m[1][3] * m[3][2] - m[0][2] * m[1][1] * m[3][3] + m[0][1] * m[1][2] * m[3][3]) * rd,
-				(m[0][3] * m[1][2] * m[2][1] - m[0][2] * m[1][3] * m[2][1] - m[0][3] * m[1][1] * m[2][2] + m[0][1] * m[1][3] * m[2][2] + m[0][2] * m[1][1] * m[2][3] - m[0][1] * m[1][2] * m[2][3]) * rd
-			},
-
-			{
-				(m[1][3] * m[2][2] * m[3][0] - m[1][2] * m[2][3] * m[3][0] - m[1][3] * m[2][0] * m[3][2] + m[1][0] * m[2][3] * m[3][2] + m[1][2] * m[2][0] * m[3][3] - m[1][0] * m[2][2] * m[3][3]) * rd,
-				(m[0][2] * m[2][3] * m[3][0] - m[0][3] * m[2][2] * m[3][0] + m[0][3] * m[2][0] * m[3][2] - m[0][0] * m[2][3] * m[3][2] - m[0][2] * m[2][0] * m[3][3] + m[0][0] * m[2][2] * m[3][3]) * rd,
-				(m[0][3] * m[1][2] * m[3][0] - m[0][2] * m[1][3] * m[3][0] - m[0][3] * m[1][0] * m[3][2] + m[0][0] * m[1][3] * m[3][2] + m[0][2] * m[1][0] * m[3][3] - m[0][0] * m[1][2] * m[3][3]) * rd,
-				(m[0][2] * m[1][3] * m[2][0] - m[0][3] * m[1][2] * m[2][0] + m[0][3] * m[1][0] * m[2][2] - m[0][0] * m[1][3] * m[2][2] - m[0][2] * m[1][0] * m[2][3] + m[0][0] * m[1][2] * m[2][3]) * rd
-			},
-
-			{
-				(m[1][1] * m[2][3] * m[3][0] - m[1][3] * m[2][1] * m[3][0] + m[1][3] * m[2][0] * m[3][1] - m[1][0] * m[2][3] * m[3][1] - m[1][1] * m[2][0] * m[3][3] + m[1][0] * m[2][1] * m[3][3]) * rd,
-				(m[0][3] * m[2][1] * m[3][0] - m[0][1] * m[2][3] * m[3][0] - m[0][3] * m[2][0] * m[3][1] + m[0][0] * m[2][3] * m[3][1] + m[0][1] * m[2][0] * m[3][3] - m[0][0] * m[2][1] * m[3][3]) * rd,
-				(m[0][1] * m[1][3] * m[3][0] - m[0][3] * m[1][1] * m[3][0] + m[0][3] * m[1][0] * m[3][1] - m[0][0] * m[1][3] * m[3][1] - m[0][1] * m[1][0] * m[3][3] + m[0][0] * m[1][1] * m[3][3]) * rd,
-				(m[0][3] * m[1][1] * m[2][0] - m[0][1] * m[1][3] * m[2][0] - m[0][3] * m[1][0] * m[2][1] + m[0][0] * m[1][3] * m[2][1] + m[0][1] * m[1][0] * m[2][3] - m[0][0] * m[1][1] * m[2][3]) * rd
-			},
-
-			{
-				(m[1][2] * m[2][1] * m[3][0] - m[1][1] * m[2][2] * m[3][0] - m[1][2] * m[2][0] * m[3][1] + m[1][0] * m[2][2] * m[3][1] + m[1][1] * m[2][0] * m[3][2] - m[1][0] * m[2][1] * m[3][2]) * rd,
-				(m[0][1] * m[2][2] * m[3][0] - m[0][2] * m[2][1] * m[3][0] + m[0][2] * m[2][0] * m[3][1] - m[0][0] * m[2][2] * m[3][1] - m[0][1] * m[2][0] * m[3][2] + m[0][0] * m[2][1] * m[3][2]) * rd,
-				(m[0][2] * m[1][1] * m[3][0] - m[0][1] * m[1][2] * m[3][0] - m[0][2] * m[1][0] * m[3][1] + m[0][0] * m[1][2] * m[3][1] + m[0][1] * m[1][0] * m[3][2] - m[0][0] * m[1][1] * m[3][2]) * rd,
-				(m[0][1] * m[1][2] * m[2][0] - m[0][2] * m[1][1] * m[2][0] + m[0][2] * m[1][0] * m[2][1] - m[0][0] * m[1][2] * m[2][1] - m[0][1] * m[1][0] * m[2][2] + m[0][0] * m[1][1] * m[2][2]) * rd
-			}
-		}};
-}
+	const T inv_det = T{1} / determinant(m);
 
-template <class T>
-constexpr matrix<T, 2, 2> componentwise_mul(const matrix<T, 2, 2>& x, const matrix<T, 2, 2>& y)
-{
 	return
-		{{
-			 {x[0][0] * y[0][0], x[0][1] * y[0][1]},
-			 {x[1][0] * y[1][0], x[1][1] * y[1][1]}
-		}};
-}
+	{
+		(m[1][1] * m[2][2] - m[1][2] * m[2][1]) * inv_det,
+		(m[0][2] * m[2][1] - m[0][1] * m[2][2]) * inv_det,
+		(m[0][1] * m[1][2] - m[0][2] * m[1][1]) * inv_det,
+		
+		(m[1][2] * m[2][0] - m[1][0] * m[2][2]) * inv_det,
+		(m[0][0] * m[2][2] - m[0][2] * m[2][0]) * inv_det,
+		(m[0][2] * m[1][0] - m[0][0] * m[1][2]) * inv_det,
 
-template <class T>
-constexpr matrix<T, 3, 3> componentwise_mul(const matrix<T, 3, 3>& x, const matrix<T, 3, 3>& y)
-{
-	return
-		{{
-			 {x[0][0] * y[0][0], x[0][1] * y[0][1], x[0][2] * y[0][2]},
-			 {x[1][0] * y[1][0], x[1][1] * y[1][1], x[1][2] * y[1][2]},
-			 {x[2][0] * y[2][0], x[2][1] * y[2][1], x[2][2] * y[2][2]}
-		}};
+		(m[1][0] * m[2][1] - m[1][1] * m[2][0]) * inv_det,
+		(m[0][1] * m[2][0] - m[0][0] * m[2][1]) * inv_det,
+		(m[0][0] * m[1][1] - m[0][1] * m[1][0]) * inv_det
+	};
 }
 
 template <class T>
-constexpr matrix<T, 4, 4> componentwise_mul(const matrix<T, 4, 4>& x, const matrix<T, 4, 4>& y)
+constexpr matrix<T, 4, 4> inverse(const matrix<T, 4, 4>& m) noexcept
 {
+	const T inv_det = T{1} / determinant(m);
+	
 	return
-		{{
-			 {x[0][0] * y[0][0], x[0][1] * y[0][1], x[0][2] * y[0][2], x[0][3] * y[0][3]},
-			 {x[1][0] * y[1][0], x[1][1] * y[1][1], x[1][2] * y[1][2], x[1][3] * y[1][3]},
-			 {x[2][0] * y[2][0], x[2][1] * y[2][1], x[2][2] * y[2][2], x[2][3] * y[2][3]},
-			 {x[3][0] * y[3][0], x[3][1] * y[3][1], x[3][2] * y[3][2], x[3][3] * y[3][3]}
-		}};
+	{
+		(m[1][2] * m[2][3] * m[3][1] - m[1][3] * m[2][2] * m[3][1] + m[1][3] * m[2][1] * m[3][2] - m[1][1] * m[2][3] * m[3][2] - m[1][2] * m[2][1] * m[3][3] + m[1][1] * m[2][2] * m[3][3]) * inv_det,
+		(m[0][3] * m[2][2] * m[3][1] - m[0][2] * m[2][3] * m[3][1] - m[0][3] * m[2][1] * m[3][2] + m[0][1] * m[2][3] * m[3][2] + m[0][2] * m[2][1] * m[3][3] - m[0][1] * m[2][2] * m[3][3]) * inv_det,
+		(m[0][2] * m[1][3] * m[3][1] - m[0][3] * m[1][2] * m[3][1] + m[0][3] * m[1][1] * m[3][2] - m[0][1] * m[1][3] * m[3][2] - m[0][2] * m[1][1] * m[3][3] + m[0][1] * m[1][2] * m[3][3]) * inv_det,
+		(m[0][3] * m[1][2] * m[2][1] - m[0][2] * m[1][3] * m[2][1] - m[0][3] * m[1][1] * m[2][2] + m[0][1] * m[1][3] * m[2][2] + m[0][2] * m[1][1] * m[2][3] - m[0][1] * m[1][2] * m[2][3]) * inv_det,
+
+		(m[1][3] * m[2][2] * m[3][0] - m[1][2] * m[2][3] * m[3][0] - m[1][3] * m[2][0] * m[3][2] + m[1][0] * m[2][3] * m[3][2] + m[1][2] * m[2][0] * m[3][3] - m[1][0] * m[2][2] * m[3][3]) * inv_det,
+		(m[0][2] * m[2][3] * m[3][0] - m[0][3] * m[2][2] * m[3][0] + m[0][3] * m[2][0] * m[3][2] - m[0][0] * m[2][3] * m[3][2] - m[0][2] * m[2][0] * m[3][3] + m[0][0] * m[2][2] * m[3][3]) * inv_det,
+		(m[0][3] * m[1][2] * m[3][0] - m[0][2] * m[1][3] * m[3][0] - m[0][3] * m[1][0] * m[3][2] + m[0][0] * m[1][3] * m[3][2] + m[0][2] * m[1][0] * m[3][3] - m[0][0] * m[1][2] * m[3][3]) * inv_det,
+		(m[0][2] * m[1][3] * m[2][0] - m[0][3] * m[1][2] * m[2][0] + m[0][3] * m[1][0] * m[2][2] - m[0][0] * m[1][3] * m[2][2] - m[0][2] * m[1][0] * m[2][3] + m[0][0] * m[1][2] * m[2][3]) * inv_det,
+
+		(m[1][1] * m[2][3] * m[3][0] - m[1][3] * m[2][1] * m[3][0] + m[1][3] * m[2][0] * m[3][1] - m[1][0] * m[2][3] * m[3][1] - m[1][1] * m[2][0] * m[3][3] + m[1][0] * m[2][1] * m[3][3]) * inv_det,
+		(m[0][3] * m[2][1] * m[3][0] - m[0][1] * m[2][3] * m[3][0] - m[0][3] * m[2][0] * m[3][1] + m[0][0] * m[2][3] * m[3][1] + m[0][1] * m[2][0] * m[3][3] - m[0][0] * m[2][1] * m[3][3]) * inv_det,
+		(m[0][1] * m[1][3] * m[3][0] - m[0][3] * m[1][1] * m[3][0] + m[0][3] * m[1][0] * m[3][1] - m[0][0] * m[1][3] * m[3][1] - m[0][1] * m[1][0] * m[3][3] + m[0][0] * m[1][1] * m[3][3]) * inv_det,
+		(m[0][3] * m[1][1] * m[2][0] - m[0][1] * m[1][3] * m[2][0] - m[0][3] * m[1][0] * m[2][1] + m[0][0] * m[1][3] * m[2][1] + m[0][1] * m[1][0] * m[2][3] - m[0][0] * m[1][1] * m[2][3]) * inv_det,
+
+		(m[1][2] * m[2][1] * m[3][0] - m[1][1] * m[2][2] * m[3][0] - m[1][2] * m[2][0] * m[3][1] + m[1][0] * m[2][2] * m[3][1] + m[1][1] * m[2][0] * m[3][2] - m[1][0] * m[2][1] * m[3][2]) * inv_det,
+		(m[0][1] * m[2][2] * m[3][0] - m[0][2] * m[2][1] * m[3][0] + m[0][2] * m[2][0] * m[3][1] - m[0][0] * m[2][2] * m[3][1] - m[0][1] * m[2][0] * m[3][2] + m[0][0] * m[2][1] * m[3][2]) * inv_det,
+		(m[0][2] * m[1][1] * m[3][0] - m[0][1] * m[1][2] * m[3][0] - m[0][2] * m[1][0] * m[3][1] + m[0][0] * m[1][2] * m[3][1] + m[0][1] * m[1][0] * m[3][2] - m[0][0] * m[1][1] * m[3][2]) * inv_det,
+		(m[0][1] * m[1][2] * m[2][0] - m[0][2] * m[1][1] * m[2][0] + m[0][2] * m[1][0] * m[2][1] - m[0][0] * m[1][2] * m[2][1] - m[0][1] * m[1][0] * m[2][2] + m[0][0] * m[1][1] * m[2][2]) * inv_det
+	};
 }
 
-template <class T>
-constexpr matrix<T, 4, 4> look_at(const vector<T, 3>& position, const vector<T, 3>& target, vector<T, 3> up)
+/// @private
+template <class T, std::size_t N, std::size_t M, std::size_t... I>
+constexpr inline matrix<T, N, M> componentwise_mul(const matrix<T, N, M>& a, const matrix<T, N, M>& b, std::index_sequence<I...>) noexcept
 {
-	vector<T, 3> forward = normalize(sub(target, position));
-	vector<T, 3> right = normalize(cross(forward, up));
-	up = cross(right, forward);
-
-	matrix<T, 4, 4> m = 
-		{{
-			{right[0], up[0], -forward[0], T(0)},
-			{right[1], up[1], -forward[1], T(0)},
-			{right[2], up[2], -forward[2], T(0)},
-			{T(0), T(0), T(0), T(1)}
-		}};
-	
-	return translate(m, negate(position));
+	return {(a[I] * b[I]) ...};
 }
 
-template <class T>
-constexpr matrix<T, 2, 2> mul(const matrix<T, 2, 2>& x, const matrix<T, 2, 2>& y)
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> componentwise_mul(const matrix<T, N, M>& a, const matrix<T, N, M>& b) noexcept
 {
-	return
-		{{
-			x[0] * y[0][0] + x[1] * y[0][1],
-			x[0] * y[1][0] + x[1] * y[1][1]
-		}};
+	return componentwise_mul(a, b, std::make_index_sequence<N>{});
 }
 
-template <class T>
-constexpr matrix<T, 3, 3> mul(const matrix<T, 3, 3>& x, const matrix<T, 3, 3>& y)
+/// @private
+template <class T, std::size_t N, std::size_t M, std::size_t... I>
+constexpr inline matrix<T, N, M> div(const matrix<T, N, M>& a, const matrix<T, N, M>& b, std::index_sequence<I...>) noexcept
 {
-	return
-		{{
-			x[0] * y[0][0] + x[1] * y[0][1] + x[2] * y[0][2],
-			x[0] * y[1][0] + x[1] * y[1][1] + x[2] * y[1][2],
-			x[0] * y[2][0] + x[1] * y[2][1] + x[2] * y[2][2]
-		}};
+	return {(a[I] / b[I]) ...};
 }
 
-template <class T>
-constexpr matrix<T, 4, 4> mul(const matrix<T, 4, 4>& x, const matrix<T, 4, 4>& y)
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> div(const matrix<T, N, M>& a, const matrix<T, N, M>& b) noexcept
 {
-	return
-		{{
-			x[0] * y[0][0] + x[1] * y[0][1] + x[2] * y[0][2] + x[3] * y[0][3],
-			x[0] * y[1][0] + x[1] * y[1][1] + x[2] * y[1][2] + x[3] * y[1][3],
-			x[0] * y[2][0] + x[1] * y[2][1] + x[2] * y[2][2] + x[3] * y[2][3],
-			x[0] * y[3][0] + x[1] * y[3][1] + x[2] * y[3][2] + x[3] * y[3][3]
-		}};
+	return div(a, b, std::make_index_sequence<N>{});
 }
 
 /// @private
 template <class T, std::size_t N, std::size_t M, std::size_t... I>
-constexpr inline matrix<T, N, M> mul(const matrix<T, N, M>& m, T s, std::index_sequence<I...>)
+constexpr inline matrix<T, N, M> div(const matrix<T, N, M>& a, T b, std::index_sequence<I...>) noexcept
 {
-	return {{(m[I] * s)...}};
+	return {(a[I] / b) ...};
 }
 
 template <class T, std::size_t N, std::size_t M>
-constexpr inline matrix<T, N, M> mul(const matrix<T, N, M>& m, T s)
+constexpr matrix<T, N, M> div(const matrix<T, N, M>& a, T b) noexcept
 {
-	return mul(m, s, std::make_index_sequence<N>{});
+	return div(a, b, std::make_index_sequence<N>{});
 }
 
-template <class T>
-constexpr vector<T, 2> mul(const matrix<T, 2, 2>& m, const vector<T, 2>& v)
+/// @private
+template <class T, std::size_t N, std::size_t M, std::size_t... I>
+constexpr inline matrix<T, N, M> div(T a, const matrix<T, N, M>& b, std::index_sequence<I...>) noexcept
 {
-	return
-		{
-			m[0][0] * v[0] + m[1][0] * v[1],
-			m[0][1] * v[0] + m[1][1] * v[1]
-		};
+	return {(a / b[I]) ...};
 }
 
-template <class T>
-constexpr vector<T, 3> mul(const matrix<T, 3, 3>& m, const vector<T, 3>& v)
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> div(T a, const matrix<T, N, M>& b) noexcept
 {
-	return
-		{
-			m[0][0] * v[0] + m[1][0] * v[1] + m[2][0] * v[2],
-			m[0][1] * v[0] + m[1][1] * v[1] + m[2][1] * v[2],
-			m[0][2] * v[0] + m[1][2] * v[1] + m[2][2] * v[2]
-		};
+	return div(a, b, std::make_index_sequence<N>{});
 }
 
 template <class T>
-constexpr vector<T, 4> mul(const matrix<T, 4, 4>& m, const vector<T, 4>& v)
+constexpr matrix<T, 4, 4> look_at(const vector<T, 3>& position, const vector<T, 3>& target, vector<T, 3> up)
 {
-	return
-		{
-			m[0][0] * v[0] + m[1][0] * v[1] + m[2][0] * v[2] + m[3][0] * v[3],
-			m[0][1] * v[0] + m[1][1] * v[1] + m[2][1] * v[2] + m[3][1] * v[3],
-			m[0][2] * v[0] + m[1][2] * v[1] + m[2][2] * v[2] + m[3][2] * v[3],
-			m[0][3] * v[0] + m[1][3] * v[1] + m[2][3] * v[2] + m[3][3] * v[3]
-		};
-}
+	vector<T, 3> forward = normalize(sub(target, position));
+	vector<T, 3> right = normalize(cross(forward, up));
+	up = cross(right, forward);
 
-template <class T>
-constexpr matrix<T, 4, 4> ortho(T left, T right, T bottom, T top, T z_near, T z_far)
-{
-	return
+	matrix<T, 4, 4> m = 
 		{{
-			 {T(2) / (right - left), T(0), T(0), T(0)},
-			 {T(0), T(2) / (top - bottom), T(0), T(0)},
-			 {T(0), T(0), T(-2) / (z_far - z_near), T(0)},
-			 {-((right + left) / (right - left)), -((top + bottom) / (top - bottom)), -((z_far + z_near) / (z_far - z_near)), T(1)}
+			{right[0], up[0], -forward[0], T(0)},
+			{right[1], up[1], -forward[1], T(0)},
+			{right[2], up[2], -forward[2], T(0)},
+			{T(0), T(0), T(0), T(1)}
 		}};
+	
+	return translate(m, negate(position));
 }
 
-template <class T>
-constexpr matrix<T, 4, 4> ortho_half_z(T left, T right, T bottom, T top, T z_near, T z_far)
+template <typename T, std::size_t N, std::size_t M, std::size_t P>
+constexpr matrix<T, P, M> mul(const matrix<T, N, M>& a, const matrix<T, P, N>& b) noexcept
 {
-	return
-		{{
-			 {T(2) / (right - left), T(0), T(0), T(0)},
-			 {T(0), T(2) / (top - bottom), T(0), T(0)},
-			 {T(0), T(0), T(-1) / (z_far - z_near), T(0)},
-			 {-((right + left) / (right - left)), -((top + bottom) / (top - bottom)), -z_near / (z_far - z_near), T(1)}
-		}};
+	matrix<T, P, M> c = matrix<T, P, M>::zero();
+	
+	for (std::size_t i = 0; i < P; ++i)
+	{
+		for (std::size_t j = 0; j < M; ++j)
+		{
+			for (std::size_t k = 0; k < N; ++k)
+			{
+				c[i][j] += a[k][j] * b[i][k];
+			}
+		}
+	}
+	
+	return c;
 }
 
-template <class T>
-constexpr matrix<T, 2, 2> outer_product(const vector<T, 2>& c, const vector<T, 2>& r)
+/// @private
+template <class T, std::size_t N, std::size_t M, std::size_t... I>
+constexpr inline matrix<T, N, M> mul(const matrix<T, N, M>& a, T b, std::index_sequence<I...>) noexcept
 {
-	return
-		{{
-			 {c[0] * r[0], c[1] * r[0]},
-			 {c[0] * r[1], c[1] * r[1]}
-		}};
+	return {(a[I] * b) ...};
 }
 
-template <class T>
-constexpr matrix<T, 3, 3> outer_product(const vector<T, 3>& c, const vector<T, 3>& r)
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> mul(const matrix<T, N, M>& a, T b) noexcept
 {
-	return
-		{{
-			 {c[0] * r[0], c[1] * r[0], c[2] * r[0]},
-			 {c[0] * r[1], c[1] * r[1], c[2] * r[1]},
-			 {c[0] * r[2], c[1] * r[2], c[2] * r[2]}
-		}};
+	return mul(a, b, std::make_index_sequence<N>{});
 }
 
-template <class T>
-constexpr matrix<T, 4, 4> outer_product(const vector<T, 4>& c, const vector<T, 4> r)
+/// @private
+template <typename T, std::size_t N, std::size_t M, std::size_t... I>
+constexpr inline typename matrix<T, N, M>::column_vector_type mul(const matrix<T, N, M>& a, const typename matrix<T, N, M>::row_vector_type& b, std::index_sequence<I...>) noexcept
 {
-	return
-		{{
-			 {c[0] * r[0], c[1] * r[0], c[2] * r[0], c[3] * r[0]},
-			 {c[0] * r[1], c[1] * r[1], c[2] * r[1], c[3] * r[1]},
-			 {c[0] * r[2], c[1] * r[2], c[2] * r[2], c[3] * r[2]},
-			 {c[0] * r[3], c[1] * r[3], c[2] * r[3], c[3] * r[3]}
-		}};
+	return ((a[I] * b[I]) + ...);
 }
 
-template <class T>
-matrix<T, 4, 4> perspective(T vertical_fov, T aspect_ratio, T z_near, T z_far)
+template <typename T, std::size_t N, std::size_t M>
+constexpr typename matrix<T, N, M>::column_vector_type mul(const matrix<T, N, M>& a, const typename matrix<T, N, M>::row_vector_type& b) noexcept
 {
-	T half_fov = vertical_fov * T(0.5);
-	T f = std::cos(half_fov) / std::sin(half_fov);
-
-	return
-		{{
-			 {f / aspect_ratio, T(0), T(0), T(0)},
-			 {T(0), f, T(0), T(0)},
-			 {T(0), T(0), (z_far + z_near) / (z_near - z_far), T(-1)},
-			 {T(0), T(0), (T(2) * z_far * z_near) / (z_near - z_far), T(0)}
-		}};
+	return mul(a, b, std::make_index_sequence<N>{});
 }
 
-template <class T>
-matrix<T, 4, 4> perspective_half_z(T vertical_fov, T aspect_ratio, T z_near, T z_far)
+/// @private
+template <typename T, std::size_t N, std::size_t M, std::size_t... I>
+constexpr inline typename matrix<T, N, M>::row_vector_type mul(const typename matrix<T, N, M>::column_vector_type& a, const matrix<T, N, M>& b, std::index_sequence<I...>) noexcept
 {
-	T half_fov = vertical_fov * T(0.5);
-	T f = std::cos(half_fov) / std::sin(half_fov);
-
-	return
-		{{
-			 {f / aspect_ratio, T(0), T(0), T(0)},
-			 {T(0), f, T(0), T(0)},
-			 {T(0), T(0), z_far / (z_near - z_far), T(-1)},
-			 {T(0), T(0), -(z_far * z_near) / (z_far - z_near), T(0)}
-		}};
+	return {dot(a, b[I]) ...};
 }
 
-template <std::size_t N1, std::size_t M1, class T, std::size_t N0, std::size_t M0>
-constexpr matrix<T, N1, M1> resize(const matrix<T, N0, M0>& m)
+template <typename T, std::size_t N, std::size_t M>
+constexpr typename matrix<T, N, M>::row_vector_type mul(const typename matrix<T, N, M>::column_vector_type& a, const matrix<T, N, M>& b) noexcept
 {
-	matrix<T, N1, M1> resized;
-
-	for (std::size_t i = 0; i < N1; ++i)
-	{
-		for (std::size_t j = 0; j < M1; ++j)
-		{
-			resized[i][j] = (i < N0 && j < M0) ? m[i][j] : ((i == j) ? T(1) : T(0));
-		}
-	}
-
-	return resized;
+	return mul(a, b, std::make_index_sequence<N>{});
 }
 
 template <class T>
@@ -1108,37 +903,43 @@ constexpr matrix<T, 4, 4> scale(const matrix<T, 4, 4>& m, const vector<T, 3>& v)
 		}});
 }
 
-template <class T>
-constexpr matrix<T, 2, 2> sub(const matrix<T, 2, 2>& x, const matrix<T, 2, 2>& y)
+/// @private
+template <class T, std::size_t N, std::size_t M, std::size_t... I>
+constexpr inline matrix<T, N, M> sub(const matrix<T, N, M>& a, const matrix<T, N, M>& b, std::index_sequence<I...>) noexcept
 {
-	return
-		{{
-			x[0] - y[0],
-			x[1] - y[1]
-		}};
+	return {(a[I] - b[I]) ...};
 }
 
-template <class T>
-constexpr matrix<T, 3, 3> sub(const matrix<T, 3, 3>& x, const matrix<T, 3, 3>& y)
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> sub(const matrix<T, N, M>& a, const matrix<T, N, M>& b) noexcept
 {
-	return
-		{{
-			x[0] - y[0],
-			x[1] - y[1],
-			x[2] - y[2]
-		}};
+	return sub(a, b, std::make_index_sequence<N>{});
 }
 
-template <class T>
-constexpr matrix<T, 4, 4> sub(const matrix<T, 4, 4>& x, const matrix<T, 4, 4>& y)
+/// @private
+template <class T, std::size_t N, std::size_t M, std::size_t... I>
+constexpr inline matrix<T, N, M> sub(const matrix<T, N, M>& a, T b, std::index_sequence<I...>) noexcept
 {
-	return
-		{{
-			x[0] - y[0],
-			x[1] - y[1],
-			x[2] - y[2],
-			x[3] - y[3]
-		}};
+	return {(a[I] - b) ...};
+}
+
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> sub(const matrix<T, N, M>& a, T b) noexcept
+{
+	return sub(a, b, std::make_index_sequence<N>{});
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t M, std::size_t... I>
+constexpr inline matrix<T, N, M> sub(T a, const matrix<T, N, M>& b, std::index_sequence<I...>) noexcept
+{
+	return {(a - b[I]) ...};
+}
+
+template <class T, std::size_t N, std::size_t M>
+constexpr matrix<T, N, M> sub(T a, const matrix<T, N, M>& b) noexcept
+{
+	return sub(a, b, std::make_index_sequence<N>{});
 }
 
 template <class T>
@@ -1153,161 +954,213 @@ constexpr matrix<T, 4, 4> translate(const matrix<T, 4, 4>& m, const vector<T, 3>
 		}});
 }
 
-template <class T>
-constexpr matrix<T, 2, 2> transpose(const matrix<T, 2, 2>& m)
+/// @private
+template <typename T, std::size_t N, std::size_t M, std::size_t... I>
+constexpr inline typename matrix<T, M, N>::column_vector_type transpose_column(const matrix<T, N, M>& m, std::size_t i, std::index_sequence<I...>) noexcept
 {
-
-	return
-		{{
-			{
-				m[0][0], m[1][0]
-			},
-
-			{
-				m[0][1], m[1][1]
-			}
-		}};
+	return {m[I][i] ...};
 }
 
-template <class T>
-constexpr matrix<T, 3, 3> transpose(const matrix<T, 3, 3>& m)
+/// @private
+template <typename T, std::size_t N, std::size_t M, std::size_t... I>
+constexpr inline matrix<T, M, N> transpose(const matrix<T, N, M>& m, std::index_sequence<I...>) noexcept
 {
-
-	return
-		{{
-			{
-				m[0][0], m[1][0], m[2][0]
-			},
-
-			{
-				m[0][1], m[1][1], m[2][1]
-			},
-
-			{
-				m[0][2], m[1][2], m[2][2]
-			}
-		}};
+	return {transpose_column(m, I, std::make_index_sequence<N>{}) ...};
 }
 
-template <class T>
-constexpr matrix<T, 4, 4> transpose(const matrix<T, 4, 4>& m)
+template <typename T, std::size_t N, std::size_t M>
+constexpr matrix<T, M, N> transpose(const matrix<T, N, M>& m) noexcept
 {
+	return transpose(m, std::make_index_sequence<M>{});
+}
 
-	return
-		{{
-			{
-				m[0][0], m[1][0], m[2][0], m[3][0]
-			},
-
-			{
-				m[0][1], m[1][1], m[2][1], m[3][1]
-			},
-
-			{
-				m[0][2], m[1][2], m[2][2], m[3][2]
-			},
+} // namespace math
 
-			{
-				m[0][3], m[1][3], m[2][3], m[3][3]
-			}
-		}};
+/// @copydoc math::add(const math::matrix<T, N, M>&, const math::matrix<T, N, M>&)
+template <class T, std::size_t N, std::size_t M>
+constexpr inline math::matrix<T, N, M> operator+(const math::matrix<T, N, M>& a, const math::matrix<T, N, M>& b) noexcept
+{
+	return math::add(a, b);
 }
 
-/// @private
-template <class T2, class T1, std::size_t N, std::size_t M, std::size_t... I>
-constexpr inline matrix<T2, N, M> type_cast(const matrix<T1, N, M>& m, std::index_sequence<I...>)
+/// @copydoc math::add(const math::matrix<T, N, M>&, T)
+/// @{
+template <class T, std::size_t N, std::size_t M>
+constexpr inline math::matrix<T, N, M> operator+(const math::matrix<T, N, M>& a, T b) noexcept
 {
-	return {type_cast<T2>(m[I])...};
+	return math::add(a, b);
 }
-
-template <class T2, class T1, std::size_t N, std::size_t M>
-constexpr matrix<T2, N, M> type_cast(const matrix<T1, N, M>& m)
+template <class T, std::size_t N, std::size_t M>
+constexpr inline math::matrix<T, N, M> operator+(T a, const math::matrix<T, N, M>& b) noexcept
 {
-	return type_cast<T2>(m, std::make_index_sequence<N>{}); 
+	return math::add(b, a);
 }
+/// @}
 
-} // namespace math
+/// @copydoc math::div(const math::matrix<T, N, M>&, const math::matrix<T, N, M>&)
+template <class T, std::size_t N, std::size_t M>
+constexpr inline math::matrix<T, N, M> operator/(const math::matrix<T, N, M>& a, const math::matrix<T, N, M>& b) noexcept
+{
+	return math::div(a, b);
+}
 
-/// @copydoc math::add(const math::matrix<T, 2, 2>&, const math::matrix<T, 2, 2>&)
+/// @copydoc math::div(const math::matrix<T, N, M>&, T)
 template <class T, std::size_t N, std::size_t M>
-constexpr math::matrix<T, N, M> operator+(const math::matrix<T, N, M>& x, const math::matrix<T, N, M>& y);
+constexpr inline math::matrix<T, N, M> operator/(const math::matrix<T, N, M>& a, T b) noexcept
+{
+	return math::div(a, b);
+}
 
-/// @copydoc math::mul(const math::matrix<T, 2, 2>&, const math::matrix<T, 2, 2>&)
+/// @copydoc math::div(T, const math::matrix<T, N, M>&)
 template <class T, std::size_t N, std::size_t M>
-constexpr math::matrix<T, N, M> operator*(const math::matrix<T, N, M>& x, const math::matrix<T, N, M>& y);
+constexpr inline math::matrix<T, N, M> operator/(T a, const math::matrix<T, N, M>& b) noexcept
+{
+	return math::div(a, b);
+}
+
+/// @copydoc math::mul(const math::matrix<T, N, M>&, const math::matrix<T, P, N>&)
+template <typename T, std::size_t N, std::size_t M, std::size_t P>
+constexpr inline math::matrix<T, P, M> operator*(const math::matrix<T, N, M>& a, const math::matrix<T, P, N>& b) noexcept
+{
+	return math::mul(a, b);
+}
 
 /// @copydoc math::mul(const math::matrix<T, N, M>&, T)
+/// @{
 template <class T, std::size_t N, std::size_t M>
-constexpr math::matrix<T, N, M> operator*(const math::matrix<T, N, M>& m, T s);
+constexpr inline math::matrix<T, N, M> operator*(const math::matrix<T, N, M>& a, T b) noexcept
+{
+	return math::mul(a, b);
+}
+template <class T, std::size_t N, std::size_t M>
+constexpr inline math::matrix<T, N, M> operator*(T a, const math::matrix<T, N, M>& b) noexcept
+{
+	return math::mul(b, a);
+}
+/// @}
 
-/// @copydoc math::mul(const math::matrix<T, N, M>&, T)
+/// @copydoc math::mul(const math::matrix<T, N, M>&, const typename math::matrix<T, N, M>::row_vector_type&)
+template <typename T, std::size_t N, std::size_t M>
+constexpr inline typename math::matrix<T, N, M>::column_vector_type operator*(const math::matrix<T, N, M>& a, const typename math::matrix<T, N, M>::row_vector_type& b) noexcept
+{
+	return math::mul(a, b);
+}
+
+/// @copydoc math::mul(const typename math::matrix<T, N, M>::column_vector_type&, const math::matrix<T, N, M>&)
+template <typename T, std::size_t N, std::size_t M>
+constexpr inline typename math::matrix<T, N, M>::row_vector_type operator*(const typename math::matrix<T, N, M>::column_vector_type& a, const math::matrix<T, N, M>& b) noexcept
+{
+	return math::mul(a, b);
+}
+
+/// @copydoc math::sub(const math::matrix<T, N, M>&, const math::matrix<T, N, M>&)
 template <class T, std::size_t N, std::size_t M>
-constexpr math::matrix<T, N, M> operator*(T s, const math::matrix<T, N, M>& m);
+constexpr inline math::matrix<T, N, M> operator-(const math::matrix<T, N, M>& a, const math::matrix<T, N, M>& b) noexcept
+{
+	return math::sub(a, b);
+}
 
-/// @copydoc math::mul(const math::matrix<T, 2, 2>&, const math::vector<T, 2>&)
-template <class T, std::size_t N>
-constexpr math::vector<T, N> operator*(const math::matrix<T, N, N>& m, const math::vector<T, N>& v);
+/// @copydoc math::sub(const math::matrix<T, N, M>&, T)
+template <class T, std::size_t N, std::size_t M>
+constexpr inline math::matrix<T, N, M> operator-(const math::matrix<T, N, M>& a, T b) noexcept
+{
+	return math::sub(a, b);
+}
 
-/// @copydoc math::sub(const math::matrix<T, 2, 2>&, const math::matrix<T, 2, 2>&)
+/// @copydoc math::sub(T, const math::matrix<T, N, M>&)
 template <class T, std::size_t N, std::size_t M>
-constexpr math::matrix<T, N, M> operator-(const math::matrix<T, N, M>& x, const math::matrix<T, N, M>& y);
+constexpr inline math::matrix<T, N, M> operator-(T a, const math::matrix<T, N, M>& b) noexcept
+{
+	return math::sub(a, b);
+}
 
 /**
- * Writes the elements of a matrix to an output stream, with each element delimeted by a space.
+ * Adds two values and stores the result in the first value.
  *
- * @param os Output stream.
- * @param m Matrix.
- * @return Output stream.
+ * @param a First value.
+ * @param b Second value.
+ * @return Reference to the first value.
  */
+/// @{
 template <class T, std::size_t N, std::size_t M>
-std::ostream& operator<<(std::ostream& os, const math::matrix<T, N, M>& m);
+constexpr inline math::matrix<T, N, M>& operator+=(math::matrix<T, N, M>& a, const math::matrix<T, N, M>& b) noexcept
+{
+	return (a = a + b);
+}
+template <class T, std::size_t N, std::size_t M>
+constexpr inline math::matrix<T, N, M>& operator+=(math::matrix<T, N, M>& a, T b) noexcept
+{
+	return (a = a + b);
+}
+/// @}
 
 /**
- * Reads the elements of a matrix from an input stream, with each element delimeted by a space.
+ * Subtracts the first value by the second value and stores the result in the first value.
  *
- * @param is Input stream.
- * @param m Matrix.
- * @return Input stream.
+ * @param a First value.
+ * @param b Second value.
+ * @return Reference to the first value.
  */
+/// @{
 template <class T, std::size_t N, std::size_t M>
-std::istream& operator>>(std::istream& is, math::matrix<T, N, M>& m);
-
-template <class T, std::size_t N, std::size_t M>
-constexpr inline math::matrix<T, N, M> operator+(const math::matrix<T, N, M>& x, const math::matrix<T, N, M>& y)
+constexpr inline math::matrix<T, N, M>& operator-=(math::matrix<T, N, M>& a, const math::matrix<T, N, M>& b) noexcept
 {
-	return math::add(x, y);
+	return (a = a - b);
 }
-
 template <class T, std::size_t N, std::size_t M>
-constexpr inline math::matrix<T, N, M> operator*(const math::matrix<T, N, M>& x, const math::matrix<T, N, M>& y)
+constexpr inline math::matrix<T, N, M>& operator-=(math::matrix<T, N, M>& a, T b) noexcept
 {
-	return math::mul(x, y);
+	return (a = a - b);
 }
+/// @}
 
-template <class T, std::size_t N, std::size_t M>
-constexpr inline math::matrix<T, N, M> operator*(const math::matrix<T, N, M>& m, T s)
+/**
+ * Multiplies two values and stores the result in the first value.
+ *
+ * @param a First value.
+ * @param b Second value.
+ * @return Reference to the first value.
+ */
+/// @{
+template <class T, std::size_t N>
+constexpr inline math::matrix<T, N, N>& operator*=(math::matrix<T, N, N>& a, math::matrix<T, N, N>& b) noexcept
 {
-	return math::mul(m, s);
+	return (a = a * b);
 }
-
 template <class T, std::size_t N, std::size_t M>
-constexpr inline math::matrix<T, N, M> operator*(T s, const math::matrix<T, N, M>& m)
+constexpr inline math::matrix<T, N, M>& operator*=(math::matrix<T, N, M>& a, T b) noexcept
 {
-	return math::mul(m, s);
+	return (a = a * b);
 }
+/// @}
 
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N> operator*(const math::matrix<T, N, N>& m, const math::vector<T, N>& v)
+/**
+ * Divides the first value by the second value and stores the result in the first value.
+ *
+ * @param a First value.
+ * @param b Second value.
+ * @return Reference to the first value.
+ */
+/// @{
+template <class T, std::size_t N, std::size_t M>
+constexpr inline math::matrix<T, N, M>& operator/=(math::matrix<T, N, M>& a, const math::matrix<T, N, M>& b) noexcept
 {
-	return math::mul(m, v);
+	return (a = a / b);
 }
-
 template <class T, std::size_t N, std::size_t M>
-constexpr inline math::matrix<T, N, M> operator-(const math::matrix<T, N, M>& x, const math::matrix<T, N, M>& y)
+constexpr inline math::matrix<T, N, M>& operator/=(math::matrix<T, N, M>& a, T b) noexcept
 {
-	return math::sub(x, y);
+	return (a = a / b);
 }
+/// @}
 
+/**
+ * Writes the elements of a matrix to an output stream, with each element delimeted by a space.
+ *
+ * @param os Output stream.
+ * @param m Matrix.
+ * @return Output stream.
+ */
 template <class T, std::size_t N, std::size_t M>
 std::ostream& operator<<(std::ostream& os, const math::matrix<T, N, M>& m)
 {
@@ -1321,6 +1174,13 @@ std::ostream& operator<<(std::ostream& os, const math::matrix<T, N, M>& m)
 	return os;
 }
 
+/**
+ * Reads the elements of a matrix from an input stream, with each element delimeted by a space.
+ *
+ * @param is Input stream.
+ * @param m Matrix.
+ * @return Input stream.
+ */
 template <class T, std::size_t N, std::size_t M>
 std::istream& operator>>(std::istream& is, math::matrix<T, N, M>& m)
 {
diff --git a/src/math/projection.hpp b/src/math/projection.hpp
index a72ab69..5e9446b 100644
--- a/src/math/projection.hpp
+++ b/src/math/projection.hpp
@@ -20,6 +20,7 @@
 #ifndef ANTKEEPER_MATH_PROJECTION_HPP
 #define ANTKEEPER_MATH_PROJECTION_HPP
 
+#include "math/matrix.hpp"
 #include <cmath>
 
 namespace math {
@@ -56,6 +57,104 @@ T vertical_fov(T h, T r)
 	return T{2} * std::atan(std::tan(h * T{0.5}) / r);
 }
 
+/**
+ * Creates an orthographic projection matrix which will transform the near and far clipping planes to `[-1, 1]`, respectively.
+ *
+ * @param left Signed distance to the left clipping plane.
+ * @param right Signed distance to the right clipping plane.
+ * @param bottom Signed distance to the bottom clipping plane.
+ * @param top Signed distance to the top clipping plane.
+ * @param z_near Signed distance to the near clipping plane.
+ * @param z_far Signed distance to the far clipping plane.
+ *
+ * @return Orthographic projection matrix.
+ */
+template <class T>
+constexpr matrix<T, 4, 4> ortho(T left, T right, T bottom, T top, T z_near, T z_far) noexcept
+{
+	return
+		{{
+			 {T(2) / (right - left), T(0), T(0), T(0)},
+			 {T(0), T(2) / (top - bottom), T(0), T(0)},
+			 {T(0), T(0), T(-2) / (z_far - z_near), T(0)},
+			 {-((right + left) / (right - left)), -((top + bottom) / (top - bottom)), -((z_far + z_near) / (z_far - z_near)), T(1)}
+		}};
+}
+
+/**
+ * Creates an orthographic projection matrix which will transform the near and far clipping planes to `[0, 1]`, respectively.
+ *
+ * @param left Signed distance to the left clipping plane.
+ * @param right Signed distance to the right clipping plane.
+ * @param bottom Signed distance to the bottom clipping plane.
+ * @param top Signed distance to the top clipping plane.
+ * @param z_near Signed distance to the near clipping plane.
+ * @param z_far Signed distance to the far clipping plane.
+ *
+ * @return Orthographic projection matrix.
+ */
+template <class T>
+constexpr matrix<T, 4, 4> ortho_half_z(T left, T right, T bottom, T top, T z_near, T z_far) noexcept
+{
+	return
+		{{
+			 {T(2) / (right - left), T(0), T(0), T(0)},
+			 {T(0), T(2) / (top - bottom), T(0), T(0)},
+			 {T(0), T(0), T(-1) / (z_far - z_near), T(0)},
+			 {-((right + left) / (right - left)), -((top + bottom) / (top - bottom)), -z_near / (z_far - z_near), T(1)}
+		}};
+}
+
+/**
+ * Creates a perspective projection matrix which will transform the near and far clipping planes to `[-1, 1]`, respectively.
+ *
+ * @param vertical_fov Vertical field of view angle, in radians.
+ * @param aspect_ratio Aspect ratio which determines the horizontal field of view.
+ * @param z_near Distance to the near clipping plane.
+ * @param z_far Distance to the far clipping plane.
+ *
+ * @return Perspective projection matrix.
+ */
+template <class T>
+matrix<T, 4, 4> perspective(T vertical_fov, T aspect_ratio, T z_near, T z_far)
+{
+	T half_fov = vertical_fov * T(0.5);
+	T f = std::cos(half_fov) / std::sin(half_fov);
+
+	return
+		{{
+			 {f / aspect_ratio, T(0), T(0), T(0)},
+			 {T(0), f, T(0), T(0)},
+			 {T(0), T(0), (z_far + z_near) / (z_near - z_far), T(-1)},
+			 {T(0), T(0), (T(2) * z_far * z_near) / (z_near - z_far), T(0)}
+		}};
+}
+
+/**
+ * Creates a perspective projection matrix which will transform the near and far clipping planes to `[0, 1]`, respectively.
+ *
+ * @param vertical_fov Vertical field of view angle, in radians.
+ * @param aspect_ratio Aspect ratio which determines the horizontal field of view.
+ * @param z_near Distance to the near clipping plane.
+ * @param z_far Distance to the far clipping plane.
+ *
+ * @return Perspective projection matrix.
+ */
+template <class T>
+matrix<T, 4, 4> perspective_half_z(T vertical_fov, T aspect_ratio, T z_near, T z_far)
+{
+	T half_fov = vertical_fov * T(0.5);
+	T f = std::cos(half_fov) / std::sin(half_fov);
+
+	return
+		{{
+			 {f / aspect_ratio, T(0), T(0), T(0)},
+			 {T(0), f, T(0), T(0)},
+			 {T(0), T(0), z_far / (z_near - z_far), T(-1)},
+			 {T(0), T(0), -(z_far * z_near) / (z_far - z_near), T(0)}
+		}};
+}
+
 } // namespace math
 
 #endif // ANTKEEPER_MATH_PROJECTION_HPP
diff --git a/src/math/se3.hpp b/src/math/se3.hpp
index 91693a5..7b23c2d 100644
--- a/src/math/se3.hpp
+++ b/src/math/se3.hpp
@@ -95,7 +95,7 @@ se3<T> se3<T>::inverse() const
 template <class T>
 typename se3<T>::matrix_type se3<T>::matrix() const
 {
-	matrix_type m = math::resize<4, 4>(math::matrix_cast<T>(r));
+	matrix_type m = math::matrix<T, 4, 4>(math::matrix_cast<T>(r));
 	
 	m[3].x() = t.x();
 	m[3].y() = t.y();
diff --git a/src/math/transform-functions.hpp b/src/math/transform-functions.hpp
index 18fcfbf..639f770 100644
--- a/src/math/transform-functions.hpp
+++ b/src/math/transform-functions.hpp
@@ -75,7 +75,7 @@ transform<T> inverse(const transform<T>& t)
 template <class T>
 matrix<T, 4, 4> matrix_cast(const transform<T>& t)
 {
-	matrix<T, 4, 4> transformation = resize<4, 4>(matrix_cast(t.rotation));
+	matrix<T, 4, 4> transformation = matrix<T, 4, 4>(matrix_cast(t.rotation));
 	transformation[3] = {t.translation[0], t.translation[1], t.translation[2], T(1)};
 	return scale(transformation, t.scale);
 }
diff --git a/src/math/vector.hpp b/src/math/vector.hpp
index 67ee3ca..62fd7a2 100644
--- a/src/math/vector.hpp
+++ b/src/math/vector.hpp
@@ -308,9 +308,9 @@ constexpr vector<T, N> abs(const vector<T, N>& x);
  */
 /// @{
 template <class T, std::size_t N>
-constexpr vector<T, N> add(const vector<T, N>& x, const vector<T, N>& y);
+constexpr vector<T, N> add(const vector<T, N>& x, const vector<T, N>& y) noexcept;
 template <class T, std::size_t N>
-constexpr vector<T, N> add(const vector<T, N>& x, T y);
+constexpr vector<T, N> add(const vector<T, N>& x, T y) noexcept;
 /// @}
 
 /**
@@ -320,7 +320,7 @@ constexpr vector<T, N> add(const vector<T, N>& x, T y);
  * @return `true` if all elements are `true`, `false` otherwise.
  */
 template <std::size_t N>
-constexpr bool all(const vector<bool, N>& x);
+constexpr bool all(const vector<bool, N>& x) noexcept;
 
 /**
  * Checks if any elements of a boolean vector are `true`.
@@ -329,15 +329,7 @@ constexpr bool all(const vector<bool, N>& x);
  * @return `true` if any elements are `true`, `false` otherwise.
  */
 template <std::size_t N>
-constexpr bool any(const vector<bool, N>& x);
-
-/**
- * Reinterprets data as an *n*-dimensional vector of type `T`.
- *
- * @param data Data to reinterpret.
- */
-template <std::size_t N, typename T>
-constexpr vector<T, N>& as_vector(T& data);
+constexpr bool any(const vector<bool, N>& x) noexcept;
 
 /**
  * Performs a element-wise ceil operation.
@@ -352,15 +344,15 @@ constexpr vector<T, N> ceil(const vector<T, N>& x);
  * Clamps the values of a vector's elements.
  *
  * @param x Vector to clamp.
- * @param min_val Minimum value.
- * @param max_val Maximum value.
+ * @param min Minimum value.
+ * @param max Maximum value.
  * @return Clamped vector.
  */
 /// @{
 template <class T, std::size_t N>
-constexpr vector<T, N> clamp(const vector<T, N>& x, const vector<T, N>& min_val, const vector<T, N>& max_val);
+constexpr vector<T, N> clamp(const vector<T, N>& x, const vector<T, N>& min, const vector<T, N>& max);
 template <class T, std::size_t N>
-constexpr vector<T, N> clamp(const vector<T, N>& x, T min_val, T max_val);
+constexpr vector<T, N> clamp(const vector<T, N>& x, T min, T max);
 /// @}
 
 /**
@@ -381,7 +373,7 @@ vector<T, N> clamp_length(const vector<T, N>& x, T max_length);
  * @return Cross product of the two vectors.
  */
 template <class T>
-constexpr vector<T, 3> cross(const vector<T, 3>& x, const vector<T, 3>& y);
+constexpr vector<T, 3> cross(const vector<T, 3>& x, const vector<T, 3>& y) noexcept;
 
 /**
  * Calculates the distance between two points.
@@ -401,7 +393,7 @@ T distance(const vector<T, N>& p0, const vector<T, N>& p1);
  * @return Squared distance between the two points.
  */
 template <class T, std::size_t N>
-constexpr T distance_squared(const vector<T, N>& p0, const vector<T, N>& p1);
+constexpr T distance_squared(const vector<T, N>& p0, const vector<T, N>& p1) noexcept;
 
 /**
  * Divides a vector by a value.
@@ -412,9 +404,11 @@ constexpr T distance_squared(const vector<T, N>& p0, const vector<T, N>& p1);
  */
 /// @{
 template <class T, std::size_t N>
-constexpr vector<T, N> div(const vector<T, N>& x, const vector<T, N>& y);
+constexpr vector<T, N> div(const vector<T, N>& x, const vector<T, N>& y) noexcept;
 template <class T, std::size_t N>
-constexpr vector<T, N> div(const vector<T, N>& x, T y);
+constexpr vector<T, N> div(const vector<T, N>& x, T y) noexcept;
+template <class T, std::size_t N>
+constexpr vector<T, N> div(T x, const vector<T, N>& y) noexcept;
 /// @}
 
 /**
@@ -425,7 +419,7 @@ constexpr vector<T, N> div(const vector<T, N>& x, T y);
  * @return Dot product of the two vectors.
  */
 template <class T, std::size_t N>
-constexpr T dot(const vector<T, N>& x, const vector<T, N>& y);
+constexpr T dot(const vector<T, N>& x, const vector<T, N>& y) noexcept;
 
 /**
  * Compares two vectors for equality
@@ -435,7 +429,7 @@ constexpr T dot(const vector<T, N>& x, const vector<T, N>& y);
  * @return Boolean vector containing the result of the element comparisons.
  */
 template <class T, std::size_t N>
-constexpr vector<bool, N> equal(const vector<T, N>& x, const vector<T, N>& y);
+constexpr vector<bool, N> equal(const vector<T, N>& x, const vector<T, N>& y) noexcept;
 
 /**
  * Performs a element-wise floor operation.
@@ -479,7 +473,7 @@ constexpr vector<T, N> fract(const vector<T, N>& x);
  * @return Boolean vector containing the result of the element comparisons.
  */
 template <class T, std::size_t N>
-constexpr vector<bool, N> greater_than(const vector<T, N>& x, const vector<T, N>& y);
+constexpr vector<bool, N> greater_than(const vector<T, N>& x, const vector<T, N>& y) noexcept;
 
 /**
  * Performs a element-wise greater-than or equal-to comparison of two vectors.
@@ -489,7 +483,7 @@ constexpr vector<bool, N> greater_than(const vector<T, N>& x, const vector<T, N>
  * @return Boolean vector containing the result of the element comparisons.
  */
 template <class T, std::size_t N>
-constexpr vector<bool, N> greater_than_equal(const vector<T, N>& x, const vector<T, N>& y);
+constexpr vector<bool, N> greater_than_equal(const vector<T, N>& x, const vector<T, N>& y) noexcept;
 
 /**
  * Calculates the length of a vector.
@@ -507,7 +501,7 @@ T length(const vector<T, N>& x);
  * @return Squared length of the vector.
  */
 template <class T, std::size_t N>
-constexpr T length_squared(const vector<T, N>& x);
+constexpr T length_squared(const vector<T, N>& x) noexcept;
 
 /**
  * Performs a element-wise less-than comparison of two vectors.
@@ -517,7 +511,7 @@ constexpr T length_squared(const vector<T, N>& x);
  * @return Boolean vector containing the result of the element comparisons.
  */
 template <class T, std::size_t N>
-constexpr vector<bool, N> less_than(const vector<T, N>& x, const vector<T, N>& y);
+constexpr vector<bool, N> less_than(const vector<T, N>& x, const vector<T, N>& y) noexcept;
 
 /**
  * Performs a element-wise less-than or equal-to comparison of two vectors.
@@ -527,7 +521,7 @@ constexpr vector<bool, N> less_than(const vector<T, N>& x, const vector<T, N>& y
  * @return Boolean vector containing the result of the element comparisons.
  */
 template <class T, std::size_t N>
-constexpr vector<bool, N> less_than_equal(const vector<T, N>& x, const vector<T, N>& y);
+constexpr vector<bool, N> less_than_equal(const vector<T, N>& x, const vector<T, N>& y) noexcept;
 
 /**
  * Returns a vector containing the maximum elements of two vectors.
@@ -594,9 +588,9 @@ constexpr vector<T, N> mod(const vector<T, N>& x, T y);
  */
 /// @{
 template <class T, std::size_t N>
-constexpr vector<T, N> mul(const vector<T, N>& x, const vector<T, N>& y);
+constexpr vector<T, N> mul(const vector<T, N>& x, const vector<T, N>& y) noexcept;
 template <class T, std::size_t N>
-constexpr vector<T, N> mul(const vector<T, N>& x, T y);
+constexpr vector<T, N> mul(const vector<T, N>& x, T y) noexcept;
 /// @}
 
 /**
@@ -606,7 +600,7 @@ constexpr vector<T, N> mul(const vector<T, N>& x, T y);
  * @return Negated vector.
  */
 template <class T, std::size_t N>
-constexpr vector<T, N> negate(const vector<T, N>& x);
+constexpr vector<T, N> negate(const vector<T, N>& x) noexcept;
 
 /**
  * Calculates the unit vector in the same direction as the original vector.
@@ -624,7 +618,7 @@ vector<T, N> normalize(const vector<T, N>& x);
  * @return Logically inverted vector.
  */
 template <class T, std::size_t N>
-constexpr vector<bool, N> not(const vector<T, N>& x);
+constexpr vector<bool, N> not(const vector<T, N>& x) noexcept;
 
 /**
  * Compares two vectors for inequality
@@ -634,7 +628,7 @@ constexpr vector<bool, N> not(const vector<T, N>& x);
  * @return Boolean vector containing the result of the element comparisons.
  */
 template <class T, std::size_t N>
-constexpr vector<bool, N> not_equal(const vector<T, N>& x, const vector<T, N>& y);
+constexpr vector<bool, N> not_equal(const vector<T, N>& x, const vector<T, N>& y) noexcept;
 
 /**
  * Raises each element to a power.
@@ -651,18 +645,6 @@ template <class T, std::size_t N>
 vector<T, N> pow(const vector<T, N>& x, T y);
 /// @}
 
-/**
- * Resizes a vector. Any new elements will be set to `0`.
- *
- * @tparam M Target number of elements.
- *
- * @param x Vector to resize.
- *
- * @return Resized vector.
- */
-template <std::size_t M, class T, std::size_t N>
-constexpr vector<T, M> resize(const vector<T, N>& x);
-
 /**
  * Performs a element-wise round operation.
  *
@@ -700,9 +682,11 @@ vector<T, N> sqrt(const vector<T, N>& x);
  */
 /// @{
 template <class T, std::size_t N>
-constexpr vector<T, N> sub(const vector<T, N>& x, const vector<T, N>& y);
+constexpr vector<T, N> sub(const vector<T, N>& x, const vector<T, N>& y) noexcept;
+template <class T, std::size_t N>
+constexpr vector<T, N> sub(const vector<T, N>& x, T y) noexcept;
 template <class T, std::size_t N>
-constexpr vector<T, N> sub(const vector<T, N>& x, T y);
+constexpr vector<T, N> sub(T x, const vector<T, N>& y) noexcept;
 /// @}
 
 /**
@@ -712,7 +696,7 @@ constexpr vector<T, N> sub(const vector<T, N>& x, T y);
  * @return Sum of the vector's elements.
  */
 template <class T, std::size_t N>
-constexpr T sum(const vector<T, N>& x);
+constexpr T sum(const vector<T, N>& x) noexcept;
 
 /**
  * Makes an *m*-dimensional vector by rearranging and/or duplicating elements of an *n*-dimensional vector.
@@ -726,7 +710,7 @@ constexpr T sum(const vector<T, N>& x);
  * @return Vector containing elements from @p x in the order specified by @p Indices. The size of the returned vector is equivalent to the number of indices in @p Indices.
  */
 template <std::size_t... Indices, class T, std::size_t N>
-constexpr vector<T, sizeof...(Indices)> swizzle(const vector<T, N>& x);
+constexpr vector<T, sizeof...(Indices)> swizzle(const vector<T, N>& x) noexcept;
 
 /**
  * Performs a element-wise trunc operation.
@@ -737,19 +721,6 @@ constexpr vector<T, sizeof...(Indices)> swizzle(const vector<T, N>& x);
 template <class T, std::size_t N>
 constexpr vector<T, N> trunc(const vector<T, N>& x);
 
-/**
- * Types casts each vector element and returns a vector of the casted type.
- *
- * @tparam U Target vector element type.
- * @tparam T Source vector element type.
- * @tparam N Number of dimensions.
- * @param x Vector to type cast.
- *
- * @return Type-casted vector.
- */
-template <class U, class T, std::size_t N>
-constexpr vector<U, N> type_cast(const vector<T, N>& x);
-
 /// @private
 template <class T, std::size_t N, std::size_t... I>
 constexpr inline vector<T, N> abs(const vector<T, N>& x, std::index_sequence<I...>)
@@ -765,63 +736,56 @@ constexpr inline vector<T, N> abs(const vector<T, N>& x)
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<T, N> add(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+constexpr inline vector<T, N> add(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>) noexcept
 {
 	return {(x[I] + y[I])...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<T, N> add(const vector<T, N>& x, const vector<T, N>& y)
+constexpr inline vector<T, N> add(const vector<T, N>& x, const vector<T, N>& y) noexcept
 {
 	return add(x, y, std::make_index_sequence<N>{});
 }
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<T, N> add(const vector<T, N>& x, T y, std::index_sequence<I...>)
+constexpr inline vector<T, N> add(const vector<T, N>& x, T y, std::index_sequence<I...>) noexcept
 {
 	return {(x[I] + y)...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<T, N> add(const vector<T, N>& x, T y)
+constexpr inline vector<T, N> add(const vector<T, N>& x, T y) noexcept
 {
 	return add(x, y, std::make_index_sequence<N>{});
 }
 
 /// @private
 template <std::size_t N, std::size_t... I>
-constexpr inline bool all(const vector<bool, N>& x, std::index_sequence<I...>)
+constexpr inline bool all(const vector<bool, N>& x, std::index_sequence<I...>) noexcept
 {
 	return (x[I] && ...);
 }
 
 template <std::size_t N>
-constexpr inline bool all(const vector<bool, N>& x)
+constexpr inline bool all(const vector<bool, N>& x) noexcept
 {
 	return all(x, std::make_index_sequence<N>{});
 }
 
 /// @private
 template <std::size_t N, std::size_t... I>
-constexpr inline bool any(const vector<bool, N>& x, std::index_sequence<I...>)
+constexpr inline bool any(const vector<bool, N>& x, std::index_sequence<I...>) noexcept
 {
 	return (x[I] || ...);
 }
 
 template <std::size_t N>
-constexpr inline bool any(const vector<bool, N>& x)
+constexpr inline bool any(const vector<bool, N>& x) noexcept
 {
 	return any(x, std::make_index_sequence<N>{});
 }
 
-template <std::size_t N, typename T>
-constexpr inline vector<T, N>& as_vector(T& data)
-{
-	static_assert(std::is_pod<vector<T, N>>::value);
-	return reinterpret_cast<vector<T, N>&>(data);
-}
-
 /// @private
 template <class T, std::size_t N, std::size_t... I>
 constexpr inline vector<T, N> ceil(const vector<T, N>& x, std::index_sequence<I...>)
@@ -844,22 +808,22 @@ constexpr inline vector<T, N> clamp(const vector<T, N>& x, const vector<T, N>& m
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<T, N> clamp(const vector<T, N>& x, const vector<T, N>& min_val, const vector<T, N>& max_val)
+constexpr inline vector<T, N> clamp(const vector<T, N>& x, const vector<T, N>& min, const vector<T, N>& max)
 {
-	return clamp(x, min_val, max_val, std::make_index_sequence<N>{});
+	return clamp(x, min, max, std::make_index_sequence<N>{});
 }
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<T, N> clamp(const vector<T, N>& x, T min_val, T max_val, std::index_sequence<I...>)
+constexpr inline vector<T, N> clamp(const vector<T, N>& x, T min, T max, std::index_sequence<I...>)
 {
-	return {std::min<T>(max_val, std::max<T>(min_val, x[I]))...};
+	return {std::min<T>(max, std::max<T>(min, x[I]))...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<T, N> clamp(const vector<T, N>& x, T min_val, T max_val)
+constexpr inline vector<T, N> clamp(const vector<T, N>& x, T min, T max)
 {
-	return clamp(x, min_val, max_val, std::make_index_sequence<N>{});
+	return clamp(x, min, max, std::make_index_sequence<N>{});
 }
 
 template <class T, std::size_t N>
@@ -871,7 +835,7 @@ vector<T, N> clamp_length(const vector<T, N>& x, T max_length)
 }
 
 template <class T>
-constexpr inline vector<T, 3> cross(const vector<T, 3>& x, const vector<T, 3>& y)
+constexpr inline vector<T, 3> cross(const vector<T, 3>& x, const vector<T, 3>& y) noexcept
 {
 	return
 	{
@@ -889,60 +853,72 @@ inline T distance(const vector<T, N>& p0, const vector<T, N>& p1)
 }
 
 template <class T, std::size_t N>
-constexpr inline T distance_squared(const vector<T, N>& p0, const vector<T, N>& p1)
+constexpr inline T distance_squared(const vector<T, N>& p0, const vector<T, N>& p1) noexcept
 {
-	static_assert(std::is_floating_point<T>::value);
 	return length_squared(sub(p0, p1));
 }
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<T, N> div(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+constexpr inline vector<T, N> div(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>) noexcept
 {
 	return {(x[I] / y[I])...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<T, N> div(const vector<T, N>& x, const vector<T, N>& y)
+constexpr inline vector<T, N> div(const vector<T, N>& x, const vector<T, N>& y) noexcept
 {
 	return div(x, y, std::make_index_sequence<N>{});
 }
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<T, N> div(const vector<T, N>& x, T y, std::index_sequence<I...>)
+constexpr inline vector<T, N> div(const vector<T, N>& x, T y, std::index_sequence<I...>) noexcept
 {
 	return {(x[I] / y)...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<T, N> div(const vector<T, N>& x, T y)
+constexpr inline vector<T, N> div(const vector<T, N>& x, T y) noexcept
+{
+	return div(x, y, std::make_index_sequence<N>{});
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+constexpr inline vector<T, N> div(T x, const vector<T, N>& y, std::index_sequence<I...>) noexcept
+{
+	return {(x / y[I])...};
+}
+
+template <class T, std::size_t N>
+constexpr inline vector<T, N> div(T x, const vector<T, N>& y) noexcept
 {
 	return div(x, y, std::make_index_sequence<N>{});
 }
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline T dot(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+constexpr inline T dot(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>) noexcept
 {
 	return ((x[I] * y[I]) + ...);
 }
 
 template <class T, std::size_t N>
-constexpr inline T dot(const vector<T, N>& x, const vector<T, N>& y)
+constexpr inline T dot(const vector<T, N>& x, const vector<T, N>& y) noexcept
 {
 	return dot(x, y, std::make_index_sequence<N>{});
 }
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<bool, N> equal(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+constexpr inline vector<bool, N> equal(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>) noexcept
 {
 	return {(x[I] == y[I])...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<bool, N> equal(const vector<T, N>& x, const vector<T, N>& y)
+constexpr inline vector<bool, N> equal(const vector<T, N>& x, const vector<T, N>& y) noexcept
 {
 	return equal(x, y, std::make_index_sequence<N>{});
 }
@@ -1005,26 +981,26 @@ constexpr inline vector<T, N> fract(const vector<T, N>& x)
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<bool, N> greater_than(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+constexpr inline vector<bool, N> greater_than(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>) noexcept
 {
 	return {(x[I] > y[I])...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<bool, N> greater_than(const vector<T, N>& x, const vector<T, N>& y)
+constexpr inline vector<bool, N> greater_than(const vector<T, N>& x, const vector<T, N>& y) noexcept
 {
 	return greater_than(x, y, std::make_index_sequence<N>{});
 }
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<bool, N> greater_than_equal(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+constexpr inline vector<bool, N> greater_than_equal(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>) noexcept
 {
 	return {(x[I] >= y[I])...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<bool, N> greater_than_equal(const vector<T, N>& x, const vector<T, N>& y)
+constexpr inline vector<bool, N> greater_than_equal(const vector<T, N>& x, const vector<T, N>& y) noexcept
 {
 	return greater_than_equal(x, y, std::make_index_sequence<N>{});
 }
@@ -1037,33 +1013,33 @@ inline T length(const vector<T, N>& x)
 }
 
 template <class T, std::size_t N>
-constexpr inline T length_squared(const vector<T, N>& x)
+constexpr inline T length_squared(const vector<T, N>& x) noexcept
 {
 	return dot(x, x);
 }
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<bool, N> less_than(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+constexpr inline vector<bool, N> less_than(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>) noexcept
 {
 	return {(x[I] < y[I])...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<bool, N> less_than(const vector<T, N>& x, const vector<T, N>& y)
+constexpr inline vector<bool, N> less_than(const vector<T, N>& x, const vector<T, N>& y) noexcept
 {
 	return less_than(x, y, std::make_index_sequence<N>{});
 }
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<bool, N> less_than_equal(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+constexpr inline vector<bool, N> less_than_equal(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>) noexcept
 {
 	return {(x[I] <= y[I])...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<bool, N> less_than_equal(const vector<T, N>& x, const vector<T, N>& y)
+constexpr inline vector<bool, N> less_than_equal(const vector<T, N>& x, const vector<T, N>& y) noexcept
 {
 	return less_than_equal(x, y, std::make_index_sequence<N>{});
 }
@@ -1136,39 +1112,39 @@ constexpr inline vector<T, N> mod(const vector<T, N>& x, T y)
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<T, N> mul(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+constexpr inline vector<T, N> mul(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>) noexcept
 {
 	return {(x[I] * y[I])...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<T, N> mul(const vector<T, N>& x, const vector<T, N>& y)
+constexpr inline vector<T, N> mul(const vector<T, N>& x, const vector<T, N>& y) noexcept
 {
 	return mul(x, y, std::make_index_sequence<N>{});
 }
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<T, N> mul(const vector<T, N>& x, T y, std::index_sequence<I...>)
+constexpr inline vector<T, N> mul(const vector<T, N>& x, T y, std::index_sequence<I...>) noexcept
 {
 	return {(x[I] * y)...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<T, N> mul(const vector<T, N>& x, T y)
+constexpr inline vector<T, N> mul(const vector<T, N>& x, T y) noexcept
 {
 	return mul(x, y, std::make_index_sequence<N>{});
 }
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<T, N> negate(const vector<T, N>& x, std::index_sequence<I...>)
+constexpr inline vector<T, N> negate(const vector<T, N>& x, std::index_sequence<I...>) noexcept
 {
 	return {(-x[I])...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<T, N> negate(const vector<T, N>& x)
+constexpr inline vector<T, N> negate(const vector<T, N>& x) noexcept
 {
 	return negate(x, std::make_index_sequence<N>{});
 }
@@ -1177,31 +1153,31 @@ template <class T, std::size_t N>
 inline vector<T, N> normalize(const vector<T, N>& x)
 {
 	static_assert(std::is_floating_point<T>::value);
-	return mul(x, T(1) / length(x));
+	return mul(x, T{1} / length(x));
 }
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<bool, N> not(const vector<T, N>& x, std::index_sequence<I...>)
+constexpr inline vector<bool, N> not(const vector<T, N>& x, std::index_sequence<I...>) noexcept
 {
 	return {!x[I]...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<bool, N> not(const vector<T, N>& x)
+constexpr inline vector<bool, N> not(const vector<T, N>& x) noexcept
 {
 	return not(x, std::make_index_sequence<N>{});
 }
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<bool, N> not_equal(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+constexpr inline vector<bool, N> not_equal(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>) noexcept
 {
 	return {(x[I] != y[I])...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<bool, N> not_equal(const vector<T, N>& x, const vector<T, N>& y)
+constexpr inline vector<bool, N> not_equal(const vector<T, N>& x, const vector<T, N>& y) noexcept
 {
 	return not_equal(x, y, std::make_index_sequence<N>{});
 }
@@ -1234,19 +1210,6 @@ inline vector<T, N> pow(const vector<T, N>& x, T y)
 	return pow(x, y, std::make_index_sequence<N>{});
 }
 
-/// @private
-template <std::size_t M, class T, std::size_t N, std::size_t... I>
-constexpr inline vector<T, M> resize(const vector<T, N>& x, std::index_sequence<I...>)
-{
-	return {(I < N) ? x[I] : T{0} ...};
-}
-
-template <std::size_t M, class T, std::size_t N>
-constexpr inline vector<T, M> resize(const vector<T, N>& x)
-{
-	return resize<M>(x, std::make_index_sequence<M>{});
-}
-
 /// @private
 template <class T, std::size_t N, std::size_t... I>
 constexpr inline vector<T, N> round(const vector<T, N>& x, std::index_sequence<I...>)
@@ -1291,45 +1254,58 @@ inline vector<T, N> sqrt(const vector<T, N>& x, const vector<T, N>& y)
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<T, N> sub(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>)
+constexpr inline vector<T, N> sub(const vector<T, N>& x, const vector<T, N>& y, std::index_sequence<I...>) noexcept
 {
 	return {(x[I] - y[I])...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<T, N> sub(const vector<T, N>& x, const vector<T, N>& y)
+constexpr inline vector<T, N> sub(const vector<T, N>& x, const vector<T, N>& y) noexcept
 {
 	return sub(x, y, std::make_index_sequence<N>{});
 }
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline vector<T, N> sub(const vector<T, N>& x, T y, std::index_sequence<I...>)
+constexpr inline vector<T, N> sub(const vector<T, N>& x, T y, std::index_sequence<I...>) noexcept
 {
 	return {(x[I] - y)...};
 }
 
 template <class T, std::size_t N>
-constexpr inline vector<T, N> sub(const vector<T, N>& x, T y)
+constexpr inline vector<T, N> sub(const vector<T, N>& x, T y) noexcept
+{
+	return sub(x, y, std::make_index_sequence<N>{});
+}
+
+/// @private
+template <class T, std::size_t N, std::size_t... I>
+constexpr inline vector<T, N> sub(T x, const vector<T, N>& y, std::index_sequence<I...>) noexcept
+{
+	return {(x - y[I])...};
+}
+
+template <class T, std::size_t N>
+constexpr inline vector<T, N> sub(T x, const vector<T, N>& y) noexcept
 {
 	return sub(x, y, std::make_index_sequence<N>{});
 }
 
 /// @private
 template <class T, std::size_t N, std::size_t... I>
-constexpr inline T sum(const vector<T, N>& x, std::index_sequence<I...>)
+constexpr inline T sum(const vector<T, N>& x, std::index_sequence<I...>) noexcept
 {
 	return (x[I] + ...);
 }
 
 template <class T, std::size_t N>
-constexpr inline T sum(const vector<T, N>& x)
+constexpr inline T sum(const vector<T, N>& x) noexcept
 {
 	return sum(x, std::make_index_sequence<N>{});
 }
 
 template <std::size_t... Indices, class T, std::size_t N>
-constexpr inline vector<T, sizeof...(Indices)> swizzle(const vector<T, N>& x)
+constexpr inline vector<T, sizeof...(Indices)> swizzle(const vector<T, N>& x) noexcept
 {
 	return {x[Indices]...};
 }
@@ -1348,64 +1324,98 @@ constexpr inline vector<T, N> trunc(const vector<T, N>& x)
 	return trunc(x, std::make_index_sequence<N>{});
 }
 
-/// @private
-template <class U, class T, std::size_t N, std::size_t... I>
-constexpr inline vector<U, N> type_cast(const vector<T, N>& x, std::index_sequence<I...>)
-{
-	return {static_cast<U>(x[I])...};
-}
-
-template <class U, class T, std::size_t N>
-constexpr inline vector<U, N> type_cast(const vector<T, N>& x)
-{
-	return type_cast<U>(x, std::make_index_sequence<N>{}); 
-}
-
 } // namespace math
 
 /// @copydoc math::add(const math::vector<T, N>&, const math::vector<T, N>&)
 template <class T, std::size_t N>
-constexpr math::vector<T, N> operator+(const math::vector<T, N>& x, const math::vector<T, N>& y);
+constexpr inline math::vector<T, N> operator+(const math::vector<T, N>& x, const math::vector<T, N>& y) noexcept
+{
+	return math::add(x, y);
+}
 
 /// @copydoc math::add(const math::vector<T, N>&, T)
 /// @{
 template <class T, std::size_t N>
-constexpr math::vector<T, N> operator+(const math::vector<T, N>& x, T y);
+constexpr inline math::vector<T, N> operator+(const math::vector<T, N>& x, T y) noexcept
+{
+	return math::add(x, y);
+}
 template <class T, std::size_t N>
-constexpr math::vector<T, N> operator+(T x, const math::vector<T, N>& y);
+constexpr inline math::vector<T, N> operator+(T x, const math::vector<T, N>& y) noexcept
+{
+	return math::add(y, x);
+}
 /// @}
 
 /// @copydoc math::div(const math::vector<T, N>&, const math::vector<T, N>&)
 template <class T, std::size_t N>
-constexpr math::vector<T, N> operator/(const math::vector<T, N>& x, const math::vector<T, N>& y);
+constexpr inline math::vector<T, N> operator/(const math::vector<T, N>& x, const math::vector<T, N>& y) noexcept
+{
+	return math::div(x, y);
+}
 
-/// @copydoc math::div(const math::vector<T, N>&, T y)
+/// @copydoc math::div(const math::vector<T, N>&, T)
 template <class T, std::size_t N>
-constexpr math::vector<T, N> operator/(const math::vector<T, N>& x, T y);
+constexpr inline math::vector<T, N> operator/(const math::vector<T, N>& x, T y) noexcept
+{
+	return math::div(x, y);
+}
+
+/// @copydoc math::div(T, const math::vector<T, N>&)
+template <class T, std::size_t N>
+constexpr inline math::vector<T, N> operator/(T x, const math::vector<T, N>& y) noexcept
+{
+	return math::div(x, y);
+}
 
 /// @copydoc math::mul(const math::vector<T, N>&, const math::vector<T, N>&)
 template <class T, std::size_t N>
-constexpr math::vector<T, N> operator*(const math::vector<T, N>& x, const math::vector<T, N>& y);
+constexpr inline math::vector<T, N> operator*(const math::vector<T, N>& x, const math::vector<T, N>& y) noexcept
+{
+	return math::mul(x, y);
+}
 
 /// @copydoc math::mul(const math::vector<T, N>&, T)
 /// @{
 template <class T, std::size_t N>
-constexpr math::vector<T, N> operator*(const math::vector<T, N>& x, T y);
+constexpr inline math::vector<T, N> operator*(const math::vector<T, N>& x, T y) noexcept
+{
+	return math::mul(x, y);
+}
 template <class T, std::size_t N>
-constexpr math::vector<T, N> operator*(T x, const math::vector<T, N>& y);
+constexpr inline math::vector<T, N> operator*(T x, const math::vector<T, N>& y) noexcept
+{
+	return math::mul(y, x);
+}
 /// @}
 
 /// @copydoc math::negate(const math::vector<T, N>&)
 template <class T, std::size_t N>
-constexpr math::vector<T, N> operator-(const math::vector<T, N>& x);
+constexpr inline math::vector<T, N> operator-(const math::vector<T, N>& x) noexcept
+{
+	return math::negate(x);
+}
 
 /// @copydoc math::sub(const math::vector<T, N>&, const math::vector<T, N>&)
 template <class T, std::size_t N>
-constexpr math::vector<T, N> operator-(const math::vector<T, N>& x, const math::vector<T, N>& y);
+constexpr inline math::vector<T, N> operator-(const math::vector<T, N>& x, const math::vector<T, N>& y) noexcept
+{
+	return math::sub(x, y);
+}
 
 /// @copydoc math::sub(const math::vector<T, N>&, T)
 template <class T, std::size_t N>
-constexpr math::vector<T, N> operator-(const math::vector<T, N>& x, T y);
+constexpr inline math::vector<T, N> operator-(const math::vector<T, N>& x, T y) noexcept
+{
+	return math::sub(x, y);
+}
+
+/// @copydoc math::sub(T, const math::vector<T, N>&)
+template <class T, std::size_t N>
+constexpr inline math::vector<T, N> operator-(T x, const math::vector<T, N>& y) noexcept
+{
+	return math::sub(x, y);
+}
 
 /**
  * Adds two values and stores the result in the first value.
@@ -1416,9 +1426,15 @@ constexpr math::vector<T, N> operator-(const math::vector<T, N>& x, T y);
  */
 /// @{
 template <class T, std::size_t N>
-constexpr math::vector<T, N>& operator+=(math::vector<T, N>& x, const math::vector<T, N>& y);
+constexpr inline math::vector<T, N>& operator+=(math::vector<T, N>& x, const math::vector<T, N>& y) noexcept
+{
+	return (x = x + y);
+}
 template <class T, std::size_t N>
-constexpr math::vector<T, N>& operator+=(math::vector<T, N>& x, T y);
+constexpr inline math::vector<T, N>& operator+=(math::vector<T, N>& x, T y) noexcept
+{
+	return (x = x + y);
+}
 /// @}
 
 /**
@@ -1430,9 +1446,15 @@ constexpr math::vector<T, N>& operator+=(math::vector<T, N>& x, T y);
  */
 /// @{
 template <class T, std::size_t N>
-constexpr math::vector<T, N>& operator-=(math::vector<T, N>& x, const math::vector<T, N>& y);
+constexpr inline math::vector<T, N>& operator-=(math::vector<T, N>& x, const math::vector<T, N>& y) noexcept
+{
+	return (x = x - y);
+}
 template <class T, std::size_t N>
-constexpr math::vector<T, N>& operator-=(math::vector<T, N>& x, T y);
+constexpr inline math::vector<T, N>& operator-=(math::vector<T, N>& x, T y) noexcept
+{
+	return (x = x - y);
+}
 /// @}
 
 /**
@@ -1444,9 +1466,15 @@ constexpr math::vector<T, N>& operator-=(math::vector<T, N>& x, T y);
  */
 /// @{
 template <class T, std::size_t N>
-constexpr math::vector<T, N>& operator*=(math::vector<T, N>& x, const math::vector<T, N>& y);
+constexpr inline math::vector<T, N>& operator*=(math::vector<T, N>& x, const math::vector<T, N>& y) noexcept
+{
+	return (x = x * y);
+}
 template <class T, std::size_t N>
-constexpr math::vector<T, N>& operator*=(math::vector<T, N>& x, T y);
+constexpr inline math::vector<T, N>& operator*=(math::vector<T, N>& x, T y) noexcept
+{
+	return (x = x * y);
+}
 /// @}
 
 /**
@@ -1458,9 +1486,15 @@ constexpr math::vector<T, N>& operator*=(math::vector<T, N>& x, T y);
  */
 /// @{
 template <class T, std::size_t N>
-constexpr math::vector<T, N>& operator/=(math::vector<T, N>& x, const math::vector<T, N>& y);
+constexpr inline math::vector<T, N>& operator/=(math::vector<T, N>& x, const math::vector<T, N>& y) noexcept
+{
+	return (x = x / y);
+}
 template <class T, std::size_t N>
-constexpr math::vector<T, N>& operator/=(math::vector<T, N>& x, T y);
+constexpr inline math::vector<T, N>& operator/=(math::vector<T, N>& x, T y) noexcept
+{
+	return (x = x / y);
+}
 /// @}
 
 /**
@@ -1471,134 +1505,6 @@ constexpr math::vector<T, N>& operator/=(math::vector<T, N>& x, T y);
  *
  * @return Output stream.
  */
-template <class T, std::size_t N>
-std::ostream& operator<<(std::ostream& os, const math::vector<T, N>& x);
-
-/**
- * Reads the elements of a vector from an input stream, with each element delimeted by a space.
- *
- * @param is Input stream.
- * @param x Vector.
- *
- * @return Input stream.
- */
-template <class T, std::size_t N>
-std::istream& operator>>(std::istream& is, math::vector<T, N>& v);
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N> operator+(const math::vector<T, N>& x, const math::vector<T, N>& y)
-{
-	return math::add(x, y);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N> operator+(const math::vector<T, N>& x, T y)
-{
-	return math::add(x, y);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N> operator+(T x, const math::vector<T, N>& y)
-{
-	return math::add(y, x);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N> operator-(const math::vector<T, N>& x)
-{
-	return math::negate(x);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N> operator-(const math::vector<T, N>& x, const math::vector<T, N>& y)
-{
-	return math::sub(x, y);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N> operator-(const math::vector<T, N>& x, T y)
-{
-	return math::sub(x, y);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N> operator*(const math::vector<T, N>& x, const math::vector<T, N>& y)
-{
-	return math::mul(x, y);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N> operator*(const math::vector<T, N>& x, T y)
-{
-	return math::mul(x, y);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N> operator*(T x, const math::vector<T, N>& y)
-{
-	return math::mul(y, x);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N> operator/(const math::vector<T, N>& x, const math::vector<T, N>& y)
-{
-	return math::div(x, y);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N> operator/(const math::vector<T, N>& x, T y)
-{
-	return math::div(x, y);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N>& operator+=(math::vector<T, N>& x, const math::vector<T, N>& y)
-{
-	return (x = x + y);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N>& operator+=(math::vector<T, N>& x, T y)
-{
-	return (x = x + y);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N>& operator-=(math::vector<T, N>& x, const math::vector<T, N>& y)
-{
-	return (x = x - y);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N>& operator-=(math::vector<T, N>& x, T y)
-{
-	return (x = x - y);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N>& operator*=(math::vector<T, N>& x, const math::vector<T, N>& y)
-{
-	return (x = x * y);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N>& operator*=(math::vector<T, N>& x, T y)
-{
-	return (x = x * y);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N>& operator/=(math::vector<T, N>& x, const math::vector<T, N>& y)
-{
-	return (x = x / y);
-}
-
-template <class T, std::size_t N>
-constexpr inline math::vector<T, N>& operator/=(math::vector<T, N>& x, T y)
-{
-	return (x = x / y);
-}
-
 template <class T, std::size_t N>
 std::ostream& operator<<(std::ostream& os, const math::vector<T, N>& x)
 {
@@ -1612,6 +1518,14 @@ std::ostream& operator<<(std::ostream& os, const math::vector<T, N>& x)
 	return os;
 }
 
+/**
+ * Reads the elements of a vector from an input stream, with each element delimeted by a space.
+ *
+ * @param is Input stream.
+ * @param x Vector.
+ *
+ * @return Input stream.
+ */
 template <class T, std::size_t N>
 std::istream& operator>>(std::istream& is, math::vector<T, N>& x)
 {
diff --git a/src/render/passes/ground-pass.cpp b/src/render/passes/ground-pass.cpp
index 1feeece..f74e996 100644
--- a/src/render/passes/ground-pass.cpp
+++ b/src/render/passes/ground-pass.cpp
@@ -86,7 +86,7 @@ void ground_pass::render(const render::context& ctx, render::queue& queue) const
 	float clip_far = camera.get_clip_far_tween().interpolate(ctx.alpha);
 	float3 model_scale = float3{1.0f, 1.0f, 1.0f} * (clip_near + clip_far) * 0.5f;
 	float4x4 model = math::scale(math::matrix4<float>::identity(), model_scale);
-	float4x4 view = math::resize<4, 4>(math::resize<3, 3>(ctx.view));
+	float4x4 view = float4x4(float3x3(ctx.view));
 	float4x4 model_view = view * model;
 	const float4x4& projection = ctx.projection;
 	const float4x4& view_projection = ctx.view_projection;
diff --git a/src/render/passes/material-pass.cpp b/src/render/passes/material-pass.cpp
index d8092bc..452b1c3 100644
--- a/src/render/passes/material-pass.cpp
+++ b/src/render/passes/material-pass.cpp
@@ -43,6 +43,7 @@
 #include "scene/spot-light.hpp"
 #include "config.hpp"
 #include "math/quaternion-operators.hpp"
+#include "math/projection.hpp"
 #include <cmath>
 #include <glad/glad.h>
 
@@ -504,8 +505,8 @@ void material_pass::render(const render::context& ctx, render::queue& queue) con
 		model = operation.transform;
 		model_view_projection = view_projection * model;
 		model_view = view * model;
-		normal_model = math::transpose(math::inverse(math::resize<3, 3>(model)));
-		normal_model_view = math::transpose(math::inverse(math::resize<3, 3>(model_view)));
+		normal_model = math::transpose(math::inverse(math::matrix<float, 3, 3>(model)));
+		normal_model_view = math::transpose(math::inverse(math::matrix<float, 3, 3>(model_view)));
 		
 		// Skinning palette
 		if (operation.bone_count && parameters->skinning_palette)
diff --git a/src/render/passes/shadow-map-pass.cpp b/src/render/passes/shadow-map-pass.cpp
index f844167..01f93fb 100644
--- a/src/render/passes/shadow-map-pass.cpp
+++ b/src/render/passes/shadow-map-pass.cpp
@@ -35,6 +35,7 @@
 #include "math/vector.hpp"
 #include "math/matrix.hpp"
 #include "math/quaternion-operators.hpp"
+#include "math/projection.hpp"
 #include <cmath>
 #include <glad/glad.h>
 
diff --git a/src/render/passes/sky-pass.cpp b/src/render/passes/sky-pass.cpp
index e675210..8cdd733 100644
--- a/src/render/passes/sky-pass.cpp
+++ b/src/render/passes/sky-pass.cpp
@@ -194,7 +194,7 @@ void sky_pass::render(const render::context& ctx, render::queue& queue) const
 	float clip_far = camera.get_clip_far_tween().interpolate(ctx.alpha);
 	float3 model_scale = float3{1.0f, 1.0f, 1.0f} * (clip_near + clip_far) * 0.5f;
 	float4x4 model = math::scale(math::matrix4<float>::identity(), model_scale);
-	float4x4 view = math::resize<4, 4>(math::resize<3, 3>(ctx.view));
+	float4x4 view = float4x4(float3x3(ctx.view));
 	float4x4 model_view = view * model;
 	const float4x4& projection = ctx.projection;
 	float4x4 view_projection = projection * view;
@@ -300,7 +300,7 @@ void sky_pass::render(const render::context& ctx, render::queue& queue) const
 	{
 		float star_distance = (clip_near + clip_far) * 0.5f;
 		
-		model = math::resize<4, 4>(math::matrix_cast<float>(icrf_to_eus.r));
+		model = float4x4(math::matrix_cast<float>(icrf_to_eus.r));
 		model = math::scale(model, {star_distance, star_distance, star_distance});
 		
 		model_view = view * model;
@@ -332,7 +332,7 @@ void sky_pass::render(const render::context& ctx, render::queue& queue) const
 		moon_transform.scale = {moon_radius, moon_radius, moon_radius};
 		
 		model = math::matrix_cast(moon_transform);		
-		float3x3 normal_model = math::transpose(math::inverse(math::resize<3, 3>(model)));
+		float3x3 normal_model = math::transpose(math::inverse(float3x3(model)));
 		
 		rasterizer->use_program(*moon_shader_program);
 		if (moon_model_input)
diff --git a/src/render/renderer.cpp b/src/render/renderer.cpp
index ff54494..cf0be48 100644
--- a/src/render/renderer.cpp
+++ b/src/render/renderer.cpp
@@ -176,7 +176,7 @@ void renderer::process_model_instance(const render::context& ctx, render::queue&
 
 	render::operation operation;
 	operation.transform = math::matrix_cast(model_instance->get_transform_tween().interpolate(ctx.alpha));
-	operation.depth = ctx.clip_near.signed_distance(math::resize<3>(operation.transform[3]));
+	operation.depth = ctx.clip_near.signed_distance(float3(operation.transform[3]));
 	operation.vertex_array = model->get_vertex_array();
 	operation.instance_count = model_instance->get_instance_count();
 	
@@ -223,7 +223,7 @@ void renderer::process_billboard(const render::context& ctx, render::queue& queu
 	
 	math::transform<float> billboard_transform = billboard->get_transform_tween().interpolate(ctx.alpha);
 	billboard_op.material = billboard->get_material();
-	billboard_op.depth = ctx.clip_near.signed_distance(math::resize<3>(billboard_transform.translation));
+	billboard_op.depth = ctx.clip_near.signed_distance(float3(billboard_transform.translation));
 	
 	// Align billboard
 	if (billboard->get_billboard_type() == scene::billboard_type::spherical)
diff --git a/src/render/renderer.hpp b/src/render/renderer.hpp
index ddf17e1..79cf1db 100644
--- a/src/render/renderer.hpp
+++ b/src/render/renderer.hpp
@@ -52,7 +52,7 @@ public:
 	 * @param t Current time, in seconds.
 	 * @param dt Timestep, in seconds.
 	 * @param alpha Subframe interpolation factor.
-	 * @parma collection Collection of scene objects to render.
+	 * @param collection Collection of scene objects to render.
 	 */
 	void render(float t, float dt, float alpha, const scene::collection& collection) const;
 	
diff --git a/src/scene/camera.cpp b/src/scene/camera.cpp
index a6c7327..4c82fe1 100644
--- a/src/scene/camera.cpp
+++ b/src/scene/camera.cpp
@@ -22,6 +22,7 @@
 #include "math/constants.hpp"
 #include "math/interpolation.hpp"
 #include "math/quaternion-operators.hpp"
+#include "math/projection.hpp"
 
 namespace scene {
 
@@ -101,7 +102,7 @@ float3 camera::project(const float3& object, const float4& viewport) const
 	result[0] = result[0] * viewport[2] + viewport[0];
 	result[1] = result[1] * viewport[3] + viewport[1];
 	
-	return math::resize<3>(result);
+	return math::vector<float, 3>(result);
 }
 
 float3 camera::unproject(const float3& window, const float4& viewport) const
@@ -116,7 +117,7 @@ float3 camera::unproject(const float3& window, const float4& viewport) const
 	
 	result = math::inverse(view_projection[1]) * result;
 
-	return math::resize<3>(result) * (1.0f / result[3]);
+	return math::vector<float, 3>(result) * (1.0f / result[3]);
 }
 
 void camera::set_perspective(float fov, float aspect_ratio, float clip_near, float clip_far)
diff --git a/src/scene/text.cpp b/src/scene/text.cpp
index 2f5b8ee..92839df 100644
--- a/src/scene/text.cpp
+++ b/src/scene/text.cpp
@@ -104,7 +104,7 @@ void text::render(const render::context& ctx, render::queue& queue) const
 		return;
 	
 	render_op.transform = math::matrix_cast(get_transform_tween().interpolate(ctx.alpha));
-	render_op.depth = ctx.clip_near.signed_distance(math::resize<3>(render_op.transform[3]));
+	render_op.depth = ctx.clip_near.signed_distance(math::vector<float, 3>(render_op.transform[3]));
 	queue.push_back(render_op);
 }