RandomLib
1.10

The Random class generates two classes of uniform real results: fixed point numbers (where the spacing between possible results is a constant) and floating point numbers (where the spacing between possible results varies).
The results returned by Fixed(), FixedU(), FixedN(), FixedW(), FixedO(), and FixedC() are "fixedpoint reals" with precision p. These are of the form i / 2^{p} where i is an integer. If, for real data type RealType, we restrict p > 0 and p ≤ std::numeric_limits<RealType>::digits, then all such numbers in [1,1] are representable. For p = 3, the set of numbers in [0,1] is {0, 1/8, 1/4, 3/8, 1/2, 5/8, 3/4, 7/8, 1}.
The results returned FixedS() are of "offset fixedpoint reals" with precision p. These are of the form (i + 1/2) / 2^{p} where i is an integer. If, for real data type RealType, we restrict p > 0 and p ≤ std::numeric_limits<RealType>::digits, then all such numbers in (−1/2,1/2) are representable. (This only "works" for radix 2 systems; with larger bases, the results can't be represented exactly.) Note that possibly "exceptional" numbers, −1/2, 0, and 1/2, are not included in this set. For p = 3, the set of numbers in (−1/2,1/2) is {−7/16, −5/16, −3/16, −1/16, 1/16, 3/16, 5/16, 7/16}.
The results returned by Float(), FloatU(), FloatN(), FloatW() are "floatingpoint reals" with precision p and exponent range e. The possible results for such floating numbers in [−1,1] consist of
For real data type RealType, such numbers are representable if
For e = 0, the numbers become the fixed point numbers with precision p.
Here is an example of the floating point points with p = 3 and e = 2 together with the probabilities yielded by Float(), FloatU(), FloatN().
X  prob(Float() = X)  prob(FloatU() = X)  prob(FloatN() = X) 

0  1/32  0  0.5/32 
1/32  1/32  1/32  1/32 
2/32  1/32  1/32  1/32 
3/32  1/32  1/32  1/32 
4/32  1/32  1/32  1/32 
5/32  1/32  1/32  1/32 
6/32  1/32  1/32  1/32 
7/32  1/32  1/32  1/32 
4/16  1/16  1/32  1.5/32 
5/16  1/16  1/16  1/16 
6/16  1/16  1/16  1/16 
7/16  1/16  1/16  1/16 
4/8  1/8  1/16  1.5/16 
5/8  1/8  1/8  1/8 
6/8  1/8  1/8  1/8 
7/8  1/8  1/8  1/8 
1  0  1/8  0.5/8 
The description of floatingpoint numbers assumes that the underlying hardware supports denormalized numbers. This is the case with most modern computers. The code attempts to deal also with older hardware where there's a gap between 0 and 1/2^{e + 1}, but this has not been tested.
The following table provides a succinct definition of each of the member functions of Random routines that return a real result. Here u is a uniformly distributed random number in (0,1). This is drawn from a continuous distribution; i.e., it may be thought of as consisting of a binary point followed by an infinite sequence of random binary bits. (This is just a useful conceptual framework. None of the implementations of these functions require an explicit realization of u.)
The term "fixed" means a fixedpoint real with precision p, and we have h = 1/2^{p} (the smallest positive fixedpoint number). The term "float" means a floatingpoint real with precision p and exponent range e.
routine  mnemonic  definition 

Fixed()  default (down)  round u down to previous fixed 
FixedU()  upper  round u up to next fixed 
FixedN()  nearest  round u to nearest fixed 
FixedW()  wide  round 2u − 1 to nearest fixed 
FixedS()  symmetric  round u − 1/2 to nearest offset fixed 
FixedO()  open  round (1 − h)u up to nearest fixed 
FixedC()  closed  round (1 + h)u down to nearest fixed 
Float()  default (down)  round u down to previous float 
FloatU()  upper  round u up to next float 
FloatN()  nearest  round u to nearest float 
FloatW()  wide  round 2u − 1 to nearest float 
The precision and exponent range are determined as follows. Each of the fixedpoint routines comes in 3 variants, for example,
Similarly each of the floatingpoint routines comes in 3 variants, for example,
Typical values of digits and min_exponent are given by
type  digits  min_exponent 

float  24  125 
double  53  1021 
long double  64  16381 
long double (Power PC)  106  968 
long double (Sun)  113  16494 
In the following tables, the columns have the following meanings
next(X) is the next representable float following X. prev(X) is the previous representable float preceding X.
routine  min  max  num  prob 

Fixed()  0  1−h  2^{p}  h 
FixedU()  h  1  2^{p}  h 
FixedN()  0  1  2^{p} + 1  h (h/2 at endpoints) 
FixedW()  −1  1  2^{p+1} + 1  h/2 (h/4 at endpoints) 
FixedS()  −(1−h)/2  (1−h)/2  2^{p}  h 
FixedO()  h  1−h  2^{p}−1  h/(1−h) 
FixedC()  0  1  2^{p}+1  h/(1+h) 
Float()  0  1−h  2^{p}(1+e/2)  min(1,next(X))−X 
FloatU()  1/2^{p+e}  1  2^{p}(1+e/2)  X−max(0,prev(X)) 
FloatN()  0  1  2^{p}(1+e/2)+1  (min(1,next(X))−max(0,prev(X)))/2 
FloatW()  −1  1  2^{p+1}(1+e/2)+1  (min(1,next(X))−max(−1,prev(X)))/4 
From these definitions, we can show that:
We can easily show that certain distributions are equivalent:
Function  Equivalent 

FixedU<RealType,p>()  Fixed<RealType,p>() + h 
FixedU<RealType,p>()  1 − Fixed<RealType,p>() 
FixedN<RealType,p>()  Boolean() ? Fixed<RealType,p>() : FixedU<RealType,p>() 
FixedN<RealType,p>()  1 − FixedN<RealType,p>() 
FixedO<RealType,p>()  1 − FixedO<RealType,p>() 
FixedC<RealType,p>()  1 − FixedC<RealType,p>() 
FixedS<RealType,p>()  Fixed<RealType,p>() − (1−h)/2 
FixedS<RealType,p>()  − FixedS<RealType,p>() 
FixedW<RealType,p>()  (Boolean() ? 1 : −1) FixedN<RealType,p>() 
FixedW<RealType,p>()  − FixedW<RealType,p>() 
FixedW<RealType,p−1>()  2 FixedN<RealType,p>() − 1 
FloatN<RealType,p,e>()  Boolean() ? Float<RealType,p,e>() : FloatU<RealType,p,e>() 
FloatW<RealType,p,e>()  (Boolean() ? 1 : −1) FloatN<RealType,p,e>() 
FloatW<RealType,p,e>()  − FloatW<RealType,p,e>() 
Float<RealType,p,0>()  Fixed<RealType,p>() 
FloatU<RealType,p,0>()  FixedU<RealType,p>() 
FloatN<RealType,p,0>()  FixedN<RealType,p>() 
FloatW<RealType,p,0>()  FixedW<RealType,p>() 
A caution about FixedO() and FixedC(). All four of FixedS(), FixedN() − 0.5, FixedO() − 0.5, and FixedC() − 0.5 produce results which are strictly symmetric about 0 and would thus be suitable for an unbiased random walk. The variances of the first two distributions are 1/12 + O(h^{2}), close to the ideal value. On the other hand the variances of the distributions of FixedO() and FixedC() are 1/12 + O(h), significantly further from the ideal value. (In practice, using a strictly symmetric normal distribution is preferable for simulating a random walk.)
Similarly, consider estimating the value of π by randomly selecting points in a unit square and determining what fraction lie in a circle of diameter 1. Sampling in the square using FixedO() or FixedC() gives poorer results than FixedN(), and FixedS() gives slightly better results.
Because of their definitions in terms of u, any of Fixed(), FixedU(), FixedN(), or FixedS() can be used to cover the periodic unit interval in an unbiased way. Thus to sample an angle uniformly, use 2 * π * Fixed() or 2 * π * FixedS(). The latter has the advantage that it is strictly symmetric about zero. In addition, angles which are multiples of π/2 are avoided (which may obviate the need to check for special cases).
Usually, these real routines would be invoked by specifying the type and allowing the precision to be determined from the type, e.g., Real<double>(). However, in some cases it might be useful to specify a lower precision: