+-- | Compute the determinant of the 4x4 matrix,
+--
+-- [1]
+-- [x]
+-- [y]
+-- [z]
+--
+-- where [1] = [1, 1, 1, 1],
+-- [x] = [x1,x2,x3,x4],
+--
+-- et cetera.
+--
+-- The termX nonsense is an attempt to prevent Double overflow.
+-- which has been observed to happen with large coordinates.
+--
det :: Point -> Point -> Point -> Point -> Double
det p0 p1 p2 p3 =
--- Both of these results are just copy/pasted from Sage. One of them
--- might be more numerically stable, faster, or both.
---
--- x1*y2*z4 - x1*y2*z3 + x1*y3*z2 - x1*y3*z4 - x1*y4*z2 + x1*y4*z3 +
--- x2*y1*z3 - x2*y1*z4 - x2*y3*z1 + x2*y3*z4 +
--- x2*y4*z1 - x2*y4*z3 - x3*y1*z2 + x3*y1*z4 + x3*y2*z1 - x3*y2*z4 - x3*y4*z1 +
--- x3*y4*z2 + x4*y1*z2 - x4*y1*z3 - x4*y2*z1 + x4*y2*z3 + x4*y3*z1 - x4*y3*z2
- -((x2 - x3)*y1 - (x1 - x3)*y2 + (x1 - x2)*y3)*z4 + ((x2 - x4)*y1 - (x1 - x4)*y2 + (x1 - x2)*y4)*z3 + ((x3 - x4)*y2 - (x2 - x4)*y3 + (x2 - x3)*y4)*z1 - ((x3 - x4)*y1 - (x1 - x4)*y3 + (x1 - x3)*y4)*z2
+ term5 + term6
where
(x1, y1, z1) = p0
(x2, y2, z2) = p1
(x3, y3, z3) = p2
(x4, y4, z4) = p3
+ term1 = ((x2 - x4)*y1 - (x1 - x4)*y2 + (x1 - x2)*y4)*z3
+ term2 = ((x2 - x3)*y1 - (x1 - x3)*y2 + (x1 - x2)*y3)*z4
+ term3 = ((x3 - x4)*y2 - (x2 - x4)*y3 + (x2 - x3)*y4)*z1
+ term4 = ((x3 - x4)*y1 - (x1 - x4)*y3 + (x1 - x3)*y4)*z2
+ term5 = term1 - term2
+ term6 = term3 - term4
-- | Computed using the formula from Lai & Schumaker, Definition 15.4,