Apply different scaling computation

bench/Main.hs

@@ -4,6 +4,7 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 module Main (main) where
 
+import Data.Bits
 import Control.Monad
 import Control.Monad.State.Strict
 import Data.Int

src/System/Random.hs

@@ -357,7 +357,7 @@ class Random a where
   -- independently:
   --
   -- >>> fst $ randomR (('a', 5.0), ('z', 10.0)) $ mkStdGen 26
-  -- ('z',7.27305019146949)
+  -- ('z',7.72694980853051)
   --
   -- In case when a lawful range is desired `uniformR` should be used
   -- instead.

src/System/Random/Internal.hs

@@ -68,6 +68,7 @@ module System.Random.Internal
   , shuffleListM
   , isInRangeOrd
   , isInRangeEnum
+  , scaleFloating
 
   -- * Generators for sequences of pseudo-random bytes
   , uniformByteStringM
@@ -1394,13 +1395,16 @@ instance UniformRange Double where
     | l == h = return l
     | isInfinite l || isInfinite h =
       -- Optimisation exploiting absorption:
-      --   (-Infinity) + (anything but +Infinity) = -Infinity
-      --   (anything but -Infinity) + (+Infinity) = +Infinity
-      --                (-Infinity) + (+Infinity) = NaN
+      --    (+Infinity) + (-Infinity) = NaN
+      --    (-Infinity) + (+Infinity) = NaN
+      --    (+Infinity) + _           = +Infinity
+      --    (-Infinity) + _           = -Infinity
+      --              _ + (+Infinity) = +Infinity
+      --              _ + (-Infinity) = -Infinity
       return $! h + l
     | otherwise = do
       x <- uniformDouble01M g
-      return $! x * l + (1 - x) * h
+      pure $! scaleFloating l h x
   {-# INLINE uniformRM #-}
   isInRange = isInRangeOrd
 
@@ -1437,13 +1441,48 @@ instance UniformRange Float where
   uniformRM (l, h) g
     | l == h = return l
     | isInfinite l || isInfinite h =
+      -- Optimisation exploiting absorption:
+      --    (+Infinity) + (-Infinity) = NaN
+      --    (-Infinity) + (+Infinity) = NaN
+      --    (+Infinity) + _           = +Infinity
+      --    (-Infinity) + _           = -Infinity
+      --              _ + (+Infinity) = +Infinity
+      --              _ + (-Infinity) = -Infinity
       return $! h + l
     | otherwise = do
       x <- uniformFloat01M g
-      return $! x * l + (1 - x) * h
+      pure $! scaleFloating l h x
   {-# INLINE uniformRM #-}
   isInRange = isInRangeOrd
 
+-- | This is the function that is used to scale a floating point values from @[0, 1]@ range
+-- to the custom @[low, high]@ range.
+--
+-- @since 1.3.0
+scaleFloating ::
+     (Ord a, Num a)
+  => a
+  -- ^ Low
+  -> a
+  -- ^ High
+  -> a
+  -- ^ Value in [0, 1] range that needs to be scaled.
+  -> a
+scaleFloating l h x =
+  if abs l < abs h
+    -- Without clamping value furthest from zero we can get into a situation where scaled value can
+    -- get outside of it. This is not a problem for the value closest to zero, due to the
+    -- conditional above.
+    then let !y = l + x * (h - l)
+          in if l < h -- `l` is closer to zero
+               then min y h -- `l` is the true low, ensure `y` is not higher than `h`
+               else max y h -- `l` is the high, ensure `y` is not lower than `h`
+    else let !y = h + x * (l - h)
+          in if l < h -- `h` is closer to zero
+               then max y l -- `h` is the true high, ensure `y` is not lower than `l`
+               else min y l -- `h` is the low, ensure `y` is not higher than `l`
+{-# INLINE scaleFloating #-}
+
 -- | Generates uniformly distributed 'Float' in the range \([0, 1]\).
 --   Numbers are generated by generating uniform 'Word32' and dividing
 --   it by \(2^{32}\). It's used to implement 'UniformRange' instance for 'Float'.

src/System/Random/Stateful.hs

@@ -127,6 +127,7 @@ module System.Random.Stateful
   -- $implemenstatefulegen
 
   -- ** Floating point number caveats #fpcaveats#
+  , scaleFloating
   -- $floating
 
   -- * References
@@ -771,79 +772,66 @@ applyTGen f (TGenM tvar) = do
 -- $floating
 --
 -- The 'UniformRange' instances for 'Float' and 'Double' use the following
--- procedure to generate a random value in a range for @uniformRM (a, b) g@:
+-- procedure to generate a random value in a range for @uniformRM (l, h) g@:
 --
--- If \(a = b\), return \(a\). Otherwise:
+-- * If @l == h@, return: @l@.
+-- * If @`isInfinite` l == True@ or @`isInfinite` h == True@, return: @l + h@
+-- * Otherwise:
 --
--- 1.  Generate \(x\) uniformly such that \(0 \leq x \leq 1\).
+--     1.  Generate @x@ uniformly such that \(0 \leq x \leq 1\).
 --
---     The method by which \(x\) is sampled does not cover all representable
---     floating point numbers in the unit interval. The method never generates
---     denormal floating point numbers, for example.
+--         The method by which @x@ is sampled does not cover all representable
+--         floating point numbers in the unit interval. The method never generates
+--         denormal floating point numbers, for example.
 --
--- 2.  Return \(x \cdot a + (1 - x) \cdot b\).
+--     2.  Return: @if `abs` l < `abs` h then l + x * (h - l) else h + x * (l - h)@
 --
---     Due to rounding errors, floating point operations are neither
---     associative nor distributive the way the corresponding operations on
---     real numbers are. Additionally, floating point numbers admit special
---     values @NaN@ as well as negative and positive infinity.
+--         It is defined in the library as `scaleFloating` function.
 --
--- For pathological values, step 2 can yield surprising results.
---
--- *   The result may be greater than @max a b@.
---
---     >>> :{
---     let (a, b, x) = (-2.13238e-29, -2.1323799e-29, 0.27736077)
---         result = x * a + (1 - x) * b :: Float
---     in (result, result > max a b)
---     :}
---     (-2.1323797e-29,True)
---
--- *   The result may be smaller than @min a b@.
---
---     >>> :{
---     let (a, b, x) = (-1.9087862, -1.908786, 0.4228573)
---         result = x * a + (1 - x) * b :: Float
---     in (result, result < min a b)
---     :}
---     (-1.9087863,True)
+--         Due to rounding errors, floating point operations are neither
+--         associative nor distributive the way the corresponding operations on
+--         real numbers are. Additionally, floating point numbers admit special
+--         values @NaN@ as well as negative and positive infinity.
 --
 -- What happens when @NaN@ or @Infinity@ are given to 'uniformRM'? We first
 -- define them as constants:
 --
 -- >>> nan = read "NaN" :: Float
 -- >>> inf = read "Infinity" :: Float
+-- >>> g <- newIOGenM (mkStdGen 2024)
 --
--- *   If at least one of \(a\) or \(b\) is @NaN@, the result is @NaN@.
+-- *   If at least one of \(l\) or \(h\) is @NaN@, the result is @NaN@.
 --
---     >>> let (a, b, x) = (nan, 1, 0.5) in x * a + (1 - x) * b
+--     >>> scaleFloating nan 1 <$> uniformFloat01M g
 --     NaN
---     >>> let (a, b, x) = (-1, nan, 0.5) in x * a + (1 - x) * b
+--     >>> scaleFloating (-1) nan <$> uniformFloat01M g
 --     NaN
 --
--- *   If \(a\) is @-Infinity@ and \(b\) is @Infinity@, the result is @NaN@.
+-- *   If \(l\) and \(h\) are both @Infinity@ with opposing signes, then the result is @NaN@.
 --
---     >>> let (a, b, x) = (-inf, inf, 0.5) in x * a + (1 - x) * b
+--     >>> scaleFloating (-inf) inf <$> uniformFloat01M g
+--     NaN
+--     >>> scaleFloating inf (-inf) <$> uniformFloat01M g
 --     NaN
 --
--- *   Otherwise, if \(a\) is @Infinity@ or @-Infinity@, the result is \(a\).
+-- *   Otherwise, if \(l\) is @Infinity@ or @-Infinity@, the result is \(l\).
 --
---     >>> let (a, b, x) = (inf, 1, 0.5) in x * a + (1 - x) * b
+--     >>> scaleFloating inf 1 <$> uniformFloat01M g
 --     Infinity
---     >>> let (a, b, x) = (-inf, 1, 0.5) in x * a + (1 - x) * b
+--     >>> scaleFloating (-inf) 1 <$> uniformFloat01M g
 --     -Infinity
 --
--- *   Otherwise, if \(b\) is @Infinity@ or @-Infinity@, the result is \(b\).
+-- *   Otherwise, if \(h\) is @Infinity@ or @-Infinity@, the result is \(h\).
 --
---     >>> let (a, b, x) = (1, inf, 0.5) in x * a + (1 - x) * b
+--     >>> scaleFloating 1 inf <$> uniformFloat01M g
 --     Infinity
---     >>> let (a, b, x) = (1, -inf, 0.5) in x * a + (1 - x) * b
+--     >>> scaleFloating 1 (-inf) <$> uniformFloat01M g
 --     -Infinity
 --
 -- Note that the [GCC 10.1.0 C++ standard library](https://gcc.gnu.org/git/?p=gcc.git;a=blob;f=libstdc%2B%2B-v3/include/bits/random.h;h=19307fbc3ca401976ef6823e8fda893e4a263751;hb=63fa67847628e5f358e7e2e7edb8314f0ee31f30#l1859),
 -- the [Java 10 standard library](https://docs.oracle.com/javase/10/docs/api/java/util/Random.html#doubles%28double,double%29)
 -- and [CPython 3.8](https://github.com/python/cpython/blob/3.8/Lib/random.py#L417)
--- use the same procedure to generate floating point values in a range.
+-- use a similar procedure to generate floating point values in a range.
 --
 -- $implemenstatefulegen
 --