Forked from
Core Modules / dune-common
4175 commits behind, 90 commits ahead of the upstream repository.
-
Nicolas Schwenck authored
floating comparisons (cherry picked from commit 94a327ac) Signed-off-by:
Steffen Müthing <muething@dune-project.org>
Nicolas Schwenck authoredfloating comparisons (cherry picked from commit 94a327ac) Signed-off-by:
float_cmp.hh 14.85 KiB
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_COMMON_FLOAT_CMP_HH
#define DUNE_COMMON_FLOAT_CMP_HH
/** \file
* \brief Various ways to compare floating-point numbers
*/
/**
@addtogroup FloatCmp FloatCmp
@ingroup Common
@section How_to_compare How to compare floats
When comparing floating point numbers for equality, one often faces the
problem that floating point operations are not always exact. For example on
i386 the expression
@code
0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 + 0.2 == 2.0
@endcode
evaluates to
@code
1.99999999999999977796 == 2.00000000000000000000
@endcode
which is false. One solution is to compare approximately, using an epsilon
which says how much deviation to accept.
The most straightforward way of comparing is using an @em absolute epsilon.
This means comparison for equality is replaced by
@code
abs(first-second) <= epsilon
@endcode
This has a severe disadvantage: if you have an epsilon like 1e-10 but first
and second are of the magnitude 1e-15 everything will compare equal which is
certainly not what you want. This can be overcome by selecting an
appropriate epsilon. Nevertheless this method of comparing is not
recommended in general, and we will present a more robus method in the
next paragraph.
There is another way of comparing approximately, using a @em relative
epsilon which is then scaled with first:
@code
abs(first-second) <= epsilon * abs(first)
@endcode
Of cource the comparison should be symmetric in first and second so we
cannot arbitrarily select either first or second to scale epsilon. The are
two symmetric variants, @em relative_weak
@code
abs(first-second) <= epsilon * max(abs(first), abs(second))
@endcode
and @em relative_strong
@code
abs(first-second) <= epsilon * min(abs(first), abs(second))
@endcode
Both variants are good, but in practice the relative_weak variant is
preferred. This is also the default variant.
\note Although using a relative epsilon is better than using an absolute
epsilon, using a relative epsilon leads to problems if either first or
second equals 0. In principle the relative method can be combined
with an absolute method using an epsilon near the minimum
representable positive value, but this is not implemented here.
There is a completely different way of comparing floats. Instead of giving
an epsilon, the programmer states how many representable value are allowed
between first and second. See the "Comparing using integers" section in
http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
for more about that.
@section Interface Interface
To do the comparison, you can use the free functions @link
Dune::FloatCmp::eq eq()@endlink, @link Dune::FloatCmp::ne ne()@endlink,
@link Dune::FloatCmp::gt gt()@endlink, @link Dune::FloatCmp::lt
lt()@endlink, @link Dune::FloatCmp::ge ge()@endlink and @link
Dune::FloatCmp::le le()@endlink from the namespace Dune::FloatCmp. They
take the values to compare and optionally an epsilon, which defaults to 8
times the machine epsilon (the difference between 1.0 and the smallest
representable value > 1.0) for relative comparisons, or simply 1e-6 for
absolute comparisons. The compare style can be given as an optional second
template parameter and defaults to relative_weak.
You can also use the class Dune::FloatCmpOps which has @link
Dune::FloatCmpOps::eq eq()@endlink, @link Dune::FloatCmpOps::ne
ne()@endlink, @link Dune::FloatCmpOps::gt gt()@endlink, @link
Dune::FloatCmpOps::lt lt()@endlink, @link Dune::FloatCmpOps::ge ge()@endlink
and @link Dune::FloatCmpOps::le le()@endlink as member functions. In this
case the class encapsulates the epsilon and the comparison style (again the
defaults from the previous paragraph apply). This may be more convenient if
you write your own class utilizing floating point comparisons, and you want
the user of you class to specify epsilon and compare style.
*/
//! Dune namespace
namespace Dune {
//! FloatCmp namespace
//! @ingroup FloatCmp
namespace FloatCmp {
// basic constants
//! How to compare
//! @ingroup FloatCmp
enum CmpStyle {
//! |a-b|/|a| <= epsilon || |a-b|/|b| <= epsilon
relativeWeak,
//! |a-b|/|a| <= epsilon && |a-b|/|b| <= epsilon
relativeStrong,
//! |a-b| <= epsilon
absolute,
//! the global default compare style (relative_weak)
defaultCmpStyle = relativeWeak
};
//! How to round or truncate
//! @ingroup FloatCmp
enum RoundingStyle {
//! always round toward 0
towardZero,
//! always round away from 0
towardInf,
//! round toward \f$-\infty\f$
downward,
//! round toward \f$+\infty\f$
upward,
//! the global default rounding style (toward_zero)
defaultRoundingStyle = towardZero
};
template<class T> struct EpsilonType;
//! mapping from a value type and a compare style to a default epsilon
/**
* @ingroup FloatCmp
* @tparam T The value type to map from
* @tparam style The compare style to map from
*/
template<class T, CmpStyle style = defaultCmpStyle>
struct DefaultEpsilon {
//! Returns the default epsilon for the given value type and compare style
static typename EpsilonType<T>::Type value();
};
// operations in functional style
//! @addtogroup FloatCmp
//! @{
//! test for equality using epsilon
/**
* @tparam T Type of the values to compare
* @tparam style How to compare. This defaults to defaultCmpStyle.
* @param first left operand of equals operation
* @param second right operand of equals operation
* @param epsilon The epsilon to use in the comparison
*/
template <class T, CmpStyle style /*= defaultCmpStyle*/>
bool eq(const T &first,
const T &second,
typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, style>::value());
//! test for inequality using epsilon
/**
* @tparam T Type of the values to compare
* @tparam style How to compare. This defaults to defaultCmpStyle.
* @param first left operand of not-equal operation
* @param second right operand of not-equal operation
* @param epsilon The epsilon to use in the comparison
* @return !eq(first, second, epsilon)
*/
template <class T, CmpStyle style /*= defaultCmpStyle*/>
bool ne(const T &first,
const T &second,
typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, style>::value());
//! test if first greater than second
/**
* @tparam T Type of the values to compare
* @tparam style How to compare. This defaults to defaultCmpStyle.
* @param first left operand of greater-than operation
* @param second right operand of greater-than operation
* @param epsilon The epsilon to use in the comparison
* @return ne(first, second, epsilon) && first > second
*
* this is like first > second but the region that compares equal with an
* epsilon is excluded
*/
template <class T, CmpStyle style /*= defaultCmpStyle*/>
bool gt(const T &first,
const T &second,
typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, style>::value());
//! test if first lesser than second
/**
* @tparam T Type of the values to compare
* @tparam style How to compare. This defaults to defaultCmpStyle.
* @param first left operand of less-than operation
* @param second right operand of less-than operation
* @param epsilon The epsilon to use in the comparison
* @return ne(first, second, epsilon) && first < second
*
* this is like first < second, but the region that compares equal with an
* epsilon is excluded
*/
template <class T, CmpStyle style /*= defaultCmpStyle*/>
bool lt(const T &first,
const T &second,
typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, style>::value());
//! test if first greater or equal second
/**
* @tparam T Type of the values to compare
* @tparam style How to compare. This defaults to defaultCmpStyle.
* @param first left operand of greater-or-equals operation
* @param second right operand of greater-or-equals operation
* @param epsilon The epsilon to use in the comparison * @return eq(first, second, epsilon) || first > second
*
* this is like first > second, but the region that compares equal with an
* epsilon is also included
*/
template <class T, CmpStyle style /*= defaultCmpStyle*/>
bool ge(const T &first,
const T &second,
typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, style>::value());
//! test if first lesser or equal second
/**
* @tparam T Type of the values to compare
* @tparam style How to compare. This defaults to defaultCmpStyle.
* @param first left operand of less-or-equals operation
* @param second right operand of less-or-equals operation
* @param epsilon The epsilon to use in the comparison
* @return eq(first, second) || first < second
*
* this is like first < second, but the region that compares equal with an
* epsilon is also included
*/
template <class T, CmpStyle style /*= defaultCmpStyle*/>
bool le(const T &first,
const T &second,
typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, style>::value());
// rounding operations
//! round using epsilon
/**
* @tparam I The integral type to round to
* @tparam T Type of the value to round
* @tparam cstyle How to compare. This defaults to defaultCmpStyle.
* @tparam rstyle How to round. This defaults to defaultRoundingStyle.
* @param val The value to round
* @param epsilon The epsilon to use in comparisons
* @return The rounded value
*
* Round according to rstyle. If val is already near the mean of two
* adjacent integers in terms of epsilon, the result will be the rounded
* mean.
*/
template<class I, class T, CmpStyle cstyle /*= defaultCmpStyle*/, RoundingStyle rstyle /*= defaultRoundingStyle*/>
I round(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, cstyle>::value());
// truncation
//! truncate using epsilon
/**
* @tparam I The integral type to truncate to
* @tparam T Type of the value to truncate
* @tparam cstyle How to compare. This defaults to defaultCmpStyle.
* @tparam rstyle How to truncate. This defaults to defaultRoundingStyle.
* @param val The value to truncate
* @param epsilon The epsilon to use in comparisons
* @return The truncated value
*
* Truncate according to rstyle. If val is already near an integer in
* terms of epsilon, the result will be that integer instead of the real
* truncated value.
*/
template<class I, class T, CmpStyle cstyle /*= defaultCmpStyle*/, RoundingStyle rstyle /*= defaultRoundingStyle*/>
I trunc(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, cstyle>::value());
//! @}
// group FloatCmp
} //namespace FloatCmp
// oo interface
//! Class encapsulating a default epsilon
/**
* @ingroup FloatCmp * @tparam T Type of the values to compare
* @tparam cstyle_ How to compare
* @tparam rstyle_ How to round
*/
template<class T, FloatCmp::CmpStyle cstyle_ = FloatCmp::defaultCmpStyle,
FloatCmp::RoundingStyle rstyle_ = FloatCmp::defaultRoundingStyle>
class FloatCmpOps {
typedef FloatCmp::CmpStyle CmpStyle;
typedef FloatCmp::RoundingStyle RoundingStyle;
public:
// record template parameters
//! How comparisons are done
static const CmpStyle cstyle = cstyle_;
//! How rounding is done
static const RoundingStyle rstyle = rstyle_;
//! Type of the values to compare
typedef T ValueType;
//! Type of the epsilon.
/**
* May be different from the value type, for example for complex<double>
*/
typedef typename FloatCmp::EpsilonType<T>::Type EpsilonType;
private:
EpsilonType epsilon_;
typedef FloatCmp::DefaultEpsilon<EpsilonType, cstyle> DefaultEpsilon;
public:
//! construct an operations object
/**
* @param epsilon Use the specified epsilon for comparing
*/
FloatCmpOps(EpsilonType epsilon = DefaultEpsilon::value());
//! return the current epsilon
EpsilonType epsilon() const;
//! set new epsilon
void epsilon(EpsilonType epsilon__);
//! test for equality using epsilon
bool eq(const ValueType &first, const ValueType &second) const;
//! test for inequality using epsilon
/**
* this is exactly !eq(first, second)
*/
bool ne(const ValueType &first, const ValueType &second) const;
//! test if first greater than second
/**
* this is exactly ne(first, second) && first > second, i.e. greater but
* the region that compares equal with an epsilon is excluded
*/
bool gt(const ValueType &first, const ValueType &second) const;
//! test if first lesser than second
/**
* this is exactly ne(first, second) && first < second, i.e. lesser but
* the region that compares equal with an epsilon is excluded
*/
bool lt(const ValueType &first, const ValueType &second) const;
//! test if first greater or equal second
/**
* this is exactly eq(first, second) || first > second, i.e. greater but
* the region that compares equal with an epsilon is also included
*/
bool ge(const ValueType &first, const ValueType &second) const;
//! test if first lesser or equal second
/**
* this is exactly eq(first, second) || first > second, i.e. lesser but
* the region that compares equal with an epsilon is also included */
bool le(const ValueType &first, const ValueType &second) const;
//! round using epsilon
/**
* @tparam I The integral type to round to
*
* @param val The value to round
*
* Round according to rstyle. If val is already near the mean of two
* adjacent integers in terms of epsilon, the result will be the rounded
* mean.
*/
template<class I>
I round(const ValueType &val) const;
//! truncate using epsilon
/**
* @tparam I The integral type to truncate to
*
* @param val The value to truncate
*
* Truncate according to rstyle. If val is already near an integer in
* terms of epsilon, the result will be that integer instead of the real
* truncated value.
*/
template<class I>
I trunc(const ValueType &val) const;
};
} //namespace Dune
#include "float_cmp.cc"
#endif //DUNE_COMMON_FLOAT_CMP_HH