MPFRNormalK.hpp
5.08 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
/**
* \file MPFRNormalK.hpp
* \brief Header for MPFRNormalK
*
* Sampling exactly from the normal distribution for MPFR.
*
* Copyright (c) Charles Karney (2012) <charles@karney.com> and licensed under
* the MIT/X11 License. For more information, see
* http://randomlib.sourceforge.net/
**********************************************************************/
#if !defined(RANDOMLIB_MPFRNORMALK_HPP)
#define RANDOMLIB_MPFRNORMALK_HPP 1
#include <algorithm> // for max
#include <RandomLib/MPFRRandom.hpp>
#include <RandomLib/MPFRExponential.hpp>
#if HAVE_MPFR || defined(DOXYGEN)
namespace RandomLib {
/**
* \brief The normal distribution for MPFR (Kahn algorithm).
*
* This class is <b>DEPRECATED</b>. It is included for illustrative purposes
* only. The MPFRNormal class provides a somewhat more efficient method for
* sampling from the normal distribution.
*
* Refs:
* - H. Kahn, Rand Report RM-1237-AEC, p. 41 (1954).
* - M. Abramowitz and I. A. Stegun, p. 953, Sec. 26.8.6.a(4) (1964).
* .
* N.B. Damien Stehle' drew my attention to this algorithm as a useful way to
* compute normal deviates exactly.
*
* This class uses mutable private objects. So a single MPFRNormalK object
* cannot safely be used by multiple threads. In a multi-processing
* environment, each thread should use a thread-specific MPFRNormalK object.
*
* @tparam bits the number of bits in each digit.
**********************************************************************/
template<int bits = 32> class MPFRNormalK {
public:
/**
* Initialize the MPFRNormalK object.
**********************************************************************/
MPFRNormalK()
{ mpfr_init2(_xf, MPFR_PREC_MIN); mpfr_init2(_zf, MPFR_PREC_MIN); }
/**
* Destroy the MPFRNormalK object.
**********************************************************************/
~MPFRNormalK()
{ mpfr_clear(_zf); mpfr_clear(_xf); }
/**
* Sample from the normal distribution with mean 0 and variance 1 returning
* a MPFRRandom.
*
* @param[out] t the MPFRRandom result.
* @param[in,out] r a GMP random generator.
**********************************************************************/
void operator()(MPFRRandom<bits>& t, gmp_randstate_t r) const
{ Compute(r); _x.swap(t); }
/**
* Sample from the normal distribution with mean 0 and variance 1.
*
* @param[out] val the sample from the normal distribution
* @param[in,out] r a GMP random generator.
* @param[in] round the rounding direction.
* @return the MPFR ternary result (±1 if val is larger/smaller than
* the exact sample).
**********************************************************************/
int operator()(mpfr_t val, gmp_randstate_t r, mpfr_rnd_t round) const
{ Compute(r); return _x(val, r, round); }
private:
// disable copy constructor and assignment operator
MPFRNormalK(const MPFRNormalK&);
MPFRNormalK& operator=(const MPFRNormalK&);
void Compute(gmp_randstate_t r) const {
// The algorithm is sample x and z from the exponential distribution; if
// (x-1)^2 < 2*z, return (random sign)*x; otherwise repeat. Probability
// of acceptance is sqrt(pi/2) * exp(-1/2) = 0.7602.
while (true) {
_edist(_x, r);
_edist(_z, r);
for (mp_size_t k = 1; ; ++k) {
_x.ExpandTo(r, k - 1);
_z.ExpandTo(r, k - 1);
mpfr_prec_t prec = (std::max)(mpfr_prec_t(MPFR_PREC_MIN), k * bits);
mpfr_set_prec(_xf, prec);
mpfr_set_prec(_zf, prec);
// Try for acceptance first; so compute upper limit on (y-1)^2 and
// lower limit on 2*z.
if (_x.UInteger() == 0) {
_x(_xf, MPFR_RNDD);
mpfr_ui_sub(_xf, 1u, _xf, MPFR_RNDU);
} else {
_x(_xf, MPFR_RNDU);
mpfr_sub_ui(_xf, _xf, 1u, MPFR_RNDU);
}
mpfr_sqr(_xf, _xf, MPFR_RNDU);
_z(_zf, MPFR_RNDD);
mpfr_mul_2ui(_zf, _zf, 1u, MPFR_RNDD);
if (mpfr_cmp(_xf, _zf) < 0) { // (y-1)^2 < 2*z, so accept
if (_x.Boolean(r)) _x.Negate(); // include a random sign
return;
}
// Try for rejection; so compute lower limit on (y-1)^2 and upper
// limit on 2*z.
if (_x.UInteger() == 0) {
_x(_xf, MPFR_RNDU);
mpfr_ui_sub(_xf, 1u, _xf, MPFR_RNDD);
} else {
_x(_xf, MPFR_RNDD);
mpfr_sub_ui(_xf, _xf, 1u, MPFR_RNDD);
}
mpfr_sqr(_xf, _xf, MPFR_RNDD);
_z(_zf, MPFR_RNDU);
mpfr_mul_2ui(_zf, _zf, 1u, MPFR_RNDU);
if (mpfr_cmp(_xf, _zf) > 0) // (y-1)^2 > 2*z, so reject
break;
// Otherwise repeat with more precision
}
// Reject and start over with a new y and z
}
}
mutable MPFRRandom<bits> _x;
mutable MPFRRandom<bits> _z;
mutable mpfr_t _xf;
mutable mpfr_t _zf;
const MPFRExponential<bits> _edist;
};
} // namespace RandomLib
#endif // HAVE_MPFR
#endif // RANDOMLIB_MPFRNORMALK_HPP