Random.h

Go to the documentation of this file.
00001 
00007 /*
00008  * Copyright (c) 1993 Martin Birgmeier
00009  * All rights reserved.
00010  *
00011  * You may redistribute unmodified or modified versions of this source
00012  * code provided that the above copyright notice and this and the
00013  * following conditions are retained.
00014  *
00015  * This software is provided ``as is'', and comes with no warranties
00016  * of any kind. I shall in no event be liable for anything that happens
00017  * to anyone/anything when using this software.
00018  */
00019 
00020 #ifndef RANDOM_H
00021 #define RANDOM_H
00022 
00023 #include "common.h"
00024 #include "Vector.h"
00025 
00026 #ifdef _MSC_VER
00027 #define INT64_LITERAL(X) X ## i64
00028 #else
00029 #define INT64_LITERAL(X) X ## LL
00030 #endif
00031 
00032 #define RAND48_SEED   INT64_LITERAL(0x00001234abcd330e)
00033 #define RAND48_MULT   INT64_LITERAL(0x00000005deece66d)
00034 #define RAND48_ADD    INT64_LITERAL(0x000000000000000b)
00035 #define RAND48_MASK   INT64_LITERAL(0x0000ffffffffffff)
00036 
00037 class Random {
00038 
00039 private:
00040 
00041   double second_gaussian;
00042   int64 second_gaussian_waiting;
00043   int64 rand48_seed;
00044   int64 rand48_mult;
00045   int64 rand48_add;
00046 
00047 public:
00048 
00049   // default constructor
00050   Random(void) {
00051     init(0);
00052     rand48_seed = RAND48_SEED;
00053   }
00054 
00055   // constructor with seed
00056   Random(unsigned long seed) {
00057     init(seed);
00058   }
00059 
00060   // reinitialize with seed
00061   void init(unsigned long seed) {
00062     second_gaussian = 0;
00063     second_gaussian_waiting = 0;
00064     rand48_seed = seed & INT64_LITERAL(0x00000000ffffffff);
00065     rand48_seed = rand48_seed << 16;
00066     rand48_seed |= RAND48_SEED & INT64_LITERAL(0x0000ffff);
00067     rand48_mult = RAND48_MULT;
00068     rand48_add = RAND48_ADD;
00069   }
00070 
00071   // advance generator by one (seed = seed * mult + add, to 48 bits)
00072   void skip(void) {
00073     rand48_seed = ( rand48_seed * rand48_mult + rand48_add ) & RAND48_MASK;
00074   }
00075 
00076   // split into numStreams different steams and take stream iStream
00077   void split(int iStream, int numStreams) {
00078 
00079     int i;
00080 
00081     // make sure that numStreams is odd to ensure maximum period
00082     numStreams |= 1;
00083 
00084     // iterate to get to the correct stream
00085     for ( i = 0; i < iStream; ++i ) skip();
00086 
00087     // save seed and add so we can use skip() for our calculations
00088     int64 save_seed = rand48_seed;
00089 
00090     // calculate c *= ( 1 + a + ... + a^(numStreams-1) )
00091     rand48_seed = rand48_add;
00092     for ( i = 1; i < numStreams; ++i ) skip();
00093     int64 new_add = rand48_seed;
00094 
00095     // calculate a = a^numStreams
00096     rand48_seed = rand48_mult;
00097     rand48_add  = 0;
00098     for ( i = 1; i < numStreams; ++i ) skip();
00099     rand48_mult = rand48_seed;
00100 
00101     rand48_add  = new_add;
00102     rand48_seed = save_seed;
00103 
00104     second_gaussian = 0;
00105     second_gaussian_waiting = 0;
00106   }
00107 
00108   // return a number uniformly distributed between 0 and 1
00109   BigReal uniform(void) {
00110     skip();
00111     const double exp48 = ( 1.0 / (double)(INT64_LITERAL(1) << 48) );
00112     return ( (double) rand48_seed * exp48 );
00113   }
00114 
00115   // return a number from a standard gaussian distribution
00116   BigReal gaussian(void) {
00117     BigReal fac, r, v1, v2;
00118 
00119     if (second_gaussian_waiting) {
00120       second_gaussian_waiting = 0;
00121       return second_gaussian;
00122     } else {
00123       r = 2.;                 // r >= 1.523e-8 ensures abs result < 6
00124       while (r >=1. || r < 1.523e-8) { // make sure we are within unit circle
00125         v1 = 2.0 * uniform() - 1.0;
00126         v2 = 2.0 * uniform() - 1.0;
00127         r = v1*v1 + v2*v2;
00128       }
00129       fac = sqrt(-2.0 * log(r)/r);
00130       // now make the Box-Muller transformation to get two normally
00131       // distributed random numbers. Save one and return the other.
00132       second_gaussian_waiting = 1;
00133       second_gaussian = v1 * fac;
00134       return v2 * fac;
00135     }
00136   }
00137 
00138   // return a vector of gaussian random numbers
00139   Vector gaussian_vector(void) {
00140     return Vector( gaussian(), gaussian(), gaussian() );
00141   }
00142 
00143   // return a random long
00144   // signed int32 ranges from 0 to (2^31)-1
00145   long integer(void) {
00146     skip();
00147     return ( ( rand48_seed >> 17 ) & INT64_LITERAL(0x000000007fffffff) );
00148   }
00149 
00150   // randomly order an array of whatever
00151   template <class Elem> void reorder(Elem *a, int n) {
00152     for ( int i = 0; i < (n-1); ++i ) {
00153       int ie = i + ( integer() % (n-i) );
00154       if ( ie == i ) continue;
00155       const Elem e = a[ie];
00156       a[ie] = a[i];
00157       a[i] = e;
00158     }
00159   }
00160 
00164   void uniform_array_f(float *a, int n) {
00165     // for now just call uniform()
00166     // ultimately we will provide faster implementation
00167     for (int i=0;  i < n;  i++) {
00168       a[i] = (float)uniform();
00169     }
00170   }
00171 
00175   void gaussian_array_f(float *a, int n) {
00176     // for now just call gaussian()
00177     // ultimately we will provide faster implementation
00178     for (int i=0;  i < n;  i++) {
00179       a[i] = (float)gaussian();
00180     }
00181   }
00182 
00183 };
00184 
00185 #endif  // RANDOM_H
00186 

Generated on Sat Nov 18 01:17:15 2017 for NAMD by  doxygen 1.4.7