From ef90a0fc3362cfe092642f6a7770fd2f6d0b0f4e Mon Sep 17 00:00:00 2001
From: Leonhard Markert <curiousleo@users.noreply.github.com>
Date: Wed, 2 Dec 2020 19:20:03 +0100
Subject: [PATCH] Optimise handling of infinities

---
 src/System/Random/Internal.hs | 26 ++++++++++++--------------
 1 file changed, 12 insertions(+), 14 deletions(-)

diff --git a/src/System/Random/Internal.hs b/src/System/Random/Internal.hs
index d1916d3cd..be879949d 100644
--- a/src/System/Random/Internal.hs
+++ b/src/System/Random/Internal.hs
@@ -855,16 +855,15 @@ instance UniformRange Bool where
 instance UniformRange Double where
   uniformRM (l, h) g
     | l == h = return l
-    | isInfiniteL && isInfiniteH = return (0/0) -- NaN
-    | isInfiniteL = return l
-    | isInfiniteH = return h
+    | isInfinite l || isInfinite h =
+      -- Optimisation exploiting absorption:
+      --   (-Infinity) + (anything but +Infinity) = -Infinity
+      --   (anything but -Infinity) + (+Infinity) = +Infinity
+      --                (-Infinity) + (+Infinity) = NaN
+      return $! h + l
     | otherwise = do
       x <- uniformDouble01M g
       return $ x * l + (1 -x) * h
-    where
-      -- Optimisations
-      isInfiniteL = isInfinite l
-      isInfiniteH = isInfinite h
 
 -- | Generates uniformly distributed 'Double' in the range \([0, 1]\).
 --   Numbers are generated by generating uniform 'Word64' and dividing
@@ -896,16 +895,15 @@ uniformDoublePositive01M g = (+ d) <$> uniformDouble01M g
 instance UniformRange Float where
   uniformRM (l, h) g
     | l == h = return l
-    | isInfiniteL && isInfiniteH = return (0/0) -- NaN
-    | isInfiniteL = return l
-    | isInfiniteH = return h
+    | isInfinite l || isInfinite h =
+      -- Optimisation exploiting absorption:
+      --   (-Infinity) + (anything but +Infinity) = -Infinity
+      --   (anything but -Infinity) + (+Infinity) = +Infinity
+      --                (-Infinity) + (+Infinity) = NaN
+      return $! h + l
     | otherwise = do
       x <- uniformFloat01M g
       return $ x * l + (1 - x) * h
-    where
-      -- Optimisations
-      isInfiniteL = isInfinite l
-      isInfiniteH = isInfinite h
 
 -- | Generates uniformly distributed 'Float' in the range \([0, 1]\).
 --   Numbers are generated by generating uniform 'Word32' and dividing