feedc0de/mp-units
1
0
Fork 0
forked from mpusz/mp-units
Code Pull Requests Activity

rational powers

Browse Source
This commit is contained in:
dkavolis
2020-10-17 20:50:58 +01:00
committed by Mateusz Pusz
parent d534d2c62e
commit 23dd21465d
9 changed files with 555 additions and 108 deletions
Download Patch File Download Diff File Expand all files Collapse all files

239
src/include/units/bits/constexpr_math.h Normal file
View File

@@ -0,0 +1,239 @@
// The MIT License (MIT)
//
// Copyright (c) 2018 Mateusz Pusz
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
#include <gsl/gsl_assert>
#include <units/bits/pow.h>
#include <units/bits/ratio_maths.h>
namespace units::detail {
struct decimal_fp {
double significand;
std::intmax_t mantissa;
};
[[nodiscard]] constexpr decimal_fp to_decimal(double v) noexcept
{
if (v == 0) {
return {.significand = 0.0, .mantissa = 0};
}
double significand = abs(v);
std::intmax_t mantissa = 0;
while (significand < 1) {
significand *= 10.0;
--mantissa;
}
while (significand >= 10) {
significand /= 10.0;
++mantissa;
}
if (v < 0) {
significand = -significand;
}
return {.significand = significand, .mantissa = mantissa};
}
/* approximate natural log as https://math.stackexchange.com/a/977836
far slower than std::log but works at compile time with similar accuracy
*/
[[nodiscard]] constexpr double constexpr_log(double v) noexcept
{
Expects(v > 0);
// lookup table to speed up convergence for all significand values
// significand values of 7 and greater benefit mostly as they now converge in 5 terms compared to O(10)-O(100)
// required without the table
//
// using python:
// >>> import math
// >>> for i in range(1, 100):
// ... print(f"/* log({i:>2d}) = */ {math.log(i):.16f},")
constexpr std::array<double, 99> log_table{
/* log( 1) = */ 0.0000000000000000,
/* log( 2) = */ 0.6931471805599453,
/* log( 3) = */ 1.0986122886681098,
/* log( 4) = */ 1.3862943611198906,
/* log( 5) = */ 1.6094379124341003,
/* log( 6) = */ 1.7917594692280550,
/* log( 7) = */ 1.9459101490553132,
/* log( 8) = */ 2.0794415416798357,
/* log( 9) = */ 2.1972245773362196,
/* log(10) = */ 2.3025850929940459,
/* log(11) = */ 2.3978952727983707,
/* log(12) = */ 2.4849066497880004,
/* log(13) = */ 2.5649493574615367,
/* log(14) = */ 2.6390573296152584,
/* log(15) = */ 2.7080502011022101,
/* log(16) = */ 2.7725887222397811,
/* log(17) = */ 2.8332133440562162,
/* log(18) = */ 2.8903717578961645,
/* log(19) = */ 2.9444389791664403,
/* log(20) = */ 2.9957322735539909,
/* log(21) = */ 3.0445224377234230,
/* log(22) = */ 3.0910424533583161,
/* log(23) = */ 3.1354942159291497,
/* log(24) = */ 3.1780538303479458,
/* log(25) = */ 3.2188758248682006,
/* log(26) = */ 3.2580965380214821,
/* log(27) = */ 3.2958368660043291,
/* log(28) = */ 3.3322045101752038,
/* log(29) = */ 3.3672958299864741,
/* log(30) = */ 3.4011973816621555,
/* log(31) = */ 3.4339872044851463,
/* log(32) = */ 3.4657359027997265,
/* log(33) = */ 3.4965075614664802,
/* log(34) = */ 3.5263605246161616,
/* log(35) = */ 3.5553480614894135,
/* log(36) = */ 3.5835189384561099,
/* log(37) = */ 3.6109179126442243,
/* log(38) = */ 3.6375861597263857,
/* log(39) = */ 3.6635616461296463,
/* log(40) = */ 3.6888794541139363,
/* log(41) = */ 3.7135720667043080,
/* log(42) = */ 3.7376696182833684,
/* log(43) = */ 3.7612001156935624,
/* log(44) = */ 3.7841896339182610,
/* log(45) = */ 3.8066624897703196,
/* log(46) = */ 3.8286413964890951,
/* log(47) = */ 3.8501476017100584,
/* log(48) = */ 3.8712010109078911,
/* log(49) = */ 3.8918202981106265,
/* log(50) = */ 3.9120230054281460,
/* log(51) = */ 3.9318256327243257,
/* log(52) = */ 3.9512437185814275,
/* log(53) = */ 3.9702919135521220,
/* log(54) = */ 3.9889840465642745,
/* log(55) = */ 4.0073331852324712,
/* log(56) = */ 4.0253516907351496,
/* log(57) = */ 4.0430512678345503,
/* log(58) = */ 4.0604430105464191,
/* log(59) = */ 4.0775374439057197,
/* log(60) = */ 4.0943445622221004,
/* log(61) = */ 4.1108738641733114,
/* log(62) = */ 4.1271343850450917,
/* log(63) = */ 4.1431347263915326,
/* log(64) = */ 4.1588830833596715,
/* log(65) = */ 4.1743872698956368,
/* log(66) = */ 4.1896547420264252,
/* log(67) = */ 4.2046926193909657,
/* log(68) = */ 4.2195077051761070,
/* log(69) = */ 4.2341065045972597,
/* log(70) = */ 4.2484952420493594,
/* log(71) = */ 4.2626798770413155,
/* log(72) = */ 4.2766661190160553,
/* log(73) = */ 4.2904594411483910,
/* log(74) = */ 4.3040650932041702,
/* log(75) = */ 4.3174881135363101,
/* log(76) = */ 4.3307333402863311,
/* log(77) = */ 4.3438054218536841,
/* log(78) = */ 4.3567088266895917,
/* log(79) = */ 4.3694478524670215,
/* log(80) = */ 4.3820266346738812,
/* log(81) = */ 4.3944491546724391,
/* log(82) = */ 4.4067192472642533,
/* log(83) = */ 4.4188406077965983,
/* log(84) = */ 4.4308167988433134,
/* log(85) = */ 4.4426512564903167,
/* log(86) = */ 4.4543472962535073,
/* log(87) = */ 4.4659081186545837,
/* log(88) = */ 4.4773368144782069,
/* log(89) = */ 4.4886363697321396,
/* log(90) = */ 4.4998096703302650,
/* log(91) = */ 4.5108595065168497,
/* log(92) = */ 4.5217885770490405,
/* log(93) = */ 4.5325994931532563,
/* log(94) = */ 4.5432947822700038,
/* log(95) = */ 4.5538768916005408,
/* log(96) = */ 4.5643481914678361,
/* log(97) = */ 4.5747109785033828,
/* log(98) = */ 4.5849674786705723,
/* log(99) = */ 4.5951198501345898,
};
decimal_fp x = to_decimal(v);
// dividing the significand by nearest lower value in [1.0, 1.1, 1.2, ..., 9.9] will greatly improve convergence
x.significand *= 10;
const auto isignificand = static_cast<std::size_t>(x.significand);
x.significand /= static_cast<double>(isignificand);
const double result = static_cast<double>(x.mantissa - 1) * log_table[9] + log_table[isignificand - 1];
// 1.0 <= significand < 1.1 converges rapidly
const double y = (x.significand - 1) / (x.significand + 1);
const double y_squared = y * y;
double sum = 0;
// 5 terms are needed for convergence to machine precision in the worst case scenario
for (int k = 4; k > 0; --k) {
sum = y_squared * (1 / (2 * static_cast<double>(k) + 1) + sum);
}
sum = 2 * y * (1 + sum); // k = 0 term
return result + sum;
}
/* approximate e^x as Taylor series e^x = 1 + x/1! + x^2/2! + x^3/3! +... where N is the order of the Taylor series
use https://math.stackexchange.com/a/1988927 to improve convergence for large values
larger Factor values improve convergence for all values but reduce the precision
*/
template<std::size_t N, std::intmax_t Factor = 256>
[[nodiscard]] constexpr double constexpr_exp(double v) noexcept requires requires { Factor > 0; }
{
if constexpr (N == 0) {
return 1.0;
} else {
constexpr auto coefficients = []() {
std::array<double, N> coeffs;
std::size_t factorial = 1;
for (std::size_t i = 0; i < N; ++i) {
factorial *= i + 1;
coeffs[i] = 1.0 / static_cast<double>(factorial);
}
return coeffs;
}();
const double x = v / static_cast<double>(Factor);
double result = 0;
for (auto i = static_cast<std::intmax_t>(N - 1); i >= 0; --i) {
result = x * (coefficients[static_cast<std::size_t>(i)] + result);
}
// for factors of power of 2 this should be replaced by log2(Factor) multiplications by the compiler
return pow_impl<Factor>(1 + result);
}
}
// default template arguments provide reasonable precision even for fairly large exponents
// see constexpr_exp for template arguments
template<std::size_t ExpOrder = 10, std::intmax_t Factor = 128>
[[nodiscard]] constexpr double constexpr_pow(double v, double exponent) noexcept
{
const double x = exponent * constexpr_log(v);
return constexpr_exp<ExpOrder, Factor>(x);
}
} // namespace units::detail

64
src/include/units/bits/dimension_op.h
View File

@@ -159,71 +159,43 @@ using dimension_multiply = TYPENAME detail::dimension_multiply_impl<D1, D2>::typ
template<Dimension D1, Dimension D2> template<Dimension D1, Dimension D2>
using dimension_divide = TYPENAME detail::dimension_multiply_impl<D1, dim_invert<D2>>::type; using dimension_divide = TYPENAME detail::dimension_multiply_impl<D1, dim_invert<D2>>::type;
// dimension_sqrt
namespace detail {
template<Dimension D>
struct dimension_sqrt_impl;
template<BaseDimension D>
struct dimension_sqrt_impl<D> {
using type = downcast_dimension<derived_dimension_base<exponent<D, 1, 2>>>;
};
template<BaseDimension D>
struct dimension_sqrt_impl<derived_dimension_base<exponent<D, 2>>> {
using type = D;
};
template<DerivedDimension D>
struct dimension_sqrt_impl<D> {
using type = TYPENAME dimension_sqrt_impl<typename D::downcast_base_type>::type;
};
template<typename... Es>
struct dimension_sqrt_impl<derived_dimension_base<Es...>> {
using type = downcast_dimension<derived_dimension_base<exponent_multiply<Es, 1, 2>...>>;
};
} // namespace detail
template<Dimension D>
using dimension_sqrt = TYPENAME detail::dimension_sqrt_impl<D>::type;
// dimension_pow // dimension_pow
namespace detail { namespace detail {
template<Dimension D, std::intmax_t N> template<Dimension D, std::intmax_t Num, std::intmax_t Den = 1>
struct dimension_pow_impl; struct dimension_pow_impl;
template<BaseDimension D, std::intmax_t N> template<BaseDimension D, std::intmax_t Num, std::intmax_t Den>
struct dimension_pow_impl<D, N> { struct dimension_pow_impl<D, Num, Den> {
using type = downcast_dimension<derived_dimension_base<exponent<D, N>>>; using type = downcast_dimension<derived_dimension_base<exponent<D, Num, Den>>>;
}; };
template<BaseDimension D> template<BaseDimension D>
struct dimension_pow_impl<D, 1> { struct dimension_pow_impl<D, 1, 1> {
using type = D; using type = D;
}; };
template<BaseDimension D, std::intmax_t N> template<BaseDimension D, std::intmax_t Num, std::intmax_t Den>
struct dimension_pow_impl<derived_dimension_base<exponent<D, 1, N>>, N> { struct dimension_pow_impl<derived_dimension_base<exponent<D, Den, Num>>, Num, Den> {
using type = D; using type = D;
}; };
template<DerivedDimension D, std::intmax_t N> template<DerivedDimension D, std::intmax_t Num, std::intmax_t Den>
struct dimension_pow_impl<D, N> { struct dimension_pow_impl<D, Num, Den> {
using type = TYPENAME dimension_pow_impl<downcast_base_t<D>, N>::type; using type = TYPENAME dimension_pow_impl<downcast_base_t<D>, Num, Den>::type;
}; };
template<typename... Es, std::intmax_t N> template<typename... Es, std::intmax_t Num, std::intmax_t Den>
struct dimension_pow_impl<derived_dimension_base<Es...>, N> { struct dimension_pow_impl<derived_dimension_base<Es...>, Num, Den> {
using type = downcast_dimension<derived_dimension_base<exponent_multiply<Es, N, 1>...>>; using type = downcast_dimension<derived_dimension_base<exponent_multiply<Es, Num, Den>...>>;
}; };
} // namespace detail } // namespace detail
template<Dimension D, std::intmax_t N> template<Dimension D, std::intmax_t Num, std::intmax_t Den = 1>
using dimension_pow = TYPENAME detail::dimension_pow_impl<D, N>::type; using dimension_pow = TYPENAME detail::dimension_pow_impl<D, Num, Den>::type;
template<Dimension D>
using dimension_sqrt = TYPENAME detail::dimension_pow_impl<D, 1, 2>::type;
} // namespace units } // namespace units

16
src/include/units/bits/pow.h
View File

@@ -58,4 +58,20 @@ constexpr Rep fpow10(std::intmax_t exp)
return result; return result;
} }
template<std::intmax_t N, typename T>
constexpr T pow_impl(const T& v) noexcept
{
if constexpr (N == 0) {
return T(1);
} else if constexpr (N == 1) {
return v;
} else if constexpr (N < 0) {
return 1 / pow_impl<-N>(v);
} else if constexpr (N % 2 == 0) { // even
return pow_impl<N / 2>(v * v);
} else { // odd
return v * pow_impl<(N - 1) / 2>(v * v);
}
}
} // namespace units::detail } // namespace units::detail

94
src/include/units/bits/root.h Normal file
View File

@@ -0,0 +1,94 @@
// The MIT License (MIT)
//
// Copyright (c) 2018 Mateusz Pusz
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
#include <gsl/gsl_assert>
#include <units/bits/constexpr_math.h>
#include <units/bits/pow.h>
#include <units/bits/ratio_maths.h>
#include <cmath>
namespace units::detail {
template<std::intmax_t N, typename F>
[[nodiscard]] constexpr std::intmax_t iroot_impl(std::intmax_t v, F const& pow_function) noexcept requires requires { N > 0; }
{
if constexpr (N == 1) {
return v;
} else {
Expects(v >= 0);
if (v == 0) {
return 0;
}
constexpr double exponent = 1.0 / static_cast<double>(N);
const auto root = static_cast<std::intmax_t>(pow_function(static_cast<double>(v), exponent));
// integer roots may be truncated down by 1 or in even rarer cases up by 1 due to finite precision of pow and
// exponent, check both cases
if (v == pow_impl<N>(root + 1)) {
return root + 1;
}
if (v < pow_impl<N>(root)) {
return root - 1;
}
return root;
}
}
// maximum v is std::numeric_limits<std::intmax_t>::max() which is the worst case for exp convergence
// ExpOrder = 12 and Factor = 64 give a precision of about O(1e-15) for a wide range of 1 / N exponents
// Factor = 32 needs quite a few more terms to converge
// https://godbolt.org/z/odWq1o
template<std::intmax_t N, std::size_t ExpOrder = 12, std::intmax_t Factor = 64>
[[nodiscard]] constexpr std::intmax_t iroot_compile(std::intmax_t v) noexcept requires requires { N > 0; }
{
return iroot_impl<N>(v, [](double x, double exponent) { return constexpr_pow<ExpOrder, Factor>(x, exponent); });
}
template<std::intmax_t N>
[[nodiscard]] std::intmax_t iroot_runtime(std::intmax_t v) noexcept requires requires { N > 0; }
{
return iroot_impl<N>(v, [](double x, double exponent) {
if constexpr (N == 2) {
return std::sqrt(x);
} else if constexpr (N == 3) {
return std::cbrt(x);
} else {
return std::pow(x, exponent);
}
});
}
template<std::intmax_t N>
[[nodiscard]] constexpr std::intmax_t iroot(std::intmax_t v) noexcept requires requires { N > 0; }
{
// compile time version is much slower, use faster version at runtime
if (std::is_constant_evaluated()) {
return iroot_compile<N>(v);
}
return iroot_runtime<N>(v);
}
} // namespace units::detail

39
src/include/units/math.h
View File

@@ -34,22 +34,23 @@ namespace units {
* *
* Both the quantity value and its dimension are the base of the operation. * Both the quantity value and its dimension are the base of the operation.
* *
* @tparam N Exponent * @tparam Num Exponent numerator
* @tparam Den Exponent denominator
* @param q Quantity being the base of the operation * @param q Quantity being the base of the operation
* @return Quantity The result of computation * @return Quantity The result of computation
*/ */
template<std::intmax_t N, Quantity Q> template<std::intmax_t Num, std::intmax_t Den = 1, Quantity Q>
[[nodiscard]] inline auto pow(const Q& q) noexcept [[nodiscard]] inline auto pow(const Q& q) noexcept
requires requires { std::pow(q.count(), N); } requires requires { std::pow(q.count(), 1.0); Den != 0; }
{ {
using rep = TYPENAME Q::rep; using rep = TYPENAME Q::rep;
if constexpr(N == 0) { if constexpr (Num == 0) {
return rep(1); return rep(1);
} } else {
else { using dim = dimension_pow<typename Q::dimension, Num, Den>;
using dim = dimension_pow<typename Q::dimension, N>; using unit = downcast_unit<dim, pow<Num, Den>(Q::unit::ratio)>;
using unit = downcast_unit<dim, pow<N>(Q::unit::ratio)>; return quantity<dim, unit, rep>(
return quantity<dim, unit, rep>(static_cast<rep>(std::pow(q.count(), N))); static_cast<rep>(std::pow(q.count(), static_cast<double>(Num) / static_cast<double>(Den))));
} }
} }
@@ -65,12 +66,30 @@ template<Quantity Q>
[[nodiscard]] inline Quantity auto sqrt(const Q& q) noexcept [[nodiscard]] inline Quantity auto sqrt(const Q& q) noexcept
requires requires { std::sqrt(q.count()); } requires requires { std::sqrt(q.count()); }
{ {
using dim = dimension_sqrt<typename Q::dimension>; using dim = dimension_pow<typename Q::dimension, 1, 2>;
using unit = downcast_unit<dim, sqrt(Q::unit::ratio)>; using unit = downcast_unit<dim, sqrt(Q::unit::ratio)>;
using rep = TYPENAME Q::rep; using rep = TYPENAME Q::rep;
return quantity<dim, unit, rep>(static_cast<rep>(std::sqrt(q.count()))); return quantity<dim, unit, rep>(static_cast<rep>(std::sqrt(q.count())));
} }
/**
* @brief Computes the cubic root of a quantity
*
* Both the quantity value and its dimension are the base of the operation.
*
* @param q Quantity being the base of the operation
* @return Quantity The result of computation
*/
template<Quantity Q>
[[nodiscard]] inline Quantity auto cbrt(const Q& q) noexcept
requires requires { std::cbrt(q.count()); }
{
using dim = dimension_pow<typename Q::dimension, 1, 3>;
using unit = downcast_unit<dim, cbrt(Q::unit::ratio)>;
using rep = TYPENAME Q::rep;
return quantity<dim, unit, rep>(static_cast<rep>(std::cbrt(q.count())));
}
/** /**
* @brief Computes Euler's raised to the given power * @brief Computes Euler's raised to the given power
* *

90
src/include/units/ratio.h
View File

@@ -24,6 +24,8 @@
#include <units/bits/external/hacks.h> #include <units/bits/external/hacks.h>
#include <units/bits/ratio_maths.h> #include <units/bits/ratio_maths.h>
#include <units/bits/pow.h>
#include <units/bits/root.h>
#include <cstdint> #include <cstdint>
#include <numeric> #include <numeric>
#include <type_traits> #include <type_traits>
@@ -85,64 +87,70 @@ struct ratio {
} }
} }
template<std::intmax_t N>
[[nodiscard]] constexpr ratio pow(const ratio& r)
{
if constexpr(N == 0)
return ratio(1);
else if constexpr(N == 1)
return r;
else
return pow<N-1>(r) * r;
}
namespace detail { namespace detail {
// sqrt_impl avoids overflow and recursion [[nodiscard]] constexpr auto make_exp_align(const ratio& r, std::intmax_t alignment)
// from http://www.codecodex.com/wiki/Calculate_an_integer_square_root#C.2B.2B
// if v >= place this will fail (so we can't quite use the last bit)
[[nodiscard]] constexpr std::intmax_t sqrt_impl(std::intmax_t v)
{ {
// place = 0x4000 0000 for 32bit Expects(alignment > 0);
// place = 0x4000 0000 0000 0000 for 64bit const std::intmax_t rem = r.exp % alignment;
std::intmax_t place = static_cast<std::intmax_t>(1) << (sizeof(std::intmax_t) * 8 - 2);
while (place > v) place /= 4; // optimized by complier as place >>= 2
std::intmax_t root = 0; if (rem == 0) { // already aligned
while (place) { return std::array{r.num, r.den, r.exp};
if (v >= root + place) {
v -= root + place;
root += place * 2;
}
root /= 2;
place /= 4;
} }
return root;
if (r.exp > 0) { // remainder is positive
return std::array{r.num * ipow10(rem), r.den, r.exp - rem};
}
// remainder is negative
return std::array{r.num, r.den * ipow10(-rem), r.exp - rem};
} }
[[nodiscard]] constexpr auto make_exp_even(const ratio& r) template<std::intmax_t N>
[[nodiscard]] constexpr ratio root(const ratio& r) requires requires { N > 0; }
{ {
if(r.exp % 2 == 0) if constexpr (N == 1) {
return std::array{r.num, r.den, r.exp}; // already even (incl zero) return r;
} else {
if (r.num == 0) {
return ratio(0);
}
// safely make exp even, so it can be divided by 2 for square root const auto aligned = make_exp_align(r, N);
if(r.exp > 0) return ratio(iroot<N>(aligned[0]), iroot<N>(aligned[1]), aligned[2] / N);
return std::array{r.num * 10, r.den, r.exp - 1}; }
else
return std::array{r.num, r.den * 10, r.exp + 1};
} }
} // namespace detail } // namespace detail
[[nodiscard]] constexpr ratio sqrt(const ratio& r) template<std::intmax_t Num, std::intmax_t Den = 1>
[[nodiscard]] constexpr ratio pow(const ratio& r) requires requires { Den != 0; }
{ {
if(r.num == 0) if constexpr (Num == 0) {
return ratio(0); return ratio(1);
} else if constexpr (Num == Den) {
return r;
} else {
// simplify factors first and compute power for positive exponent
constexpr std::intmax_t gcd = std::gcd(Num, Den);
constexpr std::intmax_t num = detail::abs(Num / gcd);
constexpr std::intmax_t den = detail::abs(Den / gcd);
const auto even = detail::make_exp_even(r); // integer root loses precision so do pow first
return ratio(detail::sqrt_impl(even[0]), detail::sqrt_impl(even[1]), even[2] / 2); const ratio result = detail::root<den>(detail::pow_impl<num>(r));
if constexpr (Num * Den < 0) { // account for negative exponent
return inverse(result);
} else {
return result;
}
}
} }
[[nodiscard]] constexpr ratio sqrt(const ratio& r) { return pow<1, 2>(r); }
[[nodiscard]] constexpr ratio cbrt(const ratio& r) { return pow<1, 3>(r); }
// common_ratio // common_ratio
[[nodiscard]] constexpr ratio common_ratio(const ratio& r1, const ratio& r2) [[nodiscard]] constexpr ratio common_ratio(const ratio& r1, const ratio& r2)
{ {

81
test/unit_test/runtime/math_test.cpp
View File

@@ -54,6 +54,16 @@ TEST_CASE("'sqrt()' on quantity changes the value and the dimension accordingly"
REQUIRE(sqrt(4_q_m2) == 2_q_m); REQUIRE(sqrt(4_q_m2) == 2_q_m);
} }
TEST_CASE("'cbrt()' on quantity changes the value and the dimension accordingly", "[math][cbrt]")
{
REQUIRE(cbrt(8_q_m3) == 2_q_m);
}
TEST_CASE("'pow<Num, Den>()' on quantity changes the value and the dimension accordingly", "[math][pow]")
{
REQUIRE(pow<1, 4>(16_q_m2) == sqrt(4_q_m));
}
TEST_CASE("absolute functions on quantity returns the absolute value", "[math][abs][fabs]") TEST_CASE("absolute functions on quantity returns the absolute value", "[math][abs][fabs]")
{ {
SECTION ("'abs()' on a negative quantity returns the abs") SECTION ("'abs()' on a negative quantity returns the abs")
@@ -99,3 +109,74 @@ TEST_CASE("numeric_limits functions", "[limits]")
REQUIRE(epsilon<decltype(1_q_m)>().count() != std::numeric_limits<decltype(1._q_m)::rep>::epsilon()); REQUIRE(epsilon<decltype(1_q_m)>().count() != std::numeric_limits<decltype(1._q_m)::rep>::epsilon());
} }
} }
TEMPLATE_TEST_CASE_SIG("pow<N>() implementation exponentiates values to power N", "[math][pow][exp]",
(std::intmax_t N, N), 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 25)
{
auto v = GENERATE(range(0.5, 20.0, 0.5));
REQUIRE(detail::pow_impl<N>(v) == Approx(std::pow(v, N)).epsilon(1e-15).margin(0));
}
template<std::intmax_t N>
struct Pow {
constexpr static std::intmax_t exponent = N;
template<typename T>
[[nodiscard]] constexpr static auto pow(T const& v) noexcept
{
return detail::pow_impl<N>(v);
}
};
template<std::intmax_t N>
struct CompileRoot : Pow<N> {
[[nodiscard]] constexpr static std::intmax_t root(std::intmax_t v) noexcept { return detail::iroot_compile<N>(v); }
};
template<std::intmax_t N>
struct RuntimeRoot : Pow<N> {
[[nodiscard]] static std::intmax_t root(std::intmax_t v) noexcept { return detail::iroot_runtime<N>(v); }
};
// test to make sure precision is not lost when rounding what should be integer roots
template<typename TestType>
static void root_test()
{
SECTION ("Roots are truncated down") {
auto base = GENERATE(range(1.0, 10.0, 1.0)); // doubles to guard against overflow
if (TestType::pow(base) < static_cast<double>(std::numeric_limits<std::intmax_t>::max())) {
const std::intmax_t x = TestType::pow(static_cast<std::intmax_t>(base));
const auto expect = static_cast<std::intmax_t>(base);
REQUIRE(TestType::root(x - 1) == expect - 1);
REQUIRE(TestType::root(x) == expect);
}
}
SECTION ("Roots are truncated correctly for very large inputs") {
auto exponent = GENERATE(range(10, std::numeric_limits<std::intmax_t>::digits10, 1));
const auto large_val = static_cast<std::intmax_t>(std::pow(10, exponent));
const auto expected = static_cast<std::intmax_t>(std::pow(10, exponent / static_cast<double>(TestType::exponent)));
REQUIRE(TestType::root(large_val) == expected);
}
}
/* Catch2 uses int for indexing in TEMPLATE_PRODUCT_TEST_CASE_SIG so it does not compile with -Werror=sign-conversion
* https://github.com/catchorg/Catch2/pull/2074
TEMPLATE_PRODUCT_TEST_CASE_SIG("detail::iroot<N>()", "[math][pow][iroot]", (std::intmax_t N, N),
(CompileRoot, RuntimeRoot), (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 25))
{
root_test<TestType>();
}
*/
#define ROOT_TEST_CASE(Type) \
TEMPLATE_TEST_CASE_SIG("detail::iroot<N>() - " #Type, "[math][pow][iroot]", (std::intmax_t N, N), 1, 2, 3, 4, 5, 6, \
7, 8, 9, 10, 15, 20, 25) \
{ \
root_test<Type<N>>(); \
}
ROOT_TEST_CASE(CompileRoot)
ROOT_TEST_CASE(RuntimeRoot)

11
test/unit_test/static/math_test.cpp
View File

@@ -39,5 +39,16 @@ static_assert(compare<decltype(pow<2>(2_q_ft)), decltype(4_q_ft2)>);
static_assert(compare<decltype(sqrt(4_q_m2)), decltype(2_q_m)>); static_assert(compare<decltype(sqrt(4_q_m2)), decltype(2_q_m)>);
static_assert(compare<decltype(sqrt(4_q_km2)), decltype(2_q_km)>); static_assert(compare<decltype(sqrt(4_q_km2)), decltype(2_q_km)>);
static_assert(compare<decltype(sqrt(4_q_ft2)), decltype(2_q_ft)>); static_assert(compare<decltype(sqrt(4_q_ft2)), decltype(2_q_ft)>);
static_assert(compare<decltype(cbrt(8_q_m3)), decltype(2_q_m)>);
static_assert(compare<decltype(cbrt(8_q_km3)), decltype(2_q_km)>);
static_assert(compare<decltype(cbrt(8_q_ft3)), decltype(2_q_ft)>);
static_assert(compare<decltype(pow<1, 4>(4_q_m2 * 4_q_m2)), decltype(2_q_m)>);
static_assert(compare<decltype(pow<1, 4>(4_q_km2 * 4_q_km2)), decltype(2_q_km)>);
static_assert(compare<decltype(pow<1, 4>(4_q_ft2 * 4_q_ft2)), decltype(2_q_ft)>);
// rational dimensions
static_assert(compare<decltype(pow<1, 4>(4_q_m2)), decltype(sqrt(2_q_m))>);
static_assert(compare<decltype(pow<1, 4>(4_q_km2)), decltype(sqrt(2_q_km))>);
static_assert(compare<decltype(pow<1, 4>(4_q_ft2)), decltype(sqrt(2_q_ft))>);
} // namespace } // namespace

11
test/unit_test/static/ratio_test.cpp
View File

@@ -63,18 +63,25 @@ static_assert(pow<2>(ratio(1, 2)) == ratio(1, 4));
static_assert(pow<3>(ratio(1, 2)) == ratio(1, 8)); static_assert(pow<3>(ratio(1, 2)) == ratio(1, 8));
// pow with exponents // pow with exponents
static_assert(pow<2>(ratio(1, 2, 3)) == ratio(1, 4, 6)); static_assert(pow<2>(ratio(1, 2, 3)) == ratio(1, 4, 6));
static_assert(pow<4, 2>(ratio(1, 2, 3)) == ratio(1, 4, 6));
static_assert(pow<3>(ratio(1, 2, -6)) == ratio(1, 8, -18)); static_assert(pow<3>(ratio(1, 2, -6)) == ratio(1, 8, -18));
static_assert(sqrt(ratio(9)) == ratio(3)); static_assert(sqrt(ratio(9)) == ratio(3));
static_assert(cbrt(ratio(27)) == ratio(3));
static_assert(sqrt(ratio(4)) == ratio(2)); static_assert(sqrt(ratio(4)) == ratio(2));
static_assert(cbrt(ratio(8)) == ratio(2));
static_assert(sqrt(ratio(1)) == ratio(1)); static_assert(sqrt(ratio(1)) == ratio(1));
static_assert(cbrt(ratio(1)) == ratio(1));
static_assert(sqrt(ratio(0)) == ratio(0)); static_assert(sqrt(ratio(0)) == ratio(0));
static_assert(cbrt(ratio(0)) == ratio(0));
static_assert(sqrt(ratio(1, 4)) == ratio(1, 2)); static_assert(sqrt(ratio(1, 4)) == ratio(1, 2));
static_assert(cbrt(ratio(1, 8)) == ratio(1, 2));
// sqrt with exponents // sqrt with exponents
static_assert(sqrt(ratio(9, 1, 2)) == ratio(3, 1, 1)); static_assert(sqrt(ratio(9, 1, 2)) == ratio(3, 1, 1));
static_assert(sqrt(ratio(4)) == ratio(2)); static_assert(cbrt(ratio(27, 1, 3)) == ratio(3, 1, 1));
static_assert(cbrt(ratio(27, 1, 2)) == ratio(13, 1, 0));
// common_ratio // common_ratio
static_assert(common_ratio(ratio(1), ratio(1000)) == ratio(1)); static_assert(common_ratio(ratio(1), ratio(1000)) == ratio(1));