diff --git a/include/fggl/math/imath.hpp b/include/fggl/math/imath.hpp
index 4c3996fd5b21111cdaf09840010a7f6b11cf46f3..337352b99af9b6d030c041d8cb9cab19a0a4b3e1 100644
--- a/include/fggl/math/imath.hpp
+++ b/include/fggl/math/imath.hpp
@@ -19,10 +19,12 @@
 #ifndef FGGL_MATH_IMATH_HPP
 #define FGGL_MATH_IMATH_HPP
 
+#include <cstdint>
+
 namespace fggl::math {
 
 	// wrap value in range [0, max)
-	inline int wrap(int value, int max) {
+	constexpr int wrap(int value, int max) {
 		if (value < 0)
 			return (value - 1) - (-1 - value) % value;
 		if (value >= max)
@@ -31,7 +33,7 @@ namespace fggl::math {
 	}
 
 	// wrap value in range [min, max)
-	inline int wrap(int value, int const min, int const max) {
+	constexpr int wrap(int value, int const min, int const max) {
 		int range = max - min + 1;
 		if (value < min) {
 			value += range * ((min - value) / range + 1);
@@ -39,19 +41,47 @@ namespace fggl::math {
 		return min + (value - min) % range;
 	}
 
-	inline
-	int clamp(int value, int min, int max) {
+	constexpr int clamp(int value, int min, int max) {
 		const int i = value > max ? max : value;
 		return i < min ? min : value;
 	}
 
 	template<typename T>
-	inline
-	T clampT(T value, T min, T max) {
+	constexpr T clampT(T value, T min, T max) {
 		const T i = value > max ? max : value;
 		return i < min ? min : value;
 	}
 
+	/**
+	 * Calculate the sum of the first N positive integers.
+	 *
+	 * @param n the number to stop at
+	 * @return sum(0 ... N)
+	 */
+	constexpr uint64_t calcSum(uint64_t n) {
+		return (n * (n + 1)) / 2;
+	}
+
+	/**
+	 * Calculate the squared sum of the first N positive integers.
+	 *
+	 * @param n the number to stop at
+	 * @return sum(0 ... N) * sum(0 ... N)
+	 */
+	constexpr uint64_t calcSumSquared(uint64_t n) {
+		return (n * (n + 1) * (2*n + 1)) / 6;
+	}
+
+	/**
+	 * Calculate the cubed sum of the first N positive integers.
+	 *
+	 * @param n the number to stop at
+	 * @return sum(0 ... N) * sum(0 ... N) * sum(0 ... N)
+	 */
+	constexpr uint64_t calcSumCubed(uint64_t n) {
+		return ((n * n) * ((n + 1) * (n + 1))) / 4;
+	}
+
 } // namespace fggl::math
 
 #endif //FGGL_MATH_IMATH_HPP
diff --git a/include/fggl/math/statistics.hpp b/include/fggl/math/statistics.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..fadec87a254017055e5f4e2dce3fb59283a14fe1
--- /dev/null
+++ b/include/fggl/math/statistics.hpp
@@ -0,0 +1,148 @@
+/*
+ * This file is part of FGGL.
+ *
+ * FGGL is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any
+ * later version.
+ *
+ * FGGL is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
+ * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License along with FGGL.
+ * If not, see <https://www.gnu.org/licenses/>.
+ */
+
+//
+// Created by webpigeon on 27/11/22.
+//
+
+#ifndef FGGL_MATH_STATS_HPP
+#define FGGL_MATH_STATS_HPP
+
+#include <cstdint>
+#include <cmath>
+
+namespace fggl::math {
+
+	class CumulativeAverage {
+		public:
+			inline void add(float item) {
+				m_current = ( item + (m_count * m_current) ) / (m_count + 1);
+				m_count++;
+			}
+
+			inline float average() {
+				return m_current;
+			}
+
+		private:
+			float m_current;
+			std::size_t m_count;
+	};
+
+	/**
+	 * Useful Statistics class.
+	 * ported to C++ from Piers' Java implementation.
+	 *
+	 * @tparam T
+	 */
+	template<typename T>
+	class Statistics {
+		public:
+			void add(T observation) {
+				m_count++;
+				m_sum += observation;
+				m_sumSquared += (observation * observation);
+				m_min = std::min(m_min, observation);
+				m_max = std::max(m_max, observation);
+				m_dirty = true;
+			}
+
+			inline T average() {
+				compute();
+				return m_average;
+			}
+
+			inline T average() const {
+				if ( m_dirty ) {
+					return m_sum / m_count;
+				}
+				return m_average;
+			}
+
+			inline T standardDeviation() {
+				compute();
+				return m_standardDiv;
+			}
+
+			inline T standardDeviation() const {
+				if ( !m_dirty ) {
+					return m_standardDiv;
+				}
+				T avg = average();
+				T num = m_sumSquared - (m_count * avg * avg);
+				return num < 0 ? 0 : sqrt( num / (m_count-1));
+			}
+
+			inline T standardError() {
+				compute();
+				return m_standardDiv / sqrt(m_count);
+			}
+
+			inline T standardError() const {
+				if ( !m_dirty ) {
+					return m_standardDiv / sqrt(m_count);
+				}
+				return standardDeviation() / sqrt(m_count);
+			}
+
+			inline T squareSum () const {
+				return m_sumSquared;
+			}
+
+			inline std::size_t count() const {
+				return m_count;
+			}
+
+			inline T sum() const {
+				return m_sum;
+			}
+
+			inline T min() const {
+				return m_min;
+			}
+
+			inline T max() const {
+				return m_max;
+			}
+
+		private:
+			T m_average{0};
+			T m_sum{0};
+			T m_sumSquared{0};
+			T m_standardDiv{0};
+			T m_min{0};
+			T m_max{0};
+			std::size_t m_count{0};
+			bool m_dirty{false};
+
+			void compute() {
+				if ( !m_dirty ) {
+					return;
+				}
+				m_average = m_sum / m_count;
+				T num = m_sumSquared - (m_count * m_average * m_average);
+				if ( num < 0 )  {
+					num = 0;
+				}
+				m_standardDiv = sqrt( num / (m_count-1) );
+				m_dirty = false;
+			}
+	};
+
+	using StatisticsF = Statistics<float>;
+	using StatisticsD = Statistics<double>;
+
+} // namespace fggl::math
+
+#endif //FGGL_MATH_STATS_HPP