diff --git a/sprout/math/asin.hpp b/sprout/math/asin.hpp index 14304e06..4c3b585d 100644 --- a/sprout/math/asin.hpp +++ b/sprout/math/asin.hpp @@ -9,7 +9,6 @@ #include #include #include -#include #include #include #include @@ -20,21 +19,54 @@ namespace sprout { namespace detail { template inline SPROUT_CONSTEXPR T - asin_impl_1(T x, std::size_t n, std::size_t last) { - return last - n == 1 - ? sprout::math::unchecked_factorial(2 * n) - / sprout::detail::pow_n(T(4), n) / sprout::detail::pow2(sprout::math::unchecked_factorial(n)) / (2 * n + 1) - * sprout::detail::pow_n(x, 2 * n + 1) - : sprout::math::detail::asin_impl_1(x, n, n + (last - n) / 2) - + sprout::math::detail::asin_impl_1(x, n + (last - n) / 2, last) + asin_impl_center_1(T x, T x2) { + return (((((((((((( + + 0.0316658385792867081040808) * x2 + + -0.0158620440988475212803145) * x2 + + 0.0192942786775238654913582) * x2 + + 0.0066153165197009078340075) * x2 + + 0.0121483892822292648695383) * x2 + + 0.0138885410156894774969889) * x2 + + 0.0173593516996479249428647) * x2 + + 0.0223717830666671020710108) * x2 + + 0.0303819580081956423799529) * x2 + + 0.0446428568582815922683933) * x2 + + 0.0750000000029696112392353) * x2 + + 0.1666666666666558995379880) * x2 + * x + x + ; + } + template + inline SPROUT_CONSTEXPR T + asin_impl_center(T x) { + return sprout::math::detail::asin_impl_center_1(x, x * x); + } + template + inline SPROUT_CONSTEXPR T + asin_impl_tail(T x) { + return sprout::math::half_pi() + sprout::math::sqrt(T(1) - x) + * ((((((((((((( + + -0.0000121189820098929624806) * x + + 0.0001307564187657962919394) * x + + -0.0006702485124770180942917) * x + + 0.0021912255981979442677477) * x + + -0.0052049731575223952626203) * x + + 0.0097868293573384001221447) * x + + -0.0156746038587246716524035) * x + + 0.0229883479552557203133368) * x + + -0.0331919619444009606270380) * x + + 0.0506659694457588602631748) * x + + -0.0890259194305537131666744) * x + + 0.2145993335526539017488949) * x + + -1.5707961988153774692344105) ; } template inline SPROUT_CONSTEXPR T asin_impl(T x) { - return x > sprout::math::half_root_two() - ? sprout::math::half_pi() - sprout::math::detail::asin_impl_1(sprout::math::sqrt(1 - x * x), 0, sprout::math::factorial_limit() / 2 + 1) - : x + sprout::math::detail::asin_impl_1(x, 1, sprout::math::factorial_limit() / 2 + 1) + return x < T(-0.5) ? -sprout::math::detail::asin_impl_tail(-x) + : x > T(0.5) ? sprout::math::detail::asin_impl_tail(x) + : sprout::math::detail::asin_impl_center(x) ; } diff --git a/sprout/math/cos.hpp b/sprout/math/cos.hpp index a1389c8f..8e323200 100644 --- a/sprout/math/cos.hpp +++ b/sprout/math/cos.hpp @@ -8,8 +8,11 @@ #include #include #include +//#include +//#include #include #include +#include #include #include #include @@ -35,6 +38,31 @@ namespace sprout { ); } +// template +// inline SPROUT_CONSTEXPR T +// cos_impl_2(T x) { +// return x <= sprout::math::quarter_pi() ? -sprout::math::detail::cosp(x) +// : x <= sprout::math::half_pi() ? -sprout::math::detail::sinp(sprout::math::half_pi() - x) +// : x <= sprout::math::three_quarters_pi() ? sprout::math::detail::sinp(x - sprout::math::half_pi()) +// : sprout::math::detail::cosp(sprout::math::pi() - x) +// ; +// } +// template +// inline SPROUT_CONSTEXPR T +// cos_impl_1(T x) { +// return x > sprout::math::pi() ? sprout::math::detail::cos_impl_2(x - sprout::math::pi()) +// : x <= sprout::math::quarter_pi() ? sprout::math::detail::cosp(x) +// : x <= sprout::math::half_pi() ? sprout::math::detail::sinp(sprout::math::half_pi() - x) +// : x <= sprout::math::three_quarters_pi() ? -sprout::math::detail::sinp(x - sprout::math::half_pi()) +// : -sprout::math::detail::cosp(sprout::math::pi() - x) +// ; +// } +// template +// inline SPROUT_CONSTEXPR T +// cos_impl(T x) { +// return sprout::math::detail::cos_impl_1(sprout::math::fmod(sprout::math::fabs(x), sprout::math::two_pi())); +// } + template< typename FloatType, typename sprout::enabler_if::value>::type = sprout::enabler diff --git a/sprout/math/cosh.hpp b/sprout/math/cosh.hpp index af8e0c20..f91d0d32 100644 --- a/sprout/math/cosh.hpp +++ b/sprout/math/cosh.hpp @@ -5,11 +5,10 @@ #include #include #include -#include #include #include #include -#include +#include #include namespace sprout { @@ -17,20 +16,13 @@ namespace sprout { namespace detail { template inline SPROUT_CONSTEXPR T - cosh_impl_1(T x2, std::size_t n, std::size_t last) { - return last - n == 1 - ? sprout::detail::pow_n(x2, n) / sprout::math::unchecked_factorial(2 * n) - : sprout::math::detail::cosh_impl_1(x2, n, n + (last - n) / 2) - + sprout::math::detail::cosh_impl_1(x2, n + (last - n) / 2, last) - ; + cosh_impl_1(T t) { + return T(0.5) * t + T(0.5) / t; } template inline SPROUT_CONSTEXPR T cosh_impl(T x) { - return T(1) + sprout::math::detail::cosh_impl_1( - x * x, - 1, sprout::math::factorial_limit() / 2 + 1 - ); + return sprout::math::detail::cosh_impl_1(sprout::math::exp(x)); } template< diff --git a/sprout/math/detail/cosp.hpp b/sprout/math/detail/cosp.hpp new file mode 100644 index 00000000..b47ad6da --- /dev/null +++ b/sprout/math/detail/cosp.hpp @@ -0,0 +1,32 @@ +#ifndef SPROUT_MATH_DETAIL_COSP_HPP +#define SPROUT_MATH_DETAIL_COSP_HPP + +#include + +namespace sprout { + namespace math { + namespace detail { + template + inline SPROUT_CONSTEXPR T + cosp_impl(T x2) { + return ((((((( + -1.13585365213876817300e-11) * x2 + + 2.08757008419747316778e-9) * x2 + - 2.75573141792967388112e-7) * x2 + + 2.48015872888517045348e-5) * x2 + - 1.38888888888730564116e-3) * x2 + + 4.16666666666665929218e-2) * x2 + - 0.5) * x2 + + 1.0 + ; + } + template + inline SPROUT_CONSTEXPR T + cosp(T x) { + return sprout::math::detail::cosp_impl(x * x); + } + } // namespace detail + } // namespace math +} // namespace sprout + +#endif // #ifndef SPROUT_MATH_DETAIL_COSP_HPP diff --git a/sprout/math/detail/sinp.hpp b/sprout/math/detail/sinp.hpp new file mode 100644 index 00000000..1400ac00 --- /dev/null +++ b/sprout/math/detail/sinp.hpp @@ -0,0 +1,30 @@ +#ifndef SPROUT_MATH_DETAIL_SINP_HPP +#define SPROUT_MATH_DETAIL_SINP_HPP + +#include + +namespace sprout { + namespace math { + namespace detail { + template + inline SPROUT_CONSTEXPR T + sinp_impl(T x, T x2) { + return x + x * ((((((( + 1.58962301576546568060e-10) * x2 + - 2.50507477628578072866e-8) * x2 + + 2.75573136213857245213e-6) * x2 + - 1.98412698295895385996e-4) * x2 + + 8.33333333332211858878e-3) * x2 + - 1.66666666666666307295e-1) * x2 + ); + } + template + inline SPROUT_CONSTEXPR T + sinp(T x) { + return sprout::math::detail::sinp_impl(x, x * x); + } + } // namespace detail + } // namespace math +} // namespace sprout + +#endif // #ifndef SPROUT_MATH_DETAIL_SINP_HPP diff --git a/sprout/math/exp.hpp b/sprout/math/exp.hpp index ea560720..c01aa614 100644 --- a/sprout/math/exp.hpp +++ b/sprout/math/exp.hpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include @@ -37,7 +38,8 @@ namespace sprout { > inline SPROUT_CONSTEXPR FloatType exp(FloatType x) { - return x == -std::numeric_limits::infinity() ? FloatType(0) + return sprout::math::isnan(x) ? x + : x == -std::numeric_limits::infinity() ? FloatType(0) : x == std::numeric_limits::infinity() ? std::numeric_limits::infinity() #if SPROUT_USE_BUILTIN_CMATH_FUNCTION : std::exp(x) diff --git a/sprout/math/exp10.hpp b/sprout/math/exp10.hpp index 6cb82b12..ddfbc8b1 100644 --- a/sprout/math/exp10.hpp +++ b/sprout/math/exp10.hpp @@ -6,8 +6,9 @@ #include #include #include -#include #include +#include +#include #include namespace sprout { @@ -25,7 +26,8 @@ namespace sprout { > inline SPROUT_CONSTEXPR FloatType exp10(FloatType x) { - return x == -std::numeric_limits::infinity() ? FloatType(0) + return sprout::math::isnan(x) ? x + : x == -std::numeric_limits::infinity() ? FloatType(0) : x == std::numeric_limits::infinity() ? std::numeric_limits::infinity() : x == 0 ? FloatType(1) : static_cast(sprout::math::detail::exp10_impl(static_cast::type>(x))) diff --git a/sprout/math/exp2.hpp b/sprout/math/exp2.hpp index 67d64f2e..9b88e002 100644 --- a/sprout/math/exp2.hpp +++ b/sprout/math/exp2.hpp @@ -6,8 +6,9 @@ #include #include #include -#include #include +#include +#include #include namespace sprout { @@ -25,7 +26,8 @@ namespace sprout { > inline SPROUT_CONSTEXPR FloatType exp2(FloatType x) { - return x == -std::numeric_limits::infinity() ? FloatType(0) + return sprout::math::isnan(x) ? x + : x == -std::numeric_limits::infinity() ? FloatType(0) : x == std::numeric_limits::infinity() ? std::numeric_limits::infinity() #if SPROUT_USE_BUILTIN_CMATH_FUNCTION : std::exp2(x) diff --git a/sprout/math/expm1.hpp b/sprout/math/expm1.hpp index ba3e1aa7..6e4cef0f 100644 --- a/sprout/math/expm1.hpp +++ b/sprout/math/expm1.hpp @@ -6,9 +6,10 @@ #include #include #include +#include +#include #include #include -#include #include namespace sprout { @@ -26,7 +27,8 @@ namespace sprout { > inline SPROUT_CONSTEXPR FloatType expm1(FloatType x) { - return x == -std::numeric_limits::infinity() ? FloatType(-1) + return sprout::math::isnan(x) ? x + : x == -std::numeric_limits::infinity() ? FloatType(-1) : x == std::numeric_limits::infinity() ? std::numeric_limits::infinity() #if SPROUT_USE_BUILTIN_CMATH_FUNCTION : std::expm1(x) diff --git a/sprout/math/ldexp.hpp b/sprout/math/ldexp.hpp index fd880849..f1390108 100644 --- a/sprout/math/ldexp.hpp +++ b/sprout/math/ldexp.hpp @@ -4,9 +4,10 @@ #include #include #include +#include #include #include -#include +#include #include namespace sprout { @@ -24,7 +25,8 @@ namespace sprout { > inline SPROUT_CONSTEXPR FloatType ldexp(FloatType x, int exp) { - return x == std::numeric_limits::infinity() ? std::numeric_limits::infinity() + return sprout::math::isnan(x) ? x + : x == std::numeric_limits::infinity() ? std::numeric_limits::infinity() : x == -std::numeric_limits::infinity() ? -std::numeric_limits::infinity() : exp == 0 ? x : x == 0 ? x diff --git a/sprout/math/log.hpp b/sprout/math/log.hpp index 9cd58059..ddd7d899 100644 --- a/sprout/math/log.hpp +++ b/sprout/math/log.hpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include @@ -47,7 +48,8 @@ namespace sprout { > inline SPROUT_CONSTEXPR FloatType log(FloatType x) { - return x == 0 ? -std::numeric_limits::infinity() + return sprout::math::isnan(x) ? x + : x == 0 ? -std::numeric_limits::infinity() : x == std::numeric_limits::infinity() ? std::numeric_limits::infinity() : x < 0 ? std::numeric_limits::quiet_NaN() #if SPROUT_USE_BUILTIN_CMATH_FUNCTION diff --git a/sprout/math/log10.hpp b/sprout/math/log10.hpp index 615dc6bf..236f00f4 100644 --- a/sprout/math/log10.hpp +++ b/sprout/math/log10.hpp @@ -6,8 +6,9 @@ #include #include #include -#include #include +#include +#include #include namespace sprout { @@ -25,7 +26,8 @@ namespace sprout { > inline SPROUT_CONSTEXPR FloatType log10(FloatType x) { - return x == 0 ? -std::numeric_limits::infinity() + return sprout::math::isnan(x) ? x + : x == 0 ? -std::numeric_limits::infinity() : x == std::numeric_limits::infinity() ? std::numeric_limits::infinity() : x < 0 ? std::numeric_limits::quiet_NaN() #if SPROUT_USE_BUILTIN_CMATH_FUNCTION diff --git a/sprout/math/log1p.hpp b/sprout/math/log1p.hpp index 6182f313..9a00908b 100644 --- a/sprout/math/log1p.hpp +++ b/sprout/math/log1p.hpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include @@ -24,9 +25,10 @@ namespace sprout { > inline SPROUT_CONSTEXPR FloatType log1p(FloatType x) { - return x == -1 ? -std::numeric_limits::infinity() + return sprout::math::isnan(x) ? x + : x == -1 ? -std::numeric_limits::infinity() : x == std::numeric_limits::infinity() ? std::numeric_limits::infinity() - : x < -1 ? std::numeric_limits::quiet_NaN() + : x < -1 ? -std::numeric_limits::quiet_NaN() #if SPROUT_USE_BUILTIN_CMATH_FUNCTION : std::log1p(x) #else diff --git a/sprout/math/log2.hpp b/sprout/math/log2.hpp index f54b5304..b7fdeafd 100644 --- a/sprout/math/log2.hpp +++ b/sprout/math/log2.hpp @@ -6,8 +6,9 @@ #include #include #include -#include #include +#include +#include #include namespace sprout { @@ -25,7 +26,8 @@ namespace sprout { > inline SPROUT_CONSTEXPR FloatType log2(FloatType x) { - return x == 0 ? -std::numeric_limits::infinity() + return sprout::math::isnan(x) ? x + : x == 0 ? -std::numeric_limits::infinity() : x == std::numeric_limits::infinity() ? std::numeric_limits::infinity() : x < 0 ? std::numeric_limits::quiet_NaN() #if SPROUT_USE_BUILTIN_CMATH_FUNCTION diff --git a/sprout/math/logb.hpp b/sprout/math/logb.hpp index 27c2a166..0b70c0f4 100644 --- a/sprout/math/logb.hpp +++ b/sprout/math/logb.hpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -39,7 +40,7 @@ namespace sprout { inline SPROUT_CONSTEXPR T logb_impl_3_pos_lo(T x, T x0, T base, T exp) { return base < 1 ? sprout::math::detail::logb_impl_3_pos_lo( - x, x0 * std::numeric_limits::radix, x / (x0 / std::numeric_limits::radix), exp + 1 + x, x0 * std::numeric_limits::radix, x / (x0 / std::numeric_limits::radix), exp - 1 ) : exp ; @@ -48,7 +49,7 @@ namespace sprout { inline SPROUT_CONSTEXPR T logb_impl_3_pos_hi(T x, T x0, T base, T exp) { return !(base < std::numeric_limits::radix) ? sprout::math::detail::logb_impl_3_pos_hi( - x, x0 / std::numeric_limits::radix, x / (x0 * std::numeric_limits::radix), exp - 1 + x, x0 / std::numeric_limits::radix, x / (x0 * std::numeric_limits::radix), exp + 1 ) : exp ; @@ -65,10 +66,10 @@ namespace sprout { ) : exp : base < 1 ? sprout::math::detail::logb_impl_3_pos_lo( - x, x0 * std::numeric_limits::radix, x / (x0 / std::numeric_limits::radix), exp + 1 + x, x0 * std::numeric_limits::radix, x / (x0 / std::numeric_limits::radix), exp - 1 ) : !(base < std::numeric_limits::radix) ? sprout::math::detail::logb_impl_3_pos_hi( - x, x0 / std::numeric_limits::radix, x / (x0 * std::numeric_limits::radix), exp - 1 + x, x0 / std::numeric_limits::radix, x / (x0 * std::numeric_limits::radix), exp + 1 ) : exp ; @@ -99,7 +100,8 @@ namespace sprout { > inline SPROUT_CONSTEXPR FloatType logb(FloatType x) { - return x == 0 ? -std::numeric_limits::infinity() + return sprout::math::isnan(x) ? x + : x == 0 ? -std::numeric_limits::infinity() #if SPROUT_USE_BUILTIN_CMATH_FUNCTION # if defined(__GNUC__) : x == -std::numeric_limits::infinity() diff --git a/sprout/math/logb2.hpp b/sprout/math/logb2.hpp index 87ff5595..f49ca33d 100644 --- a/sprout/math/logb2.hpp +++ b/sprout/math/logb2.hpp @@ -13,6 +13,7 @@ # include # include # include +# include # include # include # include @@ -61,7 +62,7 @@ namespace sprout { inline SPROUT_CONSTEXPR T logb2_impl_3_pos_lo(T x, T x0, T base, T exp) { return base < 1 ? sprout::math::detail::logb2_impl_3_pos_lo( - x, x0 * 2, x / (x0 / 2), exp + 1 + x, x0 * 2, x / (x0 / 2), exp - 1 ) : exp ; @@ -70,7 +71,7 @@ namespace sprout { inline SPROUT_CONSTEXPR T logb2_impl_3_pos_hi(T x, T x0, T base, T exp) { return !(base < 2) ? sprout::math::detail::logb2_impl_3_pos_hi( - x, x0 / 2, x / (x0 * 2), exp - 1 + x, x0 / 2, x / (x0 * 2), exp + 1 ) : exp ; @@ -87,10 +88,10 @@ namespace sprout { ) : exp : base < 1 ? sprout::math::detail::logb2_impl_3_pos_lo( - x, x0 * 2, x / (x0 / 2), exp + 1 + x, x0 * 2, x / (x0 / 2), exp - 1 ) : !(base < 2) ? sprout::math::detail::logb2_impl_3_pos_hi( - x, x0 / 2, x / (x0 * 2), exp - 1 + x, x0 / 2, x / (x0 * 2), exp + 1 ) : exp ; @@ -121,7 +122,8 @@ namespace sprout { > inline SPROUT_CONSTEXPR FloatType logb2(FloatType x) { - return x == 0 ? -std::numeric_limits::infinity() + return sprout::math::isnan(x) ? x + : x == 0 ? -std::numeric_limits::infinity() : x == std::numeric_limits::infinity() || x == -std::numeric_limits::infinity() ? std::numeric_limits::infinity() : static_cast(sprout::math::detail::logb2_impl(static_cast::type>(x))) diff --git a/sprout/math/scalbln.hpp b/sprout/math/scalbln.hpp index 4a8630ff..ff20c640 100644 --- a/sprout/math/scalbln.hpp +++ b/sprout/math/scalbln.hpp @@ -4,9 +4,10 @@ #include #include #include +#include #include #include -#include +#include #include namespace sprout { @@ -24,7 +25,8 @@ namespace sprout { > inline SPROUT_CONSTEXPR FloatType scalbln(FloatType x, long exp) { - return x == std::numeric_limits::infinity() ? std::numeric_limits::infinity() + return sprout::math::isnan(x) ? x + : x == std::numeric_limits::infinity() ? std::numeric_limits::infinity() : x == -std::numeric_limits::infinity() ? -std::numeric_limits::infinity() : exp == 0 ? x : x == 0 ? x diff --git a/sprout/math/scalbn.hpp b/sprout/math/scalbn.hpp index 90afe13e..256cf892 100644 --- a/sprout/math/scalbn.hpp +++ b/sprout/math/scalbn.hpp @@ -4,9 +4,10 @@ #include #include #include +#include #include #include -#include +#include #include namespace sprout { @@ -24,7 +25,8 @@ namespace sprout { > inline SPROUT_CONSTEXPR FloatType scalbn(FloatType x, int exp) { - return x == std::numeric_limits::infinity() ? std::numeric_limits::infinity() + return sprout::math::isnan(x) ? x + : x == std::numeric_limits::infinity() ? std::numeric_limits::infinity() : x == -std::numeric_limits::infinity() ? -std::numeric_limits::infinity() : exp == 0 ? x : x == 0 ? x diff --git a/sprout/math/sinh.hpp b/sprout/math/sinh.hpp index adf355ab..82ec686d 100644 --- a/sprout/math/sinh.hpp +++ b/sprout/math/sinh.hpp @@ -5,11 +5,10 @@ #include #include #include -#include #include #include #include -#include +#include #include namespace sprout { @@ -17,20 +16,13 @@ namespace sprout { namespace detail { template inline SPROUT_CONSTEXPR T - sinh_impl_1(T x, std::size_t n, std::size_t last) { - return last - n == 1 - ? sprout::detail::pow_n(x, 2 * n + 1) / sprout::math::unchecked_factorial(2 * n + 1) - : sprout::math::detail::sinh_impl_1(x, n, n + (last - n) / 2) - + sprout::math::detail::sinh_impl_1(x, n + (last - n) / 2, last) - ; + sinh_impl_1(T t) { + return T(0.5) * t - T(0.5) / t; } template inline SPROUT_CONSTEXPR T sinh_impl(T x) { - return x + sprout::math::detail::sinh_impl_1( - x, - 1, (sprout::math::factorial_limit() - 1) / 2 + 1 - ); + return sprout::math::detail::sinh_impl_1(sprout::math::exp(x)); } template< diff --git a/sprout/math/tanh.hpp b/sprout/math/tanh.hpp index 0bd1d9db..1c2725df 100644 --- a/sprout/math/tanh.hpp +++ b/sprout/math/tanh.hpp @@ -7,17 +7,26 @@ #include #include #include -#include -#include +#include #include namespace sprout { namespace math { namespace detail { + template + inline SPROUT_CONSTEXPR T + tanh_impl_2(T t, T u) { + return (t - u) / (t + u); + } + template + inline SPROUT_CONSTEXPR T + tanh_impl_1(T t) { + return sprout::math::detail::tanh_impl_2(t, T(1) / t); + } template inline SPROUT_CONSTEXPR T tanh_impl(T x) { - return sprout::math::sinh(x) / sprout::math::cosh(x); + return sprout::math::detail::tanh_impl_1(sprout::math::exp(x)); } template<