diff --git a/CMakeLists.txt b/CMakeLists.txt
index 20085212..88d77454 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -435,6 +435,9 @@ add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/librapid/vendor/xsimd")
 target_compile_definitions(fmt PUBLIC FMT_HEADER_ONLY)
 target_link_libraries(${module_name} PUBLIC fmt xsimd)
 
+# Disable "parameter passing for argument of type ... changed to match ..." for xsimd
+target_compile_options(xsimd INTERFACE $<$<OR:$<CXX_COMPILER_ID:GNU>,$<CXX_COMPILER_ID:Clang>>:-Wno-psabi>)
+
 if (${LIBRAPID_USE_MULTIPREC})
     # Load MPIR
     find_package(MPIR QUIET)
diff --git a/librapid/include/librapid/array/arrayContainer.hpp b/librapid/include/librapid/array/arrayContainer.hpp
index 4768ce6f..5ee80f21 100644
--- a/librapid/include/librapid/array/arrayContainer.hpp
+++ b/librapid/include/librapid/array/arrayContainer.hpp
@@ -56,7 +56,13 @@ namespace librapid {
 			static constexpr bool supportsArithmetic   = TypeInfo<Scalar>::supportsArithmetic;
 			static constexpr bool supportsLogical	   = TypeInfo<Scalar>::supportsLogical;
 			static constexpr bool supportsBinary	   = TypeInfo<Scalar>::supportsBinary;
-			static constexpr bool allowVectorisation   = TypeInfo<Scalar>::packetWidth > 1;
+			static constexpr bool allowVectorisation   = []() {
+				  if constexpr (typetraits::HasAllowVectorisation<TypeInfo<Scalar>>::value) {
+					  return TypeInfo<Scalar>::allowVectorisation;
+				  } else {
+					  return TypeInfo<Scalar>::packetWidth > 1;
+				  }
+			}();
 
 #if defined(LIBRAPID_HAS_CUDA)
 			static constexpr cudaDataType_t CudaType = TypeInfo<Scalar>::CudaType;
@@ -67,15 +73,6 @@ namespace librapid {
 			static constexpr int64_t canMemcpy = false;
 		};
 
-		/// Evaluates as true if the input type is an ArrayContainer instance
-		/// \tparam T Input type
-		template<typename T>
-		struct IsArrayContainer : std::false_type {};
-
-		template<typename ShapeType, typename StorageScalar>
-		struct IsArrayContainer<array::ArrayContainer<ShapeType, StorageScalar>> : std::true_type {
-		};
-
 		LIBRAPID_DEFINE_AS_TYPE(typename StorageScalar,
 								array::ArrayContainer<Shape COMMA StorageScalar>);
 
@@ -686,7 +683,7 @@ namespace librapid {
 			  sizeof...(Indices),
 			  m_shape.ndim());
 
-			int dim = 0;
+			int dim		  = 0;
 			int64_t index = 0;
 			for (int64_t i : {indices...}) {
 				LIBRAPID_ASSERT(
diff --git a/librapid/include/librapid/array/assignOps.hpp b/librapid/include/librapid/array/assignOps.hpp
index 3ba5d15b..ba95ae6c 100644
--- a/librapid/include/librapid/array/assignOps.hpp
+++ b/librapid/include/librapid/array/assignOps.hpp
@@ -24,11 +24,17 @@ namespace librapid {
 			using Function = detail::Function<descriptor::Trivial, Functor_, Args...>;
 			using Scalar =
 			  typename array::ArrayContainer<ShapeType_, Storage<StorageScalar>>::Scalar;
-			constexpr int64_t packetWidth = typetraits::TypeInfo<Scalar>::packetWidth;
 			constexpr bool allowVectorisation =
 			  typetraits::TypeInfo<
 				detail::Function<descriptor::Trivial, Functor_, Args...>>::allowVectorisation &&
 			  Function::argsAreSameType;
+			constexpr int64_t packetWidth = []() {
+				if constexpr (allowVectorisation) {
+					return typetraits::TypeInfo<Scalar>::packetWidth;
+				} else {
+					return 1;
+				}
+			}();
 
 			const int64_t size		 = function.shape().size();
 			const int64_t vectorSize = size - (size % packetWidth);
@@ -74,13 +80,21 @@ namespace librapid {
 			using Scalar =
 			  typename array::ArrayContainer<ShapeType_,
 											 FixedStorage<StorageScalar, StorageSize...>>::Scalar;
-			constexpr int64_t packetWidth = typetraits::TypeInfo<Scalar>::packetWidth;
-			constexpr int64_t elements	  = ::librapid::product<StorageSize...>();
-			constexpr int64_t vectorSize  = elements - (elements % packetWidth);
+
 			constexpr bool allowVectorisation =
 			  typetraits::TypeInfo<
 				detail::Function<descriptor::Trivial, Functor_, Args...>>::allowVectorisation &&
 			  Function::argsAreSameType;
+			constexpr int64_t packetWidth = []() {
+				if constexpr (allowVectorisation) {
+					return typetraits::TypeInfo<Scalar>::packetWidth;
+				} else {
+					return 1;
+				}
+			}();
+
+			constexpr int64_t elements	 = ::librapid::product<StorageSize...>();
+			constexpr int64_t vectorSize = elements - (elements % packetWidth);
 
 			// Ensure the function can actually be assigned to the array container
 			static_assert(
@@ -124,12 +138,18 @@ namespace librapid {
 			using Function = detail::Function<descriptor::Trivial, Functor_, Args...>;
 			using Scalar =
 			  typename array::ArrayContainer<ShapeType_, Storage<StorageScalar>>::Scalar;
-			constexpr size_t packetWidth = typetraits::TypeInfo<Scalar>::packetWidth;
 
 			constexpr bool allowVectorisation =
 			  typetraits::TypeInfo<
 				detail::Function<descriptor::Trivial, Functor_, Args...>>::allowVectorisation &&
 			  Function::argsAreSameType;
+			constexpr int64_t packetWidth = []() {
+				if constexpr (allowVectorisation) {
+					return typetraits::TypeInfo<Scalar>::packetWidth;
+				} else {
+					return 1;
+				}
+			}();
 
 			const size_t size		= function.size();
 			const size_t vectorSize = size - (size % packetWidth);
@@ -180,12 +200,18 @@ namespace librapid {
 			using Scalar =
 			  typename array::ArrayContainer<ShapeType_,
 											 FixedStorage<StorageScalar, StorageSize...>>::Scalar;
-			constexpr int64_t packetWidth = typetraits::TypeInfo<Scalar>::packetWidth;
 
 			constexpr bool allowVectorisation =
 			  typetraits::TypeInfo<
 				detail::Function<descriptor::Trivial, Functor_, Args...>>::allowVectorisation &&
 			  Function::argsAreSameType;
+			constexpr int64_t packetWidth = []() {
+				if constexpr (allowVectorisation) {
+					return typetraits::TypeInfo<Scalar>::packetWidth;
+				} else {
+					return 1;
+				}
+			}();
 
 			constexpr int64_t size		 = ::librapid::product<StorageSize...>();
 			constexpr int64_t vectorSize = size - (size % packetWidth);
diff --git a/librapid/include/librapid/array/function.hpp b/librapid/include/librapid/array/function.hpp
index dc2c61c1..8d92ab87 100644
--- a/librapid/include/librapid/array/function.hpp
+++ b/librapid/include/librapid/array/function.hpp
@@ -35,13 +35,60 @@ namespace librapid {
 			}
 		}
 
+		// Normally, we want to use the scalar type of the input. This said, there are a few edge
+		// cases where it is necessary to use the type itself.
+
+		// Default
+		template<typename T>
+		struct ScalarTypeHelper {
+			using Type = typename TypeInfo<std::decay_t<T>>::Scalar;
+		};
+
+		// Vectors
+		template<typename T, uint64_t NumDims>
+		struct ScalarTypeHelper<Vector<T, NumDims>> {
+			using Type = Vector<T, NumDims>;
+		};
+
+		// Once we have the correct scalar types, we need to check if the result is a lazy-evaluated
+		// function. If so, we need to extract the actual return type from the function.
+
+		// Default
+		template<typename T>
+		struct ReturnTypeHelper {
+			using Type = T;
+		};
+
+		// Binary vector operation
+		template<typename LHS, typename RHS, typename Op>
+		struct ReturnTypeHelper<vectorDetail::BinaryVecOp<LHS, RHS, Op>> {
+			using IntermediateType = vectorDetail::BinaryVecOp<LHS, RHS, Op>;
+			using Type			   = decltype(std::declval<IntermediateType>().eval());
+		};
+
+		// Unary vector operation
+		template<typename Val, typename Op>
+		struct ReturnTypeHelper<vectorDetail::UnaryVecOp<Val, Op>> {
+			using IntermediateType = vectorDetail::UnaryVecOp<Val, Op>;
+			using Type			   = decltype(std::declval<IntermediateType>().eval());
+		};
+
 		template<typename desc, typename Functor_, typename... Args>
 		struct TypeInfo<::librapid::detail::Function<desc, Functor_, Args...>> {
 			static constexpr detail::LibRapidType type = detail::LibRapidType::ArrayFunction;
-			using Scalar							   = decltype(std::declval<Functor_>()(
-				std::declval<typename TypeInfo<std::decay_t<Args>>::Scalar>()...));
-			using Packet							   = typename TypeInfo<Scalar>::Packet;
-			using Backend							   = decltype(commonBackend<Args...>());
+			// using Scalar	= decltype(std::declval<Functor_>()(
+			// 		std::declval<typename TypeInfo<std::decay_t<Args>>::Scalar>()...));
+
+			// using Scalar = decltype(std::declval<Functor_>()(
+			//   std::declval<typename ScalarTypeHelper<Args>::Type>()...));
+
+			using TempScalar = decltype(std::declval<Functor_>()(
+			  std::declval<typename ScalarTypeHelper<Args>::Type>()...));
+
+			using Scalar = typename ReturnTypeHelper<TempScalar>::Type;
+
+			using Packet  = typename TypeInfo<Scalar>::Packet;
+			using Backend = decltype(commonBackend<Args...>());
 			using ShapeType =
 			  typename detail::ShapeTypeHelper<typename TypeInfo<Args>::ShapeType...>::Type;
 
@@ -81,16 +128,12 @@ namespace librapid {
 			return Packet(obj);
 		}
 
-		template<typename T, typename std::enable_if_t<typetraits::TypeInfo<T>::type !=
-														 ::librapid::detail::LibRapidType::Scalar,
-													   int> = 0>
+		template<typename T, typename std::enable_if_t<detail::IsArrayType<T>::val, int> = 0>
 		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto scalarExtractor(const T &obj, size_t index) {
 			return obj.scalar(index);
 		}
 
-		template<typename T, typename std::enable_if_t<typetraits::TypeInfo<T>::type ==
-														 ::librapid::detail::LibRapidType::Scalar,
-													   int> = 0>
+		template<typename T, typename std::enable_if_t<!detail::IsArrayType<T>::val, int> = 0>
 		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto scalarExtractor(const T &obj, size_t) {
 			return obj;
 		}
diff --git a/librapid/include/librapid/array/operations.hpp b/librapid/include/librapid/array/operations.hpp
index b15d67be..d17c2fc9 100644
--- a/librapid/include/librapid/array/operations.hpp
+++ b/librapid/include/librapid/array/operations.hpp
@@ -67,20 +67,37 @@
 		return std::get<0>(args).shape();                                                          \
 	}
 
+//#define LIBRAPID_BINARY_SHAPE_EXTRACTOR                                                            \
+//	template<typename First, typename Second>                                                      \
+//	LIBRAPID_NODISCARD static LIBRAPID_ALWAYS_INLINE auto getShapeImpl(                            \
+//	  const std::tuple<First, Second> &tup) {                                                      \
+//		if constexpr (TypeInfo<std::decay_t<First>>::type != detail::LibRapidType::Scalar &&       \
+//					  TypeInfo<std::decay_t<Second>>::type != detail::LibRapidType::Scalar) {      \
+//			LIBRAPID_ASSERT(std::get<0>(tup).shape() == std::get<1>(tup).shape(),                  \
+//							"Shapes must match for binary operations");                            \
+//			return std::get<0>(tup).shape();                                                       \
+//		} else if constexpr (TypeInfo<std::decay_t<First>>::type ==                                \
+//							 detail::LibRapidType::Scalar) {                                       \
+//			return std::get<1>(tup).shape();                                                       \
+//		} else {                                                                                   \
+//			return std::get<0>(tup).shape();                                                       \
+//		}                                                                                          \
+//	}                                                                                                 \
+
+
 #define LIBRAPID_BINARY_SHAPE_EXTRACTOR                                                            \
 	template<typename First, typename Second>                                                      \
 	LIBRAPID_NODISCARD static LIBRAPID_ALWAYS_INLINE auto getShapeImpl(                            \
 	  const std::tuple<First, Second> &tup) {                                                      \
-		if constexpr (TypeInfo<std::decay_t<First>>::type != detail::LibRapidType::Scalar &&       \
-					  TypeInfo<std::decay_t<Second>>::type != detail::LibRapidType::Scalar) {      \
-			LIBRAPID_ASSERT(std::get<0>(tup).shape() == std::get<1>(tup).shape(),                  \
-							"Shapes must match for binary operations");                            \
+		if constexpr (IsArrayType<std::decay_t<First>>::value) {                                   \
+			if constexpr (IsArrayType<std::decay_t<Second>>::value) {                              \
+				LIBRAPID_ASSERT(std::get<0>(tup).shape() == std::get<1>(tup).shape(),              \
+								"Shapes must match for binary operations");                        \
+				return std::get<0>(tup).shape();                                                   \
+			}                                                                                      \
 			return std::get<0>(tup).shape();                                                       \
-		} else if constexpr (TypeInfo<std::decay_t<First>>::type ==                                \
-							 detail::LibRapidType::Scalar) {                                       \
+		} else if constexpr (IsArrayType<std::decay_t<Second>>::value) {                           \
 			return std::get<1>(tup).shape();                                                       \
-		} else {                                                                                   \
-			return std::get<0>(tup).shape();                                                       \
 		}                                                                                          \
 	}                                                                                              \
                                                                                                    \
@@ -523,6 +540,12 @@ namespace librapid {
 	} // namespace typetraits
 
 	namespace detail {
+		template<typename T, LibRapidType... ValidTypes>
+		constexpr bool isType() {
+			constexpr LibRapidType t = typetraits::TypeInfo<std::decay_t<T>>::type;
+			return ((t == ValidTypes) || ...);
+		}
+
 		template<typename VAL>
 		constexpr bool isArrayOp() {
 			return (typetraits::IsArrayContainer<std::decay_t<VAL>>::value ||
@@ -531,18 +554,59 @@ namespace librapid {
 
 		template<typename LHS, typename RHS>
 		constexpr bool isArrayOpArray() {
-			return (typetraits::TypeInfo<std::decay_t<LHS>>::type != LibRapidType::Scalar) &&
-				   (typetraits::TypeInfo<std::decay_t<RHS>>::type != LibRapidType::Scalar) &&
-				   typetraits::IsLibRapidType<std::decay_t<LHS>>::value &&
-				   typetraits::IsLibRapidType<std::decay_t<RHS>>::value;
+			// return (typetraits::TypeInfo<std::decay_t<LHS>>::type != LibRapidType::Scalar) &&
+			// 	   (typetraits::TypeInfo<std::decay_t<RHS>>::type != LibRapidType::Scalar) &&
+			// 	   typetraits::IsLibRapidType<std::decay_t<LHS>>::value &&
+			// 	   typetraits::IsLibRapidType<std::decay_t<RHS>>::value;
+
+			// ArrayContainer,
+			// ArrayFunction,
+			// GeneralArrayView
+
+			constexpr bool lhsIsValid = isType<LHS,
+											   LibRapidType::ArrayContainer,
+											   LibRapidType::ArrayFunction,
+											   LibRapidType::GeneralArrayView>();
+
+			constexpr bool rhsIsValid = isType<RHS,
+											   LibRapidType::ArrayContainer,
+											   LibRapidType::ArrayFunction,
+											   LibRapidType::GeneralArrayView>();
+
+			return lhsIsValid && rhsIsValid;
 		}
 
 		template<typename LHS, typename RHS>
 		constexpr bool isArrayOpWithScalar() {
-			return (typetraits::IsLibRapidType<std::decay_t<LHS>>::value &&
-					typetraits::TypeInfo<std::decay_t<RHS>>::type == LibRapidType::Scalar) ||
-				   (typetraits::TypeInfo<std::decay_t<LHS>>::type == LibRapidType::Scalar &&
-					typetraits::IsLibRapidType<std::decay_t<RHS>>::value);
+			// // We allow operations with vectors here
+			// return (typetraits::IsLibRapidType<std::decay_t<LHS>>::value &&
+			// 		(
+			// 			(typetraits::TypeInfo<std::decay_t<RHS>>::type == LibRapidType::Scalar) ||
+			// 			(typetraits::TypeInfo<std::decay_t<RHS>>::type == LibRapidType::Vector)
+			// 		)
+			// 	   ) ||
+			// 	   (
+			// 	   	(
+			// 	   		(typetraits::TypeInfo<std::decay_t<LHS>>::type == LibRapidType::Scalar) ||
+			// 	   		(typetraits::TypeInfo<std::decay_t<LHS>>::type == LibRapidType::Vector)
+			// 	   	) &&
+			// 		typetraits::IsLibRapidType<std::decay_t<RHS>>::value);
+
+			constexpr bool lhsIsArray = isType<LHS,
+											   LibRapidType::ArrayContainer,
+											   LibRapidType::ArrayFunction,
+											   LibRapidType::GeneralArrayView>();
+
+			constexpr bool rhsIsArray = isType<RHS,
+											   LibRapidType::ArrayContainer,
+											   LibRapidType::ArrayFunction,
+											   LibRapidType::GeneralArrayView>();
+
+			constexpr bool lhsIsScalar = isType<LHS, LibRapidType::Scalar, LibRapidType::Vector>();
+
+			constexpr bool rhsIsScalar = isType<RHS, LibRapidType::Scalar, LibRapidType::Vector>();
+
+			return (lhsIsArray ^ rhsIsArray) && (lhsIsScalar ^ rhsIsScalar);
 		}
 	} // namespace detail
 
diff --git a/librapid/include/librapid/array/shape.hpp b/librapid/include/librapid/array/shape.hpp
index 6e7c5903..061e296b 100644
--- a/librapid/include/librapid/array/shape.hpp
+++ b/librapid/include/librapid/array/shape.hpp
@@ -776,34 +776,19 @@ namespace librapid {
 			using Type = VectorShape;
 		};
 
-		template<>
-		struct ShapeTypeHelperImpl<Shape, std::false_type> {
-			using Type = Shape;
+		template<typename NonFalseType>
+		struct ShapeTypeHelperImpl<NonFalseType, std::false_type> {
+			using Type = NonFalseType;
 		};
 
-		template<>
-		struct ShapeTypeHelperImpl<std::false_type, Shape> {
+		template<typename NonFalseType>
+		struct ShapeTypeHelperImpl<std::false_type, NonFalseType> {
 			using Type = Shape;
 		};
 
 		template<>
-		struct ShapeTypeHelperImpl<MatrixShape, std::false_type> {
-			using Type = MatrixShape;
-		};
-
-		template<>
-		struct ShapeTypeHelperImpl<std::false_type, MatrixShape> {
-			using Type = MatrixShape;
-		};
-
-		template<>
-		struct ShapeTypeHelperImpl<VectorShape, std::false_type> {
-			using Type = VectorShape;
-		};
-
-		template<>
-		struct ShapeTypeHelperImpl<std::false_type, VectorShape> {
-			using Type = VectorShape;
+		struct ShapeTypeHelperImpl<std::false_type, std::false_type> {
+			using Type = VectorShape; // Fastest
 		};
 
 		template<typename... Args>
diff --git a/librapid/include/librapid/core/forward.hpp b/librapid/include/librapid/core/forward.hpp
index f7594eb6..1ea1b423 100644
--- a/librapid/include/librapid/core/forward.hpp
+++ b/librapid/include/librapid/core/forward.hpp
@@ -28,6 +28,17 @@ namespace librapid {
         class ArrayContainer;
     }
 
+    namespace typetraits {
+        /// Evaluates as true if the input type is an ArrayContainer instance
+        /// \tparam T Input type
+        template<typename T>
+        struct IsArrayContainer : std::false_type {};
+
+        template<typename ShapeType, typename StorageScalar>
+        struct IsArrayContainer<array::ArrayContainer<ShapeType, StorageScalar>> : std::true_type {
+        };
+    }
+
     namespace detail {
         /// \brief Identifies which type of function is being used
         namespace descriptor {
diff --git a/librapid/include/librapid/core/typetraits.hpp b/librapid/include/librapid/core/typetraits.hpp
index c4cd8962..28964131 100644
--- a/librapid/include/librapid/core/typetraits.hpp
+++ b/librapid/include/librapid/core/typetraits.hpp
@@ -46,6 +46,12 @@ namespace librapid::typetraits {
         std::true_type testCast(int);
         template<typename From, typename To>
         std::false_type testCast(float);
+
+		// Test for T::allowVectorisation (static constexpr bool)
+		template<typename T, typename = decltype(T::allowVectorisation)>
+		std::true_type testAllowVectorisation(int);
+		template<typename T>
+		std::false_type testAllowVectorisation(float);
     } // namespace impl
 
     template<typename T, typename Index = int64_t>
@@ -63,6 +69,10 @@ namespace librapid::typetraits {
     // Detect whether a class can be default constructed
     template<class T>
     using TriviallyDefaultConstructible = std::is_trivially_default_constructible<T>;
+
+	// Detect whether a class has a static constexpr bool member called allowVectorization
+	template<typename T>
+	struct HasAllowVectorisation : public decltype(impl::testAllowVectorisation<T>(1)) {};
 } // namespace librapid::typetraits
 
 #endif // LIBRAPID_CORE_TYPETRAITS_HPP
\ No newline at end of file
diff --git a/librapid/include/librapid/datastructures/bitset.hpp b/librapid/include/librapid/datastructures/bitset.hpp
index 1dcab841..2991eda5 100644
--- a/librapid/include/librapid/datastructures/bitset.hpp
+++ b/librapid/include/librapid/datastructures/bitset.hpp
@@ -394,10 +394,10 @@ namespace librapid {
 	}
 } // namespace librapid
 
-template<size_t numBits, typename Char>
-struct fmt::formatter<librapid::BitSet<numBits>, Char> {
+template<uint64_t numBits, bool stackAlloc, typename Char>
+struct fmt::formatter<librapid::BitSet<numBits, stackAlloc>, Char> {
 private:
-	using Type = librapid::BitSet<numBits>;
+	using Type = librapid::BitSet<numBits, stackAlloc>;
 	using Base = fmt::formatter<int, Char>;
 	Base m_base;
 
@@ -414,11 +414,11 @@ struct fmt::formatter<librapid::BitSet<numBits>, Char> {
 	}
 };
 
-LIBRAPID_SIMPLE_IO_NORANGE(uint64_t numBits, librapid::BitSet<numBits>)
+LIBRAPID_SIMPLE_IO_NORANGE(uint64_t numBits COMMA bool stackAlloc, librapid::BitSet<numBits COMMA stackAlloc>)
 
 // std ostream support
-template<uint64_t numElements>
-std::ostream &operator<<(std::ostream &os, const librapid::BitSet<numElements> &bitset) {
+template<uint64_t numElements, bool stackAlloc>
+std::ostream &operator<<(std::ostream &os, const librapid::BitSet<numElements, stackAlloc> &bitset) {
 	return os << fmt::format("{}", bitset);
 }
 
diff --git a/librapid/include/librapid/math/vector.hpp b/librapid/include/librapid/math/vector.hpp
index 8dc25aac..d578470a 100644
--- a/librapid/include/librapid/math/vector.hpp
+++ b/librapid/include/librapid/math/vector.hpp
@@ -18,6 +18,41 @@ namespace librapid {
     using Vec2 = Vector<double, 2>;
     using Vec3 = Vector<double, 3>;
     using Vec4 = Vector<double, 4>;
+
+    template<typename T = float, uint64_t numDims = 3>
+    Vector<T, numDims> randomPointOnSphere() {
+        // Given X = [x_1, x_2, ..., x_n], where x_n ~ N(1, 0),
+        // X / |X| will be uniformly distributed on the unit sphere
+
+        Vector<T, numDims> result;
+        for (uint64_t i = 0; i < numDims; ++i) {
+            result[i] = ::librapid::randomGaussian();
+        }
+
+        return norm(result);
+    }
+
+    /// \brief Generate a random point within the unit sphere
+    /// \tparam T Scalar type of the vector, defaulting to ``float``
+    /// \return 3D vector, \f$\vec{v}\f$, with \f$0 \leq |\vec{v}| \leq 1\f$
+    template<typename T = float>
+    Vector<T, 3> randomPointInSphere() {
+        // Adapted from https://karthikkaranth.me/blog/generating-random-points-in-a-sphere/
+
+        auto u = ::librapid::random<T>(); // [0, 1)
+        auto v = ::librapid::random<T>(); // [0, 1)
+        auto theta = u * 2 * PI;
+        auto phi = ::librapid::acos(2 * v - 1);
+        auto r = ::librapid::cbrt(random<T>());
+        auto sinTheta = ::librapid::sin(theta);
+        auto cosTheta = ::librapid::cos(theta);
+        auto sinPhi = ::librapid::sin(phi);
+        auto cosPhi = ::librapid::cos(phi);
+        auto x = r * sinPhi * cosTheta;
+        auto y = r * sinPhi * sinTheta;
+        auto z = r * cosPhi;
+        return {x, y, z};
+    }
 } // namespace librapid
 
 #endif // LIBRAPID_MATH_VECTOR_OLD_HPP
\ No newline at end of file
diff --git a/librapid/include/librapid/math/vectorForward.hpp b/librapid/include/librapid/math/vectorForward.hpp
index a825195f..0d53203f 100644
--- a/librapid/include/librapid/math/vectorForward.hpp
+++ b/librapid/include/librapid/math/vectorForward.hpp
@@ -3,13 +3,13 @@
 
 namespace librapid {
 	namespace vectorDetail {
-		template<typename T, size_t N>
+		template<typename T, uint64_t N>
 		struct GenericVectorStorage;
 
-		template<typename T, size_t N>
+		template<typename T, uint64_t N>
 		struct SimdVectorStorage;
 
-		template<typename T, size_t N>
+		template<typename T, uint64_t N>
 		struct VectorStorageType {
 			using type = std::conditional_t<(typetraits::TypeInfo<T>::packetWidth > 1),
 											SimdVectorStorage<T, N>, GenericVectorStorage<T, N>>;
@@ -17,18 +17,23 @@ namespace librapid {
 
 		template<typename Storage0, typename Storage1>
 		auto vectorStorageTypeMerger() {
-			using Scalar0						 = typename typetraits::TypeInfo<Storage0>::Scalar;
-			using Scalar1						 = typename typetraits::TypeInfo<Storage1>::Scalar;
-			static constexpr size_t packetWidth0 = typetraits::TypeInfo<Scalar0>::packetWidth;
-			static constexpr size_t packetWidth1 = typetraits::TypeInfo<Scalar1>::packetWidth;
-			if constexpr (packetWidth0 > 1 && packetWidth1 > 1) {
+			using Scalar0 = typename typetraits::TypeInfo<Storage0>::Scalar;
+			using Scalar1 = typename typetraits::TypeInfo<Storage1>::Scalar;
+			static constexpr uint64_t packetWidth0 = typetraits::TypeInfo<Scalar0>::packetWidth;
+			static constexpr uint64_t packetWidth1 = typetraits::TypeInfo<Scalar1>::packetWidth;
+			if constexpr (typetraits::TypeInfo<Storage0>::type == detail::LibRapidType::Scalar) {
+				return Storage1 {};
+			} else if constexpr (typetraits::TypeInfo<Storage1>::type ==
+								 detail::LibRapidType::Scalar) {
+				return Storage0 {};
+			} else if constexpr (packetWidth0 > 1 && packetWidth1 > 1) {
 				return SimdVectorStorage<typename Storage0::Scalar, Storage0::dims> {};
 			} else {
 				return GenericVectorStorage<typename Storage0::Scalar, Storage0::dims> {};
 			}
 		}
 
-		template<typename T, size_t N>
+		template<typename T, uint64_t N>
 		using VectorStorage = typename VectorStorageType<T, N>::type;
 
 		template<typename Storage0, typename Storage1>
@@ -50,7 +55,9 @@ namespace librapid {
 				return static_cast<Derived &>(*this);
 			}
 
-			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto eval() const { return derived(); }
+//			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE virtual Derived eval() const {
+//				return derived();
+//			}
 
 			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE virtual IndexTypeConst
 			operator[](int64_t index) const {
@@ -84,13 +91,13 @@ namespace librapid {
 				derived().str(formatter, ctx);
 			}
 
-			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE virtual GetType _get(size_t index) const {
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE virtual GetType _get(uint64_t index) const {
 				return derived()._get(index);
 			}
 		};
 	} // namespace vectorDetail
 
-	template<typename ScalarType, size_t NumDims>
+	template<typename ScalarType, uint64_t NumDims>
 	class Vector;
 
 	namespace vectorDetail {
@@ -100,27 +107,35 @@ namespace librapid {
 		template<typename Val, typename Op>
 		struct UnaryVecOp;
 
-		template<typename Scalar, size_t N, typename LHS, typename RHS, typename Op,
+		template<typename Scalar, uint64_t N, typename LHS, typename RHS, typename Op,
 				 size_t... Indices>
 		LIBRAPID_ALWAYS_INLINE void assignImpl(Vector<Scalar, N> &dst,
 											   const BinaryVecOp<LHS, RHS, Op> &src,
 											   std::index_sequence<Indices...>);
 
-		template<typename Scalar, size_t N, typename Val, typename Op, size_t... Indices>
+		template<typename Scalar, uint64_t N, typename Val, typename Op, size_t... Indices>
 		LIBRAPID_ALWAYS_INLINE void assignImpl(Vector<Scalar, N> &dst,
 											   const UnaryVecOp<Val, Op> &src,
 											   std::index_sequence<Indices...>);
 
-		template<typename Scalar, size_t N, typename LHS, typename RHS, typename Op>
+		template<typename Scalar, uint64_t N, typename LHS, typename RHS, typename Op>
 		LIBRAPID_ALWAYS_INLINE void assign(Vector<Scalar, N> &dst,
 										   const BinaryVecOp<LHS, RHS, Op> &src);
 
-		template<typename Scalar, size_t N, typename Val, typename Op>
+		template<typename Scalar, uint64_t N, typename Val, typename Op>
 		LIBRAPID_ALWAYS_INLINE void assign(Vector<Scalar, N> &dst, const UnaryVecOp<Val, Op> &src);
 	} // namespace vectorDetail
 
-	template<typename ScalarType, size_t NumDims>
-	class Vector;
+	namespace typetraits {
+		LIBRAPID_DEFINE_AS_TYPE(typename ScalarType COMMA uint64_t NumDims,
+								Vector<ScalarType COMMA NumDims>);
+
+		LIBRAPID_DEFINE_AS_TYPE(typename LHS COMMA typename RHS COMMA typename Op,
+								vectorDetail::BinaryVecOp<LHS COMMA RHS COMMA Op>);
+
+		LIBRAPID_DEFINE_AS_TYPE(typename Val COMMA typename Op,
+								vectorDetail::UnaryVecOp<Val COMMA Op>);
+	} // namespace typetraits
 } // namespace librapid
 
 #endif // LIBRAPID_MATH_VECTOR_FORWARD_HPP
\ No newline at end of file
diff --git a/librapid/include/librapid/math/vectorImpl.hpp b/librapid/include/librapid/math/vectorImpl.hpp
index c805c3f4..bc8ff529 100644
--- a/librapid/include/librapid/math/vectorImpl.hpp
+++ b/librapid/include/librapid/math/vectorImpl.hpp
@@ -8,13 +8,13 @@ namespace librapid {
 		template<typename T>
 		struct IsVector : std::false_type {};
 
-		template<typename T, size_t N>
+		template<typename T, uint64_t N>
 		struct IsVector<vectorDetail::GenericVectorStorage<T, N>> : std::true_type {};
 
-		template<typename T, size_t N>
+		template<typename T, uint64_t N>
 		struct IsVector<vectorDetail::SimdVectorStorage<T, N>> : std::true_type {};
 
-		template<typename T, size_t N>
+		template<typename T, uint64_t N>
 		struct IsVector<Vector<T, N>> : std::true_type {};
 
 		template<typename LHS, typename RHS, typename Op>
@@ -23,20 +23,20 @@ namespace librapid {
 		template<typename VAL, typename Op>
 		struct IsVector<vectorDetail::UnaryVecOp<VAL, Op>> : std::true_type {};
 
-		template<typename T, size_t N>
+		template<typename T, uint64_t N>
 		struct TypeInfo<vectorDetail::GenericVectorStorage<T, N>> {
 			static constexpr detail::LibRapidType type = detail::LibRapidType::Vector;
 			using Scalar							   = T;
 			using IndexType							   = T &;
 			using IndexTypeConst					   = const T &;
 			using GetType							   = const T &;
+			using StorageType						   = vectorDetail::GenericVectorStorage<T, N>;
+			using ShapeType							   = std::false_type;
 
-			using StorageType = vectorDetail::GenericVectorStorage<T, N>;
-
-			static constexpr size_t length = N;
+			static constexpr uint64_t length = N;
 		};
 
-		template<typename T, size_t N>
+		template<typename T, uint64_t N>
 		struct TypeInfo<vectorDetail::SimdVectorStorage<T, N>> {
 			static constexpr detail::LibRapidType type = detail::LibRapidType::Vector;
 			using Scalar							   = T;
@@ -44,24 +44,30 @@ namespace librapid {
 			using IndexType							   = T &;
 			using IndexTypeConst					   = const T &;
 			using GetType							   = Packet;
+			using Backend							   = backend::CPU;
+			using StorageType						   = vectorDetail::SimdVectorStorage<T, N>;
+			using ShapeType							   = std::false_type;
 
-			using StorageType = vectorDetail::SimdVectorStorage<T, N>;
-
-			static constexpr size_t packetWidth = TypeInfo<T>::packetWidth;
-			static constexpr size_t length =
+			static constexpr uint64_t packetWidth = TypeInfo<T>::packetWidth;
+			static constexpr uint64_t length =
 			  (N + TypeInfo<T>::packetWidth - 1) / TypeInfo<T>::packetWidth;
 		};
 
-		template<typename ScalarType, size_t NumDims>
+		template<typename ScalarType, uint64_t NumDims>
 		struct TypeInfo<Vector<ScalarType, NumDims>> {
 			static constexpr detail::LibRapidType type = detail::LibRapidType::Vector;
 			using Scalar							   = ScalarType;
-			static constexpr size_t dims			   = NumDims;
-			using StorageType			   = vectorDetail::VectorStorage<Scalar, NumDims>;
-			static constexpr size_t length = StorageType::length;
-			using IndexTypeConst		   = typename StorageType::IndexTypeConst;
-			using IndexType				   = typename StorageType::IndexType;
-			using GetType				   = typename StorageType::GetType;
+			using Packet							   = std::false_type;
+			static constexpr uint64_t dims			   = NumDims;
+			using StorageType				 = vectorDetail::VectorStorage<Scalar, NumDims>;
+			static constexpr uint64_t length = StorageType::length;
+			using IndexTypeConst			 = typename StorageType::IndexTypeConst;
+			using IndexType					 = typename StorageType::IndexType;
+			using GetType					 = typename StorageType::GetType;
+			using Backend					 = backend::CPU;
+			using ShapeType					 = std::false_type;
+
+			static constexpr bool allowVectorisation = false;
 		};
 
 		template<typename LHS, typename RHS, typename Op>
@@ -70,71 +76,81 @@ namespace librapid {
 			using ScalarLHS							   = typename typetraits::TypeInfo<LHS>::Scalar;
 			using ScalarRHS							   = typename typetraits::TypeInfo<RHS>::Scalar;
 			using Scalar = decltype(Op()(std::declval<ScalarLHS>(), std::declval<ScalarRHS>()));
-			using IndexTypeConst		   = Scalar;
-			using IndexType				   = Scalar;
-			using StorageType			   = typename vectorDetail::VectorStorageMerger<LHS, RHS>;
-			static constexpr size_t dims   = StorageType::dims;
-			static constexpr size_t length = StorageType::length;
-			using GetType				   = typename std::decay_t<typename StorageType::GetType>;
+			using Packet = std::false_type;
+			using IndexTypeConst			 = Scalar;
+			using IndexType					 = Scalar;
+			using StorageType				 = typename vectorDetail::VectorStorageMerger<LHS, RHS>;
+			static constexpr uint64_t dims	 = StorageType::dims;
+			static constexpr uint64_t length = StorageType::length;
+			using GetType					 = typename std::decay_t<typename StorageType::GetType>;
+			using Backend					 = backend::CPU;
+			using ShapeType					 = std::false_type;
+
+			static constexpr bool allowVectorisation = false;
 		};
 
 		template<typename Val, typename Op>
 		struct TypeInfo<vectorDetail::UnaryVecOp<Val, Op>> {
 			static constexpr detail::LibRapidType type = detail::LibRapidType::Vector;
 			using Scalar							   = typename typetraits::TypeInfo<Val>::Scalar;
+			using Packet							   = std::false_type;
 			using IndexTypeConst					   = Scalar;
 			using IndexType							   = Scalar;
-			using StorageType			 = typename vectorDetail::VectorStorage<Scalar, Val::dims>;
-			static constexpr size_t dims = StorageType::dims;
-			static constexpr size_t length = StorageType::length;
-			using GetType				   = typename std::decay_t<typename StorageType::GetType>;
+			using StorageType = typename vectorDetail::VectorStorage<Scalar, Val::dims>;
+			static constexpr uint64_t dims	 = StorageType::dims;
+			static constexpr uint64_t length = StorageType::length;
+			using GetType					 = typename std::decay_t<typename StorageType::GetType>;
+			using Backend					 = backend::CPU;
+
+			static constexpr bool allowVectorisation = false;
 		};
 	} // namespace typetraits
 
 	namespace vectorDetail {
-		template<typename T, size_t N, typename... Args, size_t... Indices>
+		template<typename T, uint64_t N, typename... Args, size_t... Indices>
 		LIBRAPID_ALWAYS_INLINE void vectorStorageAssigner(std::index_sequence<Indices...>,
 														  GenericVectorStorage<T, N> &dst,
 														  const Args &...args) {
 			((dst[Indices] = args), ...);
 		}
 
-		template<typename T, size_t N, typename... Args, size_t... Indices>
+		template<typename T, uint64_t N, typename... Args, size_t... Indices>
 		LIBRAPID_ALWAYS_INLINE void vectorStorageAssigner(std::index_sequence<Indices...>,
 														  SimdVectorStorage<T, N> &dst,
 														  const Args &...args) {
 			((dst.data.scalar[Indices] = args), ...);
 		}
 
-		template<typename T, size_t N, typename T2, size_t N2, size_t... Indices>
+		template<typename T, uint64_t N, typename T2, uint64_t N2, size_t... Indices>
 		LIBRAPID_ALWAYS_INLINE void vectorStorageAssigner(std::index_sequence<Indices...>,
 														  GenericVectorStorage<T, N> &dst,
 														  const GenericVectorStorage<T2, N2> &src) {
 			((dst[Indices] = src[Indices]), ...);
 		}
 
-		template<typename T, size_t N, typename T2, size_t N2, size_t... Indices>
+		template<typename T, uint64_t N, typename T2, uint64_t N2, size_t... Indices>
 		LIBRAPID_ALWAYS_INLINE void vectorStorageAssigner(std::index_sequence<Indices...>,
 														  GenericVectorStorage<T, N> &dst,
 														  const SimdVectorStorage<T2, N2> &src) {
 			((dst[Indices] = src.data.scalar[Indices]), ...);
 		}
 
-		template<typename T, size_t N, typename T2, size_t N2, size_t... Indices>
+		template<typename T, uint64_t N, typename T2, uint64_t N2, size_t... Indices>
 		LIBRAPID_ALWAYS_INLINE void vectorStorageAssigner(std::index_sequence<Indices...>,
 														  SimdVectorStorage<T, N> &dst,
 														  const GenericVectorStorage<T2, N2> &src) {
 			((dst.data.scalar[Indices] = src[Indices]), ...);
 		}
 
-		template<typename T, size_t N, typename T2, size_t N2, size_t... Indices>
+		template<typename T, uint64_t N, typename T2, uint64_t N2, size_t... Indices>
 		LIBRAPID_ALWAYS_INLINE void
-		vectorStorageAssigner_simdHelper(std::index_sequence<Indices...>, SimdVectorStorage<T, N> &dst,
-								   const SimdVectorStorage<T2, N2> &src) {
+		vectorStorageAssigner_simdHelper(std::index_sequence<Indices...>,
+										 SimdVectorStorage<T, N> &dst,
+										 const SimdVectorStorage<T2, N2> &src) {
 			((dst.data.simd[Indices] = src.data.simd[Indices]), ...);
 		}
 
-		template<typename T, size_t N, typename T2, size_t N2, size_t... Indices>
+		template<typename T, uint64_t N, typename T2, uint64_t N2, size_t... Indices>
 		LIBRAPID_ALWAYS_INLINE void vectorStorageAssigner(std::index_sequence<Indices...>,
 														  SimdVectorStorage<T, N> &dst,
 														  const SimdVectorStorage<T2, N2> &src) {
@@ -144,19 +160,20 @@ namespace librapid {
 			// packets.
 
 			if constexpr (std::is_same_v<T, T2>) {
-				constexpr size_t packetWidth = typetraits::TypeInfo<T>::packetWidth;
-				constexpr size_t length		 = (N + packetWidth - 1) / packetWidth;
-				vectorStorageAssigner_simdHelper(std::make_index_sequence<length>(), dst, src);
+				constexpr uint64_t packetWidth = typetraits::TypeInfo<T>::packetWidth;
+				constexpr uint64_t length	   = (N + packetWidth - 1) / packetWidth;
+				vectorStorageAssigner_simdHelper(
+				  std::make_index_sequence<static_cast<size_t>(length)>(), dst, src);
 			} else {
 				((dst.data.scalar[Indices] = src.data.scalar[Indices]), ...);
 			}
 		}
 
-		template<typename ScalarType, size_t NumDims>
+		template<typename ScalarType, uint64_t NumDims>
 		struct GenericVectorStorage {
-			using Scalar				   = ScalarType;
-			static constexpr size_t dims   = NumDims;
-			static constexpr size_t length = typetraits::TypeInfo<GenericVectorStorage>::length;
+			using Scalar					 = ScalarType;
+			static constexpr uint64_t dims	 = NumDims;
+			static constexpr uint64_t length = typetraits::TypeInfo<GenericVectorStorage>::length;
 			using IndexType = typename typetraits::TypeInfo<GenericVectorStorage>::IndexType;
 			using IndexTypeConst =
 			  typename typetraits::TypeInfo<GenericVectorStorage>::IndexTypeConst;
@@ -177,7 +194,7 @@ namespace librapid {
 
 			template<typename T>
 			GenericVectorStorage(const T &other) {
-				for (size_t i = 0; i < length; ++i) { data[i] = other[i]; }
+				for (uint64_t i = 0; i < length; ++i) { data[i] = other[i]; }
 			}
 
 			template<typename T>
@@ -186,8 +203,10 @@ namespace librapid {
 								"Initializer list for Vector is too long ({} > {})",
 								other.size(),
 								dims);
-				const size_t minDims = (other.size() < dims) ? other.size() : dims;
-				for (size_t i = 0; i < minDims; ++i) { this->operator[](i) = *(other.begin() + i); }
+				const uint64_t minDims = (other.size() < dims) ? other.size() : dims;
+				for (uint64_t i = 0; i < minDims; ++i) {
+					this->operator[](i) = *(other.begin() + i);
+				}
 			}
 
 			template<typename T>
@@ -196,8 +215,8 @@ namespace librapid {
 								"Initializer list for Vector is too long ({} > {})",
 								other.size(),
 								dims);
-				const size_t minDims = (other.size() < dims) ? other.size() : dims;
-				for (size_t i = 0; i < minDims; ++i) { this->operator[](i) = other[i]; }
+				const uint64_t minDims = (other.size() < dims) ? other.size() : dims;
+				for (uint64_t i = 0; i < minDims; ++i) { this->operator[](i) = other[i]; }
 			}
 
 			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE IndexTypeConst
@@ -219,17 +238,17 @@ namespace librapid {
 
 			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Scalar sum() const {
 				Scalar sum = Scalar(0);
-				for (size_t i = 0; i < dims; ++i) { sum += data[i]; }
+				for (uint64_t i = 0; i < dims; ++i) { sum += data[i]; }
 				return sum;
 			}
 
 			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Scalar sum2() const {
 				Scalar sum = Scalar(0);
-				for (size_t i = 0; i < dims; ++i) { sum += data[i] * data[i]; }
+				for (uint64_t i = 0; i < dims; ++i) { sum += data[i] * data[i]; }
 				return sum;
 			}
 
-			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE const Scalar &_get(size_t index) const {
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE const Scalar &_get(uint64_t index) const {
 				LIBRAPID_ASSERT(index >= 0 && index < dims,
 								"Index {} out of bounds for Vector of length {}",
 								index,
@@ -237,7 +256,7 @@ namespace librapid {
 				return data[index];
 			}
 
-			LIBRAPID_ALWAYS_INLINE void _set(size_t index, const Scalar &value) {
+			LIBRAPID_ALWAYS_INLINE void _set(uint64_t index, const Scalar &value) {
 				LIBRAPID_ASSERT(index >= 0 && index < dims,
 								"Index {} out of bounds for Vector of length {}",
 								index,
@@ -250,12 +269,12 @@ namespace librapid {
 		/// equal to \f$(N + L - 1) / L\f$, where \f$L\f$ is the packet width of the scalar type.
 		/// \tparam Scalar_ The scalar type to store
 		/// \tparam N The number of scalars to store
-		template<typename Scalar_, size_t N>
+		template<typename Scalar_, uint64_t N>
 		union ScalarToSimd {
-			using Scalar						= Scalar_;
-			using Packet						= typename typetraits::TypeInfo<Scalar>::Packet;
-			static constexpr size_t packetWidth = typetraits::TypeInfo<Scalar>::packetWidth;
-			static constexpr size_t length		= (N + packetWidth - 1) / packetWidth;
+			using Scalar						  = Scalar_;
+			using Packet						  = typename typetraits::TypeInfo<Scalar>::Packet;
+			static constexpr uint64_t packetWidth = typetraits::TypeInfo<Scalar>::packetWidth;
+			static constexpr uint64_t length	  = (N + packetWidth - 1) / packetWidth;
 
 #if defined(LIBRAPID_NATIVE_ARCH)
 			alignas(LIBRAPID_MEM_ALIGN) std::array<Scalar, N> scalar;
@@ -266,13 +285,13 @@ namespace librapid {
 #endif
 		};
 
-		template<typename ScalarType, size_t NumDims>
+		template<typename ScalarType, uint64_t NumDims>
 		struct SimdVectorStorage {
-			using Scalar						= ScalarType;
-			static constexpr size_t dims		= NumDims;
-			using Packet						= typename typetraits::TypeInfo<Scalar>::Packet;
-			static constexpr size_t packetWidth = typetraits::TypeInfo<Scalar>::packetWidth;
-			static constexpr size_t length		= (dims + packetWidth - 1) / packetWidth;
+			using Scalar						  = ScalarType;
+			static constexpr uint64_t dims		  = NumDims;
+			using Packet						  = typename typetraits::TypeInfo<Scalar>::Packet;
+			static constexpr uint64_t packetWidth = typetraits::TypeInfo<Scalar>::packetWidth;
+			static constexpr uint64_t length	  = (dims + packetWidth - 1) / packetWidth;
 
 			using IndexType		 = typename typetraits::TypeInfo<SimdVectorStorage>::IndexType;
 			using IndexTypeConst = typename typetraits::TypeInfo<SimdVectorStorage>::IndexTypeConst;
@@ -285,9 +304,9 @@ namespace librapid {
 
 			template<typename... Args>
 			explicit SimdVectorStorage(Args... args) {
-				constexpr size_t minLength = (sizeof...(Args) < dims) ? sizeof...(Args) : dims;
+				constexpr uint64_t minLength = (sizeof...(Args) < dims) ? sizeof...(Args) : dims;
 				vectorDetail::vectorStorageAssigner(
-				  std::make_index_sequence<minLength>(), *this, args...);
+				  std::make_index_sequence<static_cast<size_t>(minLength)>(), *this, args...);
 			}
 
 			template<typename T>
@@ -296,8 +315,8 @@ namespace librapid {
 								"Initializer list for Vector is too long ({} > {})",
 								other.size(),
 								dims);
-				const size_t minDims = (other.size() < dims) ? other.size() : dims;
-				for (size_t i = 0; i < minDims; ++i) { data.scalar[i] = *(other.begin() + i); }
+				const uint64_t minDims = (other.size() < dims) ? other.size() : dims;
+				for (uint64_t i = 0; i < minDims; ++i) { data.scalar[i] = *(other.begin() + i); }
 			}
 
 			template<typename T>
@@ -306,8 +325,8 @@ namespace librapid {
 								"Initializer list for Vector is too long ({} > {})",
 								other.size(),
 								dims);
-				const size_t minDims = (other.size() < dims) ? other.size() : dims;
-				for (size_t i = 0; i < minDims; ++i) { data.scalar[i] = other[i]; }
+				const uint64_t minDims = (other.size() < dims) ? other.size() : dims;
+				for (uint64_t i = 0; i < minDims; ++i) { data.scalar[i] = other[i]; }
 			}
 
 			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE IndexTypeConst
@@ -329,19 +348,19 @@ namespace librapid {
 
 			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto sum() const -> Scalar {
 				Packet sum = Packet(0);
-				for (size_t i = 0; i < length; ++i) { sum += data.simd[i]; }
+				for (uint64_t i = 0; i < length; ++i) { sum += data.simd[i]; }
 				// return sum.sum();
 				return xsimd::reduce_add(sum);
 			}
 
 			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto sum2() const -> Scalar {
 				Packet sum = Packet(0);
-				for (size_t i = 0; i < length; ++i) { sum += data.simd[i] * data.simd[i]; }
+				for (uint64_t i = 0; i < length; ++i) { sum += data.simd[i] * data.simd[i]; }
 				// return sum.sum();
 				return xsimd::reduce_add(sum);
 			}
 
-			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE const Packet &_get(size_t index) const {
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE const Packet &_get(uint64_t index) const {
 				LIBRAPID_ASSERT(index >= 0 && index < dims,
 								"Index {} out of bounds for Vector of length {}",
 								index,
@@ -349,7 +368,7 @@ namespace librapid {
 				return data.simd[index];
 			}
 
-			LIBRAPID_ALWAYS_INLINE void _set(size_t index, const Packet &value) {
+			LIBRAPID_ALWAYS_INLINE void _set(uint64_t index, const Packet &value) {
 				LIBRAPID_ASSERT(index >= 0 && index < dims,
 								"Index {} out of bounds for Vector of length {}",
 								index,
@@ -360,25 +379,25 @@ namespace librapid {
 
 		template<typename T>
 		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE constexpr auto
-		scalarSubscriptHelper(const T &val, size_t index) {
+		scalarSubscriptHelper(const T &val, uint64_t index) {
 			return val;
 		}
 
-		template<typename ScalarType, size_t NumDims>
+		template<typename ScalarType, uint64_t NumDims>
 		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE constexpr auto
-		scalarSubscriptHelper(const Vector<ScalarType, NumDims> &val, size_t index) {
+		scalarSubscriptHelper(const Vector<ScalarType, NumDims> &val, uint64_t index) {
 			return val[index];
 		}
 
 		template<typename LHS, typename RHS, typename Op>
 		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE constexpr auto
-		scalarSubscriptHelper(const BinaryVecOp<LHS, RHS, Op> &val, size_t index) {
+		scalarSubscriptHelper(const BinaryVecOp<LHS, RHS, Op> &val, uint64_t index) {
 			return val[index];
 		}
 
 		template<typename Val, typename Op>
 		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE constexpr auto
-		scalarSubscriptHelper(const UnaryVecOp<Val, Op> &val, size_t index) {
+		scalarSubscriptHelper(const UnaryVecOp<Val, Op> &val, uint64_t index) {
 			return val[index];
 		}
 
@@ -388,13 +407,13 @@ namespace librapid {
 			return val;
 		}
 
-		template<typename ScalarType, size_t NumDims>
+		template<typename ScalarType, uint64_t NumDims>
 		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto constexpr scalarGetHelper(
 		  const Vector<ScalarType, NumDims> &val, size_t index) {
 			return val._get(index);
 		}
 
-		template<typename ScalarType, size_t NumDims>
+		template<typename ScalarType, uint64_t NumDims>
 		LIBRAPID_NODISCARD
 		  LIBRAPID_ALWAYS_INLINE auto constexpr scalarGetHelper(Vector<ScalarType, NumDims> &val,
 																size_t index) {
@@ -419,7 +438,7 @@ namespace librapid {
 			using type = std::false_type;
 		};
 
-		template<typename ScalarType, size_t NumDims>
+		template<typename ScalarType, uint64_t NumDims>
 		struct VectorScalarStorageExtractor<Vector<ScalarType, NumDims>> {
 			using type = typename typetraits::TypeInfo<Vector<ScalarType, NumDims>>::StorageType;
 		};
@@ -436,35 +455,35 @@ namespace librapid {
 
 		template<typename T>
 		struct VectorScalarDimensionExtractor {
-			static constexpr size_t value = 0;
+			static constexpr uint64_t value = 0;
 		};
 
-		template<typename ScalarType, size_t NumDims>
+		template<typename ScalarType, uint64_t NumDims>
 		struct VectorScalarDimensionExtractor<Vector<ScalarType, NumDims>> {
-			static constexpr size_t value = NumDims;
+			static constexpr uint64_t value = NumDims;
 		};
 
 		template<typename LHS, typename RHS, typename Op>
 		struct VectorScalarDimensionExtractor<BinaryVecOp<LHS, RHS, Op>> {
-			static constexpr size_t value = BinaryVecOp<LHS, RHS, Op>::dims;
+			static constexpr uint64_t value = BinaryVecOp<LHS, RHS, Op>::dims;
 		};
 
 		template<typename Val, typename Op>
 		struct VectorScalarDimensionExtractor<UnaryVecOp<Val, Op>> {
-			static constexpr size_t value = UnaryVecOp<Val, Op>::dims;
+			static constexpr uint64_t value = UnaryVecOp<Val, Op>::dims;
 		};
 	} // namespace vectorDetail
 
-	template<typename ScalarType, size_t NumDims>
+	template<typename ScalarType, uint64_t NumDims>
 	class Vector : public vectorDetail::VectorBase<Vector<ScalarType, NumDims>> {
 	public:
-		using Scalar				   = ScalarType;
-		static constexpr size_t dims   = NumDims;
-		using StorageType			   = vectorDetail::VectorStorage<Scalar, NumDims>;
-		static constexpr size_t length = StorageType::length;
-		using IndexTypeConst		   = typename StorageType::IndexTypeConst;
-		using IndexType				   = typename StorageType::IndexType;
-		using GetType				   = typename StorageType::GetType;
+		using Scalar					 = ScalarType;
+		static constexpr uint64_t dims	 = NumDims;
+		using StorageType				 = vectorDetail::VectorStorage<Scalar, NumDims>;
+		static constexpr uint64_t length = StorageType::length;
+		using IndexTypeConst			 = typename StorageType::IndexTypeConst;
+		using IndexType					 = typename StorageType::IndexType;
+		using GetType					 = typename StorageType::GetType;
 
 		Vector()						= default;
 		Vector(const Vector &other)		= default;
@@ -479,7 +498,7 @@ namespace librapid {
 		template<typename T>
 		explicit Vector(const std::vector<T> &args) : m_data(args) {}
 
-		template<typename OtherScalar, size_t OtherDims>
+		template<typename OtherScalar, uint64_t OtherDims>
 		explicit Vector(const Vector<OtherScalar, OtherDims> &other) {
 			*this = other.template cast<Scalar, dims>();
 		}
@@ -498,27 +517,29 @@ namespace librapid {
 
 		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE static auto one() -> Vector {
 			Vector ret;
-			for (size_t i = 0; i < dims; ++i) { ret[i] = Scalar(1); }
+			for (uint64_t i = 0; i < dims; ++i) { ret[i] = Scalar(1); }
 			return ret;
 		}
 
 		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE static auto full(Scalar val) -> Vector {
 			Vector ret;
-			for (size_t i = 0; i < dims; ++i) { ret[i] = val; }
+			for (uint64_t i = 0; i < dims; ++i) { ret[i] = val; }
 			return ret;
 		}
 
 		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE static auto random(Scalar lower = 0,
 																	 Scalar upper = 1) {
 			Vector ret;
-			for (size_t i = 0; i < dims; ++i) { ret[i] = ::librapid::random<Scalar>(lower, upper); }
+			for (uint64_t i = 0; i < dims; ++i) {
+				ret[i] = ::librapid::random<Scalar>(lower, upper);
+			}
 			return ret;
 		}
 
 		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE static auto
 		random(const Vector &lower = Vector::zero(), const Vector &upper = Vector::one()) {
 			Vector ret;
-			for (size_t i = 0; i < dims; ++i) {
+			for (uint64_t i = 0; i < dims; ++i) {
 				ret[i] = ::librapid::random<Scalar>(lower[i], upper[i]);
 			}
 			return ret;
@@ -538,7 +559,7 @@ namespace librapid {
 		auto operator=(const Vector &other) -> Vector	  & = default;
 		auto operator=(Vector &&other) noexcept -> Vector & = default;
 
-		template<typename OtherScalar, size_t OtherDims>
+		template<typename OtherScalar, uint64_t OtherDims>
 		auto operator=(const Vector<OtherScalar, OtherDims> &other) -> Vector & {
 			*this = other.template cast<Scalar, dims>();
 			return *this;
@@ -567,17 +588,17 @@ namespace librapid {
 
 		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Vector eval() const { return *this; }
 
-		template<typename NewScalar, size_t NewDims>
+		template<typename NewScalar, uint64_t NewDims>
 		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto cast() const {
-			using NewVectorType		 = Vector<NewScalar, NewDims>;
-			constexpr size_t minDims = (NewVectorType::dims < dims) ? NewVectorType::dims : dims;
+			using NewVectorType		   = Vector<NewScalar, NewDims>;
+			constexpr uint64_t minDims = (NewVectorType::dims < dims) ? NewVectorType::dims : dims;
 			NewVectorType ret;
 			vectorDetail::vectorStorageAssigner(
-			  std::make_index_sequence<minDims>(), ret.storage(), m_data);
+			  std::make_index_sequence<static_cast<size_t>(minDims)>(), ret.storage(), m_data);
 			return ret;
 		}
 
-		template<typename NewScalar, size_t NewDims>
+		template<typename NewScalar, uint64_t NewDims>
 		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE operator Vector<NewScalar, NewDims>() const {
 			return cast<NewScalar, NewDims>();
 		}
@@ -601,7 +622,7 @@ namespace librapid {
 
 		//		LIBRAPID_NODISCARD std::string str(const std::string &format) const override {
 		//			std::string ret = "(";
-		//			for (size_t i = 0; i < dims; ++i) {
+		//			for (uint64_t i = 0; i < dims; ++i) {
 		//				ret += fmt::format(format, m_data[i]);
 		//				if (i != dims - 1) { ret += ", "; }
 		//			}
@@ -612,7 +633,7 @@ namespace librapid {
 		template<typename T_, typename Char, typename Ctx>
 		void str(const fmt::formatter<T_, Char> &formatter, Ctx &ctx) const {
 			fmt::format_to(ctx.out(), "(");
-			for (size_t i = 0; i < dims; ++i) {
+			for (uint64_t i = 0; i < dims; ++i) {
 				formatter.format(m_data[i], ctx);
 				if (i != dims - 1) { fmt::format_to(ctx.out(), ", "); }
 			}
@@ -624,11 +645,11 @@ namespace librapid {
 		}
 		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE StorageType &storage() { return m_data; }
 
-		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE GetType _get(size_t index) const override {
+		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE GetType _get(uint64_t index) const override {
 			return m_data._get(index);
 		}
 
-		LIBRAPID_ALWAYS_INLINE void _set(size_t index, const GetType &value) {
+		LIBRAPID_ALWAYS_INLINE void _set(uint64_t index, const GetType &value) {
 			m_data._set(index, value);
 		}
 
@@ -639,12 +660,12 @@ namespace librapid {
 	namespace vectorDetail {
 		template<typename LHS, typename RHS, typename Op>
 		struct BinaryVecOp : public VectorBase<BinaryVecOp<LHS, RHS, Op>> {
-			using Scalar				   = typename typetraits::TypeInfo<BinaryVecOp>::Scalar;
-			using StorageLHS			   = typename VectorScalarStorageExtractor<LHS>::type;
-			using StorageRHS			   = typename VectorScalarStorageExtractor<RHS>::type;
-			using StorageType			   = VectorStorageMerger<StorageLHS, StorageRHS>;
-			static constexpr size_t dims   = StorageType::dims;
-			static constexpr size_t length = StorageType::length;
+			using Scalar					 = typename typetraits::TypeInfo<BinaryVecOp>::Scalar;
+			using StorageLHS				 = typename VectorScalarStorageExtractor<LHS>::type;
+			using StorageRHS				 = typename VectorScalarStorageExtractor<RHS>::type;
+			using StorageType				 = VectorStorageMerger<StorageLHS, StorageRHS>;
+			static constexpr uint64_t dims	 = StorageType::dims;
+			static constexpr uint64_t length = StorageType::length;
 			using IndexTypeConst = typename typetraits::TypeInfo<BinaryVecOp>::IndexTypeConst;
 			using IndexType		 = typename typetraits::TypeInfo<BinaryVecOp>::IndexType;
 			using GetType		 = typename typetraits::TypeInfo<BinaryVecOp>::GetType;
@@ -668,45 +689,40 @@ namespace librapid {
 				return result;
 			}
 
-			template<typename NewScalar, size_t NewDims>
+			template<typename NewScalar = Scalar, uint64_t NewDims = dims>
 			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto cast() const {
 				return eval().template cast<NewScalar, NewDims>();
 			}
 
-			template<typename NewScalar, size_t NewDims>
+			template<typename NewScalar = Scalar, uint64_t NewDims = dims>
 			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE operator Vector<NewScalar, NewDims>() const {
 				return cast<NewScalar, NewDims>();
 			}
 
-			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE IndexType
-			operator[](int64_t index) const override {
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE IndexType operator[](int64_t index) const {
 				return op(scalarSubscriptHelper(left, index), scalarSubscriptHelper(right, index));
 			}
 
-			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE IndexType operator[](int64_t index) override {
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE IndexType operator[](int64_t index) {
 				return op(scalarSubscriptHelper(left, index), scalarSubscriptHelper(right, index));
 			}
 
-			//			LIBRAPID_NODISCARD std::string str(const std::string &format) const override
-			//{ 				return eval().str(format);
-			//			}
-
 			template<typename T_, typename Char, typename Ctx>
 			void str(const fmt::formatter<T_, Char> &formatter, Ctx &ctx) const {
 				eval().str(formatter, ctx);
 			}
 
-			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE GetType _get(size_t index) const override {
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE GetType _get(uint64_t index) const {
 				return op(scalarGetHelper(left, index), scalarGetHelper(right, index));
 			}
 		};
 
 		template<typename Val, typename Op>
 		struct UnaryVecOp : public VectorBase<UnaryVecOp<Val, Op>> {
-			using Scalar				   = typename typetraits::TypeInfo<UnaryVecOp>::Scalar;
-			using StorageType			   = typename VectorScalarStorageExtractor<Val>::type;
-			static constexpr size_t dims   = StorageType::dims;
-			static constexpr size_t length = StorageType::length;
+			using Scalar					 = typename typetraits::TypeInfo<UnaryVecOp>::Scalar;
+			using StorageType				 = typename VectorScalarStorageExtractor<Val>::type;
+			static constexpr uint64_t dims	 = StorageType::dims;
+			static constexpr uint64_t length = StorageType::length;
 			using IndexTypeConst = typename typetraits::TypeInfo<UnaryVecOp>::IndexTypeConst;
 			using IndexType		 = typename typetraits::TypeInfo<UnaryVecOp>::IndexType;
 			using GetType		 = typename typetraits::TypeInfo<UnaryVecOp>::GetType;
@@ -728,12 +744,12 @@ namespace librapid {
 				return result;
 			}
 
-			template<typename NewScalar, size_t NewDims>
+			template<typename NewScalar, uint64_t NewDims>
 			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto cast() const {
 				return eval().template cast<NewScalar, NewDims>();
 			}
 
-			template<typename NewScalar, size_t NewDims>
+			template<typename NewScalar, uint64_t NewDims>
 			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE operator Vector<NewScalar, NewDims>() const {
 				return cast<NewScalar, NewDims>();
 			}
@@ -756,12 +772,12 @@ namespace librapid {
 				eval().str(formatter, ctx);
 			}
 
-			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE GetType _get(size_t index) const override {
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE GetType _get(uint64_t index) const override {
 				return op(scalarGetHelper(val, index));
 			}
 		};
 
-		template<typename Scalar, size_t N, typename LHS, typename RHS, typename Op,
+		template<typename Scalar, uint64_t N, typename LHS, typename RHS, typename Op,
 				 size_t... Indices>
 		LIBRAPID_ALWAYS_INLINE void assignImpl(Vector<Scalar, N> &dst,
 											   const BinaryVecOp<LHS, RHS, Op> &src,
@@ -772,37 +788,37 @@ namespace librapid {
 			 ...);
 		}
 
-		template<typename Scalar, size_t N, typename Val, typename Op, size_t... Indices>
+		template<typename Scalar, uint64_t N, typename Val, typename Op, size_t... Indices>
 		LIBRAPID_ALWAYS_INLINE void assignImpl(Vector<Scalar, N> &dst,
 											   const UnaryVecOp<Val, Op> &src,
 											   std::index_sequence<Indices...>) {
 			((dst._set(Indices, src.op(scalarGetHelper(src.val, Indices)))), ...);
 		}
 
-		template<typename Scalar, size_t N, typename LHS, typename RHS, typename Op>
+		template<typename Scalar, uint64_t N, typename LHS, typename RHS, typename Op>
 		LIBRAPID_ALWAYS_INLINE void assign(Vector<Scalar, N> &dst,
 										   const BinaryVecOp<LHS, RHS, Op> &src) {
 			using ScalarDst = typename typetraits::TypeInfo<Vector<Scalar, N>>::Scalar;
 			using ScalarSrc = typename typetraits::TypeInfo<BinaryVecOp<LHS, RHS, Op>>::Scalar;
 			if constexpr (std::is_same_v<ScalarDst, ScalarSrc>) {
-				constexpr size_t lengthDst = Vector<Scalar, N>::length;
-				constexpr size_t lengthSrc = BinaryVecOp<LHS, RHS, Op>::length;
-				constexpr size_t minLength = (lengthDst < lengthSrc) ? lengthDst : lengthSrc;
-				assignImpl(dst, src, std::make_index_sequence<minLength>());
+				constexpr uint64_t lengthDst = Vector<Scalar, N>::length;
+				constexpr uint64_t lengthSrc = BinaryVecOp<LHS, RHS, Op>::length;
+				constexpr uint64_t minLength = (lengthDst < lengthSrc) ? lengthDst : lengthSrc;
+				assignImpl(dst, src, std::make_index_sequence<static_cast<size_t>(minLength)>());
 			} else {
 				dst = src.template cast<Scalar, N>();
 			}
 		}
 
-		template<typename Scalar, size_t N, typename Val, typename Op>
+		template<typename Scalar, uint64_t N, typename Val, typename Op>
 		LIBRAPID_ALWAYS_INLINE void assign(Vector<Scalar, N> &dst, const UnaryVecOp<Val, Op> &src) {
 			using ScalarDst = typename typetraits::TypeInfo<Vector<Scalar, N>>::Scalar;
 			using ScalarSrc = typename typetraits::TypeInfo<UnaryVecOp<Val, Op>>::Scalar;
 			if constexpr (std::is_same_v<ScalarDst, ScalarSrc>) {
-				constexpr size_t lengthDst = Vector<Scalar, N>::length;
-				constexpr size_t lengthSrc = UnaryVecOp<Val, Op>::length;
-				constexpr size_t minLength = (lengthDst < lengthSrc) ? lengthDst : lengthSrc;
-				assignImpl(dst, src, std::make_index_sequence<minLength>());
+				constexpr uint64_t lengthDst = Vector<Scalar, N>::length;
+				constexpr uint64_t lengthSrc = UnaryVecOp<Val, Op>::length;
+				constexpr uint64_t minLength = (lengthDst < lengthSrc) ? lengthDst : lengthSrc;
+				assignImpl(dst, src, std::make_index_sequence<static_cast<size_t>(minLength)>());
 			} else {
 				dst = src.template cast<Scalar, N>();
 			}
@@ -819,22 +835,22 @@ namespace librapid {
 		// 	return static_cast<Scalar>(val);
 		// }
 
-		template<typename NewScalar, size_t NewDims, typename T>
+		template<typename NewScalar, uint64_t NewDims, typename T>
 		constexpr auto scalarVectorCaster(const T &val) {
 			return static_cast<NewScalar>(val);
 		}
 
-		template<typename NewScalar, size_t NewDims, typename ScalarType, size_t NumDims>
+		template<typename NewScalar, uint64_t NewDims, typename ScalarType, uint64_t NumDims>
 		constexpr auto scalarVectorCaster(const Vector<ScalarType, NumDims> &val) {
 			return val.template cast<NewScalar, NewDims>();
 		}
 
-		template<typename NewScalar, size_t NewDims, typename LHS, typename RHS, typename Op>
+		template<typename NewScalar, uint64_t NewDims, typename LHS, typename RHS, typename Op>
 		constexpr auto scalarVectorCaster(const BinaryVecOp<LHS, RHS, Op> &val) {
 			return val.template cast<NewScalar, NewDims>();
 		}
 
-		template<typename NewScalar, size_t NewDims, typename Val, typename Op>
+		template<typename NewScalar, uint64_t NewDims, typename Val, typename Op>
 		constexpr auto scalarVectorCaster(const UnaryVecOp<Val, Op> &val) {
 			return val.template cast<NewScalar, NewDims>();
 		}
@@ -850,8 +866,10 @@ namespace librapid {
                                                                                                    \
 	template<typename LHS,                                                                         \
 			 typename RHS,                                                                         \
-			 typename std::enable_if_t<typetraits::IsVector<LHS>::value ||                         \
-										 typetraits::IsVector<RHS>::value,                         \
+			 typename std::enable_if_t<(typetraits::IsVector<LHS>::value ||                        \
+										typetraits::IsVector<RHS>::value) &&                       \
+										 (!typetraits::IsArrayContainer<LHS>::value &&             \
+										  !typetraits::IsArrayContainer<RHS>::value),              \
 									   int> = 0>                                                   \
 	LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto operator OP_(const LHS &lhs, const RHS &rhs) {  \
 		using namespace ::librapid::vectorDetail;                                                  \
@@ -863,9 +881,9 @@ namespace librapid {
 			return BinaryVecOp<LHS, RHS, Op> {lhs, rhs, Op {}};                                    \
 		} else {                                                                                   \
 			using Scalar = decltype(std::declval<ScalarLeft>() + std::declval<ScalarRight>());     \
-			constexpr size_t dimsLhs = VectorScalarDimensionExtractor<LHS>::value;                 \
-			constexpr size_t dimsRhs = VectorScalarDimensionExtractor<RHS>::value;                 \
-			constexpr size_t maxDims = (dimsLhs > dimsRhs) ? dimsLhs : dimsRhs;                    \
+			constexpr uint64_t dimsLhs = VectorScalarDimensionExtractor<LHS>::value;               \
+			constexpr uint64_t dimsRhs = VectorScalarDimensionExtractor<RHS>::value;               \
+			constexpr uint64_t maxDims = (dimsLhs > dimsRhs) ? dimsLhs : dimsRhs;                  \
 			return BinaryVecOp {scalarVectorCaster<Scalar, maxDims>(lhs),                          \
 								scalarVectorCaster<Scalar, maxDims>(rhs),                          \
 								Op {}};                                                            \
@@ -884,7 +902,6 @@ namespace librapid {
                                                                                                    \
 	template<typename Val, typename std::enable_if_t<typetraits::IsVector<Val>::value, int> = 0>   \
 	LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto OP_NAME_(const Val &val) {                      \
-		using namespace ::librapid::vectorDetail;                                                  \
 		using Op = NAME_;                                                                          \
 		return ::librapid::vectorDetail::UnaryVecOp<Val, Op> {val, Op {}};                         \
 	}
@@ -944,7 +961,7 @@ namespace librapid {
 	} // namespace vectorDetail
 
 #define VECTOR_FUNC_IMPL_DEF(NAME_)                                                                \
-	template<typename ScalarType, size_t NumDims>                                                  \
+	template<typename ScalarType, uint64_t NumDims>                                                \
 	auto NAME_(const Vector<ScalarType, NumDims> &vec) {                                           \
 		return vectorDetail::UnaryVecOp {vec, vectorDetail::Vector_##NAME_ {}};                    \
 	}                                                                                              \
@@ -1047,7 +1064,7 @@ struct fmt::formatter<librapid::vectorDetail::VectorBase<Derived>, Char> {
 	}
 };
 
-template<typename Scalar, size_t NumDims, typename Char>
+template<typename Scalar, uint64_t NumDims, typename Char>
 struct fmt::formatter<librapid::Vector<Scalar, NumDims>, Char> {
 	using Base = fmt::formatter<Scalar, Char>;
 	Base m_base;
@@ -1058,7 +1075,7 @@ struct fmt::formatter<librapid::Vector<Scalar, NumDims>, Char> {
 	}
 
 	template<typename FormatContext>
-	auto format(const librapid::Vector<Scalar, NumDims> &vec, FormatContext &ctx) {
+	auto format(const librapid::Vector<Scalar, NumDims> &vec, FormatContext &ctx) const {
 		vec.str(m_base, ctx);
 		return ctx.out();
 	}
diff --git a/librapid/vendor/fmt b/librapid/vendor/fmt
index 130cf54c..d9063baf 160000
--- a/librapid/vendor/fmt
+++ b/librapid/vendor/fmt
@@ -1 +1 @@
-Subproject commit 130cf54cbc77abf0c8200a1206638cb70c9e26f1
+Subproject commit d9063baf227882da0f48c761abcbb08247eb1296