Include float_cmp in dune-common. Fixes FS#445

[[Imported from SVN: r5418]]

common/Makefile.am

@@ -22,7 +22,8 @@ commoninclude_HEADERS = dlist.cc alignment.hh \
   collectivecommunication.hh mpihelper.hh singleton.hh \
   mpicollectivecommunication.hh geometrytype.hh utility.hh \
   bartonnackmanifcheck.hh binaryfunctions.hh lru.hh fassign.hh \
-  static_assert.hh smallobject.hh version.hh
+  static_assert.hh smallobject.hh version.hh \
+  float_cmp.cc float_cmp.hh
 
 if EXPRESSIONTEMPLATES
 commoninclude_HEADERS += exprtmpl.hh exprtmpl/scalar.inc exprtmpl/exprexpr.inc

common/float_cmp.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include "float_cmp.hh"
+
+#include <vector>
+#include <limits>
+#include <algorithm>
+#include <cstdlib>
+#include <dune/common/fvector.hh>
+
+namespace Dune {
+
+
+  namespace FloatCmp {
+    // traits
+    //! Mapping of value type to epsilon type
+    /**
+     * @ingroup FloatCmp
+     * @tparam T The value type
+     */
+    template<class T> struct EpsilonType {
+      //! The epsilon type corresponding to value type T
+      typedef T Type;
+    };
+    //! Specialization of EpsilonType for std::vector
+    /**
+     * @ingroup FloatCmp
+     * @tparam T The value_type of the std::vector
+     * @tparam A The Allocator of the std::vector
+     */
+    template<class T, typename A>
+    struct EpsilonType<std::vector<T, A> > {
+      //! The epsilon type corresponding to value type std::vector<T, A>
+      typedef EpsilonType<T> Type;
+    };
+    //! Specialization of EpsilonType for Dune::FieldVector
+    /**
+     * @ingroup FloatCmp
+     * @tparam T The field_type of the Dune::FieldVector
+     * @tparam n The size of the Dune::FieldVector
+     */
+    template<class T, int n>
+    struct EpsilonType<FieldVector<T, n> > {
+      //! The epsilon type corresponding to value type Dune::FieldVector<T, n>
+      typedef EpsilonType<T> Type;
+    };
+
+    // default epsilon
+    template<class T>
+    struct DefaultEpsilon<T, relativeWeak> {
+      static typename EpsilonType<T>::Type value()
+      { return std::numeric_limits<typename EpsilonType<T>::Type>::epsilon()*8; }
+    };
+    template<class T>
+    struct DefaultEpsilon<T, relativeStrong> {
+      static typename EpsilonType<T>::Type value()
+      { return std::numeric_limits<typename EpsilonType<T>::Type>::epsilon()*8; }
+    };
+    template<class T>
+    struct DefaultEpsilon<T, absolute> {
+      static typename EpsilonType<T>::Type value()
+      { return std::max(std::numeric_limits<typename EpsilonType<T>::Type>::epsilon(), 1e-6); }
+    };
+
+    namespace Detail {
+      // basic comparison
+      template<class T, CmpStyle style = defaultCmpStyle>
+      struct eq_t;
+      template<class T>
+      struct eq_t<T, relativeWeak> {
+        static bool eq(const T &first,
+                       const T &second,
+                       typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T>::value())
+        { return std::abs(first - second) <= epsilon*std::max(std::abs(first), std::abs(second)); }
+      };
+      template<class T>
+      struct eq_t<T, relativeStrong> {
+        static bool eq(const T &first,
+                       const T &second,
+                       typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T>::value())
+        { return std::abs(first - second) <= epsilon*std::min(std::abs(first), std::abs(second)); }
+      };
+      template<class T>
+      struct eq_t<T, absolute> {
+        static bool eq(const T &first,
+                       const T &second,
+                       typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T>::value())
+        { return std::abs(first-second) <= epsilon; }
+      };
+      template<class T, CmpStyle cstyle>
+      struct eq_t<std::vector<T>, cstyle> {
+        static bool eq(const std::vector<T> &first,
+                       const std::vector<T> &second,
+                       typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T>::value()) {
+          unsigned int size = first.size();
+          if(size != second.size()) return false;
+          for(unsigned int i = 0; i < size; ++i)
+            if(!eq_t<T, cstyle>(first[i], second[i], epsilon))
+              return false;
+          return true;
+        }
+      };
+      template<class T, int n, CmpStyle cstyle>
+      struct eq_t<Dune::FieldVector<T, n>, cstyle> {
+        static bool eq(const Dune::FieldVector<T, n> &first,
+                       const Dune::FieldVector<T, n> &second,
+                       typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T>::value()) {
+          for(int i = 0; i < n; ++i)
+            if(!eq_t<T, cstyle>(first[i], second[i], epsilon))
+              return false;
+          return true;
+        }
+      };
+    } // namespace Detail
+
+    // operations in functional style
+    template <class T, CmpStyle style>
+    bool eq(const T &first,
+            const T &second,
+            typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, style>::value())
+    {
+      return Detail::eq_t<T, style>::eq(first, second, epsilon);
+    }
+    template <class T, CmpStyle style>
+    bool ne(const T &first,
+            const T &second,
+            typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, style>::value())
+    {
+      return !eq<T, style>(first, second, epsilon);
+    }
+    template <class T, CmpStyle style>
+    bool gt(const T &first,
+            const T &second,
+            typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, style>::value())
+    {
+      return first > second && ne<T, style>(first, second, epsilon);
+    }
+    template <class T, CmpStyle style>
+    bool lt(const T &first,
+            const T &second,
+            typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, style>::value())
+    {
+      return first < second && ne<T, style>(first, second, epsilon);
+    }
+    template <class T, CmpStyle style>
+    bool ge(const T &first,
+            const T &second,
+            typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, style>::value())
+    {
+      return first > second || eq<T, style>(first, second, epsilon);
+    }
+    template <class T, CmpStyle style>
+    bool le(const T &first,
+            const T &second,
+            typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, style>::value())
+    {
+      return first < second || eq<T, style>(first, second, epsilon);
+    }
+
+    // default template arguments
+    template <class T>
+    bool eq(const T &first,
+            const T &second,
+            typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, defaultCmpStyle>::value())
+    {
+      return eq<T, defaultCmpStyle>(first, second, epsilon);
+    }
+    template <class T>
+    bool ne(const T &first,
+            const T &second,
+            typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, defaultCmpStyle>::value())
+    {
+      return ne<T, defaultCmpStyle>(first, second, epsilon);
+    }
+    template <class T>
+    bool gt(const T &first,
+            const T &second,
+            typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, defaultCmpStyle>::value())
+    {
+      return gt<T, defaultCmpStyle>(first, second, epsilon);
+    }
+    template <class T>
+    bool lt(const T &first,
+            const T &second,
+            typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, defaultCmpStyle>::value())
+    {
+      return lt<T, defaultCmpStyle>(first, second, epsilon);
+    }
+    template <class T>
+    bool ge(const T &first,
+            const T &second,
+            typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, defaultCmpStyle>::value())
+    {
+      return ge<T, defaultCmpStyle>(first, second, epsilon);
+    }
+    template <class T>
+    bool le(const T &first,
+            const T &second,
+            typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, defaultCmpStyle>::value())
+    {
+      return le<T, defaultCmpStyle>(first, second, epsilon);
+    }
+
+    // rounding operations
+    namespace Detail {
+      template<class I, class T, CmpStyle cstyle = defaultCmpStyle, RoundingStyle rstyle = defaultRoundingStyle>
+      struct round_t;
+      template<class I, class T, CmpStyle cstyle>
+      struct round_t<I, T, cstyle, downward> {
+        static I
+        round(const T &val,
+              typename EpsilonType<T>::Type epsilon = (DefaultEpsilon<T, cstyle>::value())) {
+          // first get an approximation
+          I lower = I(val);
+          I upper;
+          if(eq<T, cstyle>(T(lower), val, epsilon)) return lower;
+          if(T(lower) > val) { upper = lower; lower--; }
+          else upper = lower+1;
+          if(le<T, cstyle>(val - T(lower), T(upper) - val, epsilon))
+            return lower;
+          else return upper;
+        }
+      };
+      template<class I, class T, CmpStyle cstyle>
+      struct round_t<I, T, cstyle, upward> {
+        static I
+        round(const T &val,
+              typename EpsilonType<T>::Type epsilon = (DefaultEpsilon<T, cstyle>::value())) {
+          // first get an approximation
+          I lower = I(val);
+          I upper;
+          if(eq<T, cstyle>(T(lower), val, epsilon)) return lower;
+          if(T(lower) > val) { upper = lower; lower--; }
+          else upper = lower+1;
+          if(lt<T, cstyle>(val - T(lower), T(upper) - val, epsilon))
+            return lower;
+          else return upper;
+        }
+      };
+      template<class I, class T, CmpStyle cstyle>
+      struct round_t<I, T, cstyle, towardZero> {
+        static I
+        round(const T &val,
+              typename EpsilonType<T>::Type epsilon = (DefaultEpsilon<T, cstyle>::value())) {
+          if(val > T(0))
+            return round_t<I, T, cstyle, downward>::round(val, epsilon);
+          else return round_t<I, T, cstyle, upward>::round(val, epsilon);
+        }
+      };
+      template<class I, class T, CmpStyle cstyle>
+      struct round_t<I, T, cstyle, towardInf> {
+        static I
+        round(const T &val,
+              typename EpsilonType<T>::Type epsilon = (DefaultEpsilon<T, cstyle>::value())) {
+          if(val > T(0))
+            return round_t<I, T, cstyle, upward>::round(val, epsilon);
+          else return round_t<I, T, cstyle, downward>::round(val, epsilon);
+        }
+      };
+      template<class I, class T, CmpStyle cstyle, RoundingStyle rstyle>
+      struct round_t<std::vector<I>, std::vector<T>, cstyle, rstyle> {
+        static std::vector<I>
+        round(const T &val,
+              typename EpsilonType<T>::Type epsilon = (DefaultEpsilon<T, cstyle>::value())) {
+          unsigned int size = val.size();
+          std::vector<I> res(size);
+          for(unsigned int i = 0; i < size; ++i)
+            res[i] = round_t<I, T, cstyle, rstyle>::round(val[i], epsilon);
+          return res;
+        }
+      };
+      template<class I, class T, int n, CmpStyle cstyle, RoundingStyle rstyle>
+      struct round_t<Dune::FieldVector<I, n>, Dune::FieldVector<T, n>, cstyle, rstyle> {
+        static Dune::FieldVector<I, n>
+        round(const T &val,
+              typename EpsilonType<T>::Type epsilon = (DefaultEpsilon<T, cstyle>::value())) {
+          Dune::FieldVector<I, n> res;
+          for(int i = 0; i < n; ++i)
+            res[i] = round_t<I, T, cstyle, rstyle>::round(val[i], epsilon);
+          return res;
+        }
+      };
+    } // namespace Detail
+    template<class I, class T, CmpStyle cstyle, RoundingStyle rstyle>
+    I round(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, cstyle>::value())
+    {
+      return Detail::round_t<I, T, cstyle, rstyle>::round(val, epsilon);
+    }
+    template<class I, class T, CmpStyle cstyle>
+    I round(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, cstyle>::value())
+    {
+      return round<I, T, cstyle, defaultRoundingStyle>(val, epsilon);
+    }
+    template<class I, class T, RoundingStyle rstyle>
+    I round(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, defaultCmpStyle>::value())
+    {
+      return round<I, T, defaultCmpStyle, rstyle>(val, epsilon);
+    }
+    template<class I, class T>
+    I round(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, defaultCmpStyle>::value())
+    {
+      return round<I, T, defaultCmpStyle>(val, epsilon);
+    }
+
+    // truncation
+    namespace Detail {
+      template<class I, class T, CmpStyle cstyle = defaultCmpStyle, RoundingStyle rstyle = defaultRoundingStyle>
+      struct trunc_t;
+      template<class I, class T, CmpStyle cstyle>
+      struct trunc_t<I, T, cstyle, downward> {
+        static I
+        trunc(const T &val,
+              typename EpsilonType<T>::Type epsilon = (DefaultEpsilon<T, cstyle>::value())) {
+          // this sould be optimized away unless needed
+          if(!std::numeric_limits<I>::is_signed)
+            // make sure this works for all useful cases even if I is an unsigned type
+            if(eq<T, cstyle>(val, T(0), epsilon)) return I(0);
+          // first get an approximation
+          I lower = I(val); // now |val-lower| < 1
+          // make sure we're really lower in case the cast truncated to an unexpected direction
+          if(T(lower) > val) lower--; // now val-lower < 1
+          // check whether lower + 1 is approximately val
+          if(eq<T, cstyle>(T(lower+1), val, epsilon))
+            return lower+1;
+          else return lower;
+        }
+      };
+      template<class I, class T, CmpStyle cstyle>
+      struct trunc_t<I, T, cstyle, upward> {
+        static I
+        trunc(const T &val,
+              typename EpsilonType<T>::Type epsilon = (DefaultEpsilon<T, cstyle>::value())) {
+          I upper = trunc_t<I, T, cstyle, downward>::trunc(val, epsilon);
+          if(ne<T, cstyle>(T(upper), val, epsilon)) ++upper;
+          return upper;
+        }
+      };
+      template<class I, class T, CmpStyle cstyle>
+      struct trunc_t<I, T, cstyle, towardZero> {
+        static I
+        trunc(const T &val,
+              typename EpsilonType<T>::Type epsilon = (DefaultEpsilon<T, cstyle>::value())) {
+          if(val > T(0)) return trunc_t<I, T, cstyle, downward>::trunc(val, epsilon);
+          else return trunc_t<I, T, cstyle, upward>::trunc(val, epsilon);
+        }
+      };
+      template<class I, class T, CmpStyle cstyle>
+      struct trunc_t<I, T, cstyle, towardInf> {
+        static I
+        trunc(const T &val,
+              typename EpsilonType<T>::Type epsilon = (DefaultEpsilon<T, cstyle>::value())) {
+          if(val > T(0)) return trunc_t<I, T, cstyle, upward>::trunc(val, epsilon);
+          else return trunc_t<I, T, cstyle, downward>::trunc(val, epsilon);
+        }
+      };
+      template<class I, class T, CmpStyle cstyle, RoundingStyle rstyle>
+      struct trunc_t<std::vector<I>, std::vector<T>, cstyle, rstyle> {
+        static std::vector<I>
+        trunc(const std::vector<T> &val,
+              typename EpsilonType<T>::Type epsilon = (DefaultEpsilon<T, cstyle>::value())) {
+          unsigned int size = val.size();
+          std::vector<I> res(size);
+          for(unsigned int i = 0; i < size; ++i)
+            res[i] = trunc_t<I, T, cstyle, rstyle>::trunc(val[i], epsilon);
+          return res;
+        }
+      };
+      template<class I, class T, int n, CmpStyle cstyle, RoundingStyle rstyle>
+      struct trunc_t<Dune::FieldVector<I, n>, Dune::FieldVector<T, n>, cstyle, rstyle> {
+        static Dune::FieldVector<I, n>
+        trunc(const Dune::FieldVector<T, n> &val,
+              typename EpsilonType<T>::Type epsilon = (DefaultEpsilon<T, cstyle>::value())) {
+          Dune::FieldVector<I, n> res;
+          for(int i = 0; i < n; ++i)
+            res[i] = trunc_t<I, T, cstyle, rstyle>::trunc(val[i], epsilon);
+          return res;
+        }
+      };
+    } // namespace Detail
+    template<class I, class T, CmpStyle cstyle, RoundingStyle rstyle>
+    I trunc(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, cstyle>::value())
+    {
+      return Detail::trunc_t<I, T, cstyle, rstyle>::trunc(val, epsilon);
+    }
+    template<class I, class T, CmpStyle cstyle>
+    I trunc(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, cstyle>::value())
+    {
+      return trunc<I, T, cstyle, defaultRoundingStyle>(val, epsilon);
+    }
+    template<class I, class T, RoundingStyle rstyle>
+    I trunc(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, defaultCmpStyle>::value())
+    {
+      return trunc<I, T, defaultCmpStyle, rstyle>(val, epsilon);
+    }
+    template<class I, class T>
+    I trunc(const T &val, typename EpsilonType<T>::Type epsilon = DefaultEpsilon<T, defaultCmpStyle>::value())
+    {
+      return trunc<I, T, defaultCmpStyle>(val, epsilon);
+    }
+  } //namespace Dune
+
+  // oo interface
+  template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
+  FloatCmpOps<T, cstyle_, rstyle_>::
+  FloatCmpOps(EpsilonType epsilon) : epsilon_(epsilon) {}
+
+
+  template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
+  typename FloatCmpOps<T, cstyle_, rstyle_>::EpsilonType
+  FloatCmpOps<T, cstyle_, rstyle_>::epsilon() const
+  {
+    return epsilon_;
+  }
+
+  template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
+  void
+  FloatCmpOps<T, cstyle_, rstyle_>::epsilon(EpsilonType epsilon__)
+  {
+    epsilon_ = epsilon__;
+  }
+
+
+  template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
+  bool FloatCmpOps<T, cstyle_, rstyle_>::
+  eq(const ValueType &first, const ValueType &second) const
+  {
+    return Dune::FloatCmp::eq<ValueType, cstyle>(first, second, epsilon_);
+  }
+
+  template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
+  bool FloatCmpOps<T, cstyle_, rstyle_>::
+  ne(const ValueType &first, const ValueType &second) const
+  {
+    return Dune::FloatCmp::ne<ValueType, cstyle>(first, second, epsilon_);
+  }
+
+  template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
+  bool FloatCmpOps<T, cstyle_, rstyle_>::
+  gt(const ValueType &first, const ValueType &second) const
+  {
+    return Dune::FloatCmp::gt<ValueType, cstyle>(first, second, epsilon_);
+  }
+
+  template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
+  bool FloatCmpOps<T, cstyle_, rstyle_>::
+  lt(const ValueType &first, const ValueType &second) const
+  {
+    return Dune::FloatCmp::lt<ValueType, cstyle>(first, second, epsilon_);
+  }
+
+  template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
+  bool FloatCmpOps<T, cstyle_, rstyle_>::
+  ge(const ValueType &first, const ValueType &second) const
+  {
+    return Dune::FloatCmp::ge<ValueType, cstyle>(first, second, epsilon_);
+  }
+
+  template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
+  bool FloatCmpOps<T, cstyle_, rstyle_>::
+  le(const ValueType &first, const ValueType &second) const
+  {
+    return Dune::FloatCmp::le<ValueType, cstyle>(first, second, epsilon_);
+  }
+
+
+  template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
+  template<class I>
+  I FloatCmpOps<T, cstyle_, rstyle_>::
+  round(const ValueType &val) const
+  {
+    return Dune::FloatCmp::round<I, ValueType, cstyle, rstyle_>(val, epsilon_);
+  }
+
+  template<class T, FloatCmp::CmpStyle cstyle_, FloatCmp::RoundingStyle rstyle_>
+  template<class I>
+  I FloatCmpOps<T, cstyle_, rstyle_>::
+  trunc(const ValueType &val) const
+  {
+    return Dune::FloatCmp::trunc<I, ValueType, cstyle, rstyle_>(val, epsilon_);
+  }
+
+} //namespace Dune

common/float_cmp.hh

+// -*- 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
+
+/**
+   @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.  The solution is to compare approximately, using an epsilon
+   with says how much deviation to accept.
+
+   The most straightforward way of comparing is using an @em absolute epsilon.
+   The comparison is done like
+   @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.  Nethertheless this method of comparing is not
+   recommended in general, since 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 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 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 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 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 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 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

common/test/Makefile.am

+# -*- tab-width: 4; indent-tabs-mode: nil -*-
 # $Id$
 
 TESTPROGS = parsetest test-stack arraylisttest smartpointertest \
@@ -10,7 +11,8 @@ TESTPROGS = parsetest test-stack arraylisttest smartpointertest \
 	testfassign4 \
     testfassign_fail1 testfassign_fail2 testfassign_fail3 \
     testfassign_fail4 testfassign_fail5 testfassign_fail6 \
-    conversiontest bitsetvectortest deprecatedtuplestest
+    conversiontest bitsetvectortest deprecatedtuplestest \
+    float_cmp
 
 # which tests to run
 COMPILE_XFAIL=$(DUNE_COMMON_ROOT)/bin/xfail-compile-tests
@@ -131,4 +133,6 @@ conversiontest_SOURCES = conversiontest.cc
 
 sourcescheck_NOSOURCES = exprtmpl.cc timing.cc
 
+float_cmp_SOURCES = float_cmp.cc
+
 include $(top_srcdir)/am/global-rules

common/test/float_cmp.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+// Test the new (Dune) interface of float_cmp
+
+#include <iostream>
+
+#include <dune/common/float_cmp.hh>
+
+using std::cout;
+using std::endl;
+using std::flush;
+
+/////////////////////////
+//
+//  compile time checks
+//
+
+// check that we can access the functions as FloatCmp::function from within the Dune namespace
+namespace Dune {
+  void checkNamespaceAccess() {
+    FloatCmp::eq(0.0, 0.0);
+  }
+} // namespace Dune
+  // check that we can access the functions as FloatCmp::function with using namespace Dune
+void checkUsingAccess() {
+  using namespace Dune;
+  FloatCmp::eq(0.0, 0.0);
+}
+
+// run time checks
+const char* repr(bool b) {
+  if(b) return "true ";
+  else return "false";
+}
+
+int passed = 0;
+int failed = 0;
+
+void count(bool pass) {
+  if(pass) { cout << "passed"; ++passed; }
+  else     { cout << "failed"; ++failed; }
+}
+
+template<typename F>
+void tests(F f1, F f2, typename Dune::FloatCmp::EpsilonType<F>::Type eps, bool inside)
+{
+  bool result;
+
+  cout << "eq(" << f1 << ", " << f2 << ", " << eps << ") = " << flush;
+  result = Dune::FloatCmp::eq(f1, f2, eps);
+  cout << repr(result) << "\t";
+  count(result == inside);
+  cout << endl;
+
+  cout << "ge(" << f1 << ", " << f2 << ", " << eps << ") = " << flush;
+  result = Dune::FloatCmp::ge(f1, f2, eps);
+  cout << repr(result) << "\t";
+  count(result == (inside || f1 > f2));
+  cout << endl;
+
+  cout << "le(" << f1 << ", " << f2 << ", " << eps << ") = " << flush;
+  result = Dune::FloatCmp::le(f1, f2, eps);
+  cout << repr(result) << "\t";
+  count(result == (inside || f1 < f2));
+  cout << endl;
+
+  cout << "ne(" << f1 << ", " << f2 << ", " << eps << ") = " << flush;
+  result = Dune::FloatCmp::ne(f1, f2, eps);
+  cout << repr(result) << "\t";
+  count(result == !inside);
+  cout << endl;
+
+  cout << "gt(" << f1 << ", " << f2 << ", " << eps << ") = " << flush;
+  result = Dune::FloatCmp::gt(f1, f2, eps);
+  cout << repr(result) << "\t";
+  count(result == (!inside && f1 > f2));
+  cout << endl;
+
+  cout << "lt(" << f1 << ", " << f2 << ", " << eps << ") = " << flush;
+  result = Dune::FloatCmp::lt(f1, f2, eps);
+  cout << repr(result) << "\t";
+  count(result == (!inside && f1 < f2));
+  cout << endl;
+}
+
+template<typename F>
+void tests(F f1, F f2, const typename Dune::FloatCmpOps<F> &ops, bool inside)
+{
+  bool result;
+  cout << "ops = operations(" << ops.epsilon() << ")" << endl;
+
+  cout << "ops.eq(" << f1 << ", " << f2 << ") = " << flush;
+  result = ops.eq(f1, f2);
+  cout << repr(result) << "\t";
+  count(result == inside);
+  cout << endl;
+
+  cout << "ops.ge(" << f1 << ", " << f2 << ") = " << flush;
+  result = ops.ge(f1, f2);
+  cout << repr(result) << "\t";
+  count(result == (inside || f1 > f2));
+  cout << endl;
+
+  cout << "ops.le(" << f1 << ", " << f2 << ") = " << flush;
+  result = ops.le(f1, f2);
+  cout << repr(result) << "\t";
+  count(result == (inside || f1 < f2));
+  cout << endl;
+
+  cout << "ops.ne(" << f1 << ", " << f2 << ") = " << flush;
+  result = ops.ne(f1, f2);
+  cout << repr(result) << "\t";
+  count(result == !inside);
+  cout << endl;
+
+  cout << "ops.gt(" << f1 << ", " << f2 << ") = " << flush;
+  result = ops.gt(f1, f2);
+  cout << repr(result) << "\t";
+  count(result == (!inside && f1 > f2));
+  cout << endl;
+
+  cout << "ops.lt(" << f1 << ", " << f2 << ") = " << flush;
+  result = ops.lt(f1, f2);
+  cout << repr(result) << "\t";
+  count(result == (!inside && f1 < f2));
+  cout << endl;
+}
+
+int main() {
+  cout.setf(std::ios_base::scientific, std::ios_base::floatfield);
+  cout.precision(16);
+  Dune::FloatCmpOps<double> ops(1e-7);
+
+  cout << "Tests inside the epsilon environment" << endl;
+  tests<double>(1.0, 1.00000001, 1e-7, true);
+  tests<double>(1.0, 1.00000001, ops,  true);
+
+  cout << "Tests outside the epsilon environment, f1 < f2" << endl;
+  tests<double>(1.0, 1.000001, 1e-7, false);
+  tests<double>(1.0, 1.000001, ops,  false);
+
+  cout << "Tests outside the epsilon environment, f1 > f2" << endl;
+  tests<double>(1.000001, 1.0, 1e-7, false);
+  tests<double>(1.000001, 1.0, ops,  false);
+
+  cout << "Tests with f1 = f2 = 0" << endl;
+  tests<double>(0, 0, 1e-7, true);
+  tests<double>(0, 0, ops,  true);
+
+  int total = passed + failed;
+  cout << passed << "/" << total << " tests passed; " << failed << "/" << total << " tests failed" << endl;
+  if(failed > 0) return 1;
+  else return 0;
+}