Add quadratic solver

libs/geom/include/psemek/geom/math.hpp

@@ -2,6 +2,9 @@
 
 #include <psemek/geom/constants.hpp>
 
+#include <cmath>
+#include <optional>
+
 namespace psemek::geom
 {
 
@@ -29,4 +32,20 @@ namespace psemek::geom
 		return x0 * (T(1) - t) + x1 * t;
 	}
 
+	// roots of ax^2 + bx + c = 0
+	// in increasing order
+	template <typename T>
+	std::optional<std::pair<T, T>> solve_quadratic(T const & a, T const & b, T const & c)
+	{
+		auto const D = b * b - 4 * a * c;
+		if (D < 0) return std::nullopt;
+
+		auto t = - (b + (b > 0 ? 1 : -1) * std::sqrt(D));
+		auto x1 = t / a;
+		auto x2 = c / t;
+
+		if (x1 > x2) std::swap(x1, x2);
+		return std::pair{x1, x2};
+	}
+
 }