GCC Code Coverage Report

Directory: ./
File: lib/geogram/numerics/predicates.cpp
Date: 2026-09-07 02:36:43
Exec Total Coverage
Lines: 507 973 52.1%
Functions: 51 74 68.9%
Branches: 501 2908 17.2%
Line Branch Exec Source
1 /*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40 #include <geogram/basic/common.h>
41 #include <geogram/basic/numeric.h>
42
43 // This makes sure the compiler will not optimize y = a*x+b
44 // with fused multiply-add, this would break the exact
45 // predicates.
46 GEO_FP_CONTRACT_OFF
47
48 #include <geogram/numerics/predicates.h>
49 #include <geogram/numerics/multi_precision.h>
50 #include <geogram/basic/assert.h>
51 #include <geogram/basic/logger.h>
52 #include <geogram/basic/command_line.h>
53 #include <geogram/basic/matrix.h>
54 #include <algorithm>
55
56 #define FPG_UNCERTAIN_VALUE 0
57
58 #include <geogram/numerics/predicates/side1.h>
59 #include <geogram/numerics/predicates/side2.h>
60 #include <geogram/numerics/predicates/side3.h>
61 #include <geogram/numerics/predicates/side3h.h>
62 #include <geogram/numerics/predicates/side3_2dlifted.h>
63 #include <geogram/numerics/predicates/side4.h>
64 #include <geogram/numerics/predicates/side4h.h>
65 #include <geogram/numerics/predicates/orient2d.h>
66 #include <geogram/numerics/predicates/orient3d.h>
67 #include <geogram/numerics/predicates/det3d.h>
68 #include <geogram/numerics/predicates/det4d.h>
69 #include <geogram/numerics/predicates/dot3d.h>
70 #include <geogram/numerics/predicates/dot_compare_3d.h>
71 #include <geogram/numerics/predicates/det_compare_4d.h>
72 #include <geogram/numerics/predicates/aligned3d.h>
73
74 #ifdef __SSE2__
75 #include <emmintrin.h>
76 #endif
77
78
79 namespace {
80
81 using namespace GEO;
82
83 GEO::PCK::SOSMode SOS_mode_ = GEO::PCK::SOS_ADDRESS;
84
85 /**
86 * \brief Comparator class for nD points using lexicographic order.
87 * \details Used by symbolic perturbations.
88 */
89 class LexicoCompare {
90 public:
91
92 /**
93 * \brief LexicoCompare constructor.
94 * \param[in] dim dimension of the points to compare.
95 */
96 ✗ LexicoCompare(index_t dim) : dim_(dim) {
97 ✗ }
98
99 /**
100 * \brief Compares two points with respect to the lexicographic
101 * order.
102 * \param[in] x , y pointers to the coordinates of the two points.
103 * \retval true if x is strictly before y in the lexicographic order.
104 * \retval false otherwise.
105 */
106 ✗ bool operator()(const double* x, const double* y) const {
107 ✗ for(index_t i=0; i<dim_-1; ++i) {
108 ✗ if(x[i] < y[i]) {
109 ✗ return true;
110 }
111 ✗ if(x[i] > y[i]) {
112 ✗ return false;
113 }
114 }
115 ✗ return (x[dim_-1] < y[dim_-1]);
116 }
117 private:
118 index_t dim_;
119 };
120
121 /**
122 * \brief Compares two 3D points with respect to the lexicographic
123 * order.
124 * \param[in] x , y pointers to the coordinates of the two 3D points.
125 * \retval true if x is strictly before y in the lexicographic order.
126 * \retval false otherwise.
127 */
128 3670430 bool lexico_compare_3d(const double* x, const double* y) {
129
2/2
✓ Branch 0 taken 830578 times.
✓ Branch 1 taken 2839852 times.
 3670430 if(x[0] < y[0]) {
130 830578 return true;
131 }
132
2/2
✓ Branch 0 taken 1339714 times.
✓ Branch 1 taken 1500138 times.
 2839852 if(x[0] > y[0]) {
133 1339714 return false;
134 }
135
2/2
✓ Branch 0 taken 374518 times.
✓ Branch 1 taken 1125620 times.
 1500138 if(x[1] < y[1]) {
136 374518 return true;
137 }
138
2/2
✓ Branch 0 taken 633604 times.
✓ Branch 1 taken 492016 times.
 1125620 if(x[1] > y[1]) {
139 633604 return false;
140 }
141 492016 return x[2] < y[2];
142 }
143
144 /**
145 * \brief Sorts an array of pointers to points.
146 * \details set_SOS_mode() alters the behavior of this function.
147 * If set to PCK::SOS_ADDRESS, then just the addresses of the points
148 * are sorted. If set to PCK::SOS_LEXICO, then the points are sorted
149 * in function of the lexicographic order of their coordinates.
150 * \param[in] begin a pointer to the first point.
151 * \param[in] end one position past the pointer to the last point.
152 * \param[in] dim the dimension of the points.
153 */
154 994465 void GEOGRAM_API SOS_sort(
155 const double** begin, const double** end, index_t dim
156 ) {
157
2/2
✓ Branch 0 taken 439217 times.
✓ Branch 1 taken 555248 times.
 994465 if(SOS_mode_ == PCK::SOS_ADDRESS) {
158 439217 std::sort(begin, end);
159 } else {
160
1/2
✓ Branch 0 taken 555248 times.
✗ Branch 1 not taken.
 555248 if(dim == 3) {
161 555248 std::sort(begin, end, lexico_compare_3d);
162 } else {
163 ✗ std::sort(begin, end, LexicoCompare(dim));
164 }
165 }
166 994465 }
167
168
169 /**
170 * \brief Gets the maximum of 4 double precision numbers.
171 * \param[in] x1 , x2 , x3 , x4 the four numbers.
172 * \return the maximum.
173 */
174 1893975 inline double max4(double x1, double x2, double x3, double x4) {
175 #ifdef __SSE2__
176 double result;
177 1893975 __m128d X1 =_mm_load_sd(&x1);
178 1893975 __m128d X2 =_mm_load_sd(&x2);
179 1893975 __m128d X3 =_mm_load_sd(&x3);
180 1893975 __m128d X4 =_mm_load_sd(&x4);
181 1893975 X1 = _mm_max_sd(X1,X2);
182 1893975 X3 = _mm_max_sd(X3,X4);
183 1893975 X1 = _mm_max_sd(X1,X3);
184 _mm_store_sd(&result, X1);
185 1893975 return result;
186 #else
187 return std::max(std::max(x1,x2),std::max(x3,x4));
188 #endif
189 }
190
191
192 /**
193 * \brief Gets the minimum and maximum of 3 double precision numbers.
194 * \param[in] x1 , x2 , x3 the three numbers.
195 * \param[out] m the minimum
196 * \param[out] M the maximum
197 */
198 631325 inline void get_minmax3(
199 double& m, double& M, double x1, double x2, double x3
200 ) {
201 #ifdef __SSE2__
202 631325 __m128d X1 =_mm_load_sd(&x1);
203 631325 __m128d X2 =_mm_load_sd(&x2);
204 631325 __m128d X3 =_mm_load_sd(&x3);
205 631325 __m128d MIN12 = _mm_min_sd(X1,X2);
206 631325 __m128d MAX12 = _mm_max_sd(X1,X2);
207 631325 X1 = _mm_min_sd(MIN12, X3);
208 631325 X3 = _mm_max_sd(MAX12, X3);
209 _mm_store_sd(&m, X1);
210 _mm_store_sd(&M, X3);
211 #else
212 m = std::min(std::min(x1,x2), x3);
213 M = std::max(std::max(x1,x2), x3);
214 #endif
215 631325 }
216
217 /**
218 * \brief Arithmetic filter for the in_sphere_3d_SOS() predicate.
219 * \details This filter was optimized by hand by Sylvain Pion
220 * (may be faster than FPG/PCK-generated filter).
221 * Since it is used massively by Delaunay_3d, using the
222 * optimized version may be worth it.
223 * \param[in] p first vertex of the tetrahedron
224 * \param[in] q second vertex of the tetrahedron
225 * \param[in] r third vertex of the tetrahedron
226 * \param[in] s fourth vertex of the tetrahedron
227 * \param[in] t point to be tested
228 * \retval +1 if \p t was determined to be outside
229 * the circumsphere of \p p,\p q,\p r,\p s
230 * \retval -1 if \p t was determined to be inside
231 * the circumsphere of \p p,\p q,\p r,\p s
232 * \retval 0 if the position of \p t could be be determined
233 */
234 631325 inline int in_sphere_3d_filter_optim(
235 const double* p, const double* q,
236 const double* r, const double* s, const double* t
237 ) {
238 631325 double ptx = p[0] - t[0];
239 631325 double pty = p[1] - t[1];
240 631325 double ptz = p[2] - t[2];
241 631325 double pt2 = geo_sqr(ptx) + geo_sqr(pty) + geo_sqr(ptz);
242
243 631325 double qtx = q[0] - t[0];
244 631325 double qty = q[1] - t[1];
245 631325 double qtz = q[2] - t[2];
246 631325 double qt2 = geo_sqr(qtx) + geo_sqr(qty) + geo_sqr(qtz);
247
248 631325 double rtx = r[0] - t[0];
249 631325 double rty = r[1] - t[1];
250 631325 double rtz = r[2] - t[2];
251 631325 double rt2 = geo_sqr(rtx) + geo_sqr(rty) + geo_sqr(rtz);
252
253 631325 double stx = s[0] - t[0];
254 631325 double sty = s[1] - t[1];
255 631325 double stz = s[2] - t[2];
256 631325 double st2 = geo_sqr(stx) + geo_sqr(sty) + geo_sqr(stz);
257
258 // Compute the semi-static bound.
259 631325 double maxx = ::fabs(ptx);
260 631325 double maxy = ::fabs(pty);
261 631325 double maxz = ::fabs(ptz);
262
263 631325 double aqtx = ::fabs(qtx);
264 631325 double artx = ::fabs(rtx);
265 631325 double astx = ::fabs(stx);
266
267 631325 double aqty = ::fabs(qty);
268 631325 double arty = ::fabs(rty);
269 631325 double asty = ::fabs(sty);
270
271 631325 double aqtz = ::fabs(qtz);
272 631325 double artz = ::fabs(rtz);
273 631325 double astz = ::fabs(stz);
274
275 631325 maxx = max4(maxx, aqtx, artx, astx);
276 631325 maxy = max4(maxy, aqty, arty, asty);
277 631325 maxz = max4(maxz, aqtz, artz, astz);
278
279 631325 double eps = 1.2466136531027298e-13 * maxx * maxy * maxz;
280
281 double min_max;
282 double max_max;
283 631325 get_minmax3(min_max, max_max, maxx, maxy, maxz);
284
285 631325 double det = det4x4(
286 ptx,pty,ptz,pt2,
287 rtx,rty,rtz,rt2,
288 qtx,qty,qtz,qt2,
289 stx,sty,stz,st2
290 );
291
292
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 631325 times.
 631325 if (min_max < 1e-58) { /* sqrt^5(min_double/eps) */
293 // Protect against underflow in the computation of eps.
294 ✗ return FPG_UNCERTAIN_VALUE;
295
1/2
✓ Branch 0 taken 631325 times.
✗ Branch 1 not taken.
 631325 } else if (max_max < 1e61) { /* sqrt^5(max_double/4 [hadamard]) */
296 // Protect against overflow in the computation of det.
297 631325 eps *= (max_max * max_max);
298 // Note: inverted as compared to CGAL
299 // CGAL: in_sphere_3d (called side_of_oriented_sphere())
300 // positive side is outside the sphere.
301 // PCK: in_sphere_3d : positive side is inside the sphere
302
2/2
✓ Branch 0 taken 97551 times.
✓ Branch 1 taken 533774 times.
 631325 if (det > eps) return -1;
303
2/2
✓ Branch 0 taken 379737 times.
✓ Branch 1 taken 154037 times.
 533774 if (det < -eps) return 1;
304 }
305
306 154037 return FPG_UNCERTAIN_VALUE;
307 }
308
309 using namespace GEO;
310
311 PCK::PredicateStats stats_side1("side1");
312 PCK::PredicateStats stats_side2("side2");
313 PCK::PredicateStats stats_side3("side3");
314 PCK::PredicateStats stats_side3h("side3h");
315 PCK::PredicateStats stats_side4("side4/insphere");
316 PCK::PredicateStats stats_orient2d("orient2d");
317 PCK::PredicateStats stats_orient3d("orient3d");
318 PCK::PredicateStats stats_orient3dh("orient3dh");
319 PCK::PredicateStats stats_det3d("det3d");
320 PCK::PredicateStats stats_det4d("det4d");
321
322 // ================= side1 =========================================
323
324 /**
325 * \brief Exact implementation of the side1() predicate using low-level
326 * exact arithmetics API (expansion class).
327 */
328 207030 Sign side1_exact_SOS(
329 const double* p0, const double* p1,
330 const double* q0,
331 coord_index_t dim
332 ) {
333 207030 stats_side1.log_exact();
334
2/6
✓ Branch 6 taken 207030 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 207030 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 207030 expansion& l = expansion_sq_dist(p0, p1, dim);
335
2/6
✓ Branch 6 taken 207030 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 207030 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 207030 expansion& a = expansion_dot_at(p1, q0, p0, dim).scale_fast(2.0);
336
2/6
✓ Branch 6 taken 207030 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 207030 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 207030 expansion& r = expansion_diff(l, a);
337 207030 Sign r_sign = r.sign();
338 // Symbolic perturbation, Simulation of Simplicity
339
2/2
✓ Branch 0 taken 204580 times.
✓ Branch 1 taken 2450 times.
 207030 if(r_sign == ZERO) {
340 204580 stats_side1.log_SOS();
341
2/2
✓ Branch 0 taken 154574 times.
✓ Branch 1 taken 50006 times.
 204580 return (p0 < p1) ? POSITIVE : NEGATIVE;
342 }
343 2450 return r_sign;
344 }
345
346 /**
347 * \brief Implements side1() in 3d.
348 */
349 3914542 Sign side1_3d_SOS(
350 const double* p0, const double* p1, const double* q0
351 ) {
352 3914542 Sign result = Sign(side1_3d_filter(p0, p1, q0));
353
2/2
✓ Branch 0 taken 207030 times.
✓ Branch 1 taken 3707512 times.
 3914542 if(result == ZERO) {
354 207030 result = side1_exact_SOS(p0, p1, q0, 3);
355 }
356 3914542 return result;
357 }
358
359 /**
360 * \brief Implements side1() in 4d.
361 */
362 50792 Sign side1_4d_SOS(
363 const double* p0, const double* p1, const double* q0
364 ) {
365 50792 Sign result = Sign(side1_4d_filter(p0, p1, q0));
366
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 50792 times.
 50792 if(result == ZERO) {
367 ✗ result = side1_exact_SOS(p0, p1, q0, 4);
368 }
369 50792 return result;
370 }
371
372 /**
373 * \brief Implements side1() in 6d.
374 */
375 59448 Sign side1_6d_SOS(
376 const double* p0, const double* p1, const double* q0
377 ) {
378 59448 Sign result = Sign(side1_6d_filter(p0, p1, q0));
379
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 59448 times.
 59448 if(result == ZERO) {
380 ✗ result = side1_exact_SOS(p0, p1, q0, 6);
381 }
382 59448 return result;
383 }
384
385 /**
386 * \brief Implements side1() in 7d.
387 */
388 ✗ Sign side1_7d_SOS(
389 const double* p0, const double* p1, const double* q0
390 ) {
391 ✗ Sign result = Sign(side1_7d_filter(p0, p1, q0));
392 ✗ if(result == ZERO) {
393 ✗ result = side1_exact_SOS(p0, p1, q0, 7);
394 }
395 ✗ return result;
396 }
397
398 /**
399 * \brief Implements side1() in 8d.
400 */
401 57825 Sign side1_8d_SOS(
402 const double* p0, const double* p1, const double* q0
403 ) {
404 57825 Sign result = Sign(side1_8d_filter(p0, p1, q0));
405
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 57825 times.
 57825 if(result == ZERO) {
406 ✗ result = side1_exact_SOS(p0, p1, q0, 8);
407 }
408 57825 return result;
409 }
410
411 // ================= side2 =========================================
412
413 /**
414 * \brief Exact implementation of the side2() predicate using low-level
415 * exact arithmetics API (expansion class).
416 */
417 279520 Sign side2_exact_SOS(
418 const double* p0, const double* p1, const double* p2,
419 const double* q0, const double* q1,
420 coord_index_t dim
421 ) {
422 279520 stats_side2.log_exact();
423
424
2/6
✓ Branch 6 taken 279520 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 279520 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 279520 const expansion& l1 = expansion_sq_dist(p1, p0, dim);
425
2/6
✓ Branch 6 taken 279520 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 279520 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 279520 const expansion& l2 = expansion_sq_dist(p2, p0, dim);
426
427
2/6
✓ Branch 6 taken 279520 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 279520 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 279520 const expansion& a10 = expansion_dot_at(p1,q0,p0, dim).scale_fast(2.0);
428
2/6
✓ Branch 6 taken 279520 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 279520 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 279520 const expansion& a11 = expansion_dot_at(p1,q1,p0, dim).scale_fast(2.0);
429
2/6
✓ Branch 6 taken 279520 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 279520 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 279520 const expansion& a20 = expansion_dot_at(p2,q0,p0, dim).scale_fast(2.0);
430
2/6
✓ Branch 6 taken 279520 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 279520 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 279520 const expansion& a21 = expansion_dot_at(p2,q1,p0, dim).scale_fast(2.0);
431
432
2/6
✓ Branch 6 taken 279520 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 279520 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 279520 const expansion& Delta = expansion_diff(a11, a10);
433
434 279520 Sign Delta_sign = Delta.sign();
435 // Should not occur with symbolic
436 // perturbation done at previous steps.
437
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 279520 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 279520 geo_assert(Delta_sign != ZERO);
438
439 // [ Lambda0 ] [ -1 ] [ a11 ]
440 // Delta [ ] = [ ] * l1 + [ ]
441 // [ Lambda1 ] [ 1 ] [ -a10 ]
442
443
2/6
✓ Branch 6 taken 279520 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 279520 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 279520 const expansion& DeltaLambda0 = expansion_diff(a11, l1);
444
2/6
✓ Branch 6 taken 279520 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 279520 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 279520 const expansion& DeltaLambda1 = expansion_diff(l1, a10);
445
446 // r = Delta*l2 - ( a20*DeltaLambda0 + a21*DeltaLambda1 )
447
448
2/6
✓ Branch 6 taken 279520 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 279520 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 279520 const expansion& r0 = expansion_product(Delta, l2);
449
2/6
✓ Branch 6 taken 279520 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 279520 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 279520 const expansion& r1 = expansion_product(a20, DeltaLambda0).negate();
450
2/6
✓ Branch 6 taken 279520 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 279520 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 279520 const expansion& r2 = expansion_product(a21, DeltaLambda1).negate();
451
2/6
✓ Branch 6 taken 279520 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 279520 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 279520 const expansion& r = expansion_sum3(r0, r1, r2);
452
453 279520 Sign r_sign = r.sign();
454
455 // Simulation of Simplicity (symbolic perturbation)
456
2/2
✓ Branch 0 taken 271431 times.
✓ Branch 1 taken 8089 times.
 279520 if(r_sign == ZERO) {
457 271431 stats_side2.log_SOS();
458 271431 const double* p_sort[3] = {p0, p1, p2};
459
1/2
✓ Branch 1 taken 271431 times.
✗ Branch 2 not taken.
 271431 SOS_sort(p_sort, p_sort + 3, dim);
460
1/2
✓ Branch 0 taken 324648 times.
✗ Branch 1 not taken.
 324648 for(index_t i = 0; i < 3; ++i) {
461
2/2
✓ Branch 0 taken 129865 times.
✓ Branch 1 taken 194783 times.
 324648 if(p_sort[i] == p0) {
462
2/6
✓ Branch 6 taken 129865 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 129865 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 129865 const expansion& z1 = expansion_diff(Delta, a21);
463
2/6
✓ Branch 6 taken 129865 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 129865 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 129865 const expansion& z = expansion_sum(z1, a20);
464
1/2
✓ Branch 1 taken 129865 times.
✗ Branch 2 not taken.
 129865 Sign z_sign = z.sign();
465
2/2
✓ Branch 0 taken 96912 times.
✓ Branch 1 taken 32953 times.
 129865 if(z_sign != ZERO) {
466 96912 return Sign(Delta_sign * z_sign);
467 }
468 }
469
2/2
✓ Branch 0 taken 134337 times.
✓ Branch 1 taken 93399 times.
 227736 if(p_sort[i] == p1) {
470
2/6
✓ Branch 6 taken 134337 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 134337 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 134337 const expansion& z = expansion_diff(a21, a20);
471
1/2
✓ Branch 1 taken 134337 times.
✗ Branch 2 not taken.
 134337 Sign z_sign = z.sign();
472
2/2
✓ Branch 0 taken 114073 times.
✓ Branch 1 taken 20264 times.
 134337 if(z_sign != ZERO) {
473 114073 return Sign(Delta_sign * z_sign);
474 }
475 }
476
2/2
✓ Branch 0 taken 60446 times.
✓ Branch 1 taken 53217 times.
 113663 if(p_sort[i] == p2) {
477 60446 return NEGATIVE;
478 }
479 }
480 ✗ geo_assert_not_reached;
481 }
482
483 8089 return Sign(Delta_sign * r_sign);
484 }
485
486 /**
487 * \brief Implements side2() in 3d.
488 */
489 5090168 Sign side2_3d_SOS(
490 const double* p0, const double* p1, const double* p2,
491 const double* q0, const double* q1
492 ) {
493 5090168 Sign result = Sign(side2_3d_filter(p0, p1, p2, q0, q1));
494
2/2
✓ Branch 0 taken 279520 times.
✓ Branch 1 taken 4810648 times.
 5090168 if(result == ZERO) {
495 279520 result = side2_exact_SOS(p0, p1, p2, q0, q1, 3);
496 }
497 5090168 return result;
498 }
499
500 /**
501 * \brief Implements side2() in 4d.
502 */
503 189360 Sign side2_4d_SOS(
504 const double* p0, const double* p1, const double* p2,
505 const double* q0, const double* q1
506 ) {
507 189360 Sign result = Sign(side2_4d_filter(p0, p1, p2, q0, q1));
508
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 189360 times.
 189360 if(result == ZERO) {
509 ✗ result = side2_exact_SOS(p0, p1, p2, q0, q1, 4);
510 }
511 189360 return result;
512 }
513
514 /**
515 * \brief Implements side2() in 6d.
516 */
517 260704 Sign side2_6d_SOS(
518 const double* p0, const double* p1, const double* p2,
519 const double* q0, const double* q1
520 ) {
521 260704 Sign result = Sign(side2_6d_filter(p0, p1, p2, q0, q1));
522
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 260704 times.
 260704 if(result == ZERO) {
523 ✗ result = side2_exact_SOS(p0, p1, p2, q0, q1, 6);
524 }
525 260704 return result;
526 }
527
528 /**
529 * \brief Implements side2() in 7d.
530 */
531 ✗ Sign side2_7d_SOS(
532 const double* p0, const double* p1, const double* p2,
533 const double* q0, const double* q1
534 ) {
535 ✗ Sign result = Sign(side2_7d_filter(p0, p1, p2, q0, q1));
536 ✗ if(result == ZERO) {
537 ✗ result = side2_exact_SOS(p0, p1, p2, q0, q1, 7);
538 }
539 ✗ return result;
540 }
541
542 /**
543 * \brief Implements side2() in 8d.
544 */
545 258047 Sign side2_8d_SOS(
546 const double* p0, const double* p1, const double* p2,
547 const double* q0, const double* q1
548 ) {
549 258047 Sign result = Sign(side2_8d_filter(p0, p1, p2, q0, q1));
550
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 258047 times.
 258047 if(result == ZERO) {
551 ✗ result = side2_exact_SOS(p0, p1, p2, q0, q1, 8);
552 }
553 258047 return result;
554 }
555
556 // ================= side3 =========================================
557
558 /**
559 * \brief Exact implementation of the side3() predicate using low-level
560 * exact arithmetics API (expansion class).
561 */
562 80656 Sign side3_exact_SOS(
563 const double* p0, const double* p1, const double* p2, const double* p3,
564 const double* q0, const double* q1, const double* q2,
565 coord_index_t dim
566 ) {
567 80656 stats_side3.log_exact();
568
569
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& l1 = expansion_sq_dist(p1, p0, dim);
570
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& l2 = expansion_sq_dist(p2, p0, dim);
571
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& l3 = expansion_sq_dist(p3, p0, dim);
572
573
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 80656 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 80656 const expansion& a10 = expansion_dot_at(p1,q0,p0, dim).scale_fast(2.0);
574
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 80656 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 80656 const expansion& a11 = expansion_dot_at(p1,q1,p0, dim).scale_fast(2.0);
575
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 80656 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 80656 const expansion& a12 = expansion_dot_at(p1,q2,p0, dim).scale_fast(2.0);
576
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 80656 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 80656 const expansion& a20 = expansion_dot_at(p2,q0,p0, dim).scale_fast(2.0);
577
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 80656 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 80656 const expansion& a21 = expansion_dot_at(p2,q1,p0, dim).scale_fast(2.0);
578
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 80656 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 80656 const expansion& a22 = expansion_dot_at(p2,q2,p0, dim).scale_fast(2.0);
579
580
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 80656 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 80656 const expansion& a30 = expansion_dot_at(p3,q0,p0, dim).scale_fast(2.0);
581
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 80656 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 80656 const expansion& a31 = expansion_dot_at(p3,q1,p0, dim).scale_fast(2.0);
582
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 80656 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 80656 const expansion& a32 = expansion_dot_at(p3,q2,p0, dim).scale_fast(2.0);
583
584 // [ b00 b01 b02 ] [ 1 1 1 ]-1
585 // [ b10 b11 b12 ] = Delta * [ a10 a11 a12 ]
586 // [ b20 b21 b22 ] [ a20 a21 a22 ]
587
588
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& b00 = expansion_det2x2(a11, a12, a21, a22);
589
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& b01 = expansion_diff(a21, a22);
590
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& b02 = expansion_diff(a12, a11);
591
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& b10 = expansion_det2x2(a12, a10, a22, a20);
592
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& b11 = expansion_diff(a22, a20);
593
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& b12 = expansion_diff(a10, a12);
594
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& b20 = expansion_det2x2(a10, a11, a20, a21);
595
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& b21 = expansion_diff(a20, a21);
596
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& b22 = expansion_diff(a11, a10);
597
598
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& Delta = expansion_sum3(b00, b10, b20);
599 80656 Sign Delta_sign = Delta.sign();
600 // Should not occur with symbolic
601 // perturbation done at previous steps.
602
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 80656 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 80656 geo_assert(Delta_sign != ZERO);
603
604 // [ Lambda0 ] [ b01 b02 ] [ l1 ] [ b00 ]
605 // Delta [ Lambda1 ] = [ b11 b12 ] * [ ] + [ b10 ]
606 // [ Lambda2 ] [ b21 b22 ] [ l2 ] [ b20 ]
607
608
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& b01_l1 = expansion_product(b01, l1);
609
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& b02_l2 = expansion_product(b02, l2);
610
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& DeltaLambda0 = expansion_sum3(b01_l1, b02_l2, b00);
611
612
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& b11_l1 = expansion_product(b11, l1);
613
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& b12_l2 = expansion_product(b12, l2);
614
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& DeltaLambda1 = expansion_sum3(b11_l1, b12_l2, b10);
615
616
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& b21_l1 = expansion_product(b21, l1);
617
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& b22_l2 = expansion_product(b22, l2);
618
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& DeltaLambda2 = expansion_sum3(b21_l1, b22_l2, b20);
619
620 // r = Delta*l3-(a30*DeltaLambda0+a31*DeltaLambda1+a32*DeltaLambda2)
621
622
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& r0 = expansion_product(Delta, l3);
623
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 80656 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 80656 const expansion& r1 = expansion_product(a30, DeltaLambda0).negate();
624
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 80656 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 80656 const expansion& r2 = expansion_product(a31, DeltaLambda1).negate();
625
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 80656 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 80656 const expansion& r3 = expansion_product(a32, DeltaLambda2).negate();
626
2/6
✓ Branch 6 taken 80656 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 80656 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 80656 const expansion& r = expansion_sum4(r0, r1, r2, r3);
627 80656 Sign r_sign = r.sign();
628
629 // Simulation of Simplicity (symbolic perturbation)
630
2/2
✓ Branch 0 taken 72918 times.
✓ Branch 1 taken 7738 times.
 80656 if(r_sign == ZERO) {
631 72918 stats_side3.log_SOS();
632 72918 const double* p_sort[4] = {p0, p1, p2, p3};
633
1/2
✓ Branch 1 taken 72918 times.
✗ Branch 2 not taken.
 72918 SOS_sort(p_sort, p_sort + 4, dim);
634
1/2
✓ Branch 0 taken 104050 times.
✗ Branch 1 not taken.
 104050 for(index_t i = 0; i < 4; ++i) {
635
2/2
✓ Branch 0 taken 35915 times.
✓ Branch 1 taken 68135 times.
 104050 if(p_sort[i] == p0) {
636
2/6
✓ Branch 6 taken 35915 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 35915 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 35915 const expansion& z1_0 = expansion_sum(b01, b02);
637
2/6
✓ Branch 6 taken 35915 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 35915 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 35915 const expansion& z1 = expansion_product(a30, z1_0).negate();
638
2/6
✓ Branch 6 taken 35915 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 35915 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 35915 const expansion& z2_0 = expansion_sum(b11, b12);
639
2/6
✓ Branch 6 taken 35915 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 35915 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 35915 const expansion& z2 = expansion_product(a31, z2_0).negate();
640
2/6
✓ Branch 6 taken 35915 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 35915 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 35915 const expansion& z3_0 = expansion_sum(b21, b22);
641
2/6
✓ Branch 6 taken 35915 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 35915 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 35915 const expansion& z3 = expansion_product(a32, z3_0).negate();
642
2/6
✓ Branch 6 taken 35915 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 35915 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 35915 const expansion& z = expansion_sum4(Delta, z1, z2, z3);
643
1/2
✓ Branch 1 taken 35915 times.
✗ Branch 2 not taken.
 35915 Sign z_sign = z.sign();
644
2/2
✓ Branch 0 taken 15601 times.
✓ Branch 1 taken 20314 times.
 35915 if(z_sign != ZERO) {
645 15601 return Sign(Delta_sign * z_sign);
646 }
647
2/2
✓ Branch 0 taken 10534 times.
✓ Branch 1 taken 57601 times.
 68135 } else if(p_sort[i] == p1) {
648
2/6
✓ Branch 6 taken 10534 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 10534 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 10534 const expansion& z1 = expansion_product(a30, b01);
649
2/6
✓ Branch 6 taken 10534 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 10534 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 10534 const expansion& z2 = expansion_product(a31, b11);
650
2/6
✓ Branch 6 taken 10534 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 10534 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 10534 const expansion& z3 = expansion_product(a32, b21);
651
2/6
✓ Branch 6 taken 10534 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 10534 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 10534 const expansion& z = expansion_sum3(z1, z2, z3);
652
1/2
✓ Branch 1 taken 10534 times.
✗ Branch 2 not taken.
 10534 Sign z_sign = z.sign();
653
2/2
✓ Branch 0 taken 10393 times.
✓ Branch 1 taken 141 times.
 10534 if(z_sign != ZERO) {
654 10393 return Sign(Delta_sign * z_sign);
655 }
656
2/2
✓ Branch 0 taken 45832 times.
✓ Branch 1 taken 11769 times.
 57601 } else if(p_sort[i] == p2) {
657
2/6
✓ Branch 6 taken 45832 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 45832 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 45832 const expansion& z1 = expansion_product(a30, b02);
658
2/6
✓ Branch 6 taken 45832 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 45832 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 45832 const expansion& z2 = expansion_product(a31, b12);
659
2/6
✓ Branch 6 taken 45832 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 45832 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 45832 const expansion& z3 = expansion_product(a32, b22);
660
2/6
✓ Branch 6 taken 45832 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 45832 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 45832 const expansion& z = expansion_sum3(z1, z2, z3);
661
1/2
✓ Branch 1 taken 45832 times.
✗ Branch 2 not taken.
 45832 Sign z_sign = z.sign();
662
2/2
✓ Branch 0 taken 35155 times.
✓ Branch 1 taken 10677 times.
 45832 if(z_sign != ZERO) {
663 35155 return Sign(Delta_sign * z_sign);
664 }
665
1/2
✓ Branch 0 taken 11769 times.
✗ Branch 1 not taken.
 11769 } else if(p_sort[i] == p3) {
666 11769 return NEGATIVE;
667 }
668 }
669 ✗ geo_assert_not_reached;
670 }
671 7738 return Sign(Delta_sign * r_sign);
672 }
673
674
675 /**
676 * \brief Exact implementation of the side3_3dlifted() predicate
677 * using low-level exact arithmetics API (expansion class).
678 */
679 ✗ Sign side3h_exact_SOS(
680 const double* p0, const double* p1, const double* p2, const double* p3,
681 double h0, double h1, double h2, double h3,
682 const double* q0, const double* q1, const double* q2
683 ) {
684 ✗ stats_side3h.log_exact();
685
686 ✗ const expansion& l1 = expansion_diff(h1,h0);
687 ✗ const expansion& l2 = expansion_diff(h2,h0);
688 ✗ const expansion& l3 = expansion_diff(h3,h0);
689
690 ✗ const expansion& a10 = expansion_dot_at(p1, q0, p0, 3).scale_fast(2.0);
691 ✗ const expansion& a11 = expansion_dot_at(p1, q1, p0, 3).scale_fast(2.0);
692 ✗ const expansion& a12 = expansion_dot_at(p1, q2, p0, 3).scale_fast(2.0);
693 ✗ const expansion& a20 = expansion_dot_at(p2, q0, p0, 3).scale_fast(2.0);
694 ✗ const expansion& a21 = expansion_dot_at(p2, q1, p0, 3).scale_fast(2.0);
695 ✗ const expansion& a22 = expansion_dot_at(p2, q2, p0, 3).scale_fast(2.0);
696
697 ✗ const expansion& a30 = expansion_dot_at(p3, q0, p0, 3).scale_fast(2.0);
698 ✗ const expansion& a31 = expansion_dot_at(p3, q1, p0, 3).scale_fast(2.0);
699 ✗ const expansion& a32 = expansion_dot_at(p3, q2, p0, 3).scale_fast(2.0);
700
701 // [ b00 b01 b02 ] [ 1 1 1 ]-1
702 // [ b10 b11 b12 ] = Delta * [ a10 a11 a12 ]
703 // [ b20 b21 b22 ] [ a20 a21 a22 ]
704
705 ✗ const expansion& b00 = expansion_det2x2(a11, a12, a21, a22);
706 ✗ const expansion& b01 = expansion_diff(a21, a22);
707 ✗ const expansion& b02 = expansion_diff(a12, a11);
708 ✗ const expansion& b10 = expansion_det2x2(a12, a10, a22, a20);
709 ✗ const expansion& b11 = expansion_diff(a22, a20);
710 ✗ const expansion& b12 = expansion_diff(a10, a12);
711 ✗ const expansion& b20 = expansion_det2x2(a10, a11, a20, a21);
712 ✗ const expansion& b21 = expansion_diff(a20, a21);
713 ✗ const expansion& b22 = expansion_diff(a11, a10);
714
715 ✗ const expansion& Delta = expansion_sum3(b00, b10, b20);
716 ✗ Sign Delta_sign = Delta.sign();
717 // Should not occur with symbolic
718 // perturbation done at previous steps.
719 ✗ geo_assert(Delta_sign != ZERO);
720
721 // [ Lambda0 ] [ b01 b02 ] [ l1 ] [ b00 ]
722 // Delta [ Lambda1 ] = [ b11 b12 ] * [ ] + [ b10 ]
723 // [ Lambda2 ] [ b21 b22 ] [ l2 ] [ b20 ]
724
725 ✗ const expansion& b01_l1 = expansion_product(b01, l1);
726 ✗ const expansion& b02_l2 = expansion_product(b02, l2);
727 ✗ const expansion& DeltaLambda0 = expansion_sum3(b01_l1, b02_l2, b00);
728
729 ✗ const expansion& b11_l1 = expansion_product(b11, l1);
730 ✗ const expansion& b12_l2 = expansion_product(b12, l2);
731 ✗ const expansion& DeltaLambda1 = expansion_sum3(b11_l1, b12_l2, b10);
732
733 ✗ const expansion& b21_l1 = expansion_product(b21, l1);
734 ✗ const expansion& b22_l2 = expansion_product(b22, l2);
735 ✗ const expansion& DeltaLambda2 = expansion_sum3(b21_l1, b22_l2, b20);
736
737 // r = Delta*l3-(a30*DeltaLambda0+a31*DeltaLambda1+a32*DeltaLambda2)
738
739 ✗ const expansion& r0 = expansion_product(Delta, l3);
740 ✗ const expansion& r1 = expansion_product(a30, DeltaLambda0).negate();
741 ✗ const expansion& r2 = expansion_product(a31, DeltaLambda1).negate();
742 ✗ const expansion& r3 = expansion_product(a32, DeltaLambda2).negate();
743 ✗ const expansion& r = expansion_sum4(r0, r1, r2, r3);
744 ✗ Sign r_sign = r.sign();
745
746 // Simulation of Simplicity (symbolic perturbation)
747 ✗ if(r_sign == ZERO) {
748 ✗ stats_side3h.log_SOS();
749 ✗ const double* p_sort[4] = {p0, p1, p2, p3};
750 ✗ SOS_sort(p_sort, p_sort + 4, 3);
751 ✗ for(index_t i = 0; i < 4; ++i) {
752 ✗ if(p_sort[i] == p0) {
753 ✗ const expansion& z1_0 = expansion_sum(b01, b02);
754 ✗ const expansion& z1 = expansion_product(a30, z1_0).negate();
755 ✗ const expansion& z2_0 = expansion_sum(b11, b12);
756 ✗ const expansion& z2 = expansion_product(a31, z2_0).negate();
757 ✗ const expansion& z3_0 = expansion_sum(b21, b22);
758 ✗ const expansion& z3 = expansion_product(a32, z3_0).negate();
759 ✗ const expansion& z = expansion_sum4(Delta, z1, z2, z3);
760 ✗ Sign z_sign = z.sign();
761 ✗ if(z_sign != ZERO) {
762 ✗ return Sign(Delta_sign * z_sign);
763 }
764 ✗ } else if(p_sort[i] == p1) {
765 ✗ const expansion& z1 = expansion_product(a30, b01);
766 ✗ const expansion& z2 = expansion_product(a31, b11);
767 ✗ const expansion& z3 = expansion_product(a32, b21);
768 ✗ const expansion& z = expansion_sum3(z1, z2, z3);
769 ✗ Sign z_sign = z.sign();
770 ✗ if(z_sign != ZERO) {
771 ✗ return Sign(Delta_sign * z_sign);
772 }
773 ✗ } else if(p_sort[i] == p2) {
774 ✗ const expansion& z1 = expansion_product(a30, b02);
775 ✗ const expansion& z2 = expansion_product(a31, b12);
776 ✗ const expansion& z3 = expansion_product(a32, b22);
777 ✗ const expansion& z = expansion_sum3(z1, z2, z3);
778 ✗ Sign z_sign = z.sign();
779 ✗ if(z_sign != ZERO) {
780 ✗ return Sign(Delta_sign * z_sign);
781 }
782 ✗ } else if(p_sort[i] == p3) {
783 ✗ return NEGATIVE;
784 }
785 }
786 ✗ geo_assert_not_reached;
787 }
788 ✗ return Sign(Delta_sign * r_sign);
789 }
790
791
792 /**
793 * \brief Implements side3() in 3d.
794 */
795 11053772 Sign side3_3d_SOS(
796 const double* p0, const double* p1, const double* p2, const double* p3,
797 const double* q0, const double* q1, const double* q2
798 ) {
799 11053772 Sign result = Sign(side3_3d_filter(p0, p1, p2, p3, q0, q1, q2));
800
2/2
✓ Branch 0 taken 80628 times.
✓ Branch 1 taken 10973144 times.
 11053772 if(result == ZERO) {
801 80628 result = side3_exact_SOS(p0, p1, p2, p3, q0, q1, q2, 3);
802 }
803 11053772 return result;
804 }
805
806
807 /**
808 * \brief Implements side3() in 4d.
809 */
810 149639 Sign side3_4d_SOS(
811 const double* p0, const double* p1, const double* p2, const double* p3,
812 const double* q0, const double* q1, const double* q2
813 ) {
814 149639 Sign result = Sign(side3_4d_filter(p0, p1, p2, p3, q0, q1, q2));
815
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 149639 times.
 149639 if(result == ZERO) {
816 ✗ result = side3_exact_SOS(p0, p1, p2, p3, q0, q1, q2, 4);
817 }
818 149639 return result;
819 }
820
821 /**
822 * \brief Implements side3() in 6d.
823 */
824 165662 Sign side3_6d_SOS(
825 const double* p0, const double* p1, const double* p2, const double* p3,
826 const double* q0, const double* q1, const double* q2
827 ) {
828 165662 Sign result = Sign(side3_6d_filter(p0, p1, p2, p3, q0, q1, q2));
829
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 165662 times.
 165662 if(result == ZERO) {
830 ✗ result = side3_exact_SOS(p0, p1, p2, p3, q0, q1, q2, 6);
831 }
832 165662 return result;
833 }
834
835 /**
836 * \brief Implements side3() in 7d.
837 */
838 ✗ Sign side3_7d_SOS(
839 const double* p0, const double* p1, const double* p2, const double* p3,
840 const double* q0, const double* q1, const double* q2
841 ) {
842 ✗ Sign result = Sign(side3_7d_filter(p0, p1, p2, p3, q0, q1, q2));
843 ✗ if(result == ZERO) {
844 ✗ result = side3_exact_SOS(p0, p1, p2, p3, q0, q1, q2, 7);
845 }
846 ✗ return result;
847 }
848
849 /**
850 * \brief Implements side3() in 7d.
851 */
852 154806 Sign side3_8d_SOS(
853 const double* p0, const double* p1, const double* p2, const double* p3,
854 const double* q0, const double* q1, const double* q2
855 ) {
856 154806 Sign result = Sign(side3_8d_filter(p0, p1, p2, p3, q0, q1, q2));
857
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 154806 times.
 154806 if(result == ZERO) {
858 ✗ result = side3_exact_SOS(p0, p1, p2, p3, q0, q1, q2, 8);
859 }
860 154806 return result;
861 }
862
863 // ================= side4 =========================================
864
865 /**
866 * \brief Exact implementation of the side4_3d_SOS() predicate
867 * using low-level exact arithmetics API (expansion class).
868 * \param[in] sos if true, applies symbolic perturbation when
869 * result is zero, else returns zero
870 */
871 252392 Sign side4_3d_exact_SOS(
872 const double* p0, const double* p1, const double* p2, const double* p3,
873 const double* p4, bool sos = true
874 ) {
875 252392 stats_side4.log_exact();
876
877
3/8
✓ Branch 4 taken 252392 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 252392 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 252392 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 252392 const expansion& a11 = expansion_diff(p1[0], p0[0]);
878
3/8
✓ Branch 4 taken 252392 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 252392 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 252392 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 252392 const expansion& a12 = expansion_diff(p1[1], p0[1]);
879
3/8
✓ Branch 4 taken 252392 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 252392 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 252392 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 252392 const expansion& a13 = expansion_diff(p1[2], p0[2]);
880
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 252392 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 252392 const expansion& a14 = expansion_sq_dist(p1, p0, 3).negate();
881
882
3/8
✓ Branch 4 taken 252392 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 252392 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 252392 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 252392 const expansion& a21 = expansion_diff(p2[0], p0[0]);
883
3/8
✓ Branch 4 taken 252392 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 252392 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 252392 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 252392 const expansion& a22 = expansion_diff(p2[1], p0[1]);
884
3/8
✓ Branch 4 taken 252392 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 252392 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 252392 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 252392 const expansion& a23 = expansion_diff(p2[2], p0[2]);
885
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 252392 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 252392 const expansion& a24 = expansion_sq_dist(p2, p0, 3).negate();
886
887
3/8
✓ Branch 4 taken 252392 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 252392 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 252392 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 252392 const expansion& a31 = expansion_diff(p3[0], p0[0]);
888
3/8
✓ Branch 4 taken 252392 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 252392 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 252392 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 252392 const expansion& a32 = expansion_diff(p3[1], p0[1]);
889
3/8
✓ Branch 4 taken 252392 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 252392 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 252392 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 252392 const expansion& a33 = expansion_diff(p3[2], p0[2]);
890
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 252392 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 252392 const expansion& a34 = expansion_sq_dist(p3, p0, 3).negate();
891
892
3/8
✓ Branch 4 taken 252392 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 252392 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 252392 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 252392 const expansion& a41 = expansion_diff(p4[0], p0[0]);
893
3/8
✓ Branch 4 taken 252392 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 252392 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 252392 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 252392 const expansion& a42 = expansion_diff(p4[1], p0[1]);
894
3/8
✓ Branch 4 taken 252392 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 252392 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 252392 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 252392 const expansion& a43 = expansion_diff(p4[2], p0[2]);
895
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 252392 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 252392 const expansion& a44 = expansion_sq_dist(p4, p0, 3).negate();
896
897 // This commented-out version does not reuse
898 // the 2x2 minors.
899 /*
900 const expansion& Delta1 = expansion_det3x3(
901 a21, a22, a23,
902 a31, a32, a33,
903 a41, a42, a43
904 );
905 const expansion& Delta2 = expansion_det3x3(
906 a11, a12, a13,
907 a31, a32, a33,
908 a41, a42, a43
909 );
910 const expansion& Delta3 = expansion_det3x3(
911 a11, a12, a13,
912 a21, a22, a23,
913 a41, a42, a43
914 );
915 const expansion& Delta4 = expansion_det3x3(
916 a11, a12, a13,
917 a21, a22, a23,
918 a31, a32, a33
919 );
920 */
921
922 // Optimized version that reuses the 2x2 minors
923
924
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& m12 = expansion_det2x2(a12,a13,a22,a23);
925
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& m13 = expansion_det2x2(a12,a13,a32,a33);
926
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& m14 = expansion_det2x2(a12,a13,a42,a43);
927
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& m23 = expansion_det2x2(a22,a23,a32,a33);
928
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& m24 = expansion_det2x2(a22,a23,a42,a43);
929
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& m34 = expansion_det2x2(a32,a33,a42,a43);
930
931
932
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& z11 = expansion_product(a21,m34);
933
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 252392 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 252392 const expansion& z12 = expansion_product(a31,m24).negate();
934
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& z13 = expansion_product(a41,m23);
935
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& Delta1 = expansion_sum3(z11,z12,z13);
936
937
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& z21 = expansion_product(a11,m34);
938
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 252392 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 252392 const expansion& z22 = expansion_product(a31,m14).negate();
939
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& z23 = expansion_product(a41,m13);
940
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& Delta2 = expansion_sum3(z21,z22,z23);
941
942
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& z31 = expansion_product(a11,m24);
943
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 252392 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 252392 const expansion& z32 = expansion_product(a21,m14).negate();
944
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& z33 = expansion_product(a41,m12);
945
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& Delta3 = expansion_sum3(z31,z32,z33);
946
947
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& z41 = expansion_product(a11,m23);
948
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 252392 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 252392 const expansion& z42 = expansion_product(a21,m13).negate();
949
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& z43 = expansion_product(a31,m12);
950
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& Delta4 = expansion_sum3(z41,z42,z43);
951
952
953 252392 Sign Delta4_sign = Delta4.sign();
954
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 252392 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 252392 geo_assert(Delta4_sign != ZERO);
955
956
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& r_1 = expansion_product(Delta1, a14);
957
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 252392 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 252392 const expansion& r_2 = expansion_product(Delta2, a24).negate();
958
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& r_3 = expansion_product(Delta3, a34);
959
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 252392 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 252392 const expansion& r_4 = expansion_product(Delta4, a44).negate();
960
2/6
✓ Branch 6 taken 252392 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 252392 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 252392 const expansion& r = expansion_sum4(r_1, r_2, r_3, r_4);
961 252392 Sign r_sign = r.sign();
962
963 // Simulation of Simplicity (symbolic perturbation)
964
3/4
✓ Branch 0 taken 252392 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 94751 times.
✓ Branch 3 taken 157641 times.
 252392 if(sos && r_sign == ZERO) {
965 94751 stats_side4.log_SOS();
966 94751 const double* p_sort[5] = {p0, p1, p2, p3, p4};
967
1/2
✓ Branch 1 taken 94751 times.
✗ Branch 2 not taken.
 94751 SOS_sort(p_sort, p_sort + 5, 3);
968
1/2
✓ Branch 0 taken 274387 times.
✗ Branch 1 not taken.
 274387 for(index_t i = 0; i < 5; ++i) {
969
2/2
✓ Branch 0 taken 90383 times.
✓ Branch 1 taken 184004 times.
 274387 if(p_sort[i] == p0) {
970
2/6
✓ Branch 6 taken 90383 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 90383 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 90383 const expansion& z1 = expansion_diff(Delta2, Delta1);
971
2/6
✓ Branch 6 taken 90383 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 90383 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 90383 const expansion& z2 = expansion_diff(Delta4, Delta3);
972
2/6
✓ Branch 6 taken 90383 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 90383 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 90383 const expansion& z = expansion_sum(z1, z2);
973
1/2
✓ Branch 1 taken 90383 times.
✗ Branch 2 not taken.
 90383 Sign z_sign = z.sign();
974
2/2
✓ Branch 0 taken 965 times.
✓ Branch 1 taken 89418 times.
 90383 if(z_sign != ZERO) {
975 94751 return Sign(Delta4_sign * z_sign);
976 }
977
2/2
✓ Branch 0 taken 30921 times.
✓ Branch 1 taken 153083 times.
 184004 } else if(p_sort[i] == p1) {
978
1/2
✓ Branch 1 taken 30921 times.
✗ Branch 2 not taken.
 30921 Sign Delta1_sign = Delta1.sign();
979
2/2
✓ Branch 0 taken 30509 times.
✓ Branch 1 taken 412 times.
 30921 if(Delta1_sign != ZERO) {
980 30509 return Sign(Delta4_sign * Delta1_sign);
981 }
982
2/2
✓ Branch 0 taken 60662 times.
✓ Branch 1 taken 92421 times.
 153083 } else if(p_sort[i] == p2) {
983
1/2
✓ Branch 1 taken 60662 times.
✗ Branch 2 not taken.
 60662 Sign Delta2_sign = Delta2.sign();
984
2/2
✓ Branch 0 taken 30577 times.
✓ Branch 1 taken 30085 times.
 60662 if(Delta2_sign != ZERO) {
985 30577 return Sign(-Delta4_sign * Delta2_sign);
986 }
987
2/2
✓ Branch 0 taken 90327 times.
✓ Branch 1 taken 2094 times.
 92421 } else if(p_sort[i] == p3) {
988
1/2
✓ Branch 1 taken 90327 times.
✗ Branch 2 not taken.
 90327 Sign Delta3_sign = Delta3.sign();
989
2/2
✓ Branch 0 taken 30606 times.
✓ Branch 1 taken 59721 times.
 90327 if(Delta3_sign != ZERO) {
990 30606 return Sign(Delta4_sign * Delta3_sign);
991 }
992
1/2
✓ Branch 0 taken 2094 times.
✗ Branch 1 not taken.
 2094 } else if(p_sort[i] == p4) {
993 2094 return NEGATIVE;
994 }
995 }
996 }
997 157641 return Sign(Delta4_sign * r_sign);
998 }
999
1000 /**
1001 * \brief Exact implementation of the side4() predicate using low-level
1002 * exact arithmetics API (expansion class).
1003 */
1004 ✗ Sign side4_exact_SOS(
1005 const double* p0, const double* p1, const double* p2, const double* p3,
1006 const double* p4,
1007 const double* q0, const double* q1, const double* q2, const double* q3,
1008 coord_index_t dim
1009 ) {
1010 ✗ stats_side4.log_exact();
1011
1012 ✗ const expansion& l1 = expansion_sq_dist(p1, p0, dim);
1013 ✗ const expansion& l2 = expansion_sq_dist(p2, p0, dim);
1014 ✗ const expansion& l3 = expansion_sq_dist(p3, p0, dim);
1015 ✗ const expansion& l4 = expansion_sq_dist(p4, p0, dim);
1016
1017 ✗ const expansion& a10 = expansion_dot_at(p1,q0,p0, dim).scale_fast(2.0);
1018 ✗ const expansion& a11 = expansion_dot_at(p1,q1,p0, dim).scale_fast(2.0);
1019 ✗ const expansion& a12 = expansion_dot_at(p1,q2,p0, dim).scale_fast(2.0);
1020 ✗ const expansion& a13 = expansion_dot_at(p1,q3,p0, dim).scale_fast(2.0);
1021
1022 ✗ const expansion& a20 = expansion_dot_at(p2,q0,p0, dim).scale_fast(2.0);
1023 ✗ const expansion& a21 = expansion_dot_at(p2,q1,p0, dim).scale_fast(2.0);
1024 ✗ const expansion& a22 = expansion_dot_at(p2,q2,p0, dim).scale_fast(2.0);
1025 ✗ const expansion& a23 = expansion_dot_at(p2,q3,p0, dim).scale_fast(2.0);
1026
1027 ✗ const expansion& a30 = expansion_dot_at(p3,q0,p0, dim).scale_fast(2.0);
1028 ✗ const expansion& a31 = expansion_dot_at(p3,q1,p0, dim).scale_fast(2.0);
1029 ✗ const expansion& a32 = expansion_dot_at(p3,q2,p0, dim).scale_fast(2.0);
1030 ✗ const expansion& a33 = expansion_dot_at(p3,q3,p0, dim).scale_fast(2.0);
1031
1032 ✗ const expansion& a40 = expansion_dot_at(p4,q0,p0, dim).scale_fast(2.0);
1033 ✗ const expansion& a41 = expansion_dot_at(p4,q1,p0, dim).scale_fast(2.0);
1034 ✗ const expansion& a42 = expansion_dot_at(p4,q2,p0, dim).scale_fast(2.0);
1035 ✗ const expansion& a43 = expansion_dot_at(p4,q3,p0, dim).scale_fast(2.0);
1036
1037 // [ b00 b01 b02 b03 ] [ 1 1 1 1 ]-1
1038 // [ b10 b11 b12 b13 ] [ a10 a11 a12 a13 ]
1039 // [ b20 b21 b22 b23 ] = Delta * [ a20 a21 a22 a23 ]
1040 // [ b30 b31 b32 b33 ] [ a30 a31 a32 a33 ]
1041
1042 // Note: we could probably reuse some of the co-factors
1043 // (but for now I'd rather keep this form that is easier to
1044 // read ... and to debug if need be !)
1045
1046 ✗ const expansion& b00 = expansion_det3x3(a11, a12, a13, a21, a22, a23, a31, a32, a33);
1047 ✗ const expansion& b01 = expansion_det_111_2x3(a21, a22, a23, a31, a32, a33).negate();
1048 ✗ const expansion& b02 = expansion_det_111_2x3(a11, a12, a13, a31, a32, a33);
1049 ✗ const expansion& b03 = expansion_det_111_2x3(a11, a12, a13, a21, a22, a23).negate();
1050
1051 ✗ const expansion& b10 = expansion_det3x3(a10, a12, a13, a20, a22, a23, a30, a32, a33).negate();
1052 ✗ const expansion& b11 = expansion_det_111_2x3(a20, a22, a23, a30, a32, a33);
1053 ✗ const expansion& b12 = expansion_det_111_2x3(a10, a12, a13, a30, a32, a33).negate();
1054 ✗ const expansion& b13 = expansion_det_111_2x3(a10, a12, a13, a20, a22, a23);
1055
1056 ✗ const expansion& b20 = expansion_det3x3(a10, a11, a13, a20, a21, a23, a30, a31, a33);
1057 ✗ const expansion& b21 = expansion_det_111_2x3(a20, a21, a23, a30, a31, a33).negate();
1058 ✗ const expansion& b22 = expansion_det_111_2x3(a10, a11, a13, a30, a31, a33);
1059 ✗ const expansion& b23 = expansion_det_111_2x3(a10, a11, a13, a20, a21, a23).negate();
1060
1061 ✗ const expansion& b30 = expansion_det3x3(a10, a11, a12, a20, a21, a22, a30, a31, a32).negate();
1062 ✗ const expansion& b31 = expansion_det_111_2x3(a20, a21, a22, a30, a31, a32);
1063 ✗ const expansion& b32 = expansion_det_111_2x3(a10, a11, a12, a30, a31, a32).negate();
1064 ✗ const expansion& b33 = expansion_det_111_2x3(a10, a11, a12, a20, a21, a22);
1065
1066 ✗ const expansion& Delta = expansion_sum4(b00, b10, b20, b30);
1067 ✗ Sign Delta_sign = Delta.sign();
1068 ✗ geo_assert(Delta_sign != ZERO);
1069
1070 // [ Lambda0 ] [ b01 b02 b03 ] [ l1 ] [ b00 ]
1071 // [ Lambda1 ] [ b11 b12 b13 ] [ l2 ] [ b10 ]
1072 // Delta [ Lambda2 ] = [ b21 b22 b23 ] * [ l3 ] + [ b20 ]
1073 // [ Lambda3 ] [ b31 b32 b33 ] [ l4 ] [ b30 ]
1074
1075 ✗ const expansion& b01_l1 = expansion_product(b01, l1);
1076 ✗ const expansion& b02_l2 = expansion_product(b02, l2);
1077 ✗ const expansion& b03_l3 = expansion_product(b03, l3);
1078 ✗ const expansion& DeltaLambda0 = expansion_sum4(b01_l1, b02_l2, b03_l3, b00);
1079
1080 ✗ const expansion& b11_l1 = expansion_product(b11, l1);
1081 ✗ const expansion& b12_l2 = expansion_product(b12, l2);
1082 ✗ const expansion& b13_l3 = expansion_product(b13, l3);
1083 ✗ const expansion& DeltaLambda1 = expansion_sum4(b11_l1, b12_l2, b13_l3, b10);
1084
1085 ✗ const expansion& b21_l1 = expansion_product(b21, l1);
1086 ✗ const expansion& b22_l2 = expansion_product(b22, l2);
1087 ✗ const expansion& b23_l3 = expansion_product(b23, l3);
1088 ✗ const expansion& DeltaLambda2 = expansion_sum4(b21_l1, b22_l2, b23_l3, b20);
1089
1090 ✗ const expansion& b31_l1 = expansion_product(b31, l1);
1091 ✗ const expansion& b32_l2 = expansion_product(b32, l2);
1092 ✗ const expansion& b33_l3 = expansion_product(b33, l3);
1093 ✗ const expansion& DeltaLambda3 = expansion_sum4(b31_l1, b32_l2, b33_l3, b30);
1094
1095 // r = Delta*l4 - (
1096 // a40*DeltaLambda0+
1097 // a41*DeltaLambda1+
1098 // a42*DeltaLambda2+
1099 // a43*DeltaLambda3
1100 // )
1101
1102 ✗ const expansion& r0 = expansion_product(Delta, l4);
1103 ✗ const expansion& r1 = expansion_product(a40, DeltaLambda0);
1104 ✗ const expansion& r2 = expansion_product(a41, DeltaLambda1);
1105 ✗ const expansion& r3 = expansion_product(a42, DeltaLambda2);
1106 ✗ const expansion& r4 = expansion_product(a43, DeltaLambda3);
1107 ✗ const expansion& r1234 = expansion_sum4(r1, r2, r3, r4);
1108 ✗ const expansion& r = expansion_diff(r0, r1234);
1109 ✗ Sign r_sign = r.sign();
1110
1111 // Simulation of Simplicity (symbolic perturbation)
1112 ✗ if(r_sign == ZERO) {
1113 ✗ stats_side4.log_SOS();
1114 ✗ const double* p_sort[5] = {p0, p1, p2, p3, p4};
1115 ✗ SOS_sort(p_sort, p_sort + 5, dim);
1116 ✗ for(index_t i = 0; i < 5; ++i) {
1117 ✗ if(p_sort[i] == p0) {
1118 ✗ const expansion& z1_0 = expansion_sum3(b01, b02, b03);
1119 ✗ const expansion& z1 = expansion_product(a30, z1_0);
1120 ✗ const expansion& z2_0 = expansion_sum3(b11, b12, b13);
1121 ✗ const expansion& z2 = expansion_product(a31, z2_0);
1122 ✗ const expansion& z3_0 = expansion_sum3(b21, b22, b23);
1123 ✗ const expansion& z3 = expansion_product(a32, z3_0);
1124 ✗ const expansion& z4_0 = expansion_sum3(b31, b32, b33);
1125 ✗ const expansion& z4 = expansion_product(a33, z4_0);
1126 ✗ const expansion& z1234 = expansion_sum4(z1, z2, z3, z4);
1127 ✗ const expansion& z = expansion_diff(Delta, z1234);
1128 ✗ Sign z_sign = z.sign();
1129 ✗ if(z_sign != ZERO) {
1130 ✗ return Sign(Delta_sign * z_sign);
1131 }
1132 ✗ } else if(p_sort[i] == p1) {
1133 ✗ const expansion& z1 = expansion_product(a30, b01);
1134 ✗ const expansion& z2 = expansion_product(a31, b11);
1135 ✗ const expansion& z3 = expansion_product(a32, b21);
1136 ✗ const expansion& z4 = expansion_product(a33, b31);
1137 ✗ const expansion& z = expansion_sum4(z1, z2, z3, z4);
1138 ✗ Sign z_sign = z.sign();
1139 ✗ if(z_sign != ZERO) {
1140 ✗ return Sign(Delta_sign * z_sign);
1141 }
1142 ✗ } else if(p_sort[i] == p2) {
1143 ✗ const expansion& z1 = expansion_product(a30, b02);
1144 ✗ const expansion& z2 = expansion_product(a31, b12);
1145 ✗ const expansion& z3 = expansion_product(a32, b22);
1146 ✗ const expansion& z4 = expansion_product(a33, b32);
1147 ✗ const expansion& z = expansion_sum4(z1, z2, z3, z4);
1148 ✗ Sign z_sign = z.sign();
1149 ✗ if(z_sign != ZERO) {
1150 ✗ return Sign(Delta_sign * z_sign);
1151 }
1152 ✗ } else if(p_sort[i] == p3) {
1153 ✗ const expansion& z1 = expansion_product(a30, b03);
1154 ✗ const expansion& z2 = expansion_product(a31, b13);
1155 ✗ const expansion& z3 = expansion_product(a32, b23);
1156 ✗ const expansion& z4 = expansion_product(a33, b33);
1157 ✗ const expansion& z = expansion_sum4(z1, z2, z3, z4);
1158 ✗ Sign z_sign = z.sign();
1159 ✗ if(z_sign != ZERO) {
1160 ✗ return Sign(Delta_sign * z_sign);
1161 }
1162 ✗ } else if(p_sort[i] == p4) {
1163 ✗ return NEGATIVE;
1164 }
1165 }
1166 ✗ geo_assert_not_reached;
1167 }
1168 ✗ return Sign(r_sign * Delta_sign);
1169 }
1170
1171 /**
1172 * \brief Implements side4() in 4d.
1173 */
1174 25769 Sign side4_4d_SOS(
1175 const double* p0,
1176 const double* p1, const double* p2, const double* p3, const double* p4,
1177 const double* q0, const double* q1, const double* q2, const double* q3
1178 ) {
1179 25769 Sign result = Sign(side4_4d_filter(p0, p1, p2, p3, p4, q0, q1, q2, q3));
1180
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 25769 times.
 25769 if(result == ZERO) {
1181 ✗ result = side4_exact_SOS(p0, p1, p2, p3, p4, q0, q1, q2, q3, 4);
1182 }
1183 25769 return result;
1184 }
1185
1186 /**
1187 * \brief Implements side4() in 6d.
1188 */
1189 20970 Sign side4_6d_SOS(
1190 const double* p0,
1191 const double* p1, const double* p2, const double* p3, const double* p4,
1192 const double* q0, const double* q1, const double* q2, const double* q3
1193 ) {
1194 20970 Sign result = Sign(side4_6d_filter(p0, p1, p2, p3, p4, q0, q1, q2, q3));
1195
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 20970 times.
 20970 if(result == ZERO) {
1196 ✗ result = side4_exact_SOS(p0, p1, p2, p3, p4, q0, q1, q2, q3, 6);
1197 }
1198 20970 return result;
1199 }
1200
1201 /**
1202 * \brief Implements side4() in 7d.
1203 */
1204 ✗ Sign side4_7d_SOS(
1205 const double* p0,
1206 const double* p1, const double* p2, const double* p3, const double* p4,
1207 const double* q0, const double* q1, const double* q2, const double* q3
1208 ) {
1209 ✗ Sign result = Sign(side4_7d_filter(p0, p1, p2, p3, p4, q0, q1, q2, q3));
1210 ✗ if(result == ZERO) {
1211 ✗ result = side4_exact_SOS(p0, p1, p2, p3, p4, q0, q1, q2, q3, 7);
1212 }
1213 ✗ return result;
1214 }
1215
1216 /**
1217 * \brief Implements side4() in 8d.
1218 */
1219 18812 Sign side4_8d_SOS(
1220 const double* p0,
1221 const double* p1, const double* p2, const double* p3, const double* p4,
1222 const double* q0, const double* q1, const double* q2, const double* q3
1223 ) {
1224 18812 Sign result = Sign(side4_8d_filter(p0, p1, p2, p3, p4, q0, q1, q2, q3));
1225
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 18812 times.
 18812 if(result == ZERO) {
1226 ✗ result = side4_exact_SOS(p0, p1, p2, p3, p4, q0, q1, q2, q3, 8);
1227 }
1228 18812 return result;
1229 }
1230
1231 // ============ orient2d ====================================================
1232
1233 1184972 Sign orient_2d_exact(const double* p0, const double* p1, const double* p2) {
1234 1184972 stats_orient2d.log_exact();
1235
3/8
✓ Branch 4 taken 1184972 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1184972 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1184972 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 1184972 const expansion& a11 = expansion_diff(p1[0], p0[0]);
1236
3/8
✓ Branch 4 taken 1184972 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1184972 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1184972 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 1184972 const expansion& a12 = expansion_diff(p1[1], p0[1]);
1237
3/8
✓ Branch 4 taken 1184972 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1184972 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1184972 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 1184972 const expansion& a21 = expansion_diff(p2[0], p0[0]);
1238
3/8
✓ Branch 4 taken 1184972 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1184972 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1184972 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 1184972 const expansion& a22 = expansion_diff(p2[1], p0[1]);
1239
2/6
✓ Branch 6 taken 1184972 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1184972 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 1184972 const expansion& Delta = expansion_det2x2(a11, a12, a21, a22);
1240 1184972 return Delta.sign();
1241 }
1242
1243 // ============ orient3d ===================================================
1244
1245 4931619 Sign orient_3d_exact(
1246 const double* p0, const double* p1, const double* p2, const double* p3
1247 ) {
1248 4931619 stats_orient3d.log_exact();
1249
1250
3/8
✓ Branch 4 taken 4931619 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 4931619 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 4931619 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 4931619 const expansion& a11 = expansion_diff(p1[0], p0[0]);
1251
3/8
✓ Branch 4 taken 4931619 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 4931619 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 4931619 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 4931619 const expansion& a12 = expansion_diff(p1[1], p0[1]);
1252
3/8
✓ Branch 4 taken 4931619 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 4931619 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 4931619 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 4931619 const expansion& a13 = expansion_diff(p1[2], p0[2]);
1253
1254
3/8
✓ Branch 4 taken 4931619 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 4931619 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 4931619 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 4931619 const expansion& a21 = expansion_diff(p2[0], p0[0]);
1255
3/8
✓ Branch 4 taken 4931619 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 4931619 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 4931619 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 4931619 const expansion& a22 = expansion_diff(p2[1], p0[1]);
1256
3/8
✓ Branch 4 taken 4931619 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 4931619 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 4931619 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 4931619 const expansion& a23 = expansion_diff(p2[2], p0[2]);
1257
1258
3/8
✓ Branch 4 taken 4931619 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 4931619 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 4931619 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 4931619 const expansion& a31 = expansion_diff(p3[0], p0[0]);
1259
3/8
✓ Branch 4 taken 4931619 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 4931619 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 4931619 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 4931619 const expansion& a32 = expansion_diff(p3[1], p0[1]);
1260
3/8
✓ Branch 4 taken 4931619 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 4931619 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 4931619 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 4931619 const expansion& a33 = expansion_diff(p3[2], p0[2]);
1261
1262
2/6
✓ Branch 6 taken 4931619 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 4931619 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 4931619 const expansion& Delta = expansion_det3x3(
1263 a11, a12, a13, a21, a22, a23, a31, a32, a33
1264 );
1265
1266 4931619 return Delta.sign();
1267 }
1268
1269 ✗ Sign side4h_3d_exact_SOS(
1270 const double* p0, const double* p1,
1271 const double* p2, const double* p3, const double* p4,
1272 double h0, double h1, double h2, double h3, double h4,
1273 bool sos = true
1274 ) {
1275 ✗ stats_orient3dh.log_exact();
1276
1277 ✗ const expansion& a11 = expansion_diff(p1[0], p0[0]);
1278 ✗ const expansion& a12 = expansion_diff(p1[1], p0[1]);
1279 ✗ const expansion& a13 = expansion_diff(p1[2], p0[2]);
1280 ✗ const expansion& a14 = expansion_diff(h0,h1);
1281
1282 ✗ const expansion& a21 = expansion_diff(p2[0], p0[0]);
1283 ✗ const expansion& a22 = expansion_diff(p2[1], p0[1]);
1284 ✗ const expansion& a23 = expansion_diff(p2[2], p0[2]);
1285 ✗ const expansion& a24 = expansion_diff(h0,h2);
1286
1287 ✗ const expansion& a31 = expansion_diff(p3[0], p0[0]);
1288 ✗ const expansion& a32 = expansion_diff(p3[1], p0[1]);
1289 ✗ const expansion& a33 = expansion_diff(p3[2], p0[2]);
1290 ✗ const expansion& a34 = expansion_diff(h0,h3);
1291
1292 ✗ const expansion& a41 = expansion_diff(p4[0], p0[0]);
1293 ✗ const expansion& a42 = expansion_diff(p4[1], p0[1]);
1294 ✗ const expansion& a43 = expansion_diff(p4[2], p0[2]);
1295 ✗ const expansion& a44 = expansion_diff(h0,h4);
1296
1297 // Note: we could probably reuse some of the 2x2 co-factors
1298 // (but for now I'd rather keep this form that is easier to
1299 // read ... and to debug if need be !)
1300 ✗ const expansion& Delta1 = expansion_det3x3(
1301 a21, a22, a23,
1302 a31, a32, a33,
1303 a41, a42, a43
1304 );
1305 ✗ const expansion& Delta2 = expansion_det3x3(
1306 a11, a12, a13,
1307 a31, a32, a33,
1308 a41, a42, a43
1309 );
1310 ✗ const expansion& Delta3 = expansion_det3x3(
1311 a11, a12, a13,
1312 a21, a22, a23,
1313 a41, a42, a43
1314 );
1315 ✗ const expansion& Delta4 = expansion_det3x3(
1316 a11, a12, a13,
1317 a21, a22, a23,
1318 a31, a32, a33
1319 );
1320
1321 ✗ Sign Delta4_sign = Delta4.sign();
1322 ✗ geo_assert(Delta4_sign != ZERO);
1323
1324 ✗ const expansion& r_1 = expansion_product(Delta1, a14);
1325 ✗ const expansion& r_2 = expansion_product(Delta2, a24).negate();
1326 ✗ const expansion& r_3 = expansion_product(Delta3, a34);
1327 ✗ const expansion& r_4 = expansion_product(Delta4, a44).negate();
1328 ✗ const expansion& r = expansion_sum4(r_1, r_2, r_3, r_4);
1329
1330 ✗ Sign r_sign = r.sign();
1331
1332 // Simulation of Simplicity (symbolic perturbation)
1333 ✗ if(sos && r_sign == ZERO) {
1334 ✗ stats_orient3dh.log_SOS();
1335 ✗ const double* p_sort[5] = {p0, p1, p2, p3, p4};
1336 ✗ SOS_sort(p_sort, p_sort + 5, 3);
1337 ✗ for(index_t i = 0; i < 5; ++i) {
1338 ✗ if(p_sort[i] == p0) {
1339 ✗ const expansion& z1 = expansion_diff(Delta2, Delta1);
1340 ✗ const expansion& z2 = expansion_diff(Delta4, Delta3);
1341 ✗ const expansion& z = expansion_sum(z1, z2);
1342 ✗ Sign z_sign = z.sign();
1343 ✗ if(z_sign != ZERO) {
1344 ✗ return Sign(Delta4_sign * z_sign);
1345 }
1346 ✗ } else if(p_sort[i] == p1) {
1347 ✗ Sign Delta1_sign = Delta1.sign();
1348 ✗ if(Delta1_sign != ZERO) {
1349 ✗ return Sign(Delta4_sign * Delta1_sign);
1350 }
1351 ✗ } else if(p_sort[i] == p2) {
1352 ✗ Sign Delta2_sign = Delta2.sign();
1353 ✗ if(Delta2_sign != ZERO) {
1354 ✗ return Sign(-Delta4_sign * Delta2_sign);
1355 }
1356 ✗ } else if(p_sort[i] == p3) {
1357 ✗ Sign Delta3_sign = Delta3.sign();
1358 ✗ if(Delta3_sign != ZERO) {
1359 ✗ return Sign(Delta4_sign * Delta3_sign);
1360 }
1361 ✗ } else if(p_sort[i] == p4) {
1362 ✗ return NEGATIVE;
1363 }
1364 }
1365 }
1366 ✗ return Sign(Delta4_sign * r_sign);
1367 }
1368
1369
1370 ✗ Sign side3h_2d_exact_SOS(
1371 const double* p0, const double* p1,
1372 const double* p2, const double* p3,
1373 double h0, double h1, double h2, double h3,
1374 bool sos = true
1375 ) {
1376
1377 ✗ const expansion& a11 = expansion_diff(p1[0], p0[0]);
1378 ✗ const expansion& a12 = expansion_diff(p1[1], p0[1]);
1379 ✗ const expansion& a13 = expansion_diff(h0,h1);
1380
1381 ✗ const expansion& a21 = expansion_diff(p2[0], p0[0]);
1382 ✗ const expansion& a22 = expansion_diff(p2[1], p0[1]);
1383 ✗ const expansion& a23 = expansion_diff(h0,h2);
1384
1385 ✗ const expansion& a31 = expansion_diff(p3[0], p0[0]);
1386 ✗ const expansion& a32 = expansion_diff(p3[1], p0[1]);
1387 ✗ const expansion& a33 = expansion_diff(h0,h3);
1388
1389 ✗ const expansion& Delta1 = expansion_det2x2(
1390 a21, a22,
1391 a31, a32
1392 );
1393 ✗ const expansion& Delta2 = expansion_det2x2(
1394 a11, a12,
1395 a31, a32
1396 );
1397 ✗ const expansion& Delta3 = expansion_det2x2(
1398 a11, a12,
1399 a21, a22
1400 );
1401
1402 ✗ Sign Delta3_sign = Delta3.sign();
1403 ✗ geo_assert(Delta3_sign != ZERO);
1404
1405 ✗ const expansion& r_1 = expansion_product(Delta1, a13);
1406 ✗ const expansion& r_2 = expansion_product(Delta2, a23).negate();
1407 ✗ const expansion& r_3 = expansion_product(Delta3, a33);
1408 ✗ const expansion& r = expansion_sum3(r_1, r_2, r_3);
1409
1410 ✗ Sign r_sign = r.sign();
1411
1412 // Simulation of Simplicity (symbolic perturbation)
1413 ✗ if(sos && r_sign == ZERO) {
1414 ✗ const double* p_sort[4] = {p0, p1, p2, p3};
1415 ✗ SOS_sort(p_sort, p_sort + 4, 2);
1416 ✗ for(index_t i = 0; i < 4; ++i) {
1417 ✗ if(p_sort[i] == p0) {
1418 ✗ const expansion& z1 = expansion_diff(Delta2, Delta1);
1419 ✗ const expansion& z = expansion_sum(z1, Delta3);
1420 ✗ Sign z_sign = z.sign();
1421 ✗ if(z_sign != ZERO) {
1422 ✗ return Sign(Delta3_sign * z_sign);
1423 }
1424 ✗ } else if(p_sort[i] == p1) {
1425 ✗ Sign Delta1_sign = Delta1.sign();
1426 ✗ if(Delta1_sign != ZERO) {
1427 ✗ return Sign(Delta3_sign * Delta1_sign);
1428 }
1429 ✗ } else if(p_sort[i] == p2) {
1430 ✗ Sign Delta2_sign = Delta2.sign();
1431 ✗ if(Delta2_sign != ZERO) {
1432 ✗ return Sign(-Delta3_sign * Delta2_sign);
1433 }
1434 ✗ } else if(p_sort[i] == p3) {
1435 ✗ return NEGATIVE;
1436 }
1437 }
1438 }
1439 ✗ return Sign(Delta3_sign * r_sign);
1440 }
1441
1442
1443 // ================================ det and dot =======================
1444
1445 /**
1446 * \brief Computes the sign of the determinant of a 3x3
1447 * matrix formed by three 3d points using exact arithmetics.
1448 * \param[in] p0 , p1 , p2 the three points
1449 * \return the sign of the determinant of the matrix.
1450 */
1451 ✗ Sign det_3d_exact(
1452 const double* p0, const double* p1, const double* p2
1453 ) {
1454 ✗ stats_det3d.log_exact();
1455
1456 ✗ const expansion& p0_0 = expansion_create(p0[0]);
1457 ✗ const expansion& p0_1 = expansion_create(p0[1]);
1458 ✗ const expansion& p0_2 = expansion_create(p0[2]);
1459
1460 ✗ const expansion& p1_0 = expansion_create(p1[0]);
1461 ✗ const expansion& p1_1 = expansion_create(p1[1]);
1462 ✗ const expansion& p1_2 = expansion_create(p1[2]);
1463
1464 ✗ const expansion& p2_0 = expansion_create(p2[0]);
1465 ✗ const expansion& p2_1 = expansion_create(p2[1]);
1466 ✗ const expansion& p2_2 = expansion_create(p2[2]);
1467
1468 ✗ const expansion& Delta = expansion_det3x3(
1469 p0_0, p0_1, p0_2,
1470 p1_0, p1_1, p1_2,
1471 p2_0, p2_1, p2_2
1472 );
1473
1474 ✗ return Delta.sign();
1475 }
1476
1477
1478 /**
1479 * \brief Tests whether three points are aligned using
1480 * exact arithmetics.
1481 * \param[in] p0 , p1 , p2 the three points
1482 * \retval true if the three points are aligned.
1483 * \retval false otherwise.
1484 */
1485 2092743 bool aligned_3d_exact(
1486 const double* p0, const double* p1, const double* p2
1487 ) {
1488
3/8
✓ Branch 4 taken 2092743 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2092743 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2092743 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 2092743 const expansion& U_0 = expansion_diff(p1[0],p0[0]);
1489
3/8
✓ Branch 4 taken 2092743 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2092743 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2092743 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 2092743 const expansion& U_1 = expansion_diff(p1[1],p0[1]);
1490
3/8
✓ Branch 4 taken 2092743 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2092743 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2092743 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 2092743 const expansion& U_2 = expansion_diff(p1[2],p0[2]);
1491
1492
3/8
✓ Branch 4 taken 2092743 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2092743 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2092743 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 2092743 const expansion& V_0 = expansion_diff(p2[0],p0[0]);
1493
3/8
✓ Branch 4 taken 2092743 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2092743 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2092743 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 2092743 const expansion& V_1 = expansion_diff(p2[1],p0[1]);
1494
3/8
✓ Branch 4 taken 2092743 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2092743 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 2092743 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 2092743 const expansion& V_2 = expansion_diff(p2[2],p0[2]);
1495
1496
2/6
✓ Branch 6 taken 2092743 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2092743 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 2092743 const expansion& N_0 = expansion_det2x2(U_1, V_1, U_2, V_2);
1497
2/6
✓ Branch 6 taken 2092743 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2092743 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 2092743 const expansion& N_1 = expansion_det2x2(U_2, V_2, U_0, V_0);
1498
2/6
✓ Branch 6 taken 2092743 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2092743 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 2092743 const expansion& N_2 = expansion_det2x2(U_0, V_0, U_1, V_1);
1499
1500 return(
1501 2092743 N_0.sign() == 0 &&
1502
4/4
✓ Branch 0 taken 386069 times.
✓ Branch 1 taken 1706674 times.
✓ Branch 3 taken 317610 times.
✓ Branch 4 taken 68459 times.
 2410353 N_1.sign() == 0 &&
1503
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 317610 times.
 317610 N_2.sign() == 0
1504 2092743 );
1505 }
1506
1507 /**
1508 * \brief Computes the sign of the dot product between two
1509 * vectors using exact arithmetics.
1510 * \param[in] p0 , p1 , p2 three 3d points.
1511 * \return the sign of the dot product between the vectors
1512 * p0p1 and p0p2.
1513 */
1514 145176 Sign dot_3d_exact(
1515 const double* p0, const double* p1, const double* p2
1516 ) {
1517
3/8
✓ Branch 4 taken 145176 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 145176 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 145176 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 145176 const expansion& U_0 = expansion_diff(p1[0],p0[0]);
1518
3/8
✓ Branch 4 taken 145176 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 145176 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 145176 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 145176 const expansion& U_1 = expansion_diff(p1[1],p0[1]);
1519
3/8
✓ Branch 4 taken 145176 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 145176 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 145176 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 145176 const expansion& U_2 = expansion_diff(p1[2],p0[2]);
1520
1521
3/8
✓ Branch 4 taken 145176 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 145176 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 145176 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 145176 const expansion& V_0 = expansion_diff(p2[0],p0[0]);
1522
3/8
✓ Branch 4 taken 145176 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 145176 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 145176 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 145176 const expansion& V_1 = expansion_diff(p2[1],p0[1]);
1523
3/8
✓ Branch 4 taken 145176 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 145176 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 145176 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 145176 const expansion& V_2 = expansion_diff(p2[2],p0[2]);
1524
1525
2/6
✓ Branch 6 taken 145176 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 145176 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 145176 const expansion& UV_0 = expansion_product(U_0, V_0);
1526
2/6
✓ Branch 6 taken 145176 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 145176 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 145176 const expansion& UV_1 = expansion_product(U_1, V_1);
1527
2/6
✓ Branch 6 taken 145176 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 145176 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 145176 const expansion& UV_2 = expansion_product(U_2, V_2);
1528
1529
2/6
✓ Branch 6 taken 145176 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 145176 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 145176 const expansion& Delta = expansion_sum3(UV_0, UV_1, UV_2);
1530
1531 145176 return Delta.sign();
1532 }
1533
1534 /**
1535 * \brief Compares two dot products using exact arithmetics.
1536 * \param[in] v0 , v1 , v2 three vectors
1537 * \return the sign of v0.v1 - v0.v2
1538 */
1539 ✗ Sign dot_compare_3d_exact(
1540 const double* v0, const double* v1, const double* v2
1541 ) {
1542 ✗ const expansion& d01_0 = expansion_product(v0[0], v1[0]);
1543 ✗ const expansion& d01_1 = expansion_product(v0[1], v1[1]);
1544 ✗ const expansion& d01_2 = expansion_product(v0[2], v1[2]);
1545 ✗ const expansion& d01_12 = expansion_sum(d01_1, d01_2);
1546 ✗ const expansion& d01 = expansion_sum(d01_0, d01_12);
1547
1548 ✗ const expansion& d02_0 = expansion_product(v0[0], v2[0]);
1549 ✗ const expansion& d02_1 = expansion_product(v0[1], v2[1]);
1550 ✗ const expansion& d02_2 = expansion_product(v0[2], v2[2]);
1551 ✗ const expansion& d02_12 = expansion_sum(d02_1, d02_2);
1552 ✗ const expansion& d02 = expansion_sum(d02_0, d02_12);
1553
1554 ✗ const expansion& result = expansion_diff(d01, d02);
1555
1556 ✗ return result.sign();
1557 }
1558 }
1559
1560 /****************************************************************************/
1561
1562 namespace GEO {
1563
1564 namespace PCK {
1565
1566 122 void set_SOS_mode(SOSMode m) {
1567 122 SOS_mode_ = m;
1568 122 }
1569
1570 60 SOSMode get_SOS_mode() {
1571 60 return SOS_mode_;
1572 }
1573
1574
1575 4082607 Sign side1_SOS(
1576 const double* p0, const double* p1,
1577 const double* q0,
1578 coord_index_t DIM
1579 ) {
1580 4082607 stats_side1.log_invoke();
1581
4/6
✓ Branch 0 taken 3914542 times.
✓ Branch 1 taken 50792 times.
✓ Branch 2 taken 59448 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 57825 times.
✗ Branch 5 not taken.
 4082607 switch(DIM) {
1582 3914542 case 3:
1583 3914542 return side1_3d_SOS(p0, p1, q0);
1584 50792 case 4:
1585 50792 return side1_4d_SOS(p0, p1, q0);
1586 59448 case 6:
1587 59448 return side1_6d_SOS(p0, p1, q0);
1588 ✗ case 7:
1589 ✗ return side1_7d_SOS(p0, p1, q0);
1590 57825 case 8:
1591 57825 return side1_8d_SOS(p0, p1, q0);
1592 }
1593 ✗ geo_assert_not_reached;
1594 }
1595
1596 5798279 Sign side2_SOS(
1597 const double* p0, const double* p1, const double* p2,
1598 const double* q0, const double* q1,
1599 coord_index_t DIM
1600 ) {
1601 5798279 stats_side2.log_invoke();
1602
4/6
✓ Branch 0 taken 5090168 times.
✓ Branch 1 taken 189360 times.
✓ Branch 2 taken 260704 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 258047 times.
✗ Branch 5 not taken.
 5798279 switch(DIM) {
1603 5090168 case 3:
1604 5090168 return side2_3d_SOS(p0, p1, p2, q0, q1);
1605 189360 case 4:
1606 189360 return side2_4d_SOS(p0, p1, p2, q0, q1);
1607 260704 case 6:
1608 260704 return side2_6d_SOS(p0, p1, p2, q0, q1);
1609 ✗ case 7:
1610 ✗ return side2_7d_SOS(p0, p1, p2, q0, q1);
1611 258047 case 8:
1612 258047 return side2_8d_SOS(p0, p1, p2, q0, q1);
1613 }
1614 ✗ geo_assert_not_reached;
1615 }
1616
1617 11483641 Sign side3_SOS(
1618 const double* p0, const double* p1,
1619 const double* p2, const double* p3,
1620 const double* q0, const double* q1, const double* q2,
1621 coord_index_t DIM
1622 ) {
1623 11483641 stats_side3.log_invoke();
1624
4/6
✓ Branch 0 taken 11013534 times.
✓ Branch 1 taken 149639 times.
✓ Branch 2 taken 165662 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 154806 times.
✗ Branch 5 not taken.
 11483641 switch(DIM) {
1625 11013534 case 3:
1626 11013534 return side3_3d_SOS(p0, p1, p2, p3, q0, q1, q2);
1627 149639 case 4:
1628 149639 return side3_4d_SOS(p0, p1, p2, p3, q0, q1, q2);
1629 165662 case 6:
1630 165662 return side3_6d_SOS(p0, p1, p2, p3, q0, q1, q2);
1631 ✗ case 7:
1632 ✗ return side3_7d_SOS(p0, p1, p2, p3, q0, q1, q2);
1633 154806 case 8:
1634 154806 return side3_8d_SOS(p0, p1, p2, p3, q0, q1, q2);
1635 }
1636 ✗ geo_assert_not_reached;
1637 }
1638
1639
1640 ✗ Sign side3_3dlifted_SOS(
1641 const double* p0, const double* p1,
1642 const double* p2, const double* p3,
1643 double h0, double h1, double h2, double h3,
1644 const double* q0, const double* q1, const double* q2,
1645 bool SOS
1646 ) {
1647 Sign result = Sign(
1648 ✗ side3h_3d_filter(p0, p1, p2, p3, h0, h1, h2, h3, q0, q1, q2)
1649 );
1650 ✗ if(SOS && result == ZERO) {
1651 ✗ result = side3h_exact_SOS(
1652 p0, p1, p2, p3, h0, h1, h2, h3, q0, q1, q2
1653 );
1654 }
1655 ✗ return result;
1656 }
1657
1658 65551 Sign side4_SOS(
1659 const double* p0,
1660 const double* p1, const double* p2,
1661 const double* p3, const double* p4,
1662 const double* q0, const double* q1,
1663 const double* q2, const double* q3,
1664 coord_index_t DIM
1665 ) {
1666
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 25769 times.
✓ Branch 2 taken 20970 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 18812 times.
✗ Branch 5 not taken.
 65551 switch(DIM) {
1667 ✗ case 3:
1668 // 3d is a special case for side4()
1669 // (intrinsic dim == ambient dim)
1670 // therefore embedding tet q0,q1,q2,q3 is not needed.
1671 // WARNING: cnt_side4_total is not incremented here,
1672 // because it is
1673 // incremented in side4_3d_SOS().
1674 ✗ return side4_3d_SOS(p0, p1, p2, p3, p4);
1675 25769 case 4:
1676 25769 stats_side4.log_invoke();
1677 25769 return side4_4d_SOS(p0, p1, p2, p3, p4, q0, q1, q2, q3);
1678 20970 case 6:
1679 20970 stats_side4.log_invoke();
1680 20970 return side4_6d_SOS(p0, p1, p2, p3, p4, q0, q1, q2, q3);
1681 ✗ case 7:
1682 ✗ stats_side4.log_invoke();
1683 ✗ return side4_7d_SOS(p0, p1, p2, p3, p4, q0, q1, q2, q3);
1684 18812 case 8:
1685 18812 stats_side4.log_invoke();
1686 18812 return side4_8d_SOS(p0, p1, p2, p3, p4, q0, q1, q2, q3);
1687 }
1688 ✗ geo_assert_not_reached;
1689 }
1690
1691
1692 ✗ Sign side4_3d(
1693 const double* p0, const double* p1, const double* p2,
1694 const double* p3, const double* p4
1695 ) {
1696 ✗ stats_side4.log_invoke();
1697 ✗ Sign result = Sign(side4_3d_filter(p0, p1, p2, p3, p4));
1698 ✗ if(result == 0) {
1699 // last argument is false: do not apply symbolic perturbation
1700 ✗ result = side4_3d_exact_SOS(p0, p1, p2, p3, p4, false);
1701 }
1702 ✗ return result;
1703 }
1704
1705 30218722 Sign side4_3d_SOS(
1706 const double* p0, const double* p1,
1707 const double* p2, const double* p3,
1708 const double* p4
1709 ) {
1710 30218722 stats_side4.log_invoke();
1711 30218722 Sign result = Sign(side4_3d_filter(p0, p1, p2, p3, p4));
1712
2/2
✓ Branch 0 taken 98355 times.
✓ Branch 1 taken 30120367 times.
 30218722 if(result == 0) {
1713 98355 result = side4_3d_exact_SOS(p0, p1, p2, p3, p4);
1714 }
1715 30218722 return result;
1716 }
1717
1718
1719 631325 Sign in_sphere_3d_SOS(
1720 const double* p0, const double* p1,
1721 const double* p2, const double* p3,
1722 const double* p4
1723 ) {
1724 // in_sphere_3d is simply implemented using side4_3d.
1725 // Both predicates are equivalent through duality as can
1726 // be easily seen:
1727 // side4_3d(p0,p1,p2,p3,p4) returns POSITIVE if
1728 // d(q,p0) < d(q,p4)
1729 // where q denotes the circumcenter of (p0,p1,p2,p3)
1730 // Note that d(q,p0) = R (radius of circumscribed sphere)
1731 // In other words, side4_3d(p0,p1,p2,p3,p4) returns POSITIVE if
1732 // d(q,p4) > R which means whenever p4 is not in the
1733 // circumscribed sphere of (p0,p1,p2,p3).
1734 // Therefore:
1735 // in_sphere_3d(p0,p1,p2,p3,p4) = -side4_3d(p0,p1,p2,p3,p4)
1736
1737 631325 stats_side4.log_invoke();
1738
1739 // This specialized filter supposes that orient_3d(p0,p1,p2,p3) > 0
1740
1741 631325 Sign result = Sign(in_sphere_3d_filter_optim(p0, p1, p2, p3, p4));
1742
1743
2/2
✓ Branch 0 taken 154037 times.
✓ Branch 1 taken 477288 times.
 631325 if(result == 0) {
1744 154037 result = side4_3d_exact_SOS(p0, p1, p2, p3, p4);
1745 }
1746 631325 return Sign(-result);
1747 }
1748
1749 936185 Sign GEOGRAM_API in_circle_2d_SOS(
1750 const double* p0, const double* p1, const double* p2,
1751 const double* p3
1752 ) {
1753 // in_circle_2d is simply implemented using side3_2d.
1754 // Both predicates are equivalent through duality as can
1755 // be easily seen:
1756 // side3_2d(p0,p1,p2,p3,p0,p1,p2) returns POSITIVE if
1757 // d(q,p0) < d(q,p3)
1758 // where q denotes the circumcenter of (p0,p1,p2)
1759 // Note that d(q,p0) = R (radius of circumscribed circle)
1760 // In other words, side3_2d(p0,p1,p2,p3,p4) returns POSITIVE if
1761 // d(q,p3) > R which means whenever p3 is not in the
1762 // circumscribed circle of (p0,p1,p2).
1763 // Therefore:
1764 // in_circle_2d(p0,p1,p2,p3) = -side3_2d(p0,p1,p2,p3)
1765
1766 // TODO: implement specialized filter like the one used
1767 // by "in-sphere".
1768 936185 Sign s = Sign(-side3_2d_filter(p0, p1, p2, p3, p0, p1, p2));
1769
2/2
✓ Branch 0 taken 936157 times.
✓ Branch 1 taken 28 times.
 936185 if(s != ZERO) {
1770 936157 return s;
1771 }
1772 28 return Sign(-side3_exact_SOS(p0, p1, p2, p3, p0, p1, p2, 2));
1773 }
1774
1775 40238 Sign GEOGRAM_API in_circle_3d_SOS(
1776 const double* p0, const double* p1, const double* p2,
1777 const double* p3
1778 ) {
1779 // in_circle_3d is simply implemented using side3_3d.
1780 // Both predicates are equivalent through duality as can
1781 // be easily seen:
1782 // side3_3d(p0,p1,p2,p3,p0,p1,p2) returns POSITIVE if
1783 // d(q,p0) < d(q,p3)
1784 // where q denotes the circumcenter of (p0,p1,p2)
1785 // Note that d(q,p0) = R (radius of circumscribed circle)
1786 // In other words, side3_3d(p0,p1,p2,p3,p4) returns POSITIVE if
1787 // d(q,p3) > R which means whenever p3 is not in the
1788 // circumscribed circle of (p0,p1,p2).
1789 // Therefore:
1790 // in_circle_3d(p0,p1,p2,p3) = -side3_3d(p0,p1,p2,p3)
1791 40238 return Sign(-side3_3d_SOS(p0,p1,p2,p3,p0,p1,p2));
1792 }
1793
1794 ✗ Sign GEOGRAM_API in_circle_3dlifted_SOS(
1795 const double* p0, const double* p1, const double* p2,
1796 const double* p3,
1797 double h0, double h1, double h2, double h3,
1798 bool SOS
1799 ) {
1800 // in_circle_3dlifted is simply implemented using side3_3dlifted.
1801 // Both predicates are equivalent through duality
1802 // (see comment in in_circle_3d_SOS(), the same
1803 // remark applies).
1804 return Sign(
1805 ✗ -side3_3dlifted_SOS(p0,p1,p2,p3,h0,h1,h2,h3,p0,p1,p2,SOS)
1806 ✗ );
1807 }
1808
1809
1810 19283551 Sign orient_2d(
1811 const double* p0, const double* p1, const double* p2
1812 ) {
1813 19283551 stats_orient2d.log_invoke();
1814 19283551 Sign result = Sign(orient_2d_filter(p0, p1, p2));
1815
2/2
✓ Branch 0 taken 1184972 times.
✓ Branch 1 taken 18098579 times.
 19283551 if(result == 0) {
1816 1184972 result = orient_2d_exact(p0, p1, p2);
1817 }
1818 19283551 return result;
1819 }
1820
1821 ✗ Sign orient_2dlifted_SOS(
1822 const double* p0, const double* p1,
1823 const double* p2, const double* p3,
1824 double h0, double h1, double h2, double h3
1825 ) {
1826 Sign result = Sign(
1827 ✗ side3_2dlifted_2d_filter(
1828 p0, p1, p2, p3, h0, h1, h2, h3
1829 )
1830 );
1831 ✗ if(result == 0) {
1832 ✗ result = side3h_2d_exact_SOS(
1833 p0, p1, p2, p3, h0, h1, h2, h3
1834 );
1835 }
1836 // orient_3d() is opposite to side3h()
1837 // (like in_sphere() that is opposite to side3())
1838 ✗ return result;
1839 }
1840
1841 19048953 Sign orient_3d(
1842 const double* p0, const double* p1,
1843 const double* p2, const double* p3
1844 ) {
1845 19048953 stats_orient3d.log_invoke();
1846 19048953 Sign result = Sign(orient_3d_filter(p0, p1, p2, p3));
1847
2/2
✓ Branch 0 taken 4931619 times.
✓ Branch 1 taken 14117334 times.
 19048953 if(result == 0) {
1848 4931619 result = orient_3d_exact(p0, p1, p2, p3);
1849 }
1850 19048953 return result;
1851 }
1852
1853 6044649 Sign orient_3d_SOS(
1854 const double* p0, const double* p1,
1855 const double* p2, const double* p3
1856 ) {
1857 struct SOS {
1858 555365 SOS(
1859 const double* p0, const double* p1,
1860 const double* p2, const double* p3
1861 555365 ) : p_orig{p0,p1,p2,p3}, p_sort{p0, p1, p2, p3} {
1862 555365 SOS_sort(p_sort, p_sort+4, 3);
1863 555365 parity = Permutation::permutation_is_odd(p_orig, p_sort, 4)
1864
2/2
✓ Branch 0 taken 277923 times.
✓ Branch 1 taken 277442 times.
 555365 ? NEGATIVE : POSITIVE;
1865 555365 }
1866 77530 Sign orient_1d(index_t i, index_t j, index_t ax) const {
1867 77530 return Sign(parity * geo_cmp(p_sort[i][ax], p_sort[j][ax]));
1868 }
1869 1177590 Sign orient_2d(
1870 index_t i, index_t j, index_t k, index_t ax1, index_t ax2
1871 ) const {
1872 1177590 double x0 = p_sort[i][ax1]; double y0 = p_sort[i][ax2];
1873 1177590 double x1 = p_sort[j][ax1]; double y1 = p_sort[j][ax2];
1874 1177590 double x2 = p_sort[k][ax1]; double y2 = p_sort[k][ax2];
1875
3/8
✓ Branch 4 taken 1177590 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1177590 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1177590 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 1177590 const expansion& a11 = expansion_diff(x1, x0);
1876
3/8
✓ Branch 4 taken 1177590 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1177590 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1177590 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 1177590 const expansion& a12 = expansion_diff(y1, y0);
1877
3/8
✓ Branch 4 taken 1177590 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1177590 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1177590 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 1177590 const expansion& a21 = expansion_diff(x2, x0);
1878
3/8
✓ Branch 4 taken 1177590 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1177590 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1177590 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 1177590 const expansion& a22 = expansion_diff(y2, y0);
1879
2/6
✓ Branch 6 taken 1177590 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1177590 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 1177590 const expansion& D = expansion_det2x2(a11, a12, a21, a22);
1880 1177590 return Sign(parity * D.sign());
1881 }
1882 const double* p_orig[4];
1883 const double* p_sort[4];
1884 Sign parity;
1885 };
1886
1887 6044649 return orient_3d_SOS_impl<const double*, SOS>(p0,p1,p2,p3);
1888 }
1889
1890 ✗ Sign orient_3dlifted(
1891 const double* p0, const double* p1,
1892 const double* p2, const double* p3, const double* p4,
1893 double h0, double h1, double h2, double h3, double h4
1894 ) {
1895 ✗ stats_orient3dh.log_invoke();
1896 Sign result = Sign(
1897 ✗ side4h_3d_filter(
1898 p0, p1, p2, p3, p4, h0, h1, h2, h3, h4
1899 )
1900 );
1901 ✗ if(result == 0) {
1902 // last argument is false -> do not perturb.
1903 ✗ result = side4h_3d_exact_SOS(
1904 p0, p1, p2, p3, p4, h0, h1, h2, h3, h4, false
1905 );
1906 }
1907 // orient_4d() is opposite to side4h()
1908 // (like in_sphere() that is opposite to side4())
1909 ✗ return Sign(-result);
1910 }
1911
1912 ✗ Sign orient_3dlifted_SOS(
1913 const double* p0, const double* p1,
1914 const double* p2, const double* p3, const double* p4,
1915 double h0, double h1, double h2, double h3, double h4
1916 ) {
1917 ✗ stats_orient3dh.log_invoke();
1918 Sign result = Sign(
1919 ✗ side4h_3d_filter(
1920 p0, p1, p2, p3, p4, h0, h1, h2, h3, h4
1921 )
1922 );
1923 ✗ if(result == 0) {
1924 ✗ result = side4h_3d_exact_SOS(
1925 p0, p1, p2, p3, p4, h0, h1, h2, h3, h4
1926 );
1927 }
1928 // orient_4d() is opposite to side4h()
1929 // (like in_sphere() that is opposite to side4())
1930 ✗ return Sign(-result);
1931 }
1932
1933 ✗ Sign det_3d(
1934 const double* p0, const double* p1, const double* p2
1935 ) {
1936 ✗ stats_det3d.log_invoke();
1937 Sign result = Sign(
1938 ✗ det_3d_filter(p0, p1, p2)
1939 );
1940 ✗ if(result == 0) {
1941 ✗ result = det_3d_exact(p0, p1, p2);
1942 }
1943 ✗ return result;
1944 }
1945
1946
1947 ✗ Sign det_4d(
1948 const double* p0, const double* p1,
1949 const double* p2, const double* p3
1950 ) {
1951 ✗ stats_det4d.log_invoke();
1952 Sign result = Sign(
1953 ✗ det_4d_filter(p0, p1, p2, p3)
1954 );
1955
1956 ✗ if(result == 0) {
1957 ✗ stats_det4d.log_exact();
1958
1959 ✗ const expansion& p0_0 = expansion_create(p0[0]);
1960 ✗ const expansion& p0_1 = expansion_create(p0[1]);
1961 ✗ const expansion& p0_2 = expansion_create(p0[2]);
1962 ✗ const expansion& p0_3 = expansion_create(p0[3]);
1963
1964 ✗ const expansion& p1_0 = expansion_create(p1[0]);
1965 ✗ const expansion& p1_1 = expansion_create(p1[1]);
1966 ✗ const expansion& p1_2 = expansion_create(p1[2]);
1967 ✗ const expansion& p1_3 = expansion_create(p1[3]);
1968
1969 ✗ const expansion& p2_0 = expansion_create(p2[0]);
1970 ✗ const expansion& p2_1 = expansion_create(p2[1]);
1971 ✗ const expansion& p2_2 = expansion_create(p2[2]);
1972 ✗ const expansion& p2_3 = expansion_create(p2[3]);
1973
1974 ✗ const expansion& p3_0 = expansion_create(p3[0]);
1975 ✗ const expansion& p3_1 = expansion_create(p3[1]);
1976 ✗ const expansion& p3_2 = expansion_create(p3[2]);
1977 ✗ const expansion& p3_3 = expansion_create(p3[3]);
1978
1979 ✗ result = sign_of_expansion_determinant(
1980 p0_0, p0_1, p0_2, p0_3,
1981 p1_0, p1_1, p1_2, p1_3,
1982 p2_0, p2_1, p2_2, p2_3,
1983 p3_0, p3_1, p3_2, p3_3
1984 );
1985 }
1986 ✗ return result;
1987 }
1988
1989
1990 ✗ Sign det_compare_4d(
1991 const double* p0, const double* p1,
1992 const double* p2, const double* p3,
1993 const double* p4
1994 ) {
1995 Sign result = Sign(
1996 ✗ det_compare_4d_filter(p0, p1, p2, p3, p4)
1997 );
1998 ✗ if(result == 0) {
1999 ✗ const expansion& p0_0 = expansion_create(p0[0]);
2000 ✗ const expansion& p0_1 = expansion_create(p0[1]);
2001 ✗ const expansion& p0_2 = expansion_create(p0[2]);
2002 ✗ const expansion& p0_3 = expansion_create(p0[3]);
2003
2004 ✗ const expansion& p1_0 = expansion_create(p1[0]);
2005 ✗ const expansion& p1_1 = expansion_create(p1[1]);
2006 ✗ const expansion& p1_2 = expansion_create(p1[2]);
2007 ✗ const expansion& p1_3 = expansion_create(p1[3]);
2008
2009 ✗ const expansion& p2_0 = expansion_create(p2[0]);
2010 ✗ const expansion& p2_1 = expansion_create(p2[1]);
2011 ✗ const expansion& p2_2 = expansion_create(p2[2]);
2012 ✗ const expansion& p2_3 = expansion_create(p2[3]);
2013
2014 ✗ const expansion& a3_0 = expansion_diff(p4[0],p3[0]);
2015 ✗ const expansion& a3_1 = expansion_diff(p4[1],p3[1]);
2016 ✗ const expansion& a3_2 = expansion_diff(p4[2],p3[2]);
2017 ✗ const expansion& a3_3 = expansion_diff(p4[3],p3[3]);
2018
2019 ✗ result = sign_of_expansion_determinant(
2020 p0_0, p0_1, p0_2, p0_3,
2021 p1_0, p1_1, p1_2, p1_3,
2022 p2_0, p2_1, p2_2, p2_3,
2023 a3_0, a3_1, a3_2, a3_3
2024 );
2025 }
2026 ✗ return result;
2027 }
2028
2029
2030 2092743 bool aligned_3d(
2031 const double* p0, const double* p1, const double* p2
2032 ) {
2033 /*
2034 Sign result = Sign(
2035 aligned_3d_filter(p0,p1,p2)
2036 );
2037 if(result != 0) {
2038 return false;
2039 }
2040 */
2041 2092743 return aligned_3d_exact(p0, p1, p2);
2042 }
2043
2044 145176 Sign dot_3d(
2045 const double* p0, const double* p1, const double* p2
2046 ) {
2047 145176 Sign result = Sign(det_3d_filter(p0, p1, p2));
2048
1/2
✓ Branch 0 taken 145176 times.
✗ Branch 1 not taken.
 145176 if(result == 0) {
2049 145176 result = dot_3d_exact(p0, p1, p2);
2050 }
2051 145176 return result;
2052 }
2053
2054 ✗ Sign dot_compare_3d(
2055 const double* v0, const double* v1, const double* v2
2056 ) {
2057 ✗ Sign result = Sign(dot_compare_3d_filter(v0, v1, v2));
2058 ✗ if(result == 0) {
2059 ✗ result = dot_compare_3d_exact(v0, v1, v2);
2060 }
2061 ✗ return result;
2062 }
2063
2064
2065 3 bool points_are_identical_2d(
2066 const double* p1,
2067 const double* p2
2068 ) {
2069 return
2070
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 6 (p1[0] == p2[0]) &&
2071
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
 6 (p1[1] == p2[1])
2072 ;
2073 }
2074
2075 7 bool points_are_identical_3d(
2076 const double* p1,
2077 const double* p2
2078 ) {
2079 return
2080 13 (p1[0] == p2[0]) &&
2081
4/4
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 3 times.
 10 (p1[1] == p2[1]) &&
2082
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
 10 (p1[2] == p2[2])
2083 ;
2084 }
2085
2086 7 bool points_are_colinear_3d(
2087 const double* p1,
2088 const double* p2,
2089 const double* p3
2090 ) {
2091 // Colinearity is tested by using four coplanarity
2092 // tests with four points that are not coplanar.
2093 // TODO: use PCK::aligned_3d() instead (to be tested)
2094 static const double q000[3] = {0.0, 0.0, 0.0};
2095 static const double q001[3] = {0.0, 0.0, 1.0};
2096 static const double q010[3] = {0.0, 1.0, 0.0};
2097 static const double q100[3] = {1.0, 0.0, 0.0};
2098 return
2099 7 PCK::orient_3d(p1, p2, p3, q000) == ZERO &&
2100
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 PCK::orient_3d(p1, p2, p3, q001) == ZERO &&
2101
3/4
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 4 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
 13 PCK::orient_3d(p1, p2, p3, q010) == ZERO &&
2102
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
 10 PCK::orient_3d(p1, p2, p3, q100) == ZERO
2103 ;
2104 }
2105
2106 249 void initialize() {
2107 249 expansion::initialize();
2108 249 }
2109
2110 249 void terminate() {
2111 // Nothing to do.
2112 249 }
2113
2114 ✗ void show_stats() {
2115 ✗ PredicateStats::show_all_stats();
2116 ✗ expansion::show_all_stats();
2117 ✗ }
2118 }
2119 }
2120