RandomLib  1.10
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
MPFRExponential.hpp
Go to the documentation of this file.
1 /**
2  * \file MPFRExponential.hpp
3  * \brief Header for MPFRExponential
4  *
5  * Sampling exactly from the normal distribution for MPFR.
6  *
7  * Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed under
8  * the MIT/X11 License. For more information, see
9  * http://randomlib.sourceforge.net/
10  **********************************************************************/
11 
12 #if !defined(RANDOMLIB_MPFREXPONENTIAL_HPP)
13 #define RANDOMLIB_MPFREXPONENTIAL_HPP 1
14 
15 #include <RandomLib/MPFRRandom.hpp>
16 
17 #if HAVE_MPFR || defined(DOXYGEN)
18 
19 namespace RandomLib {
20 
21  /**
22  * \brief The exponential distribution for MPFR.
23  *
24  * This is a transcription of ExactExponential (version 1.4) for use with
25  * MPFR.
26  *
27  * This class uses mutable private objects. So a single MPFRExponential
28  * object cannot safely be used by multiple threads. In a multi-processing
29  * environment, each thread should use a thread-specific MPFRExponential
30  * object.
31  *
32  * @tparam bits the number of bits in each digit.
33  **********************************************************************/
34  template<int bits = 32> class MPFRExponential {
35  public:
36 
37  /**
38  * Initialize the MPFRExponential object.
39  **********************************************************************/
41  /**
42  * Sample from the exponential distribution with mean 1 returning a
43  * MPFRRandom.
44  *
45  * @param[out] t the MPFRRandom result.
46  * @param[in,out] r a GMP random generator.
47  **********************************************************************/
48  void operator()(MPFRRandom<bits>& t, gmp_randstate_t r) const
49  { Compute(r); _x.swap(t); }
50  /**
51  * Sample from the exponential distribution with mean 1.
52  *
53  * @param[out] val the sample from the exponential distribution
54  * @param[in,out] r a GMP random generator.
55  * @param[in] round the rounding direction.
56  * @return the MPFR ternary result (&plusmn; if val is larger/smaller than
57  * the exact sample).
58  **********************************************************************/
59  int operator()(mpfr_t val, gmp_randstate_t r, mpfr_rnd_t round) const
60  { Compute(r); return _x(val, r, round); }
61  private:
62  // Disable copy constructor and assignment operator
64  MPFRExponential& operator=(const MPFRExponential&);
65  int ExpFraction(gmp_randstate_t r, MPFRRandom<bits>& p) const {
66  // The early bale out
67  if (p.TestHighBit(r)) return 0;
68  // Implement the von Neumann algorithm
69  _w.Init();
70  if (!_w.LessThan(r, p)) return 1;
71  while (true) {
72  _v.Init(); if (!_v.LessThan(r, _w)) return 0;
73  _w.Init(); if (!_w.LessThan(r, _v)) return 1;
74  }
75  }
76  void Compute(gmp_randstate_t r) const {
77  _x.Init();
78  unsigned k = 0;
79  while (!ExpFraction(r, _x)) { ++k; _x.Init(); }
80  if (k & 1) _x.SetHighBit(r);
81  _x.AddInteger(k >> 1);
82  return;
83  }
84  mutable MPFRRandom<bits> _x;
85  mutable MPFRRandom<bits> _v;
86  mutable MPFRRandom<bits> _w;
87  };
88 
89 } // namespace RandomLib
90 
91 #endif // HAVE_MPFR
92 #endif // RANDOMLIB_MPFREXPONENTIAL_HPP
void operator()(MPFRRandom< bits > &t, gmp_randstate_t r) const
int TestHighBit(gmp_randstate_t r)
Definition: MPFRRandom.hpp:227
Handling random numbers in MPFR.
Definition: MPFRRandom.hpp:52
The exponential distribution for MPFR.
Header for MPFRRandom.
int operator()(mpfr_t val, gmp_randstate_t r, mpfr_rnd_t round) const