Add Simpson's 1/3 rule quadrature

src/math/quadrature.hpp

@@ -28,6 +28,19 @@ namespace math {
 /// Numerical integration functions.
 namespace quadrature {
 
+/**
+ * Approximates the definite integral of a function using Simpson's 1/3 rule.
+ *
+ * @param f Unary function object to integrate.
+ * @param first,last Range of sample points on `[first, last)`.
+ * @return Approximated integral of @p f.
+ *
+ * @see https://en.wikipedia.org/wiki/Simpson%27s_rule
+ */
+template<class UnaryOp, class InputIt>
+typename std::invoke_result<UnaryOp, typename std::iterator_traits<InputIt>::value_type>::type
+	simpson(UnaryOp f, InputIt first, InputIt last);
+
 /**
  * Approximates the definite integral of a function using the trapezoidal rule.
  *
@@ -38,14 +51,59 @@ namespace quadrature {
  * @see https://en.wikipedia.org/wiki/Trapezoidal_rule
  */
 template<class UnaryOp, class InputIt>
-typename std::invoke_result<UnaryOp, typename std::iterator_traits<InputIt>::value_type>::type trapezoid(UnaryOp f, InputIt first, InputIt last)
+typename std::invoke_result<UnaryOp, typename std::iterator_traits<InputIt>::value_type>::type
+	trapezoid(UnaryOp f, InputIt first, InputIt last);
+
+template<class UnaryOp, class InputIt>
+typename std::invoke_result<UnaryOp, typename std::iterator_traits<InputIt>::value_type>::type
+	simpson(UnaryOp f, InputIt first, InputIt last)
+{
+	typedef typename std::iterator_traits<InputIt>::value_type input_type;
+	typedef typename std::invoke_result<UnaryOp, input_type>::type output_type;
+	
+	if (first == last)
+		return output_type(0);
+	
+	output_type f_a = f(*first);
+	
+	InputIt second = first;
+	++second;
+	
+	if (second == last)
+		return f_a;
+	
+	output_type f_c = f(*second);
+	
+	input_type h = *second - *first;
+	output_type f_b = f(*first + h / input_type(2));
+	
+	output_type sum = ((f_a + f_b * input_type(4) + f_c) / input_type(6)) * h;
+	
+	for (first = second++; second != last; first = second++)
+	{
+		h = *second - *first;
+		
+		f_a = f_c;
+		f_c = f(*second);
+		f_b = f(*first + h / input_type(2));
+		
+		sum += ((f_a + f_b * input_type(4) + f_c) / input_type(6)) * h;
+	}
+	
+	return sum;
+}
+
+template<class UnaryOp, class InputIt>
+typename std::invoke_result<UnaryOp, typename std::iterator_traits<InputIt>::value_type>::type
+	trapezoid(UnaryOp f, InputIt first, InputIt last)
 {
-	typedef typename std::invoke_result<UnaryOp, typename std::iterator_traits<InputIt>::value_type>::type result_type;
+	typedef typename std::iterator_traits<InputIt>::value_type input_type;
+	typedef typename std::invoke_result<UnaryOp, input_type>::type output_type;
 	
 	if (first == last)
-		return result_type(0);
+		return output_type(0);
 	
-	result_type f_a = f(*first);
+	output_type f_a = f(*first);
 	
 	InputIt second = first;
 	++second;
@@ -53,14 +111,14 @@ typename std::invoke_result::val
 	if (second == last)
 		return f_a;
 	
-	result_type f_b = f(*second);
-	result_type sum = (f_a + f_b) / result_type(2) * (*second - *first);
+	output_type f_b = f(*second);
+	output_type sum = ((f_a + f_b) / input_type(2)) * (*second - *first);
 	
 	for (first = second++; second != last; first = second++)
 	{
 		f_a = f_b;
 		f_b = f(*second);
-		sum += (f_a + f_b) / result_type(2) * (*second - *first);
+		sum += ((f_a + f_b) / input_type(2)) * (*second - *first);
 	}
 	
 	return sum;

src/physics/constants.hpp

@@ -65,7 +65,7 @@ constexpr T stefan_boltzmann = T(5.670374419184429453970996731889230875840122970
 template <class T>
 constexpr T speed_of_light = T(299792458);
 
-} // namespace constant
+} // namespace constants
 
 } // namespace physics