diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorBase.h b/unsupported/Eigen/CXX11/src/Tensor/TensorBase.h
index 22fbe9b2c..d44060258 100644
--- a/unsupported/Eigen/CXX11/src/Tensor/TensorBase.h
+++ b/unsupported/Eigen/CXX11/src/Tensor/TensorBase.h
@@ -53,8 +53,8 @@ class TensorBase<Derived, ReadOnlyAccessors>
     }
     template <typename RandomGenerator> EIGEN_DEVICE_FUNC
     EIGEN_STRONG_INLINE const TensorCwiseNullaryOp<RandomGenerator, const Derived>
-    random() const {
-      return nullaryExpr(RandomGenerator());
+    random(const RandomGenerator& gen = RandomGenerator()) const {
+      return nullaryExpr(gen);
     }
 
     // Generic unary operation support.
diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorFunctors.h b/unsupported/Eigen/CXX11/src/Tensor/TensorFunctors.h
index 25f085a59..cb95755a7 100644
--- a/unsupported/Eigen/CXX11/src/Tensor/TensorFunctors.h
+++ b/unsupported/Eigen/CXX11/src/Tensor/TensorFunctors.h
@@ -182,18 +182,45 @@ template <typename T> struct ProdReducer
   }
 };
 
+
+// Random number generation
+namespace {
+int get_random_seed() {
+#if defined _WIN32
+    SYSTEMTIME st;
+    GetSystemTime(&st);
+    return st.wSecond + 1000 * st.wMilliseconds;
+#elif defined __APPLE__
+    return mach_absolute_time();
+#elif defined __CUDA_ARCH__
+    return clock();
+#else
+    timespec ts;
+    clock_gettime(CLOCK_REALTIME, &ts);
+    return ts.tv_nsec;
+#endif
+}
+}
+
 #if !defined (EIGEN_USE_GPU) || !defined(__CUDACC__) || !defined(__CUDA_ARCH__)
 // We're not compiling a cuda kernel
-template <typename T> struct UniformRandomGenerator {
+template <typename T> class UniformRandomGenerator {
 
+ public:
   static const bool PacketAccess = true;
 
+  UniformRandomGenerator(bool deterministic = true) {
+    if (!deterministic) {
+      srand(get_random_seed());
+    }
+  }
+
   template<typename Index>
   T operator()(Index, Index = 0) const {
     return random<T>();
   }
   template<typename Index>
-  typename internal::packet_traits<T>::type packetOp(Index, Index = 0) const {
+  typename internal::packet_traits<T>::type packetOp(Index i, Index j = 0) const {
     const int packetSize = internal::packet_traits<T>::size;
     EIGEN_ALIGN_DEFAULT T values[packetSize];
     for (int i = 0; i < packetSize; ++i) {
@@ -203,26 +230,95 @@ template <typename T> struct UniformRandomGenerator {
   }
 };
 
+#if __cplusplus > 199711
+template <> class UniformRandomGenerator<float> {
+ public:
+  static const bool PacketAccess = true;
+
+  UniformRandomGenerator(bool deterministic = true) {
+    if (!deterministic) {
+      m_generator.seed(get_random_seed());
+    }
+  }
+  UniformRandomGenerator(const UniformRandomGenerator<float>& other) {
+    m_generator.seed(other(0, 0) * UINT_MAX);
+  }
+
+  template<typename Index>
+  float operator()(Index, Index = 0) const {
+    return m_distribution(m_generator);
+  }
+  template<typename Index>
+  typename internal::packet_traits<float>::type packetOp(Index i, Index j = 0) const {
+    const int packetSize = internal::packet_traits<float>::size;
+    EIGEN_ALIGN_DEFAULT float values[packetSize];
+    for (int i = 0; i < packetSize; ++i) {
+      values[i] = this->operator()(i, j);
+    }
+    return internal::pload<typename internal::packet_traits<float>::type>(values);
+  }
+
+ private:
+  UniformRandomGenerator& operator = (const UniformRandomGenerator&);
+  mutable std::mt19937 m_generator;
+  mutable std::uniform_real_distribution<float> m_distribution;
+};
+
+template <> class UniformRandomGenerator<double> {
+ public:
+  static const bool PacketAccess = true;
+
+  UniformRandomGenerator(bool deterministic = true) {
+    if (!deterministic) {
+      m_generator.seed(get_random_seed());
+    }
+  }
+  UniformRandomGenerator(const UniformRandomGenerator<double>& other) {
+    m_generator.seed(other(0, 0) * UINT_MAX);
+  }
+
+  template<typename Index>
+  double operator()(Index, Index = 0) const {
+    return m_distribution(m_generator);
+  }
+  template<typename Index>
+  typename internal::packet_traits<double>::type packetOp(Index i, Index j = 0) const {
+    const int packetSize = internal::packet_traits<double>::size;
+    EIGEN_ALIGN_DEFAULT double values[packetSize];
+    for (int i = 0; i < packetSize; ++i) {
+      values[i] = this->operator()(i, j);
+    }
+    return internal::pload<typename internal::packet_traits<double>::type>(values);
+  }
+
+ private:
+  UniformRandomGenerator& operator = (const UniformRandomGenerator&);
+  mutable std::mt19937 m_generator;
+  mutable std::uniform_real_distribution<double> m_distribution;
+};
+#endif
+
 #else
 
 // We're compiling a cuda kernel
-template <typename T> struct UniformRandomGenerator;
-
-template <> struct UniformRandomGenerator<float> {
+template <typename T> class UniformRandomGenerator;
 
+template <> class UniformRandomGenerator<float> {
+ public:
   static const bool PacketAccess = true;
 
-  EIGEN_DEVICE_FUNC UniformRandomGenerator() {
+  EIGEN_DEVICE_FUNC UniformRandomGenerator(bool deterministic = true) {
     const int tid = blockIdx.x * blockDim.x + threadIdx.x;
-    curand_init(0, tid, 0, &m_state);
+    const int seed = deterministic ? 0 : get_random_seed();
+    curand_init(seed, tid, 0, &m_state);
   }
 
-  template<typename Index> EIGEN_DEVICE_FUNC
-  float operator()(Index, Index = 0) const {
+  template<typename Index>
+  EIGEN_DEVICE_FUNC float operator()(Index, Index = 0) const {
     return curand_uniform(&m_state);
   }
-  template<typename Index> EIGEN_DEVICE_FUNC
-  float4 packetOp(Index, Index = 0) const {
+  template<typename Index>
+  EIGEN_DEVICE_FUNC float4 packetOp(Index, Index = 0) const {
     return curand_uniform4(&m_state);
   }
 
@@ -230,20 +326,21 @@ template <> struct UniformRandomGenerator<float> {
   mutable curandStatePhilox4_32_10_t m_state;
 };
 
-template <> struct UniformRandomGenerator<double> {
-
+template <> class UniformRandomGenerator<double> {
+ public:
   static const bool PacketAccess = true;
 
-  EIGEN_DEVICE_FUNC UniformRandomGenerator() {
+  EIGEN_DEVICE_FUNC UniformRandomGenerator(bool deterministic = true) {
     const int tid = blockIdx.x * blockDim.x + threadIdx.x;
-    curand_init(0, tid, 0, &m_state);
+    const int seed = deterministic ? 0 : get_random_seed();
+    curand_init(seed, tid, 0, &m_state);
   }
-  template<typename Index> EIGEN_DEVICE_FUNC
-  double operator()(Index, Index = 0) const {
+  template<typename Index>
+  EIGEN_DEVICE_FUNC double operator()(Index, Index = 0) const {
     return curand_uniform_double(&m_state);
   }
-  template<typename Index> EIGEN_DEVICE_FUNC
-  double2 packetOp(Index, Index = 0) const {
+  template<typename Index>
+  EIGEN_DEVICE_FUNC double2 packetOp(Index, Index = 0) const {
     return curand_uniform2_double(&m_state);
   }
 
@@ -256,12 +353,18 @@ template <> struct UniformRandomGenerator<double> {
 
 #if (!defined (EIGEN_USE_GPU) || !defined(__CUDACC__) || !defined(__CUDA_ARCH__)) && __cplusplus > 199711
 // We're not compiling a cuda kernel
-template <typename T> struct NormalRandomGenerator {
-
+template <typename T> class NormalRandomGenerator {
+ public:
   static const bool PacketAccess = true;
 
-  NormalRandomGenerator() : m_distribution(0, 1) {}
-  NormalRandomGenerator(const NormalRandomGenerator& other) : m_distribution(other.m_distribution) { }
+  NormalRandomGenerator(bool deterministic = true) : m_distribution(0, 1) {
+    if (!deterministic) {
+      m_generator.seed(get_random_seed());
+    }
+  }
+  NormalRandomGenerator(const NormalRandomGenerator& other) : m_distribution(other.m_distribution) {
+    m_generator.seed(other(0, 0) * UINT_MAX);
+  }
 
   template<typename Index>
   T operator()(Index, Index = 0) const {
@@ -278,29 +381,30 @@ template <typename T> struct NormalRandomGenerator {
   }
 
   mutable std::normal_distribution<T> m_distribution;
-  mutable std::default_random_engine m_generator;
+  mutable std::mt19937 m_generator;
 };
 
 #elif defined (EIGEN_USE_GPU) && defined(__CUDACC__) && defined(__CUDA_ARCH__)
 
 // We're compiling a cuda kernel
-template <typename T> struct NormalRandomGenerator;
-
-template <> struct NormalRandomGenerator<float> {
+template <typename T> class NormalRandomGenerator;
 
+template <> class NormalRandomGenerator<float> {
+ public:
   static const bool PacketAccess = true;
 
-  EIGEN_DEVICE_FUNC NormalRandomGenerator() {
+  EIGEN_DEVICE_FUNC NormalRandomGenerator(bool deterministic = true) {
     const int tid = blockIdx.x * blockDim.x + threadIdx.x;
-    curand_init(0, tid, 0, &m_state);
+    const int seed = deterministic ? 0 : get_random_seed();
+    curand_init(seed, tid, 0, &m_state);
   }
 
-  template<typename Index> EIGEN_DEVICE_FUNC
-  float operator()(Index, Index = 0) const {
+  template<typename Index>
+  EIGEN_DEVICE_FUNC float operator()(Index, Index = 0) const {
     return curand_normal(&m_state);
   }
-  template<typename Index> EIGEN_DEVICE_FUNC
-  float4 packetOp(Index, Index = 0) const {
+  template<typename Index>
+  EIGEN_DEVICE_FUNC float4 packetOp(Index, Index = 0) const {
     return curand_normal4(&m_state);
   }
 
@@ -308,20 +412,21 @@ template <> struct NormalRandomGenerator<float> {
   mutable curandStatePhilox4_32_10_t m_state;
 };
 
-template <> struct NormalRandomGenerator<double> {
-
+template <> class NormalRandomGenerator<double> {
+ public:
   static const bool PacketAccess = true;
 
-  EIGEN_DEVICE_FUNC NormalRandomGenerator() {
+  EIGEN_DEVICE_FUNC NormalRandomGenerator(bool deterministic = true) {
     const int tid = blockIdx.x * blockDim.x + threadIdx.x;
-    curand_init(0, tid, 0, &m_state);
+    const int seed = deterministic ? 0 : get_random_seed();
+    curand_init(seed, tid, 0, &m_state);
   }
-  template<typename Index> EIGEN_DEVICE_FUNC
-  double operator()(Index, Index = 0) const {
+  template<typename Index>
+  EIGEN_DEVICE_FUNC double operator()(Index, Index = 0) const {
     return curand_normal_double(&m_state);
   }
-  template<typename Index> EIGEN_DEVICE_FUNC
-  double2 packetOp(Index, Index = 0) const {
+  template<typename Index>
+  EIGEN_DEVICE_FUNC double2 packetOp(Index, Index = 0) const {
     return curand_normal2_double(&m_state);
   }
 
@@ -329,6 +434,13 @@ template <> struct NormalRandomGenerator<double> {
   mutable curandStatePhilox4_32_10_t m_state;
 };
 
+#else
+
+template <typename T> class NormalRandomGenerator {
+ public:
+  NormalRandomGenerator(bool = true) {}
+};
+
 #endif