planck_rng.h

00001 /*
00002  *  This file is part of Healpix_cxx.
00003  *
00004  *  Healpix_cxx is free software; you can redistribute it and/or modify
00005  *  it under the terms of the GNU General Public License as published by
00006  *  the Free Software Foundation; either version 2 of the License, or
00007  *  (at your option) any later version.
00008  *
00009  *  Healpix_cxx is distributed in the hope that it will be useful,
00010  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012  *  GNU General Public License for more details.
00013  *
00014  *  You should have received a copy of the GNU General Public License
00015  *  along with Healpix_cxx; if not, write to the Free Software
00016  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00017  *
00018  *  For more information about HEALPix, see http://healpix.jpl.nasa.gov
00019  */
00020 
00021 /*
00022  *  Healpix_cxx is being developed at the Max-Planck-Institut fuer Astrophysik
00023  *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
00024  *  (DLR).
00025  */
00026 
00027 /*! \file planck_rng.h
00028  *  This file contains the random number generator
00029  *  used by the Planck LevelS package.
00030  *  The generator is a C++ port of the xorshift generator xor128() described
00031  *  in Marsaglia, Journal of Statistical Software 2003, vol 8.
00032  *  It has a period of 2^128 - 1.
00033  *
00034  *  Copyright (C) 2003, 2004, 2005 Max-Planck-Society
00035  *  \author Martin Reinecke
00036  */
00037 
00038 #ifndef PLANCK_RNG_H
00039 #define PLANCK_RNG_H
00040 
00041 #include <cmath>
00042 #include "cxxutils.h"
00043 
00044 /*! C++ port of the xorshift generator xor128() described in Marsaglia,
00045     Journal of Statistical Software 2003, vol 8.
00046     It has a period of 2^128 - 1. */
00047 class planck_rng
00048   {
00049   private:
00050     unsigned int x,y,z,w;
00051     double small, gset;
00052     bool empty;
00053 
00054     void twiddle (unsigned int &v)
00055       {
00056       for (int i=0; i<9; ++i)
00057         {
00058         v ^= v<<13;
00059         v ^= v>>17;
00060         v ^= v<<5;
00061         }
00062       }
00063 
00064     void init_rng ()
00065       {
00066       // avoid zero seeds
00067       if (x==0) x = 123456789;
00068       if (y==0) y = 362436069;
00069       if (z==0) z = 521288629;
00070       if (w==0) w = 88675123;
00071 
00072       // shuffle the bits of the seeds
00073       twiddle(x); twiddle(y); twiddle(z); twiddle(w);
00074 
00075       // burn in the RNG
00076       for (int i=0; i<16; ++i)
00077         int_rand_uni();
00078       }
00079 
00080   public:
00081     /*! Initializes the generator with 0 to 4 seed values. */
00082     planck_rng (unsigned int x1=123456789, unsigned int y1=362436069,
00083                 unsigned int z1=521288629, unsigned int w1=88675123)
00084       : x(x1), y(y1), z(z1), w(w1),
00085         small(1./(1.+double(0xFFFFFFFF))), empty(true)
00086       {
00087       planck_assert (sizeof(unsigned int)==4, "wrong integer size for RNG");
00088       init_rng();
00089       }
00090 
00091     /*! Re-initializes the generator with 0 to 4 seed values. */
00092     void seed (unsigned int x1=123456789, unsigned int y1=362436069,
00093       unsigned int z1=521288629, unsigned int w1=88675123)
00094       {
00095       x = x1; y = y1; z = z1; w = w1;
00096       empty = true;
00097       init_rng();
00098       }
00099 
00100     /*! Returns uniformly distributed random integer numbers from the
00101         interval [0;0xFFFFFFFF]. */
00102     unsigned int int_rand_uni()
00103       {
00104       unsigned int t = x^(x<<11);
00105       x = y;
00106       y = z;
00107       z = w;
00108 
00109       return w=(w^(w>>19))^(t^(t>>8));
00110       }
00111 
00112     //! Returns uniformly distributed random numbers from the interval [0;1[.
00113     double rand_uni()
00114       {
00115       return small*int_rand_uni();
00116       }
00117 
00118     //! Returns random numbers with Gaussian distribution (mean=0, sigma=1).
00119     /*! Uses rand_uni() internally. */
00120     double rand_gauss()
00121       {
00122       using namespace std;
00123       if (empty)
00124         {
00125         double v1,v2,rsq;
00126         do
00127           {
00128           v1=2*rand_uni()-1.;
00129           v2=2*rand_uni()-1.;
00130           rsq=v1*v1+v2*v2;
00131           }
00132         while ((rsq>=1) || (rsq==0));
00133         double fac=sqrt(-2*log(rsq)/rsq);
00134         gset=v1*fac;
00135         empty=false;
00136         return v2*fac;
00137         }
00138       else
00139         {
00140         empty=true;
00141         return gset;
00142         }
00143       }
00144 
00145     //! Returns exponentially distributed random numbers (mean=1, nonnegative)
00146     /*! Uses rand_uni() internally. */
00147     double rand_exp()
00148       {
00149       using namespace std;
00150       double val=rand_uni();
00151       if (val==0.) val=1.;
00152       return -log(val);
00153       }
00154   };
00155 
00156 #endif