00001 /***************************************************************************
00002                           random.cpp  -  description
00003                              -------------------
00004     begin                : Mon Jul 24 2000
00005     copyright            : (C) 2000 by Michael Peeters
00006     email                : Michael.Peeters@vub.ac.be
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  *                                                                         *
00011  *   This program is free software; you can redistribute it and/or modify  *
00012  *   it under the terms of the GNU General Public License as published by  *
00013  *   the Free Software Foundation; either version 2 of the License, or     *
00014  *   (at your option) any later version.                                   *
00015  *                                                                         *
00016  ***************************************************************************/
00017 
00018 #include "random.h"
00019 
00020 namespace MODEL {
00021 
00022   // Uniform
00023   //------------------------------------------------------------
00024   Uniform::Uniform(counter initial_seed) : Random(initial_seed)
00025   {
00026         // Load the shuffle table after 8 warmups
00027         for (counter j=tabelsize+7;j>=0;j--)
00028           if (j < tabelsize) tabel[j]=seed_update();
00029         lastseed=tabel[0];
00030   }
00031 
00032   const counter& Uniform::seed_update(void)
00033   {
00034         counter k=_seed/IQ;
00035         _seed=IA*(_seed-k*IQ)-IR*k;
00036         if(_seed<0) _seed+=IM;
00037         return _seed;
00038   }
00039 
00040   number& Uniform::generate(number& storage)
00041   {
00042         seed_update();
00043         counter shuffle=lastseed/NDIV;
00044         lastseed=tabel[shuffle];
00045         tabel[shuffle]=_seed;
00046 
00047         if ((storage=AM*lastseed)>RNMX) return storage=RNMX; // no end values
00048         else return storage;
00049   }
00050 
00051   // Normal
00052   //------------------------------------------------------------
00053         
00054   Normal::Normal(counter seed) 
00055         : regenerate2(true)
00056   {
00057         rnd = new Uniform(seed);
00058   }
00059 
00060   Normal::~Normal()
00061   {
00062         delete rnd;
00063   }
00064 
00065   number& Normal::generate(number& storage)
00066   {
00067         number v1,v2,rsq,fac;
00068         if(regenerate2) // we need new numbers
00069           {
00070                 do { // until we find a pair inside the unit circle
00071                   v1=2.0*(*rnd)()-1.0;
00072                   v2=2.0*(*rnd)()-1.0;
00073                   rsq=v1*v1+v2*v2;
00074                 } while (rsq >=1.0 || rsq==0.0);
00075                 fac=sqrt(-2.0*log(rsq)/rsq);
00076                 saved=v1*fac; regenerate2=false;
00077                 return storage=v2*fac;
00078           }
00079         else 
00080           {
00081                 regenerate2=true;
00082                 return storage=saved;
00083           }
00084   }
00085 
00086 } // end namespace
00087 
00088 /*********************************************************************
00089 $Id: random.cpp,v 1.1 2001/05/22 10:54:55 mpeeters Exp $
00090 **********************************************************************
00091 
00092 $Log: random.cpp,v $
00093 Revision 1.1  2001/05/22 10:54:55  mpeeters
00094 Moved sources and headers for libModel to model/ subdirectory, in an attempt to rationalize the source tree. This should make things "netter".
00095 
00096 Revision 1.2  2000/09/15 10:26:31  mpeeters
00097 Cleaned out header and added CVS tails to files, removed superfuous
00098 @author comments, inserted dates
00099 
00100 *********************************************************************/