Add make_ratio() helper and prime factorization

src/core/include/units/magnitude.h

@@ -80,8 +80,10 @@ template <std::intmax_t N>
 using prime_factorization_t = typename prime_factorization<N>::type;
 } // namespace detail
 
-template <std::intmax_t num = 1, std::intmax_t den = 1>
-using ratio_t = quotient_t<detail::prime_factorization_t<num>, detail::prime_factorization_t<den>>;
+template <std::intmax_t N, std::intmax_t D = 1>
+auto make_ratio() {
+  return quotient_t<detail::prime_factorization_t<N>, detail::prime_factorization_t<D>>{};
+}
 
 struct pi {
   static inline constexpr long double value = std::numbers::pi_v<long double>;
@@ -207,4 +209,62 @@ struct product<T, U, Tail...>
   using type = product_t<product_t<T, U>, Tail...>;
 };
 
+namespace detail
+{
+
+////////////////////////////////////////////////////////////////////////////////////////////////////
+// Prime factorization implementation.
+
+// Find the smallest prime factor of `n`.
+constexpr std::intmax_t find_first_factor(std::intmax_t n)
+{
+  for (std::intmax_t f = 2; f * f <= n; f += 1 + (f % 2))
+  {
+    if (n % f == 0) { return f; }
+  }
+  return n;
+}
+
+// The exponent of `factor` in the prime factorization of `n`.
+constexpr std::intmax_t multiplicity(std::intmax_t factor, std::intmax_t n)
+{
+  std::intmax_t m = 0;
+  while (n % factor == 0)
+  {
+    n /= factor;
+    ++m;
+  }
+  return m;
+}
+
+// Divide a number by a given base raised to some power.
+//
+// Undefined unless base > 1, pow >= 0, and (base ^ pow) evenly divides n.
+constexpr std::intmax_t remove_power(std::intmax_t base, std::intmax_t pow, std::intmax_t n)
+{
+  while (pow-- > 0)
+  {
+    n /= base;
+  }
+  return n;
+}
+
+// Specialization for the prime factorization of 1 (base case).
+template <>
+struct prime_factorization<1> { using type = magnitude<>; };
+
+// Specialization for the prime factorization of larger numbers (recursive case).
+template <std::intmax_t N>
+struct prime_factorization {
+  static_assert(N > 1, "Can only factor positive integers.");
+
+  static inline constexpr std::intmax_t first_base = find_first_factor(N);
+  static inline constexpr std::intmax_t first_power = multiplicity(first_base, N);
+  static inline constexpr std::intmax_t remainder = remove_power(first_base, first_power, N);
+
+  using type = product_t<
+    magnitude<int_base_power<first_base, first_power>>, prime_factorization_t<remainder>>;
+};
+
+} // namespace detail
 } // namespace units::mag

test/unit_test/runtime/magnitude_test.cpp

@@ -25,8 +25,8 @@
 #include <catch2/catch.hpp>
 #include <type_traits>
 
-using namespace units;
-using namespace units::mag;
+namespace units::mag
+{
 
 TEST_CASE("Magnitude is invertible")
 {
@@ -211,16 +211,74 @@ TEST_CASE("strictly_increasing")
   }
 }
 
-// TEST_CASE("Ratio shortcut performs prime factorization")
-// {
-//   // CHECK(std::is_same_v<ratio<>, magnitude<>>);
-//   // CHECK(std::is_same_v<ratio<2>, magnitude<int_base_power<2>>>);
-// }
-// 
-// TEST_CASE("Equality works for ratios")
-// {
-//   CHECK(ratio<>{} == ratio<>{});
-//   CHECK(ratio<3>{} == ratio<3>{});
-// 
-//   // CHECK(ratio<3>{} != ratio<4>{});
-// }
+TEST_CASE("make_ratio performs prime factorization correctly")
+{
+  SECTION("Performs prime factorization when denominator is 1")
+  {
+    CHECK(std::is_same_v<decltype(make_ratio<1>()), magnitude<>>);
+    CHECK(std::is_same_v<decltype(make_ratio<2>()), magnitude<int_base_power<2>>>);
+    CHECK(std::is_same_v<decltype(make_ratio<3>()), magnitude<int_base_power<3>>>);
+    CHECK(std::is_same_v<decltype(make_ratio<4>()), magnitude<int_base_power<2, 2>>>);
+
+    CHECK(std::is_same_v<
+        decltype(make_ratio<792>()),
+        magnitude<int_base_power<2, 3>, int_base_power<3, 2>, int_base_power<11>>>);
+  }
+
+  SECTION("Reduces fractions to lowest terms")
+  {
+    CHECK(std::is_same_v<decltype(make_ratio<8, 8>()), magnitude<>>);
+    CHECK(std::is_same_v<
+        decltype(make_ratio<50, 80>()),
+        magnitude<int_base_power<2, -3>, int_base_power<5>>>);
+  }
+}
+
+TEST_CASE("Equality works for magnitudes")
+{
+  SECTION("Equivalent ratios are equal")
+  {
+    CHECK(make_ratio<1>() == make_ratio<1>());
+    CHECK(make_ratio<3>() == make_ratio<3>());
+    CHECK(make_ratio<3, 4>() == make_ratio<9, 12>());
+  }
+
+  SECTION("Different ratios are unequal")
+  {
+    CHECK(make_ratio<3>() != make_ratio<5>());
+    CHECK(make_ratio<3>() != make_ratio<3, 2>());
+  }
+}
+
+namespace detail
+{
+
+TEST_CASE("Prime factorization")
+{
+  SECTION ("1 factors into the null magnitude")
+  {
+    CHECK(std::is_same_v<prime_factorization_t<1>, magnitude<>>);
+  }
+
+  SECTION ("Prime numbers factor into themselves")
+  {
+    CHECK(std::is_same_v<prime_factorization_t<2>, magnitude<int_base_power<2>>>);
+    CHECK(std::is_same_v<prime_factorization_t<3>, magnitude<int_base_power<3>>>);
+    CHECK(std::is_same_v<prime_factorization_t<5>, magnitude<int_base_power<5>>>);
+    CHECK(std::is_same_v<prime_factorization_t<7>, magnitude<int_base_power<7>>>);
+    CHECK(std::is_same_v<prime_factorization_t<11>, magnitude<int_base_power<11>>>);
+
+    CHECK(std::is_same_v<prime_factorization_t<41>, magnitude<int_base_power<41>>>);
+  }
+
+  SECTION("Prime factorization finds factors and multiplicities")
+  {
+    CHECK(std::is_same_v<
+        prime_factorization_t<792>,
+        magnitude<int_base_power<2, 3>, int_base_power<3, 2>, int_base_power<11>>>);
+  }
+}
+
+} // namespace detail
+
+} // namespace units::mag