diff --git a/src/tg/detail/functions.hh b/src/tg/detail/functions.hh
index 48f006fd3389876f718ebd3401fc4f2f882888f4..d48c1ebd1f57b61f00c4442ce4bdc7239b875990 100644
--- a/src/tg/detail/functions.hh
+++ b/src/tg/detail/functions.hh
@@ -15,6 +15,7 @@
 #include "functions/distance.hh"
 #include "functions/dot.hh"
 #include "functions/interpolate.hh"
+#include "functions/inverse.hh"
 #include "functions/length.hh"
 #include "functions/minmax.hh"
 #include "functions/mix.hh"
diff --git a/src/tg/detail/functions/inverse.hh b/src/tg/detail/functions/inverse.hh
new file mode 100644
index 0000000000000000000000000000000000000000..632c84616f8c7f94bdd5ad0261e1828ac3ade308
--- /dev/null
+++ b/src/tg/detail/functions/inverse.hh
@@ -0,0 +1,53 @@
+#pragma once
+
+#include "../../types/mat.hh"
+#include "../../types/vec.hh"
+#include "../scalar_traits.hh"
+#include "determinant.hh"
+
+namespace tg
+{
+template <class ScalarT>
+mat<1, 1, fractional_result<ScalarT>> inverse(mat<1, 1, ScalarT> const& m)
+{
+    return {ScalarT(1.0) / m[0].x};
+}
+template <class ScalarT>
+mat<2, 2, fractional_result<ScalarT>> inverse(mat<2, 2, ScalarT> const& m)
+{
+    auto invdet = ScalarT(1.0) / determinant(m);
+
+    mat<2, 2, fractional_result<ScalarT>> res;
+    res[0][0] = +m[1][1] * invdet;
+    res[0][1] = -m[0][1] * invdet;
+    res[1][0] = -m[1][0] * invdet;
+    res[1][1] = +m[0][0] * invdet;
+    return res;
+}
+template <class ScalarT>
+mat<3, 3, fractional_result<ScalarT>> inverse(mat<3, 3, ScalarT> const& m)
+{
+    auto invdet = ScalarT(1.0) / determinant(m);
+
+    mat<3, 3, fractional_result<ScalarT>> res;
+    res[0][0] = +(m[1][1] * m[2][2] - m[2][1] * m[1][2]) * invdet;
+    res[0][1] = -(m[0][1] * m[2][2] - m[2][1] * m[0][2]) * invdet;
+    res[0][2] = +(m[0][1] * m[1][2] - m[1][1] * m[0][2]) * invdet;
+    res[1][0] = -(m[1][0] * m[2][2] - m[2][0] * m[1][2]) * invdet;
+    res[1][1] = +(m[0][0] * m[2][2] - m[2][0] * m[0][2]) * invdet;
+    res[1][2] = -(m[0][0] * m[1][2] - m[1][0] * m[0][2]) * invdet;
+    res[2][0] = +(m[1][0] * m[2][1] - m[2][0] * m[1][1]) * invdet;
+    res[2][1] = -(m[0][0] * m[2][1] - m[2][0] * m[0][1]) * invdet;
+    res[2][2] = +(m[0][0] * m[1][1] - m[1][0] * m[0][1]) * invdet;
+    return res;
+}
+// template <class ScalarT>
+// mat<4, 4, fractional_result<ScalarT>> inverse(mat<4, 4, ScalarT> const& m)
+// {
+//     auto invdet = ScalarT(1.0) / determinant(m);
+//
+//     mat<4, 4, fractional_result<ScalarT>> res;
+//     ... TODO ...
+//     return res;
+// }
+} // namespace tg