Much more convenient, less over-engineered NumTraits. Done during this KDE-Edu weekend.

src/Core/Conjugate.h

@@ -52,7 +52,7 @@ template<typename MatrixType> class Conjugate
     
     Scalar _read(int row, int col) const
     {
-      return NumTraits<Scalar>::conj(m_matrix.read(row, col));
+      return conj(m_matrix.read(row, col));
     }
     
   protected:

src/Core/Dot.h

@@ -32,7 +32,7 @@ struct DotUnroller
   static void run(const Derived1 &v1, const Derived2& v2, typename Derived1::Scalar &dot)
   {
     DotUnroller<Index-1, Size, Derived1, Derived2>::run(v1, v2, dot);
-    dot += v1[Index] * NumTraits<typename Derived1::Scalar>::conj(v2[Index]);
+    dot += v1[Index] * conj(v2[Index]);
   }
 };
 
@@ -41,7 +41,7 @@ struct DotUnroller<0, Size, Derived1, Derived2>
 {
   static void run(const Derived1 &v1, const Derived2& v2, typename Derived1::Scalar &dot)
   {
-    dot = v1[0] * NumTraits<typename Derived1::Scalar>::conj(v2[0]);
+    dot = v1[0] * conj(v2[0]);
   }
 };
 
@@ -67,9 +67,9 @@ Scalar MatrixBase<Scalar, Derived>::dot(const OtherDerived& other) const
       ::run(*static_cast<const Derived*>(this), other, res);
   else
   {
-    res = (*this)[0] * NumTraits<Scalar>::conj(other[0]);
+    res = (*this)[0] * conj(other[0]);
     for(int i = 1; i < size(); i++)
-      res += (*this)[i]* NumTraits<Scalar>::conj(other[i]);
+      res += (*this)[i]* conj(other[i]);
   }
   return res;
 }
@@ -77,13 +77,13 @@ Scalar MatrixBase<Scalar, Derived>::dot(const OtherDerived& other) const
 template<typename Scalar, typename Derived>
 typename NumTraits<Scalar>::Real MatrixBase<Scalar, Derived>::norm2() const
 {
-  return NumTraits<Scalar>::real(dot(*this));
+  return real(dot(*this));
 }
 
 template<typename Scalar, typename Derived>
 typename NumTraits<Scalar>::Real MatrixBase<Scalar, Derived>::norm() const
 {
-  return std::sqrt(norm2());
+  return sqrt(norm2());
 }
 
 template<typename Scalar, typename Derived>

src/Core/Fuzzy.h

@@ -56,12 +56,12 @@ bool MatrixBase<Scalar, Derived>::isMuchSmallerThan(
 {
   if(IsVector)
   {
-    return(norm2() <= NumTraits<Scalar>::abs2(other) * prec * prec);
+    return(norm2() <= abs2(other) * prec * prec);
   }
   else
   {
     for(int i = 0; i < cols(); i++)
-      if(col(i).norm2() > NumTraits<Scalar>::abs2(other) * prec * prec)
+      if(col(i).norm2() > abs2(other) * prec * prec)
         return false;
     return true;
   }

src/Core/MatrixBase.h

@@ -37,16 +37,16 @@ template<typename Scalar, typename Derived> class MatrixBase
     template<typename OtherDerived>
     bool _isApprox_helper(
       const OtherDerived& other,
-      const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::precision()
+      const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()
     ) const;
     bool _isMuchSmallerThan_helper(
       const Scalar& other,
-      const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::precision()
+      const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()
     ) const;
     template<typename OtherDerived>
     bool _isMuchSmallerThan_helper(
       const MatrixBase<Scalar, OtherDerived>& other,
-      const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::precision()
+      const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()
     ) const;
     
   public:
@@ -121,16 +121,16 @@ template<typename Scalar, typename Derived> class MatrixBase
     template<typename OtherDerived>
     bool isApprox(
       const OtherDerived& other,
-      const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::precision()
+      const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()
     ) const;
     bool isMuchSmallerThan(
       const Scalar& other,
-      const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::precision()
+      const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()
     ) const;
     template<typename OtherDerived>
     bool isMuchSmallerThan(
       const MatrixBase<Scalar, OtherDerived>& other,
-      const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::precision()
+      const typename NumTraits<Scalar>::Real& prec = precision<Scalar>()
     ) const;
     
     template<typename OtherDerived>

src/Core/NumTraits.h

@@ -23,144 +23,176 @@
 // License. This exception does not invalidate any other reasons why a work
 // based on this file might be covered by the GNU General Public License.
 
-#ifndef EIGEN_NUMERIC_H
+#ifndef EIGEN_NUMTRAITS_H
-#define EIGEN_NUMERIC_H
+#define EIGEN_NUMTRAITS_H
 
 template<typename T> struct NumTraits;
 
 template<> struct NumTraits<int>
 {
   typedef int Real;
-  typedef double FloatingPoint;
-  typedef double RealFloatingPoint;
-  
   static const bool IsComplex = false;
   static const bool HasFloatingPoint = false;
-  
-  static int precision() { return 0; }
-  static int real(const int& x) { return x; }
-  static int imag(const int& x) { EIGEN_UNUSED(x); return 0; }
-  static int conj(const int& x) { return x; }
-  static int abs2(const int& x) { return x*x; }
-  static int random()
-  {
-    // "random() % n" is bad, they say, because the low-order bits are not random enough.
-    // However here, 21 is odd, so random() % 21 uses the high-order bits
-    // as well, so there's no problem.
-    return (std::rand() % 21) - 10;
-  }
-  static bool isMuchSmallerThan(const int& a, const int& b, const int& prec = precision())
-  {
-    EIGEN_UNUSED(b);
-    EIGEN_UNUSED(prec);
-    return a == 0;
-  }
-  static bool isApprox(const int& a, const int& b, const int& prec = precision())
-  {
-    EIGEN_UNUSED(prec);
-    return a == b;
-  }
-  static bool isApproxOrLessThan(const int& a, const int& b, const int& prec = precision())
-  {
-    EIGEN_UNUSED(prec);
-    return a <= b;
-  }
 };
 
 template<> struct NumTraits<float>
 {
   typedef float Real;
-  typedef float FloatingPoint;
-  typedef float RealFloatingPoint;
-  
   static const bool IsComplex = false;
   static const bool HasFloatingPoint = true;
-  
-  static float precision() { return 1e-5f; }
-  static float real(const float& x) { return x; }
-  static float imag(const float& x) { EIGEN_UNUSED(x); return 0; }
-  static float conj(const float& x) { return x; }
-  static float abs2(const float& x) { return x*x; }
-  static float random()
-  {
-    return std::rand() / (RAND_MAX/20.0f) - 10.0f;
-  }
-  static bool isMuchSmallerThan(const float& a, const float& b, const float& prec = precision())
-  {
-    return std::abs(a) <= std::abs(b) * prec;
-  }
-  static bool isApprox(const float& a, const float& b, const float& prec = precision())
-  {
-    return std::abs(a - b) <= std::min(std::abs(a), std::abs(b)) * prec;
-  }
-  static bool isApproxOrLessThan(const float& a, const float& b, const float& prec = precision())
-  {
-    return a <= b || isApprox(a, b, prec);
-  }
 };
 
 template<> struct NumTraits<double>
 {
   typedef double Real;
-  typedef double FloatingPoint;
-  typedef double RealFloatingPoint;
-  
   static const bool IsComplex = false;
   static const bool HasFloatingPoint = true;
-  
-  static double precision() { return 1e-11; }
-  static double real(const double& x) { return x; }
-  static double imag(const double& x) { EIGEN_UNUSED(x); return 0; }
-  static double conj(const double& x) { return x; }
-  static double abs2(const double& x) { return x*x; }
-  static double random()
-  {
-    return std::rand() / (RAND_MAX/20.0) - 10.0;
-  }
-  static bool isMuchSmallerThan(const double& a, const double& b, const double& prec = precision())
-  {
-    return std::abs(a) <= std::abs(b) * prec;
-  }
-  static bool isApprox(const double& a, const double& b, const double& prec = precision())
-  {
-    return std::abs(a - b) <= std::min(std::abs(a), std::abs(b)) * prec;
-  }
-  static bool isApproxOrLessThan(const double& a, const double& b, const double& prec = precision())
-  {
-    return a <= b || isApprox(a, b, prec);
-  }
 };
 
 template<typename _Real> struct NumTraits<std::complex<_Real> >
 {
   typedef _Real Real;
-  typedef std::complex<Real> Complex;
-  typedef std::complex<double> FloatingPoint;
-  typedef typename NumTraits<Real>::FloatingPoint RealFloatingPoint;
-  
   static const bool IsComplex = true;
   static const bool HasFloatingPoint = NumTraits<Real>::HasFloatingPoint;
+};
 
-  static Real precision() { return NumTraits<Real>::precision(); }
+template<typename T> inline typename NumTraits<T>::Real precision();
-  static Real real(const Complex& x) { return std::real(x); }
+template<typename T> inline T random();
-  static Real imag(const Complex& x) { return std::imag(x); }
+
-  static Complex conj(const Complex& x) { return std::conj(x); }
+template<> inline int precision<int>() { return 0; }
-  static Real abs2(const Complex& x)
+inline int real(const int& x)  { return x; }
-  { return std::norm(x); }
+inline int imag(const int& x)  { EIGEN_UNUSED(x); return 0; }
-  static Complex random()
+inline int conj(const int& x)  { return x; }
+inline int abs(const int& x)   { return std::abs(x); }
+inline int abs2(const int& x)  { return x*x; }
+inline int sqrt(const int& x)
 {
-    return Complex(NumTraits<Real>::random(), NumTraits<Real>::random());
+  EIGEN_UNUSED(x);
+  // Taking the square root of integers is not allowed
+  // (the square root does not always exist within the integers).
+  // Please cast to a floating-point type.
+  assert(false);
 }
-  static bool isMuchSmallerThan(const Complex& a, const Complex& b, const Real& prec = precision())
+template<> inline int random()
+{
+  // "rand() % n" is bad, they say, because the low-order bits are not random enough.
+  // However here, 21 is odd, so random() % 21 uses the high-order bits
+  // as well, so there's no problem.
+  return (std::rand() % 21) - 10;
+}
+inline bool isMuchSmallerThan(const int& a, const int& b, const int& prec = precision<int>())
+{
+  EIGEN_UNUSED(b);
+  EIGEN_UNUSED(prec);
+  return a == 0;
+}
+inline bool isApprox(const int& a, const int& b, const int& prec = precision<int>())
+{
+  EIGEN_UNUSED(prec);
+  return a == b;
+}
+inline bool isApproxOrLessThan(const int& a, const int& b, const int& prec = precision<int>())
+{
+  EIGEN_UNUSED(prec);
+  return a <= b;
+}
+
+template<> inline float precision<float>() { return 1e-5f; }
+inline float real(const float& x)  { return x; }
+inline float imag(const float& x)  { EIGEN_UNUSED(x); return 0.f; }
+inline float conj(const float& x)  { return x; }
+inline float abs(const float& x)   { return std::abs(x); }
+inline float abs2(const float& x)  { return x*x; }
+inline float sqrt(const float& x)  { return std::sqrt(x); }
+template<> inline float random()
+{
+  return std::rand() / (RAND_MAX/20.0f) - 10.0f;
+}
+inline bool isMuchSmallerThan(const float& a, const float& b, const float& prec = precision<float>())
+{
+  return std::abs(a) <= std::abs(b) * prec;
+}
+inline bool isApprox(const float& a, const float& b, const float& prec = precision<float>())
+{
+  return std::abs(a - b) <= std::min(std::abs(a), std::abs(b)) * prec;
+}
+inline bool isApproxOrLessThan(const float& a, const float& b, const float& prec = precision<float>())
+{
+  return a <= b || isApprox(a, b, prec);
+}
+
+template<> inline double precision<double>() { return 1e-11; }
+inline double real(const double& x)  { return x; }
+inline double imag(const double& x)  { EIGEN_UNUSED(x); return 0.; }
+inline double conj(const double& x)  { return x; }
+inline double abs(const double& x)   { return std::abs(x); }
+inline double abs2(const double& x)  { return x*x; }
+inline double sqrt(const double& x)  { return std::sqrt(x); }
+template<> inline double random()
+{
+  return std::rand() / (RAND_MAX/20.0) - 10.0;
+}
+inline bool isMuchSmallerThan(const double& a, const double& b, const double& prec = precision<double>())
+{
+  return std::abs(a) <= std::abs(b) * prec;
+}
+inline bool isApprox(const double& a, const double& b, const double& prec = precision<double>())
+{
+  return std::abs(a - b) <= std::min(std::abs(a), std::abs(b)) * prec;
+}
+inline bool isApproxOrLessThan(const double& a, const double& b, const double& prec = precision<double>())
+{
+  return a <= b || isApprox(a, b, prec);
+}
+
+template<> inline float precision<std::complex<float> >() { return precision<float>(); }
+inline float real(const std::complex<float>& x) { return std::real(x); }
+inline float imag(const std::complex<float>& x) { return std::imag(x); }
+inline std::complex<float> conj(const std::complex<float>& x) { return std::conj(x); }
+inline float abs(const std::complex<float>& x) { return std::abs(x); }
+inline float abs2(const std::complex<float>& x) { return std::norm(x); }
+inline std::complex<float> sqrt(const std::complex<float>& x)
+{
+  EIGEN_UNUSED(x);
+  // Taking the square roots of complex numbers is not allowed,
+  // as this is ambiguous (there are two square roots).
+  // What were you trying to do?
+  assert(false);
+}
+template<> inline std::complex<float> random()
+{
+  return std::complex<float>(random<float>(), random<float>());
+}
+inline bool isMuchSmallerThan(const std::complex<float>& a, const std::complex<float>& b, const float& prec = precision<float>())
 {
   return abs2(a) <= abs2(b) * prec * prec;
 }
-  static bool isApprox(const Complex& a, const Complex& b, const Real& prec = precision())
+inline bool isApprox(const std::complex<float>& a, const std::complex<float>& b, const float& prec = precision<float>())
 {
-    return NumTraits<Real>::isApprox(std::real(a), std::real(b), prec)
+  return isApprox(std::real(a), std::real(b), prec)
-        && NumTraits<Real>::isApprox(std::imag(a), std::imag(b), prec);
+      && isApprox(std::imag(a), std::imag(b), prec);
 }
 // isApproxOrLessThan wouldn't make sense for complex numbers
-};
 
-#endif // EIGEN_NUMERIC_H
+template<> inline double precision<std::complex<double> >() { return precision<double>(); }
+inline double real(const std::complex<double>& x) { return std::real(x); }
+inline double imag(const std::complex<double>& x) { return std::imag(x); }
+inline std::complex<double> conj(const std::complex<double>& x) { return std::conj(x); }
+inline double abs(const std::complex<double>& x) { return std::abs(x); }
+inline double abs2(const std::complex<double>& x) { return std::norm(x); }
+template<> inline std::complex<double> random()
+{
+  return std::complex<double>(random<double>(), random<double>());
+}
+inline bool isMuchSmallerThan(const std::complex<double>& a, const std::complex<double>& b, const double& prec = precision<double>())
+{
+  return abs2(a) <= abs2(b) * prec * prec;
+}
+inline bool isApprox(const std::complex<double>& a, const std::complex<double>& b, const double& prec = precision<double>())
+{
+  return isApprox(std::real(a), std::real(b), prec)
+      && isApprox(std::imag(a), std::imag(b), prec);
+}
+// isApproxOrLessThan wouldn't make sense for complex numbers
+
+#endif // EIGEN_NUMTRAITS_H

src/Core/Random.h

@@ -53,7 +53,7 @@ template<typename MatrixType> class Random
     {
       EIGEN_UNUSED(row);
       EIGEN_UNUSED(col);
-      return NumTraits<Scalar>::random();
+      return random<Scalar>();
     }
     
   protected:

test/adjoint.cpp

@@ -25,6 +25,8 @@
 
 #include "main.h"
 
+namespace Eigen {
+
 template<typename MatrixType> void adjoint(const MatrixType& m)
 {
   /* this test covers the following files:
@@ -49,8 +51,8 @@ template<typename MatrixType> void adjoint(const MatrixType& m)
              v3 = VectorType::random(rows),
              vzero = VectorType::zero(rows);
 
-  Scalar s1 = NumTraits<Scalar>::random(),
+  Scalar s1 = random<Scalar>(),
-         s2 = NumTraits<Scalar>::random();
+         s2 = random<Scalar>();
   
   // check involutivity of adjoint, transpose, conjugate
   QVERIFY(m1.transpose().transpose().isApprox(m1));
@@ -67,28 +69,32 @@ template<typename MatrixType> void adjoint(const MatrixType& m)
   QVERIFY((m1.adjoint() * m2).adjoint().isApprox(m2.adjoint() * m1));
   QVERIFY((m1.transpose() * m2).conjugate().isApprox(m1.adjoint() * m2.conjugate()));
   QVERIFY((s1 * m1).transpose().isApprox(s1 * m1.transpose()));
-  QVERIFY((s1 * m1).conjugate().isApprox(NumTraits<Scalar>::conj(s1) * m1.conjugate()));
+  QVERIFY((s1 * m1).conjugate().isApprox(conj(s1) * m1.conjugate()));
-  QVERIFY((s1 * m1).adjoint().isApprox(NumTraits<Scalar>::conj(s1) * m1.adjoint()));
+  QVERIFY((s1 * m1).adjoint().isApprox(conj(s1) * m1.adjoint()));
   
   // check basic properties of dot, norm, norm2
   typedef typename NumTraits<Scalar>::Real RealScalar;
-  QVERIFY(NumTraits<Scalar>::isApprox((s1 * v1 + s2 * v2).dot(v3), s1 * v1.dot(v3) + s2 * v2.dot(v3)));
+  QVERIFY(isApprox((s1 * v1 + s2 * v2).dot(v3), s1 * v1.dot(v3) + s2 * v2.dot(v3)));
-  QVERIFY(NumTraits<Scalar>::isApprox(v3.dot(s1 * v1 + s2 * v2), NumTraits<Scalar>::conj(s1) * v3.dot(v1) + NumTraits<Scalar>::conj(s2) * v3.dot(v2)));
+  QVERIFY(isApprox(v3.dot(s1 * v1 + s2 * v2), conj(s1) * v3.dot(v1) + conj(s2) * v3.dot(v2)));
-  QVERIFY(NumTraits<Scalar>::isApprox(NumTraits<Scalar>::conj(v1.dot(v2)), v2.dot(v1)));
+  QVERIFY(isApprox(conj(v1.dot(v2)), v2.dot(v1)));
-  QVERIFY(NumTraits<RealScalar>::isApprox(abs(v1.dot(v1)), v1.norm2()));
+  QVERIFY(isApprox(abs(v1.dot(v1)), v1.norm2()));
-  if(NumTraits<Scalar>::HasFloatingPoint) QVERIFY(NumTraits<RealScalar>::isApprox(v1.norm2(), v1.norm() * v1.norm()));
+  if(NumTraits<Scalar>::HasFloatingPoint)
-  QVERIFY(NumTraits<RealScalar>::isMuchSmallerThan(abs(vzero.dot(v1)), 1));
+    QVERIFY(isApprox(v1.norm2(), v1.norm() * v1.norm()));
-  QVERIFY(NumTraits<RealScalar>::isMuchSmallerThan(vzero.norm(), 1));
+  QVERIFY(isMuchSmallerThan(abs(vzero.dot(v1)), static_cast<RealScalar>(1)));
+  if(NumTraits<Scalar>::HasFloatingPoint)
+    QVERIFY(isMuchSmallerThan(vzero.norm(), static_cast<RealScalar>(1)));
   
   // check compatibility of dot and adjoint
-  QVERIFY(NumTraits<Scalar>::isApprox(v1.dot(square * v2), (square.adjoint() * v1).dot(v2)));
+  QVERIFY(isApprox(v1.dot(square * v2), (square.adjoint() * v1).dot(v2)));
 }
 
 void EigenTest::testAdjoint()
 {
   adjoint(Matrix<float, 1, 1>());
-  adjoint(Matrix<complex<double>, 4, 4>());
+  adjoint(Matrix4cd());
   adjoint(MatrixXcf(3, 3));
   adjoint(MatrixXi(8, 12));
   adjoint(MatrixXd(20, 20));
 }
+
+} // namespace Eigen

test/basicstuff.cpp

@@ -25,6 +25,8 @@
 
 #include "main.h"
 
+namespace Eigen {
+
 template<typename MatrixType> void basicStuff(const MatrixType& m)
 {
   /* this test covers the following files:
@@ -56,13 +58,14 @@ template<typename MatrixType> void basicStuff(const MatrixType& m)
              v2 = VectorType::random(rows),
              vzero = VectorType::zero(rows);
 
-  Scalar s1 = NumTraits<Scalar>::random(),
+  Scalar s1 = random<Scalar>(),
-         s2 = NumTraits<Scalar>::random();
+         s2 = random<Scalar>();
   
   // test Fuzzy.h and Zero.h.
   QVERIFY(v1.isApprox(v1));
   QVERIFY(!v1.isApprox(2*v1));
   QVERIFY(vzero.isMuchSmallerThan(v1));
+  if(NumTraits<Scalar>::HasFloatingPoint)
     QVERIFY(vzero.isMuchSmallerThan(v1.norm()));
   QVERIFY(!v1.isMuchSmallerThan(v1));
   QVERIFY(vzero.isApprox(v1-v1));
@@ -134,8 +137,10 @@ template<typename MatrixType> void basicStuff(const MatrixType& m)
 void EigenTest::testBasicStuff()
 {
   basicStuff(Matrix<float, 1, 1>());
-  basicStuff(Matrix<complex<double>, 4, 4>());
+  basicStuff(Matrix4cd());
   basicStuff(MatrixXcf(3, 3));
   basicStuff(MatrixXi(8, 12));
   basicStuff(MatrixXd(20, 20));
 }
+
+} // namespace Eigen

test/main.cpp

@@ -25,7 +25,7 @@
 
 #include "main.h"
 
-EigenTest::EigenTest()
+Eigen::EigenTest::EigenTest()
 {
     unsigned int t = (unsigned int) time(NULL);
     qDebug() << "Initializing random number generator with seed"
@@ -33,5 +33,6 @@ EigenTest::EigenTest()
     srand(t);
 }
 
-QTEST_APPLESS_MAIN( EigenTest )
+QTEST_APPLESS_MAIN(Eigen::EigenTest)
+
 #include "main.moc"

test/main.h

@@ -29,12 +29,10 @@
 #include <QtTest/QtTest>
 #include "../src/Core.h"
 
-using namespace Eigen;
-
 #include <cstdlib>
 #include <ctime>
 
-using namespace std;
+namespace Eigen {
 
 class EigenTest : public QObject
 {
@@ -48,4 +46,6 @@ class EigenTest : public QObject
     void testAdjoint();
 };
 
+} // end namespace Eigen
+
 #endif // EIGEN_TEST_MAIN_H