diff --git a/extensions/example/generic_robust_predicates/static_side_2d.cpp b/extensions/example/generic_robust_predicates/static_side_2d.cpp index c18ef3c4f16..04ba5b05153 100644 --- a/extensions/example/generic_robust_predicates/static_side_2d.cpp +++ b/extensions/example/generic_robust_predicates/static_side_2d.cpp @@ -56,7 +56,7 @@ struct side_robust_with_static_filter bg::get<1>(p2), bg::get<0>(p), bg::get<1>(p)); - if(sign != bg::detail::generic_robust_predicates::sign_uncertain) + if (sign != bg::detail::generic_robust_predicates::sign_uncertain) { return sign; } @@ -68,7 +68,7 @@ struct side_robust_with_static_filter } }; -int main(int argc, char** argv) +int main() { point p1(0.0, 0.0); point p2(1.0, 1.0); diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/almost_static_filter.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/almost_static_filter.hpp index 6dcfd615dc6..7727c1c4bc7 100644 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/almost_static_filter.hpp +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/almost_static_filter.hpp @@ -28,7 +28,7 @@ template struct make_filter_impl { template - static Filter apply(const ExtremaArray& extrema, const Reals&... args) + static Filter apply(ExtremaArray const& extrema, Reals const&... args) { return make_filter_impl ::apply(extrema, args..., extrema[N]); @@ -39,13 +39,27 @@ template struct make_filter_impl { template - static Filter apply(const ExtremaArray& extrema, const Reals&... args) + static Filter apply(ExtremaArray const& extrema, Reals const&... args) { Filter f(args...); return f; } }; +//The almost static filter holds an instance of a static filter and applies it +//when it is called. Its static filter can be updated and applied like this: +// +//almost_static_filter<...> f; +// +//f.update(max1, max2, ..., min1, min2, ...); +// +//f.apply(arg1, arg2, ...); +// +//Unlike with the static filter, global bounds do not need to be known at +//construction time and with incremental algorithms where inputs with higher +//magnitude are added later, the earlier applications of the filter may benefit +//from more optimistic error bounds. + template < typename Expression, @@ -72,13 +86,13 @@ class almost_static_filter std::numeric_limits::infinity()); } template - int apply(const Reals&... args) const + int apply(Reals const&... args) const { return m_filter.apply(args...); } template - inline void update_extrema(const Reals&... args) + inline void update_extrema(Reals const&... args) { std::array input {{ static_cast(args)... }}; for(int i = 0; i < m_extrema.size() / 2; ++i) @@ -92,7 +106,7 @@ class almost_static_filter } template - inline bool update_extrema_check(const Reals&... args) + inline bool update_extrema_check(Reals const&... args) { bool changed = false; std::array input {{ static_cast(args)... }}; diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp index 65d4f100bf9..53dabca59ea 100644 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp @@ -28,6 +28,10 @@ namespace boost { namespace geometry namespace detail { namespace generic_robust_predicates { +//The templates in this file are meant to be used for the evaluation of +//expressions with floating-point precision and floating-point rounding. +//The most important template in this file is approximate_interim. + template < typename Node, @@ -37,7 +41,7 @@ template > struct get_nth_real_impl { - static inline Real apply(const InputArr& input) + static inline Real apply(InputArr const& input) { return input[N - 1]; } @@ -51,7 +55,7 @@ template > struct get_nth_real_impl { - static inline Real apply(const InputArr&) + static inline Real apply(InputArr const&) { return Node::value; } @@ -64,7 +68,7 @@ template typename Real, typename InputArr > -inline Real get_nth_real(const InputArr& input) +inline Real get_nth_real(InputArr const& input) { return get_nth_real_impl::apply(input); } @@ -80,7 +84,7 @@ template > struct get_approx_impl { - static inline Real apply(Arr& interim_results, const InputArr& input) + static inline Real apply(Arr& interim_results, InputArr const&) { return interim_results[boost::mp11::mp_find::value]; } @@ -96,7 +100,7 @@ template > struct get_approx_impl { - static inline Real apply(Arr& interim_results, const InputArr& input) + static inline Real apply(Arr&, InputArr const& input) { return get_nth_real(input); } @@ -110,7 +114,7 @@ template typename Arr, typename InputArr > -inline Real get_approx(Arr& interim_results, const InputArr& input) +inline Real get_approx(Arr& interim_results, InputArr const& input) { return get_approx_impl < @@ -145,7 +149,7 @@ template > struct approximate_remainder_impl { - static inline void apply(Arr& interim_results, const InputArr& input) + static inline void apply(Arr& interim_results, InputArr const& input) { using node = boost::mp11::mp_front; approximate_interim_impl @@ -178,7 +182,7 @@ struct approximate_remainder_impl InputArr > { - static inline void apply(Arr& interim_results, const InputArr& input) {} + static inline void apply(Arr&, InputArr const&) {} }; template @@ -189,7 +193,7 @@ template typename Arr, typename InputArr > -inline void approximate_remainder(Arr& interim_results, const InputArr& input) +inline void approximate_remainder(Arr& interim_results, InputArr const& input) { approximate_remainder_impl < @@ -220,7 +224,7 @@ struct approximate_interim_impl InputArr > { - static inline void apply(Arr& interim_results, const InputArr& input) + static inline void apply(Arr& interim_results, InputArr const& input) { using node = boost::mp11::mp_front; interim_results[boost::mp11::mp_find::value] = @@ -255,7 +259,7 @@ struct approximate_interim_impl InputArr > { - static inline void apply(Arr& interim_results, const InputArr& input) + static inline void apply(Arr& interim_results, InputArr const& input) { using node = boost::mp11::mp_front; interim_results[boost::mp11::mp_find::value] = std::max( @@ -290,7 +294,7 @@ struct approximate_interim_impl InputArr > { - static inline void apply(Arr& interim_results, const InputArr& input) + static inline void apply(Arr& interim_results, InputArr const& input) { using node = boost::mp11::mp_front; interim_results[boost::mp11::mp_find::value] = std::min( @@ -325,7 +329,7 @@ struct approximate_interim_impl InputArr > { - static inline void apply(Arr& interim_results, const InputArr& input) + static inline void apply(Arr& interim_results, InputArr const& input) { using node = boost::mp11::mp_front; interim_results[boost::mp11::mp_find::value] = @@ -360,7 +364,7 @@ struct approximate_interim_impl InputArr > { - static inline void apply(Arr& interim_results, const InputArr& input) + static inline void apply(Arr& interim_results, InputArr const& input) { using node = boost::mp11::mp_front; interim_results[boost::mp11::mp_find::value] = @@ -395,7 +399,7 @@ struct approximate_interim_impl InputArr > { - static inline void apply(Arr& interim_results, const InputArr& input) + static inline void apply(Arr& interim_results, InputArr const& input) { using node = boost::mp11::mp_front; interim_results[boost::mp11::mp_find::value] = @@ -432,7 +436,7 @@ struct approximate_interim_impl InputArr > { - static inline void apply(Arr& interim_results, const InputArr& input) + static inline void apply(Arr& interim_results, InputArr const& input) { approximate_remainder < @@ -443,6 +447,13 @@ struct approximate_interim_impl } }; +//All expects an boost::mp11::mp_list of all expressions that need to be +//evaluated. Remaining expects an boost::mp11::mp_list of the expressions +//that are left to be evaluated. In the first call it is expected to be +//equal to All and it serves as an anchor for template recursion. Real is a +//floating-point type. The remaining template arguments are deduced from +//parameters. + template < typename All, @@ -451,7 +462,7 @@ template typename Arr, typename InputArr > -inline void approximate_interim(Arr& interim_results, const InputArr& input) +inline void approximate_interim(Arr& interim_results, InputArr const& input) { approximate_remainder < @@ -462,7 +473,7 @@ inline void approximate_interim(Arr& interim_results, const InputArr& input) } template -inline Real approximate_value(const InputArr& input) +inline Real approximate_value(InputArr const& input) { using stack = typename boost::mp11::mp_unique>; using evals = typename boost::mp11::mp_remove_if; diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp index 46e39a0acf9..331f983e2f3 100644 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/coefficient_list.hpp @@ -25,6 +25,18 @@ namespace boost { namespace geometry namespace detail { namespace generic_robust_predicates { +//The templates in this file are helpers for the computation of error maps and +//error expressions (see also error_bound.hpp). A coefficient list is a +//boost::mp11::mp_list that contains integers that represent the coefficients +//of an polynomial. This is used for the polynomials in epsilon that occur in +//forward error analysis, e.g. epsilon + 2 * epsilon^2 + epsilon^3 would be +//represented by the type mp_list, mp_int<2>, mp_int<1>>. +// +//See p. 39 in "Adaptive Precision Floating-Point Arithmetic and Fast Robust +//Geometric Predicates" by Jonathan Richard Shewchuk (can be downloaded from +//https://www.cs.cmu.edu/~quake/robust.html) for the inspiration for the +//methods implemented below. + //only meant to be used in assertions template using is_mp_int = boost::mp11::mp_same>; @@ -68,12 +80,21 @@ struct coeff_truncate_impl static_assert(is_coeff_list::value, "type should be a coefficient list"); }; +//Because of the way floating-point numbers are rounded, for a given epsilon +//polynomial we only care about the first two non-zero coefficients. +//e.g. if a*eps + b*eps^2 + c*eps^3 is considered then c * eps^3 is always +//neglible unless c is so large that it must have resulted from the analysis +//of an expression that is so complex that it cannot be feasibly processed +//by the methods implemented in this feature. template using coeff_truncate = typename coeff_truncate_impl::type; template using app_zero_b = boost::mp11::mp_push_front>; template using app_zero_f = boost::mp11::mp_push_back>; + +//For a coefficient list representing the polynomial p, the following template +//produces the coefficient list representing p * (1 + eps). template using mult_by_1_p_eps = coeff_truncate < boost::mp11::mp_transform @@ -110,6 +131,9 @@ struct mult_by_1_p_eps_pow_impl using type = L; }; + +//For a coefficient list representing the polynomial p, the following template +//produces the cofficient list representing p * (1 + eps)^N template using mult_by_1_p_eps_pow = coeff_truncate < typename mult_by_1_p_eps_pow_impl::type @@ -123,6 +147,9 @@ template using div_by_1_m_eps_helper = boost::mp11::mp_plus >; +//For a coefficient list representing a polynomial p in eps, the following +//template computes a coefficient list representing a polynomial q such that +//q * (1 - p) >= p. template using div_by_1_m_eps = coeff_truncate < boost::mp11::mp_push_back @@ -406,6 +433,11 @@ struct coeff_round_impl { using type = L; }; template struct coeff_round_impl { using type = L; }; +//The following template rounds a coefficient list up such that the +//corresponding polynomial evaluated for epsilon is a representable +//floating-point number. +// +//e.g. 3 * eps + 12 * eps^2 + 1 * eps^3 is rounded to 3 * eps + 16 * eps^2 template using coeff_round = typename coeff_round_impl::type; template using indexed_value_product = diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/error_bound.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/error_bound.hpp index 6df2bafda69..d42a519aa93 100644 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/error_bound.hpp +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/error_bound.hpp @@ -26,6 +26,24 @@ namespace boost { namespace geometry namespace detail { namespace generic_robust_predicates { +//The templates in this file are utilities for the generation of error +//expressions for given arithmetic expressions. This is done using forward +//error analysis methods at compile time. This analysis produces an error map +//that is folded into an error expression. +// +//An Error Map is a template fulfilling the boost::mp11 map concept. The keys K +//are expressions (as in expression trees and the values V are lists of +//integers. Each key-value pair represents the contribution of an error term to +//the total error, i.e. the error expression is something like the sum of +//abs * (V[0] * epsilon + V[1] * epsilon^2 + ...) +// +//over all key-value pairs KV in the error map. +// +//See p. 39 in "Adaptive Precision Floating-Point Arithmetic and Fast Robust +//Geometric Predicates" by Jonathan Richard Shewchuk (can be downloaded from +//https://www.cs.cmu.edu/~quake/robust.html) for the inspiration for the +//methods implemented below. + template using second_is_coeff_list = is_coeff_list< boost::mp11::mp_second >; diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp index 4ac3e558b07..91bbd6d3d8b 100644 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp @@ -25,6 +25,19 @@ namespace boost { namespace geometry namespace detail { namespace generic_robust_predicates { +//The templates in this file are designed to represent arithmetic expressions +//In the type system so that they can be processed at compile-time for the +//automatic creation of floating-point filters. +// +//E.g. the arithmetic expression a * b + c can be represented by the type +// +//sum< product < argument<1>, argument<2> >, argument<3> > +// +//The arguments represent the input variables and are the leaves of the +//expression tree (other possible leaves are constants). Because they represent +//placeholders their indexing is one-based rather than zero-based, similar to +//the placeholders for std::bind. + enum class operator_types { sum, difference, product, abs, no_op, max, min }; @@ -144,6 +157,12 @@ struct static_constant_interface : public leaf template using is_leaf = boost::mp11::mp_bool; +//post_order_impl and post_order are templates for compile-time traversion of +//expression trees. The post_order traversal was chosen because it guarantees +//that children (subexpressions) are evaluated before there parents as is +//necessary for the evaluation of the arithmetic expressions that these +//expression trees represent. + template < typename In, diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expressions.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expressions.hpp index 6042b1834ec..b586605eedd 100644 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expressions.hpp +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expressions.hpp @@ -32,10 +32,10 @@ using det2x2 = difference >; using orient2d = det2x2 - < - difference <_1, _5>, difference<_2, _6>, - difference <_3, _5>, difference<_4, _6> - >; + < + difference <_1, _5>, difference<_2, _6>, + difference <_3, _5>, difference<_4, _6> + >; template < diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/interval_error_bound.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/interval_error_bound.hpp index 59e4d3f47b9..9a9de3d1ff0 100644 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/interval_error_bound.hpp +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/interval_error_bound.hpp @@ -20,6 +20,13 @@ namespace boost { namespace geometry namespace detail { namespace generic_robust_predicates { +//This file contains methods to manipulate expression trees based on the ideas +//found in https://en.wikipedia.org/wiki/Interval_arithmetic. The purpose of +//this is to transform error expressions for semi-static filters, i.e. error +//expressions that compute error bounds for specific inputs, to error +//expressions for static filters that compute error bounds for upper and lower +//bounds on input values. + template struct interval_min_impl { diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/semi_static_filter.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/semi_static_filter.hpp index 26ca2884fa9..43095caab68 100644 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/semi_static_filter.hpp +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/semi_static_filter.hpp @@ -27,6 +27,22 @@ namespace boost { namespace geometry namespace detail { namespace generic_robust_predicates { +//The semi static filter evaluates an Expression and an ErrorExpression. On +//application it evaluates both expressions and verifies the sign of the +//Expression by the following rules: +//If the absolute value of the number resulting from the evaluation of the +//expression is larger than or equal to the result for the error expression, +//then the sign can be returned. If not, then a constant is returned to +//represent an uncertain sign. +// +//The filters that are build from this template are meant to be semi static in +//the sense that the error expression including constants that depend only the +//epsilon of the calculation type are known statically at compile-time but the +//final value of the error bound depends on the specific inputs for each call +//to the filter and parts of it need to be computed at each call. +// +//This filter is stateless. + template < typename Expression, @@ -55,7 +71,7 @@ struct semi_static_filter using ct = CalculationType; public: template - static inline int apply(const Reals&... args) + static inline int apply(Reals const&... args) { std::array input {{ static_cast(args)... }}; diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a.hpp index 96725db461a..9fafd9f13ba 100644 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a.hpp +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a.hpp @@ -31,6 +31,9 @@ namespace boost { namespace geometry namespace detail { namespace generic_robust_predicates { +//The following templates apply the ideas implemented in error_bound.hpp and +//interval_error_bound.hpp to generate error expressions from expressions. + template struct stage_a_error_bound_impl { diff --git a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_filter.hpp b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_filter.hpp index d93eb6497f9..6f88418a7b9 100644 --- a/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_filter.hpp +++ b/include/boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/static_filter.hpp @@ -27,6 +27,16 @@ namespace boost { namespace geometry namespace detail { namespace generic_robust_predicates { +//The static filter works similar to the semi static filter with the exception +//that the error bound is only computed once at construction time. For this +//purpose, the filter is stateful and it requires upper and lower bounds on the +//inputs at compile time like this: +// +//static_filter<...> f(max1, max2, ..., min1, min2, ...); +// +//for each call to the filter, i.e. f.apply(arg1, arg2, ...); it must hold that +//min1 <= arg1 <= max1 and so on. + template < typename Expression, @@ -50,7 +60,7 @@ class static_filter ct error_bound() const { return m_error_bound; } template - inline static_filter(const Reals&... args) + inline static_filter(Reals const&... args) : m_error_bound(approximate_value( std::array {static_cast(args)...})) @@ -60,7 +70,7 @@ class static_filter } template - inline int apply(const Reals&... args) const + inline int apply(Reals const&... args) const { std::array input {static_cast(args)...}; std::array::value> results; @@ -85,7 +95,7 @@ class static_filter } template - inline int operator()(const Reals&... args) const + inline int operator()(Reals const&... args) const { return apply(args...); }