blitz Version 1.0.2
Loading...
Searching...
No Matches
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
41namespace ranlib {
42
43/*****************************************************************************
44 * UniformClosedOpen generator: uniform random numbers in [0,1).
45 *****************************************************************************/
46
47template<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
52const long double norm32open = .2328306436538696289062500000000000000000E-9L;
53const long double norm64open = .5421010862427522170037264004349708557129E-19L;
54const long double norm96open = .1262177448353618888658765704452457967477E-28L;
55const long double norm128open = .2938735877055718769921841343055614194547E-38L;
56
57
58template<typename IRNG, typename stateTag>
59class UniformClosedOpen<float,IRNG,stateTag>
60 : public IRNGWrapper<IRNG,stateTag>
61{
62
63public:
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
101template<typename IRNG, typename stateTag>
102class UniformClosedOpen<double,IRNG,stateTag>
103 : public IRNGWrapper<IRNG,stateTag>
104{
105public:
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
144template<typename IRNG, typename stateTag>
145class UniformClosedOpen<long double,IRNG,stateTag>
146 : public IRNGWrapper<IRNG,stateTag>
147{
148public:
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
189template<class T, typename IRNG = defaultIRNG,
190 typename stateTag = defaultState>
191class Uniform : public UniformClosedOpen<T,IRNG,stateTag>
192{
193public:
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)
204const long double norm32closed = .2328306437080797375431469961868475648078E-9L;
205
206// These constants are 2^32/(2^64-1) and 1/(2^64-1), respectively.
207
208const long double norm64closed1 =
209 .23283064365386962891887177448353618888727188481031E-9L;
210const 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)
214const long double norm96closed1 = .2328306436538696289062500000029387358771E-9L;
215const long double norm96closed2 =
216 .5421010862427522170037264004418131333707E-19L;
217const 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).
222const long double norm128clos1 = .2328306436538696289062500000000000000007E-9L;
223const long double norm128clos2 = .5421010862427522170037264004349708557145E-19L;
224const long double norm128clos3 = .1262177448353618888658765704452457967481E-28L;
225const long double norm128clos4 = .2938735877055718769921841343055614194555E-38L;
226
227
228template<typename T = double, typename IRNG = defaultIRNG,
229 typename stateTag = defaultState>
231
232template<typename IRNG, typename stateTag>
233class UniformClosed<float,IRNG,stateTag> : public IRNGWrapper<IRNG,stateTag> {
234
235public:
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
275 { return random(); }
276};
277
278template<typename IRNG, typename stateTag>
279class UniformClosed<double,IRNG,stateTag> : public IRNGWrapper<IRNG,stateTag> {
280
281public:
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
324template<typename IRNG, typename stateTag>
325class UniformClosed<long double,IRNG,stateTag>
326 : public IRNGWrapper<IRNG,stateTag> {
327
328public:
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
375template<typename T = double, typename IRNG = defaultIRNG,
376 typename stateTag = defaultState>
377class 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
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
406template<typename T = double, typename IRNG = defaultIRNG,
407 typename stateTag = defaultState>
408class UniformOpenClosed : public UniformClosedOpen<T,IRNG,stateTag> {
409
410public:
411 typedef T T_numtype;
412
414 UniformOpenClosed(unsigned int i) :
415 UniformClosedOpen<T,IRNG,stateTag>(i) {};
416
417
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
Definition default.h:68
double getUniform()
Definition uniform.h:141
UniformClosedOpen(unsigned int i)
Definition uniform.h:109
UniformClosedOpen(unsigned int i)
Definition uniform.h:67
long double T_numtype
Definition uniform.h:149
UniformClosedOpen(unsigned int i)
Definition uniform.h:152
long double random()
Definition uniform.h:155
long double getUniform()
Definition uniform.h:183
Definition uniform.h:49
double T_numtype
Definition uniform.h:282
double random()
Definition uniform.h:288
double getUniform()
Definition uniform.h:320
UniformClosed(unsigned int i)
Definition uniform.h:285
float random()
Definition uniform.h:242
UniformClosed(unsigned int i)
Definition uniform.h:239
float T_numtype
Definition uniform.h:236
float getUniform()
Definition uniform.h:274
long double random()
Definition uniform.h:335
long double T_numtype
Definition uniform.h:329
long double getUniform()
Definition uniform.h:366
UniformClosed(unsigned int i)
Definition uniform.h:332
Definition uniform.h:230
Definition uniform.h:408
T getUniform()
Definition uniform.h:425
T random()
Definition uniform.h:418
T T_numtype
Definition uniform.h:411
UniformOpenClosed(unsigned int i)
Definition uniform.h:414
UniformOpenClosed()
Definition uniform.h:413
Definition uniform.h:377
T T_numtype
Definition uniform.h:379
UniformOpen()
Definition uniform.h:381
UniformOpen(unsigned int i)
Definition uniform.h:382
T random()
Definition uniform.h:385
T getUniform()
Definition uniform.h:397
Definition uniform.h:192
Uniform()
Definition uniform.h:194
Uniform(unsigned int i)
Definition uniform.h:195
Definition beta.h:50
const long double norm128clos1
Definition uniform.h:222
const long double norm128clos2
Definition uniform.h:223
const long double norm64open
Definition uniform.h:53
const long double norm32open
Definition uniform.h:52
float defaultType
Definition default.h:46
const long double norm96closed2
Definition uniform.h:215
const long double norm128clos4
Definition uniform.h:225
const long double norm96open
Definition uniform.h:54
const long double norm96closed1
Definition uniform.h:214
const long double norm64closed1
Definition uniform.h:208
const long double norm64closed2
Definition uniform.h:210
sharedState defaultState
Definition default.h:55
unsigned int IRNG_int
Definition default.h:57
const long double norm128open
Definition uniform.h:55
MersenneTwister defaultIRNG
Definition default.h:120
const long double norm32closed
Definition uniform.h:204
const long double norm96closed3
Definition uniform.h:217
const long double norm128clos3
Definition uniform.h:224
Definition default.h:53