RandomLib  1.8
MPFR interface
Back to Other random distributions. Forward to Saving and restoring the state. Up to Contents.

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 header-only 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:

Times (us) for sampling from the exponential distribution
prec(1)(2)(3)(4)
4 0.3010
8 0.307.7
16 0.307.0
24 0.0560.210.306.4
32 0.338.6
48 0.358.5
53 0.0790.290.358.5
64 0.0690.300.358.4
128 0.3713
256 0.4123
1024 0.6482
4096 1.5 520
16384 4.9 5200
65536 19 53000
262144 81 490000
1048576300 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:

Times (us) for sampling from the normal distribution
prec(1)(2)(3)(3')(4)(4')
4 0.911.20.62 2.7
8 0.911.20.61 2.4
16 0.911.20.61 2.6
24 0.0340.540.911.20.65 2.8
32 0.951.20.92 3.0
48 0.961.30.92 3.1
53 0.0470.610.961.30.93 3.2
64 0.0430.640.971.31.0 3.4
128 0.991.31.2 4.8
256 1.0 1.31.5 7.7
1024 1.3 1.53.7 26
4096 2.1 2.418 160
16384 5.4 5.7150 1600
65536 19 19 1300 15000
262144 77 76 11000140000
1048576300 300870001000000

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 timing tests used RandomLib 1.4, MPFR 3.1.0, g++ 4.6.1 on a 64-bit Intel Xeon 2.66GHz running Linux Fedora Release 15.
• Methods (1) use standard algorithms for sampling from these distributions. These are reasonably accurate but do not yield correctly rounded results (in the sense of MPFR).
• Methods (2) use adaptions of von Neumann's method. These give correctly rounded results.
• Methods (3) adapt methods (2) to MPFR. For high precision, the time scales linearly with precision. Method (3'), which uses RandomLib::MPFRExponential internally, is an alternate method of sampling from the normal distribution using Kahn's algorithm. Damien Stehlé suggested the use of Kahn's algorithm to me. For high precision, the time is proportional to the precision for these methods and they all have the same constant of proportionality.
• Methods (4) adapt methods (1) to MPFR. For high precision, the time scales approximately as the 3/2 power of the precision. Method (4') uses the grandom function.
• RandomLib::MPFRExponential and RandomLib::MPFRNormal have meager internal documentation. Refer to RandomLib::ExactExponential and RandomLib::ExactNormal (which implement the same algorithms) for more information.

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:

/**
* \file MPFRExample.cpp
* \brief An example of calling RandomLib::MPFRNormal.
*
* Compile, link, and run with, e.g.,
* g++ -I../include -O2 -o MPFRExample MPFRExample.cpp -lmpfr -lgmp
* ./MPFRExample
*
* http://randomlib.sourceforge.net/
**********************************************************************/
#include <cstdio>
#include <iostream>
#include <ctime> // for time()
int main() {
#if HAVE_MPFR
gmp_randstate_t r;
gmp_randinit_mt(r);
time_t t0 = std::time(0);
gmp_randseed_ui(r, t0);
mpfr_t z;
{
mpfr_prec_t prec = 240; mpfr_init2(z, prec);
std::cout << "Sample from the unit normal distribution at precision "
<< prec << "\n";
RandomLib::MPFRNormal<> norm; // bits = 32, by default
for (int k = 0; k < 10; ++k) {
norm(z, r, MPFR_RNDN); // Obtain a normal deviate
mpfr_out_str(stdout, 10, 0, z, MPFR_RNDN); std::cout << "\n";
}
}
{
mpfr_prec_t prec = 20; mpfr_set_prec(z, prec);
std::cout << "Sample ranges from the normal distribution at precision "
<< prec << "\n";
RandomLib::MPFRNormal<1> norm; // choose bits = 1 so that the ranges
RandomLib::MPFRRandom<1> x; // are not too narrow
for (int k = 0; k < 10; ++k) {
norm(x, r); // Obtain an MPFRRandom range
x(z, MPFR_RNDD); // Lower bound of range
std::cout << "[" << mpfr_get_d(z, MPFR_RNDD) << ",";
x(z, MPFR_RNDU); // Upper bound of range
std::cout << mpfr_get_d(z, MPFR_RNDU) << "] -> ";
x(z, r, MPFR_RNDN); // Realize the normal deviate
mpfr_out_str(stdout, 10, 0, z, MPFR_RNDN); std::cout << "\n";
}
}
// Clean up
mpfr_clear(z); mpfr_free_cache(); gmp_randclear(r);
return 0;
#else
std::cerr << "Need MPFR version 3.0 or later to run MPFRExample\n";
return 1;
#endif // HAVE_MPFR
}

Typical output from this program is:

Sample from the unit normal distribution at precision 240
-5.6692556912807683301494494874784943926662835210924611284269456310097367350e-2
8.5362069032310026750031502283620341952982028910093176612208572465159728631e-1
...
Sample ranges from the normal distribution at precision 20
[0,1] -> 7.2809029e-1
[-2.25,-2.28125] -> -2.2580643
...

Back to Other random distributions. Forward to Saving and restoring the state. Up to Contents.