GeomPredicates-SSE-0.2: Geometric predicates (Intel SSE)



Hardware based, exact computation using the ESSA algorithm in double precision (1)

  • We're using Float inputs on Double precision ESSA at the moment. Hopefully later we can add support for Double inputs on Quadruple ESSA
  • Line intersection is based on the algorithm presented in (4)
  • ccw and incircle based on (5)
  • See (2) and (3) for more information on the splitDouble operation
  • We assume realToFrac is broken and that CFloat == Float and CDouble == Double
  1. Helmut Ratschek, Jon Rokne. "Exact computation of the sign of a finite sum". Applied Mathematics and Computation, Volume 99, Issue 2-3, Pages: 99-127, ISSN:0096-3003, 1999.
  2. Siegfried M. Rump. "High precision evaluation of nonlinear functions"
  3. T.J. Dekker. "A Floating-Point Technique for Extending the Available Precision". Numerische Mathematik, 18:224-242, 1971.
  4. Marina Gavrilova, Jon Rokne. "Reliable line segment intersection testing"
  5. Marina Gavrilova, Helmut Ratschek and Jon Rokne. "Exact computation of Voronoi diagram and Delaunay triangulation"



cinttESSA :: Float -> Float -> Float -> BoolSource

Test if p3 is within the closed interval specified by [p1,p2]

ccwESSA :: Vector2 Float -> Vector2 Float -> Vector2 Float -> OrderingSource

Counter-clockwise orientation test. Classifies p3 in relation to the line formed by p1 and p2.

Result: LT=Right, GT=Left, EQ=Coincident

incircleESSA :: (Vector2 Float, Vector2 Float, Vector2 Float) -> Vector2 Float -> OrderingSource

Test the relation of a point to the circle formed by (p1..p3). (p1..p3) must be in counterclockwise order.

Result: GT=inside, EQ=border, LT=outside.

Note: this is the sum of 192 multiplications.

essa :: (Functor t, Foldable t) => t Double -> OrderingSource

Compute the exact sign of the sum of the input sequence. It is the caller's responsibility to ensure that the inputs have not suffered a loss of precision.

splitDouble :: Double -> (Double, Double)Source

Split a 53 bit double into two 26 bit halves so that: let (lo,hi) = splitDouble x in x == lo + hi

The trick is that the sign is used as the additional bit.

Note that the multiplication of two 26-bit floating point numbers is exact in double precision.

If you're new to this function you may want to read paper (5), as using this function properly may be trickier than it seems.