Main Page   Modules   Class Hierarchy   Compound List   File List   Compound Members   Related Pages   Examples  

Rnd_no_gen.cpp

00001 /* A C-program for MT19937: Real number version                */
00002 /*   GenRand() generates one pseudorandom real number (double) */
00003 /* which is uniformly distributed on [0,1]-interval, for each  */
00004 /* call. SetGenRand(seed) set initial values to the working area */
00005 /* of 624 words. Before GenRand(), SetGenRand(seed) must be      */
00006 /* called once. (seed is any 32-bit integer except for 0).     */
00007 /* Integer generator is obtained by modifying two lines.       */
00008 /*   Coded by Takuji Nishimura, considering the suggestions by */
00009 /* Topher Cooper and Marc Rieffel in July-Aug. 1997.           */
00010 
00011 /* This library is free software; you can redistribute it and/or   */
00012 /* modify it under the terms of the GNU Library General Public     */
00013 /* License as published by the Free Software Foundation; either    */
00014 /* version 2 of the License, or (at your option) any later         */
00015 /* version.                                                        */
00016 /* This library is distributed in the hope that it will be useful, */
00017 /* but WITHOUT ANY WARRANTY; without even the implied warranty of  */
00018 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.            */
00019 /* See the GNU Library General Public License for more details.    */
00020 /* You should have received a copy of the GNU Library General      */
00021 /* Public License along with this library; if not, write to the    */
00022 /* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA   */ 
00023 /* 02111-1307  USA                                                 */
00024 
00025 /* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.       */
00026 /* Any feedback is very welcome. For any question, comments,       */
00027 /* see http://www.math.keio.ac.jp/matumoto/emt.html or email       */
00028 /* matumoto@math.keio.ac.jp                                        */
00029 
00030 #include<stdio.h>
00031 #include <math.h>
00032 #include "Rnd_no_gen.h"
00033 
00034 #define N 624
00035 #define M 397
00036 #define MATRIX_A 0x9908b0df   
00037 #define UPPER_MASK 0x80000000 
00038 #define LOWER_MASK 0x7fffffff 
00039 
00040 #define TEMPERING_MASK_B 0x9d2c5680
00041 #define TEMPERING_MASK_C 0xefc60000
00042 #define TEMPERING_SHIFT_U(y)  (y >> 11)
00043 #define TEMPERING_SHIFT_S(y)  (y << 7)
00044 #define TEMPERING_SHIFT_T(y)  (y << 15)
00045 #define TEMPERING_SHIFT_L(y)  (y >> 18)
00046 
00047 static unsigned long mt[N]; 
00048 static int mti=N+1; 
00049 
00050 void SetGenRand(unsigned long seed)
00051 {
00052   mt[0]= seed & 0xffffffff;
00053   for (mti=1; mti<N; mti++)
00054     mt[mti] = (69069 * mt[mti-1]) & 0xffffffff;
00055 }
00056 
00057 double GenRand()
00058 {
00059   unsigned long y;
00060   static unsigned long mag01[2]={0x0, MATRIX_A};
00061 
00062   if (mti >= N) { 
00063     int kk;
00064 
00065     if (mti == N+1)   
00066       SetGenRand(4357); 
00067 
00068     for (kk=0;kk<N-M;kk++) {
00069       y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
00070       mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
00071     }
00072     for (;kk<N-1;kk++) {
00073       y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
00074       mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
00075     }
00076     y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
00077     mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
00078     
00079     mti = 0;
00080   }
00081   
00082   y = mt[mti++];
00083   y ^= TEMPERING_SHIFT_U(y);
00084   y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
00085   y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
00086   y ^= TEMPERING_SHIFT_L(y);
00087   
00088   return ( (double)y / (unsigned long)0xffffffff ); 
00089 }
00090 
00091 double GenRand(double lo, double hi)
00092 {
00093   double res;
00094   res = lo + (hi - lo)*GenRand();
00095   return res;
00096 }
00097 
00098 #define IA 16807
00099 #define IM 2147483647
00100 #define AM (1.0/IM)
00101 #define IQ 127773
00102 #define IR 2836
00103 #define NTAB
00104 
00105 double GaussRand()
00106 {
00107   static int iset = 0;
00108   static double gset;
00109   double fac, rsq, v1, v2;
00110 
00111   if (iset == 0) {
00112     do {
00113       v1 = 2.0*GenRand() - 1.0;
00114       v2 = 2.0*GenRand() - 1.0;
00115       rsq = v1*v1 +v2*v2;
00116   } while (rsq >= 1.0 || rsq == 0.0);
00117     fac = sqrt(-2.0*log(rsq)/rsq);
00118     gset = v1*fac;
00119     iset = 1;
00120     return v2*fac;
00121   } else {
00122     iset = 0;
00123     return gset;
00124   }
00125 }
00126 
00127 double GaussRand(double mean, double stdev)
00128 {
00129   double res;
00130   res = mean + stdev*GaussRand();
00131   return res;
00132 }
00133 
00134 
00135 
00136 
00137 
00138