diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3f3126f..b9647dd 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,9 @@
 This changelog is generated by `changelog.py` using manually added semantic version tags to classify commits as breaking changes, additions or fixes.
 
 
+#### [v15.3.1](https://github.com/nophead/NopSCADlib/releases/tag/v15.3.1 "show release") Fixes [...](https://github.com/nophead/NopSCADlib/compare/v15.3.0...v15.3.1 "diff with v15.3.0")
+*   2021-04-02 [`f3376ed`](https://github.com/nophead/NopSCADlib/commit/f3376edaf186b32f442b94d6d0b42f1ba0c7612c "show commit") [C.P.](# "Chris Palmer") Documented `xor()` function.
+
 ### [v15.3.0](https://github.com/nophead/NopSCADlib/releases/tag/v15.3.0 "show release") Additions [...](https://github.com/nophead/NopSCADlib/compare/v15.2.0...v15.3.0 "diff with v15.2.0")
 *   2021-04-02 [`c073419`](https://github.com/nophead/NopSCADlib/commit/c073419c0b4eddcda4cda5bd0f8d48268b6e58ec "show commit") [C.P.](# "Chris Palmer") Added `opengrab_screw_depth()` function.
 
diff --git a/readme.md b/readme.md
index 33acd90..96dbd5d 100644
--- a/readme.md
+++ b/readme.md
@@ -5834,6 +5834,7 @@ Maths utilities for manipulating vectors and matrices.
 | `circle_intersect(c1, r1, c2, r2)` | Calculate one point where two circles in the X-Z plane intersect, clockwise around c1 |
 | `cosh(x)` | hyperbolic cosine |
 | `coth(x)` | hyperbolic cotangent |
+| `cubic_real_roots(a, b, c, d)` | Returns real roots of cubic equation |
 | `degrees(radians)` | Convert degrees to radians |
 | `euler(R)` | Convert a rotation matrix to a Euler rotation vector. |
 | `identity(n, x = 1)` | Construct an arbitrary size identity matrix |
@@ -5841,6 +5842,7 @@ Maths utilities for manipulating vectors and matrices.
 | `map(v, func)` | make a new vector where the func function argument is applied to each element of the vector v |
 | `mapi(v, func)` | make a new vector where the func function argument is applied to each element of the vector v. The func will get the index number as first argument, and the element as second argument. |
 | `nearly_zero(x)` | True if x is close to zero |
+| `quadratic_real_roots(a, b, c)` | Returns real roots of a quadratic equation, biggest first. Returns empty list if no real roots |
 | `radians(degrees)` | Convert radians to degrees |
 | `reduce(v, func, unity)` | reduce a vector v to a single entity by applying the func function recursively to the reduced value so far and the next element, starting with unity as the initial reduced value |
 | `reverse(v)` | Reverse a vector |
diff --git a/utils/maths.scad b/utils/maths.scad
index 6d878ab..7ae3208 100644
--- a/utils/maths.scad
+++ b/utils/maths.scad
@@ -162,3 +162,28 @@ function sumv(v) = reduce(v, function(a, b) a + b, 0); //! sum a vector of value
 
 function xor(a,b) = (a && !b) || (!a && b);     //! Logical exclusive OR
 function cuberoot(x)= sign(x)*abs(x)^(1/3);
+
+function quadratic_real_roots(a, b, c) =        //! Returns real roots of a quadratic equation, biggest first. Returns empty list if no real roots
+    let(2a = 2 * a,
+        2c = 2 * c,
+        det = b^2 - 2a * 2c
+    ) det < 0 ? [] :
+        let(r = sqrt(det),
+            x1 = b < 0 ? 2c / (-b + r) : (-b - r) / 2a,
+            x2 = b < 0 ? (-b + r) / 2a : 2c / (-b - r)
+        ) [x2, x1];
+
+function cubic_real_roots(a, b, c, d) = //! Returns real roots of cubic equation
+    let(b = b / a,
+        c = c / a,
+        d = d / a,
+        inflection = -b / 3,
+        p = c - b^2 / 3,
+        q = 2 * b^3 / 27 - b * c / 3 + d,
+        det = q^2 / 4 + p^3 / 27,
+        roots = !p && !q ? 1 : nearly_zero(det) ? 2 :  det < 0 ? 3 : 1,
+        r = sqrt(det),
+        x = cuberoot(-q / 2 - r) + cuberoot(-q / 2 + r)
+    ) roots == 1 ? [x] :
+      roots == 2 ? [3 * q /p + inflection, -3 * q / p / 2 + inflection] :
+      [for(i = [0 : roots - 1]) 2 * sqrt(-p / 3) * cos(acos(3 * q * sqrt(-3 / p) / p / 2) - i * 120) + inflection];