Compiler Explorer

Source code

#include <algorithm>
#include <concepts>
#include <cstdint>
#include <exception>
#include <limits>
#include <numbers>
#include <numeric>
#include <optional>
#include <span>
#include <string_view>
#include <type_traits>
#include <utility>
#include <variant>
#include <vector>

#define BOOST_MP_STANDALONE 1
#include <boost/multiprecision/cpp_int.hpp>

#include <boost/hana/core/make.hpp>

namespace waarudo {

#if WAARUDO_DOCUMENTING and false
// Overview.
// - Number types.
struct exponentiation; // $bⁿ$.
struct rational;       // $p/q$.
struct radical;        // $ⁿ√x$.
union number;          // \expos, sum of `std::intmax_t` + the above, converts from `std::intmax_t` and the above.
union exponent;        // \expos, sum of `std::intmax_t` + `rational`, converts from `std::intmax_t` and `rational`.
struct unit_prefix;
// - Symbol types.
struct transcendental_number;
struct unit;
// - Reference types.
union ref;        // NTTP, converts from `reference`.
struct reference; // Product of `number` × ``range of symbol types'', converts/constructs from ``symbol type''/`number`.
constexpr reference operator*(unit_prefix, unit) noexcept;
constexpr reference operator*(unit_prefix, reference) noexcept;
constexpr reference operator*(number, reference) noexcept;
constexpr reference operator*(reference, reference) noexcept;
constexpr reference operator/(reference, number) noexcept;
constexpr reference operator/(number, reference) noexcept;
constexpr reference operator/(reference, reference) noexcept;
constexpr reference operator^(reference, exponent) noexcept;
// - Symbol to number.
template<const transcendental_number&> struct to_number;            // Not defined, program-specializable.
template<const transcendental_number& T, class N> consteval N to(); // _Returns_: `to_number<T>::template value<N>()`.
// - Variables.
// SI and binary prefixes, π, SI units ("the core of the set", SI Brochure Tables 2 and 4).
#endif

namespace details {

template<auto F> concept constant_invocable = requires { typename std::integral_constant<int, (F(), 0)>; };

constexpr void expects(const bool truth) noexcept {
    if (not truth) std::terminate();
  }

[[nodiscard]] constexpr bool is_negative(const std::totally_ordered_with<std::intmax_t> auto x) noexcept {
    return x < 0;
  }

[[nodiscard]] constexpr bool is_positive(const std::totally_ordered_with<std::intmax_t> auto x) noexcept {
    return x > 0;
  }

[[nodiscard]] constexpr bool is_even(const std::intmax_t x) noexcept { return x % 2 == 0; }

[[nodiscard]] constexpr int sign(const auto& x) noexcept { return x < 0 ? -1 : 1; };

[[nodiscard]] constexpr std::intmax_t abs(const std::intmax_t x) noexcept { return x < 0 ? -x : x; }

template<class T> concept ratio = //
    requires(T r) {
      r.numerator;
      r.denominator;
      T{r.numerator, r.denominator};
    };

constexpr void to_canonical_form(ratio auto& r) noexcept {
    const auto abs_num = abs(r.numerator);
    const auto abs_den = abs(r.denominator);

using std::gcd;
    const auto gcd_ = gcd(abs_num, abs_den);

r.numerator   = sign(r.numerator) * sign(r.denominator) * abs_num / gcd_;
    r.denominator = abs_den / gcd_;
  }

// Workaround LWG3436.
  template<class V> [[nodiscard]] constexpr V construct_at(
    const std::integral auto i,
    std::invocable<std::decay_t<decltype(std::declval<const V&>().a1)>> auto alt_init) noexcept {
    V res = {.a1 = {}};
    switch (i) {
    case 1: alt_init((res = V{.a1 = {}}).a1); break;
    case 2: alt_init((res = V{.a2 = {}}).a2); break;
    case 3: alt_init((res = V{.a3 = {}}).a3); break;
    case 4: alt_init((res = V{.a4 = {}}).a4); break;
    case 5: alt_init((res = V{.a5 = {}}).a5); break;
    case 6: alt_init((res = V{.a6 = {}}).a6); break;
    case 7: alt_init((res = V{.a7 = {}}).a7); break;
    case 8: alt_init((res = V{.a8 = {}}).a8); break;
    case 9: alt_init((res = V{.a9 = {}}).a9); break;
    default: std::unreachable();
    }
    return res;
  }

} // namespace details

struct exponentiation {
  std::intmax_t base;
  std::intmax_t power;
};

struct rational {
  std::intmax_t numerator;
  std::intmax_t denominator;

[[nodiscard]] constexpr explicit(false) rational(const std::intmax_t n, const std::intmax_t d = 1) noexcept
    : numerator{n}, denominator{d} {
    details::expects(d != 0);
    details::to_canonical_form(*this);
  }

[[nodiscard]] bool operator==(const rational&) const = default;

[[nodiscard]] friend constexpr std::strong_ordering operator<=>(const rational l, const rational r) noexcept {
    return l.numerator * r.denominator <=> r.numerator * l.denominator;
  }

[[nodiscard]] constexpr rational operator+() const noexcept { return *this; }
  [[nodiscard]] constexpr rational operator-() const noexcept { return {-numerator, denominator}; }

[[nodiscard]] friend constexpr rational operator+(const rational l, const rational r) noexcept {
    return {l.numerator * r.denominator + r.numerator * l.denominator, l.denominator * r.denominator};
  }
  [[nodiscard]] friend constexpr rational operator-(const rational l, const rational r) noexcept { return l + -r; }
  [[nodiscard]] friend constexpr rational operator*(const rational l, const rational r) noexcept {
    return {l.numerator * r.numerator, l.denominator * r.denominator};
  }
  [[nodiscard]] friend constexpr rational operator/(const rational l, const rational r) noexcept {
    return l * inverse(r);
  }

constexpr rational& operator+=(const rational r) & noexcept { return *this = *this + r; }
  constexpr rational& operator-=(const rational r) & noexcept { return *this = *this - r; }
  constexpr rational& operator*=(const rational r) & noexcept { return *this = *this * r; }
  constexpr rational& operator/=(const rational r) & noexcept { return *this = *this / r; }

[[nodiscard]] friend constexpr rational inverse(const rational r) noexcept { return {r.denominator, r.numerator}; }
};

// $ⁿ√x$, an unresolved nth root.
struct radical {
  std::intmax_t index = 1;
  rational radicand;

[[nodiscard]] constexpr radical(const std::intmax_t i, const rational r) noexcept
    : index{r == 0 ? 1 : i}, radicand{r} {
    details::expects(details::is_positive(i));
  }

[[nodiscard]] constexpr explicit(false) radical(const std::convertible_to<rational> auto r) noexcept : radicand{r} { }

[[nodiscard]] bool operator==(const radical&) const = default;
};

namespace details {

using mp_int = boost::multiprecision::checked_int1024_t;

[[nodiscard]] constexpr mp_int pow(const mp_int& base, const mp_int& power) noexcept {
    return powm(base, abs(power), std::numeric_limits<mp_int>::max());
  }

template<ratio R> [[nodiscard]] constexpr R pow(R base, const mp_int& power) noexcept {
    return {pow(base.numerator, power), pow(base.denominator, power)};
  }

[[nodiscard]] constexpr std::optional<std::intmax_t> try_narrow(const mp_int& x) noexcept {
    if (std::numeric_limits<std::intmax_t>::min() <= x and x <= std::numeric_limits<std::intmax_t>::max())
      return std::intmax_t{x};
    return {};
  }

struct mp_ratio {
    mp_int numerator;
    mp_int denominator = 1;

[[nodiscard]] constexpr explicit(false)
      mp_ratio(const std::convertible_to<mp_int> auto& n, const mp_int& d = 1) noexcept
      : numerator{n}, denominator{d} {
      details::expects(d != 0);
      details::to_canonical_form(*this);
    }

[[nodiscard]] constexpr explicit(false) mp_ratio(const rational r) noexcept
      : numerator{r.numerator}, denominator{r.denominator} { }

[[nodiscard]] constexpr explicit(false) mp_ratio(const exponentiation e) noexcept : numerator{1} {
      (e.power >= 0 ? numerator : denominator) = pow(e.base, e.power);
    }

[[nodiscard]] bool operator==(const mp_ratio& r) const = default;

[[nodiscard]] friend constexpr bool operator<(const mp_ratio& l, const mp_ratio& r) noexcept {
      return l.numerator * r.denominator < r.numerator * l.denominator;
    }
    [[nodiscard]] friend constexpr bool operator>(const mp_ratio& l, const mp_ratio& r) noexcept {
      return l.numerator * r.denominator > r.numerator * l.denominator;
    }
    [[nodiscard]] friend constexpr bool operator<=(const mp_ratio& l, const mp_ratio& r) noexcept {
      return l.numerator * r.denominator <= r.numerator * l.denominator;
    }
    [[nodiscard]] friend constexpr bool operator>=(const mp_ratio& l, const mp_ratio& r) noexcept {
      return l.numerator * r.denominator >= r.numerator * l.denominator;
    }

[[nodiscard]] constexpr mp_ratio operator+() const noexcept { return *this; }
    [[nodiscard]] constexpr mp_ratio operator-() const noexcept { return {-numerator, denominator}; }

[[nodiscard]] friend constexpr mp_ratio operator+(const mp_ratio& l, const mp_ratio& r) noexcept {
      return {l.numerator * r.denominator + r.numerator * l.denominator, l.denominator * r.denominator};
    }
    [[nodiscard]] friend constexpr mp_ratio operator-(const mp_ratio& l, const mp_ratio& r) noexcept { return l + -r; }
    [[nodiscard]] friend constexpr mp_ratio operator*(const mp_ratio& l, const mp_ratio& r) noexcept {
      return {l.numerator * r.numerator, l.denominator * r.denominator};
    }
    [[nodiscard]] friend constexpr mp_ratio operator/(const mp_ratio& l, const mp_ratio& r) noexcept {
      return l * inverse(r);
    }
    [[nodiscard]] friend constexpr mp_ratio operator^(const mp_ratio& l, const mp_int& r) noexcept { return pow(l, r); }

constexpr mp_ratio& operator+=(const mp_ratio& r) & noexcept { return *this = *this + r; }
    constexpr mp_ratio& operator-=(const mp_ratio& r) & noexcept { return *this = *this - r; }
    constexpr mp_ratio& operator*=(const mp_ratio& r) & noexcept { return *this = *this * r; }
    constexpr mp_ratio& operator/=(const mp_ratio& r) & noexcept { return *this = *this / r; }
    constexpr mp_ratio& operator^=(const mp_int& r) & noexcept { return *this = *this ^ r; }

[[nodiscard]] friend constexpr mp_ratio inverse(const mp_ratio& r) noexcept { return {r.denominator, r.numerator}; }
  };

struct mp_radical;

constexpr mp_radical operator^(const mp_ratio&, const mp_ratio&) noexcept;

} // namespace details

union exp { // \expos
  struct arg {
    rational value;

[[nodiscard]] constexpr arg(const std::intmax_t x) noexcept : value{x} { }
    [[nodiscard]] constexpr arg(const rational x) noexcept : value{x} { }

[[nodiscard]] constexpr operator rational() const noexcept { return value; }
    [[nodiscard]] constexpr operator details::mp_ratio() const noexcept { return value; }
  };

std::intmax_t exp_ = 1;
  rational exp_rational;

[[nodiscard]] constexpr exp() noexcept { }

[[nodiscard]] constexpr explicit(false) exp(const std::intmax_t x) noexcept : exp_{x} { }

[[nodiscard]] constexpr explicit(false) exp(const rational x) noexcept : exp_rational{x} {
    if (x.numerator == 0 or x.denominator == 1) *this = x.numerator;
  }
};

namespace details {

// $ⁿ√x$, an unresolved nth root.
  struct mp_radical {
    mp_int index = 1;
    mp_ratio radicand;

[[nodiscard]] constexpr mp_radical(const mp_int& i, const mp_ratio& r) noexcept
      : index{r == 0 ? 1 : i}, radicand{r} {
      details::expects(details::is_positive(i));
    }

[[nodiscard]] constexpr explicit(false) mp_radical(const radical r) noexcept
      : index{r.index}, radicand{r.radicand} { }

[[nodiscard]] constexpr explicit(false) mp_radical(const std::convertible_to<mp_ratio> auto& r) noexcept
      : radicand{r} { }

[[nodiscard]] constexpr bool operator==(const mp_radical& r) const noexcept {
      const auto common_index = lcm(index, r.index);
      return common_radicand(common_index, *this) == common_radicand(common_index, r);
    }

[[nodiscard]] friend constexpr mp_radical operator*(const mp_radical& l, const mp_radical& r) noexcept {
      const auto common_index = lcm(l.index, r.index);
      return {common_index, common_radicand(common_index, l) * common_radicand(common_index, r)};
    }

[[nodiscard]] friend constexpr mp_radical operator^(const mp_radical& l, const mp_ratio& r) noexcept {
      return (not is_negative(r) ? l.radicand : inverse(l.radicand)) ^ (r / l.index);
    }

constexpr mp_radical& operator*=(const mp_radical& r) noexcept { return *this = *this * r; }
    constexpr mp_radical& operator^=(const mp_ratio& r) noexcept { return *this = *this ^ r; }

constexpr void invert() & noexcept { radicand = inverse(radicand); }

[[nodiscard]] friend constexpr mp_ratio common_radicand(const mp_int& common_index, const mp_radical& r) noexcept {
      expects(common_index % r.index == 0);
      return pow(r.radicand, common_index / r.index);
    }
  };

[[nodiscard]] constexpr std::optional<radical> try_narrow(const mp_radical& x) noexcept {
    const auto i = try_narrow(x.index), n = try_narrow(x.radicand.numerator), d = try_narrow(x.radicand.denominator);
    if (i and n and d) return {{*i, {*n, *d}}};
    return {};
  }

[[nodiscard]] constexpr mp_radical operator^(const mp_ratio& l, const mp_ratio& r) noexcept {
    return {r.denominator, l ^ r.numerator};
  }

} // namespace details

struct transcendental_number {
  std::string_view symbol;
  const transcendental_number* name = this;
};

union exp_ {
  std::int_least32_t exp;
};

struct times { // \expos
  std::int_least32_t base;
  exp_ power;
};

struct term { // \expos
  std::intmax_t digit;
  times exponentiation;
};

union sum { // \expos, an arbitrary-precision integer.
  // Its representation is a positional notation of base:
  static constexpr exponentiation base = {10, std::numeric_limits<std::intmax_t>::digits10};
  // [ _Example_: $``yotta'' + 1$ is $1₁₀ × 10²⁴ + 1₁₀ × 10⁰$, which as an expression in an error is
  // `sum{term [2]{term{1, times{10, exp_{24}}}, term{1, times{10, exp_{0}}}}}` for GCC and
  // `sum{.a2 = {{1, {10, {.exp_ = 24}}}, {1, {10, {.exp_ = 0}}}}}` for Clang. -- _end example_ ]
  term a1[1];
  term a2[2];
  term a3[3];
  term a4[4];
  term a5[5];
  term a6[6];
  term a7[7];
  term a8[8];
  term a9[9];

[[nodiscard]] static constexpr sum from(details::mp_int x) noexcept {
    constexpr auto divisor = details::mp_ratio{base}.numerator;
    std::vector<term> terms;
    details::mp_int remainder;
    times exponentiation = {base.base, {0}};
    while (x != 0) {
      while (x % base.base == 0) {
        x /= base.base;
        ++exponentiation.power.exp;
      }
      divide_qr(x, divisor, x, remainder);
      terms.push_back(
        {.digit = std::intmax_t{remainder},
         // Workaround ⁰ being omitted from the mangled value; $10⁰ ≠ 10$:
         .exponentiation = exponentiation.power.exp != 0 ? exponentiation : times{1, {0}}});
      exponentiation.power.exp += base.power;
    }
    return details::construct_at<sum>(terms.size(), [&](term* o) { std::copy(terms.rbegin(), terms.rend(), o); });
  }
};

union num { // \expos
  std::intmax_t i;
  struct exponentiation exponentiation;
  struct rational rational;
  struct radical radical;
  union sum sum;

struct arg {
    details::mp_radical value;

[[nodiscard]] constexpr arg(const std::intmax_t x) noexcept : value{x} { }
    [[nodiscard]] constexpr arg(const struct exponentiation x) noexcept : value{x} { }
    [[nodiscard]] constexpr arg(const struct rational x) noexcept : value{x} { }
    [[nodiscard]] constexpr arg(const struct radical x) noexcept : value{x} { }

[[nodiscard]] constexpr operator const details::mp_radical&() const noexcept { return value; }
  };

using variant_base = std::variant<decltype(i), decltype(exponentiation), decltype(rational), decltype(radical)>;

struct variant : variant_base {
    [[nodiscard]] constexpr explicit(false) variant(const std::intmax_t x) noexcept : variant_base{x} { }

[[nodiscard]] constexpr explicit(false) variant(const struct exponentiation x) noexcept : variant_base{x} {
      if (x.base == 1 or x.power == 0) *this = 1;
    }

[[nodiscard]] constexpr explicit(false) variant(const struct rational x) noexcept : variant_base{x} {
      if (x.numerator == 0 or x.denominator == 1) *this = x.numerator;
    }

[[nodiscard]] constexpr explicit(false) variant(const struct radical x) noexcept : variant_base{x} {
      if (x.index == 1) *this = x.radicand;
    }
  };

[[nodiscard]] static constexpr num from(const variant x) noexcept {
    return visit(
      []<class X>(/*const*/ X x) -> num { // Workaround GCC bug 107222.
        if constexpr (std::is_same_v<X, std::intmax_t>) return {.i = x};
        else if constexpr (std::is_same_v<X, struct exponentiation>) return {.exponentiation = x};
        else if constexpr (std::is_same_v<X, struct rational>) return {.rational = x};
        else if constexpr (std::is_same_v<X, struct radical>) return {.radical = x};
      },
      x);
  }

[[nodiscard]] static constexpr num from(const details::mp_radical x) noexcept {
    if (auto r = try_narrow(x)) return from(variant{*r});

details::expects(x.index == 1 and x.radicand.denominator == 1);
    return {.sum = sum::from(x.radicand.numerator)};
  }
};

struct unit_prefix {
  num::variant factor;
  std::string_view symbol;
};

class reference;

struct unit {
  using prefixes_t = std::span<const unit_prefix* const>;

std::string_view symbol = {};
  prefixes_t prefixes     = {};
  reference (*base)()     = nullptr;
  const unit* name        = this;
};

union ref {
  // \expos nested classes.

struct t {
    num number;
    const struct transcendental_number* transcendental_number;
    exp exponent;
  };

struct u {
    num number;
    const struct unit* unit;
    exp exponent;
  };

struct p {
    const unit_prefix* prefix;
    const struct unit* unit;
    exp exponent;
  };

struct te {
    const transcendental_number* name;
    exp exponent;
  };

struct ue {
    const unit* name;
    exp exponent;
  };

union vt {
    te a1[1];
    te a2[2];
    te a3[3];
    te a4[4];
    te a5[5];
    te a6[6];
    te a7[7];
    te a8[8];
    te a9[9];
  };

union vu {
    ue a1[1];
    ue a2[2];
    ue a3[3];
    ue a4[4];
    ue a5[5];
    ue a6[6];
    ue a7[7];
    ue a8[8];
    ue a9[9];
  };

struct us {
    num number;
    vu units;
  };

struct ss {
    num number;
    vt transcendental_numbers;
    vu units;
  };

struct sp {
    num number;
    const te* transcendental_numbers;
    const ue* units;
  };

// \expos data members.

// Expression in error, GCC vs. Clang. For how omitting 1 might look like, see https://gcc.godbolt.org/z/q6oWzn45K.
  num n;      // `ref{num{1}}`                              `{.n = {.i = 1}}`
  ref::t t;   // `ref{ref::t{num{1}, (& pi), exp{1}}}`      `{.t = {{.i = 1}, &pi, {.exp_ = 1}}}`
  ref::u u;   // `ref{ref::u{num{1}, (& metre), exp{1}}}`   `{.u = {{.i = 1}, &metre, {.exp_ = 1}}}`
  ref::p p;   // `ref{ref::p{(& kilo), (& metre), exp{1}}}` `{.p = {&kilo, &metre, {.exp_ = 1}}}`
  ref::us us; // `ref{ref::us{num{1}, ref::vu{ref::ue [2]{ref::ue{(& metre), exp{1}}, ref::ue{(& second), exp{-2}}}}}}`
              // `{.us = {{.i = 1}, {.a2 = {{&metre, {.exp_ = 1}}, {&second, {.exp_ = -2}}}}}}`
  ref::ss ss; // `ref{ref::ss{num{1}, ref::vt{ref::te [1]{ref::te{(& pi), exp{1}}}},
              //                      ref::vu{ref::ue [2]{ref::ue{(& radian), exp{1}}, ref::ue{(& metre), exp{1}}}}}}`
              // `{.ss = {{.i = 1}, {.a1 = {{&pi, {.exp_ = 1}}}},
              //                    {.a2 = {{&radian, {.exp_ = 1}}, {&metre, {.exp_ = 1}}}}}}`
  ref::sp sp;
};

class reference {
private:
  template<class Symbol, class Se> struct symbol {
    const Symbol* name;
    rational exponent = 1;

[[nodiscard]] constexpr operator Se() const noexcept { return {name, exponent}; }
  };

details::mp_radical number = 1;
  std::vector<symbol<transcendental_number, ref::te>> transcendental_numbers;
  std::vector<symbol<unit, ref::ue>> units;

public:
  [[nodiscard]] constexpr explicit reference(const std::intmax_t x) noexcept : number{x} {
    details::expects(details::is_positive(x));
  }

[[nodiscard]] constexpr explicit reference(const exponentiation x) noexcept : number{x} {
    details::expects(not details::is_negative(x.base) or details::is_even(x.power));
  }

[[nodiscard]] constexpr explicit reference(const rational x) noexcept : number{x} {
    details::expects(details::is_positive(x));
  }

[[nodiscard]] constexpr explicit reference(const radical x) noexcept : number{x} {
    details::expects(details::is_positive(x.radicand));
  }

[[nodiscard]] constexpr explicit reference(const num::arg x) noexcept // \expos
    : reference{*details::try_narrow(x)} { }

[[nodiscard]] constexpr explicit(false) reference(const transcendental_number& x) noexcept
    : transcendental_numbers{{&x}} { }

[[nodiscard]] constexpr explicit(false) reference(const unit& x) noexcept : units{{&x}} { }

private:
  [[nodiscard]] constexpr const unit_prefix* number_to_prefix() const noexcept {
    details::expects(transcendental_numbers.empty() and units.size() == 1);

const auto unit_prefixes = units[0].name->prefixes;
    auto p                   = std::ranges::find(unit_prefixes, number, [](const unit_prefix* const p) {
      return std::visit(boost::hana::make<decltype(number)>, p->factor);
    });
    if (p != unit_prefixes.end()) return *p;
    return nullptr;
  };

template<class V> [[nodiscard]] constexpr V make_multiple_symbols(auto& symbols) const noexcept {
    return details::construct_at<V>(symbols.size(), [&](auto* a) { std::ranges::copy(symbols, a); });
  }

template<class V> static constexpr auto max_extent = sizeof(V) / sizeof(std::declval<V&>().a1[0]);

public:
  [[nodiscard]] constexpr explicit(false) operator ref() const noexcept {
    details::expects(
      ((void)"Dynamic size exceeds `ref`'s static capacity. Use `nttp_fit<[] { return R; }>` to convert `R` to `ref`.",
       transcendental_numbers.size() <= max_extent<ref::vt> and units.size() <= max_extent<ref::vu>));

/*const*/ num n = num::from(number); // Workaround GCC bug 107222.

if (transcendental_numbers.empty() and units.empty()) return {.n = n};

if (transcendental_numbers.size() == 1 and units.empty())
      return {.t = {n, transcendental_numbers[0].name, transcendental_numbers[0].exponent}};

if (transcendental_numbers.empty()) {
      if (units.size() == 1) {
        if (auto p = number_to_prefix()) return {.p = {p, units[0].name, units[0].exponent}};
        return {.u = {n, units[0].name, units[0].exponent}};
      }
      return {.us = {n, make_multiple_symbols<ref::vu>(units)}};
    }

return {.ss = {n, make_multiple_symbols<ref::vt>(transcendental_numbers), make_multiple_symbols<ref::vu>(units)}};
  }

friend constexpr reference operator*(unit_prefix, reference) noexcept;
  friend constexpr reference operator*(reference, reference) noexcept;
  friend constexpr reference operator/(reference, reference) noexcept;
  friend constexpr reference operator^(reference, exp::arg) noexcept;

// With P2641 (`std::is_active_member`), each `from` overload becomes non-templated.
  // This is simulated at https://gcc.godbolt.org/z/q6oWzn45K.

template<num N> [[nodiscard]] static consteval reference from() noexcept { // \expos
    if constexpr (details::constant_invocable<[] { (void)auto{N.i}; }>) return reference{N.i};
    else if constexpr (details::constant_invocable<[] { (void)auto{N.exponentiation}; }>)
      return reference{N.exponentiation};
    else if constexpr (details::constant_invocable<[] { (void)auto{N.rational}; }>) return reference{N.rational};
    else if constexpr (details::constant_invocable<[] { (void)auto{N.radical}; }>) return reference{N.radical};
  }

template<exp E> [[nodiscard]] static consteval rational from() noexcept { // \expos
    if constexpr (details::constant_invocable<[] { (void)auto{E.exp_}; }>) return rational{E.exp_};
    else if constexpr (details::constant_invocable<[] { (void)auto{E.exp_rational}; }>) return E.exp_rational;
  }

template<auto M, auto V> [[nodiscard]] static constexpr reference from() noexcept { // \expos
    constexpr auto& r = V.*M;
    return []<int... I>(std::integer_sequence<int, I...>) {
      return (... * (reference{*r[I].name} ^ from<r[I].exponent>()));
    }(std::make_integer_sequence<int, std::extent_v<std::remove_reference_t<decltype(r)>>>{});
  }

template<auto V> requires requires { V.a1[0]; } [[nodiscard]] static constexpr reference from() noexcept { // \expos
    using Vt = decltype(V);
    if constexpr (details::constant_invocable<[] { (void)auto{V.a1[0]}; }>) return from<&Vt::a1, V>();
    else if constexpr (details::constant_invocable<[] { (void)auto{V.a2[0]}; }>) return from<&Vt::a2, V>();
    else if constexpr (details::constant_invocable<[] { (void)auto{V.a3[0]}; }>) return from<&Vt::a3, V>();
    else if constexpr (details::constant_invocable<[] { (void)auto{V.a4[0]}; }>) return from<&Vt::a4, V>();
    else if constexpr (details::constant_invocable<[] { (void)auto{V.a5[0]}; }>) return from<&Vt::a5, V>();
    else if constexpr (details::constant_invocable<[] { (void)auto{V.a6[0]}; }>) return from<&Vt::a6, V>();
    else if constexpr (details::constant_invocable<[] { (void)auto{V.a7[0]}; }>) return from<&Vt::a7, V>();
    else if constexpr (details::constant_invocable<[] { (void)auto{V.a8[0]}; }>) return from<&Vt::a8, V>();
    else if constexpr (details::constant_invocable<[] { (void)auto{V.a9[0]}; }>) return from<&Vt::a9, V>();
  }

template<ref R> [[nodiscard]] static constexpr reference from() noexcept {
    if constexpr (details::constant_invocable<[] { (void)auto{R.n}; }>) return from<R.n>();
    else if constexpr (details::constant_invocable<[] { (void)auto{R.t}; }>)
      return from<R.t.number>() * (*R.t.transcendental_number ^ from<R.t.exponent>());
    else if constexpr (details::constant_invocable<[] { (void)auto{R.u}; }>)
      return from<R.u.number>() * (*R.u.unit ^ from<R.u.exponent>());
    else if constexpr (details::constant_invocable<[] { (void)auto{R.us}; }>)
      return from<R.us.number>() * from<R.us.units>();
    else if constexpr (details::constant_invocable<[] { (void)auto{R.ss}; }>)
      return from<R.ss.number>() * from<R.ss.transcendental_numbers>() * from<R.ss.units>();
  }
};

[[nodiscard]] constexpr reference operator*(const unit_prefix p, const unit& u) noexcept {
  return visit(boost::hana::make<reference>, p.factor) * u;
}

[[nodiscard]] constexpr reference operator*(const unit_prefix p, reference r) noexcept {
  details::expects(r.number == 1 and r.transcendental_numbers.empty() and r.units.size() == 1);
  return visit(boost::hana::make<reference>, p.factor) * std::move(r);
}

[[nodiscard]] constexpr reference operator*(const num::arg n, reference r) noexcept { return reference{n} * r; }

[[nodiscard]] constexpr reference operator*(reference l, const reference r) noexcept {
  l.number *= r.number;

constexpr auto merge = [symbol_order_proj = [](auto symbol) { return symbol.name->symbol; },
                          constant_erase =
                            [](auto& set, auto it) {
                              std::swap(*it, set.back());
                              set.pop_back();
                            }](auto& set, const auto& r) {
    for (auto symbol : r) {
      auto pos = std::ranges::lower_bound(set, symbol_order_proj(symbol), {}, symbol_order_proj);

if (pos == set.end() or pos->name != symbol.name) set.insert(pos, symbol);
      else if ((pos->exponent += symbol.exponent) == 0) constant_erase(set, pos);
    }
  };

merge(l.transcendental_numbers, r.transcendental_numbers);
  merge(l.units, r.units);

return l;
}

[[nodiscard]] constexpr reference operator/(reference l, reference r) noexcept {
  r.number.invert();
  for (auto& t : r.transcendental_numbers) t.exponent *= -1;
  for (auto& u : r.units) u.exponent *= -1;
  return std::move(l) * std::move(r);
}

[[nodiscard]] constexpr reference operator^(reference r, const exp::arg e) noexcept {
  r.number ^= e;
  for (auto& t : r.transcendental_numbers) t.exponent *= e;
  for (auto& u : r.units) u.exponent *= e;
  return r;
}

[[nodiscard]] constexpr reference operator/(reference r, const num::arg n) noexcept {
  return std::move(r) / reference{n};
}

[[nodiscard]] constexpr reference operator/(const num::arg n, reference r) noexcept {
  return reference{n} / std::move(r);
}

// Prefix variables.

inline constexpr unit_prefix yocto = {exponentiation{10, -24}, "y"};
inline constexpr unit_prefix zepto = {exponentiation{10, -21}, "z"};
inline constexpr unit_prefix atto  = {exponentiation{10, -18}, "a"};
inline constexpr unit_prefix femto = {exponentiation{10, -15}, "f"};
inline constexpr unit_prefix pico  = {exponentiation{10, -12}, "p"};
inline constexpr unit_prefix nano  = {exponentiation{10, -9}, "n"};
inline constexpr unit_prefix micro = {exponentiation{10, -6}, "µ"};
inline constexpr unit_prefix milli = {exponentiation{10, -3}, "m"};
inline constexpr unit_prefix centi = {exponentiation{10, -2}, "c"};
inline constexpr unit_prefix deci  = {exponentiation{10, -1}, "d"};
inline constexpr unit_prefix deca  = {exponentiation{10, 1}, "da"};
inline constexpr unit_prefix hecto = {exponentiation{10, 2}, "h"};
inline constexpr unit_prefix kilo  = {exponentiation{10, 3}, "k"};
inline constexpr unit_prefix mega  = {exponentiation{10, 6}, "M"};
inline constexpr unit_prefix giga  = {exponentiation{10, 9}, "G"};
inline constexpr unit_prefix tera  = {exponentiation{10, 12}, "T"};
inline constexpr unit_prefix peta  = {exponentiation{10, 15}, "P"};
inline constexpr unit_prefix exa   = {exponentiation{10, 18}, "E"};
inline constexpr unit_prefix zetta = {exponentiation{10, 21}, "Z"};
inline constexpr unit_prefix yotta = {exponentiation{10, 24}, "Y"};

inline constexpr unit_prefix kibi = {exponentiation{2, 10}, "Ki"};
inline constexpr unit_prefix mebi = {exponentiation{2, 20}, "Mi"};
inline constexpr unit_prefix gibi = {exponentiation{2, 30}, "Gi"};
inline constexpr unit_prefix tebi = {exponentiation{2, 40}, "Ti"};
inline constexpr unit_prefix pebi = {exponentiation{2, 50}, "Pi"};
inline constexpr unit_prefix exbi = {exponentiation{2, 60}, "Ei"};
inline constexpr unit_prefix zebi = {exponentiation{2, 70}, "Zi"};
inline constexpr unit_prefix yobi = {exponentiation{2, 80}, "Yi"};

inline constexpr const unit_prefix* si_and_binary_prefixes[] = {
  &yocto, &zepto, &atto, &femto, &pico,  &nano,  &micro, &milli, &centi, &deci, &deca, &hecto, &kilo, &mega,
  &giga,  &tera,  &peta, &exa,   &zetta, &yotta, &kibi,  &mebi,  &gibi,  &tebi, &pebi, &exbi,  &zebi, &yobi};

inline constexpr unit::prefixes_t si_prefixes = {si_and_binary_prefixes, 20};

static_assert(si_prefixes.back() == &yotta);

// Transcendental number.

template<const transcendental_number&> struct to_number;

template<const transcendental_number& T, class N> requires requires { to_number<T>::template value<N>(); }
consteval N to() {
  return to_number<T>::template value<N>();
}

inline constexpr transcendental_number pi = {"π"};

template<> struct to_number<pi> {
  template<class N> requires requires { std::numbers::pi_v<N>; } static consteval N value() {
    return std::numbers::pi_v<N>;
  }
};

// Unit variables.
// // Base units.
// // // SI.
inline constexpr unit ampere  = {.symbol = "A", .prefixes = si_prefixes};
inline constexpr unit candela = {.symbol = "cd", .prefixes = si_prefixes};
inline constexpr unit gram    = {.symbol = "g", .prefixes = si_prefixes}; // Not ``kilogram''.
inline constexpr unit kelvin  = {.symbol = "K", .prefixes = si_prefixes};
inline constexpr unit metre   = {.symbol = "m", .prefixes = si_prefixes};
inline constexpr unit mole    = {.symbol = "mol", .prefixes = si_prefixes};
inline constexpr unit radian  = {.symbol = "rad", .prefixes = si_prefixes}; // SI extension.
inline constexpr unit second  = {.symbol = "s", .prefixes = si_prefixes};
// // // Others.
inline constexpr unit byte  = {.symbol = "b", .prefixes = si_and_binary_prefixes};
inline constexpr unit pixel = {.symbol = "px", .prefixes = si_and_binary_prefixes};
// // Derived units.
// // // Unit one.
inline constexpr unit percent = {.symbol = "%", .base = [] { return reference{rational{1, 100}}; }};
// // // SI.
inline constexpr unit steradian = {.symbol = "sr", .prefixes = si_prefixes, .base = [] { return radian ^ 2; }};
inline constexpr unit hertz     = {.symbol = "Hz", .prefixes = si_prefixes, .base = [] { return second ^ -1; }};
inline constexpr unit newton    = {
     .symbol = "N", .prefixes = si_prefixes, .base = [] { return (kilo * gram) * metre * (second ^ -2); }};
inline constexpr unit pascal  = {.symbol = "Pa", .prefixes = si_prefixes, .base = [] { return newton * (metre ^ -2); }};
inline constexpr unit joule   = {.symbol = "J", .prefixes = si_prefixes, .base = [] { return newton * metre; }};
inline constexpr unit watt    = {.symbol = "W", .prefixes = si_prefixes, .base = [] { return joule / second; }};
inline constexpr unit coulomb = {.symbol = "C", .prefixes = si_prefixes, .base = [] { return ampere * second; }};
inline constexpr unit volt    = {.symbol = "V", .prefixes = si_prefixes, .base = [] { return watt / ampere; }};
inline constexpr unit farad   = {.symbol = "F", .prefixes = si_prefixes, .base = [] { return coulomb / volt; }};
inline constexpr unit ohm     = {.symbol = "Ω", .prefixes = si_prefixes, .base = [] { return volt / ampere; }};
inline constexpr unit siemens = {.symbol = "S", .prefixes = si_prefixes, .base = [] { return ampere / volt; }};
inline constexpr unit weber   = {.symbol = "Wb", .prefixes = si_prefixes, .base = [] { return volt * second; }};
inline constexpr unit tesla   = {.symbol = "T", .prefixes = si_prefixes, .base = [] { return weber / (metre ^ 2); }};
inline constexpr unit henry   = {.symbol = "H", .prefixes = si_prefixes, .base = [] { return weber / ampere; }};
inline constexpr unit degree_Celsius = {
  .symbol = "°C", .prefixes = si_prefixes, .base = [] { return reference{kelvin}; }};
inline constexpr unit lumen     = {.symbol = "lm", .prefixes = si_prefixes, .base = [] { return candela * steradian; }};
inline constexpr unit lux       = {.symbol = "lx", .prefixes = si_prefixes, .base = [] { return lumen * (metre ^ 2); }};
inline constexpr unit becquerel = {.symbol = "Bq", .prefixes = si_prefixes, .base = [] { return second ^ -1; }};
inline constexpr unit gray    = {.symbol = "Gy", .prefixes = si_prefixes, .base = [] { return joule / (kilo * gram); }};
inline constexpr unit sievert = {.symbol = "Sv", .prefixes = si_prefixes, .base = [] { return joule / (kilo * gram); }};
inline constexpr unit katal   = {.symbol = "kat", .prefixes = si_prefixes, .base = [] { return mole * (second ^ -1); }};
// // // Time.
inline constexpr unit minute = {.symbol = "min", .base = [] { return 60 * second; }};
inline constexpr unit hour   = {.symbol = "h", .base = [] { return 60 * minute; }};
inline constexpr unit day    = {.symbol = "d", .base = [] { return 24 * hour; }};
inline constexpr unit week   = {.base = [] { return 7 * day; }};
inline constexpr unit year   = {.base = [] { return 146'097 * day / 400; }};
inline constexpr unit month  = {.base = [] { return year / 12; }};

} // namespace waarudo

// From https://mpusz.github.io/units/use_cases/extensions.html.
namespace desk_user {

inline constexpr waarudo::unit_prefix package = {6, "pkg"};
inline constexpr waarudo::unit_prefix lorry   = {6 * 40, "lorry"};

inline constexpr const waarudo::unit_prefix* prefixes[] = {&package, &lorry};

inline constexpr waarudo::unit desk = {
  .symbol = "desk", .prefixes = prefixes, .base = [] { return waarudo::rational(3, 10) * (waarudo::metre ^ 2); }};

} // namespace desk_user

template<auto, auto, waarudo::ref...> struct error { };

template<waarudo::ref V, waarudo::ref... U> consteval void assert_template_argument_equivalent() {
  waarudo::details::expects((... and (&V == &U)));
}

int main() {
  using waarudo::details::constant_invocable;

// Non-positive magnitudes are rejected.
  []<auto... N>() { static_assert((... and not constant_invocable<[] { (void)waarudo::reference{N}; }>)); }
    .template operator()<
      -1,                             //
      0,                              //
      waarudo::exponentiation{-1, 1}, //
      waarudo::rational{-1},          //
      waarudo::rational{0, 1},        //
      waarudo::radical{1, -1},        //
      waarudo::radical{2, -1},        // Numbers without a positive real root are also rejected.
      waarudo::radical{1, 0}>();

{
    static constexpr auto max            = std::numeric_limits<std::intmax_t>::max();
    static constexpr waarudo::rational l = {max - 0, max - 1};
    static constexpr waarudo::rational r = {max - 1, max - 2};
    static_assert(not constant_invocable<[] { void(l * r); }>, "Only `reference` uses arbitrary-precision arithmetic.");
    static_assert(constant_invocable<[] { void(waarudo::reference{l} * waarudo::reference{r}); }>);
  }

assert_template_argument_equivalent<
    waarudo::reference{1},                             //
    waarudo::reference{waarudo::exponentiation{0, 0}}, //
    waarudo::reference{waarudo::exponentiation{1, 9}}, //
    waarudo::reference{waarudo::exponentiation{9, 0}}, //
    waarudo::reference{waarudo::rational{1, 1}},       //
    waarudo::reference{waarudo::radical{1, 1}}>();

assert_template_argument_equivalent<
    waarudo::reference{waarudo::rational{1, 2}}, waarudo::reference{1} / 2, waarudo::reference{2} ^ -1>();
  assert_template_argument_equivalent<
    waarudo::reference{waarudo::rational{1, 8}}, waarudo::reference{1} / 8, waarudo::reference{2} ^ -3>();
  assert_template_argument_equivalent<
    waarudo::reference{8}, waarudo::reference{2} ^ 3, waarudo::reference{waarudo::rational{1, 2}} ^ -3>();

assert_template_argument_equivalent<waarudo::metre ^ 2, waarudo::metre * waarudo::metre>();
  assert_template_argument_equivalent<
    waarudo::reference{1}, waarudo::metre*(waarudo::metre ^ -1), waarudo::metre / waarudo::metre>();
  assert_template_argument_equivalent<
    waarudo::micro*(waarudo::metre ^ 6), waarudo::deci * waarudo::metre ^ 6, waarudo::deci * waarudo::metre ^ 2 ^ 3>();

assert_template_argument_equivalent<waarudo::second / 2, waarudo::rational{1, 2} * waarudo::second>();
  assert_template_argument_equivalent<waarudo::second ^ -1, 1 / waarudo::second>();
  assert_template_argument_equivalent<waarudo::second / 2 ^ -1, 2 / waarudo::second>();

assert_template_argument_equivalent<waarudo::reference{1}, waarudo::reference::from<waarudo::reference{1}>()>();
  assert_template_argument_equivalent<
    waarudo::reference{waarudo::pi}, waarudo::reference::from<waarudo::reference{waarudo::pi}>()>();
  assert_template_argument_equivalent<
    waarudo::reference{waarudo::metre}, waarudo::reference::from<waarudo::reference{waarudo::metre}>()>();
  assert_template_argument_equivalent<
    waarudo::pi * waarudo::radian, waarudo::reference::from<waarudo::pi * waarudo::radian>()>();
  assert_template_argument_equivalent<
    waarudo::metre / waarudo::second, waarudo::reference::from<waarudo::metre / waarudo::second>()>();

[[maybe_unused]] static constexpr error e = error<
    sizeof(waarudo::ref),
    waarudo::sum::from(waarudo::details::mp_ratio{get<waarudo::exponentiation>(waarudo::yotta.factor)}.numerator + 1),
    waarudo::reference{9} ^ 0,                                //
    waarudo::reference{waarudo::exponentiation{2, 3}},        //
    waarudo::reference{waarudo::radical{2, {2, 3}}},          //
    waarudo::reference{waarudo::metre},                       //
    waarudo::kilo * waarudo::metre ^ 2,                       //
    waarudo::kilo * waarudo::metre ^ waarudo::rational{1, 2}, //
    waarudo::radical{2, 31} * waarudo::minute ^ 2,            //
    42 * waarudo::metre / (waarudo::second ^ 2),              //
    1729 * (waarudo::pi ^ 2) * waarudo::radian                //
    >();
#ifndef WAARUDO_NOT_ERROR_TESTING
  +e;
#endif
}