blitz Version 0.10
random/uniform.h
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 /***************************************************************************
00003  * random/uniform.h       Uniform IRNG wrapper class
00004  *
00005  * $Id: uniform.h,v 1.5 2011/03/25 22:41:17 julianc Exp $
00006  *
00007  * Copyright (C) 1997-2011 Todd Veldhuizen <tveldhui@acm.org>
00008  *
00009  * This file is a part of Blitz.
00010  *
00011  * Blitz is free software: you can redistribute it and/or modify 
00012  * it under the terms of the GNU Lesser General Public License
00013  * as published by the Free Software Foundation, either version 3
00014  * of the License, or (at your option) any later version.
00015  *
00016  * Blitz is distributed in the hope that it will be useful,
00017  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00018  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00019  * GNU Lesser General Public License for more details.
00020  *
00021  * You should have received a copy of the GNU Lesser General Public 
00022  * License along with Blitz.  If not, see <http://www.gnu.org/licenses/>.
00023  * 
00024  * Suggestions:          blitz-devel@lists.sourceforge.net
00025  * Bugs:                 blitz-support@lists.sourceforge.net    
00026  *
00027  * For more information, please see the Blitz++ Home Page:
00028  *    https://sourceforge.net/projects/blitz/
00029  *
00030  ***************************************************************************/
00031 
00032 #ifndef BZ_RANDOM_UNIFORM_H
00033 #define BZ_RANDOM_UNIFORM_H
00034 
00035 #include <random/default.h>
00036 
00037 #ifndef FLT_MANT_DIG
00038  #include <float.h>
00039 #endif
00040 
00041 BZ_NAMESPACE(ranlib)
00042 
00043 /*****************************************************************************
00044  * UniformClosedOpen generator: uniform random numbers in [0,1).
00045  *****************************************************************************/
00046 
00047 template<typename T = defaultType, typename IRNG = defaultIRNG, 
00048     typename stateTag = defaultState>
00049 class UniformClosedOpen { };
00050 
00051 // These constants are 1/2^32, 1/2^64, 1/2^96, 1/2^128
00052 const long double norm32open = .2328306436538696289062500000000000000000E-9L;
00053 const long double norm64open = .5421010862427522170037264004349708557129E-19L;
00054 const long double norm96open = .1262177448353618888658765704452457967477E-28L;
00055 const long double norm128open = .2938735877055718769921841343055614194547E-38L;
00056 
00057 
00058 template<typename IRNG, typename stateTag>
00059 class UniformClosedOpen<float,IRNG,stateTag> 
00060   : public IRNGWrapper<IRNG,stateTag> 
00061 {
00062 
00063 public:
00064     typedef float T_numtype;
00065 
00066   UniformClosedOpen() {};
00067   UniformClosedOpen(unsigned int i) : 
00068     IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
00069 
00070     float random()
00071     {
00072 #if FLT_MANT_DIG > 96
00073  #if (FLT_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
00074   #error RNG code assumes float mantissa is <= 128 bits (not true for your platform).  Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
00075  #endif
00076         IRNG_int i1 = this->irng_.random();
00077         IRNG_int i2 = this->irng_.random();
00078         IRNG_int i3 = this->irng_.random();
00079         IRNG_int i4 = this->irng_.random();
00080         return i1 * norm128open + i2 * norm96open + i3 * norm64open
00081             + i4 * norm32open;
00082 #elif FLT_MANT_DIG > 64
00083         IRNG_int i1 = this->irng_.random();
00084         IRNG_int i2 = this->irng_.random();
00085         IRNG_int i3 = this->irng_.random();
00086         return i1 * norm96open + i2 * norm64open + i3 * norm32open;
00087 #elif FLT_MANT_DIG > 32
00088         IRNG_int i1 = this->irng_.random();
00089         IRNG_int i2 = this->irng_.random();
00090         return i1 * norm64open + i2 * norm32open;
00091 #else
00092         IRNG_int i1 = this->irng_.random();
00093         return i1 * norm32open;
00094 #endif
00095     }    
00096 
00097     float getUniform()
00098     { return random(); }
00099 };
00100 
00101 template<typename IRNG, typename stateTag>
00102 class UniformClosedOpen<double,IRNG,stateTag> 
00103   : public IRNGWrapper<IRNG,stateTag> 
00104 {
00105 public:
00106     typedef double T_numtype;
00107 
00108   UniformClosedOpen() {}
00109   UniformClosedOpen(unsigned int i) : 
00110     IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
00111 
00112     double random()
00113     {    
00114 #if DBL_MANT_DIG > 96
00115  #if (DBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
00116   #error RNG code assumes double mantissa is <= 128 bits (not true for your platform).  Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
00117  #endif
00118 
00119         IRNG_int i1 = this->irng_.random();
00120         IRNG_int i2 = this->irng_.random();
00121         IRNG_int i3 = this->irng_.random();
00122         IRNG_int i4 = this->irng_.random();
00123         return i1 * norm128open + i2 * norm96open + i3 * norm64open
00124             + i4 * norm32open;
00125 #elif DBL_MANT_DIG > 64
00126         IRNG_int i1 = this->irng_.random();
00127         IRNG_int i2 = this->irng_.random();
00128         IRNG_int i3 = this->irng_.random();
00129         return i1 * norm96open + i2 * norm64open + i3 * norm32open;
00130 #elif DBL_MANT_DIG > 32
00131         IRNG_int i1 = this->irng_.random();
00132         IRNG_int i2 = this->irng_.random();
00133         return i1 * norm64open + i2 * norm32open;
00134 #else
00135         IRNG_int i1 = this->irng_.random();
00136         return i1 * norm32open;
00137 #endif
00138 
00139     }
00140 
00141     double getUniform() { return random(); }
00142 };
00143 
00144 template<typename IRNG, typename stateTag>
00145 class UniformClosedOpen<long double,IRNG,stateTag>
00146   : public IRNGWrapper<IRNG,stateTag> 
00147 {
00148 public:
00149     typedef long double T_numtype;
00150 
00151   UniformClosedOpen() {};
00152   UniformClosedOpen(unsigned int i) : 
00153     IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
00154 
00155     long double random()
00156     {
00157 #if LDBL_MANT_DIG > 96
00158  #if (LDBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
00159   #error RNG code assumes long double mantissa is <= 128 bits (not true for your platform).  Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.  
00160 #endif
00161 
00162         IRNG_int i1 = this->irng_.random();
00163         IRNG_int i2 = this->irng_.random();
00164         IRNG_int i3 = this->irng_.random();
00165         IRNG_int i4 = this->irng_.random();
00166         return i1 * norm128open + i2 * norm96open + i3 * norm64open
00167             + i4 * norm32open;
00168 #elif LDBL_MANT_DIG > 64
00169         IRNG_int i1 = this->irng_.random();
00170         IRNG_int i2 = this->irng_.random();
00171         IRNG_int i3 = this->irng_.random();
00172         return i1 * norm96open + i2 * norm64open + i3 * norm32open;
00173 #elif LDBL_MANT_DIG > 32
00174         IRNG_int i1 = this->irng_.random();
00175         IRNG_int i2 = this->irng_.random();
00176         return i1 * norm64open + i2 * norm32open;
00177 #else
00178         IRNG_int i1 = this->irng_.random();
00179         return i1 * norm32open;
00180 #endif
00181     }
00182 
00183     long double getUniform() { return random(); }
00184 };
00185 
00186 // For people who don't care about open or closed: just give them
00187 // ClosedOpen (this is the fastest).
00188 
00189 template<class T, typename IRNG = defaultIRNG, 
00190     typename stateTag = defaultState>
00191 class Uniform : public UniformClosedOpen<T,IRNG,stateTag> 
00192 {
00193 public:
00194   Uniform() {};
00195   Uniform(unsigned int i) : 
00196     UniformClosedOpen<T,IRNG,stateTag>::UniformClosedOpen(i) {}
00197  };
00198 
00199 /*****************************************************************************
00200  * UniformClosed generator: uniform random numbers in [0,1].
00201  *****************************************************************************/
00202 
00203 // This constant is 1/(2^32-1)
00204 const long double norm32closed = .2328306437080797375431469961868475648078E-9L;
00205 
00206 // These constants are 2^32/(2^64-1) and 1/(2^64-1), respectively.
00207 
00208 const long double norm64closed1 =
00209     .23283064365386962891887177448353618888727188481031E-9L;
00210 const long double norm64closed2 =
00211     .54210108624275221703311375920552804341370213034169E-19L;
00212 
00213 // These constants are 2^64/(2^96-1), 2^32/(2^96-1), and 1/(2^96-1)
00214 const long double norm96closed1 = .2328306436538696289062500000029387358771E-9L;
00215 const long double norm96closed2 =
00216     .5421010862427522170037264004418131333707E-19L;
00217 const long double norm96closed3 =
00218     .1262177448353618888658765704468388886588E-28L;
00219 
00220 // These constants are 2^96/(2^128-1), 2^64/(2^128-1), 2^32/(2^128-1) and
00221 // 1/(2^128-1).
00222 const long double norm128clos1 = .2328306436538696289062500000000000000007E-9L;
00223 const long double norm128clos2 = .5421010862427522170037264004349708557145E-19L;
00224 const long double norm128clos3 = .1262177448353618888658765704452457967481E-28L;
00225 const long double norm128clos4 = .2938735877055718769921841343055614194555E-38L;
00226 
00227 
00228 template<typename T = double, typename IRNG = defaultIRNG, 
00229     typename stateTag = defaultState>
00230 class UniformClosed { };
00231 
00232 template<typename IRNG, typename stateTag>
00233 class UniformClosed<float,IRNG,stateTag> : public IRNGWrapper<IRNG,stateTag> {
00234 
00235 public:
00236     typedef float T_numtype;
00237 
00238   UniformClosed() {};
00239   UniformClosed(unsigned int i) : 
00240     IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
00241 
00242     float random()
00243     {
00244 #if FLTMANT_DIG > 96
00245  #if (FLT_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
00246   #error RNG code assumes float mantissa is <= 128 bits (not true for your platform).  Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
00247  #endif
00248         IRNG_int i1 = this->irng_.random();
00249         IRNG_int i2 = this->irng_.random();
00250         IRNG_int i3 = this->irng_.random();
00251         IRNG_int i4 = this->irng_.random();
00252 
00253         return i1 * norm128clos1 + i2 * norm128clos2
00254           + i3 * norm128clos3 + i4 * norm128clos4;
00255 #elif FLT_MANT_DIG > 64
00256         IRNG_int i1 = this->irng_.random();
00257         IRNG_int i2 = this->irng_.random();
00258         IRNG_int i3 = this->irng_.random();
00259 
00260         return i1 * norm96closed1 + i2 * norm96closed2
00261           + i3 * norm96closed3;
00262 #elif FLT_MANT_DIG > 32
00263         IRNG_int i1 = this->irng_.random();
00264         IRNG_int i2 = this->irng_.random();
00265 
00266         return i1 * norm64closed1 + i2 * norm64closed2;
00267 #else
00268         IRNG_int i = this->irng_.random();
00269         return i * norm32closed;
00270 #endif
00271 
00272     }
00273 
00274     float getUniform()
00275     { return random(); }
00276 };
00277 
00278 template<typename IRNG, typename stateTag>
00279 class UniformClosed<double,IRNG,stateTag> : public IRNGWrapper<IRNG,stateTag> {
00280 
00281 public:
00282     typedef double T_numtype;
00283 
00284   UniformClosed() {};
00285   UniformClosed(unsigned int i) : 
00286     IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
00287 
00288     double random()
00289     {
00290 #if DBL_MANT_DIG > 96
00291  #if (DBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
00292   #error RNG code assumes double mantissa is <= 128 bits (not true for your platform).  Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
00293  #endif
00294         IRNG_int i1 = this->irng_.random();
00295         IRNG_int i2 = this->irng_.random();
00296         IRNG_int i3 = this->irng_.random();
00297         IRNG_int i4 = this->irng_.random();
00298 
00299         return i1 * norm128clos1 + i2 * norm128clos2
00300           + i3 * norm128clos3 + i4 * norm128clos4;
00301 #elif DBL_MANT_DIG > 64
00302         IRNG_int i1 = this->irng_.random();
00303         IRNG_int i2 = this->irng_.random();
00304         IRNG_int i3 = this->irng_.random();
00305 
00306         return i1 * norm96closed1 + i2 * norm96closed2
00307           + i3 * norm96closed3;
00308 #elif DBL_MANT_DIG > 32
00309         IRNG_int i1 = this->irng_.random();
00310         IRNG_int i2 = this->irng_.random();
00311 
00312         return i1 * norm64closed1 + i2 * norm64closed2;
00313 #else
00314         IRNG_int i = this->irng_.random();
00315         return i * norm32closed;
00316 #endif
00317 
00318     }
00319 
00320     double getUniform()
00321     { return random(); }
00322 };
00323 
00324 template<typename IRNG, typename stateTag>
00325 class UniformClosed<long double,IRNG,stateTag>
00326   : public IRNGWrapper<IRNG,stateTag> {
00327 
00328 public:
00329     typedef long double T_numtype;
00330 
00331   UniformClosed() {};
00332   UniformClosed(unsigned int i) : 
00333     IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {};
00334 
00335     long double random()
00336     {
00337 #if LDBL_MANT_DIG > 96
00338  #if (LDBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS)
00339   #error RNG code assumes long double mantissa is <= 128 bits (not true for your platform).  Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning.
00340  #endif
00341         IRNG_int i1 = this->irng_.random();
00342         IRNG_int i2 = this->irng_.random();
00343         IRNG_int i3 = this->irng_.random();
00344         IRNG_int i4 = this->irng_.random();
00345 
00346         return i1 * norm128clos1 + i2 * norm128clos2 
00347           + i3 * norm128clos3 + i4 * norm128clos4;
00348 #elif LDBL_MANT_DIG > 64
00349         IRNG_int i1 = this->irng_.random();
00350         IRNG_int i2 = this->irng_.random();
00351         IRNG_int i3 = this->irng_.random();
00352 
00353         return i1 * norm96closed1 + i2 * norm96closed2
00354           + i3 * norm96closed3;
00355 #elif LDBL_MANT_DIG > 32
00356         IRNG_int i1 = this->irng_.random();
00357         IRNG_int i2 = this->irng_.random();
00358 
00359         return i1 * norm64closed1 + i2 * norm64closed2;
00360 #else
00361         IRNG_int i = this->irng_.random();
00362         return i * norm32closed;
00363 #endif
00364     }
00365 
00366     long double getUniform()
00367     { return random(); }
00368   
00369 };
00370 
00371 /*****************************************************************************
00372  * UniformOpen generator: uniform random numbers in (0,1).
00373  *****************************************************************************/
00374 
00375 template<typename T = double, typename IRNG = defaultIRNG,
00376     typename stateTag = defaultState>
00377 class UniformOpen : public UniformClosedOpen<T,IRNG,stateTag> {
00378   public:
00379     typedef T T_numtype;
00380 
00381   UniformOpen() {};
00382   UniformOpen(unsigned int i) : 
00383     UniformClosedOpen<T,IRNG,stateTag>(i) {};
00384 
00385     T random()
00386     {
00387         // Turn a [0,1) into a (0,1) interval by weeding out
00388         // any zeros.
00389         T x;
00390         do {
00391           x = UniformClosedOpen<T,IRNG,stateTag>::random();
00392         } while (x == 0.0L);
00393 
00394         return x;
00395     }
00396 
00397     T getUniform()
00398     { return random(); }
00399 
00400 };
00401 
00402 /*****************************************************************************
00403  * UniformOpenClosed generator: uniform random numbers in (0,1]
00404  *****************************************************************************/
00405 
00406 template<typename T = double, typename IRNG = defaultIRNG,
00407     typename stateTag = defaultState>
00408 class UniformOpenClosed : public UniformClosedOpen<T,IRNG,stateTag> {
00409 
00410 public:
00411     typedef T T_numtype;
00412 
00413   UniformOpenClosed() {};
00414   UniformOpenClosed(unsigned int i) : 
00415     UniformClosedOpen<T,IRNG,stateTag>(i) {};
00416 
00417 
00418     T random()
00419     {
00420         // Antithetic value: taking 1-X where X is [0,1) results
00421         // in a (0,1] distribution.
00422         return 1.0 - UniformClosedOpen<T,IRNG,stateTag>::random();
00423     }
00424 
00425     T getUniform()
00426     { return random(); }
00427 };
00428 
00429 BZ_NAMESPACE_END
00430 
00431 #endif // BZ_RANDOM_UNIFORM_H
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines