mirror/Sprout
1
0
Fork 0
mirror of https://github.com/bolero-MURAKAMI/Sprout.git synced 2024-11-14 10:39:05 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

fix erf, erfc

Browse source 
fix for special values: trigonometric functions
This commit is contained in:
bolero-MURAKAMI 2013-04-23 23:07:20 +09:00
parent 5f40808f75
commit a27c83e939
12 changed files with 313 additions and 87 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

33
sprout/math/acos.hpp
View file

@ -5,25 +5,37 @@
#include <type_traits>
#include <sprout/config.hpp>
#include <sprout/math/detail/config.hpp>
#include <sprout/math/detail/float_compute.hpp>
#include <sprout/math/asin.hpp>
#include <sprout/math/fabs.hpp>
#include <sprout/math/constants.hpp>
#include <sprout/type_traits/enabler_if.hpp>
#include <sprout/type_traits/float_promote.hpp>
namespace sprout {
namespace math {
namespace detail {
template<typename T>
inline SPROUT_CONSTEXPR T
acos_impl(T x) {
return sprout::math::half_pi<T>() - sprout::math::asin(x);
}
template<
typename FloatType,
typename sprout::enabler_if<std::is_floating_point<FloatType>::value>::type = sprout::enabler
>
inline SPROUT_CONSTEXPR FloatType
acos(FloatType x) {
return x == 1 ? FloatType(0)
: x > 1 || x < -1 ? std::numeric_limits<FloatType>::quiet_NaN()
: sprout::math::half_pi<FloatType>() - sprout::asin(x)
return sprout::math::fabs(x) > 1 ? std::numeric_limits<FloatType>::quiet_NaN()
#if SPROUT_USE_BUILTIN_CMATH_FUNCTION
: std::acos(x)
#else
: x == 1 ? FloatType(0)
: static_cast<FloatType>(sprout::math::detail::acos_impl(static_cast<typename sprout::math::detail::float_compute<FloatType>::type>(x)))
#endif
;
}
template<
typename IntType,
typename sprout::enabler_if<std::is_integral<IntType>::value>::type = sprout::enabler
@ -33,8 +45,17 @@ namespace sprout {
return sprout::math::detail::acos(static_cast<double>(x));
}
} // namespace detail
using NS_SPROUT_MATH_DETAIL::acos;
//
// acos
//
template<
typename ArithmeticType,
typename sprout::enabler_if<std::is_arithmetic<ArithmeticType>::value>::type = sprout::enabler
>
inline SPROUT_CONSTEXPR typename sprout::float_promote<ArithmeticType>::type
acos(ArithmeticType x) {
return sprout::math::detail::acos(x);
}
} // namespace math
using sprout::math::acos;

35
sprout/math/asin.hpp
View file

@ -10,8 +10,10 @@
#include <sprout/math/detail/float_compute.hpp>
#include <sprout/math/constants.hpp>
#include <sprout/math/factorial.hpp>
#include <sprout/math/fabs.hpp>
#include <sprout/math/sqrt.hpp>
#include <sprout/type_traits/enabler_if.hpp>
#include <sprout/type_traits/float_promote.hpp>
namespace sprout {
namespace math {
@ -20,8 +22,8 @@ namespace sprout {
inline SPROUT_CONSTEXPR T
asin_impl_1(T x, std::size_t n, std::size_t last) {
return last - n == 1
? sprout::math::factorial<T>(2 * n)
/ sprout::detail::pow_n(T(4), n) / sprout::detail::pow2(sprout::math::factorial<T>(n)) / (2 * n + 1)
? sprout::math::unchecked_factorial<T>(2 * n)
/ sprout::detail::pow_n(T(4), n) / sprout::detail::pow2(sprout::math::unchecked_factorial<T>(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)
@ -42,14 +44,18 @@ namespace sprout {
>
inline SPROUT_CONSTEXPR FloatType
asin(FloatType x) {
typedef typename sprout::math::detail::float_compute<FloatType>::type type;
return x == 0 ? FloatType(0)
: x > 1 || x < -1 ? std::numeric_limits<FloatType>::quiet_NaN()
: x < 0 ? -static_cast<FloatType>(sprout::math::detail::asin_impl(static_cast<type>(-x)))
: static_cast<FloatType>(sprout::math::detail::asin_impl(static_cast<type>(x)))
return sprout::math::fabs(x) > 1 ? std::numeric_limits<FloatType>::quiet_NaN()
#if SPROUT_USE_BUILTIN_CMATH_FUNCTION
: std::asin(x)
#else
: x == 0 ? FloatType(0)
: static_cast<FloatType>(
x < 0 ? -sprout::math::detail::asin_impl(static_cast<typename sprout::math::detail::float_compute<FloatType>::type>(-x))
: sprout::math::detail::asin_impl(static_cast<typename sprout::math::detail::float_compute<FloatType>::type>(x))
)
#endif
;
}
template<
typename IntType,
typename sprout::enabler_if<std::is_integral<IntType>::value>::type = sprout::enabler
@ -59,8 +65,17 @@ namespace sprout {
return sprout::math::detail::asin(static_cast<double>(x));
}
} // namespace detail
using NS_SPROUT_MATH_DETAIL::asin;
//
// asin
//
template<
typename ArithmeticType,
typename sprout::enabler_if<std::is_arithmetic<ArithmeticType>::value>::type = sprout::enabler
>
inline SPROUT_CONSTEXPR typename sprout::float_promote<ArithmeticType>::type
asin(ArithmeticType x) {
return sprout::math::detail::asin(x);
}
} // namespace math
using sprout::math::asin;

24
sprout/math/atan.hpp
View file

@ -11,6 +11,7 @@
#include <sprout/math/constants.hpp>
#include <sprout/math/factorial.hpp>
#include <sprout/type_traits/enabler_if.hpp>
#include <sprout/type_traits/float_promote.hpp>
namespace sprout {
namespace math {
@ -41,17 +42,19 @@ namespace sprout {
>
inline SPROUT_CONSTEXPR FloatType
atan(FloatType x) {
typedef typename sprout::math::detail::float_compute<FloatType>::type type;
return x == 0 ? FloatType(0)
: x == std::numeric_limits<FloatType>::infinity() ? sprout::math::half_pi<FloatType>()
: x == -std::numeric_limits<FloatType>::infinity() ? -sprout::math::half_pi<FloatType>()
#if SPROUT_USE_BUILTIN_CMATH_FUNCTION
: std::atan(x)
#else
: static_cast<FloatType>(
x < 0 ? -sprout::math::detail::atan_impl(static_cast<type>(-x))
: sprout::math::detail::atan_impl(static_cast<type>(x))
x < 0 ? -sprout::math::detail::atan_impl(static_cast<typename sprout::math::detail::float_compute<FloatType>::type>(-x))
: sprout::math::detail::atan_impl(static_cast<typename sprout::math::detail::float_compute<FloatType>::type>(x))
)
#endif
;
}
template<
typename IntType,
typename sprout::enabler_if<std::is_integral<IntType>::value>::type = sprout::enabler
@ -61,8 +64,17 @@ namespace sprout {
return sprout::math::detail::atan(static_cast<double>(x));
}
} // namespace detail
using NS_SPROUT_MATH_DETAIL::atan;
//
// atan
//
template<
typename ArithmeticType,
typename sprout::enabler_if<std::is_arithmetic<ArithmeticType>::value>::type = sprout::enabler
>
inline SPROUT_CONSTEXPR typename sprout::float_promote<ArithmeticType>::type
atan(ArithmeticType x) {
return sprout::math::detail::atan(x);
}
} // namespace math
using sprout::math::atan;

60
sprout/math/atan2.hpp
View file

@ -5,40 +5,62 @@
#include <type_traits>
#include <sprout/config.hpp>
#include <sprout/math/detail/config.hpp>
#include <sprout/math/atan.hpp>
#include <sprout/math/detail/float_compute.hpp>
#include <sprout/math/constants.hpp>
#include <sprout/math/atan.hpp>
#include <sprout/type_traits/float_promote.hpp>
#include <sprout/type_traits/enabler_if.hpp>
#include <sprout/type_traits/float_promote.hpp>
namespace sprout {
namespace math {
namespace detail {
template<typename T>
inline SPROUT_CONSTEXPR T
atan2_impl(T y, T x) {
return x < 0 ? sprout::math::atan(y / x) + (
y < 0 ? -sprout::math::pi<T>()
: sprout::math::pi<T>()
)
: sprout::math::atan(y / x)
;
}
template<
typename FloatType,
typename sprout::enabler_if<std::is_floating_point<FloatType>::value>::type = sprout::enabler
>
inline SPROUT_CONSTEXPR FloatType
atan2(FloatType y, FloatType x) {
return y == 0
? x == 0 ? FloatType(0)
: x < 0 ? sprout::math::pi<FloatType>()
: FloatType(0)
: x == 0 ? (y < 0 ? -1 : 1) * sprout::math::half_pi<FloatType>()
: x == -std::numeric_limits<FloatType>::infinity()
return x == -std::numeric_limits<FloatType>::infinity()
? y == std::numeric_limits<FloatType>::infinity() ? sprout::math::three_quarters_pi<FloatType>()
: y == -std::numeric_limits<FloatType>::infinity() ? -sprout::math::three_quarters_pi<FloatType>()
: (y < 0 ? -1 : 1) * sprout::math::half_pi<FloatType>()
: y < 0 ? -sprout::math::half_pi<FloatType>()
: sprout::math::half_pi<FloatType>()
: x == std::numeric_limits<FloatType>::infinity()
? y == std::numeric_limits<FloatType>::infinity() ? sprout::math::quarter_pi<FloatType>()
: y == -std::numeric_limits<FloatType>::infinity() ? -sprout::math::quarter_pi<FloatType>()
: FloatType(0)
: y == std::numeric_limits<FloatType>::infinity() ? sprout::math::half_pi<FloatType>()
: y == -std::numeric_limits<FloatType>::infinity() ? -sprout::math::half_pi<FloatType>()
: x < 0 ? sprout::atan(y / x) + (y < 0 ? -1 : 1) * sprout::math::pi<FloatType>()
: sprout::atan(y / x)
#if SPROUT_USE_BUILTIN_CMATH_FUNCTION
: std::atan2(y, x)
#else
: y == 0
? x < 0 ? sprout::math::pi<FloatType>()
: FloatType(0)
: x == 0
? y < 0 ? -sprout::math::half_pi<FloatType>()
: sprout::math::half_pi<FloatType>()
: static_cast<FloatType>(
sprout::math::detail::atan2_impl(
static_cast<typename sprout::math::detail::float_compute<FloatType>::type>(y),
static_cast<typename sprout::math::detail::float_compute<FloatType>::type>(x)
)
)
#endif
;
}
template<
typename ArithmeticType1,
typename ArithmeticType2,
@ -52,8 +74,20 @@ namespace sprout {
return sprout::math::detail::atan2(static_cast<type>(y), static_cast<type>(x));
}
} // namespace detail
using NS_SPROUT_MATH_DETAIL::atan2;
//
// atan2
//
template<
typename ArithmeticType1,
typename ArithmeticType2,
typename sprout::enabler_if<
std::is_arithmetic<ArithmeticType1>::value && std::is_arithmetic<ArithmeticType2>::value
>::type = sprout::enabler
>
inline SPROUT_CONSTEXPR typename sprout::float_promote<ArithmeticType1, ArithmeticType2>::type
atan2(ArithmeticType1 y, ArithmeticType2 x) {
return sprout::math::detail::atan2(y, x);
}
} // namespace math
using sprout::math::atan2;

41
sprout/math/cos.hpp
View file

@ -12,6 +12,7 @@
#include <sprout/math/constants.hpp>
#include <sprout/math/fmod.hpp>
#include <sprout/type_traits/enabler_if.hpp>
#include <sprout/type_traits/float_promote.hpp>
namespace sprout {
namespace math {
@ -25,15 +26,12 @@ namespace sprout {
+ sprout::math::detail::cos_impl_1(x2, n + (last - n) / 2, last)
;
}
template<typename FloatType>
inline SPROUT_CONSTEXPR FloatType
cos_impl(FloatType x) {
typedef typename sprout::math::detail::float_compute<FloatType>::type type;
return static_cast<FloatType>(
type(1) + sprout::math::detail::cos_impl_1(
sprout::detail::pow2(static_cast<type>(x)),
1, sprout::math::factorial_limit<type>() / 2 + 1
)
template<typename T>
inline SPROUT_CONSTEXPR T
cos_impl(T x) {
return T(1) + sprout::math::detail::cos_impl_1(
sprout::detail::pow2(sprout::math::fmod(x, sprout::math::two_pi<T>())),
1, sprout::math::factorial_limit<T>() / 2 + 1
);
}
@ -43,13 +41,17 @@ namespace sprout {
>
inline SPROUT_CONSTEXPR FloatType
cos(FloatType x) {
return x == 0 ? FloatType(1)
: x == std::numeric_limits<FloatType>::infinity() || x == -std::numeric_limits<FloatType>::infinity()
return x == std::numeric_limits<FloatType>::infinity()
|| x == -std::numeric_limits<FloatType>::infinity()
? std::numeric_limits<FloatType>::quiet_NaN()
: sprout::math::detail::cos_impl(sprout::math::fmod(x, sprout::math::two_pi<FloatType>()))
#if SPROUT_USE_BUILTIN_CMATH_FUNCTION
: std::cos(x)
#else
: x == 0 ? FloatType(1)
: static_cast<FloatType>(sprout::math::detail::cos_impl(static_cast<typename sprout::math::detail::float_compute<FloatType>::type>(x)))
#endif
;
}
template<
typename IntType,
typename sprout::enabler_if<std::is_integral<IntType>::value>::type = sprout::enabler
@ -59,8 +61,17 @@ namespace sprout {
return sprout::math::detail::cos(static_cast<double>(x));
}
} // namespace detail
using NS_SPROUT_MATH_DETAIL::cos;
//
// cos
//
template<
typename ArithmeticType,
typename sprout::enabler_if<std::is_arithmetic<ArithmeticType>::value>::type = sprout::enabler
>
inline SPROUT_CONSTEXPR typename sprout::float_promote<ArithmeticType>::type
cos(ArithmeticType x) {
return sprout::math::detail::cos(x);
}
} // namespace math
using sprout::math::cos;

116
sprout/math/erf.hpp
View file

@ -1,15 +1,14 @@
#ifndef SPROUT_MATH_ERF_HPP
#define SPROUT_MATH_ERF_HPP
#include <cstddef>
#include <limits>
#include <type_traits>
#include <sprout/config.hpp>
#include <sprout/detail/pow.hpp>
#include <sprout/math/detail/config.hpp>
#include <sprout/math/detail/float_compute.hpp>
#include <sprout/math/factorial.hpp>
#include <sprout/math/constants.hpp>
#include <sprout/math/itrunc.hpp>
#include <sprout/math/fractional_part.hpp>
#include <sprout/detail/pow.hpp>
#include <sprout/type_traits/enabler_if.hpp>
namespace sprout {
@ -17,25 +16,108 @@ namespace sprout {
namespace detail {
template<typename T>
inline SPROUT_CONSTEXPR T
erf_impl_1(T x, std::size_t n, std::size_t last) {
return last - n == 1
? (n % 2 ? -1 : 1) / sprout::math::factorial<T>(n) / (2 * n + 1) * sprout::detail::pow_n(x, 2 * n + 1)
: sprout::math::detail::erf_impl_1(x, n, n + (last - n) / 2)
+ sprout::math::detail::erf_impl_1(x, n + (last - n) / 2, last)
erf_impl_3(T x, T y) {
return x < 0 ? -y
: y
;
}
template<typename T>
inline SPROUT_CONSTEXPR T
erf_impl_2_b(T x, T w, T t, int k) {
return sprout::math::detail::erf_impl_3(
x,
T(1) - sprout::detail::pow_n(
k == 0 ? (((((((((((
-2.9734388465e-10 * t + 2.69776334046e-9) * t
+ -6.40788827665e-9) * t + -1.6678201321e-8) * t + -2.1854388148686e-7) * t
+ 2.66246030457984e-6) * t + 1.612722157047886e-5) * t + -2.5616361025506629e-4) * t
+ 1.5380842432375365e-4) * t + 0.00815533022524927908) * t + -0.01402283663896319337) * t
+ -0.19746892495383021487) * t + 0.71511720328842845913
: k == 1 ? (((((((((((
-1.951073787e-11 * t + -3.2302692214e-10) * t
+ 5.22461866919e-9) * t + 3.42940918551e-9) * t + -3.5772874310272e-7) * t
+ 1.9999935792654e-7) * t + 2.687044575042908e-5) * t + -1.1843240273775776e-4) * t
+ -8.0991728956032271e-4) * t + 0.00661062970502241174) * t + 0.00909530922354827295) * t
+ -0.2016007277849101314) * t + 0.51169696718727644908
: k == 2 ? (((((((((((
3.147682272e-11 * t + -4.8465972408e-10) * t
+ 6.3675740242e-10) * t + 3.377623323271e-8) * t + -1.5451139637086e-7) * t
+ -2.03340624738438e-6) * t + 1.947204525295057e-5) * t + 2.854147231653228e-5) * t
+ -0.00101565063152200272) * t + 0.00271187003520095655) * t + 0.02328095035422810727) * t
+ -0.16725021123116877197) * t + 0.32490054966649436974
: k == 3 ? (((((((((((
2.31936337e-11 * t + -6.303206648e-11) * t
+ -2.64888267434e-9) * t + 2.050708040581e-8) * t + 1.1371857327578e-7) * t
+ -2.11211337219663e-6) * t + 3.68797328322935e-6) * t + 9.823686253424796e-5) * t
+ -6.5860243990455368e-4) * t + -7.5285814895230877e-4) * t + 0.02585434424202960464) * t
+ -0.11637092784486193258) * t + 0.18267336775296612024
: (((((((((((
-3.67789363e-12 * t + 2.0876046746e-10) * t
+ -1.93319027226e-9) * t + -4.35953392472e-9) * t + 1.8006992266137e-7) * t
+ -7.8441223763969e-7) * t + -6.75407647949153e-6) * t + 8.428418334440096e-5) * t
+ -1.7604388937031815e-4) * t + -0.0023972961143507161) * t + 0.0206412902387602297) * t
+ -0.06905562880005864105) * t + 0.09084526782065478489
,
16
)
);
}
template<typename T>
inline SPROUT_CONSTEXPR T
erf_impl_2_a_1(T x, T w, T t, int k) {
return sprout::math::detail::erf_impl_3(
x,
k == 0 ? ((((((((((((
5.958930743e-11 * t + -1.13739022964e-9) * t
+ 1.466005199839e-8) * t + -1.635035446196e-7) * t + 1.6461004480962e-6) * t
+ -1.492559551950604e-5) * t + 1.2055331122299265e-4) * t + -8.548326981129666e-4) * t
+ 0.00522397762482322257) * t + -0.0268661706450773342) * t + 0.11283791670954881569) * t
+ -0.37612638903183748117) * t + 1.12837916709551257377) * w
: k == 1 ? ((((((((((((
2.372510631e-11 * t + -4.5493253732e-10) * t
+ 5.90362766598e-9) * t + -6.642090827576e-8) * t + 6.7595634268133e-7) * t
+ -6.21188515924e-6) * t + 5.10388300970969e-5) * t + -3.7015410692956173e-4) * t
+ 0.00233307631218880978) * t + -0.0125498847718219221) * t + 0.05657061146827041994) * t
+ -0.2137966477645600658) * t + 0.84270079294971486929) * w
: k == 2 ? ((((((((((((
9.49905026e-12 * t + -1.8310229805e-10) * t
+ 2.39463074e-9) * t + -2.721444369609e-8) * t + 2.8045522331686e-7) * t
+ -2.61830022482897e-6) * t + 2.195455056768781e-5) * t + -1.6358986921372656e-4) * t
+ 0.00107052153564110318) * t + -0.00608284718113590151) * t + 0.02986978465246258244) * t
+ -0.13055593046562267625) * t + 0.67493323603965504676) * w
: k == 3 ? ((((((((((((
3.82722073e-12 * t + -7.421598602e-11) * t
+ 9.793057408e-10) * t + -1.126008898854e-8) * t + 1.1775134830784e-7) * t
+ -1.1199275838265e-6) * t + 9.62023443095201e-6) * t + -7.404402135070773e-5) * t
+ 5.0689993654144881e-4) * t + -0.00307553051439272889) * t + 0.01668977892553165586) * t
+ -0.08548534594781312114) * t + 0.56909076642393639985) * w
: ((((((((((((
1.55296588e-12 * t + -3.032205868e-11) * t
+ 4.0424830707e-10) * t + -4.71135111493e-9) * t + 5.011915876293e-8) * t
+ -4.8722516178974e-7) * t + 4.30683284629395e-6) * t + -3.445026145385764e-5) * t
+ 2.4879276133931664e-4) * t + -0.00162940941748079288) * t + 0.00988786373932350462) * t
+ -0.05962426839442303805) * t + 0.49766113250947636708) * w
);
}
template<typename T>
inline SPROUT_CONSTEXPR T
erf_impl_2_a(T x, T w, T t) {
return sprout::math::detail::erf_impl_2_a_1(x, w, sprout::fractional_part(t), sprout::itrunc(t));
}
template<typename T>
inline SPROUT_CONSTEXPR T
erf_impl_1(T x, T w) {
return w < T(2.2) ? sprout::math::detail::erf_impl_2_a(x, w, w * w)
: w < T(6.9) ? sprout::math::detail::erf_impl_2_b(x, w, sprout::fractional_part(w), sprout::itrunc(w) - 2)
: sprout::math::detail::erf_impl_3(x, T(1))
;
}
template<typename T>
inline SPROUT_CONSTEXPR T
erf_impl(T x) {
return sprout::math::two_div_root_pi<T>()
* (x + sprout::math::detail::erf_impl_1(x, 1, sprout::math::factorial_limit<T>() + 1))
;
}
template<typename T>
inline SPROUT_CONSTEXPR T
erf_impl(T x2, T a) {
return T(1) - sprout::exp(-x2 / (T(1) + a * x2) * (sprout::math::quarter_pi<T>() + a * x2));
return sprout::math::detail::erf_impl_1(x, x < 0 ? -x : x);
}
template<
typename FloatType,
typename sprout::enabler_if<std::is_floating_point<FloatType>::value>::type = sprout::enabler

12
sprout/math/factorial.hpp
View file

@ -763,17 +763,27 @@ namespace sprout {
return sprout::math::detail::factorials<type>::limit;
}
//
// unchecked_factorial
//
template<typename T, typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
inline SPROUT_CONSTEXPR T unchecked_factorial(std::size_t x) {
typedef typename std::remove_cv<T>::type type;
return sprout::math::detail::factorials<type>::table[x];
}
//
// factorial
//
template<typename T, typename = typename std::enable_if<std::is_arithmetic<T>::value>::type>
inline SPROUT_CONSTEXPR T factorial(std::size_t x) {
typedef typename std::remove_cv<T>::type type;
return SPROUT_ASSERT(x <= sprout::math::factorial_limit<type>()),
sprout::math::detail::factorials<type>::table[x]
sprout::math::unchecked_factorial<T>(x)
;
}
} // namespace math
using sprout::math::factorial_limit;
using sprout::math::unchecked_factorial;
using sprout::math::factorial;
} // namespace sprout

4
sprout/math/hypot.hpp
View file

@ -19,8 +19,8 @@ namespace sprout {
>
inline SPROUT_CONSTEXPR FloatType
hypot(FloatType x, FloatType y) {
return y == 0 ? sprout::fabs(x)
: x == 0 ? sprout::fabs(y)
return y == 0 ? sprout::math::fabs(x)
: x == 0 ? sprout::math::fabs(y)
: y == std::numeric_limits<FloatType>::infinity() || y == -std::numeric_limits<FloatType>::infinity()
? std::numeric_limits<FloatType>::infinity()
: x == std::numeric_limits<FloatType>::infinity() || x == -std::numeric_limits<FloatType>::infinity()

4
sprout/math/lgamma.hpp
View file

@ -10,7 +10,7 @@
#include <sprout/math/factorial.hpp>
#include <sprout/math/log.hpp>
#include <sprout/math/sin.hpp>
#include <sprout/math/abs.hpp>
#include <sprout/math/fabs.hpp>
#include <sprout/math/itrunc.hpp>
#include <sprout/type_traits/enabler_if.hpp>
@ -20,7 +20,7 @@ namespace sprout {
template<typename T>
inline SPROUT_CONSTEXPR T
lgamma_impl_3(T x, T y) {
return x < 0 ? sprout::log(sprout::math::pi<T>() / sprout::abs(x * sprout::sin(x * sprout::math::pi<T>()))) - y
return x < 0 ? sprout::log(sprout::math::pi<T>() / sprout::math::fabs(x * sprout::math::sin(x * sprout::math::pi<T>()))) - y
: y
;
}

35
sprout/math/sin.hpp
View file

@ -5,26 +5,38 @@
#include <type_traits>
#include <sprout/config.hpp>
#include <sprout/math/detail/config.hpp>
#include <sprout/math/cos.hpp>
#include <sprout/math/detail/float_compute.hpp>
#include <sprout/math/constants.hpp>
#include <sprout/math/cos.hpp>
#include <sprout/type_traits/enabler_if.hpp>
#include <sprout/type_traits/float_promote.hpp>
namespace sprout {
namespace math {
namespace detail {
template<typename T>
inline SPROUT_CONSTEXPR T
sin_impl(T x) {
return -sprout::math::cos(x + sprout::math::half_pi<T>());
}
template<
typename FloatType,
typename sprout::enabler_if<std::is_floating_point<FloatType>::value>::type = sprout::enabler
>
inline SPROUT_CONSTEXPR FloatType
sin(FloatType x) {
return x == 0 ? FloatType(0)
: x == std::numeric_limits<FloatType>::infinity() || x == -std::numeric_limits<FloatType>::infinity()
return x == std::numeric_limits<FloatType>::infinity()
|| x == -std::numeric_limits<FloatType>::infinity()
? std::numeric_limits<FloatType>::quiet_NaN()
: -sprout::cos(x + sprout::math::half_pi<FloatType>())
#if SPROUT_USE_BUILTIN_CMATH_FUNCTION
: std::sin(x)
#else
: x == 0 ? FloatType(0)
: static_cast<FloatType>(sprout::math::detail::sin_impl(static_cast<typename sprout::math::detail::float_compute<FloatType>::type>(x)))
#endif
;
}
template<
typename IntType,
typename sprout::enabler_if<std::is_integral<IntType>::value>::type = sprout::enabler
@ -34,8 +46,17 @@ namespace sprout {
return sprout::math::detail::sin(static_cast<double>(x));
}
} // namespace detail
using NS_SPROUT_MATH_DETAIL::sin;
//
// sin
//
template<
typename ArithmeticType,
typename sprout::enabler_if<std::is_arithmetic<ArithmeticType>::value>::type = sprout::enabler
>
inline SPROUT_CONSTEXPR typename sprout::float_promote<ArithmeticType>::type
sin(ArithmeticType x) {
return sprout::math::detail::sin(x);
}
} // namespace math
using sprout::math::sin;

34
sprout/math/tan.hpp
View file

@ -5,26 +5,37 @@
#include <type_traits>
#include <sprout/config.hpp>
#include <sprout/math/detail/config.hpp>
#include <sprout/math/sin.hpp>
#include <sprout/math/detail/float_compute.hpp>
#include <sprout/math/cos.hpp>
#include <sprout/math/sin.hpp>
#include <sprout/type_traits/enabler_if.hpp>
#include <sprout/type_traits/float_promote.hpp>
namespace sprout {
namespace math {
namespace detail {
template<typename T>
inline SPROUT_CONSTEXPR T
tan_impl(T x) {
return sprout::math::sin(x) / sprout::math::cos(x);
}
template<
typename FloatType,
typename sprout::enabler_if<std::is_floating_point<FloatType>::value>::type = sprout::enabler
>
inline SPROUT_CONSTEXPR FloatType
tan(FloatType x) {
return x == 0 ? FloatType(0)
: x == std::numeric_limits<FloatType>::infinity() || x == -std::numeric_limits<FloatType>::infinity()
return x == std::numeric_limits<FloatType>::infinity() || x == -std::numeric_limits<FloatType>::infinity()
? std::numeric_limits<FloatType>::quiet_NaN()
: sprout::sin(x) / sprout::cos(x)
#if SPROUT_USE_BUILTIN_CMATH_FUNCTION
: std::sin(x)
#else
: x == 0 ? FloatType(0)
: static_cast<FloatType>(sprout::math::detail::tan_impl(static_cast<typename sprout::math::detail::float_compute<FloatType>::type>(x)))
#endif
;
}
template<
typename IntType,
typename sprout::enabler_if<std::is_integral<IntType>::value>::type = sprout::enabler
@ -34,8 +45,17 @@ namespace sprout {
return sprout::math::detail::tan(static_cast<double>(x));
}
} // namespace detail
using NS_SPROUT_MATH_DETAIL::tan;
//
// tan
//
template<
typename ArithmeticType,
typename sprout::enabler_if<std::is_arithmetic<ArithmeticType>::value>::type = sprout::enabler
>
inline SPROUT_CONSTEXPR typename sprout::float_promote<ArithmeticType>::type
tan(ArithmeticType x) {
return sprout::math::detail::tan(x);
}
} // namespace math
using sprout::math::tan;

2
sprout/math/tanh.hpp
View file

@ -21,7 +21,7 @@ namespace sprout {
return x == 0 ? FloatType(0)
: x == std::numeric_limits<FloatType>::infinity() ? FloatType(1)
: x == -std::numeric_limits<FloatType>::infinity() ? FloatType(-1)
: sprout::sinh(x) / sprout::cosh(x)
: sprout::math::sinh(x) / sprout::math::cosh(x)
;
}