Go back to the "simplified" determinant formula in an attempt to avoid overflows.

diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs
index 2342ba139d2b79cf64c2b5dbac60c3fa0de55b3a..a7b4047398d5a12d5102618c03aad60bddbcd0d6 100644 (file)
--- a/src/Tetrahedron.hs
+++ b/src/Tetrahedron.hs
@@ -293,17 +293,23 @@ c _ _ _ _ _ = error "coefficient index out of bounds"
 --
 --   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 =
-  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 + x3*y1*z4 +
-  x3*y2*z1 - x3*y2*z4 - x3*y4*z1 - x2*y4*z3 - x3*y1*z2 + x3*y4*z2 +
-  x4*y1*z2 - x4*y1*z3 - x4*y2*z1 + x4*y2*z3 + x4*y3*z1 - x4*y3*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,