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