[Haskell-cafe] Working around floating-point inaccuracies

[Michail Pevnev]

Hello everyone!

For a library I write at the moment, I made a class for approximate
geometric equality comparisons (and dubbed it `ApproxEq`). It's not a
very lawful class, but there is one invariant that should hold
absolutely always - if two entities coincide but have different
representations, they should compare "equal". 

All instances for my geometric entities seem to uphold this invariant
(btw, QuickCheck is awesome). Except for the `Plane` type, which gives
me a false negative once in about 50000 runs. I've chosen to represent
planes by two vectors: an arbitrary point on a plane, and a normal. To
test this invariant I translate the plane by a vector parallel to the
plane, and then compare the translated plane with the original. 

The general pattern to these failures is that one of the coordinates of
the origin point of the untranslated plane is huge compared to the other
two - by 4-6 orders of magnitude, while the translation vector is not
particularly big and has about the same magnitude as the two smaller
coordinates.

I suspect that this difference throws off the function that tests if a
point belongs to a plane (which is in turn used in the equality
comparison) by filling lower digits of participating numbers with
floating-point garbage.

I can live with this false negative, but I'd prefer not to. There
probably is a clever way to work around it, I just can't find it. Any
ideas? Thanks in advance. The relevant code is below.

```
-- | A plane is represented by an arbitrary point on it and a normal vector.
data Plane a = Plane (V3 a) (V3 a)
           deriving (Show)

instance (Epsilon a, Floating a) => ApproxEq (Plane a) where
    p1@(Plane o1 n1) ~== p2@(Plane o2 n2) = res
        where n1' = normalize n1
              n2' = normalize n2
              parallel = n1' ~== n2'
              antiparallel = (-n1') ~== n2'
              originsAreTheSame = o1 ~== o2
              originsBelong = pointOnPlane p1 o2 || pointOnPlane p2 o1
              res = (parallel || antiparallel) && (originsAreTheSame || originsBelong)

-- | Return True if a point lies on a plane.
pointOnPlane :: (Epsilon a, Floating a) => Plane a -> V3 a -> Bool
pointOnPlane (Plane o n) v = nearZero $ delta `dot` n'
    where delta = normalize $ o - v
          n' = normalize n
```

-- 
Michail.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 488 bytes
Desc: not available
URL: <http://mail.haskell.org/pipermail/haskell-cafe/attachments/20180805/45080659/attachment.sig>