forked from mpusz/mp-units
After review
This commit is contained in:
@@ -23,40 +23,41 @@
|
||||
#pragma once
|
||||
|
||||
#include <gsl/gsl_assert>
|
||||
#include <units/bits/math_concepts.h>
|
||||
#include <units/bits/pow.h>
|
||||
#include <units/bits/ratio_maths.h>
|
||||
|
||||
namespace units::detail {
|
||||
|
||||
struct decimal_fp {
|
||||
double significand;
|
||||
std::intmax_t mantissa;
|
||||
double significant;
|
||||
std::intmax_t exponent;
|
||||
};
|
||||
|
||||
[[nodiscard]] constexpr decimal_fp to_decimal(double v) noexcept
|
||||
{
|
||||
if (v == 0) {
|
||||
return {.significand = 0.0, .mantissa = 0};
|
||||
return {.significant = 0.0, .exponent = 0};
|
||||
}
|
||||
|
||||
double significand = abs(v);
|
||||
std::intmax_t mantissa = 0;
|
||||
double significant = abs(v);
|
||||
std::intmax_t exponent = 0;
|
||||
|
||||
while (significand < 1) {
|
||||
significand *= 10.0;
|
||||
--mantissa;
|
||||
while (significant < 1) {
|
||||
significant *= 10.0;
|
||||
--exponent;
|
||||
}
|
||||
|
||||
while (significand >= 10) {
|
||||
significand /= 10.0;
|
||||
++mantissa;
|
||||
while (significant >= 10) {
|
||||
significant /= 10.0;
|
||||
++exponent;
|
||||
}
|
||||
|
||||
if (v < 0) {
|
||||
significand = -significand;
|
||||
significant = -significant;
|
||||
}
|
||||
|
||||
return {.significand = significand, .mantissa = mantissa};
|
||||
return {.significant = significant, .exponent = exponent};
|
||||
}
|
||||
|
||||
/* approximate natural log as https://math.stackexchange.com/a/977836
|
||||
@@ -66,8 +67,8 @@ struct decimal_fp {
|
||||
{
|
||||
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)
|
||||
// lookup table to speed up convergence for all significant values
|
||||
// significant 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:
|
||||
@@ -177,14 +178,14 @@ struct decimal_fp {
|
||||
};
|
||||
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];
|
||||
// dividing the significant by nearest lower value in [1.0, 1.1, 1.2, ..., 9.9] will greatly improve convergence
|
||||
x.significant *= 10;
|
||||
const auto isignificant = static_cast<std::size_t>(x.significant);
|
||||
x.significant /= static_cast<double>(isignificant);
|
||||
const double result = static_cast<double>(x.exponent - 1) * log_table[9] + log_table[isignificant - 1];
|
||||
|
||||
// 1.0 <= significand < 1.1 converges rapidly
|
||||
const double y = (x.significand - 1) / (x.significand + 1);
|
||||
// 1.0 <= significant < 1.1 converges rapidly
|
||||
const double y = (x.significant - 1) / (x.significant + 1);
|
||||
const double y_squared = y * y;
|
||||
double sum = 0;
|
||||
// 5 terms are needed for convergence to machine precision in the worst case scenario
|
||||
@@ -201,7 +202,8 @@ struct decimal_fp {
|
||||
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; }
|
||||
requires gt_zero<Factor>
|
||||
[[nodiscard]] constexpr double constexpr_exp(double v) noexcept
|
||||
{
|
||||
if constexpr (N == 0) {
|
||||
return 1.0;
|
||||
|
35
src/include/units/bits/math_concepts.h
Normal file
35
src/include/units/bits/math_concepts.h
Normal file
@@ -0,0 +1,35 @@
|
||||
// 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 <cstdint>
|
||||
|
||||
namespace units::detail {
|
||||
|
||||
template<std::intmax_t N>
|
||||
concept gt_zero = (N > 0);
|
||||
|
||||
template<std::intmax_t N>
|
||||
concept non_zero = (N != 0);
|
||||
|
||||
} // namespace units::detail
|
@@ -24,6 +24,7 @@
|
||||
|
||||
#include <gsl/gsl_assert>
|
||||
#include <units/bits/constexpr_math.h>
|
||||
#include <units/bits/math_concepts.h>
|
||||
#include <units/bits/pow.h>
|
||||
#include <units/bits/ratio_maths.h>
|
||||
#include <cmath>
|
||||
@@ -31,7 +32,8 @@
|
||||
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; }
|
||||
requires gt_zero<N>
|
||||
[[nodiscard]] constexpr std::intmax_t iroot_impl(std::intmax_t v, F const& pow_function) noexcept
|
||||
{
|
||||
if constexpr (N == 1) {
|
||||
return v;
|
||||
@@ -61,13 +63,15 @@ template<std::intmax_t N, typename F>
|
||||
// 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; }
|
||||
requires gt_zero<N>
|
||||
[[nodiscard]] constexpr std::intmax_t iroot_compile(std::intmax_t v) noexcept
|
||||
{
|
||||
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; }
|
||||
requires gt_zero<N>
|
||||
[[nodiscard]] std::intmax_t iroot_runtime(std::intmax_t v) noexcept
|
||||
{
|
||||
return iroot_impl<N>(v, [](double x, double exponent) {
|
||||
if constexpr (N == 2) {
|
||||
@@ -81,7 +85,8 @@ template<std::intmax_t N>
|
||||
}
|
||||
|
||||
template<std::intmax_t N>
|
||||
[[nodiscard]] constexpr std::intmax_t iroot(std::intmax_t v) noexcept requires requires { N > 0; }
|
||||
requires gt_zero<N>
|
||||
[[nodiscard]] constexpr std::intmax_t iroot(std::intmax_t v) noexcept
|
||||
{
|
||||
// compile time version is much slower, use faster version at runtime
|
||||
if (std::is_constant_evaluated()) {
|
||||
|
@@ -40,8 +40,8 @@ namespace units {
|
||||
* @return Quantity The result of computation
|
||||
*/
|
||||
template<std::intmax_t Num, std::intmax_t Den = 1, Quantity Q>
|
||||
[[nodiscard]] inline auto pow(const Q& q) noexcept
|
||||
requires requires { std::pow(q.count(), 1.0); Den != 0; }
|
||||
requires detail::non_zero<Den>
|
||||
[[nodiscard]] inline auto pow(const Q& q) noexcept requires requires { std::pow(q.count(), 1.0); }
|
||||
{
|
||||
using rep = TYPENAME Q::rep;
|
||||
if constexpr (Num == 0) {
|
||||
|
@@ -23,6 +23,7 @@
|
||||
#pragma once
|
||||
|
||||
#include <units/bits/external/hacks.h>
|
||||
#include <units/bits/math_concepts.h>
|
||||
#include <units/bits/ratio_maths.h>
|
||||
#include <units/bits/pow.h>
|
||||
#include <units/bits/root.h>
|
||||
@@ -107,7 +108,8 @@ namespace detail {
|
||||
}
|
||||
|
||||
template<std::intmax_t N>
|
||||
[[nodiscard]] constexpr ratio root(const ratio& r) requires requires { N > 0; }
|
||||
requires gt_zero<N>
|
||||
[[nodiscard]] constexpr ratio root(const ratio& r)
|
||||
{
|
||||
if constexpr (N == 1) {
|
||||
return r;
|
||||
@@ -124,7 +126,8 @@ template<std::intmax_t N>
|
||||
} // namespace detail
|
||||
|
||||
template<std::intmax_t Num, std::intmax_t Den = 1>
|
||||
[[nodiscard]] constexpr ratio pow(const ratio& r) requires requires { Den != 0; }
|
||||
requires detail::non_zero<Den>
|
||||
[[nodiscard]] constexpr ratio pow(const ratio& r)
|
||||
{
|
||||
if constexpr (Num == 0) {
|
||||
return ratio(1);
|
||||
|
Reference in New Issue
Block a user