Changeset 209:765619b7cbb2 in lemon-1.2 for lemon/random.h

Timestamp:

07/13/08 20:51:02 (17 years ago)

Author:

Alpar Juttner <alpar@…>

Branch:

default

Phase:

public

Message:

Apply unify-sources.sh to the source tree

File:

• lemon/random.h (modified) (51 diffs)

lemon/random.h

@@ -1 @@
-/* -*- C++ -*-
- *
- * This file is a part of LEMON, a generic C++ optimization library
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
  *
  * Copyright (C) 2003-2008
@@ -22 +22 @@
  *
  * See the appropriate copyright notice below.
- * 
+ *
  * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
- * All rights reserved.                          
+ * All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
@@ -37 +37 @@
  *    documentation and/or other materials provided with the distribution.
  *
- * 3. The names of its contributors may not be used to endorse or promote 
- *    products derived from this software without specific prior written 
+ * 3. The names of its contributors may not be used to endorse or promote
+ *    products derived from this software without specific prior written
  *    permission.
  *
@@ -88 +88 @@
 
   namespace _random_bits {
-    
+
     template <typename _Word, int _bits = std::numeric_limits<_Word>::digits>
     struct RandomTraits {};
@@ -100 +100 @@
       static const int length = 624;
       static const int shift = 397;
-      
+
       static const Word mul = 0x6c078965u;
       static const Word arrayInit = 0x012BD6AAu;
@@ -168 +168 @@
           0x12345u, 0x23456u, 0x34567u, 0x45678u
         };
-    
+
         initState(seedArray, seedArray + 4);
       }
@@ -176 +176 @@
         static const Word mul = RandomTraits<Word>::mul;
 
-        current = state; 
+        current = state;
 
         Word *curr = state + length - 1;
@@ -202 +202 @@
         num = length > end - begin ? length : end - begin;
         while (num--) {
-          curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1)) 
+          curr[0] = (curr[0] ^ ((curr[1] ^ (curr[1] >> (bits - 2))) * mul1))
             + *it + cnt;
           ++it; ++cnt;
@@ -224 +224 @@
           }
         }
-        
+
         state[length - 1] = Word(1) << (bits - 1);
       }
-      
+
       void copyState(const RandomCore& other) {
         std::copy(other.state, other.state + length, state);
@@ -242 +242 @@
     private:
 
-  
+
       void fillState() {
         static const Word mask[2] = { 0x0ul, RandomTraits<Word>::mask };
@@ -248 +248 @@
         static const Word hiMask = RandomTraits<Word>::hiMask;
 
-        current = state + length; 
+        current = state + length;
 
         register Word *curr = state + length - 1;
         register long num;
-      
+
         num = length - shift;
         while (num--) {
@@ -270 +270 @@
       }
 
-  
+
       Word *current;
       Word state[length];
-      
-    };
-
-
-    template <typename Result, 
+
+    };
+
+
+    template <typename Result,
               int shift = (std::numeric_limits<Result>::digits + 1) / 2>
     struct Masker {
@@ -285 +285 @@
       }
     };
-    
+
     template <typename Result>
     struct Masker<Result, 1> {
@@ -293 +293 @@
     };
 
-    template <typename Result, typename Word, 
-              int rest = std::numeric_limits<Result>::digits, int shift = 0, 
+    template <typename Result, typename Word,
+              int rest = std::numeric_limits<Result>::digits, int shift = 0,
               bool last = rest <= std::numeric_limits<Word>::digits>
     struct IntConversion {
       static const int bits = std::numeric_limits<Word>::digits;
-    
+
       static Result convert(RandomCore<Word>& rnd) {
         return static_cast<Result>(rnd() >> (bits - rest)) << shift;
       }
-      
-    }; 
-
-    template <typename Result, typename Word, int rest, int shift> 
+
+    };
+
+    template <typename Result, typename Word, int rest, int shift>
     struct IntConversion<Result, Word, rest, shift, false> {
       static const int bits = std::numeric_limits<Word>::digits;
 
       static Result convert(RandomCore<Word>& rnd) {
-        return (static_cast<Result>(rnd()) << shift) | 
+        return (static_cast<Result>(rnd()) << shift) |
           IntConversion<Result, Word, rest - bits, shift + bits>::convert(rnd);
       }
@@ -317 +317 @@
 
     template <typename Result, typename Word,
-              bool one_word = (std::numeric_limits<Word>::digits < 
-                               std::numeric_limits<Result>::digits) >
+              bool one_word = (std::numeric_limits<Word>::digits <
+                               std::numeric_limits<Result>::digits) >
     struct Mapping {
       static Result map(RandomCore<Word>& rnd, const Result& bound) {
@@ -325 +325 @@
         Result num;
         do {
-          num = IntConversion<Result, Word>::convert(rnd) & mask; 
+          num = IntConversion<Result, Word>::convert(rnd) & mask;
         } while (num > max);
         return num;
@@ -351 +351 @@
         res *= res;
         if ((exp & 1) == 1) res *= static_cast<Result>(2.0);
-        return res; 
+        return res;
       }
     };
@@ -361 +361 @@
         res *= res;
         if ((exp & 1) == 1) res *= static_cast<Result>(0.5);
-        return res; 
+        return res;
       }
     };
@@ -368 +368 @@
     struct ShiftMultiplier<Result, 0, true> {
       static const Result multiplier() {
-        return static_cast<Result>(1.0); 
+        return static_cast<Result>(1.0);
       }
     };
@@ -375 +375 @@
     struct ShiftMultiplier<Result, -20, true> {
       static const Result multiplier() {
-        return static_cast<Result>(1.0/1048576.0); 
-      }
-    };
-    
+        return static_cast<Result>(1.0/1048576.0);
+      }
+    };
+
     template <typename Result>
     struct ShiftMultiplier<Result, -32, true> {
       static const Result multiplier() {
-        return static_cast<Result>(1.0/424967296.0); 
+        return static_cast<Result>(1.0/424967296.0);
       }
     };
@@ -389 +389 @@
     struct ShiftMultiplier<Result, -53, true> {
       static const Result multiplier() {
-        return static_cast<Result>(1.0/9007199254740992.0); 
+        return static_cast<Result>(1.0/9007199254740992.0);
       }
     };
@@ -396 +396 @@
     struct ShiftMultiplier<Result, -64, true> {
       static const Result multiplier() {
-        return static_cast<Result>(1.0/18446744073709551616.0); 
+        return static_cast<Result>(1.0/18446744073709551616.0);
       }
     };
@@ -408 +408 @@
 
     template <typename Result, typename Word,
-              int rest = std::numeric_limits<Result>::digits, int shift = 0, 
+              int rest = std::numeric_limits<Result>::digits, int shift = 0,
               bool last = rest <= std::numeric_limits<Word>::digits>
-    struct RealConversion{ 
+    struct RealConversion{
       static const int bits = std::numeric_limits<Word>::digits;
 
@@ -420 +420 @@
 
     template <typename Result, typename Word, int rest, int shift>
-    struct RealConversion<Result, Word, rest, shift, false> { 
+    struct RealConversion<Result, Word, rest, shift, false> {
       static const int bits = std::numeric_limits<Word>::digits;
 
@@ -459 +459 @@
       Word buffer;
       int num;
-      
+
       BoolProducer() : num(0) {}
 
@@ -530 +530 @@
     // Architecture word
     typedef unsigned long Word;
-    
+
     _random_bits::RandomCore<Word> core;
     _random_bits::BoolProducer<Word> bool_producer;
-    
+
 
   public:
@@ -555 +555 @@
     /// to the architecture word type.
     template <typename Number>
-    Random(Number seed) { 
+    Random(Number seed) {
       _random_bits::Initializer<Number, Word>::init(core, seed);
     }
@@ -565 +565 @@
     /// architecture word type.
     template <typename Iterator>
-    Random(Iterator begin, Iterator end) { 
+    Random(Iterator begin, Iterator end) {
       typedef typename std::iterator_traits<Iterator>::value_type Number;
       _random_bits::Initializer<Number, Word>::init(core, begin, end);
@@ -598 +598 @@
     /// converted to the architecture word type.
     template <typename Number>
-    void seed(Number seed) { 
+    void seed(Number seed) {
       _random_bits::Initializer<Number, Word>::init(core, seed);
     }
@@ -608 +608 @@
     /// architecture word type.
     template <typename Iterator>
-    void seed(Iterator begin, Iterator end) { 
+    void seed(Iterator begin, Iterator end) {
       typedef typename std::iterator_traits<Iterator>::value_type Number;
       _random_bits::Initializer<Number, Word>::init(core, begin, end);
@@ -626 +626 @@
       return false;
     }
-    
+
     /// \brief Seeding from file
     ///
@@ -641 +641 @@
     /// \return True when the seeding successes.
 #ifndef WIN32
-    bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0) 
+    bool seedFromFile(const std::string& file = "/dev/urandom", int offset = 0)
 #else
-    bool seedFromFile(const std::string& file = "", int offset = 0) 
+    bool seedFromFile(const std::string& file = "", int offset = 0)
 #endif
     {
@@ -661 +661 @@
     /// random sequence.
     /// \return Currently always true.
-    bool seedFromTime() {       
+    bool seedFromTime() {
 #ifndef WIN32
       timeval tv;
@@ -697 +697 @@
     /// It returns a random real number from the range [0, b).
     template <typename Number>
-    Number real(Number b) { 
-      return real<Number>() * b; 
+    Number real(Number b) {
+      return real<Number>() * b;
     }
 
@@ -705 +705 @@
     /// It returns a random real number from the range [a, b).
     template <typename Number>
-    Number real(Number a, Number b) { 
-      return real<Number>() * (b - a) + a; 
+    Number real(Number a, Number b) {
+      return real<Number>() * (b - a) + a;
     }
 
@@ -726 +726 @@
     /// It returns a random real number from the range [0, b).
     template <typename Number>
-    Number operator()(Number b) { 
-      return real<Number>() * b; 
+    Number operator()(Number b) {
+      return real<Number>() * b;
     }
 
@@ -734 +734 @@
     /// It returns a random real number from the range [a, b).
     template <typename Number>
-    Number operator()(Number a, Number b) { 
-      return real<Number>() * (b - a) + a; 
+    Number operator()(Number a, Number b) {
+      return real<Number>() * (b - a) + a;
     }
 
@@ -785 +785 @@
     template <typename Number>
     Number integer() {
-      static const int nb = std::numeric_limits<Number>::digits + 
+      static const int nb = std::numeric_limits<Number>::digits +
         (std::numeric_limits<Number>::is_signed ? 1 : 0);
       return _random_bits::IntConversion<Number, Word, nb>::convert(core);
@@ -793 +793 @@
       return integer<int>();
     }
-    
+
     /// \brief Returns a random bool
     ///
@@ -807 +807 @@
     ///\name Non-uniform distributions
     ///
-    
+
     ///@{
-    
+
     /// \brief Returns a random bool
     ///
@@ -823 +823 @@
     /// transformation is used to generate a random normal distribution.
     /// \todo Consider using the "ziggurat" method instead.
-    double gauss() 
+    double gauss()
     {
       double V1,V2,S;
       do {
-        V1=2*real<double>()-1;
-        V2=2*real<double>()-1;
-        S=V1*V1+V2*V2;
+        V1=2*real<double>()-1;
+        V2=2*real<double>()-1;
+        S=V1*V1+V2*V2;
       } while(S>=1);
       return std::sqrt(-2*std::log(S)/S)*V1;
@@ -855 +855 @@
 
     /// This function generates a gamma distribution random number.
-    /// 
+    ///
     ///\param k shape parameter (<tt>k>0</tt> integer)
-    double gamma(int k) 
+    double gamma(int k)
     {
       double s = 0;
@@ -863 +863 @@
       return s;
     }
-    
+
     /// Gamma distribution with given shape and scale parameter
 
     /// This function generates a gamma distribution random number.
-    /// 
+    ///
     ///\param k shape parameter (<tt>k>0</tt>)
     ///\param theta scale parameter
@@ -877 +877 @@
       const double v0=E/(E-delta);
       do {
-        double V0=1.0-real<double>();
-        double V1=1.0-real<double>();
-        double V2=1.0-real<double>();
-        if(V2<=v0) 
-          {
-            xi=std::pow(V1,1.0/delta);
-            nu=V0*std::pow(xi,delta-1.0);
-          }
-        else 
-          {
-            xi=1.0-std::log(V1);
-            nu=V0*std::exp(-xi);
-          }
+        double V0=1.0-real<double>();
+        double V1=1.0-real<double>();
+        double V2=1.0-real<double>();
+        if(V2<=v0)
+          {
+            xi=std::pow(V1,1.0/delta);
+            nu=V0*std::pow(xi,delta-1.0);
+          }
+        else
+          {
+            xi=1.0-std::log(V1);
+            nu=V0*std::exp(-xi);
+          }
       } while(nu>std::pow(xi,delta-1.0)*std::exp(-xi));
       return theta*(xi+gamma(int(std::floor(k))));
     }
-    
+
     /// Weibull distribution
 
     /// This function generates a Weibull distribution random number.
-    /// 
+    ///
     ///\param k shape parameter (<tt>k>0</tt>)
     ///\param lambda scale parameter (<tt>lambda>0</tt>)
@@ -904 +904 @@
     {
       return lambda*pow(-std::log(1.0-real<double>()),1.0/k);
-    }  
-      
+    }
+
     /// Pareto distribution
 
     /// This function generates a Pareto distribution random number.
-    /// 
+    ///
     ///\param k shape parameter (<tt>k>0</tt>)
     ///\param x_min location parameter (<tt>x_min>0</tt>)
@@ -916 +916 @@
     {
       return exponential(gamma(k,1.0/x_min))+x_min;
-    }  
-      
+    }
+
     /// Poisson distribution
 
     /// This function generates a Poisson distribution random number with
     /// parameter \c lambda.
-    /// 
+    ///
     /// The probability mass function of this distribusion is
     /// \f[ \frac{e^{-\lambda}\lambda^k}{k!} \f]
@@ -928 +928 @@
     /// ''Seminumerical Algorithms'' (1969). Its running time is linear in the
     /// return value.
-    
+
     int poisson(double lambda)
     {
@@ -935 +935 @@
       double p = 1.0;
       do {
-        k++;
-        p*=real<double>();
+        k++;
+        p*=real<double>();
       } while (p>=l);
       return k-1;
-    }  
-      
+    }
+
     ///@}
-    
+
     ///\name Two dimensional distributions
     ///
 
     ///@{
-    
+
     /// Uniform distribution on the full unit circle
 
     /// Uniform distribution on the full unit circle.
     ///
-    dim2::Point<double> disc() 
+    dim2::Point<double> disc()
     {
       double V1,V2;
       do {
-        V1=2*real<double>()-1;
-        V2=2*real<double>()-1;
-        
+        V1=2*real<double>()-1;
+        V2=2*real<double>()-1;
+
       } while(V1*V1+V2*V2>=1);
       return dim2::Point<double>(V1,V2);
@@ -974 +974 @@
       double V1,V2,S;
       do {
-        V1=2*real<double>()-1;
-        V2=2*real<double>()-1;
-        S=V1*V1+V2*V2;
+        V1=2*real<double>()-1;
+        V2=2*real<double>()-1;
+        S=V1*V1+V2*V2;
       } while(S>=1);
       double W=std::sqrt(-2*std::log(S)/S);
@@ -985 +985 @@
     /// This function provides a turning symmetric two-dimensional distribution.
     /// The x-coordinate is of conditionally exponential distribution
-    /// with the condition that x is positive and y=0. If x is negative and 
+    /// with the condition that x is positive and y=0. If x is negative and
     /// y=0 then, -x is of exponential distribution. The same is true for the
     /// y-coordinate.
-    dim2::Point<double> exponential2() 
+    dim2::Point<double> exponential2()
     {
       double V1,V2,S;
       do {
-        V1=2*real<double>()-1;
-        V2=2*real<double>()-1;
-        S=V1*V1+V2*V2;
+        V1=2*real<double>()-1;
+        V2=2*real<double>()-1;
+        S=V1*V1+V2*V2;
       } while(S>=1);
       double W=-std::log(S)/S;
@@ -1000 +1000 @@
     }
 
-    ///@}    
+    ///@}
   };