source: XIOS/dev/dev_olga/src/extern/blitz/include/random/uniform.h @ 1022

Last change on this file since 1022 was 1022, checked in by mhnguyen, 7 years ago
File size: 13.5 KB
Line 
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
41BZ_NAMESPACE(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>
49class UniformClosedOpen { };
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
66  UniformClosedOpen() {};
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
108  UniformClosedOpen() {}
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
151  UniformClosedOpen() {};
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:
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)
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>
230class UniformClosed { };
231
232template<typename IRNG, typename stateTag>
233class UniformClosed<float,IRNG,stateTag> : public IRNGWrapper<IRNG,stateTag> {
234
235public:
236    typedef float T_numtype;
237
238  UniformClosed() {};
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
278template<typename IRNG, typename stateTag>
279class UniformClosed<double,IRNG,stateTag> : public IRNGWrapper<IRNG,stateTag> {
280
281public:
282    typedef double T_numtype;
283
284  UniformClosed() {};
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
331  UniformClosed() {};
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
381  UniformOpen() {};
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 {
391          x = UniformClosedOpen<T,IRNG,stateTag>::random();
392        } while (x == 0.0L);
393
394        return x;
395    }
396
397    T getUniform()
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
413  UniformOpenClosed() {};
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.
422        return 1.0 - UniformClosedOpen<T,IRNG,stateTag>::random();
423    }
424
425    T getUniform()
426    { return random(); }
427};
428
429BZ_NAMESPACE_END
430
431#endif // BZ_RANDOM_UNIFORM_H
Note: See TracBrowser for help on using the repository browser.