Skip to content

make uniformRM (0, 1) for floats equivalent to uniformDouble01M / uniformFloat01M #166

@Flupp

Description

@Flupp

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:

  1. x * l + (1 - x) * h = h + x * (l - h)
  2. (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.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions