GCC Code Coverage Report

Directory: ./
File: lib/geogram/numerics/exact_geometry.cpp
Date: 2026-09-07 02:36:43
Exec Total Coverage
Lines: 272 338 80.5%
Functions: 17 26 65.4%
Branches: 368 898 41.0%
Line Branch Exec Source
1 /*
2 * Copyright (c) 2000-2023 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/numerics/exact_geometry.h>
41 #include <geogram/numerics/interval_nt.h>
42 #include <geogram/numerics/predicates.h>
43 #include <geogram/numerics/PCK.h>
44 #include <geogram/basic/logger.h>
45
46 namespace GEO {
47
48
49 namespace PCK {
50
51 670938 Sign orient_2d(
52 const vec2HE& p0, const vec2HE& p1, const vec2HE& p2
53 ) {
54
5/8
✓ Branch 0 taken 37 times.
✓ Branch 1 taken 670901 times.
✓ Branch 3 taken 34 times.
✓ Branch 4 taken 3 times.
✓ Branch 6 taken 34 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 670938 static PredicateStats stats("orient_2d(vec2HE)");
55 670938 stats.log_invoke();
56 // Filter, using interval arithmetics
57 {
58 670938 interval_nt::Rounding rounding;
59
1/2
✓ Branch 1 taken 670938 times.
✗ Branch 2 not taken.
 670938 interval_nt Delta = det3x3(
60
3/6
✓ Branch 1 taken 670938 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 670938 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 670938 times.
✗ Branch 8 not taken.
 670938 interval_nt(p0.x),interval_nt(p0.y),interval_nt(p0.w),
61
3/6
✓ Branch 1 taken 670938 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 670938 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 670938 times.
✗ Branch 8 not taken.
 670938 interval_nt(p1.x),interval_nt(p1.y),interval_nt(p1.w),
62
3/6
✓ Branch 1 taken 670938 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 670938 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 670938 times.
✗ Branch 8 not taken.
 670938 interval_nt(p2.x),interval_nt(p2.y),interval_nt(p2.w)
63 );
64 670938 interval_nt::Sign2 s = Delta.sign();
65
2/2
✓ Branch 1 taken 660902 times.
✓ Branch 2 taken 10036 times.
 670938 if(interval_nt::sign_is_determined(s)) {
66 return Sign(
67
1/2
✓ Branch 1 taken 660902 times.
✗ Branch 2 not taken.
 660902 interval_nt::convert_sign(s)*
68
3/6
✓ Branch 1 taken 660902 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 660902 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 660902 times.
✗ Branch 8 not taken.
 660902 p0.w.sign()*p1.w.sign()*p2.w.sign()
69 660902 );
70 }
71
2/2
✓ Branch 1 taken 10036 times.
✓ Branch 2 taken 660902 times.
 670938 }
72 10036 stats.log_exact();
73 #ifdef GEO_HAS_BIG_STACK
74
2/6
✓ Branch 33 taken 10036 times.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✓ Branch 36 taken 10036 times.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
 10036 const expansion& Delta = expansion_det3x3(
75 p0.x.rep(), p0.y.rep(), p0.w.rep(),
76 p1.x.rep(), p1.y.rep(), p1.w.rep(),
77 p2.x.rep(), p2.y.rep(), p2.w.rep()
78 );
79 #else
80 expansion_nt Delta = det3x3(
81 p0.x, p0.y, p0.w,
82 p1.x, p1.y, p1.w,
83 p2.x, p2.y, p2.w
84 );
85 #endif
86 return Sign(
87 10036 Delta.sign()*
88 10036 p0.w.rep().sign()*
89 10036 p1.w.rep().sign()*
90 10036 p2.w.rep().sign()
91 10036 );
92 }
93
94 23715 Sign orient_3d(
95 const vec3HE& p0, const vec3HE& p1,
96 const vec3HE& p2, const vec3HE& p3
97 ) {
98
5/8
✓ Branch 0 taken 31 times.
✓ Branch 1 taken 23684 times.
✓ Branch 3 taken 25 times.
✓ Branch 4 taken 6 times.
✓ Branch 6 taken 25 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 23715 static PredicateStats stats("orient_3d(vec3HE)");
99 23715 stats.log_invoke();
100 // Filter
101 {
102 23715 interval_nt::Rounding rounding;
103
1/2
✓ Branch 1 taken 23715 times.
✗ Branch 2 not taken.
 23715 vec3HI p0I(p0);
104
2/4
✓ Branch 1 taken 23715 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23715 times.
✗ Branch 5 not taken.
 23715 vec3HI U = vec3HI(p1)-p0I;
105
2/4
✓ Branch 1 taken 23715 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23715 times.
✗ Branch 5 not taken.
 23715 vec3HI V = vec3HI(p2)-p0I;
106
2/4
✓ Branch 1 taken 23715 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23715 times.
✗ Branch 5 not taken.
 23715 vec3HI W = vec3HI(p3)-p0I;
107 23715 interval_nt::Sign2 s1 = U.w.sign();
108 23715 interval_nt::Sign2 s2 = V.w.sign();
109 23715 interval_nt::Sign2 s3 = W.w.sign();
110 23715 if(
111
1/2
✓ Branch 1 taken 23715 times.
✗ Branch 2 not taken.
 47430 interval_nt::sign_is_non_zero(s1) &&
112
3/6
✓ Branch 0 taken 23715 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 23715 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 23715 times.
✗ Branch 6 not taken.
 47430 interval_nt::sign_is_non_zero(s2) &&
113 23715 interval_nt::sign_is_non_zero(s3)
114 ) {
115
1/2
✓ Branch 1 taken 23715 times.
✗ Branch 2 not taken.
 23715 interval_nt Delta = det3x3(
116 U.x, U.y, U.z,
117 V.x, V.y, V.z,
118 W.x, W.y, W.z
119 );
120 23715 interval_nt::Sign2 s = Delta.sign();
121
2/2
✓ Branch 1 taken 16167 times.
✓ Branch 2 taken 7548 times.
 23715 if(interval_nt::sign_is_non_zero(s)) {
122 return Sign(
123
1/2
✓ Branch 1 taken 16167 times.
✗ Branch 2 not taken.
 16167 interval_nt::convert_sign(s)*
124
1/2
✓ Branch 1 taken 16167 times.
✗ Branch 2 not taken.
 16167 interval_nt::convert_sign(s1)*
125
1/2
✓ Branch 1 taken 16167 times.
✗ Branch 2 not taken.
 16167 interval_nt::convert_sign(s2)*
126
1/2
✓ Branch 1 taken 16167 times.
✗ Branch 2 not taken.
 16167 interval_nt::convert_sign(s3)
127 16167 );
128 }
129 }
130
2/2
✓ Branch 1 taken 7548 times.
✓ Branch 2 taken 16167 times.
 23715 }
131
132 7548 stats.log_exact();
133
134
1/2
✓ Branch 1 taken 7548 times.
✗ Branch 2 not taken.
 7548 vec3HE U = p1-p0;
135
1/2
✓ Branch 1 taken 7548 times.
✗ Branch 2 not taken.
 7548 vec3HE V = p2-p0;
136
1/2
✓ Branch 1 taken 7548 times.
✗ Branch 2 not taken.
 7548 vec3HE W = p3-p0;
137
138 // Here we do not use expansion_det3x3() (that
139 // allocates on the stack), because
140 // RadialSort uses generted points that
141 // can have very looonng expansions that
142 // can cause stack overflow.
143 expansion_nt Delta = det3x3(
144 U.x, U.y, U.z,
145 V.x, V.y, V.z,
146 W.x, W.y, W.z
147
1/2
✓ Branch 1 taken 7548 times.
✗ Branch 2 not taken.
 7548 );
148
149 Sign result = Sign(
150
1/2
✓ Branch 1 taken 7548 times.
✗ Branch 2 not taken.
 7548 Delta.sign()*
151
1/2
✓ Branch 2 taken 7548 times.
✗ Branch 3 not taken.
 7548 U.w.rep().sign()*
152
1/2
✓ Branch 2 taken 7548 times.
✗ Branch 3 not taken.
 7548 V.w.rep().sign()*
153
1/2
✓ Branch 2 taken 7548 times.
✗ Branch 3 not taken.
 7548 W.w.rep().sign()
154 7548 );
155
156 7548 return result;
157 7548 }
158
159 4321212 Sign orient_2d_projected(
160 const vec3HE& p0, const vec3HE& p1, const vec3HE& p2,
161 coord_index_t axis
162 ) {
163
5/8
✓ Branch 0 taken 111 times.
✓ Branch 1 taken 4321101 times.
✓ Branch 3 taken 34 times.
✓ Branch 4 taken 77 times.
✓ Branch 6 taken 34 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 4321212 static PredicateStats stats("orient_2d_projected(vec3HE)");
164 4321212 stats.log_invoke();
165
166 4321212 coord_index_t u = coord_index_t((axis+1)%3);
167 4321212 coord_index_t v = coord_index_t((axis+2)%3);
168
169 // Filter, using interval arithmetics
170 {
171 4321212 interval_nt::Rounding rounding;
172
173
1/2
✓ Branch 1 taken 4321212 times.
✗ Branch 2 not taken.
 4321212 interval_nt Delta = det3x3(
174
5/10
✓ Branch 1 taken 4321212 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4321212 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4321212 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4321212 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 4321212 times.
✗ Branch 14 not taken.
 4321212 interval_nt(p0[u]),interval_nt(p0[v]),interval_nt(p0.w),
175
5/10
✓ Branch 1 taken 4321212 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4321212 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4321212 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4321212 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 4321212 times.
✗ Branch 14 not taken.
 4321212 interval_nt(p1[u]),interval_nt(p1[v]),interval_nt(p1.w),
176
5/10
✓ Branch 1 taken 4321212 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4321212 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4321212 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4321212 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 4321212 times.
✗ Branch 14 not taken.
 4321212 interval_nt(p2[u]),interval_nt(p2[v]),interval_nt(p2.w)
177 );
178 4321212 interval_nt::Sign2 s = Delta.sign();
179
2/2
✓ Branch 1 taken 3419330 times.
✓ Branch 2 taken 901882 times.
 4321212 if(interval_nt::sign_is_determined(s)) {
180 return Sign(
181
1/2
✓ Branch 1 taken 3419330 times.
✗ Branch 2 not taken.
 3419330 interval_nt::convert_sign(s)*
182
3/6
✓ Branch 1 taken 3419330 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3419330 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3419330 times.
✗ Branch 8 not taken.
 3419330 p0.w.sign()*p1.w.sign()*p2.w.sign()
183 3419330 );
184 }
185
2/2
✓ Branch 1 taken 901882 times.
✓ Branch 2 taken 3419330 times.
 4321212 }
186
187 901882 stats.log_exact();
188
189 901882 Sign result = ZERO;
190 {
191 #ifdef GEO_HAS_BIG_STACK
192
14/30
✓ Branch 20 taken 901882 times.
✗ Branch 21 not taken.
✓ Branch 24 taken 901882 times.
✗ Branch 25 not taken.
✓ Branch 29 taken 901882 times.
✗ Branch 30 not taken.
✓ Branch 33 taken 901882 times.
✗ Branch 34 not taken.
✓ Branch 38 taken 901882 times.
✗ Branch 39 not taken.
✓ Branch 42 taken 901882 times.
✗ Branch 43 not taken.
✓ Branch 49 taken 901882 times.
✗ Branch 50 not taken.
✓ Branch 53 taken 901882 times.
✗ Branch 54 not taken.
✓ Branch 58 taken 901882 times.
✗ Branch 59 not taken.
✓ Branch 62 taken 901882 times.
✗ Branch 63 not taken.
✓ Branch 67 taken 901882 times.
✗ Branch 68 not taken.
✓ Branch 71 taken 901882 times.
✗ Branch 72 not taken.
✓ Branch 75 taken 901882 times.
✗ Branch 76 not taken.
✗ Branch 77 not taken.
✓ Branch 78 taken 901882 times.
✗ Branch 80 not taken.
✗ Branch 81 not taken.
 901882 const expansion& Delta = expansion_det3x3(
193 p0[u].rep(), p0[v].rep(), p0.w.rep(),
194 p1[u].rep(), p1[v].rep(), p1.w.rep(),
195 p2[u].rep(), p2[v].rep(), p2.w.rep()
196 );
197 #else
198 expansion_nt Delta = det3x3(
199 p0[u], p0[v], p0.w,
200 p1[u], p1[v], p1.w,
201 p2[u], p2[v], p2.w
202 );
203 #endif
204 901882 result = Sign(
205 901882 Delta.sign()*
206 901882 p0.w.rep().sign()*
207 901882 p1.w.rep().sign()*
208 901882 p2.w.rep().sign()
209 );
210 }
211 901882 return result;
212 }
213
214 ✗ Sign dot_2d(const vec2HE& p0, const vec2HE& p1, const vec2HE& p2) {
215 ✗ static PredicateStats stats("dot_2d(vec2HE)");
216 ✗ stats.log_invoke();
217
218 // TODO: filter
219
220 ✗ vec2HE U = p1 - p0;
221 ✗ vec2HE V = p2 - p0;
222 #ifdef GEO_HAS_BIG_STACK
223 ✗ const expansion& x1x2 = expansion_product(U.x.rep(), V.x.rep());
224 ✗ const expansion& y1y2 = expansion_product(U.y.rep(), V.y.rep());
225 ✗ const expansion& S = expansion_sum(x1x2, y1y2);
226 #else
227 expansion_nt S = U.x*V.x+U.y*V.y;
228 #endif
229 ✗ return Sign(S.sign()*U.w.sign()*V.w.sign());
230 ✗ }
231
232 /******************************************************************************/
233
234 /**
235 * \brief Computes the sign of
236 * det3x3(x1,y1,1,x2,y2,1,x3,y3,1)
237 * \param[in] p1 , p2 , p3 the three points in
238 * homogeneous exact coordiates
239 * \return the sign of the determinant
240 */
241 12696 static inline Sign det3_111_sign(
242 const vec2HE& p1,
243 const vec2HE& p2,
244 const vec2HE& p3
245 ) {
246
1/2
✓ Branch 1 taken 12696 times.
✗ Branch 2 not taken.
 12696 expansion_nt m1 = det2x2(p2.x, p2.y, p3.x, p3.y);
247
1/2
✓ Branch 1 taken 12696 times.
✗ Branch 2 not taken.
 12696 expansion_nt m2 = det2x2(p1.x, p1.y, p3.x, p3.y);
248
1/2
✓ Branch 1 taken 12696 times.
✗ Branch 2 not taken.
 12696 expansion_nt m3 = det2x2(p1.x, p1.y, p2.x, p2.y);
249
3/6
✓ Branch 1 taken 12696 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12696 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12696 times.
✗ Branch 8 not taken.
 12696 m1.optimize(); m2.optimize(); m3.optimize();
250
5/10
✓ Branch 1 taken 12696 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12696 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12696 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 12696 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 12696 times.
✗ Branch 14 not taken.
 12696 expansion_nt D = p1.w*m1-p2.w*m2+p3.w*m3 ;
251
4/8
✓ Branch 1 taken 12696 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12696 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12696 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 12696 times.
✗ Branch 11 not taken.
 25392 return Sign(p1.w.sign()*p2.w.sign()*p3.w.sign()*D.sign());
252 12696 }
253
254 342079 Sign incircle_2d_SOS_with_lengths(
255 const vec2HE& p0, const vec2HE& p1,
256 const vec2HE& p2, const vec2HE& p3,
257 double l0, double l1, double l2, double l3
258 ) {
259
4/8
✓ Branch 0 taken 35 times.
✓ Branch 1 taken 342044 times.
✓ Branch 3 taken 35 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 35 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 342079 static PredicateStats stats("incircle_2d_SOS_with_lengths(vec2HE)");
260 342079 stats.log_invoke();
261
262 342079 Sign result = ZERO;
263
264 // "Documentation is a love letter that you write to your
265 // future self." - Damian Conway (or that you write to the
266 // poor guy who will have one day to dive again in this stuff,
267 // but it's probably me anyway)
268 //
269 // Determinant to compute:
270 // | x0 y0 l0 1 |
271 // | x1 y1 l1 1 |
272 // | x2 y2 l2 1 |
273 // | x3 y3 l3 1 |
274 // where li = xi^2 + yi^2
275 // (positive if (p0,p1,p2) counterclockwise and p3 in circumcircle
276 // of (p0,p1,p2)). Sign changes if (p0,p1,p2) is clockwise).
277 // We suppose that the li's are *given* numbers (it is like
278 // perturbating a regular (weighted) triangulation instead of a
279 // Delaunay triangulation).
280 // It allows to use arithmetic expansions without
281 // overflowing/underflowing too soon.
282 //
283 // Subtract last row to first three rows
284 // (does not change determinant):
285 //
286 // | x0-x3 y0-y3 l0-l3 0 |
287 // | x1-x3 y1-y3 l1-l3 0 |
288 // | x2-x3 y2-y3 l2-l3 0 |
289 // | x3 y3 l3 1 |
290 //
291 // Develop along last column:
292 // | x0-x3 y0-y3 l0-l3 |
293 // | x1-x3 y1-y3 l1-l3 |
294 // | x2-x3 y2-y3 l2-l3 |
295 //
296 // let (Xi+1,Yi+1,Wi+1) = (xi,yi)-(x3,y3) in homogeneous coordinates
297 // let Li+1 = li-l3:
298 // | X1/W1 Y1/W1 L1 |
299 // | X2/W2 Y2/W2 L2 |
300 // | X3/W3 Y3/W3 L3 |
301 //
302 // Develop along last column, factor-out the Wi's
303 //
304 // | X2 Y2 | | X1 Y1 | | X1 Y1 |
305 // (L1/W2W3) | X3 Y3 | - (L2/W1W3) | X3 Y3 | + (L3/W1W2) | X2 Y2 |
306 //
307 // Multiply everything by W1W2W3:
308 //
309 // | X2 Y2 | | X1 Y1 | | X1 Y1 |
310 // sign( L1W1 | X3 Y3 | - L2W2 | X3 Y3 | + L3W3 | X2 Y2 | ) *
311 // sign(W1) * sign(W2) * sign(W3)
312
313 // The four approximated li's. It is OK since they will
314 // always have the same value for the same vertex.
315 // We do it like that because computing them exactly (and
316 // properly propagating the wi^2's) makes expansions
317 // overflow/underflow.
318 // It is like perturbating a regular (weighted) triangulation
319 // instead of a Delaunay triangulation.
320 // However, if incircle(p1,p2,p3,p4) is lower than 0, it does
321 // not imply that (p1,p2,p3,p4) forms a convex quadrilateral
322 // (needs to be tested in addition, it is what CDT2d does
323 // when exact_incircle_ is set to false).
324
325 // We could also compute them each time, as in incircle_2d_SOS()
326 // (but we are caching them in MeshSurfaceIntersection's temporary
327 // Vertex objects), this gains 20-25% performance so it is
328 // worth it.
329
330 // Filter
331 {
332 342079 interval_nt::Rounding rounding;
333
1/2
✓ Branch 1 taken 342079 times.
✗ Branch 2 not taken.
 342079 interval_nt l3I(l3);
334
2/4
✓ Branch 1 taken 342079 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 342079 times.
✗ Branch 5 not taken.
 342079 interval_nt L1 = interval_nt(l0) - l3I;
335
2/4
✓ Branch 1 taken 342079 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 342079 times.
✗ Branch 5 not taken.
 342079 interval_nt L2 = interval_nt(l1) - l3I;
336
2/4
✓ Branch 1 taken 342079 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 342079 times.
✗ Branch 5 not taken.
 342079 interval_nt L3 = interval_nt(l2) - l3I;
337
338
1/2
✓ Branch 1 taken 342079 times.
✗ Branch 2 not taken.
 342079 vec2HI p3I(p3);
339
2/4
✓ Branch 1 taken 342079 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 342079 times.
✗ Branch 5 not taken.
 342079 vec2HI P1 = vec2HI(p0) - p3I;
340
2/4
✓ Branch 1 taken 342079 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 342079 times.
✗ Branch 5 not taken.
 342079 vec2HI P2 = vec2HI(p1) - p3I;
341
2/4
✓ Branch 1 taken 342079 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 342079 times.
✗ Branch 5 not taken.
 342079 vec2HI P3 = vec2HI(p2) - p3I;
342
343 342079 interval_nt::Sign2 s1 = P1.w.sign();
344 342079 interval_nt::Sign2 s2 = P2.w.sign();
345 342079 interval_nt::Sign2 s3 = P3.w.sign();
346
347 342079 if(
348
1/2
✓ Branch 1 taken 342079 times.
✗ Branch 2 not taken.
 684158 interval_nt::sign_is_non_zero(s1) &&
349
3/6
✓ Branch 0 taken 342079 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 342079 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 342079 times.
✗ Branch 6 not taken.
 684158 interval_nt::sign_is_non_zero(s2) &&
350 342079 interval_nt::sign_is_non_zero(s3)
351 ) {
352
353
1/2
✓ Branch 1 taken 342079 times.
✗ Branch 2 not taken.
 342079 interval_nt M1 = det2x2(P2.x, P2.y, P3.x, P3.y);
354
1/2
✓ Branch 1 taken 342079 times.
✗ Branch 2 not taken.
 342079 interval_nt M2 = det2x2(P1.x, P1.y, P3.x, P3.y);
355
1/2
✓ Branch 1 taken 342079 times.
✗ Branch 2 not taken.
 342079 interval_nt M3 = det2x2(P1.x, P1.y, P2.x, P2.y);
356
357
2/4
✓ Branch 1 taken 342079 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 342079 times.
✗ Branch 5 not taken.
 342079 interval_nt D = L1*P1.w*M1
358
3/6
✓ Branch 1 taken 342079 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 342079 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 342079 times.
✗ Branch 8 not taken.
 684158 - L2*P2.w*M2
359
3/6
✓ Branch 1 taken 342079 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 342079 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 342079 times.
✗ Branch 8 not taken.
 684158 + L3*P3.w*M3 ;
360
361 342079 interval_nt::Sign2 s = D.sign();
362
2/2
✓ Branch 1 taken 306125 times.
✓ Branch 2 taken 35954 times.
 342079 if(interval_nt::sign_is_non_zero(s)) {
363 return Sign(
364
1/2
✓ Branch 1 taken 306125 times.
✗ Branch 2 not taken.
 306125 interval_nt::convert_sign(s) *
365
1/2
✓ Branch 1 taken 306125 times.
✗ Branch 2 not taken.
 306125 interval_nt::convert_sign(s1) *
366
1/2
✓ Branch 1 taken 306125 times.
✗ Branch 2 not taken.
 306125 interval_nt::convert_sign(s2) *
367
1/2
✓ Branch 1 taken 306125 times.
✗ Branch 2 not taken.
 306125 interval_nt::convert_sign(s3)
368 306125 );
369 }
370 }
371
2/2
✓ Branch 1 taken 35954 times.
✓ Branch 2 taken 306125 times.
 342079 }
372
373 // Exact
374 35954 stats.log_exact();
375 {
376 // These ones can be computed on the stack even
377 // under MacOSX since they are at most of length 2
378
1/2
✓ Branch 1 taken 35954 times.
✗ Branch 2 not taken.
 35954 expansion_nt L1(expansion_nt::DIFF, l0, l3);
379
1/2
✓ Branch 1 taken 35954 times.
✗ Branch 2 not taken.
 35954 expansion_nt L2(expansion_nt::DIFF, l1, l3);
380
1/2
✓ Branch 1 taken 35954 times.
✗ Branch 2 not taken.
 35954 expansion_nt L3(expansion_nt::DIFF, l2, l3);
381
3/6
✓ Branch 1 taken 35954 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 35954 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 35954 times.
✗ Branch 8 not taken.
 35954 L1.optimize(); L2.optimize(); L3.optimize();
382
383
1/2
✓ Branch 1 taken 35954 times.
✗ Branch 2 not taken.
 35954 vec2HE P1 = p0 - p3;
384
1/2
✓ Branch 1 taken 35954 times.
✗ Branch 2 not taken.
 35954 vec2HE P2 = p1 - p3;
385
1/2
✓ Branch 1 taken 35954 times.
✗ Branch 2 not taken.
 35954 vec2HE P3 = p2 - p3;
386
3/6
✓ Branch 1 taken 35954 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 35954 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 35954 times.
✗ Branch 8 not taken.
 35954 P1.optimize(); P2.optimize(); P3.optimize();
387
388
389
1/2
✓ Branch 1 taken 35954 times.
✗ Branch 2 not taken.
 35954 expansion_nt M1 = det2x2(P2.x, P2.y, P3.x, P3.y);
390
1/2
✓ Branch 1 taken 35954 times.
✗ Branch 2 not taken.
 35954 expansion_nt M2 = det2x2(P1.x, P1.y, P3.x, P3.y);
391
1/2
✓ Branch 1 taken 35954 times.
✗ Branch 2 not taken.
 35954 expansion_nt M3 = det2x2(P1.x, P1.y, P2.x, P2.y);
392
3/6
✓ Branch 1 taken 35954 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 35954 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 35954 times.
✗ Branch 8 not taken.
 35954 M1.optimize(); M2.optimize(); M3.optimize();
393
394
2/4
✓ Branch 1 taken 35954 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 35954 times.
✗ Branch 5 not taken.
 71908 expansion_nt D = L1*P1.w*M1
395
3/6
✓ Branch 1 taken 35954 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 35954 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 35954 times.
✗ Branch 8 not taken.
 143816 - L2*P2.w*M2
396
3/6
✓ Branch 1 taken 35954 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 35954 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 35954 times.
✗ Branch 8 not taken.
 107862 + L3*P3.w*M3 ;
397
398
4/8
✓ Branch 1 taken 35954 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 35954 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 35954 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 35954 times.
✗ Branch 11 not taken.
 35954 result = Sign(D.sign()*P1.w.sign()*P2.w.sign()*P3.w.sign());
399 35954 }
400
401
2/2
✓ Branch 0 taken 23258 times.
✓ Branch 1 taken 12696 times.
 35954 if(result != ZERO) {
402 23258 return result;
403 }
404
405 12696 stats.log_SOS();
406
407 // Symbolic perturbation.
408 //
409 // We use the simple form of the predicate:
410 // | x0 y0 (x0^2+y0^2+eps^i0) 1 |
411 // | x1 y1 (x1^2+y1^2+eps^i1) 1 |
412 // | x2 y2 (x2^2+y2^2+eps^i2) 1 |
413 // | x3 y3 (x3^2+y3^2+eps^i3) 1 |
414 // where i0,i1,i2,i3 denote the indices of the points (here, they
415 // are local indices, coming from geometric sorting, lexico order)
416 // Develop along the third row (keeping only the terms in epsilon):
417 // | x1 y1 1 | | x0 y0 1 |
418 // eps^i0 | x2 y2 1 | - eps^i1 | x2 y2 1 |
419 // | x3 y3 1 | | x3 y3 1 |
420 //
421 // | x0 y0 1 | | x0 y0 1 |
422 // + eps^i2 | x1 y1 1 | - eps^i3 | x1 y1 1 |
423 // | x3 y3 1 | | x2 y2 1 |
424 12696 return SOS(
425 vec2HgLexicoCompare<expansion_nt>(),
426 5660 p0, SOS_result( det3_111_sign(p1,p2,p3)),
427 2252 p1, SOS_result(-det3_111_sign(p0,p2,p3)),
428 3422 p2, SOS_result( det3_111_sign(p0,p1,p3)),
429 1362 p3, SOS_result(-det3_111_sign(p0,p1,p2))
430 12696 );
431 }
432
433 724326 coord_index_t triangle_normal_axis(
434 const vec3& p1, const vec3& p2, const vec3& p3
435 ) {
436
5/8
✓ Branch 0 taken 35 times.
✓ Branch 1 taken 724291 times.
✓ Branch 3 taken 34 times.
✓ Branch 4 taken 1 times.
✓ Branch 6 taken 34 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 724326 static PredicateStats stats("triangle_normal_axis");
437 724326 stats.log_invoke();
438
439 // Filter using interval arithmetics
440 {
441 724326 interval_nt::Rounding rounding;
442
1/2
✓ Branch 1 taken 724326 times.
✗ Branch 2 not taken.
 724326 vec3I p1I(p1);
443
2/4
✓ Branch 1 taken 724326 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 724326 times.
✗ Branch 5 not taken.
 724326 vec3I U = vec3I(p2) - p1I;
444
2/4
✓ Branch 1 taken 724326 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 724326 times.
✗ Branch 5 not taken.
 724326 vec3I V = vec3I(p3) - p1I;
445
1/2
✓ Branch 1 taken 724326 times.
✗ Branch 2 not taken.
 724326 vec3I N = cross(U,V);
446 724326 interval_nt::Sign2 sx = N.x.sign();
447 724326 interval_nt::Sign2 sy = N.y.sign();
448 724326 interval_nt::Sign2 sz = N.z.sign();
449 724326 if(
450 724326 !interval_nt::sign_is_determined(sx) ||
451
6/6
✓ Branch 0 taken 196996 times.
✓ Branch 1 taken 527330 times.
✓ Branch 3 taken 165999 times.
✓ Branch 4 taken 30997 times.
✓ Branch 5 taken 577059 times.
✓ Branch 6 taken 147267 times.
 890325 !interval_nt::sign_is_determined(sy) ||
452
2/2
✓ Branch 1 taken 18732 times.
✓ Branch 2 taken 147267 times.
 165999 !interval_nt::sign_is_determined(sz)
453 ) {
454 577059 goto exact; // Yes, goto, why not ?
455 }
456
457
3/4
✓ Branch 1 taken 147267 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 74849 times.
✓ Branch 4 taken 72418 times.
 147267 if(interval_nt::convert_sign(sx) != POSITIVE) {
458
1/2
✓ Branch 1 taken 74849 times.
✗ Branch 2 not taken.
 74849 N.x.negate();
459 }
460
3/4
✓ Branch 1 taken 147267 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 72306 times.
✓ Branch 4 taken 74961 times.
 147267 if(interval_nt::convert_sign(sy) != POSITIVE) {
461
1/2
✓ Branch 1 taken 72306 times.
✗ Branch 2 not taken.
 72306 N.y.negate();
462 }
463
3/4
✓ Branch 1 taken 147267 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 65655 times.
✓ Branch 4 taken 81612 times.
 147267 if(interval_nt::convert_sign(sz) != POSITIVE) {
464
1/2
✓ Branch 1 taken 65655 times.
✗ Branch 2 not taken.
 65655 N.z.negate();
465 }
466
1/2
✓ Branch 1 taken 147267 times.
✗ Branch 2 not taken.
 147267 interval_nt::Sign2 sxy = (N.x - N.y).sign();
467
1/2
✓ Branch 1 taken 147267 times.
✗ Branch 2 not taken.
 147267 interval_nt::Sign2 sxz = (N.x - N.z).sign();
468 147267 if(
469
4/4
✓ Branch 1 taken 147243 times.
✓ Branch 2 taken 24 times.
✓ Branch 3 taken 146 times.
✓ Branch 4 taken 147121 times.
 294510 !interval_nt::sign_is_determined(sxy) ||
470
2/2
✓ Branch 1 taken 122 times.
✓ Branch 2 taken 147121 times.
 147243 !interval_nt::sign_is_determined(sxz)
471 ) {
472 146 goto exact; // Ahaha, another one !!!
473 }
474 147121 if(
475
5/6
✓ Branch 1 taken 147121 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 73640 times.
✓ Branch 4 taken 73481 times.
✓ Branch 5 taken 49926 times.
✓ Branch 6 taken 97195 times.
 220761 interval_nt::convert_sign(sxy) >= 0 &&
476
3/4
✓ Branch 1 taken 73640 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 49926 times.
✓ Branch 4 taken 23714 times.
 73640 interval_nt::convert_sign(sxz) >= 0
477 ) {
478 49926 return 0;
479 }
480
1/2
✓ Branch 1 taken 97195 times.
✗ Branch 2 not taken.
 97195 interval_nt::Sign2 syz = (N.y - N.z).sign();
481
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 97195 times.
 97195 if(!interval_nt::sign_is_determined(syz)) {
482 ✗ goto exact; // The last one (for now !)
483 }
484
3/4
✓ Branch 1 taken 97195 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 50210 times.
✓ Branch 4 taken 46985 times.
 97195 if(interval_nt::convert_sign(syz) >=0 ) {
485 50210 return 1;
486 }
487 46985 return 2;
488
2/2
✓ Branch 1 taken 577205 times.
✓ Branch 2 taken 147121 times.
 724326 }
489
490 // Exact computation, using low-level expansion API
491 // (expansions allocated on the stack, better for
492 // multithreading)
493 577205 exact:
494 577205 stats.log_exact();
495
496 // These ones can be computed on the stack even
497 // under MacOSX since they are at most of length 2
498
499
3/8
✓ Branch 4 taken 577205 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 577205 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 577205 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 577205 const expansion& Ux = expansion_diff(p2.x, p1.x);
500
3/8
✓ Branch 4 taken 577205 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 577205 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 577205 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 577205 const expansion& Uy = expansion_diff(p2.y, p1.y);
501
3/8
✓ Branch 4 taken 577205 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 577205 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 577205 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 577205 const expansion& Uz = expansion_diff(p2.z, p1.z);
502
503
3/8
✓ Branch 4 taken 577205 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 577205 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 577205 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 577205 const expansion& Vx = expansion_diff(p3.x, p1.x);
504
3/8
✓ Branch 4 taken 577205 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 577205 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 577205 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 577205 const expansion& Vy = expansion_diff(p3.y, p1.y);
505
3/8
✓ Branch 4 taken 577205 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 577205 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 577205 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 577205 const expansion& Vz = expansion_diff(p3.z, p1.z);
506
507
2/6
✓ Branch 6 taken 577205 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 577205 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 577205 expansion& Nx = expansion_det2x2(Uy,Vy,Uz,Vz);
508
2/6
✓ Branch 6 taken 577205 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 577205 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 577205 expansion& Ny = expansion_det2x2(Uz,Vz,Ux,Vx);
509
2/6
✓ Branch 6 taken 577205 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 577205 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 577205 expansion& Nz = expansion_det2x2(Ux,Vx,Uy,Vy);
510
511
2/2
✓ Branch 1 taken 552536 times.
✓ Branch 2 taken 24669 times.
 577205 if(Nx.sign() != POSITIVE) {
512 552536 Nx.negate();
513 }
514
515
2/2
✓ Branch 1 taken 550580 times.
✓ Branch 2 taken 26625 times.
 577205 if(Ny.sign() != POSITIVE) {
516 550580 Ny.negate();
517 }
518
519
2/2
✓ Branch 1 taken 422659 times.
✓ Branch 2 taken 154546 times.
 577205 if(Nz.sign() != POSITIVE) {
520 422659 Nz.negate();
521 }
522
523
6/6
✓ Branch 1 taken 535472 times.
✓ Branch 2 taken 41733 times.
✓ Branch 4 taken 41274 times.
✓ Branch 5 taken 494198 times.
✓ Branch 6 taken 41274 times.
✓ Branch 7 taken 535931 times.
 577205 if(Nx.compare(Ny) >= 0 && Nx.compare(Nz) >= 0) {
524
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 41274 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 41274 geo_debug_assert(Nx.sign() != ZERO);
525 41274 return 0;
526 }
527
2/2
✓ Branch 1 taken 37945 times.
✓ Branch 2 taken 497986 times.
 535931 if(Ny.compare(Nz) >= 0) {
528
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 37945 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 37945 geo_debug_assert(Ny.sign() != ZERO);
529 37945 return 1;
530 }
531
1/6
✗ Branch 1 not taken.
✓ Branch 2 taken 497986 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 497986 geo_assert(Nz.sign() != ZERO);
532 497986 return 2;
533 }
534
535 ✗ bool aligned_3d(
536 const vec3HE& p0, const vec3HE& p1, const vec3HE& p2
537 ) {
538 // TODO: filter if need be
539 ✗ vec3HE U = p1-p0;
540 ✗ vec3HE V = p2-p0;
541 return (
542 ✗ det2x2(U.x,V.x,U.y,V.y).sign() == ZERO &&
543 ✗ det2x2(U.y,V.y,U.z,V.z).sign() == ZERO &&
544 ✗ det2x2(U.z,V.z,U.x,V.x).sign() == ZERO
545 ✗ );
546 ✗ }
547
548 246914 bool on_segment_3d(
549 const vec3HE& p, const vec3HE& q1, const vec3HE& q2
550 ) {
551 // TODO: filter if need be
552
1/2
✓ Branch 1 taken 246914 times.
✗ Branch 2 not taken.
 246914 vec3HE U = p-q1;
553
1/2
✓ Branch 1 taken 246914 times.
✗ Branch 2 not taken.
 246914 vec3HE V = p-q2;
554 246914 if (
555
2/4
✓ Branch 1 taken 238471 times.
✓ Branch 2 taken 8443 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 246914 det2x2(U.x,V.x,U.y,V.y).sign() != ZERO ||
556
9/16
✓ Branch 1 taken 246914 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 246914 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 19962 times.
✓ Branch 7 taken 226952 times.
✓ Branch 9 taken 19962 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 19962 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 11583 times.
✓ Branch 15 taken 8379 times.
✓ Branch 17 taken 246914 times.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 505411 det2x2(U.y,V.y,U.z,V.z).sign() != ZERO ||
557
8/12
✓ Branch 1 taken 11583 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11583 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 3140 times.
✓ Branch 7 taken 8443 times.
✓ Branch 8 taken 11583 times.
✓ Branch 9 taken 235331 times.
✓ Branch 11 taken 19962 times.
✓ Branch 12 taken 226952 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 258497 det2x2(U.z,V.z,U.x,V.x).sign() != ZERO
558 ) {
559 238471 return false;
560 }
561
562 return (
563 (
564
6/12
✓ Branch 1 taken 8443 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8443 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8443 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 8443 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 8443 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 8443 times.
✗ Branch 17 not taken.
 16886 (U.x*V.x + U.y*V.y + U.z*V.z).sign() *
565
2/4
✓ Branch 1 taken 8443 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8443 times.
✗ Branch 5 not taken.
 8443 U.w.sign() * V.w.sign()
566 ) <= ZERO
567 8443 );
568 246914 }
569
570
571 24895 vec3 approximate(const vec3HE& p) {
572 // TODO: find a way of computing the round to nearest approxomation.
573 // see division operation for expansions,
574 // here:
575 // https://www.jucs.org/jucs_5_6/division_of_floating_point/Daumas_M.pdf
576 24895 double w = p.w.estimate();
577
3/6
✓ Branch 1 taken 24895 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24895 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 24895 times.
✗ Branch 8 not taken.
 24895 return vec3(p.x.estimate()/w, p.y.estimate()/w, p.z.estimate()/w);
578 }
579
580 81 vec2 approximate(const vec2HE& p) {
581 // TODO: find a way of computing the round to nearest approxomation.
582 // see division operation for expansions,
583 // here:
584 // https://www.jucs.org/jucs_5_6/division_of_floating_point/Daumas_M.pdf
585 81 double w = p.w.estimate();
586
2/4
✓ Branch 1 taken 81 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 81 times.
✗ Branch 5 not taken.
 81 return vec2(p.x.estimate()/w, p.y.estimate()/w);
587 }
588
589 }
590
591 /*****************************************************************/
592
593 // Under Linux we got 10 Mb of stack (!) Then some operations can be
594 // made faster by using the low-level expansion API (that allocates
595 // intermediary multiprecision values on stack rather than in the heap).
596 // These optimized functions are written as template specializations
597 // (used automatically).
598
599 #ifdef GEO_HAS_BIG_STACK
600
601 310518 template<> expansion_nt det(const vec2E& v1, const vec2E& v2) {
602 310518 expansion* result = expansion::new_expansion_on_heap(
603 expansion::det2x2_capacity(
604 v1.x.rep(), v1.y.rep(),
605 v2.x.rep(), v2.y.rep()
606 )
607 );
608 310518 result->assign_det2x2(
609 v1.x.rep(), v1.y.rep(),
610 v2.x.rep(), v2.y.rep()
611 );
612 310518 return expansion_nt(result);
613 }
614
615 ✗ template<> expansion_nt dot(const vec2E& v1, const vec2E& v2) {
616 ✗ const expansion& m1 = expansion_product(v1.x.rep(), v2.x.rep());
617 ✗ const expansion& m2 = expansion_product(v1.y.rep(), v2.y.rep());
618 ✗ return expansion_nt(expansion_nt::SUM, m1, m2);
619 }
620
621 180391 template<> expansion_nt dot(const vec3E& v1, const vec3E& v2) {
622
2/6
✓ Branch 12 taken 180391 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 180391 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
 180391 const expansion& m1 = expansion_product(v1.x.rep(), v2.x.rep());
623
2/6
✓ Branch 12 taken 180391 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 180391 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
 180391 const expansion& m2 = expansion_product(v1.y.rep(), v2.y.rep());
624
2/6
✓ Branch 12 taken 180391 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 180391 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
 180391 const expansion& m3 = expansion_product(v1.z.rep(), v2.z.rep());
625 180391 return expansion_nt(expansion_nt::SUM,m1,m2,m3);
626 }
627
628 /*********************************************/
629
630 222821 template<> vec3Hg<expansion_nt> mix(
631 const rationalg<expansion_nt>& t, const vec3& p1, const vec3& p2
632 ) {
633 222821 expansion& st_d = const_cast<expansion&>(t.denom().rep());
634 222821 st_d.optimize();
635 222821 expansion& t_n = const_cast<expansion&>(t.num().rep());
636 222821 t_n.optimize();
637
2/6
✓ Branch 6 taken 222821 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 222821 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 222821 const expansion& s_n = expansion_diff(st_d, t_n);
638
3/8
✓ Branch 4 taken 222821 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 222821 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 222821 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 222821 const expansion& sx = expansion_product(s_n, p1.x);
639
3/8
✓ Branch 4 taken 222821 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 222821 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 222821 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 222821 const expansion& tx = expansion_product(t_n, p2.x);
640
3/8
✓ Branch 4 taken 222821 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 222821 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 222821 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 222821 const expansion& sy = expansion_product(s_n, p1.y);
641
3/8
✓ Branch 4 taken 222821 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 222821 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 222821 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 222821 const expansion& ty = expansion_product(t_n, p2.y);
642
3/8
✓ Branch 4 taken 222821 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 222821 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 222821 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 222821 const expansion& sz = expansion_product(s_n, p1.z);
643
3/8
✓ Branch 4 taken 222821 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 222821 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 222821 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 222821 const expansion& tz = expansion_product(t_n, p2.z);
644 return vec3HE(
645
1/2
✓ Branch 1 taken 222821 times.
✗ Branch 2 not taken.
 445642 expansion_nt(expansion_nt::SUM, sx,tx),
646
1/2
✓ Branch 1 taken 222821 times.
✗ Branch 2 not taken.
 445642 expansion_nt(expansion_nt::SUM, sy,ty),
647
1/2
✓ Branch 1 taken 222821 times.
✗ Branch 2 not taken.
 445642 expansion_nt(expansion_nt::SUM, sz,tz),
648
1/2
✓ Branch 1 taken 222821 times.
✗ Branch 2 not taken.
 445642 expansion_nt(st_d)
649
1/2
✓ Branch 1 taken 222821 times.
✗ Branch 2 not taken.
 445642 );
650 }
651
652 ✗ template<> vec2Hg<expansion_nt> mix(
653 const rationalg<expansion_nt>& t, const vec2& p1, const vec2& p2
654 ) {
655 ✗ expansion& st_d = const_cast<expansion&>(t.denom().rep());
656 ✗ st_d.optimize();
657 ✗ expansion& t_n = const_cast<expansion&>(t.num().rep());
658 ✗ t_n.optimize();
659 ✗ const expansion& s_n = expansion_diff(st_d, t_n);
660 ✗ const expansion& sx = expansion_product(s_n, p1.x);
661 ✗ const expansion& tx = expansion_product(t_n, p2.x);
662 ✗ const expansion& sy = expansion_product(s_n, p1.y);
663 ✗ const expansion& ty = expansion_product(t_n, p2.y);
664 return vec2HE(
665 ✗ expansion_nt(expansion_nt::SUM, sx,tx),
666 ✗ expansion_nt(expansion_nt::SUM, sy,ty),
667 ✗ expansion_nt(st_d)
668 ✗ );
669 }
670
671 /*********************************************/
672
673 17052 template <> vec3E triangle_normal<vec3E>(
674 const vec3& p1, const vec3& p2, const vec3& p3
675 ) {
676
3/8
✓ Branch 4 taken 17052 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 17052 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 17052 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 17052 const expansion& Ux = expansion_diff(p2.x,p1.x);
677
3/8
✓ Branch 4 taken 17052 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 17052 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 17052 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 17052 const expansion& Uy = expansion_diff(p2.y,p1.y);
678
3/8
✓ Branch 4 taken 17052 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 17052 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 17052 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 17052 const expansion& Uz = expansion_diff(p2.z,p1.z);
679
3/8
✓ Branch 4 taken 17052 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 17052 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 17052 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 17052 const expansion& Vx = expansion_diff(p3.x,p1.x);
680
3/8
✓ Branch 4 taken 17052 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 17052 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 17052 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 17052 const expansion& Vy = expansion_diff(p3.y,p1.y);
681
3/8
✓ Branch 4 taken 17052 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 17052 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 17052 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
 17052 const expansion& Vz = expansion_diff(p3.z,p1.z);
682 17052 expansion* Nx = expansion::new_expansion_on_heap(
683 expansion::det2x2_capacity(Uy,Uz,Vy,Vz)
684 );
685 17052 Nx->assign_det2x2(Uy,Uz,Vy,Vz);
686 17052 expansion* Ny = expansion::new_expansion_on_heap(
687 expansion::det2x2_capacity(Uz,Ux,Vz,Vx)
688 );
689 17052 Ny->assign_det2x2(Uz,Ux,Vz,Vx);
690 17052 expansion* Nz = expansion::new_expansion_on_heap(
691 expansion::det2x2_capacity(Ux,Uy,Vx,Vy)
692 );
693 17052 Nz->assign_det2x2(Ux,Uy,Vx,Vy);
694
1/2
✓ Branch 4 taken 17052 times.
✗ Branch 5 not taken.
 17052 return vec3E(expansion_nt(Nx), expansion_nt(Ny), expansion_nt(Nz));
695 }
696 #endif
697 }
698
699 /**************************************************************************/
700
701 #ifndef GEOGRAM_PSM
702 namespace GEO {
703 namespace PCK {
704 ✗ Sign orient_3d_SOS(
705 const exact::vec3h& p0, const exact::vec3h& p1,
706 const exact::vec3h& p2, const exact::vec3h& p3
707 ) {
708
709 struct SOS {
710 enum {W=3};
711
712 ✗ SOS(
713 const exact::vec3h& p0, const exact::vec3h& p1,
714 const exact::vec3h& p2, const exact::vec3h& p3
715 ✗ ) : p_orig{&p0, &p1, &p2, &p3}, p_sort{&p0, &p1, &p2, &p3} {
716 ✗ std::sort(
717 ✗ p_sort, p_sort+4,
718 ✗ [](
719 const exact::vec3h* pp1, const exact::vec3h* pp2
720 )->bool {
721 vec3HgLexicoCompare<exact::scalar> cmp;
722 ✗ return cmp(*pp1,*pp2);
723 }
724 );
725 ✗ parity = Permutation::permutation_is_odd(p_orig, p_sort, 4)
726 ✗ ? NEGATIVE : POSITIVE;
727 ✗ }
728
729 ✗ Sign orient_1d(index_t i, index_t j, index_t axis) const {
730 ✗ coord_index_t ax = coord_index_t(axis);
731 ✗ Sign s = geo_cmp(
732 ✗ exact::rational((*p_sort[i])[ax], (*p_sort[i])[W]),
733 ✗ exact::rational((*p_sort[j])[ax], (*p_sort[j])[W])
734 );
735 ✗ return Sign(s*parity);
736 }
737
738 ✗ Sign orient_2d(
739 index_t i, index_t j, index_t k, index_t axis1, index_t axis2
740 ) const {
741 ✗ coord_index_t ax1 = coord_index_t(axis1);
742 ✗ coord_index_t ax2 = coord_index_t(axis2);
743 exact::vec2h pi{
744 ✗ (*p_sort[i])[ax1], (*p_sort[i])[ax2], (*p_sort[i])[W]
745 ✗ };
746 exact::vec2h pj{
747 ✗ (*p_sort[j])[ax1], (*p_sort[j])[ax2], (*p_sort[j])[W]
748 ✗ };
749 exact::vec2h pk{
750 ✗ (*p_sort[k])[ax1], (*p_sort[k])[ax2], (*p_sort[k])[W]
751 ✗ };
752 ✗ Sign s = PCK::orient_2d(pi,pj,pk);
753 ✗ return Sign(s*parity);
754 ✗ }
755
756 const exact::vec3h* p_orig[4];
757 const exact::vec3h* p_sort[4];
758 Sign parity;
759 };
760
761 ✗ return orient_3d_SOS_impl<exact::vec3h, SOS>(p0,p1,p2,p3);
762 }
763 }
764 }
765 #endif
766