From: Michael Orlitzky Date: Mon, 3 Oct 2011 18:50:38 +0000 (-0400) Subject: Go back to the "simplified" determinant formula in an attempt to avoid overflows. X-Git-Tag: 0.0.1~113 X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=f66108655b9d4c3ba5ab8493c38882d00098509f;p=spline3.git Go back to the "simplified" determinant formula in an attempt to avoid overflows. --- diff --git a/src/Tetrahedron.hs b/src/Tetrahedron.hs index 2342ba1..a7b4047 100644 --- 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,