diff -r 86a97665dcb0 src/core/random-variable.cc
--- a/src/core/random-variable.cc	Tue Jan 20 16:36:08 2009 -0500
+++ b/src/core/random-variable.cc	Thu Jan 22 12:21:57 2009 -0500
@@ -577,9 +577,13 @@
     m_generator->InitializeStream();
     m_generator->ResetNthSubstream(SeedManager::GetRun());
   }
-  double r = -m_mean*log(m_generator->RandU01());
-  if (m_bound != 0 && r > m_bound) return m_bound;
-  return r;
+
+  while(1)
+  {
+      double r = -m_mean*log(m_generator->RandU01());
+      if (m_bound == 0 || r <= m_bound) return r;
+      //otherwise, try again
+  }
 }
 
 RandomVariableBase* ExponentialVariableImpl::Copy() const
@@ -679,10 +683,14 @@
     m_generator->InitializeStream();
     m_generator->ResetNthSubstream(SeedManager::GetRun());
   }
+
   double scale = m_mean * ( m_shape - 1.0) / m_shape;
-  double r = (scale * ( 1.0 / pow(m_generator->RandU01(), 1.0 / m_shape)));
-  if (m_bound != 0 && r > m_bound) return m_bound;
-  return r;
+  while(1)
+  {
+      double r = (scale * ( 1.0 / pow(m_generator->RandU01(), 1.0 / m_shape)));
+      if (m_bound == 0 || r <= m_bound) return r;
+      //otherwise, try again
+  }
 }
 
 RandomVariableBase* ParetoVariableImpl::Copy() const
@@ -781,10 +789,14 @@
     m_generator->InitializeStream();
     m_generator->ResetNthSubstream(SeedManager::GetRun());
   }
+
   double exponent = 1.0 / m_alpha;
-  double r = m_mean * pow( -log(m_generator->RandU01()), exponent);
-  if (m_bound != 0 && r > m_bound) return m_bound;
-  return r;
+  while(1)
+  {
+      double r = m_mean * pow( -log(m_generator->RandU01()), exponent);
+      if (m_bound == 0 || r <= m_bound) return r;
+      //otherwise, try again
+  }
 }
 
 RandomVariableBase* WeibullVariableImpl::Copy() const
@@ -823,7 +835,7 @@
    * \brief Construct a normal random variable with specified mean and variance
    * \param m Mean value
    * \param v Variance
-   * \param b Bound.  The NormalVariableImpl is bounded within +-bound.
+   * \param b Bound.  The NormalVariableImpl is bounded within +-bound of the mean.
    */ 
   NormalVariableImpl(double m, double v, double b = INFINITE_VALUE);
 
@@ -838,7 +850,7 @@
 private:
   double m_mean;      // Mean value of RV
   double m_variance;  // Mean value of RV
-  double m_bound;     // Bound on value (absolute value)
+  double m_bound;     // Bound on value's difference from the mean (absolute value)
   bool   m_nextValid; // True if next valid
   double m_next;      // The algorithm produces two values at a time
   static bool   m_static_nextValid;
@@ -889,11 +901,21 @@
         { // Got good pair
           double y = sqrt((-2 * log(w))/w);
           m_next = m_mean + v2 * y * sqrt(m_variance);
-          if (fabs(m_next) > m_bound) m_next = m_bound * (m_next)/fabs(m_next);
-          m_nextValid = true;
+          //if next is in bounds, it is valid
+          m_nextValid = fabs(m_next-m_mean) <= m_bound;
           double x1 = m_mean + v1 * y * sqrt(m_variance);
-          if (fabs(x1) > m_bound) x1 = m_bound * (x1)/fabs(x1);
-          return x1;
+          //if x1 is in bounds, return it
+          if (fabs(x1-m_mean) <= m_bound)
+            {
+              return x1;
+            }
+          //otherwise try and return m_next if it is valid
+          else if (m_nextValid)
+            {
+              m_nextValid = false;
+              return m_next;
+            }
+          //otherwise, just run this loop again
         }
     }
 }