GCC Code Coverage Report

Directory: ./
File: numerics/multi_precision.cpp
Date: 2026-09-27 03:22:43
Exec Total Coverage
Lines: 352 475 74.1%
Functions: 23 36 63.9%
Branches: 264 802 32.9%
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/multi_precision.h>
49 #include <geogram/numerics/PCK.h>
50 #include <geogram/basic/process.h>
51 #include <geogram/basic/logger.h>
52
53 namespace {
54
55 using namespace GEO;
56
57 /************************************************************************/
58
59 /**
60 * \brief Computes the sum of a length 2 expansion and a double
61 * into a length 3 expansion.
62 * \param[in] a1 high-magnitude component of first argument
63 * \param[in] a0 low-magnitude component of first argument
64 * \param[in] b second argument
65 * \param[in] x2 high-magnitude component of the result
66 * \param[in] x1 component of the result
67 * \param[in] x0 low-magnitude component of the result
68 * \details By Jonathan Shewchuk.
69 */
70 18158202 inline void two_one_sum(
71 double a1, double a0, double b, double& x2, double& x1, double& x0
72 ) {
73 double _i;
74 18158202 two_sum(a0, b, _i, x0);
75 18158202 two_sum(a1, _i, x2, x1);
76 18158202 }
77
78 /**
79 * \brief Computes the sum of two length 2 expansions
80 * into a length 4 expansion.
81 * \param[in] a1 high-magnitude component of first argument
82 * \param[in] a0 low-magnitude component of first argument
83 * \param[in] b1 high-magnitude component of second argument
84 * \param[in] b0 high-magnitude component of second argument
85 * \param[in] x3 high-magnitude component of the result
86 * \param[in] x2 component of the result
87 * \param[in] x1 component of the result
88 * \param[in] x0 low-magnitude component of the result
89 * \details By Jonathan Shewchuk.
90 */
91 6052734 inline void two_two_sum(
92 double a1, double a0, double b1, double b0,
93 double& x3, double& x2, double& x1, double& x0
94 ) {
95 double _j, _0;
96 6052734 two_one_sum(a1, a0, b0, _j, _0, x0);
97 6052734 two_one_sum(_j, _0, b1, x3, x2, x1);
98 6052734 }
99
100 #ifndef FP_FAST_FMA
101
102 /**
103 * \brief Computes the product between two doubles where
104 * the second one have already been split.
105 * \param[in] a first argument
106 * \param[in] b second argument
107 * \param[in] bhi high-magnitude part of second argument
108 * \param[in] blo low-magnitude part of second argument
109 * \param[out] x high-magnitude component of the result
110 * \param[out] y low-magnitude component of the result
111 * \details By Jonathan Shewchuk.
112 */
113 inline void two_product_presplit(
114 double a, double b, double bhi, double blo, double& x, double& y
115 ) {
116 x = a * b;
117 double ahi;
118 double alo;
119 split(a, ahi, alo);
120 double err1 = x - (ahi * bhi);
121 double err2 = err1 - (alo * bhi);
122 double err3 = err2 - (ahi * blo);
123 y = (alo * blo) - err3;
124 }
125
126 /**
127 * \brief Computes the product between two doubles
128 * where both have already been split.
129 * \param[in] a first argument
130 * \param[in] ahi high-magnitude part of first argument
131 * \param[in] alo low-magnitude part of first argument
132 * \param[in] b second argument
133 * \param[in] bhi high-magnitude part of second argument
134 * \param[in] blo low-magnitude part of second argument
135 * \param[out] x high-magnitude component of the result
136 * \param[out] y low-magnitude component of the result
137 * \details By Jonathan Shewchuk.
138 */
139 inline void two_product_2presplit(
140 double a, double ahi, double alo,
141 double b, double bhi, double blo,
142 double& x, double& y
143 ) {
144 x = a * b;
145 double err1 = x - (ahi * bhi);
146 double err2 = err1 - (alo * bhi);
147 double err3 = err2 - (ahi * blo);
148 y = (alo * blo) - err3;
149 }
150
151 #endif
152
153 /**
154 * \brief Computes the square of an expansion of length 2.
155 * \param[in] a1 high-magnitude component of the argument
156 * \param[in] a0 low-magnitude component of the argument
157 * \param[out] x an array of six doubles to store the result.
158 * \details By Jonathan Shewchuk.
159 * An expansion of length two can be squared more quickly than finding the
160 * product of two different expansions of length two, and the result is
161 * guaranteed to have no more than six (rather than eight) components.
162 */
163 6052734 inline void two_square(
164 double a1, double a0,
165 double* x
166 ) {
167 double _0, _1, _2;
168 double _j, _k, _l;
169 6052734 square(a0, _j, x[0]);
170 6052734 _0 = a0 + a0;
171 6052734 two_product(a1, _0, _k, _1);
172 6052734 two_one_sum(_k, _1, _j, _l, _2, x[1]);
173 6052734 square(a1, _j, _1);
174 6052734 two_two_sum(_j, _1, _l, _2, x[5], x[4], x[3], x[2]);
175 6052734 }
176
177 /**
178 * \brief Computes the product of two expansions of length 2.
179 * \param[in] a first argument (array of 2 doubles)
180 * \param[in] b second argument (array of 2 doubles)
181 * \param[out] x an array of 8 doubles to store the result
182 * \details By Jonathan Shewchuk.
183 */
184 69593685 void two_two_product(
185 const double* a,
186 const double* b,
187 double* x
188 ) {
189 double _0, _1, _2;
190 double _i, _j, _k, _l, _m, _n;
191
192 // If the target processor supports the FMA (Fused Multiply Add)
193 // instruction, then the product of two doubles into a length-2
194 // expansion can be implemented as follows. Thanks to Marc Glisse
195 // for the information.
196 // Note: under gcc, automatic generations of fma() for a*b+c needs
197 // to be deactivated, using -ffp-contract=off, else it may break
198 // other functions such as fast_expansion_sum_zeroelim().
199 #ifdef FP_FAST_FMA
200 69593685 two_product(a[0],b[0],_i,x[0]);
201 69593685 two_product(a[1],b[0],_j,_0);
202 69593685 two_sum(_i, _0, _k, _1);
203 69593685 fast_two_sum(_j, _k, _l, _2);
204 69593685 two_product(a[0], b[1], _i, _0);
205 69593685 two_sum(_1, _0, _k, x[1]);
206 69593685 two_sum(_2, _k, _j, _1);
207 69593685 two_sum(_l, _j, _m, _2);
208 69593685 two_product(a[1], b[1], _j, _0);
209 69593685 two_sum(_i, _0, _n, _0);
210 69593685 two_sum(_1, _0, _i, x[2]);
211 69593685 two_sum(_2, _i, _k, _1);
212 69593685 two_sum(_m, _k, _l, _2);
213 69593685 two_sum(_j, _n, _k, _0);
214 69593685 two_sum(_1, _0, _j, x[3]);
215 69593685 two_sum(_2, _j, _i, _1);
216 69593685 two_sum(_l, _i, _m, _2);
217 69593685 two_sum(_1, _k, _i, x[4]);
218 69593685 two_sum(_2, _i, _k, x[5]);
219 69593685 two_sum(_m, _k, x[7], x[6]);
220 #else
221 double a0hi, a0lo;
222 split(a[0], a0hi, a0lo);
223 double bhi, blo;
224 split(b[0], bhi, blo);
225 two_product_2presplit(
226 a[0], a0hi, a0lo, b[0], bhi, blo, _i, x[0]
227 );
228 double a1hi, a1lo;
229 split(a[1], a1hi, a1lo);
230 two_product_2presplit(
231 a[1], a1hi, a1lo, b[0], bhi, blo, _j, _0
232 );
233 two_sum(_i, _0, _k, _1);
234 fast_two_sum(_j, _k, _l, _2);
235 split(b[1], bhi, blo);
236 two_product_2presplit(
237 a[0], a0hi, a0lo, b[1], bhi, blo, _i, _0
238 );
239 two_sum(_1, _0, _k, x[1]);
240 two_sum(_2, _k, _j, _1);
241 two_sum(_l, _j, _m, _2);
242 two_product_2presplit(
243 a[1], a1hi, a1lo, b[1], bhi, blo, _j, _0
244 );
245 two_sum(_i, _0, _n, _0);
246 two_sum(_1, _0, _i, x[2]);
247 two_sum(_2, _i, _k, _1);
248 two_sum(_m, _k, _l, _2);
249 two_sum(_j, _n, _k, _0);
250 two_sum(_1, _0, _j, x[3]);
251 two_sum(_2, _j, _i, _1);
252 two_sum(_l, _i, _m, _2);
253 two_sum(_1, _k, _i, x[4]);
254 two_sum(_2, _i, _k, x[5]);
255 two_sum(_m, _k, x[7], x[6]);
256 #endif
257 69593685 }
258
259 // [Shewchuk 97]
260 // (https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf)
261 // Section 2.8: other operations
262 // Compression
263 // Note: when converting the algorithms in Shewchuk's article
264 // into code, indices in the article go from 1 to m, and in the
265 // code they go from 0 to m-1 !!!
266 // /!\ there is a bug in the original article,
267 // line 14 of the algorithm should be h_top <= q (small q and not capital Q)
268
269 /**
270 * \brief Compresses an expansion
271 * \details Modifies in-place an expansion in such a way that it
272 * is shorter. The represented value is not modified.
273 * \param[in,out] e a reference to the expansion to be compressed
274 */
275 2290108 void compress_expansion(expansion& e) {
276 2290108 expansion& h = e;
277
278 2290108 index_t m = e.length();
279 double Qnew,q;
280
281 2290108 index_t bottom = m-1;
282
1/2
✓ Branch 1 taken 2290108 times.
✗ Branch 2 not taken.
 2290108 double Q = e[bottom];
283
284
2/2
✓ Branch 0 taken 3841700 times.
✓ Branch 1 taken 2290108 times.
 6131808 for(int i=int(m)-2; i>=0; --i) {
285
1/2
✓ Branch 1 taken 3841700 times.
✗ Branch 2 not taken.
 3841700 fast_two_sum(Q, e[index_t(i)], Qnew, q);
286 3841700 Q = Qnew;
287
2/2
✓ Branch 0 taken 2428276 times.
✓ Branch 1 taken 1413424 times.
 3841700 if(q != 0.0) {
288
1/2
✓ Branch 1 taken 2428276 times.
✗ Branch 2 not taken.
 2428276 h[bottom] = Q;
289 2428276 --bottom;
290 2428276 Q = q;
291 }
292 }
293
1/2
✓ Branch 1 taken 2290108 times.
✗ Branch 2 not taken.
 2290108 h[bottom] = Q;
294
295 2290108 index_t top = 0;
296
2/2
✓ Branch 0 taken 2428276 times.
✓ Branch 1 taken 2290108 times.
 4718384 for(index_t i=bottom+1; i<m; ++i) {
297
1/2
✓ Branch 1 taken 2428276 times.
✗ Branch 2 not taken.
 2428276 fast_two_sum(h[i],Q,Qnew,q);
298 2428276 Q = Qnew;
299
1/2
✓ Branch 0 taken 2428276 times.
✗ Branch 1 not taken.
 2428276 if(q != 0) {
300
1/2
✓ Branch 1 taken 2428276 times.
✗ Branch 2 not taken.
 2428276 h[top] = q;
301 2428276 ++top;
302 }
303 }
304
1/2
✓ Branch 1 taken 2290108 times.
✗ Branch 2 not taken.
 2290108 h[top] = Q;
305
1/2
✓ Branch 1 taken 2290108 times.
✗ Branch 2 not taken.
 2290108 h.set_length(top+1);
306 2290108 }
307 }
308
309 namespace GEO {
310
311 ✗ void grow_expansion_zeroelim(
312 const expansion& e, double b, expansion& h
313 ) {
314 double Q, hh;
315 double Qnew;
316 index_t eindex, hindex;
317 ✗ index_t elen = e.length();
318
319 ✗ hindex = 0;
320 ✗ Q = b;
321 ✗ for(eindex = 0; eindex < elen; eindex++) {
322 ✗ double enow = e[eindex];
323 ✗ two_sum(Q, enow, Qnew, hh);
324 ✗ Q = Qnew;
325 ✗ if(hh != 0.0) {
326 ✗ h[hindex++] = hh;
327 }
328 }
329 ✗ if((Q != 0.0) || (hindex == 0)) {
330 ✗ h[hindex++] = Q;
331 }
332 ✗ h.set_length(hindex);
333 ✗ }
334
335 45136452 void scale_expansion_zeroelim(
336 const expansion& e, double b, expansion& h
337 ) {
338 double Q, sum;
339 double hh;
340 double product1;
341 double product0;
342 index_t eindex, hindex;
343
344 // If the target processor supports the FMA (Fused Multiply Add)
345 // instruction, then the product of two doubles into a length-2
346 // expansion can be implemented as follows. Thanks to Marc Glisse
347 // for the information.
348 // Note: under gcc, automatic generations of fma() for a*b+c needs
349 // to be deactivated, using -ffp-contract=off, else it may break
350 // other functions such as fast_expansion_sum_zeroelim().
351 #ifndef FP_FAST_FMA
352 double bhi, blo;
353 #endif
354 45136452 index_t elen = e.length();
355
356 // Sanity check: e and h cannot be the same.
357
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 45136452 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 45136452 geo_debug_assert(&e != &h);
358
359 #ifdef FP_FAST_FMA
360
1/2
✓ Branch 1 taken 45136452 times.
✗ Branch 2 not taken.
 45136452 two_product(e[0], b, Q, hh);
361 #else
362 split(b, bhi, blo);
363 two_product_presplit(e[0], b, bhi, blo, Q, hh);
364 #endif
365
366 45136452 hindex = 0;
367
2/2
✓ Branch 0 taken 18624459 times.
✓ Branch 1 taken 26511993 times.
 45136452 if(hh != 0) {
368
1/2
✓ Branch 1 taken 18624459 times.
✗ Branch 2 not taken.
 18624459 h[hindex++] = hh;
369 }
370
2/2
✓ Branch 0 taken 113077063 times.
✓ Branch 1 taken 45136452 times.
 158213515 for(eindex = 1; eindex < elen; eindex++) {
371
1/2
✓ Branch 1 taken 113077063 times.
✗ Branch 2 not taken.
 113077063 double enow = e[eindex];
372 #ifdef FP_FAST_FMA
373 113077063 two_product(enow, b, product1, product0);
374 #else
375 two_product_presplit(enow, b, bhi, blo, product1, product0);
376 #endif
377 113077063 two_sum(Q, product0, sum, hh);
378
2/2
✓ Branch 0 taken 13320085 times.
✓ Branch 1 taken 99756978 times.
 113077063 if(hh != 0) {
379
1/2
✓ Branch 1 taken 13320085 times.
✗ Branch 2 not taken.
 13320085 h[hindex++] = hh;
380 }
381 113077063 fast_two_sum(product1, sum, Q, hh);
382
2/2
✓ Branch 0 taken 85127172 times.
✓ Branch 1 taken 27949891 times.
 113077063 if(hh != 0) {
383
1/2
✓ Branch 1 taken 85127172 times.
✗ Branch 2 not taken.
 85127172 h[hindex++] = hh;
384 }
385 }
386
3/4
✓ Branch 0 taken 15634045 times.
✓ Branch 1 taken 29502407 times.
✓ Branch 2 taken 15634045 times.
✗ Branch 3 not taken.
 45136452 if((Q != 0.0) || (hindex == 0)) {
387
1/2
✓ Branch 1 taken 45136452 times.
✗ Branch 2 not taken.
 45136452 h[hindex++] = Q;
388 }
389
1/2
✓ Branch 1 taken 45136452 times.
✗ Branch 2 not taken.
 45136452 h.set_length(hindex);
390 45136452 }
391
392 43362207 void fast_expansion_sum_zeroelim(
393 const expansion& e, const expansion& f, expansion& h
394 ) {
395 double Q;
396 double Qnew;
397 double hh;
398 index_t eindex, findex, hindex;
399 double enow, fnow;
400 43362207 index_t elen = e.length();
401 43362207 index_t flen = f.length();
402
403 // sanity check: h cannot be e or f
404
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 43362207 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 43362207 geo_debug_assert(&h != &e);
405
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 43362207 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 43362207 geo_debug_assert(&h != &f);
406
407
1/2
✓ Branch 1 taken 43362207 times.
✗ Branch 2 not taken.
 43362207 enow = e[0];
408
1/2
✓ Branch 1 taken 43362207 times.
✗ Branch 2 not taken.
 43362207 fnow = f[0];
409 43362207 eindex = findex = 0;
410
2/2
✓ Branch 0 taken 38520289 times.
✓ Branch 1 taken 4841918 times.
 43362207 if((fnow > enow) == (fnow > -enow)) {
411 38520289 Q = enow;
412
1/2
✓ Branch 1 taken 38520289 times.
✗ Branch 2 not taken.
 38520289 enow = e[++eindex];
413 } else {
414 4841918 Q = fnow;
415
1/2
✓ Branch 1 taken 4841918 times.
✗ Branch 2 not taken.
 4841918 fnow = f[++findex];
416 }
417 43362207 hindex = 0;
418
4/4
✓ Branch 0 taken 31107497 times.
✓ Branch 1 taken 12254710 times.
✓ Branch 2 taken 30287494 times.
✓ Branch 3 taken 820003 times.
 43362207 if((eindex < elen) && (findex < flen)) {
419
2/2
✓ Branch 0 taken 23877455 times.
✓ Branch 1 taken 6410039 times.
 30287494 if((fnow > enow) == (fnow > -enow)) {
420 23877455 fast_two_sum(enow, Q, Qnew, hh);
421
1/2
✓ Branch 1 taken 23877455 times.
✗ Branch 2 not taken.
 23877455 enow = e[++eindex];
422 } else {
423 6410039 fast_two_sum(fnow, Q, Qnew, hh);
424
1/2
✓ Branch 1 taken 6410039 times.
✗ Branch 2 not taken.
 6410039 fnow = f[++findex];
425 }
426 30287494 Q = Qnew;
427
2/2
✓ Branch 0 taken 4753199 times.
✓ Branch 1 taken 25534295 times.
 30287494 if(hh != 0.0) {
428
1/2
✓ Branch 1 taken 4753199 times.
✗ Branch 2 not taken.
 4753199 h[hindex++] = hh;
429 }
430
4/4
✓ Branch 0 taken 268557682 times.
✓ Branch 1 taken 22511190 times.
✓ Branch 2 taken 260781378 times.
✓ Branch 3 taken 7776304 times.
 291068872 while((eindex < elen) && (findex < flen)) {
431
2/2
✓ Branch 0 taken 144914011 times.
✓ Branch 1 taken 115867367 times.
 260781378 if((fnow > enow) == (fnow > -enow)) {
432 144914011 two_sum(Q, enow, Qnew, hh);
433
1/2
✓ Branch 1 taken 144914011 times.
✗ Branch 2 not taken.
 144914011 enow = e[++eindex];
434 } else {
435 115867367 two_sum(Q, fnow, Qnew, hh);
436
1/2
✓ Branch 1 taken 115867367 times.
✗ Branch 2 not taken.
 115867367 fnow = f[++findex];
437 }
438 260781378 Q = Qnew;
439
2/2
✓ Branch 0 taken 124391699 times.
✓ Branch 1 taken 136389679 times.
 260781378 if(hh != 0.0) {
440
1/2
✓ Branch 1 taken 124391699 times.
✗ Branch 2 not taken.
 124391699 h[hindex++] = hh;
441 }
442 }
443 }
444
2/2
✓ Branch 0 taken 10983109 times.
✓ Branch 1 taken 43362207 times.
 54345316 while(eindex < elen) {
445 10983109 two_sum(Q, enow, Qnew, hh);
446
1/2
✓ Branch 1 taken 10983109 times.
✗ Branch 2 not taken.
 10983109 enow = e[++eindex];
447 10983109 Q = Qnew;
448
2/2
✓ Branch 0 taken 3659691 times.
✓ Branch 1 taken 7323418 times.
 10983109 if(hh != 0.0) {
449
1/2
✓ Branch 1 taken 3659691 times.
✗ Branch 2 not taken.
 3659691 h[hindex++] = hh;
450 }
451 }
452
2/2
✓ Branch 0 taken 63294707 times.
✓ Branch 1 taken 43362207 times.
 106656914 while(findex < flen) {
453 63294707 two_sum(Q, fnow, Qnew, hh);
454
1/2
✓ Branch 1 taken 63294707 times.
✗ Branch 2 not taken.
 63294707 fnow = f[++findex];
455 63294707 Q = Qnew;
456
2/2
✓ Branch 0 taken 23295827 times.
✓ Branch 1 taken 39998880 times.
 63294707 if(hh != 0.0) {
457
1/2
✓ Branch 1 taken 23295827 times.
✗ Branch 2 not taken.
 23295827 h[hindex++] = hh;
458 }
459 }
460
4/4
✓ Branch 0 taken 14602294 times.
✓ Branch 1 taken 28759913 times.
✓ Branch 2 taken 14556073 times.
✓ Branch 3 taken 46221 times.
 43362207 if((Q != 0.0) || (hindex == 0)) {
461
1/2
✓ Branch 1 taken 43315986 times.
✗ Branch 2 not taken.
 43315986 h[hindex++] = Q;
462 }
463
1/2
✓ Branch 1 taken 43362207 times.
✗ Branch 2 not taken.
 43362207 h.set_length(hindex);
464 43362207 }
465
466 39119417 void fast_expansion_diff_zeroelim(
467 const expansion& e, const expansion& f, expansion& h
468 ) {
469 double Q;
470 double Qnew;
471 double hh;
472 index_t eindex, findex, hindex;
473 double enow, fnow;
474 39119417 index_t elen = e.length();
475 39119417 index_t flen = f.length();
476
477 // sanity check: h cannot be e or f
478
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 39119417 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 39119417 geo_debug_assert(&h != &e);
479
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 39119417 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 39119417 geo_debug_assert(&h != &f);
480
481
1/2
✓ Branch 1 taken 39119417 times.
✗ Branch 2 not taken.
 39119417 enow = e[0];
482
1/2
✓ Branch 1 taken 39119417 times.
✗ Branch 2 not taken.
 39119417 fnow = -f[0];
483 39119417 eindex = findex = 0;
484
2/2
✓ Branch 0 taken 34238550 times.
✓ Branch 1 taken 4880867 times.
 39119417 if((fnow > enow) == (fnow > -enow)) {
485 34238550 Q = enow;
486
1/2
✓ Branch 1 taken 34238550 times.
✗ Branch 2 not taken.
 34238550 enow = e[++eindex];
487 } else {
488 4880867 Q = fnow;
489
1/2
✓ Branch 1 taken 4880867 times.
✗ Branch 2 not taken.
 4880867 fnow = -f[++findex];
490 }
491 39119417 hindex = 0;
492
4/4
✓ Branch 0 taken 37123581 times.
✓ Branch 1 taken 1995836 times.
✓ Branch 2 taken 35515546 times.
✓ Branch 3 taken 1608035 times.
 39119417 if((eindex < elen) && (findex < flen)) {
493
2/2
✓ Branch 0 taken 32071348 times.
✓ Branch 1 taken 3444198 times.
 35515546 if((fnow > enow) == (fnow > -enow)) {
494 32071348 fast_two_sum(enow, Q, Qnew, hh);
495
1/2
✓ Branch 1 taken 32071348 times.
✗ Branch 2 not taken.
 32071348 enow = e[++eindex];
496 } else {
497 3444198 fast_two_sum(fnow, Q, Qnew, hh);
498
1/2
✓ Branch 1 taken 3444198 times.
✗ Branch 2 not taken.
 3444198 fnow = -f[++findex];
499 }
500 35515546 Q = Qnew;
501
2/2
✓ Branch 0 taken 556154 times.
✓ Branch 1 taken 34959392 times.
 35515546 if(hh != 0.0) {
502
1/2
✓ Branch 1 taken 556154 times.
✗ Branch 2 not taken.
 556154 h[hindex++] = hh;
503 }
504
4/4
✓ Branch 0 taken 328061914 times.
✓ Branch 1 taken 24050662 times.
✓ Branch 2 taken 316597030 times.
✓ Branch 3 taken 11464884 times.
 352112576 while((eindex < elen) && (findex < flen)) {
505
2/2
✓ Branch 0 taken 188841784 times.
✓ Branch 1 taken 127755246 times.
 316597030 if((fnow > enow) == (fnow > -enow)) {
506 188841784 two_sum(Q, enow, Qnew, hh);
507
1/2
✓ Branch 1 taken 188841784 times.
✗ Branch 2 not taken.
 188841784 enow = e[++eindex];
508 } else {
509 127755246 two_sum(Q, fnow, Qnew, hh);
510
1/2
✓ Branch 1 taken 127755246 times.
✗ Branch 2 not taken.
 127755246 fnow = -f[++findex];
511 }
512 316597030 Q = Qnew;
513
2/2
✓ Branch 0 taken 33290379 times.
✓ Branch 1 taken 283306651 times.
 316597030 if(hh != 0.0) {
514
1/2
✓ Branch 1 taken 33290379 times.
✗ Branch 2 not taken.
 33290379 h[hindex++] = hh;
515 }
516 }
517 }
518
2/2
✓ Branch 0 taken 15114516 times.
✓ Branch 1 taken 39119417 times.
 54233933 while(eindex < elen) {
519 15114516 two_sum(Q, enow, Qnew, hh);
520
1/2
✓ Branch 1 taken 15114516 times.
✗ Branch 2 not taken.
 15114516 enow = e[++eindex];
521 15114516 Q = Qnew;
522
2/2
✓ Branch 0 taken 4884591 times.
✓ Branch 1 taken 10229925 times.
 15114516 if(hh != 0.0) {
523
1/2
✓ Branch 1 taken 4884591 times.
✗ Branch 2 not taken.
 4884591 h[hindex++] = hh;
524 }
525 }
526
2/2
✓ Branch 0 taken 134118319 times.
✓ Branch 1 taken 39119417 times.
 173237736 while(findex < flen) {
527 134118319 two_sum(Q, fnow, Qnew, hh);
528
1/2
✓ Branch 1 taken 134118319 times.
✗ Branch 2 not taken.
 134118319 fnow = -f[++findex];
529 134118319 Q = Qnew;
530
2/2
✓ Branch 0 taken 4665264 times.
✓ Branch 1 taken 129453055 times.
 134118319 if(hh != 0.0) {
531
1/2
✓ Branch 1 taken 4665264 times.
✗ Branch 2 not taken.
 4665264 h[hindex++] = hh;
532 }
533 }
534
4/4
✓ Branch 0 taken 16356479 times.
✓ Branch 1 taken 22762938 times.
✓ Branch 2 taken 16327216 times.
✓ Branch 3 taken 29263 times.
 39119417 if((Q != 0.0) || (hindex == 0)) {
535
1/2
✓ Branch 1 taken 39090154 times.
✗ Branch 2 not taken.
 39090154 h[hindex++] = Q;
536 }
537
1/2
✓ Branch 1 taken 39119417 times.
✗ Branch 2 not taken.
 39119417 h.set_length(hindex);
538 39119417 }
539
540 }
541
542 /****************************************************************************/
543
544 namespace GEO {
545
546 double expansion_splitter_;
547 double expansion_epsilon_;
548 bool expansion_initialized_ = false;
549
550 254 void expansion::initialize() {
551 // Taken from Jonathan Shewchuk's exactinit.
552 double half;
553 double check, lastcheck;
554 int every_other;
555
556 254 every_other = 1;
557 254 half = 0.5;
558 254 expansion_epsilon_ = 1.0;
559 254 expansion_splitter_ = 1.0;
560 254 check = 1.0;
561 // Repeatedly divide `epsilon' by two until it is too small to add to
562 // one without causing roundoff. (Also check if the sum is equal to
563 // the previous sum, for machines that round up instead of using exact
564 // rounding. Not that this library will work on such machines anyway.
565 do {
566 13462 lastcheck = check;
567 13462 expansion_epsilon_ *= half;
568
2/2
✓ Branch 0 taken 6858 times.
✓ Branch 1 taken 6604 times.
 13462 if(every_other) {
569 6858 expansion_splitter_ *= 2.0;
570 }
571 13462 every_other = !every_other;
572 13462 check = 1.0 + expansion_epsilon_;
573
3/4
✓ Branch 0 taken 13208 times.
✓ Branch 1 taken 254 times.
✓ Branch 2 taken 13208 times.
✗ Branch 3 not taken.
 13462 } while((check != 1.0) && (check != lastcheck));
574 254 expansion_splitter_ += 1.0;
575 254 expansion_initialized_ = true;
576 254 }
577
578 // ====== Initialization from expansion and double ===============
579
580 ✗ expansion& expansion::assign_sum(const expansion& a, double b) {
581 ✗ geo_debug_assert(capacity() >= sum_capacity(a, b));
582 ✗ grow_expansion_zeroelim(a, b, *this);
583 ✗ return *this;
584 }
585
586 ✗ expansion& expansion::assign_diff(const expansion& a, double b) {
587 ✗ geo_debug_assert(capacity() >= diff_capacity(a, b));
588 ✗ grow_expansion_zeroelim(a, -b, *this);
589 ✗ return *this;
590 }
591
592 29090386 expansion& expansion::assign_product(const expansion& a, double b) {
593 // TODO: implement special case where the double argument
594 // is a power of two.
595
1/6
✗ Branch 2 not taken.
✓ Branch 3 taken 29090386 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 29090386 geo_debug_assert(capacity() >= product_capacity(a, b));
596 29090386 scale_expansion_zeroelim(a, b, *this);
597 29090386 return *this;
598 }
599
600 // ============= expansion sum and difference =========================
601
602 43362207 expansion& expansion::assign_sum(
603 const expansion& a, const expansion& b
604 ) {
605
1/6
✗ Branch 2 not taken.
✓ Branch 3 taken 43362207 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 43362207 geo_debug_assert(capacity() >= sum_capacity(a, b));
606 43362207 fast_expansion_sum_zeroelim(a, b, *this);
607 43362207 return *this;
608 }
609
610 7906229 expansion& expansion::assign_sum(
611 const expansion& a, const expansion& b, const expansion& c
612 ) {
613
1/6
✗ Branch 2 not taken.
✓ Branch 3 taken 7906229 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 7906229 geo_debug_assert(capacity() >= sum_capacity(a, b, c));
614
2/6
✓ Branch 6 taken 7906229 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 7906229 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 7906229 expansion& ab = expansion_sum(a, b);
615 7906229 this->assign_sum(ab, c);
616 7906229 return *this;
617 }
618
619 374371 expansion& expansion::assign_sum(
620 const expansion& a, const expansion& b,
621 const expansion& c, const expansion& d
622 ) {
623
1/6
✗ Branch 2 not taken.
✓ Branch 3 taken 374371 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 374371 geo_debug_assert(capacity() >= sum_capacity(a, b, c));
624
2/6
✓ Branch 6 taken 374371 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 374371 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 374371 expansion& ab = expansion_sum(a, b);
625
2/6
✓ Branch 6 taken 374371 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 374371 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 374371 expansion& cd = expansion_sum(c, d);
626 374371 this->assign_sum(ab, cd);
627 374371 return *this;
628 }
629
630 39119417 expansion& expansion::assign_diff(const expansion& a, const expansion& b) {
631
1/6
✗ Branch 2 not taken.
✓ Branch 3 taken 39119417 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 39119417 geo_debug_assert(capacity() >= diff_capacity(a, b));
632 39119417 fast_expansion_diff_zeroelim(a, b, *this);
633 39119417 return *this;
634 }
635
636 // ============= expansion product ==================================
637
638 // Recursive helper function for product implementation
639 66871 expansion& expansion::assign_sub_product(
640 const double* a, index_t a_length, const expansion& b
641 ) {
642
1/6
✗ Branch 3 not taken.
✓ Branch 4 taken 66871 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 66871 geo_debug_assert(
643 capacity() >= sub_product_capacity(a_length, b.length())
644 );
645
2/2
✓ Branch 0 taken 34317 times.
✓ Branch 1 taken 32554 times.
 66871 if(a_length == 1) {
646 34317 scale_expansion_zeroelim(b, a[0], *this);
647 } else {
648 // "Distillation" (see Shewchuk's paper) is computed recursively,
649 // by splitting the list of expansions to sum into two halves.
650
651 32554 const double* a1 = a;
652 32554 index_t a1_length = a_length / 2;
653 32554 const double* a2 = a1 + a1_length;
654 32554 index_t a2_length = a_length - a1_length;
655
656 // Allocate both halves on the stack or on the heap if too large
657 // (some platformes, e.g. MacOSX, have a small stack)
658
659 32554 index_t a1b_capa = sub_product_capacity(a1_length, b.length());
660 32554 index_t a2b_capa = sub_product_capacity(a2_length, b.length());
661
662 32554 bool a1b_on_heap = (a1b_capa > MAX_CAPACITY_ON_STACK);
663 32554 bool a2b_on_heap = (a2b_capa > MAX_CAPACITY_ON_STACK);
664
665 32554 expansion* a1b = a1b_on_heap ?
666 18 new_expansion_on_heap(a1b_capa) :
667
5/6
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 32536 times.
✓ Branch 5 taken 32536 times.
✓ Branch 6 taken 18 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 32536 times.
 32572 new_expansion_on_stack(a1b_capa);
668
669 32554 a1b->assign_sub_product(a1, a1_length, b);
670
671 32554 expansion* a2b = a2b_on_heap ?
672 23 new_expansion_on_heap(a2b_capa) :
673
5/6
✓ Branch 0 taken 23 times.
✓ Branch 1 taken 32531 times.
✓ Branch 5 taken 32531 times.
✓ Branch 6 taken 23 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 32531 times.
 32577 new_expansion_on_stack(a2b_capa);
674
675 32554 a2b->assign_sub_product(a2, a2_length, b);
676
677 32554 this->assign_sum(*a1b, *a2b);
678
679
2/2
✓ Branch 0 taken 18 times.
✓ Branch 1 taken 32536 times.
 32554 if(a1b_on_heap) {
680 18 delete_expansion_on_heap(a1b);
681 }
682
683
2/2
✓ Branch 0 taken 23 times.
✓ Branch 1 taken 32531 times.
 32554 if(a2b_on_heap) {
684 23 delete_expansion_on_heap(a2b);
685 }
686 }
687 66871 return *this;
688 }
689
690 94953111 expansion& expansion::assign_product(
691 const expansion& a, const expansion& b
692 ) {
693
1/6
✗ Branch 2 not taken.
✓ Branch 3 taken 94953111 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 94953111 geo_debug_assert(capacity() >= product_capacity(a, b));
694
3/6
✓ Branch 1 taken 94953111 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 94953111 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 94953111 times.
 94953111 if(a.length() == 0 || b.length() == 0) {
695 ✗ x_[0] = 0.0;
696 ✗ set_length(0);
697
6/6
✓ Branch 1 taken 6552034 times.
✓ Branch 2 taken 88401077 times.
✓ Branch 4 taken 4568993 times.
✓ Branch 5 taken 1983041 times.
✓ Branch 6 taken 4568993 times.
✓ Branch 7 taken 90384118 times.
 94953111 } else if(a.length() == 1 && b.length() == 1) {
698 4568993 two_product(a[0], b[0], x_[1], x_[0]);
699 4568993 set_length(2);
700
2/2
✓ Branch 1 taken 1983041 times.
✓ Branch 2 taken 88401077 times.
 90384118 } else if(a.length() == 1) {
701 1983041 scale_expansion_zeroelim(b, a[0], *this);
702
2/2
✓ Branch 1 taken 14028708 times.
✓ Branch 2 taken 74372369 times.
 88401077 } else if(b.length() == 1) {
703 14028708 scale_expansion_zeroelim(a, b[0], *this);
704
6/6
✓ Branch 1 taken 67517639 times.
✓ Branch 2 taken 6854730 times.
✓ Branch 4 taken 63437358 times.
✓ Branch 5 taken 4080281 times.
✓ Branch 6 taken 63437358 times.
✓ Branch 7 taken 10935011 times.
 74372369 } else if(a.length() == 2 && b.length() == 2) {
705 63437358 two_two_product(a.data(), b.data(), x_);
706 63437358 set_length(8);
707 } else {
708
709
710 10935011 const expansion* pa = &a;
711 10935011 const expansion* pb = &b;
712
713
2/2
✓ Branch 2 taken 3264956 times.
✓ Branch 3 taken 7670055 times.
 10935011 if(pa->length() > pb->length()) {
714 3264956 std::swap(pa, pb);
715 }
716
717 // [Shewchuk 97]
718 // (https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf)
719 // Section 2.8: other operations
720 // Distillation: sum of k values.
721 // Worst case: 1/2*k*(k-1)
722 // But O(k log(k)) if the "summing tree" is well balanced
723 // and using fast_expansion_sum().
724 // Recommended way of computing a product:
725 // compute a1*b, a2*b ... ak*b using scale_expansion_zeroelim()
726 // sum them using a well-balanced tree
727 // However, there is an extra cost for the recursion (and more
728 // importantly, for allocating the intermediary sums, especially
729 // when they do not fit on the stack). So when there are less than
730 // 16 values to add, we simply accumulate them.
731
732 10935011 bool use_balanced_distillation = (pa->length() >= 16);
733
734
2/2
✓ Branch 0 taken 1763 times.
✓ Branch 1 taken 10933248 times.
 10935011 if(use_balanced_distillation) {
735 // assign_sub_product() is a recursive function that
736 // creates a balanced distillation tree on the stack.
737
1/2
✓ Branch 3 taken 1763 times.
✗ Branch 4 not taken.
 1763 assign_sub_product(pa->data(), pa->length(),*pb);
738 } else {
739 // trivial implementation: compute all the products
740 // P = ak*b and accumulate them into S
741
742
1/2
✓ Branch 1 taken 10933248 times.
✗ Branch 2 not taken.
 10933248 index_t P_capa = product_capacity(*pb, 3.0); // 3.0, or any
743 // number that is
744 // not a power of 2
745
746 10933248 index_t S_capa = capacity(); // same capacity as this,
747 // enough to store sum.
748
749 10933248 bool P_on_heap = (P_capa > MAX_CAPACITY_ON_STACK);
750 10933248 bool S_on_heap = (S_capa > MAX_CAPACITY_ON_STACK);
751
752 10933248 expansion* P = P_on_heap ?
753 ✗ new_expansion_on_heap(P_capa) :
754
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 10933248 times.
✓ Branch 5 taken 10933248 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 10933248 times.
 10933248 new_expansion_on_stack(P_capa);
755
756 10933248 expansion* S = S_on_heap ?
757 1 new_expansion_on_heap(S_capa) :
758
5/6
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 10933247 times.
✓ Branch 5 taken 10933247 times.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 10933247 times.
 10933249 new_expansion_on_stack(S_capa);
759
760 10933248 expansion* S1 = S;
761 10933248 expansion* S2 = this;
762
763
2/2
✓ Branch 1 taken 6784860 times.
✓ Branch 2 taken 4148388 times.
 10933248 if((pa->length()%2) == 0) {
764 6784860 std::swap(S1,S2);
765 }
766
767
2/2
✓ Branch 1 taken 27753466 times.
✓ Branch 2 taken 10933248 times.
 38686714 for(index_t i=0; i<pa->length(); ++i) {
768
2/2
✓ Branch 0 taken 10933248 times.
✓ Branch 1 taken 16820218 times.
 27753466 if(i == 0) {
769
2/4
✓ Branch 1 taken 10933248 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10933248 times.
✗ Branch 5 not taken.
 10933248 S2->assign_product(*pb, (*pa)[i]);
770 } else {
771
2/4
✓ Branch 1 taken 16820218 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16820218 times.
✗ Branch 5 not taken.
 16820218 P->assign_product(*pb, (*pa)[i]);
772
1/2
✓ Branch 1 taken 16820218 times.
✗ Branch 2 not taken.
 16820218 S2->assign_sum(*S1,*P);
773 }
774 27753466 std::swap(S1,S2);
775 }
776
777
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 10933248 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 10933248 geo_assert(S1 == this);
778
779
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 10933247 times.
 10933248 if(S_on_heap) {
780 1 delete_expansion_on_heap(S);
781 }
782
783
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10933248 times.
 10933248 if(P_on_heap) {
784 ✗ delete_expansion_on_heap(P);
785 }
786 }
787 }
788 94953111 return *this;
789 }
790
791 ✗ expansion& expansion::assign_product(
792 const expansion& a, const expansion& b, const expansion& c
793 ) {
794 ✗ const expansion& bc = expansion_product(b, c);
795 ✗ this->assign_product(a, bc);
796 ✗ return *this;
797 }
798
799 ✗ expansion& expansion::assign_square(const expansion& a) {
800 ✗ geo_debug_assert(capacity() >= square_capacity(a));
801 ✗ if(a.length() == 1) {
802 ✗ square(a[0], x_[1], x_[0]);
803 ✗ set_length(2);
804 ✗ } else if(a.length() == 2) {
805 ✗ two_square(a[1], a[0], x_);
806 ✗ set_length(6);
807 } else {
808 ✗ this->assign_product(a, a);
809 }
810 ✗ return *this;
811 }
812
813 // ============= determinants ==========================================
814
815 32668732 expansion& expansion::assign_det2x2(
816 const expansion& a11, const expansion& a12,
817 const expansion& a21, const expansion& a22
818 ) {
819
2/6
✓ Branch 6 taken 32668732 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 32668732 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 32668732 const expansion& a11a22 = expansion_product(a11, a22);
820
2/6
✓ Branch 6 taken 32668732 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 32668732 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 32668732 const expansion& a12a21 = expansion_product(a12, a21);
821 32668732 return this->assign_diff(a11a22, a12a21);
822 }
823
824 5894622 expansion& expansion::assign_det3x3(
825 const expansion& a11, const expansion& a12, const expansion& a13,
826 const expansion& a21, const expansion& a22, const expansion& a23,
827 const expansion& a31, const expansion& a32, const expansion& a33
828 ) {
829 // Development w.r.t. first row
830
2/6
✓ Branch 6 taken 5894622 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 5894622 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 5894622 const expansion& c11 = expansion_det2x2(a22, a23, a32, a33);
831
2/6
✓ Branch 6 taken 5894622 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 5894622 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 5894622 const expansion& c12 = expansion_det2x2(a23, a21, a33, a31);
832
2/6
✓ Branch 6 taken 5894622 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 5894622 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 5894622 const expansion& c13 = expansion_det2x2(a21, a22, a31, a32);
833
2/6
✓ Branch 6 taken 5894622 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 5894622 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 5894622 const expansion& a11c11 = expansion_product(a11, c11);
834
2/6
✓ Branch 6 taken 5894622 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 5894622 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 5894622 const expansion& a12c12 = expansion_product(a12, c12);
835
2/6
✓ Branch 6 taken 5894622 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 5894622 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 5894622 const expansion& a13c13 = expansion_product(a13, c13);
836 5894622 return this->assign_sum(a11c11, a12c12, a13c13);
837 }
838
839 ✗ expansion& expansion::assign_det_111_2x3(
840 const expansion& a21, const expansion& a22, const expansion& a23,
841 const expansion& a31, const expansion& a32, const expansion& a33
842 ) {
843 ✗ const expansion& c11 = expansion_det2x2(a22, a23, a32, a33);
844 ✗ const expansion& c12 = expansion_det2x2(a23, a21, a33, a31);
845 ✗ const expansion& c13 = expansion_det2x2(a21, a22, a31, a32);
846 ✗ return this->assign_sum(c11, c12, c13);
847 }
848
849 // ============= geometric operations ==================================
850
851 10087862 expansion& expansion::assign_sq_dist(
852 const double* p1, const double* p2, coord_index_t dim
853 ) {
854
1/6
✗ Branch 2 not taken.
✓ Branch 3 taken 10087862 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 10087862 geo_debug_assert(capacity() >= sq_dist_capacity(dim));
855
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 10087862 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 10087862 geo_debug_assert(dim > 0);
856
2/2
✓ Branch 0 taken 6052734 times.
✓ Branch 1 taken 4035128 times.
 10087862 if(dim == 1) {
857 double d0, d1;
858 6052734 two_diff(p1[0], p2[0], d1, d0);
859 6052734 two_square(d1, d0, x_);
860
1/2
✓ Branch 1 taken 6052734 times.
✗ Branch 2 not taken.
 6052734 set_length(6);
861 } else {
862 // "Distillation" (see Shewchuk's paper) is computed recursively,
863 // by splitting the list of expansions to sum into two halves.
864 4035128 coord_index_t dim1 = dim / 2;
865 4035128 coord_index_t dim2 = coord_index_t(dim - dim1);
866 4035128 const double* p1_2 = p1 + dim1;
867 4035128 const double* p2_2 = p2 + dim1;
868
2/6
✓ Branch 6 taken 4035128 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 4035128 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 4035128 expansion& d1 = expansion_sq_dist(p1, p2, dim1);
869
2/6
✓ Branch 6 taken 4035128 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 4035128 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 4035128 expansion& d2 = expansion_sq_dist(p1_2, p2_2, dim2);
870 4035128 this->assign_sum(d1, d2);
871 }
872 10087862 return *this;
873 }
874
875 10260461 expansion& expansion::assign_dot_at(
876 const double* p1, const double* p2, const double* p0,
877 coord_index_t dim
878 ) {
879
1/6
✗ Branch 2 not taken.
✓ Branch 3 taken 10260461 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 10260461 geo_debug_assert(capacity() >= dot_at_capacity(dim));
880
2/2
✓ Branch 0 taken 6156327 times.
✓ Branch 1 taken 4104134 times.
 10260461 if(dim == 1) {
881
882 double v[2];
883 6156327 two_diff(p1[0], p0[0], v[1], v[0]);
884 double w[2];
885 6156327 two_diff(p2[0], p0[0], w[1], w[0]);
886 6156327 two_two_product(v, w, x_);
887
1/2
✓ Branch 1 taken 6156327 times.
✗ Branch 2 not taken.
 6156327 set_length(8);
888 } else {
889 // "Distillation" (see Shewchuk's paper) is computed recursively,
890 // by splitting the list of expansions to sum into two halves.
891 4104134 coord_index_t dim1 = dim / 2;
892 4104134 coord_index_t dim2 = coord_index_t(dim - dim1);
893 4104134 const double* p1_2 = p1 + dim1;
894 4104134 const double* p2_2 = p2 + dim1;
895 4104134 const double* p0_2 = p0 + dim1;
896
2/6
✓ Branch 6 taken 4104134 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 4104134 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 4104134 expansion& d1 = expansion_dot_at(p1, p2, p0, dim1);
897
2/6
✓ Branch 6 taken 4104134 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 4104134 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 4104134 expansion& d2 = expansion_dot_at(p1_2, p2_2, p0_2, dim2);
898 4104134 this->assign_sum(d1, d2);
899 }
900 10260461 return *this;
901 }
902
903 ✗ expansion& expansion::assign_length2(
904 const expansion& x, const expansion& y, const expansion& z
905 ) {
906 ✗ const expansion& x2 = expansion_square(x);
907 ✗ const expansion& y2 = expansion_square(y);
908 ✗ const expansion& z2 = expansion_square(z);
909 ✗ this->assign_sum(x2,y2,z2);
910 ✗ return *this;
911 }
912
913 /************************************************************************/
914
915 1575968 bool expansion::is_same_as(const expansion& rhs) const {
916
2/2
✓ Branch 2 taken 406222 times.
✓ Branch 3 taken 1169746 times.
 1575968 if(length() != rhs.length()) {
917 406222 return false;
918 }
919
2/2
✓ Branch 1 taken 1405988 times.
✓ Branch 2 taken 747164 times.
 2153152 for(index_t i=0; i<length(); ++i) {
920
2/2
✓ Branch 0 taken 422582 times.
✓ Branch 1 taken 983406 times.
 1405988 if(x_[i] != rhs.x_[i]) {
921 422582 return false;
922 }
923 }
924 747164 return true;
925 }
926
927 ✗ bool expansion::is_same_as(double rhs) const {
928 ✗ if(length() != 1) {
929 ✗ return false;
930 }
931 ✗ return (x_[0] == rhs);
932 }
933
934 4214299 Sign expansion::compare(const expansion& rhs) const {
935 // Fast path: different signs or both zero
936 4214299 Sign s1 = sign();
937 4214299 Sign s2 = rhs.sign();
938
4/4
✓ Branch 0 taken 1514137 times.
✓ Branch 1 taken 2700162 times.
✓ Branch 2 taken 492609 times.
✓ Branch 3 taken 1021528 times.
 4214299 if(s1 == ZERO && s2 == ZERO) {
939 492609 return ZERO;
940 }
941
2/2
✓ Branch 0 taken 2145722 times.
✓ Branch 1 taken 1575968 times.
 3721690 if(s1 != s2) {
942
2/2
✓ Branch 0 taken 694983 times.
✓ Branch 1 taken 1450739 times.
 2145722 return (int(s1) > int(s2) ? POSITIVE : NEGATIVE);
943 }
944
945 // Fast path: same internal representation
946
2/2
✓ Branch 1 taken 747164 times.
✓ Branch 2 taken 828804 times.
 1575968 if(is_same_as(rhs)) {
947 747164 return ZERO;
948 }
949
950 // Compute difference and return sign of difference
951 828804 index_t capa = diff_capacity(*this, rhs);
952
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 828804 times.
 828804 if(capa > MAX_CAPACITY_ON_STACK) {
953 ✗ expansion* d = new_expansion_on_heap(capa);
954 ✗ d->assign_diff(*this, rhs);
955 ✗ Sign result = d->sign();
956 ✗ delete_expansion_on_heap(d);
957 ✗ return result;
958 }
959
2/6
✓ Branch 6 taken 828804 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 828804 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 828804 const expansion& d = expansion_diff(*this, rhs);
960 828804 return d.sign();
961 }
962
963 ✗ Sign expansion::compare(double rhs) const {
964 // Fast path: different signs or both zero
965 ✗ Sign s1 = sign();
966 ✗ Sign s2 = geo_sgn(rhs);
967 ✗ if(s1 == ZERO && s2 == ZERO) {
968 ✗ return ZERO;
969 }
970 ✗ if(s1 != s2) {
971 ✗ return (int(s1) > int(s2) ? POSITIVE : NEGATIVE);
972 }
973
974 // Fast path: same internal representation
975 ✗ if(is_same_as(rhs)) {
976 ✗ return ZERO;
977 }
978
979 // Compute difference and return sign of difference
980 ✗ index_t capa = diff_capacity(*this, rhs);
981 ✗ if(capa > MAX_CAPACITY_ON_STACK) {
982 ✗ expansion* d = new_expansion_on_heap(capa);
983 ✗ d->assign_diff(*this, rhs);
984 ✗ Sign result = d->sign();
985 ✗ delete_expansion_on_heap(d);
986 ✗ return result;
987 }
988 ✗ const expansion& d = expansion_diff(*this, rhs);
989 ✗ return d.sign();
990 }
991
992
993 /************************************************************************/
994
995 ✗ void expansion::show_all_stats() {
996 #ifdef PCK_STATS
997 // Place holder: if we compute statistics for expansions,
998 // the code here will be called if sys:stats is specified
999 // on command line.
1000 #endif
1001 ✗ }
1002
1003 /************************************************************************/
1004
1005 ✗ Sign sign_of_expansion_determinant(
1006 const expansion& a00,const expansion& a01,
1007 const expansion& a10,const expansion& a11
1008 ) {
1009 ✗ const expansion& result = expansion_det2x2(a00, a01, a10, a11);
1010 ✗ return result.sign();
1011 }
1012
1013 ✗ Sign sign_of_expansion_determinant(
1014 const expansion& a00,const expansion& a01,const expansion& a02,
1015 const expansion& a10,const expansion& a11,const expansion& a12,
1016 const expansion& a20,const expansion& a21,const expansion& a22
1017 ) {
1018 // First compute the det2x2
1019 const expansion& m01 =
1020 ✗ expansion_det2x2(a00, a10, a01, a11);
1021 const expansion& m02 =
1022 ✗ expansion_det2x2(a00, a20, a01, a21);
1023 const expansion& m12 =
1024 ✗ expansion_det2x2(a10, a20, a11, a21);
1025
1026 // Now compute the minors of rank 3
1027 ✗ const expansion& z1 = expansion_product(m01,a22);
1028 ✗ const expansion& z2 = expansion_product(m02,a12).negate();
1029 ✗ const expansion& z3 = expansion_product(m12,a02);
1030
1031 ✗ const expansion& result = expansion_sum3(z1,z2,z3);
1032 ✗ return result.sign();
1033 }
1034
1035 ✗ Sign sign_of_expansion_determinant(
1036 const expansion& a00,const expansion& a01,
1037 const expansion& a02,const expansion& a03,
1038 const expansion& a10,const expansion& a11,
1039 const expansion& a12,const expansion& a13,
1040 const expansion& a20,const expansion& a21,
1041 const expansion& a22,const expansion& a23,
1042 const expansion& a30,const expansion& a31,
1043 const expansion& a32,const expansion& a33
1044 ) {
1045
1046 // First compute the det2x2
1047 const expansion& m01 =
1048 ✗ expansion_det2x2(a10,a00,a11,a01);
1049 const expansion& m02 =
1050 ✗ expansion_det2x2(a20,a00,a21,a01);
1051 const expansion& m03 =
1052 ✗ expansion_det2x2(a30,a00,a31,a01);
1053 const expansion& m12 =
1054 ✗ expansion_det2x2(a20,a10,a21,a11);
1055 const expansion& m13 =
1056 ✗ expansion_det2x2(a30,a10,a31,a11);
1057 const expansion& m23 =
1058 ✗ expansion_det2x2(a30,a20,a31,a21);
1059
1060 // Now compute the minors of rank 3
1061 ✗ const expansion& m012_1 = expansion_product(m12,a02);
1062 ✗ expansion& m012_2 = expansion_product(m02,a12); m012_2.negate();
1063 ✗ const expansion& m012_3 = expansion_product(m01,a22);
1064 ✗ const expansion& m012 = expansion_sum3(m012_1, m012_2, m012_3);
1065
1066 ✗ const expansion& m013_1 = expansion_product(m13,a02);
1067 ✗ expansion& m013_2 = expansion_product(m03,a12); m013_2.negate();
1068
1069 ✗ const expansion& m013_3 = expansion_product(m01,a32);
1070 ✗ const expansion& m013 = expansion_sum3(m013_1, m013_2, m013_3);
1071
1072 ✗ const expansion& m023_1 = expansion_product(m23,a02);
1073 ✗ expansion& m023_2 = expansion_product(m03,a22); m023_2.negate();
1074 ✗ const expansion& m023_3 = expansion_product(m02,a32);
1075 ✗ const expansion& m023 = expansion_sum3(m023_1, m023_2, m023_3);
1076
1077 ✗ const expansion& m123_1 = expansion_product(m23,a12);
1078 ✗ expansion& m123_2 = expansion_product(m13,a22); m123_2.negate();
1079 ✗ const expansion& m123_3 = expansion_product(m12,a32);
1080 ✗ const expansion& m123 = expansion_sum3(m123_1, m123_2, m123_3);
1081
1082 // Now compute the minors of rank 4
1083 ✗ const expansion& m0123_1 = expansion_product(m123,a03);
1084 ✗ const expansion& m0123_2 = expansion_product(m023,a13);
1085 ✗ const expansion& m0123_3 = expansion_product(m013,a23);
1086 ✗ const expansion& m0123_4 = expansion_product(m012,a33);
1087
1088 ✗ const expansion& z1 = expansion_sum(m0123_1, m0123_3);
1089 ✗ const expansion& z2 = expansion_sum(m0123_2, m0123_4);
1090
1091 ✗ const expansion& result = expansion_diff(z1,z2);
1092 ✗ return result.sign();
1093 }
1094
1095 /************************************************************************/
1096
1097 2290108 void expansion::optimize() {
1098 2290108 compress_expansion(*this);
1099 2290108 }
1100
1101 /************************************************************************/
1102
1103 }
1104