blitz/rand-uniform.h

Go to the documentation of this file.
00001 /***************************************************************************
00002  * blitz/rand-uniform.h    Uniform class, which provides uniformly
00003  *                         distributed random numbers.
00004  *
00005  * $Id: rand-uniform.h,v 1.3 2003/01/14 11:29:18 patricg Exp $
00006  *
00007  * Copyright (C) 1997-2001 Todd Veldhuizen <tveldhui@oonumerics.org>
00008  *
00009  * This program is free software; you can redistribute it and/or
00010  * modify it under the terms of the GNU General Public License
00011  * as published by the Free Software Foundation; either version 2
00012  * of the License, or (at your option) any later version.
00013  *
00014  * This program is distributed in the hope that it will be useful,
00015  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00016  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00017  * GNU General Public License for more details.
00018  *
00019  * Suggestions:          blitz-dev@oonumerics.org
00020  * Bugs:                 blitz-bugs@oonumerics.org
00021  *
00022  * For more information, please see the Blitz++ Home Page:
00023  *    http://oonumerics.org/blitz/
00024  *
00025  ***************************************************************************
00026  *
00027  * This random number generator is based on the LAPACK auxilliary
00028  * routine DLARAN by Jack Dongarra.  It's a multiplicative congruential 
00029  * generator with modulus 2^48 and multiplier 33952834046453. 
00030  *
00031  * See also: G. S. Fishman, Multiplicative congruential random number
00032  * generators with modulus 2^b: an exhaustive analysis for b=32 and
00033  * a partial analysis for b=48, Math. Comp. 189, pp 331-344, 1990.
00034  * 
00035  * This routine requires 32-bit integers.
00036  *
00037  * The generated number lies in the open interval (low,high).  i.e. low and
00038  * high themselves will never be generated.
00039  *
00040  ***************************************************************************/
00041 
00042 #ifndef BZ_RAND_UNIFORM_H
00043 #define BZ_RAND_UNIFORM_H
00044 
00045 #ifndef BZ_RANDOM_H
00046  #include <blitz/random.h>
00047 #endif
00048 
00049 BZ_NAMESPACE(blitz)
00050 
00051 class Uniform {
00052 
00053 public:
00054     typedef double T_numtype;
00055 
00056     Uniform(double low = 0.0, double high = 1.0, double = 0.0)
00057         : low_(low), length_(high-low)
00058     { 
00059         BZPRECONDITION(sizeof(int) >= 4);   // Need 32 bit integers!
00060 
00061         seed[0] = 24;       // All seeds in the range [0,4095]
00062         seed[1] = 711;
00063         seed[2] = 3;
00064         seed[3] = 3721;     // The last seed must be odd
00065     }
00066 
00067     void randomize() 
00068     { 
00069         BZ_NOT_IMPLEMENTED();            // NEEDS_WORK
00070 
00071         BZPOSTCONDITION(seed[3] % 2 == 1);
00072     }
00073   
00074     // I'm trying to avoid having a compiled 
00075     // portion of the library, so this is inline until I
00076     // figure out a better way to do this or I change my mind.
00077     // -- TV
00078     // NEEDS_WORK
00079     double random()
00080     { 
00081         BZPRECONDITION(seed[3] % 2 == 1);
00082 
00083         int it0, it1, it2, it3;
00084         it3 = seed[3] * 2549;
00085         it2 = it3 / 4096;
00086         it3 -= it2 << 12;
00087         it2 += seed[2] * 2549 + seed[3] * 2508;
00088         it1 = it2 / 4096;
00089         it2 -= it1 << 12;
00090         it1 += seed[1] * 2549 + seed[2] * 2508 + seed[3] * 322;
00091         it0 = it1 / 4096;
00092         it1 -= it0 << 12;
00093         it0 += seed[0] * 2549 + seed[1] * 2508 + seed[2] * 322 + seed[3] * 494;
00094         it0 %= 4096;
00095         seed[0] = it0;
00096         seed[1] = it1;
00097         seed[2] = it2;
00098         seed[3] = it3;
00099       
00100         const double z = 1 / 4096.;
00101         return low_ + length_ * (it0 + (it1 + (it2 + it3 * z) * z) * z) * z;
00102     } 
00103 
00104     operator double() 
00105     { return random(); }
00106 
00107 private:
00108     double low_, length_;
00109 
00110     int seed[4];
00111 };
00112 
00113 BZ_NAMESPACE_END
00114 
00115 #endif // BZ_RAND_UNIFORM_H
00116 

Generated on Mon Oct 22 09:27:56 2007 for blitz by  doxygen 1.4.7