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>
49 class UniformClosedOpen { };
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 
UniformClosedOpen()66   UniformClosedOpen() {};
UniformClosedOpen(unsigned int i)67   UniformClosedOpen(unsigned int i) :
68     IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
69 
random()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 
getUniform()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 
UniformClosedOpen()108   UniformClosedOpen() {}
UniformClosedOpen(unsigned int i)109   UniformClosedOpen(unsigned int i) :
110     IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
111 
random()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 
getUniform()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 
UniformClosedOpen()151   UniformClosedOpen() {};
UniformClosedOpen(unsigned int i)152   UniformClosedOpen(unsigned int i) :
153     IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
154 
random()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 
getUniform()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:
Uniform()194   Uniform() {};
Uniform(unsigned int i)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 
UniformClosed()238   UniformClosed() {};
UniformClosed(unsigned int i)239   UniformClosed(unsigned int i) :
240     IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
241 
random()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 
getUniform()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 
UniformClosed()284   UniformClosed() {};
UniformClosed(unsigned int i)285   UniformClosed(unsigned int i) :
286     IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
287 
random()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 
getUniform()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 
UniformClosed()331   UniformClosed() {};
UniformClosed(unsigned int i)332   UniformClosed(unsigned int i) :
333     IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
334 
random()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 
getUniform()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 
UniformOpen()381   UniformOpen() {};
UniformOpen(unsigned int i)382   UniformOpen(unsigned int i) :
383     UniformClosedOpen<T,IRNG,stateTag>(i) {};
384 
random()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 
getUniform()397     T getUniform()
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 
UniformOpenClosed()413   UniformOpenClosed() {};
UniformOpenClosed(unsigned int i)414   UniformOpenClosed(unsigned int i) :
415     UniformClosedOpen<T,IRNG,stateTag>(i) {};
416 
417 
random()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 
getUniform()425     T getUniform()
426     { return random(); }
427 };
428 
429 }
430 
431 #endif // BZ_RANDOM_UNIFORM_H
432