Scopira 20080306

distrandom.h

00001 
00002 /*
00003  *  Copyright (c) 2002    National Research Council
00004  *
00005  *  All rights reserved.
00006  *
00007  *  This material is confidential and proprietary information of
00008  *  National Research Council Canada ("Confidential Information").
00009  *  This Confidential Information may only be used and reproduced
00010  *  in accordance with the terms of the license agreement.
00011  *
00012  */
00013 
00014 #ifndef __INCLUDED__SCOPIRA_TOOL_DISTRANDOM_H__
00015 #define __INCLUDED__SCOPIRA_TOOL_DISTRANDOM_H__
00016 
00017 #include <scopira/tool/limits.h>
00018 #include <scopira/tool/linconrandom.h>
00019 
00020 namespace scopira
00021 {
00022   namespace tool
00023   {
00024     template <class UniRanGen, class RealType = double>
00025       class uni01_dist;
00026     template <class UniRanGen, class RealType = double>
00027       class unireal_dist;
00028     template <class UniRanGen, class IntType = int>
00029       class unismallint_dist;
00030     template<class UniRanGen, class IntType = int>
00031       class uniint_dist;
00032 
00033     // "typical" distrbutions
00034     typedef uni01_dist<lincon_core> lincon_01;
00035     typedef unireal_dist<lincon_core> lincon_real;
00036     typedef unismallint_dist<lincon_core> lincon_smallint;  // how small is small? who knows
00037     typedef uniint_dist<lincon_core> lincon_int;
00038 
00039     typedef uni01_dist<scopira::tool::rand48> rand48_01;
00040     typedef unireal_dist<scopira::tool::rand48> rand48_real;
00041     typedef unismallint_dist<scopira::tool::rand48> rand48_smallint;  // how small is small? who knows
00042     typedef uniint_dist<scopira::tool::rand48> rand48_int;
00043   }
00044 }
00045 
00046 // various random distributions
00047 
00056 template <class UniRanGen, class RealType>
00057 class scopira::tool::uni01_dist
00058 {
00059   public:
00060     typedef UniRanGen base_type;
00061     typedef RealType result_type;
00062     static const bool has_fixed_range = true;
00063 
00065     explicit uni01_dist(base_type &rng) : m_rng(rng) { }
00066 
00068     result_type min() const { return 0.0; }
00070     result_type max() const { return 1.0; }
00071 
00072     result_type operator()(void) {
00073       return static_cast<result_type>(m_rng() - m_rng.min()) /
00074         (static_cast<result_type>(m_rng.max()-m_rng.min()) +
00075          (std::numeric_limits<base_result>::is_integer ? 1.0 : 0.0));
00076     }
00077 
00078   private:
00079     typedef typename base_type::result_type base_result;
00080     base_type & m_rng;
00081 };
00082 
00088 template <class UniRanGen, class RealType>
00089 class scopira::tool::unireal_dist
00090 {
00091   public:
00092     typedef UniRanGen base_type;
00093     typedef RealType result_type;
00094     static const bool has_fixed_range = false;
00095 
00097     unireal_dist(base_type &rng, RealType min, RealType max)
00098       : m_rng(rng), m_min(min), m_max(max) { }
00099 
00101     result_type min(void) const { return m_min; }
00103     result_type max(void) const { return m_max; }
00104 
00105     result_type operator()() { return m_rng() * (m_max - m_min) + m_min; }
00106 
00107   private:
00108     uni01_dist<base_type, result_type> m_rng;
00109     RealType m_min, m_max;
00110 };
00111 
00117 template <class UniRanGen, class IntType>
00118 class scopira::tool::unismallint_dist
00119 {
00120   public:
00121     typedef UniRanGen base_type;
00122     typedef IntType result_type;
00123     static const bool has_fixed_range = false;
00124 
00125   public:
00127     unismallint_dist(base_type &rng, IntType min, IntType max);
00128 
00130     result_type min(void) const { return m_min; }
00132     result_type max(void) const { return m_max; }
00133 
00135     result_type operator()(void) {
00136       // we must not use the low bits here, because LCGs get very bad then
00137       return long((m_rng() - m_rng.min()) / m_factor) % long(m_range) + m_min;  // MDA  casts used to prevent bad types
00138     }
00139 
00140   private:
00141     typedef typename base_type::result_type base_result;
00142     base_type & m_rng;
00143     IntType m_min, m_max;
00144     base_result m_range;
00145     int m_factor;
00146 };
00147 
00149 template <class UniRanGen, class IntType>
00150   scopira::tool::unismallint_dist<UniRanGen, IntType>::unismallint_dist
00151   (base_type & rng, IntType min, IntType max) 
00152   : m_rng(rng), m_min(min), m_max(max),
00153     m_range(static_cast<base_result>(m_max-m_min)+1), m_factor(1)
00154 {
00155   //assert(min < max);
00156   
00157   // LCGs get bad when only taking the low bits.
00158   // (probably put this logic into a partial template specialization)
00159   // Check how many low bits we can ignore before we get too much
00160   // quantization error.
00161   base_result r_base = m_rng.max() - m_rng.min();
00162   if (r_base == std::numeric_limits<base_result>::max()) {
00163     m_factor = 2;
00164     r_base /= 2;
00165   }
00166   r_base += 1;
00167   if(long(r_base) % long(m_range) == 0) {   // MDA invalid double operators to %, added casts
00168     // No quantization effects, good
00169     m_factor = int(double(r_base) / double(m_range));  
00170   } else {
00171     //const base_result r = 32 * m_range * m_range;
00172     //for(; r_base >= r; m_factor *= 2)
00173       //r_base /= 2;
00174     // carefully avoid overflow; pessimizing heree
00175     for ( ; r_base/m_range/32 >= m_range; m_factor *= 2)
00176       r_base /= 2;
00177   }
00178 }
00179 
00180 template<class UniRanGen, class IntType>
00181 class scopira::tool::uniint_dist
00182 {
00183   public:
00184     typedef UniRanGen base_type;
00185     typedef IntType input_type;
00186     typedef IntType result_type;
00187     static const bool has_fixed_range = false;
00188 
00189     explicit uniint_dist(UniRanGen &rng, IntType min, IntType max)
00190       : m_rng(rng), m_min(min), m_max(max) { m_range = m_max - m_min; }
00191 
00193     result_type min(void) const { return m_min; }
00195     result_type max(void) const { return m_max; }
00196   
00198     result_type operator()(void);
00199 
00200   private:
00201     typedef typename base_type::result_type base_result;
00202     base_type & m_rng;
00203     IntType m_min, m_max, m_range;
00204 };
00205 
00206 template<class UniRanGen, class IntType>
00207 IntType scopira::tool::uniint_dist<UniRanGen, IntType>::operator ()(void)
00208 {
00209   base_result bmin = (m_rng.min)();
00210   base_result brange = (m_rng.max)() - (m_rng.min)();
00211 
00212   if (m_range == 0) {
00213     return m_min;    
00214   } else if (brange == m_range) {
00215     // this will probably never happen in real life
00216     // basically nothing to do; just take care we don't overflow / underflow
00217     return static_cast<result_type>(m_rng() - bmin) + m_min;
00218   } else if (brange < m_range) {
00219     // use rejection method to handle things like 0..3 --> 0..4
00220     for(;;) {
00221       // concatenate several invocations of the base RNG
00222       // take extra care to avoid overflows
00223       result_type limit;
00224       if(m_range == (std::numeric_limits<result_type>::max)()) {
00225         limit = m_range/(result_type(brange)+1);
00226         if(m_range % result_type(brange)+1 == result_type(brange))
00227           ++limit;
00228       } else {
00229         limit = (m_range+1)/(result_type(brange)+1);
00230       }
00231       // We consider "result" as expressed to base (brange+1):
00232       // For every power of (brange+1), we determine a random factor
00233       result_type result = result_type(0);
00234       result_type mult = result_type(1);
00235       while(mult <= limit) {
00236         result += (m_rng() - bmin) * mult;
00237         mult *= result_type(brange)+result_type(1);
00238       }
00239       if(mult == limit)
00240         // m_range+1 is an integer power of brange+1: no rejections required
00241         return result;
00242       // m_range/mult < brange+1  -> no endless loop
00243       result += uniint_dist<UniRanGen, result_type>(m_rng, 0, m_range/mult)() * mult;
00244       if(result <= m_range)
00245         return result + m_min;
00246     }
00247   } else {                   // brange > range
00248     if (brange / m_range > 4 /* quantization_cutoff */ ) {
00249       // the new range is vastly smaller than the source range,
00250       // so quantization effects are not relevant
00251       return unismallint_dist<UniRanGen, result_type>(m_rng, m_min, m_max)();
00252     } else {
00253       // use rejection method to handle cases like 0..5 -> 0..4
00254       for(;;) {
00255         base_result result = m_rng() - bmin;
00256         // result and range are non-negative, and result is possibly larger
00257         // than range, so the cast is safe
00258         if(result <= static_cast<base_result>(m_range))
00259           return result + m_min;
00260       }
00261     }
00262   }
00263 }
00264 
00265 #endif
00266