From a27c83e9390075af5da3885d470e1a9594a913a8 Mon Sep 17 00:00:00 2001
From: bolero-MURAKAMI <bolero.murakami@gmail.com>
Date: Tue, 23 Apr 2013 23:07:20 +0900
Subject: [PATCH] fix erf, erfc fix for special values: trigonometric functions

---
 sprout/math/acos.hpp      |  33 +++++++++--
 sprout/math/asin.hpp      |  35 ++++++++----
 sprout/math/atan.hpp      |  24 ++++++--
 sprout/math/atan2.hpp     |  60 +++++++++++++++-----
 sprout/math/cos.hpp       |  41 +++++++++-----
 sprout/math/erf.hpp       | 116 ++++++++++++++++++++++++++++++++------
 sprout/math/factorial.hpp |  12 +++-
 sprout/math/hypot.hpp     |   4 +-
 sprout/math/lgamma.hpp    |   4 +-
 sprout/math/sin.hpp       |  35 +++++++++---
 sprout/math/tan.hpp       |  34 ++++++++---
 sprout/math/tanh.hpp      |   2 +-
 12 files changed, 313 insertions(+), 87 deletions(-)

diff --git a/sprout/math/acos.hpp b/sprout/math/acos.hpp
index 511134db..2342e8c3 100644
--- a/sprout/math/acos.hpp
+++ b/sprout/math/acos.hpp
@@ -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;
diff --git a/sprout/math/asin.hpp b/sprout/math/asin.hpp
index 7ceeb50e..69242657 100644
--- a/sprout/math/asin.hpp
+++ b/sprout/math/asin.hpp
@@ -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;
diff --git a/sprout/math/atan.hpp b/sprout/math/atan.hpp
index bea9e5ff..57230652 100644
--- a/sprout/math/atan.hpp
+++ b/sprout/math/atan.hpp
@@ -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;
diff --git a/sprout/math/atan2.hpp b/sprout/math/atan2.hpp
index 6a5215ff..e64ae0aa 100644
--- a/sprout/math/atan2.hpp
+++ b/sprout/math/atan2.hpp
@@ -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;
diff --git a/sprout/math/cos.hpp b/sprout/math/cos.hpp
index 85ee35a7..d1eb19d5 100644
--- a/sprout/math/cos.hpp
+++ b/sprout/math/cos.hpp
@@ -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;
diff --git a/sprout/math/erf.hpp b/sprout/math/erf.hpp
index 1d387cd0..04c4c0bf 100644
--- a/sprout/math/erf.hpp
+++ b/sprout/math/erf.hpp
@@ -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
diff --git a/sprout/math/factorial.hpp b/sprout/math/factorial.hpp
index fc664e91..8a0f7266 100644
--- a/sprout/math/factorial.hpp
+++ b/sprout/math/factorial.hpp
@@ -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
 
diff --git a/sprout/math/hypot.hpp b/sprout/math/hypot.hpp
index a6c90c51..950a89ce 100644
--- a/sprout/math/hypot.hpp
+++ b/sprout/math/hypot.hpp
@@ -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()
diff --git a/sprout/math/lgamma.hpp b/sprout/math/lgamma.hpp
index ebacd50d..6ed5c5e6 100644
--- a/sprout/math/lgamma.hpp
+++ b/sprout/math/lgamma.hpp
@@ -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
 					;
 			}
diff --git a/sprout/math/sin.hpp b/sprout/math/sin.hpp
index 6bfbcdf5..7ca86499 100644
--- a/sprout/math/sin.hpp
+++ b/sprout/math/sin.hpp
@@ -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;
diff --git a/sprout/math/tan.hpp b/sprout/math/tan.hpp
index f744e6a7..2e393059 100644
--- a/sprout/math/tan.hpp
+++ b/sprout/math/tan.hpp
@@ -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;
diff --git a/sprout/math/tanh.hpp b/sprout/math/tanh.hpp
index 238cc3f8..007b2225 100644
--- a/sprout/math/tanh.hpp
+++ b/sprout/math/tanh.hpp
@@ -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)
 					;
 			}