-
Notifications
You must be signed in to change notification settings - Fork 56
Description
uniformRM (l, h)
for Float
and Double
is defined by drawing x
from [0, 1] using uniformDouble01M
/ uniformFloat01M
and then returning x * l + (1 - x) * h
.
One might expect that uniformRM (0, 1)
may produce the same values as uniformDouble01M
/ uniformFloat01M
. However, this is not the case because of rounding errors when calculating 1 - x
.
For visualization of the rounding problem try:
mapM_ print [ (x == x', x', x)
| let m = fromIntegral (maxBound :: Word64) :: Double
, i <- [0, 2 ^ 8 .. 2 ^ 14]
, let x = fromIntegral i / m -- resembles `uniformDouble01M`
, let x' = 1 - (1 - x)
]
Note that the x'
values have a lot of repetitions while each x
is different.
An easy fix would be changing the calculation to (1 - x) * l + x * h
. However, then the possible results of uniformRM (-1, 0)
are reduced instead.
A possible solution could be to first (mathematically) transform the two calculations as follows:
x * l + (1 - x) * h = h + x * (l - h)
(1 - x) * l + x * h = l + x * (h - l)
Then use calculation 2 when l
is closer to zero than h
, and calculation 2 otherwise. This may be implemented as follows:
instance UniformRange Float where -- analogously for `Double`
uniformRM (l, h)
| l == h = const $ return l
| isInfinite l || isInfinite h = const $ return $! h + l
-- Optimisation exploiting absorption:
-- (-Infinity) + (anything but +Infinity) = -Infinity
-- (anything but -Infinity) + (+Infinity) = +Infinity
-- (-Infinity) + (+Infinity) = NaN
| let f = if abs l < abs h then (let d = h - l in \ x -> l + x * d)
else (let d = l - h in \ x -> h + x * d)
= fmap f . uniformFloat01M
{-# INLINE uniformRM #-}
(This also drops the binding for the second parameter of uniformRM
in order to allow memoization of the bound preprocessing.)
This approach also slightly reduces the documented floating point number caveats, because this implementation guarantees that the bound that is closer to 0 is not exceeded.
Note: Obviously, this would break backwards compatibility as it changes the result of uniformRM
for a given state of the RNG.