diff --git a/fggl/math/shapes.cpp b/fggl/math/shapes.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b67f6f52318a66fe680e7a5faf6c35f1418be909
--- /dev/null
+++ b/fggl/math/shapes.cpp
@@ -0,0 +1,178 @@
+/*
+ * ${license.title}
+ * Copyright (C) 2022 ${license.owner}
+ * ${license.mailto}
+ *
+ * This program 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.
+ *
+ * This program 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 this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ */
+
+#include "fggl/math/shapes.hpp"
+#include <vector>
+
+namespace fggl::math::phs3d {
+
+	void AABB::emtpy() {
+		min = {FLT_MAX, FLT_MAX, FLT_MAX};
+		max = { - FLT_MAX, - FLT_MAX, -FLT_MAX};
+	}
+
+	void AABB::add(const math::vec3 &p) {
+		if ( p.x < min.x ) min.x = p.x;
+		if ( p.x > max.x ) max.x = p.x;
+		if ( p.y < min.y ) min.y = p.y;
+		if ( p.y > min.y ) max.y = p.y;
+		if ( p.z < min.z ) min.z = p.z;
+		if ( p.z > max.z ) max.z = p.z;
+	}
+
+	AABB AABB::fromPoints(const std::vector<math::vec3> &points) {
+		AABB box;
+		for (const auto& point : points) {
+			box.add(point);
+		}
+		return box;
+	}
+
+	void AABB::set(const AABB &other, const math::mat4 &m) {
+		// TODO this needs testing, I'm not sure if our vectors are oriented the same as 9.4.4
+		min = max = m[3]; // should be the translation component of the matrix
+
+		// this feels like something that should be vectorizable...
+		for ( int i = 0; i < 3; i++) {
+			for (int j = 0; j < 3; j++) {
+				if ( m[i][j] > 0.0f ) {
+					min[j] += m[i][j] * other.min[j];
+					max[j] += m[i][j] * other.max[j];
+				} else {
+					min[j] += m[i][j] * other.max[j];
+					max[j] += m[i][j] * other.min[j];
+				}
+			}
+		}
+	}
+
+	Plane Plane::fromPoints(const math::vec3 p1, const math::vec3 p2, const math::vec3 p3) {
+		const auto e3 = p2 - p1;
+		const auto e1 = p3 - p2;
+		auto normal = glm::normalize(glm::cross(e3, e1));
+		auto d = glm::dot(p2, normal);
+		return {normal, d};
+	}
+
+	static math::vec3 bestFitNormal(const std::vector<math::vec3>& points) {
+		assert(!points.empty());
+
+		math::vec3 result;
+
+		math::vec3 p = points.back();
+		for ( std::size_t i = 0; i < points.size(); ++i ){
+			math::vec3 c = points[i];
+
+			result.x += (p.z + c.z) * (p.y - c.y);
+			result.y += (p.x + c.x) * (p.z - c.z);
+			result.z += (p.y + c.y) - (p.x - c.x);
+
+			p = c;
+		}
+
+		return glm::normalize(result);
+	};
+
+	static float bestFitD(const std::vector<math::vec3>& points, glm::vec3 normal) {
+		math::vec3 sum;
+		for (auto& point : points) {
+			sum += point;
+		}
+		sum *= 1.0F/points.size();
+		return glm::dot( sum, normal );
+	}
+
+	const char X = 0;
+	const char Y = 1;
+	const char Z = 2;
+
+	struct bary_axis {
+		const char a;
+		const char b;
+	};
+
+	static bary_axis baryCalcAxis(const math::vec3& normal) {
+		if ( (fabs(normal.x) >= fabs(normal.y)) && (fabs(normal.x) >= fabs(normal.z))) {
+			// discard x
+			return {Y, Z};
+		} else if ( fabs(normal.y) >= fabs(normal.z) ) {
+			// discard y
+			return {Z, X};
+		} else {
+			// discard z
+			return {X, Y};
+		}
+	}
+
+	bool Triangle::CartToBarycentric(const math::vec3& cart, Barycentric& outVal) {
+		// everything is const because I'm paying the compiler is smarter than me...
+
+		const auto d1 = v[1] - v[0];
+		const auto d2 = v[2] - v[1];
+		const auto n = glm::cross(d1, d2);
+
+		// 0 = x, 1 = y, 2 = z
+		const auto ax = baryCalcAxis(n);
+
+		// first axis
+		const float u1 = v[0][ax.a] - v[2][ax.a];
+		const float u2 = v[1][ax.a] - v[2][ax.a];
+		const float u3 = cart[ax.a] - v[0][ax.a];
+		const float u4 = cart[ax.a] - v[2][ax.a];
+
+		// second axis
+		const float v1 = v[0][ax.b] - v[2][ax.b];
+		const float v2 = v[1][ax.b] - v[2][ax.b];
+		const float v3 = cart[ax.b] - v[0][ax.b];
+		const float v4 = cart[ax.b] - v[2][ax.b];
+
+		const float denom = v1*u2 - v2*u1;
+		if ( denom == 0.0f) {
+			return false;
+		}
+
+		// finally, we can work it out
+		const float oneOverDenom = 1.0f / denom;
+		outVal.b[0] = (v4*u2 - v2*u4) * oneOverDenom;
+		outVal.b[1] = (v1*u3 - v3*u1) * oneOverDenom;
+		outVal.b[2] = 1.0f - outVal.b[0] - outVal.b[1];
+		return true;
+	}
+
+	bool Triangle::CartToBarycentric2(const math::vec3& cart, Barycentric& outVal) {
+		const auto e1 = v[2] - v[1];
+		const auto e2 = v[0] - v[2];
+		const auto e3 = v[1] - v[0];
+
+		const auto d1 = cart - v[0];
+		const auto d2 = cart - v[1];
+		const auto d3 = cart - v[2];
+
+		const auto normal = glm::normalize(glm::cross(e1, e2));
+		const auto denom = glm::dot(glm::cross(e1, e2), normal);
+		assert( denom != 0.0f);
+
+		outVal.b[0] = glm::dot(glm::cross(e1, d3), normal) / denom;
+		outVal.b[1] = glm::dot(glm::cross(e2, d1), normal) / denom;
+		outVal.b[2] = glm::dot(glm::cross(e3, d2), normal) / denom;
+		return true;
+	}
+
+} // namespace fggl::math::phys3d
\ No newline at end of file
diff --git a/include/fggl/math/fmath.hpp b/include/fggl/math/fmath.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..265e58b17b4dee8444345e054333def41f434113
--- /dev/null
+++ b/include/fggl/math/fmath.hpp
@@ -0,0 +1,49 @@
+/*
+ * ${license.title}
+ * Copyright (C) 2022 ${license.owner}
+ * ${license.mailto}
+ *
+ * This program 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.
+ *
+ * This program 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 this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ */
+
+#ifndef FGGL_MATH_FMATH_H
+#define FGGL_MATH_FMATH_H
+
+#include <cmath>
+
+namespace fggl::math {
+
+	// wrap value in range [0, max)
+	inline float wrap(float value, float max) {
+		return fmod(max + fmod(value, max), max);
+	}
+
+	// wrap value in range [min, max)
+	inline float wrap(float value, float min, float max) {
+		if ( min > max ){ std::swap(min, max); };
+		value -= min;
+		float rangeSize = max - min;
+		return value - (rangeSize* std::floor(value/rangeSize)) + min;
+	}
+
+	// if the value is out of bounds, return that bound
+	inline float clamp(float value, float min, float max) {
+		const float t = value < min ? min : value;
+		return t > max ? max : t;
+	}
+
+} // namespace fggl::math
+
+#endif //FGGL_MATH_FMATH_H
diff --git a/include/fggl/math/imath.hpp b/include/fggl/math/imath.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ae0fe4f5be1ab5bf89dd3d3b51bc958ce0e7ef46
--- /dev/null
+++ b/include/fggl/math/imath.hpp
@@ -0,0 +1,61 @@
+/*
+ * ${license.title}
+ * Copyright (C) 2022 ${license.owner}
+ * ${license.mailto}
+ *
+ * This program 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.
+ *
+ * This program 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 this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ */
+
+//
+// Created by webpigeon on 20/04/22.
+//
+
+#ifndef FGGL_MATH_IMATH_HPP
+#define FGGL_MATH_IMATH_HPP
+
+namespace fggl::math {
+
+	// wrap value in range [0, max)
+	inline int wrap(int value, int max) {
+		if ( value < 0 ) return (n-1)-(-1-value) % n;
+		if ( value >= max ) return value % max;
+		return value;
+	}
+
+	// wrap value in range [min, max)
+	inline int wrap(int value, int const min, int const max) {
+		int range = max - min + 1;
+		if ( value < min) {
+			value += range * ( (min - value) / range + 1 );
+		}
+		return min + (value - min) % range;
+	}
+
+	inline
+	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) {
+		const T i = value > max ? max : value;
+		return i < min ? min : value;
+	}
+
+} // namespace fggl::math
+
+#endif //FGGL_MATH_IMATH_HPP
diff --git a/include/fggl/math/shapes.hpp b/include/fggl/math/shapes.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1a5c7124be8d33eee89c880f233edca38cc6c085
--- /dev/null
+++ b/include/fggl/math/shapes.hpp
@@ -0,0 +1,210 @@
+/*
+ * ${license.title}
+ * Copyright (C) 2022 ${license.owner}
+ * ${license.mailto}
+ *
+ * This program 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.
+ *
+ * This program 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 this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ */
+
+#ifndef FGGL_MATH_PHYS3D_SHAPES_HPP
+#define FGGL_MATH_PHYS3D_SHAPES_HPP
+
+#include "fggl/math/types.hpp"
+#include <vector>
+#include <array>
+
+// based on https://gamemath.com/book/geomprims.html
+namespace fggl::math::phs3d {
+
+	struct Line {
+		math::vec2 normal;
+		float d;
+	};
+
+	struct LineSegment {
+		math::vec3 start;
+		math::vec3 end;
+	};
+
+	struct Ray {
+		math::vec3 start;
+		math::vec3 length;
+
+		inline math::vec3 end() {
+			return start + length;
+		}
+
+		static Ray fromPos(math::vec3 start, math::vec3 end);
+	};
+
+	// nice bounding shapes
+
+	struct Sphere {
+		math::vec3 center;
+		float radius;
+
+		float volume() const;
+		float area() const;
+	};
+
+	struct AABB {
+		math::vec3 min;
+		math::vec3 max;
+
+		void add(const math::vec3& p);
+		void emtpy();
+		void set(const AABB& other, const math::mat4& m);
+
+		inline math::vec3 center() const {
+			return (min + max) / 2.0F;
+		}
+
+		inline math::vec3 size() const {
+			return max - min;
+		}
+
+		inline math::vec3 radius() const {
+			return max - center();
+		}
+
+		static AABB fromPoints(const std::vector<math::vec3>& points);
+	};
+
+	struct Plane {
+		glm::vec3 normal;
+		float d; // distance to origin
+
+		inline bool contains(const math::vec3& point) {
+			return glm::dot(point, normal) == d;
+		}
+
+		inline float distance(const math::vec3& q) {
+			return glm::dot(q, normal) - d;
+		}
+
+		static Plane fromPoints(const math::vec3, const math::vec3, const math::vec3);
+		static Plane bestFit(const std::vector<math::vec3>& points);
+	};
+
+	struct Barycentric {
+
+		std::array<float, 3> b;
+
+		bool inTriangle() {
+			return 0 <= b[0] && b[0] <= 1 &&
+				   0 <= b[1] && b[1] <= 1 &&
+				   0 <= b[2] && b[2] <= 1;
+		}
+
+		inline bool isLegal() const {
+			return (b[0] + b[1] + b[2]) == 1.0f;
+		}
+	};
+
+	constexpr float THIRD = 1.0F/3.0F;
+	const Barycentric cGrav = { {THIRD, THIRD, THIRD}};
+
+	struct Triangle {
+		std::array<math::vec3, 3> v;
+
+		inline float l1() const {
+			return length(2, 1);
+		}
+
+		inline float l2() const {
+			return length(0, 2);
+		}
+
+		inline float l3() const {
+			return length(1, 0);
+		}
+
+		inline math::vec3 edge(int v1, int v2) const {
+			return v[v1] - v[v2];
+		}
+
+		inline float length(int a, int b) const {
+			return glm::length( edge(a, b));
+		}
+
+		inline float perimeter() const {
+			return length(3, 2) + length(1,3) + length(2, 1);
+		}
+
+		inline float area() const {
+			return glm::length(glm::cross(edge(3, 2), edge(1, 3))) / 2;
+		}
+
+		inline math::vec3 cGrav() const {
+			return (v[0] + v[1] + v[2]) / 3.0f;
+		}
+
+		inline math::vec3 inCenter() const {
+			auto p = perimeter();
+			Barycentric bs{
+				l1() / p,
+				l2() / p,
+				l3() / p
+			};
+			return BarycentricToCart(bs);
+		}
+
+		inline float centerRadius() const {
+			return area() / perimeter();
+		}
+
+		// circle that passes through all three points = cirumcenter + circumraidus
+		inline math::vec3 circumcenter() {
+			// d1 = - e2 DOT e3
+			// d2 = - e3 DOT e1
+			// d3 = - e1 DOT e2
+			// C1 = d2 * d3
+			// C2 = d3 * d1
+			// C3 = d1 * d2
+			// c = c1 + c2 + c3
+			// c_2 = 2 * c
+			// bary{
+			//   (c2 + c3 ) / c_2
+			//   (c3 + c1 ) / c_2
+			//   (c1 + c2 ) / c_2
+			// return BarycentricToCart(bary)
+		}
+
+		inline float circumradius() {
+			// d1 = - e2 DOT e3
+			// d2 = - e3 DOT e1
+			// d3 = - e1 DOT e2
+			// C1 = d2 * d3
+			// C2 = d3 * d1
+			// C3 = d1 * d2
+			// c = c1 + c2 + c3
+			// c_2 = 2 * c
+			// return sqrt( (d1+d2) * (d2+d3) * ( d3 + d1 ) / c ) / 2
+		}
+
+		bool CartToBarycentric(const math::vec3& cart, Barycentric& outVal);
+		bool CartToBarycentric2(const math::vec3& cart, Barycentric& outVal);
+
+		math::vec3 BarycentricToCart(const Barycentric& p) const {
+			return v[0] * p.b[0] +
+				   v[1] * p.b[1] +
+				   v[2] * p.b[2];
+		}
+	};
+
+} // namespace fggl::math::phys3d
+
+
+#endif //FGGL_MATH_PHYS3D_SHAPES_HPP