diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index 54bd2d58a..b794b0c43 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -206,7 +206,7 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
       // in "Analysis of the Cholesky Decomposition of a Semi-definite Matrix" by
       // Nicholas J. Higham. Also see "Accuracy and Stability of Numerical
       // Algorithms" page 217, also by Higham.
-      cutoff = ei_abs(epsilon<Scalar>() * RealScalar(size) * biggest_in_corner);
+      cutoff = ei_abs(NumTraits<Scalar>::epsilon() * RealScalar(size) * biggest_in_corner);
 
       m_sign = ei_real(m_matrix.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
     }
diff --git a/Eigen/src/Core/DenseBase.h b/Eigen/src/Core/DenseBase.h
index 9ef9108e4..700c11929 100644
--- a/Eigen/src/Core/DenseBase.h
+++ b/Eigen/src/Core/DenseBase.h
@@ -381,17 +381,17 @@ template<typename Derived> class DenseBase
 
     template<typename OtherDerived>
     bool isApprox(const DenseBase<OtherDerived>& other,
-                  RealScalar prec = dummy_precision<Scalar>()) const;
+                  RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
     bool isMuchSmallerThan(const RealScalar& other,
-                           RealScalar prec = dummy_precision<Scalar>()) const;
+                           RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
     template<typename OtherDerived>
     bool isMuchSmallerThan(const DenseBase<OtherDerived>& other,
-                           RealScalar prec = dummy_precision<Scalar>()) const;
+                           RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
 
-    bool isApproxToConstant(const Scalar& value, RealScalar prec = dummy_precision<Scalar>()) const;
-    bool isConstant(const Scalar& value, RealScalar prec = dummy_precision<Scalar>()) const;
-    bool isZero(RealScalar prec = dummy_precision<Scalar>()) const;
-    bool isOnes(RealScalar prec = dummy_precision<Scalar>()) const;
+    bool isApproxToConstant(const Scalar& value, RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
+    bool isConstant(const Scalar& value, RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
+    bool isZero(RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
+    bool isOnes(RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
 
     inline Derived& operator*=(const Scalar& other);
     inline Derived& operator/=(const Scalar& other);
diff --git a/Eigen/src/Core/IO.h b/Eigen/src/Core/IO.h
index 8f91e272e..3e8d2bc66 100644
--- a/Eigen/src/Core/IO.h
+++ b/Eigen/src/Core/IO.h
@@ -146,7 +146,7 @@ std::ostream & ei_print_matrix(std::ostream & s, const Derived& _m, const IOForm
     if (NumTraits<Scalar>::HasFloatingPoint)
     {
       typedef typename NumTraits<Scalar>::Real RealScalar;
-      RealScalar explicit_precision_fp = std::ceil(-ei_log(epsilon<Scalar>())/ei_log(10.0));
+      RealScalar explicit_precision_fp = std::ceil(-ei_log(NumTraits<Scalar>::epsilon())/ei_log(10.0));
       explicit_precision = static_cast<std::streamsize>(explicit_precision_fp);
     }
     else
diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h
index 1ea0116b4..1a4561555 100644
--- a/Eigen/src/Core/MathFunctions.h
+++ b/Eigen/src/Core/MathFunctions.h
@@ -25,13 +25,6 @@
 #ifndef EIGEN_MATHFUNCTIONS_H
 #define EIGEN_MATHFUNCTIONS_H
 
-template<typename T> inline typename NumTraits<T>::Real epsilon()
-{
- return std::numeric_limits<typename NumTraits<T>::Real>::epsilon();
-}
-
-template<typename T> inline typename NumTraits<T>::Real dummy_precision();
-
 template<typename T> inline T ei_random(T a, T b);
 template<typename T> inline T ei_random();
 template<typename T> inline T ei_random_amplitude()
@@ -55,7 +48,6 @@ template<typename T> inline typename NumTraits<T>::Real ei_hypot(T x, T y)
 ***   int   ***
 **************/
 
-template<> inline int dummy_precision<int>() { return 0; }
 inline int ei_real(int x)  { return x; }
 inline int& ei_real_ref(int& x)  { return x; }
 inline int ei_imag(int)    { return 0; }
@@ -92,15 +84,15 @@ template<> inline int ei_random()
 {
   return ei_random<int>(-ei_random_amplitude<int>(), ei_random_amplitude<int>());
 }
-inline bool ei_isMuchSmallerThan(int a, int, int = dummy_precision<int>())
+inline bool ei_isMuchSmallerThan(int a, int, int = NumTraits<int>::dummy_precision())
 {
   return a == 0;
 }
-inline bool ei_isApprox(int a, int b, int = dummy_precision<int>())
+inline bool ei_isApprox(int a, int b, int = NumTraits<int>::dummy_precision())
 {
   return a == b;
 }
-inline bool ei_isApproxOrLessThan(int a, int b, int = dummy_precision<int>())
+inline bool ei_isApproxOrLessThan(int a, int b, int = NumTraits<int>::dummy_precision())
 {
   return a <= b;
 }
@@ -109,7 +101,6 @@ inline bool ei_isApproxOrLessThan(int a, int b, int = dummy_precision<int>())
 *** float   ***
 **************/
 
-template<> inline float dummy_precision<float>() { return 1e-5f; }
 inline float ei_real(float x)  { return x; }
 inline float& ei_real_ref(float& x)  { return x; }
 inline float ei_imag(float)    { return 0.f; }
@@ -140,15 +131,15 @@ template<> inline float ei_random()
 {
   return ei_random<float>(-ei_random_amplitude<float>(), ei_random_amplitude<float>());
 }
-inline bool ei_isMuchSmallerThan(float a, float b, float prec = dummy_precision<float>())
+inline bool ei_isMuchSmallerThan(float a, float b, float prec = NumTraits<float>::dummy_precision())
 {
   return ei_abs(a) <= ei_abs(b) * prec;
 }
-inline bool ei_isApprox(float a, float b, float prec = dummy_precision<float>())
+inline bool ei_isApprox(float a, float b, float prec = NumTraits<float>::dummy_precision())
 {
   return ei_abs(a - b) <= std::min(ei_abs(a), ei_abs(b)) * prec;
 }
-inline bool ei_isApproxOrLessThan(float a, float b, float prec = dummy_precision<float>())
+inline bool ei_isApproxOrLessThan(float a, float b, float prec = NumTraits<float>::dummy_precision())
 {
   return a <= b || ei_isApprox(a, b, prec);
 }
@@ -157,8 +148,6 @@ inline bool ei_isApproxOrLessThan(float a, float b, float prec = dummy_precision
 *** double  ***
 **************/
 
-template<> inline double dummy_precision<double>() { return 1e-12; }
-
 inline double ei_real(double x)  { return x; }
 inline double& ei_real_ref(double& x)  { return x; }
 inline double ei_imag(double)    { return 0.; }
@@ -189,15 +178,15 @@ template<> inline double ei_random()
 {
   return ei_random<double>(-ei_random_amplitude<double>(), ei_random_amplitude<double>());
 }
-inline bool ei_isMuchSmallerThan(double a, double b, double prec = dummy_precision<double>())
+inline bool ei_isMuchSmallerThan(double a, double b, double prec = NumTraits<double>::dummy_precision())
 {
   return ei_abs(a) <= ei_abs(b) * prec;
 }
-inline bool ei_isApprox(double a, double b, double prec = dummy_precision<double>())
+inline bool ei_isApprox(double a, double b, double prec = NumTraits<double>::dummy_precision())
 {
   return ei_abs(a - b) <= std::min(ei_abs(a), ei_abs(b)) * prec;
 }
-inline bool ei_isApproxOrLessThan(double a, double b, double prec = dummy_precision<double>())
+inline bool ei_isApproxOrLessThan(double a, double b, double prec = NumTraits<double>::dummy_precision())
 {
   return a <= b || ei_isApprox(a, b, prec);
 }
@@ -206,7 +195,6 @@ inline bool ei_isApproxOrLessThan(double a, double b, double prec = dummy_precis
 *** complex<float> ***
 *********************/
 
-template<> inline float dummy_precision<std::complex<float> >() { return dummy_precision<float>(); }
 inline float ei_real(const std::complex<float>& x) { return std::real(x); }
 inline float ei_imag(const std::complex<float>& x) { return std::imag(x); }
 inline float& ei_real_ref(std::complex<float>& x) { return reinterpret_cast<float*>(&x)[0]; }
@@ -225,15 +213,15 @@ template<> inline std::complex<float> ei_random()
 {
   return std::complex<float>(ei_random<float>(), ei_random<float>());
 }
-inline bool ei_isMuchSmallerThan(const std::complex<float>& a, const std::complex<float>& b, float prec = dummy_precision<float>())
+inline bool ei_isMuchSmallerThan(const std::complex<float>& a, const std::complex<float>& b, float prec = NumTraits<float>::dummy_precision())
 {
   return ei_abs2(a) <= ei_abs2(b) * prec * prec;
 }
-inline bool ei_isMuchSmallerThan(const std::complex<float>& a, float b, float prec = dummy_precision<float>())
+inline bool ei_isMuchSmallerThan(const std::complex<float>& a, float b, float prec = NumTraits<float>::dummy_precision())
 {
   return ei_abs2(a) <= ei_abs2(b) * prec * prec;
 }
-inline bool ei_isApprox(const std::complex<float>& a, const std::complex<float>& b, float prec = dummy_precision<float>())
+inline bool ei_isApprox(const std::complex<float>& a, const std::complex<float>& b, float prec = NumTraits<float>::dummy_precision())
 {
   return ei_isApprox(ei_real(a), ei_real(b), prec)
       && ei_isApprox(ei_imag(a), ei_imag(b), prec);
@@ -244,7 +232,6 @@ inline bool ei_isApprox(const std::complex<float>& a, const std::complex<float>&
 *** complex<double> ***
 **********************/
 
-template<> inline double dummy_precision<std::complex<double> >() { return dummy_precision<double>(); }
 inline double ei_real(const std::complex<double>& x) { return std::real(x); }
 inline double ei_imag(const std::complex<double>& x) { return std::imag(x); }
 inline double& ei_real_ref(std::complex<double>& x) { return reinterpret_cast<double*>(&x)[0]; }
@@ -263,15 +250,15 @@ template<> inline std::complex<double> ei_random()
 {
   return std::complex<double>(ei_random<double>(), ei_random<double>());
 }
-inline bool ei_isMuchSmallerThan(const std::complex<double>& a, const std::complex<double>& b, double prec = dummy_precision<double>())
+inline bool ei_isMuchSmallerThan(const std::complex<double>& a, const std::complex<double>& b, double prec = NumTraits<double>::dummy_precision())
 {
   return ei_abs2(a) <= ei_abs2(b) * prec * prec;
 }
-inline bool ei_isMuchSmallerThan(const std::complex<double>& a, double b, double prec = dummy_precision<double>())
+inline bool ei_isMuchSmallerThan(const std::complex<double>& a, double b, double prec = NumTraits<double>::dummy_precision())
 {
   return ei_abs2(a) <= ei_abs2(b) * prec * prec;
 }
-inline bool ei_isApprox(const std::complex<double>& a, const std::complex<double>& b, double prec = dummy_precision<double>())
+inline bool ei_isApprox(const std::complex<double>& a, const std::complex<double>& b, double prec = NumTraits<double>::dummy_precision())
 {
   return ei_isApprox(ei_real(a), ei_real(b), prec)
       && ei_isApprox(ei_imag(a), ei_imag(b), prec);
@@ -283,7 +270,6 @@ inline bool ei_isApprox(const std::complex<double>& a, const std::complex<double
 *** long double ***
 ******************/
 
-template<> inline long double dummy_precision<long double>() { return dummy_precision<double>(); }
 inline long double ei_real(long double x)  { return x; }
 inline long double& ei_real_ref(long double& x)  { return x; }
 inline long double ei_imag(long double)    { return 0.; }
@@ -306,15 +292,15 @@ template<> inline long double ei_random()
 {
   return ei_random<double>(-ei_random_amplitude<double>(), ei_random_amplitude<double>());
 }
-inline bool ei_isMuchSmallerThan(long double a, long double b, long double prec = dummy_precision<long double>())
+inline bool ei_isMuchSmallerThan(long double a, long double b, long double prec = NumTraits<long double>::dummy_precision())
 {
   return ei_abs(a) <= ei_abs(b) * prec;
 }
-inline bool ei_isApprox(long double a, long double b, long double prec = dummy_precision<long double>())
+inline bool ei_isApprox(long double a, long double b, long double prec = NumTraits<long double>::dummy_precision())
 {
   return ei_abs(a - b) <= std::min(ei_abs(a), ei_abs(b)) * prec;
 }
-inline bool ei_isApproxOrLessThan(long double a, long double b, long double prec = dummy_precision<long double>())
+inline bool ei_isApproxOrLessThan(long double a, long double b, long double prec = NumTraits<long double>::dummy_precision())
 {
   return a <= b || ei_isApprox(a, b, prec);
 }
@@ -323,7 +309,6 @@ inline bool ei_isApproxOrLessThan(long double a, long double b, long double prec
 ***  bool  ***
 **************/
 
-template<> inline bool dummy_precision<bool>() { return 0; }
 inline bool ei_real(bool x)  { return x; }
 inline bool& ei_real_ref(bool& x)  { return x; }
 inline bool ei_imag(bool)    { return 0; }
@@ -336,15 +321,15 @@ template<> inline bool ei_random()
 {
   return (ei_random<int>(0,1) == 1);
 }
-inline bool ei_isMuchSmallerThan(bool a, bool, bool = dummy_precision<bool>())
+inline bool ei_isMuchSmallerThan(bool a, bool, bool = NumTraits<bool>::dummy_precision())
 {
   return !a;
 }
-inline bool ei_isApprox(bool a, bool b, bool = dummy_precision<bool>())
+inline bool ei_isApprox(bool a, bool b, bool = NumTraits<bool>::dummy_precision())
 {
   return a == b;
 }
-inline bool ei_isApproxOrLessThan(bool a, bool b, bool = dummy_precision<bool>())
+inline bool ei_isApproxOrLessThan(bool a, bool b, bool = NumTraits<bool>::dummy_precision())
 {
   return int(a) <= int(b);
 }
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index 65b3bfb7c..229195046 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -253,16 +253,16 @@ template<typename Derived> class MatrixBase
     Derived& setIdentity();
     Derived& setIdentity(int rows, int cols);
 
-    bool isIdentity(RealScalar prec = dummy_precision<Scalar>()) const;
-    bool isDiagonal(RealScalar prec = dummy_precision<Scalar>()) const;
+    bool isIdentity(RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
+    bool isDiagonal(RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
 
-    bool isUpperTriangular(RealScalar prec = dummy_precision<Scalar>()) const;
-    bool isLowerTriangular(RealScalar prec = dummy_precision<Scalar>()) const;
+    bool isUpperTriangular(RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
+    bool isLowerTriangular(RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
 
     template<typename OtherDerived>
     bool isOrthogonal(const MatrixBase<OtherDerived>& other,
-                      RealScalar prec = dummy_precision<Scalar>()) const;
-    bool isUnitary(RealScalar prec = dummy_precision<Scalar>()) const;
+                      RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
+    bool isUnitary(RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
 
     /** \returns true if each coefficients of \c *this and \a other are all exactly equal.
       * \warning When using floating point scalar values you probably should rather use a
@@ -310,13 +310,13 @@ template<typename Derived> class MatrixBase
       ResultType& inverse,
       typename ResultType::Scalar& determinant,
       bool& invertible,
-      const RealScalar& absDeterminantThreshold = dummy_precision<Scalar>()
+      const RealScalar& absDeterminantThreshold = NumTraits<Scalar>::dummy_precision()
     ) const;
     template<typename ResultType>
     void computeInverseWithCheck(
       ResultType& inverse,
       bool& invertible,
-      const RealScalar& absDeterminantThreshold = dummy_precision<Scalar>()
+      const RealScalar& absDeterminantThreshold = NumTraits<Scalar>::dummy_precision()
     ) const;
     Scalar determinant() const;
 
diff --git a/Eigen/src/Core/NumTraits.h b/Eigen/src/Core/NumTraits.h
index 304e2c1d6..7279e156d 100644
--- a/Eigen/src/Core/NumTraits.h
+++ b/Eigen/src/Core/NumTraits.h
@@ -49,6 +49,7 @@
 template<typename T> struct NumTraits;
 
 template<> struct NumTraits<int>
+  : std::numeric_limits<int>
 {
   typedef int Real;
   typedef double FloatingPoint;
@@ -60,9 +61,12 @@ template<> struct NumTraits<int>
     AddCost = 1,
     MulCost = 1
   };
+
+  inline static int dummy_precision() { return 0; }
 };
 
 template<> struct NumTraits<float>
+  : std::numeric_limits<float>
 {
   typedef float Real;
   typedef float FloatingPoint;
@@ -74,9 +78,12 @@ template<> struct NumTraits<float>
     AddCost = 1,
     MulCost = 1
   };
+
+  inline static float dummy_precision() { return 1e-5f; }
 };
 
 template<> struct NumTraits<double>
+  : std::numeric_limits<double>
 {
   typedef double Real;
   typedef double FloatingPoint;
@@ -88,9 +95,12 @@ template<> struct NumTraits<double>
     AddCost = 1,
     MulCost = 1
   };
+
+  inline static double dummy_precision() { return 1e-12; }
 };
 
 template<typename _Real> struct NumTraits<std::complex<_Real> >
+  : std::numeric_limits<std::complex<_Real> >
 {
   typedef _Real Real;
   typedef std::complex<_Real> FloatingPoint;
@@ -102,9 +112,13 @@ template<typename _Real> struct NumTraits<std::complex<_Real> >
     AddCost = 2 * NumTraits<Real>::AddCost,
     MulCost = 4 * NumTraits<Real>::MulCost + 2 * NumTraits<Real>::AddCost
   };
+
+  inline static Real epsilon() { return std::numeric_limits<Real>::epsilon(); }
+  inline static Real dummy_precision() { return NumTraits<Real>::dummy_precision(); }
 };
 
 template<> struct NumTraits<long long int>
+  : std::numeric_limits<long long int>
 {
   typedef long long int Real;
   typedef long double FloatingPoint;
@@ -119,6 +133,7 @@ template<> struct NumTraits<long long int>
 };
 
 template<> struct NumTraits<long double>
+  : std::numeric_limits<long double>
 {
   typedef long double Real;
   typedef long double FloatingPoint;
@@ -130,9 +145,12 @@ template<> struct NumTraits<long double>
     AddCost = 1,
     MulCost = 1
   };
+
+  static inline long double dummy_precision() { return NumTraits<double>::dummy_precision(); }
 };
 
 template<> struct NumTraits<bool>
+  : std::numeric_limits<bool>
 {
   typedef bool Real;
   typedef float FloatingPoint;
@@ -144,6 +162,8 @@ template<> struct NumTraits<bool>
     AddCost = 1,
     MulCost = 1
   };
+
+  inline static bool dummy_precision() { return 0; }
 };
 
 #endif // EIGEN_NUMTRAITS_H
diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
index d55dc2a96..27c52b4dc 100644
--- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h
+++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
@@ -94,7 +94,7 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
   m_eivalues.resize(n,1);
   m_eivec.resize(n,n);
 
-  RealScalar eps = epsilon<RealScalar>();
+  RealScalar eps = NumTraits<RealScalar>::epsilon();
 
   // Reduce to complex Schur form
   ComplexSchur<MatrixType> schur(matrix);
diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h
index 38a1f56dc..5deac3247 100644
--- a/Eigen/src/Eigenvalues/ComplexSchur.h
+++ b/Eigen/src/Eigenvalues/ComplexSchur.h
@@ -159,7 +159,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
   RealScalar d,sd,sf;
   Complex c,b,disc,r1,r2,kappa;
 
-  RealScalar eps = epsilon<RealScalar>();
+  RealScalar eps = NumTraits<RealScalar>::epsilon();
 
   int iter = 0;
   while(true)
diff --git a/Eigen/src/Geometry/AlignedBox.h b/Eigen/src/Geometry/AlignedBox.h
index 1865b179b..5c1f992c6 100644
--- a/Eigen/src/Geometry/AlignedBox.h
+++ b/Eigen/src/Geometry/AlignedBox.h
@@ -171,7 +171,7 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim)
     * determined by \a prec.
     *
     * \sa MatrixBase::isApprox() */
-  bool isApprox(const AlignedBox& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const
+  bool isApprox(const AlignedBox& other, typename NumTraits<Scalar>::Real prec = NumTraits<Scalar>::dummy_precision()) const
   { return m_min.isApprox(other.m_min, prec) && m_max.isApprox(other.m_max, prec); }
 
 protected:
diff --git a/Eigen/src/Geometry/AngleAxis.h b/Eigen/src/Geometry/AngleAxis.h
index 284c60a74..d02ecf4c5 100644
--- a/Eigen/src/Geometry/AngleAxis.h
+++ b/Eigen/src/Geometry/AngleAxis.h
@@ -146,7 +146,7 @@ public:
     * determined by \a prec.
     *
     * \sa MatrixBase::isApprox() */
-  bool isApprox(const AngleAxis& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const
+  bool isApprox(const AngleAxis& other, typename NumTraits<Scalar>::Real prec = NumTraits<Scalar>::dummy_precision()) const
   { return m_axis.isApprox(other.m_axis, prec) && ei_isApprox(m_angle,other.m_angle, prec); }
 };
 
@@ -165,7 +165,7 @@ template<typename QuatDerived>
 AngleAxis<Scalar>& AngleAxis<Scalar>::operator=(const QuaternionBase<QuatDerived>& q)
 {
   Scalar n2 = q.vec().squaredNorm();
-  if (n2 < dummy_precision<Scalar>()*dummy_precision<Scalar>())
+  if (n2 < NumTraits<Scalar>::dummy_precision()*NumTraits<Scalar>::dummy_precision())
   {
     m_angle = 0;
     m_axis << 1, 0, 0;
diff --git a/Eigen/src/Geometry/EulerAngles.h b/Eigen/src/Geometry/EulerAngles.h
index b6dbf8ae9..13d23761a 100644
--- a/Eigen/src/Geometry/EulerAngles.h
+++ b/Eigen/src/Geometry/EulerAngles.h
@@ -50,7 +50,7 @@ MatrixBase<Derived>::eulerAngles(int a0, int a1, int a2) const
 
   Matrix<Scalar,3,1> res;
   typedef Matrix<typename Derived::Scalar,2,1> Vector2;
-  const Scalar epsilon = dummy_precision<Scalar>();
+  const Scalar epsilon = NumTraits<Scalar>::dummy_precision();
 
   const int odd = ((a0+1)%3 == a1) ? 0 : 1;
   const int i = a0;
diff --git a/Eigen/src/Geometry/Hyperplane.h b/Eigen/src/Geometry/Hyperplane.h
index aab3d5b35..1d0b299ba 100644
--- a/Eigen/src/Geometry/Hyperplane.h
+++ b/Eigen/src/Geometry/Hyperplane.h
@@ -257,7 +257,7 @@ public:
     * determined by \a prec.
     *
     * \sa MatrixBase::isApprox() */
-  bool isApprox(const Hyperplane& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const
+  bool isApprox(const Hyperplane& other, typename NumTraits<Scalar>::Real prec = NumTraits<Scalar>::dummy_precision()) const
   { return m_coeffs.isApprox(other.m_coeffs, prec); }
 
 protected:
diff --git a/Eigen/src/Geometry/ParametrizedLine.h b/Eigen/src/Geometry/ParametrizedLine.h
index 21a5595b9..1846a440a 100644
--- a/Eigen/src/Geometry/ParametrizedLine.h
+++ b/Eigen/src/Geometry/ParametrizedLine.h
@@ -123,7 +123,7 @@ public:
     * determined by \a prec.
     *
     * \sa MatrixBase::isApprox() */
-  bool isApprox(const ParametrizedLine& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const
+  bool isApprox(const ParametrizedLine& other, typename NumTraits<Scalar>::Real prec = NumTraits<Scalar>::dummy_precision()) const
   { return m_origin.isApprox(other.m_origin, prec) && m_direction.isApprox(other.m_direction, prec); }
 
 protected:
diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h
index 98b91e2de..67319d15b 100644
--- a/Eigen/src/Geometry/Quaternion.h
+++ b/Eigen/src/Geometry/Quaternion.h
@@ -163,7 +163,7 @@ public:
     *
     * \sa MatrixBase::isApprox() */
   template<class OtherDerived>
-  bool isApprox(const QuaternionBase<OtherDerived>& other, RealScalar prec = dummy_precision<Scalar>()) const
+  bool isApprox(const QuaternionBase<OtherDerived>& other, RealScalar prec = NumTraits<Scalar>::dummy_precision()) const
   { return coeffs().isApprox(other.coeffs(), prec); }
 
 	/** return the result vector of \a v through the rotation*/
@@ -508,7 +508,7 @@ inline Derived& QuaternionBase<Derived>::setFromTwoVectors(const MatrixBase<Deri
   //    under the constraint:
   //       ||x|| = 1
   //    which yields a singular value problem
-  if (c < Scalar(-1)+dummy_precision<Scalar>())
+  if (c < Scalar(-1)+NumTraits<Scalar>::dummy_precision())
   {
     c = std::max<Scalar>(c,-1);
     Matrix<Scalar,2,3> m; m << v0.transpose(), v1.transpose();
@@ -584,7 +584,7 @@ template <class OtherDerived>
 Quaternion<typename ei_traits<Derived>::Scalar>
 QuaternionBase<Derived>::slerp(Scalar t, const QuaternionBase<OtherDerived>& other) const
 {
-  static const Scalar one = Scalar(1) - epsilon<Scalar>();
+  static const Scalar one = Scalar(1) - NumTraits<Scalar>::epsilon();
   Scalar d = this->dot(other);
   Scalar absD = ei_abs(d);
 
diff --git a/Eigen/src/Geometry/Rotation2D.h b/Eigen/src/Geometry/Rotation2D.h
index 7f24a3eae..73f731f76 100644
--- a/Eigen/src/Geometry/Rotation2D.h
+++ b/Eigen/src/Geometry/Rotation2D.h
@@ -121,7 +121,7 @@ public:
     * determined by \a prec.
     *
     * \sa MatrixBase::isApprox() */
-  bool isApprox(const Rotation2D& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const
+  bool isApprox(const Rotation2D& other, typename NumTraits<Scalar>::Real prec = NumTraits<Scalar>::dummy_precision()) const
   { return ei_isApprox(m_angle,other.m_angle, prec); }
 };
 
diff --git a/Eigen/src/Geometry/Scaling.h b/Eigen/src/Geometry/Scaling.h
index 4695914fd..2ff8eaba3 100644
--- a/Eigen/src/Geometry/Scaling.h
+++ b/Eigen/src/Geometry/Scaling.h
@@ -107,7 +107,7 @@ public:
     * determined by \a prec.
     *
     * \sa MatrixBase::isApprox() */
-  bool isApprox(const UniformScaling& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const
+  bool isApprox(const UniformScaling& other, typename NumTraits<Scalar>::Real prec = NumTraits<Scalar>::dummy_precision()) const
   { return ei_isApprox(m_factor, other.factor(), prec); }
 
 };
diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h
index 309c1e1fd..40bc7033a 100644
--- a/Eigen/src/Geometry/Transform.h
+++ b/Eigen/src/Geometry/Transform.h
@@ -429,7 +429,7 @@ public:
     * determined by \a prec.
     *
     * \sa MatrixBase::isApprox() */
-  bool isApprox(const Transform& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const
+  bool isApprox(const Transform& other, typename NumTraits<Scalar>::Real prec = NumTraits<Scalar>::dummy_precision()) const
   { return m_matrix.isApprox(other.m_matrix, prec); }
 
   /** Sets the last row to [0 ... 0 1]
diff --git a/Eigen/src/Geometry/Translation.h b/Eigen/src/Geometry/Translation.h
index b7477df9f..1e3906bde 100644
--- a/Eigen/src/Geometry/Translation.h
+++ b/Eigen/src/Geometry/Translation.h
@@ -154,7 +154,7 @@ public:
     * determined by \a prec.
     *
     * \sa MatrixBase::isApprox() */
-  bool isApprox(const Translation& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const
+  bool isApprox(const Translation& other, typename NumTraits<Scalar>::Real prec = NumTraits<Scalar>::dummy_precision()) const
   { return m_coeffs.isApprox(other.m_coeffs, prec); }
 
 };
diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index 657cbbf1c..72e878223 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -276,7 +276,7 @@ template<typename _MatrixType> class FullPivLU
       return m_usePrescribedThreshold ? m_prescribedThreshold
       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
       // and turns out to be identical to Higham's formula used already in LDLt.
-                                      : epsilon<Scalar>() * m_lu.diagonalSize();
+                                      : NumTraits<Scalar>::epsilon() * m_lu.diagonalSize();
     }
 
     /** \returns the rank of the matrix of which *this is the LU decomposition.
diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h
index 63d069081..bf28605dd 100644
--- a/Eigen/src/QR/ColPivHouseholderQR.h
+++ b/Eigen/src/QR/ColPivHouseholderQR.h
@@ -282,7 +282,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
       return m_usePrescribedThreshold ? m_prescribedThreshold
       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
       // and turns out to be identical to Higham's formula used already in LDLt.
-                                      : epsilon<Scalar>() * m_qr.diagonalSize();
+                                      : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
     }
 
     /** \returns the number of nonzero pivots in the QR decomposition.
@@ -350,7 +350,7 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
   for(int k = 0; k < cols; ++k)
     colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
 
-  RealScalar threshold_helper = colSqNorms.maxCoeff() * ei_abs2(epsilon<Scalar>()) / rows;
+  RealScalar threshold_helper = colSqNorms.maxCoeff() * ei_abs2(NumTraits<Scalar>::epsilon()) / rows;
 
   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
   m_maxpivot = RealScalar(0);
diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h
index 8e8dbbbf8..1ec60aeaf 100644
--- a/Eigen/src/QR/FullPivHouseholderQR.h
+++ b/Eigen/src/QR/FullPivHouseholderQR.h
@@ -270,7 +270,7 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
 
   RowVectorType temp(cols);
 
-  m_precision = epsilon<Scalar>() * size;
+  m_precision = NumTraits<Scalar>::epsilon() * size;
 
   m_rows_transpositions.resize(matrix.rows());
   IntRowVectorType cols_transpositions(matrix.cols());
@@ -370,7 +370,7 @@ struct ei_solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
       RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, dec().rank(), c.cols()).cwiseAbs().maxCoeff();
       RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-dec().rank(), c.cols()).cwiseAbs().maxCoeff();
       // FIXME brain dead
-      const RealScalar m_precision = epsilon<Scalar>() * std::min(rows,cols);
+      const RealScalar m_precision = NumTraits<Scalar>::epsilon() * std::min(rows,cols);
       if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
         return;
     }
diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h
index 55fac3d12..94375725f 100644
--- a/Eigen/src/SVD/JacobiSVD.h
+++ b/Eigen/src/SVD/JacobiSVD.h
@@ -295,7 +295,7 @@ JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const Ma
   int cols = matrix.cols();
   int diagSize = std::min(rows, cols);
   m_singularValues.resize(diagSize);
-  const RealScalar precision = 2 * epsilon<Scalar>();
+  const RealScalar precision = 2 * NumTraits<Scalar>::epsilon();
 
   if(!ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options>::run(matrix, work_matrix, *this)
   && !ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options>::run(matrix, work_matrix, *this))
diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h
index 1d4bbe417..c308ff3ee 100644
--- a/Eigen/src/SVD/SVD.h
+++ b/Eigen/src/SVD/SVD.h
@@ -193,7 +193,7 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
   int i=0,its=0,j=0,k=0,l=0,nm=0;
   Scalar anorm, c, f, g, h, s, scale, x, y, z;
   bool convergence = true;
-  Scalar eps = dummy_precision<Scalar>();
+  Scalar eps = NumTraits<Scalar>::dummy_precision();
 
   Matrix<Scalar,Dynamic,1> rv1(n);
   g = scale = anorm = 0;
diff --git a/Eigen/src/Sparse/AmbiVector.h b/Eigen/src/Sparse/AmbiVector.h
index 6bb6ee3e4..1ac28272b 100644
--- a/Eigen/src/Sparse/AmbiVector.h
+++ b/Eigen/src/Sparse/AmbiVector.h
@@ -296,7 +296,7 @@ class AmbiVector<_Scalar>::Iterator
       * In practice, all coefficients having a magnitude smaller than \a epsilon
       * are skipped.
       */
-    Iterator(const AmbiVector& vec, RealScalar epsilon = RealScalar(0.1)*dummy_precision<RealScalar>())
+    Iterator(const AmbiVector& vec, RealScalar epsilon = RealScalar(0.1)*NumTraits<RealScalar>::dummy_precision())
       : m_vector(vec)
     {
       m_epsilon = epsilon;
diff --git a/Eigen/src/Sparse/CompressedStorage.h b/Eigen/src/Sparse/CompressedStorage.h
index b25b05e91..4fc1797d1 100644
--- a/Eigen/src/Sparse/CompressedStorage.h
+++ b/Eigen/src/Sparse/CompressedStorage.h
@@ -185,7 +185,7 @@ class CompressedStorage
       return m_values[id];
     }
     
-    void prune(Scalar reference, RealScalar epsilon = dummy_precision<RealScalar>())
+    void prune(Scalar reference, RealScalar epsilon = NumTraits<RealScalar>::dummy_precision())
     {
       size_t k = 0;
       size_t n = size();
diff --git a/Eigen/src/Sparse/DynamicSparseMatrix.h b/Eigen/src/Sparse/DynamicSparseMatrix.h
index 96385df50..2594ffebc 100644
--- a/Eigen/src/Sparse/DynamicSparseMatrix.h
+++ b/Eigen/src/Sparse/DynamicSparseMatrix.h
@@ -209,7 +209,7 @@ class DynamicSparseMatrix
 
     inline void finalize() {}
 
-    void prune(Scalar reference, RealScalar epsilon = dummy_precision<RealScalar>())
+    void prune(Scalar reference, RealScalar epsilon = NumTraits<RealScalar>::dummy_precision())
     {
       for (int j=0; j<outerSize(); ++j)
         m_data[j].prune(reference,epsilon);
diff --git a/Eigen/src/Sparse/SparseLDLT.h b/Eigen/src/Sparse/SparseLDLT.h
index 7dfcf329c..28797a6c4 100644
--- a/Eigen/src/Sparse/SparseLDLT.h
+++ b/Eigen/src/Sparse/SparseLDLT.h
@@ -94,7 +94,7 @@ class SparseLDLT
       : m_flags(flags), m_status(0)
     {
       ei_assert((MatrixType::Flags&RowMajorBit)==0);
-      m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar>();
+      m_precision = RealScalar(0.1) * Eigen::NumTraits<RealScalar>::dummy_precision();
     }
 
     /** Creates a LDLT object and compute the respective factorization of \a matrix using
@@ -103,7 +103,7 @@ class SparseLDLT
       : m_matrix(matrix.rows(), matrix.cols()), m_flags(flags), m_status(0)
     {
       ei_assert((MatrixType::Flags&RowMajorBit)==0);
-      m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar>();
+      m_precision = RealScalar(0.1) * Eigen::NumTraits<RealScalar>::dummy_precision();
       compute(matrix);
     }
 
diff --git a/Eigen/src/Sparse/SparseLLT.h b/Eigen/src/Sparse/SparseLLT.h
index 9018fe0df..a1c10ba13 100644
--- a/Eigen/src/Sparse/SparseLLT.h
+++ b/Eigen/src/Sparse/SparseLLT.h
@@ -54,7 +54,7 @@ class SparseLLT
     SparseLLT(int flags = 0)
       : m_flags(flags), m_status(0)
     {
-      m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar>();
+      m_precision = RealScalar(0.1) * Eigen::NumTraits<RealScalar>::dummy_precision();
     }
 
     /** Creates a LLT object and compute the respective factorization of \a matrix using
@@ -62,7 +62,7 @@ class SparseLLT
     SparseLLT(const MatrixType& matrix, int flags = 0)
       : m_matrix(matrix.rows(), matrix.cols()), m_flags(flags), m_status(0)
     {
-      m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar>();
+      m_precision = RealScalar(0.1) * Eigen::NumTraits<RealScalar>::dummy_precision();
       compute(matrix);
     }
 
diff --git a/Eigen/src/Sparse/SparseLU.h b/Eigen/src/Sparse/SparseLU.h
index 0802b0807..3fe9df1ef 100644
--- a/Eigen/src/Sparse/SparseLU.h
+++ b/Eigen/src/Sparse/SparseLU.h
@@ -59,7 +59,7 @@ class SparseLU
     SparseLU(int flags = 0)
       : m_flags(flags), m_status(0)
     {
-      m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar>();
+      m_precision = RealScalar(0.1) * Eigen::NumTraits<RealScalar>::dummy_precision();
     }
 
     /** Creates a LU object and compute the respective factorization of \a matrix using
@@ -67,7 +67,7 @@ class SparseLU
     SparseLU(const MatrixType& matrix, int flags = 0)
       : /*m_matrix(matrix.rows(), matrix.cols()),*/ m_flags(flags), m_status(0)
     {
-      m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar>();
+      m_precision = RealScalar(0.1) * Eigen::NumTraits<RealScalar>::dummy_precision();
       compute(matrix);
     }
 
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h
index bba0404c3..68385f564 100644
--- a/Eigen/src/Sparse/SparseMatrix.h
+++ b/Eigen/src/Sparse/SparseMatrix.h
@@ -350,7 +350,7 @@ class SparseMatrix
       }
     }
 
-    void prune(Scalar reference, RealScalar epsilon = dummy_precision<RealScalar>())
+    void prune(Scalar reference, RealScalar epsilon = NumTraits<RealScalar>::dummy_precision())
     {
       int k = 0;
       for (int j=0; j<m_outerSize; ++j)
diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h
index 9e426c2b9..205325939 100644
--- a/Eigen/src/Sparse/SparseMatrixBase.h
+++ b/Eigen/src/Sparse/SparseMatrixBase.h
@@ -513,32 +513,32 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived
 
     template<typename OtherDerived>
     bool isApprox(const SparseMatrixBase<OtherDerived>& other,
-                  RealScalar prec = dummy_precision<Scalar>()) const
+                  RealScalar prec = NumTraits<Scalar>::dummy_precision()) const
     { return toDense().isApprox(other.toDense(),prec); }
 
     template<typename OtherDerived>
     bool isApprox(const MatrixBase<OtherDerived>& other,
-                  RealScalar prec = dummy_precision<Scalar>()) const
+                  RealScalar prec = NumTraits<Scalar>::dummy_precision()) const
     { return toDense().isApprox(other,prec); }
 //     bool isMuchSmallerThan(const RealScalar& other,
-//                            RealScalar prec = dummy_precision<Scalar>()) const;
+//                            RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
 //     template<typename OtherDerived>
 //     bool isMuchSmallerThan(const MatrixBase<OtherDerived>& other,
-//                            RealScalar prec = dummy_precision<Scalar>()) const;
+//                            RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
 
-//     bool isApproxToConstant(const Scalar& value, RealScalar prec = dummy_precision<Scalar>()) const;
-//     bool isZero(RealScalar prec = dummy_precision<Scalar>()) const;
-//     bool isOnes(RealScalar prec = dummy_precision<Scalar>()) const;
-//     bool isIdentity(RealScalar prec = dummy_precision<Scalar>()) const;
-//     bool isDiagonal(RealScalar prec = dummy_precision<Scalar>()) const;
+//     bool isApproxToConstant(const Scalar& value, RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
+//     bool isZero(RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
+//     bool isOnes(RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
+//     bool isIdentity(RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
+//     bool isDiagonal(RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
 
-//     bool isUpper(RealScalar prec = dummy_precision<Scalar>()) const;
-//     bool isLower(RealScalar prec = dummy_precision<Scalar>()) const;
+//     bool isUpper(RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
+//     bool isLower(RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
 
 //     template<typename OtherDerived>
 //     bool isOrthogonal(const MatrixBase<OtherDerived>& other,
-//                       RealScalar prec = dummy_precision<Scalar>()) const;
-//     bool isUnitary(RealScalar prec = dummy_precision<Scalar>()) const;
+//                       RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
+//     bool isUnitary(RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
 
 //     template<typename OtherDerived>
 //     inline bool operator==(const MatrixBase<OtherDerived>& other) const
diff --git a/Eigen/src/Sparse/SparseVector.h b/Eigen/src/Sparse/SparseVector.h
index 13b4f75a6..58f3ec342 100644
--- a/Eigen/src/Sparse/SparseVector.h
+++ b/Eigen/src/Sparse/SparseVector.h
@@ -202,7 +202,7 @@ class SparseVector
     EIGEN_DEPRECATED void endFill() {}
     inline void finalize() {}
 
-    void prune(Scalar reference, RealScalar epsilon = dummy_precision<RealScalar>())
+    void prune(Scalar reference, RealScalar epsilon = NumTraits<RealScalar>::dummy_precision())
     {
       m_data.prune(reference,epsilon);
     }
diff --git a/test/prec_inverse_4x4.cpp b/test/prec_inverse_4x4.cpp
index e1b05aa0d..59ccbbcaf 100644
--- a/test/prec_inverse_4x4.cpp
+++ b/test/prec_inverse_4x4.cpp
@@ -35,7 +35,7 @@ template<typename MatrixType> void inverse_permutation_4x4()
   {
     MatrixType m = PermutationMatrix<4>(indices);
     MatrixType inv = m.inverse();
-    double error = double( (m*inv-MatrixType::Identity()).norm() / epsilon<Scalar>() );
+    double error = double( (m*inv-MatrixType::Identity()).norm() / NumTraits<Scalar>::epsilon() );
     VERIFY(error == 0.0);
     std::next_permutation(indices.data(),indices.data()+4);
   }
@@ -53,9 +53,9 @@ template<typename MatrixType> void inverse_general_4x4(int repeat)
     do {
       m = MatrixType::Random();
       absdet = ei_abs(m.determinant());
-    } while(absdet < epsilon<Scalar>());
+    } while(absdet < NumTraits<Scalar>::epsilon());
     MatrixType inv = m.inverse();
-    double error = double( (m*inv-MatrixType::Identity()).norm() * absdet / epsilon<Scalar>() );
+    double error = double( (m*inv-MatrixType::Identity()).norm() * absdet / NumTraits<Scalar>::epsilon() );
     error_sum += error;
     error_max = std::max(error_max, error);
   }
diff --git a/test/product.h b/test/product.h
index 85afbe1e4..f6109fae4 100644
--- a/test/product.h
+++ b/test/product.h
@@ -26,7 +26,7 @@
 #include <Eigen/QR>
 
 template<typename Derived1, typename Derived2>
-bool areNotApprox(const MatrixBase<Derived1>& m1, const MatrixBase<Derived2>& m2, typename Derived1::RealScalar epsilon = dummy_precision<typename Derived1::RealScalar>())
+bool areNotApprox(const MatrixBase<Derived1>& m1, const MatrixBase<Derived2>& m2, typename Derived1::RealScalar epsilon = NumTraits<typename Derived1::RealScalar>::dummy_precision())
 {
   return !((m1-m2).cwiseAbs2().maxCoeff() < epsilon * epsilon
                           * std::max(m1.cwiseAbs2().maxCoeff(), m2.cwiseAbs2().maxCoeff()));