Implement my own 4x4 determinant.

diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs
index 613b863a49953b57120884b15c83d7ebd733a187..e84b09103d7ad1bb850bed4b9380aeb47d6bb3b3 100644 (file)
--- a/src/Tetrahedron.hs
+++ b/src/Tetrahedron.hs
@@ -17,7 +17,7 @@ import qualified Data.Vector as V (
   snoc,
   sum
   )
-import Numeric.LinearAlgebra hiding (i, scale)
+
 import Prelude hiding (LT)
 import Test.Framework (Test, testGroup)
 import Test.Framework.Providers.HUnit (testCase)
@@ -281,31 +281,39 @@ c _ _ _ _ _ = error "coefficient index out of bounds"
 
 
 
--- | The matrix used in the tetrahedron volume calculation as given in
---   Lai & Schumaker, Definition 15.4, page 436.
-vol_matrix :: Tetrahedron -> Matrix Double
-vol_matrix t = (4><4)
-               [1,  1,  1,  1,
-                x1, x2, x3, x4,
-                y1, y2, y3, y4,
-                z1, z2, z3, z4 ]
-    where
-      (x1, y1, z1) = v0 t
-      (x2, y2, z2) = v1 t
-      (x3, y3, z3) = v2 t
-      (x4, y4, z4) = v3 t
+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
+  where
+    (x1, y1, z1) = p0
+    (x2, y2, z2) = p1
+    (x3, y3, z3) = p2
+    (x4, y4, z4) = p3
+
 
 -- | Computed using the formula from Lai & Schumaker, Definition 15.4,
 --   page 436.
 volume :: Tetrahedron -> Double
 volume t
-       | (v0 t) == (v1 t) = 0
-       | (v0 t) == (v2 t) = 0
-       | (v0 t) == (v3 t) = 0
-       | (v1 t) == (v2 t) = 0
-       | (v1 t) == (v3 t) = 0
-       | (v2 t) == (v3 t) = 0
-       | otherwise = (1/6)*(det (vol_matrix t))
+       | v0' == v1' = 0
+       | v0' == v2' = 0
+       | v0' == v3' = 0
+       | v1' == v2' = 0
+       | v1' == v3' = 0
+       | v2' == v3' = 0
+       | otherwise = (1/6)*(det v0' v1' v2' v3')
+  where
+    v0' = v0 t
+    v1' = v1 t
+    v2' = v2 t
+    v3' = v3 t
 
 
 -- | The barycentric coordinates of a point with respect to v0.