RandomLib  1.10
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
MPFRNormalR.hpp
Go to the documentation of this file.
1 /**
2  * \file MPFRNormalR.hpp
3  * \brief Header for MPFRNormalR
4  *
5  * Sampling exactly from the normal distribution for MPFR using the ratio
6  * method.
7  *
8  * Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed under
9  * the MIT/X11 License. For more information, see
10  * http://randomlib.sourceforge.net/
11  **********************************************************************/
12 
13 #if !defined(RANDOMLIB_MPFRNORMALR_HPP)
14 #define RANDOMLIB_MPFRNORMALR_HPP 1
15 
16 #include <algorithm> // for max/min
17 #include <cmath> // for pow
18 #include <mpfr.h>
19 
20 #define HAVE_MPFR (MPFR_VERSION_MAJOR >= 3)
21 
22 #if HAVE_MPFR || defined(DOXYGEN)
23 
24 namespace RandomLib {
25 
26  /**
27  * \brief The normal distribution for MPFR (ratio method).
28  *
29  * This class is <b>DEPRECATED</b>. It is included for illustrative purposes
30  * only. The MPFRNormal class provides a much more efficient method for
31  * sampling from the normal distribution.
32  *
33  * This is an adaption of NormalDistribution to MPFR. The changes are
34  * - Use MPFR's random number generator
35  * - Use sufficient precision internally to ensure that a correctly rounded
36  * result is returned.
37  *
38  * This class uses a mutable private object. So a single MPFRNormalR
39  * object cannot safely be used by multiple threads. In a multi-processing
40  * environment, each thread should use a thread-specific MPFRNormalR
41  * object.
42  **********************************************************************/
43  class MPFRNormalR {
44  private:
45  // The number of bits of randomness to add at a time. Require that Leva's
46  // bounds "work" at a precision of 2^-chunk and that an unsigned long can
47  // hold this many bits.
48  static const long chunk_ = 32;
49  static const unsigned long m = 3684067834; // ceil(2^chunk*sqrt(2/e))
50 
51  public:
52  /**
53  * Initialize the MPFRNormalR object.
54  **********************************************************************/
56  mpz_init(_ui);
57  mpz_init(_vi);
58  mpfr_init2(_eps, chunk_);
59  mpfr_init2(_u, chunk_);
60  mpfr_init2(_v, chunk_);
61  mpfr_init2(_up, chunk_);
62  mpfr_init2(_vp, chunk_);
63  mpfr_init2(_vx, chunk_);
64  mpfr_init2(_x1, chunk_);
65  mpfr_init2(_x2, chunk_);
66  }
67  /**
68  * Destroy the MPFRNormalR object.
69  **********************************************************************/
71  mpfr_clear(_x2);
72  mpfr_clear(_x1);
73  mpfr_clear(_vx);
74  mpfr_clear(_vp);
75  mpfr_clear(_up);
76  mpfr_clear(_v);
77  mpfr_clear(_u);
78  mpfr_clear(_eps);
79  mpz_clear(_vi);
80  mpz_clear(_ui);
81  }
82  /**
83  * Sample from the normal distribution with mean 0 and variance 1.
84  *
85  * @param[out] val the sample from the normal distribution
86  * @param[in,out] r a GMP random generator.
87  * @param[in] round the rounding direction.
88  * @return the MPFR ternary result (&plusmn;1 if val is larger/smaller than
89  * the exact sample).
90  **********************************************************************/
91  int operator()(mpfr_t val, gmp_randstate_t r, mpfr_rnd_t round) const {
92  const double
93  s = 0.449871, // Constants from Leva
94  t = -0.386595,
95  a = 0.19600 ,
96  b = 0.25472 ,
97  r1 = 0.27597 ,
98  r2 = 0.27846 ,
99  u1 = 0.606530, // sqrt(1/e) rounded down and up
100  u2 = 0.606531,
101  scale = std::pow(2.0, -chunk_); // for turning randoms into doubles
102 
103  while (true) {
104  mpz_urandomb(_vi, r, chunk_);
105  if (mpz_cmp_ui(_vi, m) >= 0) continue; // Very early reject
106  double vf = (mpz_get_ui(_vi) + 0.5) * scale;
107  mpz_urandomb(_ui, r, chunk_);
108  double uf = (mpz_get_ui(_ui) + 0.5) * scale;
109  double
110  x = uf - s,
111  y = vf - t,
112  Q = x*x + y * (a*y - b*x);
113  if (Q >= r2) continue; // Early reject
114  mpfr_set_ui_2exp(_eps, 1u, -chunk_, MPFR_RNDN);
115  mpfr_prec_t prec = chunk_;
116  mpfr_set_prec(_u, prec);
117  mpfr_set_prec(_v, prec);
118  // (u,v) = sw corner of range
119  mpfr_set_z_2exp(_u, _ui, -prec, MPFR_RNDN);
120  mpfr_set_z_2exp(_v, _vi, -prec, MPFR_RNDN);
121  mpfr_set_prec(_up, prec);
122  mpfr_set_prec(_vp, prec);
123  // (up,vp) = ne corner of range
124  mpfr_add(_up, _u, _eps, MPFR_RNDN);
125  mpfr_add(_vp, _v, _eps, MPFR_RNDN);
126  // Estimate how many extra bits will be needed to achieve the desired
127  // precision.
128  mpfr_prec_t prec_guard = 3 + chunk_ -
129  (std::max)(mpz_sizeinbase(_ui, 2), mpz_sizeinbase(_vi, 2));
130  if (Q > r1) {
131  int reject;
132  while (true) {
133  // Rejection curve v^2 + 4 * u^2 * log(u) < 0 has a peak at u =
134  // exp(-1/2) = 0.60653066. So treat uf in (0.606530, 0.606531) =
135  // (u1, u2) specially
136 
137  // Try for rejection first
138  if (uf <= u1)
139  reject = Reject(_u, _vp, prec, MPFR_RNDU);
140  else if (uf >= u2)
141  reject = Reject(_up, _vp, prec, MPFR_RNDU);
142  else { // u in (u1, u2)
143  mpfr_set_prec(_vx, prec);
144  mpfr_add(_vx, _vp, _eps, MPFR_RNDN);
145  reject = Reject(_u, _vx, prec, MPFR_RNDU); // Could use _up too
146  }
147  if (reject < 0) break; // tried to reject but failed, so accept
148 
149  // Try for acceptance
150  if (uf <= u1)
151  reject = Reject(_up, _v, prec, MPFR_RNDD);
152  else if (uf >= u2)
153  reject = Reject(_u, _v, prec, MPFR_RNDD);
154  else { // u in (u2, u2)
155  mpfr_sub(_vx, _v, _eps, MPFR_RNDN);
156  reject = Reject(_u, _vx, prec, MPFR_RNDD); // Could use _up too
157  }
158  if (reject > 0) break; // tried to accept but failed, so reject
159 
160  prec = Refine(r, prec); // still can't decide, to refine
161  }
162  if (reject > 0) continue; // reject, back to outer loop
163  }
164  // Now evaluate v/u to the necessary precision
165  mpfr_prec_t prec0 = mpfr_get_prec (val);
166  // while (prec < prec0 + prec_guard) prec = Refine(r, prec);
167  if (prec < prec0 + prec_guard)
168  prec = Refine(r, prec,
169  (prec0 + prec_guard - prec + chunk_ - 1) / chunk_);
170  mpfr_set_prec(_x1, prec0);
171  mpfr_set_prec(_x2, prec0);
172  int flag;
173  while (true) {
174  int
175  f1 = mpfr_div(_x1, _v, _up, round), // min slope
176  f2 = mpfr_div(_x2, _vp, _u, round); // max slope
177  if (f1 == f2 && mpfr_equal_p(_x1, _x2)) {
178  flag = f1;
179  break;
180  }
181  prec = Refine(r, prec);
182  }
183  mpz_urandomb(_ui, r, 1);
184  if (mpz_tstbit(_ui, 0)) {
185  flag = -flag;
186  mpfr_neg(val, _x1, MPFR_RNDN);
187  } else
188  mpfr_set(val, _x1, MPFR_RNDN);
189  // std::cerr << uf << " " << vf << " " << Q << "\n";
190  return flag;
191  }
192  }
193  private:
194  // disable copy constructor and assignment operator
195  MPFRNormalR(const MPFRNormalR&);
196  MPFRNormalR& operator=(const MPFRNormalR&);
197  // Refine the random square
198  mpfr_prec_t Refine(gmp_randstate_t r, mpfr_prec_t prec, long num = 1)
199  const {
200  if (num <= 0) return prec;
201  // Use _vx as scratch
202  prec += num * chunk_;
203  mpfr_div_2ui(_eps, _eps, num * chunk_, MPFR_RNDN);
204 
205  mpz_urandomb(_ui, r, num * chunk_);
206  mpfr_set_prec(_up, prec);
207  mpfr_set_z_2exp(_up, _ui, -prec, MPFR_RNDN);
208  mpfr_set_prec(_vx, prec);
209  mpfr_add(_vx, _u, _up, MPFR_RNDN);
210  mpfr_swap(_u, _vx); // u = vx
211  mpfr_add(_up, _u, _eps, MPFR_RNDN);
212 
213  mpz_urandomb(_vi, r, num * chunk_);
214  mpfr_set_prec(_vp, prec);
215  mpfr_set_z_2exp(_vp, _vi, -prec, MPFR_RNDN);
216  mpfr_set_prec(_vx, prec);
217  mpfr_add(_vx, _v, _vp, MPFR_RNDN);
218  mpfr_swap(_v, _vx); // v = vx
219  mpfr_add(_vp, _v, _eps, MPFR_RNDN);
220 
221  return prec;
222  }
223  // Evaluate the sign of the rejection condition v^2 + 4*u^2*log(u)
224  int Reject(mpfr_t u, mpfr_t v, mpfr_prec_t prec, mpfr_rnd_t round) const {
225  // Use x1, x2 as scratch
226  mpfr_set_prec(_x1, prec);
227 
228  mpfr_log(_x1, u, round);
229  mpfr_mul(_x1, _x1, u, round); // Important to do the multiplications in
230  mpfr_mul(_x1, _x1, u, round); // this order so that rounding works right.
231  mpfr_mul_2ui(_x1, _x1, 2u, round); // 4*u^2*log(u)
232 
233  mpfr_set_prec(_x2, prec);
234  mpfr_mul(_x2, v, v, round); // v^2
235 
236  mpfr_add(_x1, _x1, _x2, round); // v^2 + 4*u^2*log(u)
237 
238  return mpfr_sgn(_x1);
239  }
240  mutable mpz_t _ui;
241  mutable mpz_t _vi;
242  mutable mpfr_t _eps;
243  mutable mpfr_t _u;
244  mutable mpfr_t _v;
245  mutable mpfr_t _up;
246  mutable mpfr_t _vp;
247  mutable mpfr_t _vx;
248  mutable mpfr_t _x1;
249  mutable mpfr_t _x2;
250  };
251 
252 } // namespace RandomLib
253 
254 #endif // HAVE_MPFR
255 #endif // RANDOMLIB_MPFRNORMALR_HPP
int operator()(mpfr_t val, gmp_randstate_t r, mpfr_rnd_t round) const
Definition: MPFRNormalR.hpp:91
The normal distribution for MPFR (ratio method).
Definition: MPFRNormalR.hpp:43