diff --git a/test/unit_test/runtime/CMakeLists.txt b/test/unit_test/runtime/CMakeLists.txt
index b5cea543..ee86ef23 100644
--- a/test/unit_test/runtime/CMakeLists.txt
+++ b/test/unit_test/runtime/CMakeLists.txt
@@ -23,14 +23,12 @@
 cmake_minimum_required(VERSION 3.2)
 
 find_package(Catch2 3 CONFIG REQUIRED)
+find_package(wg21_linear_algebra CONFIG REQUIRED)
 
-add_executable(
-    unit_tests_runtime
-    distribution_test.cpp fmt_test.cpp
-    # fmt_units_test.cpp
-    math_test.cpp
+add_executable(unit_tests_runtime distribution_test.cpp fmt_test.cpp linear_algebra_test.cpp math_test.cpp)
+target_link_libraries(
+    unit_tests_runtime PRIVATE mp-units::mp-units Catch2::Catch2WithMain wg21_linear_algebra::wg21_linear_algebra
 )
-target_link_libraries(unit_tests_runtime PRIVATE mp-units::mp-units Catch2::Catch2WithMain)
 
 if(CMAKE_CXX_COMPILER_ID STREQUAL "MSVC")
     target_compile_options(
diff --git a/test/unit_test/runtime/fmt_units_test.cpp b/test/unit_test/runtime/fmt_units_test.cpp
deleted file mode 100644
index 03067f3a..00000000
--- a/test/unit_test/runtime/fmt_units_test.cpp
+++ /dev/null
@@ -1,334 +0,0 @@
-// The MIT License (MIT)
-//
-// Copyright (c) 2018 Mateusz Pusz
-//
-// Permission is hereby granted, free of charge, to any person obtaining a copy
-// of this software and associated documentation files (the "Software"), to deal
-// in the Software without restriction, including without limitation the rights
-// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-// copies of the Software, and to permit persons to whom the Software is
-// furnished to do so, subject to the following conditions:
-//
-// The above copyright notice and this permission notice shall be included in all
-// copies or substantial portions of the Software.
-//
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-// SOFTWARE.
-
-#include <catch2/catch_test_macros.hpp>
-#include <units/bits/external/hacks.h>  // IWYU pragma: keep
-#include <units/format.h>
-#include <units/isq/iec80000/iec80000.h>
-#include <units/isq/si/iau/iau.h>
-#include <units/isq/si/imperial/imperial.h>
-#include <units/isq/si/international/international.h>
-#include <units/isq/si/si.h>
-#include <units/isq/si/typographic/typographic.h>
-#include <units/isq/si/uscs/uscs.h>
-#include <units/magnitude.h>
-
-using namespace units::isq::si;
-using namespace units::isq::si::references;
-using namespace units::isq::si::international;
-using namespace units::isq::si::uscs;
-using namespace units::isq::si::iau;
-using namespace units::isq::si::imperial;
-using namespace units::isq::si::imperial::references;
-using namespace units::isq::si::typographic;
-using namespace units::isq::iec80000::references;
-
-TEST_CASE("std::format on synthesized unit symbols", "[text][fmt]")
-{
-  SECTION("time")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_ns) == "1 ns");
-    CHECK(STD_FMT::format("{}", 1_q_us) == "1 µs");
-    CHECK(STD_FMT::format("{}", 1_q_ms) == "1 ms");
-
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_us) == "1 us");
-  }
-
-  SECTION("length")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_mm) == "1 mm");
-    CHECK(STD_FMT::format("{}", 1_q_cm) == "1 cm");
-    CHECK(STD_FMT::format("{}", 1_q_km) == "1 km");
-    CHECK(STD_FMT::format("{}", 1 * ft) == "1 ft");
-    CHECK(STD_FMT::format("{}", 1_q_ft_us) == "1 ft(us)");
-    CHECK(STD_FMT::format("{}", 1_q_yd) == "1 yd");
-    CHECK(STD_FMT::format("{}", 1_q_in) == "1 in");
-    CHECK(STD_FMT::format("{}", 1_q_fathom) == "1 fathom");
-    CHECK(STD_FMT::format("{}", 1_q_fathom_us) == "1 fathom(us)");
-    CHECK(STD_FMT::format("{}", 1_q_mi) == "1 mi");
-    CHECK(STD_FMT::format("{}", 1_q_mi_us) == "1 mi(us)");
-    CHECK(STD_FMT::format("{}", 1_q_naut_mi) == "1 mi(naut)");
-    CHECK(STD_FMT::format("{}", 1_q_ch) == "1 ch");
-    CHECK(STD_FMT::format("{}", 1_q_rd) == "1 rd");
-    CHECK(STD_FMT::format("{}", 1_q_thou) == "1 thou");
-    CHECK(STD_FMT::format("{}", 1_q_pc) == "1 pc");
-    CHECK(STD_FMT::format("{}", 1_q_ly) == "1 ly");
-    CHECK(STD_FMT::format("{}", 1_q_pc) == "1 pc");
-    CHECK(STD_FMT::format("{}", 1_q_angstrom) == "1 angstrom");
-    CHECK(STD_FMT::format("{}", 1_q_au) == "1 au");
-    CHECK(STD_FMT::format("{}", 1_q_pica_comp) == "1 pica(comp)");
-    CHECK(STD_FMT::format("{}", 1_q_pica_prn) == "1 pica(prn)");
-    CHECK(STD_FMT::format("{}", 1_q_point_comp) == "1 point(comp)");
-    CHECK(STD_FMT::format("{}", 1_q_point_prn) == "1 point(prn)");
-  }
-
-  SECTION("mass") { CHECK(STD_FMT::format("{}", 1_q_kg) == "1 kg"); }
-
-  SECTION("area")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_m2) == "1 m²");
-    CHECK(STD_FMT::format("{}", 1_q_mm2) == "1 mm²");
-    CHECK(STD_FMT::format("{}", 1_q_cm2) == "1 cm²");
-    CHECK(STD_FMT::format("{}", 1_q_km2) == "1 km²");
-    CHECK(STD_FMT::format("{}", 1_q_ft2) == "1 ft²");
-
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_m2) == "1 m^2");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_mm2) == "1 mm^2");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_cm2) == "1 cm^2");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_km2) == "1 km^2");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_ft2) == "1 ft^2");
-  }
-
-  SECTION("density")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_kg_per_m3) == "1 kg/m³");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_kg_per_m3) == "1 kg/m^3");
-  }
-
-  SECTION("resistance")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_R) == "1 Ω");
-    CHECK(STD_FMT::format("{}", 1_q_kR) == "1 kΩ");
-    CHECK(STD_FMT::format("{}", 1_q_mR) == "1 mΩ");
-    CHECK(STD_FMT::format("{}", 1_q_MR) == "1 MΩ");
-
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_R) == "1 ohm");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_kR) == "1 kohm");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_mR) == "1 mohm");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_MR) == "1 Mohm");
-  }
-
-  SECTION("voltage")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_V) == "1 V");
-    CHECK(STD_FMT::format("{}", 1_q_mV) == "1 mV");
-    CHECK(STD_FMT::format("{}", 1_q_uV) == "1 µV");
-    CHECK(STD_FMT::format("{}", 1_q_nV) == "1 nV");
-    CHECK(STD_FMT::format("{}", 1_q_pV) == "1 pV");
-
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_uV) == "1 uV");
-  }
-
-  SECTION("volume")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_m3) == "1 m³");
-    CHECK(STD_FMT::format("{}", 1_q_mm3) == "1 mm³");
-    CHECK(STD_FMT::format("{}", 1_q_cm3) == "1 cm³");
-    CHECK(STD_FMT::format("{}", 1_q_km3) == "1 km³");
-    CHECK(STD_FMT::format("{}", 1_q_ft3) == "1 ft³");
-
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_m3) == "1 m^3");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_mm3) == "1 mm^3");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_cm3) == "1 cm^3");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_km3) == "1 km^3");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_ft3) == "1 ft^3");
-  }
-
-  SECTION("frequency")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_mHz) == "1 mHz");
-    CHECK(STD_FMT::format("{}", 1_q_kHz) == "1 kHz");
-    CHECK(STD_FMT::format("{}", 1_q_MHz) == "1 MHz");
-    CHECK(STD_FMT::format("{}", 1_q_GHz) == "1 GHz");
-    CHECK(STD_FMT::format("{}", 1_q_THz) == "1 THz");
-  }
-
-  SECTION("speed")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_m_per_s) == "1 m/s");
-    CHECK(STD_FMT::format("{}", 1_q_km_per_h) == "1 km/h");
-    CHECK(STD_FMT::format("{}", 1_q_mi_per_h) == "1 mi/h");
-  }
-
-  SECTION("acceleration")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_m_per_s2) == "1 m/s²");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_m_per_s2) == "1 m/s^2");
-  }
-
-  SECTION("momentum")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_kg_m_per_s) == "1 kg⋅m/s");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_kg_m_per_s) == "1 kg m/s");
-  }
-
-  SECTION("energy")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_mJ) == "1 mJ");
-    CHECK(STD_FMT::format("{}", 1_q_kJ) == "1 kJ");
-    CHECK(STD_FMT::format("{}", 1_q_MJ) == "1 MJ");
-    CHECK(STD_FMT::format("{}", 1_q_GJ) == "1 GJ");
-  }
-
-  SECTION("power")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_mW) == "1 mW");
-    CHECK(STD_FMT::format("{}", 1_q_kW) == "1 kW");
-    CHECK(STD_FMT::format("{}", 1_q_MW) == "1 MW");
-    CHECK(STD_FMT::format("{}", 1_q_GW) == "1 GW");
-  }
-
-  SECTION("surface tension") { CHECK(STD_FMT::format("{}", 1_q_N_per_m) == "1 N/m"); }
-
-  SECTION("magnetic induction") { CHECK(STD_FMT::format("{}", 1_q_T) == "1 T"); }
-
-  SECTION("magnetic flux")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_Wb) == "1 Wb");
-    CHECK(STD_FMT::format("{}", 1_q_G) == "1 G");
-  }
-
-  SECTION("inductance")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_H) == "1 H");
-    CHECK(STD_FMT::format("{}", 1_q_mH) == "1 mH");
-  }
-
-  SECTION("conductance")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_S) == "1 S");
-    CHECK(STD_FMT::format("{}", 1_q_nS) == "1 nS");
-  }
-
-  SECTION("catalytic activity")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_kat) == "1 kat");
-    CHECK(STD_FMT::format("{}", 1_q_U) == "1 U");
-  }
-
-  SECTION("absorbed dose")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_Gy) == "1 Gy");
-    CHECK(STD_FMT::format("{}", 1_q_kGy) == "1 kGy");
-    CHECK(STD_FMT::format("{}", 1_q_mGy) == "1 mGy");
-  }
-
-  SECTION("addition with common ratio") { CHECK(STD_FMT::format("{}", 1_q_in + 1_q_yd) == "37 in"); }
-
-  SECTION("current density")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_A_per_m2) == "1 A/m²");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_A_per_m2) == "1 A/m^2");
-  }
-
-  SECTION("concentration")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_mol_per_m3) == "1 mol/m³");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_mol_per_m3) == "1 mol/m^3");
-  }
-
-  SECTION("luminance")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_cd_per_m2) == "1 cd/m²");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_cd_per_m2) == "1 cd/m^2");
-  }
-
-  SECTION("dynamic viscosity")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_Pa_s) == "1 Pa⋅s");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_Pa_s) == "1 Pa s");
-  }
-
-  SECTION("heat capacity") { CHECK(STD_FMT::format("{}", 1_q_J_per_K) == "1 J/K"); }
-
-  SECTION("specific heat capacity")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_J_per_kg_K) == "1 J⋅K⁻¹⋅kg⁻¹");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_J_per_kg_K) == "1 J K^-1 kg^-1");
-  }
-
-  SECTION("molar heath capacity")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_J_per_mol_K) == "1 J⋅K⁻¹⋅mol⁻¹");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_J_per_mol_K) == "1 J K^-1 mol^-1");
-  }
-
-  SECTION("thermal conductivity")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_W_per_m_K) == "1 W⋅m⁻¹⋅K⁻¹");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_W_per_m_K) == "1 W m^-1 K^-1");
-  }
-
-  SECTION("electric field strength") { CHECK(STD_FMT::format("{}", 1_q_V_per_m) == "1 V/m"); }
-
-  SECTION("charge density")
-  {
-    CHECK(STD_FMT::format("{}", 1_q_C_per_m3) == "1 C/m³");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_C_per_m3) == "1 C/m^3");
-    CHECK(STD_FMT::format("{}", 1_q_C_per_m2) == "1 C/m²");
-    CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_C_per_m2) == "1 C/m^2");
-  }
-
-  SECTION("permittivity") { CHECK(STD_FMT::format("{}", 1_q_F_per_m) == "1 F/m"); }
-
-  SECTION("permeability") { CHECK(STD_FMT::format("{}", 1_q_H_per_m) == "1 H/m"); }
-
-  SECTION("molar energy") { CHECK(STD_FMT::format("{}", 1_q_J_per_mol) == "1 J/mol"); }
-
-  SECTION("torque") { CHECK(STD_FMT::format("{}", 1_q_N_m_per_rad) == "1 N⋅m/rad"); }
-
-  SECTION("storage_capacity")
-  {
-    CHECK(STD_FMT::format("{}", 1 * bit) == "1 bit");
-    CHECK(STD_FMT::format("{}", 1 * kbit) == "1 kbit");
-    CHECK(STD_FMT::format("{}", 1 * Tibit) == "1 Tibit");
-    CHECK(STD_FMT::format("{}", 1 * B) == "1 B");
-    CHECK(STD_FMT::format("{}", 1 * kB) == "1 kB");
-    CHECK(STD_FMT::format("{}", 1 * TiB) == "1 TiB");
-  }
-
-  SECTION("transfer_rate")
-  {
-    CHECK(STD_FMT::format("{}", 1 * (B / s)) == "1 B/s");
-    CHECK(STD_FMT::format("{}", 1 * (kB / s)) == "1 kB/s");
-    CHECK(STD_FMT::format("{}", 1 * (TB / s)) == "1 TB/s");
-  }
-
-  SECTION("traffic_intesity") { CHECK(STD_FMT::format("{}", 1 * E) == "1 E"); }
-
-  SECTION("modulation_rate")
-  {
-    using namespace units::isq::iec80000;
-    CHECK(STD_FMT::format("{}", 1 * Bd) == "1 Bd");
-    CHECK(STD_FMT::format("{}", 1 * kBd) == "1 kBd");
-    CHECK(STD_FMT::format("{}", 1 * TBd) == "1 TBd");
-    CHECK(STD_FMT::format("{}", quantity_cast<baud>(4 / (2 * s))) == "2 Bd");
-  }
-
-  SECTION("incoherent units with powers")
-  {
-    // TODO(chogg): Reinstate after format/Magnitude redesign.
-    // CHECK(STD_FMT::format("{}", 1_q_mi * 1_q_mi * 1_q_mi) == "1 [15900351812136/3814697265625 × 10⁹] m³");
-    // CHECK(STD_FMT::format("{}", 1_q_au * 1_q_au) == "1 [2237952291797391849 × 10⁴] m²");
-    //
-    // CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_mi * 1_q_mi * 1_q_mi) == "1 [15900351812136/3814697265625 x 10^9] m^3");
-    // CHECK(STD_FMT::format("{:%Q %Aq}", 1_q_au * 1_q_au) == "1 [2237952291797391849 x 10^4] m^2");
-  }
-
-  SECTION("unknown scaled unit with reference different than the dimension's coherent unit")
-  {
-    // TODO(chogg): Reinstate after format/Magnitude redesign.
-    // constexpr auto mag = units::mag<units::ratio{2, 3}>();
-    // CHECK(STD_FMT::format("{}", mass<units::scaled_unit<mag, gram>>(1)) == "1 [2/3 × 10⁻³] kg");
-    // CHECK(STD_FMT::format("{:%Q %Aq}", mass<units::scaled_unit<mag, gram>>(1)) == "1 [2/3 x 10^-3] kg");
-  }
-}
diff --git a/test/unit_test/runtime/linear_algebra_test.cpp b/test/unit_test/runtime/linear_algebra_test.cpp
new file mode 100644
index 00000000..e64ec2af
--- /dev/null
+++ b/test/unit_test/runtime/linear_algebra_test.cpp
@@ -0,0 +1,289 @@
+// The MIT License (MIT)
+//
+// Copyright (c) 2018 Mateusz Pusz
+//
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is
+// furnished to do so, subject to the following conditions:
+//
+// The above copyright notice and this permission notice shall be included in all
+// copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+// SOFTWARE.
+
+#include <linear_algebra.hpp>
+// linear_algebra.hpp has to be included first otherwise the header will fail to compile!
+#include <catch2/catch_test_macros.hpp>
+#include <mp_units/customization_points.h>
+#include <mp_units/format.h>
+#include <mp_units/quantity_io.h>
+#include <mp_units/systems/isq/mechanics.h>
+#include <mp_units/systems/isq/space_and_time.h>
+#include <mp_units/systems/si/si.h>
+#include <cmath>
+
+namespace STD_LA {
+
+template<class ET, class OT>
+std::ostream& operator<<(std::ostream& os, const vector<ET, OT>& v)
+{
+  os << "|";
+  for (auto i = 0U; i < v.size(); ++i) {
+    os << STD_FMT::format(" {:>9}", v(i));
+  }
+  os << " |";
+  return os;
+}
+
+template<class ET, class OT>
+std::ostream& operator<<(std::ostream& os, const matrix<ET, OT>& v)
+{
+  for (auto i = 0U; i < v.rows(); ++i) {
+    os << "|";
+    for (auto j = 0U; j < v.columns(); ++j) {
+      os << STD_FMT::format(" {:>9}", v(i, j));
+    }
+    os << (i != v.rows() - 1U ? " |\n" : " |");
+  }
+  return os;
+}
+
+}  // namespace STD_LA
+
+using namespace mp_units;
+using namespace mp_units::si::unit_symbols;
+
+template<typename Rep = double>
+using vector = std::math::fs_vector<Rep, 3>;
+
+template<typename Rep>
+inline constexpr bool mp_units::is_vector<vector<Rep>> = true;
+
+namespace {
+
+template<typename T>
+[[nodiscard]] auto get_magnitude(const vector<T>& v)
+{
+  return std::hypot(v[0], v[1], v[2]);
+}
+
+template<typename T, typename U>
+[[nodiscard]] vector<decltype(T{} * U{})> cross_product(const vector<T>& a, const vector<U>& b)
+{
+  return {a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]};
+}
+
+template<Quantity Q1, Quantity Q2>
+  requires is_vector<typename Q1::rep> && is_vector<typename Q2::rep> &&
+           requires(typename Q1::rep v1, typename Q2::rep v2) { cross_product(v1, v2); }
+[[nodiscard]] quantity_of<Q1::reference * Q2::reference> auto cross_product(const Q1& q1, const Q2& q2)
+{
+  return (Q1::reference * Q2::reference)(cross_product(q1.number(), q2.number()));
+}
+
+}  // namespace
+
+TEST_CASE("vector quantity", "[la]")
+{
+  SECTION("cast of unit")
+  {
+    SECTION("non-truncating")
+    {
+      const quantity<isq::position_vector[km], vector<int>> v{vector<int>{3, 2, 1}};
+      CHECK(v[m].number() == vector<int>{3000, 2000, 1000});
+    }
+
+    SECTION("truncating")
+    {
+      const quantity<isq::position_vector[m], vector<int>> v{vector<int>{1001, 1002, 1003}};
+      CHECK(quantity_cast<km>(v).number() == vector<int>{1, 1, 1});
+    }
+  }
+
+  SECTION("to scalar magnitude")
+  {
+    const quantity<isq::velocity[km / h], vector<int>> v{vector<int>{2, 3, 6}};
+    const auto speed = get_magnitude(v.number()) * isq::speed[v.unit];  // TODO can we do better here?
+    CHECK(speed.number() == 7);
+  }
+
+  SECTION("multiply by scalar value")
+  {
+    const quantity<isq::position_vector[m], vector<int>> v{vector<int>{1, 2, 3}};
+
+    SECTION("integral")
+    {
+      SECTION("scalar on LHS") { CHECK((2 * v).number() == vector<int>{2, 4, 6}); }
+      SECTION("scalar on RHS") { CHECK((v * 2).number() == vector<int>{2, 4, 6}); }
+    }
+
+    SECTION("floating-point")
+    {
+      SECTION("scalar on LHS") { CHECK((0.5 * v).number() == vector<double>{0.5, 1., 1.5}); }
+      SECTION("scalar on RHS") { CHECK((v * 0.5).number() == vector<double>{0.5, 1., 1.5}); }
+    }
+  }
+
+  SECTION("divide by scalar value")
+  {
+    const quantity<isq::position_vector[m], vector<int>> v{vector<int>{2, 4, 6}};
+
+    SECTION("integral") { CHECK((v / 2).number() == vector<int>{1, 2, 3}); }
+    SECTION("floating-point") { CHECK((v / 0.5).number() == vector<double>{4., 8., 12.}); }
+  }
+
+  SECTION("add")
+  {
+    const quantity<isq::position_vector[m], vector<int>> v{vector<int>{1, 2, 3}};
+
+    SECTION("same unit")
+    {
+      const quantity<isq::position_vector[m], vector<int>> u{vector<int>{3, 2, 1}};
+      CHECK((v + u).number() == vector<int>{4, 4, 4});
+    }
+    SECTION("different units")
+    {
+      const quantity<isq::position_vector[km], vector<int>> u{vector<int>{3, 2, 1}};
+      CHECK((v + u).number() == vector<int>{3001, 2002, 1003});
+    }
+  }
+
+  SECTION("subtract")
+  {
+    const quantity<isq::position_vector[m], vector<int>> v{vector<int>{1, 2, 3}};
+
+    SECTION("same unit")
+    {
+      const quantity<isq::position_vector[m], vector<int>> u{vector<int>{3, 2, 1}};
+      CHECK((v - u).number() == vector<int>{-2, 0, 2});
+    }
+    SECTION("different units")
+    {
+      const quantity<isq::position_vector[km], vector<int>> u{vector<int>{3, 2, 1}};
+      CHECK((v - u).number() == vector<int>{-2999, -1998, -997});
+    }
+  }
+
+  SECTION("multiply by scalar quantity")
+  {
+    const quantity<isq::velocity[m / s], vector<int>> v{vector<int>{1, 2, 3}};
+
+    SECTION("integral")
+    {
+      const auto mass = 2 * isq::mass[kg];
+
+      SECTION("derived_quantity_spec")
+      {
+        SECTION("scalar on LHS") { CHECK((mass * v).number() == vector<int>{2, 4, 6}); }
+        SECTION("scalar on RHS") { CHECK((v * mass).number() == vector<int>{2, 4, 6}); }
+      }
+      SECTION("quantity_cast to momentum")
+      {
+        SECTION("scalar on LHS") { CHECK(quantity_cast<isq::momentum>(mass * v).number() == vector<int>{2, 4, 6}); }
+        SECTION("scalar on RHS") { CHECK(quantity_cast<isq::momentum>(v * mass).number() == vector<int>{2, 4, 6}); }
+      }
+      SECTION("quantity of momentum")
+      {
+        SECTION("scalar on LHS")
+        {
+          const quantity<isq::momentum[kg * m / s], vector<int>> momentum = mass * v;
+          CHECK(momentum.number() == vector<int>{2, 4, 6});
+        }
+        SECTION("scalar on RHS")
+        {
+          const quantity<isq::momentum[kg * m / s], vector<int>> momentum = v * mass;
+          CHECK(momentum.number() == vector<int>{2, 4, 6});
+        }
+      }
+    }
+
+    SECTION("floating-point")
+    {
+      const auto mass = 0.5 * isq::mass[kg];
+
+      SECTION("derived_quantity_spec")
+      {
+        SECTION("scalar on LHS") { CHECK((mass * v).number() == vector<double>{0.5, 1., 1.5}); }
+        SECTION("scalar on RHS") { CHECK((v * mass).number() == vector<double>{0.5, 1., 1.5}); }
+      }
+      SECTION("quantity_cast to momentum")
+      {
+        SECTION("scalar on LHS")
+        {
+          CHECK(quantity_cast<isq::momentum>(mass * v).number() == vector<double>{0.5, 1., 1.5});
+        }
+        SECTION("scalar on RHS")
+        {
+          CHECK(quantity_cast<isq::momentum>(v * mass).number() == vector<double>{0.5, 1., 1.5});
+        }
+      }
+      SECTION("quantity of momentum")
+      {
+        SECTION("scalar on LHS")
+        {
+          const quantity<isq::momentum[kg * m / s], vector<double>> momentum = mass * v;
+          CHECK(momentum.number() == vector<double>{0.5, 1., 1.5});
+        }
+        SECTION("scalar on RHS")
+        {
+          const quantity<isq::momentum[kg * m / s], vector<double>> momentum = v * mass;
+          CHECK(momentum.number() == vector<double>{0.5, 1., 1.5});
+        }
+      }
+    }
+  }
+
+  SECTION("divide by scalar quantity")
+  {
+    const quantity<isq::position_vector[km], vector<int>> pos{vector<int>{30, 20, 10}};
+
+    SECTION("integral")
+    {
+      const auto dur = 2 * isq::duration[h];
+
+      SECTION("derived_quantity_spec") { CHECK((pos / dur).number() == vector<int>{15, 10, 5}); }
+      SECTION("quantity_cast to velocity")
+      {
+        CHECK(quantity_cast<isq::velocity>(pos / dur).number() == vector<int>{15, 10, 5});
+      }
+      SECTION("quantity of velocity")
+      {
+        const quantity<isq::velocity[km / h], vector<int>> v = pos / dur;
+        CHECK(v.number() == vector<int>{15, 10, 5});
+      }
+    }
+
+    SECTION("floating-point")
+    {
+      const auto dur = 0.5 * isq::duration[h];
+
+      SECTION("derived_quantity_spec") { CHECK((pos / dur).number() == vector<double>{60, 40, 20}); }
+      SECTION("quantity_cast to velocity")
+      {
+        CHECK(quantity_cast<isq::velocity>(pos / dur).number() == vector<double>{60, 40, 20});
+      }
+      SECTION("quantity of velocity")
+      {
+        const quantity<isq::velocity[km / h], vector<double>> v = pos / dur;
+        CHECK(v.number() == vector<double>{60, 40, 20});
+      }
+    }
+  }
+
+  SECTION("cross product with a vector quantity")
+  {
+    const quantity<isq::position_vector[m], vector<int>> r{vector<int>{3, 0, 0}};
+    const quantity<isq::force[N], vector<int>> f{vector<int>{0, 10, 0}};
+
+    CHECK(cross_product(r, f) == isq::moment_of_force[N * m](vector<int>{0, 0, 30}));
+  }
+}