GCC Code Coverage Report

Directory: ./
File: numerics/predicates.cpp
Date: 2026-09-27 03:10:11
Exec Total Coverage
Lines: 582 851 68.4%
Functions: 56 70 80.0%
Branches: 185 334 55.4%
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 3681406 bool lexico_compare_3d(const double* x, const double* y) {
129
2/2
✓ Branch 0 taken 2848370 times.
✓ Branch 1 taken 833036 times.
 3681406 if(x[0] < y[0]) {
130 return true;
131 }
132
2/2
✓ Branch 0 taken 1504293 times.
✓ Branch 1 taken 1344077 times.
 2848370 if(x[0] > y[0]) {
133 return false;
134 }
135
2/2
✓ Branch 0 taken 1128249 times.
✓ Branch 1 taken 376044 times.
 1504293 if(x[1] < y[1]) {
136 return true;
137 }
138
2/2
✓ Branch 0 taken 492699 times.
✓ Branch 1 taken 635550 times.
 1128249 if(x[1] > y[1]) {
139 return false;
140 }
141 492699 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 993937 void GEOGRAM_API SOS_sort(
155 const double** begin, const double** end, index_t dim
156 ) {
157
2/2
✓ Branch 0 taken 437467 times.
✓ Branch 1 taken 556470 times.
 993937 if(SOS_mode_ == PCK::SOS_ADDRESS) {
158 std::sort(begin, end);
159 } else {
160
1/2
✓ Branch 0 taken 556470 times.
✗ Branch 1 not taken.
 556470 if(dim == 3) {
161 std::sort(begin, end, lexico_compare_3d);
162 } else {
163 std::sort(begin, end, LexicoCompare(dim));
164 }
165 }
166 993937 }
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 inline double max4(double x1, double x2, double x3, double x4) {
175 #ifdef __SSE2__
176 double result;
177 __m128d X1 =_mm_load_sd(&x1);
178 __m128d X2 =_mm_load_sd(&x2);
179 __m128d X3 =_mm_load_sd(&x3);
180 __m128d X4 =_mm_load_sd(&x4);
181 X1 = _mm_max_sd(X1,X2);
182 X3 = _mm_max_sd(X3,X4);
183 X1 = _mm_max_sd(X1,X3);
184 _mm_store_sd(&result, X1);
185 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 inline void get_minmax3(
199 double& m, double& M, double x1, double x2, double x3
200 ) {
201 #ifdef __SSE2__
202 __m128d X1 =_mm_load_sd(&x1);
203 __m128d X2 =_mm_load_sd(&x2);
204 __m128d X3 =_mm_load_sd(&x3);
205 __m128d MIN12 = _mm_min_sd(X1,X2);
206 __m128d MAX12 = _mm_max_sd(X1,X2);
207 X1 = _mm_min_sd(MIN12, X3);
208 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 }
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 865999 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 865999 double ptx = p[0] - t[0];
239 865999 double pty = p[1] - t[1];
240
1/2
✓ Branch 0 taken 865999 times.
✗ Branch 1 not taken.
 865999 double ptz = p[2] - t[2];
241 865999 double pt2 = geo_sqr(ptx) + geo_sqr(pty) + geo_sqr(ptz);
242
243 865999 double qtx = q[0] - t[0];
244 865999 double qty = q[1] - t[1];
245 865999 double qtz = q[2] - t[2];
246 865999 double qt2 = geo_sqr(qtx) + geo_sqr(qty) + geo_sqr(qtz);
247
248 865999 double rtx = r[0] - t[0];
249 865999 double rty = r[1] - t[1];
250 865999 double rtz = r[2] - t[2];
251 865999 double rt2 = geo_sqr(rtx) + geo_sqr(rty) + geo_sqr(rtz);
252
253 865999 double stx = s[0] - t[0];
254 865999 double sty = s[1] - t[1];
255 865999 double stz = s[2] - t[2];
256 865999 double st2 = geo_sqr(stx) + geo_sqr(sty) + geo_sqr(stz);
257
258 // Compute the semi-static bound.
259 865999 double maxx = ::fabs(ptx);
260 865999 double maxy = ::fabs(pty);
261 865999 double maxz = ::fabs(ptz);
262
263 865999 double aqtx = ::fabs(qtx);
264 865999 double artx = ::fabs(rtx);
265 865999 double astx = ::fabs(stx);
266
267 865999 double aqty = ::fabs(qty);
268 865999 double arty = ::fabs(rty);
269 865999 double asty = ::fabs(sty);
270
271 865999 double aqtz = ::fabs(qtz);
272 865999 double artz = ::fabs(rtz);
273
1/2
✓ Branch 0 taken 865999 times.
✗ Branch 1 not taken.
 865999 double astz = ::fabs(stz);
274
275 maxx = max4(maxx, aqtx, artx, astx);
276 maxy = max4(maxy, aqty, arty, asty);
277 maxz = max4(maxz, aqtz, artz, astz);
278
279
1/2
✓ Branch 0 taken 865999 times.
✗ Branch 1 not taken.
 865999 double eps = 1.2466136531027298e-13 * maxx * maxy * maxz;
280
281 double min_max;
282 double max_max;
283 get_minmax3(min_max, max_max, maxx, maxy, maxz);
284
285 865999 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 taken 865999 times.
✗ Branch 1 not taken.
 865999 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 865999 times.
✗ Branch 1 not taken.
 865999 } else if (max_max < 1e61) { /* sqrt^5(max_double/4 [hadamard]) */
296 // Protect against overflow in the computation of det.
297 865999 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 655350 times.
✓ Branch 1 taken 210649 times.
 865999 if (det > eps) return -1;
303
2/2
✓ Branch 0 taken 572913 times.
✓ Branch 1 taken 82437 times.
 655350 if (det < -eps) return 1;
304 }
305
306 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 stats_side1.log_exact();
334 207030 expansion& l = expansion_sq_dist(p0, p1, dim);
335 207030 expansion& a = expansion_dot_at(p1, q0, p0, dim).scale_fast(2.0);
336 207030 expansion& r = expansion_diff(l, a);
337 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 stats_side1.log_SOS();
341
2/2
✓ Branch 0 taken 50006 times.
✓ Branch 1 taken 154574 times.
 254586 return (p0 < p1) ? POSITIVE : NEGATIVE;
342 }
343 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 stats_side2.log_exact();
423
424 279520 const expansion& l1 = expansion_sq_dist(p1, p0, dim);
425 279520 const expansion& l2 = expansion_sq_dist(p2, p0, dim);
426
427 279520 const expansion& a10 = expansion_dot_at(p1,q0,p0, dim).scale_fast(2.0);
428 279520 const expansion& a11 = expansion_dot_at(p1,q1,p0, dim).scale_fast(2.0);
429 279520 const expansion& a20 = expansion_dot_at(p2,q0,p0, dim).scale_fast(2.0);
430 279520 const expansion& a21 = expansion_dot_at(p2,q1,p0, dim).scale_fast(2.0);
431
432 279520 const expansion& Delta = expansion_diff(a11, a10);
433
434 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 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 279520 geo_assert(Delta_sign != ZERO);
438
439 // [ Lambda0 ] [ -1 ] [ a11 ]
440 // Delta [ ] = [ ] * l1 + [ ]
441 // [ Lambda1 ] [ 1 ] [ -a10 ]
442
443 279520 const expansion& DeltaLambda0 = expansion_diff(a11, l1);
444 279520 const expansion& DeltaLambda1 = expansion_diff(l1, a10);
445
446 // r = Delta*l2 - ( a20*DeltaLambda0 + a21*DeltaLambda1 )
447
448 279520 const expansion& r0 = expansion_product(Delta, l2);
449 279520 const expansion& r1 = expansion_product(a20, DeltaLambda0).negate();
450 279520 const expansion& r2 = expansion_product(a21, DeltaLambda1).negate();
451 279520 const expansion& r = expansion_sum3(r0, r1, r2);
452
453 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 stats_side2.log_SOS();
458 271431 const double* p_sort[3] = {p0, p1, p2};
459 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 129865 const expansion& z1 = expansion_diff(Delta, a21);
463 129865 const expansion& z = expansion_sum(z1, a20);
464 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 134337 const expansion& z = expansion_diff(a21, a20);
471 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 53217 times.
✓ Branch 1 taken 60446 times.
 113663 if(p_sort[i] == p2) {
477 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 80655 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 stats_side3.log_exact();
568
569 80655 const expansion& l1 = expansion_sq_dist(p1, p0, dim);
570 80655 const expansion& l2 = expansion_sq_dist(p2, p0, dim);
571 80655 const expansion& l3 = expansion_sq_dist(p3, p0, dim);
572
573 80655 const expansion& a10 = expansion_dot_at(p1,q0,p0, dim).scale_fast(2.0);
574 80655 const expansion& a11 = expansion_dot_at(p1,q1,p0, dim).scale_fast(2.0);
575 80655 const expansion& a12 = expansion_dot_at(p1,q2,p0, dim).scale_fast(2.0);
576 80655 const expansion& a20 = expansion_dot_at(p2,q0,p0, dim).scale_fast(2.0);
577 80655 const expansion& a21 = expansion_dot_at(p2,q1,p0, dim).scale_fast(2.0);
578 80655 const expansion& a22 = expansion_dot_at(p2,q2,p0, dim).scale_fast(2.0);
579
580 80655 const expansion& a30 = expansion_dot_at(p3,q0,p0, dim).scale_fast(2.0);
581 80655 const expansion& a31 = expansion_dot_at(p3,q1,p0, dim).scale_fast(2.0);
582 80655 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 80655 const expansion& b00 = expansion_det2x2(a11, a12, a21, a22);
589 80655 const expansion& b01 = expansion_diff(a21, a22);
590 80655 const expansion& b02 = expansion_diff(a12, a11);
591 80655 const expansion& b10 = expansion_det2x2(a12, a10, a22, a20);
592 80655 const expansion& b11 = expansion_diff(a22, a20);
593 80655 const expansion& b12 = expansion_diff(a10, a12);
594 80655 const expansion& b20 = expansion_det2x2(a10, a11, a20, a21);
595 80655 const expansion& b21 = expansion_diff(a20, a21);
596 80655 const expansion& b22 = expansion_diff(a11, a10);
597
598 80655 const expansion& Delta = expansion_sum3(b00, b10, b20);
599 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 80655 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 80655 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 80655 const expansion& b01_l1 = expansion_product(b01, l1);
609 80655 const expansion& b02_l2 = expansion_product(b02, l2);
610 80655 const expansion& DeltaLambda0 = expansion_sum3(b01_l1, b02_l2, b00);
611
612 80655 const expansion& b11_l1 = expansion_product(b11, l1);
613 80655 const expansion& b12_l2 = expansion_product(b12, l2);
614 80655 const expansion& DeltaLambda1 = expansion_sum3(b11_l1, b12_l2, b10);
615
616 80655 const expansion& b21_l1 = expansion_product(b21, l1);
617 80655 const expansion& b22_l2 = expansion_product(b22, l2);
618 80655 const expansion& DeltaLambda2 = expansion_sum3(b21_l1, b22_l2, b20);
619
620 // r = Delta*l3-(a30*DeltaLambda0+a31*DeltaLambda1+a32*DeltaLambda2)
621
622 80655 const expansion& r0 = expansion_product(Delta, l3);
623 80655 const expansion& r1 = expansion_product(a30, DeltaLambda0).negate();
624 80655 const expansion& r2 = expansion_product(a31, DeltaLambda1).negate();
625 80655 const expansion& r3 = expansion_product(a32, DeltaLambda2).negate();
626 80655 const expansion& r = expansion_sum4(r0, r1, r2, r3);
627 Sign r_sign = r.sign();
628
629 // Simulation of Simplicity (symbolic perturbation)
630
2/2
✓ Branch 0 taken 72918 times.
✓ Branch 1 taken 7737 times.
 80655 if(r_sign == ZERO) {
631 stats_side3.log_SOS();
632 72918 const double* p_sort[4] = {p0, p1, p2, p3};
633 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 35915 const expansion& z1_0 = expansion_sum(b01, b02);
637 35915 const expansion& z1 = expansion_product(a30, z1_0).negate();
638 35915 const expansion& z2_0 = expansion_sum(b11, b12);
639 35915 const expansion& z2 = expansion_product(a31, z2_0).negate();
640 35915 const expansion& z3_0 = expansion_sum(b21, b22);
641 35915 const expansion& z3 = expansion_product(a32, z3_0).negate();
642 35915 const expansion& z = expansion_sum4(Delta, z1, z2, z3);
643 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 10534 const expansion& z1 = expansion_product(a30, b01);
649 10534 const expansion& z2 = expansion_product(a31, b11);
650 10534 const expansion& z3 = expansion_product(a32, b21);
651 10534 const expansion& z = expansion_sum3(z1, z2, z3);
652 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 45832 const expansion& z1 = expansion_product(a30, b02);
658 45832 const expansion& z2 = expansion_product(a31, b12);
659 45832 const expansion& z3 = expansion_product(a32, b22);
660 45832 const expansion& z = expansion_sum3(z1, z2, z3);
661 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 not taken.
✓ Branch 1 taken 11769 times.
 11769 } else if(p_sort[i] == p3) {
666 return NEGATIVE;
667 }
668 }
669 ✗ geo_assert_not_reached;
670 }
671 7737 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 137 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 137 const expansion& l1 = expansion_diff(h1,h0);
687 137 const expansion& l2 = expansion_diff(h2,h0);
688 137 const expansion& l3 = expansion_diff(h3,h0);
689
690 137 const expansion& a10 = expansion_dot_at(p1, q0, p0, 3).scale_fast(2.0);
691 137 const expansion& a11 = expansion_dot_at(p1, q1, p0, 3).scale_fast(2.0);
692 137 const expansion& a12 = expansion_dot_at(p1, q2, p0, 3).scale_fast(2.0);
693 137 const expansion& a20 = expansion_dot_at(p2, q0, p0, 3).scale_fast(2.0);
694 137 const expansion& a21 = expansion_dot_at(p2, q1, p0, 3).scale_fast(2.0);
695 137 const expansion& a22 = expansion_dot_at(p2, q2, p0, 3).scale_fast(2.0);
696
697 137 const expansion& a30 = expansion_dot_at(p3, q0, p0, 3).scale_fast(2.0);
698 137 const expansion& a31 = expansion_dot_at(p3, q1, p0, 3).scale_fast(2.0);
699 137 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 137 const expansion& b00 = expansion_det2x2(a11, a12, a21, a22);
706 137 const expansion& b01 = expansion_diff(a21, a22);
707 137 const expansion& b02 = expansion_diff(a12, a11);
708 137 const expansion& b10 = expansion_det2x2(a12, a10, a22, a20);
709 137 const expansion& b11 = expansion_diff(a22, a20);
710 137 const expansion& b12 = expansion_diff(a10, a12);
711 137 const expansion& b20 = expansion_det2x2(a10, a11, a20, a21);
712 137 const expansion& b21 = expansion_diff(a20, a21);
713 137 const expansion& b22 = expansion_diff(a11, a10);
714
715 137 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
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 137 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 137 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 137 const expansion& b01_l1 = expansion_product(b01, l1);
726 137 const expansion& b02_l2 = expansion_product(b02, l2);
727 137 const expansion& DeltaLambda0 = expansion_sum3(b01_l1, b02_l2, b00);
728
729 137 const expansion& b11_l1 = expansion_product(b11, l1);
730 137 const expansion& b12_l2 = expansion_product(b12, l2);
731 137 const expansion& DeltaLambda1 = expansion_sum3(b11_l1, b12_l2, b10);
732
733 137 const expansion& b21_l1 = expansion_product(b21, l1);
734 137 const expansion& b22_l2 = expansion_product(b22, l2);
735 137 const expansion& DeltaLambda2 = expansion_sum3(b21_l1, b22_l2, b20);
736
737 // r = Delta*l3-(a30*DeltaLambda0+a31*DeltaLambda1+a32*DeltaLambda2)
738
739 137 const expansion& r0 = expansion_product(Delta, l3);
740 137 const expansion& r1 = expansion_product(a30, DeltaLambda0).negate();
741 137 const expansion& r2 = expansion_product(a31, DeltaLambda1).negate();
742 137 const expansion& r3 = expansion_product(a32, DeltaLambda2).negate();
743 137 const expansion& r = expansion_sum4(r0, r1, r2, r3);
744 Sign r_sign = r.sign();
745
746 // Simulation of Simplicity (symbolic perturbation)
747
2/2
✓ Branch 0 taken 52 times.
✓ Branch 1 taken 85 times.
 137 if(r_sign == ZERO) {
748 stats_side3h.log_SOS();
749 52 const double* p_sort[4] = {p0, p1, p2, p3};
750 52 SOS_sort(p_sort, p_sort + 4, 3);
751
1/2
✓ Branch 0 taken 52 times.
✗ Branch 1 not taken.
 52 for(index_t i = 0; i < 4; ++i) {
752
2/2
✓ Branch 0 taken 17 times.
✓ Branch 1 taken 35 times.
 52 if(p_sort[i] == p0) {
753 17 const expansion& z1_0 = expansion_sum(b01, b02);
754 17 const expansion& z1 = expansion_product(a30, z1_0).negate();
755 17 const expansion& z2_0 = expansion_sum(b11, b12);
756 17 const expansion& z2 = expansion_product(a31, z2_0).negate();
757 17 const expansion& z3_0 = expansion_sum(b21, b22);
758 17 const expansion& z3 = expansion_product(a32, z3_0).negate();
759 17 const expansion& z = expansion_sum4(Delta, z1, z2, z3);
760 Sign z_sign = z.sign();
761
1/2
✓ Branch 0 taken 17 times.
✗ Branch 1 not taken.
 17 if(z_sign != ZERO) {
762 17 return Sign(Delta_sign * z_sign);
763 }
764
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 29 times.
 35 } else if(p_sort[i] == p1) {
765 6 const expansion& z1 = expansion_product(a30, b01);
766 6 const expansion& z2 = expansion_product(a31, b11);
767 6 const expansion& z3 = expansion_product(a32, b21);
768 6 const expansion& z = expansion_sum3(z1, z2, z3);
769 Sign z_sign = z.sign();
770
1/2
✓ Branch 0 taken 6 times.
✗ Branch 1 not taken.
 6 if(z_sign != ZERO) {
771 6 return Sign(Delta_sign * z_sign);
772 }
773
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 15 times.
 29 } else if(p_sort[i] == p2) {
774 14 const expansion& z1 = expansion_product(a30, b02);
775 14 const expansion& z2 = expansion_product(a31, b12);
776 14 const expansion& z3 = expansion_product(a32, b22);
777 14 const expansion& z = expansion_sum3(z1, z2, z3);
778 Sign z_sign = z.sign();
779
1/2
✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
 14 if(z_sign != ZERO) {
780 14 return Sign(Delta_sign * z_sign);
781 }
782
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
 15 } else if(p_sort[i] == p3) {
783 return NEGATIVE;
784 }
785 }
786 ✗ geo_assert_not_reached;
787 }
788 85 return Sign(Delta_sign * r_sign);
789 }
790
791
792 /**
793 * \brief Implements side3() in 3d.
794 */
795 11053634 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 11053634 Sign result = Sign(side3_3d_filter(p0, p1, p2, p3, q0, q1, q2));
800
2/2
✓ Branch 0 taken 80627 times.
✓ Branch 1 taken 10973007 times.
 11053634 if(result == ZERO) {
801 80627 result = side3_exact_SOS(p0, p1, p2, p3, q0, q1, q2, 3);
802 }
803 11053634 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 180792 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 stats_side4.log_exact();
876
877 180792 const expansion& a11 = expansion_diff(p1[0], p0[0]);
878 180792 const expansion& a12 = expansion_diff(p1[1], p0[1]);
879 180792 const expansion& a13 = expansion_diff(p1[2], p0[2]);
880 180792 const expansion& a14 = expansion_sq_dist(p1, p0, 3).negate();
881
882 180792 const expansion& a21 = expansion_diff(p2[0], p0[0]);
883 180792 const expansion& a22 = expansion_diff(p2[1], p0[1]);
884 180792 const expansion& a23 = expansion_diff(p2[2], p0[2]);
885 180792 const expansion& a24 = expansion_sq_dist(p2, p0, 3).negate();
886
887 180792 const expansion& a31 = expansion_diff(p3[0], p0[0]);
888 180792 const expansion& a32 = expansion_diff(p3[1], p0[1]);
889 180792 const expansion& a33 = expansion_diff(p3[2], p0[2]);
890 180792 const expansion& a34 = expansion_sq_dist(p3, p0, 3).negate();
891
892 180792 const expansion& a41 = expansion_diff(p4[0], p0[0]);
893 180792 const expansion& a42 = expansion_diff(p4[1], p0[1]);
894 180792 const expansion& a43 = expansion_diff(p4[2], p0[2]);
895 180792 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 180792 const expansion& m12 = expansion_det2x2(a12,a13,a22,a23);
925 180792 const expansion& m13 = expansion_det2x2(a12,a13,a32,a33);
926 180792 const expansion& m14 = expansion_det2x2(a12,a13,a42,a43);
927 180792 const expansion& m23 = expansion_det2x2(a22,a23,a32,a33);
928 180792 const expansion& m24 = expansion_det2x2(a22,a23,a42,a43);
929 180792 const expansion& m34 = expansion_det2x2(a32,a33,a42,a43);
930
931
932 180792 const expansion& z11 = expansion_product(a21,m34);
933 180792 const expansion& z12 = expansion_product(a31,m24).negate();
934 180792 const expansion& z13 = expansion_product(a41,m23);
935 180792 const expansion& Delta1 = expansion_sum3(z11,z12,z13);
936
937 180792 const expansion& z21 = expansion_product(a11,m34);
938 180792 const expansion& z22 = expansion_product(a31,m14).negate();
939 180792 const expansion& z23 = expansion_product(a41,m13);
940 180792 const expansion& Delta2 = expansion_sum3(z21,z22,z23);
941
942 180792 const expansion& z31 = expansion_product(a11,m24);
943 180792 const expansion& z32 = expansion_product(a21,m14).negate();
944 180792 const expansion& z33 = expansion_product(a41,m12);
945 180792 const expansion& Delta3 = expansion_sum3(z31,z32,z33);
946
947 180792 const expansion& z41 = expansion_product(a11,m23);
948 180792 const expansion& z42 = expansion_product(a21,m13).negate();
949 180792 const expansion& z43 = expansion_product(a31,m12);
950 180792 const expansion& Delta4 = expansion_sum3(z41,z42,z43);
951
952
953 Sign Delta4_sign = Delta4.sign();
954
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 180792 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 180792 geo_assert(Delta4_sign != ZERO);
955
956 180792 const expansion& r_1 = expansion_product(Delta1, a14);
957 180792 const expansion& r_2 = expansion_product(Delta2, a24).negate();
958 180792 const expansion& r_3 = expansion_product(Delta3, a34);
959 180792 const expansion& r_4 = expansion_product(Delta4, a44).negate();
960 180792 const expansion& r = expansion_sum4(r_1, r_2, r_3, r_4);
961 Sign r_sign = r.sign();
962
963 // Simulation of Simplicity (symbolic perturbation)
964
2/2
✓ Branch 0 taken 92941 times.
✓ Branch 1 taken 87851 times.
 180792 if(sos && r_sign == ZERO) {
965 stats_side4.log_SOS();
966 92941 const double* p_sort[5] = {p0, p1, p2, p3, p4};
967 92941 SOS_sort(p_sort, p_sort + 5, 3);
968
1/2
✓ Branch 0 taken 272195 times.
✗ Branch 1 not taken.
 272195 for(index_t i = 0; i < 5; ++i) {
969
2/2
✓ Branch 0 taken 90257 times.
✓ Branch 1 taken 181938 times.
 272195 if(p_sort[i] == p0) {
970 90257 const expansion& z1 = expansion_diff(Delta2, Delta1);
971 90257 const expansion& z2 = expansion_diff(Delta4, Delta3);
972 90257 const expansion& z = expansion_sum(z1, z2);
973 Sign z_sign = z.sign();
974
2/2
✓ Branch 0 taken 921 times.
✓ Branch 1 taken 89336 times.
 90257 if(z_sign != ZERO) {
975 92941 return Sign(Delta4_sign * z_sign);
976 }
977
2/2
✓ Branch 0 taken 30601 times.
✓ Branch 1 taken 151337 times.
 181938 } else if(p_sort[i] == p1) {
978 Sign Delta1_sign = Delta1.sign();
979
2/2
✓ Branch 0 taken 30299 times.
✓ Branch 1 taken 302 times.
 30601 if(Delta1_sign != ZERO) {
980 30299 return Sign(Delta4_sign * Delta1_sign);
981 }
982
2/2
✓ Branch 0 taken 60318 times.
✓ Branch 1 taken 91019 times.
 151337 } else if(p_sort[i] == p2) {
983 Sign Delta2_sign = Delta2.sign();
984
2/2
✓ Branch 0 taken 30338 times.
✓ Branch 1 taken 29980 times.
 60318 if(Delta2_sign != ZERO) {
985 30338 return Sign(-Delta4_sign * Delta2_sign);
986 }
987
2/2
✓ Branch 0 taken 89998 times.
✓ Branch 1 taken 1021 times.
 91019 } else if(p_sort[i] == p3) {
988 Sign Delta3_sign = Delta3.sign();
989
2/2
✓ Branch 0 taken 30362 times.
✓ Branch 1 taken 59636 times.
 89998 if(Delta3_sign != ZERO) {
990 30362 return Sign(Delta4_sign * Delta3_sign);
991 }
992
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1021 times.
 1021 } else if(p_sort[i] == p4) {
993 return NEGATIVE;
994 }
995 }
996 }
997 87851 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 1178251 Sign orient_2d_exact(const double* p0, const double* p1, const double* p2) {
1234 stats_orient2d.log_exact();
1235 1178251 const expansion& a11 = expansion_diff(p1[0], p0[0]);
1236 1178251 const expansion& a12 = expansion_diff(p1[1], p0[1]);
1237 1178251 const expansion& a21 = expansion_diff(p2[0], p0[0]);
1238 1178251 const expansion& a22 = expansion_diff(p2[1], p0[1]);
1239 1178251 const expansion& Delta = expansion_det2x2(a11, a12, a21, a22);
1240 1178251 return Delta.sign();
1241 }
1242
1243 // ============ orient3d ===================================================
1244
1245 2161322 Sign orient_3d_exact(
1246 const double* p0, const double* p1, const double* p2, const double* p3
1247 ) {
1248 stats_orient3d.log_exact();
1249
1250 2161322 const expansion& a11 = expansion_diff(p1[0], p0[0]);
1251 2161322 const expansion& a12 = expansion_diff(p1[1], p0[1]);
1252 2161322 const expansion& a13 = expansion_diff(p1[2], p0[2]);
1253
1254 2161322 const expansion& a21 = expansion_diff(p2[0], p0[0]);
1255 2161322 const expansion& a22 = expansion_diff(p2[1], p0[1]);
1256 2161322 const expansion& a23 = expansion_diff(p2[2], p0[2]);
1257
1258 2161322 const expansion& a31 = expansion_diff(p3[0], p0[0]);
1259 2161322 const expansion& a32 = expansion_diff(p3[1], p0[1]);
1260 2161322 const expansion& a33 = expansion_diff(p3[2], p0[2]);
1261
1262 2161322 const expansion& Delta = expansion_det3x3(
1263 a11, a12, a13, a21, a22, a23, a31, a32, a33
1264 );
1265
1266 2161322 return Delta.sign();
1267 }
1268
1269 3237 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 3237 const expansion& a11 = expansion_diff(p1[0], p0[0]);
1278 3237 const expansion& a12 = expansion_diff(p1[1], p0[1]);
1279 3237 const expansion& a13 = expansion_diff(p1[2], p0[2]);
1280 3237 const expansion& a14 = expansion_diff(h0,h1);
1281
1282 3237 const expansion& a21 = expansion_diff(p2[0], p0[0]);
1283 3237 const expansion& a22 = expansion_diff(p2[1], p0[1]);
1284 3237 const expansion& a23 = expansion_diff(p2[2], p0[2]);
1285 3237 const expansion& a24 = expansion_diff(h0,h2);
1286
1287 3237 const expansion& a31 = expansion_diff(p3[0], p0[0]);
1288 3237 const expansion& a32 = expansion_diff(p3[1], p0[1]);
1289 3237 const expansion& a33 = expansion_diff(p3[2], p0[2]);
1290 3237 const expansion& a34 = expansion_diff(h0,h3);
1291
1292 3237 const expansion& a41 = expansion_diff(p4[0], p0[0]);
1293 3237 const expansion& a42 = expansion_diff(p4[1], p0[1]);
1294 3237 const expansion& a43 = expansion_diff(p4[2], p0[2]);
1295 3237 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 3237 const expansion& Delta1 = expansion_det3x3(
1301 a21, a22, a23,
1302 a31, a32, a33,
1303 a41, a42, a43
1304 );
1305 3237 const expansion& Delta2 = expansion_det3x3(
1306 a11, a12, a13,
1307 a31, a32, a33,
1308 a41, a42, a43
1309 );
1310 3237 const expansion& Delta3 = expansion_det3x3(
1311 a11, a12, a13,
1312 a21, a22, a23,
1313 a41, a42, a43
1314 );
1315 3237 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
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 3237 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 3237 geo_assert(Delta4_sign != ZERO);
1323
1324 3237 const expansion& r_1 = expansion_product(Delta1, a14);
1325 3237 const expansion& r_2 = expansion_product(Delta2, a24).negate();
1326 3237 const expansion& r_3 = expansion_product(Delta3, a34);
1327 3237 const expansion& r_4 = expansion_product(Delta4, a44).negate();
1328 3237 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
2/2
✓ Branch 0 taken 1234 times.
✓ Branch 1 taken 2003 times.
 3237 if(sos && r_sign == ZERO) {
1334 stats_orient3dh.log_SOS();
1335 1234 const double* p_sort[5] = {p0, p1, p2, p3, p4};
1336 1234 SOS_sort(p_sort, p_sort + 5, 3);
1337
1/2
✓ Branch 0 taken 1660 times.
✗ Branch 1 not taken.
 1660 for(index_t i = 0; i < 5; ++i) {
1338
2/2
✓ Branch 0 taken 402 times.
✓ Branch 1 taken 1258 times.
 1660 if(p_sort[i] == p0) {
1339 402 const expansion& z1 = expansion_diff(Delta2, Delta1);
1340 402 const expansion& z2 = expansion_diff(Delta4, Delta3);
1341 402 const expansion& z = expansion_sum(z1, z2);
1342 Sign z_sign = z.sign();
1343
2/2
✓ Branch 0 taken 303 times.
✓ Branch 1 taken 99 times.
 402 if(z_sign != ZERO) {
1344 1234 return Sign(Delta4_sign * z_sign);
1345 }
1346
2/2
✓ Branch 0 taken 255 times.
✓ Branch 1 taken 1003 times.
 1258 } else if(p_sort[i] == p1) {
1347 Sign Delta1_sign = Delta1.sign();
1348
2/2
✓ Branch 0 taken 137 times.
✓ Branch 1 taken 118 times.
 255 if(Delta1_sign != ZERO) {
1349 137 return Sign(Delta4_sign * Delta1_sign);
1350 }
1351
2/2
✓ Branch 0 taken 284 times.
✓ Branch 1 taken 719 times.
 1003 } else if(p_sort[i] == p2) {
1352 Sign Delta2_sign = Delta2.sign();
1353
2/2
✓ Branch 0 taken 168 times.
✓ Branch 1 taken 116 times.
 284 if(Delta2_sign != ZERO) {
1354 168 return Sign(-Delta4_sign * Delta2_sign);
1355 }
1356
2/2
✓ Branch 0 taken 284 times.
✓ Branch 1 taken 435 times.
 719 } else if(p_sort[i] == p3) {
1357 Sign Delta3_sign = Delta3.sign();
1358
2/2
✓ Branch 0 taken 191 times.
✓ Branch 1 taken 93 times.
 284 if(Delta3_sign != ZERO) {
1359 191 return Sign(Delta4_sign * Delta3_sign);
1360 }
1361
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 435 times.
 435 } else if(p_sort[i] == p4) {
1362 return NEGATIVE;
1363 }
1364 }
1365 }
1366 2003 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 1348 Sign det_3d_exact(
1452 const double* p0, const double* p1, const double* p2
1453 ) {
1454 stats_det3d.log_exact();
1455
1456 1348 const expansion& p0_0 = expansion_create(p0[0]);
1457 1348 const expansion& p0_1 = expansion_create(p0[1]);
1458 1348 const expansion& p0_2 = expansion_create(p0[2]);
1459
1460 1348 const expansion& p1_0 = expansion_create(p1[0]);
1461 1348 const expansion& p1_1 = expansion_create(p1[1]);
1462 1348 const expansion& p1_2 = expansion_create(p1[2]);
1463
1464 1348 const expansion& p2_0 = expansion_create(p2[0]);
1465 1348 const expansion& p2_1 = expansion_create(p2[1]);
1466 1348 const expansion& p2_2 = expansion_create(p2[2]);
1467
1468 1348 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 1348 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 2077270 bool aligned_3d_exact(
1486 const double* p0, const double* p1, const double* p2
1487 ) {
1488 2077270 const expansion& U_0 = expansion_diff(p1[0],p0[0]);
1489 2077270 const expansion& U_1 = expansion_diff(p1[1],p0[1]);
1490 2077270 const expansion& U_2 = expansion_diff(p1[2],p0[2]);
1491
1492 2077270 const expansion& V_0 = expansion_diff(p2[0],p0[0]);
1493 2077270 const expansion& V_1 = expansion_diff(p2[1],p0[1]);
1494 2077270 const expansion& V_2 = expansion_diff(p2[2],p0[2]);
1495
1496 2077270 const expansion& N_0 = expansion_det2x2(U_1, V_1, U_2, V_2);
1497 2077270 const expansion& N_1 = expansion_det2x2(U_2, V_2, U_0, V_0);
1498 2077270 const expansion& N_2 = expansion_det2x2(U_0, V_0, U_1, V_1);
1499
1500 return(
1501
2/2
✓ Branch 0 taken 314184 times.
✓ Branch 1 taken 67279 times.
 381463 N_0.sign() == 0 &&
1502
3/4
✓ Branch 0 taken 381463 times.
✓ Branch 1 taken 1695807 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 314184 times.
 2391454 N_1.sign() == 0 &&
1503 N_2.sign() == 0
1504 2077270 );
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 145152 Sign dot_3d_exact(
1515 const double* p0, const double* p1, const double* p2
1516 ) {
1517 145152 const expansion& U_0 = expansion_diff(p1[0],p0[0]);
1518 145152 const expansion& U_1 = expansion_diff(p1[1],p0[1]);
1519 145152 const expansion& U_2 = expansion_diff(p1[2],p0[2]);
1520
1521 145152 const expansion& V_0 = expansion_diff(p2[0],p0[0]);
1522 145152 const expansion& V_1 = expansion_diff(p2[1],p0[1]);
1523 145152 const expansion& V_2 = expansion_diff(p2[2],p0[2]);
1524
1525 145152 const expansion& UV_0 = expansion_product(U_0, V_0);
1526 145152 const expansion& UV_1 = expansion_product(U_1, V_1);
1527 145152 const expansion& UV_2 = expansion_product(U_2, V_2);
1528
1529 145152 const expansion& Delta = expansion_sum3(UV_0, UV_1, UV_2);
1530
1531 145152 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 123 void set_SOS_mode(SOSMode m) {
1567 123 SOS_mode_ = m;
1568 123 }
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 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 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 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 3951 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 3951 side3h_3d_filter(p0, p1, p2, p3, h0, h1, h2, h3, q0, q1, q2)
1649 );
1650
2/2
✓ Branch 0 taken 137 times.
✓ Branch 1 taken 3814 times.
 3951 if(SOS && result == ZERO) {
1651 137 result = side3h_exact_SOS(
1652 p0, p1, p2, p3, h0, h1, h2, h3, q0, q1, q2
1653 );
1654 }
1655 3951 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 case 4:
1676 stats_side4.log_invoke();
1677 25769 return side4_4d_SOS(p0, p1, p2, p3, p4, q0, q1, q2, q3);
1678 case 6:
1679 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 case 8:
1685 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 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 865999 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 stats_side4.log_invoke();
1738
1739 // This specialized filter supposes that orient_3d(p0,p1,p2,p3) > 0
1740
1741 865999 Sign result = Sign(in_sphere_3d_filter_optim(p0, p1, p2, p3, p4));
1742
1743
2/2
✓ Branch 0 taken 82437 times.
✓ Branch 1 taken 783562 times.
 865999 if(result == 0) {
1744 82437 result = side4_3d_exact_SOS(p0, p1, p2, p3, p4);
1745 }
1746 865999 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 28 times.
✓ Branch 1 taken 936157 times.
 936185 if(s != ZERO) {
1770 return s;
1771 }
1772 28 return Sign(-side3_exact_SOS(p0, p1, p2, p3, p0, p1, p2, 2));
1773 }
1774
1775 40100 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 40100 return Sign(-side3_3d_SOS(p0,p1,p2,p3,p0,p1,p2));
1792 }
1793
1794 3951 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 3951 -side3_3dlifted_SOS(p0,p1,p2,p3,h0,h1,h2,h3,p0,p1,p2,SOS)
1806 3951 );
1807 }
1808
1809
1810 17501481 Sign orient_2d(
1811 const double* p0, const double* p1, const double* p2
1812 ) {
1813 stats_orient2d.log_invoke();
1814 17501481 Sign result = Sign(orient_2d_filter(p0, p1, p2));
1815
2/2
✓ Branch 0 taken 1178251 times.
✓ Branch 1 taken 16323230 times.
 17501481 if(result == 0) {
1816 1178251 result = orient_2d_exact(p0, p1, p2);
1817 }
1818 17501481 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 11263781 Sign orient_3d(
1842 const double* p0, const double* p1,
1843 const double* p2, const double* p3
1844 ) {
1845 stats_orient3d.log_invoke();
1846 11263781 Sign result = Sign(orient_3d_filter(p0, p1, p2, p3));
1847
2/2
✓ Branch 0 taken 2161322 times.
✓ Branch 1 taken 9102459 times.
 11263781 if(result == 0) {
1848 2161322 result = orient_3d_exact(p0, p1, p2, p3);
1849 }
1850 11263781 return result;
1851 }
1852
1853 6044552 Sign orient_3d_SOS(
1854 const double* p0, const double* p1,
1855 const double* p2, const double* p3
1856 ) {
1857 struct SOS {
1858 555361 SOS(
1859 const double* p0, const double* p1,
1860 const double* p2, const double* p3
1861 555361 ) : p_orig{p0,p1,p2,p3}, p_sort{p0, p1, p2, p3} {
1862 555361 SOS_sort(p_sort, p_sort+4, 3);
1863 1110722 parity = Permutation::permutation_is_odd(p_orig, p_sort, 4)
1864
2/2
✓ Branch 0 taken 277443 times.
✓ Branch 1 taken 277918 times.
 555361 ? NEGATIVE : POSITIVE;
1865 555361 }
1866 Sign orient_1d(index_t i, index_t j, index_t ax) const {
1867
2/2
✓ Branch 0 taken 31672 times.
✓ Branch 1 taken 45872 times.
 77544 return Sign(parity * geo_cmp(p_sort[i][ax], p_sort[j][ax]));
1868 }
1869 1177574 Sign orient_2d(
1870 index_t i, index_t j, index_t k, index_t ax1, index_t ax2
1871 ) const {
1872 1177574 double x0 = p_sort[i][ax1]; double y0 = p_sort[i][ax2];
1873 1177574 double x1 = p_sort[j][ax1]; double y1 = p_sort[j][ax2];
1874 1177574 double x2 = p_sort[k][ax1]; double y2 = p_sort[k][ax2];
1875 1177574 const expansion& a11 = expansion_diff(x1, x0);
1876 1177574 const expansion& a12 = expansion_diff(y1, y0);
1877 1177574 const expansion& a21 = expansion_diff(x2, x0);
1878 1177574 const expansion& a22 = expansion_diff(y2, y0);
1879 1177574 const expansion& D = expansion_det2x2(a11, a12, a21, a22);
1880
1/2
✓ Branch 0 taken 1177574 times.
✗ Branch 1 not taken.
 1177574 return Sign(parity * D.sign());
1881 }
1882 const double* p_orig[4];
1883 const double* p_sort[4];
1884 Sign parity;
1885 };
1886
1887 6044552 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 113342 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 113342 side4h_3d_filter(
1920 p0, p1, p2, p3, p4, h0, h1, h2, h3, h4
1921 )
1922 );
1923
2/2
✓ Branch 0 taken 3237 times.
✓ Branch 1 taken 110105 times.
 113342 if(result == 0) {
1924 3237 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 113342 return Sign(-result);
1931 }
1932
1933 2124 Sign det_3d(
1934 const double* p0, const double* p1, const double* p2
1935 ) {
1936 stats_det3d.log_invoke();
1937 Sign result = Sign(
1938 2124 det_3d_filter(p0, p1, p2)
1939 );
1940
2/2
✓ Branch 0 taken 1348 times.
✓ Branch 1 taken 776 times.
 2124 if(result == 0) {
1941 1348 result = det_3d_exact(p0, p1, p2);
1942 }
1943 2124 return result;
1944 }
1945
1946
1947 4386 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 4386 det_4d_filter(p0, p1, p2, p3)
1954 );
1955
1956
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4386 times.
 4386 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 4386 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 2077270 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 2077270 return aligned_3d_exact(p0, p1, p2);
2042 }
2043
2044 145152 Sign dot_3d(
2045 const double* p0, const double* p1, const double* p2
2046 ) {
2047 145152 Sign result = Sign(det_3d_filter(p0, p1, p2));
2048
1/2
✓ Branch 0 taken 145152 times.
✗ Branch 1 not taken.
 145152 if(result == 0) {
2049 145152 result = dot_3d_exact(p0, p1, p2);
2050 }
2051 145152 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.
 3 (p1[0] == p2[0]) &&
2071
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
 3 (p1[1] == p2[1])
2072 ;
2073 }
2074
2075 10 bool points_are_identical_3d(
2076 const double* p1,
2077 const double* p2
2078 ) {
2079 return
2080 18 (p1[0] == p2[0]) &&
2081
4/4
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 5 times.
 10 (p1[1] == p2[1]) &&
2082
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 (p1[2] == p2[2])
2083 ;
2084 }
2085
2086 10 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
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
 13 PCK::orient_3d(p1, p2, p3, q000) == ZERO &&
2100
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
 6 PCK::orient_3d(p1, p2, p3, q001) == ZERO &&
2101
3/4
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 7 times.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
 16 PCK::orient_3d(p1, p2, p3, q010) == ZERO &&
2102 3 PCK::orient_3d(p1, p2, p3, q100) == ZERO
2103 ;
2104 }
2105
2106 253 void initialize() {
2107 253 expansion::initialize();
2108 253 }
2109
2110 253 void terminate() {
2111 // Nothing to do.
2112 253 }
2113
2114 ✗ void show_stats() {
2115 ✗ PredicateStats::show_all_stats();
2116 ✗ expansion::show_all_stats();
2117 ✗ }
2118 }
2119 }
2120