diff --git a/corsika/framework/units/quantity.hpp b/corsika/framework/units/quantity.hpp index 3e897cfc80151beebdd2fb700af35fd27ae94e47..7b5bd35877b755fe1d751965975c148d42a8f20e 100644 --- a/corsika/framework/units/quantity.hpp +++ b/corsika/framework/units/quantity.hpp @@ -468,6 +468,9 @@ namespace phys { template <typename D, typename X> friend detail::Root<D, 2, X> sqrt(quantity<D, X> const& x); + template <typename D, typename X> + friend detail::Root<D, 3, X> cbrt(quantity<D, X> const& x); + // comparison template <typename D, typename X, typename Y> @@ -676,7 +679,17 @@ namespace phys { static_assert(detail::root<D, 2, X>::all_even_multiples, "root result dimensions must be integral"); - return detail::Root<D, 2, X>(std::pow(x.m_value, X(1.0) / 2)); + return detail::Root<D, 2, X>(std::sqrt(x.m_value)); + } + + /// cubic root. + + template <typename D, typename X> + detail::Root<D, 3, X> cbrt(quantity<D, X> const& x) { + static_assert(detail::root<D, 3, X>::all_even_multiples, + "root result dimensions must be integral"); + + return detail::Root<D, 3, X>(std::cbrt(x.m_value)); } // Comparison operators diff --git a/tests/framework/testUnits.cpp b/tests/framework/testUnits.cpp index 67db321d8734fc488a833348dca02c120a831b52..8d15544aca0d5c2dc1923af5853e6de7cafeb4f4 100644 --- a/tests/framework/testUnits.cpp +++ b/tests/framework/testUnits.cpp @@ -92,6 +92,9 @@ TEST_CASE("PhysicalUnits", "[Units]") { const auto E3 = E2 + 100_GeV + pow(10, lgE) * 1_GeV; CHECK(E3 == 180_GeV); + + CHECK(sqrt(5_GeV * 5_GeV) / 5_GeV == Approx(1)); + CHECK(cbrt(static_pow<3>(5_GeV)) / 5_GeV == Approx(1)); } SECTION("Output") {