diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index ef1c11f62..1d0369b65 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -440,7 +440,7 @@ template<> struct ldlt_inplace<Lower>
       // Update the terms of L
       Index rs = size-j-1;
       w.tail(rs) -= wj * mat.col(j).tail(rs);
-      if(gamma != 0)
+      if(!numext::is_exactly_zero(gamma))
         mat.col(j).tail(rs) += (sigma*numext::conj(wj)/gamma)*w.tail(rs);
     }
     return true;
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index c223da266..ee590ebe7 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -297,7 +297,7 @@ static Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const
       if(rs)
       {
         temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
-        if(gamma != 0)
+        if(!numext::is_exactly_zero(gamma))
           mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*numext::conj(wj)/gamma)*temp.tail(rs);
       }
     }
diff --git a/Eigen/src/Core/ConditionEstimator.h b/Eigen/src/Core/ConditionEstimator.h
index bd4455f41..694be8ba7 100644
--- a/Eigen/src/Core/ConditionEstimator.h
+++ b/Eigen/src/Core/ConditionEstimator.h
@@ -162,12 +162,12 @@ rcond_estimate_helper(typename Decomposition::RealScalar matrix_norm, const Deco
 {
   typedef typename Decomposition::RealScalar RealScalar;
   eigen_assert(dec.rows() == dec.cols());
-  if (dec.rows() == 0)              return NumTraits<RealScalar>::infinity();
-  if (matrix_norm == RealScalar(0)) return RealScalar(0);
-  if (dec.rows() == 1)              return RealScalar(1);
+  if (dec.rows() == 0)                        return NumTraits<RealScalar>::infinity();
+  if (numext::is_exactly_zero(matrix_norm)) return RealScalar(0);
+  if (dec.rows() == 1)                        return RealScalar(1);
   const RealScalar inverse_matrix_norm = rcond_invmatrix_L1_norm_estimate(dec);
-  return (inverse_matrix_norm == RealScalar(0) ? RealScalar(0)
-                                               : (RealScalar(1) / inverse_matrix_norm) / matrix_norm);
+  return (numext::is_exactly_zero(inverse_matrix_norm) ? RealScalar(0)
+                                                       : (RealScalar(1) / inverse_matrix_norm) / matrix_norm);
 }
 
 }  // namespace internal
diff --git a/Eigen/src/Core/GeneralProduct.h b/Eigen/src/Core/GeneralProduct.h
index 0d24a00af..783a3b6c4 100644
--- a/Eigen/src/Core/GeneralProduct.h
+++ b/Eigen/src/Core/GeneralProduct.h
@@ -219,7 +219,6 @@ template<> struct gemv_dense_selector<OnTheRight,ColMajor,true>
     typedef typename Lhs::Scalar   LhsScalar;
     typedef typename Rhs::Scalar   RhsScalar;
     typedef typename Dest::Scalar  ResScalar;
-    typedef typename Dest::RealScalar  RealScalar;
     
     typedef internal::blas_traits<Lhs> LhsBlasTraits;
     typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType;
@@ -264,7 +263,7 @@ template<> struct gemv_dense_selector<OnTheRight,ColMajor,true>
     {
       gemv_static_vector_if<ResScalar,ActualDest::SizeAtCompileTime,ActualDest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest;
 
-      const bool alphaIsCompatible = (!ComplexByReal) || (numext::imag(actualAlpha)==RealScalar(0));
+      const bool alphaIsCompatible = (!ComplexByReal) || (numext::is_exactly_zero(numext::imag(actualAlpha)));
       const bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible;
 
       ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(),
diff --git a/Eigen/src/Core/MathFunctionsImpl.h b/Eigen/src/Core/MathFunctionsImpl.h
index 182dd371e..6b248d59c 100644
--- a/Eigen/src/Core/MathFunctionsImpl.h
+++ b/Eigen/src/Core/MathFunctionsImpl.h
@@ -119,7 +119,7 @@ RealScalar positive_real_hypot(const RealScalar& x, const RealScalar& y)
   EIGEN_USING_STD(sqrt);
   RealScalar p, qp;
   p = numext::maxi(x,y);
-  if(p==RealScalar(0)) return RealScalar(0);
+  if(numext::is_exactly_zero(p)) return RealScalar(0);
   qp = numext::mini(y,x) / p;
   return p * sqrt(RealScalar(1) + qp*qp);
 }
@@ -169,8 +169,8 @@ EIGEN_DEVICE_FUNC std::complex<T> complex_sqrt(const std::complex<T>& z) {
 
   return
     (numext::isinf)(y) ? std::complex<T>(NumTraits<T>::infinity(), y)
-      : x == zero ? std::complex<T>(w, y < zero ? -w : w)
-      : x > zero ? std::complex<T>(w, y / (2 * w))
+      : numext::is_exactly_zero(x) ? std::complex<T>(w, y < zero ? -w : w)
+                                   : x > zero ? std::complex<T>(w, y / (2 * w))
       : std::complex<T>(numext::abs(y) / (2 * w), y < zero ? -w : w );
 }
 
@@ -208,10 +208,10 @@ EIGEN_DEVICE_FUNC std::complex<T> complex_rsqrt(const std::complex<T>& z) {
   const T woz = w / abs_z;
   // Corner cases consistent with 1/sqrt(z) on gcc/clang.
   return
-    abs_z == zero ? std::complex<T>(NumTraits<T>::infinity(), NumTraits<T>::quiet_NaN())
-      : ((numext::isinf)(x) || (numext::isinf)(y)) ? std::complex<T>(zero, zero)
-      : x == zero ? std::complex<T>(woz, y < zero ? woz : -woz)
-      : x > zero ? std::complex<T>(woz, -y / (2 * w * abs_z))
+          numext::is_exactly_zero(abs_z) ? std::complex<T>(NumTraits<T>::infinity(), NumTraits<T>::quiet_NaN())
+                                         : ((numext::isinf)(x) || (numext::isinf)(y)) ? std::complex<T>(zero, zero)
+      : numext::is_exactly_zero(x) ? std::complex<T>(woz, y < zero ? woz : -woz)
+                                   : x > zero ? std::complex<T>(woz, -y / (2 * w * abs_z))
       : std::complex<T>(numext::abs(y) / (2 * w * abs_z), y < zero ? woz : -woz );
 }
 
diff --git a/Eigen/src/Core/ProductEvaluators.h b/Eigen/src/Core/ProductEvaluators.h
index 5f9ffe86b..a874ee326 100644
--- a/Eigen/src/Core/ProductEvaluators.h
+++ b/Eigen/src/Core/ProductEvaluators.h
@@ -460,7 +460,7 @@ protected:
   void eval_dynamic_impl(Dst& dst, const LhsT& lhs, const RhsT& rhs, const Func &func, const Scalar&  s /* == 1 */, false_type)
   {
     EIGEN_UNUSED_VARIABLE(s);
-    eigen_internal_assert(s==Scalar(1));
+    eigen_internal_assert(numext::is_exactly_one(s));
     call_restricted_packet_assignment_no_alias(dst, lhs.lazyProduct(rhs), func);
   }
 
diff --git a/Eigen/src/Core/products/TriangularMatrixMatrix.h b/Eigen/src/Core/products/TriangularMatrixMatrix.h
index 5b8ca1232..5f65798ba 100644
--- a/Eigen/src/Core/products/TriangularMatrixMatrix.h
+++ b/Eigen/src/Core/products/TriangularMatrixMatrix.h
@@ -453,12 +453,12 @@ struct triangular_product_impl<Mode,LhsIsTriangular,Lhs,false,Rhs,false>
     // Apply correction if the diagonal is unit and a scalar factor was nested:
     if ((Mode&UnitDiag)==UnitDiag)
     {
-      if (LhsIsTriangular && lhs_alpha!=LhsScalar(1))
+      if (LhsIsTriangular && !numext::is_exactly_one(lhs_alpha))
       {
         Index diagSize = (std::min)(lhs.rows(),lhs.cols());
         dst.topRows(diagSize) -= ((lhs_alpha-LhsScalar(1))*a_rhs).topRows(diagSize);
       }
-      else if ((!LhsIsTriangular) && rhs_alpha!=RhsScalar(1))
+      else if ((!LhsIsTriangular) && !numext::is_exactly_one(rhs_alpha))
       {
         Index diagSize = (std::min)(rhs.rows(),rhs.cols());
         dst.leftCols(diagSize) -= (rhs_alpha-RhsScalar(1))*a_lhs.leftCols(diagSize);
diff --git a/Eigen/src/Core/products/TriangularMatrixVector.h b/Eigen/src/Core/products/TriangularMatrixVector.h
index c6d5afa10..72b3c13a0 100644
--- a/Eigen/src/Core/products/TriangularMatrixVector.h
+++ b/Eigen/src/Core/products/TriangularMatrixVector.h
@@ -211,7 +211,6 @@ template<int Mode> struct trmv_selector<Mode,ColMajor>
     typedef typename Lhs::Scalar      LhsScalar;
     typedef typename Rhs::Scalar      RhsScalar;
     typedef typename Dest::Scalar     ResScalar;
-    typedef typename Dest::RealScalar RealScalar;
     
     typedef internal::blas_traits<Lhs> LhsBlasTraits;
     typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType;
@@ -237,7 +236,7 @@ template<int Mode> struct trmv_selector<Mode,ColMajor>
 
     gemv_static_vector_if<ResScalar,Dest::SizeAtCompileTime,Dest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest;
 
-    bool alphaIsCompatible = (!ComplexByReal) || (numext::imag(actualAlpha)==RealScalar(0));
+    bool alphaIsCompatible = (!ComplexByReal) || numext::is_exactly_zero(numext::imag(actualAlpha));
     bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible;
 
     RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha);
@@ -278,7 +277,7 @@ template<int Mode> struct trmv_selector<Mode,ColMajor>
         dest = MappedDest(actualDestPtr, dest.size());
     }
 
-    if ( ((Mode&UnitDiag)==UnitDiag) && (lhs_alpha!=LhsScalar(1)) )
+    if ( ((Mode&UnitDiag)==UnitDiag) && !numext::is_exactly_one(lhs_alpha) )
     {
       Index diagSize = (std::min)(lhs.rows(),lhs.cols());
       dest.head(diagSize) -= (lhs_alpha-LhsScalar(1))*rhs.head(diagSize);
@@ -337,7 +336,7 @@ template<int Mode> struct trmv_selector<Mode,RowMajor>
             dest.data(),dest.innerStride(),
             actualAlpha);
 
-    if ( ((Mode&UnitDiag)==UnitDiag) && (lhs_alpha!=LhsScalar(1)) )
+    if ( ((Mode&UnitDiag)==UnitDiag) && !numext::is_exactly_one(lhs_alpha) )
     {
       Index diagSize = (std::min)(lhs.rows(),lhs.cols());
       dest.head(diagSize) -= (lhs_alpha-LhsScalar(1))*rhs.head(diagSize);
diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h
index 67ec8c939..0646a9a91 100755
--- a/Eigen/src/Core/util/Meta.h
+++ b/Eigen/src/Core/util/Meta.h
@@ -397,6 +397,8 @@ struct aligned_storage {
 
 } // end namespace internal
 
+template<typename T> struct NumTraits;
+
 namespace numext {
 
 #if defined(EIGEN_GPU_COMPILE_PHASE)
@@ -429,6 +431,20 @@ template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC
 bool equal_strict(const double& x,const double& y) { return std::equal_to<double>()(x,y); }
 #endif
 
+/**
+ * \internal Performs an exact comparison of x to zero, e.g. to decide whether a term can be ignored.
+ * Use this to to bypass -Wfloat-equal warnings when exact zero is what needs to be tested.
+*/
+template<typename X> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC
+bool is_exactly_zero(const X& x) { return equal_strict(x, typename NumTraits<X>::Literal{0}); }
+
+/**
+ * \internal Performs an exact comparison of x to one, e.g. to decide whether a factor needs to be multiplied.
+ * Use this to to bypass -Wfloat-equal warnings when exact one is what needs to be tested.
+*/
+template<typename X> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC
+bool is_exactly_one(const X& x) { return equal_strict(x, typename NumTraits<X>::Literal{1}); }
+
 template<typename X, typename Y> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC
 bool not_equal_strict(const X& x,const Y& y) { return x != y; }
 
diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h
index b4f82492e..80a28fb39 100644
--- a/Eigen/src/Eigenvalues/ComplexSchur.h
+++ b/Eigen/src/Eigenvalues/ComplexSchur.h
@@ -308,7 +308,7 @@ typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::compu
   // In this case, det==0, and all we have to do is checking that eival2_norm!=0
   if(eival1_norm > eival2_norm)
     eival2 = det / eival1;
-  else if(eival2_norm!=RealScalar(0))
+  else if(!numext::is_exactly_zero(eival2_norm))
     eival1 = det / eival2;
 
   // choose the eigenvalue closest to the bottom entry of the diagonal
diff --git a/Eigen/src/Eigenvalues/RealQZ.h b/Eigen/src/Eigenvalues/RealQZ.h
index 5564f7fc4..545918fa8 100644
--- a/Eigen/src/Eigenvalues/RealQZ.h
+++ b/Eigen/src/Eigenvalues/RealQZ.h
@@ -239,7 +239,7 @@ namespace Eigen {
         for (Index i=dim-1; i>=j+2; i--) {
           JRs G;
           // kill S(i,j)
-          if(m_S.coeff(i,j) != 0)
+          if(!numext::is_exactly_zero(m_S.coeff(i, j)))
           {
             G.makeGivens(m_S.coeff(i-1,j), m_S.coeff(i,j), &m_S.coeffRef(i-1, j));
             m_S.coeffRef(i,j) = Scalar(0.0);
@@ -250,7 +250,7 @@ namespace Eigen {
               m_Q.applyOnTheRight(i-1,i,G);
           }
           // kill T(i,i-1)
-          if(m_T.coeff(i,i-1)!=Scalar(0))
+          if(!numext::is_exactly_zero(m_T.coeff(i, i - 1)))
           {
             G.makeGivens(m_T.coeff(i,i), m_T.coeff(i,i-1), &m_T.coeffRef(i,i));
             m_T.coeffRef(i,i-1) = Scalar(0.0);
@@ -288,7 +288,7 @@ namespace Eigen {
       while (res > 0)
       {
         Scalar s = abs(m_S.coeff(res-1,res-1)) + abs(m_S.coeff(res,res));
-        if (s == Scalar(0.0))
+        if (numext::is_exactly_zero(s))
           s = m_normOfS;
         if (abs(m_S.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
           break;
@@ -318,7 +318,7 @@ namespace Eigen {
       using std::abs;
       using std::sqrt;
       const Index dim=m_S.cols();
-      if (abs(m_S.coeff(i+1,i))==Scalar(0))
+      if (numext::is_exactly_zero(abs(m_S.coeff(i + 1, i))))
         return;
       Index j = findSmallDiagEntry(i,i+1);
       if (j==i-1)
@@ -629,7 +629,7 @@ namespace Eigen {
       {
         for(Index i=0; i<dim-1; ++i)
         {
-          if(m_S.coeff(i+1, i) != Scalar(0))
+          if(!numext::is_exactly_zero(m_S.coeff(i + 1, i)))
           {
             JacobiRotation<Scalar> j_left, j_right;
             internal::real_2x2_jacobi_svd(m_T, i, i+1, &j_left, &j_right);
diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h
index 96372dcb5..98176668a 100644
--- a/Eigen/src/Eigenvalues/RealSchur.h
+++ b/Eigen/src/Eigenvalues/RealSchur.h
@@ -314,7 +314,7 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMa
   Scalar considerAsZero = numext::maxi<Scalar>( norm * numext::abs2(NumTraits<Scalar>::epsilon()),
                                                 (std::numeric_limits<Scalar>::min)() );
 
-  if(norm!=Scalar(0))
+  if(!numext::is_exactly_zero(norm))
   {
     while (iu >= 0)
     {
@@ -517,7 +517,7 @@ inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Inde
     Matrix<Scalar, 2, 1> ess;
     v.makeHouseholder(ess, tau, beta);
     
-    if (beta != Scalar(0)) // if v is not zero
+    if (!numext::is_exactly_zero(beta)) // if v is not zero
     {
       if (firstIteration && k > il)
         m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1);
@@ -537,7 +537,7 @@ inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Inde
   Matrix<Scalar, 1, 1> ess;
   v.makeHouseholder(ess, tau, beta);
 
-  if (beta != Scalar(0)) // if v is not zero
+  if (!numext::is_exactly_zero(beta)) // if v is not zero
   {
     m_matT.coeffRef(iu-1, iu-2) = beta;
     m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
index 70d370c1b..d196ec054 100644
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -447,7 +447,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
   // map the matrix coefficients to [-1:1] to avoid over- and underflow.
   mat = matrix.template triangularView<Lower>();
   RealScalar scale = mat.cwiseAbs().maxCoeff();
-  if(scale==RealScalar(0)) scale = RealScalar(1);
+  if(numext::is_exactly_zero(scale)) scale = RealScalar(1);
   mat.template triangularView<Lower>() /= scale;
   m_subdiag.resize(n-1);
   m_hcoeffs.resize(n-1);
@@ -526,7 +526,7 @@ ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag
     }
 
     // find the largest unreduced block at the end of the matrix.
-    while (end>0 && subdiag[end-1]==RealScalar(0))
+    while (end>0 && numext::is_exactly_zero(subdiag[end - 1]))
     {
       end--;
     }
@@ -538,7 +538,7 @@ ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag
     if(iter > maxIterations * n) break;
 
     start = end - 1;
-    while (start>0 && subdiag[start-1]!=0)
+    while (start>0 && !numext::is_exactly_zero(subdiag[start - 1]))
       start--;
 
     internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), subdiag.data(), start, end, computeEigenvectors ? eivec.data() : (Scalar*)0, n);
@@ -843,12 +843,12 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta
   //   RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
   // This explain the following, somewhat more complicated, version:
   RealScalar mu = diag[end];
-  if(td==RealScalar(0)) {
+  if(numext::is_exactly_zero(td)) {
     mu -= numext::abs(e);
-  } else if (e != RealScalar(0)) {
+  } else if (!numext::is_exactly_zero(e)) {
     const RealScalar e2 = numext::abs2(e);
     const RealScalar h = numext::hypot(td,e);
-    if(e2 == RealScalar(0)) {
+    if(numext::is_exactly_zero(e2)) {
       mu -= e / ((td + (td>RealScalar(0) ? h : -h)) / e);
     } else {
       mu -= e2 / (td + (td>RealScalar(0) ? h : -h)); 
@@ -859,7 +859,7 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta
   RealScalar z = subdiag[start];
   // If z ever becomes zero, the Givens rotation will be the identity and
   // z will stay zero for all future iterations.
-  for (Index k = start; k < end && z != RealScalar(0); ++k)
+  for (Index k = start; k < end && !numext::is_exactly_zero(z); ++k)
   {
     JacobiRotation<RealScalar> rot;
     rot.makeGivens(x, z);
diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h
index dc6bf3e14..bf3f45667 100644
--- a/Eigen/src/Householder/Householder.h
+++ b/Eigen/src/Householder/Householder.h
@@ -124,7 +124,7 @@ void MatrixBase<Derived>::applyHouseholderOnTheLeft(
   {
     *this *= Scalar(1)-tau;
   }
-  else if(tau!=Scalar(0))
+  else if(!numext::is_exactly_zero(tau))
   {
     Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace,cols());
     Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows()-1, cols());
@@ -162,7 +162,7 @@ void MatrixBase<Derived>::applyHouseholderOnTheRight(
   {
     *this *= Scalar(1)-tau;
   }
-  else if(tau!=Scalar(0))
+  else if(!numext::is_exactly_zero(tau))
   {
     Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace,rows());
     Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols()-1);
diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h
index 6a533a06b..d515a1f3f 100644
--- a/Eigen/src/Jacobi/Jacobi.h
+++ b/Eigen/src/Jacobi/Jacobi.h
@@ -234,13 +234,13 @@ void JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar
 {
   using std::sqrt;
   using std::abs;
-  if(q==Scalar(0))
+  if(numext::is_exactly_zero(q))
   {
     m_c = p<Scalar(0) ? Scalar(-1) : Scalar(1);
     m_s = Scalar(0);
     if(r) *r = abs(p);
   }
-  else if(p==Scalar(0))
+  else if(numext::is_exactly_zero(p))
   {
     m_c = Scalar(0);
     m_s = q<Scalar(0) ? Scalar(1) : Scalar(-1);
@@ -468,7 +468,7 @@ void /*EIGEN_DONT_INLINE*/ apply_rotation_in_the_plane(DenseBase<VectorX>& xpr_x
 
   OtherScalar c = j.c();
   OtherScalar s = j.s();
-  if (c==OtherScalar(1) && s==OtherScalar(0))
+  if (numext::is_exactly_one(c) && numext::is_exactly_zero(s))
     return;
 
   apply_rotation_in_the_plane_selector<
diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index fce7c34e1..259b549b9 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -519,7 +519,7 @@ void FullPivLU<MatrixType>::computeInPlace()
     row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
     col_of_biggest_in_corner += k; // need to add k to them.
 
-    if(biggest_in_corner==Score(0))
+    if(numext::is_exactly_zero(biggest_in_corner))
     {
       // before exiting, make sure to initialize the still uninitialized transpositions
       // in a sane state without destroying what we already have.
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h
index aba4a67d1..3a32f19e5 100644
--- a/Eigen/src/LU/PartialPivLU.h
+++ b/Eigen/src/LU/PartialPivLU.h
@@ -378,7 +378,7 @@ struct partial_lu_impl
 
       row_transpositions[k] = PivIndex(row_of_biggest_in_col);
 
-      if(biggest_in_corner != Score(0))
+      if(!numext::is_exactly_zero(biggest_in_corner))
       {
         if(k != row_of_biggest_in_col)
         {
@@ -404,7 +404,7 @@ struct partial_lu_impl
     {
       Index k = endk;
       row_transpositions[k] = PivIndex(k);
-      if (Scoring()(lu(k, k)) == Score(0) && first_zero_pivot == -1)
+      if (numext::is_exactly_zero(Scoring()(lu(k, k))) && first_zero_pivot == -1)
         first_zero_pivot = k;
     }
 
diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h
index b9500c865..dd5db97e1 100644
--- a/Eigen/src/QR/ColPivHouseholderQR.h
+++ b/Eigen/src/QR/ColPivHouseholderQR.h
@@ -552,7 +552,7 @@ void ColPivHouseholderQR<MatrixType>::computeInPlace()
       // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf
       // and used in LAPACK routines xGEQPF and xGEQP3.
       // See lines 278-297 in http://www.netlib.org/lapack/explore-html/dc/df4/sgeqpf_8f_source.html
-      if (m_colNormsUpdated.coeffRef(j) != RealScalar(0)) {
+      if (!numext::is_exactly_zero(m_colNormsUpdated.coeffRef(j))) {
         RealScalar temp = abs(m_qr.coeffRef(k, j)) / m_colNormsUpdated.coeffRef(j);
         temp = (RealScalar(1) + temp) * (RealScalar(1) - temp);
         temp = temp <  RealScalar(0) ? RealScalar(0) : temp;
diff --git a/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h b/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h
index f1c29dd4f..9803f1272 100644
--- a/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h
+++ b/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h
@@ -139,7 +139,7 @@ class SPQR : public SparseSolverBase<SPQR<MatrixType_> >
       {
         RealScalar max2Norm = 0.0;
         for (int j = 0; j < mat.cols(); j++) max2Norm = numext::maxi(max2Norm, mat.col(j).norm());
-        if(max2Norm==RealScalar(0))
+        if(numext::is_exactly_zero(max2Norm))
           max2Norm = RealScalar(1);
         pivotThreshold = 20 * (mat.rows() + mat.cols()) * max2Norm * NumTraits<RealScalar>::epsilon();
       }
diff --git a/Eigen/src/SVD/BDCSVD.h b/Eigen/src/SVD/BDCSVD.h
index 0ad453f58..7b2064980 100644
--- a/Eigen/src/SVD/BDCSVD.h
+++ b/Eigen/src/SVD/BDCSVD.h
@@ -282,7 +282,7 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
     return *this;
   }
 
-  if(scale==Literal(0)) scale = Literal(1);
+  if(numext::is_exactly_zero(scale)) scale = Literal(1);
   MatrixX copy;
   if (m_isTranspose) copy = matrix.adjoint()/scale;
   else               copy = matrix/scale;
@@ -621,7 +621,10 @@ void BDCSVD<MatrixType>::computeSVDofM(Eigen::Index firstCol, Eigen::Index n, Ma
   // but others are interleaved and we must ignore them at this stage.
   // To this end, let's compute a permutation skipping them:
   Index actual_n = n;
-  while(actual_n>1 && diag(actual_n-1)==Literal(0)) {--actual_n; eigen_internal_assert(col0(actual_n)==Literal(0)); }
+  while(actual_n>1 && numext::is_exactly_zero(diag(actual_n - 1))) {
+    --actual_n;
+    eigen_internal_assert(numext::is_exactly_zero(col0(actual_n)));
+  }
   Index m = 0; // size of the deflated problem
   for(Index k=0;k<actual_n;++k)
     if(abs(col0(k))>considerZero)
@@ -753,11 +756,11 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
   Index actual_n = n;
   // Note that here actual_n is computed based on col0(i)==0 instead of diag(i)==0 as above
   // because 1) we have diag(i)==0 => col0(i)==0 and 2) if col0(i)==0, then diag(i) is already a singular value.
-  while(actual_n>1 && col0(actual_n-1)==Literal(0)) --actual_n;
+  while(actual_n>1 && numext::is_exactly_zero(col0(actual_n - 1))) --actual_n;
 
   for (Index k = 0; k < n; ++k)
   {
-    if (col0(k) == Literal(0) || actual_n==1)
+    if (numext::is_exactly_zero(col0(k)) || actual_n == 1)
     {
       // if col0(k) == 0, then entry is deflated, so singular value is on diagonal
       // if actual_n==1, then the deflated problem is already diagonalized
@@ -778,7 +781,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
       // recall that at this stage we assume that z[j]!=0 and all entries for which z[j]==0 have been put aside.
       // This should be equivalent to using perm[]
       Index l = k+1;
-      while(col0(l)==Literal(0)) { ++l; eigen_internal_assert(l<actual_n); }
+      while(numext::is_exactly_zero(col0(l))) { ++l; eigen_internal_assert(l < actual_n); }
       right = diag(l);
     }
 
@@ -813,7 +816,8 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
     {
       // check that after the shift, f(mid) is still negative:
       RealScalar midShifted = (right - left) / RealScalar(2);
-      if(shift==right)
+      // we can test exact equality here, because shift comes from `... ? left : right`
+      if(numext::equal_strict(shift, right))
         midShifted = -midShifted;
       RealScalar fMidShifted = secularEq(midShifted, col0, diag, perm, diagShifted, shift);
       if(fMidShifted>0)
@@ -826,7 +830,8 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
 
     // initial guess
     RealScalar muPrev, muCur;
-    if (shift == left)
+    // we can test exact equality here, because shift comes from `... ? left : right`
+    if (numext::equal_strict(shift, left))
     {
       muPrev = (right - left) * RealScalar(0.1);
       if (k == actual_n-1) muCur = right - left;
@@ -849,7 +854,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
     // rational interpolation: fit a function of the form a / mu + b through the two previous
     // iterates and use its zero to compute the next iterate
     bool useBisection = fPrev*fCur>Literal(0);
-    while (fCur!=Literal(0) && abs(muCur - muPrev) > Literal(8) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(abs(muCur), abs(muPrev)) && abs(fCur - fPrev)>NumTraits<RealScalar>::epsilon() && !useBisection)
+    while (!numext::is_exactly_zero(fCur) && abs(muCur - muPrev) > Literal(8) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(abs(muCur), abs(muPrev)) && abs(fCur - fPrev) > NumTraits<RealScalar>::epsilon() && !useBisection)
     {
       ++m_numIters;
 
@@ -869,8 +874,9 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
       muCur = muZero;
       fCur = fZero;
 
-      if (shift == left  && (muCur < Literal(0) || muCur > right - left)) useBisection = true;
-      if (shift == right && (muCur < -(right - left) || muCur > Literal(0))) useBisection = true;
+      // we can test exact equality here, because shift comes from `... ? left : right`
+      if (numext::equal_strict(shift, left)  && (muCur < Literal(0) || muCur > right - left)) useBisection = true;
+      if (numext::equal_strict(shift, right) && (muCur < -(right - left) || muCur > Literal(0))) useBisection = true;
       if (abs(fCur)>abs(fPrev)) useBisection = true;
     }
 
@@ -881,7 +887,8 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
       std::cout << "useBisection for k = " << k << ", actual_n = " << actual_n << "\n";
 #endif
       RealScalar leftShifted, rightShifted;
-      if (shift == left)
+      // we can test exact equality here, because shift comes from `... ? left : right`
+      if (numext::equal_strict(shift, left))
       {
         // to avoid overflow, we must have mu > max(real_min, |z(k)|/sqrt(real_max)),
         // the factor 2 is to be more conservative
@@ -959,7 +966,8 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
         // Instead fo abbording or entering an infinite loop,
         // let's just use the middle as the estimated zero-crossing:
         muCur = (right - left) * RealScalar(0.5);
-        if(shift == right)
+        // we can test exact equality here, because shift comes from `... ? left : right`
+        if(numext::equal_strict(shift, right))
           muCur = -muCur;
       }
     }
@@ -1004,7 +1012,7 @@ void BDCSVD<MatrixType>::perturbCol0
   // The offset permits to skip deflated entries while computing zhat
   for (Index k = 0; k < n; ++k)
   {
-    if (col0(k) == Literal(0)) // deflated
+    if (numext::is_exactly_zero(col0(k))) // deflated
       zhat(k) = Literal(0);
     else
     {
@@ -1077,7 +1085,7 @@ void BDCSVD<MatrixType>::computeSingVecs
 
   for (Index k = 0; k < n; ++k)
   {
-    if (zhat(k) == Literal(0))
+    if (numext::is_exactly_zero(zhat(k)))
     {
       U.col(k) = VectorType::Unit(n+1, k);
       if (m_compV) V.col(k) = VectorType::Unit(n, k);
@@ -1123,7 +1131,7 @@ void BDCSVD<MatrixType>::deflation43(Eigen::Index firstCol, Eigen::Index shift,
   RealScalar c = m_computed(start, start);
   RealScalar s = m_computed(start+i, start);
   RealScalar r = numext::hypot(c,s);
-  if (r == Literal(0))
+  if (numext::is_exactly_zero(r))
   {
     m_computed(start+i, start+i) = Literal(0);
     return;
@@ -1163,7 +1171,7 @@ void BDCSVD<MatrixType>::deflation44(Eigen::Index firstColu , Eigen::Index first
     << m_computed(firstColm + i+1, firstColm+i+1) << " "
     << m_computed(firstColm + i+2, firstColm+i+2) << "\n";
 #endif
-  if (r==Literal(0))
+  if (numext::is_exactly_zero(r))
   {
     m_computed(firstColm + i, firstColm + i) = m_computed(firstColm + j, firstColm + j);
     return;
diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h
index e69d13aa6..c57b20199 100644
--- a/Eigen/src/SVD/JacobiSVD.h
+++ b/Eigen/src/SVD/JacobiSVD.h
@@ -377,7 +377,7 @@ struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
     const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
     const RealScalar precision = NumTraits<Scalar>::epsilon();
 
-    if(n==0)
+    if(numext::is_exactly_zero(n))
     {
       // make sure first column is zero
       work_matrix.coeffRef(p,p) = work_matrix.coeffRef(q,p) = Scalar(0);
@@ -684,7 +684,7 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
     m_info = InvalidInput;
     return *this;
   }
-  if(scale==RealScalar(0)) scale = RealScalar(1);
+  if(numext::is_exactly_zero(scale)) scale = RealScalar(1);
   
   /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
 
@@ -777,7 +777,7 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
   {
     Index pos;
     RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
-    if(maxRemainingSingularValue == RealScalar(0))
+    if(numext::is_exactly_zero(maxRemainingSingularValue))
     {
       m_nonzeroSingularValues = i;
       break;
diff --git a/Eigen/src/SparseCore/TriangularSolver.h b/Eigen/src/SparseCore/TriangularSolver.h
index 76c32f207..73e2a7987 100644
--- a/Eigen/src/SparseCore/TriangularSolver.h
+++ b/Eigen/src/SparseCore/TriangularSolver.h
@@ -116,7 +116,7 @@ struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,ColMajor>
       for(Index i=0; i<lhs.cols(); ++i)
       {
         Scalar& tmp = other.coeffRef(i,col);
-        if (tmp!=Scalar(0)) // optimization when other is actually sparse
+        if (!numext::is_exactly_zero(tmp)) // optimization when other is actually sparse
         {
           LhsIterator it(lhsEval, i);
           while(it && it.index()<i)
@@ -151,7 +151,7 @@ struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,ColMajor>
       for(Index i=lhs.cols()-1; i>=0; --i)
       {
         Scalar& tmp = other.coeffRef(i,col);
-        if (tmp!=Scalar(0)) // optimization when other is actually sparse
+        if (!numext::is_exactly_zero(tmp)) // optimization when other is actually sparse
         {
           if(!(Mode & UnitDiag))
           {
@@ -241,7 +241,7 @@ struct sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
       {
         tempVector.restart();
         Scalar& ci = tempVector.coeffRef(i);
-        if (ci!=Scalar(0))
+        if (!numext::is_exactly_zero(ci))
         {
           // find
           typename Lhs::InnerIterator it(lhs, i);
diff --git a/test/AnnoyingScalar.h b/test/AnnoyingScalar.h
index d4cca7963..4362de263 100644
--- a/test/AnnoyingScalar.h
+++ b/test/AnnoyingScalar.h
@@ -32,12 +32,12 @@ class AnnoyingScalar
 {
   public:
     AnnoyingScalar()                { init(); *v = 0;  }
-    AnnoyingScalar(long double _v)  { init(); *v = _v; }
-    AnnoyingScalar(double _v)       { init(); *v = _v; }
+    AnnoyingScalar(long double _v)  { init(); *v = static_cast<float>(_v); }
+    AnnoyingScalar(double _v)       { init(); *v = static_cast<float>(_v); }
     AnnoyingScalar(float _v)        { init(); *v = _v; }
-    AnnoyingScalar(int _v)          { init(); *v = _v; }
-    AnnoyingScalar(long _v)         { init(); *v = _v; }
-    AnnoyingScalar(long long _v)    { init(); *v = _v; }
+    AnnoyingScalar(int _v)          { init(); *v = static_cast<float>(_v); }
+    AnnoyingScalar(long _v)         { init(); *v = static_cast<float>(_v); }
+    AnnoyingScalar(long long _v)    { init(); *v = static_cast<float>(_v); }
     AnnoyingScalar(const AnnoyingScalar& other) { init(); *v = *(other.v); }
     ~AnnoyingScalar() {
       if(v!=&data)
@@ -81,8 +81,8 @@ class AnnoyingScalar
     AnnoyingScalar& operator/=(const AnnoyingScalar& other) { *v /= *other.v; return *this; }
     AnnoyingScalar& operator= (const AnnoyingScalar& other) { *v  = *other.v; return *this; }
 
-    bool operator==(const AnnoyingScalar& other) const { return *v == *other.v; }
-    bool operator!=(const AnnoyingScalar& other) const { return *v != *other.v; }
+    bool operator==(const AnnoyingScalar& other) const { return numext::equal_strict(*v, *other.v); }
+    bool operator!=(const AnnoyingScalar& other) const { return numext::not_equal_strict(*v, *other.v); }
     bool operator<=(const AnnoyingScalar& other) const { return *v <= *other.v; }
     bool operator< (const AnnoyingScalar& other) const { return *v <  *other.v; }
     bool operator>=(const AnnoyingScalar& other) const { return *v >= *other.v; }
diff --git a/test/adjoint.cpp b/test/adjoint.cpp
index 84f430cc5..da8f9589c 100644
--- a/test/adjoint.cpp
+++ b/test/adjoint.cpp
@@ -45,7 +45,7 @@ template<> struct adjoint_specific<false> {
     VERIFY_IS_APPROX((v1*0).normalized(), (v1*0));
 #if (!EIGEN_ARCH_i386) || defined(EIGEN_VECTORIZE)
     RealScalar very_small = (std::numeric_limits<RealScalar>::min)();
-    VERIFY( (v1*very_small).norm() == 0 );
+    VERIFY( numext::is_exactly_zero((v1*very_small).norm()) );
     VERIFY_IS_APPROX((v1*very_small).normalized(), (v1*very_small));
     v3 = v1*very_small;
     v3.normalize();
diff --git a/test/block.cpp b/test/block.cpp
index 43849adb6..a2396d1a9 100644
--- a/test/block.cpp
+++ b/test/block.cpp
@@ -149,11 +149,11 @@ template<typename MatrixType> void block(const MatrixType& m)
   }
 
   // stress some basic stuffs with block matrices
-  VERIFY(numext::real(ones.col(c1).sum()) == RealScalar(rows));
-  VERIFY(numext::real(ones.row(r1).sum()) == RealScalar(cols));
+  VERIFY_IS_EQUAL(numext::real(ones.col(c1).sum()), RealScalar(rows));
+  VERIFY_IS_EQUAL(numext::real(ones.row(r1).sum()), RealScalar(cols));
 
-  VERIFY(numext::real(ones.col(c1).dot(ones.col(c2))) == RealScalar(rows));
-  VERIFY(numext::real(ones.row(r1).dot(ones.row(r2))) == RealScalar(cols));
+  VERIFY_IS_EQUAL(numext::real(ones.col(c1).dot(ones.col(c2))), RealScalar(rows));
+  VERIFY_IS_EQUAL(numext::real(ones.row(r1).dot(ones.row(r2))), RealScalar(cols));
   
   // check that linear acccessors works on blocks
   m1 = m1_copy;
diff --git a/test/geo_eulerangles.cpp b/test/geo_eulerangles.cpp
index 693c627a9..bea241965 100644
--- a/test/geo_eulerangles.cpp
+++ b/test/geo_eulerangles.cpp
@@ -26,7 +26,7 @@ void verify_euler(const Matrix<Scalar,3,1>& ea, int i, int j, int k)
   VERIFY_IS_APPROX(m,  mbis); 
   /* If I==K, and ea[1]==0, then there no unique solution. */ 
   /* The remark apply in the case where I!=K, and |ea[1]| is close to pi/2. */ 
-  if( (i!=k || ea[1]!=0) && (i==k || !internal::isApprox(abs(ea[1]),Scalar(EIGEN_PI/2),test_precision<Scalar>())) ) 
+  if((i!=k || !numext::is_exactly_zero(ea[1])) && (i == k || !internal::isApprox(abs(ea[1]), Scalar(EIGEN_PI / 2), test_precision<Scalar>())) )
     VERIFY((ea-eabis).norm() <= test_precision<Scalar>());
   
   // approx_or_less_than does not work for 0
diff --git a/test/mapstride.cpp b/test/mapstride.cpp
index fde73f2ec..93e880acc 100644
--- a/test/mapstride.cpp
+++ b/test/mapstride.cpp
@@ -29,8 +29,8 @@ template<int Alignment,typename VectorType> void map_class_vector(const VectorTy
     map = v;
     for(int i = 0; i < size; ++i)
     {
-      VERIFY(array[3*i] == v[i]);
-      VERIFY(map[i] == v[i]);
+      VERIFY_IS_EQUAL(array[3*i], v[i]);
+      VERIFY_IS_EQUAL(map[i], v[i]);
     }
   }
 
@@ -39,8 +39,8 @@ template<int Alignment,typename VectorType> void map_class_vector(const VectorTy
     map = v;
     for(int i = 0; i < size; ++i)
     {
-      VERIFY(array[2*i] == v[i]);
-      VERIFY(map[i] == v[i]);
+      VERIFY_IS_EQUAL(array[2*i], v[i]);
+      VERIFY_IS_EQUAL(map[i], v[i]);
     }
   }
 
@@ -84,8 +84,8 @@ template<int Alignment,typename MatrixType> void map_class_matrix(const MatrixTy
     for(int i = 0; i < m.outerSize(); ++i)
       for(int j = 0; j < m.innerSize(); ++j)
       {
-        VERIFY(array[map.outerStride()*i+j] == m.coeffByOuterInner(i,j));
-        VERIFY(map.coeffByOuterInner(i,j) == m.coeffByOuterInner(i,j));
+        VERIFY_IS_EQUAL(array[map.outerStride()*i+j], m.coeffByOuterInner(i,j));
+        VERIFY_IS_EQUAL(map.coeffByOuterInner(i,j), m.coeffByOuterInner(i,j));
       }
     VERIFY_IS_APPROX(s1*map,s1*m);
     map *= s1;
@@ -111,8 +111,8 @@ template<int Alignment,typename MatrixType> void map_class_matrix(const MatrixTy
     for(int i = 0; i < m.outerSize(); ++i)
       for(int j = 0; j < m.innerSize(); ++j)
       {
-        VERIFY(array[map.outerStride()*i+j] == m.coeffByOuterInner(i,j));
-        VERIFY(map.coeffByOuterInner(i,j) == m.coeffByOuterInner(i,j));
+        VERIFY_IS_EQUAL(array[map.outerStride()*i+j], m.coeffByOuterInner(i,j));
+        VERIFY_IS_EQUAL(map.coeffByOuterInner(i,j), m.coeffByOuterInner(i,j));
       }
     VERIFY_IS_APPROX(s1*map,s1*m);
     map *= s1;
@@ -133,8 +133,8 @@ template<int Alignment,typename MatrixType> void map_class_matrix(const MatrixTy
     for(int i = 0; i < m.outerSize(); ++i)
       for(int j = 0; j < m.innerSize(); ++j)
       {
-        VERIFY(array[map.outerStride()*i+map.innerStride()*j] == m.coeffByOuterInner(i,j));
-        VERIFY(map.coeffByOuterInner(i,j) == m.coeffByOuterInner(i,j));
+        VERIFY_IS_EQUAL(array[map.outerStride()*i+map.innerStride()*j], m.coeffByOuterInner(i,j));
+        VERIFY_IS_EQUAL(map.coeffByOuterInner(i,j), m.coeffByOuterInner(i,j));
       }
     VERIFY_IS_APPROX(s1*map,s1*m);
     map *= s1;
@@ -154,8 +154,8 @@ template<int Alignment,typename MatrixType> void map_class_matrix(const MatrixTy
     for(int i = 0; i < m.outerSize(); ++i)
       for(int j = 0; j < m.innerSize(); ++j)
       {
-        VERIFY(array[map.innerSize()*i*2+j*2] == m.coeffByOuterInner(i,j));
-        VERIFY(map.coeffByOuterInner(i,j) == m.coeffByOuterInner(i,j));
+        VERIFY_IS_EQUAL(array[map.innerSize()*i*2+j*2], m.coeffByOuterInner(i,j));
+        VERIFY_IS_EQUAL(map.coeffByOuterInner(i,j), m.coeffByOuterInner(i,j));
       }
     VERIFY_IS_APPROX(s1*map,s1*m);
     map *= s1;
diff --git a/test/nullary.cpp b/test/nullary.cpp
index 9b25ea4f3..2c4d93806 100644
--- a/test/nullary.cpp
+++ b/test/nullary.cpp
@@ -13,24 +13,20 @@
 template<typename MatrixType>
 bool equalsIdentity(const MatrixType& A)
 {
-  typedef typename MatrixType::Scalar Scalar;
-  Scalar zero = static_cast<Scalar>(0);
-
   bool offDiagOK = true;
   for (Index i = 0; i < A.rows(); ++i) {
     for (Index j = i+1; j < A.cols(); ++j) {
-      offDiagOK = offDiagOK && (A(i,j) == zero);
+      offDiagOK = offDiagOK && numext::is_exactly_zero(A(i, j));
     }
   }
   for (Index i = 0; i < A.rows(); ++i) {
     for (Index j = 0; j < (std::min)(i, A.cols()); ++j) {
-      offDiagOK = offDiagOK && (A(i,j) == zero);
+      offDiagOK = offDiagOK && numext::is_exactly_zero(A(i, j));
     }
   }
 
   bool diagOK = (A.diagonal().array() == 1).all();
   return offDiagOK && diagOK;
-
 }
 
 template<typename VectorType>
diff --git a/test/numext.cpp b/test/numext.cpp
index 8a2fde501..ee879c9ac 100644
--- a/test/numext.cpp
+++ b/test/numext.cpp
@@ -11,7 +11,7 @@
 
 template<typename T, typename U>
 bool check_if_equal_or_nans(const T& actual, const U& expected) {
-  return ((actual == expected) || ((numext::isnan)(actual) && (numext::isnan)(expected)));
+  return (numext::equal_strict(actual, expected) || ((numext::isnan)(actual) && (numext::isnan)(expected)));
 }
 
 template<typename T, typename U>
diff --git a/test/packetmath_test_shared.h b/test/packetmath_test_shared.h
index 8624fe2fe..709c4ae43 100644
--- a/test/packetmath_test_shared.h
+++ b/test/packetmath_test_shared.h
@@ -100,7 +100,7 @@ template<typename Scalar> bool areApprox(const Scalar* a, const Scalar* b, int s
 {
   for (int i=0; i<size; ++i)
   {
-    if ( a[i]!=b[i] && !internal::isApprox(a[i],b[i]) 
+    if ( numext::not_equal_strict(a[i], b[i]) && !internal::isApprox(a[i],b[i])
          && !((numext::isnan)(a[i]) && (numext::isnan)(b[i])) )
     {
       print_mismatch(a, b, size);
@@ -114,7 +114,7 @@ template<typename Scalar> bool areEqual(const Scalar* a, const Scalar* b, int si
 {
   for (int i=0; i<size; ++i)
   {
-    if ( (a[i] != b[i]) && !((numext::isnan)(a[i]) && (numext::isnan)(b[i])) )
+    if ( numext::not_equal_strict(a[i], b[i]) && !((numext::isnan)(a[i]) && (numext::isnan)(b[i])) )
     {
       print_mismatch(a, b, size);
       return false;
diff --git a/test/real_qz.cpp b/test/real_qz.cpp
index 1cf7aba2d..ea4a2701a 100644
--- a/test/real_qz.cpp
+++ b/test/real_qz.cpp
@@ -18,7 +18,6 @@ template<typename MatrixType> void real_qz(const MatrixType& m)
      RealQZ.h
   */
   using std::abs;
-  typedef typename MatrixType::Scalar Scalar;
   
   Index dim = m.cols();
   
@@ -52,17 +51,18 @@ template<typename MatrixType> void real_qz(const MatrixType& m)
   bool all_zeros = true;
   for (Index i=0; i<A.cols(); i++)
     for (Index j=0; j<i; j++) {
-      if (abs(qz.matrixT()(i,j))!=Scalar(0.0))
+      if (!numext::is_exactly_zero(abs(qz.matrixT()(i, j))))
       {
         std::cerr << "Error: T(" << i << "," << j << ") = " << qz.matrixT()(i,j) << std::endl;
         all_zeros = false;
       }
-      if (j<i-1 && abs(qz.matrixS()(i,j))!=Scalar(0.0))
+      if (j<i-1 && !numext::is_exactly_zero(abs(qz.matrixS()(i, j))))
       {
         std::cerr << "Error: S(" << i << "," << j << ") = " << qz.matrixS()(i,j) << std::endl;
         all_zeros = false;
       }
-      if (j==i-1 && j>0 && abs(qz.matrixS()(i,j))!=Scalar(0.0) && abs(qz.matrixS()(i-1,j-1))!=Scalar(0.0))
+      if (j==i-1 && j>0 && !numext::is_exactly_zero(abs(qz.matrixS()(i, j))) &&
+              !numext::is_exactly_zero(abs(qz.matrixS()(i - 1, j - 1))))
       {
         std::cerr << "Error: S(" << i << "," << j << ") = " << qz.matrixS()(i,j)  << " && S(" << i-1 << "," << j-1 << ") = " << qz.matrixS()(i-1,j-1) << std::endl;
         all_zeros = false;
diff --git a/test/schur_real.cpp b/test/schur_real.cpp
index 945461027..8b40ddd7e 100644
--- a/test/schur_real.cpp
+++ b/test/schur_real.cpp
@@ -19,15 +19,15 @@ template<typename MatrixType> void verifyIsQuasiTriangular(const MatrixType& T)
   // Check T is lower Hessenberg
   for(int row = 2; row < size; ++row) {
     for(int col = 0; col < row - 1; ++col) {
-      VERIFY(T(row,col) == Scalar(0));
+      VERIFY_IS_EQUAL(T(row,col), Scalar(0));
     }
   }
 
   // Check that any non-zero on the subdiagonal is followed by a zero and is
   // part of a 2x2 diagonal block with imaginary eigenvalues.
   for(int row = 1; row < size; ++row) {
-    if (T(row,row-1) != Scalar(0)) {
-      VERIFY(row == size-1 || T(row+1,row) == 0);
+    if (!numext::is_exactly_zero(T(row, row - 1))) {
+      VERIFY(row == size-1 || numext::is_exactly_zero(T(row + 1, row)));
       Scalar tr = T(row-1,row-1) + T(row,row);
       Scalar det = T(row-1,row-1) * T(row,row) - T(row-1,row) * T(row,row-1);
       VERIFY(4 * det > tr * tr);
diff --git a/test/sparse.h b/test/sparse.h
index 9a63e0deb..f3e697d5e 100644
--- a/test/sparse.h
+++ b/test/sparse.h
@@ -54,7 +54,8 @@ initSparse(double density,
   enum { IsRowMajor = SparseMatrix<Scalar,Opt2,StorageIndex>::IsRowMajor };
   sparseMat.setZero();
   //sparseMat.reserve(int(refMat.rows()*refMat.cols()*density));
-  sparseMat.reserve(VectorXi::Constant(IsRowMajor ? refMat.rows() : refMat.cols(), int((1.5*density)*(IsRowMajor?refMat.cols():refMat.rows()))));
+  int nnz = static_cast<int>((1.5 * density) * static_cast<double>(IsRowMajor ? refMat.cols() : refMat.rows()));
+  sparseMat.reserve(VectorXi::Constant(IsRowMajor ? refMat.rows() : refMat.cols(), nnz));
 
   Index insert_count = 0;
   for(Index j=0; j<sparseMat.outerSize(); j++)
@@ -82,7 +83,7 @@ initSparse(double density,
       if ((flags&ForceRealDiag) && (i==j))
         v = numext::real(v);
 
-      if (v!=Scalar(0))
+      if (!numext::is_exactly_zero(v))
       {
         //sparseMat.insertBackByOuterInner(j,i) = v;
         sparseMat.insertByOuterInner(j,i) = v;
@@ -115,7 +116,7 @@ initSparse(double density,
   for(int i=0; i<refVec.size(); i++)
   {
     Scalar v = (internal::random<double>(0,1) < density) ? internal::random<Scalar>() : Scalar(0);
-    if (v!=Scalar(0))
+    if (!numext::is_exactly_zero(v))
     {
       sparseVec.insertBack(i) = v;
       if (nonzeroCoords)
diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp
index 85a607700..86944900b 100644
--- a/test/sparse_basic.cpp
+++ b/test/sparse_basic.cpp
@@ -679,7 +679,7 @@ void big_sparse_triplet(Index rows, Index cols, double density) {
   typedef typename SparseMatrixType::Scalar Scalar;
   typedef Triplet<Scalar,Index> TripletType;
   std::vector<TripletType> triplets;
-  double nelements = density * rows*cols;
+  double nelements = density * static_cast<double>(rows*cols);
   VERIFY(nelements>=0 && nelements < static_cast<double>(NumTraits<StorageIndex>::highest()));
   Index ntriplets = Index(nelements);
   triplets.reserve(ntriplets);
diff --git a/test/sparse_block.cpp b/test/sparse_block.cpp
index b4905b053..452926f91 100644
--- a/test/sparse_block.cpp
+++ b/test/sparse_block.cpp
@@ -90,11 +90,11 @@ template<typename SparseMatrixType> void sparse_block(const SparseMatrixType& re
           
           VERIFY_IS_APPROX(m.middleCols(j,w).coeff(r,c), refMat.middleCols(j,w).coeff(r,c));
           VERIFY_IS_APPROX(m.middleRows(i,h).coeff(r,c), refMat.middleRows(i,h).coeff(r,c));
-          if(m.middleCols(j,w).coeff(r,c) != Scalar(0))
+          if(!numext::is_exactly_zero(m.middleCols(j, w).coeff(r, c)))
           {
             VERIFY_IS_APPROX(m.middleCols(j,w).coeffRef(r,c), refMat.middleCols(j,w).coeff(r,c));
           }
-          if(m.middleRows(i,h).coeff(r,c) != Scalar(0))
+          if(!numext::is_exactly_zero(m.middleRows(i, h).coeff(r, c)))
           {
             VERIFY_IS_APPROX(m.middleRows(i,h).coeff(r,c), refMat.middleRows(i,h).coeff(r,c));
           }
@@ -166,14 +166,14 @@ template<typename SparseMatrixType> void sparse_block(const SparseMatrixType& re
     {
       VERIFY(j==numext::real(m3.innerVector(j).nonZeros()));
       if(j>0)
-        VERIFY(RealScalar(j)==numext::real(m3.innerVector(j).lastCoeff()));
+        VERIFY_IS_EQUAL(RealScalar(j), numext::real(m3.innerVector(j).lastCoeff()));
     }
     m3.makeCompressed();
     for(Index j=0; j<(std::min)(outer, inner); ++j)
     {
       VERIFY(j==numext::real(m3.innerVector(j).nonZeros()));
       if(j>0)
-        VERIFY(RealScalar(j)==numext::real(m3.innerVector(j).lastCoeff()));
+        VERIFY_IS_EQUAL(RealScalar(j), numext::real(m3.innerVector(j).lastCoeff()));
     }
 
     VERIFY(m3.innerVector(j0).nonZeros() == m3.transpose().innerVector(j0).nonZeros());
diff --git a/test/sparse_permutations.cpp b/test/sparse_permutations.cpp
index e93493c39..5974c74f8 100644
--- a/test/sparse_permutations.cpp
+++ b/test/sparse_permutations.cpp
@@ -53,7 +53,7 @@ template<int OtherStorage, typename SparseMatrixType> void sparse_permutations(c
 //   bool IsRowMajor1 = SparseMatrixType::IsRowMajor;
 //   bool IsRowMajor2 = OtherSparseMatrixType::IsRowMajor;
   
-  double density = (std::max)(8./(rows*cols), 0.01);
+  double density = (std::max)(8./static_cast<double>(rows*cols), 0.01);
   
   SparseMatrixType mat(rows, cols), up(rows,cols), lo(rows,cols);
   OtherSparseMatrixType res;
diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp
index 488a3920c..dbf549a8c 100644
--- a/test/sparse_product.cpp
+++ b/test/sparse_product.cpp
@@ -390,7 +390,7 @@ void test_mixing_types()
   typedef Matrix<Cplx,Dynamic,Dynamic> DenseMatCplx;
 
   Index n = internal::random<Index>(1,100);
-  double density = (std::max)(8./(n*n), 0.2);
+  double density = (std::max)(8./static_cast<double>(n*n), 0.2);
 
   SpMatReal sR1(n,n);
   SpMatCplx sC1(n,n), sC2(n,n), sC3(n,n);
diff --git a/test/sparse_solver.h b/test/sparse_solver.h
index 2b3b403ce..a9b18b803 100644
--- a/test/sparse_solver.h
+++ b/test/sparse_solver.h
@@ -350,7 +350,7 @@ int generate_sparse_spd_problem(Solver& , typename Solver::MatrixType& A, typena
   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
 
   int size = internal::random<int>(1,maxSize);
-  double density = (std::max)(8./(size*size), 0.01);
+  double density = (std::max)(8./static_cast<double>(size*size), 0.01);
 
   Mat M(size, size);
   DenseMatrix dM(size, size);
@@ -419,7 +419,7 @@ template<typename Solver> void check_sparse_spd_solving(Solver& solver, int maxS
 
     // generate the right hand sides
     int rhsCols = internal::random<int>(1,16);
-    double density = (std::max)(8./(size*rhsCols), 0.1);
+    double density = (std::max)(8./static_cast<double>(size*rhsCols), 0.1);
     SpMat B(size,rhsCols);
     DenseVector b = DenseVector::Random(size);
     DenseMatrix dB(size,rhsCols);
@@ -510,7 +510,7 @@ Index generate_sparse_square_problem(Solver&, typename Solver::MatrixType& A, De
   typedef typename Mat::Scalar Scalar;
 
   Index size = internal::random<int>(1,maxSize);
-  double density = (std::max)(8./(size*size), 0.01);
+  double density = (std::max)(8./static_cast<double>(size*size), 0.01);
   
   A.resize(size,size);
   dA.resize(size,size);
@@ -551,7 +551,7 @@ template<typename Solver> void check_sparse_square_solving(Solver& solver, int m
     DenseVector b = DenseVector::Random(size);
     DenseMatrix dB(size,rhsCols);
     SpMat B(size,rhsCols);
-    double density = (std::max)(8./(size*rhsCols), 0.1);
+    double density = (std::max)(8./double(size*rhsCols), 0.1);
     initSparse<Scalar>(density, dB, B, ForceNonZeroDiag);
     B.makeCompressed();
     SpVec c = B.col(0);
diff --git a/test/sparse_vector.cpp b/test/sparse_vector.cpp
index 7bd57cdbe..e1e74f340 100644
--- a/test/sparse_vector.cpp
+++ b/test/sparse_vector.cpp
@@ -47,8 +47,8 @@ template<typename Scalar,typename StorageIndex> void sparse_vector(int rows, int
     for (typename SparseVectorType::InnerIterator it(v1); it; ++it,++j)
     {
       VERIFY(nonzerocoords[j]==it.index());
-      VERIFY(it.value()==v1.coeff(it.index()));
-      VERIFY(it.value()==refV1.coeff(it.index()));
+      VERIFY_IS_EQUAL(it.value(), v1.coeff(it.index()));
+      VERIFY_IS_EQUAL(it.value(), refV1.coeff(it.index()));
     }
   }
   VERIFY_IS_APPROX(v1, refV1);
diff --git a/test/stl_iterators.cpp b/test/stl_iterators.cpp
index aab3be9e8..121eb865e 100644
--- a/test/stl_iterators.cpp
+++ b/test/stl_iterators.cpp
@@ -438,14 +438,14 @@ void test_stl_iterators(int rows=Rows, int cols=Cols)
     i = internal::random<Index>(0,A.rows()-1);
     A.setRandom();
     A.row(i).setZero();
-    VERIFY_IS_EQUAL( std::find_if(A.rowwise().begin(),  A.rowwise().end(),  [](typename ColMatrixType::RowXpr x) { return x.squaredNorm() == Scalar(0); })-A.rowwise().begin(),  i );
-    VERIFY_IS_EQUAL( std::find_if(A.rowwise().rbegin(), A.rowwise().rend(), [](typename ColMatrixType::RowXpr x) { return x.squaredNorm() == Scalar(0); })-A.rowwise().rbegin(), (A.rows()-1) - i );
+    VERIFY_IS_EQUAL(std::find_if(A.rowwise().begin(),  A.rowwise().end(),  [](typename ColMatrixType::RowXpr x) { return numext::is_exactly_zero(x.squaredNorm()); }) - A.rowwise().begin(), i );
+    VERIFY_IS_EQUAL(std::find_if(A.rowwise().rbegin(), A.rowwise().rend(), [](typename ColMatrixType::RowXpr x) { return numext::is_exactly_zero(x.squaredNorm()); }) - A.rowwise().rbegin(), (A.rows() - 1) - i );
 
     j = internal::random<Index>(0,A.cols()-1);
     A.setRandom();
     A.col(j).setZero();
-    VERIFY_IS_EQUAL( std::find_if(A.colwise().begin(),  A.colwise().end(),  [](typename ColMatrixType::ColXpr x) { return x.squaredNorm() == Scalar(0); })-A.colwise().begin(),  j );
-    VERIFY_IS_EQUAL( std::find_if(A.colwise().rbegin(), A.colwise().rend(), [](typename ColMatrixType::ColXpr x) { return x.squaredNorm() == Scalar(0); })-A.colwise().rbegin(), (A.cols()-1) - j );
+    VERIFY_IS_EQUAL(std::find_if(A.colwise().begin(),  A.colwise().end(),  [](typename ColMatrixType::ColXpr x) { return numext::is_exactly_zero(x.squaredNorm()); }) - A.colwise().begin(), j );
+    VERIFY_IS_EQUAL(std::find_if(A.colwise().rbegin(), A.colwise().rend(), [](typename ColMatrixType::ColXpr x) { return numext::is_exactly_zero(x.squaredNorm()); }) - A.colwise().rbegin(), (A.cols() - 1) - j );
   }
 
   {
diff --git a/test/visitor.cpp b/test/visitor.cpp
index 20fb2c3ed..05c2a4838 100644
--- a/test/visitor.cpp
+++ b/test/visitor.cpp
@@ -21,7 +21,7 @@ template<typename MatrixType> void matrixVisitor(const MatrixType& p)
   m = MatrixType::Random(rows, cols);
   for(Index i = 0; i < m.size(); i++)
     for(Index i2 = 0; i2 < i; i2++)
-      while(m(i) == m(i2)) // yes, ==
+      while(numext::equal_strict(m(i), m(i2))) // yes, strict equality
         m(i) = internal::random<Scalar>();
   
   Scalar minc = Scalar(1000), maxc = Scalar(-1000);