diff --git a/src/Core.h b/src/Core.h index 874a8f5e1..658f7fe08 100644 --- a/src/Core.h +++ b/src/Core.h @@ -13,3 +13,4 @@ #include "Core/Conjugate.h" #include "Core/Trace.h" #include "Core/Dot.h" +#include "Core/Random.h" diff --git a/src/Core/Dot.h b/src/Core/Dot.h index 175c3e57d..0d000ea70 100644 --- a/src/Core/Dot.h +++ b/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::run(v1, v2, dot); - if(i == Size - 1) - dot = v1[i] * EiConj(v2[i]); - else - dot += v1[i] * EiConj(v2[i]); + dot += v1[Index-1] * EiConj(v2[Index-1]); } }; template -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); - EI_UNUSED(v2); - EI_UNUSED(dot); + dot = v1[0] * EiConj(v2[0]); } }; @@ -80,4 +74,23 @@ Scalar EiObject::dot(const OtherDerived& other) const return res; } +template +typename EiNumTraits::Real EiObject::norm2() const +{ + assert(IsVector); + return EiReal(dot(*this)); +} + +template +typename EiNumTraits::Real EiObject::norm() const +{ + return EiSqrt(norm2()); +} + +template +EiScalarProduct EiObject::normalized() const +{ + return (*this) / norm(); +} + #endif // EI_DOT_H diff --git a/src/Core/Fuzzy.h b/src/Core/Fuzzy.h new file mode 100644 index 000000000..44b20615d --- /dev/null +++ b/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 +// +// 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 +template +bool EiObject::isApprox(const OtherDerived& other) const +{ + if(IsVector) + { + return((*this - other).norm2() + <= std::min(norm2(), other.norm2()) + * EiAbs2(EiNumTraits::epsilon())); + } + else + { + for(int i = 0; i < cols(); i++) + if(!col(i).isApprox(other.col(i))) + return false; + return true; + } +} + +template +bool EiObject::isNegligble(const Scalar& other) const +{ + if(IsVector) + { + return(norm2() <= EiAbs2(other) * EiAbs2(EiNumTraits::epsilon())); + } + else + { + for(int i = 0; i < cols(); i++) + if(!col(i).isNegligible(other)) + return false; + return true; + } +} + +template +template +bool EiObject::isNegligble(const EiObject& other) const +{ + if(IsVector) + { + return(norm2() <= other.norm2() * EiAbs2(EiNumTraits::epsilon())); + } + else + { + for(int i = 0; i < cols(); i++) + if(!col(i).isNegligible(other.col(i))) + return false; + return true; + } +} + +#endif // EI_FUZZY_H \ No newline at end of file diff --git a/src/Core/Numeric.h b/src/Core/Numeric.h index 48a36f236..23a6c03e6 100644 --- a/src/Core/Numeric.h +++ b/src/Core/Numeric.h @@ -27,9 +27,9 @@ #ifndef EI_NUMERIC_H #define EI_NUMERIC_H -template struct EiTraits; +template struct EiNumTraits; -template<> struct EiTraits +template<> struct EiNumTraits { typedef int Real; typedef double FloatingPoint; @@ -45,16 +45,16 @@ template<> struct EiTraits static double sqrt(const int& x) { return std::sqrt(static_cast(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 seems to round the nearest int. // I don't know if that's part of the standard. - return -10 + static_cast(rand() / ((RAND_MAX + 1.0)/20.0)); + return -10 + static_cast(std::rand() / ((RAND_MAX + 1.0)/20.0)); } }; -template<> struct EiTraits +template<> struct EiNumTraits { typedef float Real; typedef float FloatingPoint; @@ -70,13 +70,13 @@ template<> struct EiTraits 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 +template<> struct EiNumTraits { typedef double Real; typedef double FloatingPoint; @@ -92,23 +92,23 @@ template<> struct EiTraits 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 struct EiTraits > +template struct EiNumTraits > { typedef _Real Real; typedef std::complex Complex; typedef std::complex FloatingPoint; - typedef typename EiTraits::FloatingPoint RealFloatingPoint; + typedef typename EiNumTraits::FloatingPoint RealFloatingPoint; static const bool IsComplex = true; - static const bool HasFloatingPoint = EiTraits::HasFloatingPoint; + static const bool HasFloatingPoint = EiNumTraits::HasFloatingPoint; - static Real epsilon() { return EiTraits::epsilon(); } + static Real epsilon() { return EiNumTraits::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 struct EiTraits > { return std::abs(static_cast(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::random(), EiTraits::random()); + return Complex(EiNumTraits::rand(), EiNumTraits::rand()); } }; -template typename EiTraits::Real EiReal(const T& x) -{ return EiTraits::real(x); } +template typename EiNumTraits::Real EiReal(const T& x) +{ return EiNumTraits::real(x); } -template typename EiTraits::Real EiImag(const T& x) -{ return EiTraits::imag(x); } +template typename EiNumTraits::Real EiImag(const T& x) +{ return EiNumTraits::imag(x); } template T EiConj(const T& x) -{ return EiTraits::conj(x); } +{ return EiNumTraits::conj(x); } -template typename EiTraits::FloatingPoint EiSqrt(const T& x) -{ return EiTraits::sqrt(x); } +template typename EiNumTraits::FloatingPoint EiSqrt(const T& x) +{ return EiNumTraits::sqrt(x); } -template typename EiTraits::RealFloatingPoint EiAbs(const T& x) -{ return EiTraits::abs(x); } +template typename EiNumTraits::RealFloatingPoint EiAbs(const T& x) +{ return EiNumTraits::abs(x); } -template typename EiTraits::Real EiAbs2(const T& x) -{ return EiTraits::abs2(x); } +template typename EiNumTraits::Real EiAbs2(const T& x) +{ return EiNumTraits::abs2(x); } -template T EiRandom() -{ return EiTraits::random(); } +template T EiRand() +{ return EiNumTraits::rand(); } template bool EiNegligible(const T& a, const T& b) { - return(EiAbs(a) <= EiAbs(b) * EiTraits::epsilon()); + return(EiAbs(a) <= EiAbs(b) * EiNumTraits::epsilon()); } template bool EiApprox(const T& a, const T& b) { - if(EiTraits::IsFloat) - return(EiAbs(a - b) <= std::min(EiAbs(a), EiAbs(b)) * EiTraits::epsilon()); + if(EiNumTraits::IsFloat) + return(EiAbs(a - b) <= std::min(EiAbs(a), EiAbs(b)) * EiNumTraits::epsilon()); else return(a == b); } template bool EiLessThanOrApprox(const T& a, const T& b) { - if(EiTraits::IsFloat) + if(EiNumTraits::IsFloat) return(a < b || EiApprox(a, b)); else return(a <= b); diff --git a/src/Core/Object.h b/src/Core/Object.h index 89c62b7b5..d1a310a20 100644 --- a/src/Core/Object.h +++ b/src/Core/Object.h @@ -42,6 +42,7 @@ template class EiObject typedef typename EiForwardDecl::Ref Ref; typedef typename EiForwardDecl::ConstRef ConstRef; + typedef typename EiNumTraits::Real RealScalar; int rows() const { return static_cast(this)->_rows(); } int cols() const { return static_cast(this)->_cols(); } @@ -92,8 +93,17 @@ template class EiObject template Scalar dot(const OtherDerived& other) const; - Scalar norm2() const { assert(IsVector); return dot(*this); } - Scalar norm() const { assert(IsVector); return EiSqrt(dot(*this)); } + RealScalar norm2() const; + RealScalar norm() const; + EiScalarProduct normalized() const; + + static EiEval > + random(int rows = RowsAtCompileTime, int cols = ColsAtCompileTime); + + template + bool isApprox(const OtherDerived& other) const; + template + bool isNegligible(const OtherDerived& other) const; template EiMatrixProduct diff --git a/src/Core/Random.h b/src/Core/Random.h new file mode 100644 index 000000000..70485028f --- /dev/null +++ b/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 +// +// 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 class EiRandom + : public EiObject > +{ + public: + typedef typename MatrixType::Scalar Scalar; + friend class EiObject >; + + 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(); + } + + protected: + int m_rows, m_cols; +}; + +template +EiEval > EiObject::random(int rows, int cols) +{ + return EiRandom(rows, cols).eval(); +} + +#endif // EI_RANDOM_H diff --git a/src/Core/ScalarOps.h b/src/Core/ScalarOps.h index 2f56ff3b4..c058978e6 100644 --- a/src/Core/ScalarOps.h +++ b/src/Core/ScalarOps.h @@ -83,7 +83,7 @@ EiScalarProduct \ operator/(const EiObject& matrix, \ OtherScalar scalar) \ { \ - assert(EiTraits::HasFloatingPoint); \ + assert(EiNumTraits::HasFloatingPoint); \ return matrix * (static_cast(1) / scalar); \ } \ \ diff --git a/src/Core/Util.h b/src/Core/Util.h index 859913908..460eecb1a 100644 --- a/src/Core/Util.h +++ b/src/Core/Util.h @@ -54,6 +54,7 @@ template class EiSum; template class EiDifference; template class EiMatrixProduct; template class EiScalarProduct; +template class EiRandom; template class EiEval; template struct EiForwardDecl diff --git a/test/main.h b/test/main.h index 91d63c7f8..05853481c 100644 --- a/test/main.h +++ b/test/main.h @@ -47,7 +47,7 @@ class EigenTest : public QObject void testMatrixManip(); }; -template inline typename EiTraits::Real TestEpsilon(); +template inline typename EiNumTraits::Real TestEpsilon(); template<> inline int TestEpsilon() { return 0; } template<> inline float TestEpsilon() { return 1e-2f; } template<> inline double TestEpsilon() { return 1e-4; } @@ -62,7 +62,7 @@ template bool TestNegligible(const T& a, const T& b) template bool TestApprox(const T& a, const T& b) { - if(EiTraits::IsFloat) + if(EiNumTraits::IsFloat) return(EiAbs(a - b) <= std::min(EiAbs(a), EiAbs(b)) * TestEpsilon()); else return(a == b); @@ -70,7 +70,7 @@ template bool TestApprox(const T& a, const T& b) template bool TestLessThanOrApprox(const T& a, const T& b) { - if(EiTraits::IsFloat) + if(EiNumTraits::IsFloat) return(a < b || EiApprox(a, b)); else return(a <= b);