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

Note: Adaptions of MPFRExponential and MPFRNormal will be included in the next (post 3.1.2) release of MPFR as functions mpfr_erandom and mpfr_nrandom.

The classes MPFRUniform, MPFRExponential, and MPFRNormal, provide a efficient methods of computing random deviates and returning the results in the form of MPFR numbers. These classes (together with 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.,

#include <RandomLib/MPFRNormal.hpp>

(There are 3 additional classes included in RandomLib: MPFRExponentialL, MPFRNormalK, and 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 MPFRUniform, MPFRExponential, and 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.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) ExponentialDistribution;
(2) ExactExponential, with bits = 32;
(3) MPFRExponential, with bits = 32, which uses the same algorithm as ExactExponential;
(4) MPFRExponentialL, which uses the same algorithm as 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.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) NormalDistribution;
(2) ExactNormal, with bits = 32;
(3) MPFRNormal, with bits = 32, which uses the same algorithm as ExactNormal;
(3') MPFRNormalK, with bits = 32, which uses the Kahn algorithm;
(4) MPFRNormalR, which uses the same algorithm as 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 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.
MPFRExponential and MPFRNormal have meager internal documentation. Refer to ExactExponential and 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 MPFRRandom class holds this intermediate object. The utility of MPFRRandom is seen by considering Kahn's algorithm for normal sampling: pick y and z from the exponential distribution until (y − 1)² < 2z and then return y. If y and z are returned as 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 MPFRNormalK for details. MPFRUniform is a thin wrapper for MPFRRandom which returns samples in [0,1).

The classes MPFRRandom, MPFRUniform, MPFRExponential, MPFRNormal, and MPFRNormalK all take a template parameter bits (default value 32) which gives the number of bits in each "digit" of the 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 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
 *
 * Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed under
 * the MIT/X11 License.  For more information, see
 * http://randomlib.sourceforge.net/
 **********************************************************************/
#include <cstdio>
#include <iostream>
#include <ctime>                // for time()
#include <RandomLib/MPFRRandom.hpp>
#include <RandomLib/MPFRNormal.hpp>
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.