blitz  Version 1.0.2
uniform.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 /***************************************************************************
3  * random/uniform.h Uniform IRNG wrapper class
4  *
5  * $Id$
6  *
7  * Copyright (C) 1997-2011 Todd Veldhuizen <tveldhui@acm.org>
8  *
9  * This file is a part of Blitz.
10  *
11  * Blitz is free software: you can redistribute it and/or modify
12  * it under the terms of the GNU Lesser General Public License
13  * as published by the Free Software Foundation, either version 3
14  * of the License, or (at your option) any later version.
15  *
16  * Blitz is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  * GNU Lesser General Public License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with Blitz. If not, see <http://www.gnu.org/licenses/>.
23  *
24  * Suggestions: blitz-devel@lists.sourceforge.net
25  * Bugs: blitz-support@lists.sourceforge.net
26  *
27  * For more information, please see the Blitz++ Home Page:
28  * https://sourceforge.net/projects/blitz/
29  *
30  ***************************************************************************/
31 
32 #ifndef BZ_RANDOM_UNIFORM_H
33 #define BZ_RANDOM_UNIFORM_H
34 
35 #include <random/default.h>
36 
37 #ifndef FLT_MANT_DIG
38  #include <float.h>
39 #endif
40 
41 namespace ranlib {
42 
43 /*****************************************************************************
44  * UniformClosedOpen generator: uniform random numbers in [0,1).
45  *****************************************************************************/
46 
47 template<typename T = defaultType, typename IRNG = defaultIRNG,
48  typename stateTag = defaultState>
50 
51 // These constants are 1/2^32, 1/2^64, 1/2^96, 1/2^128
52 const long double norm32open = .2328306436538696289062500000000000000000E-9L;
53 const long double norm64open = .5421010862427522170037264004349708557129E-19L;
54 const long double norm96open = .1262177448353618888658765704452457967477E-28L;
55 const long double norm128open = .2938735877055718769921841343055614194547E-38L;
56 
57 
58 template<typename IRNG, typename stateTag>
59 class UniformClosedOpen<float,IRNG,stateTag>
60  : public IRNGWrapper<IRNG,stateTag>
61 {
62 
63 public:
64  typedef float T_numtype;
65 
67  UniformClosedOpen(unsigned int i) :
68  IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
69 
70  float random()
71  {
72 #if FLT_MANT_DIG > 96
73  #if (FLT_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
74  #error RNG code assumes float mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
75  #endif
76  IRNG_int i1 = this->irng_.random();
77  IRNG_int i2 = this->irng_.random();
78  IRNG_int i3 = this->irng_.random();
79  IRNG_int i4 = this->irng_.random();
80  return i1 * norm128open + i2 * norm96open + i3 * norm64open
81  + i4 * norm32open;
82 #elif FLT_MANT_DIG > 64
83  IRNG_int i1 = this->irng_.random();
84  IRNG_int i2 = this->irng_.random();
85  IRNG_int i3 = this->irng_.random();
86  return i1 * norm96open + i2 * norm64open + i3 * norm32open;
87 #elif FLT_MANT_DIG > 32
88  IRNG_int i1 = this->irng_.random();
89  IRNG_int i2 = this->irng_.random();
90  return i1 * norm64open + i2 * norm32open;
91 #else
92  IRNG_int i1 = this->irng_.random();
93  return i1 * norm32open;
94 #endif
95  }
96 
97  float getUniform()
98  { return random(); }
99 };
100 
101 template<typename IRNG, typename stateTag>
102 class UniformClosedOpen<double,IRNG,stateTag>
103  : public IRNGWrapper<IRNG,stateTag>
104 {
105 public:
106  typedef double T_numtype;
107 
109  UniformClosedOpen(unsigned int i) :
110  IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
111 
112  double random()
113  {
114 #if DBL_MANT_DIG > 96
115  #if (DBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
116  #error RNG code assumes 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 DBL_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 DBL_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 
141  double getUniform() { return random(); }
142 };
143 
144 template<typename IRNG, typename stateTag>
145 class UniformClosedOpen<long double,IRNG,stateTag>
146  : public IRNGWrapper<IRNG,stateTag>
147 {
148 public:
149  typedef long double T_numtype;
150 
152  UniformClosedOpen(unsigned int i) :
153  IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
154 
155  long double random()
156  {
157 #if LDBL_MANT_DIG > 96
158  #if (LDBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
159  #error RNG code assumes long double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
160 #endif
161 
162  IRNG_int i1 = this->irng_.random();
163  IRNG_int i2 = this->irng_.random();
164  IRNG_int i3 = this->irng_.random();
165  IRNG_int i4 = this->irng_.random();
166  return i1 * norm128open + i2 * norm96open + i3 * norm64open
167  + i4 * norm32open;
168 #elif LDBL_MANT_DIG > 64
169  IRNG_int i1 = this->irng_.random();
170  IRNG_int i2 = this->irng_.random();
171  IRNG_int i3 = this->irng_.random();
172  return i1 * norm96open + i2 * norm64open + i3 * norm32open;
173 #elif LDBL_MANT_DIG > 32
174  IRNG_int i1 = this->irng_.random();
175  IRNG_int i2 = this->irng_.random();
176  return i1 * norm64open + i2 * norm32open;
177 #else
178  IRNG_int i1 = this->irng_.random();
179  return i1 * norm32open;
180 #endif
181  }
182 
183  long double getUniform() { return random(); }
184 };
185 
186 // For people who don't care about open or closed: just give them
187 // ClosedOpen (this is the fastest).
188 
189 template<class T, typename IRNG = defaultIRNG,
190  typename stateTag = defaultState>
191 class Uniform : public UniformClosedOpen<T,IRNG,stateTag>
192 {
193 public:
194  Uniform() {};
195  Uniform(unsigned int i) :
196  UniformClosedOpen<T,IRNG,stateTag>::UniformClosedOpen(i) {}
197  };
198 
199 /*****************************************************************************
200  * UniformClosed generator: uniform random numbers in [0,1].
201  *****************************************************************************/
202 
203 // This constant is 1/(2^32-1)
204 const long double norm32closed = .2328306437080797375431469961868475648078E-9L;
205 
206 // These constants are 2^32/(2^64-1) and 1/(2^64-1), respectively.
207 
208 const long double norm64closed1 =
209  .23283064365386962891887177448353618888727188481031E-9L;
210 const long double norm64closed2 =
211  .54210108624275221703311375920552804341370213034169E-19L;
212 
213 // These constants are 2^64/(2^96-1), 2^32/(2^96-1), and 1/(2^96-1)
214 const long double norm96closed1 = .2328306436538696289062500000029387358771E-9L;
215 const long double norm96closed2 =
216  .5421010862427522170037264004418131333707E-19L;
217 const long double norm96closed3 =
218  .1262177448353618888658765704468388886588E-28L;
219 
220 // These constants are 2^96/(2^128-1), 2^64/(2^128-1), 2^32/(2^128-1) and
221 // 1/(2^128-1).
222 const long double norm128clos1 = .2328306436538696289062500000000000000007E-9L;
223 const long double norm128clos2 = .5421010862427522170037264004349708557145E-19L;
224 const long double norm128clos3 = .1262177448353618888658765704452457967481E-28L;
225 const long double norm128clos4 = .2938735877055718769921841343055614194555E-38L;
226 
227 
228 template<typename T = double, typename IRNG = defaultIRNG,
229  typename stateTag = defaultState>
230 class UniformClosed { };
231 
232 template<typename IRNG, typename stateTag>
233 class UniformClosed<float,IRNG,stateTag> : public IRNGWrapper<IRNG,stateTag> {
234 
235 public:
236  typedef float T_numtype;
237 
239  UniformClosed(unsigned int i) :
240  IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
241 
242  float random()
243  {
244 #if FLTMANT_DIG > 96
245  #if (FLT_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
246  #error RNG code assumes float mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
247  #endif
248  IRNG_int i1 = this->irng_.random();
249  IRNG_int i2 = this->irng_.random();
250  IRNG_int i3 = this->irng_.random();
251  IRNG_int i4 = this->irng_.random();
252 
253  return i1 * norm128clos1 + i2 * norm128clos2
254  + i3 * norm128clos3 + i4 * norm128clos4;
255 #elif FLT_MANT_DIG > 64
256  IRNG_int i1 = this->irng_.random();
257  IRNG_int i2 = this->irng_.random();
258  IRNG_int i3 = this->irng_.random();
259 
260  return i1 * norm96closed1 + i2 * norm96closed2
261  + i3 * norm96closed3;
262 #elif FLT_MANT_DIG > 32
263  IRNG_int i1 = this->irng_.random();
264  IRNG_int i2 = this->irng_.random();
265 
266  return i1 * norm64closed1 + i2 * norm64closed2;
267 #else
268  IRNG_int i = this->irng_.random();
269  return i * norm32closed;
270 #endif
271 
272  }
273 
274  float getUniform()
275  { return random(); }
276 };
277 
278 template<typename IRNG, typename stateTag>
279 class UniformClosed<double,IRNG,stateTag> : public IRNGWrapper<IRNG,stateTag> {
280 
281 public:
282  typedef double T_numtype;
283 
285  UniformClosed(unsigned int i) :
286  IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
287 
288  double random()
289  {
290 #if DBL_MANT_DIG > 96
291  #if (DBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
292  #error RNG code assumes double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
293  #endif
294  IRNG_int i1 = this->irng_.random();
295  IRNG_int i2 = this->irng_.random();
296  IRNG_int i3 = this->irng_.random();
297  IRNG_int i4 = this->irng_.random();
298 
299  return i1 * norm128clos1 + i2 * norm128clos2
300  + i3 * norm128clos3 + i4 * norm128clos4;
301 #elif DBL_MANT_DIG > 64
302  IRNG_int i1 = this->irng_.random();
303  IRNG_int i2 = this->irng_.random();
304  IRNG_int i3 = this->irng_.random();
305 
306  return i1 * norm96closed1 + i2 * norm96closed2
307  + i3 * norm96closed3;
308 #elif DBL_MANT_DIG > 32
309  IRNG_int i1 = this->irng_.random();
310  IRNG_int i2 = this->irng_.random();
311 
312  return i1 * norm64closed1 + i2 * norm64closed2;
313 #else
314  IRNG_int i = this->irng_.random();
315  return i * norm32closed;
316 #endif
317 
318  }
319 
320  double getUniform()
321  { return random(); }
322 };
323 
324 template<typename IRNG, typename stateTag>
325 class UniformClosed<long double,IRNG,stateTag>
326  : public IRNGWrapper<IRNG,stateTag> {
327 
328 public:
329  typedef long double T_numtype;
330 
332  UniformClosed(unsigned int i) :
333  IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
334 
335  long double random()
336  {
337 #if LDBL_MANT_DIG > 96
338  #if (LDBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
339  #error RNG code assumes long double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
340  #endif
341  IRNG_int i1 = this->irng_.random();
342  IRNG_int i2 = this->irng_.random();
343  IRNG_int i3 = this->irng_.random();
344  IRNG_int i4 = this->irng_.random();
345 
346  return i1 * norm128clos1 + i2 * norm128clos2
347  + i3 * norm128clos3 + i4 * norm128clos4;
348 #elif LDBL_MANT_DIG > 64
349  IRNG_int i1 = this->irng_.random();
350  IRNG_int i2 = this->irng_.random();
351  IRNG_int i3 = this->irng_.random();
352 
353  return i1 * norm96closed1 + i2 * norm96closed2
354  + i3 * norm96closed3;
355 #elif LDBL_MANT_DIG > 32
356  IRNG_int i1 = this->irng_.random();
357  IRNG_int i2 = this->irng_.random();
358 
359  return i1 * norm64closed1 + i2 * norm64closed2;
360 #else
361  IRNG_int i = this->irng_.random();
362  return i * norm32closed;
363 #endif
364  }
365 
366  long double getUniform()
367  { return random(); }
368 
369 };
370 
371 /*****************************************************************************
372  * UniformOpen generator: uniform random numbers in (0,1).
373  *****************************************************************************/
374 
375 template<typename T = double, typename IRNG = defaultIRNG,
376  typename stateTag = defaultState>
377 class UniformOpen : public UniformClosedOpen<T,IRNG,stateTag> {
378  public:
379  typedef T T_numtype;
380 
382  UniformOpen(unsigned int i) :
383  UniformClosedOpen<T,IRNG,stateTag>(i) {};
384 
385  T random()
386  {
387  // Turn a [0,1) into a (0,1) interval by weeding out
388  // any zeros.
389  T x;
390  do {
392  } while (x == 0.0L);
393 
394  return x;
395  }
396 
398  { return random(); }
399 
400 };
401 
402 /*****************************************************************************
403  * UniformOpenClosed generator: uniform random numbers in (0,1]
404  *****************************************************************************/
405 
406 template<typename T = double, typename IRNG = defaultIRNG,
407  typename stateTag = defaultState>
408 class UniformOpenClosed : public UniformClosedOpen<T,IRNG,stateTag> {
409 
410 public:
411  typedef T T_numtype;
412 
414  UniformOpenClosed(unsigned int i) :
415  UniformClosedOpen<T,IRNG,stateTag>(i) {};
416 
417 
418  T random()
419  {
420  // Antithetic value: taking 1-X where X is [0,1) results
421  // in a (0,1] distribution.
423  }
424 
426  { return random(); }
427 };
428 
429 }
430 
431 #endif // BZ_RANDOM_UNIFORM_H
float random()
Definition: uniform.h:242
Definition: beta.h:50
double T_numtype
Definition: uniform.h:282
UniformClosedOpen(unsigned int i)
Definition: uniform.h:109
float getUniform()
Definition: uniform.h:97
T random()
Definition: uniform.h:385
Definition: uniform.h:230
_bz_global blitz::IndexPlaceholder< 0 > i
Definition: indexexpr.h:256
const long double norm128clos1
Definition: uniform.h:222
const long double norm96closed3
Definition: uniform.h:217
long double getUniform()
Definition: uniform.h:366
long double getUniform()
Definition: uniform.h:183
const long double norm32closed
Definition: uniform.h:204
const long double norm64closed1
Definition: uniform.h:208
Definition: uniform.h:191
T getUniform()
Definition: uniform.h:397
UniformOpenClosed(unsigned int i)
Definition: uniform.h:414
Definition: default.h:68
const long double norm64closed2
Definition: uniform.h:210
const long double norm32open
Definition: uniform.h:52
float getUniform()
Definition: uniform.h:274
const long double norm96closed2
Definition: uniform.h:215
float defaultType
Definition: default.h:46
UniformOpen(unsigned int i)
Definition: uniform.h:382
UniformClosed(unsigned int i)
Definition: uniform.h:239
UniformOpen()
Definition: uniform.h:381
T getUniform()
Definition: uniform.h:425
long double random()
Definition: uniform.h:335
const long double norm128clos2
Definition: uniform.h:223
UniformOpenClosed()
Definition: uniform.h:413
const long double norm128clos3
Definition: uniform.h:224
T T_numtype
Definition: uniform.h:411
unsigned int IRNG_int
Definition: default.h:57
long double T_numtype
Definition: uniform.h:329
const long double norm96closed1
Definition: uniform.h:214
Definition: default.h:53
MersenneTwister defaultIRNG
Definition: default.h:120
long double T_numtype
Definition: uniform.h:149
double getUniform()
Definition: uniform.h:320
Definition: uniform.h:377
Definition: uniform.h:408
float T_numtype
Definition: uniform.h:236
const long double norm64open
Definition: uniform.h:53
const long double norm96open
Definition: uniform.h:54
UniformClosed(unsigned int i)
Definition: uniform.h:285
double getUniform()
Definition: uniform.h:141
Uniform()
Definition: uniform.h:194
const long double norm128clos4
Definition: uniform.h:225
T T_numtype
Definition: uniform.h:379
UniformClosedOpen(unsigned int i)
Definition: uniform.h:152
UniformClosedOpen(unsigned int i)
Definition: uniform.h:67
UniformClosed(unsigned int i)
Definition: uniform.h:332
double random()
Definition: uniform.h:288
sharedState defaultState
Definition: default.h:55
const long double norm128open
Definition: uniform.h:55
Uniform(unsigned int i)
Definition: uniform.h:195
T random()
Definition: uniform.h:418
double random()
Definition: uniform.h:112
Definition: uniform.h:49
long double random()
Definition: uniform.h:155