src/core/random-variable.cc
changeset 3292 bd8d2601af21
parent 3267 d0c70dbe918e
child 3325 11719a408a0c
--- a/src/core/random-variable.cc	Wed Jun 18 21:29:21 2008 -0700
+++ b/src/core/random-variable.cc	Thu Jun 19 14:46:27 2008 -0700
@@ -1033,7 +1033,8 @@
     }
   while(1)
     { // See Simulation Modeling and Analysis p. 466 (Averill Law)
-      // for algorithm
+      // for algorithm; basically a Box-Muller transform:
+      // http://en.wikipedia.org/wiki/Box-Muller_transform
       double u1 = m_generator->RandU01();
       double u2 = m_generator->RandU01();;
       double v1 = 2 * u1 - 1;
@@ -1082,7 +1083,8 @@
     }
   while(1)
     { // See Simulation Modeling and Analysis p. 466 (Averill Law)
-      // for algorithm
+      // for algorithm; basically a Box-Muller transform:
+      // http://en.wikipedia.org/wiki/Box-Muller_transform
       double u1 = m_static_generator->RandU01();
       double u2 = m_static_generator->RandU01();;
       double v1 = 2 * u1 - 1;
@@ -1092,11 +1094,21 @@
         { // Got good pair
           double y = sqrt((-2 * log(w))/w);
           m_static_next = m + v2 * y * sqrt(v);
-          if (fabs(m_static_next) > b) m_static_next = b * (m_static_next)/fabs(m_static_next);
-          m_static_nextValid = true;
+          //if next is in bounds, it is valid
+          m_static_nextValid = fabs(m_static_next-m) <= b;;
           double x1 = m + v1 * y * sqrt(v);
-          if (fabs(x1) > b) x1 = b * (x1)/fabs(x1);
-          return x1;
+          //if x1 is in bounds, return it
+          if (fabs(x1-m) <= b)
+          {
+            return x1;
+          }
+          //otherwise try and return m_next if it is valid
+          else if (m_static_nextValid)
+          {
+            m_static_nextValid = false;
+            return m_static_next;
+          }
+          //otherwise, just run this loop again
         }
     }
 }