RandomLib
1.9

Note: Adaptions of RandomLib::MPFRExponential and RandomLib::MPFRNormal will be included in the next (post 3.1.2) release of MPFR as functions mpfr_erandom
and mpfr_nrandom
.
The classes RandomLib::MPFRUniform, RandomLib::MPFRExponential, and RandomLib::MPFRNormal, provide a efficient methods of computing random deviates and returning the results in the form of MPFR numbers. These classes (together with RandomLib::MPFRRandom) are headeronly implementations which do not depend on the rest of RandomLib. In order to use these headers, you can copy them into the directory tree where you are developing an MPFR application (make sure they stay within a directory called RandomLib) and include them with e.g.,
(There are 3 additional classes included in RandomLib: RandomLib::MPFRExponentialL, RandomLib::MPFRNormalK, and RandomLib::MPFRNormalR. But these are included for illustrative purposes only and therefore they are marked as deprecated in the documentation.)
For high precision, the time per sample for RandomLib::MPFRUniform, RandomLib::MPFRExponential, and RandomLib::MPFRNormal is proportional to the precision where the constant of proportionality is governed by generating the required number of random bits and copying these into the result. Thus, using these methods, the generation of random numbers in MPFR is amongst the cheapest of operations (cheaper, for example, than multiplication). The timing data is giving in the tables below.
This table gives the time required to generate a random sample from an exponential distribution:
prec  (1)  (2)  (3)  (4) 

4  0.30  10  
8  0.30  7.7  
16  0.30  7.0  
24  0.056  0.21  0.30  6.4 
32  0.33  8.6  
48  0.35  8.5  
53  0.079  0.29  0.35  8.5 
64  0.069  0.30  0.35  8.4 
128  0.37  13  
256  0.41  23  
1024  0.64  82  
4096  1.5  520  
16384  4.9  5200  
65536  19  53000  
262144  81  490000  
1048576  300  4100000 
Key to methods for sampling from the exponential distribution:
(1) RandomLib::ExponentialDistribution;
(2) RandomLib::ExactExponential, with bits = 32;
(3) RandomLib::MPFRExponential, with bits = 32, which uses the same algorithm as RandomLib::ExactExponential;
(4) RandomLib::MPFRExponentialL, which uses the same algorithm as RandomLib::ExponentialDistribution (taking the log of a uniform deviate).
This table gives the time required to generate a random sample from an normal distribution:
prec  (1)  (2)  (3)  (3')  (4)  (4') 

4  0.91  1.2  0.62  2.7  
8  0.91  1.2  0.61  2.4  
16  0.91  1.2  0.61  2.6  
24  0.034  0.54  0.91  1.2  0.65  2.8 
32  0.95  1.2  0.92  3.0  
48  0.96  1.3  0.92  3.1  
53  0.047  0.61  0.96  1.3  0.93  3.2 
64  0.043  0.64  0.97  1.3  1.0  3.4 
128  0.99  1.3  1.2  4.8  
256  1.0  1.3  1.5  7.7  
1024  1.3  1.5  3.7  26  
4096  2.1  2.4  18  160  
16384  5.4  5.7  150  1600  
65536  19  19  1300  15000  
262144  77  76  11000  140000  
1048576  300  300  87000  1000000 
Key to methods for sampling from the normal distribution:
(1) RandomLib::NormalDistribution;
(2) RandomLib::ExactNormal, with bits = 32;
(3) RandomLib::MPFRNormal, with bits = 32, which uses the same algorithm as RandomLib::ExactNormal;
(3') RandomLib::MPFRNormalK, with bits = 32, which uses the Kahn algorithm;
(4) RandomLib::MPFRNormalR, which uses the same algorithm as RandomLib::NormalDistribution (the ratio method);
(4') MPFR's grandom, which uses the polar method.
Notes:
The von Neumann methods work in two phases. (1) Some number of uniform random numbers are consumed to form the initial digits of the random sample. (2) Additional digits are copied directly from the random number generator to the random sample. It is convenient to be able to interrupt the process after phase (1) and the RandomLib::MPFRRandom class holds this intermediate object. The utility of RandomLib::MPFRRandom is seen by considering Kahn's algorithm for normal sampling: pick y and z from the exponential distribution until (y − 1)^{2} < 2z and then return y. If y and z are returned as RandomLib::MPFRRandom objects then only sufficient digits need to be generated to determine the acceptance test. If the test passes, then digits only need to be added to y (and not to z). Since the implementation of RandomLib::MPFRNormalK for details. RandomLib::MPFRUniform is a thin wrapper for RandomLib::MPFRRandom which returns samples in [0,1).
The classes RandomLib::MPFRRandom, RandomLib::MPFRUniform, RandomLib::MPFRExponential, RandomLib::MPFRNormal, and RandomLib::MPFRNormalK all take a template parameter bits (default value 32) which gives the number of bits in each "digit" of the RandomLib::MPFRRandom object. This must evenly divide GMP_LIMB_BITS. For portability you should not set this to 64. Typically 32 will provide the best efficiency. Smaller values are useful for debugging, for studying algorithmic complexity (see Knuth and Yao, 1976), or if random bits are expensive (e.g., they are being generated by a hardware generator).
Here is an example of using RandomLib::MPFRNormal:
Typical output from this program is:
Sample from the unit normal distribution at precision 240 5.6692556912807683301494494874784943926662835210924611284269456310097367350e2 8.5362069032310026750031502283620341952982028910093176612208572465159728631e1 ... Sample ranges from the normal distribution at precision 20 [0,1] > 7.2809029e1 [2.25,2.28125] > 2.2580643 ...