2007-10-07 20:44:42 +08:00
|
|
|
|
|
|
|
// 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.
|
|
|
|
|
2007-10-07 23:12:41 +08:00
|
|
|
#ifndef EI_NUMERIC_H
|
|
|
|
#define EI_NUMERIC_H
|
2007-10-07 20:44:42 +08:00
|
|
|
|
2007-10-10 14:09:56 +08:00
|
|
|
template<typename T> struct EiNumTraits;
|
2007-10-07 20:44:42 +08:00
|
|
|
|
2007-10-10 14:09:56 +08:00
|
|
|
template<> struct EiNumTraits<int>
|
2007-10-07 20:44:42 +08:00
|
|
|
{
|
|
|
|
typedef int Real;
|
|
|
|
typedef double FloatingPoint;
|
|
|
|
typedef double RealFloatingPoint;
|
|
|
|
|
|
|
|
static const bool IsComplex = false;
|
|
|
|
static const bool HasFloatingPoint = false;
|
|
|
|
|
2007-10-07 23:12:41 +08:00
|
|
|
static int epsilon() { return 0; }
|
2007-10-07 20:44:42 +08:00
|
|
|
static int real(const int& x) { return x; }
|
|
|
|
static int imag(const int& x) { EI_UNUSED(x); return 0; }
|
|
|
|
static int conj(const int& x) { return x; }
|
|
|
|
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; }
|
2007-10-10 14:09:56 +08:00
|
|
|
static int rand()
|
2007-10-07 20:44:42 +08:00
|
|
|
{
|
|
|
|
// "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.
|
2007-10-10 14:09:56 +08:00
|
|
|
return -10 + static_cast<int>(std::rand() / ((RAND_MAX + 1.0)/20.0));
|
2007-10-07 20:44:42 +08:00
|
|
|
}
|
|
|
|
};
|
|
|
|
|
2007-10-10 14:09:56 +08:00
|
|
|
template<> struct EiNumTraits<float>
|
2007-10-07 20:44:42 +08:00
|
|
|
{
|
|
|
|
typedef float Real;
|
|
|
|
typedef float FloatingPoint;
|
|
|
|
typedef float RealFloatingPoint;
|
|
|
|
|
|
|
|
static const bool IsComplex = false;
|
|
|
|
static const bool HasFloatingPoint = true;
|
|
|
|
|
2007-10-07 23:12:41 +08:00
|
|
|
static float epsilon() { return 1e-5f; }
|
2007-10-07 20:44:42 +08:00
|
|
|
static float real(const float& x) { return x; }
|
|
|
|
static float imag(const float& x) { EI_UNUSED(x); return 0; }
|
|
|
|
static float conj(const float& x) { return x; }
|
|
|
|
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; }
|
2007-10-10 14:09:56 +08:00
|
|
|
static float rand()
|
2007-10-07 20:44:42 +08:00
|
|
|
{
|
2007-10-10 14:09:56 +08:00
|
|
|
return std::rand() / (RAND_MAX/20.0f) - 10.0f;
|
2007-10-07 20:44:42 +08:00
|
|
|
}
|
|
|
|
};
|
|
|
|
|
2007-10-10 14:09:56 +08:00
|
|
|
template<> struct EiNumTraits<double>
|
2007-10-07 20:44:42 +08:00
|
|
|
{
|
|
|
|
typedef double Real;
|
|
|
|
typedef double FloatingPoint;
|
|
|
|
typedef double RealFloatingPoint;
|
|
|
|
|
|
|
|
static const bool IsComplex = false;
|
|
|
|
static const bool HasFloatingPoint = true;
|
|
|
|
|
2007-10-07 23:12:41 +08:00
|
|
|
static double epsilon() { return 1e-11; }
|
2007-10-07 20:44:42 +08:00
|
|
|
static double real(const double& x) { return x; }
|
|
|
|
static double imag(const double& x) { EI_UNUSED(x); return 0; }
|
|
|
|
static double conj(const double& x) { return x; }
|
|
|
|
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; }
|
2007-10-10 14:09:56 +08:00
|
|
|
static double rand()
|
2007-10-07 20:44:42 +08:00
|
|
|
{
|
2007-10-10 14:09:56 +08:00
|
|
|
return std::rand() / (RAND_MAX/20.0) - 10.0;
|
2007-10-07 20:44:42 +08:00
|
|
|
}
|
|
|
|
};
|
|
|
|
|
2007-10-10 14:09:56 +08:00
|
|
|
template<typename _Real> struct EiNumTraits<std::complex<_Real> >
|
2007-10-07 20:44:42 +08:00
|
|
|
{
|
|
|
|
typedef _Real Real;
|
|
|
|
typedef std::complex<Real> Complex;
|
|
|
|
typedef std::complex<double> FloatingPoint;
|
2007-10-10 14:09:56 +08:00
|
|
|
typedef typename EiNumTraits<Real>::FloatingPoint RealFloatingPoint;
|
2007-10-07 20:44:42 +08:00
|
|
|
|
|
|
|
static const bool IsComplex = true;
|
2007-10-10 14:09:56 +08:00
|
|
|
static const bool HasFloatingPoint = EiNumTraits<Real>::HasFloatingPoint;
|
2007-10-07 20:44:42 +08:00
|
|
|
|
2007-10-10 14:09:56 +08:00
|
|
|
static Real epsilon() { return EiNumTraits<Real>::epsilon(); }
|
2007-10-07 20:44:42 +08:00
|
|
|
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); }
|
|
|
|
static FloatingPoint sqrt(const Complex& x)
|
|
|
|
{ return std::sqrt(static_cast<FloatingPoint>(x)); }
|
|
|
|
static RealFloatingPoint abs(const Complex& x)
|
|
|
|
{ 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); }
|
2007-10-10 14:09:56 +08:00
|
|
|
static Complex rand()
|
2007-10-07 20:44:42 +08:00
|
|
|
{
|
2007-10-10 14:09:56 +08:00
|
|
|
return Complex(EiNumTraits<Real>::rand(), EiNumTraits<Real>::rand());
|
2007-10-07 20:44:42 +08:00
|
|
|
}
|
|
|
|
};
|
|
|
|
|
2007-10-10 14:09:56 +08:00
|
|
|
template<typename T> typename EiNumTraits<T>::Real EiReal(const T& x)
|
|
|
|
{ return EiNumTraits<T>::real(x); }
|
2007-10-07 20:44:42 +08:00
|
|
|
|
2007-10-10 14:09:56 +08:00
|
|
|
template<typename T> typename EiNumTraits<T>::Real EiImag(const T& x)
|
|
|
|
{ return EiNumTraits<T>::imag(x); }
|
2007-10-07 20:44:42 +08:00
|
|
|
|
|
|
|
template<typename T> T EiConj(const T& x)
|
2007-10-10 14:09:56 +08:00
|
|
|
{ return EiNumTraits<T>::conj(x); }
|
2007-10-07 20:44:42 +08:00
|
|
|
|
2007-10-10 14:09:56 +08:00
|
|
|
template<typename T> typename EiNumTraits<T>::FloatingPoint EiSqrt(const T& x)
|
|
|
|
{ return EiNumTraits<T>::sqrt(x); }
|
2007-10-07 20:44:42 +08:00
|
|
|
|
2007-10-10 14:09:56 +08:00
|
|
|
template<typename T> typename EiNumTraits<T>::RealFloatingPoint EiAbs(const T& x)
|
|
|
|
{ return EiNumTraits<T>::abs(x); }
|
2007-10-07 20:44:42 +08:00
|
|
|
|
2007-10-10 14:09:56 +08:00
|
|
|
template<typename T> typename EiNumTraits<T>::Real EiAbs2(const T& x)
|
|
|
|
{ return EiNumTraits<T>::abs2(x); }
|
2007-10-07 20:44:42 +08:00
|
|
|
|
2007-10-10 14:09:56 +08:00
|
|
|
template<typename T> T EiRand()
|
|
|
|
{ return EiNumTraits<T>::rand(); }
|
2007-10-07 20:44:42 +08:00
|
|
|
|
|
|
|
template<typename T> bool EiNegligible(const T& a, const T& b)
|
|
|
|
{
|
2007-10-10 14:09:56 +08:00
|
|
|
return(EiAbs(a) <= EiAbs(b) * EiNumTraits<T>::epsilon());
|
2007-10-07 20:44:42 +08:00
|
|
|
}
|
|
|
|
|
|
|
|
template<typename T> bool EiApprox(const T& a, const T& b)
|
|
|
|
{
|
2007-10-10 14:09:56 +08:00
|
|
|
if(EiNumTraits<T>::IsFloat)
|
|
|
|
return(EiAbs(a - b) <= std::min(EiAbs(a), EiAbs(b)) * EiNumTraits<T>::epsilon());
|
2007-10-07 20:44:42 +08:00
|
|
|
else
|
|
|
|
return(a == b);
|
|
|
|
}
|
|
|
|
|
|
|
|
template<typename T> bool EiLessThanOrApprox(const T& a, const T& b)
|
|
|
|
{
|
2007-10-10 14:09:56 +08:00
|
|
|
if(EiNumTraits<T>::IsFloat)
|
2007-10-07 20:44:42 +08:00
|
|
|
return(a < b || EiApprox(a, b));
|
|
|
|
else
|
|
|
|
return(a <= b);
|
|
|
|
}
|
|
|
|
|
2007-10-07 23:12:41 +08:00
|
|
|
#endif // EI_NUMERIC_H
|