boost/geometry/util/precise_math.hpp
// Boost.Geometry (aka GGL, Generic Geometry Library) // Copyright (c) 2019 Tinko Bartels, Berlin, Germany. // Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland. // Contributed and/or modified by Tinko Bartels, // as part of Google Summer of Code 2019 program. // This file was modified by Oracle on 2021. // Modifications copyright (c) 2021, Oracle and/or its affiliates. // Contributed and/or modified by Vissarion Fisikopoulos, on behalf of Oracle // Use, modification and distribution is subject to the Boost Software License, // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) #ifndef BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_CARTESIAN_DETAIL_PRECISE_MATH_HPP #define BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_CARTESIAN_DETAIL_PRECISE_MATH_HPP #include<numeric> #include<cmath> #include<limits> #include<array> #include <boost/geometry/core/access.hpp> #include <boost/geometry/util/condition.hpp> // The following code is based on "Adaptive Precision Floating-Point Arithmetic // and Fast Robust Geometric Predicates" by Richard Shewchuk, // J. Discrete Comput Geom (1997) 18: 305. https://doi.org/10.1007/PL00009321 namespace boost { namespace geometry { namespace detail { namespace precise_math { // See Theorem 6, page 6 template < typename RealNumber > inline std::array<RealNumber, 2> fast_two_sum(RealNumber const a, RealNumber const b) { RealNumber x = a + b; RealNumber b_virtual = x - a; return {{x, b - b_virtual}}; } // See Theorem 7, page 7 - 8 template < typename RealNumber > inline std::array<RealNumber, 2> two_sum(RealNumber const a, RealNumber const b) { RealNumber x = a + b; RealNumber b_virtual = x - a; RealNumber a_virtual = x - b_virtual; RealNumber b_roundoff = b - b_virtual; RealNumber a_roundoff = a - a_virtual; RealNumber y = a_roundoff + b_roundoff; return {{ x, y }}; } // See bottom of page 8 template < typename RealNumber > inline RealNumber two_diff_tail(RealNumber const a, RealNumber const b, RealNumber const x) { RealNumber b_virtual = a - x; RealNumber a_virtual = x + b_virtual; RealNumber b_roundoff = b_virtual - b; RealNumber a_roundoff = a - a_virtual; return a_roundoff + b_roundoff; } // see bottom of page 8 template < typename RealNumber > inline std::array<RealNumber, 2> two_diff(RealNumber const a, RealNumber const b) { RealNumber x = a - b; RealNumber y = two_diff_tail(a, b, x); return {{ x, y }}; } // see theorem 18, page 19 template < typename RealNumber > inline RealNumber two_product_tail(RealNumber const a, RealNumber const b, RealNumber const x) { return std::fma(a, b, -x); } // see theorem 18, page 19 template < typename RealNumber > inline std::array<RealNumber, 2> two_product(RealNumber const a, RealNumber const b) { RealNumber x = a * b; RealNumber y = two_product_tail(a, b, x); return {{ x , y }}; } // see theorem 12, figure 7, page 11 - 12, // this is the 2 by 2 case for the corresponding diff-method // note that this method takes input in descending order of magnitude and // returns components in ascending order of magnitude template < typename RealNumber > inline std::array<RealNumber, 4> two_two_expansion_diff( std::array<RealNumber, 2> const a, std::array<RealNumber, 2> const b) { std::array<RealNumber, 4> h; std::array<RealNumber, 2> Qh = two_diff(a[1], b[1]); h[0] = Qh[1]; Qh = two_sum( a[0], Qh[0] ); RealNumber _j = Qh[0]; Qh = two_diff(Qh[1], b[0]); h[1] = Qh[1]; Qh = two_sum( _j, Qh[0] ); h[2] = Qh[1]; h[3] = Qh[0]; return h; } // see theorem 13, figure 8. This implementation uses zero elimination as // suggested on page 17, second to last paragraph. Returns the number of // non-zero components in the result and writes the result to h. // the merger into a single sequence g is done implicitly template < typename RealNumber, std::size_t InSize1, std::size_t InSize2, std::size_t OutSize > inline int fast_expansion_sum_zeroelim( std::array<RealNumber, InSize1> const& e, std::array<RealNumber, InSize2> const& f, std::array<RealNumber, OutSize> & h, int m = InSize1, int n = InSize2) { std::array<RealNumber, 2> Qh; int i_e = 0; int i_f = 0; int i_h = 0; if (std::abs(f[0]) > std::abs(e[0])) { Qh[0] = e[i_e++]; } else { Qh[0] = f[i_f++]; } i_h = 0; if ((i_e < m) && (i_f < n)) { if (std::abs(f[i_f]) > std::abs(e[i_e])) { Qh = fast_two_sum(e[i_e++], Qh[0]); } else { Qh = fast_two_sum(f[i_f++], Qh[0]); } if (Qh[1] != 0.0) { h[i_h++] = Qh[1]; } while ((i_e < m) && (i_f < n)) { if (std::abs(f[i_f]) > std::abs(e[i_e])) { Qh = two_sum(Qh[0], e[i_e++]); } else { Qh = two_sum(Qh[0], f[i_f++]); } if (Qh[1] != 0.0) { h[i_h++] = Qh[1]; } } } while (i_e < m) { Qh = two_sum(Qh[0], e[i_e++]); if (Qh[1] != 0.0) { h[i_h++] = Qh[1]; } } while (i_f < n) { Qh = two_sum(Qh[0], f[i_f++]); if (Qh[1] != 0.0) { h[i_h++] = Qh[1]; } } if ((Qh[0] != 0.0) || (i_h == 0)) { h[i_h++] = Qh[0]; } return i_h; } // see theorem 19, figure 13, page 20 - 21. This implementation uses zero // elimination as suggested on page 17, second to last paragraph. Returns the // number of non-zero components in the result and writes the result to h. template < typename RealNumber, std::size_t InSize > inline int scale_expansion_zeroelim( std::array<RealNumber, InSize> const& e, RealNumber const b, std::array<RealNumber, 2 * InSize> & h, int e_non_zeros = InSize) { std::array<RealNumber, 2> Qh = two_product(e[0], b); int i_h = 0; if (Qh[1] != 0) { h[i_h++] = Qh[1]; } for (int i_e = 1; i_e < e_non_zeros; i_e++) { std::array<RealNumber, 2> Tt = two_product(e[i_e], b); Qh = two_sum(Qh[0], Tt[1]); if (Qh[1] != 0) { h[i_h++] = Qh[1]; } Qh = fast_two_sum(Tt[0], Qh[0]); if (Qh[1] != 0) { h[i_h++] = Qh[1]; } } if ((Qh[0] != 0.0) || (i_h == 0)) { h[i_h++] = Qh[0]; } return i_h; } template<typename RealNumber> struct vec2d { RealNumber x; RealNumber y; }; template < typename RealNumber, std::size_t Robustness > inline RealNumber orient2dtail(vec2d<RealNumber> const& p1, vec2d<RealNumber> const& p2, vec2d<RealNumber> const& p3, std::array<RealNumber, 2>& t1, std::array<RealNumber, 2>& t2, std::array<RealNumber, 2>& t3, std::array<RealNumber, 2>& t4, std::array<RealNumber, 2>& t5_01, std::array<RealNumber, 2>& t6_01, RealNumber const& magnitude) { t5_01[1] = two_product_tail(t1[0], t2[0], t5_01[0]); t6_01[1] = two_product_tail(t3[0], t4[0], t6_01[0]); std::array<RealNumber, 4> tA_03 = two_two_expansion_diff(t5_01, t6_01); RealNumber det = std::accumulate(tA_03.begin(), tA_03.end(), static_cast<RealNumber>(0)); if (BOOST_GEOMETRY_CONDITION(Robustness == 1)) { return det; } // see p.39, mind the different definition of epsilon for error bound RealNumber B_relative_bound = (1 + 3 * std::numeric_limits<RealNumber>::epsilon()) * std::numeric_limits<RealNumber>::epsilon(); RealNumber absolute_bound = B_relative_bound * magnitude; if (std::abs(det) >= absolute_bound) { return det; //B estimate } t1[1] = two_diff_tail(p1.x, p3.x, t1[0]); t2[1] = two_diff_tail(p2.y, p3.y, t2[0]); t3[1] = two_diff_tail(p1.y, p3.y, t3[0]); t4[1] = two_diff_tail(p2.x, p3.x, t4[0]); if ((t1[1] == 0) && (t3[1] == 0) && (t2[1] == 0) && (t4[1] == 0)) { return det; //If all tails are zero, there is noething else to compute } RealNumber sub_bound = (1.5 + 2 * std::numeric_limits<RealNumber>::epsilon()) * std::numeric_limits<RealNumber>::epsilon(); // see p.39, mind the different definition of epsilon for error bound RealNumber C_relative_bound = (2.25 + 8 * std::numeric_limits<RealNumber>::epsilon()) * std::numeric_limits<RealNumber>::epsilon() * std::numeric_limits<RealNumber>::epsilon(); absolute_bound = C_relative_bound * magnitude + sub_bound * std::abs(det); det += (t1[0] * t2[1] + t2[0] * t1[1]) - (t3[0] * t4[1] + t4[0] * t3[1]); if (BOOST_GEOMETRY_CONDITION(Robustness == 2) || std::abs(det) >= absolute_bound) { return det; //C estimate } std::array<RealNumber, 8> D_left; int D_left_nz; { std::array<RealNumber, 2> t5_23 = two_product(t1[1], t2[0]); std::array<RealNumber, 2> t6_23 = two_product(t3[1], t4[0]); std::array<RealNumber, 4> tA_47 = two_two_expansion_diff(t5_23, t6_23); D_left_nz = fast_expansion_sum_zeroelim(tA_03, tA_47, D_left); } std::array<RealNumber, 8> D_right; int D_right_nz; { std::array<RealNumber, 2> t5_45 = two_product(t1[0], t2[1]); std::array<RealNumber, 2> t6_45 = two_product(t3[0], t4[1]); std::array<RealNumber, 4> tA_8_11 = two_two_expansion_diff(t5_45, t6_45); std::array<RealNumber, 2> t5_67 = two_product(t1[1], t2[1]); std::array<RealNumber, 2> t6_67 = two_product(t3[1], t4[1]); std::array<RealNumber, 4> tA_12_15 = two_two_expansion_diff(t5_67, t6_67); D_right_nz = fast_expansion_sum_zeroelim(tA_8_11, tA_12_15, D_right); } std::array<RealNumber, 16> D; int D_nz = fast_expansion_sum_zeroelim(D_left, D_right, D, D_left_nz, D_right_nz); // only return component of highest magnitude because we mostly care about the sign. return(D[D_nz - 1]); } // see page 38, Figure 21 for the calculations, notation follows the notation // in the figure. template < typename RealNumber, std::size_t Robustness = 3, typename EpsPolicy > inline RealNumber orient2d(vec2d<RealNumber> const& p1, vec2d<RealNumber> const& p2, vec2d<RealNumber> const& p3, EpsPolicy& eps_policy) { std::array<RealNumber, 2> t1, t2, t3, t4; t1[0] = p1.x - p3.x; t2[0] = p2.y - p3.y; t3[0] = p1.y - p3.y; t4[0] = p2.x - p3.x; eps_policy = EpsPolicy(t1[0], t2[0], t3[0], t4[0]); std::array<RealNumber, 2> t5_01, t6_01; t5_01[0] = t1[0] * t2[0]; t6_01[0] = t3[0] * t4[0]; RealNumber det = t5_01[0] - t6_01[0]; if (BOOST_GEOMETRY_CONDITION(Robustness == 0)) { return det; } RealNumber const magnitude = std::abs(t5_01[0]) + std::abs(t6_01[0]); // see p.39, mind the different definition of epsilon for error bound RealNumber const A_relative_bound = (1.5 + 4 * std::numeric_limits<RealNumber>::epsilon()) * std::numeric_limits<RealNumber>::epsilon(); RealNumber absolute_bound = A_relative_bound * magnitude; if ( std::abs(det) >= absolute_bound ) { return det; //A estimate } if ( (t5_01[0] > 0 && t6_01[0] <= 0) || (t5_01[0] < 0 && t6_01[0] >= 0) ) { //if diagonal and antidiagonal have different sign, the sign of det is //obvious return det; } return orient2dtail<RealNumber, Robustness>(p1, p2, p3, t1, t2, t3, t4, t5_01, t6_01, magnitude); } // This method adaptively computes increasingly precise approximations of the following // determinant using Laplace expansion along the last column. // det A = // | p1_x - p4_x p1_y - p4_y ( p1_x - p4_x ) ^ 2 + ( p1_y - p4_y ) ^ 2 | // | p2_x - p4_x p2_y - p4_y ( p2_x - p4_x ) ^ 2 + ( p1_y - p4_y ) ^ 2 | // | p3_x - p4_x p3_y - p4_y ( p3_x - p4_x ) ^ 2 + ( p3_y - p4_y ) ^ 2 | // = a_13 * C_13 + a_23 * C_23 + a_33 * C_33 // where a_ij is the i-j-entry and C_ij is the i_j Cofactor template < typename RealNumber, std::size_t Robustness = 2 > RealNumber incircle(std::array<RealNumber, 2> const& p1, std::array<RealNumber, 2> const& p2, std::array<RealNumber, 2> const& p3, std::array<RealNumber, 2> const& p4) { RealNumber A_11 = p1[0] - p4[0]; RealNumber A_21 = p2[0] - p4[0]; RealNumber A_31 = p3[0] - p4[0]; RealNumber A_12 = p1[1] - p4[1]; RealNumber A_22 = p2[1] - p4[1]; RealNumber A_32 = p3[1] - p4[1]; std::array<RealNumber, 2> A_21_x_A_32, A_31_x_A_22, A_31_x_A_12, A_11_x_A_32, A_11_x_A_22, A_21_x_A_12; A_21_x_A_32[0] = A_21 * A_32; A_31_x_A_22[0] = A_31 * A_22; RealNumber A_13 = A_11 * A_11 + A_12 * A_12; A_31_x_A_12[0] = A_31 * A_12; A_11_x_A_32[0] = A_11 * A_32; RealNumber A_23 = A_21 * A_21 + A_22 * A_22; A_11_x_A_22[0] = A_11 * A_22; A_21_x_A_12[0] = A_21 * A_12; RealNumber A_33 = A_31 * A_31 + A_32 * A_32; RealNumber det = A_13 * (A_21_x_A_32[0] - A_31_x_A_22[0]) + A_23 * (A_31_x_A_12[0] - A_11_x_A_32[0]) + A_33 * (A_11_x_A_22[0] - A_21_x_A_12[0]); if (BOOST_GEOMETRY_CONDITION(Robustness == 0)) { return det; } RealNumber magnitude = (std::abs(A_21_x_A_32[0]) + std::abs(A_31_x_A_22[0])) * A_13 + (std::abs(A_31_x_A_12[0]) + std::abs(A_11_x_A_32[0])) * A_23 + (std::abs(A_11_x_A_22[0]) + std::abs(A_21_x_A_12[0])) * A_33; RealNumber A_relative_bound = (5 + 24 * std::numeric_limits<RealNumber>::epsilon()) * std::numeric_limits<RealNumber>::epsilon(); RealNumber absolute_bound = A_relative_bound * magnitude; if (std::abs(det) > absolute_bound) { return det; } // (p2_x - p4_x) * (p3_y - p4_y) A_21_x_A_32[1] = two_product_tail(A_21, A_32, A_21_x_A_32[0]); // (p3_x - p4_x) * (p2_y - p4_y) A_31_x_A_22[1] = two_product_tail(A_31, A_22, A_31_x_A_22[0]); // (bx - dx) * (cy - dy) - (cx - dx) * (by - dy) std::array<RealNumber, 4> C_13 = two_two_expansion_diff(A_21_x_A_32, A_31_x_A_22); std::array<RealNumber, 8> C_13_x_A11; // ( (bx - dx) * (cy - dy) - (cx - dx) * (by - dy) ) * ( ax - dx ) int C_13_x_A11_nz = scale_expansion_zeroelim(C_13, A_11, C_13_x_A11); std::array<RealNumber, 16> C_13_x_A11_sq; // ( (bx - dx) * (cy - dy) - (cx - dx) * (by - dy) ) * ( ax - dx ) * (ax - dx) int C_13_x_A11_sq_nz = scale_expansion_zeroelim(C_13_x_A11, A_11, C_13_x_A11_sq, C_13_x_A11_nz); std::array<RealNumber, 8> C_13_x_A12; // ( (bx - dx) * (cy - dy) - (cx - dx) * (by - dy) ) * ( ay - dy ) int C_13_x_A12_nz = scale_expansion_zeroelim(C_13, A_12, C_13_x_A12); std::array<RealNumber, 16> C_13_x_A12_sq; // ( (bx - dx) * (cy - dy) - (cx - dx) * (by - dy) ) * ( ay - dy ) * ( ay - dy ) int C_13_x_A12_sq_nz = scale_expansion_zeroelim(C_13_x_A12, A_12, C_13_x_A12_sq, C_13_x_A12_nz); std::array<RealNumber, 32> A_13_x_C13; // ( (bx - dx) * (cy - dy) - (cx - dx) * (by - dy) ) // * ( ( ay - dy ) * ( ay - dy ) + ( ax - dx ) * (ax - dx) ) int A_13_x_C13_nz = fast_expansion_sum_zeroelim(C_13_x_A11_sq, C_13_x_A12_sq, A_13_x_C13, C_13_x_A11_sq_nz, C_13_x_A12_sq_nz); // (cx - dx) * (ay - dy) A_31_x_A_12[1] = two_product_tail(A_31, A_12, A_31_x_A_12[0]); // (ax - dx) * (cy - dy) A_11_x_A_32[1] = two_product_tail(A_11, A_32, A_11_x_A_32[0]); // (cx - dx) * (ay - dy) - (ax - dx) * (cy - dy) std::array<RealNumber, 4> C_23 = two_two_expansion_diff(A_31_x_A_12, A_11_x_A_32); std::array<RealNumber, 8> C_23_x_A_21; // ( (cx - dx) * (ay - dy) - (ax - dx) * (cy - dy) ) * ( bx - dx ) int C_23_x_A_21_nz = scale_expansion_zeroelim(C_23, A_21, C_23_x_A_21); std::array<RealNumber, 16> C_23_x_A_21_sq; // ( (cx - dx) * (ay - dy) - (ax - dx) * (cy - dy) ) * ( bx - dx ) * ( bx - dx ) int C_23_x_A_21_sq_nz = scale_expansion_zeroelim(C_23_x_A_21, A_21, C_23_x_A_21_sq, C_23_x_A_21_nz); std::array<RealNumber, 8> C_23_x_A_22; // ( (cx - dx) * (ay - dy) - (ax - dx) * (cy - dy) ) * ( by - dy ) int C_23_x_A_22_nz = scale_expansion_zeroelim(C_23, A_22, C_23_x_A_22); std::array<RealNumber, 16> C_23_x_A_22_sq; // ( (cx - dx) * (ay - dy) - (ax - dx) * (cy - dy) ) * ( by - dy ) * ( by - dy ) int C_23_x_A_22_sq_nz = scale_expansion_zeroelim(C_23_x_A_22, A_22, C_23_x_A_22_sq, C_23_x_A_22_nz); std::array<RealNumber, 32> A_23_x_C_23; // ( (cx - dx) * (ay - dy) - (ax - dx) * (cy - dy) ) // * ( ( bx - dx ) * ( bx - dx ) + ( by - dy ) * ( by - dy ) ) int A_23_x_C_23_nz = fast_expansion_sum_zeroelim(C_23_x_A_21_sq, C_23_x_A_22_sq, A_23_x_C_23, C_23_x_A_21_sq_nz, C_23_x_A_22_sq_nz); // (ax - dx) * (by - dy) A_11_x_A_22[1] = two_product_tail(A_11, A_22, A_11_x_A_22[0]); // (bx - dx) * (ay - dy) A_21_x_A_12[1] = two_product_tail(A_21, A_12, A_21_x_A_12[0]); // (ax - dx) * (by - dy) - (bx - dx) * (ay - dy) std::array<RealNumber, 4> C_33 = two_two_expansion_diff(A_11_x_A_22, A_21_x_A_12); std::array<RealNumber, 8> C_33_x_A31; // ( (ax - dx) * (by - dy) - (bx - dx) * (ay - dy) ) * ( cx - dx ) int C_33_x_A31_nz = scale_expansion_zeroelim(C_33, A_31, C_33_x_A31); std::array<RealNumber, 16> C_33_x_A31_sq; // ( (ax - dx) * (by - dy) - (bx - dx) * (ay - dy) ) * ( cx - dx ) * ( cx - dx ) int C_33_x_A31_sq_nz = scale_expansion_zeroelim(C_33_x_A31, A_31, C_33_x_A31_sq, C_33_x_A31_nz); std::array<RealNumber, 8> C_33_x_A_32; // ( (ax - dx) * (by - dy) - (bx - dx) * (ay - dy) ) * ( cy - dy ) int C_33_x_A_32_nz = scale_expansion_zeroelim(C_33, A_32, C_33_x_A_32); std::array<RealNumber, 16> C_33_x_A_32_sq; // ( (ax - dx) * (by - dy) - (bx - dx) * (ay - dy) ) * ( cy - dy ) * ( cy - dy ) int C_33_x_A_32_sq_nz = scale_expansion_zeroelim(C_33_x_A_32, A_32, C_33_x_A_32_sq, C_33_x_A_32_nz); std::array<RealNumber, 32> A_33_x_C_33; int A_33_x_C_33_nz = fast_expansion_sum_zeroelim(C_33_x_A31_sq, C_33_x_A_32_sq, A_33_x_C_33, C_33_x_A31_sq_nz, C_33_x_A_32_sq_nz); std::array<RealNumber, 64> A_13_x_C13_p_A_13_x_C13; int A_13_x_C13_p_A_13_x_C13_nz = fast_expansion_sum_zeroelim( A_13_x_C13, A_23_x_C_23, A_13_x_C13_p_A_13_x_C13, A_13_x_C13_nz, A_23_x_C_23_nz); std::array<RealNumber, 96> det_expansion; int det_expansion_nz = fast_expansion_sum_zeroelim( A_13_x_C13_p_A_13_x_C13, A_33_x_C_33, det_expansion, A_13_x_C13_p_A_13_x_C13_nz, A_33_x_C_33_nz); det = std::accumulate(det_expansion.begin(), det_expansion.begin() + det_expansion_nz, static_cast<RealNumber>(0)); if (BOOST_GEOMETRY_CONDITION(Robustness == 1)) { return det; } RealNumber B_relative_bound = (2 + 12 * std::numeric_limits<RealNumber>::epsilon()) * std::numeric_limits<RealNumber>::epsilon(); absolute_bound = B_relative_bound * magnitude; if (std::abs(det) >= absolute_bound) { return det; } RealNumber A_11tail = two_diff_tail(p1[0], p4[0], A_11); RealNumber A_12tail = two_diff_tail(p1[1], p4[1], A_12); RealNumber A_21tail = two_diff_tail(p2[0], p4[0], A_21); RealNumber A_22tail = two_diff_tail(p2[1], p4[1], A_22); RealNumber A_31tail = two_diff_tail(p3[0], p4[0], A_31); RealNumber A_32tail = two_diff_tail(p3[1], p4[1], A_32); if ((A_11tail == 0) && (A_21tail == 0) && (A_31tail == 0) && (A_12tail == 0) && (A_22tail == 0) && (A_32tail == 0)) { return det; } // RealNumber sub_bound = (1.5 + 2.0 * std::numeric_limits<RealNumber>::epsilon()) // * std::numeric_limits<RealNumber>::epsilon(); // RealNumber C_relative_bound = (11.0 + 72.0 * std::numeric_limits<RealNumber>::epsilon()) // * std::numeric_limits<RealNumber>::epsilon() // * std::numeric_limits<RealNumber>::epsilon(); //absolute_bound = C_relative_bound * magnitude + sub_bound * std::abs(det); det += ((A_11 * A_11 + A_12 * A_12) * ((A_21 * A_32tail + A_32 * A_21tail) - (A_22 * A_31tail + A_31 * A_22tail)) + 2 * (A_11 * A_11tail + A_12 * A_12tail) * (A_21 * A_32 - A_22 * A_31)) + ((A_21 * A_21 + A_22 * A_22) * ((A_31 * A_12tail + A_12 * A_31tail) - (A_32 * A_11tail + A_11 * A_32tail)) + 2 * (A_21 * A_21tail + A_22 * A_22tail) * (A_31 * A_12 - A_32 * A_11)) + ((A_31 * A_31 + A_32 * A_32) * ((A_11 * A_22tail + A_22 * A_11tail) - (A_12 * A_21tail + A_21 * A_12tail)) + 2 * (A_31 * A_31tail + A_32 * A_32tail) * (A_11 * A_22 - A_12 * A_21)); //if (std::abs(det) >= absolute_bound) //{ return det; //} } }} // namespace detail::precise_math }} // namespace boost::geometry #endif // BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_CARTESIAN_DETAIL_PRECISE_MATH_HPP