libs/CppAddons/Random.cpp

   1 // Copyright 2003 "Gilles Degottex"
   2
   3 // This file is part of "CppAddons"
   4
   5 // "CppAddons" is free software; you can redistribute it and/or modify
   6 // it under the terms of the GNU Lesser General Public License as published by
   7 // the Free Software Foundation; either version 2.1 of the License, or
   8 // (at your option) any later version.
   9 //
  10 // "CppAddons" is distributed in the hope that it will be useful,
  11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  13 // GNU Lesser General Public License for more details.
  14 //
  15 // You should have received a copy of the GNU Lesser General Public License
  16 // along with this program; if not, write to the Free Software
  17 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  18
  19
  20 #include "Random.h"
  21 #include <stdlib.h>
  22 #include <cmath>
  23
  24 const int Random::A = 48271;
  25 //const int Random::M = 2147483647;
  26 const int Random::M = RAND_MAX;
  27 const int Random::Q = M / A;
  28 const int Random::R = M % A;
  29
  30 Random Random::s_random;
  31
  32 Random::Random(long seed)
  33 {
  34         m_haveNextNextGaussian=false;
  35
  36         setSeed(seed);
  37
  38         next();
  39 }
  40
  41 void Random::setSeed(long seed)
  42 {
  43         srand(seed);
  44
  45         if(seed < 0)    seed += M+1;
  46         m_seed = seed;
  47         if(m_seed==0)   m_seed = 1;
  48
  49         m_haveNextNextGaussian = false;
  50 }
  51
  52 long Random::next()
  53 {
  54 //      int tmp = A * (m_seed % Q) - R * (m_seed / Q);
  55 //      if(tmp>=0)      m_seed = tmp;
  56 //      else            m_seed = tmp + M;
  57 //
  58 //      return m_seed;
  59
  60         return rand();
  61 }
  62
  63 double Random::nextGaussian()
  64 {
  65 //      Throw("nextGaussian", "not yet implemented");
  66
  67         // See Knuth, ACP, Section 3.4.1 Algorithm C.
  68         if(m_haveNextNextGaussian)
  69         {
  70         m_haveNextNextGaussian = false;
  71         return m_nextNextGaussian;
  72         }
  73         else
  74         {
  75                 double v1, v2, s;
  76         do
  77                 {
  78                         v1 = 2 * nextDouble() - 1; // between -1 and 1
  79                         v2 = 2 * nextDouble() - 1; // between -1 and 1
  80                         s = v1 * v1 + v2 * v2;
  81         }
  82                 while (s >= 1 || s == 0);
  83
  84                 double multiplier = sqrt(-2 * log(s)/s);
  85
  86         m_nextNextGaussian = v2 * multiplier;
  87                 m_haveNextNextGaussian = true;
  88
  89                 return v1 * multiplier;
  90         }
  91 }
  92