diff --git a/doc/example.cpp b/doc/example.cpp
index ca10d077a..b5582f006 100644
--- a/doc/example.cpp
+++ b/doc/example.cpp
@@ -11,7 +11,7 @@ void foo(const Eigen::Object<Scalar, Derived>& m)
 }
 
 template<typename Scalar, typename Derived>
-Eigen::ScalarProduct<Derived>
+Eigen::ScalarMultiple<Derived>
 twice(const Eigen::Object<Scalar, Derived>& m)
 {
 	return 2 * m;
diff --git a/src/Core/Dot.h b/src/Core/Dot.h
index 158b852eb..3bf8514ba 100644
--- a/src/Core/Dot.h
+++ b/src/Core/Dot.h
@@ -88,7 +88,7 @@ typename NumTraits<Scalar>::Real Object<Scalar, Derived>::norm() const
 }
 
 template<typename Scalar, typename Derived>
-ScalarProduct<Derived> Object<Scalar, Derived>::normalized() const
+ScalarMultiple<Derived> Object<Scalar, Derived>::normalized() const
 {
   return (*this) / norm();
 }
diff --git a/src/Core/Fuzzy.h b/src/Core/Fuzzy.h
index 0334e6d96..4980de8d1 100644
--- a/src/Core/Fuzzy.h
+++ b/src/Core/Fuzzy.h
@@ -28,34 +28,38 @@
 
 template<typename Scalar, typename Derived>
 template<typename OtherDerived>
-bool Object<Scalar, Derived>::isApprox(const OtherDerived& other) const
+bool Object<Scalar, Derived>::isApprox(
+  const OtherDerived& other,
+  const typename NumTraits<Scalar>::Real& prec
+) const
 {
   if(IsVector)
   {
-    return((*this - other).norm2()
-           <= std::min(norm2(), other.norm2())
-            * NumTraits<Scalar>::epsilon2());
+    return((*this - other).norm2() <= std::min(norm2(), other.norm2()) * prec * prec);
   }
   else
   {
     for(int i = 0; i < cols(); i++)
-      if(!col(i).isApprox(other.col(i)))
+      if(!col(i).isApprox(other.col(i), prec))
         return false;
     return true;
   }
 }
 
 template<typename Scalar, typename Derived>
-bool Object<Scalar, Derived>::isNegligble(const Scalar& other) const
+bool Object<Scalar, Derived>::isMuchSmallerThan(
+  const Scalar& other,
+  const typename NumTraits<Scalar>::Real& prec
+) const
 {
   if(IsVector)
   {
-    return(norm2() <= NumTraits<Scalar>::abs2(other) * NumTraits<Scalar>::epsilon2());
+    return(norm2() <= NumTraits<Scalar>::abs2(other) * prec * prec);
   }
   else
   {
     for(int i = 0; i < cols(); i++)
-      if(!col(i).isNegligible(other))
+      if(!col(i).isMuchSmallerThan(other, prec))
         return false;
     return true;
   }
@@ -63,16 +67,19 @@ bool Object<Scalar, Derived>::isNegligble(const Scalar& other) const
 
 template<typename Scalar, typename Derived>
 template<typename OtherDerived>
-bool Object<Scalar, Derived>::isNegligble(const Object<Scalar, OtherDerived>& other) const
+bool Object<Scalar, Derived>::isMuchSmallerThan(
+  const Object<Scalar, OtherDerived>& other,
+  const typename NumTraits<Scalar>::Real& prec
+) const
 {
   if(IsVector)
   {
-    return(norm2() <= other.norm2() * NumTraits<Scalar>::epsilon2());
+    return(norm2() <= other.norm2() * prec * prec);
   }
   else
   {
     for(int i = 0; i < cols(); i++)
-      if(!col(i).isNegligible(other.col(i)))
+      if(!col(i).isMuchSmallerThan(other.col(i), prec))
         return false;
     return true;
   }
diff --git a/src/Core/MatrixOps.h b/src/Core/MatrixOps.h
index a1a1f47a0..7cec76407 100644
--- a/src/Core/MatrixOps.h
+++ b/src/Core/MatrixOps.h
@@ -206,7 +206,7 @@ operator-(const Object<Scalar, Derived1> &mat1, const Object<Scalar, Derived2> &
 template<typename Scalar, typename Derived>
 template<typename OtherDerived>
 MatrixProduct<Derived, OtherDerived>
-Object<Scalar, Derived>::lazyMul(const Object<Scalar, OtherDerived> &other) const
+Object<Scalar, Derived>::lazyProduct(const Object<Scalar, OtherDerived> &other) const
 {
   return MatrixProduct<Derived, OtherDerived>(constRef(), other.constRef());
 }
@@ -215,7 +215,7 @@ template<typename Scalar, typename Derived1, typename Derived2>
 Eval<MatrixProduct<Derived1, Derived2> >
 operator*(const Object<Scalar, Derived1> &mat1, const Object<Scalar, Derived2> &mat2)
 {
-  return mat1.lazyMul(mat2).eval();
+  return mat1.lazyProduct(mat2).eval();
 }
 
 template<typename Scalar, typename Derived>
diff --git a/src/Core/Numeric.h b/src/Core/Numeric.h
index 8437dfb2e..04721e019 100644
--- a/src/Core/Numeric.h
+++ b/src/Core/Numeric.h
@@ -38,8 +38,7 @@ template<> struct NumTraits<int>
   static const bool IsComplex = false;
   static const bool HasFloatingPoint = false;
   
-  static int epsilon() { return 0; }
-  static int epsilon2() { return 0; }
+  static int precision() { return 0; }
   static int real(const int& x) { return x; }
   static int imag(const int& x) { EI_UNUSED(x); return 0; }
   static int conj(const int& x) { return x; }
@@ -53,18 +52,21 @@ template<> struct NumTraits<int>
     // as well, so there's no problem.
     return (std::rand() % 21) - 10;
   }
-  static bool negligible(const int& a, const int& b)
+  static bool isMuchSmallerThan(const int& a, const int& b, const int& prec = precision())
   {
     EI_UNUSED(b);
-    return(a == 0);
+    EI_UNUSED(prec);
+    return a == 0;
   }
-  static bool approx(const int& a, const int& b)
+  static bool isApprox(const int& a, const int& b, const int& prec = precision())
   {
-    return(a == b);
+    EI_UNUSED(prec);
+    return a == b;
   }
-  static bool lessThanOrApprox(const int& a, const int& b)
+  static bool isApproxOrLessThan(const int& a, const int& b, const int& prec = precision())
   {
-    return(a <= b);
+    EI_UNUSED(prec);
+    return a <= b;
   }
 };
 
@@ -77,8 +79,7 @@ template<> struct NumTraits<float>
   static const bool IsComplex = false;
   static const bool HasFloatingPoint = true;
   
-  static float epsilon() { return 1e-5f; }
-  static float epsilon2() { return epsilon() * epsilon(); }
+  static float precision() { return 1e-5f; }
   static float real(const float& x) { return x; }
   static float imag(const float& x) { EI_UNUSED(x); return 0; }
   static float conj(const float& x) { return x; }
@@ -89,17 +90,17 @@ template<> struct NumTraits<float>
   {
     return std::rand() / (RAND_MAX/20.0f) - 10.0f;
   }
-  static bool negligible(const float& a, const float& b)
+  static bool isMuchSmallerThan(const float& a, const float& b, const float& prec = precision())
   {
-    return(abs(a) <= abs(b) * epsilon());
+    return abs(a) <= abs(b) * prec;
   }
-  static bool approx(const float& a, const float& b)
+  static bool isApprox(const float& a, const float& b, const float& prec = precision())
   {
-    return(abs(a - b) <= std::min(abs(a), abs(b)) * epsilon());
+    return abs(a - b) <= std::min(abs(a), abs(b)) * prec;
   }
-  static bool lessThanOrApprox(const float& a, const float& b)
+  static bool isApproxOrLessThan(const float& a, const float& b, const float& prec = precision())
   {
-    return(a <= b || approx(a, b));
+    return a <= b || isApprox(a, b, prec);
   }
 };
 
@@ -112,8 +113,7 @@ template<> struct NumTraits<double>
   static const bool IsComplex = false;
   static const bool HasFloatingPoint = true;
   
-  static double epsilon() { return 1e-11; }
-  static double epsilon2() { return epsilon() * epsilon(); }
+  static double precision() { return 1e-11; }
   static double real(const double& x) { return x; }
   static double imag(const double& x) { EI_UNUSED(x); return 0; }
   static double conj(const double& x) { return x; }
@@ -124,17 +124,17 @@ template<> struct NumTraits<double>
   {
     return std::rand() / (RAND_MAX/20.0) - 10.0;
   }
-  static bool negligible(const double& a, const double& b)
+  static bool isMuchSmallerThan(const double& a, const double& b, const double& prec = precision())
   {
-    return(abs(a) <= abs(b) * epsilon());
+    return abs(a) <= abs(b) * prec;
   }
-  static bool approx(const double& a, const double& b)
+  static bool isApprox(const double& a, const double& b, const double& prec = precision())
   {
-    return(abs(a - b) <= std::min(abs(a), abs(b)) * epsilon());
+    return abs(a - b) <= std::min(abs(a), abs(b)) * prec;
   }
-  static bool lessThanOrApprox(const double& a, const double& b)
+  static bool isApproxOrLessThan(const double& a, const double& b, const double& prec = precision())
   {
-    return(a <= b || approx(a, b));
+    return a <= b || isApprox(a, b, prec);
   }
 };
 
@@ -148,8 +148,7 @@ template<typename _Real> struct NumTraits<std::complex<_Real> >
   static const bool IsComplex = true;
   static const bool HasFloatingPoint = NumTraits<Real>::HasFloatingPoint;
   
-  static Real epsilon() { return NumTraits<Real>::epsilon(); }
-  static Real epsilon2() { return epsilon() * epsilon(); }
+  static Real precision() { return NumTraits<Real>::precision(); }
   static Real real(const Complex& x) { return std::real(x); }
   static Real imag(const Complex& x) { return std::imag(x); }
   static Complex conj(const Complex& x) { return std::conj(x); }
@@ -163,16 +162,16 @@ template<typename _Real> struct NumTraits<std::complex<_Real> >
   {
     return Complex(NumTraits<Real>::rand(), NumTraits<Real>::rand());
   }
-  static bool negligible(const Complex& a, const Complex& b)
+  static bool isMuchSmallerThan(const Complex& a, const Complex& b, const Real& prec = precision())
   {
-    return(abs2(a) <= abs2(b) * epsilon2());
+    return abs2(a) <= abs2(b) * prec * prec;
   }
-  static bool approx(const Complex& a, const Complex& b)
+  static bool isApprox(const Complex& a, const Complex& b, const Real& prec = precision())
   {
-    return(NumTraits<Real>::approx(std::real(a), std::real(b))
-        && NumTraits<Real>::approx(std::imag(a), std::imag(b)));
+    return NumTraits<Real>::isApprox(std::real(a), std::real(b), prec)
+        && NumTraits<Real>::isApprox(std::imag(a), std::imag(b), prec);
   }
-  // lessThanOrApprox wouldn't make sense for complex numbers
+  // isApproxOrLessThan wouldn't make sense for complex numbers
 };
 
 #endif // EI_NUMERIC_H
diff --git a/src/Core/Object.h b/src/Core/Object.h
index 27e5e456a..90962c46c 100644
--- a/src/Core/Object.h
+++ b/src/Core/Object.h
@@ -95,19 +95,29 @@ template<typename Scalar, typename Derived> class Object
     
     RealScalar norm2() const;
     RealScalar norm()  const;
-    ScalarProduct<Derived> normalized() const;
+    ScalarMultiple<Derived> normalized() const;
     
     static Eval<Random<Derived> >
     random(int rows = RowsAtCompileTime, int cols = ColsAtCompileTime);
     
     template<typename OtherDerived>
-    bool isApprox(const OtherDerived& other) const;
+    bool isApprox(
+      const OtherDerived& other,
+      const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::precision()
+    ) const;
+    bool isMuchSmallerThan(
+      const Scalar& other,
+      const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::precision()
+    ) const;
     template<typename OtherDerived>
-    bool isNegligible(const OtherDerived& other) const;
+    bool isMuchSmallerThan(
+      const OtherDerived& other,
+      const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::precision()
+    ) const;
     
     template<typename OtherDerived>
     MatrixProduct<Derived, OtherDerived>
-    lazyMul(const Object<Scalar, OtherDerived>& other) const EI_ALWAYS_INLINE;
+    lazyProduct(const Object<Scalar, OtherDerived>& other) const EI_ALWAYS_INLINE;
     
     template<typename OtherDerived>
     Derived& operator+=(const Object<Scalar, OtherDerived>& other);
diff --git a/src/Core/ScalarOps.h b/src/Core/ScalarOps.h
index 4a71b403e..77c9b9342 100644
--- a/src/Core/ScalarOps.h
+++ b/src/Core/ScalarOps.h
@@ -26,28 +26,28 @@
 #ifndef EI_SCALAROPS_H
 #define EI_SCALAROPS_H
 
-template<typename MatrixType> class ScalarProduct
-  : public Object<typename MatrixType::Scalar, ScalarProduct<MatrixType> >
+template<typename MatrixType> class ScalarMultiple
+  : public Object<typename MatrixType::Scalar, ScalarMultiple<MatrixType> >
 {
   public:
     typedef typename MatrixType::Scalar Scalar;
     typedef typename MatrixType::ConstRef MatRef;
-    friend class Object<typename MatrixType::Scalar, ScalarProduct<MatrixType> >;
+    friend class Object<typename MatrixType::Scalar, ScalarMultiple<MatrixType> >;
 
     static const int RowsAtCompileTime = MatrixType::RowsAtCompileTime,
                      ColsAtCompileTime = MatrixType::ColsAtCompileTime;
 
-    ScalarProduct(const MatRef& matrix, Scalar scalar)
+    ScalarMultiple(const MatRef& matrix, Scalar scalar)
       : m_matrix(matrix), m_scalar(scalar) {}
 
-    ScalarProduct(const ScalarProduct& other)
+    ScalarMultiple(const ScalarMultiple& other)
       : m_matrix(other.m_matrix), m_scalar(other.m_scalar) {}
 
-    EI_INHERIT_ASSIGNMENT_OPERATORS(ScalarProduct)
+    EI_INHERIT_ASSIGNMENT_OPERATORS(ScalarMultiple)
 
   private:
-    const ScalarProduct& _ref() const { return *this; }
-    const ScalarProduct& _constRef() const { return *this; }
+    const ScalarMultiple& _ref() const { return *this; }
+    const ScalarMultiple& _constRef() const { return *this; }
     int _rows() const { return m_matrix.rows(); }
     int _cols() const { return m_matrix.cols(); }
 
@@ -63,23 +63,23 @@ template<typename MatrixType> class ScalarProduct
 
 #define EI_MAKE_SCALAR_OPS(OtherScalar)                                \
 template<typename Scalar, typename Derived>                            \
-ScalarProduct<Derived>                                               \
+ScalarMultiple<Derived>                                               \
 operator*(const Object<Scalar, Derived>& matrix,                     \
           OtherScalar scalar)                                          \
 {                                                                      \
-  return ScalarProduct<Derived>(matrix.constRef(), scalar);          \
+  return ScalarMultiple<Derived>(matrix.constRef(), scalar);          \
 }                                                                      \
                                                                        \
 template<typename Scalar, typename Derived>                            \
-ScalarProduct<Derived>                                               \
+ScalarMultiple<Derived>                                               \
 operator*(OtherScalar scalar,                                          \
           const Object<Scalar, Derived>& matrix)                     \
 {                                                                      \
-  return ScalarProduct<Derived>(matrix.constRef(), scalar);          \
+  return ScalarMultiple<Derived>(matrix.constRef(), scalar);          \
 }                                                                      \
                                                                        \
 template<typename Scalar, typename Derived>                            \
-ScalarProduct<Derived>                                               \
+ScalarMultiple<Derived>                                               \
 operator/(const Object<Scalar, Derived>& matrix,                     \
           OtherScalar scalar)                                          \
 {                                                                      \
diff --git a/src/Core/Util.h b/src/Core/Util.h
index 54edeed50..431c008ed 100644
--- a/src/Core/Util.h
+++ b/src/Core/Util.h
@@ -53,7 +53,7 @@ template<typename MatrixType> class Conjugate;
 template<typename Lhs, typename Rhs> class Sum;
 template<typename Lhs, typename Rhs> class Difference;
 template<typename Lhs, typename Rhs> class MatrixProduct;
-template<typename MatrixType> class ScalarProduct;
+template<typename MatrixType> class ScalarMultiple;
 template<typename MatrixType> class Random;
 template<typename ExpressionType> class Eval;
 
diff --git a/test/main.h b/test/main.h
index 85968134c..fcb1b23ab 100644
--- a/test/main.h
+++ b/test/main.h
@@ -57,32 +57,47 @@ template<> inline int TestEpsilon<std::complex<int> >() { return TestEpsilon<int
 template<> inline float TestEpsilon<std::complex<float> >() { return TestEpsilon<float>(); }
 template<> inline double TestEpsilon<std::complex<double> >() { return TestEpsilon<double>(); }
 
-template<typename T> bool TestNegligible(const T& a, const T& b)
+template<typename T> bool TestMuchSmallerThan(const T& a, const T& b)
 {
-  return(Abs(a) <= Abs(b) * TestEpsilon<T>());
+  return NumTraits<T>::isMuchSmallerThan(a, b, TestEpsilon<T>());
 }
 
-//template<typename Scalar, typename Derived, typename OtherDerived>
-//bool TestNegligible
+template<typename Scalar, typename Derived, typename OtherDerived>
+bool TestMuchSmallerThan(
+  const Object<Scalar, Derived>& a,
+  const Object<Scalar, OtherDerived>& b)
+{
+  return a.isMuchSmallerThan(b, TestEpsilon<Scalar>());
+}
 
 template<typename T> bool TestApprox(const T& a, const T& b)
 {
-  if(Eigen::NumTraits<T>::IsFloat)
-    return(Abs(a - b) <= std::min(Abs(a), Abs(b)) * TestEpsilon<T>());
-  else
-    return(a == b);
+  return NumTraits<T>::isApprox(a, b, TestEpsilon<T>());
 }
 
-template<typename T> bool TestLessThanOrApprox(const T& a, const T& b)
+template<typename Scalar, typename Derived, typename OtherDerived>
+bool TestApprox(
+  const Object<Scalar, Derived>& a,
+  const Object<Scalar, OtherDerived>& b)
 {
-  if(Eigen::NumTraits<T>::IsFloat)
-    return(a < b || Approx(a, b));
-  else
-    return(a <= b);
+  return a.isApprox(b, TestEpsilon<Scalar>());
 }
 
-#define QVERIFY_NEGLIGIBLE(a, b) QVERIFY(TestNegligible(a, b))
+template<typename T> bool TestApproxOrLessThan(const T& a, const T& b)
+{
+  return NumTraits<T>::isApproxOrLessThan(a, b, TestEpsilon<T>());
+}
+
+template<typename Scalar, typename Derived, typename OtherDerived>
+bool TestApproxOrLessThan(
+  const Object<Scalar, Derived>& a,
+  const Object<Scalar, OtherDerived>& b)
+{
+  return a.isApproxOrLessThan(b, TestEpsilon<Scalar>());
+}
+
+#define QVERIFY_MUCH_SMALLER_THAN(a, b) QVERIFY(TestMuchSmallerThan(a, b))
 #define QVERIFY_APPROX(a, b) QVERIFY(TestApprox(a, b))
-#define QVERIFY_LESS_THAN_OR_APPROX(a, b) QVERIFY(TestLessThanOrApprox(a, b))
+#define QVERIFY_APPROX_OR_LESS_THAN(a, b) QVERIFY(TestApproxOrLessThan(a, b))
 
 #endif // EI_TEST_MAIN_H
diff --git a/test/matrixops.cpp b/test/matrixops.cpp
index c763ff687..9c0933d89 100644
--- a/test/matrixops.cpp
+++ b/test/matrixops.cpp
@@ -51,7 +51,7 @@ template<typename MatrixType1,
   if(rows1 == cols1)
   {
     a *= b;
-    a.lazyMul(b);
+    a.lazyProduct(b);
   }
   
   MatrixType1 d(rows1, cols1);