diff --git a/cmake/modules/AddQuadMathFlags.cmake b/cmake/modules/AddQuadMathFlags.cmake index cae9bbfc9d6d75941c339cd7988f12fbecfbf1d5..85fc3181bbc02628ab10e8940487a1c4381f9766 100644 --- a/cmake/modules/AddQuadMathFlags.cmake +++ b/cmake/modules/AddQuadMathFlags.cmake @@ -21,7 +21,7 @@ function(add_dune_quadmath_flags _targets) if(${CMAKE_CXX_COMPILER_ID} STREQUAL GNU) set_property(TARGET ${_target} APPEND_STRING - PROPERTY COMPILE_FLAGS "-fext-numeric-literals ") + PROPERTY COMPILE_FLAGS "-fext-numeric-literals -Wno-pedantic ") endif() endforeach(_target ${_targets}) endif(QUADMATH_FOUND) diff --git a/cmake/modules/DuneCommonExtensionsMacros.cmake b/cmake/modules/DuneCommonExtensionsMacros.cmake index 613dfb664b75999b008f0003a4d7c409cbf409ad..f35d463d1205c5940b681e9532c60a3ee96a93a2 100644 --- a/cmake/modules/DuneCommonExtensionsMacros.cmake +++ b/cmake/modules/DuneCommonExtensionsMacros.cmake @@ -1 +1,4 @@ # File for module specific CMake tests. + +include(AddQuadMathFlags) +find_package(QuadMath) \ No newline at end of file diff --git a/cmake/modules/FindQuadMath.cmake b/cmake/modules/FindQuadMath.cmake index ff70b337941ff611571df0dc4af78c391a05686a..08e455aa674922cf93635f0bb003d1d7c0658adb 100644 --- a/cmake/modules/FindQuadMath.cmake +++ b/cmake/modules/FindQuadMath.cmake @@ -48,7 +48,7 @@ set(HAVE_QUADMATH ${QUADMATH_FOUND}) # -fext-numeric-literals is a GCC extension not available in other compilers like clang if(${CMAKE_CXX_COMPILER_ID} STREQUAL GNU) - set(_QUADMATH_EXT_NUMERIC_LITERALS "-fext-numeric-literals") + set(_QUADMATH_EXT_NUMERIC_LITERALS "-fext-numeric-literals" "-Wno-pedantic") endif() # register all QuadMath related flags diff --git a/config.h.cmake b/config.h.cmake index 8d7e6c0215f150942c02299ab2210ec6797c15c6..06c6ea79069bc25d1b5d4ef962ea1f8bf805dfc5 100644 --- a/config.h.cmake +++ b/config.h.cmake @@ -4,6 +4,10 @@ overwritten */ +/* Define if you have the GCC Quad-Precision library. The value should be ENABLE_QUADMATH + to facilitate activating and deactivating QuadMath using compile flags. */ +#cmakedefine HAVE_QUADMATH ENABLE_QUADMATH + /* begin private */ /* Name of package */ #define PACKAGE "@DUNE_MOD_NAME@" diff --git a/dune/common/concurrentcache.hh b/dune/common/concurrentcache.hh index d8095fd623ab62218db11868bf563f6cbdd96c49..7edddc84948430cde5d8030c7eb0966d701bb875 100644 --- a/dune/common/concurrentcache.hh +++ b/dune/common/concurrentcache.hh @@ -116,6 +116,7 @@ namespace Dune mutable container_type cachedData_; }; + // implementation of the ThreadLocal policy. Data is stored in thread_local variable. template struct ThreadLocalPolicy diff --git a/dune/common/float_cmp.cc b/dune/common/float_cmp.cc new file mode 100644 index 0000000000000000000000000000000000000000..1dc18a89b6103b48c5d7bd97b00ac89220297ad5 --- /dev/null +++ b/dune/common/float_cmp.cc @@ -0,0 +1,506 @@ +// -*- 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 +#include +#include +#include +#include + +namespace Dune { + + + namespace FloatCmp { + // traits + //! Mapping of value type to epsilon type + /** + * @ingroup FloatCmp + * @tparam T The value type + */ + template 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 + struct EpsilonType > { + //! The epsilon type corresponding to value type std::vector + typedef typename EpsilonType::Type 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 + struct EpsilonType > { + //! The epsilon type corresponding to value type Dune::FieldVector + typedef typename EpsilonType::Type Type; + }; + + // default epsilon + template + struct DefaultEpsilon { + static typename EpsilonType::Type value() + { return std::numeric_limits::Type>::epsilon()*8.; } + }; + template + struct DefaultEpsilon { + static typename EpsilonType::Type value() + { return std::numeric_limits::Type>::epsilon()*8.; } + }; + template + struct DefaultEpsilon { + static typename EpsilonType::Type value() + { return std::max(std::numeric_limits::Type>::epsilon(), 1e-6); } + }; + + namespace Impl { + // basic comparison + template + struct eq_t; + template + struct eq_t { + static bool eq(const T &first, + const T &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::value()) + { + using std::abs; + return abs(first - second) <= epsilon*std::max(abs(first), abs(second)); + } + }; + template + struct eq_t { + static bool eq(const T &first, + const T &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::value()) + { + using std::abs; + return abs(first - second) <= epsilon*std::min(abs(first), abs(second)); + } + }; + template + struct eq_t { + static bool eq(const T &first, + const T &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::value()) + { + using std::abs; + return abs(first-second) <= epsilon; + } + }; + template + struct eq_t_std_vec { + typedef std::vector V; + static bool eq(const V &first, + const V &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::value()) { + auto size = first.size(); + if(size != second.size()) return false; + for(unsigned int i = 0; i < size; ++i) + if(!eq_t::eq(first[i], second[i], epsilon)) + return false; + return true; + } + }; + template< class T> + struct eq_t, relativeWeak> : eq_t_std_vec {}; + template< class T> + struct eq_t, relativeStrong> : eq_t_std_vec {}; + template< class T> + struct eq_t, absolute> : eq_t_std_vec {}; + + template + struct eq_t_fvec { + typedef Dune::FieldVector V; + static bool eq(const V &first, + const V &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::value()) { + for(int i = 0; i < n; ++i) + if(!eq_t::eq(first[i], second[i], epsilon)) + return false; + return true; + } + }; + template< class T, int n > + struct eq_t< Dune::FieldVector, relativeWeak> : eq_t_fvec {}; + template< class T, int n > + struct eq_t< Dune::FieldVector, relativeStrong> : eq_t_fvec {}; + template< class T, int n > + struct eq_t< Dune::FieldVector, absolute> : eq_t_fvec {}; + } // namespace Impl + + // operations in functional style + template + bool eq(const T &first, + const T &second, + typename EpsilonType::Type epsilon) + { + return Impl::eq_t::eq(first, second, epsilon); + } + template + bool ne(const T &first, + const T &second, + typename EpsilonType::Type epsilon) + { + return !eq(first, second, epsilon); + } + template + bool gt(const T &first, + const T &second, + typename EpsilonType::Type epsilon) + { + return first > second && ne(first, second, epsilon); + } + template + bool lt(const T &first, + const T &second, + typename EpsilonType::Type epsilon) + { + return first < second && ne(first, second, epsilon); + } + template + bool ge(const T &first, + const T &second, + typename EpsilonType::Type epsilon) + { + return first > second || eq(first, second, epsilon); + } + template + bool le(const T &first, + const T &second, + typename EpsilonType::Type epsilon) + { + return first < second || eq(first, second, epsilon); + } + + // default template arguments + template + bool eq(const T &first, + const T &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::value()) + { + return eq(first, second, epsilon); + } + template + bool ne(const T &first, + const T &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::value()) + { + return ne(first, second, epsilon); + } + template + bool gt(const T &first, + const T &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::value()) + { + return gt(first, second, epsilon); + } + template + bool lt(const T &first, + const T &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::value()) + { + return lt(first, second, epsilon); + } + template + bool ge(const T &first, + const T &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::value()) + { + return ge(first, second, epsilon); + } + template + bool le(const T &first, + const T &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::value()) + { + return le(first, second, epsilon); + } + + // rounding operations + namespace Impl { + template + struct round_t; + template + struct round_t { + static I + round(const T &val, + typename EpsilonType::Type epsilon = (DefaultEpsilon::value())) { + // first get an approximation + I lower = I(val); + I upper; + if(eq(T(lower), val, epsilon)) return lower; + if(T(lower) > val) { upper = lower; lower--; } + else upper = lower+1; + if(le(val - T(lower), T(upper) - val, epsilon)) + return lower; + else return upper; + } + }; + template + struct round_t { + static I + round(const T &val, + typename EpsilonType::Type epsilon = (DefaultEpsilon::value())) { + // first get an approximation + I lower = I(val); + I upper; + if(eq(T(lower), val, epsilon)) return lower; + if(T(lower) > val) { upper = lower; lower--; } + else upper = lower+1; + if(lt(val - T(lower), T(upper) - val, epsilon)) + return lower; + else return upper; + } + }; + template + struct round_t { + static I + round(const T &val, + typename EpsilonType::Type epsilon = (DefaultEpsilon::value())) { + if(val > T(0)) + return round_t::round(val, epsilon); + else return round_t::round(val, epsilon); + } + }; + template + struct round_t { + static I + round(const T &val, + typename EpsilonType::Type epsilon = (DefaultEpsilon::value())) { + if(val > T(0)) + return round_t::round(val, epsilon); + else return round_t::round(val, epsilon); + } + }; + template + struct round_t, std::vector, cstyle, rstyle> { + static std::vector + round(const T &val, + typename EpsilonType::Type epsilon = (DefaultEpsilon::value())) { + unsigned int size = val.size(); + std::vector res(size); + for(unsigned int i = 0; i < size; ++i) + res[i] = round_t::round(val[i], epsilon); + return res; + } + }; + template + struct round_t, Dune::FieldVector, cstyle, rstyle> { + static Dune::FieldVector + round(const T &val, + typename EpsilonType::Type epsilon = (DefaultEpsilon::value())) { + Dune::FieldVector res; + for(int i = 0; i < n; ++i) + res[i] = round_t::round(val[i], epsilon); + return res; + } + }; + } // end namespace Impl + template + I round(const T &val, typename EpsilonType::Type epsilon /*= DefaultEpsilon::value()*/) + { + return Impl::round_t::round(val, epsilon); + } + template + I round(const T &val, typename EpsilonType::Type epsilon = DefaultEpsilon::value()) + { + return round(val, epsilon); + } + template + I round(const T &val, typename EpsilonType::Type epsilon = DefaultEpsilon::value()) + { + return round(val, epsilon); + } + template + I round(const T &val, typename EpsilonType::Type epsilon = DefaultEpsilon::value()) + { + return round(val, epsilon); + } + + // truncation + namespace Impl { + template + struct trunc_t; + template + struct trunc_t { + static I + trunc(const T &val, + typename EpsilonType::Type epsilon = (DefaultEpsilon::value())) { + // this should be optimized away unless needed + if(!std::numeric_limits::is_signed) + // make sure this works for all useful cases even if I is an unsigned type + if(eq(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(lower+1), val, epsilon)) + return lower+1; + else return lower; + } + }; + template + struct trunc_t { + static I + trunc(const T &val, + typename EpsilonType::Type epsilon = (DefaultEpsilon::value())) { + I upper = trunc_t::trunc(val, epsilon); + if(ne(T(upper), val, epsilon)) ++upper; + return upper; + } + }; + template + struct trunc_t { + static I + trunc(const T &val, + typename EpsilonType::Type epsilon = (DefaultEpsilon::value())) { + if(val > T(0)) return trunc_t::trunc(val, epsilon); + else return trunc_t::trunc(val, epsilon); + } + }; + template + struct trunc_t { + static I + trunc(const T &val, + typename EpsilonType::Type epsilon = (DefaultEpsilon::value())) { + if(val > T(0)) return trunc_t::trunc(val, epsilon); + else return trunc_t::trunc(val, epsilon); + } + }; + template + struct trunc_t, std::vector, cstyle, rstyle> { + static std::vector + trunc(const std::vector &val, + typename EpsilonType::Type epsilon = (DefaultEpsilon::value())) { + unsigned int size = val.size(); + std::vector res(size); + for(unsigned int i = 0; i < size; ++i) + res[i] = trunc_t::trunc(val[i], epsilon); + return res; + } + }; + template + struct trunc_t, Dune::FieldVector, cstyle, rstyle> { + static Dune::FieldVector + trunc(const Dune::FieldVector &val, + typename EpsilonType::Type epsilon = (DefaultEpsilon::value())) { + Dune::FieldVector res; + for(int i = 0; i < n; ++i) + res[i] = trunc_t::trunc(val[i], epsilon); + return res; + } + }; + } // namespace Impl + template + I trunc(const T &val, typename EpsilonType::Type epsilon /*= DefaultEpsilon::value()*/) + { + return Impl::trunc_t::trunc(val, epsilon); + } + template + I trunc(const T &val, typename EpsilonType::Type epsilon = DefaultEpsilon::value()) + { + return trunc(val, epsilon); + } + template + I trunc(const T &val, typename EpsilonType::Type epsilon = DefaultEpsilon::value()) + { + return trunc(val, epsilon); + } + template + I trunc(const T &val, typename EpsilonType::Type epsilon = DefaultEpsilon::value()) + { + return trunc(val, epsilon); + } + } //namespace Dune + + // oo interface + template + FloatCmpOps:: + FloatCmpOps(EpsilonType epsilon) : epsilon_(epsilon) {} + + + template + typename FloatCmpOps::EpsilonType + FloatCmpOps::epsilon() const + { + return epsilon_; + } + + template + void + FloatCmpOps::epsilon(EpsilonType epsilon__) + { + epsilon_ = epsilon__; + } + + + template + bool FloatCmpOps:: + eq(const ValueType &first, const ValueType &second) const + { + return Dune::FloatCmp::eq(first, second, epsilon_); + } + + template + bool FloatCmpOps:: + ne(const ValueType &first, const ValueType &second) const + { + return Dune::FloatCmp::ne(first, second, epsilon_); + } + + template + bool FloatCmpOps:: + gt(const ValueType &first, const ValueType &second) const + { + return Dune::FloatCmp::gt(first, second, epsilon_); + } + + template + bool FloatCmpOps:: + lt(const ValueType &first, const ValueType &second) const + { + return Dune::FloatCmp::lt(first, second, epsilon_); + } + + template + bool FloatCmpOps:: + ge(const ValueType &first, const ValueType &second) const + { + return Dune::FloatCmp::ge(first, second, epsilon_); + } + + template + bool FloatCmpOps:: + le(const ValueType &first, const ValueType &second) const + { + return Dune::FloatCmp::le(first, second, epsilon_); + } + + + template + template + I FloatCmpOps:: + round(const ValueType &val) const + { + return Dune::FloatCmp::round(val, epsilon_); + } + + template + template + I FloatCmpOps:: + trunc(const ValueType &val) const + { + return Dune::FloatCmp::trunc(val, epsilon_); + } + +} //namespace Dune diff --git a/dune/common/float_cmp.hh b/dune/common/float_cmp.hh new file mode 100644 index 0000000000000000000000000000000000000000..3ac90672454cbeaf8682b9e13b51fe214d7a184b --- /dev/null +++ b/dune/common/float_cmp.hh @@ -0,0 +1,385 @@ +// -*- 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 + + @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 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 + struct DefaultEpsilon { + //! Returns the default epsilon for the given value type and compare style + static typename EpsilonType::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 + bool eq(const T &first, + const T &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::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 + bool ne(const T &first, + const T &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::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 + bool gt(const T &first, + const T &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::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 + bool lt(const T &first, + const T &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::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 + bool ge(const T &first, + const T &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::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 + bool le(const T &first, + const T &second, + typename EpsilonType::Type epsilon = DefaultEpsilon::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 + I round(const T &val, typename EpsilonType::Type epsilon = DefaultEpsilon::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 + I trunc(const T &val, typename EpsilonType::Type epsilon = DefaultEpsilon::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 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 + */ + typedef typename FloatCmp::EpsilonType::Type EpsilonType; + + private: + EpsilonType epsilon_; + + typedef FloatCmp::DefaultEpsilon 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 + 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 + I trunc(const ValueType &val) const; + + }; + +} //namespace Dune + +#include "float_cmp.cc" + +#endif //DUNE_COMMON_FLOAT_CMP_HH diff --git a/dune/common/quadmath.hh b/dune/common/quadmath.hh index ccbe88cd164ac63b07f3b149523faf0fa1b79cca..3aa1880344226c90384e237fa38822061318f1bc 100644 --- a/dune/common/quadmath.hh +++ b/dune/common/quadmath.hh @@ -15,6 +15,7 @@ #include #include +#include #include namespace Dune @@ -53,7 +54,7 @@ namespace Dune {} // constructor from pointer to null-terminated byte string - Float128(const char* str) noexcept + explicit Float128(const char* str) noexcept : value_(strtoflt128(str, NULL)) {} @@ -389,6 +390,21 @@ namespace Dune struct IsNumber : public std::true_type {}; + + namespace Impl + { + // Specialization of StrToNumber converter for Float128 type + template<> + struct StrToNumber + { + static Float128 eval (const char* first, const char* last) + { + auto parser = [](const char* str, char** end) { return strtoflt128(str, end); }; + return StrToNumberParser::eval(first,last,parser); + } + }; + + } // end namespace Impl } // end namespace Dune namespace std diff --git a/dune/common/strtonumber.hh b/dune/common/strtonumber.hh index 62bc4ac941c6be41e2d443fad9c015921482a4c4..09949e1db4fbaad6ee21d437e710351a47c5fa2c 100644 --- a/dune/common/strtonumber.hh +++ b/dune/common/strtonumber.hh @@ -33,7 +33,7 @@ namespace Dune { static T eval (const char* first, const char* /*last*/) { - return T(first);; + return T(first); } }; @@ -60,10 +60,16 @@ namespace Dune T value; auto result = Std::from_chars(first, last, value); + // require that all characters are consumed during conversion + char const* end = result.ptr; + bool all_consumed = (end != first); + while (all_consumed && (end != last) && (*end != '\0')) + all_consumed = std::isspace(*end++); + if (result.ec == std::errc::result_out_of_range) { DUNE_THROW(RangeError, std::error_condition(result.ec).message()); } - else if (result.ec == std::errc::invalid_argument) { + else if (result.ec == std::errc::invalid_argument || !all_consumed) { DUNE_THROW(InvalidArgument, "Conversion of '" << first << "' to number failed. Possible reason: invalid string or locale format."); } @@ -90,7 +96,7 @@ namespace Dune // The parser has the signature `T(const char*, char**)` and may set the errno // in case of a range error. template - T eval (const char* first, const char* /*last*/, Parser parser) + static T eval (const char* first, const char* last, Parser parser) { char* end; auto old_errno = errno; @@ -104,7 +110,7 @@ namespace Dune // test whether all non-space characters are consumed during conversion bool all_consumed = (end != first); - while (all_consumed && (*end != '\0')) + while (all_consumed && (end != last) && (*end != '\0')) all_consumed = std::isspace(*end++); if (!all_consumed) { @@ -115,7 +121,6 @@ namespace Dune return convertToRange(x); } - private: // Check whether a numeric conversion is safe template static T convertToRange (U const& u) diff --git a/dune/common/test/quadmathtest.cc b/dune/common/test/quadmathtest.cc index cd9381df24e807cc5e94c586bce26cee1c19466f..52d380a67cfd268fc7e7816e2c6e982a1558ec60 100644 --- a/dune/common/test/quadmathtest.cc +++ b/dune/common/test/quadmathtest.cc @@ -44,10 +44,10 @@ int main() Float128 x3 = 1.0; Float128 x4 = 1.0l; - int z1 = x1; - float z2 = x2; - double z3 = x3; - long double z4 = x4; + DUNE_UNUSED int z1 = x1; + DUNE_UNUSED float z2 = x2; + DUNE_UNUSED double z3 = x3; + DUNE_UNUSED long double z4 = x4; // field-vector FieldVector v{1,2,3}, x; @@ -78,8 +78,8 @@ int main() M.solve(v, x); // x = M^(-1)*v - auto M3 = M.leftmultiplyany(M2); - auto M4 = M.rightmultiplyany(M2); + DUNE_UNUSED auto M3 = M.leftmultiplyany(M2); + DUNE_UNUSED auto M4 = M.rightmultiplyany(M2); using namespace FMatrixHelp; diff --git a/dune/common/test/strtonumbertest.cc b/dune/common/test/strtonumbertest.cc index 38f97e397d4b3f505a6043e2eb2383f8ce275704..c98b933aa6737f170c1a519114f4d17dab8ae5de 100644 --- a/dune/common/test/strtonumbertest.cc +++ b/dune/common/test/strtonumbertest.cc @@ -121,19 +121,11 @@ int main() std::srand(std::time(nullptr)); using Types = std::tuple; TestSuite test; - test.check(strTo("0") == '\0'); - test.check(strTo("0") == '\0'); - test.check(strTo("0") == '\0'); - Hybrid::forEach(Types{}, [&test](auto value) { using namespace Dune; @@ -161,7 +153,7 @@ int main() bool exception = false; try { - T val3 = strTo(("1" + oss2.str()).c_str()); // add one more digit in front of number + DUNE_UNUSED T val3 = strTo(("1" + oss2.str()).c_str()); // add one more digit in front of number } catch(Dune::RangeError const&) { exception = true; } @@ -189,8 +181,8 @@ int main() test.check(cmp(value, val0), "Locale_C"); char* new_loc = setlocale(LC_NUMERIC, "de_DE.UTF-8"); - if (new_loc) { - T val1 = strTo("1,5"); + if (new_loc && Std::is_detected::value) { + T val1 = strTo("1.5"); test.check(cmp(value, val1), "Locale_de_DE"); } else { std::cout << "### skipped locale de_DE.UTF-8 test\n"; @@ -199,7 +191,7 @@ int main() setlocale(LC_NUMERIC, "C"); bool exception = false; try { - T val2 = strTo("1,5"); + DUNE_UNUSED T val2 = strTo("1,5"); } catch(Dune::InvalidArgument const&) { exception = true; } @@ -207,7 +199,7 @@ int main() exception = false; try { - T val2 = strTo("1.5__"); + DUNE_UNUSED T val2 = strTo("1.5__"); } catch(Dune::InvalidArgument const&) { exception = true; } @@ -215,7 +207,7 @@ int main() exception = false; try { - T val2 = strTo("1.5 "); + DUNE_UNUSED T val2 = strTo("1.5 "); } catch(Dune::InvalidArgument const&) { exception = true; } @@ -223,7 +215,7 @@ int main() exception = false; try { - T val2 = strTo(" 1.5"); + DUNE_UNUSED T val2 = strTo(" 1.5"); } catch(Dune::InvalidArgument const&) { exception = true; }