diff --git a/sprout/cmath.hpp b/sprout/cmath.hpp
index 11a04603..a4a20dcd 100644
--- a/sprout/cmath.hpp
+++ b/sprout/cmath.hpp
@@ -7,6 +7,7 @@
 #include <sprout/math/hyperbolic.hpp>
 #include <sprout/math/exponential.hpp>
 #include <sprout/math/power.hpp>
+#include <sprout/math/error.hpp>
 #include <sprout/math/nearest.hpp>
 #include <sprout/math/remainders.hpp>
 #include <sprout/math/minmax.hpp>
diff --git a/sprout/compost/effects/auto_pan.hpp b/sprout/compost/effects/auto_pan.hpp
index 13ae262c..365a247d 100644
--- a/sprout/compost/effects/auto_pan.hpp
+++ b/sprout/compost/effects/auto_pan.hpp
@@ -33,7 +33,7 @@ namespace sprout {
 				typename std::iterator_traits<Outdirected>::value_type
 			>::type
 			calc(Outdirected const& x) const {
-				return (1 + depth_ * sprout::sin(2 * sprout::math::pi<Value>() * rate_ * x.index() / samples_per_sec_)) * *x;
+				return (1 + depth_ * sprout::sin(sprout::math::two_pi<Value>() * rate_ * x.index() / samples_per_sec_)) * *x;
 			}
 			template<bool Left, typename Outdirected>
 			SPROUT_CONSTEXPR typename std::enable_if<
@@ -41,7 +41,7 @@ namespace sprout {
 				typename std::iterator_traits<Outdirected>::value_type
 			>::type
 			calc(Outdirected const& x) const {
-				return (1 + depth_ * sprout::sin(2 * sprout::math::pi<Value>() * rate_ * x.index() / samples_per_sec_ + sprout::math::pi<Value>())) * *x;
+				return (1 + depth_ * sprout::sin(sprout::math::two_pi<Value>() * rate_ * x.index() / samples_per_sec_ + sprout::math::pi<Value>())) * *x;
 			}
 		public:
 			SPROUT_CONSTEXPR auto_pan_outdirected_value(
diff --git a/sprout/compost/effects/chorus.hpp b/sprout/compost/effects/chorus.hpp
index 6a51da74..e0acd5f5 100644
--- a/sprout/compost/effects/chorus.hpp
+++ b/sprout/compost/effects/chorus.hpp
@@ -57,7 +57,7 @@ namespace sprout {
 			template<typename Outdirected>
 			SPROUT_CONSTEXPR typename std::iterator_traits<Outdirected>::value_type
 			operator()(Outdirected const& x) const {
-				return calc(x, d_ + depth_ * sprout::math::sin(2 * sprout::math::pi<Value>() * rate_ * x.index() / samples_per_sec_));
+				return calc(x, d_ + depth_ * sprout::math::sin(sprout::math::two_pi<Value>() * rate_ * x.index() / samples_per_sec_));
 			}
 		};
 
diff --git a/sprout/compost/effects/tremolo.hpp b/sprout/compost/effects/tremolo.hpp
index 4220f5e8..3f37a04c 100644
--- a/sprout/compost/effects/tremolo.hpp
+++ b/sprout/compost/effects/tremolo.hpp
@@ -31,7 +31,7 @@ namespace sprout {
 			template<typename Outdirected>
 			SPROUT_CONSTEXPR typename std::iterator_traits<Outdirected>::value_type
 			operator()(Outdirected const& x) const {
-				return (1 + depth_ * sprout::math::sin(2 * sprout::math::pi<Value>() * rate_ * x.index() / samples_per_sec_)) * *x;
+				return (1 + depth_ * sprout::math::sin(sprout::math::two_pi<Value>() * rate_ * x.index() / samples_per_sec_)) * *x;
 			}
 		};
 
diff --git a/sprout/compost/effects/vibrato.hpp b/sprout/compost/effects/vibrato.hpp
index adda38e5..4f171b44 100644
--- a/sprout/compost/effects/vibrato.hpp
+++ b/sprout/compost/effects/vibrato.hpp
@@ -56,7 +56,7 @@ namespace sprout {
 			template<typename Outdirected>
 			SPROUT_CONSTEXPR typename std::iterator_traits<Outdirected>::value_type
 			operator()(Outdirected const& x) const {
-				return calc(x, d_ + depth_ * sprout::math::sin(2 * sprout::math::pi<Value>() * rate_ * x.index() / samples_per_sec_));
+				return calc(x, d_ + depth_ * sprout::math::sin(sprout::math::two_pi<Value>() * rate_ * x.index() / samples_per_sec_));
 			}
 		};
 
diff --git a/sprout/compost/utility/iir_filter.hpp b/sprout/compost/utility/iir_filter.hpp
index 4a352c8b..c1797e7a 100644
--- a/sprout/compost/utility/iir_filter.hpp
+++ b/sprout/compost/utility/iir_filter.hpp
@@ -18,7 +18,7 @@ namespace sprout {
 			iir_fc(T const& fc) {
 				typedef typename sprout::float_promote<T>::type type;
 				using sprout::tan;
-				return tan(sprout::math::pi<type>() * fc) / (2 * sprout::math::pi<type>());
+				return tan(sprout::math::pi<type>() * fc) / sprout::math::two_pi<type>();
 			}
 			template<typename T>
 			inline SPROUT_CONSTEXPR typename sprout::float_promote<T>::type
@@ -60,7 +60,7 @@ namespace sprout {
 			iir_lpf_impl(T const& fc, T const& q, A const& a, B const& b) {
 				return sprout::compost::detail::iir_lpf_impl_1<Result>(
 					fc, q, a, b,
-					2 * sprout::math::pi<T>() * fc
+					sprout::math::two_pi<T>() * fc
 					);
 			}
 		}	// namespace detail
@@ -112,7 +112,7 @@ namespace sprout {
 			iir_hpf_impl(T const& fc, T const& q, A const& a, B const& b) {
 				return sprout::compost::detail::iir_hpf_impl_1<Result>(
 					fc, q, a, b,
-					2 * sprout::math::pi<T>() * fc
+					sprout::math::two_pi<T>() * fc
 					);
 			}
 		}	// namespace detail
@@ -164,8 +164,8 @@ namespace sprout {
 			iir_bpf_impl(T const& fc1, T const& fc2, A const& a, B const& b) {
 				return sprout::compost::detail::iir_bpf_impl_1<Result>(
 					fc1, fc2, a, b,
-					2 * sprout::math::pi<T>() * (fc2 - fc1),
-					4 * sprout::math::pi<T>() * fc1 * fc2
+					sprout::math::two_pi<T>() * (fc2 - fc1),
+					sprout::math::four_pi<T>() * fc1 * fc2
 					);
 			}
 		}	// namespace detail
@@ -220,8 +220,8 @@ namespace sprout {
 			iir_bef_impl(T const& fc1, T const& fc2, A const& a, B const& b) {
 				return sprout::compost::detail::iir_bef_impl_1<Result>(
 					fc1, fc2, a, b,
-					2 * sprout::math::pi<T>() * (fc2 - fc1),
-					4 * sprout::math::pi<T>() * fc1 * fc2
+					sprout::math::two_pi<T>() * (fc2 - fc1),
+					sprout::math::four_pi<T>() * fc1 * fc2
 					);
 			}
 		}	// namespace detail
@@ -276,7 +276,7 @@ namespace sprout {
 			iir_resonator_impl(T const& fc, T const& q, A const& a, B const& b) {
 				return sprout::compost::detail::iir_resonator_impl_1<Result>(
 					fc, q, a, b,
-					2 * sprout::math::pi<T>() * fc
+					sprout::math::two_pi<T>() * fc
 					);
 			}
 		}	// namespace detail
@@ -328,7 +328,7 @@ namespace sprout {
 			iir_notch_impl(T const& fc, T const& q, A const& a, B const& b) {
 				return sprout::compost::detail::iir_notch_impl_1<Result>(
 					fc, q, a, b,
-					2 * sprout::math::pi<T>() * fc
+					sprout::math::two_pi<T>() * fc
 					);
 			}
 		}	// namespace detail
@@ -381,7 +381,7 @@ namespace sprout {
 			iir_low_shelving_impl(T const& fc, T const& q, T const& g, A const& a, B const& b) {
 				return sprout::compost::detail::iir_low_shelving_impl_1<Result>(
 					fc, q, a, b,
-					2 * sprout::math::pi<T>() * fc
+					sprout::math::two_pi<T>() * fc
 					);
 			}
 		}	// namespace detail
@@ -437,7 +437,7 @@ namespace sprout {
 			iir_high_shelving_impl(T const& fc, T const& q, T const& g, A const& a, B const& b) {
 				return sprout::compost::detail::iir_high_shelving_impl_1<Result>(
 					fc, q, a, b,
-					2 * sprout::math::pi<T>() * fc
+					sprout::math::two_pi<T>() * fc
 					);
 			}
 		}	// namespace detail
@@ -492,7 +492,7 @@ namespace sprout {
 			iir_peaking_impl(T const& fc, T const& q, T const& g, A const& a, B const& b) {
 				return sprout::compost::detail::iir_peaking_impl_1<Result>(
 					fc, q, a, b,
-					2 * sprout::math::pi<T>() * fc
+					sprout::math::two_pi<T>() * fc
 					);
 			}
 		}	// namespace detail
diff --git a/sprout/iterator/sinusoid_iterator.hpp b/sprout/iterator/sinusoid_iterator.hpp
index 78435cb0..7a4a5873 100644
--- a/sprout/iterator/sinusoid_iterator.hpp
+++ b/sprout/iterator/sinusoid_iterator.hpp
@@ -53,7 +53,7 @@ namespace sprout {
 			, frequency_(1)
 			, amplitude_(1)
 			, phase_(0)
-			, d_(value_type(2) * sprout::math::pi<value_type>())
+			, d_(sprout::math::two_pi<value_type>())
 		{}
 		sinusoid_iterator(sinusoid_iterator const&) = default;
 		explicit SPROUT_CONSTEXPR sinusoid_iterator(
@@ -66,7 +66,7 @@ namespace sprout {
 			, frequency_(frequency)
 			, amplitude_(amplitude)
 			, phase_(phase)
-			, d_(value_type(2) * sprout::math::pi<value_type>() * frequency)
+			, d_(sprout::math::two_pi<value_type>() * frequency)
 		{}
 		template<typename U>
 		SPROUT_CONSTEXPR sinusoid_iterator(sinusoid_iterator<U> const& it)
diff --git a/sprout/math/acos.hpp b/sprout/math/acos.hpp
index c0b7d5c4..d3cd21c6 100644
--- a/sprout/math/acos.hpp
+++ b/sprout/math/acos.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_ACOS_HPP
 #define SPROUT_MATH_ACOS_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -17,7 +18,10 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			acos(FloatType x) {
-				return sprout::math::half_pi<FloatType>() - sprout::math::asin(x);
+				return x == 1 ? FloatType(0)
+					: x > 1 || x < -1 ? std::numeric_limits<FloatType>::quiet_NaN()
+					: sprout::math::half_pi<FloatType>() - sprout::math::asin(x)
+					;
 			}
 
 			template<
diff --git a/sprout/math/acosh.hpp b/sprout/math/acosh.hpp
index 5f6996a2..9c46c502 100644
--- a/sprout/math/acosh.hpp
+++ b/sprout/math/acosh.hpp
@@ -18,7 +18,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			acosh(FloatType x) {
-				return x < 1 ? std::numeric_limits<FloatType>::quiet_NaN()
+				return x == 1 ? FloatType(0)
+					: x < 1 ? std::numeric_limits<FloatType>::quiet_NaN()
 					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
 					: sprout::math::log(x + sprout::math::sqrt(x * x - 1))
 					;
diff --git a/sprout/math/asin.hpp b/sprout/math/asin.hpp
index 2e88708f..7ceeb50e 100644
--- a/sprout/math/asin.hpp
+++ b/sprout/math/asin.hpp
@@ -7,6 +7,7 @@
 #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/constants.hpp>
 #include <sprout/math/factorial.hpp>
 #include <sprout/math/sqrt.hpp>
@@ -41,8 +42,9 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			asin(FloatType x) {
-				typedef double type;
-				return x > 1 || x < -1 ? std::numeric_limits<FloatType>::quiet_NaN()
+				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)))
 					;
diff --git a/sprout/math/asinh.hpp b/sprout/math/asinh.hpp
index e86d1367..abfb5d86 100644
--- a/sprout/math/asinh.hpp
+++ b/sprout/math/asinh.hpp
@@ -18,7 +18,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			asinh(FloatType x) {
-				return x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+				return x == 0 ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
 					: x == -std::numeric_limits<FloatType>::infinity() ? -std::numeric_limits<FloatType>::infinity()
 					: sprout::math::log(x + sprout::math::sqrt(x * x + 1))
 					;
diff --git a/sprout/math/atan.hpp b/sprout/math/atan.hpp
index ad4fa039..bea9e5ff 100644
--- a/sprout/math/atan.hpp
+++ b/sprout/math/atan.hpp
@@ -2,10 +2,12 @@
 #define SPROUT_MATH_ATAN_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/constants.hpp>
 #include <sprout/math/factorial.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
@@ -39,11 +41,15 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			atan(FloatType x) {
-				typedef double type;
-				return static_cast<FloatType>(
-					x < 0 ? -sprout::math::detail::atan_impl(static_cast<type>(-x))
-						: sprout::math::detail::atan_impl(static_cast<type>(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>()
+					: static_cast<FloatType>(
+						x < 0 ? -sprout::math::detail::atan_impl(static_cast<type>(-x))
+							: sprout::math::detail::atan_impl(static_cast<type>(x))
+						)
+					;
 			}
 
 			template<
diff --git a/sprout/math/atan2.hpp b/sprout/math/atan2.hpp
index 928da281..94ab68e9 100644
--- a/sprout/math/atan2.hpp
+++ b/sprout/math/atan2.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_ATAN2_HPP
 #define SPROUT_MATH_ATAN2_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -18,8 +19,22 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			atan2(FloatType y, FloatType x) {
-				return x < 0
-					? sprout::math::atan(y / x) + (y < 0 ? -1 : 1) * sprout::math::pi<FloatType>()
+				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()
+						? 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>()
+					: 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::math::atan(y / x) + (y < 0 ? -1 : 1) * sprout::math::pi<FloatType>()
 					: sprout::math::atan(y / x)
 					;
 			}
diff --git a/sprout/math/atanh.hpp b/sprout/math/atanh.hpp
index 28737df9..8e1da4fe 100644
--- a/sprout/math/atanh.hpp
+++ b/sprout/math/atanh.hpp
@@ -17,9 +17,10 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			atanh(FloatType x) {
-				return x < -1 || x > 1 ? std::numeric_limits<FloatType>::quiet_NaN()
-					: x == -1 ? -std::numeric_limits<FloatType>::infinity()
+				return x == 0 ? FloatType(0)
 					: x == 1 ? std::numeric_limits<FloatType>::infinity()
+					: x == -1 ? -std::numeric_limits<FloatType>::infinity()
+					: x > 1 || x < -1 ? std::numeric_limits<FloatType>::quiet_NaN()
 					: sprout::math::log((1 + x) / (1 - x)) / 2
 					;
 			}
diff --git a/sprout/math/cbrt.hpp b/sprout/math/cbrt.hpp
index 48963c54..d5493a04 100644
--- a/sprout/math/cbrt.hpp
+++ b/sprout/math/cbrt.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_CBRT_HPP
 #define SPROUT_MATH_CBRT_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -17,7 +18,10 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			cbrt(FloatType x) {
-				return x < 0 ? -sprout::math::pow(-x, sprout::math::third<FloatType>())
+				return x == 0 ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: x == -std::numeric_limits<FloatType>::infinity() ? -std::numeric_limits<FloatType>::infinity()
+					: x < 0 ? -sprout::math::pow(-x, sprout::math::third<FloatType>())
 					: sprout::math::pow(x, sprout::math::third<FloatType>())
 					;
 			}
diff --git a/sprout/math/ceil.hpp b/sprout/math/ceil.hpp
index 4cd9cd3d..2ee870c4 100644
--- a/sprout/math/ceil.hpp
+++ b/sprout/math/ceil.hpp
@@ -7,7 +7,6 @@
 #include <stdexcept>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
-#include <sprout/math/isinf.hpp>
 #include <sprout/math/equal_to.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
@@ -27,7 +26,9 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			ceil(FloatType x) {
-				return sprout::math::isinf(x) ? x
+				return x == 0 ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: x == -std::numeric_limits<FloatType>::infinity() ? -std::numeric_limits<FloatType>::infinity()
 					: std::numeric_limits<std::uintmax_t>::max() < x || std::numeric_limits<std::uintmax_t>::max() < -x
 						? SPROUT_MATH_THROW_LARGE_FLOAT_ROUNDING(std::domain_error("ceil: large float rounding."), x)
 					: x < 0 ? -static_cast<FloatType>(static_cast<std::uintmax_t>(-x))
diff --git a/sprout/math/constants.hpp b/sprout/math/constants.hpp
index 1a258c64..77af9ebb 100644
--- a/sprout/math/constants.hpp
+++ b/sprout/math/constants.hpp
@@ -7,20 +7,87 @@ namespace sprout {
 	namespace math {
 		//
 		// pi
-		// half_pi
-		// quarter_pi
 		//
 		template<typename T>
 		inline SPROUT_CONSTEXPR T pi() {
-			return 3.141592653589793238462643383279502884197169399375105820974944L;
+			return 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651L;
 		}
+		//
+		// half_pi
+		//
 		template<typename T>
 		inline SPROUT_CONSTEXPR T half_pi() {
-			return 1.570796326794896619231321691639751442098584699687552910487472L;
+			return 1.57079632679489661923132169163975144209858469968755291048747229615390820314310449931401741267105853399107404326L;
+		}
+		//
+		// third_pi
+		// two_thirds_pi
+		//
+		template<typename T>
+		inline SPROUT_CONSTEXPR T third_pi() {
+			return 1.04719755119659774615421446109316762806572313312503527365831486410260546876206966620934494178070568932738269550L;
 		}
 		template<typename T>
+		inline SPROUT_CONSTEXPR T two_thirds_pi() {
+			return 2.09439510239319549230842892218633525613144626625007054731662972820521093752413933241868988356141137865476539101L;
+		}
+		//
+		// quarter_pi
+		// three_quarters_pi
+		//
+		template<typename T>
 		inline SPROUT_CONSTEXPR T quarter_pi() {
-			return 0.785398163397448309615660845819875721049292349843776455243736L;
+			return 0.78539816339744830961566084581987572104929234984377645524373614807695410157155224965700870633552926699553702163L;
+		}
+		template<typename T>
+		inline SPROUT_CONSTEXPR T three_quarters_pi() {
+			return 2.35619449019234492884698253745962716314787704953132936573120844423086230471465674897102611900658780098661106488L;
+		}
+		//
+		// two_pi
+		// four_pi
+		//
+		template<typename T>
+		inline SPROUT_CONSTEXPR T two_pi() {
+			return 6.28318530717958647692528676655900576839433879875021164194988918461563281257241799725606965068423413596429617303L;
+		}
+		template<typename T>
+		inline SPROUT_CONSTEXPR T four_pi() {
+			return 12.56637061435917295385057353311801153678867759750042328389977836923126562514483599451213930136846827192859234606L;
+		}
+		//
+		// two_div_pi
+		// root_two_div_pi
+		//
+		template<typename T>
+		inline SPROUT_CONSTEXPR T two_div_pi() {
+			return 0.636619772367581343075535053490057448137838582961825794990669376235587190536906140360455211065012343824291370907L;
+		}
+		template<typename T>
+		inline SPROUT_CONSTEXPR T root_two_div_pi() {
+			return 0.797884560802865355879892119868763736951717262329869315331851659341315851798603677002504667814613872860605117725L;
+		}
+		//
+		// root_pi
+		// one_div_root_pi
+		// root_one_div_pi
+		// two_div_root_pi
+		//
+		template<typename T>
+		inline SPROUT_CONSTEXPR T root_pi() {
+			return 1.77245385090551602729816748334114518279754945612238712821380778985291128459103218137495065673854466541622682362L;
+		}
+		template<typename T>
+		inline SPROUT_CONSTEXPR T one_div_root_pi() {
+			return 0.564189583547756286948079451560772585844050629328998856844085721710642468441493414486743660202107363443028347906L;
+		}
+		template<typename T>
+		inline SPROUT_CONSTEXPR T root_one_div_pi() {
+			return 0.564189583547756286948079451560772585844050629328998856844085721710642468441493414486743660202107363443028347906L;
+		}
+		template<typename T>
+		inline SPROUT_CONSTEXPR T two_div_root_pi() {
+			return 1.128379167095512573896158903121545171688101258657997713688171443421284936882986828973487320404214726886056695812L;
 		}
 		//
 		// half
@@ -33,11 +100,11 @@ namespace sprout {
 		}
 		template<typename T>
 		inline SPROUT_CONSTEXPR T third() {
-			return 0.3333333333333333333333333333333333333333333333333333333333333333333333L;
+			return 0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333L;
 		}
 		template<typename T>
 		inline SPROUT_CONSTEXPR T twothirds() {
-			return 0.6666666666666666666666666666666666666666666666666666666666666666666666L;
+			return 0.66666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666L;
 		}
 		//
 		// root_two
@@ -45,18 +112,18 @@ namespace sprout {
 		//
 		template<typename T>
 		inline SPROUT_CONSTEXPR T root_two() {
-			return 1.414213562373095048801688724209698078569671875376948073L;
+			return 1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623L;
 		}
 		template<typename T>
 		inline SPROUT_CONSTEXPR T half_root_two() {
-			return 0.70710678118654752440084436210484903928483593756084L;
+			return 0.707106781186547524400844362104849039284835937688474036588339868995366239231053519425193767163820786367506923115L;
 		}
 		//
 		// e
 		//
 		template<typename T>
 		inline SPROUT_CONSTEXPR T e() {
-			return 2.7182818284590452353602874713526624977572470936999595749669676L;
+			return 2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193L;
 		}
 		//
 		// ln_ten
@@ -64,11 +131,11 @@ namespace sprout {
 		//
 		template<typename T>
 		inline SPROUT_CONSTEXPR T ln_ten() {
-			return 2.302585092994045684017991454684364207601101488628772976L;
+			return 2.30258509299404568401799145468436420760110148862877297603332790096757260967735248023599720508959829834196778404L;
 		}
 		template<typename T>
 		inline SPROUT_CONSTEXPR T ln_two() {
-			return 0.693147180559945309417232121458176568075500134360255254L;
+			return 0.693147180559945309417232121458176568075500134360255254120680009493393621969694715605863326996418687542001481021L;
 		}
 	}	// namespace math
 }	// namespace sprout
diff --git a/sprout/math/cos.hpp b/sprout/math/cos.hpp
index 001e4822..85ee35a7 100644
--- a/sprout/math/cos.hpp
+++ b/sprout/math/cos.hpp
@@ -7,6 +7,7 @@
 #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/fmod.hpp>
@@ -27,10 +28,10 @@ namespace sprout {
 			template<typename FloatType>
 			inline SPROUT_CONSTEXPR FloatType
 			cos_impl(FloatType x) {
-				typedef double type;
+				typedef typename sprout::math::detail::float_compute<FloatType>::type type;
 				return static_cast<FloatType>(
 					type(1) + sprout::math::detail::cos_impl_1(
-						static_cast<type>(x) * static_cast<type>(x),
+						sprout::detail::pow2(static_cast<type>(x)),
 						1, sprout::math::factorial_limit<type>() / 2 + 1
 						)
 					);
@@ -42,11 +43,10 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			cos(FloatType x) {
-				typedef double type;
-				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, 2 * sprout::math::pi<FloatType>()))
+				return x == 0 ? FloatType(1)
+					: 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>()))
 					;
 			}
 
diff --git a/sprout/math/cosh.hpp b/sprout/math/cosh.hpp
index c3c9d185..17564cb2 100644
--- a/sprout/math/cosh.hpp
+++ b/sprout/math/cosh.hpp
@@ -2,10 +2,12 @@
 #define SPROUT_MATH_COSH_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/type_traits/enabler_if.hpp>
 
@@ -28,13 +30,17 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			cosh(FloatType x) {
-				typedef double type;
-				return static_cast<FloatType>(
-					type(1) + sprout::math::detail::cosh_impl(
-						static_cast<type>(x) * static_cast<type>(x),
-						1, sprout::math::factorial_limit<type>() / 2 + 1
+				typedef typename sprout::math::detail::float_compute<FloatType>::type type;
+				return x == 0 ? FloatType(1)
+					: x == std::numeric_limits<FloatType>::infinity() || x == -std::numeric_limits<FloatType>::infinity()
+						? std::numeric_limits<FloatType>::infinity()
+					: static_cast<FloatType>(
+						type(1) + sprout::math::detail::cosh_impl(
+							static_cast<type>(x) * static_cast<type>(x),
+							1, sprout::math::factorial_limit<type>() / 2 + 1
+							)
 						)
-					);
+					;
 			}
 
 			template<
diff --git a/sprout/math/detail/float_compute.hpp b/sprout/math/detail/float_compute.hpp
new file mode 100644
index 00000000..1f444461
--- /dev/null
+++ b/sprout/math/detail/float_compute.hpp
@@ -0,0 +1,18 @@
+#ifndef SPROUT_MATH_DETAIL_FLOAT_COMPUTE_HPP
+#define SPROUT_MATH_DETAIL_FLOAT_COMPUTE_HPP
+
+#include <sprout/config.hpp>
+#include <sprout/type_traits/float_promote.hpp>
+
+namespace sprout {
+	namespace math {
+		namespace detail {
+			template<typename... Types>
+			struct float_compute
+				: public sprout::float_promote<double, Types...>
+			{};
+		}	// namespace detail
+	}	// namespace math
+}	// namespace sprout
+
+#endif	// #ifndef SPROUT_MATH_DETAIL_FLOAT_COMPUTE_HPP
diff --git a/sprout/math/erf.hpp b/sprout/math/erf.hpp
new file mode 100644
index 00000000..062dbcc2
--- /dev/null
+++ b/sprout/math/erf.hpp
@@ -0,0 +1,69 @@
+#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/type_traits/enabler_if.hpp>
+
+namespace sprout {
+	namespace math {
+		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)
+					;
+			}
+			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::math::exp(-x2 / (T(1) + a * x2) * (sprout::math::quarter_pi<T>() + a * x2));
+			}
+			template<
+				typename FloatType,
+				typename sprout::enabler_if<std::is_floating_point<FloatType>::value>::type = sprout::enabler
+			>
+			inline SPROUT_CONSTEXPR FloatType
+			erf(FloatType x) {
+				typedef typename sprout::math::detail::float_compute<FloatType>::type type;
+				return x == 0 ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? FloatType(1)
+					: x == -std::numeric_limits<FloatType>::infinity() ? FloatType(-1)
+					: static_cast<FloatType>(sprout::math::detail::erf_impl(static_cast<type>(x)))
+					;
+			}
+
+			template<
+				typename IntType,
+				typename sprout::enabler_if<std::is_integral<IntType>::value>::type = sprout::enabler
+			>
+			inline SPROUT_CONSTEXPR double
+			erf(IntType x) {
+				return sprout::math::detail::erf(static_cast<double>(x));
+			}
+		}	// namespace detail
+
+		using NS_SPROUT_MATH_DETAIL::erf;
+	}	// namespace math
+
+	using sprout::math::erf;
+}	// namespace sprout
+
+#endif	// #ifndef SPROUT_MATH_ERF_HPP
diff --git a/sprout/math/erfc.hpp b/sprout/math/erfc.hpp
new file mode 100644
index 00000000..a0d2f201
--- /dev/null
+++ b/sprout/math/erfc.hpp
@@ -0,0 +1,41 @@
+#ifndef SPROUT_MATH_ERFC_HPP
+#define SPROUT_MATH_ERFC_HPP
+
+#include <type_traits>
+#include <sprout/config.hpp>
+#include <sprout/math/detail/config.hpp>
+#include <sprout/math/erf.hpp>
+#include <sprout/type_traits/enabler_if.hpp>
+
+namespace sprout {
+	namespace math {
+		namespace detail {
+			template<
+				typename FloatType,
+				typename sprout::enabler_if<std::is_floating_point<FloatType>::value>::type = sprout::enabler
+			>
+			inline SPROUT_CONSTEXPR FloatType
+			erfc(FloatType x) {
+				return x == std::numeric_limits<FloatType>::infinity() ? FloatType(0)
+					: x == -std::numeric_limits<FloatType>::infinity() ? FloatType(2)
+					: FloatType(1) - sprout::math::erf(x)
+					;
+			}
+
+			template<
+				typename IntType,
+				typename sprout::enabler_if<std::is_integral<IntType>::value>::type = sprout::enabler
+			>
+			inline SPROUT_CONSTEXPR double
+			erfc(IntType x) {
+				return sprout::math::detail::erfc(static_cast<double>(x));
+			}
+		}	// namespace detail
+
+		using NS_SPROUT_MATH_DETAIL::erfc;
+	}	// namespace math
+
+	using sprout::math::erfc;
+}	// namespace sprout
+
+#endif	// #ifndef SPROUT_MATH_ERFC_HPP
diff --git a/sprout/math/error.hpp b/sprout/math/error.hpp
new file mode 100644
index 00000000..c930020c
--- /dev/null
+++ b/sprout/math/error.hpp
@@ -0,0 +1,8 @@
+#ifndef SPROUT_MATH_ERROR_HPP
+#define SPROUT_MATH_ERROR_HPP
+
+#include <sprout/config.hpp>
+#include <sprout/math/erf.hpp>
+#include <sprout/math/erfc.hpp>
+
+#endif	// #ifndef SPROUT_MATH_ERROR_HPP
diff --git a/sprout/math/exp.hpp b/sprout/math/exp.hpp
index 1468a38d..9d84c9b1 100644
--- a/sprout/math/exp.hpp
+++ b/sprout/math/exp.hpp
@@ -2,10 +2,12 @@
 #define SPROUT_MATH_EXP_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/type_traits/enabler_if.hpp>
 
@@ -24,7 +26,7 @@ namespace sprout {
 			template<typename T>
 			inline SPROUT_CONSTEXPR T
 			exp_impl(T x) {
-				return 1 + sprout::math::detail::exp_impl_1(x, 1, sprout::math::factorial_limit<T>() + 1);
+				return T(1) + sprout::math::detail::exp_impl_1(x, 1, sprout::math::factorial_limit<T>() / 2 + 1);
 			}
 
 			template<
@@ -33,8 +35,11 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			exp(FloatType x) {
-				typedef double type;
-				return !(x > -1) ? static_cast<FloatType>(1 / sprout::math::detail::exp_impl(-static_cast<type>(x)))
+				typedef typename sprout::math::detail::float_compute<FloatType>::type type;
+				return x == 0 ? FloatType(1)
+					: x == -std::numeric_limits<FloatType>::infinity() ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: !(x > -1) ? static_cast<FloatType>(type(1) / sprout::math::detail::exp_impl(-static_cast<type>(x)))
 					: static_cast<FloatType>(sprout::math::detail::exp_impl(static_cast<type>(x)))
 					;
 			}
diff --git a/sprout/math/exp10.hpp b/sprout/math/exp10.hpp
index e5a27d0d..ff0214d3 100644
--- a/sprout/math/exp10.hpp
+++ b/sprout/math/exp10.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_EXP10_HPP
 #define SPROUT_MATH_EXP10_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/exp.hpp>
@@ -16,7 +17,11 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			exp10(FloatType x) {
-				return sprout::math::exp(x * sprout::math::ln_ten<FloatType>());
+				return x == 0 ? FloatType(1)
+					: x == -std::numeric_limits<FloatType>::infinity() ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: sprout::math::exp(x * sprout::math::ln_ten<FloatType>())
+					;
 			}
 
 			template<
diff --git a/sprout/math/exp2.hpp b/sprout/math/exp2.hpp
index 72fbcd3e..794ec434 100644
--- a/sprout/math/exp2.hpp
+++ b/sprout/math/exp2.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_EXP2_HPP
 #define SPROUT_MATH_EXP2_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -17,7 +18,11 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			exp2(FloatType x) {
-				return sprout::math::exp(x * sprout::math::ln_two<FloatType>());
+				return x == 0 ? FloatType(1)
+					: x == -std::numeric_limits<FloatType>::infinity() ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: sprout::math::exp(x * sprout::math::ln_two<FloatType>())
+					;
 			}
 
 			template<
diff --git a/sprout/math/expm1.hpp b/sprout/math/expm1.hpp
index cccbde77..b68dc124 100644
--- a/sprout/math/expm1.hpp
+++ b/sprout/math/expm1.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_EXPM1_HPP
 #define SPROUT_MATH_EXPM1_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -17,7 +18,11 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			expm1(FloatType x) {
-				return sprout::math::exp(x) - 1;
+				return x == 0 ? FloatType(0)
+					: x == -std::numeric_limits<FloatType>::infinity() ? FloatType(-1)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: sprout::math::exp(x) - FloatType(1)
+					;
 			}
 
 			template<
diff --git a/sprout/math/exponential.hpp b/sprout/math/exponential.hpp
index 40628bd1..1f901acc 100644
--- a/sprout/math/exponential.hpp
+++ b/sprout/math/exponential.hpp
@@ -8,9 +8,8 @@
 #include <sprout/math/expm1.hpp>
 #include <sprout/math/log.hpp>
 #include <sprout/math/log10.hpp>
+#include <sprout/math/log1p.hpp>
 #include <sprout/math/log2.hpp>
-#include <sprout/math/log1p.hpp>
-#include <sprout/math/log1p.hpp>
 #include <sprout/math/floating_point.hpp>
 
 #endif	// #ifndef SPROUT_MATH_EXPONENTIAL_HPP
diff --git a/sprout/math/float2_exponent.hpp b/sprout/math/float2_exponent.hpp
index 9e84c3a8..70bbf5af 100644
--- a/sprout/math/float2_exponent.hpp
+++ b/sprout/math/float2_exponent.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_FLOAT2_EXPONENT_HPP
 #define SPROUT_MATH_FLOAT2_EXPONENT_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -16,7 +17,9 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR int
 			float2_exponent(FloatType x) {
-				return sprout::math::ilogb2(x) + 1;
+				return x == 0 ? 0
+					: sprout::math::ilogb2(x) + 1
+					;
 			}
 
 			template<
diff --git a/sprout/math/float2_sig_exp.hpp b/sprout/math/float2_sig_exp.hpp
index a79a16a4..a54317e7 100644
--- a/sprout/math/float2_sig_exp.hpp
+++ b/sprout/math/float2_sig_exp.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_FLOAT2_SIG_EXP_HPP
 #define SPROUT_MATH_FLOAT2_SIG_EXP_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -15,7 +16,13 @@ namespace sprout {
 			template<typename T>
 			inline SPROUT_CONSTEXPR sprout::pair<T, int>
 			float2_sig_exp_impl(T x, int exp) {
-				return sprout::pair<T, int>(x / sprout::detail::pow_n(T(2), exp), exp);
+				typedef sprout::pair<T, int> type;
+				return x == 0 ? type(T(0), exp)
+					: x == std::numeric_limits<T>::infinity() ? type(std::numeric_limits<T>::infinity(), exp)
+					: x == -std::numeric_limits<T>::infinity() ? type(-std::numeric_limits<T>::infinity(), exp)
+					: x == std::numeric_limits<T>::quiet_NaN() ? type(std::numeric_limits<T>::quiet_NaN(), exp)
+					: type(x / sprout::detail::pow_n(T(2), exp), exp)
+					;
 			}
 
 			template<
diff --git a/sprout/math/float2_significand.hpp b/sprout/math/float2_significand.hpp
index 5ae3ce77..5ca0f6ac 100644
--- a/sprout/math/float2_significand.hpp
+++ b/sprout/math/float2_significand.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_FLOAT2_SIGNIFICAND_HPP
 #define SPROUT_MATH_FLOAT2_SIGNIFICAND_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -17,7 +18,12 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			float2_significand(FloatType x) {
-				return x / sprout::detail::pow_n(FloatType(2), sprout::math::float2_exponent(x));
+				return x == 0 ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: x == -std::numeric_limits<FloatType>::infinity() ? -std::numeric_limits<FloatType>::infinity()
+					: x == std::numeric_limits<FloatType>::quiet_NaN() ? std::numeric_limits<FloatType>::quiet_NaN()
+					: x / sprout::detail::pow_n(FloatType(2), sprout::math::float2_exponent(x))
+					;
 			}
 
 			template<
diff --git a/sprout/math/float_exponent.hpp b/sprout/math/float_exponent.hpp
index 5dbe9455..661dac6e 100644
--- a/sprout/math/float_exponent.hpp
+++ b/sprout/math/float_exponent.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_FLOAT_EXPONENT_HPP
 #define SPROUT_MATH_FLOAT_EXPONENT_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -16,7 +17,9 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR int
 			float_exponent(FloatType x) {
-				return sprout::math::ilogb(x) + 1;
+				return x == 0 ? 0
+					: sprout::math::ilogb(x) + 1
+					;
 			}
 
 			template<
diff --git a/sprout/math/float_sig_exp.hpp b/sprout/math/float_sig_exp.hpp
index 425e5484..74fc1bdf 100644
--- a/sprout/math/float_sig_exp.hpp
+++ b/sprout/math/float_sig_exp.hpp
@@ -16,7 +16,13 @@ namespace sprout {
 			template<typename T>
 			inline SPROUT_CONSTEXPR sprout::pair<T, int>
 			float_sig_exp_impl(T x, int exp) {
-				return sprout::pair<T, int>(x / sprout::detail::pow_n(T(std::numeric_limits<T>::radix), exp), exp);
+				typedef sprout::pair<T, int> type;
+				return x == 0 ? type(T(0), exp)
+					: x == std::numeric_limits<T>::infinity() ? type(std::numeric_limits<T>::infinity(), exp)
+					: x == -std::numeric_limits<T>::infinity() ? type(-std::numeric_limits<T>::infinity(), exp)
+					: x == std::numeric_limits<T>::quiet_NaN() ? type(std::numeric_limits<T>::quiet_NaN(), exp)
+					: type(x / sprout::detail::pow_n(T(std::numeric_limits<T>::radix), exp), exp)
+					;
 			}
 
 			template<
diff --git a/sprout/math/float_significand.hpp b/sprout/math/float_significand.hpp
index c582b5b1..6cee2bca 100644
--- a/sprout/math/float_significand.hpp
+++ b/sprout/math/float_significand.hpp
@@ -18,7 +18,12 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			float_significand(FloatType x) {
-				return x / sprout::detail::pow_n(FloatType(std::numeric_limits<FloatType>::radix), sprout::math::float_exponent(x));
+				return x == 0 ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: x == -std::numeric_limits<FloatType>::infinity() ? -std::numeric_limits<FloatType>::infinity()
+					: x == std::numeric_limits<FloatType>::quiet_NaN() ? std::numeric_limits<FloatType>::quiet_NaN()
+					: x / sprout::detail::pow_n(FloatType(std::numeric_limits<FloatType>::radix), sprout::math::float_exponent(x))
+					;
 			}
 
 			template<
diff --git a/sprout/math/floor.hpp b/sprout/math/floor.hpp
index decda316..6f56ba14 100644
--- a/sprout/math/floor.hpp
+++ b/sprout/math/floor.hpp
@@ -7,7 +7,6 @@
 #include <stdexcept>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
-#include <sprout/math/isinf.hpp>
 #include <sprout/math/equal_to.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
@@ -27,7 +26,9 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			floor(FloatType x) {
-				return sprout::math::isinf(x) ? x
+				return x == 0 ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: x == -std::numeric_limits<FloatType>::infinity() ? -std::numeric_limits<FloatType>::infinity()
 					: std::numeric_limits<std::uintmax_t>::max() < x || std::numeric_limits<std::uintmax_t>::max() < -x
 						? SPROUT_MATH_THROW_LARGE_FLOAT_ROUNDING(std::domain_error("floor: large float rounding."), x)
 					: x < 0 ? sprout::math::detail::floor_impl(x, -static_cast<FloatType>(static_cast<std::uintmax_t>(-x)))
diff --git a/sprout/math/fmax.hpp b/sprout/math/fmax.hpp
index e7768aa1..3ec0bec5 100644
--- a/sprout/math/fmax.hpp
+++ b/sprout/math/fmax.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_FMAX_HPP
 #define SPROUT_MATH_FMAX_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -16,7 +17,7 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			fmax(FloatType x, FloatType y) {
-				return x < y ? y : x;
+				return x < y && !y == std::numeric_limits<FloatType>::quiet_NaN() ? y : x;
 			}
 
 			template<
diff --git a/sprout/math/fmin.hpp b/sprout/math/fmin.hpp
index a151bf81..4bf9cf4e 100644
--- a/sprout/math/fmin.hpp
+++ b/sprout/math/fmin.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_FMIN_HPP
 #define SPROUT_MATH_FMIN_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -16,7 +17,7 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			fmin(FloatType x, FloatType y) {
-				return x > y ? y : x;
+				return x > y && !y == std::numeric_limits<FloatType>::quiet_NaN() ? y : x;
 			}
 
 			template<
diff --git a/sprout/math/fmod.hpp b/sprout/math/fmod.hpp
index 93ce53d5..2b8fd467 100644
--- a/sprout/math/fmod.hpp
+++ b/sprout/math/fmod.hpp
@@ -1,7 +1,7 @@
 #ifndef SPROUT_MATH_FMOD_HPP
 #define SPROUT_MATH_FMOD_HPP
 
-#include <cstdint>
+#include <limits>
 #include <type_traits>
 #include <stdexcept>
 #include <sprout/config.hpp>
@@ -18,7 +18,10 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			fmod(FloatType x, FloatType y) {
-				return y == 0 ? throw std::domain_error("fmod: divide by zero.")
+				return x == std::numeric_limits<FloatType>::infinity() || x == -std::numeric_limits<FloatType>::infinity() || y == 0
+						? std::numeric_limits<FloatType>::quiet_NaN()
+					: x == 0 ? FloatType(0)
+					: y == std::numeric_limits<FloatType>::infinity() || y == -std::numeric_limits<FloatType>::infinity() ? x
 					: x - sprout::math::trunc(x / y) * y
 					;
 			}
diff --git a/sprout/math/frac_int.hpp b/sprout/math/frac_int.hpp
index b923d196..50eb1f49 100644
--- a/sprout/math/frac_int.hpp
+++ b/sprout/math/frac_int.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_FRAC_INT_HPP
 #define SPROUT_MATH_FRAC_INT_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -14,7 +15,11 @@ namespace sprout {
 			template<typename T>
 			inline SPROUT_CONSTEXPR sprout::pair<T, T>
 			frac_int_impl(T x, T ipart) {
-				return sprout::pair<T, T>(x - ipart, ipart);
+				typedef sprout::pair<T, T> type;
+				return  x == std::numeric_limits<T>::infinity() || x == -std::numeric_limits<T>::infinity() ? type(T(0), ipart)
+					: x == std::numeric_limits<T>::quiet_NaN() ? type(std::numeric_limits<T>::quiet_NaN(), ipart)
+					: type(x - ipart, ipart)
+					;
 			}
 
 			template<
diff --git a/sprout/math/fractional_part.hpp b/sprout/math/fractional_part.hpp
index f966e0e3..7a584476 100644
--- a/sprout/math/fractional_part.hpp
+++ b/sprout/math/fractional_part.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_FRACTIONAL_PART_HPP
 #define SPROUT_MATH_FRACTIONAL_PART_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -16,7 +17,10 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			fractional_part(FloatType x) {
-				return x - sprout::math::integer_part(x);
+				return x == std::numeric_limits<FloatType>::infinity() || x == -std::numeric_limits<FloatType>::infinity() ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::quiet_NaN() ? std::numeric_limits<FloatType>::quiet_NaN()
+					: x - sprout::math::integer_part(x)
+					;
 			}
 
 			template<
diff --git a/sprout/math/functions.hpp b/sprout/math/functions.hpp
index 011cbfca..baedfc21 100644
--- a/sprout/math/functions.hpp
+++ b/sprout/math/functions.hpp
@@ -7,6 +7,7 @@
 #include <sprout/math/hyperbolic.hpp>
 #include <sprout/math/exponential.hpp>
 #include <sprout/math/power.hpp>
+#include <sprout/math/error.hpp>
 #include <sprout/math/nearest.hpp>
 #include <sprout/math/remainders.hpp>
 #include <sprout/math/minmax.hpp>
diff --git a/sprout/math/hyperbolic.hpp b/sprout/math/hyperbolic.hpp
index a6d2df9c..31919310 100644
--- a/sprout/math/hyperbolic.hpp
+++ b/sprout/math/hyperbolic.hpp
@@ -2,11 +2,11 @@
 #define SPROUT_MATH_HYPERBOLIC_HPP
 
 #include <sprout/config.hpp>
-#include <sprout/math/sinh.hpp>
-#include <sprout/math/cosh.hpp>
-#include <sprout/math/tanh.hpp>
-#include <sprout/math/asinh.hpp>
 #include <sprout/math/acosh.hpp>
+#include <sprout/math/asinh.hpp>
 #include <sprout/math/atanh.hpp>
+#include <sprout/math/cosh.hpp>
+#include <sprout/math/sinh.hpp>
+#include <sprout/math/tanh.hpp>
 
 #endif	// #ifndef SPROUT_MATH_HYPERBOLIC_HPP
diff --git a/sprout/math/hypot.hpp b/sprout/math/hypot.hpp
index 5b2645bc..4334d580 100644
--- a/sprout/math/hypot.hpp
+++ b/sprout/math/hypot.hpp
@@ -1,11 +1,13 @@
 #ifndef SPROUT_MATH_HYPOT_HPP
 #define SPROUT_MATH_HYPOT_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/type_traits/float_promote.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 #include <sprout/math/detail/config.hpp>
+#include <sprout/math/fabs.hpp>
 #include <sprout/math/sqrt.hpp>
 
 namespace sprout {
@@ -17,7 +19,14 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			hypot(FloatType x, FloatType y) {
-				return sprout::math::sqrt(x * x + y * 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()
+						? std::numeric_limits<FloatType>::infinity()
+					: sprout::math::sqrt(x * x + y * y)
+					;
 			}
 
 			template<
diff --git a/sprout/math/iceil.hpp b/sprout/math/iceil.hpp
index 4ad04ef9..b0fa5d69 100644
--- a/sprout/math/iceil.hpp
+++ b/sprout/math/iceil.hpp
@@ -49,7 +49,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR To
 			iceil(FloatType x) {
-				return std::numeric_limits<To>::max() < x || std::numeric_limits<To>::min() > x
+				return x == 0 ? To(0)
+					: std::numeric_limits<To>::max() < x || std::numeric_limits<To>::min() > x
 						? SPROUT_MATH_THROW_LARGE_FLOAT_ROUNDING(std::domain_error("iceil: large float rounding."), static_cast<To>(x))
 					: sprout::math::detail::iceil_impl(x, static_cast<To>(x))
 					;
diff --git a/sprout/math/ifloor.hpp b/sprout/math/ifloor.hpp
index 964e6b06..45d6e0b3 100644
--- a/sprout/math/ifloor.hpp
+++ b/sprout/math/ifloor.hpp
@@ -49,7 +49,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR To
 			ifloor(FloatType x) {
-				return std::numeric_limits<To>::max() < x || std::numeric_limits<To>::min() > x
+				return x == 0 ? To(0)
+					: std::numeric_limits<To>::max() < x || std::numeric_limits<To>::min() > x
 						? SPROUT_MATH_THROW_LARGE_FLOAT_ROUNDING(std::domain_error("ifloor: large float rounding."), static_cast<To>(x))
 					: sprout::math::detail::ifloor_impl(x, static_cast<To>(x))
 					;
diff --git a/sprout/math/integer_part.hpp b/sprout/math/integer_part.hpp
index f8bf9962..338fa39e 100644
--- a/sprout/math/integer_part.hpp
+++ b/sprout/math/integer_part.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_INTEGER_PART_HPP
 #define SPROUT_MATH_INTEGER_PART_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -16,7 +17,11 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			integer_part(FloatType x) {
-				return sprout::math::trunc(x);
+				return x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: x == -std::numeric_limits<FloatType>::infinity() ? -std::numeric_limits<FloatType>::infinity()
+					: x == std::numeric_limits<FloatType>::quiet_NaN() ? std::numeric_limits<FloatType>::quiet_NaN()
+					: sprout::math::trunc(x)
+					;
 			}
 
 			template<
diff --git a/sprout/math/iround.hpp b/sprout/math/iround.hpp
index c7920da8..1ae6fa9e 100644
--- a/sprout/math/iround.hpp
+++ b/sprout/math/iround.hpp
@@ -56,7 +56,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR To
 			iround(FloatType x) {
-				return std::numeric_limits<To>::max() < x || std::numeric_limits<To>::min() > x
+				return x == 0 ? To(0)
+					: std::numeric_limits<To>::max() < x || std::numeric_limits<To>::min() > x
 						? SPROUT_MATH_THROW_LARGE_FLOAT_ROUNDING(std::domain_error("iround: large float irounding."), x)
 					: x < 0 ? sprout::math::detail::iround_impl_nagative(x, static_cast<To>(x))
 					: sprout::math::detail::iround_impl_positive(x, static_cast<To>(x))
diff --git a/sprout/math/itrunc.hpp b/sprout/math/itrunc.hpp
index d1f2ff30..9f1e44fb 100644
--- a/sprout/math/itrunc.hpp
+++ b/sprout/math/itrunc.hpp
@@ -40,7 +40,8 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR To
 			itrunc(FloatType x) {
-				return std::numeric_limits<To>::max() < x || std::numeric_limits<To>::min() > x
+				return x == 0 ? To(0)
+					: std::numeric_limits<To>::max() < x || std::numeric_limits<To>::min() > x
 						? SPROUT_MATH_THROW_LARGE_FLOAT_ROUNDING(std::domain_error("itrunc: large float rounding."), static_cast<To>(x))
 					: static_cast<To>(x)
 					;
diff --git a/sprout/math/ldexp.hpp b/sprout/math/ldexp.hpp
index 8029f30a..9804ee67 100644
--- a/sprout/math/ldexp.hpp
+++ b/sprout/math/ldexp.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_LDEXP_HPP
 #define SPROUT_MATH_LDEXP_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -16,7 +17,12 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			ldexp(FloatType x, int exp) {
-				return x * sprout::detail::pow_n(FloatType(2), exp);
+				return x == 0 ? FloatType(0)
+					: 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 * sprout::detail::pow_n(FloatType(2), exp)
+					;
 			}
 
 			template<
diff --git a/sprout/math/log.hpp b/sprout/math/log.hpp
index b265e7e2..4309d08c 100644
--- a/sprout/math/log.hpp
+++ b/sprout/math/log.hpp
@@ -7,6 +7,7 @@
 #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/constants.hpp>
 #include <sprout/math/factorial.hpp>
 #include <sprout/math/sqrt.hpp>
@@ -39,9 +40,11 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			log(FloatType x) {
-				typedef double type;
-				return x == 0 ? std::numeric_limits<FloatType>::quiet_NaN()
-					: !(x > 0) ? -std::numeric_limits<FloatType>::infinity()
+				typedef typename sprout::math::detail::float_compute<FloatType>::type type;
+				return x == 0 ? -std::numeric_limits<FloatType>::infinity()
+					: x == 1 ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: x < 0 ? std::numeric_limits<FloatType>::quiet_NaN()
 					: x < 1 ? static_cast<FloatType>(-sprout::math::detail::log_impl(1 / static_cast<type>(x)))
 					: static_cast<FloatType>(sprout::math::detail::log_impl(static_cast<type>(x)))
 					;
diff --git a/sprout/math/log10.hpp b/sprout/math/log10.hpp
index bfb98933..819ab27f 100644
--- a/sprout/math/log10.hpp
+++ b/sprout/math/log10.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_LOG10_HPP
 #define SPROUT_MATH_LOG10_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -17,7 +18,12 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			log10(FloatType x) {
-				return sprout::math::log(x) / sprout::math::ln_ten<FloatType>();
+				return x == 0 ? -std::numeric_limits<FloatType>::infinity()
+					: x == 1 ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: x < 0 ? std::numeric_limits<FloatType>::quiet_NaN()
+					: sprout::math::log(x) / sprout::math::ln_ten<FloatType>()
+					;
 			}
 
 			template<
diff --git a/sprout/math/log1p.hpp b/sprout/math/log1p.hpp
index 605f6dd9..c38525bc 100644
--- a/sprout/math/log1p.hpp
+++ b/sprout/math/log1p.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_LOG1P_HPP
 #define SPROUT_MATH_LOG1P_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -16,7 +17,12 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			log1p(FloatType x) {
-				return sprout::math::log(1 + x);
+				return x == 0 ? FloatType(0)
+					: 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()
+					: sprout::math::log(1 + x)
+					;
 			}
 
 			template<
diff --git a/sprout/math/log2.hpp b/sprout/math/log2.hpp
index 385e2921..6524bce8 100644
--- a/sprout/math/log2.hpp
+++ b/sprout/math/log2.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_LOG2_HPP
 #define SPROUT_MATH_LOG2_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -17,7 +18,12 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			log2(FloatType x) {
-				return sprout::math::log(x) / sprout::math::ln_two<FloatType>();
+				return x == 0 ? -std::numeric_limits<FloatType>::infinity()
+					: x == 1 ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: x < 0 ? std::numeric_limits<FloatType>::quiet_NaN()
+					: sprout::math::log(x) / sprout::math::ln_two<FloatType>()
+					;
 			}
 
 			template<
diff --git a/sprout/math/logb.hpp b/sprout/math/logb.hpp
index 6d59d515..a47619e4 100644
--- a/sprout/math/logb.hpp
+++ b/sprout/math/logb.hpp
@@ -90,12 +90,11 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			logb(FloatType x) {
-				return x < 0 ? sprout::math::detail::logb_impl(
-						-x, sprout::math::trunc(sprout::math::log_a(FloatType(std::numeric_limits<FloatType>::radix), -x))
-						)
-					: sprout::math::detail::logb_impl(
-						x, sprout::math::trunc(sprout::math::log_a(FloatType(std::numeric_limits<FloatType>::radix), x))
-						)
+				return x == 0 ? -std::numeric_limits<FloatType>::infinity()
+					: x == std::numeric_limits<FloatType>::infinity() || x == -std::numeric_limits<FloatType>::infinity()
+						? std::numeric_limits<FloatType>::infinity()
+					: x < 0 ? sprout::math::detail::logb_impl(-x, sprout::math::trunc(sprout::math::log_a(FloatType(std::numeric_limits<FloatType>::radix), -x)))
+					: sprout::math::detail::logb_impl(x, sprout::math::trunc(sprout::math::log_a(FloatType(std::numeric_limits<FloatType>::radix), x)))
 					;
 			}
 
diff --git a/sprout/math/logb2.hpp b/sprout/math/logb2.hpp
index 3f5ba26a..2648c975 100644
--- a/sprout/math/logb2.hpp
+++ b/sprout/math/logb2.hpp
@@ -10,6 +10,7 @@
 #	include <sprout/math/logb.hpp>
 #else
 #	include <cstdint>
+#	include <limits>
 #	include <sprout/detail/pow.hpp>
 #	include <sprout/math/log_a.hpp>
 #	include <sprout/math/trunc.hpp>
@@ -113,12 +114,11 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			logb2(FloatType x) {
-				return x < 0 ? sprout::math::detail::logb2_impl(
-						-x, sprout::math::trunc(sprout::math::log_a(FloatType(2), -x))
-						)
-					: sprout::math::detail::logb2_impl(
-						x, sprout::math::trunc(sprout::math::log_a(FloatType(2), x))
-						)
+				return x == 0 ? -std::numeric_limits<FloatType>::infinity()
+					: x == std::numeric_limits<FloatType>::infinity() || x == -std::numeric_limits<FloatType>::infinity()
+						? std::numeric_limits<FloatType>::infinity()
+					: x < 0 ? sprout::math::detail::logb2_impl(-x, sprout::math::trunc(sprout::math::log_a(FloatType(2), -x)))
+					: sprout::math::detail::logb2_impl(x, sprout::math::trunc(sprout::math::log_a(FloatType(2), x)))
 					;
 			}
 
diff --git a/sprout/math/nearest.hpp b/sprout/math/nearest.hpp
index 176b4b55..c2cd9c9a 100644
--- a/sprout/math/nearest.hpp
+++ b/sprout/math/nearest.hpp
@@ -4,12 +4,12 @@
 #include <sprout/config.hpp>
 #include <sprout/math/ceil.hpp>
 #include <sprout/math/floor.hpp>
-#include <sprout/math/trunc.hpp>
 #include <sprout/math/round.hpp>
+#include <sprout/math/trunc.hpp>
 #include <sprout/math/iceil.hpp>
 #include <sprout/math/ifloor.hpp>
-#include <sprout/math/itrunc.hpp>
 #include <sprout/math/iround.hpp>
+#include <sprout/math/itrunc.hpp>
 #include <sprout/math/lround.hpp>
 #include <sprout/math/llround.hpp>
 
diff --git a/sprout/math/pow.hpp b/sprout/math/pow.hpp
index 216871ac..415595ce 100644
--- a/sprout/math/pow.hpp
+++ b/sprout/math/pow.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_POW_HPP
 #define SPROUT_MATH_POW_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -19,7 +20,29 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			pow(FloatType x, FloatType y) {
-				return x == 0 && y > 0 ? FloatType(0)
+				return x == 0
+						? y < 0 ? std::numeric_limits<FloatType>::infinity()
+						: y > 0 ? FloatType(0)
+						: sprout::math::exp(y * sprout::math::log(x))
+					: x == -1 && (y == std::numeric_limits<FloatType>::infinity() || y == -std::numeric_limits<FloatType>::infinity()) ? FloatType(1)
+					: x == 1 ? FloatType(1)
+					: y == 0 ? FloatType(1)
+					: y == -std::numeric_limits<FloatType>::infinity()
+						? x < 1 && x > -1 ? std::numeric_limits<FloatType>::infinity()
+						: x > 1 || x < -1 ? FloatType(0)
+						: sprout::math::exp(y * sprout::math::log(x))
+					: y == std::numeric_limits<FloatType>::infinity()
+						? x < 1 && x > -1 ? FloatType(0)
+						: x > 1 || x < -1 ? std::numeric_limits<FloatType>::infinity()
+						: sprout::math::exp(y * sprout::math::log(x))
+					: x == -std::numeric_limits<FloatType>::infinity()
+						? y < 0 ? FloatType(0)
+						: y > 0 ? std::numeric_limits<FloatType>::infinity()
+						: sprout::math::exp(y * sprout::math::log(x))
+					: x == std::numeric_limits<FloatType>::infinity()
+						? y < 0 ? FloatType(0)
+						: y > 0 ? std::numeric_limits<FloatType>::infinity()
+						: sprout::math::exp(y * sprout::math::log(x))
 					: sprout::math::exp(y * sprout::math::log(x))
 					;
 			}
diff --git a/sprout/math/power.hpp b/sprout/math/power.hpp
index dfe9854d..31c4a427 100644
--- a/sprout/math/power.hpp
+++ b/sprout/math/power.hpp
@@ -2,12 +2,12 @@
 #define SPROUT_MATH_POWER_HPP
 
 #include <sprout/config.hpp>
-#include <sprout/math/sqrt.hpp>
 #include <sprout/math/cbrt.hpp>
-#include <sprout/math/pow.hpp>
-#include <sprout/math/log_a.hpp>
-#include <sprout/math/hypot.hpp>
-#include <sprout/math/abs.hpp>
 #include <sprout/math/fabs.hpp>
+#include <sprout/math/abs.hpp>
+#include <sprout/math/hypot.hpp>
+#include <sprout/math/pow.hpp>
+#include <sprout/math/sqrt.hpp>
+#include <sprout/math/log_a.hpp>
 
 #endif	// #ifndef SPROUT_MATH_POWER_HPP
diff --git a/sprout/math/quotient.hpp b/sprout/math/quotient.hpp
index 5cfe4358..8fb8db54 100644
--- a/sprout/math/quotient.hpp
+++ b/sprout/math/quotient.hpp
@@ -1,7 +1,7 @@
 #ifndef SPROUT_MATH_QUOTIENT_HPP
 #define SPROUT_MATH_QUOTIENT_HPP
 
-#include <cstdint>
+#include <limits>
 #include <type_traits>
 #include <stdexcept>
 #include <sprout/config.hpp>
@@ -20,7 +20,10 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR R
 			quotient(FloatType x, FloatType y) {
-				return y == 0 ? throw std::domain_error("quotient: divide by zero.")
+				return x == std::numeric_limits<FloatType>::infinity() || x == -std::numeric_limits<FloatType>::infinity() || y == 0
+						? std::numeric_limits<FloatType>::quiet_NaN()
+					: x == 0 ? FloatType(0)
+					: y == std::numeric_limits<FloatType>::infinity() || y == -std::numeric_limits<FloatType>::infinity() ? FloatType(0)
 					: sprout::math::iround<R>(x / y)
 					;
 			}
diff --git a/sprout/math/rem_quo.hpp b/sprout/math/rem_quo.hpp
index 63381ccd..37691d40 100644
--- a/sprout/math/rem_quo.hpp
+++ b/sprout/math/rem_quo.hpp
@@ -1,7 +1,7 @@
 #ifndef SPROUT_MATH_REM_QUO_HPP
 #define SPROUT_MATH_REM_QUO_HPP
 
-#include <cstdint>
+#include <limits>
 #include <type_traits>
 #include <stdexcept>
 #include <sprout/config.hpp>
@@ -17,7 +17,13 @@ namespace sprout {
 			template<typename R, typename T>
 			inline SPROUT_CONSTEXPR sprout::pair<T, R>
 			rem_quo_impl(T x, T y, R quo) {
-				return sprout::pair<T, T>(x - quo * y, quo);
+				typedef sprout::pair<T, R> type;
+				return x == std::numeric_limits<T>::infinity() || x == -std::numeric_limits<T>::infinity() || y == 0
+						? type(std::numeric_limits<T>::quiet_NaN(), quo)
+					: x == 0 ? type(T(0), quo)
+					: y == std::numeric_limits<T>::infinity() || y == -std::numeric_limits<T>::infinity() ? type(x, quo)
+					: type(x - quo * y, quo)
+					;
 			}
 			template<
 				typename R = int,
diff --git a/sprout/math/remainder.hpp b/sprout/math/remainder.hpp
index 56f177f4..af4fd035 100644
--- a/sprout/math/remainder.hpp
+++ b/sprout/math/remainder.hpp
@@ -1,7 +1,7 @@
 #ifndef SPROUT_MATH_REMAINDER_HPP
 #define SPROUT_MATH_REMAINDER_HPP
 
-#include <cstdint>
+#include <limits>
 #include <type_traits>
 #include <stdexcept>
 #include <sprout/config.hpp>
@@ -19,7 +19,10 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			remainder(FloatType x, FloatType y) {
-				return y == 0 ? throw std::domain_error("remainder: divide by zero.")
+				return x == std::numeric_limits<FloatType>::infinity() || x == -std::numeric_limits<FloatType>::infinity() || y == 0
+						? std::numeric_limits<FloatType>::quiet_NaN()
+					: x == 0 ? FloatType(0)
+					: y == std::numeric_limits<FloatType>::infinity() || y == -std::numeric_limits<FloatType>::infinity() ? x
 					: x - sprout::math::round(x / y) * y
 					;
 			}
diff --git a/sprout/math/round.hpp b/sprout/math/round.hpp
index f2f76630..ea539f19 100644
--- a/sprout/math/round.hpp
+++ b/sprout/math/round.hpp
@@ -7,7 +7,6 @@
 #include <stdexcept>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
-#include <sprout/math/isinf.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
 namespace sprout {
@@ -33,7 +32,9 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			round(FloatType x) {
-				return sprout::math::isinf(x) ? x
+				return x == 0 ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: x == -std::numeric_limits<FloatType>::infinity() ? -std::numeric_limits<FloatType>::infinity()
 					: std::numeric_limits<std::uintmax_t>::max() < x || std::numeric_limits<std::uintmax_t>::max() < -x
 						? SPROUT_MATH_THROW_LARGE_FLOAT_ROUNDING(std::domain_error("round: large float rounding."), x)
 					: x < 0 ? sprout::math::detail::round_impl_nagative(x, -static_cast<FloatType>(static_cast<std::uintmax_t>(-x)))
diff --git a/sprout/math/scalbln.hpp b/sprout/math/scalbln.hpp
index 3a86f06f..a21e63a6 100644
--- a/sprout/math/scalbln.hpp
+++ b/sprout/math/scalbln.hpp
@@ -17,7 +17,12 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			scalbln(FloatType x, long exp) {
-				return x * sprout::detail::pow_n(FloatType(std::numeric_limits<FloatType>::radix), exp);
+				return x == 0 ? FloatType(0)
+					: 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 * sprout::detail::pow_n(FloatType(std::numeric_limits<FloatType>::radix), exp)
+					;
 			}
 
 			template<
diff --git a/sprout/math/scalbn.hpp b/sprout/math/scalbn.hpp
index 4290e99f..f8f7f9b8 100644
--- a/sprout/math/scalbn.hpp
+++ b/sprout/math/scalbn.hpp
@@ -17,7 +17,12 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			scalbn(FloatType x, int exp) {
-				return x * sprout::detail::pow_n(FloatType(std::numeric_limits<FloatType>::radix), exp);
+				return x == 0 ? FloatType(0)
+					: 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 * sprout::detail::pow_n(FloatType(std::numeric_limits<FloatType>::radix), exp)
+					;
 			}
 
 			template<
diff --git a/sprout/math/sin.hpp b/sprout/math/sin.hpp
index e145807f..9f5003ad 100644
--- a/sprout/math/sin.hpp
+++ b/sprout/math/sin.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_SIN_HPP
 #define SPROUT_MATH_SIN_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -17,7 +18,11 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			sin(FloatType x) {
-				return -sprout::math::cos(x + sprout::math::half_pi<FloatType>());
+				return x == 0 ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() || x == -std::numeric_limits<FloatType>::infinity()
+						? std::numeric_limits<FloatType>::quiet_NaN()
+					: -sprout::math::cos(x + sprout::math::half_pi<FloatType>())
+					;
 			}
 
 			template<
diff --git a/sprout/math/sinh.hpp b/sprout/math/sinh.hpp
index 29211b66..3248a87f 100644
--- a/sprout/math/sinh.hpp
+++ b/sprout/math/sinh.hpp
@@ -2,10 +2,12 @@
 #define SPROUT_MATH_SINH_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/type_traits/enabler_if.hpp>
 
@@ -28,13 +30,17 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			sinh(FloatType x) {
-				typedef double type;
-				return static_cast<FloatType>(
-					static_cast<type>(x) + sprout::math::detail::sinh_impl(
-						static_cast<type>(x),
-						1, (sprout::math::factorial_limit<type>() - 1) / 2 + 1
+				typedef typename sprout::math::detail::float_compute<FloatType>::type type;
+				return x == 0 ? FloatType(1)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: x == -std::numeric_limits<FloatType>::infinity() ? -std::numeric_limits<FloatType>::infinity()
+					: static_cast<FloatType>(
+						static_cast<type>(x) + sprout::math::detail::sinh_impl(
+							static_cast<type>(x),
+							1, (sprout::math::factorial_limit<type>() - 1) / 2 + 1
+							)
 						)
-					);
+					;
 			}
 
 			template<
diff --git a/sprout/math/sqrt.hpp b/sprout/math/sqrt.hpp
index 1d67a91e..ed4d50c4 100644
--- a/sprout/math/sqrt.hpp
+++ b/sprout/math/sqrt.hpp
@@ -5,6 +5,7 @@
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
+#include <sprout/math/detail/float_compute.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
 namespace sprout {
@@ -14,22 +15,13 @@ namespace sprout {
 			inline SPROUT_CONSTEXPR T
 			sqrt_impl_1(T x, T s, T s2) {
 				return !(s < s2) ? s2
-					: sprout::math::detail::sqrt_impl_1(
-						x,
-						(x / s + s) / 2,
-						s
-						)
+					: sprout::math::detail::sqrt_impl_1(x, (x / s + s) / 2, s)
 					;
 			}
 			template<typename T>
 			inline SPROUT_CONSTEXPR T
 			sqrt_impl(T x, T s) {
-				return sprout::math::detail::sqrt_impl_1(
-						x,
-						(x / s + s) / 2,
-						s
-						)
-					;
+				return sprout::math::detail::sqrt_impl_1(x, (x / s + s) / 2, s);
 			}
 
 			template<
@@ -38,13 +30,12 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			sqrt(FloatType x) {
-				typedef double type;
-				return x < 0 ? std::numeric_limits<FloatType>::quiet_NaN()
-					: x == 0 ? type(0)
-					: static_cast<FloatType>(sprout::math::detail::sqrt_impl(
-						static_cast<type>(x),
-						x > 1 ? static_cast<type>(x) : type(1)
-						));
+				typedef typename sprout::math::detail::float_compute<FloatType>::type type;
+				return x == 0 ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: x == std::numeric_limits<FloatType>::quiet_NaN() ? std::numeric_limits<FloatType>::quiet_NaN()
+					: x < 0 ? std::numeric_limits<FloatType>::quiet_NaN()
+					: static_cast<FloatType>(sprout::math::detail::sqrt_impl(static_cast<type>(x), x > 1 ? static_cast<type>(x) : type(1)));
 			}
 
 			template<
diff --git a/sprout/math/tan.hpp b/sprout/math/tan.hpp
index 198be179..c23bcd4a 100644
--- a/sprout/math/tan.hpp
+++ b/sprout/math/tan.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_TAN_HPP
 #define SPROUT_MATH_TAN_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -17,7 +18,11 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			tan(FloatType x) {
-				return sprout::math::sin(x) / sprout::math::cos(x);
+				return x == 0 ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() || x == -std::numeric_limits<FloatType>::infinity()
+						? std::numeric_limits<FloatType>::quiet_NaN()
+					: sprout::math::sin(x) / sprout::math::cos(x)
+					;
 			}
 
 			template<
diff --git a/sprout/math/tanh.hpp b/sprout/math/tanh.hpp
index 918e7975..007b2225 100644
--- a/sprout/math/tanh.hpp
+++ b/sprout/math/tanh.hpp
@@ -1,6 +1,7 @@
 #ifndef SPROUT_MATH_TANH_HPP
 #define SPROUT_MATH_TANH_HPP
 
+#include <limits>
 #include <type_traits>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
@@ -17,7 +18,11 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			tanh(FloatType x) {
-				return sprout::math::sinh(x) / sprout::math::cosh(x);
+				return x == 0 ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? FloatType(1)
+					: x == -std::numeric_limits<FloatType>::infinity() ? FloatType(-1)
+					: sprout::math::sinh(x) / sprout::math::cosh(x)
+					;
 			}
 
 			template<
diff --git a/sprout/math/trigonometric.hpp b/sprout/math/trigonometric.hpp
index e6fa6cac..ee821771 100644
--- a/sprout/math/trigonometric.hpp
+++ b/sprout/math/trigonometric.hpp
@@ -2,12 +2,12 @@
 #define SPROUT_MATH_TRIGONOMETRIC_HPP
 
 #include <sprout/config.hpp>
-#include <sprout/math/sin.hpp>
-#include <sprout/math/cos.hpp>
-#include <sprout/math/tan.hpp>
-#include <sprout/math/asin.hpp>
 #include <sprout/math/acos.hpp>
+#include <sprout/math/asin.hpp>
 #include <sprout/math/atan.hpp>
 #include <sprout/math/atan2.hpp>
+#include <sprout/math/cos.hpp>
+#include <sprout/math/sin.hpp>
+#include <sprout/math/tan.hpp>
 
 #endif	// #ifndef SPROUT_MATH_TRIGONOMETRIC_HPP
diff --git a/sprout/math/trunc.hpp b/sprout/math/trunc.hpp
index 14b1b31e..8cfce95c 100644
--- a/sprout/math/trunc.hpp
+++ b/sprout/math/trunc.hpp
@@ -7,7 +7,6 @@
 #include <stdexcept>
 #include <sprout/config.hpp>
 #include <sprout/math/detail/config.hpp>
-#include <sprout/math/isinf.hpp>
 #include <sprout/type_traits/enabler_if.hpp>
 
 namespace sprout {
@@ -19,7 +18,9 @@ namespace sprout {
 			>
 			inline SPROUT_CONSTEXPR FloatType
 			trunc(FloatType x) {
-				return sprout::math::isinf(x) ? x
+				return x == 0 ? FloatType(0)
+					: x == std::numeric_limits<FloatType>::infinity() ? std::numeric_limits<FloatType>::infinity()
+					: x == -std::numeric_limits<FloatType>::infinity() ? -std::numeric_limits<FloatType>::infinity()
 					: std::numeric_limits<std::uintmax_t>::max() < x || std::numeric_limits<std::uintmax_t>::max() < -x
 						? SPROUT_MATH_THROW_LARGE_FLOAT_ROUNDING(std::domain_error("trunc: large float rounding."), x)
 					: x < 0 ? -static_cast<FloatType>(static_cast<std::uintmax_t>(-x))
diff --git a/sprout/numeric/dft/dft_element.hpp b/sprout/numeric/dft/dft_element.hpp
index 702182e2..24b5ee57 100644
--- a/sprout/numeric/dft/dft_element.hpp
+++ b/sprout/numeric/dft/dft_element.hpp
@@ -21,7 +21,7 @@ namespace sprout {
 			return sprout::detail::dft_element_gen(
 				first,
 				last,
-				-(2 * sprout::math::pi<elem_type>() * i / size)
+				-(sprout::math::two_pi<elem_type>() * i / size)
 				);
 		}
 	}	// namespace detail
diff --git a/sprout/numeric/dft/fixed/sinusoid.hpp b/sprout/numeric/dft/fixed/sinusoid.hpp
index 058c1082..92ab2901 100644
--- a/sprout/numeric/dft/fixed/sinusoid.hpp
+++ b/sprout/numeric/dft/fixed/sinusoid.hpp
@@ -50,7 +50,7 @@ namespace sprout {
 				typedef typename sprout::container_traits<Container>::value_type value_type;
 				return sprout::fixed::detail::sinusoid_impl(
 					cont,
-					value_type(2) * sprout::math::pi<value_type>() * frequency / value_type(sprout::size(cont)),
+					sprout::math::two_pi<value_type>() * frequency / value_type(sprout::size(cont)),
 					amplitude,
 					phase,
 					sprout::container_indexes<Container>::make(),
diff --git a/sprout/numeric/dft/fixed/triangle.hpp b/sprout/numeric/dft/fixed/triangle.hpp
index 84f3f3d7..bb84b98c 100644
--- a/sprout/numeric/dft/fixed/triangle.hpp
+++ b/sprout/numeric/dft/fixed/triangle.hpp
@@ -22,7 +22,7 @@ namespace sprout {
 			triangle_value(T const& t) {
 				using sprout::sin;
 				using sprout::asin;
-				return T(2) / T(sprout::math::pi<double>()) * asin(sin(T(2) * T(sprout::math::pi<double>()) * t));
+				return T(sprout::math::two_div_pi<double>()) * asin(sin(T(sprout::math::two_pi<double>()) * t));
 			}
 
 			template<typename Container, sprout::index_t... Indexes>
diff --git a/sprout/numeric/dft/idft_element.hpp b/sprout/numeric/dft/idft_element.hpp
index 9972a2d9..6acfa299 100644
--- a/sprout/numeric/dft/idft_element.hpp
+++ b/sprout/numeric/dft/idft_element.hpp
@@ -21,7 +21,7 @@ namespace sprout {
 			return sprout::detail::dft_element_gen(
 				first,
 				last,
-				2 * sprout::math::pi<elem_type>() * i / size
+				sprout::math::two_pi<elem_type>() * i / size
 				)
 				/ static_cast<elem_type>(size)
 				;
diff --git a/sprout/random/normal_distribution.hpp b/sprout/random/normal_distribution.hpp
index eb5150e7..3f95fe9e 100644
--- a/sprout/random/normal_distribution.hpp
+++ b/sprout/random/normal_distribution.hpp
@@ -127,8 +127,8 @@ namespace sprout {
 				return sprout::random::random_result<Engine, normal_distribution>(
 					cached_rho
 						* (valid
-							? sprout::cos(result_type(2) * sprout::math::pi<result_type>() * r1)
-							: sprout::sin(result_type(2) * sprout::math::pi<result_type>() * r1)
+							? sprout::cos(sprout::math::two_pi<result_type>() * r1)
+							: sprout::sin(sprout::math::two_pi<result_type>() * r1)
 							)
 						* sigma_ + mean_,
 					eng,