From fccb16687b67ea21c248850d7410203147cfc1bb Mon Sep 17 00:00:00 2001
From: bolero-MURAKAMI <bolero.murakami@gmail.com>
Date: Sun, 5 May 2013 00:06:30 +0900
Subject: [PATCH] fix math functions

---
 sprout/math/asin.hpp        | 54 +++++++++++++++++++++++++++++--------
 sprout/math/cos.hpp         | 28 +++++++++++++++++++
 sprout/math/cosh.hpp        | 16 +++--------
 sprout/math/detail/cosp.hpp | 32 ++++++++++++++++++++++
 sprout/math/detail/sinp.hpp | 30 +++++++++++++++++++++
 sprout/math/exp.hpp         |  4 ++-
 sprout/math/exp10.hpp       |  6 +++--
 sprout/math/exp2.hpp        |  6 +++--
 sprout/math/expm1.hpp       |  6 +++--
 sprout/math/ldexp.hpp       |  6 +++--
 sprout/math/log.hpp         |  4 ++-
 sprout/math/log10.hpp       |  6 +++--
 sprout/math/log1p.hpp       |  6 +++--
 sprout/math/log2.hpp        |  6 +++--
 sprout/math/logb.hpp        | 12 +++++----
 sprout/math/logb2.hpp       | 12 +++++----
 sprout/math/scalbln.hpp     |  6 +++--
 sprout/math/scalbn.hpp      |  6 +++--
 sprout/math/sinh.hpp        | 16 +++--------
 sprout/math/tanh.hpp        | 15 ++++++++---
 20 files changed, 209 insertions(+), 68 deletions(-)
 create mode 100644 sprout/math/detail/cosp.hpp
 create mode 100644 sprout/math/detail/sinp.hpp

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 <sprout/math/detail/config.hpp>
 #include <sprout/math/detail/float_compute.hpp>
 #include <sprout/math/constants.hpp>
-#include <sprout/math/factorial.hpp>
 #include <sprout/math/isnan.hpp>
 #include <sprout/math/fabs.hpp>
 #include <sprout/math/sqrt.hpp>
@@ -20,21 +19,54 @@ namespace sprout {
 		namespace detail {
 			template<typename T>
 			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<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)
+			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<typename T>
+			inline SPROUT_CONSTEXPR T
+			asin_impl_center(T x) {
+				return sprout::math::detail::asin_impl_center_1(x, x * x);
+			}
+			template<typename T>
+			inline SPROUT_CONSTEXPR T
+			asin_impl_tail(T x) {
+				return sprout::math::half_pi<T>() + 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<typename T>
 			inline SPROUT_CONSTEXPR T
 			asin_impl(T x) {
-				return x > sprout::math::half_root_two<T>()
-					? sprout::math::half_pi<T>() - sprout::math::detail::asin_impl_1(sprout::math::sqrt(1 - x * x), 0, sprout::math::factorial_limit<T>() / 2 + 1)
-					: x + sprout::math::detail::asin_impl_1(x, 1, sprout::math::factorial_limit<T>() / 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 <sprout/detail/pow.hpp>
 #include <sprout/math/detail/config.hpp>
 #include <sprout/math/detail/float_compute.hpp>
+//#include <sprout/math/detail/cosp.hpp>
+//#include <sprout/math/detail/sinp.hpp>
 #include <sprout/math/constants.hpp>
 #include <sprout/math/isnan.hpp>
+#include <sprout/math/fabs.hpp>
 #include <sprout/math/factorial.hpp>
 #include <sprout/math/fmod.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
@@ -35,6 +38,31 @@ namespace sprout {
 					);
 			}
 
+//			template<typename T>
+//			inline SPROUT_CONSTEXPR T
+//			cos_impl_2(T x) {
+//				return x <= sprout::math::quarter_pi<T>() ? -sprout::math::detail::cosp(x)
+//					: x <= sprout::math::half_pi<T>() ? -sprout::math::detail::sinp(sprout::math::half_pi<T>() - x)
+//					: x <= sprout::math::three_quarters_pi<T>() ? sprout::math::detail::sinp(x - sprout::math::half_pi<T>())
+//					: sprout::math::detail::cosp(sprout::math::pi<T>() - x)
+//					;
+//			}
+//			template<typename T>
+//			inline SPROUT_CONSTEXPR T
+//			cos_impl_1(T x) {
+//				return x > sprout::math::pi<T>() ? sprout::math::detail::cos_impl_2(x - sprout::math::pi<T>())
+//					: x <= sprout::math::quarter_pi<T>() ? sprout::math::detail::cosp(x)
+//					: x <= sprout::math::half_pi<T>() ? sprout::math::detail::sinp(sprout::math::half_pi<T>() - x)
+//					: x <= sprout::math::three_quarters_pi<T>() ? -sprout::math::detail::sinp(x - sprout::math::half_pi<T>())
+//					: -sprout::math::detail::cosp(sprout::math::pi<T>() - x)
+//					;
+//			}
+//			template<typename T>
+//			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<T>()));
+//			}
+
 			template<
 				typename FloatType,
 				typename sprout::enabler_if<std::is_floating_point<FloatType>::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 <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/isnan.hpp>
-#include <sprout/math/factorial.hpp>
+#include <sprout/math/exp.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
 namespace sprout {
@@ -17,20 +16,13 @@ namespace sprout {
 		namespace detail {
 			template<typename T>
 			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<T>(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<typename T>
 			inline SPROUT_CONSTEXPR T
 			cosh_impl(T x) {
-				return T(1) + sprout::math::detail::cosh_impl_1(
-					x * x,
-					1, sprout::math::factorial_limit<T>() / 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 <sprout/config.hpp>
+
+namespace sprout {
+	namespace math {
+		namespace detail {
+			template<typename T>
+			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<typename T>
+			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 <sprout/config.hpp>
+
+namespace sprout {
+	namespace math {
+		namespace detail {
+			template<typename T>
+			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<typename T>
+			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 <sprout/detail/pow.hpp>
 #include <sprout/math/detail/config.hpp>
 #include <sprout/math/detail/float_compute.hpp>
+#include <sprout/math/isnan.hpp>
 #include <sprout/math/factorial.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
@@ -37,7 +38,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			exp(FloatType x) {
-				return x == -std::numeric_limits<FloatType>::infinity() ? FloatType(0)
+				return sprout::math::isnan(x) ? x
+					: x == -std::numeric_limits<FloatType>::infinity() ? FloatType(0)
 					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::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 <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
 #include <sprout/math/detail/float_compute.hpp>
-#include <sprout/math/exp.hpp>
 #include <sprout/math/constants.hpp>
+#include <sprout/math/isnan.hpp>
+#include <sprout/math/exp.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
 namespace sprout {
@@ -25,7 +26,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			exp10(FloatType x) {
-				return x == -std::numeric_limits<FloatType>::infinity() ? FloatType(0)
+				return sprout::math::isnan(x) ? x
+					: x == -std::numeric_limits<FloatType>::infinity() ? FloatType(0)
 					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
 					: x == 0 ? FloatType(1)
 					: static_cast<FloatType>(sprout::math::detail::exp10_impl(static_cast<typename sprout::math::detail::float_compute<FloatType>::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 <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
 #include <sprout/math/detail/float_compute.hpp>
-#include <sprout/math/exp.hpp>
 #include <sprout/math/constants.hpp>
+#include <sprout/math/isnan.hpp>
+#include <sprout/math/exp.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
 namespace sprout {
@@ -25,7 +26,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			exp2(FloatType x) {
-				return x == -std::numeric_limits<FloatType>::infinity() ? FloatType(0)
+				return sprout::math::isnan(x) ? x
+					: x == -std::numeric_limits<FloatType>::infinity() ? FloatType(0)
 					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::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 <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
 #include <sprout/math/detail/float_compute.hpp>
+#include <sprout/math/constants.hpp>
+#include <sprout/math/isnan.hpp>
 #include <sprout/math/copysign.hpp>
 #include <sprout/math/exp.hpp>
-#include <sprout/math/constants.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
 namespace sprout {
@@ -26,7 +27,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			expm1(FloatType x) {
-				return x == -std::numeric_limits<FloatType>::infinity() ? FloatType(-1)
+				return sprout::math::isnan(x) ? x
+					: x == -std::numeric_limits<FloatType>::infinity() ? FloatType(-1)
 					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::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 <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/detail/pow.hpp>
+#include <sprout/math/isnan.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
 namespace sprout {
@@ -24,7 +25,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			ldexp(FloatType x, int exp) {
-				return x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+				return sprout::math::isnan(x) ? x
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
 					: x == -std::numeric_limits<FloatType>::infinity() ? -std::numeric_limits<FloatType>::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 <sprout/math/detail/config.hpp>
 #include <sprout/math/detail/float_compute.hpp>
 #include <sprout/math/constants.hpp>
+#include <sprout/math/isnan.hpp>
 #include <sprout/math/factorial.hpp>
 #include <sprout/math/sqrt.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
@@ -47,7 +48,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			log(FloatType x) {
-				return x == 0 ? -std::numeric_limits<FloatType>::infinity()
+				return sprout::math::isnan(x) ? x
+					: x == 0 ? -std::numeric_limits<FloatType>::infinity()
 					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
 					: x < 0 ? std::numeric_limits<FloatType>::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 <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
 #include <sprout/math/detail/float_compute.hpp>
-#include <sprout/math/log.hpp>
 #include <sprout/math/constants.hpp>
+#include <sprout/math/isnan.hpp>
+#include <sprout/math/log.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
 namespace sprout {
@@ -25,7 +26,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			log10(FloatType x) {
-				return x == 0 ? -std::numeric_limits<FloatType>::infinity()
+				return sprout::math::isnan(x) ? x
+					: x == 0 ? -std::numeric_limits<FloatType>::infinity()
 					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
 					: x < 0 ? std::numeric_limits<FloatType>::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 <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
 #include <sprout/math/detail/float_compute.hpp>
+#include <sprout/math/isnan.hpp>
 #include <sprout/math/log.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
@@ -24,9 +25,10 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			log1p(FloatType x) {
-				return x == -1 ? -std::numeric_limits<FloatType>::infinity()
+				return sprout::math::isnan(x) ? x
+					: x == -1 ? -std::numeric_limits<FloatType>::infinity()
 					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
-					: x < -1 ? std::numeric_limits<FloatType>::quiet_NaN()
+					: x < -1 ? -std::numeric_limits<FloatType>::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 <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
 #include <sprout/math/detail/float_compute.hpp>
-#include <sprout/math/log.hpp>
 #include <sprout/math/constants.hpp>
+#include <sprout/math/isnan.hpp>
+#include <sprout/math/log.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
 namespace sprout {
@@ -25,7 +26,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			log2(FloatType x) {
-				return x == 0 ? -std::numeric_limits<FloatType>::infinity()
+				return sprout::math::isnan(x) ? x
+					: x == 0 ? -std::numeric_limits<FloatType>::infinity()
 					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
 					: x < 0 ? std::numeric_limits<FloatType>::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 <sprout/detail/pow.hpp>
 #include <sprout/math/detail/config.hpp>
 #include <sprout/math/detail/float_compute.hpp>
+#include <sprout/math/isnan.hpp>
 #include <sprout/math/log_a.hpp>
 #include <sprout/math/trunc.hpp>
 #include <sprout/math/itrunc.hpp>
@@ -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<T>::radix, x / (x0 / std::numeric_limits<T>::radix), exp + 1
+						x, x0 * std::numeric_limits<T>::radix, x / (x0 / std::numeric_limits<T>::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<T>::radix) ? sprout::math::detail::logb_impl_3_pos_hi(
-						x, x0 / std::numeric_limits<T>::radix, x / (x0 * std::numeric_limits<T>::radix), exp - 1
+						x, x0 / std::numeric_limits<T>::radix, x / (x0 * std::numeric_limits<T>::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<T>::radix, x / (x0 / std::numeric_limits<T>::radix), exp + 1
+							x, x0 * std::numeric_limits<T>::radix, x / (x0 / std::numeric_limits<T>::radix), exp - 1
 							)
 						: !(base < std::numeric_limits<T>::radix) ? sprout::math::detail::logb_impl_3_pos_hi(
-								x, x0 / std::numeric_limits<T>::radix, x / (x0 * std::numeric_limits<T>::radix), exp - 1
+								x, x0 / std::numeric_limits<T>::radix, x / (x0 * std::numeric_limits<T>::radix), exp + 1
 								)
 						: exp
 						;
@@ -99,7 +100,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			logb(FloatType x) {
-				return x == 0 ? -std::numeric_limits<FloatType>::infinity()
+				return sprout::math::isnan(x) ? x
+					: x == 0 ? -std::numeric_limits<FloatType>::infinity()
 #if SPROUT_USE_BUILTIN_CMATH_FUNCTION
 #	if defined(__GNUC__)
 					: x == -std::numeric_limits<FloatType>::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 <limits>
 #	include <sprout/detail/pow.hpp>
 #	include <sprout/math/detail/float_compute.hpp>
+#	include <sprout/math/isnan.hpp>
 #	include <sprout/math/log_a.hpp>
 #	include <sprout/math/trunc.hpp>
 #	include <sprout/math/itrunc.hpp>
@@ -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<FloatType>::infinity()
+				return sprout::math::isnan(x) ? x
+					: x == 0 ? -std::numeric_limits<FloatType>::infinity()
 					: x == std::numeric_limits<FloatType>::infinity() || x == -std::numeric_limits<FloatType>::infinity()
 						? std::numeric_limits<FloatType>::infinity()
 					: static_cast<FloatType>(sprout::math::detail::logb2_impl(static_cast<typename sprout::math::detail::float_compute<FloatType>::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 <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/detail/pow.hpp>
+#include <sprout/math/isnan.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
 namespace sprout {
@@ -24,7 +25,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			scalbln(FloatType x, long exp) {
-				return x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+				return sprout::math::isnan(x) ? x
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
 					: x == -std::numeric_limits<FloatType>::infinity() ? -std::numeric_limits<FloatType>::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 <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/detail/pow.hpp>
+#include <sprout/math/isnan.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
 namespace sprout {
@@ -24,7 +25,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			scalbn(FloatType x, int exp) {
-				return x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+				return sprout::math::isnan(x) ? x
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
 					: x == -std::numeric_limits<FloatType>::infinity() ? -std::numeric_limits<FloatType>::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 <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/isnan.hpp>
-#include <sprout/math/factorial.hpp>
+#include <sprout/math/exp.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
 namespace sprout {
@@ -17,20 +16,13 @@ namespace sprout {
 		namespace detail {
 			template<typename T>
 			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<T>(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<typename T>
 			inline SPROUT_CONSTEXPR T
 			sinh_impl(T x) {
-				return x + sprout::math::detail::sinh_impl_1(
-					x,
-					1, (sprout::math::factorial_limit<T>() - 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 <sprout/math/detail/config.hpp>
 #include <sprout/math/detail/float_compute.hpp>
 #include <sprout/math/isnan.hpp>
-#include <sprout/math/sinh.hpp>
-#include <sprout/math/cosh.hpp>
+#include <sprout/math/exp.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
 namespace sprout {
 	namespace math {
 		namespace detail {
+			template<typename T>
+			inline SPROUT_CONSTEXPR T
+			tanh_impl_2(T t, T u) {
+				return (t - u) / (t + u);
+			}
+			template<typename T>
+			inline SPROUT_CONSTEXPR T
+			tanh_impl_1(T t) {
+				return sprout::math::detail::tanh_impl_2(t, T(1) / t);
+			}
 			template<typename T>
 			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<