STARLIGHT/starlight/src/randomgenerator.cpp

   1 ///////////////////////////////////////////////////////////////////////////
   2 //
   3 //    Copyright 2010
   4 //
   5 //    This file is part of starlight.
   6 //
   7 //    starlight is free software: you can redistribute it and/or modify
   8 //    it under the terms of the GNU General Public License as published by
   9 //    the Free Software Foundation, either version 3 of the License, or
  10 //    (at your option) any later version.
  11 //
  12 //    starlight is distributed in the hope that it will be useful,
  13 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 //    GNU General Public License for more details.
  16 //
  17 //    You should have received a copy of the GNU General Public License
  18 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
  19 //
  20 ///////////////////////////////////////////////////////////////////////////
  21 //
  22 // File and Version Information:
  23 // $Rev:: 156                         $: revision of last commit
  24 // $Author:: odjuvsla                 $: author of last commit
  25 // $Date:: 2013-10-06 16:17:52 +0200 #$: date of last commit
  26 //
  27 // Description:
  28 //
  29 //
  30 //
  31 ///////////////////////////////////////////////////////////////////////////
  32
  33
  34 #include <iostream>
  35 #include <fstream>
  36 #include <cmath>
  37
  38 #include "randomgenerator.h"
  39
  40
  41 using namespace std;
  42
  43
  44 //USED IN ROOT under TRANDOM3
  45 // Random number generator class based on
  46 //   M. Matsumoto and T. Nishimura,
  47 //   Mersenne Twistor: A 623-diminsionally equidistributed
  48 //   uniform pseudorandom number generator
  49 //   ACM Transactions on Modeling and Computer Simulation,
  50 //   Vol. 8, No. 1, January 1998, pp 3--30.
  51 //
  52 // For more information see the Mersenne Twistor homepage
  53 //   http://www.math.keio.ac.jp/~matumoto/emt.html
  54 //
  55 // Advantage: large period 2**19937-1
  56 //            relativly fast
  57 //              (only two times slower than TRandom, but
  58 //               two times faster than TRandom2)
  59 // Drawback:  a relative large internal state of 624 integers
  60 //
  61 //
  62 // Aug.99 ROOT implementation based on CLHEP by P.Malzacher
  63 //
  64 // the original code contains the following copyright notice:
  65 /* This library is free software; you can redistribute it and/or   */
  66 /* modify it under the terms of the GNU Library General Public     */
  67 /* License as published by the Free Software Foundation; either    */
  68 /* version 2 of the License, or (at your option) any later         */
  69 /* version.                                                        */
  70 /* This library is distributed in the hope that it will be useful, */
  71 /* but WITHOUT ANY WARRANTY; without even the implied warranty of  */
  72 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.            */
  73 /* See the GNU Library General Public License for more details.    */
  74 /* You should have received a copy of the GNU Library General      */
  75 /* Public License along with this library; if not, write to the    */
  76 /* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA   */
  77 /* 02111-1307  USA                                                 */
  78 /* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.       */
  79 /* When you use this, send an email to: matumoto@math.keio.ac.jp   */
  80 /* with an appropriate reference to your work.                     */
  81 /////////////////////////////////////////////////////////////////////
  82
  83
  84 void randomGenerator::SetSeed(unsigned int seed)
  85 {
  86 //  Set the random generator sequence
  87 // if seed is 0 (default value) a TUUID is generated and used to fill
  88 // the first 8 integers of the seed array.
  89 // In this case the seed is guaranteed to be unique in space and time.
  90 // Use upgraded seeding procedure to fix a known problem when seeding with values
  91 // with many zero in the bit pattern (like 2**28).
  92 // see http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
  93
  94   Lockguard<MutexPosix> guard(&_mutex);
  95
  96   _count624 = 624;
  97   int i,j;
  98
  99   _Mt[0] = seed;
 100   j = 1;
 101   // use multipliers from  Knuth's "Art of Computer Programming" Vol. 2, 3rd Ed. p.106
 102   for(i=j; i<624; i++) {
 103     _Mt[i] = (1812433253 * ( _Mt[i-1]  ^ ( _Mt[i-1] >> 30)) + i);
 104   }
 105 }
 106
 107
 108 double randomGenerator::Rndom(int)
 109 {
 110
 111 //  Machine independent random number generator.
 112 //  Produces uniformly-distributed floating points in ]0,1]
 113 //  Method: Mersenne Twistor
 114
 115    Lockguard<MutexPosix> guard(&_mutex);
 116
 117    unsigned int y;
 118
 119    const int  kM = 397;
 120    const int  kN = 624;
 121    const unsigned int kTemperingMaskB =  0x9d2c5680;
 122    const unsigned int kTemperingMaskC =  0xefc60000;
 123    const unsigned int kUpperMask =       0x80000000;
 124    const unsigned int kLowerMask =       0x7fffffff;
 125    const unsigned int kMatrixA =         0x9908b0df;
 126
 127    if (_count624 >= kN) {
 128       register int i;
 129
 130       for (i=0; i < kN-kM; i++) {
 131          y = (_Mt[i] & kUpperMask) | (_Mt[i+1] & kLowerMask);
 132          _Mt[i] = _Mt[i+kM] ^ (y >> 1) ^ ((y & 0x1) ? kMatrixA : 0x0);
 133       }
 134
 135       for (   ; i < kN-1    ; i++) {
 136          y = (_Mt[i] & kUpperMask) | (_Mt[i+1] & kLowerMask);
 137          _Mt[i] = _Mt[i+kM-kN] ^ (y >> 1) ^ ((y & 0x1) ? kMatrixA : 0x0);
 138       }
 139
 140       y = (_Mt[kN-1] & kUpperMask) | (_Mt[0] & kLowerMask);
 141       _Mt[kN-1] = _Mt[kM-1] ^ (y >> 1) ^ ((y & 0x1) ? kMatrixA : 0x0);
 142       _count624 = 0;
 143    }
 144
 145    y = _Mt[_count624++];
 146    y ^=  (y >> 11);
 147    y ^= ((y << 7 ) & kTemperingMaskB );
 148    y ^= ((y << 15) & kTemperingMaskC );
 149    y ^=  (y >> 18);
 150
 151    if (y) return ( (double) y * 2.3283064365386963e-10); // * Power(2,-32)
 152    return Rndom();
 153 }