blitz  Version 0.9
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
rand-uniform.h
Go to the documentation of this file.
1 /***************************************************************************
2  * blitz/rand-uniform.h Uniform class, which provides uniformly
3  * distributed random numbers.
4  *
5  * $Id: rand-uniform.h,v 1.3 2003/01/14 11:29:18 patricg Exp $
6  *
7  * Copyright (C) 1997-2001 Todd Veldhuizen <tveldhui@oonumerics.org>
8  *
9  * This program is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU General Public License
11  * as published by the Free Software Foundation; either version 2
12  * of the License, or (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * Suggestions: blitz-dev@oonumerics.org
20  * Bugs: blitz-bugs@oonumerics.org
21  *
22  * For more information, please see the Blitz++ Home Page:
23  * http://oonumerics.org/blitz/
24  *
25  ***************************************************************************
26  *
27  * This random number generator is based on the LAPACK auxilliary
28  * routine DLARAN by Jack Dongarra. It's a multiplicative congruential
29  * generator with modulus 2^48 and multiplier 33952834046453.
30  *
31  * See also: G. S. Fishman, Multiplicative congruential random number
32  * generators with modulus 2^b: an exhaustive analysis for b=32 and
33  * a partial analysis for b=48, Math. Comp. 189, pp 331-344, 1990.
34  *
35  * This routine requires 32-bit integers.
36  *
37  * The generated number lies in the open interval (low,high). i.e. low and
38  * high themselves will never be generated.
39  *
40  ***************************************************************************/
41 
42 #ifndef BZ_RAND_UNIFORM_H
43 #define BZ_RAND_UNIFORM_H
44 
45 #ifndef BZ_RANDOM_H
46  #include <blitz/random.h>
47 #endif
48 
49 BZ_NAMESPACE(blitz)
50 
51 class Uniform {
52 
53 public:
54  typedef double T_numtype;
55 
56  Uniform(double low = 0.0, double high = 1.0, double = 0.0)
57  : low_(low), length_(high-low)
58  {
59  BZPRECONDITION(sizeof(int) >= 4); // Need 32 bit integers!
60 
61  seed[0] = 24; // All seeds in the range [0,4095]
62  seed[1] = 711;
63  seed[2] = 3;
64  seed[3] = 3721; // The last seed must be odd
65  }
66 
67  void randomize()
68  {
69  BZ_NOT_IMPLEMENTED(); // NEEDS_WORK
70 
71  BZPOSTCONDITION(seed[3] % 2 == 1);
72  }
73 
74  // I'm trying to avoid having a compiled
75  // portion of the library, so this is inline until I
76  // figure out a better way to do this or I change my mind.
77  // -- TV
78  // NEEDS_WORK
79  double random()
80  {
81  BZPRECONDITION(seed[3] % 2 == 1);
82 
83  int it0, it1, it2, it3;
84  it3 = seed[3] * 2549;
85  it2 = it3 / 4096;
86  it3 -= it2 << 12;
87  it2 += seed[2] * 2549 + seed[3] * 2508;
88  it1 = it2 / 4096;
89  it2 -= it1 << 12;
90  it1 += seed[1] * 2549 + seed[2] * 2508 + seed[3] * 322;
91  it0 = it1 / 4096;
92  it1 -= it0 << 12;
93  it0 += seed[0] * 2549 + seed[1] * 2508 + seed[2] * 322 + seed[3] * 494;
94  it0 %= 4096;
95  seed[0] = it0;
96  seed[1] = it1;
97  seed[2] = it2;
98  seed[3] = it3;
99 
100  const double z = 1 / 4096.;
101  return low_ + length_ * (it0 + (it1 + (it2 + it3 * z) * z) * z) * z;
102  }
103 
104  operator double()
105  { return random(); }
106 
107 private:
108  double low_, length_;
109 
110  int seed[4];
111 };
112 
114 
115 #endif // BZ_RAND_UNIFORM_H
116