diff --git a/docs/users_guide/framework_basics/generic_interfaces.md b/docs/users_guide/framework_basics/generic_interfaces.md
index 3dfbdfc39..df4cc98fc 100644
--- a/docs/users_guide/framework_basics/generic_interfaces.md
+++ b/docs/users_guide/framework_basics/generic_interfaces.md
@@ -65,11 +65,11 @@ some issues start to be clearly visible:
    Trying to use an integral type in this scenario will work only for `s1`, while `s2` and `s3`
    will fail to compile. Failing to compile is a good thing here as the library tries to prevent
    the user from doing a clearly wrong thing. To make the code compile, the user needs to use
-   a dedicated [`value_cast`](value_conversions.md#value-truncating-conversions) like this:
+   dedicated [`value_cast` or `force_in`](value_conversions.md#value-truncating-conversions) like this:
 
     ```cpp
     quantity<isq::speed[mi / h]> s2 = avg_speed(value_cast<km>(140 * mi), 2 * h);
-    quantity<isq::speed[m / s]> s3 = avg_speed(value_cast<km>(20 * m), value_cast<h>(2 * s));
+    quantity<isq::speed[m / s]> s3 = avg_speed((20 * m).force_in(km), (2 * s).force_in(h));
     ```
 
     but the above will obviously provide an incorrect behavior (i.e. division by `0` in the evaluation
diff --git a/docs/users_guide/framework_basics/text_output.md b/docs/users_guide/framework_basics/text_output.md
index d78996680..3f961b4ad 100644
--- a/docs/users_guide/framework_basics/text_output.md
+++ b/docs/users_guide/framework_basics/text_output.md
@@ -73,8 +73,8 @@ associated with this quantity.
     before passing it to the text output:
 
     ```cpp
-    std::cout << v1.in(km / h) << '\n';          // 110 km/h
-    std::cout << value_cast<m / s>(v1) << '\n';  // 30.5556 m/s
+    std::cout << v1.in(km / h) << '\n';       // 110 km/h
+    std::cout << v1.force_in(m / s) << '\n';  // 30.5556 m/s
     ```
 
 
diff --git a/docs/users_guide/framework_basics/value_conversions.md b/docs/users_guide/framework_basics/value_conversions.md
index 5dd8b111b..978f7d294 100644
--- a/docs/users_guide/framework_basics/value_conversions.md
+++ b/docs/users_guide/framework_basics/value_conversions.md
@@ -59,12 +59,18 @@ The second solution is to force a truncating conversion:
 ```cpp
 auto q1 = 5 * m;
 std::cout << value_cast<km>(q1) << '\n';
-quantity<si::kilo<si::metre>, int> q2 = value_cast<km>(q1);
+quantity<si::kilo<si::metre>, int> q2 = q1.force_in(km);
 ```
 
 This explicit cast makes it clear that something unsafe is going on. It is easy to spot in code
 reviews or while chasing a bug in the source code.
 
+!!! note
+
+    `q.force_in(U)` is just a shortcut to run `value_cast<U>(q)`. There is no difference in behavior
+    between those two interfaces. `q.force_in(U)` was added for consistency with `q.in(U)` and
+    `q.force_numerical_value_in(U)`.
+
 Another place where this cast is useful is when a user wants to convert a quantity with
 a floating-point representation to the one using an integral one. Again this is a truncating
 conversion, so an explicit cast is needed:
diff --git a/example/avg_speed.cpp b/example/avg_speed.cpp
index e9ef5e012..c59bacb7f 100644
--- a/example/avg_speed.cpp
+++ b/example/avg_speed.cpp
@@ -56,7 +56,7 @@ constexpr QuantityOf<isq::speed> auto avg_speed(QuantityOf<isq::length> auto d,
 template<QuantityOf<isq::length> D, QuantityOf<isq::time> T, QuantityOf<isq::speed> V>
 void print_result(D distance, T duration, V speed)
 {
-  const auto result_in_kmph = value_cast<si::kilo<si::metre> / non_si::hour>(speed);
+  const auto result_in_kmph = speed.force_in(si::kilo<si::metre> / non_si::hour);
   std::cout << "Average speed of a car that makes " << distance << " in " << duration << " is " << result_in_kmph
             << ".\n";
 }
@@ -101,7 +101,7 @@ void example()
 
     // it is not possible to make a lossless conversion of miles to meters on an integral type
     // (explicit cast needed)
-    print_result(distance, duration, fixed_int_si_avg_speed(value_cast<si::metre>(distance), duration));
+    print_result(distance, duration, fixed_int_si_avg_speed(distance.force_in(m), duration));
     print_result(distance, duration, fixed_double_si_avg_speed(distance, duration));
     print_result(distance, duration, avg_speed(distance, duration));
   }
@@ -119,7 +119,7 @@ void example()
     // also it is not possible to make a lossless conversion of miles to meters on an integral type
     // (explicit cast needed)
     print_result(distance, duration,
-                 fixed_int_si_avg_speed(value_cast<int>(value_cast<si::metre>(distance)), value_cast<int>(duration)));
+                 fixed_int_si_avg_speed(value_cast<int>(distance.force_in(m)), value_cast<int>(duration)));
     print_result(distance, duration, fixed_double_si_avg_speed(distance, duration));
     print_result(distance, duration, avg_speed(distance, duration));
   }
@@ -133,7 +133,7 @@ void example()
 
     // it is not possible to make a lossless conversion of centimeters to meters on an integral type
     // (explicit cast needed)
-    print_result(distance, duration, fixed_int_si_avg_speed(value_cast<si::metre>(distance), duration));
+    print_result(distance, duration, fixed_int_si_avg_speed(distance.force_in(m), duration));
     print_result(distance, duration, fixed_double_si_avg_speed(distance, duration));
     print_result(distance, duration, avg_speed(distance, duration));
   }
@@ -149,7 +149,7 @@ void example()
     // it is not possible to make a lossless conversion of centimeters to meters on an integral type
     // (explicit cast needed)
     print_result(distance, duration,
-                 fixed_int_si_avg_speed(value_cast<int>(value_cast<m>(distance)), value_cast<int>(duration)));
+                 fixed_int_si_avg_speed(value_cast<int>(distance.force_in(m)), value_cast<int>(duration)));
 
     print_result(distance, duration, fixed_double_si_avg_speed(distance, duration));
     print_result(distance, duration, avg_speed(distance, duration));
diff --git a/example/conversion_factor.cpp b/example/conversion_factor.cpp
index a3b810512..5a4eaa97a 100644
--- a/example/conversion_factor.cpp
+++ b/example/conversion_factor.cpp
@@ -33,7 +33,7 @@ template<mp_units::Quantity Target, mp_units::Quantity Source>
   requires std::constructible_from<Target, Source>
 inline constexpr double conversion_factor(Target, Source)
 {
-  return mp_units::value_cast<Target::unit>(1. * Source::reference).numerical_value();
+  return (1. * Source::reference).force_numerical_value_in(Target::unit);
 }
 
 }  // namespace
@@ -56,6 +56,7 @@ int main()
   std::cout << MP_UNITS_STD_FMT::format("conversion factor from lengthA::unit of {:%q} to lengthB::unit of {:%q}:\n\n",
                                         lengthA, lengthB)
             << MP_UNITS_STD_FMT::format("lengthB.value( {} ) == lengthA.value( {} ) * conversion_factor( {} )\n",
-                                        lengthB.numerical_value(), lengthA.numerical_value(),
+                                        lengthB.numerical_value_ref_in(lengthB.unit),
+                                        lengthA.numerical_value_ref_in(lengthA.unit),
                                         conversion_factor(lengthB, lengthA));
 }
diff --git a/example/currency.cpp b/example/currency.cpp
index 587f1f7da..40892ad2a 100644
--- a/example/currency.cpp
+++ b/example/currency.cpp
@@ -81,9 +81,9 @@ quantity<To, Rep> exchange_to(quantity<From, Rep> q)
 template<ReferenceOf<currency> auto To, ReferenceOf<currency> auto From, auto PO, typename Rep>
 quantity_point<To, PO, Rep> exchange_to(quantity_point<From, PO, Rep> q)
 {
-  return quantity_point{
-    zero +
-    static_cast<Rep>(exchange_rate<q.unit, get_unit(To)>() * (q - q.absolute_point_origin).numerical_value()) * To};
+  return quantity_point{zero + static_cast<Rep>(exchange_rate<q.unit, get_unit(To)>() *
+                                                (q - q.absolute_point_origin).numerical_value_in(q.unit)) *
+                                 To};
 }
 
 int main()
@@ -91,6 +91,6 @@ int main()
   quantity_point price_usd = zero + 100 * us_dollar;
   quantity_point price_euro = exchange_to<euro>(price_usd);
 
-  std::cout << price_usd.quantity_from_origin() << " -> " << price_euro.quantity_from_origin() << "\n";
-  // std::cout << price_usd.quantity_from_origin() + price_euro.quantity_from_origin() << "\n";  // does not compile
+  std::cout << price_usd.quantity_from(zero) << " -> " << price_euro.quantity_from(zero) << "\n";
+  // std::cout << price_usd.quantity_from(zero) + price_euro.quantity_from(zero) << "\n";  // does not compile
 }
diff --git a/example/glide_computer.cpp b/example/glide_computer.cpp
index 51486f8b6..69231828b 100644
--- a/example/glide_computer.cpp
+++ b/example/glide_computer.cpp
@@ -83,12 +83,12 @@ void print(const R& gliders)
     std::cout << "- Name: " << g.name << "\n";
     std::cout << "- Polar:\n";
     for (const auto& p : g.polar) {
-      const auto ratio = value_cast<one>(glide_ratio(g.polar[0]));
+      const auto ratio = glide_ratio(g.polar[0]).force_in(one);
       std::cout << MP_UNITS_STD_FMT::format("  * {:%.4Q %q} @ {:%.1Q %q} -> {:%.1Q %q} ({:%.1Q %q})\n", p.climb, p.v,
                                             ratio,
                                             // TODO is it possible to make ADL work below (we need another set of trig
                                             // functions for strong angle in a different namespace)
-                                            value_cast<si::degree>(isq::asin(1 / ratio)));
+                                            isq::asin(1 / ratio).force_in(si::degree));
     }
     std::cout << "\n";
   }
diff --git a/example/include/geographic.h b/example/include/geographic.h
index e24a0aaba..c422fa760 100644
--- a/example/include/geographic.h
+++ b/example/include/geographic.h
@@ -138,7 +138,7 @@ struct MP_UNITS_STD_FMT::formatter<geographic::latitude<T>> :
   auto format(geographic::latitude<T> lat, FormatContext& ctx)
   {
     formatter<typename geographic::latitude<T>::quantity_type>::format(
-      is_gt_zero(lat) ? lat.quantity_from_origin() : -lat.quantity_from_origin(), ctx);
+      is_gt_zero(lat) ? lat.quantity_from(geographic::equator) : -lat.quantity_from(geographic::equator), ctx);
     MP_UNITS_STD_FMT::format_to(ctx.out(), "{}", is_gt_zero(lat) ? " N" : "S");
     return ctx.out();
   }
@@ -151,7 +151,8 @@ struct MP_UNITS_STD_FMT::formatter<geographic::longitude<T>> :
   auto format(geographic::longitude<T> lon, FormatContext& ctx)
   {
     formatter<typename geographic::longitude<T>::quantity_type>::format(
-      is_gt_zero(lon) ? lon.quantity_from_origin() : -lon.quantity_from_origin(), ctx);
+      is_gt_zero(lon) ? lon.quantity_from(geographic::prime_meridian) : -lon.quantity_from(geographic::prime_meridian),
+      ctx);
     MP_UNITS_STD_FMT::format_to(ctx.out(), "{}", is_gt_zero(lon) ? " E" : " W");
     return ctx.out();
   }
@@ -171,28 +172,29 @@ template<typename T>
 distance spherical_distance(position<T> from, position<T> to)
 {
   using namespace mp_units;
-  constexpr auto earth_radius = 6'371 * isq::radius[si::kilo<si::metre>];
+  constexpr quantity earth_radius = 6'371 * isq::radius[si::kilo<si::metre>];
 
   using isq::sin, isq::cos, isq::asin, isq::acos;
 
-  const auto& from_lat = from.lat.quantity_from_origin();
-  const auto& from_lon = from.lon.quantity_from_origin();
-  const auto& to_lat = to.lat.quantity_from_origin();
-  const auto& to_lon = to.lon.quantity_from_origin();
+  const quantity from_lat = from.lat.quantity_from(equator);
+  const quantity from_lon = from.lon.quantity_from(prime_meridian);
+  const quantity to_lat = to.lat.quantity_from(equator);
+  const quantity to_lon = to.lon.quantity_from(prime_meridian);
 
   // https://en.wikipedia.org/wiki/Great-circle_distance#Formulae
   if constexpr (sizeof(T) >= 8) {
     // spherical law of cosines
-    const auto central_angle = acos(sin(from_lat) * sin(to_lat) + cos(from_lat) * cos(to_lat) * cos(to_lon - from_lon));
+    const quantity central_angle =
+      acos(sin(from_lat) * sin(to_lat) + cos(from_lat) * cos(to_lat) * cos(to_lon - from_lon));
     // const auto central_angle = 2 * asin(sqrt(0.5 - cos(to_lat - from_lat) / 2 + cos(from_lat) * cos(to_lat) * (1
     // - cos(lon2_rad - from_lon)) / 2));
 
     return quantity_cast<isq::distance>(earth_radius * central_angle);
   } else {
     // the haversine formula
-    const auto sin_lat = sin((to_lat - from_lat) / 2);
-    const auto sin_lon = sin((to_lon - from_lon) / 2);
-    const auto central_angle = 2 * asin(sqrt(sin_lat * sin_lat + cos(from_lat) * cos(to_lat) * sin_lon * sin_lon));
+    const quantity sin_lat = sin((to_lat - from_lat) / 2);
+    const quantity sin_lon = sin((to_lon - from_lon) / 2);
+    const quantity central_angle = 2 * asin(sqrt(sin_lat * sin_lat + cos(from_lat) * cos(to_lat) * sin_lon * sin_lon));
 
     return quantity_cast<isq::distance>(earth_radius * central_angle);
   }
diff --git a/example/kalman_filter/kalman.h b/example/kalman_filter/kalman.h
index 0d38ebaa7..50621236a 100644
--- a/example/kalman_filter/kalman.h
+++ b/example/kalman_filter/kalman.h
@@ -212,7 +212,7 @@ struct MP_UNITS_STD_FMT::formatter<kalman::estimation<Q>> {
       if constexpr (mp_units::Quantity<Q>)
         return t;
       else
-        return t.quantity_from_origin();
+        return t.quantity_ref_from(t.point_origin);
     }(kalman::get<0>(e.state));
 
     std::string value_buffer;
diff --git a/example/kalman_filter/kalman_filter-example_6.cpp b/example/kalman_filter/kalman_filter-example_6.cpp
index ec2954088..23660370e 100644
--- a/example/kalman_filter/kalman_filter-example_6.cpp
+++ b/example/kalman_filter/kalman_filter-example_6.cpp
@@ -45,7 +45,7 @@ template<QuantityPoint QP, QuantityOf<dimensionless> K>
 void print(auto iteration, K gain, QP measured, kalman::estimation<QP> current, kalman::estimation<QP> next)
 {
   std::cout << MP_UNITS_STD_FMT::format("{:2} | {:7%.4Q} | {:10%.3Q %q} | {:>18.3} | {:>18.3}\n", iteration, gain,
-                                        measured.quantity_from_origin(), current, next);
+                                        measured.quantity_ref_from(QP::point_origin), current, next);
 }
 
 int main()
diff --git a/example/kalman_filter/kalman_filter-example_7.cpp b/example/kalman_filter/kalman_filter-example_7.cpp
index 613a6da7c..d82d73231 100644
--- a/example/kalman_filter/kalman_filter-example_7.cpp
+++ b/example/kalman_filter/kalman_filter-example_7.cpp
@@ -45,7 +45,7 @@ template<QuantityPoint QP, QuantityOf<dimensionless> K>
 void print(auto iteration, K gain, QP measured, kalman::estimation<QP> current, kalman::estimation<QP> next)
 {
   std::cout << MP_UNITS_STD_FMT::format("{:2} | {:7%.4Q} | {:10%.3Q %q} | {:>18.3} | {:>18.3}\n", iteration, gain,
-                                        measured.quantity_from_origin(), current, next);
+                                        measured.quantity_ref_from(QP::point_origin), current, next);
 }
 
 int main()
diff --git a/example/kalman_filter/kalman_filter-example_8.cpp b/example/kalman_filter/kalman_filter-example_8.cpp
index 081b8507b..ca02b7bdd 100644
--- a/example/kalman_filter/kalman_filter-example_8.cpp
+++ b/example/kalman_filter/kalman_filter-example_8.cpp
@@ -45,7 +45,7 @@ template<QuantityPoint QP, QuantityOf<dimensionless> K>
 void print(auto iteration, K gain, QP measured, kalman::estimation<QP> current, kalman::estimation<QP> next)
 {
   std::cout << MP_UNITS_STD_FMT::format("{:2} | {:7%.3Q} | {:10%.3Q %q} | {:>16.2} | {:>16.2}\n", iteration, gain,
-                                        measured.quantity_from_origin(), current, next);
+                                        measured.quantity_ref_from(QP::point_origin), current, next);
 }
 
 int main()
diff --git a/example/total_energy.cpp b/example/total_energy.cpp
index 9d6a9f8b6..fa83a6cd7 100644
--- a/example/total_energy.cpp
+++ b/example/total_energy.cpp
@@ -80,7 +80,7 @@ void si_example()
             << "E = " << E3 << "\n";
 
   std::cout << "\n[converted from SI units back to GeV]\n"
-            << "E = " << value_cast<GeV>(E3) << "\n";
+            << "E = " << E3.force_in(GeV) << "\n";
 }
 
 void natural_example()
diff --git a/example/unmanned_aerial_vehicle.cpp b/example/unmanned_aerial_vehicle.cpp
index 0aca9c832..f880ba567 100644
--- a/example/unmanned_aerial_vehicle.cpp
+++ b/example/unmanned_aerial_vehicle.cpp
@@ -97,8 +97,8 @@ template<earth_gravity_model M>
 hae_altitude<M> to_hae(msl_altitude msl, position<long double> pos)
 {
   const auto geoid_undulation =
-    isq::height(GeographicLibWhatsMyOffset(pos.lat.quantity_from_origin().numerical_value_in(si::degree),
-                                           pos.lon.quantity_from_origin().numerical_value_in(si::degree)) *
+    isq::height(GeographicLibWhatsMyOffset(pos.lat.quantity_from(equator).numerical_value_in(si::degree),
+                                           pos.lon.quantity_from(prime_meridian).numerical_value_in(si::degree)) *
                 si::metre);
   return height_above_ellipsoid<M> + (msl - mean_sea_level - geoid_undulation);
 }
@@ -115,7 +115,7 @@ using hal_altitude = quantity_point<isq::altitude[si::metre], height_above_launc
 template<class CharT, class Traits>
 std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& os, const hal_altitude& a)
 {
-  return os << a.quantity_from_origin() << " HAL";
+  return os << a.quantity_from(height_above_launch) << " HAL";
 }
 
 template<>
@@ -123,7 +123,7 @@ struct MP_UNITS_STD_FMT::formatter<hal_altitude> : formatter<hal_altitude::quant
   template<typename FormatContext>
   auto format(const hal_altitude& a, FormatContext& ctx)
   {
-    formatter<hal_altitude::quantity_type>::format(a.quantity_from_origin(), ctx);
+    formatter<hal_altitude::quantity_type>::format(a.quantity_from(height_above_launch), ctx);
     return MP_UNITS_STD_FMT::format_to(ctx.out(), " HAL");
   }
 };
diff --git a/src/core-fmt/include/mp-units/format.h b/src/core-fmt/include/mp-units/format.h
index 4b0a3ec4c..7929bb2f7 100644
--- a/src/core-fmt/include/mp-units/format.h
+++ b/src/core-fmt/include/mp-units/format.h
@@ -224,9 +224,9 @@ struct quantity_formatter {
   const quantity_format_specs<CharT>& specs;
   Locale loc;
 
-  explicit quantity_formatter(OutputIt o, quantity<Reference, Rep> q, const quantity_format_specs<CharT>& fspecs,
+  explicit quantity_formatter(OutputIt o, const quantity<Reference, Rep>& q, const quantity_format_specs<CharT>& fspecs,
                               Locale lc) :
-      out(o), val(std::move(q).numerical_value()), specs(fspecs), loc(std::move(lc))
+      out(o), val(q.numerical_value_ref_in(q.unit)), specs(fspecs), loc(std::move(lc))
   {
   }
 
@@ -396,7 +396,8 @@ struct MP_UNITS_STD_FMT::formatter<mp_units::quantity<Reference, Rep>, CharT> {
 
     if (begin == end || *begin == '}') {
       // default format should print value followed by the unit separated with 1 space
-      out = mp_units::detail::format_units_quantity_value<CharT>(out, q.numerical_value(), specs.rep, ctx.locale());
+      out = mp_units::detail::format_units_quantity_value<CharT>(out, q.numerical_value_ref_in(q.unit), specs.rep,
+                                                                 ctx.locale());
       if constexpr (mp_units::detail::has_unit_symbol(get_unit(Reference))) {
         if constexpr (mp_units::space_before_unit_symbol<get_unit(Reference)>) *out++ = CharT(' ');
         out = unit_symbol_to<CharT>(out, get_unit(Reference));
diff --git a/src/core-io/include/mp-units/ostream.h b/src/core-io/include/mp-units/ostream.h
index 5fc8489c1..c648a8100 100644
--- a/src/core-io/include/mp-units/ostream.h
+++ b/src/core-io/include/mp-units/ostream.h
@@ -36,9 +36,9 @@ void to_stream(std::basic_ostream<CharT, Traits>& os, const quantity<R, Rep>& q)
 {
   if constexpr (is_same_v<Rep, std::uint8_t> || is_same_v<Rep, std::int8_t>)
     // promote the value to int
-    os << +q.numerical_value();
+    os << +q.numerical_value_ref_in(q.unit);
   else
-    os << q.numerical_value();
+    os << q.numerical_value_ref_in(q.unit);
   if constexpr (has_unit_symbol(get_unit(R))) {
     if constexpr (space_before_unit_symbol<get_unit(R)>) os << " ";
     unit_symbol_to<CharT>(std::ostream_iterator<CharT>(os), get_unit(R));
@@ -49,7 +49,7 @@ void to_stream(std::basic_ostream<CharT, Traits>& os, const quantity<R, Rep>& q)
 
 template<typename CharT, typename Traits, auto R, typename Rep>
 std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& os, const quantity<R, Rep>& q)
-  requires requires { os << q.numerical_value(); }
+  requires requires { os << q.numerical_value_ref_in(q.unit); }
 {
   if (os.width()) {
     // std::setw() applies to the whole quantity output so it has to be first put into std::string
diff --git a/src/core/include/mp-units/bits/quantity_cast.h b/src/core/include/mp-units/bits/quantity_cast.h
index fea62a498..4fa7d08b9 100644
--- a/src/core/include/mp-units/bits/quantity_cast.h
+++ b/src/core/include/mp-units/bits/quantity_cast.h
@@ -51,9 +51,9 @@ template<QuantitySpec auto ToQS, typename Q>
 {
   if constexpr (detail::QuantityKindSpec<std::remove_const_t<decltype(ToQS)>> &&
                 AssociatedUnit<std::remove_const_t<decltype(Q::unit)>>)
-    return make_quantity<Q::unit>(std::forward<Q>(q).numerical_value());
+    return make_quantity<Q::unit>(std::forward<Q>(q).numerical_value_);
   else
-    return make_quantity<reference<ToQS, Q::unit>{}>(std::forward<Q>(q).numerical_value());
+    return make_quantity<reference<ToQS, Q::unit>{}>(std::forward<Q>(q).numerical_value_);
 }
 
 }  // namespace mp_units
diff --git a/src/core/include/mp-units/bits/sudo_cast.h b/src/core/include/mp-units/bits/sudo_cast.h
index e8a2c1f02..49f6dc52e 100644
--- a/src/core/include/mp-units/bits/sudo_cast.h
+++ b/src/core/include/mp-units/bits/sudo_cast.h
@@ -63,10 +63,10 @@ template<Quantity To, typename From>
   if constexpr (q_unit == To::unit) {
     // no scaling of the number needed
     return make_quantity<To::reference>(static_cast<MP_UNITS_TYPENAME To::rep>(
-      std::forward<From>(q).numerical_value()));  // this is the only (and recommended) way to do
-                                                  // a truncating conversion on a number, so we are
-                                                  // using static_cast to suppress all the compiler
-                                                  // warnings on conversions
+      std::forward<From>(q).numerical_value_));  // this is the only (and recommended) way to do
+                                                 // a truncating conversion on a number, so we are
+                                                 // using static_cast to suppress all the compiler
+                                                 // warnings on conversions
   } else {
     // scale the number
     constexpr Magnitude auto c_mag = get_canonical_unit(q_unit).mag / get_canonical_unit(To::unit).mag;
@@ -78,7 +78,7 @@ template<Quantity To, typename From>
     using multiplier_type =
       conditional<treat_as_floating_point<c_rep_type>, std::common_type_t<c_mag_type, long double>, c_mag_type>;
     constexpr auto val = [](Magnitude auto m) { return get_value<multiplier_type>(m); };
-    return static_cast<MP_UNITS_TYPENAME To::rep>(static_cast<c_rep_type>(std::forward<From>(q).numerical_value()) *
+    return static_cast<MP_UNITS_TYPENAME To::rep>(static_cast<c_rep_type>(std::forward<From>(q).numerical_value_) *
                                                   val(num) / val(den) * val(irr)) *
            To::reference;
   }
diff --git a/src/core/include/mp-units/quantity.h b/src/core/include/mp-units/quantity.h
index 674f44c77..c8d1d3196 100644
--- a/src/core/include/mp-units/quantity.h
+++ b/src/core/include/mp-units/quantity.h
@@ -85,7 +85,7 @@ using common_quantity_for = quantity<common_reference(Q1::reference, Q2::referen
 template<Reference auto R, RepresentationOf<get_quantity_spec(R).character> Rep = double>
 class quantity {
 public:
-  Rep value_;  // needs to be public for a structural type
+  Rep numerical_value_;  // needs to be public for a structural type
 
   // member types and values
   static constexpr Reference auto reference = R;
@@ -126,7 +126,7 @@ class quantity {
 
   template<detail::QuantityConvertibleTo<quantity> Q>
   constexpr explicit(!std::convertible_to<typename Q::rep, Rep>) quantity(const Q& q) :
-      value_(detail::sudo_cast<quantity>(q).numerical_value())
+      numerical_value_(detail::sudo_cast<quantity>(q).numerical_value_)
   {
   }
 
@@ -141,32 +141,52 @@ class quantity {
   quantity& operator=(const quantity&) = default;
   quantity& operator=(quantity&&) = default;
 
+  // conversions
+  template<Unit U>
+    requires detail::QuantityConvertibleTo<quantity, quantity<quantity_spec[U{}], Rep>>
+  [[nodiscard]] constexpr quantity<quantity_spec[U{}], Rep> in(U) const
+  {
+    return quantity<quantity_spec[U{}], Rep>{*this};
+  }
+
+  template<Unit U>
+    requires requires(quantity q) { value_cast<U{}>(q); }
+  [[nodiscard]] constexpr quantity<quantity_spec[U{}], Rep> force_in(U) const
+  {
+    return value_cast<U{}>(*this);
+  }
+
   // data access
-#ifdef __cpp_explicit_this_parameter
-  template<typename Self>
-  [[nodiscard]] constexpr auto&& value(this Self&& self) noexcept
+  template<Unit U>
+    requires(U{} == unit)
+  [[nodiscard]] constexpr rep& numerical_value_ref_in(U) & noexcept
   {
-    return std::forward<Self>(self).value_;
+    return numerical_value_;
   }
-#else
-  [[nodiscard]] constexpr rep& numerical_value() & noexcept { return value_; }
-  [[nodiscard]] constexpr const rep& numerical_value() const& noexcept { return value_; }
-  [[nodiscard]] constexpr rep&& numerical_value() && noexcept { return std::move(value_); }
-  [[nodiscard]] constexpr const rep&& numerical_value() const&& noexcept { return std::move(value_); }
-#endif
+
+  template<Unit U>
+    requires(U{} == unit)
+  [[nodiscard]] constexpr const rep& numerical_value_ref_in(U) const& noexcept
+  {
+    return numerical_value_;
+  }
+
+  template<Unit U>
+    requires(U{} == unit)
+  constexpr const rep&& numerical_value_ref_in(U) const&& noexcept = delete;
 
   template<Unit U>
     requires requires(quantity q) { q.in(U{}); }
   [[nodiscard]] constexpr rep numerical_value_in(U) const noexcept
   {
-    return (*this).in(U{}).numerical_value();
+    return (*this).in(U{}).numerical_value_;
   }
 
   template<Unit U>
-    requires detail::QuantityConvertibleTo<quantity, quantity<quantity_spec[U{}], Rep>>
-  [[nodiscard]] constexpr quantity<quantity_spec[U{}], Rep> in(U) const
+    requires requires(quantity q) { q.force_in(U{}); }
+  [[nodiscard]] constexpr rep force_numerical_value_in(U) const noexcept
   {
-    return quantity<quantity_spec[U{}], Rep>{*this};
+    return (*this).force_in(U{}).numerical_value_;
   }
 
   // member unary operators
@@ -177,7 +197,7 @@ class quantity {
       } -> std::common_with<rep>;
     }
   {
-    return make_quantity<reference>(+numerical_value());
+    return make_quantity<reference>(+numerical_value_);
   }
 
   [[nodiscard]] constexpr Quantity auto operator-() const
@@ -187,18 +207,19 @@ class quantity {
       } -> std::common_with<rep>;
     }
   {
-    return make_quantity<reference>(-numerical_value());
+    return make_quantity<reference>(-numerical_value_);
   }
 
-  constexpr quantity& operator++()
-    requires requires(rep v) {
+  template<typename Q>
+  friend constexpr decltype(auto) operator++(Q&& q)
+    requires std::derived_from<std::remove_cvref_t<Q>, quantity> && requires(rep v) {
       {
         ++v
       } -> std::same_as<rep&>;
     }
   {
-    ++value_;
-    return *this;
+    ++q.numerical_value_;
+    return std::forward<Q>(q);
   }
 
   [[nodiscard]] constexpr Quantity auto operator++(int)
@@ -208,18 +229,19 @@ class quantity {
       } -> std::common_with<rep>;
     }
   {
-    return make_quantity<reference>(value_++);
+    return make_quantity<reference>(numerical_value_++);
   }
 
-  constexpr quantity& operator--()
-    requires requires(rep v) {
+  template<typename Q>
+  friend constexpr decltype(auto) operator--(Q&& q)
+    requires std::derived_from<std::remove_cvref_t<Q>, quantity> && requires(rep v) {
       {
         --v
       } -> std::same_as<rep&>;
     }
   {
-    --value_;
-    return *this;
+    --q.numerical_value_;
+    return std::forward<Q>(q);
   }
 
   [[nodiscard]] constexpr Quantity auto operator--(int)
@@ -229,92 +251,101 @@ class quantity {
       } -> std::common_with<rep>;
     }
   {
-    return make_quantity<reference>(value_--);
+    return make_quantity<reference>(numerical_value_--);
   }
 
   // compound assignment operators
-  constexpr quantity& operator+=(const quantity& q)
-    requires requires(rep a, rep b) {
+  template<typename Q>
+    requires std::derived_from<std::remove_cvref_t<Q>, quantity> && requires(rep a, rep b) {
       {
         a += b
       } -> std::same_as<rep&>;
     }
+  friend constexpr decltype(auto) operator+=(Q&& lhs, const quantity& rhs)
   {
-    value_ += q.numerical_value();
-    return *this;
+    lhs.numerical_value_ += rhs.numerical_value_;
+    return std::forward<Q>(lhs);
   }
 
-  constexpr quantity& operator-=(const quantity& q)
-    requires requires(rep a, rep b) {
+  template<typename Q>
+    requires std::derived_from<std::remove_cvref_t<Q>, quantity> && requires(rep a, rep b) {
       {
         a -= b
       } -> std::same_as<rep&>;
     }
+  friend constexpr decltype(auto) operator-=(Q&& lhs, const quantity& rhs)
   {
-    value_ -= q.numerical_value();
-    return *this;
+    lhs.numerical_value_ -= rhs.numerical_value_;
+    return std::forward<Q>(lhs);
   }
 
-  constexpr quantity& operator%=(const quantity& q)
-    requires(!treat_as_floating_point<rep>) && requires(rep a, rep b) {
-      {
-        a %= b
-      } -> std::same_as<rep&>;
-    }
+  template<typename Q>
+    requires std::derived_from<std::remove_cvref_t<Q>, quantity> && (!treat_as_floating_point<rep>) &&
+             requires(rep a, rep b) {
+               {
+                 a %= b
+               } -> std::same_as<rep&>;
+             }
+  friend constexpr decltype(auto) operator%=(Q&& lhs, const quantity& rhs)
+
   {
-    gsl_ExpectsAudit(q != zero());
-    value_ %= q.numerical_value();
-    return *this;
+    gsl_ExpectsAudit(rhs != zero());
+    lhs.numerical_value_ %= rhs.numerical_value_;
+    return std::forward<Q>(lhs);
   }
 
-  template<typename Value>
-    requires(!Quantity<Value>) && requires(rep a, const Value b) {
-      {
-        a *= b
-      } -> std::same_as<rep&>;
-    }
-  constexpr quantity& operator*=(const Value& v)
+  template<typename Q, typename Value>
+    requires std::derived_from<std::remove_cvref_t<Q>, quantity> && (!Quantity<Value>) &&
+             requires(rep a, const Value b) {
+               {
+                 a *= b
+               } -> std::same_as<rep&>;
+             }
+  friend constexpr decltype(auto) operator*=(Q&& lhs, const Value& v)
   {
-    value_ *= v;
-    return *this;
+    lhs.numerical_value_ *= v;
+    return std::forward<Q>(lhs);
   }
 
-  template<QuantityOf<dimension_one> Q>
-    requires(Q::unit == ::mp_units::one) && requires(rep a, const typename Q::rep b) {
-      {
-        a *= b
-      } -> std::same_as<rep&>;
-    }
-  constexpr quantity& operator*=(const Q& rhs)
+  template<typename Q1, QuantityOf<dimension_one> Q2>
+    requires std::derived_from<std::remove_cvref_t<Q1>, quantity> && (Q2::unit == ::mp_units::one) &&
+             requires(rep a, const typename Q2::rep b) {
+               {
+                 a *= b
+               } -> std::same_as<rep&>;
+             }
+  friend constexpr decltype(auto) operator*=(Q1&& lhs, const Q2& rhs)
   {
-    value_ *= rhs.numerical_value();
-    return *this;
+    lhs.numerical_value_ *= rhs.numerical_value_;
+    return std::forward<Q1>(lhs);
   }
 
-  template<typename Value>
-    requires(!Quantity<Value>) && requires(rep a, const Value b) {
-      {
-        a /= b
-      } -> std::same_as<rep&>;
-    }
-  constexpr quantity& operator/=(const Value& v)
+  template<typename Q, typename Value>
+    requires std::derived_from<std::remove_cvref_t<Q>, quantity> && (!Quantity<Value>) &&
+             requires(rep a, const Value b) {
+               {
+                 a /= b
+               } -> std::same_as<rep&>;
+             }
+  friend constexpr decltype(auto) operator/=(Q&& lhs, const Value& v)
   {
     gsl_ExpectsAudit(v != quantity_values<Value>::zero());
-    value_ /= v;
-    return *this;
+    lhs.numerical_value_ /= v;
+    return std::forward<Q>(lhs);
   }
 
-  template<QuantityOf<dimension_one> Q>
-    requires(Q::unit == ::mp_units::one) && requires(rep a, const typename Q::rep b) {
-      {
-        a /= b
-      } -> std::same_as<rep&>;
-    }
-  constexpr quantity& operator/=(const Q& rhs)
+  template<typename Q1, QuantityOf<dimension_one> Q2>
+    requires std::derived_from<std::remove_cvref_t<Q1>, quantity> && (Q2::unit == ::mp_units::one) &&
+             requires(rep a, const typename Q2::rep b) {
+               {
+                 a /= b
+               } -> std::same_as<rep&>;
+             }
+  friend constexpr decltype(auto) operator/=(Q1&& lhs, const Q2& rhs)
   {
     gsl_ExpectsAudit(rhs != rhs.zero());
-    value_ /= rhs.numerical_value();
-    return *this;
+    lhs.numerical_value_ /= rhs.numerical_value_;
+    return std::forward<Q1>(lhs);
   }
 
 private:
@@ -327,7 +358,7 @@ class quantity {
 
   template<typename Value>
     requires std::constructible_from<rep, Value&&>
-  constexpr explicit quantity(Value&& v) : value_(std::forward<Value>(v))
+  constexpr explicit quantity(Value&& v) : numerical_value_(std::forward<Value>(v))
   {
   }
 };
@@ -342,7 +373,7 @@ template<auto R1, typename Rep1, auto R2, typename Rep2>
 [[nodiscard]] constexpr Quantity auto operator+(const quantity<R1, Rep1>& lhs, const quantity<R2, Rep2>& rhs)
 {
   using ret = detail::common_quantity_for<std::plus<>, quantity<R1, Rep1>, quantity<R2, Rep2>>;
-  return make_quantity<ret::reference>(ret(lhs).numerical_value() + ret(rhs).numerical_value());
+  return make_quantity<ret::reference>(ret(lhs).numerical_value_ + ret(rhs).numerical_value_);
 }
 
 template<auto R1, typename Rep1, auto R2, typename Rep2>
@@ -350,7 +381,7 @@ template<auto R1, typename Rep1, auto R2, typename Rep2>
 [[nodiscard]] constexpr Quantity auto operator-(const quantity<R1, Rep1>& lhs, const quantity<R2, Rep2>& rhs)
 {
   using ret = detail::common_quantity_for<std::minus<>, quantity<R1, Rep1>, quantity<R2, Rep2>>;
-  return make_quantity<ret::reference>(ret(lhs).numerical_value() - ret(rhs).numerical_value());
+  return make_quantity<ret::reference>(ret(lhs).numerical_value_ - ret(rhs).numerical_value_);
 }
 
 template<auto R1, typename Rep1, auto R2, typename Rep2>
@@ -360,7 +391,7 @@ template<auto R1, typename Rep1, auto R2, typename Rep2>
 {
   gsl_ExpectsAudit(rhs != rhs.zero());
   using ret = detail::common_quantity_for<std::modulus<>, quantity<R1, Rep1>, quantity<R2, Rep2>>;
-  return make_quantity<ret::reference>(ret(lhs).numerical_value() % ret(rhs).numerical_value());
+  return make_quantity<ret::reference>(ret(lhs).numerical_value_ % ret(rhs).numerical_value_);
 }
 
 template<auto R1, typename Rep1, auto R2, typename Rep2>
@@ -368,7 +399,7 @@ template<auto R1, typename Rep1, auto R2, typename Rep2>
                                   Rep2>
 [[nodiscard]] constexpr Quantity auto operator*(const quantity<R1, Rep1>& lhs, const quantity<R2, Rep2>& rhs)
 {
-  return make_quantity<R1 * R2>(lhs.numerical_value() * rhs.numerical_value());
+  return make_quantity<R1 * R2>(lhs.numerical_value_ * rhs.numerical_value_);
 }
 
 template<auto R, typename Rep, typename Value>
@@ -376,7 +407,7 @@ template<auto R, typename Rep, typename Value>
           detail::InvokeResultOf<get_quantity_spec(R).character, std::multiplies<>, Rep, const Value&>
 [[nodiscard]] constexpr Quantity auto operator*(const quantity<R, Rep>& q, const Value& v)
 {
-  return make_quantity<R>(q.numerical_value() * v);
+  return make_quantity<R>(q.numerical_value_ * v);
 }
 
 template<typename Value, auto R, typename Rep>
@@ -384,7 +415,7 @@ template<typename Value, auto R, typename Rep>
           detail::InvokeResultOf<get_quantity_spec(R).character, std::multiplies<>, const Value&, Rep>
 [[nodiscard]] constexpr Quantity auto operator*(const Value& v, const quantity<R, Rep>& q)
 {
-  return make_quantity<R>(v * q.numerical_value());
+  return make_quantity<R>(v * q.numerical_value_);
 }
 
 template<auto R1, typename Rep1, auto R2, typename Rep2>
@@ -392,7 +423,7 @@ template<auto R1, typename Rep1, auto R2, typename Rep2>
 [[nodiscard]] constexpr Quantity auto operator/(const quantity<R1, Rep1>& lhs, const quantity<R2, Rep2>& rhs)
 {
   gsl_ExpectsAudit(rhs != rhs.zero());
-  return make_quantity<R1 / R2>(lhs.numerical_value() / rhs.numerical_value());
+  return make_quantity<R1 / R2>(lhs.numerical_value_ / rhs.numerical_value_);
 }
 
 template<auto R, typename Rep, typename Value>
@@ -401,7 +432,7 @@ template<auto R, typename Rep, typename Value>
 [[nodiscard]] constexpr Quantity auto operator/(const quantity<R, Rep>& q, const Value& v)
 {
   gsl_ExpectsAudit(v != quantity_values<Value>::zero());
-  return make_quantity<R>(q.numerical_value() / v);
+  return make_quantity<R>(q.numerical_value_ / v);
 }
 
 template<typename Value, auto R, typename Rep>
@@ -409,7 +440,7 @@ template<typename Value, auto R, typename Rep>
           detail::InvokeResultOf<get_quantity_spec(R).character, std::divides<>, const Value&, Rep>
 [[nodiscard]] constexpr Quantity auto operator/(const Value& v, const quantity<R, Rep>& q)
 {
-  return make_quantity<::mp_units::one / R>(v / q.numerical_value());
+  return make_quantity<::mp_units::one / R>(v / q.numerical_value_);
 }
 
 template<auto R1, typename Rep1, auto R2, typename Rep2>
@@ -418,7 +449,7 @@ template<auto R1, typename Rep1, auto R2, typename Rep2>
 [[nodiscard]] constexpr bool operator==(const quantity<R1, Rep1>& lhs, const quantity<R2, Rep2>& rhs)
 {
   using ct = std::common_type_t<quantity<R1, Rep1>, quantity<R2, Rep2>>;
-  return ct(lhs).numerical_value() == ct(rhs).numerical_value();
+  return ct(lhs).numerical_value_ == ct(rhs).numerical_value_;
 }
 
 template<auto R1, typename Rep1, auto R2, typename Rep2>
@@ -427,7 +458,7 @@ template<auto R1, typename Rep1, auto R2, typename Rep2>
 [[nodiscard]] constexpr auto operator<=>(const quantity<R1, Rep1>& lhs, const quantity<R2, Rep2>& rhs)
 {
   using ct = std::common_type_t<quantity<R1, Rep1>, quantity<R2, Rep2>>;
-  return ct(lhs).numerical_value() <=> ct(rhs).numerical_value();
+  return ct(lhs).numerical_value_ <=> ct(rhs).numerical_value_;
 }
 
 // make_quantity
diff --git a/src/core/include/mp-units/quantity_point.h b/src/core/include/mp-units/quantity_point.h
index d736d2790..910620233 100644
--- a/src/core/include/mp-units/quantity_point.h
+++ b/src/core/include/mp-units/quantity_point.h
@@ -82,7 +82,7 @@ class quantity_point {
   using rep = Rep;
   using quantity_type = quantity<reference, Rep>;
 
-  quantity_type q_;  // needs to be public for a structural type
+  quantity_type quantity_from_origin_;  // needs to be public for a structural type
 
   // static member functions
   [[nodiscard]] static constexpr quantity_point zero() noexcept
@@ -112,10 +112,10 @@ class quantity_point {
     requires std::constructible_from<quantity_type, typename QP::quantity_type>
   // TODO add perfect forwarding
   constexpr explicit(!std::convertible_to<typename QP::quantity_type, quantity_type>) quantity_point(const QP& qp) :
-      q_([&] {
+      quantity_from_origin_([&] {
         if constexpr (is_same_v<std::remove_const_t<decltype(point_origin)>,
                                 std::remove_const_t<decltype(QP::point_origin)>>)
-          return qp.quantity_from_origin();
+          return qp.quantity_ref_from(point_origin);
         else
           return qp - point_origin;
       }())
@@ -128,7 +128,8 @@ class quantity_point {
              std::convertible_to<
                quantity<quantity_point_like_traits<QP>::reference, typename quantity_point_like_traits<QP>::rep>,
                quantity_type>
-  constexpr explicit quantity_point(const QP& qp) : q_(quantity_point_like_traits<QP>::quantity_from_origin(qp))
+  constexpr explicit quantity_point(const QP& qp) :
+      quantity_from_origin_(quantity_point_like_traits<QP>::quantity_from_origin(qp))
   {
   }
 
@@ -145,66 +146,88 @@ class quantity_point {
   }
 
   // data access
-#ifdef __cpp_explicit_this_parameter
-  template<typename Self>
-  [[nodiscard]] constexpr auto&& quantity_from_origin(this Self&& self) noexcept
+  template<std::same_as<std::remove_const_t<decltype(PO)>> PO2>
+  [[nodiscard]] constexpr quantity_type& quantity_ref_from(PO2) & noexcept
   {
-    return std::forward<Self>(self).q_;
+    return quantity_from_origin_;
+  }
+
+  template<std::same_as<std::remove_const_t<decltype(PO)>> PO2>
+  [[nodiscard]] constexpr const quantity_type& quantity_ref_from(PO2) const& noexcept
+  {
+    return quantity_from_origin_;
+  }
+
+  template<std::same_as<std::remove_const_t<decltype(PO)>> PO2>
+  constexpr const quantity_type&& quantity_ref_from(PO2) const&& noexcept = delete;
+
+  template<PointOrigin PO2>
+    requires requires { quantity_point{} - PO2{}; }
+  [[nodiscard]] constexpr Quantity auto quantity_from(PO2) const
+  {
+    return *this - PO2{};
+  }
+
+  template<Unit U>
+    requires detail::QuantityConvertibleTo<quantity_type, quantity<quantity_spec[U{}], Rep>>
+  [[nodiscard]] constexpr quantity_point<quantity_spec[U{}], PO, Rep> in(U) const
+  {
+    return make_quantity_point<PO>(quantity_ref_from(PO).in(U{}));
   }
-#else
-  [[nodiscard]] constexpr quantity_type& quantity_from_origin() & noexcept { return q_; }
-  [[nodiscard]] constexpr const quantity_type& quantity_from_origin() const& noexcept { return q_; }
-  [[nodiscard]] constexpr quantity_type&& quantity_from_origin() && noexcept { return std::move(q_); }
-  [[nodiscard]] constexpr const quantity_type&& quantity_from_origin() const&& noexcept { return std::move(q_); }
-#endif
 
   template<Unit U>
-    requires detail::QuantityConvertibleTo<quantity_type, quantity<::mp_units::reference<quantity_spec, U{}>{}, Rep>>
-  [[nodiscard]] constexpr quantity_point<::mp_units::reference<quantity_spec, U{}>{}, PO, Rep> in(U) const
+    requires requires(quantity_type q) { value_cast<U{}>(q); }
+  [[nodiscard]] constexpr quantity_point<quantity_spec[U{}], PO, Rep> force_in(U) const
   {
-    return make_quantity_point<PO>(quantity_from_origin().in(U{}));
+    return make_quantity_point<PO>(quantity_ref_from(PO).force_in(U{}));
   }
 
   // member unary operators
-  constexpr quantity_point& operator++()
-    requires requires { ++q_; }
+  template<typename QP>
+  friend constexpr decltype(auto) operator++(QP&& qp)
+    requires std::derived_from<std::remove_cvref_t<QP>, quantity_point> && requires { ++qp.quantity_from_origin_; }
   {
-    ++q_;
-    return *this;
+    ++qp.quantity_from_origin_;
+    return std::forward<QP>(qp);
   }
 
   [[nodiscard]] constexpr quantity_point operator++(int)
-    requires requires { q_++; }
+    requires requires { quantity_from_origin_++; }
   {
-    return quantity_point(q_++);
+    return quantity_point(quantity_from_origin_++);
   }
 
-  constexpr quantity_point& operator--()
-    requires requires { --q_; }
+  template<typename QP>
+  friend constexpr decltype(auto) operator--(QP&& qp)
+    requires std::derived_from<std::remove_cvref_t<QP>, quantity_point> && requires { --qp.quantity_from_origin_; }
   {
-    --q_;
-    return *this;
+    --qp.quantity_from_origin_;
+    return std::forward<QP>(qp);
   }
 
   [[nodiscard]] constexpr quantity_point operator--(int)
-    requires requires { q_--; }
+    requires requires { quantity_from_origin_--; }
   {
-    return quantity_point(q_--);
+    return quantity_point(quantity_from_origin_--);
   }
 
   // compound assignment operators
-  constexpr quantity_point& operator+=(const quantity_type& q)
-    requires requires { q_ += q; }
+  template<typename QP>
+    requires std::derived_from<std::remove_cvref_t<QP>, quantity_point> &&
+             requires(quantity_type q) { quantity_from_origin_ += q; }
+  friend constexpr decltype(auto) operator+=(QP&& qp, const quantity_type& q)
   {
-    q_ += q;
-    return *this;
+    qp.quantity_from_origin_ += q;
+    return std::forward<QP>(qp);
   }
 
-  constexpr quantity_point& operator-=(const quantity_type& q)
-    requires requires { q_ -= q; }
+  template<typename QP>
+    requires std::derived_from<std::remove_cvref_t<QP>, quantity_point> &&
+             requires(quantity_type q) { quantity_from_origin_ -= q; }
+  friend constexpr decltype(auto) operator-=(QP&& qp, const quantity_type& q)
   {
-    q_ -= q;
-    return *this;
+    qp.quantity_from_origin_ -= q;
+    return std::forward<QP>(qp);
   }
 
 private:
@@ -219,7 +242,7 @@ class quantity_point {
   template<Quantity Q>
     requires std::constructible_from<quantity_type, Q> &&
              ReferenceOf<std::remove_const_t<decltype(Q::reference)>, PO.quantity_spec>
-  constexpr explicit quantity_point(Q&& q) : q_(std::forward<Q>(q))
+  constexpr explicit quantity_point(Q&& q) : quantity_from_origin_(std::forward<Q>(q))
   {
   }
 };
@@ -235,9 +258,9 @@ template<auto R1, auto PO1, typename Rep1, auto R2, typename Rep2>
   requires ReferenceOf<std::remove_const_t<decltype(R2)>, PO1.quantity_spec>
 [[nodiscard]] constexpr QuantityPoint auto operator+(const quantity_point<R1, PO1, Rep1>& qp,
                                                      const quantity<R2, Rep2>& q)
-  requires requires { qp.quantity_from_origin() + q; }
+  requires requires { qp.quantity_ref_from(PO1) + q; }
 {
-  return make_quantity_point<PO1>(qp.quantity_from_origin() + q);
+  return make_quantity_point<PO1>(qp.quantity_ref_from(PO1) + q);
 }
 
 template<auto R1, typename Rep1, auto R2, auto PO2, typename Rep2>
@@ -245,7 +268,7 @@ template<auto R1, typename Rep1, auto R2, auto PO2, typename Rep2>
   requires ReferenceOf<std::remove_const_t<decltype(R1)>, PO2.quantity_spec>
 [[nodiscard]] constexpr QuantityPoint auto operator+(const quantity<R1, Rep1>& q,
                                                      const quantity_point<R2, PO2, Rep2>& qp)
-  requires requires { q + qp.quantity_from_origin(); }
+  requires requires { q + qp.quantity_ref_from(PO2); }
 {
   return qp + q;
 }
@@ -269,9 +292,9 @@ template<auto R1, auto PO1, typename Rep1, auto R2, typename Rep2>
   requires ReferenceOf<std::remove_const_t<decltype(R2)>, PO1.quantity_spec>
 [[nodiscard]] constexpr QuantityPoint auto operator-(const quantity_point<R1, PO1, Rep1>& qp,
                                                      const quantity<R2, Rep2>& q)
-  requires requires { qp.quantity_from_origin() - q; }
+  requires requires { qp.quantity_ref_from(PO1) - q; }
 {
-  return make_quantity_point<PO1>(qp.quantity_from_origin() - q);
+  return make_quantity_point<PO1>(qp.quantity_ref_from(PO1) - q);
 }
 
 template<PointOrigin PO, Quantity Q>
@@ -285,13 +308,14 @@ template<PointOrigin PO, Quantity Q>
 template<QuantityPoint QP1, QuantityPointOf<QP1::absolute_point_origin> QP2>
 [[nodiscard]] constexpr Quantity auto operator-(const QP1& lhs, const QP2& rhs)
   // TODO consider constraining it for both branches
-  requires requires { lhs.quantity_from_origin() - rhs.quantity_from_origin(); }
+  requires requires { lhs.quantity_ref_from(QP1::point_origin) - rhs.quantity_ref_from(QP2::point_origin); }
 {
   if constexpr (is_same_v<std::remove_const_t<decltype(QP1::point_origin)>,
                           std::remove_const_t<decltype(QP2::point_origin)>>)
-    return lhs.quantity_from_origin() - rhs.quantity_from_origin();
+    return lhs.quantity_ref_from(QP1::point_origin) - rhs.quantity_ref_from(QP2::point_origin);
   else
-    return lhs.quantity_from_origin() - rhs.quantity_from_origin() + (lhs.point_origin - rhs.point_origin);
+    return lhs.quantity_ref_from(QP1::point_origin) - rhs.quantity_ref_from(QP2::point_origin) +
+           (lhs.point_origin - rhs.point_origin);
 }
 
 template<PointOrigin PO, QuantityPointOf<PO{}> QP>
@@ -299,19 +323,21 @@ template<PointOrigin PO, QuantityPointOf<PO{}> QP>
 [[nodiscard]] constexpr Quantity auto operator-(const QP& qp, PO po)
 {
   if constexpr (is_same_v<std::remove_const_t<decltype(QP::point_origin)>, std::remove_const_t<PO>>)
-    return qp.quantity_from_origin();
+    return qp.quantity_ref_from(QP::point_origin);
   else if constexpr (detail::is_derived_from_specialization_of_absolute_point_origin<PO>) {
     if constexpr (is_same_v<std::remove_const_t<decltype(QP::point_origin)>,
                             std::remove_const_t<decltype(QP::absolute_point_origin)>>)
-      return qp.quantity_from_origin();
+      return qp.quantity_ref_from(QP::point_origin);
     else
-      return qp.quantity_from_origin() + (qp.point_origin - qp.absolute_point_origin);
+      return qp.quantity_ref_from(QP::point_origin) + (qp.point_origin - qp.absolute_point_origin);
   } else {
     if constexpr (is_same_v<std::remove_const_t<decltype(QP::point_origin)>,
                             std::remove_const_t<decltype(po.quantity_point.point_origin)>>)
-      return qp.quantity_from_origin() - po.quantity_point.quantity_from_origin();
+      return qp.quantity_ref_from(QP::point_origin) -
+             po.quantity_point.quantity_ref_from(po.quantity_point.point_origin);
     else
-      return qp.quantity_from_origin() - po.quantity_point.quantity_from_origin() +
+      return qp.quantity_ref_from(QP::point_origin) -
+             po.quantity_point.quantity_ref_from(po.quantity_point.point_origin) +
              (qp.point_origin - po.quantity_point.point_origin);
   }
 }
@@ -344,7 +370,7 @@ template<QuantityPoint QP1, QuantityPointOf<QP1::absolute_point_origin> QP2>
 {
   if constexpr (is_same_v<std::remove_const_t<decltype(QP1::point_origin)>,
                           std::remove_const_t<decltype(QP2::point_origin)>>)
-    return lhs.quantity_from_origin() <=> rhs.quantity_from_origin();
+    return lhs.quantity_ref_from(QP1::point_origin) <=> rhs.quantity_ref_from(QP2::point_origin);
   else
     return lhs - lhs.absolute_point_origin <=> rhs - rhs.absolute_point_origin;
 }
@@ -355,7 +381,7 @@ template<QuantityPoint QP1, QuantityPointOf<QP1::absolute_point_origin> QP2>
 {
   if constexpr (is_same_v<std::remove_const_t<decltype(QP1::point_origin)>,
                           std::remove_const_t<decltype(QP2::point_origin)>>)
-    return lhs.quantity_from_origin() == rhs.quantity_from_origin();
+    return lhs.quantity_ref_from(QP1::point_origin) == rhs.quantity_ref_from(QP2::point_origin);
   else
     return lhs - lhs.absolute_point_origin == rhs - rhs.absolute_point_origin;
 }
diff --git a/src/core/include/mp-units/quantity_spec.h b/src/core/include/mp-units/quantity_spec.h
index 0bd2fe95e..b6e2e2a1a 100644
--- a/src/core/include/mp-units/quantity_spec.h
+++ b/src/core/include/mp-units/quantity_spec.h
@@ -111,7 +111,7 @@ struct quantity_spec_interface {
     requires Quantity<std::remove_cvref_t<Q>> &&
              (explicitly_convertible(std::remove_reference_t<Q>::quantity_spec, self))
   {
-    return make_quantity<reference<self, std::remove_cvref_t<Q>::unit>{}>(std::forward<Q>(q).numerical_value());
+    return make_quantity<reference<self, std::remove_cvref_t<Q>::unit>{}>(std::forward<Q>(q).numerical_value_);
   }
 #else
   template<typename Self_ = Self, UnitOf<Self_{}> U>
@@ -128,7 +128,7 @@ struct quantity_spec_interface {
              (explicitly_convertible(std::remove_reference_t<Q>::quantity_spec, Self_{}))
   [[nodiscard]] constexpr Quantity auto operator()(Q&& q) const
   {
-    return make_quantity<reference<Self{}, std::remove_cvref_t<Q>::unit>{}>(std::forward<Q>(q).numerical_value());
+    return make_quantity<reference<Self{}, std::remove_cvref_t<Q>::unit>{}>(std::forward<Q>(q).numerical_value_);
   }
 #endif
 };
@@ -307,7 +307,7 @@ struct quantity_spec<Self, QS, Args...> : std::remove_const_t<decltype(QS)> {
              (explicitly_convertible(std::remove_reference_t<Q>::quantity_spec, Self_{}))
   [[nodiscard]] constexpr Quantity auto operator()(Q&& q) const
   {
-    return make_quantity<reference<Self{}, std::remove_cvref_t<Q>::unit>{}>(std::forward<Q>(q).numerical_value());
+    return make_quantity<reference<Self{}, std::remove_cvref_t<Q>::unit>{}>(std::forward<Q>(q).numerical_value_);
   }
 #endif
 };
diff --git a/src/core/include/mp-units/unit.h b/src/core/include/mp-units/unit.h
index 98adfa068..a15b0f7ca 100644
--- a/src/core/include/mp-units/unit.h
+++ b/src/core/include/mp-units/unit.h
@@ -402,7 +402,7 @@ using type_list_of_unit_less = expr_less<T1, T2, unit_less>;
  * Multiplication by `1` returns the same unit, otherwise `scaled_unit` is being returned.
  */
 template<Magnitude M, Unit U>
-[[nodiscard]] consteval Unit auto operator*(M mag, const U u)
+[[nodiscard]] MP_UNITS_CONSTEVAL Unit auto operator*(M mag, const U u)
 {
   if constexpr (mag == mp_units::mag<1>)
     return u;
@@ -418,7 +418,7 @@ template<Magnitude M, Unit U>
  * to the derived unit and the magnitude remains outside forming another scaled unit as a result of the operation.
  */
 template<Unit Lhs, Unit Rhs>
-[[nodiscard]] consteval Unit auto operator*(Lhs lhs, Rhs rhs)
+[[nodiscard]] MP_UNITS_CONSTEVAL Unit auto operator*(Lhs lhs, Rhs rhs)
 {
   if constexpr (detail::is_specialization_of_scaled_unit<Lhs> && detail::is_specialization_of_scaled_unit<Rhs>)
     return (Lhs::mag * Rhs::mag) * (Lhs::reference_unit * Rhs::reference_unit);
@@ -436,7 +436,7 @@ template<Unit Lhs, Unit Rhs>
  * to the derived unit and the magnitude remains outside forming another scaled unit as a result of the operation.
  */
 template<Unit Lhs, Unit Rhs>
-[[nodiscard]] consteval Unit auto operator/(Lhs lhs, Rhs rhs)
+[[nodiscard]] MP_UNITS_CONSTEVAL Unit auto operator/(Lhs lhs, Rhs rhs)
 {
   if constexpr (detail::is_specialization_of_scaled_unit<Lhs> && detail::is_specialization_of_scaled_unit<Rhs>)
     return (Lhs::mag / Rhs::mag) * (Lhs::reference_unit / Rhs::reference_unit);
diff --git a/src/utility/include/mp-units/chrono.h b/src/utility/include/mp-units/chrono.h
index aca9ad487..f868d84b0 100644
--- a/src/utility/include/mp-units/chrono.h
+++ b/src/utility/include/mp-units/chrono.h
@@ -92,7 +92,7 @@ template<QuantityOf<isq::time> Q>
 {
   constexpr auto canonical = detail::get_canonical_unit(Q::unit);
   constexpr ratio r = as_ratio(canonical.mag);
-  return std::chrono::duration<typename Q::rep, std::ratio<r.num, r.den>>{q.numerical_value()};
+  return std::chrono::duration<typename Q::rep, std::ratio<r.num, r.den>>{q.numerical_value_ref_in(Q::unit)};
 }
 
 template<QuantityPointOf<isq::time> QP>
diff --git a/src/utility/include/mp-units/math.h b/src/utility/include/mp-units/math.h
index e367ce4b2..7bdf71893 100644
--- a/src/utility/include/mp-units/math.h
+++ b/src/utility/include/mp-units/math.h
@@ -54,7 +54,8 @@ template<std::intmax_t Num, std::intmax_t Den = 1, auto R, typename Rep>
   requires detail::non_zero<Den> &&
            requires { quantity_values<Rep>::one(); }
            [[nodiscard]] constexpr quantity<pow<Num, Den>(R), Rep> pow(const quantity<R, Rep>& q) noexcept
-             requires requires { pow(q.numerical_value(), 1.0); } || requires { std::pow(q.numerical_value(), 1.0); }
+             requires requires { pow(q.numerical_value_ref_in(q.unit), 1.0); } ||
+                      requires { std::pow(q.numerical_value_ref_in(q.unit), 1.0); }
 {
   if constexpr (Num == 0) {
     return quantity<pow<Num, Den>(R), Rep>::one();
@@ -63,7 +64,7 @@ template<std::intmax_t Num, std::intmax_t Den = 1, auto R, typename Rep>
   } else {
     using std::pow;
     return make_quantity<pow<Num, Den>(R)>(
-      static_cast<Rep>(pow(q.numerical_value(), static_cast<double>(Num) / static_cast<double>(Den))));
+      static_cast<Rep>(pow(q.numerical_value_ref_in(q.unit), static_cast<double>(Num) / static_cast<double>(Den))));
   }
 }
 
@@ -77,10 +78,11 @@ template<std::intmax_t Num, std::intmax_t Den = 1, auto R, typename Rep>
  */
 template<auto R, typename Rep>
 [[nodiscard]] constexpr quantity<sqrt(R), Rep> sqrt(const quantity<R, Rep>& q) noexcept
-  requires requires { sqrt(q.numerical_value()); } || requires { std::sqrt(q.numerical_value()); }
+  requires requires { sqrt(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::sqrt(q.numerical_value_ref_in(q.unit)); }
 {
   using std::sqrt;
-  return make_quantity<sqrt(R)>(static_cast<Rep>(sqrt(q.numerical_value())));
+  return make_quantity<sqrt(R)>(static_cast<Rep>(sqrt(q.numerical_value_ref_in(q.unit))));
 }
 
 /**
@@ -93,10 +95,11 @@ template<auto R, typename Rep>
  */
 template<auto R, typename Rep>
 [[nodiscard]] constexpr quantity<cbrt(R), Rep> cbrt(const quantity<R, Rep>& q) noexcept
-  requires requires { cbrt(q.numerical_value()); } || requires { std::cbrt(q.numerical_value()); }
+  requires requires { cbrt(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::cbrt(q.numerical_value_ref_in(q.unit)); }
 {
   using std::cbrt;
-  return make_quantity<cbrt(R)>(static_cast<Rep>(cbrt(q.numerical_value())));
+  return make_quantity<cbrt(R)>(static_cast<Rep>(cbrt(q.numerical_value_ref_in(q.unit))));
 }
 
 /**
@@ -109,11 +112,12 @@ template<auto R, typename Rep>
  */
 template<ReferenceOf<dimensionless> auto R, typename Rep>
 [[nodiscard]] constexpr quantity<R, Rep> exp(const quantity<R, Rep>& q)
-  requires requires { exp(q.numerical_value()); } || requires { std::exp(q.numerical_value()); }
+  requires requires { exp(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::exp(q.numerical_value_ref_in(q.unit)); }
 {
   using std::exp;
   return value_cast<get_unit(R)>(
-    make_quantity<detail::clone_reference_with<one>(R)>(static_cast<Rep>(exp(value_cast<one>(q).numerical_value()))));
+    make_quantity<detail::clone_reference_with<one>(R)>(static_cast<Rep>(exp(q.force_numerical_value_in(q.unit)))));
 }
 
 /**
@@ -124,10 +128,11 @@ template<ReferenceOf<dimensionless> auto R, typename Rep>
  */
 template<auto R, typename Rep>
 [[nodiscard]] constexpr quantity<R, Rep> abs(const quantity<R, Rep>& q) noexcept
-  requires requires { abs(q.numerical_value()); } || requires { std::abs(q.numerical_value()); }
+  requires requires { abs(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::abs(q.numerical_value_ref_in(q.unit)); }
 {
   using std::abs;
-  return make_quantity<R>(static_cast<Rep>(abs(q.numerical_value())));
+  return make_quantity<R>(static_cast<Rep>(abs(q.numerical_value_ref_in(q.unit))));
 }
 
 /**
@@ -153,10 +158,10 @@ template<Representation Rep, Reference R>
  */
 template<Unit auto To, auto R, typename Rep>
 [[nodiscard]] constexpr quantity<detail::clone_reference_with<To>(R), Rep> floor(const quantity<R, Rep>& q) noexcept
-  requires((!treat_as_floating_point<Rep>) || requires { floor(q.numerical_value()); } ||
-           requires { std::floor(q.numerical_value()); }) &&
+  requires((!treat_as_floating_point<Rep>) || requires { floor(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::floor(q.numerical_value_ref_in(q.unit)); }) &&
           (To == get_unit(R) || requires {
-            ::mp_units::value_cast<To>(q);
+            q.force_in(To);
             quantity_values<Rep>::one();
           })
 {
@@ -169,16 +174,17 @@ template<Unit auto To, auto R, typename Rep>
   if constexpr (treat_as_floating_point<Rep>) {
     using std::floor;
     if constexpr (To == get_unit(R)) {
-      return make_quantity<detail::clone_reference_with<To>(R)>(static_cast<Rep>(floor(q.numerical_value())));
+      return make_quantity<detail::clone_reference_with<To>(R)>(
+        static_cast<Rep>(floor(q.numerical_value_ref_in(q.unit))));
     } else {
-      return handle_signed_results(make_quantity<detail::clone_reference_with<To>(R)>(
-        static_cast<Rep>(floor(value_cast<To>(q).numerical_value()))));
+      return handle_signed_results(
+        make_quantity<detail::clone_reference_with<To>(R)>(static_cast<Rep>(floor(q.force_numerical_value_in(To)))));
     }
   } else {
     if constexpr (To == get_unit(R)) {
-      return value_cast<To>(q);
+      return q.force_in(To);
     } else {
-      return handle_signed_results(value_cast<To>(q));
+      return handle_signed_results(q.force_in(To));
     }
   }
 }
@@ -191,10 +197,10 @@ template<Unit auto To, auto R, typename Rep>
  */
 template<Unit auto To, auto R, typename Rep>
 [[nodiscard]] constexpr quantity<detail::clone_reference_with<To>(R), Rep> ceil(const quantity<R, Rep>& q) noexcept
-  requires((!treat_as_floating_point<Rep>) || requires { ceil(q.numerical_value()); } ||
-           requires { std::ceil(q.numerical_value()); }) &&
+  requires((!treat_as_floating_point<Rep>) || requires { ceil(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::ceil(q.numerical_value_ref_in(q.unit)); }) &&
           (To == get_unit(R) || requires {
-            ::mp_units::value_cast<To>(q);
+            q.force_in(To);
             quantity_values<Rep>::one();
           })
 {
@@ -207,16 +213,17 @@ template<Unit auto To, auto R, typename Rep>
   if constexpr (treat_as_floating_point<Rep>) {
     using std::ceil;
     if constexpr (To == get_unit(R)) {
-      return make_quantity<detail::clone_reference_with<To>(R)>(static_cast<Rep>(ceil(q.numerical_value())));
+      return make_quantity<detail::clone_reference_with<To>(R)>(
+        static_cast<Rep>(ceil(q.numerical_value_ref_in(q.unit))));
     } else {
-      return handle_signed_results(make_quantity<detail::clone_reference_with<To>(R)>(
-        static_cast<Rep>(ceil(value_cast<To>(q).numerical_value()))));
+      return handle_signed_results(
+        make_quantity<detail::clone_reference_with<To>(R)>(static_cast<Rep>(ceil(q.force_numerical_value_in(To)))));
     }
   } else {
     if constexpr (To == get_unit(R)) {
-      return value_cast<To>(q);
+      return q.force_in(To);
     } else {
-      return handle_signed_results(value_cast<To>(q));
+      return handle_signed_results(q.force_in(To));
     }
   }
 }
@@ -231,8 +238,8 @@ template<Unit auto To, auto R, typename Rep>
  */
 template<Unit auto To, auto R, typename Rep>
 [[nodiscard]] constexpr quantity<detail::clone_reference_with<To>(R), Rep> round(const quantity<R, Rep>& q) noexcept
-  requires((!treat_as_floating_point<Rep>) || requires { round(q.numerical_value()); } ||
-           requires { std::round(q.numerical_value()); }) &&
+  requires((!treat_as_floating_point<Rep>) || requires { round(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::round(q.numerical_value_ref_in(q.unit)); }) &&
           (To == get_unit(R) || requires {
             ::mp_units::floor<To>(q);
             quantity_values<Rep>::one();
@@ -241,9 +248,10 @@ template<Unit auto To, auto R, typename Rep>
   if constexpr (To == get_unit(R)) {
     if constexpr (treat_as_floating_point<Rep>) {
       using std::round;
-      return make_quantity<detail::clone_reference_with<To>(R)>(static_cast<Rep>(round(q.numerical_value())));
+      return make_quantity<detail::clone_reference_with<To>(R)>(
+        static_cast<Rep>(round(q.numerical_value_ref_in(q.unit))));
     } else {
-      return value_cast<To>(q);
+      return q.force_in(To);
     }
   } else {
     const auto res_low = mp_units::floor<To>(q);
@@ -251,7 +259,7 @@ template<Unit auto To, auto R, typename Rep>
     const auto diff0 = q - res_low;
     const auto diff1 = res_high - q;
     if (diff0 == diff1) {
-      if (static_cast<int>(res_low.numerical_value()) & 1) {
+      if (static_cast<int>(res_low.numerical_value_ref_in(To)) & 1) {
         return res_high;
       }
       return res_low;
@@ -272,8 +280,8 @@ template<auto R1, typename Rep1, auto R2, typename Rep2>
   const quantity<R1, Rep1>& x, const quantity<R2, Rep2>& y) noexcept
   requires requires { common_reference(R1, R2); } &&
            (
-             requires { hypot(x.numerical_value(), y.numerical_value()); } ||
-             requires { std::hypot(x.numerical_value(), y.numerical_value()); })
+             requires { hypot(x.numerical_value_ref_in(x.unit), y.numerical_value_ref_in(y.unit)); } ||
+             requires { std::hypot(x.numerical_value_ref_in(x.unit), y.numerical_value_ref_in(y.unit)); })
 {
   constexpr auto ref = common_reference(R1, R2);
   constexpr auto unit = get_unit(ref);
@@ -290,8 +298,14 @@ template<auto R1, typename Rep1, auto R2, typename Rep2, auto R3, typename Rep3>
   const quantity<R1, Rep1>& x, const quantity<R2, Rep2>& y, const quantity<R3, Rep3>& z) noexcept
   requires requires { common_reference(R1, R2, R3); } &&
            (
-             requires { hypot(x.numerical_value(), y.numerical_value(), z.numerical_value()); } ||
-             requires { std::hypot(x.numerical_value(), y.numerical_value(), z.numerical_value()); })
+             requires {
+               hypot(x.numerical_value_ref_in(x.unit), y.numerical_value_ref_in(y.unit),
+                     z.numerical_value_ref_in(z.unit));
+             } ||
+             requires {
+               std::hypot(x.numerical_value_ref_in(x.unit), y.numerical_value_ref_in(y.unit),
+                          z.numerical_value_ref_in(z.unit));
+             })
 {
   constexpr auto ref = common_reference(R1, R2);
   constexpr auto unit = get_unit(ref);
@@ -303,12 +317,13 @@ namespace isq {
 
 template<ReferenceOf<angular_measure> auto R, typename Rep>
 [[nodiscard]] inline QuantityOf<dimensionless> auto sin(const quantity<R, Rep>& q) noexcept
-  requires requires { sin(q.numerical_value()); } || requires { std::sin(q.numerical_value()); }
+  requires requires { sin(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::sin(q.numerical_value_ref_in(q.unit)); }
 {
   using std::sin;
   if constexpr (!treat_as_floating_point<Rep>) {
     // check what is the return type when called with the integral value
-    using rep = decltype(sin(value_cast<si::radian>(q).numerical_value()));
+    using rep = decltype(sin(q.force_numerical_value_in(si::radian)));
     // use this type ahead of calling the function to prevent narrowing if a unit conversion is needed
     return make_quantity<one>(sin(value_cast<rep>(q).numerical_value_in(si::radian)));
   } else
@@ -317,12 +332,13 @@ template<ReferenceOf<angular_measure> auto R, typename Rep>
 
 template<ReferenceOf<angular_measure> auto R, typename Rep>
 [[nodiscard]] inline QuantityOf<dimensionless> auto cos(const quantity<R, Rep>& q) noexcept
-  requires requires { cos(q.numerical_value()); } || requires { std::cos(q.numerical_value()); }
+  requires requires { cos(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::cos(q.numerical_value_ref_in(q.unit)); }
 {
   using std::cos;
   if constexpr (!treat_as_floating_point<Rep>) {
     // check what is the return type when called with the integral value
-    using rep = decltype(cos(value_cast<si::radian>(q).numerical_value()));
+    using rep = decltype(cos(q.force_numerical_value_in(si::radian)));
     // use this type ahead of calling the function to prevent narrowing if a unit conversion is needed
     return make_quantity<one>(cos(value_cast<rep>(q).numerical_value_in(si::radian)));
   } else
@@ -331,12 +347,13 @@ template<ReferenceOf<angular_measure> auto R, typename Rep>
 
 template<ReferenceOf<angular_measure> auto R, typename Rep>
 [[nodiscard]] inline QuantityOf<dimensionless> auto tan(const quantity<R, Rep>& q) noexcept
-  requires requires { tan(q.numerical_value()); } || requires { std::tan(q.numerical_value()); }
+  requires requires { tan(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::tan(q.numerical_value_ref_in(q.unit)); }
 {
   using std::tan;
   if constexpr (!treat_as_floating_point<Rep>) {
     // check what is the return type when called with the integral value
-    using rep = decltype(tan(value_cast<si::radian>(q).numerical_value()));
+    using rep = decltype(tan(q.force_numerical_value_in(si::radian)));
     // use this type ahead of calling the function to prevent narrowing if a unit conversion is needed
     return make_quantity<one>(tan(value_cast<rep>(q).numerical_value_in(si::radian)));
   } else
@@ -345,12 +362,13 @@ template<ReferenceOf<angular_measure> auto R, typename Rep>
 
 template<ReferenceOf<dimensionless> auto R, typename Rep>
 [[nodiscard]] inline QuantityOf<isq::angular_measure> auto asin(const quantity<R, Rep>& q) noexcept
-  requires requires { asin(q.numerical_value()); } || requires { std::asin(q.numerical_value()); }
+  requires requires { asin(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::asin(q.numerical_value_ref_in(q.unit)); }
 {
   using std::asin;
   if constexpr (!treat_as_floating_point<Rep>) {
     // check what is the return type when called with the integral value
-    using rep = decltype(asin(value_cast<one>(q).numerical_value()));
+    using rep = decltype(asin(q.force_numerical_value_in(one)));
     // use this type ahead of calling the function to prevent narrowing if a unit conversion is needed
     return make_quantity<si::radian>(asin(value_cast<rep>(q).numerical_value_in(one)));
   } else
@@ -359,12 +377,13 @@ template<ReferenceOf<dimensionless> auto R, typename Rep>
 
 template<ReferenceOf<dimensionless> auto R, typename Rep>
 [[nodiscard]] inline QuantityOf<isq::angular_measure> auto acos(const quantity<R, Rep>& q) noexcept
-  requires requires { acos(q.numerical_value()); } || requires { std::acos(q.numerical_value()); }
+  requires requires { acos(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::acos(q.numerical_value_ref_in(q.unit)); }
 {
   using std::acos;
   if constexpr (!treat_as_floating_point<Rep>) {
     // check what is the return type when called with the integral value
-    using rep = decltype(acos(value_cast<one>(q).numerical_value()));
+    using rep = decltype(acos(q.force_numerical_value_in(one)));
     // use this type ahead of calling the function to prevent narrowing if a unit conversion is needed
     return make_quantity<si::radian>(acos(value_cast<rep>(q).numerical_value_in(one)));
   } else
@@ -373,12 +392,13 @@ template<ReferenceOf<dimensionless> auto R, typename Rep>
 
 template<ReferenceOf<dimensionless> auto R, typename Rep>
 [[nodiscard]] inline QuantityOf<isq::angular_measure> auto atan(const quantity<R, Rep>& q) noexcept
-  requires requires { atan(q.numerical_value()); } || requires { std::atan(q.numerical_value()); }
+  requires requires { atan(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::atan(q.numerical_value_ref_in(q.unit)); }
 {
   using std::atan;
   if constexpr (!treat_as_floating_point<Rep>) {
     // check what is the return type when called with the integral value
-    using rep = decltype(atan(value_cast<one>(q).numerical_value()));
+    using rep = decltype(atan(q.force_numerical_value_in(one)));
     // use this type ahead of calling the function to prevent narrowing if a unit conversion is needed
     return make_quantity<si::radian>(atan(value_cast<rep>(q).numerical_value_in(one)));
   } else
@@ -391,12 +411,13 @@ namespace angular {
 
 template<ReferenceOf<angle> auto R, typename Rep>
 [[nodiscard]] inline QuantityOf<dimensionless> auto sin(const quantity<R, Rep>& q) noexcept
-  requires requires { sin(q.numerical_value()); } || requires { std::sin(q.numerical_value()); }
+  requires requires { sin(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::sin(q.numerical_value_ref_in(q.unit)); }
 {
   using std::sin;
   if constexpr (!treat_as_floating_point<Rep>) {
     // check what is the return type when called with the integral value
-    using rep = decltype(sin(value_cast<radian>(q).numerical_value()));
+    using rep = decltype(sin(q.force_numerical_value_in(radian)));
     // use this type ahead of calling the function to prevent narrowing if a unit conversion is needed
     return make_quantity<one>(sin(value_cast<rep>(q).numerical_value_in(radian)));
   } else
@@ -405,12 +426,13 @@ template<ReferenceOf<angle> auto R, typename Rep>
 
 template<ReferenceOf<angle> auto R, typename Rep>
 [[nodiscard]] inline QuantityOf<dimensionless> auto cos(const quantity<R, Rep>& q) noexcept
-  requires requires { cos(q.numerical_value()); } || requires { std::cos(q.numerical_value()); }
+  requires requires { cos(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::cos(q.numerical_value_ref_in(q.unit)); }
 {
   using std::cos;
   if constexpr (!treat_as_floating_point<Rep>) {
     // check what is the return type when called with the integral value
-    using rep = decltype(cos(value_cast<radian>(q).numerical_value()));
+    using rep = decltype(cos(q.force_numerical_value_in(radian)));
     // use this type ahead of calling the function to prevent narrowing if a unit conversion is needed
     return make_quantity<one>(cos(value_cast<rep>(q).numerical_value_in(radian)));
   } else
@@ -419,12 +441,13 @@ template<ReferenceOf<angle> auto R, typename Rep>
 
 template<ReferenceOf<angle> auto R, typename Rep>
 [[nodiscard]] inline QuantityOf<dimensionless> auto tan(const quantity<R, Rep>& q) noexcept
-  requires requires { tan(q.numerical_value()); } || requires { std::tan(q.numerical_value()); }
+  requires requires { tan(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::tan(q.numerical_value_ref_in(q.unit)); }
 {
   using std::tan;
   if constexpr (!treat_as_floating_point<Rep>) {
     // check what is the return type when called with the integral value
-    using rep = decltype(tan(value_cast<radian>(q).numerical_value()));
+    using rep = decltype(tan(q.force_numerical_value_in(radian)));
     // use this type ahead of calling the function to prevent narrowing if a unit conversion is needed
     return make_quantity<one>(tan(value_cast<rep>(q).numerical_value_in(radian)));
   } else
@@ -433,12 +456,13 @@ template<ReferenceOf<angle> auto R, typename Rep>
 
 template<ReferenceOf<dimensionless> auto R, typename Rep>
 [[nodiscard]] inline QuantityOf<angle> auto asin(const quantity<R, Rep>& q) noexcept
-  requires requires { asin(q.numerical_value()); } || requires { std::asin(q.numerical_value()); }
+  requires requires { asin(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::asin(q.numerical_value_ref_in(q.unit)); }
 {
   using std::asin;
   if constexpr (!treat_as_floating_point<Rep>) {
     // check what is the return type when called with the integral value
-    using rep = decltype(asin(value_cast<one>(q).numerical_value()));
+    using rep = decltype(asin(q.force_numerical_value_in(one)));
     // use this type ahead of calling the function to prevent narrowing if a unit conversion is needed
     return make_quantity<radian>(asin(value_cast<rep>(q).numerical_value_in(one)));
   } else
@@ -447,12 +471,13 @@ template<ReferenceOf<dimensionless> auto R, typename Rep>
 
 template<ReferenceOf<dimensionless> auto R, typename Rep>
 [[nodiscard]] inline QuantityOf<angle> auto acos(const quantity<R, Rep>& q) noexcept
-  requires requires { acos(q.numerical_value()); } || requires { std::acos(q.numerical_value()); }
+  requires requires { acos(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::acos(q.numerical_value_ref_in(q.unit)); }
 {
   using std::acos;
   if constexpr (!treat_as_floating_point<Rep>) {
     // check what is the return type when called with the integral value
-    using rep = decltype(acos(value_cast<one>(q).numerical_value()));
+    using rep = decltype(acos(q.force_numerical_value_in(one)));
     // use this type ahead of calling the function to prevent narrowing if a unit conversion is needed
     return make_quantity<radian>(acos(value_cast<rep>(q).numerical_value_in(one)));
   } else
@@ -461,12 +486,13 @@ template<ReferenceOf<dimensionless> auto R, typename Rep>
 
 template<ReferenceOf<dimensionless> auto R, typename Rep>
 [[nodiscard]] inline QuantityOf<angle> auto atan(const quantity<R, Rep>& q) noexcept
-  requires requires { atan(q.numerical_value()); } || requires { std::atan(q.numerical_value()); }
+  requires requires { atan(q.numerical_value_ref_in(q.unit)); } ||
+           requires { std::atan(q.numerical_value_ref_in(q.unit)); }
 {
   using std::atan;
   if constexpr (!treat_as_floating_point<Rep>) {
     // check what is the return type when called with the integral value
-    using rep = decltype(atan(value_cast<one>(q).numerical_value()));
+    using rep = decltype(atan(q.force_numerical_value_in(one)));
     // use this type ahead of calling the function to prevent narrowing if a unit conversion is needed
     return make_quantity<radian>(atan(value_cast<rep>(q).numerical_value_in(one)));
   } else
diff --git a/src/utility/include/mp-units/random.h b/src/utility/include/mp-units/random.h
index 2c26b97a9..fd2ceda9a 100644
--- a/src/utility/include/mp-units/random.h
+++ b/src/utility/include/mp-units/random.h
@@ -35,7 +35,7 @@ static std::vector<typename Q::rep> i_qty_to_rep(InputIt first, InputIt last)
   std::vector<typename Q::rep> intervals_rep;
   intervals_rep.reserve(static_cast<size_t>(std::distance(first, last)));
   for (auto itr = first; itr != last; ++itr) {
-    intervals_rep.push_back(itr->numerical_value());
+    intervals_rep.push_back(itr->numerical_value_ref_in(Q::unit));
   }
   return intervals_rep;
 }
@@ -46,7 +46,7 @@ static std::vector<typename Q::rep> bl_qty_to_rep(std::initializer_list<Q>& bl)
   std::vector<typename Q::rep> bl_rep;
   bl_rep.reserve(bl.size());
   for (const Q& qty : bl) {
-    bl_rep.push_back(qty.numerical_value());
+    bl_rep.push_back(qty.numerical_value_ref_in(Q::unit));
   }
   return bl_rep;
 }
@@ -88,7 +88,10 @@ struct uniform_int_distribution : public std::uniform_int_distribution<typename
   using base = MP_UNITS_TYPENAME std::uniform_int_distribution<rep>;
 
   uniform_int_distribution() : base() {}
-  uniform_int_distribution(const Q& a, const Q& b) : base(a.numerical_value(), b.numerical_value()) {}
+  uniform_int_distribution(const Q& a, const Q& b) :
+      base(a.numerical_value_ref_in(Q::unit), b.numerical_value_ref_in(Q::unit))
+  {
+  }
 
   template<typename Generator>
   Q operator()(Generator& g)
@@ -110,7 +113,10 @@ struct uniform_real_distribution : public std::uniform_real_distribution<typenam
   using base = MP_UNITS_TYPENAME std::uniform_real_distribution<rep>;
 
   uniform_real_distribution() : base() {}
-  uniform_real_distribution(const Q& a, const Q& b) : base(a.numerical_value(), b.numerical_value()) {}
+  uniform_real_distribution(const Q& a, const Q& b) :
+      base(a.numerical_value_ref_in(Q::unit), b.numerical_value_ref_in(Q::unit))
+  {
+  }
 
   template<typename Generator>
   Q operator()(Generator& g)
@@ -132,7 +138,7 @@ struct binomial_distribution : public std::binomial_distribution<typename Q::rep
   using base = MP_UNITS_TYPENAME std::binomial_distribution<rep>;
 
   binomial_distribution() : base() {}
-  binomial_distribution(const Q& t, double p) : base(t.numerical_value(), p) {}
+  binomial_distribution(const Q& t, double p) : base(t.numerical_value_ref_in(Q::unit), p) {}
 
   template<typename Generator>
   Q operator()(Generator& g)
@@ -153,7 +159,7 @@ struct negative_binomial_distribution : public std::negative_binomial_distributi
   using base = MP_UNITS_TYPENAME std::negative_binomial_distribution<rep>;
 
   negative_binomial_distribution() : base() {}
-  negative_binomial_distribution(const Q& k, double p) : base(k.numerical_value(), p) {}
+  negative_binomial_distribution(const Q& k, double p) : base(k.numerical_value_ref_in(Q::unit), p) {}
 
   template<typename Generator>
   Q operator()(Generator& g)
@@ -269,7 +275,7 @@ struct extreme_value_distribution : public std::extreme_value_distribution<typen
   using base = MP_UNITS_TYPENAME std::extreme_value_distribution<rep>;
 
   extreme_value_distribution() : base() {}
-  extreme_value_distribution(const Q& a, const rep& b) : base(a.numerical_value(), b) {}
+  extreme_value_distribution(const Q& a, const rep& b) : base(a.numerical_value_ref_in(Q::unit), b) {}
 
   template<typename Generator>
   Q operator()(Generator& g)
@@ -290,7 +296,10 @@ struct normal_distribution : public std::normal_distribution<typename Q::rep> {
   using base = MP_UNITS_TYPENAME std::normal_distribution<rep>;
 
   normal_distribution() : base() {}
-  normal_distribution(const Q& mean, const Q& stddev) : base(mean.numerical_value(), stddev.numerical_value()) {}
+  normal_distribution(const Q& mean, const Q& stddev) :
+      base(mean.numerical_value_ref_in(Q::unit), stddev.numerical_value_ref_in(Q::unit))
+  {
+  }
 
   template<typename Generator>
   Q operator()(Generator& g)
@@ -312,7 +321,10 @@ struct lognormal_distribution : public std::lognormal_distribution<typename Q::r
   using base = MP_UNITS_TYPENAME std::lognormal_distribution<rep>;
 
   lognormal_distribution() : base() {}
-  lognormal_distribution(const Q& m, const Q& s) : base(m.numerical_value(), s.numerical_value()) {}
+  lognormal_distribution(const Q& m, const Q& s) :
+      base(m.numerical_value_ref_in(Q::unit), s.numerical_value_ref_in(Q::unit))
+  {
+  }
 
   template<typename Generator>
   Q operator()(Generator& g)
@@ -353,7 +365,10 @@ struct cauchy_distribution : public std::cauchy_distribution<typename Q::rep> {
   using base = MP_UNITS_TYPENAME std::cauchy_distribution<rep>;
 
   cauchy_distribution() : base() {}
-  cauchy_distribution(const Q& a, const Q& b) : base(a.numerical_value(), b.numerical_value()) {}
+  cauchy_distribution(const Q& a, const Q& b) :
+      base(a.numerical_value_ref_in(Q::unit), b.numerical_value_ref_in(Q::unit))
+  {
+  }
 
   template<typename Generator>
   Q operator()(Generator& g)
@@ -470,7 +485,8 @@ class piecewise_constant_distribution : public std::piecewise_constant_distribut
 
   template<typename UnaryOperation>
   piecewise_constant_distribution(std::size_t nw, const Q& xmin, const Q& xmax, UnaryOperation fw) :
-      base(nw, xmin.numerical_value(), xmax.numerical_value(), [fw](rep val) { return fw(val * Q::reference); })
+      base(nw, xmin.numerical_value_ref_in(Q::unit), xmax.numerical_value_ref_in(Q::unit),
+           [fw](rep val) { return fw(val * Q::reference); })
   {
   }
 
@@ -528,7 +544,8 @@ class piecewise_linear_distribution : public std::piecewise_linear_distribution<
 
   template<typename UnaryOperation>
   piecewise_linear_distribution(std::size_t nw, const Q& xmin, const Q& xmax, UnaryOperation fw) :
-      base(nw, xmin.numerical_value(), xmax.numerical_value(), [fw](rep val) { return fw(val * Q::reference); })
+      base(nw, xmin.numerical_value_ref_in(Q::unit), xmax.numerical_value_ref_in(Q::unit),
+           [fw](rep val) { return fw(val * Q::reference); })
   {
   }
 
diff --git a/test/unit_test/runtime/almost_equals.h b/test/unit_test/runtime/almost_equals.h
index 69b972c1e..b3af6689b 100644
--- a/test/unit_test/runtime/almost_equals.h
+++ b/test/unit_test/runtime/almost_equals.h
@@ -36,8 +36,8 @@ struct AlmostEqualsMatcher : Catch::Matchers::MatcherGenericBase {
   {
     using std::abs;
     using common = std::common_type_t<T, U>;
-    const auto x = common(target_).numerical_value();
-    const auto y = common(other).numerical_value();
+    const auto x = common(target_).numerical_value_in(common::unit);
+    const auto y = common(other).numerical_value_in(common::unit);
     const auto maxXYOne = std::max({typename T::rep{1}, abs(x), abs(y)});
     return abs(x - y) <= std::numeric_limits<typename T::rep>::epsilon() * maxXYOne;
   }
diff --git a/test/unit_test/runtime/distribution_test.cpp b/test/unit_test/runtime/distribution_test.cpp
index b7f459713..5b14ae7f1 100644
--- a/test/unit_test/runtime/distribution_test.cpp
+++ b/test/unit_test/runtime/distribution_test.cpp
@@ -559,8 +559,8 @@ TEST_CASE("piecewise_constant_distribution")
                                               3.0 * isq::length[si::metre]};
 
     auto stl_dist = std::piecewise_constant_distribution<rep>(intervals_rep, [](rep val) { return val; });
-    auto units_dist =
-      mp_units::piecewise_constant_distribution<q>(intervals_qty, [](q qty) { return qty.numerical_value(); });
+    auto units_dist = mp_units::piecewise_constant_distribution<q>(
+      intervals_qty, [](q qty) { return qty.numerical_value_ref_in(qty.unit); });
 
     CHECK(units_dist.intervals() == intervals_qty_vec);
     CHECK(units_dist.densities() == stl_dist.densities());
@@ -573,8 +573,8 @@ TEST_CASE("piecewise_constant_distribution")
     constexpr q xmin_qty = 1.0 * isq::length[si::metre], xmax_qty = 3.0 * isq::length[si::metre];
 
     auto stl_dist = std::piecewise_constant_distribution<rep>(nw, xmin_rep, xmax_rep, [](rep val) { return val; });
-    auto units_dist =
-      mp_units::piecewise_constant_distribution<q>(nw, xmin_qty, xmax_qty, [](q qty) { return qty.numerical_value(); });
+    auto units_dist = mp_units::piecewise_constant_distribution<q>(
+      nw, xmin_qty, xmax_qty, [](q qty) { return qty.numerical_value_ref_in(qty.unit); });
 
     CHECK(units_dist.intervals() == intervals_qty_vec);
     CHECK(units_dist.densities() == stl_dist.densities());
@@ -627,8 +627,8 @@ TEST_CASE("piecewise_linear_distribution")
                                               3.0 * isq::length[si::metre]};
 
     auto stl_dist = std::piecewise_linear_distribution<rep>(intervals_rep, [](rep val) { return val; });
-    auto units_dist =
-      mp_units::piecewise_linear_distribution<q>(intervals_qty, [](q qty) { return qty.numerical_value(); });
+    auto units_dist = mp_units::piecewise_linear_distribution<q>(
+      intervals_qty, [](q qty) { return qty.numerical_value_ref_in(qty.unit); });
 
     CHECK(units_dist.intervals() == intervals_qty_vec);
     CHECK(units_dist.densities() == stl_dist.densities());
@@ -641,8 +641,8 @@ TEST_CASE("piecewise_linear_distribution")
     constexpr q xmin_qty = 1.0 * isq::length[si::metre], xmax_qty = 3.0 * isq::length[si::metre];
 
     auto stl_dist = std::piecewise_linear_distribution<rep>(nw, xmin_rep, xmax_rep, [](rep val) { return val; });
-    auto units_dist =
-      mp_units::piecewise_linear_distribution<q>(nw, xmin_qty, xmax_qty, [](q qty) { return qty.numerical_value(); });
+    auto units_dist = mp_units::piecewise_linear_distribution<q>(
+      nw, xmin_qty, xmax_qty, [](q qty) { return qty.numerical_value_ref_in(qty.unit); });
 
     CHECK(units_dist.intervals() == intervals_qty_vec);
     CHECK(units_dist.densities() == stl_dist.densities());
diff --git a/test/unit_test/runtime/linear_algebra_test.cpp b/test/unit_test/runtime/linear_algebra_test.cpp
index 818031f4b..ad9fb31f4 100644
--- a/test/unit_test/runtime/linear_algebra_test.cpp
+++ b/test/unit_test/runtime/linear_algebra_test.cpp
@@ -56,11 +56,19 @@ namespace {
 using namespace mp_units;
 using namespace mp_units::si::unit_symbols;
 
-template<typename T>
-[[nodiscard]] auto get_magnitude(const vector<T>& v)
+template<QuantitySpec auto QS, QuantityOf<QS> Q>
+  requires(Q::quantity_spec.character == quantity_character::vector) && (QS.character == quantity_character::scalar)
+[[nodiscard]] QuantityOf<QS> auto get_magnitude(const Q& q)
 {
-  using namespace std;
-  return hypot(v(0), v(1), v(2));
+  const auto& v = q.numerical_value_ref_in(q.unit);
+  return hypot(v(0) * QS[q.unit], v(1) * QS[q.unit], v(2) * QS[q.unit]);
+}
+
+template<QuantitySpec auto QS, QuantityOf<QS> T>
+  requires(T::quantity_spec.character == quantity_character::vector) && (QS.character == quantity_character::scalar)
+[[nodiscard]] QuantityOf<QS> auto get_magnitude(const vector<T>& v)
+{
+  return hypot(QS(v(0)), QS(v(1)), QS(v(2)));
 }
 
 template<typename T, typename U>
@@ -74,7 +82,8 @@ template<Quantity Q1, Quantity Q2>
            requires(typename Q1::rep v1, typename Q2::rep v2) { cross_product(v1, v2); }
 [[nodiscard]] QuantityOf<Q1::quantity_spec * Q2::quantity_spec> auto cross_product(const Q1& q1, const Q2& q2)
 {
-  return cross_product(q1.numerical_value(), q2.numerical_value()) * (Q1::reference * Q2::reference);
+  return cross_product(q1.numerical_value_ref_in(q1.unit), q2.numerical_value_ref_in(q2.unit)) *
+         (Q1::reference * Q2::reference);
 }
 
 }  // namespace
@@ -86,21 +95,21 @@ TEST_CASE("vector quantity", "[la]")
     SECTION("non-truncating")
     {
       const auto v = vector<int>{3, 2, 1} * isq::position_vector[km];
-      CHECK(v.in(m).numerical_value() == vector<int>{3000, 2000, 1000});
+      CHECK(v.numerical_value_in(m) == vector<int>{3000, 2000, 1000});
     }
 
     SECTION("truncating")
     {
       const auto v = vector<int>{1001, 1002, 1003} * isq::position_vector[m];
-      CHECK(value_cast<km>(v).numerical_value() == vector<int>{1, 1, 1});
+      CHECK(v.force_numerical_value_in(km) == vector<int>{1, 1, 1});
     }
   }
 
   SECTION("to scalar magnitude")
   {
     const auto v = vector<int>{2, 3, 6} * isq::velocity[km / h];
-    const auto speed = get_magnitude(v.numerical_value()) * isq::speed[v.unit];  // TODO can we do better here?
-    CHECK(speed.numerical_value() == 7);
+    const auto speed = get_magnitude<isq::speed>(v);
+    CHECK(speed.numerical_value_ref_in(km / h) == 7);
   }
 
   SECTION("multiply by scalar value")
@@ -109,14 +118,14 @@ TEST_CASE("vector quantity", "[la]")
 
     SECTION("integral")
     {
-      SECTION("scalar on LHS") { CHECK((2 * v).numerical_value() == vector<int>{2, 4, 6}); }
-      SECTION("scalar on RHS") { CHECK((v * 2).numerical_value() == vector<int>{2, 4, 6}); }
+      SECTION("scalar on LHS") { CHECK((2 * v).numerical_value_ == vector<int>{2, 4, 6}); }
+      SECTION("scalar on RHS") { CHECK((v * 2).numerical_value_ == vector<int>{2, 4, 6}); }
     }
 
     SECTION("floating-point")
     {
-      SECTION("scalar on LHS") { CHECK((0.5 * v).numerical_value() == vector<double>{0.5, 1., 1.5}); }
-      SECTION("scalar on RHS") { CHECK((v * 0.5).numerical_value() == vector<double>{0.5, 1., 1.5}); }
+      SECTION("scalar on LHS") { CHECK((0.5 * v).numerical_value_ == vector<double>{0.5, 1., 1.5}); }
+      SECTION("scalar on RHS") { CHECK((v * 0.5).numerical_value_ == vector<double>{0.5, 1., 1.5}); }
     }
   }
 
@@ -124,8 +133,8 @@ TEST_CASE("vector quantity", "[la]")
   {
     const auto v = vector<int>{2, 4, 6} * isq::position_vector[m];
 
-    SECTION("integral") { CHECK((v / 2).numerical_value() == vector<int>{1, 2, 3}); }
-    SECTION("floating-point") { CHECK((v / 0.5).numerical_value() == vector<double>{4., 8., 12.}); }
+    SECTION("integral") { CHECK((v / 2).numerical_value_ == vector<int>{1, 2, 3}); }
+    SECTION("floating-point") { CHECK((v / 0.5).numerical_value_ == vector<double>{4., 8., 12.}); }
   }
 
   SECTION("add")
@@ -135,12 +144,12 @@ TEST_CASE("vector quantity", "[la]")
     SECTION("same unit")
     {
       const auto u = vector<int>{3, 2, 1} * isq::position_vector[m];
-      CHECK((v + u).numerical_value() == vector<int>{4, 4, 4});
+      CHECK((v + u).numerical_value_ == vector<int>{4, 4, 4});
     }
     SECTION("different units")
     {
       const auto u = vector<int>{3, 2, 1} * isq::position_vector[km];
-      CHECK((v + u).numerical_value() == vector<int>{3001, 2002, 1003});
+      CHECK((v + u).numerical_value_ == vector<int>{3001, 2002, 1003});
     }
   }
 
@@ -151,12 +160,12 @@ TEST_CASE("vector quantity", "[la]")
     SECTION("same unit")
     {
       const auto u = vector<int>{3, 2, 1} * isq::position_vector[m];
-      CHECK((v - u).numerical_value() == vector<int>{-2, 0, 2});
+      CHECK((v - u).numerical_value_ == vector<int>{-2, 0, 2});
     }
     SECTION("different units")
     {
       const auto u = vector<int>{3, 2, 1} * isq::position_vector[km];
-      CHECK((v - u).numerical_value() == vector<int>{-2999, -1998, -997});
+      CHECK((v - u).numerical_value_ == vector<int>{-2999, -1998, -997});
     }
   }
 
@@ -170,18 +179,18 @@ TEST_CASE("vector quantity", "[la]")
 
       SECTION("derived_quantity_spec")
       {
-        SECTION("scalar on LHS") { CHECK((mass * v).numerical_value() == vector<int>{2, 4, 6}); }
-        SECTION("scalar on RHS") { CHECK((v * mass).numerical_value() == vector<int>{2, 4, 6}); }
+        SECTION("scalar on LHS") { CHECK((mass * v).numerical_value_ == vector<int>{2, 4, 6}); }
+        SECTION("scalar on RHS") { CHECK((v * mass).numerical_value_ == vector<int>{2, 4, 6}); }
       }
       SECTION("quantity_cast to momentum")
       {
         SECTION("scalar on LHS")
         {
-          CHECK(quantity_cast<isq::momentum>(mass * v).numerical_value() == vector<int>{2, 4, 6});
+          CHECK(quantity_cast<isq::momentum>(mass * v).numerical_value_ == vector<int>{2, 4, 6});
         }
         SECTION("scalar on RHS")
         {
-          CHECK(quantity_cast<isq::momentum>(v * mass).numerical_value() == vector<int>{2, 4, 6});
+          CHECK(quantity_cast<isq::momentum>(v * mass).numerical_value_ == vector<int>{2, 4, 6});
         }
       }
       SECTION("quantity of momentum")
@@ -189,12 +198,12 @@ TEST_CASE("vector quantity", "[la]")
         SECTION("scalar on LHS")
         {
           const quantity<isq::momentum[N * s], vector<int>> momentum = mass * v;
-          CHECK(momentum.numerical_value() == vector<int>{2, 4, 6});
+          CHECK(momentum.numerical_value_ == vector<int>{2, 4, 6});
         }
         SECTION("scalar on RHS")
         {
           const quantity<isq::momentum[N * s], vector<int>> momentum = v * mass;
-          CHECK(momentum.numerical_value() == vector<int>{2, 4, 6});
+          CHECK(momentum.numerical_value_ == vector<int>{2, 4, 6});
         }
       }
     }
@@ -205,18 +214,18 @@ TEST_CASE("vector quantity", "[la]")
 
       SECTION("derived_quantity_spec")
       {
-        SECTION("scalar on LHS") { CHECK((mass * v).numerical_value() == vector<double>{0.5, 1., 1.5}); }
-        SECTION("scalar on RHS") { CHECK((v * mass).numerical_value() == vector<double>{0.5, 1., 1.5}); }
+        SECTION("scalar on LHS") { CHECK((mass * v).numerical_value_ == vector<double>{0.5, 1., 1.5}); }
+        SECTION("scalar on RHS") { CHECK((v * mass).numerical_value_ == vector<double>{0.5, 1., 1.5}); }
       }
       SECTION("quantity_cast to momentum")
       {
         SECTION("scalar on LHS")
         {
-          CHECK(quantity_cast<isq::momentum>(mass * v).numerical_value() == vector<double>{0.5, 1., 1.5});
+          CHECK(quantity_cast<isq::momentum>(mass * v).numerical_value_ == vector<double>{0.5, 1., 1.5});
         }
         SECTION("scalar on RHS")
         {
-          CHECK(quantity_cast<isq::momentum>(v * mass).numerical_value() == vector<double>{0.5, 1., 1.5});
+          CHECK(quantity_cast<isq::momentum>(v * mass).numerical_value_ == vector<double>{0.5, 1., 1.5});
         }
       }
       SECTION("quantity of momentum")
@@ -224,12 +233,12 @@ TEST_CASE("vector quantity", "[la]")
         SECTION("scalar on LHS")
         {
           const quantity<isq::momentum[N * s], vector<double>> momentum = mass * v;
-          CHECK(momentum.numerical_value() == vector<double>{0.5, 1., 1.5});
+          CHECK(momentum.numerical_value_ == vector<double>{0.5, 1., 1.5});
         }
         SECTION("scalar on RHS")
         {
           const quantity<isq::momentum[N * s], vector<double>> momentum = v * mass;
-          CHECK(momentum.numerical_value() == vector<double>{0.5, 1., 1.5});
+          CHECK(momentum.numerical_value_ == vector<double>{0.5, 1., 1.5});
         }
       }
     }
@@ -243,15 +252,15 @@ TEST_CASE("vector quantity", "[la]")
     {
       const auto dur = 2 * isq::duration[h];
 
-      SECTION("derived_quantity_spec") { CHECK((pos / dur).numerical_value() == vector<int>{15, 10, 5}); }
+      SECTION("derived_quantity_spec") { CHECK((pos / dur).numerical_value_ == vector<int>{15, 10, 5}); }
       SECTION("quantity_cast to velocity")
       {
-        CHECK(quantity_cast<isq::velocity>(pos / dur).numerical_value() == vector<int>{15, 10, 5});
+        CHECK(quantity_cast<isq::velocity>(pos / dur).numerical_value_ == vector<int>{15, 10, 5});
       }
       SECTION("quantity of velocity")
       {
         const quantity<isq::velocity[km / h], vector<int>> v = pos / dur;
-        CHECK(v.numerical_value() == vector<int>{15, 10, 5});
+        CHECK(v.numerical_value_ == vector<int>{15, 10, 5});
       }
     }
 
@@ -259,15 +268,15 @@ TEST_CASE("vector quantity", "[la]")
     {
       const auto dur = 0.5 * isq::duration[h];
 
-      SECTION("derived_quantity_spec") { CHECK((pos / dur).numerical_value() == vector<double>{60, 40, 20}); }
+      SECTION("derived_quantity_spec") { CHECK((pos / dur).numerical_value_ == vector<double>{60, 40, 20}); }
       SECTION("quantity_cast to velocity")
       {
-        CHECK(quantity_cast<isq::velocity>(pos / dur).numerical_value() == vector<double>{60, 40, 20});
+        CHECK(quantity_cast<isq::velocity>(pos / dur).numerical_value_ == vector<double>{60, 40, 20});
       }
       SECTION("quantity of velocity")
       {
         const quantity<isq::velocity[km / h], vector<double>> v = pos / dur;
-        CHECK(v.numerical_value() == vector<double>{60, 40, 20});
+        CHECK(v.numerical_value_ == vector<double>{60, 40, 20});
       }
     }
   }
@@ -303,8 +312,8 @@ TEST_CASE("vector of quantities", "[la]")
   SECTION("to scalar magnitude")
   {
     const vector<quantity<isq::velocity[km / h], int>> v = {2 * (km / h), 3 * (km / h), 6 * (km / h)};
-    const auto speed = get_magnitude(v).numerical_value() * isq::speed[v(0).unit];  // TODO can we do better here?
-    CHECK(speed.numerical_value() == 7);
+    const auto speed = get_magnitude<isq::speed>(v);
+    CHECK(speed.numerical_value_ref_in(km / h) == 7);
   }
 
   SECTION("multiply by scalar value")
diff --git a/test/unit_test/runtime/math_test.cpp b/test/unit_test/runtime/math_test.cpp
index 4e838273e..fdf5e1dc3 100644
--- a/test/unit_test/runtime/math_test.cpp
+++ b/test/unit_test/runtime/math_test.cpp
@@ -90,12 +90,12 @@ TEST_CASE("numeric_limits functions", "[limits]")
 {
   SECTION("'epsilon' works as expected using default floating type")
   {
-    REQUIRE(epsilon<double>(isq::length[m]).numerical_value() ==
+    REQUIRE(epsilon<double>(isq::length[m]).numerical_value_ ==
             std::numeric_limits<decltype(1. * isq::length[m])::rep>::epsilon());
   }
   SECTION("'epsilon' works as expected using integers")
   {
-    REQUIRE(epsilon<int>(isq::length[m]).numerical_value() ==
+    REQUIRE(epsilon<int>(isq::length[m]).numerical_value_ ==
             std::numeric_limits<decltype(1 * isq::length[m])::rep>::epsilon());
   }
 }
diff --git a/test/unit_test/static/quantity_point_test.cpp b/test/unit_test/static/quantity_point_test.cpp
index e93ebf09f..feb633d65 100644
--- a/test/unit_test/static/quantity_point_test.cpp
+++ b/test/unit_test/static/quantity_point_test.cpp
@@ -241,16 +241,18 @@ static_assert(
 // static member functions
 ////////////////////////////
 
-static_assert(quantity_point<isq::height[m], mean_sea_level>::zero().quantity_from_origin().numerical_value() == 0);
-static_assert(quantity_point<isq::height[m], mean_sea_level>::min().quantity_from_origin().numerical_value() ==
+static_assert(quantity_point<isq::height[m], mean_sea_level>::zero().quantity_from_origin_.numerical_value_ == 0);
+static_assert(quantity_point<isq::height[m], mean_sea_level>::min().quantity_from_origin_.numerical_value_ ==
               std::numeric_limits<double>::lowest());
-static_assert(quantity_point<isq::height[m], mean_sea_level>::max().quantity_from_origin().numerical_value() ==
+static_assert(quantity_point<isq::height[m], mean_sea_level>::max().quantity_from_origin_.numerical_value_ ==
               std::numeric_limits<double>::max());
 
-static_assert(quantity_point<isq::height[m], ground_level, int>::zero().quantity_from_origin().numerical_value() == 0);
-static_assert(quantity_point<isq::height[m], ground_level, int>::min().quantity_from_origin().numerical_value() ==
+static_assert(
+
+  quantity_point<isq::height[m], ground_level, int>::zero().quantity_from_origin_.numerical_value_ == 0);
+static_assert(quantity_point<isq::height[m], ground_level, int>::min().quantity_from_origin_.numerical_value_ ==
               std::numeric_limits<int>::lowest());
-static_assert(quantity_point<isq::height[m], ground_level, int>::max().quantity_from_origin().numerical_value() ==
+static_assert(quantity_point<isq::height[m], ground_level, int>::max().quantity_from_origin_.numerical_value_ ==
               std::numeric_limits<int>::max());
 
 
@@ -537,36 +539,36 @@ static_assert(
 // obtaining a relative quantity
 //////////////////////////////////
 
-static_assert((mean_sea_level + 42 * m).quantity_from_origin() == 42 * m);
-static_assert((mean_sea_level + isq::height(42 * m)).quantity_from_origin() == 42 * m);
+static_assert((mean_sea_level + 42 * m).quantity_from_origin_ == 42 * m);
+static_assert((mean_sea_level + isq::height(42 * m)).quantity_from_origin_ == 42 * m);
 
-static_assert((zero + 1 * one).quantity_from_origin() == 1 * one);
-static_assert((zero + dimensionless(1 * one)).quantity_from_origin() == 1 * one);
+static_assert((zero + 1 * one).quantity_from_origin_ == 1 * one);
+static_assert((zero + dimensionless(1 * one)).quantity_from_origin_ == 1 * one);
 
-static_assert((mean_sea_level + 42 * m).quantity_from_origin() == 42 * m);
-static_assert((ground_level + 42 * m).quantity_from_origin() == 42 * m);
-static_assert((tower_peak + 42 * m).quantity_from_origin() == 42 * m);
+static_assert((mean_sea_level + 42 * m).quantity_from_origin_ == 42 * m);
+static_assert((ground_level + 42 * m).quantity_from_origin_ == 42 * m);
+static_assert((tower_peak + 42 * m).quantity_from_origin_ == 42 * m);
 
-static_assert(quantity_point<isq::height[m], mean_sea_level>(ground_level + 42 * m).quantity_from_origin() == 84 * m);
-static_assert(quantity_point<isq::height[m], mean_sea_level>(tower_peak + 42 * m).quantity_from_origin() == 126 * m);
+static_assert(quantity_point<isq::height[m], mean_sea_level>(ground_level + 42 * m).quantity_from_origin_ == 84 * m);
+static_assert(quantity_point<isq::height[m], mean_sea_level>(tower_peak + 42 * m).quantity_from_origin_ == 126 * m);
 
-static_assert(quantity_point<isq::height[m], ground_level>(mean_sea_level + 84 * m).quantity_from_origin() == 42 * m);
-static_assert(quantity_point<isq::height[m], ground_level>(tower_peak + 42 * m).quantity_from_origin() == 84 * m);
+static_assert(quantity_point<isq::height[m], ground_level>(mean_sea_level + 84 * m).quantity_from_origin_ == 42 * m);
+static_assert(quantity_point<isq::height[m], ground_level>(tower_peak + 42 * m).quantity_from_origin_ == 84 * m);
 
-static_assert(quantity_point<isq::height[m], tower_peak>(mean_sea_level + 42 * m).quantity_from_origin() == -42 * m);
-static_assert(quantity_point<isq::height[m], tower_peak>(ground_level + 84 * m).quantity_from_origin() == 42 * m);
+static_assert(quantity_point<isq::height[m], tower_peak>(mean_sea_level + 42 * m).quantity_from_origin_ == -42 * m);
+static_assert(quantity_point<isq::height[m], tower_peak>(ground_level + 84 * m).quantity_from_origin_ == 42 * m);
 
-static_assert((mean_sea_level + 42 * m).point_for(mean_sea_level).quantity_from_origin() == 42 * m);
-static_assert((ground_level + 42 * m).point_for(mean_sea_level).quantity_from_origin() == 84 * m);
-static_assert((tower_peak + 42 * m).point_for(mean_sea_level).quantity_from_origin() == 126 * m);
+static_assert((mean_sea_level + 42 * m).point_for(mean_sea_level).quantity_from_origin_ == 42 * m);
+static_assert((ground_level + 42 * m).point_for(mean_sea_level).quantity_from_origin_ == 84 * m);
+static_assert((tower_peak + 42 * m).point_for(mean_sea_level).quantity_from_origin_ == 126 * m);
 
-static_assert((ground_level + 84 * m).point_for(ground_level).quantity_from_origin() == 84 * m);
-static_assert((mean_sea_level + 84 * m).point_for(ground_level).quantity_from_origin() == 42 * m);
-static_assert((tower_peak + 42 * m).point_for(ground_level).quantity_from_origin() == 84 * m);
+static_assert((ground_level + 84 * m).point_for(ground_level).quantity_from_origin_ == 84 * m);
+static_assert((mean_sea_level + 84 * m).point_for(ground_level).quantity_from_origin_ == 42 * m);
+static_assert((tower_peak + 42 * m).point_for(ground_level).quantity_from_origin_ == 84 * m);
 
-static_assert((tower_peak + 42 * m).point_for(tower_peak).quantity_from_origin() == 42 * m);
-static_assert((mean_sea_level + 42 * m).point_for(tower_peak).quantity_from_origin() == -42 * m);
-static_assert((ground_level + 84 * m).point_for(tower_peak).quantity_from_origin() == 42 * m);
+static_assert((tower_peak + 42 * m).point_for(tower_peak).quantity_from_origin_ == 42 * m);
+static_assert((mean_sea_level + 42 * m).point_for(tower_peak).quantity_from_origin_ == -42 * m);
+static_assert((ground_level + 84 * m).point_for(tower_peak).quantity_from_origin_ == 42 * m);
 
 static_assert(is_of_type<(ground_level + isq::height(short(42) * m)).point_for(mean_sea_level),
                          quantity_point<isq::height[m], mean_sea_level, int>>);
@@ -576,15 +578,15 @@ static_assert(is_of_type<(ground_level + isq::height(short(42) * m)).point_for(m
 // converting to a different unit
 ///////////////////////////////////
 
-static_assert((mean_sea_level + 2. * km).in(km).quantity_from_origin().numerical_value() == 2.);
-static_assert((mean_sea_level + 2. * km).in(m).quantity_from_origin().numerical_value() == 2000.);
-static_assert((mean_sea_level + 2000. * m).in(km).quantity_from_origin().numerical_value() == 2.);
-static_assert((ground_level + 2. * km).in(km).quantity_from_origin().numerical_value() == 2.);
-static_assert((ground_level + 2. * km).in(m).quantity_from_origin().numerical_value() == 2000.);
-static_assert((ground_level + 2000. * m).in(km).quantity_from_origin().numerical_value() == 2.);
-static_assert((tower_peak + 2. * km).in(km).quantity_from_origin().numerical_value() == 2.);
-static_assert((tower_peak + 2. * km).in(m).quantity_from_origin().numerical_value() == 2000.);
-static_assert((tower_peak + 2000. * m).in(km).quantity_from_origin().numerical_value() == 2.);
+static_assert((mean_sea_level + 2. * km).in(km).quantity_from_origin_.numerical_value_ == 2.);
+static_assert((mean_sea_level + 2. * km).in(m).quantity_from_origin_.numerical_value_ == 2000.);
+static_assert((mean_sea_level + 2000. * m).in(km).quantity_from_origin_.numerical_value_ == 2.);
+static_assert((ground_level + 2. * km).in(km).quantity_from_origin_.numerical_value_ == 2.);
+static_assert((ground_level + 2. * km).in(m).quantity_from_origin_.numerical_value_ == 2000.);
+static_assert((ground_level + 2000. * m).in(km).quantity_from_origin_.numerical_value_ == 2.);
+static_assert((tower_peak + 2. * km).in(km).quantity_from_origin_.numerical_value_ == 2.);
+static_assert((tower_peak + 2. * km).in(m).quantity_from_origin_.numerical_value_ == 2000.);
+static_assert((tower_peak + 2000. * m).in(km).quantity_from_origin_.numerical_value_ == 2.);
 
 #if MP_UNITS_COMP_GCC != 10 || MP_UNITS_COMP_GCC_MINOR > 2
 template<template<auto, auto, typename> typename QP>
@@ -616,18 +618,18 @@ static_assert(([]() {
                 quantity_point l1{mean_sea_level + 1 * m}, l2{mean_sea_level + 2 * m};
                 return l2 = l1;
               }())
-                .quantity_from_origin() == 1 * m);
+                .quantity_from_origin_ == 1 * m);
 static_assert(([]() {
                 const quantity_point l1{mean_sea_level + 1 * m};
                 quantity_point l2{mean_sea_level + 2 * m};
                 return l2 = l1;
               }())
-                .quantity_from_origin() == 1 * m);
+                .quantity_from_origin_ == 1 * m);
 static_assert(([]() {
                 quantity_point l1{mean_sea_level + 1 * m}, l2{mean_sea_level + 2 * m};
                 return l2 = std::move(l1);
               }())
-                .quantity_from_origin() == 1 * m);
+                .quantity_from_origin_ == 1 * m);
 
 
 ////////////////////
@@ -657,14 +659,14 @@ static_assert([](auto v) {
 ////////////////////////
 
 // same type
-static_assert((mean_sea_level + 1 * m += 1 * m).quantity_from_origin().numerical_value() == 2);
-static_assert((mean_sea_level + 2 * m -= 1 * m).quantity_from_origin().numerical_value() == 1);
+static_assert((mean_sea_level + 1 * m += 1 * m).quantity_from_origin_.numerical_value_ == 2);
+static_assert((mean_sea_level + 2 * m -= 1 * m).quantity_from_origin_.numerical_value_ == 1);
 
 // different types
-static_assert((mean_sea_level + 2.5 * m += 3 * m).quantity_from_origin().numerical_value() == 5.5);
-static_assert((mean_sea_level + 123 * m += 1 * km).quantity_from_origin().numerical_value() == 1123);
-static_assert((mean_sea_level + 5.5 * m -= 3 * m).quantity_from_origin().numerical_value() == 2.5);
-static_assert((mean_sea_level + 1123 * m -= 1 * km).quantity_from_origin().numerical_value() == 123);
+static_assert((mean_sea_level + 2.5 * m += 3 * m).quantity_from_origin_.numerical_value_ == 5.5);
+static_assert((mean_sea_level + 123 * m += 1 * km).quantity_from_origin_.numerical_value_ == 1123);
+static_assert((mean_sea_level + 5.5 * m -= 3 * m).quantity_from_origin_.numerical_value_ == 2.5);
+static_assert((mean_sea_level + 1123 * m -= 1 * km).quantity_from_origin_.numerical_value_ == 123);
 
 
 template<template<auto, auto, typename> typename QP>
@@ -935,30 +937,25 @@ static_assert(is_of_type<(1 * m + tower_peak) - (1 * m + other_ground_level), qu
 
 
 // check for integral types promotion
+static_assert(is_same_v<decltype(((mean_sea_level + std::uint8_t(0) * m) + std::uint8_t(0) * m)
+                                   .quantity_from_origin_.numerical_value_),
+                        int>);
+static_assert(is_same_v<decltype((std::uint8_t(0) * m + (mean_sea_level + std::uint8_t(0) * m))
+                                   .quantity_from_origin_.numerical_value_),
+                        int>);
+static_assert(is_same_v<decltype(((mean_sea_level + std::uint8_t(0) * m) - std::uint8_t(0) * m)
+                                   .quantity_from_origin_.numerical_value_),
+                        int>);
 static_assert(
   is_same_v<
-    decltype(((mean_sea_level + std::uint8_t(0) * m) + std::uint8_t(0) * m).quantity_from_origin().numerical_value()),
-    int&&>);
-static_assert(
-  is_same_v<
-    decltype((std::uint8_t(0) * m + (mean_sea_level + std::uint8_t(0) * m)).quantity_from_origin().numerical_value()),
-    int&&>);
-static_assert(
-  is_same_v<
-    decltype(((mean_sea_level + std::uint8_t(0) * m) - std::uint8_t(0) * m).quantity_from_origin().numerical_value()),
-    int&&>);
-static_assert(is_same_v<decltype(((mean_sea_level + std::uint8_t(0) * m) - (mean_sea_level + std::uint8_t(0) * m))
-                                   .numerical_value()),
-                        int&&>);
-static_assert(
-  ((mean_sea_level + std::uint8_t(128) * m) + std::uint8_t(128) * m).quantity_from_origin().numerical_value() ==
-  std::uint8_t(128) + std::uint8_t(128));
-static_assert(
-  (std::uint8_t(128) * m + (mean_sea_level + std::uint8_t(128) * m)).quantity_from_origin().numerical_value() ==
-  std::uint8_t(128) + std::uint8_t(128));
-static_assert(((mean_sea_level + std::uint8_t(0) * m) - std::uint8_t(1) * m).quantity_from_origin().numerical_value() ==
+    decltype(((mean_sea_level + std::uint8_t(0) * m) - (mean_sea_level + std::uint8_t(0) * m)).numerical_value_), int>);
+static_assert(((mean_sea_level + std::uint8_t(128) * m) + std::uint8_t(128) * m)
+                .quantity_from_origin_.numerical_value_ == std::uint8_t(128) + std::uint8_t(128));
+static_assert((std::uint8_t(128) * m + (mean_sea_level + std::uint8_t(128) * m))
+                .quantity_from_origin_.numerical_value_ == std::uint8_t(128) + std::uint8_t(128));
+static_assert(((mean_sea_level + std::uint8_t(0) * m) - std::uint8_t(1) * m).quantity_from_origin_.numerical_value_ ==
               std::uint8_t(0) - std::uint8_t(1));
-static_assert(((mean_sea_level + std::uint8_t(0) * m) - (mean_sea_level + std::uint8_t(1) * m)).numerical_value() ==
+static_assert(((mean_sea_level + std::uint8_t(0) * m) - (mean_sea_level + std::uint8_t(1) * m)).numerical_value_ ==
               std::uint8_t(0) - std::uint8_t(1));
 
 // different representation types
@@ -1013,39 +1010,39 @@ static_assert(is_of_type<(mean_sea_level + 1 * km) - (mean_sea_level + 1. * m),
 static_assert(is_of_type<(mean_sea_level + 1. * km) - (mean_sea_level + 1. * m), quantity<si::metre, double>>);
 
 
-static_assert(((mean_sea_level + 1 * m) + 1 * m).quantity_from_origin().numerical_value() == 2);
-static_assert((1 * m + (mean_sea_level + 1 * m)).quantity_from_origin().numerical_value() == 2);
-static_assert(((mean_sea_level + 1 * m) + 1 * km).quantity_from_origin().numerical_value() == 1001);
-static_assert((1 * m + (mean_sea_level + 1 * km)).quantity_from_origin().numerical_value() == 1001);
-static_assert(((mean_sea_level + 1 * km) + 1 * m).quantity_from_origin().numerical_value() == 1001);
-static_assert((1 * km + (mean_sea_level + 1 * m)).quantity_from_origin().numerical_value() == 1001);
-static_assert(((mean_sea_level + 2 * m) - 1 * m).quantity_from_origin().numerical_value() == 1);
-static_assert(((mean_sea_level + 1 * km) - 1 * m).quantity_from_origin().numerical_value() == 999);
-
-static_assert(((mean_sea_level + 1.5 * m) + 1 * m).quantity_from_origin().numerical_value() == 2.5);
-static_assert((1.5 * m + (mean_sea_level + 1 * m)).quantity_from_origin().numerical_value() == 2.5);
-static_assert(((mean_sea_level + 1.5 * m) + 1 * km).quantity_from_origin().numerical_value() == 1001.5);
-static_assert((1.5 * m + (mean_sea_level + 1 * km)).quantity_from_origin().numerical_value() == 1001.5);
-static_assert(((mean_sea_level + 1.5 * km) + 1 * m).quantity_from_origin().numerical_value() == 1501);
-static_assert((1.5 * km + (mean_sea_level + 1 * m)).quantity_from_origin().numerical_value() == 1501);
-static_assert(((mean_sea_level + 2.5 * m) - 1 * m).quantity_from_origin().numerical_value() == 1.5);
-static_assert(((mean_sea_level + 1.5 * km) - 1 * m).quantity_from_origin().numerical_value() == 1499);
-
-static_assert(((mean_sea_level + 1 * m) + 1.5 * m).quantity_from_origin().numerical_value() == 2.5);
-static_assert((1 * m + (mean_sea_level + 1.5 * m)).quantity_from_origin().numerical_value() == 2.5);
-static_assert(((mean_sea_level + 1 * m) + 1.5 * km).quantity_from_origin().numerical_value() == 1501);
-static_assert((1 * m + (mean_sea_level + 1.5 * km)).quantity_from_origin().numerical_value() == 1501);
-static_assert(((mean_sea_level + 1 * km) + 1.5 * m).quantity_from_origin().numerical_value() == 1001.5);
-static_assert((1 * km + (mean_sea_level + 1.5 * m)).quantity_from_origin().numerical_value() == 1001.5);
-static_assert(((mean_sea_level + 2 * m) - 1.5 * m).quantity_from_origin().numerical_value() == 0.5);
-static_assert(((mean_sea_level + 1 * km) - 1.5 * m).quantity_from_origin().numerical_value() == 998.5);
-
-static_assert(((mean_sea_level + 2 * m) - (mean_sea_level + 1 * m)).numerical_value() == 1);
-static_assert(((mean_sea_level + 1 * km) - (mean_sea_level + 1 * m)).numerical_value() == 999);
-static_assert(((mean_sea_level + 2.5 * m) - (mean_sea_level + 1 * m)).numerical_value() == 1.5);
-static_assert(((mean_sea_level + 1.5 * km) - (mean_sea_level + 1 * m)).numerical_value() == 1499);
-static_assert(((mean_sea_level + 2 * m) - (mean_sea_level + 1.5 * m)).numerical_value() == 0.5);
-static_assert(((mean_sea_level + 1 * km) - (mean_sea_level + 1.5 * m)).numerical_value() == 998.5);
+static_assert(((mean_sea_level + 1 * m) + 1 * m).quantity_from_origin_.numerical_value_ == 2);
+static_assert((1 * m + (mean_sea_level + 1 * m)).quantity_from_origin_.numerical_value_ == 2);
+static_assert(((mean_sea_level + 1 * m) + 1 * km).quantity_from_origin_.numerical_value_ == 1001);
+static_assert((1 * m + (mean_sea_level + 1 * km)).quantity_from_origin_.numerical_value_ == 1001);
+static_assert(((mean_sea_level + 1 * km) + 1 * m).quantity_from_origin_.numerical_value_ == 1001);
+static_assert((1 * km + (mean_sea_level + 1 * m)).quantity_from_origin_.numerical_value_ == 1001);
+static_assert(((mean_sea_level + 2 * m) - 1 * m).quantity_from_origin_.numerical_value_ == 1);
+static_assert(((mean_sea_level + 1 * km) - 1 * m).quantity_from_origin_.numerical_value_ == 999);
+
+static_assert(((mean_sea_level + 1.5 * m) + 1 * m).quantity_from_origin_.numerical_value_ == 2.5);
+static_assert((1.5 * m + (mean_sea_level + 1 * m)).quantity_from_origin_.numerical_value_ == 2.5);
+static_assert(((mean_sea_level + 1.5 * m) + 1 * km).quantity_from_origin_.numerical_value_ == 1001.5);
+static_assert((1.5 * m + (mean_sea_level + 1 * km)).quantity_from_origin_.numerical_value_ == 1001.5);
+static_assert(((mean_sea_level + 1.5 * km) + 1 * m).quantity_from_origin_.numerical_value_ == 1501);
+static_assert((1.5 * km + (mean_sea_level + 1 * m)).quantity_from_origin_.numerical_value_ == 1501);
+static_assert(((mean_sea_level + 2.5 * m) - 1 * m).quantity_from_origin_.numerical_value_ == 1.5);
+static_assert(((mean_sea_level + 1.5 * km) - 1 * m).quantity_from_origin_.numerical_value_ == 1499);
+
+static_assert(((mean_sea_level + 1 * m) + 1.5 * m).quantity_from_origin_.numerical_value_ == 2.5);
+static_assert((1 * m + (mean_sea_level + 1.5 * m)).quantity_from_origin_.numerical_value_ == 2.5);
+static_assert(((mean_sea_level + 1 * m) + 1.5 * km).quantity_from_origin_.numerical_value_ == 1501);
+static_assert((1 * m + (mean_sea_level + 1.5 * km)).quantity_from_origin_.numerical_value_ == 1501);
+static_assert(((mean_sea_level + 1 * km) + 1.5 * m).quantity_from_origin_.numerical_value_ == 1001.5);
+static_assert((1 * km + (mean_sea_level + 1.5 * m)).quantity_from_origin_.numerical_value_ == 1001.5);
+static_assert(((mean_sea_level + 2 * m) - 1.5 * m).quantity_from_origin_.numerical_value_ == 0.5);
+static_assert(((mean_sea_level + 1 * km) - 1.5 * m).quantity_from_origin_.numerical_value_ == 998.5);
+
+static_assert(((mean_sea_level + 2 * m) - (mean_sea_level + 1 * m)).numerical_value_ == 1);
+static_assert(((mean_sea_level + 1 * km) - (mean_sea_level + 1 * m)).numerical_value_ == 999);
+static_assert(((mean_sea_level + 2.5 * m) - (mean_sea_level + 1 * m)).numerical_value_ == 1.5);
+static_assert(((mean_sea_level + 1.5 * km) - (mean_sea_level + 1 * m)).numerical_value_ == 1499);
+static_assert(((mean_sea_level + 2 * m) - (mean_sea_level + 1.5 * m)).numerical_value_ == 0.5);
+static_assert(((mean_sea_level + 1 * km) - (mean_sea_level + 1.5 * m)).numerical_value_ == 998.5);
 
 static_assert((mean_sea_level + 42 * m) - (ground_level + 42 * m) == -42 * m);
 static_assert((ground_level + 42 * m) - (mean_sea_level + 42 * m) == 42 * m);
@@ -1058,15 +1055,15 @@ static_assert((ground_level + 42 * m) - (other_ground_level + 42 * m) == -81 * m
 static_assert((other_ground_level + 42 * m) - (tower_peak + 42 * m) == 39 * m);
 static_assert((tower_peak + 42 * m) - (other_ground_level + 42 * m) == -39 * m);
 
-static_assert((mean_sea_level + 42 * m).quantity_from_origin() == 42 * m);
-static_assert((42 * m + mean_sea_level).quantity_from_origin() == 42 * m);
-static_assert((mean_sea_level - 42 * m).quantity_from_origin() == -42 * m);
-static_assert((ground_level + 42 * m).quantity_from_origin() == 42 * m);
-static_assert((42 * m + ground_level).quantity_from_origin() == 42 * m);
-static_assert((ground_level - 42 * m).quantity_from_origin() == -42 * m);
-static_assert((tower_peak + 42 * m).quantity_from_origin() == 42 * m);
-static_assert((42 * m + tower_peak).quantity_from_origin() == 42 * m);
-static_assert((tower_peak - 42 * m).quantity_from_origin() == -42 * m);
+static_assert((mean_sea_level + 42 * m).quantity_from_origin_ == 42 * m);
+static_assert((42 * m + mean_sea_level).quantity_from_origin_ == 42 * m);
+static_assert((mean_sea_level - 42 * m).quantity_from_origin_ == -42 * m);
+static_assert((ground_level + 42 * m).quantity_from_origin_ == 42 * m);
+static_assert((42 * m + ground_level).quantity_from_origin_ == 42 * m);
+static_assert((ground_level - 42 * m).quantity_from_origin_ == -42 * m);
+static_assert((tower_peak + 42 * m).quantity_from_origin_ == 42 * m);
+static_assert((42 * m + tower_peak).quantity_from_origin_ == 42 * m);
+static_assert((tower_peak - 42 * m).quantity_from_origin_ == -42 * m);
 
 static_assert((mean_sea_level + 42 * m) - ground_level == 0 * m);
 static_assert((ground_level + 42 * m) - mean_sea_level == 84 * m);
@@ -1079,6 +1076,17 @@ static_assert((ground_level + 42 * m) - other_ground_level == -39 * m);
 static_assert((other_ground_level + 42 * m) - tower_peak == 81 * m);
 static_assert((tower_peak + 42 * m) - other_ground_level == 3 * m);
 
+static_assert((mean_sea_level + 42 * m).quantity_from(ground_level) == 0 * m);
+static_assert((ground_level + 42 * m).quantity_from(mean_sea_level) == 84 * m);
+static_assert((tower_peak + 42 * m).quantity_from(ground_level) == 84 * m);
+static_assert((ground_level + 42 * m).quantity_from(tower_peak) == 0 * m);
+static_assert((tower_peak + 42 * m).quantity_from(mean_sea_level) == 126 * m);
+static_assert((mean_sea_level + 42 * m).quantity_from(tower_peak) == -42 * m);
+static_assert((other_ground_level + 42 * m).quantity_from(ground_level) == 123 * m);
+static_assert((ground_level + 42 * m).quantity_from(other_ground_level) == -39 * m);
+static_assert((other_ground_level + 42 * m).quantity_from(tower_peak) == 81 * m);
+static_assert((tower_peak + 42 * m).quantity_from(other_ground_level) == 3 * m);
+
 static_assert(mean_sea_level - (ground_level + 42 * m) == -84 * m);
 static_assert(ground_level - (mean_sea_level + 42 * m) == 0 * m);
 static_assert(tower_peak - (ground_level + 42 * m) == 0 * m);
@@ -1106,17 +1114,17 @@ inline constexpr struct zero_m_per_s : absolute_point_origin<kind_of<isq::speed>
 
 // commutativity and associativity
 static_assert(((zero_m_per_s + 10 * isq::height[m] / (2 * isq::time[s])) + 5 * isq::speed[m / s])
-                .quantity_from_origin() == 10 * isq::speed[m / s]);
+                .quantity_from_origin_ == 10 * isq::speed[m / s]);
 static_assert((10 * isq::height[m] / (2 * isq::time[s]) + (zero_m_per_s + 5 * isq::speed[m / s]))
-                .quantity_from_origin() == 10 * isq::speed[m / s]);
+                .quantity_from_origin_ == 10 * isq::speed[m / s]);
 static_assert(((zero_m_per_s + 5 * isq::speed[m / s]) + 10 * isq::height[m] / (2 * isq::time[s]))
-                .quantity_from_origin() == 10 * isq::speed[m / s]);
+                .quantity_from_origin_ == 10 * isq::speed[m / s]);
 static_assert((5 * isq::speed[m / s] + (zero_m_per_s + 10 * isq::height[m] / (2 * isq::time[s])))
-                .quantity_from_origin() == 10 * isq::speed[m / s]);
+                .quantity_from_origin_ == 10 * isq::speed[m / s]);
 static_assert(((zero_m_per_s + 10 * isq::height[m] / (2 * isq::time[s])) - 5 * isq::speed[m / s])
-                .quantity_from_origin() == 0 * isq::speed[m / s]);
+                .quantity_from_origin_ == 0 * isq::speed[m / s]);
 static_assert(((zero_m_per_s + 5 * isq::speed[m / s]) - 10 * isq::height[m] / (2 * isq::time[s]))
-                .quantity_from_origin() == 0 * isq::speed[m / s]);
+                .quantity_from_origin_ == 0 * isq::speed[m / s]);
 static_assert((zero_m_per_s + 10 * isq::height[m] / (2 * isq::time[s])) - (zero_m_per_s + 5 * isq::speed[m / s]) ==
               0 * isq::speed[m / s]);
 static_assert((zero_m_per_s + 5 * isq::speed[m / s]) - (zero_m_per_s + 10 * isq::height[m] / (2 * isq::time[s])) ==
@@ -1148,17 +1156,17 @@ static_assert(
 inline constexpr struct zero_Hz : absolute_point_origin<kind_of<isq::frequency>> {
 } zero_Hz;
 
-static_assert(((zero_Hz + 10 / (2 * isq::period_duration[s])) + 5 * isq::frequency[Hz]).quantity_from_origin() ==
+static_assert(((zero_Hz + 10 / (2 * isq::period_duration[s])) + 5 * isq::frequency[Hz]).quantity_from_origin_ ==
               10 * isq::frequency[Hz]);
-static_assert((10 / (2 * isq::period_duration[s]) + (zero_Hz + 5 * isq::frequency[Hz])).quantity_from_origin() ==
+static_assert((10 / (2 * isq::period_duration[s]) + (zero_Hz + 5 * isq::frequency[Hz])).quantity_from_origin_ ==
               10 * isq::frequency[Hz]);
-static_assert(((zero_Hz + 5 * isq::frequency[Hz]) + 10 / (2 * isq::period_duration[s])).quantity_from_origin() ==
+static_assert(((zero_Hz + 5 * isq::frequency[Hz]) + 10 / (2 * isq::period_duration[s])).quantity_from_origin_ ==
               10 * isq::frequency[Hz]);
-static_assert((5 * isq::frequency[Hz] + (zero_Hz + 10 / (2 * isq::period_duration[s]))).quantity_from_origin() ==
+static_assert((5 * isq::frequency[Hz] + (zero_Hz + 10 / (2 * isq::period_duration[s]))).quantity_from_origin_ ==
               10 * isq::frequency[Hz]);
-static_assert(((zero_Hz + 10 / (2 * isq::period_duration[s])) - 5 * isq::frequency[Hz]).quantity_from_origin() ==
+static_assert(((zero_Hz + 10 / (2 * isq::period_duration[s])) - 5 * isq::frequency[Hz]).quantity_from_origin_ ==
               0 * isq::frequency[Hz]);
-static_assert(((zero_Hz + 5 * isq::frequency[Hz]) - 10 / (2 * isq::period_duration[s])).quantity_from_origin() ==
+static_assert(((zero_Hz + 5 * isq::frequency[Hz]) - 10 / (2 * isq::period_duration[s])).quantity_from_origin_ ==
               0 * isq::frequency[Hz]);
 static_assert((zero_Hz + 10 / (2 * isq::period_duration[s])) - (zero_Hz + 5 * isq::frequency[Hz]) ==
               0 * isq::frequency[Hz]);
diff --git a/test/unit_test/static/quantity_test.cpp b/test/unit_test/static/quantity_test.cpp
index 0fdc13b3f..77091cb73 100644
--- a/test/unit_test/static/quantity_test.cpp
+++ b/test/unit_test/static/quantity_test.cpp
@@ -112,14 +112,14 @@ static_assert(is_same_v<quantity<isq::length[m], int>::rep, int>);
 // static member functions
 ////////////////////////////
 
-static_assert(quantity<isq::length[m], int>::zero().numerical_value() == 0);
-static_assert(quantity<isq::length[m], int>::one().numerical_value() == 1);
-static_assert(quantity<isq::length[m], int>::min().numerical_value() == std::numeric_limits<int>::lowest());
-static_assert(quantity<isq::length[m], int>::max().numerical_value() == std::numeric_limits<int>::max());
-static_assert(quantity<isq::length[m], double>::zero().numerical_value() == 0.0);
-static_assert(quantity<isq::length[m], double>::one().numerical_value() == 1.0);
-static_assert(quantity<isq::length[m], double>::min().numerical_value() == std::numeric_limits<double>::lowest());
-static_assert(quantity<isq::length[m], double>::max().numerical_value() == std::numeric_limits<double>::max());
+static_assert(quantity<isq::length[m], int>::zero().numerical_value_ == 0);
+static_assert(quantity<isq::length[m], int>::one().numerical_value_ == 1);
+static_assert(quantity<isq::length[m], int>::min().numerical_value_ == std::numeric_limits<int>::lowest());
+static_assert(quantity<isq::length[m], int>::max().numerical_value_ == std::numeric_limits<int>::max());
+static_assert(quantity<isq::length[m], double>::zero().numerical_value_ == 0.0);
+static_assert(quantity<isq::length[m], double>::one().numerical_value_ == 1.0);
+static_assert(quantity<isq::length[m], double>::min().numerical_value_ == std::numeric_limits<double>::lowest());
+static_assert(quantity<isq::length[m], double>::max().numerical_value_ == std::numeric_limits<double>::max());
 
 
 //////////////////////////////
@@ -190,10 +190,10 @@ static_assert(std::convertible_to<quantity<isq::length[m], int>, quantity<isq::l
 // obtaining a number
 ///////////////////////
 
-static_assert(quantity<isq::length[m], int>(123 * m).numerical_value() == 123);
-static_assert(quantity<isq::length[m], int>(2 * km).numerical_value() == 2000);
-static_assert(quantity<isq::length[km], int>(2 * km).numerical_value() == 2);
-static_assert(quantity<isq::length[km]>(1500 * m).numerical_value() == 1.5);
+static_assert(quantity<isq::length[m], int>(123 * m).numerical_value_ == 123);
+static_assert(quantity<isq::length[m], int>(2 * km).numerical_value_ == 2000);
+static_assert(quantity<isq::length[km], int>(2 * km).numerical_value_ == 2);
+static_assert(quantity<isq::length[km]>(1500 * m).numerical_value_ == 1.5);
 
 
 ///////////////////////////////////
@@ -204,11 +204,11 @@ static_assert(is_of_type<(2. * km).in(m), quantity<si::metre>>);
 static_assert(is_of_type<isq::length(2. * km).in(m), quantity<isq::length[m]>>);
 static_assert(is_of_type<isq::height(2. * km).in(m), quantity<isq::height[m]>>);
 
-static_assert(quantity<isq::length[km]>(2. * km).in(km).numerical_value() == 2.);
-static_assert(quantity<isq::length[km]>(2. * km).in(m).numerical_value() == 2000.);
-static_assert(quantity<isq::length[m]>(2000. * m).in(km).numerical_value() == 2.);
-static_assert(quantity<isq::length[km], int>(2 * km).in(km).numerical_value() == 2);
-static_assert(quantity<isq::length[km], int>(2 * km).in(m).numerical_value() == 2000);
+static_assert(quantity<isq::length[km]>(2. * km).in(km).numerical_value_ == 2.);
+static_assert(quantity<isq::length[km]>(2. * km).in(m).numerical_value_ == 2000.);
+static_assert(quantity<isq::length[m]>(2000. * m).in(km).numerical_value_ == 2.);
+static_assert(quantity<isq::length[km], int>(2 * km).in(km).numerical_value_ == 2);
+static_assert(quantity<isq::length[km], int>(2 * km).in(m).numerical_value_ == 2000);
 
 #if MP_UNITS_COMP_GCC != 10 || MP_UNITS_COMP_GCC_MINOR > 2
 template<template<auto, typename> typename Q>
@@ -308,28 +308,28 @@ static_assert([] {
   auto l1(1 * m), l2(2 * m);
   return l2 = l1;
 }()
-                .numerical_value() == 1);
+                .numerical_value_ == 1);
 static_assert([] {
   const auto l1(1 * m);
   auto l2(2 * m);
   return l2 = l1;
 }()
-                .numerical_value() == 1);
+                .numerical_value_ == 1);
 static_assert([]() {
   auto l1(1 * m), l2(2 * m);
   return l2 = std::move(l1);
 }()
-                .numerical_value() == 1);
+                .numerical_value_ == 1);
 
 
 ////////////////////
 // unary operators
 ////////////////////
 
-static_assert((+123 * m).numerical_value() == 123);
-static_assert((-123 * m).numerical_value() == -123);
-static_assert((+(-123 * m)).numerical_value() == -123);
-static_assert((-(-123 * m)).numerical_value() == 123);
+static_assert((+123 * m).numerical_value_ == 123);
+static_assert((-123 * m).numerical_value_ == -123);
+static_assert((+(-123 * m)).numerical_value_ == -123);
+static_assert((-(-123 * m)).numerical_value_ == 123);
 
 static_assert([](auto v) {
   auto vv = v++;
@@ -348,7 +348,7 @@ static_assert([](auto v) {
   return std::pair(v, vv);
 }(123 * m) == std::pair(122 * m, 122 * m));
 
-static_assert(is_same_v<decltype((+(short{0} * m)).numerical_value()), int&&>);
+static_assert(is_same_v<decltype((+(short{0} * m)).numerical_value_), int>);
 
 
 ////////////////////////
@@ -356,29 +356,31 @@ static_assert(is_same_v<decltype((+(short{0} * m)).numerical_value()), int&&>);
 ////////////////////////
 
 // same type
-static_assert((1 * m += 1 * m).numerical_value() == 2);
-static_assert((2 * m -= 1 * m).numerical_value() == 1);
-static_assert((1 * m *= 2).numerical_value() == 2);
-static_assert((2 * m /= 2).numerical_value() == 1);
-static_assert((1 * m *= 2 * one).numerical_value() == 2);
-static_assert((2 * m /= 2 * one).numerical_value() == 1);
-static_assert((7 * m %= 2 * m).numerical_value() == 1);
+static_assert((1 * m += 1 * m).numerical_value_ == 2);
+static_assert((2 * m -= 1 * m).numerical_value_ == 1);
+static_assert((1 * m *= 2).numerical_value_ == 2);
+static_assert((2 * m /= 2).numerical_value_ == 1);
+static_assert((1 * m *= 2 * one).numerical_value_ == 2);
+static_assert((2 * m /= 2 * one).numerical_value_ == 1);
+static_assert((7 * m %= 2 * m).numerical_value_ == 1);
 
 // different types
-static_assert((2.5 * m += 3 * m).numerical_value() == 5.5);
-static_assert((123 * m += 1 * km).numerical_value() == 1123);
-static_assert((5.5 * m -= 3 * m).numerical_value() == 2.5);
-static_assert((1123 * m -= 1 * km).numerical_value() == 123);
-static_assert((2.5 * m *= 3).numerical_value() == 7.5);
-static_assert((7.5 * m /= 3).numerical_value() == 2.5);
-static_assert((2.5 * m *= 3 * one).numerical_value() == 7.5);
-static_assert((7.5 * m /= 3 * one).numerical_value() == 2.5);
-static_assert((3500 * m %= 1 * km).numerical_value() == 500);
-
-// static_assert((std::uint8_t(255) * m %= 256 * m).numerical_value() != [] { std::uint8_t ui(255); return ui %= 256;
+static_assert((2.5 * m += 3 * m).numerical_value_ == 5.5);
+static_assert((123 * m += 1 * km).numerical_value_ == 1123);
+static_assert((5.5 * m -= 3 * m).numerical_value_ == 2.5);
+static_assert((1123 * m -= 1 * km).numerical_value_ == 123);
+static_assert((2.5 * m *= 3).numerical_value_ == 7.5);
+static_assert((7.5 * m /= 3).numerical_value_ == 2.5);
+static_assert((2.5 * m *= 3 * one).numerical_value_ == 7.5);
+static_assert((7.5 * m /= 3 * one).numerical_value_ == 2.5);
+static_assert((3500 * m %= 1 * km).numerical_value_ == 500);
+
+// static_assert((std::uint8_t(255) * m %= 256 * m).numerical_value_ == [] {
+//   std::uint8_t ui(255);
+//   return ui %= 256;
 // }());  // UB
 // TODO: Fix
-static_assert((std::uint8_t(255) * m %= 257 * m).numerical_value() != [] {
+static_assert((std::uint8_t(255) * m %= 257 * m).numerical_value_ != [] {
   std::uint8_t ui(255);
   return ui %= 257;
 }());
@@ -388,10 +390,10 @@ static_assert((std::uint8_t(255) * m %= 257 * m).numerical_value() != [] {
 #ifndef MP_UNITS_COMP_MSVC
 // next two lines trigger conversions warnings
 // (warning disabled in CMake for this file)
-static_assert((22 * m *= 33.33).numerical_value() == 733);
-static_assert((22 * m /= 3.33).numerical_value() == 6);
-static_assert((22 * m *= 33.33 * one).numerical_value() == 733);
-static_assert((22 * m /= 3.33 * one).numerical_value() == 6);
+static_assert((22 * m *= 33.33).numerical_value_ == 733);
+static_assert((22 * m /= 3.33).numerical_value_ == 6);
+static_assert((22 * m *= 33.33 * one).numerical_value_ == 733);
+static_assert((22 * m /= 3.33 * one).numerical_value_ == 6);
 #endif
 
 template<template<auto, typename> typename Q>
@@ -516,14 +518,14 @@ static_assert(is_of_type<1 * km % (300 * m), quantity<si::metre, int>>);
 static_assert(is_of_type<4 * one % (2 * one), quantity<one, int>>);
 
 // check for integral types promotion
-static_assert(is_same_v<decltype((std::uint8_t(0) * m + std::uint8_t(0) * m).numerical_value()), int&&>);
-static_assert(is_same_v<decltype((std::uint8_t(0) * m - std::uint8_t(0) * m).numerical_value()), int&&>);
-static_assert((std::uint8_t(128) * m + std::uint8_t(128) * m).numerical_value() ==
+static_assert(is_same_v<decltype(std::uint8_t(0) * m + std::uint8_t(0) * m)::rep, int>);
+static_assert(is_same_v<decltype(std::uint8_t(0) * m - std::uint8_t(0) * m)::rep, int>);
+static_assert((std::uint8_t(128) * m + std::uint8_t(128) * m).numerical_value_ ==
               std::uint8_t(128) + std::uint8_t(128));
-static_assert((std::uint8_t(0) * m - std::uint8_t(1) * m).numerical_value() == std::uint8_t(0) - std::uint8_t(1));
+static_assert((std::uint8_t(0) * m - std::uint8_t(1) * m).numerical_value_ == std::uint8_t(0) - std::uint8_t(1));
 
-static_assert(is_same_v<decltype(((std::uint8_t(0) * m) % (std::uint8_t(0) * m)).numerical_value()),
-                        decltype(std::uint8_t(0) % std::uint8_t(0))&&>);
+static_assert(
+  is_same_v<decltype((std::uint8_t(0) * m) % (std::uint8_t(0) * m))::rep, decltype(std::uint8_t(0) % std::uint8_t(0))>);
 
 // different representation types
 static_assert(is_of_type<1. * m + 1 * m, quantity<si::metre, double>>);
@@ -596,67 +598,67 @@ static_assert(is_of_type<1 * m / (1 * s), quantity<derived_unit<struct si::metre
 static_assert(is_of_type<1 * m / (1 * min), quantity<derived_unit<struct si::metre, per<struct si::minute>>{}, int>>);
 static_assert(is_of_type<1 * min / (1 * m), quantity<derived_unit<struct si::minute, per<struct si::metre>>{}, int>>);
 
-static_assert((1 * m + 1 * m).numerical_value() == 2);
-static_assert((1 * m + 1 * km).numerical_value() == 1001);
-static_assert((1 * km + 1 * m).numerical_value() == 1001);
-static_assert((2 * m - 1 * m).numerical_value() == 1);
-static_assert((1 * km - 1 * m).numerical_value() == 999);
-static_assert((2 * m * 2).numerical_value() == 4);
-static_assert((2 * m * (2 * one)).numerical_value() == 4);
-static_assert((2 * m * (2 * percent)).numerical_value() == 4);
-static_assert((3 * 3 * m).numerical_value() == 9);
-static_assert(((3 * one) * (3 * m)).numerical_value() == 9);
-static_assert(((3 * percent) * (3 * m)).numerical_value() == 9);
-static_assert((4 * m / 2).numerical_value() == 2);
-static_assert((4 * m / (2 * one)).numerical_value() == 2);
-static_assert((4 * m / (2 * percent)).numerical_value() == 2);
-static_assert((4 * km / (2 * m)).numerical_value() == 2);
-static_assert((4000 * m / (2 * m)).numerical_value() == 2000);
-
-static_assert((1.5 * m + 1 * m).numerical_value() == 2.5);
-static_assert((1.5 * m + 1 * km).numerical_value() == 1001.5);
-static_assert((1.5 * km + 1 * m).numerical_value() == 1501);
-static_assert((2.5 * m - 1 * m).numerical_value() == 1.5);
-static_assert((1.5 * km - 1 * m).numerical_value() == 1499);
-static_assert((2.5 * m * 2).numerical_value() == 5);
-static_assert((2.5 * m * (2 * one)).numerical_value() == 5);
-static_assert((2.5 * m * (2 * percent)).numerical_value() == 5);
-static_assert((2.5L * (2 * m)).numerical_value() == 5);
-static_assert((2.5L * one * (2 * m)).numerical_value() == 5);
-static_assert((2.5L * percent * (2 * m)).numerical_value() == 5);
-static_assert((5. * m / 2).numerical_value() == 2.5);
-static_assert((5. * m / (2 * one)).numerical_value() == 2.5);
-static_assert((5. * m / (2 * percent)).numerical_value() == 2.5);
-static_assert((5. * km / (2 * m)).numerical_value() == 2.5);
-static_assert((5000. * m / (2 * m)).numerical_value() == 2500);
-
-static_assert((1 * m + 1.5 * m).numerical_value() == 2.5);
-static_assert((1 * m + 1.5 * km).numerical_value() == 1501);
-static_assert((1 * km + 1.5 * m).numerical_value() == 1001.5);
-static_assert((2 * m - 1.5 * m).numerical_value() == 0.5);
-static_assert((1 * km - 1.5 * m).numerical_value() == 998.5);
-static_assert((2 * m * 2.5L).numerical_value() == 5);
-static_assert((2 * m * (2.5L * one)).numerical_value() == 5);
-static_assert((2 * m * (2.5L * percent)).numerical_value() == 5);
-static_assert((2 * 2.5 * m).numerical_value() == 5);
-static_assert((2 * one * (2.5 * m)).numerical_value() == 5);
-static_assert((2 * percent * (2.5 * m)).numerical_value() == 5);
-static_assert((5 * m / 2.5L).numerical_value() == 2);
-static_assert((5 * m / (2.5L * one)).numerical_value() == 2);
-static_assert((5 * m / (2.5L * percent)).numerical_value() == 2);
-static_assert((5 * km / (2.5 * m)).numerical_value() == 2);
-static_assert((5000 * m / (2.5 * m)).numerical_value() == 2000);
-
-static_assert((7 * m % (2 * m)).numerical_value() == 1);
-static_assert((7 * km % (2000 * m)).numerical_value() == 1000);
-static_assert((1300 * m % (1 * km)).numerical_value() == 300);
-static_assert((7 * one % (2 * one)).numerical_value() == 1);
+static_assert((1 * m + 1 * m).numerical_value_ == 2);
+static_assert((1 * m + 1 * km).numerical_value_ == 1001);
+static_assert((1 * km + 1 * m).numerical_value_ == 1001);
+static_assert((2 * m - 1 * m).numerical_value_ == 1);
+static_assert((1 * km - 1 * m).numerical_value_ == 999);
+static_assert((2 * m * 2).numerical_value_ == 4);
+static_assert((2 * m * (2 * one)).numerical_value_ == 4);
+static_assert((2 * m * (2 * percent)).numerical_value_ == 4);
+static_assert((3 * 3 * m).numerical_value_ == 9);
+static_assert(((3 * one) * (3 * m)).numerical_value_ == 9);
+static_assert(((3 * percent) * (3 * m)).numerical_value_ == 9);
+static_assert((4 * m / 2).numerical_value_ == 2);
+static_assert((4 * m / (2 * one)).numerical_value_ == 2);
+static_assert((4 * m / (2 * percent)).numerical_value_ == 2);
+static_assert((4 * km / (2 * m)).numerical_value_ == 2);
+static_assert((4000 * m / (2 * m)).numerical_value_ == 2000);
+
+static_assert((1.5 * m + 1 * m).numerical_value_ == 2.5);
+static_assert((1.5 * m + 1 * km).numerical_value_ == 1001.5);
+static_assert((1.5 * km + 1 * m).numerical_value_ == 1501);
+static_assert((2.5 * m - 1 * m).numerical_value_ == 1.5);
+static_assert((1.5 * km - 1 * m).numerical_value_ == 1499);
+static_assert((2.5 * m * 2).numerical_value_ == 5);
+static_assert((2.5 * m * (2 * one)).numerical_value_ == 5);
+static_assert((2.5 * m * (2 * percent)).numerical_value_ == 5);
+static_assert((2.5L * (2 * m)).numerical_value_ == 5);
+static_assert((2.5L * one * (2 * m)).numerical_value_ == 5);
+static_assert((2.5L * percent * (2 * m)).numerical_value_ == 5);
+static_assert((5. * m / 2).numerical_value_ == 2.5);
+static_assert((5. * m / (2 * one)).numerical_value_ == 2.5);
+static_assert((5. * m / (2 * percent)).numerical_value_ == 2.5);
+static_assert((5. * km / (2 * m)).numerical_value_ == 2.5);
+static_assert((5000. * m / (2 * m)).numerical_value_ == 2500);
+
+static_assert((1 * m + 1.5 * m).numerical_value_ == 2.5);
+static_assert((1 * m + 1.5 * km).numerical_value_ == 1501);
+static_assert((1 * km + 1.5 * m).numerical_value_ == 1001.5);
+static_assert((2 * m - 1.5 * m).numerical_value_ == 0.5);
+static_assert((1 * km - 1.5 * m).numerical_value_ == 998.5);
+static_assert((2 * m * 2.5L).numerical_value_ == 5);
+static_assert((2 * m * (2.5L * one)).numerical_value_ == 5);
+static_assert((2 * m * (2.5L * percent)).numerical_value_ == 5);
+static_assert((2 * 2.5 * m).numerical_value_ == 5);
+static_assert((2 * one * (2.5 * m)).numerical_value_ == 5);
+static_assert((2 * percent * (2.5 * m)).numerical_value_ == 5);
+static_assert((5 * m / 2.5L).numerical_value_ == 2);
+static_assert((5 * m / (2.5L * one)).numerical_value_ == 2);
+static_assert((5 * m / (2.5L * percent)).numerical_value_ == 2);
+static_assert((5 * km / (2.5 * m)).numerical_value_ == 2);
+static_assert((5000 * m / (2.5 * m)).numerical_value_ == 2000);
+
+static_assert((7 * m % (2 * m)).numerical_value_ == 1);
+static_assert((7 * km % (2000 * m)).numerical_value_ == 1000);
+static_assert((1300 * m % (1 * km)).numerical_value_ == 300);
+static_assert((7 * one % (2 * one)).numerical_value_ == 1);
 
 static_assert((10 * m2 * (10 * m2)) / (50 * m2) == 2 * m2);
 
-static_assert((10 * km / (5 * m)).numerical_value() == 2);
-static_assert((10 * km / (5 * m)).in(one).numerical_value() == 2000);
-static_assert((10 * s * (2 * kHz)).numerical_value() == 20);
+static_assert((10 * km / (5 * m)).numerical_value_ == 2);
+static_assert((10 * km / (5 * m)).numerical_value_in(one) == 2000);
+static_assert((10 * s * (2 * kHz)).numerical_value_ == 20);
 
 // commutativity and associativity
 static_assert(10 * isq::length[si::metre] / (2 * isq::time[s]) + 5 * isq::speed[m / s] == 10 * isq::speed[m / s]);
@@ -740,14 +742,13 @@ static_assert(is_same_v<decltype(0.0 * one - 0 * one), decltype(0.0 * one)>);
 static_assert(1 * one - 30 * percent == (100 - 30) * percent);
 static_assert(1 * one + 30 * percent == (100 + 30) * percent);
 
-static_assert(is_same_v<decltype((std::uint8_t(0) * one + std::uint8_t(0) * one).numerical_value()), int&&>);
-static_assert(is_same_v<decltype((std::uint8_t(0) * one - std::uint8_t(0) * one).numerical_value()), int&&>);
-static_assert((std::uint8_t(128) * one + std::uint8_t(128) * one).numerical_value() ==
+static_assert(is_same_v<decltype(std::uint8_t(0) * one + std::uint8_t(0) * one)::rep, int>);
+static_assert(is_same_v<decltype(std::uint8_t(0) * one - std::uint8_t(0) * one)::rep, int>);
+static_assert((std::uint8_t(128) * one + std::uint8_t(128) * one).numerical_value_ ==
               std::uint8_t(128) + std::uint8_t(128));
-static_assert((std::uint8_t(0) * one - std::uint8_t(1) * one).numerical_value() == std::uint8_t(0) - std::uint8_t(1));
-
-static_assert(is_same_v<decltype((std::uint8_t(0) * one % (std::uint8_t(0) * one)).numerical_value()),
-                        decltype(std::uint8_t(0) % std::uint8_t(0))&&>);
+static_assert((std::uint8_t(0) * one - std::uint8_t(1) * one).numerical_value_ == std::uint8_t(0) - std::uint8_t(1));
+static_assert(is_same_v<decltype(std::uint8_t(0) * one % (std::uint8_t(0) * one))::rep,
+                        decltype(std::uint8_t(0) % std::uint8_t(0))>);
 
 static_assert(2 * one * (1 * m) == 2 * m);
 static_assert(2 * one / (1 * m) == 2 / (1 * m));
@@ -880,20 +881,24 @@ static_assert(!(123 * km >= 321'000 * m));
 
 static_assert(is_of_type<10 * km / (5 * km), quantity<one, int>>);
 
-static_assert((50. * m / (100. * m)).in(percent).numerical_value() == 50);
+static_assert((50. * m / (100. * m)).numerical_value_in(percent) == 50);
 static_assert(50. * m / (100. * m) == 50 * percent);
 
-static_assert((50. * percent).in(one).numerical_value() == 0.5);
+static_assert((50. * percent).numerical_value_in(one) == 0.5);
 
 
 //////////////////
 // value_cast
 //////////////////
 
-static_assert(value_cast<m>(2 * km).numerical_value() == 2000);
-static_assert(value_cast<km>(2000 * m).numerical_value() == 2);
-static_assert(value_cast<int>(1.23 * m).numerical_value() == 1);
-static_assert(value_cast<km / h>(2000.0 * m / (3600.0 * s)).numerical_value() == 2);
+static_assert(value_cast<m>(2 * km).numerical_value_ == 2000);
+static_assert(value_cast<km>(2000 * m).numerical_value_ == 2);
+static_assert(value_cast<int>(1.23 * m).numerical_value_ == 1);
+static_assert(value_cast<km / h>(2000.0 * m / (3600.0 * s)).numerical_value_ == 2);
+
+static_assert((2 * km).force_in(m).numerical_value_ == 2000);
+static_assert((2000 * m).force_in(km).numerical_value_ == 2);
+static_assert((2000.0 * m / (3600.0 * s)).force_in(km / h).numerical_value_ == 2);
 
 //////////////////
 // quantity_cast