fix dot product, add norm/norm2/normalized

add fuzzy compares for matrices/vectors
add random matrix/vector generation

src/Core.h

@@ -13,3 +13,4 @@
 #include "Core/Conjugate.h"
 #include "Core/Trace.h"
 #include "Core/Dot.h"
+#include "Core/Random.h"

src/Core/Dot.h

@@ -31,23 +31,17 @@ struct EiDotUnroller
 {
   static void run(const Derived1 &v1, const Derived2& v2, typename Derived1::Scalar &dot)
   {
-    const int i = Index - 1;
     EiDotUnroller<Index-1, Size, Derived1, Derived2>::run(v1, v2, dot);
-    if(i == Size - 1)
+    dot += v1[Index-1] * EiConj(v2[Index-1]);
-      dot = v1[i] * EiConj(v2[i]);
-    else
-      dot += v1[i] * EiConj(v2[i]);
   }
 };
 
 template<int Size, typename Derived1, typename Derived2>
-struct EiDotUnroller<0, Size, Derived1, Derived2>
+struct EiDotUnroller<1, Size, Derived1, Derived2>
 {
   static void run(const Derived1 &v1, const Derived2& v2, typename Derived1::Scalar &dot)
   {
-    EI_UNUSED(v1);
+    dot = v1[0] * EiConj(v2[0]);
-    EI_UNUSED(v2);
-    EI_UNUSED(dot);
   }
 };
 
@@ -80,4 +74,23 @@ Scalar EiObject<Scalar, Derived>::dot(const OtherDerived& other) const
   return res;
 }
 
+template<typename Scalar, typename Derived>
+typename EiNumTraits<Scalar>::Real EiObject<Scalar, Derived>::norm2() const
+{
+  assert(IsVector);
+  return EiReal(dot(*this));
+}
+
+template<typename Scalar, typename Derived>
+typename EiNumTraits<Scalar>::Real EiObject<Scalar, Derived>::norm() const
+{
+  return EiSqrt(norm2());
+}
+
+template<typename Scalar, typename Derived>
+EiScalarProduct<Derived> EiObject<Scalar, Derived>::normalized() const
+{
+  return (*this) / norm();
+}
+
 #endif // EI_DOT_H

src/Core/Fuzzy.h

@@ -0,0 +1,81 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2006-2007 Benoit Jacob <jacob@math.jussieu.fr>
+//
+// Eigen is free software; you can redistribute it and/or modify it under the
+// terms of the GNU General Public License as published by the Free Software
+// Foundation; either version 2 or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+// details.
+//
+// You should have received a copy of the GNU General Public License along
+// with Eigen; if not, write to the Free Software Foundation, Inc., 51
+// Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
+//
+// As a special exception, if other files instantiate templates or use macros
+// or functions from this file, or you compile this file and link it
+// with other works to produce a work based on this file, this file does not
+// by itself cause the resulting work to be covered by the GNU General Public
+// 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 EI_FUZZY_H
+#define EI_FUZZY_H
+
+template<typename Scalar, typename Derived>
+template<typename OtherDerived>
+bool EiObject<Scalar, Derived>::isApprox(const OtherDerived& other) const
+{
+  if(IsVector)
+  {
+    return((*this - other).norm2()
+           <= std::min(norm2(), other.norm2())
+            * EiAbs2(EiNumTraits<Scalar>::epsilon()));
+  }
+  else
+  {
+    for(int i = 0; i < cols(); i++)
+      if(!col(i).isApprox(other.col(i)))
+        return false;
+    return true;
+  }
+}
+
+template<typename Scalar, typename Derived>
+bool EiObject<Scalar, Derived>::isNegligble(const Scalar& other) const
+{
+  if(IsVector)
+  {
+    return(norm2() <= EiAbs2(other) * EiAbs2(EiNumTraits<Scalar>::epsilon()));
+  }
+  else
+  {
+    for(int i = 0; i < cols(); i++)
+      if(!col(i).isNegligible(other))
+        return false;
+    return true;
+  }
+}
+
+template<typename Scalar, typename Derived>
+template<typename OtherDerived>
+bool EiObject<Scalar, Derived>::isNegligble(const EiObject<Scalar, OtherDerived>& other) const
+{
+  if(IsVector)
+  {
+    return(norm2() <= other.norm2() * EiAbs2(EiNumTraits<Scalar>::epsilon()));
+  }
+  else
+  {
+    for(int i = 0; i < cols(); i++)
+      if(!col(i).isNegligible(other.col(i)))
+        return false;
+    return true;
+  }
+}
+
+#endif // EI_FUZZY_H

src/Core/Numeric.h

@@ -27,9 +27,9 @@
 #ifndef EI_NUMERIC_H
 #define EI_NUMERIC_H
 
-template<typename T> struct EiTraits;
+template<typename T> struct EiNumTraits;
 
-template<> struct EiTraits<int>
+template<> struct EiNumTraits<int>
 {
   typedef int Real;
   typedef double FloatingPoint;
@@ -45,16 +45,16 @@ template<> struct EiTraits<int>
   static double sqrt(const int& x) { return std::sqrt(static_cast<double>(x)); }
   static int abs(const int& x) { return std::abs(x); }
   static int abs2(const int& x) { return x*x; }
-  static int random()
+  static int rand()
   {
     // "rand()%21" would be bad. always use the high-order bits, not the low-order bits.
     // note: here (gcc 4.1) static_cast<int> seems to round the nearest int.
     // I don't know if that's part of the standard.
-    return -10 + static_cast<int>(rand() / ((RAND_MAX + 1.0)/20.0));
+    return -10 + static_cast<int>(std::rand() / ((RAND_MAX + 1.0)/20.0));
   }
 };
 
-template<> struct EiTraits<float>
+template<> struct EiNumTraits<float>
 {
   typedef float Real;
   typedef float FloatingPoint;
@@ -70,13 +70,13 @@ template<> struct EiTraits<float>
   static float sqrt(const float& x) { return std::sqrt(x); }
   static float abs(const float& x) { return std::abs(x); }
   static float abs2(const float& x) { return x*x; }
-  static float random()
+  static float rand()
   {
-    return rand() / (RAND_MAX/20.0f) - 10.0f;
+    return std::rand() / (RAND_MAX/20.0f) - 10.0f;
   }
 };
 
-template<> struct EiTraits<double>
+template<> struct EiNumTraits<double>
 {
   typedef double Real;
   typedef double FloatingPoint;
@@ -92,23 +92,23 @@ template<> struct EiTraits<double>
   static double sqrt(const double& x) { return std::sqrt(x); }
   static double abs(const double& x) { return std::abs(x); }
   static double abs2(const double& x) { return x*x; }
-  static double random()
+  static double rand()
   {
-    return rand() / (RAND_MAX/20.0) - 10.0;
+    return std::rand() / (RAND_MAX/20.0) - 10.0;
   }
 };
 
-template<typename _Real> struct EiTraits<std::complex<_Real> >
+template<typename _Real> struct EiNumTraits<std::complex<_Real> >
 {
   typedef _Real Real;
   typedef std::complex<Real> Complex;
   typedef std::complex<double> FloatingPoint;
-  typedef typename EiTraits<Real>::FloatingPoint RealFloatingPoint;
+  typedef typename EiNumTraits<Real>::FloatingPoint RealFloatingPoint;
   
   static const bool IsComplex = true;
-  static const bool HasFloatingPoint = EiTraits<Real>::HasFloatingPoint;
+  static const bool HasFloatingPoint = EiNumTraits<Real>::HasFloatingPoint;
   
-  static Real epsilon() { return EiTraits<Real>::epsilon(); }
+  static Real epsilon() { return EiNumTraits<Real>::epsilon(); }
   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); }
@@ -118,49 +118,49 @@ template<typename _Real> struct EiTraits<std::complex<_Real> >
   { return std::abs(static_cast<FloatingPoint>(x)); }
   static Real abs2(const Complex& x)
   { return std::real(x) * std::real(x) + std::imag(x) * std::imag(x); }
-  static Complex random()
+  static Complex rand()
   {
-    return Complex(EiTraits<Real>::random(), EiTraits<Real>::random());
+    return Complex(EiNumTraits<Real>::rand(), EiNumTraits<Real>::rand());
   }
 };
 
-template<typename T> typename EiTraits<T>::Real EiReal(const T& x)
+template<typename T> typename EiNumTraits<T>::Real EiReal(const T& x)
-{ return EiTraits<T>::real(x); }
+{ return EiNumTraits<T>::real(x); }
 
-template<typename T> typename EiTraits<T>::Real EiImag(const T& x)
+template<typename T> typename EiNumTraits<T>::Real EiImag(const T& x)
-{ return EiTraits<T>::imag(x); }
+{ return EiNumTraits<T>::imag(x); }
 
 template<typename T> T EiConj(const T& x)
-{ return EiTraits<T>::conj(x); }
+{ return EiNumTraits<T>::conj(x); }
 
-template<typename T> typename EiTraits<T>::FloatingPoint EiSqrt(const T& x)
+template<typename T> typename EiNumTraits<T>::FloatingPoint EiSqrt(const T& x)
-{ return EiTraits<T>::sqrt(x); }
+{ return EiNumTraits<T>::sqrt(x); }
 
-template<typename T> typename EiTraits<T>::RealFloatingPoint EiAbs(const T& x)
+template<typename T> typename EiNumTraits<T>::RealFloatingPoint EiAbs(const T& x)
-{ return EiTraits<T>::abs(x); }
+{ return EiNumTraits<T>::abs(x); }
 
-template<typename T> typename EiTraits<T>::Real EiAbs2(const T& x)
+template<typename T> typename EiNumTraits<T>::Real EiAbs2(const T& x)
-{ return EiTraits<T>::abs2(x); }
+{ return EiNumTraits<T>::abs2(x); }
 
-template<typename T> T EiRandom()
+template<typename T> T EiRand()
-{ return EiTraits<T>::random(); }
+{ return EiNumTraits<T>::rand(); }
 
 template<typename T> bool EiNegligible(const T& a, const T& b)
 {
-  return(EiAbs(a) <= EiAbs(b) * EiTraits<T>::epsilon());
+  return(EiAbs(a) <= EiAbs(b) * EiNumTraits<T>::epsilon());
 }
 
 template<typename T> bool EiApprox(const T& a, const T& b)
 {
-  if(EiTraits<T>::IsFloat)
+  if(EiNumTraits<T>::IsFloat)
-    return(EiAbs(a - b) <= std::min(EiAbs(a), EiAbs(b)) * EiTraits<T>::epsilon());
+    return(EiAbs(a - b) <= std::min(EiAbs(a), EiAbs(b)) * EiNumTraits<T>::epsilon());
   else
     return(a == b);
 }
 
 template<typename T> bool EiLessThanOrApprox(const T& a, const T& b)
 {
-  if(EiTraits<T>::IsFloat)
+  if(EiNumTraits<T>::IsFloat)
     return(a < b || EiApprox(a, b));
   else
     return(a <= b);

src/Core/Object.h

@@ -42,6 +42,7 @@ template<typename Scalar, typename Derived> class EiObject
     
     typedef typename EiForwardDecl<Derived>::Ref Ref;
     typedef typename EiForwardDecl<Derived>::ConstRef ConstRef;
+    typedef typename EiNumTraits<Scalar>::Real RealScalar;
   
     int rows() const { return static_cast<const Derived *>(this)->_rows(); }
     int cols() const { return static_cast<const Derived *>(this)->_cols(); }
@@ -92,8 +93,17 @@ template<typename Scalar, typename Derived> class EiObject
     template<typename OtherDerived>
     Scalar dot(const OtherDerived& other) const;
     
-    Scalar norm2() const { assert(IsVector); return dot(*this); }
+    RealScalar norm2() const;
-    Scalar norm()  const { assert(IsVector); return EiSqrt(dot(*this)); }
+    RealScalar norm()  const;
+    EiScalarProduct<Derived> normalized() const;
+    
+    static EiEval<EiRandom<Derived> >
+    random(int rows = RowsAtCompileTime, int cols = ColsAtCompileTime);
+    
+    template<typename OtherDerived>
+    bool isApprox(const OtherDerived& other) const;
+    template<typename OtherDerived>
+    bool isNegligible(const OtherDerived& other) const;
     
     template<typename OtherDerived>
     EiMatrixProduct<Derived, OtherDerived>

src/Core/Random.h

@@ -0,0 +1,67 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2006-2007 Benoit Jacob <jacob@math.jussieu.fr>
+//
+// Eigen is free software; you can redistribute it and/or modify it under the
+// terms of the GNU General Public License as published by the Free Software
+// Foundation; either version 2 or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+// details.
+//
+// You should have received a copy of the GNU General Public License along
+// with Eigen; if not, write to the Free Software Foundation, Inc., 51
+// Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
+//
+// As a special exception, if other files instantiate templates or use macros
+// or functions from this file, or you compile this file and link it
+// with other works to produce a work based on this file, this file does not
+// by itself cause the resulting work to be covered by the GNU General Public
+// 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 EI_RANDOM_H
+#define EI_RANDOM_H
+
+template<typename MatrixType> class EiRandom
+  : public EiObject<typename MatrixType::Scalar, EiRandom<MatrixType> >
+{
+  public:
+    typedef typename MatrixType::Scalar Scalar;
+    friend class EiObject<Scalar, EiRandom<MatrixType> >;
+    
+    static const int RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+                     ColsAtCompileTime = MatrixType::ColsAtCompileTime;
+
+    EiRandom(int rows, int cols) : m_rows(rows), m_cols(cols)
+    {
+      assert(rows > 0 && cols > 0);
+    }
+    
+  private:
+    EiRandom& _ref() { return *this; }
+    const EiRandom& _constRef() const { return *this; }
+    int _rows() const { return m_rows; }
+    int _cols() const { return m_cols; }
+    
+    Scalar _read(int row, int col) const
+    {
+      EI_UNUSED(row);
+      EI_UNUSED(col);
+      return EiRand<Scalar>();
+    }
+    
+  protected:
+    int m_rows, m_cols;
+};
+
+template<typename Scalar, typename Derived>
+EiEval<EiRandom<Derived> > EiObject<Scalar, Derived>::random(int rows, int cols)
+{
+  return EiRandom<Derived>(rows, cols).eval();
+}
+
+#endif // EI_RANDOM_H

src/Core/ScalarOps.h

@@ -83,7 +83,7 @@ EiScalarProduct<Derived>                                               \
 operator/(const EiObject<Scalar, Derived>& matrix,                     \
           OtherScalar scalar)                                          \
 {                                                                      \
-  assert(EiTraits<Scalar>::HasFloatingPoint);                          \
+  assert(EiNumTraits<Scalar>::HasFloatingPoint);                       \
   return matrix * (static_cast<Scalar>(1) / scalar);                   \
 }                                                                      \
                                                                        \

src/Core/Util.h

@@ -54,6 +54,7 @@ template<typename Lhs, typename Rhs> class EiSum;
 template<typename Lhs, typename Rhs> class EiDifference;
 template<typename Lhs, typename Rhs> class EiMatrixProduct;
 template<typename MatrixType> class EiScalarProduct;
+template<typename MatrixType> class EiRandom;
 template<typename ExpressionType> class EiEval;
 
 template<typename T> struct EiForwardDecl

test/main.h

@@ -47,7 +47,7 @@ class EigenTest : public QObject
     void testMatrixManip();
 };
 
-template<typename T> inline typename EiTraits<T>::Real TestEpsilon();
+template<typename T> inline typename EiNumTraits<T>::Real TestEpsilon();
 template<> inline int TestEpsilon<int>() { return 0; }
 template<> inline float TestEpsilon<float>() { return 1e-2f; }
 template<> inline double TestEpsilon<double>() { return 1e-4; }
@@ -62,7 +62,7 @@ template<typename T> bool TestNegligible(const T& a, const T& b)
 
 template<typename T> bool TestApprox(const T& a, const T& b)
 {
-  if(EiTraits<T>::IsFloat)
+  if(EiNumTraits<T>::IsFloat)
     return(EiAbs(a - b) <= std::min(EiAbs(a), EiAbs(b)) * TestEpsilon<T>());
   else
     return(a == b);
@@ -70,7 +70,7 @@ template<typename T> bool TestApprox(const T& a, const T& b)
 
 template<typename T> bool TestLessThanOrApprox(const T& a, const T& b)
 {
-  if(EiTraits<T>::IsFloat)
+  if(EiNumTraits<T>::IsFloat)
     return(a < b || EiApprox(a, b));
   else
     return(a <= b);