estarlight/src/randomgenerator.cpp

0001 ///////////////////////////////////////////////////////////////////////////
0002 //
0003 //    Copyright 2010
0004 //
0005 //    This file is part of starlight.
0006 //
0007 //    starlight is free software: you can redistribute it and/or modify
0008 //    it under the terms of the GNU General Public License as published by
0009 //    the Free Software Foundation, either version 3 of the License, or
0010 //    (at your option) any later version.
0011 //
0012 //    starlight is distributed in the hope that it will be useful,
0013 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
0014 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0015 //    GNU General Public License for more details.
0016 //
0017 //    You should have received a copy of the GNU General Public License
0018 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
0019 //
0020 ///////////////////////////////////////////////////////////////////////////
0021 //
0022 // File and Version Information:
0023 // $Rev:: 213                         $: revision of last commit
0024 // $Author:: butter                   $: author of last commit
0025 // $Date:: 2015-08-15 22:08:02 +0100 #$: date of last commit
0026 //
0027 // Description:
0028 //
0029 //
0030 //
0031 ///////////////////////////////////////////////////////////////////////////
0032
0033
0034 #include <iostream>
0035 #include <fstream>
0036 #include <cmath>
0037
0038 #include "randomgenerator.h"
0039
0040
0041 using namespace std;
0042
0043
0044 //USED IN ROOT under TRANDOM3
0045 // Random number generator class based on
0046 //   M. Matsumoto and T. Nishimura,
0047 //   Mersenne Twistor: A 623-diminsionally equidistributed
0048 //   uniform pseudorandom number generator
0049 //   ACM Transactions on Modeling and Computer Simulation,
0050 //   Vol. 8, No. 1, January 1998, pp 3--30.
0051 //
0052 // For more information see the Mersenne Twistor homepage
0053 //   http://www.math.keio.ac.jp/~matumoto/emt.html
0054 //
0055 // Advantage: large period 2**19937-1
0056 //            relativly fast
0057 //              (only two times slower than TRandom, but
0058 //               two times faster than TRandom2)
0059 // Drawback:  a relative large internal state of 624 integers
0060 //
0061 //
0062 // Aug.99 ROOT implementation based on CLHEP by P.Malzacher
0063 //
0064 // the original code contains the following copyright notice:
0065 /* This library is free software; you can redistribute it and/or   */
0066 /* modify it under the terms of the GNU Library General Public     */
0067 /* License as published by the Free Software Foundation; either    */
0068 /* version 2 of the License, or (at your option) any later         */
0069 /* version.                                                        */
0070 /* This library is distributed in the hope that it will be useful, */
0071 /* but WITHOUT ANY WARRANTY; without even the implied warranty of  */
0072 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.            */
0073 /* See the GNU Library General Public License for more details.    */
0074 /* You should have received a copy of the GNU Library General      */
0075 /* Public License along with this library; if not, write to the    */
0076 /* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA   */
0077 /* 02111-1307  USA                                                 */
0078 /* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.       */
0079 /* When you use this, send an email to: matumoto@math.keio.ac.jp   */
0080 /* with an appropriate reference to your work.                     */
0081 /////////////////////////////////////////////////////////////////////
0082
0083
0084 void randomGenerator::SetSeed(unsigned int seed)
0085 {
0086 //  Set the random generator sequence
0087 // if seed is 0 (default value) a TUUID is generated and used to fill
0088 // the first 8 integers of the seed array.
0089 // In this case the seed is guaranteed to be unique in space and time.
0090 // Use upgraded seeding procedure to fix a known problem when seeding with values
0091 // with many zero in the bit pattern (like 2**28).
0092 // see http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
0093
0094
0095   _count624 = 624;
0096   int i,j;
0097
0098   _Mt[0] = seed;
0099   j = 1;
0100   // use multipliers from  Knuth's "Art of Computer Programming" Vol. 2, 3rd Ed. p.106
0101   for(i=j; i<624; i++) {
0102     _Mt[i] = (1812433253 * ( _Mt[i-1]  ^ ( _Mt[i-1] >> 30)) + i);
0103   }
0104 }
0105
0106
0107 double randomGenerator::Rndom(int)
0108 {
0109
0110 //  Machine independent random number generator.
0111 //  Produces uniformly-distributed floating points in [0,1]
0112 //  Method: Mersenne Twistor
0113
0114
0115    unsigned int y;
0116
0117    const int  kM = 397;
0118    const int  kN = 624;
0119    const unsigned int kTemperingMaskB =  0x9d2c5680;
0120    const unsigned int kTemperingMaskC =  0xefc60000;
0121    const unsigned int kUpperMask =       0x80000000;
0122    const unsigned int kLowerMask =       0x7fffffff;
0123    const unsigned int kMatrixA =         0x9908b0df;
0124
0125    if (_count624 >= kN) {
0126       int i;
0127
0128       for (i=0; i < kN-kM; i++) {
0129          y = (_Mt[i] & kUpperMask) | (_Mt[i+1] & kLowerMask);
0130          _Mt[i] = _Mt[i+kM] ^ (y >> 1) ^ ((y & 0x1) ? kMatrixA : 0x0);
0131       }
0132
0133       for (   ; i < kN-1    ; i++) {
0134          y = (_Mt[i] & kUpperMask) | (_Mt[i+1] & kLowerMask);
0135          _Mt[i] = _Mt[i+kM-kN] ^ (y >> 1) ^ ((y & 0x1) ? kMatrixA : 0x0);
0136       }
0137
0138       y = (_Mt[kN-1] & kUpperMask) | (_Mt[0] & kLowerMask);
0139       _Mt[kN-1] = _Mt[kM-1] ^ (y >> 1) ^ ((y & 0x1) ? kMatrixA : 0x0);
0140       _count624 = 0;
0141    }
0142
0143    y = _Mt[_count624++];
0144    y ^=  (y >> 11);
0145    y ^= ((y << 7 ) & kTemperingMaskB );
0146    y ^= ((y << 15) & kTemperingMaskC );
0147    y ^=  (y >> 18);
0148
0149    if (y) return ( (double) y * 2.3283064365386963e-10); // * Power(2,-32)
0150    return Rndom();
0151 }