refactor: 💥 magnitudes code cleanup + mag_pi is now mag<pi>

docs/users_guide/framework_basics/faster_than_lightspeed_constants.md

@@ -45,7 +45,7 @@ inline constexpr struct speed_of_light_in_vacuum final :
 }  // namespace si2019
 
 inline constexpr struct magnetic_constant final :
-  named_unit<{u8"μ₀", "u_0"}, mag<4> * mag_pi * mag_power<10, -7> * henry / metre> {} magnetic_constant;
+  named_unit<{u8"μ₀", "u_0"}, mag<4> * mag<pi> * mag_power<10, -7> * henry / metre> {} magnetic_constant;
 
 }  // namespace mp_units::si
 ```

docs/users_guide/framework_basics/systems_of_units.md

@@ -184,11 +184,14 @@ For some units, a magnitude might also be irrational. The best example here is a
 is defined using a floating-point magnitude having a factor of the number π (Pi):
 
 ```cpp
-inline constexpr struct mag_pi final : magnitude<std::numbers::pi_v<long double>> {} mag_pi;
+struct pi : mag_constant {
+  static constexpr auto value = std::numbers::pi_v<long double>;
+};
+inline constexpr pi pi;
 ```
 
 ```cpp
-inline constexpr struct degree final : named_unit<{u8"°", "deg"}, mag_pi / mag<180> * si::radian> {} degree;
+inline constexpr struct degree final : named_unit<{u8"°", "deg"}, mag<pi> / mag<180> * si::radian> {} degree;
 ```
 
 

src/core/CMakeLists.txt

@@ -31,7 +31,8 @@ endfunction()
 # core library definition
 add_mp_units_module(
     core mp-units-core
-    HEADERS include/mp-units/bits/core_gmf.h
+    HEADERS include/mp-units/bits/constexpr_math.h
+            include/mp-units/bits/core_gmf.h
             include/mp-units/bits/get_associated_quantity.h
             include/mp-units/bits/hacks.h
             include/mp-units/bits/math_concepts.h

src/core/include/mp-units/bits/constexpr_math.h

@@ -0,0 +1,224 @@
+// 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.
+
+#pragma once
+
+#ifndef MP_UNITS_IN_MODULE_INTERFACE
+#ifdef MP_UNITS_IMPORT_STD
+import std;
+#else
+#include <concepts>
+#include <cstdint>
+#include <cstdlib>
+#include <limits>
+#include <optional>
+#endif
+#endif
+
+namespace mp_units::detail {
+
+// Raise an arbitrary arithmetic type to a positive integer power at compile time.
+template<typename T>
+[[nodiscard]] consteval T int_power(T base, std::integral auto exp)
+{
+  // As this function should only be called at compile time, the terminations herein function as
+  // "parameter-compatible static_asserts", and should not result in terminations at runtime.
+  if (exp < 0) {
+    std::abort();  // int_power only supports positive integer powers
+  }
+
+  constexpr auto checked_multiply = [](auto a, auto b) {
+    const auto result = a * b;
+    MP_UNITS_DIAGNOSTIC_PUSH
+    MP_UNITS_DIAGNOSTIC_IGNORE_FLOAT_EQUAL
+    if (result / a != b) {
+      std::abort();  // Wraparound detected
+    }
+    MP_UNITS_DIAGNOSTIC_POP
+    return result;
+  };
+
+  constexpr auto checked_square = [checked_multiply](auto a) { return checked_multiply(a, a); };
+
+  if (exp == 0) return T{1};
+  if (exp % 2 == 1) return checked_multiply(base, int_power(base, exp - 1));
+  return checked_square(int_power(base, exp / 2));
+}
+
+template<typename T>
+[[nodiscard]] consteval std::optional<T> checked_int_pow(T base, std::uintmax_t exp)
+{
+  T result = T{1};
+  while (exp > 0u) {
+    if (exp % 2u == 1u) {
+      if (base > std::numeric_limits<T>::max() / result) {
+        return std::nullopt;
+      }
+      result *= base;
+    }
+
+    exp /= 2u;
+
+    if (base > std::numeric_limits<T>::max() / base) {
+      return (exp == 0u) ? std::make_optional(result) : std::nullopt;
+    }
+    base *= base;
+  }
+  return result;
+}
+
+template<typename T>
+[[nodiscard]] consteval std::optional<T> root(T x, std::uintmax_t n)
+{
+  // The "zeroth root" would be mathematically undefined.
+  if (n == 0) {
+    return std::nullopt;
+  }
+
+  // The "first root" is trivial.
+  if (n == 1) {
+    return x;
+  }
+
+  // We only support nontrivial roots of floating point types.
+  if (!std::is_floating_point<T>::value) {
+    return std::nullopt;
+  }
+
+  // Handle negative numbers: only odd roots are allowed.
+  if (x < 0) {
+    if (n % 2 == 0) {
+      return std::nullopt;
+    } else {
+      const auto negative_result = root(-x, n);
+      if (!negative_result.has_value()) {
+        return std::nullopt;
+      }
+      return static_cast<T>(-negative_result.value());
+    }
+  }
+
+  // Handle special cases of zero and one.
+  if (x == 0 || x == 1) {
+    return x;
+  }
+
+  // Handle numbers bewtween 0 and 1.
+  if (x < 1) {
+    const auto inverse_result = root(T{1} / x, n);
+    if (!inverse_result.has_value()) {
+      return std::nullopt;
+    }
+    return static_cast<T>(T{1} / inverse_result.value());
+  }
+
+  //
+  // At this point, error conditions are finished, and we can proceed with the "core" algorithm.
+  //
+
+  // Always use `long double` for intermediate computations.  We don't ever expect people to be
+  // calling this at runtime, so we want maximum accuracy.
+  long double xld = static_cast<long double>(x);
+  long double lo = 1.0;
+  long double hi = xld;
+
+  // Do a binary search to find the closest value such that `checked_int_pow` recovers the input.
+  //
+  // Because we know `n > 1`, and `x > 1`, and x^n is monotonically increasing, we know that
+  // `checked_int_pow(lo, n) < x < checked_int_pow(hi, n)`.  We will preserve this as an
+  // invariant.
+  while (lo < hi) {
+    long double mid = lo + (hi - lo) / 2;
+
+    auto result = checked_int_pow(mid, n);
+
+    if (!result.has_value()) {
+      return std::nullopt;
+    }
+
+    // Early return if we get lucky with an exact answer.
+    if (result.value() == xld) {
+      return static_cast<T>(mid);
+    }
+
+    // Check for stagnation.
+    if (mid == lo || mid == hi) {
+      break;
+    }
+
+    // Preserve the invariant that `checked_int_pow(lo, n) < x < checked_int_pow(hi, n)`.
+    if (result.value() < xld) {
+      lo = mid;
+    } else {
+      hi = mid;
+    }
+  }
+
+  // Pick whichever one gets closer to the target.
+  const auto lo_diff = xld - checked_int_pow(lo, n).value();
+  const auto hi_diff = checked_int_pow(hi, n).value() - xld;
+  return static_cast<T>(lo_diff < hi_diff ? lo : hi);
+}
+
+// A converter for the value member variable of magnitude (below).
+//
+// The input is the desired result, but in a (wider) intermediate type.  The point of this function
+// is to cast to the desired type, but avoid overflow in doing so.
+template<typename To, typename From>
+  requires(!std::is_integral_v<To> || std::is_integral_v<From>)
+[[nodiscard]] consteval To checked_static_cast(From x)
+{
+  // This function should only ever be called at compile time.  The purpose of these terminations is
+  // to produce compiler errors, because we cannot `static_assert` on function arguments.
+  if constexpr (std::is_integral_v<To>) {
+    if (!std::in_range<To>(x)) {
+      std::abort();  // Cannot represent magnitude in this type
+    }
+  }
+
+  return static_cast<To>(x);
+}
+
+// The exponent of `factor` in the prime factorization of `n`.
+[[nodiscard]] consteval 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.
+// NOLINTNEXTLINE(bugprone-easily-swappable-parameters)
+[[nodiscard]] consteval std::intmax_t remove_power(std::intmax_t base, std::intmax_t pow, std::intmax_t n)
+{
+  while (pow-- > 0) {
+    n /= base;
+  }
+  return n;
+}
+
+}  // namespace mp_units::detail

src/core/include/mp-units/bits/sudo_cast.h

@@ -76,8 +76,8 @@ struct conversion_type_traits {
  */
 template<Magnitude auto M, typename T>
 struct conversion_value_traits {
-  static constexpr Magnitude auto num = numerator(M);
-  static constexpr Magnitude auto den = denominator(M);
+  static constexpr Magnitude auto num = _numerator(M);
+  static constexpr Magnitude auto den = _denominator(M);
   static constexpr Magnitude auto irr = M * (den / num);
   static constexpr T num_mult = get_value<T>(num);
   static constexpr T den_mult = get_value<T>(den);
@@ -121,9 +121,9 @@ template<Quantity To, typename FwdFrom, Quantity From = std::remove_cvref_t<FwdF
 
     // scale the number
     if constexpr (is_integral(c_mag))
-      return scale([&](auto value) { return value * get_value<multiplier_type>(numerator(c_mag)); });
+      return scale([&](auto value) { return value * get_value<multiplier_type>(_numerator(c_mag)); });
     else if constexpr (is_integral(pow<-1>(c_mag)))
-      return scale([&](auto value) { return value / get_value<multiplier_type>(denominator(c_mag)); });
+      return scale([&](auto value) { return value / get_value<multiplier_type>(_denominator(c_mag)); });
     else {
       using value_traits = conversion_value_traits<c_mag, multiplier_type>;
       if constexpr (std::is_floating_point_v<multiplier_type>)

src/core/include/mp-units/framework/unit.h

@@ -654,7 +654,7 @@ template<Unit U1, Unit U2>
     else if constexpr (is_integral(canonical_rhs.mag / canonical_lhs.mag))
       return u1;
     else {
-      constexpr auto common_magnitude = detail::common_magnitude(canonical_lhs.mag, canonical_rhs.mag);
+      constexpr auto common_magnitude = _common_magnitude(canonical_lhs.mag, canonical_rhs.mag);
       return scaled_unit<common_magnitude, decltype(canonical_lhs.reference_unit)>{};
     }
   }
@@ -741,7 +741,7 @@ constexpr Out unit_symbol_impl(Out out, const scaled_unit_impl<M, U>& u, const u
     // no ratio/prefix
     return unit_symbol_impl<CharT>(out, u.reference_unit, fmt, negative_power);
   } else {
-    constexpr auto mag_txt = magnitude_text<M>();
+    constexpr auto mag_txt = _magnitude_text(M);
     out = copy<CharT>(mag_txt, fmt.encoding, out);
 
     if constexpr (space_before_unit_symbol<scaled_unit<M, U>::reference_unit>) *out++ = ' ';

src/systems/include/mp-units/systems/angular/units.h

@@ -41,7 +41,7 @@ QUANTITY_SPEC(angle, dim_angle);
 QUANTITY_SPEC(solid_angle, pow<2>(angle));
 
 inline constexpr struct radian final : named_unit<"rad", kind_of<angle>> {} radian;
-inline constexpr struct revolution final : named_unit<"rev", mag<2> * mag_pi * radian> {} revolution;
+inline constexpr struct revolution final : named_unit<"rev", mag<2> * mag<pi> * radian> {} revolution;
 inline constexpr struct degree final : named_unit<symbol_text{u8"°", "deg"}, mag_ratio<1, 360> * revolution> {} degree;
 inline constexpr struct gradian final : named_unit<symbol_text{u8"ᵍ", "grad"}, mag_ratio<1, 400> * revolution> {} gradian;
 inline constexpr struct steradian final : named_unit<"sr", square(radian)> {} steradian;

src/systems/include/mp-units/systems/si/chrono.h

@@ -122,7 +122,7 @@ namespace detail {
 [[nodiscard]] constexpr auto as_ratio(Magnitude auto m)
   requires(is_rational(m))
 {
-  return std::ratio<get_value<std::intmax_t>(numerator(m)), get_value<std::intmax_t>(denominator(m))>{};
+  return std::ratio<get_value<std::intmax_t>(_numerator(m)), get_value<std::intmax_t>(_denominator(m))>{};
 }
 
 }  // namespace detail

src/systems/include/mp-units/systems/si/constants.h

@@ -57,7 +57,7 @@ inline constexpr struct luminous_efficacy final :
 inline constexpr struct standard_gravity final :
   named_unit<symbol_text{u8"g₀", "g_0"}, mag_ratio<980'665, 100'000> * metre / square(second)> {} standard_gravity;
 inline constexpr struct magnetic_constant final :
-  named_unit<symbol_text{u8"μ₀", "u_0"}, mag<4> * mag_pi * mag_power<10, -7> * henry / metre> {} magnetic_constant;
+  named_unit<symbol_text{u8"μ₀", "u_0"}, mag<4> * mag<pi> * mag_power<10, -7> * henry / metre> {} magnetic_constant;
 // clang-format on
 
 }  // namespace mp_units::si

src/systems/include/mp-units/systems/si/units.h

@@ -100,7 +100,7 @@ inline constexpr struct minute final : named_unit<"min", mag<60> * si::second> {
 inline constexpr struct hour final : named_unit<"h", mag<60> * minute> {} hour;
 inline constexpr struct day final : named_unit<"d", mag<24> * hour> {} day;
 inline constexpr struct astronomical_unit final : named_unit<"au", mag<149'597'870'700> * si::metre> {} astronomical_unit;
-inline constexpr struct degree final : named_unit<symbol_text{u8"°", "deg"}, mag_pi / mag<180> * si::radian> {} degree;
+inline constexpr struct degree final : named_unit<symbol_text{u8"°", "deg"}, mag<pi> / mag<180> * si::radian> {} degree;
 inline constexpr struct arcminute final : named_unit<symbol_text{u8"′", "'"}, mag_ratio<1, 60> * degree> {} arcminute;
 inline constexpr struct arcsecond final : named_unit<symbol_text{u8"″", "''"}, mag_ratio<1, 60> * arcminute> {} arcsecond;
 inline constexpr struct are final : named_unit<"a", square(si::deca<si::metre>)> {} are;

test/runtime/truncation_test.cpp

@@ -42,11 +42,11 @@ using namespace mp_units;
 using namespace mp_units::angular;
 using namespace mp_units::angular::unit_symbols;
 
-inline constexpr struct half_revolution final : named_unit<"hrev", mag_pi * radian> {
+inline constexpr struct half_revolution final : named_unit<"hrev", mag<pi> * radian> {
 } half_revolution;
 inline constexpr auto hrev = half_revolution;
 
-// constexpr auto revb6 = mag_ratio<1,3> * mag_pi * rad;
+// constexpr auto revb6 = mag_ratio<1,3> * mag<pi> * rad;
 
 TEST_CASE("value_cast should not truncate for valid inputs", "[value_cast]")
 {

test/static/magnitude_test.cpp

@@ -216,9 +216,9 @@ static_assert(std::is_same_v<decltype(get_base(power_v<mag_2, 5, 8>{})), mag_2_>
 
 //   SECTION("pi to the 1 supplies correct values")
 //   {
-//     check_same_type_and_value(get_value<float>(mag_pi), std::numbers::pi_v<float>);
-//     check_same_type_and_value(get_value<double>(mag_pi), std::numbers::pi_v<double>);
-//     check_same_type_and_value(get_value<long double>(mag_pi), std::numbers::pi_v<long double>);
+//     check_same_type_and_value(get_value<float>(mag<pi>), std::numbers::pi_v<float>);
+//     check_same_type_and_value(get_value<double>(mag<pi>), std::numbers::pi_v<double>);
+//     check_same_type_and_value(get_value<long double>(mag<pi>), std::numbers::pi_v<long double>);
 //   }
 
 //   SECTION("pi to arbitrary power performs computations in most accurate type at compile time")

test/static/unit_test.cpp

@@ -74,7 +74,7 @@ inline constexpr struct degree_Celsius_ final : named_unit<symbol_text{u8"℃",
 
 inline constexpr struct minute_ final : named_unit<"min", mag<60> * second> {} minute;
 inline constexpr struct hour_ final : named_unit<"h", mag<60> * minute> {} hour;
-inline constexpr struct degree_ final : named_unit<symbol_text{u8"°", "deg"}, mag_pi / mag<180> * radian> {} degree;
+inline constexpr struct degree_ final : named_unit<symbol_text{u8"°", "deg"}, mag<pi> / mag<180> * radian> {} degree;
 
 inline constexpr struct yard_ final : named_unit<"yd", mag_ratio<9'144, 10'000> * metre> {} yard;
 inline constexpr struct mile_ final : named_unit<"mi", mag<1760> * yard> {} mile;
@@ -140,7 +140,7 @@ static_assert(get_canonical_unit(radian).mag == mag<1>);
 
 static_assert(is_of_type<degree, degree_>);
 static_assert(is_of_type<get_canonical_unit(degree).reference_unit, one_>);
-static_assert(get_canonical_unit(degree).mag == mag_pi / mag<180>);
+static_assert(get_canonical_unit(degree).mag == mag<pi> / mag<180>);
 static_assert(convertible(radian, degree));
 static_assert(radian != degree);