diff --git a/include/caffe/util/math_functions.hpp b/include/caffe/util/math_functions.hpp
index 05ed1301599..a4c63de0440 100644
--- a/include/caffe/util/math_functions.hpp
+++ b/include/caffe/util/math_functions.hpp
@@ -85,6 +85,9 @@ void caffe_div(const int N, const Dtype* a, const Dtype* b, Dtype* y);
 template <typename Dtype>
 void caffe_powx(const int n, const Dtype* a, const Dtype b, Dtype* y);
 
+template <typename Dtype>
+Dtype caffe_nextafter(const Dtype b);
+
 template <typename Dtype>
 void caffe_vRngUniform(const int n, Dtype* r, const Dtype a, const Dtype b);
 
diff --git a/src/caffe/test/test_multinomial_logistic_loss_layer.cpp b/src/caffe/test/test_multinomial_logistic_loss_layer.cpp
index adb36627606..6bd94ae24b8 100644
--- a/src/caffe/test/test_multinomial_logistic_loss_layer.cpp
+++ b/src/caffe/test/test_multinomial_logistic_loss_layer.cpp
@@ -24,6 +24,7 @@ class MultinomialLogisticLossLayerTest : public ::testing::Test {
   MultinomialLogisticLossLayerTest()
       : blob_bottom_data_(new Blob<Dtype>(10, 5, 1, 1)),
         blob_bottom_label_(new Blob<Dtype>(10, 1, 1, 1)) {
+    Caffe::set_random_seed(1701);
     // fill the values
     FillerParameter filler_param;
     PositiveUnitballFiller<Dtype> filler(filler_param);
diff --git a/src/caffe/test/test_random_number_generator.cpp b/src/caffe/test/test_random_number_generator.cpp
new file mode 100644
index 00000000000..4c3358f9f49
--- /dev/null
+++ b/src/caffe/test/test_random_number_generator.cpp
@@ -0,0 +1,67 @@
+#include <cmath>
+#include <cstring>
+#include <cuda_runtime.h>
+
+#include "gtest/gtest.h"
+#include "caffe/common.hpp"
+#include "caffe/syncedmem.hpp"
+#include "caffe/util/math_functions.hpp"
+#include "caffe/test/test_caffe_main.hpp"
+
+namespace caffe {
+
+template <typename Dtype>
+class RandomNumberGeneratorTest : public ::testing::Test {
+ public:
+  virtual ~RandomNumberGeneratorTest() {}
+
+  Dtype sample_mean(const Dtype* const seqs, const size_t sample_size)
+  {
+      double sum = 0;
+      for (int i = 0; i < sample_size; ++i) {
+          sum += seqs[i];
+      }
+      return sum / sample_size;
+  }
+
+  Dtype mean_bound(const Dtype std, const size_t sample_size)
+  {
+      return  std/sqrt((double)sample_size);
+  }
+};
+
+
+typedef ::testing::Types<float, double> Dtypes;
+TYPED_TEST_CASE(RandomNumberGeneratorTest, Dtypes);
+
+TYPED_TEST(RandomNumberGeneratorTest, TestRngGaussian) {
+  size_t sample_size = 10000;
+  SyncedMemory data_a(sample_size * sizeof(TypeParam));
+  Caffe::set_random_seed(1701);
+  TypeParam mu = 0;
+  TypeParam sigma = 1;
+  caffe_vRngGaussian(sample_size, (TypeParam*)data_a.mutable_cpu_data(), mu, sigma);
+  TypeParam true_mean = mu;
+  TypeParam true_std = sigma;
+  TypeParam bound = mean_bound(true_std, sample_size);
+  TypeParam real_mean = sample_mean((TypeParam*)data_a.cpu_data(), sample_size);
+  EXPECT_NEAR(real_mean, true_mean, bound);
+}
+
+TYPED_TEST(RandomNumberGeneratorTest, TestRngUniform) {
+  size_t sample_size = 10000;
+  SyncedMemory data_a(sample_size * sizeof(TypeParam));
+  Caffe::set_random_seed(1701);
+  TypeParam lower = 0;
+  TypeParam upper = 1;
+  caffe_vRngUniform(sample_size, (TypeParam*)data_a.mutable_cpu_data(), lower, upper);
+  TypeParam true_mean = (lower + upper) / 2;
+  TypeParam true_std = (upper - lower) / sqrt(12);
+  TypeParam bound = mean_bound(true_std, sample_size);
+  TypeParam real_mean = sample_mean((TypeParam*)data_a.cpu_data(), sample_size);
+  EXPECT_NEAR(real_mean, true_mean, bound);
+}
+
+
+
+}  // namespace caffe
diff --git a/src/caffe/util/math_functions.cpp b/src/caffe/util/math_functions.cpp
index 25de4251a00..852cf6d6d38 100644
--- a/src/caffe/util/math_functions.cpp
+++ b/src/caffe/util/math_functions.cpp
@@ -1,7 +1,9 @@
 // Copyright 2013 Yangqing Jia
 
+#include <limits>
 //#include <mkl.h>
 #include <eigen3/Eigen/Dense>
+#include <boost/math/special_functions/next.hpp>
 #include <boost/random.hpp>
 
 #include <cublas_v2.h>
@@ -280,6 +282,11 @@ void caffe_powx<double>(const int n, const double* a, const double b,
     map_vector_double_t(y, n) = const_map_vector_double_t(a, n).array().pow(b);
 }
 
+template <typename Dtype>
+Dtype caffe_nextafter(const Dtype b) {
+  return boost::math::nextafter<Dtype, Dtype>(b, std::numeric_limits<Dtype>::max());
+}
+
 template <>
 void caffe_vRngUniform<float>(const int n, float* r,
     const float a, const float b) {
@@ -287,7 +294,8 @@ void caffe_vRngUniform<float>(const int n, float* r,
   //    n, r, a, b));
 
   // FIXME check if boundaries are handled in the same way ?
-  boost::uniform_real<float> random_distribution(a, b);
+  boost::random::uniform_real_distribution<float> random_distribution(
+      a, caffe_nextafter<float>(b));
   Caffe::random_generator_t &generator = Caffe::vsl_stream();
 
   for(int i = 0; i < n; i += 1)
@@ -303,7 +311,8 @@ void caffe_vRngUniform<double>(const int n, double* r,
   //    n, r, a, b));
 
     // FIXME check if boundaries are handled in the same way ?
-    boost::uniform_real<double> random_distribution(a, b);
+    boost::random::uniform_real_distribution<double> random_distribution(
+        a, caffe_nextafter<double>(b));
     Caffe::random_generator_t &generator = Caffe::vsl_stream();
 
     for(int i = 0; i < n; i += 1)
@@ -315,6 +324,7 @@ void caffe_vRngUniform<double>(const int n, double* r,
 template <>
 void caffe_vRngGaussian<float>(const int n, float* r, const float a,
     const float sigma) {
+    DCHECK(sigma > 0);
   //VSL_CHECK(vsRngGaussian(VSL_RNG_METHOD_GAUSSIAN_BOXMULLER,
 //      Caffe::vsl_stream(), n, r, a, sigma));
 
@@ -332,6 +342,7 @@ void caffe_vRngGaussian<float>(const int n, float* r, const float a,
 template <>
 void caffe_vRngGaussian<double>(const int n, double* r, const double a,
     const double sigma) {
+    DCHECK(sigma > 0);
   //VSL_CHECK(vdRngGaussian(VSL_RNG_METHOD_GAUSSIAN_BOXMULLER,
   //    Caffe::vsl_stream(), n, r, a, sigma));