1 // -*- C++ -*- 2 /*************************************************************************** 3 * random/uniform.h Uniform IRNG wrapper class 4 * 5 * $Id$ 6 * 7 * Copyright (C) 1997-2011 Todd Veldhuizen <tveldhui@acm.org> 8 * 9 * This file is a part of Blitz. 10 * 11 * Blitz is free software: you can redistribute it and/or modify 12 * it under the terms of the GNU Lesser General Public License 13 * as published by the Free Software Foundation, either version 3 14 * of the License, or (at your option) any later version. 15 * 16 * Blitz is distributed in the hope that it will be useful, 17 * but WITHOUT ANY WARRANTY; without even the implied warranty of 18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 * GNU Lesser General Public License for more details. 20 * 21 * You should have received a copy of the GNU Lesser General Public 22 * License along with Blitz. If not, see <http://www.gnu.org/licenses/>. 23 * 24 * Suggestions: blitz-devel@lists.sourceforge.net 25 * Bugs: blitz-support@lists.sourceforge.net 26 * 27 * For more information, please see the Blitz++ Home Page: 28 * https://sourceforge.net/projects/blitz/ 29 * 30 ***************************************************************************/ 31 32 #ifndef BZ_RANDOM_UNIFORM_H 33 #define BZ_RANDOM_UNIFORM_H 34 35 #include <random/default.h> 36 37 #ifndef FLT_MANT_DIG 38 #include <float.h> 39 #endif 40 41 namespace ranlib { 42 43 /***************************************************************************** 44 * UniformClosedOpen generator: uniform random numbers in [0,1). 45 *****************************************************************************/ 46 47 template<typename T = defaultType, typename IRNG = defaultIRNG, 48 typename stateTag = defaultState> 49 class UniformClosedOpen { }; 50 51 // These constants are 1/2^32, 1/2^64, 1/2^96, 1/2^128 52 const long double norm32open = .2328306436538696289062500000000000000000E-9L; 53 const long double norm64open = .5421010862427522170037264004349708557129E-19L; 54 const long double norm96open = .1262177448353618888658765704452457967477E-28L; 55 const long double norm128open = .2938735877055718769921841343055614194547E-38L; 56 57 58 template<typename IRNG, typename stateTag> 59 class UniformClosedOpen<float,IRNG,stateTag> 60 : public IRNGWrapper<IRNG,stateTag> 61 { 62 63 public: 64 typedef float T_numtype; 65 UniformClosedOpen()66 UniformClosedOpen() {}; UniformClosedOpen(unsigned int i)67 UniformClosedOpen(unsigned int i) : 68 IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {}; 69 random()70 float random() 71 { 72 #if FLT_MANT_DIG > 96 73 #if (FLT_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS) 74 #error RNG code assumes float mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning. 75 #endif 76 IRNG_int i1 = this->irng_.random(); 77 IRNG_int i2 = this->irng_.random(); 78 IRNG_int i3 = this->irng_.random(); 79 IRNG_int i4 = this->irng_.random(); 80 return i1 * norm128open + i2 * norm96open + i3 * norm64open 81 + i4 * norm32open; 82 #elif FLT_MANT_DIG > 64 83 IRNG_int i1 = this->irng_.random(); 84 IRNG_int i2 = this->irng_.random(); 85 IRNG_int i3 = this->irng_.random(); 86 return i1 * norm96open + i2 * norm64open + i3 * norm32open; 87 #elif FLT_MANT_DIG > 32 88 IRNG_int i1 = this->irng_.random(); 89 IRNG_int i2 = this->irng_.random(); 90 return i1 * norm64open + i2 * norm32open; 91 #else 92 IRNG_int i1 = this->irng_.random(); 93 return i1 * norm32open; 94 #endif 95 } 96 getUniform()97 float getUniform() 98 { return random(); } 99 }; 100 101 template<typename IRNG, typename stateTag> 102 class UniformClosedOpen<double,IRNG,stateTag> 103 : public IRNGWrapper<IRNG,stateTag> 104 { 105 public: 106 typedef double T_numtype; 107 UniformClosedOpen()108 UniformClosedOpen() {} UniformClosedOpen(unsigned int i)109 UniformClosedOpen(unsigned int i) : 110 IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {}; 111 random()112 double random() 113 { 114 #if DBL_MANT_DIG > 96 115 #if (DBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS) 116 #error RNG code assumes double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning. 117 #endif 118 119 IRNG_int i1 = this->irng_.random(); 120 IRNG_int i2 = this->irng_.random(); 121 IRNG_int i3 = this->irng_.random(); 122 IRNG_int i4 = this->irng_.random(); 123 return i1 * norm128open + i2 * norm96open + i3 * norm64open 124 + i4 * norm32open; 125 #elif DBL_MANT_DIG > 64 126 IRNG_int i1 = this->irng_.random(); 127 IRNG_int i2 = this->irng_.random(); 128 IRNG_int i3 = this->irng_.random(); 129 return i1 * norm96open + i2 * norm64open + i3 * norm32open; 130 #elif DBL_MANT_DIG > 32 131 IRNG_int i1 = this->irng_.random(); 132 IRNG_int i2 = this->irng_.random(); 133 return i1 * norm64open + i2 * norm32open; 134 #else 135 IRNG_int i1 = this->irng_.random(); 136 return i1 * norm32open; 137 #endif 138 139 } 140 getUniform()141 double getUniform() { return random(); } 142 }; 143 144 template<typename IRNG, typename stateTag> 145 class UniformClosedOpen<long double,IRNG,stateTag> 146 : public IRNGWrapper<IRNG,stateTag> 147 { 148 public: 149 typedef long double T_numtype; 150 UniformClosedOpen()151 UniformClosedOpen() {}; UniformClosedOpen(unsigned int i)152 UniformClosedOpen(unsigned int i) : 153 IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {}; 154 random()155 long double random() 156 { 157 #if LDBL_MANT_DIG > 96 158 #if (LDBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS) 159 #error RNG code assumes long double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning. 160 #endif 161 162 IRNG_int i1 = this->irng_.random(); 163 IRNG_int i2 = this->irng_.random(); 164 IRNG_int i3 = this->irng_.random(); 165 IRNG_int i4 = this->irng_.random(); 166 return i1 * norm128open + i2 * norm96open + i3 * norm64open 167 + i4 * norm32open; 168 #elif LDBL_MANT_DIG > 64 169 IRNG_int i1 = this->irng_.random(); 170 IRNG_int i2 = this->irng_.random(); 171 IRNG_int i3 = this->irng_.random(); 172 return i1 * norm96open + i2 * norm64open + i3 * norm32open; 173 #elif LDBL_MANT_DIG > 32 174 IRNG_int i1 = this->irng_.random(); 175 IRNG_int i2 = this->irng_.random(); 176 return i1 * norm64open + i2 * norm32open; 177 #else 178 IRNG_int i1 = this->irng_.random(); 179 return i1 * norm32open; 180 #endif 181 } 182 getUniform()183 long double getUniform() { return random(); } 184 }; 185 186 // For people who don't care about open or closed: just give them 187 // ClosedOpen (this is the fastest). 188 189 template<class T, typename IRNG = defaultIRNG, 190 typename stateTag = defaultState> 191 class Uniform : public UniformClosedOpen<T,IRNG,stateTag> 192 { 193 public: Uniform()194 Uniform() {}; Uniform(unsigned int i)195 Uniform(unsigned int i) : 196 UniformClosedOpen<T,IRNG,stateTag>::UniformClosedOpen(i) {} 197 }; 198 199 /***************************************************************************** 200 * UniformClosed generator: uniform random numbers in [0,1]. 201 *****************************************************************************/ 202 203 // This constant is 1/(2^32-1) 204 const long double norm32closed = .2328306437080797375431469961868475648078E-9L; 205 206 // These constants are 2^32/(2^64-1) and 1/(2^64-1), respectively. 207 208 const long double norm64closed1 = 209 .23283064365386962891887177448353618888727188481031E-9L; 210 const long double norm64closed2 = 211 .54210108624275221703311375920552804341370213034169E-19L; 212 213 // These constants are 2^64/(2^96-1), 2^32/(2^96-1), and 1/(2^96-1) 214 const long double norm96closed1 = .2328306436538696289062500000029387358771E-9L; 215 const long double norm96closed2 = 216 .5421010862427522170037264004418131333707E-19L; 217 const long double norm96closed3 = 218 .1262177448353618888658765704468388886588E-28L; 219 220 // These constants are 2^96/(2^128-1), 2^64/(2^128-1), 2^32/(2^128-1) and 221 // 1/(2^128-1). 222 const long double norm128clos1 = .2328306436538696289062500000000000000007E-9L; 223 const long double norm128clos2 = .5421010862427522170037264004349708557145E-19L; 224 const long double norm128clos3 = .1262177448353618888658765704452457967481E-28L; 225 const long double norm128clos4 = .2938735877055718769921841343055614194555E-38L; 226 227 228 template<typename T = double, typename IRNG = defaultIRNG, 229 typename stateTag = defaultState> 230 class UniformClosed { }; 231 232 template<typename IRNG, typename stateTag> 233 class UniformClosed<float,IRNG,stateTag> : public IRNGWrapper<IRNG,stateTag> { 234 235 public: 236 typedef float T_numtype; 237 UniformClosed()238 UniformClosed() {}; UniformClosed(unsigned int i)239 UniformClosed(unsigned int i) : 240 IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {}; 241 random()242 float random() 243 { 244 #if FLTMANT_DIG > 96 245 #if (FLT_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS) 246 #error RNG code assumes float mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning. 247 #endif 248 IRNG_int i1 = this->irng_.random(); 249 IRNG_int i2 = this->irng_.random(); 250 IRNG_int i3 = this->irng_.random(); 251 IRNG_int i4 = this->irng_.random(); 252 253 return i1 * norm128clos1 + i2 * norm128clos2 254 + i3 * norm128clos3 + i4 * norm128clos4; 255 #elif FLT_MANT_DIG > 64 256 IRNG_int i1 = this->irng_.random(); 257 IRNG_int i2 = this->irng_.random(); 258 IRNG_int i3 = this->irng_.random(); 259 260 return i1 * norm96closed1 + i2 * norm96closed2 261 + i3 * norm96closed3; 262 #elif FLT_MANT_DIG > 32 263 IRNG_int i1 = this->irng_.random(); 264 IRNG_int i2 = this->irng_.random(); 265 266 return i1 * norm64closed1 + i2 * norm64closed2; 267 #else 268 IRNG_int i = this->irng_.random(); 269 return i * norm32closed; 270 #endif 271 272 } 273 getUniform()274 float getUniform() 275 { return random(); } 276 }; 277 278 template<typename IRNG, typename stateTag> 279 class UniformClosed<double,IRNG,stateTag> : public IRNGWrapper<IRNG,stateTag> { 280 281 public: 282 typedef double T_numtype; 283 UniformClosed()284 UniformClosed() {}; UniformClosed(unsigned int i)285 UniformClosed(unsigned int i) : 286 IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {}; 287 random()288 double random() 289 { 290 #if DBL_MANT_DIG > 96 291 #if (DBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS) 292 #error RNG code assumes double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning. 293 #endif 294 IRNG_int i1 = this->irng_.random(); 295 IRNG_int i2 = this->irng_.random(); 296 IRNG_int i3 = this->irng_.random(); 297 IRNG_int i4 = this->irng_.random(); 298 299 return i1 * norm128clos1 + i2 * norm128clos2 300 + i3 * norm128clos3 + i4 * norm128clos4; 301 #elif DBL_MANT_DIG > 64 302 IRNG_int i1 = this->irng_.random(); 303 IRNG_int i2 = this->irng_.random(); 304 IRNG_int i3 = this->irng_.random(); 305 306 return i1 * norm96closed1 + i2 * norm96closed2 307 + i3 * norm96closed3; 308 #elif DBL_MANT_DIG > 32 309 IRNG_int i1 = this->irng_.random(); 310 IRNG_int i2 = this->irng_.random(); 311 312 return i1 * norm64closed1 + i2 * norm64closed2; 313 #else 314 IRNG_int i = this->irng_.random(); 315 return i * norm32closed; 316 #endif 317 318 } 319 getUniform()320 double getUniform() 321 { return random(); } 322 }; 323 324 template<typename IRNG, typename stateTag> 325 class UniformClosed<long double,IRNG,stateTag> 326 : public IRNGWrapper<IRNG,stateTag> { 327 328 public: 329 typedef long double T_numtype; 330 UniformClosed()331 UniformClosed() {}; UniformClosed(unsigned int i)332 UniformClosed(unsigned int i) : 333 IRNGWrapper<IRNG,stateTag>::IRNGWrapper(i) {}; 334 random()335 long double random() 336 { 337 #if LDBL_MANT_DIG > 96 338 #if (LDBL_MANT_DIG > 128) && !defined(BZ_IGNORE_RNG_ERRORS) 339 #error RNG code assumes long double mantissa is <= 128 bits (not true for your platform). Use -DBZ_IGNORE_RNG_ERRORS to ignore this warning. 340 #endif 341 IRNG_int i1 = this->irng_.random(); 342 IRNG_int i2 = this->irng_.random(); 343 IRNG_int i3 = this->irng_.random(); 344 IRNG_int i4 = this->irng_.random(); 345 346 return i1 * norm128clos1 + i2 * norm128clos2 347 + i3 * norm128clos3 + i4 * norm128clos4; 348 #elif LDBL_MANT_DIG > 64 349 IRNG_int i1 = this->irng_.random(); 350 IRNG_int i2 = this->irng_.random(); 351 IRNG_int i3 = this->irng_.random(); 352 353 return i1 * norm96closed1 + i2 * norm96closed2 354 + i3 * norm96closed3; 355 #elif LDBL_MANT_DIG > 32 356 IRNG_int i1 = this->irng_.random(); 357 IRNG_int i2 = this->irng_.random(); 358 359 return i1 * norm64closed1 + i2 * norm64closed2; 360 #else 361 IRNG_int i = this->irng_.random(); 362 return i * norm32closed; 363 #endif 364 } 365 getUniform()366 long double getUniform() 367 { return random(); } 368 369 }; 370 371 /***************************************************************************** 372 * UniformOpen generator: uniform random numbers in (0,1). 373 *****************************************************************************/ 374 375 template<typename T = double, typename IRNG = defaultIRNG, 376 typename stateTag = defaultState> 377 class UniformOpen : public UniformClosedOpen<T,IRNG,stateTag> { 378 public: 379 typedef T T_numtype; 380 UniformOpen()381 UniformOpen() {}; UniformOpen(unsigned int i)382 UniformOpen(unsigned int i) : 383 UniformClosedOpen<T,IRNG,stateTag>(i) {}; 384 random()385 T random() 386 { 387 // Turn a [0,1) into a (0,1) interval by weeding out 388 // any zeros. 389 T x; 390 do { 391 x = UniformClosedOpen<T,IRNG,stateTag>::random(); 392 } while (x == 0.0L); 393 394 return x; 395 } 396 getUniform()397 T getUniform() 398 { return random(); } 399 400 }; 401 402 /***************************************************************************** 403 * UniformOpenClosed generator: uniform random numbers in (0,1] 404 *****************************************************************************/ 405 406 template<typename T = double, typename IRNG = defaultIRNG, 407 typename stateTag = defaultState> 408 class UniformOpenClosed : public UniformClosedOpen<T,IRNG,stateTag> { 409 410 public: 411 typedef T T_numtype; 412 UniformOpenClosed()413 UniformOpenClosed() {}; UniformOpenClosed(unsigned int i)414 UniformOpenClosed(unsigned int i) : 415 UniformClosedOpen<T,IRNG,stateTag>(i) {}; 416 417 random()418 T random() 419 { 420 // Antithetic value: taking 1-X where X is [0,1) results 421 // in a (0,1] distribution. 422 return 1.0 - UniformClosedOpen<T,IRNG,stateTag>::random(); 423 } 424 getUniform()425 T getUniform() 426 { return random(); } 427 }; 428 429 } 430 431 #endif // BZ_RANDOM_UNIFORM_H 432