blitz  Version 0.9
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
uniform.h
Go to the documentation of this file.
1 #ifndef BZ_RANDOM_UNIFORM_H
2 #define BZ_RANDOM_UNIFORM_H
3 
4 #include <random/default.h>
5 
6 #ifndef FLT_MANT_DIG
7  #include <float.h>
8 #endif
9 
10 BZ_NAMESPACE(ranlib)
11 
12 /*****************************************************************************
13  * UniformClosedOpen generator: uniform random numbers in [0,1).
14  *****************************************************************************/
15 
16 template<typename T = defaultType, typename IRNG = defaultIRNG,
17  typename stateTag = defaultState>
19 
20 // These constants are 1/2^32, 1/2^64, 1/2^96, 1/2^128
21 const long double norm32open = .2328306436538696289062500000000000000000E-9L;
22 const long double norm64open = .5421010862427522170037264004349708557129E-19L;
23 const long double norm96open = .1262177448353618888658765704452457967477E-28L;
24 const long double norm128open = .2938735877055718769921841343055614194547E-38L;
25 
26 
27 template<typename IRNG, typename stateTag>
28 class UniformClosedOpen<float,IRNG,stateTag>
29  : public IRNGWrapper<IRNG,stateTag>
30 {
31 
32 public:
33  typedef float T_numtype;
34 
35  float random()
36  {
37 #if FLT_MANT_DIG > 96
38  #if (FLT_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
39  #error RNG code assumes float mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
40  #endif
41  IRNG_int i1 = this->irng_.random();
42  IRNG_int i2 = this->irng_.random();
43  IRNG_int i3 = this->irng_.random();
44  IRNG_int i4 = this->irng_.random();
45  return i1 * norm128open + i2 * norm96open + i3 * norm64open
46  + i4 * norm32open;
47 #elif FLT_MANT_DIG > 64
48  IRNG_int i1 = this->irng_.random();
49  IRNG_int i2 = this->irng_.random();
50  IRNG_int i3 = this->irng_.random();
51  return i1 * norm96open + i2 * norm64open + i3 * norm32open;
52 #elif FLT_MANT_DIG > 32
53  IRNG_int i1 = this->irng_.random();
54  IRNG_int i2 = this->irng_.random();
55  return i1 * norm64open + i2 * norm32open;
56 #else
57  IRNG_int i1 = this->irng_.random();
58  return i1 * norm32open;
59 #endif
60  }
61 
62  float getUniform()
63  { return random(); }
64 };
65 
66 template<typename IRNG, typename stateTag>
67 class UniformClosedOpen<double,IRNG,stateTag>
68  : public IRNGWrapper<IRNG,stateTag>
69 {
70 public:
71  typedef double T_numtype;
72 
73  double random()
74  {
75 #if DBL_MANT_DIG > 96
76  #if (DBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
77  #error RNG code assumes double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
78  #endif
79 
80  IRNG_int i1 = this->irng_.random();
81  IRNG_int i2 = this->irng_.random();
82  IRNG_int i3 = this->irng_.random();
83  IRNG_int i4 = this->irng_.random();
84  return i1 * norm128open + i2 * norm96open + i3 * norm64open
85  + i4 * norm32open;
86 #elif DBL_MANT_DIG > 64
87  IRNG_int i1 = this->irng_.random();
88  IRNG_int i2 = this->irng_.random();
89  IRNG_int i3 = this->irng_.random();
90  return i1 * norm96open + i2 * norm64open + i3 * norm32open;
91 #elif DBL_MANT_DIG > 32
92  IRNG_int i1 = this->irng_.random();
93  IRNG_int i2 = this->irng_.random();
94  return i1 * norm64open + i2 * norm32open;
95 #else
96  IRNG_int i1 = this->irng_.random();
97  return i1 * norm32open;
98 #endif
99 
100  }
101 
102  double getUniform() { return random(); }
103 };
104 
105 template<typename IRNG, typename stateTag>
106 class UniformClosedOpen<long double,IRNG,stateTag>
107  : public IRNGWrapper<IRNG,stateTag>
108 {
109 public:
110  typedef long double T_numtype;
111 
112  long double random()
113  {
114 #if LDBL_MANT_DIG > 96
115  #if (LDBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
116  #error RNG code assumes long double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
117 #endif
118 
119  IRNG_int i1 = this->irng_.random();
120  IRNG_int i2 = this->irng_.random();
121  IRNG_int i3 = this->irng_.random();
122  IRNG_int i4 = this->irng_.random();
123  return i1 * norm128open + i2 * norm96open + i3 * norm64open
124  + i4 * norm32open;
125 #elif LDBL_MANT_DIG > 64
126  IRNG_int i1 = this->irng_.random();
127  IRNG_int i2 = this->irng_.random();
128  IRNG_int i3 = this->irng_.random();
129  return i1 * norm96open + i2 * norm64open + i3 * norm32open;
130 #elif LDBL_MANT_DIG > 32
131  IRNG_int i1 = this->irng_.random();
132  IRNG_int i2 = this->irng_.random();
133  return i1 * norm64open + i2 * norm32open;
134 #else
135  IRNG_int i1 = this->irng_.random();
136  return i1 * norm32open;
137 #endif
138  }
139 
140  long double getUniform() { return random(); }
141 };
142 
143 // For people who don't care about open or closed: just give them
144 // ClosedOpen (this is the fastest).
145 
146 template<class T, typename IRNG = defaultIRNG,
147  typename stateTag = defaultState>
148 class Uniform : public UniformClosedOpen<T,IRNG,stateTag>
149 { };
150 
151 /*****************************************************************************
152  * UniformClosed generator: uniform random numbers in [0,1].
153  *****************************************************************************/
154 
155 // This constant is 1/(2^32-1)
156 const long double norm32closed = .2328306437080797375431469961868475648078E-9L;
157 
158 // These constants are 2^32/(2^64-1) and 1/(2^64-1), respectively.
159 
160 const long double norm64closed1 =
161  .23283064365386962891887177448353618888727188481031E-9L;
162 const long double norm64closed2 =
163  .54210108624275221703311375920552804341370213034169E-19L;
164 
165 // These constants are 2^64/(2^96-1), 2^32/(2^96-1), and 1/(2^96-1)
166 const long double norm96closed1 = .2328306436538696289062500000029387358771E-9L;
167 const long double norm96closed2 =
168  .5421010862427522170037264004418131333707E-19L;
169 const long double norm96closed3 =
170  .1262177448353618888658765704468388886588E-28L;
171 
172 // These constants are 2^96/(2^128-1), 2^64/(2^128-1), 2^32/(2^128-1) and
173 // 1/(2^128-1).
174 const long double norm128clos1 = .2328306436538696289062500000000000000007E-9L;
175 const long double norm128clos2 = .5421010862427522170037264004349708557145E-19L;
176 const long double norm128clos3 = .1262177448353618888658765704452457967481E-28L;
177 const long double norm128clos4 = .2938735877055718769921841343055614194555E-38L;
178 
179 
180 template<typename T = double, typename IRNG = defaultIRNG,
181  typename stateTag = defaultState>
182 class UniformClosed { };
183 
184 template<typename IRNG, typename stateTag>
185 class UniformClosed<float,IRNG,stateTag> : public IRNGWrapper<IRNG,stateTag> {
186 
187 public:
188  typedef float T_numtype;
189 
190  float random()
191  {
192 #if FLTMANT_DIG > 96
193  #if (FLT_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
194  #error RNG code assumes float mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
195  #endif
196  IRNG_int i1 = this->irng_.random();
197  IRNG_int i2 = this->irng_.random();
198  IRNG_int i3 = this->irng_.random();
199  IRNG_int i4 = this->irng_.random();
200 
201  return i1 * norm128clos1 + i2 * norm128clos2
202  + i3 * norm128clos3 + i4 * norm128clos4;
203 #elif FLT_MANT_DIG > 64
204  IRNG_int i1 = this->irng_.random();
205  IRNG_int i2 = this->irng_.random();
206  IRNG_int i3 = this->irng_.random();
207 
208  return i1 * norm96closed1 + i2 * norm96closed2
209  + i3 * norm96closed3;
210 #elif FLT_MANT_DIG > 32
211  IRNG_int i1 = this->irng_.random();
212  IRNG_int i2 = this->irng_.random();
213 
214  return i1 * norm64closed1 + i2 * norm64closed2;
215 #else
216  IRNG_int i = this->irng_.random();
217  return i * norm32closed;
218 #endif
219 
220  }
221 
222  float getUniform()
223  { return random(); }
224 };
225 
226 template<typename IRNG, typename stateTag>
227 class UniformClosed<double,IRNG,stateTag> : public IRNGWrapper<IRNG,stateTag> {
228 
229 public:
230  typedef double T_numtype;
231 
232  double random()
233  {
234 #if DBL_MANT_DIG > 96
235  #if (DBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
236  #error RNG code assumes double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
237  #endif
238  IRNG_int i1 = this->irng_.random();
239  IRNG_int i2 = this->irng_.random();
240  IRNG_int i3 = this->irng_.random();
241  IRNG_int i4 = this->irng_.random();
242 
243  return i1 * norm128clos1 + i2 * norm128clos2
244  + i3 * norm128clos3 + i4 * norm128clos4;
245 #elif DBL_MANT_DIG > 64
246  IRNG_int i1 = this->irng_.random();
247  IRNG_int i2 = this->irng_.random();
248  IRNG_int i3 = this->irng_.random();
249 
250  return i1 * norm96closed1 + i2 * norm96closed2
251  + i3 * norm96closed3;
252 #elif DBL_MANT_DIG > 32
253  IRNG_int i1 = this->irng_.random();
254  IRNG_int i2 = this->irng_.random();
255 
256  return i1 * norm64closed1 + i2 * norm64closed2;
257 #else
258  IRNG_int i = this->irng_.random();
259  return i * norm32closed;
260 #endif
261 
262  }
263 
264  double getUniform()
265  { return random(); }
266 };
267 
268 template<typename IRNG, typename stateTag>
269 class UniformClosed<long double,IRNG,stateTag>
270  : public IRNGWrapper<IRNG,stateTag> {
271 
272 public:
273  typedef long double T_numtype;
274 
275  long double random()
276  {
277 #if LDBL_MANT_DIG > 96
278  #if (LDBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
279  #error RNG code assumes long double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
280  #endif
281  IRNG_int i1 = this->irng_.random();
282  IRNG_int i2 = this->irng_.random();
283  IRNG_int i3 = this->irng_.random();
284  IRNG_int i4 = this->irng_.random();
285 
286  return i1 * norm128clos1 + i2 * norm128clos2
287  + i3 * norm128clos3 + i4 * norm128clos4;
288 #elif LDBL_MANT_DIG > 64
289  IRNG_int i1 = this->irng_.random();
290  IRNG_int i2 = this->irng_.random();
291  IRNG_int i3 = this->irng_.random();
292 
293  return i1 * norm96closed1 + i2 * norm96closed2
294  + i3 * norm96closed3;
295 #elif LDBL_MANT_DIG > 32
296  IRNG_int i1 = this->irng_.random();
297  IRNG_int i2 = this->irng_.random();
298 
299  return i1 * norm64closed1 + i2 * norm64closed2;
300 #else
301  IRNG_int i = this->irng_.random();
302  return i * norm32closed;
303 #endif
304  }
305 
306  long double getUniform()
307  { return random(); }
308 
309 };
310 
311 /*****************************************************************************
312  * UniformOpen generator: uniform random numbers in (0,1).
313  *****************************************************************************/
314 
315 template<typename T = double, typename IRNG = defaultIRNG,
316  typename stateTag = defaultState>
317 class UniformOpen : public UniformClosedOpen<T,IRNG,stateTag> {
318  public:
319  typedef T T_numtype;
320 
321  T random()
322  {
323  // Turn a [0,1) into a (0,1) interval by weeding out
324  // any zeros.
325  T x;
326  do {
328  } while (x == 0.0L);
329 
330  return x;
331  }
332 
334  { return random(); }
335 
336 };
337 
338 /*****************************************************************************
339  * UniformOpenClosed generator: uniform random numbers in (0,1]
340  *****************************************************************************/
341 
342 template<typename T = double, typename IRNG = defaultIRNG,
343  typename stateTag = defaultState>
344 class UniformOpenClosed : public UniformClosedOpen<T,IRNG,stateTag> {
345 
346 public:
347  typedef T T_numtype;
348 
349  T random()
350  {
351  // Antithetic value: taking 1-X where X is [0,1) results
352  // in a (0,1] distribution.
354  }
355 
357  { return random(); }
358 };
359 
361 
362 #endif // BZ_RANDOM_UNIFORM_H