blitz Version 0.9
|
00001 #ifndef BZ_RANDOM_UNIFORM_H 00002 #define BZ_RANDOM_UNIFORM_H 00003 00004 #include <random/default.h> 00005 00006 #ifndef FLT_MANT_DIG 00007 #include <float.h> 00008 #endif 00009 00010 BZ_NAMESPACE(ranlib) 00011 00012 /***************************************************************************** 00013 * UniformClosedOpen generator: uniform random numbers in [0,1). 00014 *****************************************************************************/ 00015 00016 template<typename T = defaultType, typename IRNG = defaultIRNG, 00017 typename stateTag = defaultState> 00018 class UniformClosedOpen { }; 00019 00020 // These constants are 1/2^32, 1/2^64, 1/2^96, 1/2^128 00021 const long double norm32open = .2328306436538696289062500000000000000000E-9L; 00022 const long double norm64open = .5421010862427522170037264004349708557129E-19L; 00023 const long double norm96open = .1262177448353618888658765704452457967477E-28L; 00024 const long double norm128open = .2938735877055718769921841343055614194547E-38L; 00025 00026 00027 template<typename IRNG, typename stateTag> 00028 class UniformClosedOpen<float,IRNG,stateTag> 00029 : public IRNGWrapper<IRNG,stateTag> 00030 { 00031 00032 public: 00033 typedef float T_numtype; 00034 00035 float random() 00036 { 00037 #if FLT_MANT_DIG > 96 00038 #if (FLT_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS) 00039 #error RNG code assumes float mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning. 00040 #endif 00041 IRNG_int i1 = this->irng_.random(); 00042 IRNG_int i2 = this->irng_.random(); 00043 IRNG_int i3 = this->irng_.random(); 00044 IRNG_int i4 = this->irng_.random(); 00045 return i1 * norm128open + i2 * norm96open + i3 * norm64open 00046 + i4 * norm32open; 00047 #elif FLT_MANT_DIG > 64 00048 IRNG_int i1 = this->irng_.random(); 00049 IRNG_int i2 = this->irng_.random(); 00050 IRNG_int i3 = this->irng_.random(); 00051 return i1 * norm96open + i2 * norm64open + i3 * norm32open; 00052 #elif FLT_MANT_DIG > 32 00053 IRNG_int i1 = this->irng_.random(); 00054 IRNG_int i2 = this->irng_.random(); 00055 return i1 * norm64open + i2 * norm32open; 00056 #else 00057 IRNG_int i1 = this->irng_.random(); 00058 return i1 * norm32open; 00059 #endif 00060 } 00061 00062 float getUniform() 00063 { return random(); } 00064 }; 00065 00066 template<typename IRNG, typename stateTag> 00067 class UniformClosedOpen<double,IRNG,stateTag> 00068 : public IRNGWrapper<IRNG,stateTag> 00069 { 00070 public: 00071 typedef double T_numtype; 00072 00073 double random() 00074 { 00075 #if DBL_MANT_DIG > 96 00076 #if (DBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS) 00077 #error RNG code assumes double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning. 00078 #endif 00079 00080 IRNG_int i1 = this->irng_.random(); 00081 IRNG_int i2 = this->irng_.random(); 00082 IRNG_int i3 = this->irng_.random(); 00083 IRNG_int i4 = this->irng_.random(); 00084 return i1 * norm128open + i2 * norm96open + i3 * norm64open 00085 + i4 * norm32open; 00086 #elif DBL_MANT_DIG > 64 00087 IRNG_int i1 = this->irng_.random(); 00088 IRNG_int i2 = this->irng_.random(); 00089 IRNG_int i3 = this->irng_.random(); 00090 return i1 * norm96open + i2 * norm64open + i3 * norm32open; 00091 #elif DBL_MANT_DIG > 32 00092 IRNG_int i1 = this->irng_.random(); 00093 IRNG_int i2 = this->irng_.random(); 00094 return i1 * norm64open + i2 * norm32open; 00095 #else 00096 IRNG_int i1 = this->irng_.random(); 00097 return i1 * norm32open; 00098 #endif 00099 00100 } 00101 00102 double getUniform() { return random(); } 00103 }; 00104 00105 template<typename IRNG, typename stateTag> 00106 class UniformClosedOpen<long double,IRNG,stateTag> 00107 : public IRNGWrapper<IRNG,stateTag> 00108 { 00109 public: 00110 typedef long double T_numtype; 00111 00112 long double random() 00113 { 00114 #if LDBL_MANT_DIG > 96 00115 #if (LDBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS) 00116 #error RNG code assumes long 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 LDBL_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 LDBL_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 long double getUniform() { return random(); } 00141 }; 00142 00143 // For people who don't care about open or closed: just give them 00144 // ClosedOpen (this is the fastest). 00145 00146 template<class T, typename IRNG = defaultIRNG, 00147 typename stateTag = defaultState> 00148 class Uniform : public UniformClosedOpen<T,IRNG,stateTag> 00149 { }; 00150 00151 /***************************************************************************** 00152 * UniformClosed generator: uniform random numbers in [0,1]. 00153 *****************************************************************************/ 00154 00155 // This constant is 1/(2^32-1) 00156 const long double norm32closed = .2328306437080797375431469961868475648078E-9L; 00157 00158 // These constants are 2^32/(2^64-1) and 1/(2^64-1), respectively. 00159 00160 const long double norm64closed1 = 00161 .23283064365386962891887177448353618888727188481031E-9L; 00162 const long double norm64closed2 = 00163 .54210108624275221703311375920552804341370213034169E-19L; 00164 00165 // These constants are 2^64/(2^96-1), 2^32/(2^96-1), and 1/(2^96-1) 00166 const long double norm96closed1 = .2328306436538696289062500000029387358771E-9L; 00167 const long double norm96closed2 = 00168 .5421010862427522170037264004418131333707E-19L; 00169 const long double norm96closed3 = 00170 .1262177448353618888658765704468388886588E-28L; 00171 00172 // These constants are 2^96/(2^128-1), 2^64/(2^128-1), 2^32/(2^128-1) and 00173 // 1/(2^128-1). 00174 const long double norm128clos1 = .2328306436538696289062500000000000000007E-9L; 00175 const long double norm128clos2 = .5421010862427522170037264004349708557145E-19L; 00176 const long double norm128clos3 = .1262177448353618888658765704452457967481E-28L; 00177 const long double norm128clos4 = .2938735877055718769921841343055614194555E-38L; 00178 00179 00180 template<typename T = double, typename IRNG = defaultIRNG, 00181 typename stateTag = defaultState> 00182 class UniformClosed { }; 00183 00184 template<typename IRNG, typename stateTag> 00185 class UniformClosed<float,IRNG,stateTag> : public IRNGWrapper<IRNG,stateTag> { 00186 00187 public: 00188 typedef float T_numtype; 00189 00190 float random() 00191 { 00192 #if FLTMANT_DIG > 96 00193 #if (FLT_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS) 00194 #error RNG code assumes float mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning. 00195 #endif 00196 IRNG_int i1 = this->irng_.random(); 00197 IRNG_int i2 = this->irng_.random(); 00198 IRNG_int i3 = this->irng_.random(); 00199 IRNG_int i4 = this->irng_.random(); 00200 00201 return i1 * norm128clos1 + i2 * norm128clos2 00202 + i3 * norm128clos3 + i4 * norm128clos4; 00203 #elif FLT_MANT_DIG > 64 00204 IRNG_int i1 = this->irng_.random(); 00205 IRNG_int i2 = this->irng_.random(); 00206 IRNG_int i3 = this->irng_.random(); 00207 00208 return i1 * norm96closed1 + i2 * norm96closed2 00209 + i3 * norm96closed3; 00210 #elif FLT_MANT_DIG > 32 00211 IRNG_int i1 = this->irng_.random(); 00212 IRNG_int i2 = this->irng_.random(); 00213 00214 return i1 * norm64closed1 + i2 * norm64closed2; 00215 #else 00216 IRNG_int i = this->irng_.random(); 00217 return i * norm32closed; 00218 #endif 00219 00220 } 00221 00222 float getUniform() 00223 { return random(); } 00224 }; 00225 00226 template<typename IRNG, typename stateTag> 00227 class UniformClosed<double,IRNG,stateTag> : public IRNGWrapper<IRNG,stateTag> { 00228 00229 public: 00230 typedef double T_numtype; 00231 00232 double random() 00233 { 00234 #if DBL_MANT_DIG > 96 00235 #if (DBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS) 00236 #error RNG code assumes double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning. 00237 #endif 00238 IRNG_int i1 = this->irng_.random(); 00239 IRNG_int i2 = this->irng_.random(); 00240 IRNG_int i3 = this->irng_.random(); 00241 IRNG_int i4 = this->irng_.random(); 00242 00243 return i1 * norm128clos1 + i2 * norm128clos2 00244 + i3 * norm128clos3 + i4 * norm128clos4; 00245 #elif DBL_MANT_DIG > 64 00246 IRNG_int i1 = this->irng_.random(); 00247 IRNG_int i2 = this->irng_.random(); 00248 IRNG_int i3 = this->irng_.random(); 00249 00250 return i1 * norm96closed1 + i2 * norm96closed2 00251 + i3 * norm96closed3; 00252 #elif DBL_MANT_DIG > 32 00253 IRNG_int i1 = this->irng_.random(); 00254 IRNG_int i2 = this->irng_.random(); 00255 00256 return i1 * norm64closed1 + i2 * norm64closed2; 00257 #else 00258 IRNG_int i = this->irng_.random(); 00259 return i * norm32closed; 00260 #endif 00261 00262 } 00263 00264 double getUniform() 00265 { return random(); } 00266 }; 00267 00268 template<typename IRNG, typename stateTag> 00269 class UniformClosed<long double,IRNG,stateTag> 00270 : public IRNGWrapper<IRNG,stateTag> { 00271 00272 public: 00273 typedef long double T_numtype; 00274 00275 long double random() 00276 { 00277 #if LDBL_MANT_DIG > 96 00278 #if (LDBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS) 00279 #error RNG code assumes long double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning. 00280 #endif 00281 IRNG_int i1 = this->irng_.random(); 00282 IRNG_int i2 = this->irng_.random(); 00283 IRNG_int i3 = this->irng_.random(); 00284 IRNG_int i4 = this->irng_.random(); 00285 00286 return i1 * norm128clos1 + i2 * norm128clos2 00287 + i3 * norm128clos3 + i4 * norm128clos4; 00288 #elif LDBL_MANT_DIG > 64 00289 IRNG_int i1 = this->irng_.random(); 00290 IRNG_int i2 = this->irng_.random(); 00291 IRNG_int i3 = this->irng_.random(); 00292 00293 return i1 * norm96closed1 + i2 * norm96closed2 00294 + i3 * norm96closed3; 00295 #elif LDBL_MANT_DIG > 32 00296 IRNG_int i1 = this->irng_.random(); 00297 IRNG_int i2 = this->irng_.random(); 00298 00299 return i1 * norm64closed1 + i2 * norm64closed2; 00300 #else 00301 IRNG_int i = this->irng_.random(); 00302 return i * norm32closed; 00303 #endif 00304 } 00305 00306 long double getUniform() 00307 { return random(); } 00308 00309 }; 00310 00311 /***************************************************************************** 00312 * UniformOpen generator: uniform random numbers in (0,1). 00313 *****************************************************************************/ 00314 00315 template<typename T = double, typename IRNG = defaultIRNG, 00316 typename stateTag = defaultState> 00317 class UniformOpen : public UniformClosedOpen<T,IRNG,stateTag> { 00318 public: 00319 typedef T T_numtype; 00320 00321 T random() 00322 { 00323 // Turn a [0,1) into a (0,1) interval by weeding out 00324 // any zeros. 00325 T x; 00326 do { 00327 x = UniformClosedOpen<T,IRNG,stateTag>::random(); 00328 } while (x == 0.0L); 00329 00330 return x; 00331 } 00332 00333 T getUniform() 00334 { return random(); } 00335 00336 }; 00337 00338 /***************************************************************************** 00339 * UniformOpenClosed generator: uniform random numbers in (0,1] 00340 *****************************************************************************/ 00341 00342 template<typename T = double, typename IRNG = defaultIRNG, 00343 typename stateTag = defaultState> 00344 class UniformOpenClosed : public UniformClosedOpen<T,IRNG,stateTag> { 00345 00346 public: 00347 typedef T T_numtype; 00348 00349 T random() 00350 { 00351 // Antithetic value: taking 1-X where X is [0,1) results 00352 // in a (0,1] distribution. 00353 return 1.0 - UniformClosedOpen<T,IRNG,stateTag>::random(); 00354 } 00355 00356 T getUniform() 00357 { return random(); } 00358 }; 00359 00360 BZ_NAMESPACE_END 00361 00362 #endif // BZ_RANDOM_UNIFORM_H