diff --git a/librapid/include/librapid/array/arrayContainer.hpp b/librapid/include/librapid/array/arrayContainer.hpp
index f49b8efd..cbf6063f 100644
--- a/librapid/include/librapid/array/arrayContainer.hpp
+++ b/librapid/include/librapid/array/arrayContainer.hpp
@@ -50,6 +50,7 @@ namespace librapid {
 			using Scalar							   = typename TypeInfo<StorageType_>::Scalar;
 			using Packet							   = std::false_type;
 			using Backend							   = typename TypeInfo<StorageType_>::Backend;
+			using ShapeType							   = ShapeType_;
 			static constexpr int64_t packetWidth	   = 1;
 			static constexpr bool supportsArithmetic   = TypeInfo<Scalar>::supportsArithmetic;
 			static constexpr bool supportsLogical	   = TypeInfo<Scalar>::supportsLogical;
diff --git a/librapid/include/librapid/array/assignOps.hpp b/librapid/include/librapid/array/assignOps.hpp
index 920bfcad..747f022e 100644
--- a/librapid/include/librapid/array/assignOps.hpp
+++ b/librapid/include/librapid/array/assignOps.hpp
@@ -124,15 +124,15 @@ 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 size_t packetWidth = typetraits::TypeInfo<Scalar>::packetWidth;
 
             constexpr bool allowVectorisation =
               typetraits::TypeInfo<
                 detail::Function<descriptor::Trivial, Functor_, Args...>>::allowVectorisation &&
               Function::argsAreSameType;
 
-            const int64_t size       = function.shape().size();
-            const int64_t vectorSize = size - (size % packetWidth);
+            const size_t size       = function.shape().size();
+            const size_t vectorSize = size - (size % packetWidth);
 
 			LIBRAPID_ASSUME(vectorSize % packetWidth == 0);
 
diff --git a/librapid/include/librapid/array/function.hpp b/librapid/include/librapid/array/function.hpp
index 2bfc3928..bc52202b 100644
--- a/librapid/include/librapid/array/function.hpp
+++ b/librapid/include/librapid/array/function.hpp
@@ -2,281 +2,282 @@
 #define LIBRAPID_ARRAY_FUNCTION_HPP
 
 namespace librapid {
-    namespace typetraits {
-        // Extract allowVectorisation from the input types
-        template<typename First, typename... T>
-        constexpr bool checkAllowVectorisation() {
-            if constexpr (sizeof...(T) == 0) {
-                return TypeInfo<std::decay_t<First>>::allowVectorisation;
-            } else {
-                using T1 = typename TypeInfo<std::decay_t<First>>::Scalar;
-                return TypeInfo<std::decay_t<First>>::allowVectorisation &&
-                       checkAllowVectorisation<T...>() &&
-                       (std::is_same_v<T1, typename TypeInfo<std::decay_t<T>>::Scalar> && ...);
-            }
-        }
-
-        template<typename First, typename... Rest>
-        constexpr auto commonBackend() {
-            using FirstBackend = typename TypeInfo<std::decay_t<First>>::Backend;
-            if constexpr (sizeof...(Rest) == 0) {
-                return FirstBackend {};
-            } else {
-                using RestBackend = decltype(commonBackend<Rest...>());
-                if constexpr (std::is_same_v<FirstBackend, backend::OpenCLIfAvailable> ||
-                              std::is_same_v<RestBackend, backend::OpenCLIfAvailable>) {
-                    return backend::OpenCLIfAvailable {};
-                } else if constexpr (std::is_same_v<FirstBackend, backend::CUDAIfAvailable> ||
-                                     std::is_same_v<RestBackend, backend::CUDAIfAvailable>) {
-                    return backend::CUDAIfAvailable {};
-                } else {
-                    return backend::CPU {};
-                }
-            }
-        }
-
-        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 Backend                              = decltype(commonBackend<Args...>());
-
-            static constexpr bool allowVectorisation = checkAllowVectorisation<Args...>();
-
-            static constexpr bool supportsArithmetic = TypeInfo<Scalar>::supportsArithmetic;
-            static constexpr bool supportsLogical    = TypeInfo<Scalar>::supportsLogical;
-            static constexpr bool supportsBinary     = TypeInfo<Scalar>::supportsBinary;
-        };
-
-        LIBRAPID_DEFINE_AS_TYPE(typename desc COMMA typename Functor_ COMMA typename... Args,
-                                ::librapid::detail::Function<desc COMMA Functor_ COMMA Args...>);
-    } // namespace typetraits
-
-    namespace detail {
-        // Descriptor is defined in "forward.hpp"
-
-        template<
-          typename Packet, typename T,
-          typename std::enable_if_t<
-            typetraits::TypeInfo<T>::type != ::librapid::detail::LibRapidType::Scalar, int> = 0>
-        LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Packet packetExtractor(const T &obj,
-                                                                         size_t index) {
-            static_assert(std::is_same_v<Packet, decltype(obj.packet(index))>,
-                          "Packet types do not match");
-            return obj.packet(index);
-        }
-
-        template<
-          typename Packet, typename T,
-          typename std::enable_if_t<
-            typetraits::TypeInfo<T>::type == ::librapid::detail::LibRapidType::Scalar, int> = 0>
-        LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Packet packetExtractor(const T &obj, size_t) {
-            return Packet(obj);
-        }
-
-        template<typename T, typename std::enable_if_t<typetraits::TypeInfo<T>::type !=
-                                                         ::librapid::detail::LibRapidType::Scalar,
-                                                       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>
-        LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto scalarExtractor(const T &obj, size_t) {
-            return obj;
-        }
-
-        template<typename First, typename... Rest>
-        constexpr auto scalarTypesAreSame() {
-            if constexpr (sizeof...(Rest) == 0) {
-                using Scalar = typename typetraits::TypeInfo<std::decay_t<First>>::Scalar;
-                return Scalar {};
-            } else {
-                using RestType = decltype(scalarTypesAreSame<Rest...>());
-                if constexpr (std::is_same_v<
-                                typename typetraits::TypeInfo<std::decay_t<First>>::Scalar,
-                                RestType>) {
-                    return RestType {};
-                } else {
-                    return std::false_type {};
-                }
-            }
-        }
-
-        template<typename desc, typename Functor_, typename... Args>
-        class Function {
-        public:
-            using Type       = Function<desc, Functor_, Args...>;
-            using Functor    = Functor_;
-            using ShapeType  = Shape<LIBRAPID_MAX_ARRAY_DIMS>;
-            using StrideType = ShapeType;
-            using Scalar     = typename typetraits::TypeInfo<Type>::Scalar;
-            using Backend    = typename typetraits::TypeInfo<Type>::Backend;
-            using Packet     = typename typetraits::TypeInfo<Scalar>::Packet;
-            using Iterator   = detail::ArrayIterator<Function>;
-
-            using Descriptor = desc;
-            static constexpr bool argsAreSameType =
-              !std::is_same_v<decltype(scalarTypesAreSame<Args...>()), std::false_type>;
-
-            Function() = default;
-
-            /// Constructs a Function from a functor and arguments.
-            /// \param functor The functor to use.
-            /// \param args The arguments to use.
-            LIBRAPID_ALWAYS_INLINE explicit Function(Functor &&functor, Args &&...args);
-
-            /// Constructs a Function from another function.
-            /// \param other The Function to copy.
-            LIBRAPID_ALWAYS_INLINE Function(const Function &other) = default;
-
-            /// Construct a Function from a temporary function.
-            /// \param other The Function to move.
-            LIBRAPID_ALWAYS_INLINE Function(Function &&other) noexcept = default;
-
-            /// Assigns a Function to this function.
-            /// \param other The Function to copy.
-            /// \return A reference to this Function.
-            LIBRAPID_ALWAYS_INLINE Function &operator=(const Function &other) = default;
-
-            /// Assigns a temporary Function to this Function.
-            /// \param other The Function to move.
-            /// \return A reference to this Function.
-            LIBRAPID_ALWAYS_INLINE Function &operator=(Function &&other) noexcept = default;
-
-            /// Return the shape of the Function's result
-            /// \return The shape of the Function's result
-            LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto shape() const;
-
-            /// Return the arguments in the Function
-            /// \return The arguments in the Function
-            LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto &args() const;
-
-            /// Return an evaluated Array object
-            /// \return
-            LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto eval() const;
-
-            LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto operator[](int64_t index) const;
-
-            /// Evaluates the function at the given index, returning a Packet result.
-            /// \param index The index to evaluate at.
-            /// \return The result of the function (vectorized).
-            LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Packet packet(size_t index) const;
-
-            /// Evaluates the function at the given index, returning a Scalar result.
-            /// \param index The index to evaluate at.
-            /// \return The result of the function (scalar).
-            LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Scalar scalar(size_t index) const;
-
-            LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Iterator begin() const;
-            LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Iterator end() const;
-
-            template<typename T, typename Char, typename Ctx>
-            void str(const fmt::formatter<T, Char> &format, char bracket, char separator,
-                     Ctx &ctx) const;
-
-        private:
-            /// Implementation detail -- evaluates the function at the given index,
-            /// returning a Packet result.
-            /// \tparam I The index sequence.
-            /// \param index The index to evaluate at.
-            /// \return The result of the function (vectorized).
-            template<size_t... I>
-            LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Packet packetImpl(std::index_sequence<I...>,
-                                                                        size_t index) const;
-
-            /// Implementation detail -- evaluates the function at the given index,
-            /// returning a Scalar result.
-            /// \tparam I The index sequence.
-            /// \param index The index to evaluate at.
-            /// \return The result of the function (scalar).
-            template<size_t... I>
-            LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Scalar scalarImpl(std::index_sequence<I...>,
-                                                                        size_t index) const;
-
-            Functor m_functor;
-            std::tuple<Args...> m_args;
-        };
-
-        template<typename desc, typename Functor, typename... Args>
-        Function<desc, Functor, Args...>::Function(Functor &&functor, Args &&...args) :
-                m_functor(std::forward<Functor>(functor)), m_args(std::forward<Args>(args)...) {}
-
-        template<typename desc, typename Functor, typename... Args>
-        auto Function<desc, Functor, Args...>::shape() const {
-            return typetraits::TypeInfo<Functor>::getShape(m_args);
-        }
-
-        template<typename desc, typename Functor, typename... Args>
-        auto &Function<desc, Functor, Args...>::args() const {
-            return m_args;
-        }
-
-        template<typename desc, typename Functor, typename... Args>
-        auto Function<desc, Functor, Args...>::operator[](int64_t index) const {
-            return array::GeneralArrayView(*this)[index];
-        }
-
-        template<typename desc, typename Functor, typename... Args>
-        auto Function<desc, Functor, Args...>::eval() const {
-            auto res = Array<Scalar, Backend>(shape());
-            res      = *this;
-            return res;
-        }
-
-        template<typename desc, typename Functor, typename... Args>
-        typename Function<desc, Functor, Args...>::Packet
-        Function<desc, Functor, Args...>::packet(size_t index) const {
-            return packetImpl(std::make_index_sequence<sizeof...(Args)>(), index);
-        }
-
-        template<typename desc, typename Functor, typename... Args>
-        template<size_t... I>
-        auto Function<desc, Functor, Args...>::packetImpl(std::index_sequence<I...>,
-                                                          size_t index) const -> Packet {
-            return m_functor.packet(packetExtractor<Packet>(std::get<I>(m_args), index)...);
-        }
-
-        template<typename desc, typename Functor, typename... Args>
-        auto Function<desc, Functor, Args...>::scalar(size_t index) const -> Scalar {
-            return scalarImpl(std::make_index_sequence<sizeof...(Args)>(), index);
-        }
-
-        template<typename desc, typename Functor, typename... Args>
-        template<size_t... I>
-        auto Function<desc, Functor, Args...>::scalarImpl(std::index_sequence<I...>,
-                                                          size_t index) const -> Scalar {
-            return m_functor(scalarExtractor(std::get<I>(m_args), index)...);
-        }
-
-        template<typename desc, typename Functor, typename... Args>
-        auto Function<desc, Functor, Args...>::begin() const -> Iterator {
-            return Iterator(*this, 0);
-        }
-
-        template<typename desc, typename Functor, typename... Args>
-        auto Function<desc, Functor, Args...>::end() const -> Iterator {
-            return Iterator(*this, shape()[0]);
-        }
-
-        template<typename desc, typename Functor, typename... Args>
-        template<typename T, typename Char, typename Ctx>
-        void Function<desc, Functor, Args...>::str(const fmt::formatter<T, Char> &format,
-                                                   char bracket, char separator, Ctx &ctx) const {
-            array::GeneralArrayView(*this).str(format, bracket, separator, ctx);
-        }
-    } // namespace detail
+	namespace typetraits {
+		// Extract allowVectorisation from the input types
+		template<typename First, typename... T>
+		constexpr bool checkAllowVectorisation() {
+			if constexpr (sizeof...(T) == 0) {
+				return TypeInfo<std::decay_t<First>>::allowVectorisation;
+			} else {
+				using T1 = typename TypeInfo<std::decay_t<First>>::Scalar;
+				return TypeInfo<std::decay_t<First>>::allowVectorisation &&
+					   checkAllowVectorisation<T...>() &&
+					   (std::is_same_v<T1, typename TypeInfo<std::decay_t<T>>::Scalar> && ...);
+			}
+		}
+
+		template<typename First, typename... Rest>
+		constexpr auto commonBackend() {
+			using FirstBackend = typename TypeInfo<std::decay_t<First>>::Backend;
+			if constexpr (sizeof...(Rest) == 0) {
+				return FirstBackend {};
+			} else {
+				using RestBackend = decltype(commonBackend<Rest...>());
+				if constexpr (std::is_same_v<FirstBackend, backend::OpenCLIfAvailable> ||
+							  std::is_same_v<RestBackend, backend::OpenCLIfAvailable>) {
+					return backend::OpenCLIfAvailable {};
+				} else if constexpr (std::is_same_v<FirstBackend, backend::CUDAIfAvailable> ||
+									 std::is_same_v<RestBackend, backend::CUDAIfAvailable>) {
+					return backend::CUDAIfAvailable {};
+				} else {
+					return backend::CPU {};
+				}
+			}
+		}
+
+		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 Backend							   = decltype(commonBackend<Args...>());
+
+			static constexpr bool allowVectorisation = checkAllowVectorisation<Args...>();
+
+			static constexpr bool supportsArithmetic = TypeInfo<Scalar>::supportsArithmetic;
+			static constexpr bool supportsLogical	 = TypeInfo<Scalar>::supportsLogical;
+			static constexpr bool supportsBinary	 = TypeInfo<Scalar>::supportsBinary;
+		};
+
+		LIBRAPID_DEFINE_AS_TYPE(typename desc COMMA typename Functor_ COMMA typename... Args,
+								::librapid::detail::Function<desc COMMA Functor_ COMMA Args...>);
+	} // namespace typetraits
+
+	namespace detail {
+		// Descriptor is defined in "forward.hpp"
+
+		template<
+		  typename Packet, typename T,
+		  typename std::enable_if_t<
+			typetraits::TypeInfo<T>::type != ::librapid::detail::LibRapidType::Scalar, int> = 0>
+		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Packet packetExtractor(const T &obj,
+																		 size_t index) {
+			static_assert(std::is_same_v<Packet, decltype(obj.packet(index))>,
+						  "Packet types do not match");
+			return obj.packet(index);
+		}
+
+		template<
+		  typename Packet, typename T,
+		  typename std::enable_if_t<
+			typetraits::TypeInfo<T>::type == ::librapid::detail::LibRapidType::Scalar, int> = 0>
+		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Packet packetExtractor(const T &obj, size_t) {
+			return Packet(obj);
+		}
+
+		template<typename T, typename std::enable_if_t<typetraits::TypeInfo<T>::type !=
+														 ::librapid::detail::LibRapidType::Scalar,
+													   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>
+		LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto scalarExtractor(const T &obj, size_t) {
+			return obj;
+		}
+
+		template<typename First, typename... Rest>
+		constexpr auto scalarTypesAreSame() {
+			if constexpr (sizeof...(Rest) == 0) {
+				using Scalar = typename typetraits::TypeInfo<std::decay_t<First>>::Scalar;
+				return Scalar {};
+			} else {
+				using RestType = decltype(scalarTypesAreSame<Rest...>());
+				if constexpr (std::is_same_v<
+								typename typetraits::TypeInfo<std::decay_t<First>>::Scalar,
+								RestType>) {
+					return RestType {};
+				} else {
+					return std::false_type {};
+				}
+			}
+		}
+
+		template<typename desc, typename Functor_, typename... Args>
+		class Function {
+		public:
+			using Type	  = Function<desc, Functor_, Args...>;
+			using Functor = Functor_;
+			using ShapeType =
+			  detail::ShapeTypeHelper<typename typetraits::TypeInfo<Args>::ShapeType...>::Type;
+			using StrideType = ShapeType;
+			using Scalar	 = typename typetraits::TypeInfo<Type>::Scalar;
+			using Backend	 = typename typetraits::TypeInfo<Type>::Backend;
+			using Packet	 = typename typetraits::TypeInfo<Scalar>::Packet;
+			using Iterator	 = detail::ArrayIterator<Function>;
+
+			using Descriptor = desc;
+			static constexpr bool argsAreSameType =
+			  !std::is_same_v<decltype(scalarTypesAreSame<Args...>()), std::false_type>;
+
+			Function() = default;
+
+			/// Constructs a Function from a functor and arguments.
+			/// \param functor The functor to use.
+			/// \param args The arguments to use.
+			LIBRAPID_ALWAYS_INLINE explicit Function(Functor &&functor, Args &&...args);
+
+			/// Constructs a Function from another function.
+			/// \param other The Function to copy.
+			LIBRAPID_ALWAYS_INLINE Function(const Function &other) = default;
+
+			/// Construct a Function from a temporary function.
+			/// \param other The Function to move.
+			LIBRAPID_ALWAYS_INLINE Function(Function &&other) noexcept = default;
+
+			/// Assigns a Function to this function.
+			/// \param other The Function to copy.
+			/// \return A reference to this Function.
+			LIBRAPID_ALWAYS_INLINE Function &operator=(const Function &other) = default;
+
+			/// Assigns a temporary Function to this Function.
+			/// \param other The Function to move.
+			/// \return A reference to this Function.
+			LIBRAPID_ALWAYS_INLINE Function &operator=(Function &&other) noexcept = default;
+
+			/// Return the shape of the Function's result
+			/// \return The shape of the Function's result
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto shape() const;
+
+			/// Return the arguments in the Function
+			/// \return The arguments in the Function
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto &args() const;
+
+			/// Return an evaluated Array object
+			/// \return
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto eval() const;
+
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE auto operator[](int64_t index) const;
+
+			/// Evaluates the function at the given index, returning a Packet result.
+			/// \param index The index to evaluate at.
+			/// \return The result of the function (vectorized).
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Packet packet(size_t index) const;
+
+			/// Evaluates the function at the given index, returning a Scalar result.
+			/// \param index The index to evaluate at.
+			/// \return The result of the function (scalar).
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Scalar scalar(size_t index) const;
+
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Iterator begin() const;
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Iterator end() const;
+
+			template<typename T, typename Char, typename Ctx>
+			void str(const fmt::formatter<T, Char> &format, char bracket, char separator,
+					 Ctx &ctx) const;
+
+		private:
+			/// Implementation detail -- evaluates the function at the given index,
+			/// returning a Packet result.
+			/// \tparam I The index sequence.
+			/// \param index The index to evaluate at.
+			/// \return The result of the function (vectorized).
+			template<size_t... I>
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Packet packetImpl(std::index_sequence<I...>,
+																		size_t index) const;
+
+			/// Implementation detail -- evaluates the function at the given index,
+			/// returning a Scalar result.
+			/// \tparam I The index sequence.
+			/// \param index The index to evaluate at.
+			/// \return The result of the function (scalar).
+			template<size_t... I>
+			LIBRAPID_NODISCARD LIBRAPID_ALWAYS_INLINE Scalar scalarImpl(std::index_sequence<I...>,
+																		size_t index) const;
+
+			Functor m_functor;
+			std::tuple<Args...> m_args;
+		};
+
+		template<typename desc, typename Functor, typename... Args>
+		Function<desc, Functor, Args...>::Function(Functor &&functor, Args &&...args) :
+				m_functor(std::forward<Functor>(functor)), m_args(std::forward<Args>(args)...) {}
+
+		template<typename desc, typename Functor, typename... Args>
+		auto Function<desc, Functor, Args...>::shape() const {
+			return typetraits::TypeInfo<Functor>::getShape(m_args);
+		}
+
+		template<typename desc, typename Functor, typename... Args>
+		auto &Function<desc, Functor, Args...>::args() const {
+			return m_args;
+		}
+
+		template<typename desc, typename Functor, typename... Args>
+		auto Function<desc, Functor, Args...>::operator[](int64_t index) const {
+			return array::GeneralArrayView(*this)[index];
+		}
+
+		template<typename desc, typename Functor, typename... Args>
+		auto Function<desc, Functor, Args...>::eval() const {
+			auto res = Array<Scalar, Backend>(shape());
+			res		 = *this;
+			return res;
+		}
+
+		template<typename desc, typename Functor, typename... Args>
+		typename Function<desc, Functor, Args...>::Packet
+		Function<desc, Functor, Args...>::packet(size_t index) const {
+			return packetImpl(std::make_index_sequence<sizeof...(Args)>(), index);
+		}
+
+		template<typename desc, typename Functor, typename... Args>
+		template<size_t... I>
+		auto Function<desc, Functor, Args...>::packetImpl(std::index_sequence<I...>,
+														  size_t index) const -> Packet {
+			return m_functor.packet(packetExtractor<Packet>(std::get<I>(m_args), index)...);
+		}
+
+		template<typename desc, typename Functor, typename... Args>
+		auto Function<desc, Functor, Args...>::scalar(size_t index) const -> Scalar {
+			return scalarImpl(std::make_index_sequence<sizeof...(Args)>(), index);
+		}
+
+		template<typename desc, typename Functor, typename... Args>
+		template<size_t... I>
+		auto Function<desc, Functor, Args...>::scalarImpl(std::index_sequence<I...>,
+														  size_t index) const -> Scalar {
+			return m_functor(scalarExtractor(std::get<I>(m_args), index)...);
+		}
+
+		template<typename desc, typename Functor, typename... Args>
+		auto Function<desc, Functor, Args...>::begin() const -> Iterator {
+			return Iterator(*this, 0);
+		}
+
+		template<typename desc, typename Functor, typename... Args>
+		auto Function<desc, Functor, Args...>::end() const -> Iterator {
+			return Iterator(*this, shape()[0]);
+		}
+
+		template<typename desc, typename Functor, typename... Args>
+		template<typename T, typename Char, typename Ctx>
+		void Function<desc, Functor, Args...>::str(const fmt::formatter<T, Char> &format,
+												   char bracket, char separator, Ctx &ctx) const {
+			array::GeneralArrayView(*this).str(format, bracket, separator, ctx);
+		}
+	} // namespace detail
 } // namespace librapid
 
 // Support FMT printing
 #ifdef FMT_API
 ARRAY_TYPE_FMT_IML(typename desc COMMA typename Functor COMMA typename... Args,
-                   librapid::detail::Function<desc COMMA Functor COMMA Args...>)
+				   librapid::detail::Function<desc COMMA Functor COMMA Args...>)
 LIBRAPID_SIMPLE_IO_NORANGE(typename desc COMMA typename Functor COMMA typename... Args,
-                           librapid::detail::Function<desc COMMA Functor COMMA Args...>)
+						   librapid::detail::Function<desc COMMA Functor COMMA Args...>)
 #endif // FMT_API
 
 #endif // LIBRAPID_ARRAY_FUNCTION_HPP
\ No newline at end of file
diff --git a/librapid/include/librapid/array/shape.hpp b/librapid/include/librapid/array/shape.hpp
index 3c0009d0..39ef4f75 100644
--- a/librapid/include/librapid/array/shape.hpp
+++ b/librapid/include/librapid/array/shape.hpp
@@ -10,7 +10,7 @@ namespace librapid {
 	namespace typetraits {
 		LIBRAPID_DEFINE_AS_TYPE(size_t N, Shape<N>);
 		LIBRAPID_DEFINE_AS_TYPE_NO_TEMPLATE(MatrixShape);
-	}
+	} // namespace typetraits
 
 	template<size_t N = LIBRAPID_MAX_ARRAY_DIMS>
 	class Shape {
@@ -492,10 +492,10 @@ namespace librapid {
 	/// \param shapes Remaining (optional) inputs
 	/// \return True if all inputs have the same shape, false otherwise
 	template<typename First, typename Second, typename... Rest,
-			 std::enable_if_t<typetraits::IsSizeType<First>::value &&
-								typetraits::IsSizeType<Second>::value &&
-								(typetraits::IsSizeType<Rest>::value && ...),
-							  int> = 0>
+			 typename std::enable_if_t<typetraits::IsSizeType<First>::value &&
+										 typetraits::IsSizeType<Second>::value &&
+										 (typetraits::IsSizeType<Rest>::value && ...),
+									   int> = 0>
 	LIBRAPID_NODISCARD LIBRAPID_INLINE bool shapesMatch(const First &first, const Second &second,
 														const Rest &...shapes) {
 		if constexpr (sizeof...(Rest) == 0) {
@@ -507,10 +507,10 @@ namespace librapid {
 
 	/// \sa shapesMatch
 	template<typename First, typename Second, typename... Rest,
-			 std::enable_if_t<typetraits::IsSizeType<First>::value &&
-								typetraits::IsSizeType<Second>::value &&
-								(typetraits::IsSizeType<Rest>::value && ...),
-							  int> = 0>
+			 typename std::enable_if_t<typetraits::IsSizeType<First>::value &&
+										 typetraits::IsSizeType<Second>::value &&
+										 (typetraits::IsSizeType<Rest>::value && ...),
+									   int> = 0>
 	LIBRAPID_NODISCARD LIBRAPID_INLINE bool
 	shapesMatch(const std::tuple<First, Second, Rest...> &shapes) {
 		if constexpr (sizeof...(Rest) == 0) {
@@ -521,6 +521,44 @@ namespace librapid {
 					 [](auto, auto, auto... rest) { return std::make_tuple(rest...); }, shapes));
 		}
 	}
+
+	namespace detail {
+		template<typename First, typename Second>
+		struct ShapeTypeHelperImpl {
+			using Type = std::false_type;
+		};
+
+		template<size_t N, size_t M>
+		struct ShapeTypeHelperImpl<Shape<N>, Shape<M>> {
+			using Type = Shape<(N > M ? N : M)>;
+		};
+
+		template<size_t N>
+		struct ShapeTypeHelperImpl<Shape<N>, MatrixShape> {
+			using Type = Shape<N>;
+		};
+
+		template<size_t N>
+		struct ShapeTypeHelperImpl<MatrixShape, Shape<N>> {
+			using Type = Shape<N>;
+		};
+
+		template<>
+		struct ShapeTypeHelperImpl<MatrixShape, MatrixShape> {
+			using Type = MatrixShape;
+		};
+
+		template<typename First, typename Second, typename... Rest>
+		struct ShapeTypeHelper {
+			using FirstResult = typename ShapeTypeHelperImpl<First, Second>::Type;
+			using Type		  = typename ShapeTypeHelper<FirstResult, Rest...>::Type;
+		};
+
+		template<typename First, typename Second>
+		struct ShapeTypeHelper<First, Second> {
+			using Type = typename ShapeTypeHelperImpl<First, Second>::Type;
+		};
+	} // namespace detail
 } // namespace librapid
 
 // Support FMT printing