diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h
index 32d1499c6..babfc86a9 100644
--- a/Eigen/src/Geometry/Quaternion.h
+++ b/Eigen/src/Geometry/Quaternion.h
@@ -271,6 +271,8 @@ public:
   explicit inline Quaternion(const Quaternion<OtherScalar, OtherOptions>& other)
   { m_coeffs = other.coeffs().template cast<Scalar>(); }
 
+  static Quaternion UniformRandom();
+
   template<typename Derived1, typename Derived2>
   static Quaternion FromTwoVectors(const MatrixBase<Derived1>& a, const MatrixBase<Derived2>& b);
 
@@ -609,6 +611,26 @@ inline Derived& QuaternionBase<Derived>::setFromTwoVectors(const MatrixBase<Deri
   return derived();
 }
 
+/** Returns a random quaternion following a uniform distribution law.
+  *
+  * \returns resulting quaternion
+  *
+  * \note The implementation is based on http://planning.cs.uiuc.edu/node198.html
+  */
+template<typename Scalar, int Options>
+Quaternion<Scalar,Options> Quaternion<Scalar,Options>::UniformRandom()
+{
+  const Scalar u1 = internal::random<Scalar>(0,1),
+               u2 = internal::random<Scalar>(0, 2*M_PI),
+               u3 = internal::random<Scalar>(0, 2*M_PI);
+  const Scalar a = std::sqrt (1 - u1),
+               b = std::sqrt (u1);
+  return Quaternion (a * std::sin (u2),
+                     a * std::cos (u2),
+                     b * std::sin (u3),
+                     b * std::cos (u3));
+}
+
 
 /** Returns a quaternion representing a rotation between
   * the two arbitrary vectors \a a and \a b. In other words, the built
diff --git a/test/geo_quaternion.cpp b/test/geo_quaternion.cpp
index 25130c19a..7c06473a5 100644
--- a/test/geo_quaternion.cpp
+++ b/test/geo_quaternion.cpp
@@ -50,7 +50,6 @@ template<typename Scalar, int Options> void quaternion(void)
   using std::abs;
   typedef Matrix<Scalar,3,1> Vector3;
   typedef Matrix<Scalar,3,3> Matrix3;
-  typedef Matrix<Scalar,4,1> Vector4;
   typedef Quaternion<Scalar,Options> Quaternionx;
   typedef AngleAxis<Scalar> AngleAxisx;
 
@@ -157,8 +156,8 @@ template<typename Scalar, int Options> void quaternion(void)
   Quaternionx *q = new Quaternionx;
   delete q;
 
-  q1 = AngleAxisx(a, v0.normalized());
-  q2 = AngleAxisx(b, v1.normalized());
+  q1 = Quaternionx::UniformRandom ();
+  q2 = Quaternionx::UniformRandom ();
   check_slerp(q1,q2);
 
   q1 = AngleAxisx(b, v1.normalized());
@@ -169,7 +168,7 @@ template<typename Scalar, int Options> void quaternion(void)
   q2 = AngleAxisx(-b, -v1.normalized());
   check_slerp(q1,q2);
 
-  q1.coeffs() = Vector4::Random().normalized();
+  q1 = Quaternionx::UniformRandom ();
   q2.coeffs() = -q1.coeffs();
   check_slerp(q1,q2);
 }