geometry/util/precise_math.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2019 Tinko Bartels, Berlin, Germany.
0004 // Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
0005
0006 // Contributed and/or modified by Tinko Bartels,
0007 //   as part of Google Summer of Code 2019 program.
0008
0009 // This file was modified by Oracle on 2021.
0010 // Modifications copyright (c) 2021, Oracle and/or its affiliates.
0011 // Contributed and/or modified by Vissarion Fisikopoulos, on behalf of Oracle
0012
0013 // Use, modification and distribution is subject to the Boost Software License,
0014 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0015 // http://www.boost.org/LICENSE_1_0.txt)
0016
0017 #ifndef BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_CARTESIAN_DETAIL_PRECISE_MATH_HPP
0018 #define BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_CARTESIAN_DETAIL_PRECISE_MATH_HPP
0019
0020 #include<numeric>
0021 #include<cmath>
0022 #include<limits>
0023 #include<array>
0024
0025 #include <boost/geometry/core/access.hpp>
0026 #include <boost/geometry/util/condition.hpp>
0027
0028 // The following code is based on "Adaptive Precision Floating-Point Arithmetic
0029 // and Fast Robust Geometric Predicates" by Richard Shewchuk,
0030 // J. Discrete Comput Geom (1997) 18: 305. https://doi.org/10.1007/PL00009321
0031
0032 namespace boost { namespace geometry
0033 {
0034
0035 namespace detail { namespace precise_math
0036 {
0037
0038 // See Theorem 6, page 6
0039 template
0040 <
0041     typename RealNumber
0042 >
0043 inline std::array<RealNumber, 2> fast_two_sum(RealNumber const a,
0044                                               RealNumber const b)
0045 {
0046     RealNumber x = a + b;
0047     RealNumber b_virtual = x - a;
0048     return {{x, b - b_virtual}};
0049 }
0050
0051 // See Theorem 7, page 7 - 8
0052 template
0053 <
0054     typename RealNumber
0055 >
0056 inline std::array<RealNumber, 2> two_sum(RealNumber const a,
0057                                          RealNumber const b)
0058 {
0059     RealNumber x = a + b;
0060     RealNumber b_virtual = x - a;
0061     RealNumber a_virtual = x - b_virtual;
0062     RealNumber b_roundoff = b - b_virtual;
0063     RealNumber a_roundoff = a - a_virtual;
0064     RealNumber y = a_roundoff + b_roundoff;
0065     return {{ x,  y }};
0066 }
0067
0068 // See bottom of page 8
0069 template
0070 <
0071     typename RealNumber
0072 >
0073 inline RealNumber two_diff_tail(RealNumber const a,
0074                                 RealNumber const b,
0075                                 RealNumber const x)
0076 {
0077     RealNumber b_virtual = a - x;
0078     RealNumber a_virtual = x + b_virtual;
0079     RealNumber b_roundoff = b_virtual - b;
0080     RealNumber a_roundoff = a - a_virtual;
0081     return a_roundoff + b_roundoff;
0082 }
0083
0084 // see bottom of page 8
0085 template
0086 <
0087     typename RealNumber
0088 >
0089 inline std::array<RealNumber, 2> two_diff(RealNumber const a,
0090                                           RealNumber const b)
0091 {
0092     RealNumber x = a - b;
0093     RealNumber y = two_diff_tail(a, b, x);
0094     return {{ x, y }};
0095 }
0096
0097 // see theorem 18, page 19
0098 template
0099 <
0100     typename RealNumber
0101 >
0102 inline RealNumber two_product_tail(RealNumber const a,
0103                                    RealNumber const b,
0104                                    RealNumber const x)
0105 {
0106     return std::fma(a, b, -x);
0107 }
0108
0109 // see theorem 18, page 19
0110 template
0111 <
0112     typename RealNumber
0113 >
0114 inline std::array<RealNumber, 2> two_product(RealNumber const a,
0115                                              RealNumber const b)
0116 {
0117     RealNumber x = a * b;
0118     RealNumber y = two_product_tail(a, b, x);
0119     return {{ x , y }};
0120 }
0121
0122 // see theorem 12, figure 7, page 11 - 12,
0123 // this is the 2 by 2 case for the corresponding diff-method
0124 // note that this method takes input in descending order of magnitude and
0125 // returns components in ascending order of magnitude
0126 template
0127 <
0128     typename RealNumber
0129 >
0130 inline std::array<RealNumber, 4> two_two_expansion_diff(
0131         std::array<RealNumber, 2> const a,
0132         std::array<RealNumber, 2> const b)
0133 {
0134     std::array<RealNumber, 4> h;
0135     std::array<RealNumber, 2> Qh = two_diff(a[1], b[1]);
0136     h[0] = Qh[1];
0137     Qh = two_sum( a[0], Qh[0] );
0138     RealNumber _j = Qh[0];
0139     Qh = two_diff(Qh[1], b[0]);
0140     h[1] = Qh[1];
0141     Qh = two_sum( _j, Qh[0] );
0142     h[2] = Qh[1];
0143     h[3] = Qh[0];
0144     return h;
0145 }
0146
0147 // see theorem 13, figure 8. This implementation uses zero elimination as
0148 // suggested on page 17, second to last paragraph. Returns the number of
0149 // non-zero components in the result and writes the result to h.
0150 // the merger into a single sequence g is done implicitly
0151 template
0152 <
0153     typename RealNumber,
0154     std::size_t InSize1,
0155     std::size_t InSize2,
0156     std::size_t OutSize
0157 >
0158 inline int fast_expansion_sum_zeroelim(
0159         std::array<RealNumber, InSize1> const& e,
0160         std::array<RealNumber, InSize2> const& f,
0161         std::array<RealNumber, OutSize> & h,
0162         int m = InSize1,
0163         int n = InSize2)
0164 {
0165     std::array<RealNumber, 2> Qh;
0166     int i_e = 0;
0167     int i_f = 0;
0168     int i_h = 0;
0169     if (std::abs(f[0]) > std::abs(e[0]))
0170     {
0171         Qh[0] = e[i_e++];
0172     }
0173     else
0174     {
0175         Qh[0] = f[i_f++];
0176     }
0177     i_h = 0;
0178     if ((i_e < m) && (i_f < n))
0179     {
0180         if (std::abs(f[i_f]) > std::abs(e[i_e]))
0181         {
0182             Qh = fast_two_sum(e[i_e++], Qh[0]);
0183         }
0184         else
0185         {
0186             Qh = fast_two_sum(f[i_f++], Qh[0]);
0187         }
0188         if (Qh[1] != 0.0)
0189         {
0190             h[i_h++] = Qh[1];
0191         }
0192         while ((i_e < m) && (i_f < n))
0193         {
0194             if (std::abs(f[i_f]) > std::abs(e[i_e]))
0195             {
0196                 Qh = two_sum(Qh[0], e[i_e++]);
0197             }
0198             else
0199             {
0200                 Qh = two_sum(Qh[0], f[i_f++]);
0201             }
0202             if (Qh[1] != 0.0)
0203             {
0204                 h[i_h++] = Qh[1];
0205             }
0206         }
0207     }
0208     while (i_e < m)
0209     {
0210         Qh = two_sum(Qh[0], e[i_e++]);
0211         if (Qh[1] != 0.0)
0212         {
0213             h[i_h++] = Qh[1];
0214         }
0215     }
0216     while (i_f < n)
0217     {
0218         Qh = two_sum(Qh[0], f[i_f++]);
0219         if (Qh[1] != 0.0)
0220         {
0221             h[i_h++] = Qh[1];
0222         }
0223     }
0224     if ((Qh[0] != 0.0) || (i_h == 0))
0225     {
0226         h[i_h++] = Qh[0];
0227     }
0228     return i_h;
0229 }
0230
0231 // see theorem 19, figure 13, page 20 - 21. This implementation uses zero
0232 // elimination as suggested on page 17, second to last paragraph. Returns the
0233 // number of non-zero components in the result and writes the result to h.
0234 template
0235 <
0236     typename RealNumber,
0237     std::size_t InSize
0238 >
0239 inline int scale_expansion_zeroelim(
0240         std::array<RealNumber, InSize> const& e,
0241         RealNumber const b,
0242         std::array<RealNumber, 2 * InSize> & h,
0243         int e_non_zeros = InSize)
0244 {
0245     std::array<RealNumber, 2> Qh = two_product(e[0], b);
0246     int i_h = 0;
0247     if (Qh[1] != 0)
0248     {
0249         h[i_h++] = Qh[1];
0250     }
0251     for (int i_e = 1; i_e < e_non_zeros; i_e++)
0252     {
0253         std::array<RealNumber, 2> Tt = two_product(e[i_e], b);
0254         Qh = two_sum(Qh[0], Tt[1]);
0255         if (Qh[1] != 0)
0256         {
0257             h[i_h++] = Qh[1];
0258         }
0259         Qh = fast_two_sum(Tt[0], Qh[0]);
0260         if (Qh[1] != 0)
0261         {
0262             h[i_h++] = Qh[1];
0263         }
0264     }
0265     if ((Qh[0] != 0.0) || (i_h == 0))
0266     {
0267         h[i_h++] = Qh[0];
0268     }
0269     return i_h;
0270 }
0271
0272 template<typename RealNumber>
0273 struct vec2d
0274 {
0275     RealNumber x;
0276     RealNumber y;
0277 };
0278
0279 template
0280 <
0281     typename RealNumber,
0282     std::size_t Robustness
0283 >
0284 inline RealNumber orient2dtail(vec2d<RealNumber> const& p1,
0285                                vec2d<RealNumber> const& p2,
0286                                vec2d<RealNumber> const& p3,
0287                                std::array<RealNumber, 2>& t1,
0288                                std::array<RealNumber, 2>& t2,
0289                                std::array<RealNumber, 2>& t3,
0290                                std::array<RealNumber, 2>& t4,
0291                                std::array<RealNumber, 2>& t5_01,
0292                                std::array<RealNumber, 2>& t6_01,
0293                                RealNumber const& magnitude)
0294 {
0295     t5_01[1] = two_product_tail(t1[0], t2[0], t5_01[0]);
0296     t6_01[1] = two_product_tail(t3[0], t4[0], t6_01[0]);
0297     std::array<RealNumber, 4> tA_03 = two_two_expansion_diff(t5_01, t6_01);
0298     RealNumber det = std::accumulate(tA_03.begin(), tA_03.end(), static_cast<RealNumber>(0));
0299     if (BOOST_GEOMETRY_CONDITION(Robustness == 1))
0300     {
0301         return det;
0302     }
0303     // see p.39, mind the different definition of epsilon for error bound
0304     RealNumber B_relative_bound =
0305           (1 + 3 * std::numeric_limits<RealNumber>::epsilon())
0306         * std::numeric_limits<RealNumber>::epsilon();
0307     RealNumber absolute_bound = B_relative_bound * magnitude;
0308     if (std::abs(det) >= absolute_bound)
0309     {
0310         return det; //B estimate
0311     }
0312     t1[1] = two_diff_tail(p1.x, p3.x, t1[0]);
0313     t2[1] = two_diff_tail(p2.y, p3.y, t2[0]);
0314     t3[1] = two_diff_tail(p1.y, p3.y, t3[0]);
0315     t4[1] = two_diff_tail(p2.x, p3.x, t4[0]);
0316
0317     if ((t1[1] == 0) && (t3[1] == 0) && (t2[1] == 0) && (t4[1] == 0))
0318     {
0319         return det; //If all tails are zero, there is noething else to compute
0320     }
0321     RealNumber sub_bound =
0322           (1.5 + 2 * std::numeric_limits<RealNumber>::epsilon())
0323         * std::numeric_limits<RealNumber>::epsilon();
0324     // see p.39, mind the different definition of epsilon for error bound
0325     RealNumber C_relative_bound =
0326           (2.25 + 8 * std::numeric_limits<RealNumber>::epsilon())
0327         * std::numeric_limits<RealNumber>::epsilon()
0328         * std::numeric_limits<RealNumber>::epsilon();
0329     absolute_bound = C_relative_bound * magnitude + sub_bound * std::abs(det);
0330     det += (t1[0] * t2[1] + t2[0] * t1[1]) - (t3[0] * t4[1] + t4[0] * t3[1]);
0331     if (Robustness == 2 || std::abs(det) >= absolute_bound)
0332     {
0333         return det; //C estimate
0334     }
0335     std::array<RealNumber, 8> D_left;
0336     int D_left_nz;
0337     {
0338         std::array<RealNumber, 2> t5_23 = two_product(t1[1], t2[0]);
0339         std::array<RealNumber, 2> t6_23 = two_product(t3[1], t4[0]);
0340         std::array<RealNumber, 4> tA_47 = two_two_expansion_diff(t5_23, t6_23);
0341         D_left_nz = fast_expansion_sum_zeroelim(tA_03, tA_47, D_left);
0342     }
0343     std::array<RealNumber, 8> D_right;
0344     int D_right_nz;
0345     {
0346         std::array<RealNumber, 2> t5_45 = two_product(t1[0], t2[1]);
0347         std::array<RealNumber, 2> t6_45 = two_product(t3[0], t4[1]);
0348         std::array<RealNumber, 4> tA_8_11 = two_two_expansion_diff(t5_45, t6_45);
0349         std::array<RealNumber, 2> t5_67 = two_product(t1[1], t2[1]);
0350         std::array<RealNumber, 2> t6_67 = two_product(t3[1], t4[1]);
0351         std::array<RealNumber, 4> tA_12_15 = two_two_expansion_diff(t5_67, t6_67);
0352         D_right_nz = fast_expansion_sum_zeroelim(tA_8_11, tA_12_15, D_right);
0353     }
0354     std::array<RealNumber, 16> D;
0355     int D_nz = fast_expansion_sum_zeroelim(D_left, D_right, D, D_left_nz, D_right_nz);
0356     // only return component of highest magnitude because we mostly care about the sign.
0357     return(D[D_nz - 1]);
0358 }
0359
0360 // see page 38, Figure 21 for the calculations, notation follows the notation
0361 // in the figure.
0362 template
0363 <
0364     typename RealNumber,
0365     std::size_t Robustness = 3,
0366     typename EpsPolicy
0367 >
0368 inline RealNumber orient2d(vec2d<RealNumber> const& p1,
0369                            vec2d<RealNumber> const& p2,
0370                            vec2d<RealNumber> const& p3,
0371                            EpsPolicy& eps_policy)
0372 {
0373     std::array<RealNumber, 2> t1, t2, t3, t4;
0374     t1[0] = p1.x - p3.x;
0375     t2[0] = p2.y - p3.y;
0376     t3[0] = p1.y - p3.y;
0377     t4[0] = p2.x - p3.x;
0378
0379     eps_policy = EpsPolicy(t1[0], t2[0], t3[0], t4[0]);
0380
0381     std::array<RealNumber, 2> t5_01, t6_01;
0382     t5_01[0] = t1[0] * t2[0];
0383     t6_01[0] = t3[0] * t4[0];
0384     RealNumber det = t5_01[0] - t6_01[0];
0385
0386     if (BOOST_GEOMETRY_CONDITION(Robustness == 0))
0387     {
0388         return det;
0389     }
0390
0391     RealNumber const magnitude = std::abs(t5_01[0]) + std::abs(t6_01[0]);
0392
0393     // see p.39, mind the different definition of epsilon for error bound
0394     RealNumber const A_relative_bound =
0395           (1.5 + 4 * std::numeric_limits<RealNumber>::epsilon())
0396         * std::numeric_limits<RealNumber>::epsilon();
0397     RealNumber absolute_bound = A_relative_bound * magnitude;
0398     if ( std::abs(det) >= absolute_bound )
0399     {
0400         return det; //A estimate
0401     }
0402
0403     if ( (t5_01[0] > 0 && t6_01[0] <= 0) || (t5_01[0] < 0 && t6_01[0] >= 0) )
0404     {
0405         //if diagonal and antidiagonal have different sign, the sign of det is
0406         //obvious
0407         return det;
0408     }
0409     return orient2dtail<RealNumber, Robustness>(p1, p2, p3, t1, t2, t3, t4,
0410                                                 t5_01, t6_01, magnitude);
0411 }
0412
0413 // This method adaptively computes increasingly precise approximations of the following
0414 // determinant using Laplace expansion along the last column.
0415 // det A =
0416 //      | p1_x - p4_x    p1_y - p4_y     ( p1_x - p4_x ) ^ 2 + ( p1_y - p4_y ) ^ 2 |
0417 //      | p2_x - p4_x    p2_y - p4_y     ( p2_x - p4_x ) ^ 2 + ( p1_y - p4_y ) ^ 2 |
0418 //      | p3_x - p4_x    p3_y - p4_y     ( p3_x - p4_x ) ^ 2 + ( p3_y - p4_y ) ^ 2 |
0419 // = a_13 * C_13 + a_23 * C_23 + a_33 * C_33
0420 // where a_ij is the i-j-entry and C_ij is the i_j Cofactor
0421
0422 template
0423 <
0424     typename RealNumber,
0425     std::size_t Robustness = 2
0426 >
0427 RealNumber incircle(std::array<RealNumber, 2> const& p1,
0428                     std::array<RealNumber, 2> const& p2,
0429                     std::array<RealNumber, 2> const& p3,
0430                     std::array<RealNumber, 2> const& p4)
0431 {
0432     RealNumber A_11 = p1[0] - p4[0];
0433     RealNumber A_21 = p2[0] - p4[0];
0434     RealNumber A_31 = p3[0] - p4[0];
0435     RealNumber A_12 = p1[1] - p4[1];
0436     RealNumber A_22 = p2[1] - p4[1];
0437     RealNumber A_32 = p3[1] - p4[1];
0438
0439     std::array<RealNumber, 2> A_21_x_A_32,
0440                               A_31_x_A_22,
0441                               A_31_x_A_12,
0442                               A_11_x_A_32,
0443                               A_11_x_A_22,
0444                               A_21_x_A_12;
0445     A_21_x_A_32[0] = A_21 * A_32;
0446     A_31_x_A_22[0] = A_31 * A_22;
0447     RealNumber A_13 = A_11 * A_11 + A_12 * A_12;
0448
0449     A_31_x_A_12[0] = A_31 * A_12;
0450     A_11_x_A_32[0] = A_11 * A_32;
0451     RealNumber A_23 = A_21 * A_21 + A_22 * A_22;
0452
0453     A_11_x_A_22[0] = A_11 * A_22;
0454     A_21_x_A_12[0] = A_21 * A_12;
0455     RealNumber A_33 = A_31 * A_31 + A_32 * A_32;
0456
0457     RealNumber det = A_13 * (A_21_x_A_32[0] - A_31_x_A_22[0])
0458       + A_23 * (A_31_x_A_12[0] - A_11_x_A_32[0])
0459       + A_33 * (A_11_x_A_22[0] - A_21_x_A_12[0]);
0460     if(Robustness == 0) return det;
0461
0462     RealNumber magnitude =
0463           (std::abs(A_21_x_A_32[0]) + std::abs(A_31_x_A_22[0])) * A_13
0464         + (std::abs(A_31_x_A_12[0]) + std::abs(A_11_x_A_32[0])) * A_23
0465         + (std::abs(A_11_x_A_22[0]) + std::abs(A_21_x_A_12[0])) * A_33;
0466     RealNumber A_relative_bound =
0467           (5 + 24 * std::numeric_limits<RealNumber>::epsilon())
0468         * std::numeric_limits<RealNumber>::epsilon();
0469     RealNumber absolute_bound = A_relative_bound * magnitude;
0470     if (std::abs(det) > absolute_bound)
0471     {
0472         return det;
0473     }
0474     // (p2_x - p4_x) * (p3_y - p4_y)
0475     A_21_x_A_32[1] = two_product_tail(A_21, A_32, A_21_x_A_32[0]);
0476     // (p3_x - p4_x) * (p2_y - p4_y)
0477     A_31_x_A_22[1] = two_product_tail(A_31, A_22, A_31_x_A_22[0]);
0478     // (bx - dx) * (cy - dy) - (cx - dx) * (by - dy)
0479     std::array<RealNumber, 4> C_13 = two_two_expansion_diff(A_21_x_A_32, A_31_x_A_22);
0480     std::array<RealNumber, 8> C_13_x_A11;
0481     // ( (bx - dx) * (cy - dy) - (cx - dx) * (by - dy) ) * ( ax - dx )
0482     int C_13_x_A11_nz = scale_expansion_zeroelim(C_13, A_11, C_13_x_A11);
0483     std::array<RealNumber, 16> C_13_x_A11_sq;
0484     // ( (bx - dx) * (cy - dy) - (cx - dx) * (by - dy) ) * ( ax - dx ) * (ax - dx)
0485     int C_13_x_A11_sq_nz = scale_expansion_zeroelim(C_13_x_A11,
0486                                                     A_11,
0487                                                     C_13_x_A11_sq,
0488                                                     C_13_x_A11_nz);
0489
0490     std::array<RealNumber, 8> C_13_x_A12;
0491     // ( (bx - dx) * (cy - dy) - (cx - dx) * (by - dy) ) * ( ay - dy )
0492     int C_13_x_A12_nz = scale_expansion_zeroelim(C_13, A_12, C_13_x_A12);
0493
0494     std::array<RealNumber, 16> C_13_x_A12_sq;
0495     // ( (bx - dx) * (cy - dy) - (cx - dx) * (by - dy) ) * ( ay - dy ) * ( ay - dy )
0496     int C_13_x_A12_sq_nz = scale_expansion_zeroelim(C_13_x_A12, A_12,
0497                                                     C_13_x_A12_sq,
0498                                                     C_13_x_A12_nz);
0499
0500     std::array<RealNumber, 32> A_13_x_C13;
0501     //   ( (bx - dx) * (cy - dy) - (cx - dx) * (by - dy) )
0502     // * ( ( ay - dy ) * ( ay - dy ) + ( ax - dx ) * (ax - dx) )
0503     int A_13_x_C13_nz = fast_expansion_sum_zeroelim(C_13_x_A11_sq,
0504                                                     C_13_x_A12_sq,
0505                                                     A_13_x_C13,
0506                                                     C_13_x_A11_sq_nz,
0507                                                     C_13_x_A12_sq_nz);
0508
0509     // (cx - dx) * (ay - dy)
0510     A_31_x_A_12[1] = two_product_tail(A_31, A_12, A_31_x_A_12[0]);
0511     // (ax - dx) * (cy - dy)
0512     A_11_x_A_32[1] = two_product_tail(A_11, A_32, A_11_x_A_32[0]);
0513     // (cx - dx) * (ay - dy) - (ax - dx) * (cy - dy)
0514     std::array<RealNumber, 4> C_23 = two_two_expansion_diff(A_31_x_A_12,
0515                                                             A_11_x_A_32);
0516     std::array<RealNumber, 8> C_23_x_A_21;
0517     // ( (cx - dx) * (ay - dy) - (ax - dx) * (cy - dy) ) * ( bx - dx )
0518     int C_23_x_A_21_nz = scale_expansion_zeroelim(C_23, A_21, C_23_x_A_21);
0519     std::array<RealNumber, 16> C_23_x_A_21_sq;
0520     // ( (cx - dx) * (ay - dy) - (ax - dx) * (cy - dy) ) * ( bx - dx ) * ( bx - dx )
0521     int C_23_x_A_21_sq_nz = scale_expansion_zeroelim(C_23_x_A_21, A_21,
0522                                                      C_23_x_A_21_sq,
0523                                                      C_23_x_A_21_nz);
0524     std::array<RealNumber, 8>  C_23_x_A_22;
0525     // ( (cx - dx) * (ay - dy) - (ax - dx) * (cy - dy) ) * ( by - dy )
0526     int C_23_x_A_22_nz = scale_expansion_zeroelim(C_23, A_22, C_23_x_A_22);
0527     std::array<RealNumber, 16> C_23_x_A_22_sq;
0528     // ( (cx - dx) * (ay - dy) - (ax - dx) * (cy - dy) ) * ( by - dy ) * ( by - dy )
0529     int C_23_x_A_22_sq_nz = scale_expansion_zeroelim(C_23_x_A_22, A_22,
0530                                                      C_23_x_A_22_sq,
0531                                                      C_23_x_A_22_nz);
0532     std::array<RealNumber, 32> A_23_x_C_23;
0533     //   ( (cx - dx) * (ay - dy) - (ax - dx) * (cy - dy) )
0534     // * ( ( bx - dx ) * ( bx - dx ) + ( by - dy ) * ( by - dy ) )
0535     int A_23_x_C_23_nz = fast_expansion_sum_zeroelim(C_23_x_A_21_sq,
0536                                                      C_23_x_A_22_sq,
0537                                                      A_23_x_C_23,
0538                                                      C_23_x_A_21_sq_nz,
0539                                                      C_23_x_A_22_sq_nz);
0540
0541     // (ax - dx) * (by - dy)
0542     A_11_x_A_22[1] = two_product_tail(A_11, A_22, A_11_x_A_22[0]);
0543     // (bx - dx) * (ay - dy)
0544     A_21_x_A_12[1] = two_product_tail(A_21, A_12, A_21_x_A_12[0]);
0545     // (ax - dx) * (by - dy) - (bx - dx) * (ay - dy)
0546     std::array<RealNumber, 4> C_33 = two_two_expansion_diff(A_11_x_A_22,
0547                                                             A_21_x_A_12);
0548     std::array<RealNumber, 8>  C_33_x_A31;
0549     // ( (ax - dx) * (by - dy) - (bx - dx) * (ay - dy) ) * ( cx - dx )
0550     int C_33_x_A31_nz = scale_expansion_zeroelim(C_33, A_31, C_33_x_A31);
0551     std::array<RealNumber, 16> C_33_x_A31_sq;
0552     // ( (ax - dx) * (by - dy) - (bx - dx) * (ay - dy) ) * ( cx - dx ) * ( cx - dx )
0553     int C_33_x_A31_sq_nz = scale_expansion_zeroelim(C_33_x_A31, A_31,
0554                                                     C_33_x_A31_sq,
0555                                                     C_33_x_A31_nz);
0556     std::array<RealNumber, 8>  C_33_x_A_32;
0557     // ( (ax - dx) * (by - dy) - (bx - dx) * (ay - dy) ) * ( cy - dy )
0558     int C_33_x_A_32_nz = scale_expansion_zeroelim(C_33, A_32, C_33_x_A_32);
0559     std::array<RealNumber, 16> C_33_x_A_32_sq;
0560     // ( (ax - dx) * (by - dy) - (bx - dx) * (ay - dy) ) * ( cy - dy ) * ( cy - dy )
0561     int C_33_x_A_32_sq_nz = scale_expansion_zeroelim(C_33_x_A_32, A_32,
0562                                                      C_33_x_A_32_sq,
0563                                                      C_33_x_A_32_nz);
0564
0565     std::array<RealNumber, 32> A_33_x_C_33;
0566     int A_33_x_C_33_nz = fast_expansion_sum_zeroelim(C_33_x_A31_sq,
0567                                                      C_33_x_A_32_sq,
0568                                                      A_33_x_C_33,
0569                                                      C_33_x_A31_sq_nz,
0570                                                      C_33_x_A_32_sq_nz);
0571     std::array<RealNumber, 64> A_13_x_C13_p_A_13_x_C13;
0572     int A_13_x_C13_p_A_13_x_C13_nz = fast_expansion_sum_zeroelim(
0573             A_13_x_C13, A_23_x_C_23,
0574             A_13_x_C13_p_A_13_x_C13,
0575             A_13_x_C13_nz,
0576             A_23_x_C_23_nz);
0577     std::array<RealNumber, 96> det_expansion;
0578     int det_expansion_nz = fast_expansion_sum_zeroelim(
0579             A_13_x_C13_p_A_13_x_C13,
0580             A_33_x_C_33,
0581             det_expansion,
0582             A_13_x_C13_p_A_13_x_C13_nz,
0583             A_33_x_C_33_nz);
0584
0585     det = std::accumulate(det_expansion.begin(),
0586                           det_expansion.begin() + det_expansion_nz,
0587                           static_cast<RealNumber>(0));
0588     if(Robustness == 1) return det;
0589     RealNumber B_relative_bound =
0590           (2 + 12 * std::numeric_limits<RealNumber>::epsilon())
0591         * std::numeric_limits<RealNumber>::epsilon();
0592     absolute_bound = B_relative_bound * magnitude;
0593     if (std::abs(det) >= absolute_bound)
0594     {
0595         return det;
0596     }
0597     RealNumber A_11tail = two_diff_tail(p1[0], p4[0], A_11);
0598     RealNumber A_12tail = two_diff_tail(p1[1], p4[1], A_12);
0599     RealNumber A_21tail = two_diff_tail(p2[0], p4[0], A_21);
0600     RealNumber A_22tail = two_diff_tail(p2[1], p4[1], A_22);
0601     RealNumber A_31tail = two_diff_tail(p3[0], p4[0], A_31);
0602     RealNumber A_32tail = two_diff_tail(p3[1], p4[1], A_32);
0603     if ((A_11tail == 0) && (A_21tail == 0) && (A_31tail == 0)
0604         && (A_12tail == 0) && (A_22tail == 0) && (A_32tail == 0))
0605     {
0606         return det;
0607     }
0608     //  RealNumber sub_bound =  (1.5 + 2.0 * std::numeric_limits<RealNumber>::epsilon())
0609     //    * std::numeric_limits<RealNumber>::epsilon();
0610     //  RealNumber C_relative_bound = (11.0 + 72.0 * std::numeric_limits<RealNumber>::epsilon())
0611     //    * std::numeric_limits<RealNumber>::epsilon()
0612     //    * std::numeric_limits<RealNumber>::epsilon();
0613     //absolute_bound = C_relative_bound * magnitude + sub_bound * std::abs(det);
0614     det += ((A_11 * A_11 + A_12 * A_12) * ((A_21 * A_32tail + A_32 * A_21tail)
0615         - (A_22 * A_31tail + A_31 * A_22tail))
0616     + 2 * (A_11 * A_11tail + A_12 * A_12tail) * (A_21 * A_32 - A_22 * A_31))
0617     + ((A_21 * A_21 + A_22 * A_22) * ((A_31 * A_12tail + A_12 * A_31tail)
0618         - (A_32 * A_11tail + A_11 * A_32tail))
0619     + 2 * (A_21 * A_21tail + A_22 * A_22tail) * (A_31 * A_12 - A_32 * A_11))
0620     + ((A_31 * A_31 + A_32 * A_32) * ((A_11 * A_22tail + A_22 * A_11tail)
0621         - (A_12 * A_21tail + A_21 * A_12tail))
0622     + 2 * (A_31 * A_31tail + A_32 * A_32tail) * (A_11 * A_22 - A_12 * A_21));
0623     //if (std::abs(det) >= absolute_bound)
0624     //{
0625     return det;
0626     //}
0627 }
0628
0629 }} // namespace detail::precise_math
0630
0631 }} // namespace boost::geometry
0632
0633 #endif // BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_CARTESIAN_DETAIL_PRECISE_MATH_HPP