knuthlfg.h

Go to the documentation of this file.

/*
 *  knuthlfg.h
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  NEST is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#ifndef KNUTHLFG_H
#define KNUTHLFG_H

/*
 *  Built-in implementation of Knuth's LFG generator.
 *  This code is a C++ adaptation of the code published by Knuth on his
 *  website, http://www-cs-faculty.stanford.edu/~knuth/programs/rng.c,
 *  retrieved 8 Jan 2008. See also Knuth's header comment below.
 */
 
/* Header comment by  D E Knuth ------------------------------------------ */

/*    This program by D E Knuth is in the public domain and freely copyable.
 *    It is explained in Seminumerical Algorithms, 3rd edition, Section 3.6
 *    (or in the errata to the 2nd edition --- see
 *        http://www-cs-faculty.stanford.edu/~knuth/taocp.html
 *    in the changes to Volume 2 on pages 171 and following).              */

/*    N.B. The MODIFICATIONS introduced in the 9th printing (2002) are
      included here; there's no backwards compatibility with the original. */

/*    This version also adopts Brendan McKay's suggestion to
      accommodate naive users who forget to call ran_start(seed).          */

/*    If you find any bugs, please report them immediately to
 *                 taocp@cs.stanford.edu
 *    (and you will be rewarded if the bug is genuine). Thanks!            */

/************ see the book for explanations and caveats! *******************/
/************ in particular, you need two's complement arithmetic **********/

/* End of Header comment by  D E Knuth ----------------------------------- */

#include "randomgen.h"
#include <vector>

namespace librandom {
  
  class KnuthLFG : public RandomGen {
  public:
    
    explicit KnuthLFG(unsigned long);
    
    ~KnuthLFG() {};

    RngPtr clone(unsigned long s)
      {
    return RngPtr(new KnuthLFG(s));
      }

  private:
    void   seed_(unsigned long);
    
    double drand_();
    
  private:  
    static const long    KK_;        
    static const long    LL_;        
    static const long    MM_;        
    static const long    TT_;        
    static const long    QUALITY_;   
    static const double  I2DFactor_; 

    static long mod_diff_(long, long);  
    static bool is_odd_(long);
    
    std::vector<long> ran_x_;       
    std::vector<long> ran_buffer_;  
    const std::vector<long>::const_iterator end_;  
    std::vector<long>::const_iterator next_; 
    
    void ran_array_(std::vector<long>& rbuff);
    void ran_start_(long seed); 
    long ran_draw_();           
    
    void self_test_(); 
  };  

  inline
  void KnuthLFG::seed_(unsigned long seed)
  {  
    ran_start_(seed);
  }
  

  inline
  double KnuthLFG::drand_()
  {  
    return I2DFactor_ * ran_draw_();
  }
  
  
  inline
  long KnuthLFG::mod_diff_(long x, long y)
  {
    // modulo computation assumes two's complement
    return ( x - y ) & ( MM_ - 1 );
  }
  
  inline
  bool KnuthLFG::is_odd_(long x)
  {
    return x & 1;
  }
  
  inline long KnuthLFG::ran_draw_()
  {
    if ( next_ == end_ )
    {
      ran_array_(ran_buffer_);  // refill
      next_ = ran_buffer_.begin();
    }
    
    return *next_++;   // return next and increment
  }

}  // namespace librandom

#endif