diff --git a/dynet/CMakeLists.txt b/dynet/CMakeLists.txt
index 6b15280e1..22f9b4471 100644
--- a/dynet/CMakeLists.txt
+++ b/dynet/CMakeLists.txt
@@ -55,6 +55,7 @@ set(dynet_library_SRCS
     param-init.cc
     param-nodes.cc
     pretrain.cc
+    rand.cc
     rnn-state-machine.cc
     rnn.cc
     saxe-init.cc
@@ -135,6 +136,7 @@ set(dynet_library_HDRS
     param-init.h
     param-nodes.h
     pretrain.h
+    rand.h
     rnn-state-machine.h
     rnn.h
     saxe-init.h
diff --git a/dynet/cfsm-builder.cc b/dynet/cfsm-builder.cc
index 3032ba4f5..4acb73e09 100644
--- a/dynet/cfsm-builder.cc
+++ b/dynet/cfsm-builder.cc
@@ -1,6 +1,7 @@
 #include "dynet/cfsm-builder.h"
 #include "dynet/except.h"
 #include "dynet/param-init.h"
+#include "dynet/rand.h"
 
 #include <fstream>
 #include <iostream>
diff --git a/dynet/hsm-builder.cc b/dynet/hsm-builder.cc
index e3c8a128c..d65e61492 100644
--- a/dynet/hsm-builder.cc
+++ b/dynet/hsm-builder.cc
@@ -5,6 +5,7 @@
 #include <sstream>
 
 #include "dynet/param-init.h"
+#include "dynet/rand.h"
 
 using namespace std;
 
diff --git a/dynet/nodes-dropout.cc b/dynet/nodes-dropout.cc
index 1d052f18b..a84a69920 100644
--- a/dynet/nodes-dropout.cc
+++ b/dynet/nodes-dropout.cc
@@ -1,3 +1,4 @@
+#include "dynet/rand.h"
 #include "dynet/nodes-dropout.h"
 
 #include "dynet/nodes-macros.h"
@@ -30,7 +31,9 @@ size_t Dropout::aux_storage_size() const {
 template<class MyDevice>
 void Dropout::forward_dev_impl(const MyDevice & dev, const vector<const Tensor*>& xs, Tensor& fx) const {
   Tensor m(dim, (float*)aux_mem, fx.device, DeviceMempool::FXS);
-  TensorTools::randomize_bernoulli(m, (1.f-p), 1.f / (1.f-p));
+  Eigen::internal::UniformRandomGenerator<float> uni_rg(draw_random_seed());
+  m.tvec().device(*dev.edevice) = m.tvec().random(uni_rg);
+  m.tvec().device(*dev.edevice) = (m.tvec() < m.tvec().constant((1.f-p))).cast<float>() * 1.f / (1.f-p);
   fx.tvec().device(*dev.edevice) = xs[0]->tvec() * m.tvec();
 }
 
@@ -74,7 +77,10 @@ void DropoutDim::forward_dev_impl(const MyDevice & dev, const vector<const Tenso
   Dim mask_dim(dim);
   mask_dim.d[dimension]=1;
   Tensor m(mask_dim, (float*)aux_mem, fx.device, DeviceMempool::FXS);
-  TensorTools::randomize_bernoulli(m, (1.f-p), 1.f / (1.f-p));
+  Eigen::internal::UniformRandomGenerator<float> uni_rg(draw_random_seed());
+  m.tvec().device(*dev.edevice) = m.tvec().random(uni_rg);
+  m.tvec().device(*dev.edevice) = (m.tvec() < m.tvec().constant((1.f-p))).cast<float>() / (1.f-p);
+
   Eigen::array<ptrdiff_t, 4> bcast = {1, 1, 1, 1}; bcast[dimension] = xs[0]->d[dimension];
   fx.tb<3>().device(*dev.edevice) = xs[0]->tb<3>() * m.tb<3>().broadcast(bcast);
 }
@@ -119,7 +125,9 @@ template<class MyDevice>
 void DropoutBatch::forward_dev_impl(const MyDevice & dev, const vector<const Tensor*>& xs, Tensor& fx) const {
   Dim mask_dim({1},xs[0]->d.batch_elems());
   Tensor m(mask_dim, (float*)aux_mem, fx.device, DeviceMempool::FXS);
-  TensorTools::randomize_bernoulli(m, (1.f-p), 1.f / (1.f-p));
+  Eigen::internal::UniformRandomGenerator<float> uni_rg(draw_random_seed());
+  m.tvec().device(*dev.edevice) = m.tvec().random(uni_rg);
+  m.tvec().device(*dev.edevice) = (m.tvec() < m.tvec().constant((1.f-p))).cast<float>() * 1.f / (1.f-p);
   Eigen::array<ptrdiff_t, 2> bcast = {xs[0]->d.batch_size(), 1};
   fx.tbvec().device(*dev.edevice) = xs[0]->tbvec() * m.tbvec().broadcast(bcast);
 }
diff --git a/dynet/nodes-lstm.cc b/dynet/nodes-lstm.cc
index b1a7adcfe..9b716dc37 100644
--- a/dynet/nodes-lstm.cc
+++ b/dynet/nodes-lstm.cc
@@ -1,3 +1,4 @@
+#include "dynet/rand.h"
 #include "dynet/nodes-lstm.h"
 #include "dynet/matrix-multiply.h"
 
@@ -200,17 +201,20 @@ namespace dynet {
     if(weightnoise_std > 0.f){
       Tensor Wx_noisy(Dim({hidden_dim*4, input_dim},1), nullptr, fx.device, fx.mem_pool);
       Wx_noisy.v = static_cast<float*>(scratch_allocator->allocate(Wx_noisy.d.size() * sizeof(float)));
-      TensorTools::randomize_normal(Wx_noisy, 0, weightnoise_std);
+      Eigen::internal::NormalRandomGenerator<float> normal_rg1(draw_random_seed());
+      Wx_noisy.tvec().device(*dev.edevice) = Wx_noisy.tvec().random(normal_rg1) * weightnoise_std;
       Wx_noisy.tvec().device(*dev.edevice) += Wx->tvec();
 
       Tensor Wh_noisy(Dim({hidden_dim*4, hidden_dim},1), nullptr, fx.device, fx.mem_pool);
       Wh_noisy.v = static_cast<float*>(scratch_allocator->allocate(Wh_noisy.d.size() * sizeof(float)));
-      TensorTools::randomize_normal(Wh_noisy, 0, weightnoise_std);
+      Eigen::internal::NormalRandomGenerator<float> normal_rg2(draw_random_seed());
+      Wh_noisy.tvec().device(*dev.edevice) = Wh_noisy.tvec().random(normal_rg2) * weightnoise_std;
       Wh_noisy.tvec().device(*dev.edevice) += Wh->tvec();
 
       Tensor b_noisy(Dim({hidden_dim*4, 1},1), nullptr, fx.device, fx.mem_pool);
       b_noisy.v = static_cast<float*>(scratch_allocator->allocate(b_noisy.d.size() * sizeof(float)));
-      TensorTools::randomize_normal(b_noisy, 0, weightnoise_std);
+      Eigen::internal::NormalRandomGenerator<float> normal_rg3(draw_random_seed());
+      b_noisy.tvec().device(*dev.edevice) = b_noisy.tvec().random(normal_rg3) * weightnoise_std;
       b_noisy.tvec().device(*dev.edevice) += b->tvec();
 
     } else {
diff --git a/dynet/nodes-random.cc b/dynet/nodes-random.cc
index 1324a66ab..8b62186ba 100644
--- a/dynet/nodes-random.cc
+++ b/dynet/nodes-random.cc
@@ -1,3 +1,4 @@
+#include "dynet/rand.h"
 #include "dynet/nodes-random.h"
 
 #include "dynet/nodes-macros.h"
@@ -25,15 +26,8 @@ Dim GaussianNoise::dim_forward(const vector<Dim>& xs) const {
 
 template<class MyDevice>
 void GaussianNoise::forward_dev_impl(const MyDevice & dev, const vector<const Tensor*>& xs, Tensor& fx) const {
-
-  AlignedMemoryPool* scratch_allocator = fx.device->pools[(int)DeviceMempool::SCS];
-  Tensor noise(dim, nullptr, fx.device, fx.mem_pool);
-  noise.v = static_cast<float*>(scratch_allocator->allocate(noise.d.size() * sizeof(float)));
-  TensorTools::randomize_normal(noise, 0, stddev);
-
-  fx.tvec().device(*dev.edevice) = xs[0]->tvec() + noise.tvec();
-
-  scratch_allocator->free();
+  Eigen::internal::NormalRandomGenerator<float> normal_rg(draw_random_seed());
+  fx.tvec().device(*dev.edevice) = xs[0]->tvec() + xs[0]->tvec().random(normal_rg) * stddev;
 }
 
 template<class MyDevice>
@@ -66,7 +60,8 @@ Dim RandomNormal::dim_forward(const vector<Dim>& xs) const {
 template<class MyDevice>
 void RandomNormal::forward_dev_impl(const MyDevice & dev, const vector<const Tensor*>& xs, Tensor& fx) const {
   DYNET_ASSERT(xs.size() == 0, "Failed dimension check in RandomNormal::forward");
-  TensorTools::randomize_normal(fx);
+  Eigen::internal::NormalRandomGenerator<float> normal_rg(draw_random_seed());
+  fx.tvec().device(*dev.edevice) = fx.tvec().random(normal_rg);
 }
 
 template<class MyDevice>
@@ -99,7 +94,9 @@ Dim RandomBernoulli::dim_forward(const vector<Dim>& xs) const {
 template<class MyDevice>
 void RandomBernoulli::forward_dev_impl(const MyDevice & dev, const vector<const Tensor*>& xs, Tensor& fx) const {
   DYNET_ASSERT(xs.size() == 0, "Failed dimension check in RandomBernoulli::forward");
-  TensorTools::randomize_bernoulli(fx, p, scale);
+  Eigen::internal::UniformRandomGenerator<float> uni_rg(draw_random_seed());
+  fx.tvec().device(*dev.edevice) = fx.tvec().random(uni_rg);
+  fx.tvec().device(*dev.edevice) = (fx.tvec() < fx.tvec().constant(p)).cast<float>() * scale;
 }
 
 template<class MyDevice>
@@ -132,7 +129,8 @@ Dim RandomUniform::dim_forward(const vector<Dim>& xs) const {
 template<class MyDevice>
 void RandomUniform::forward_dev_impl(const MyDevice & dev, const vector<const Tensor*>& xs, Tensor& fx) const {
   DYNET_ASSERT(xs.size() == 0, "Failed dimension check in RandomUniform::forward");
-  TensorTools::randomize_uniform(fx, left, right);
+  Eigen::internal::UniformRandomGenerator<float> uni_rg(draw_random_seed());
+  fx.tvec().device(*dev.edevice) = (fx.tvec().random(uni_rg) * (right-left)) + left;
 }
 
 template<class MyDevice>
@@ -166,7 +164,8 @@ template<class MyDevice>
 void RandomGumbel::forward_dev_impl(const MyDevice & dev, const vector<const Tensor*>& xs, Tensor& fx) const {
   DYNET_ASSERT(xs.size() == 0, "Failed dimension check in RandomGumbel::forward");
   DYNET_ARG_CHECK(mu == 0.0 && beta == 1.0, "RandomGumbel only supports Gumbel(0,1) at the moment (pull requests welcome)");
-  TensorTools::randomize_uniform(fx, 0, 1);
+  Eigen::internal::UniformRandomGenerator<float> uni_rg(draw_random_seed());
+  fx.tvec().device(*dev.edevice) = fx.tvec().random(uni_rg);
   float eps = 1e-20;
   fx.tvec().device(*dev.edevice) = -(-fx.tvec().cwiseMax(eps).log()).cwiseMax(eps).log();
 }
diff --git a/dynet/rand.cc b/dynet/rand.cc
new file mode 100644
index 000000000..eee13ae8f
--- /dev/null
+++ b/dynet/rand.cc
@@ -0,0 +1,31 @@
+#include "dynet/rand.h"
+//#include "dynet/globals.h"
+#include <random>
+
+using namespace std;
+
+namespace dynet {
+
+  real rand01() {
+    uniform_real_distribution<real> distribution(0, 1);
+    return distribution(*rndeng);
+  }
+
+  int rand0n(int n) {
+    if (n <= 0) throw std::runtime_error("Integer upper bound is non-positive");
+    int x = rand01() * n;
+    while (n == x) { x = rand01() * n; }
+    return x;
+  }
+
+  real rand_normal() {
+    normal_distribution<real> distribution(0, 1);
+    return distribution(*rndeng);
+  }
+
+  int draw_random_seed(){
+    std::uniform_int_distribution<> seed_dist(1, 2147483647);
+    return seed_dist(*rndeng);
+  }
+} // namespace dynet
+
diff --git a/dynet/rand.h b/dynet/rand.h
new file mode 100644
index 000000000..ec86e5171
--- /dev/null
+++ b/dynet/rand.h
@@ -0,0 +1,37 @@
+#ifndef DYNET_RAND_H
+#define DYNET_RAND_H
+
+#include "dynet/tensor.h"
+
+namespace dynet {
+
+  /**
+   * \ingroup random
+   * \brief This is a helper function to sample uniformly in \f$[0,1]\f$
+   * \return \f$x\sim\mathcal U([0,1])\f$
+   */
+  real rand01();
+  /**
+   * \ingroup random
+   * \brief This is a helper function to sample uniformly in \f$\{0,\dots,n-1\}\f$
+   *
+   * \param n Upper bound (excluded)
+   * \return \f$x\sim\mathcal U(\{0,\dots,n-1\})\f$
+   */
+  int rand0n(int n);
+  /**
+   * \ingroup random
+   * \brief This is a helper function to sample from a normalized gaussian distribution
+   *
+   * \return \f$x\sim\mathcal N(0,1)\f$
+   */
+  real rand_normal();
+  /**
+   * \ingroup random
+   * \brief This returns a new random seed.
+   */
+  int draw_random_seed();
+
+} // namespace dynet
+
+#endif
diff --git a/dynet/tensor.cc b/dynet/tensor.cc
index 099c32ab7..f2353d39b 100644
--- a/dynet/tensor.cc
+++ b/dynet/tensor.cc
@@ -285,22 +285,6 @@ void TensorTools::randomize_orthonormal(Tensor& val, real scale) {
   } else { throw std::runtime_error("Bad device type"); }
 }
 
-real rand01() {
-  uniform_real_distribution<real> distribution(0, 1);
-  return distribution(*rndeng);
-}
-
-int rand0n(int n) {
-  if (n <= 0) throw std::runtime_error("Integer upper bound is non-positive");
-  int x = rand01() * n;
-  while (n == x) { x = rand01() * n; }
-  return x;
-}
-
-real rand_normal() {
-  normal_distribution<real> distribution(0, 1);
-  return distribution(*rndeng);
-}
 
 #endif
 
diff --git a/dynet/tensor.h b/dynet/tensor.h
index c9ba9fcd2..5d3bf41a0 100644
--- a/dynet/tensor.h
+++ b/dynet/tensor.h
@@ -729,28 +729,6 @@ struct TensorTools {
 
 };
 
-/**
- * \ingroup tensor
- * \brief This is a helper function to sample uniformly in \f$[0,1]\f$
- * \return \f$x\sim\mathcal U([0,1])\f$
- */
-real rand01();
-/**
- * \ingroup tensor
- * \brief This is a helper function to sample uniformly in \f$\{0,\dots,n-1\}\f$
- *
- * \param n Upper bound (excluded)
- * \return \f$x\sim\mathcal U(\{0,\dots,n-1\})\f$
- */
-int rand0n(int n);
-/**
- * \ingroup tensor
- * \brief This is a helper function to sample from a normalized gaussian distribution
- *
- * \return \f$x\sim\mathcal N(0,1)\f$
- */
-real rand_normal();
-
 } // namespace dynet
 
 #endif
diff --git a/examples/batching/rnnlm-batch.h b/examples/batching/rnnlm-batch.h
index e63255746..0e65e6236 100644
--- a/examples/batching/rnnlm-batch.h
+++ b/examples/batching/rnnlm-batch.h
@@ -19,6 +19,7 @@
 #include "dynet/lstm.h"
 #include "dynet/dict.h"
 #include "dynet/expr.h"
+#include "dynet/rand.h"
 
 #include <iostream>
 #include <fstream>
diff --git a/examples/multiprocessing/rnnlm.h b/examples/multiprocessing/rnnlm.h
index 8cc8721a6..494e0ca85 100644
--- a/examples/multiprocessing/rnnlm.h
+++ b/examples/multiprocessing/rnnlm.h
@@ -2,6 +2,7 @@
 #include "dynet/expr.h"
 #include "dynet/dict.h"
 #include "dynet/lstm.h"
+#include "dynet/rand.h"
 
 #include <vector>
 
diff --git a/examples/rnnlm/train_rnnlm.cc b/examples/rnnlm/train_rnnlm.cc
index c1bac1ec3..09cf4eb46 100644
--- a/examples/rnnlm/train_rnnlm.cc
+++ b/examples/rnnlm/train_rnnlm.cc
@@ -11,6 +11,7 @@
 #include "dynet/hsm-builder.h"
 #include "dynet/globals.h"
 #include "dynet/io.h"
+#include "dynet/rand.h"
 #include "getpid.h"
 #include "cl-args.h"
 
diff --git a/examples/sentence-embedding/train_embed-cl.cc b/examples/sentence-embedding/train_embed-cl.cc
index a2a57ec6d..aca84253f 100644
--- a/examples/sentence-embedding/train_embed-cl.cc
+++ b/examples/sentence-embedding/train_embed-cl.cc
@@ -5,6 +5,7 @@
 #include "dynet/dict.h"
 #include "dynet/expr.h"
 #include "dynet/io.h"
+#include "dynet/rand.h"
 
 #include <iostream>
 #include <fstream>
diff --git a/examples/variational-autoencoder/train_rnnlm-aevb.cc b/examples/variational-autoencoder/train_rnnlm-aevb.cc
index 6096f8fd6..b2913487c 100644
--- a/examples/variational-autoencoder/train_rnnlm-aevb.cc
+++ b/examples/variational-autoencoder/train_rnnlm-aevb.cc
@@ -9,6 +9,7 @@
 #include "dynet/expr.h"
 #include "dynet/globals.h"
 #include "dynet/io.h"
+#include "dynet/rand.h"
 #include "getpid.h"
 
 #include <iostream>