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)
    {
      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)
    {
      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)
    {
      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)
    {
      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)
    {
      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)
    {
      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