Blender  V3.3
math_boolean.cc
Go to the documentation of this file.
1 /* SPDX-License-Identifier: GPL-2.0-or-later */
2 
7 #include "BLI_hash.hh"
8 #include "BLI_math_boolean.hh"
9 #include "BLI_math_mpq.hh"
10 #include "BLI_math_vec_types.hh"
11 #include "BLI_span.hh"
12 #include "BLI_utildefines.h"
13 
14 namespace blender {
15 
16 #ifdef WITH_GMP
17 int orient2d(const mpq2 &a, const mpq2 &b, const mpq2 &c)
18 {
19  mpq_class detleft = (a[0] - c[0]) * (b[1] - c[1]);
20  mpq_class detright = (a[1] - c[1]) * (b[0] - c[0]);
21  mpq_class det = detleft - detright;
22  return sgn(det);
23 }
24 
25 int incircle(const mpq2 &a, const mpq2 &b, const mpq2 &c, const mpq2 &d)
26 {
27  mpq_class adx = a[0] - d[0];
28  mpq_class bdx = b[0] - d[0];
29  mpq_class cdx = c[0] - d[0];
30  mpq_class ady = a[1] - d[1];
31  mpq_class bdy = b[1] - d[1];
32  mpq_class cdy = c[1] - d[1];
33 
34  mpq_class bdxcdy = bdx * cdy;
35  mpq_class cdxbdy = cdx * bdy;
36  mpq_class alift = adx * adx + ady * ady;
37 
38  mpq_class cdxady = cdx * ady;
39  mpq_class adxcdy = adx * cdy;
40  mpq_class blift = bdx * bdx + bdy * bdy;
41 
42  mpq_class adxbdy = adx * bdy;
43  mpq_class bdxady = bdx * ady;
44  mpq_class clift = cdx * cdx + cdy * cdy;
45 
46  mpq_class det = alift * (bdxcdy - cdxbdy) + blift * (cdxady - adxcdy) +
47  clift * (adxbdy - bdxady);
48  return sgn(det);
49 }
50 
51 int orient3d(const mpq3 &a, const mpq3 &b, const mpq3 &c, const mpq3 &d)
52 {
53  mpq_class adx = a[0] - d[0];
54  mpq_class bdx = b[0] - d[0];
55  mpq_class cdx = c[0] - d[0];
56  mpq_class ady = a[1] - d[1];
57  mpq_class bdy = b[1] - d[1];
58  mpq_class cdy = c[1] - d[1];
59  mpq_class adz = a[2] - d[2];
60  mpq_class bdz = b[2] - d[2];
61  mpq_class cdz = c[2] - d[2];
62 
63  mpq_class bdxcdy = bdx * cdy;
64  mpq_class cdxbdy = cdx * bdy;
65 
66  mpq_class cdxady = cdx * ady;
67  mpq_class adxcdy = adx * cdy;
68 
69  mpq_class adxbdy = adx * bdy;
70  mpq_class bdxady = bdx * ady;
71 
72  mpq_class det = adz * (bdxcdy - cdxbdy) + bdz * (cdxady - adxcdy) + cdz * (adxbdy - bdxady);
73  return sgn(det);
74 }
75 #endif /* WITH_GMP */
76 
83 namespace robust_pred {
84 
85 /* Using Shewchuk's file here, edited to removed unneeded functions,
86  * change REAL to double everywhere, added const to some arguments,
87  * and to export only the following declared non-static functions.
88  *
89  * Since this is C++, an instantiated singleton class is used to make
90  * sure that #exactinit() is called once.
91  * (Because it's undefined when this is called in initialization of all modules,
92  * other modules shouldn't use these functions in initialization.)
93  */
94 
95 void exactinit();
96 double orient2dfast(const double *pa, const double *pb, const double *pc);
97 double orient2d(const double *pa, const double *pb, const double *pc);
98 double orient3dfast(const double *pa, const double *pb, const double *pc, const double *pd);
99 double orient3d(const double *pa, const double *pb, const double *pc, const double *pd);
100 double incirclefast(const double *pa, const double *pb, const double *pc, const double *pd);
101 double incircle(const double *pa, const double *pb, const double *pc, const double *pd);
102 double inspherefast(
103  const double *pa, const double *pb, const double *pc, const double *pd, const double *pe);
104 double insphere(
105  const double *pa, const double *pb, const double *pc, const double *pd, const double *pe);
106 
108  public:
110  {
111  exactinit();
112  }
113 };
114 
116 
117 /* Routines for Arbitrary Precision Floating-point Arithmetic
118  * and Fast Robust Geometric Predicates
119  * (predicates.c)
120  *
121  * May 18, 1996
122  *
123  * Placed in the public domain by
124  * Jonathan Richard Shewchuk
125  * School of Computer Science
126  * Carnegie Mellon University
127  * 5000 Forbes Avenue
128  * Pittsburgh, Pennsylvania 15213-3891
129  * <jrs@cs.cmu.edu>
130  *
131  * This file contains C implementation of algorithms for exact addition
132  * and multiplication of floating-point numbers, and predicates for
133  * robustly performing the orientation and incircle tests used in
134  * computational geometry. The algorithms and underlying theory are
135  * described in Jonathan Richard Shewchuk. "Adaptive Precision Floating-
136  * Point Arithmetic and Fast Robust Geometric Predicates." Technical
137  * Report CMU-CS-96-140, School of Computer Science, Carnegie Mellon
138  * University, Pittsburgh, Pennsylvania, May 1996. (Submitted to
139  * Discrete & Computational Geometry.)
140  *
141  * This file, the paper listed above, and other information are available
142  * from the Web page http://www.cs.cmu.edu/~quake/robust.html .
143  *
144  *
145  * Using this code:
146  *
147  * First, read the short or long version of the paper (from the Web page above).
148  *
149  * Be sure to call #exactinit() once, before calling any of the arithmetic
150  * functions or geometric predicates. Also be sure to turn on the
151  * optimizer when compiling this file.
152  */
153 
154 /* On some machines, the exact arithmetic routines might be defeated by the
155  * use of internal extended precision floating-point registers. Sometimes
156  * this problem can be fixed by defining certain values to be volatile,
157  * thus forcing them to be stored to memory and rounded off. This isn't
158  * a great solution, though, as it slows the arithmetic down.
159  *
160  * To try this out, write "#define INEXACT volatile" below. Normally,
161  * however, INEXACT should be defined to be nothing. ("#define INEXACT".)
162  */
163 
164 #define INEXACT /* Nothing */
165 /* #define INEXACT volatile */
166 
167 /* Which of the following two methods of finding the absolute values is
168  * fastest is compiler-dependent. A few compilers can inline and optimize
169  * the fabs() call; but most will incur the overhead of a function call,
170  * which is disastrously slow. A faster way on IEEE machines might be to
171  * mask the appropriate bit, but that's difficult to do in C.
172  */
173 
174 #define Absolute(a) ((a) >= 0.0 ? (a) : -(a))
175 /* #define Absolute(a) fabs(a) */
176 
177 /* Many of the operations are broken up into two pieces, a main part that
178  * performs an approximate operation, and a "tail" that computes the
179  * round-off error of that operation.
180  *
181  * The operations Fast_Two_Sum(), Fast_Two_Diff(), Two_Sum(), Two_Diff(),
182  * Split(), and Two_Product() are all implemented as described in the
183  * reference. Each of these macros requires certain variables to be
184  * defined in the calling routine. The variables `bvirt', `c', `abig',
185  * `_i', `_j', `_k', `_l', `_m', and `_n' are declared `INEXACT' because
186  * they store the result of an operation that may incur round-off error.
187  * The input parameter `x' (or the highest numbered `x_' parameter) must
188  * also be declared `INEXACT'.
189  */
190 
191 #define Fast_Two_Sum_Tail(a, b, x, y) \
192  bvirt = x - a; \
193  y = b - bvirt
194 
195 #define Fast_Two_Sum(a, b, x, y) \
196  x = (double)(a + b); \
197  Fast_Two_Sum_Tail(a, b, x, y)
198 
199 #define Fast_Two_Diff_Tail(a, b, x, y) \
200  bvirt = a - x; \
201  y = bvirt - b
202 
203 #define Fast_Two_Diff(a, b, x, y) \
204  x = (double)(a - b); \
205  Fast_Two_Diff_Tail(a, b, x, y)
206 
207 #define Two_Sum_Tail(a, b, x, y) \
208  bvirt = (double)(x - a); \
209  avirt = x - bvirt; \
210  bround = b - bvirt; \
211  around = a - avirt; \
212  y = around + bround
213 
214 #define Two_Sum(a, b, x, y) \
215  x = (double)(a + b); \
216  Two_Sum_Tail(a, b, x, y)
217 
218 #define Two_Diff_Tail(a, b, x, y) \
219  bvirt = (double)(a - x); \
220  avirt = x + bvirt; \
221  bround = bvirt - b; \
222  around = a - avirt; \
223  y = around + bround
224 
225 #define Two_Diff(a, b, x, y) \
226  x = (double)(a - b); \
227  Two_Diff_Tail(a, b, x, y)
228 
229 #define Split(a, ahi, alo) \
230  c = (double)(splitter * a); \
231  abig = (double)(c - a); \
232  ahi = c - abig; \
233  alo = a - ahi
234 
235 #define Two_Product_Tail(a, b, x, y) \
236  Split(a, ahi, alo); \
237  Split(b, bhi, blo); \
238  err1 = x - (ahi * bhi); \
239  err2 = err1 - (alo * bhi); \
240  err3 = err2 - (ahi * blo); \
241  y = (alo * blo) - err3
242 
243 #define Two_Product(a, b, x, y) \
244  x = (double)(a * b); \
245  Two_Product_Tail(a, b, x, y)
246 
247 #define Two_Product_Presplit(a, b, bhi, blo, x, y) \
248  x = (double)(a * b); \
249  Split(a, ahi, alo); \
250  err1 = x - (ahi * bhi); \
251  err2 = err1 - (alo * bhi); \
252  err3 = err2 - (ahi * blo); \
253  y = (alo * blo) - err3
254 
255 #define Two_Product_2Presplit(a, ahi, alo, b, bhi, blo, x, y) \
256  x = (double)(a * b); \
257  err1 = x - (ahi * bhi); \
258  err2 = err1 - (alo * bhi); \
259  err3 = err2 - (ahi * blo); \
260  y = (alo * blo) - err3
261 
262 #define Square_Tail(a, x, y) \
263  Split(a, ahi, alo); \
264  err1 = x - (ahi * ahi); \
265  err3 = err1 - ((ahi + ahi) * alo); \
266  y = (alo * alo) - err3
267 
268 #define Square(a, x, y) \
269  x = (double)(a * a); \
270  Square_Tail(a, x, y)
271 
272 #define Two_One_Sum(a1, a0, b, x2, x1, x0) \
273  Two_Sum(a0, b, _i, x0); \
274  Two_Sum(a1, _i, x2, x1)
275 
276 #define Two_One_Diff(a1, a0, b, x2, x1, x0) \
277  Two_Diff(a0, b, _i, x0); \
278  Two_Sum(a1, _i, x2, x1)
279 
280 #define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \
281  Two_One_Sum(a1, a0, b0, _j, _0, x0); \
282  Two_One_Sum(_j, _0, b1, x3, x2, x1)
283 
284 #define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \
285  Two_One_Diff(a1, a0, b0, _j, _0, x0); \
286  Two_One_Diff(_j, _0, b1, x3, x2, x1)
287 
288 #define Four_One_Sum(a3, a2, a1, a0, b, x4, x3, x2, x1, x0) \
289  Two_One_Sum(a1, a0, b, _j, x1, x0); \
290  Two_One_Sum(a3, a2, _j, x4, x3, x2)
291 
292 #define Four_Two_Sum(a3, a2, a1, a0, b1, b0, x5, x4, x3, x2, x1, x0) \
293  Four_One_Sum(a3, a2, a1, a0, b0, _k, _2, _1, _0, x0); \
294  Four_One_Sum(_k, _2, _1, _0, b1, x5, x4, x3, x2, x1)
295 
296 #define Four_Four_Sum(a3, a2, a1, a0, b4, b3, b1, b0, x7, x6, x5, x4, x3, x2, x1, x0) \
297  Four_Two_Sum(a3, a2, a1, a0, b1, b0, _l, _2, _1, _0, x1, x0); \
298  Four_Two_Sum(_l, _2, _1, _0, b4, b3, x7, x6, x5, x4, x3, x2)
299 
300 #define Eight_One_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b, x8, x7, x6, x5, x4, x3, x2, x1, x0) \
301  Four_One_Sum(a3, a2, a1, a0, b, _j, x3, x2, x1, x0); \
302  Four_One_Sum(a7, a6, a5, a4, _j, x8, x7, x6, x5, x4)
303 
304 #define Eight_Two_Sum( \
305  a7, a6, a5, a4, a3, a2, a1, a0, b1, b0, x9, x8, x7, x6, x5, x4, x3, x2, x1, x0) \
306  Eight_One_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b0, _k, _6, _5, _4, _3, _2, _1, _0, x0); \
307  Eight_One_Sum(_k, _6, _5, _4, _3, _2, _1, _0, b1, x9, x8, x7, x6, x5, x4, x3, x2, x1)
308 
309 #define Eight_Four_Sum(a7, \
310  a6, \
311  a5, \
312  a4, \
313  a3, \
314  a2, \
315  a1, \
316  a0, \
317  b4, \
318  b3, \
319  b1, \
320  b0, \
321  x11, \
322  x10, \
323  x9, \
324  x8, \
325  x7, \
326  x6, \
327  x5, \
328  x4, \
329  x3, \
330  x2, \
331  x1, \
332  x0) \
333  Eight_Two_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b1, b0, _l, _6, _5, _4, _3, _2, _1, _0, x1, x0); \
334  Eight_Two_Sum(_l, _6, _5, _4, _3, _2, _1, _0, b4, b3, x11, x10, x9, x8, x7, x6, x5, x4, x3, x2)
335 
336 #define Two_One_Product(a1, a0, b, x3, x2, x1, x0) \
337  Split(b, bhi, blo); \
338  Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \
339  Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
340  Two_Sum(_i, _0, _k, x1); \
341  Fast_Two_Sum(_j, _k, x3, x2)
342 
343 #define Four_One_Product(a3, a2, a1, a0, b, x7, x6, x5, x4, x3, x2, x1, x0) \
344  Split(b, bhi, blo); \
345  Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \
346  Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
347  Two_Sum(_i, _0, _k, x1); \
348  Fast_Two_Sum(_j, _k, _i, x2); \
349  Two_Product_Presplit(a2, b, bhi, blo, _j, _0); \
350  Two_Sum(_i, _0, _k, x3); \
351  Fast_Two_Sum(_j, _k, _i, x4); \
352  Two_Product_Presplit(a3, b, bhi, blo, _j, _0); \
353  Two_Sum(_i, _0, _k, x5); \
354  Fast_Two_Sum(_j, _k, x7, x6)
355 
356 #define Two_Two_Product(a1, a0, b1, b0, x7, x6, x5, x4, x3, x2, x1, x0) \
357  Split(a0, a0hi, a0lo); \
358  Split(b0, bhi, blo); \
359  Two_Product_2Presplit(a0, a0hi, a0lo, b0, bhi, blo, _i, x0); \
360  Split(a1, a1hi, a1lo); \
361  Two_Product_2Presplit(a1, a1hi, a1lo, b0, bhi, blo, _j, _0); \
362  Two_Sum(_i, _0, _k, _1); \
363  Fast_Two_Sum(_j, _k, _l, _2); \
364  Split(b1, bhi, blo); \
365  Two_Product_2Presplit(a0, a0hi, a0lo, b1, bhi, blo, _i, _0); \
366  Two_Sum(_1, _0, _k, x1); \
367  Two_Sum(_2, _k, _j, _1); \
368  Two_Sum(_l, _j, _m, _2); \
369  Two_Product_2Presplit(a1, a1hi, a1lo, b1, bhi, blo, _j, _0); \
370  Two_Sum(_i, _0, _n, _0); \
371  Two_Sum(_1, _0, _i, x2); \
372  Two_Sum(_2, _i, _k, _1); \
373  Two_Sum(_m, _k, _l, _2); \
374  Two_Sum(_j, _n, _k, _0); \
375  Two_Sum(_1, _0, _j, x3); \
376  Two_Sum(_2, _j, _i, _1); \
377  Two_Sum(_l, _i, _m, _2); \
378  Two_Sum(_1, _k, _i, x4); \
379  Two_Sum(_2, _i, _k, x5); \
380  Two_Sum(_m, _k, x7, x6)
381 
382 #define Two_Square(a1, a0, x5, x4, x3, x2, x1, x0) \
383  Square(a0, _j, x0); \
384  _0 = a0 + a0; \
385  Two_Product(a1, _0, _k, _1); \
386  Two_One_Sum(_k, _1, _j, _l, _2, x1); \
387  Square(a1, _j, _1); \
388  Two_Two_Sum(_j, _1, _l, _2, x5, x4, x3, x2)
389 
390 static double splitter; /* = 2^ceiling(p / 2) + 1. Used to split floats in half. */
391 static double epsilon; /* = 2^(-p). Used to estimate round-off errors. */
392 /* A set of coefficients used to calculate maximum round-off errors. */
393 static double resulterrbound;
398 
416 void exactinit()
417 {
418  double half;
419  double check, lastcheck;
420  int every_other;
421 
422  every_other = 1;
423  half = 0.5;
424  epsilon = 1.0;
425  splitter = 1.0;
426  check = 1.0;
427  /* Repeatedly divide `epsilon' by two until it is too small to add to
428  * one without causing round-off. (Also check if the sum is equal to
429  * the previous sum, for machines that round up instead of using exact
430  * rounding. Not that this library will work on such machines anyway. */
431  do {
432  lastcheck = check;
433  epsilon *= half;
434  if (every_other) {
435  splitter *= 2.0;
436  }
437  every_other = !every_other;
438  check = 1.0 + epsilon;
439  } while (!ELEM(check, 1.0, lastcheck));
440  splitter += 1.0;
441 
442  /* Error bounds for orientation and #incircle tests. */
443  resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
444  ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon;
445  ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon;
446  ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon;
447  o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon;
448  o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon;
449  o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;
450  iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon;
451  iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon;
452  iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon;
453  isperrboundA = (16.0 + 224.0 * epsilon) * epsilon;
454  isperrboundB = (5.0 + 72.0 * epsilon) * epsilon;
455  isperrboundC = (71.0 + 1408.0 * epsilon) * epsilon * epsilon;
456 }
457 
466  int elen, const double *e, int flen, const double *f, double *h)
467 {
468  double Q;
469  INEXACT double Qnew;
470  INEXACT double hh;
471  INEXACT double bvirt;
472  double avirt, bround, around;
473  int eindex, findex, hindex;
474  double enow, fnow;
475 
476  enow = e[0];
477  fnow = f[0];
478  eindex = findex = 0;
479  if ((fnow > enow) == (fnow > -enow)) {
480  Q = enow;
481  enow = e[++eindex];
482  }
483  else {
484  Q = fnow;
485  fnow = f[++findex];
486  }
487  hindex = 0;
488  if ((eindex < elen) && (findex < flen)) {
489  if ((fnow > enow) == (fnow > -enow)) {
490  Fast_Two_Sum(enow, Q, Qnew, hh);
491  enow = e[++eindex];
492  }
493  else {
494  Fast_Two_Sum(fnow, Q, Qnew, hh);
495  fnow = f[++findex];
496  }
497  Q = Qnew;
498  if (hh != 0.0) {
499  h[hindex++] = hh;
500  }
501  while ((eindex < elen) && (findex < flen)) {
502  if ((fnow > enow) == (fnow > -enow)) {
503  Two_Sum(Q, enow, Qnew, hh);
504  enow = e[++eindex];
505  }
506  else {
507  Two_Sum(Q, fnow, Qnew, hh);
508  fnow = f[++findex];
509  }
510  Q = Qnew;
511  if (hh != 0.0) {
512  h[hindex++] = hh;
513  }
514  }
515  }
516  while (eindex < elen) {
517  Two_Sum(Q, enow, Qnew, hh);
518  enow = e[++eindex];
519  Q = Qnew;
520  if (hh != 0.0) {
521  h[hindex++] = hh;
522  }
523  }
524  while (findex < flen) {
525  Two_Sum(Q, fnow, Qnew, hh);
526  fnow = f[++findex];
527  Q = Qnew;
528  if (hh != 0.0) {
529  h[hindex++] = hh;
530  }
531  }
532  if ((Q != 0.0) || (hindex == 0)) {
533  h[hindex++] = Q;
534  }
535  return hindex;
536 }
537 
538 /* scale_expansion_zeroelim() Multiply an expansion by a scalar,
539  * eliminating zero components from the
540  * output expansion.
541  *
542  * Sets h = be. See either version of my paper for details.
543  * e and h cannot be the same.
544  */
545 static int scale_expansion_zeroelim(int elen, const double *e, double b, double *h)
546 {
547  INEXACT double Q, sum;
548  double hh;
549  INEXACT double product1;
550  double product0;
551  int eindex, hindex;
552  double enow;
553  INEXACT double bvirt;
554  double avirt, bround, around;
555  INEXACT double c;
556  INEXACT double abig;
557  double ahi, alo, bhi, blo;
558  double err1, err2, err3;
559 
560  Split(b, bhi, blo);
561  Two_Product_Presplit(e[0], b, bhi, blo, Q, hh);
562  hindex = 0;
563  if (hh != 0) {
564  h[hindex++] = hh;
565  }
566  for (eindex = 1; eindex < elen; eindex++) {
567  enow = e[eindex];
568  Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
569  Two_Sum(Q, product0, sum, hh);
570  if (hh != 0) {
571  h[hindex++] = hh;
572  }
573  Fast_Two_Sum(product1, sum, Q, hh);
574  if (hh != 0) {
575  h[hindex++] = hh;
576  }
577  }
578  if ((Q != 0.0) || (hindex == 0)) {
579  h[hindex++] = Q;
580  }
581  return hindex;
582 }
583 
584 /* estimate() Produce a one-word estimate of an expansion's value. */
585 static double estimate(int elen, const double *e)
586 {
587  double Q;
588  int eindex;
589 
590  Q = e[0];
591  for (eindex = 1; eindex < elen; eindex++) {
592  Q += e[eindex];
593  }
594  return Q;
595 }
596 
615 double orient2dfast(const double *pa, const double *pb, const double *pc)
616 {
617  double acx, bcx, acy, bcy;
618 
619  acx = pa[0] - pc[0];
620  bcx = pb[0] - pc[0];
621  acy = pa[1] - pc[1];
622  bcy = pb[1] - pc[1];
623  return acx * bcy - acy * bcx;
624 }
625 
626 static double orient2dadapt(const double *pa, const double *pb, const double *pc, double detsum)
627 {
628  INEXACT double acx, acy, bcx, bcy;
629  double acxtail, acytail, bcxtail, bcytail;
630  INEXACT double detleft, detright;
631  double detlefttail, detrighttail;
632  double det, errbound;
633  double B[4], C1[8], C2[12], D[16];
634  INEXACT double B3;
635  int C1length, C2length, Dlength;
636  double u[4];
637  INEXACT double u3;
638  INEXACT double s1, t1;
639  double s0, t0;
640 
641  INEXACT double bvirt;
642  double avirt, bround, around;
643  INEXACT double c;
644  INEXACT double abig;
645  double ahi, alo, bhi, blo;
646  double err1, err2, err3;
647  INEXACT double _i, _j;
648  double _0;
649 
650  acx = (double)(pa[0] - pc[0]);
651  bcx = (double)(pb[0] - pc[0]);
652  acy = (double)(pa[1] - pc[1]);
653  bcy = (double)(pb[1] - pc[1]);
654 
655  Two_Product(acx, bcy, detleft, detlefttail);
656  Two_Product(acy, bcx, detright, detrighttail);
657 
658  Two_Two_Diff(detleft, detlefttail, detright, detrighttail, B3, B[2], B[1], B[0]);
659  B[3] = B3;
660 
661  det = estimate(4, B);
662  errbound = ccwerrboundB * detsum;
663  if ((det >= errbound) || (-det >= errbound)) {
664  return det;
665  }
666 
667  Two_Diff_Tail(pa[0], pc[0], acx, acxtail);
668  Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail);
669  Two_Diff_Tail(pa[1], pc[1], acy, acytail);
670  Two_Diff_Tail(pb[1], pc[1], bcy, bcytail);
671 
672  if ((acxtail == 0.0) && (acytail == 0.0) && (bcxtail == 0.0) && (bcytail == 0.0)) {
673  return det;
674  }
675 
676  errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det);
677  det += (acx * bcytail + bcy * acxtail) - (acy * bcxtail + bcx * acytail);
678  if ((det >= errbound) || (-det >= errbound)) {
679  return det;
680  }
681 
682  Two_Product(acxtail, bcy, s1, s0);
683  Two_Product(acytail, bcx, t1, t0);
684  Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
685  u[3] = u3;
686  C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1);
687 
688  Two_Product(acx, bcytail, s1, s0);
689  Two_Product(acy, bcxtail, t1, t0);
690  Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
691  u[3] = u3;
692  C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2);
693 
694  Two_Product(acxtail, bcytail, s1, s0);
695  Two_Product(acytail, bcxtail, t1, t0);
696  Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
697  u[3] = u3;
698  Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D);
699 
700  return (D[Dlength - 1]);
701 }
702 
703 double orient2d(const double *pa, const double *pb, const double *pc)
704 {
705  double detleft, detright, det;
706  double detsum, errbound;
707 
708  detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]);
709  detright = (pa[1] - pc[1]) * (pb[0] - pc[0]);
710  det = detleft - detright;
711 
712  if (detleft > 0.0) {
713  if (detright <= 0.0) {
714  return det;
715  }
716  detsum = detleft + detright;
717  }
718  else if (detleft < 0.0) {
719  if (detright >= 0.0) {
720  return det;
721  }
722  detsum = -detleft - detright;
723  }
724  else {
725  return det;
726  }
727 
728  errbound = ccwerrboundA * detsum;
729  if ((det >= errbound) || (-det >= errbound)) {
730  return det;
731  }
732 
733  return orient2dadapt(pa, pb, pc, detsum);
734 }
735 
758 double orient3dfast(const double *pa, const double *pb, const double *pc, const double *pd)
759 {
760  double adx, bdx, cdx;
761  double ady, bdy, cdy;
762  double adz, bdz, cdz;
763 
764  adx = pa[0] - pd[0];
765  bdx = pb[0] - pd[0];
766  cdx = pc[0] - pd[0];
767  ady = pa[1] - pd[1];
768  bdy = pb[1] - pd[1];
769  cdy = pc[1] - pd[1];
770  adz = pa[2] - pd[2];
771  bdz = pb[2] - pd[2];
772  cdz = pc[2] - pd[2];
773 
774  return adx * (bdy * cdz - bdz * cdy) + bdx * (cdy * adz - cdz * ady) +
775  cdx * (ady * bdz - adz * bdy);
776 }
777 
782 /* NOLINTNEXTLINE: readability-function-size */
783 static double orient3dadapt(
784  const double *pa, const double *pb, const double *pc, const double *pd, double permanent)
785 {
786  INEXACT double adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
787  double det, errbound;
788 
789  INEXACT double bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
790  double bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
791  double bc[4], ca[4], ab[4];
792  INEXACT double bc3, ca3, ab3;
793  double adet[8], bdet[8], cdet[8];
794  int alen, blen, clen;
795  double abdet[16];
796  int ablen;
797  double *finnow, *finother, *finswap;
798  double fin1[192], fin2[192];
799  int finlength;
800 
801  double adxtail, bdxtail, cdxtail;
802  double adytail, bdytail, cdytail;
803  double adztail, bdztail, cdztail;
804  INEXACT double at_blarge, at_clarge;
805  INEXACT double bt_clarge, bt_alarge;
806  INEXACT double ct_alarge, ct_blarge;
807  double at_b[4], at_c[4], bt_c[4], bt_a[4], ct_a[4], ct_b[4];
808  int at_blen, at_clen, bt_clen, bt_alen, ct_alen, ct_blen;
809  INEXACT double bdxt_cdy1, cdxt_bdy1, cdxt_ady1;
810  INEXACT double adxt_cdy1, adxt_bdy1, bdxt_ady1;
811  double bdxt_cdy0, cdxt_bdy0, cdxt_ady0;
812  double adxt_cdy0, adxt_bdy0, bdxt_ady0;
813  INEXACT double bdyt_cdx1, cdyt_bdx1, cdyt_adx1;
814  INEXACT double adyt_cdx1, adyt_bdx1, bdyt_adx1;
815  double bdyt_cdx0, cdyt_bdx0, cdyt_adx0;
816  double adyt_cdx0, adyt_bdx0, bdyt_adx0;
817  double bct[8], cat[8], abt[8];
818  int bctlen, catlen, abtlen;
819  INEXACT double bdxt_cdyt1, cdxt_bdyt1, cdxt_adyt1;
820  INEXACT double adxt_cdyt1, adxt_bdyt1, bdxt_adyt1;
821  double bdxt_cdyt0, cdxt_bdyt0, cdxt_adyt0;
822  double adxt_cdyt0, adxt_bdyt0, bdxt_adyt0;
823  double u[4], v[12], w[16];
824  INEXACT double u3;
825  int vlength, wlength;
826  double negate;
827 
828  INEXACT double bvirt;
829  double avirt, bround, around;
830  INEXACT double c;
831  INEXACT double abig;
832  double ahi, alo, bhi, blo;
833  double err1, err2, err3;
834  INEXACT double _i, _j, _k;
835  double _0;
836 
837  adx = (double)(pa[0] - pd[0]);
838  bdx = (double)(pb[0] - pd[0]);
839  cdx = (double)(pc[0] - pd[0]);
840  ady = (double)(pa[1] - pd[1]);
841  bdy = (double)(pb[1] - pd[1]);
842  cdy = (double)(pc[1] - pd[1]);
843  adz = (double)(pa[2] - pd[2]);
844  bdz = (double)(pb[2] - pd[2]);
845  cdz = (double)(pc[2] - pd[2]);
846 
847  Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
848  Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
849  Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
850  bc[3] = bc3;
851  alen = scale_expansion_zeroelim(4, bc, adz, adet);
852 
853  Two_Product(cdx, ady, cdxady1, cdxady0);
854  Two_Product(adx, cdy, adxcdy1, adxcdy0);
855  Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
856  ca[3] = ca3;
857  blen = scale_expansion_zeroelim(4, ca, bdz, bdet);
858 
859  Two_Product(adx, bdy, adxbdy1, adxbdy0);
860  Two_Product(bdx, ady, bdxady1, bdxady0);
861  Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
862  ab[3] = ab3;
863  clen = scale_expansion_zeroelim(4, ab, cdz, cdet);
864 
865  ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
866  finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
867 
868  det = estimate(finlength, fin1);
869  errbound = o3derrboundB * permanent;
870  if ((det >= errbound) || (-det >= errbound)) {
871  return det;
872  }
873 
874  Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
875  Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
876  Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
877  Two_Diff_Tail(pa[1], pd[1], ady, adytail);
878  Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
879  Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
880  Two_Diff_Tail(pa[2], pd[2], adz, adztail);
881  Two_Diff_Tail(pb[2], pd[2], bdz, bdztail);
882  Two_Diff_Tail(pc[2], pd[2], cdz, cdztail);
883 
884  if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) && (adytail == 0.0) &&
885  (bdytail == 0.0) && (cdytail == 0.0) && (adztail == 0.0) && (bdztail == 0.0) &&
886  (cdztail == 0.0)) {
887  return det;
888  }
889 
890  errbound = o3derrboundC * permanent + resulterrbound * Absolute(det);
891  det += (adz * ((bdx * cdytail + cdy * bdxtail) - (bdy * cdxtail + cdx * bdytail)) +
892  adztail * (bdx * cdy - bdy * cdx)) +
893  (bdz * ((cdx * adytail + ady * cdxtail) - (cdy * adxtail + adx * cdytail)) +
894  bdztail * (cdx * ady - cdy * adx)) +
895  (cdz * ((adx * bdytail + bdy * adxtail) - (ady * bdxtail + bdx * adytail)) +
896  cdztail * (adx * bdy - ady * bdx));
897  if ((det >= errbound) || (-det >= errbound)) {
898  return det;
899  }
900 
901  finnow = fin1;
902  finother = fin2;
903 
904  if (adxtail == 0.0) {
905  if (adytail == 0.0) {
906  at_b[0] = 0.0;
907  at_blen = 1;
908  at_c[0] = 0.0;
909  at_clen = 1;
910  }
911  else {
912  negate = -adytail;
913  Two_Product(negate, bdx, at_blarge, at_b[0]);
914  at_b[1] = at_blarge;
915  at_blen = 2;
916  Two_Product(adytail, cdx, at_clarge, at_c[0]);
917  at_c[1] = at_clarge;
918  at_clen = 2;
919  }
920  }
921  else {
922  if (adytail == 0.0) {
923  Two_Product(adxtail, bdy, at_blarge, at_b[0]);
924  at_b[1] = at_blarge;
925  at_blen = 2;
926  negate = -adxtail;
927  Two_Product(negate, cdy, at_clarge, at_c[0]);
928  at_c[1] = at_clarge;
929  at_clen = 2;
930  }
931  else {
932  Two_Product(adxtail, bdy, adxt_bdy1, adxt_bdy0);
933  Two_Product(adytail, bdx, adyt_bdx1, adyt_bdx0);
934  Two_Two_Diff(
935  adxt_bdy1, adxt_bdy0, adyt_bdx1, adyt_bdx0, at_blarge, at_b[2], at_b[1], at_b[0]);
936  at_b[3] = at_blarge;
937  at_blen = 4;
938  Two_Product(adytail, cdx, adyt_cdx1, adyt_cdx0);
939  Two_Product(adxtail, cdy, adxt_cdy1, adxt_cdy0);
940  Two_Two_Diff(
941  adyt_cdx1, adyt_cdx0, adxt_cdy1, adxt_cdy0, at_clarge, at_c[2], at_c[1], at_c[0]);
942  at_c[3] = at_clarge;
943  at_clen = 4;
944  }
945  }
946  if (bdxtail == 0.0) {
947  if (bdytail == 0.0) {
948  bt_c[0] = 0.0;
949  bt_clen = 1;
950  bt_a[0] = 0.0;
951  bt_alen = 1;
952  }
953  else {
954  negate = -bdytail;
955  Two_Product(negate, cdx, bt_clarge, bt_c[0]);
956  bt_c[1] = bt_clarge;
957  bt_clen = 2;
958  Two_Product(bdytail, adx, bt_alarge, bt_a[0]);
959  bt_a[1] = bt_alarge;
960  bt_alen = 2;
961  }
962  }
963  else {
964  if (bdytail == 0.0) {
965  Two_Product(bdxtail, cdy, bt_clarge, bt_c[0]);
966  bt_c[1] = bt_clarge;
967  bt_clen = 2;
968  negate = -bdxtail;
969  Two_Product(negate, ady, bt_alarge, bt_a[0]);
970  bt_a[1] = bt_alarge;
971  bt_alen = 2;
972  }
973  else {
974  Two_Product(bdxtail, cdy, bdxt_cdy1, bdxt_cdy0);
975  Two_Product(bdytail, cdx, bdyt_cdx1, bdyt_cdx0);
976  Two_Two_Diff(
977  bdxt_cdy1, bdxt_cdy0, bdyt_cdx1, bdyt_cdx0, bt_clarge, bt_c[2], bt_c[1], bt_c[0]);
978  bt_c[3] = bt_clarge;
979  bt_clen = 4;
980  Two_Product(bdytail, adx, bdyt_adx1, bdyt_adx0);
981  Two_Product(bdxtail, ady, bdxt_ady1, bdxt_ady0);
982  Two_Two_Diff(
983  bdyt_adx1, bdyt_adx0, bdxt_ady1, bdxt_ady0, bt_alarge, bt_a[2], bt_a[1], bt_a[0]);
984  bt_a[3] = bt_alarge;
985  bt_alen = 4;
986  }
987  }
988  if (cdxtail == 0.0) {
989  if (cdytail == 0.0) {
990  ct_a[0] = 0.0;
991  ct_alen = 1;
992  ct_b[0] = 0.0;
993  ct_blen = 1;
994  }
995  else {
996  negate = -cdytail;
997  Two_Product(negate, adx, ct_alarge, ct_a[0]);
998  ct_a[1] = ct_alarge;
999  ct_alen = 2;
1000  Two_Product(cdytail, bdx, ct_blarge, ct_b[0]);
1001  ct_b[1] = ct_blarge;
1002  ct_blen = 2;
1003  }
1004  }
1005  else {
1006  if (cdytail == 0.0) {
1007  Two_Product(cdxtail, ady, ct_alarge, ct_a[0]);
1008  ct_a[1] = ct_alarge;
1009  ct_alen = 2;
1010  negate = -cdxtail;
1011  Two_Product(negate, bdy, ct_blarge, ct_b[0]);
1012  ct_b[1] = ct_blarge;
1013  ct_blen = 2;
1014  }
1015  else {
1016  Two_Product(cdxtail, ady, cdxt_ady1, cdxt_ady0);
1017  Two_Product(cdytail, adx, cdyt_adx1, cdyt_adx0);
1018  Two_Two_Diff(
1019  cdxt_ady1, cdxt_ady0, cdyt_adx1, cdyt_adx0, ct_alarge, ct_a[2], ct_a[1], ct_a[0]);
1020  ct_a[3] = ct_alarge;
1021  ct_alen = 4;
1022  Two_Product(cdytail, bdx, cdyt_bdx1, cdyt_bdx0);
1023  Two_Product(cdxtail, bdy, cdxt_bdy1, cdxt_bdy0);
1024  Two_Two_Diff(
1025  cdyt_bdx1, cdyt_bdx0, cdxt_bdy1, cdxt_bdy0, ct_blarge, ct_b[2], ct_b[1], ct_b[0]);
1026  ct_b[3] = ct_blarge;
1027  ct_blen = 4;
1028  }
1029  }
1030 
1031  bctlen = fast_expansion_sum_zeroelim(bt_clen, bt_c, ct_blen, ct_b, bct);
1032  wlength = scale_expansion_zeroelim(bctlen, bct, adz, w);
1033  finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, finother);
1034  finswap = finnow;
1035  finnow = finother;
1036  finother = finswap;
1037 
1038  catlen = fast_expansion_sum_zeroelim(ct_alen, ct_a, at_clen, at_c, cat);
1039  wlength = scale_expansion_zeroelim(catlen, cat, bdz, w);
1040  finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, finother);
1041  finswap = finnow;
1042  finnow = finother;
1043  finother = finswap;
1044 
1045  abtlen = fast_expansion_sum_zeroelim(at_blen, at_b, bt_alen, bt_a, abt);
1046  wlength = scale_expansion_zeroelim(abtlen, abt, cdz, w);
1047  finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, finother);
1048  finswap = finnow;
1049  finnow = finother;
1050  finother = finswap;
1051 
1052  if (adztail != 0.0) {
1053  vlength = scale_expansion_zeroelim(4, bc, adztail, v);
1054  finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v, finother);
1055  finswap = finnow;
1056  finnow = finother;
1057  finother = finswap;
1058  }
1059  if (bdztail != 0.0) {
1060  vlength = scale_expansion_zeroelim(4, ca, bdztail, v);
1061  finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v, finother);
1062  finswap = finnow;
1063  finnow = finother;
1064  finother = finswap;
1065  }
1066  if (cdztail != 0.0) {
1067  vlength = scale_expansion_zeroelim(4, ab, cdztail, v);
1068  finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v, finother);
1069  finswap = finnow;
1070  finnow = finother;
1071  finother = finswap;
1072  }
1073 
1074  if (adxtail != 0.0) {
1075  if (bdytail != 0.0) {
1076  Two_Product(adxtail, bdytail, adxt_bdyt1, adxt_bdyt0);
1077  Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdz, u3, u[2], u[1], u[0]);
1078  u[3] = u3;
1079  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, finother);
1080  finswap = finnow;
1081  finnow = finother;
1082  finother = finswap;
1083  if (cdztail != 0.0) {
1084  Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdztail, u3, u[2], u[1], u[0]);
1085  u[3] = u3;
1086  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, finother);
1087  finswap = finnow;
1088  finnow = finother;
1089  finother = finswap;
1090  }
1091  }
1092  if (cdytail != 0.0) {
1093  negate = -adxtail;
1094  Two_Product(negate, cdytail, adxt_cdyt1, adxt_cdyt0);
1095  Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdz, u3, u[2], u[1], u[0]);
1096  u[3] = u3;
1097  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, finother);
1098  finswap = finnow;
1099  finnow = finother;
1100  finother = finswap;
1101  if (bdztail != 0.0) {
1102  Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdztail, u3, u[2], u[1], u[0]);
1103  u[3] = u3;
1104  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, finother);
1105  finswap = finnow;
1106  finnow = finother;
1107  finother = finswap;
1108  }
1109  }
1110  }
1111  if (bdxtail != 0.0) {
1112  if (cdytail != 0.0) {
1113  Two_Product(bdxtail, cdytail, bdxt_cdyt1, bdxt_cdyt0);
1114  Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adz, u3, u[2], u[1], u[0]);
1115  u[3] = u3;
1116  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, finother);
1117  finswap = finnow;
1118  finnow = finother;
1119  finother = finswap;
1120  if (adztail != 0.0) {
1121  Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adztail, u3, u[2], u[1], u[0]);
1122  u[3] = u3;
1123  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, finother);
1124  finswap = finnow;
1125  finnow = finother;
1126  finother = finswap;
1127  }
1128  }
1129  if (adytail != 0.0) {
1130  negate = -bdxtail;
1131  Two_Product(negate, adytail, bdxt_adyt1, bdxt_adyt0);
1132  Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdz, u3, u[2], u[1], u[0]);
1133  u[3] = u3;
1134  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, finother);
1135  finswap = finnow;
1136  finnow = finother;
1137  finother = finswap;
1138  if (cdztail != 0.0) {
1139  Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdztail, u3, u[2], u[1], u[0]);
1140  u[3] = u3;
1141  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, finother);
1142  finswap = finnow;
1143  finnow = finother;
1144  finother = finswap;
1145  }
1146  }
1147  }
1148  if (cdxtail != 0.0) {
1149  if (adytail != 0.0) {
1150  Two_Product(cdxtail, adytail, cdxt_adyt1, cdxt_adyt0);
1151  Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdz, u3, u[2], u[1], u[0]);
1152  u[3] = u3;
1153  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, finother);
1154  finswap = finnow;
1155  finnow = finother;
1156  finother = finswap;
1157  if (bdztail != 0.0) {
1158  Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdztail, u3, u[2], u[1], u[0]);
1159  u[3] = u3;
1160  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, finother);
1161  finswap = finnow;
1162  finnow = finother;
1163  finother = finswap;
1164  }
1165  }
1166  if (bdytail != 0.0) {
1167  negate = -cdxtail;
1168  Two_Product(negate, bdytail, cdxt_bdyt1, cdxt_bdyt0);
1169  Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adz, u3, u[2], u[1], u[0]);
1170  u[3] = u3;
1171  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, finother);
1172  finswap = finnow;
1173  finnow = finother;
1174  finother = finswap;
1175  if (adztail != 0.0) {
1176  Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adztail, u3, u[2], u[1], u[0]);
1177  u[3] = u3;
1178  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u, finother);
1179  finswap = finnow;
1180  finnow = finother;
1181  finother = finswap;
1182  }
1183  }
1184  }
1185 
1186  if (adztail != 0.0) {
1187  wlength = scale_expansion_zeroelim(bctlen, bct, adztail, w);
1188  finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, finother);
1189  finswap = finnow;
1190  finnow = finother;
1191  finother = finswap;
1192  }
1193  if (bdztail != 0.0) {
1194  wlength = scale_expansion_zeroelim(catlen, cat, bdztail, w);
1195  finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, finother);
1196  finswap = finnow;
1197  finnow = finother;
1198  finother = finswap;
1199  }
1200  if (cdztail != 0.0) {
1201  wlength = scale_expansion_zeroelim(abtlen, abt, cdztail, w);
1202  finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w, finother);
1203  finswap = finnow;
1204  finnow = finother;
1205  finother = finswap;
1206  }
1207 
1208  return finnow[finlength - 1];
1209 }
1210 
1211 double orient3d(const double *pa, const double *pb, const double *pc, const double *pd)
1212 {
1213  double adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
1214  double bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
1215  double det;
1216  double permanent, errbound;
1217 
1218  adx = pa[0] - pd[0];
1219  bdx = pb[0] - pd[0];
1220  cdx = pc[0] - pd[0];
1221  ady = pa[1] - pd[1];
1222  bdy = pb[1] - pd[1];
1223  cdy = pc[1] - pd[1];
1224  adz = pa[2] - pd[2];
1225  bdz = pb[2] - pd[2];
1226  cdz = pc[2] - pd[2];
1227 
1228  bdxcdy = bdx * cdy;
1229  cdxbdy = cdx * bdy;
1230 
1231  cdxady = cdx * ady;
1232  adxcdy = adx * cdy;
1233 
1234  adxbdy = adx * bdy;
1235  bdxady = bdx * ady;
1236 
1237  det = adz * (bdxcdy - cdxbdy) + bdz * (cdxady - adxcdy) + cdz * (adxbdy - bdxady);
1238 
1239  permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * Absolute(adz) +
1240  (Absolute(cdxady) + Absolute(adxcdy)) * Absolute(bdz) +
1241  (Absolute(adxbdy) + Absolute(bdxady)) * Absolute(cdz);
1242  errbound = o3derrboundA * permanent;
1243  if ((det > errbound) || (-det > errbound)) {
1244  return det;
1245  }
1246 
1247  return orient3dadapt(pa, pb, pc, pd, permanent);
1248 }
1249 
1269 double incirclefast(const double *pa, const double *pb, const double *pc, const double *pd)
1270 {
1271  double adx, ady, bdx, bdy, cdx, cdy;
1272  double abdet, bcdet, cadet;
1273  double alift, blift, clift;
1274 
1275  adx = pa[0] - pd[0];
1276  ady = pa[1] - pd[1];
1277  bdx = pb[0] - pd[0];
1278  bdy = pb[1] - pd[1];
1279  cdx = pc[0] - pd[0];
1280  cdy = pc[1] - pd[1];
1281 
1282  abdet = adx * bdy - bdx * ady;
1283  bcdet = bdx * cdy - cdx * bdy;
1284  cadet = cdx * ady - adx * cdy;
1285  alift = adx * adx + ady * ady;
1286  blift = bdx * bdx + bdy * bdy;
1287  clift = cdx * cdx + cdy * cdy;
1288 
1289  return alift * bcdet + blift * cadet + clift * abdet;
1290 }
1291 
1296 /* NOLINTNEXTLINE: readability-function-size */
1297 static double incircleadapt(
1298  const double *pa, const double *pb, const double *pc, const double *pd, double permanent)
1299 {
1300  INEXACT double adx, bdx, cdx, ady, bdy, cdy;
1301  double det, errbound;
1302 
1303  INEXACT double bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
1304  double bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
1305  double bc[4], ca[4], ab[4];
1306  INEXACT double bc3, ca3, ab3;
1307  double axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32];
1308  int axbclen, axxbclen, aybclen, ayybclen, alen;
1309  double bxca[8], bxxca[16], byca[8], byyca[16], bdet[32];
1310  int bxcalen, bxxcalen, bycalen, byycalen, blen;
1311  double cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32];
1312  int cxablen, cxxablen, cyablen, cyyablen, clen;
1313  double abdet[64];
1314  int ablen;
1315  double fin1[1152], fin2[1152];
1316  double *finnow, *finother, *finswap;
1317  int finlength;
1318 
1319  double adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
1320  INEXACT double adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1;
1321  double adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0;
1322  double aa[4], bb[4], cc[4];
1323  INEXACT double aa3, bb3, cc3;
1324  INEXACT double ti1, tj1;
1325  double ti0, tj0;
1326  double u[4], v[4];
1327  INEXACT double u3, v3;
1328  double temp8[8], temp16a[16], temp16b[16], temp16c[16];
1329  double temp32a[32], temp32b[32], temp48[48], temp64[64];
1330  int temp8len, temp16alen, temp16blen, temp16clen;
1331  int temp32alen, temp32blen, temp48len, temp64len;
1332  double axtbb[8], axtcc[8], aytbb[8], aytcc[8];
1333  int axtbblen, axtcclen, aytbblen, aytcclen;
1334  double bxtaa[8], bxtcc[8], bytaa[8], bytcc[8];
1335  int bxtaalen, bxtcclen, bytaalen, bytcclen;
1336  double cxtaa[8], cxtbb[8], cytaa[8], cytbb[8];
1337  int cxtaalen, cxtbblen, cytaalen, cytbblen;
1338  double axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8];
1339  int axtbclen, aytbclen, bxtcalen, bytcalen, cxtablen, cytablen;
1340  double axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16];
1341  int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen;
1342  double axtbctt[8], aytbctt[8], bxtcatt[8];
1343  double bytcatt[8], cxtabtt[8], cytabtt[8];
1344  int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen;
1345  double abt[8], bct[8], cat[8];
1346  int abtlen, bctlen, catlen;
1347  double abtt[4], bctt[4], catt[4];
1348  int abttlen, bcttlen, cattlen;
1349  INEXACT double abtt3, bctt3, catt3;
1350  double negate;
1351 
1352  INEXACT double bvirt;
1353  double avirt, bround, around;
1354  INEXACT double c;
1355  INEXACT double abig;
1356  double ahi, alo, bhi, blo;
1357  double err1, err2, err3;
1358  INEXACT double _i, _j;
1359  double _0;
1360 
1361  adx = (double)(pa[0] - pd[0]);
1362  bdx = (double)(pb[0] - pd[0]);
1363  cdx = (double)(pc[0] - pd[0]);
1364  ady = (double)(pa[1] - pd[1]);
1365  bdy = (double)(pb[1] - pd[1]);
1366  cdy = (double)(pc[1] - pd[1]);
1367 
1368  Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
1369  Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
1370  Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
1371  bc[3] = bc3;
1372  axbclen = scale_expansion_zeroelim(4, bc, adx, axbc);
1373  axxbclen = scale_expansion_zeroelim(axbclen, axbc, adx, axxbc);
1374  aybclen = scale_expansion_zeroelim(4, bc, ady, aybc);
1375  ayybclen = scale_expansion_zeroelim(aybclen, aybc, ady, ayybc);
1376  alen = fast_expansion_sum_zeroelim(axxbclen, axxbc, ayybclen, ayybc, adet);
1377 
1378  Two_Product(cdx, ady, cdxady1, cdxady0);
1379  Two_Product(adx, cdy, adxcdy1, adxcdy0);
1380  Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
1381  ca[3] = ca3;
1382  bxcalen = scale_expansion_zeroelim(4, ca, bdx, bxca);
1383  bxxcalen = scale_expansion_zeroelim(bxcalen, bxca, bdx, bxxca);
1384  bycalen = scale_expansion_zeroelim(4, ca, bdy, byca);
1385  byycalen = scale_expansion_zeroelim(bycalen, byca, bdy, byyca);
1386  blen = fast_expansion_sum_zeroelim(bxxcalen, bxxca, byycalen, byyca, bdet);
1387 
1388  Two_Product(adx, bdy, adxbdy1, adxbdy0);
1389  Two_Product(bdx, ady, bdxady1, bdxady0);
1390  Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
1391  ab[3] = ab3;
1392  cxablen = scale_expansion_zeroelim(4, ab, cdx, cxab);
1393  cxxablen = scale_expansion_zeroelim(cxablen, cxab, cdx, cxxab);
1394  cyablen = scale_expansion_zeroelim(4, ab, cdy, cyab);
1395  cyyablen = scale_expansion_zeroelim(cyablen, cyab, cdy, cyyab);
1396  clen = fast_expansion_sum_zeroelim(cxxablen, cxxab, cyyablen, cyyab, cdet);
1397 
1398  ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
1399  finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
1400 
1401  det = estimate(finlength, fin1);
1402  errbound = iccerrboundB * permanent;
1403  if ((det >= errbound) || (-det >= errbound)) {
1404  return det;
1405  }
1406 
1407  Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
1408  Two_Diff_Tail(pa[1], pd[1], ady, adytail);
1409  Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
1410  Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
1411  Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
1412  Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
1413  if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) && (adytail == 0.0) &&
1414  (bdytail == 0.0) && (cdytail == 0.0)) {
1415  return det;
1416  }
1417 
1418  errbound = iccerrboundC * permanent + resulterrbound * Absolute(det);
1419  det += ((adx * adx + ady * ady) *
1420  ((bdx * cdytail + cdy * bdxtail) - (bdy * cdxtail + cdx * bdytail)) +
1421  2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx)) +
1422  ((bdx * bdx + bdy * bdy) *
1423  ((cdx * adytail + ady * cdxtail) - (cdy * adxtail + adx * cdytail)) +
1424  2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx)) +
1425  ((cdx * cdx + cdy * cdy) *
1426  ((adx * bdytail + bdy * adxtail) - (ady * bdxtail + bdx * adytail)) +
1427  2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));
1428  if ((det >= errbound) || (-det >= errbound)) {
1429  return det;
1430  }
1431 
1432  finnow = fin1;
1433  finother = fin2;
1434 
1435  if ((bdxtail != 0.0) || (bdytail != 0.0) || (cdxtail != 0.0) || (cdytail != 0.0)) {
1436  Square(adx, adxadx1, adxadx0);
1437  Square(ady, adyady1, adyady0);
1438  Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]);
1439  aa[3] = aa3;
1440  }
1441  if ((cdxtail != 0.0) || (cdytail != 0.0) || (adxtail != 0.0) || (adytail != 0.0)) {
1442  Square(bdx, bdxbdx1, bdxbdx0);
1443  Square(bdy, bdybdy1, bdybdy0);
1444  Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]);
1445  bb[3] = bb3;
1446  }
1447  if ((adxtail != 0.0) || (adytail != 0.0) || (bdxtail != 0.0) || (bdytail != 0.0)) {
1448  Square(cdx, cdxcdx1, cdxcdx0);
1449  Square(cdy, cdycdy1, cdycdy0);
1450  Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]);
1451  cc[3] = cc3;
1452  }
1453 
1454  if (adxtail != 0.0) {
1455  axtbclen = scale_expansion_zeroelim(4, bc, adxtail, axtbc);
1456  temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, 2.0 * adx, temp16a);
1457 
1458  axtcclen = scale_expansion_zeroelim(4, cc, adxtail, axtcc);
1459  temp16blen = scale_expansion_zeroelim(axtcclen, axtcc, bdy, temp16b);
1460 
1461  axtbblen = scale_expansion_zeroelim(4, bb, adxtail, axtbb);
1462  temp16clen = scale_expansion_zeroelim(axtbblen, axtbb, -cdy, temp16c);
1463 
1464  temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
1465  temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
1466  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
1467  finswap = finnow;
1468  finnow = finother;
1469  finother = finswap;
1470  }
1471  if (adytail != 0.0) {
1472  aytbclen = scale_expansion_zeroelim(4, bc, adytail, aytbc);
1473  temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, 2.0 * ady, temp16a);
1474 
1475  aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb);
1476  temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b);
1477 
1478  aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc);
1479  temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c);
1480 
1481  temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
1482  temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
1483  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
1484  finswap = finnow;
1485  finnow = finother;
1486  finother = finswap;
1487  }
1488  if (bdxtail != 0.0) {
1489  bxtcalen = scale_expansion_zeroelim(4, ca, bdxtail, bxtca);
1490  temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, 2.0 * bdx, temp16a);
1491 
1492  bxtaalen = scale_expansion_zeroelim(4, aa, bdxtail, bxtaa);
1493  temp16blen = scale_expansion_zeroelim(bxtaalen, bxtaa, cdy, temp16b);
1494 
1495  bxtcclen = scale_expansion_zeroelim(4, cc, bdxtail, bxtcc);
1496  temp16clen = scale_expansion_zeroelim(bxtcclen, bxtcc, -ady, temp16c);
1497 
1498  temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
1499  temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
1500  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
1501  finswap = finnow;
1502  finnow = finother;
1503  finother = finswap;
1504  }
1505  if (bdytail != 0.0) {
1506  bytcalen = scale_expansion_zeroelim(4, ca, bdytail, bytca);
1507  temp16alen = scale_expansion_zeroelim(bytcalen, bytca, 2.0 * bdy, temp16a);
1508 
1509  bytcclen = scale_expansion_zeroelim(4, cc, bdytail, bytcc);
1510  temp16blen = scale_expansion_zeroelim(bytcclen, bytcc, adx, temp16b);
1511 
1512  bytaalen = scale_expansion_zeroelim(4, aa, bdytail, bytaa);
1513  temp16clen = scale_expansion_zeroelim(bytaalen, bytaa, -cdx, temp16c);
1514 
1515  temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
1516  temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
1517  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
1518  finswap = finnow;
1519  finnow = finother;
1520  finother = finswap;
1521  }
1522  if (cdxtail != 0.0) {
1523  cxtablen = scale_expansion_zeroelim(4, ab, cdxtail, cxtab);
1524  temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, 2.0 * cdx, temp16a);
1525 
1526  cxtbblen = scale_expansion_zeroelim(4, bb, cdxtail, cxtbb);
1527  temp16blen = scale_expansion_zeroelim(cxtbblen, cxtbb, ady, temp16b);
1528 
1529  cxtaalen = scale_expansion_zeroelim(4, aa, cdxtail, cxtaa);
1530  temp16clen = scale_expansion_zeroelim(cxtaalen, cxtaa, -bdy, temp16c);
1531 
1532  temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
1533  temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
1534  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
1535  finswap = finnow;
1536  finnow = finother;
1537  finother = finswap;
1538  }
1539  if (cdytail != 0.0) {
1540  cytablen = scale_expansion_zeroelim(4, ab, cdytail, cytab);
1541  temp16alen = scale_expansion_zeroelim(cytablen, cytab, 2.0 * cdy, temp16a);
1542 
1543  cytaalen = scale_expansion_zeroelim(4, aa, cdytail, cytaa);
1544  temp16blen = scale_expansion_zeroelim(cytaalen, cytaa, bdx, temp16b);
1545 
1546  cytbblen = scale_expansion_zeroelim(4, bb, cdytail, cytbb);
1547  temp16clen = scale_expansion_zeroelim(cytbblen, cytbb, -adx, temp16c);
1548 
1549  temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
1550  temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
1551  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
1552  finswap = finnow;
1553  finnow = finother;
1554  finother = finswap;
1555  }
1556 
1557  if ((adxtail != 0.0) || (adytail != 0.0)) {
1558  if ((bdxtail != 0.0) || (bdytail != 0.0) || (cdxtail != 0.0) || (cdytail != 0.0)) {
1559  Two_Product(bdxtail, cdy, ti1, ti0);
1560  Two_Product(bdx, cdytail, tj1, tj0);
1561  Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
1562  u[3] = u3;
1563  negate = -bdy;
1564  Two_Product(cdxtail, negate, ti1, ti0);
1565  negate = -bdytail;
1566  Two_Product(cdx, negate, tj1, tj0);
1567  Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
1568  v[3] = v3;
1569  bctlen = fast_expansion_sum_zeroelim(4, u, 4, v, bct);
1570 
1571  Two_Product(bdxtail, cdytail, ti1, ti0);
1572  Two_Product(cdxtail, bdytail, tj1, tj0);
1573  Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]);
1574  bctt[3] = bctt3;
1575  bcttlen = 4;
1576  }
1577  else {
1578  bct[0] = 0.0;
1579  bctlen = 1;
1580  bctt[0] = 0.0;
1581  bcttlen = 1;
1582  }
1583 
1584  if (adxtail != 0.0) {
1585  temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, adxtail, temp16a);
1586  axtbctlen = scale_expansion_zeroelim(bctlen, bct, adxtail, axtbct);
1587  temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, 2.0 * adx, temp32a);
1588  temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
1589  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
1590  finswap = finnow;
1591  finnow = finother;
1592  finother = finswap;
1593  if (bdytail != 0.0) {
1594  temp8len = scale_expansion_zeroelim(4, cc, adxtail, temp8);
1595  temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail, temp16a);
1596  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
1597  finswap = finnow;
1598  finnow = finother;
1599  finother = finswap;
1600  }
1601  if (cdytail != 0.0) {
1602  temp8len = scale_expansion_zeroelim(4, bb, -adxtail, temp8);
1603  temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail, temp16a);
1604  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
1605  finswap = finnow;
1606  finnow = finother;
1607  finother = finswap;
1608  }
1609 
1610  temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, adxtail, temp32a);
1611  axtbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adxtail, axtbctt);
1612  temp16alen = scale_expansion_zeroelim(axtbcttlen, axtbctt, 2.0 * adx, temp16a);
1613  temp16blen = scale_expansion_zeroelim(axtbcttlen, axtbctt, adxtail, temp16b);
1614  temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
1615  temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
1616  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
1617  finswap = finnow;
1618  finnow = finother;
1619  finother = finswap;
1620  }
1621  if (adytail != 0.0) {
1622  temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, adytail, temp16a);
1623  aytbctlen = scale_expansion_zeroelim(bctlen, bct, adytail, aytbct);
1624  temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, 2.0 * ady, temp32a);
1625  temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
1626  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
1627  finswap = finnow;
1628  finnow = finother;
1629  finother = finswap;
1630 
1631  temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, adytail, temp32a);
1632  aytbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adytail, aytbctt);
1633  temp16alen = scale_expansion_zeroelim(aytbcttlen, aytbctt, 2.0 * ady, temp16a);
1634  temp16blen = scale_expansion_zeroelim(aytbcttlen, aytbctt, adytail, temp16b);
1635  temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
1636  temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
1637  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
1638  finswap = finnow;
1639  finnow = finother;
1640  finother = finswap;
1641  }
1642  }
1643  if ((bdxtail != 0.0) || (bdytail != 0.0)) {
1644  if ((cdxtail != 0.0) || (cdytail != 0.0) || (adxtail != 0.0) || (adytail != 0.0)) {
1645  Two_Product(cdxtail, ady, ti1, ti0);
1646  Two_Product(cdx, adytail, tj1, tj0);
1647  Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
1648  u[3] = u3;
1649  negate = -cdy;
1650  Two_Product(adxtail, negate, ti1, ti0);
1651  negate = -cdytail;
1652  Two_Product(adx, negate, tj1, tj0);
1653  Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
1654  v[3] = v3;
1655  catlen = fast_expansion_sum_zeroelim(4, u, 4, v, cat);
1656 
1657  Two_Product(cdxtail, adytail, ti1, ti0);
1658  Two_Product(adxtail, cdytail, tj1, tj0);
1659  Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]);
1660  catt[3] = catt3;
1661  cattlen = 4;
1662  }
1663  else {
1664  cat[0] = 0.0;
1665  catlen = 1;
1666  catt[0] = 0.0;
1667  cattlen = 1;
1668  }
1669 
1670  if (bdxtail != 0.0) {
1671  temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, bdxtail, temp16a);
1672  bxtcatlen = scale_expansion_zeroelim(catlen, cat, bdxtail, bxtcat);
1673  temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, 2.0 * bdx, temp32a);
1674  temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
1675  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
1676  finswap = finnow;
1677  finnow = finother;
1678  finother = finswap;
1679  if (cdytail != 0.0) {
1680  temp8len = scale_expansion_zeroelim(4, aa, bdxtail, temp8);
1681  temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail, temp16a);
1682  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
1683  finswap = finnow;
1684  finnow = finother;
1685  finother = finswap;
1686  }
1687  if (adytail != 0.0) {
1688  temp8len = scale_expansion_zeroelim(4, cc, -bdxtail, temp8);
1689  temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail, temp16a);
1690  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
1691  finswap = finnow;
1692  finnow = finother;
1693  finother = finswap;
1694  }
1695 
1696  temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, bdxtail, temp32a);
1697  bxtcattlen = scale_expansion_zeroelim(cattlen, catt, bdxtail, bxtcatt);
1698  temp16alen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, 2.0 * bdx, temp16a);
1699  temp16blen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, bdxtail, temp16b);
1700  temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
1701  temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
1702  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
1703  finswap = finnow;
1704  finnow = finother;
1705  finother = finswap;
1706  }
1707  if (bdytail != 0.0) {
1708  temp16alen = scale_expansion_zeroelim(bytcalen, bytca, bdytail, temp16a);
1709  bytcatlen = scale_expansion_zeroelim(catlen, cat, bdytail, bytcat);
1710  temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, 2.0 * bdy, temp32a);
1711  temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
1712  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
1713  finswap = finnow;
1714  finnow = finother;
1715  finother = finswap;
1716 
1717  temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, bdytail, temp32a);
1718  bytcattlen = scale_expansion_zeroelim(cattlen, catt, bdytail, bytcatt);
1719  temp16alen = scale_expansion_zeroelim(bytcattlen, bytcatt, 2.0 * bdy, temp16a);
1720  temp16blen = scale_expansion_zeroelim(bytcattlen, bytcatt, bdytail, temp16b);
1721  temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
1722  temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
1723  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
1724  finswap = finnow;
1725  finnow = finother;
1726  finother = finswap;
1727  }
1728  }
1729  if ((cdxtail != 0.0) || (cdytail != 0.0)) {
1730  if ((adxtail != 0.0) || (adytail != 0.0) || (bdxtail != 0.0) || (bdytail != 0.0)) {
1731  Two_Product(adxtail, bdy, ti1, ti0);
1732  Two_Product(adx, bdytail, tj1, tj0);
1733  Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
1734  u[3] = u3;
1735  negate = -ady;
1736  Two_Product(bdxtail, negate, ti1, ti0);
1737  negate = -adytail;
1738  Two_Product(bdx, negate, tj1, tj0);
1739  Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
1740  v[3] = v3;
1741  abtlen = fast_expansion_sum_zeroelim(4, u, 4, v, abt);
1742 
1743  Two_Product(adxtail, bdytail, ti1, ti0);
1744  Two_Product(bdxtail, adytail, tj1, tj0);
1745  Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]);
1746  abtt[3] = abtt3;
1747  abttlen = 4;
1748  }
1749  else {
1750  abt[0] = 0.0;
1751  abtlen = 1;
1752  abtt[0] = 0.0;
1753  abttlen = 1;
1754  }
1755 
1756  if (cdxtail != 0.0) {
1757  temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, cdxtail, temp16a);
1758  cxtabtlen = scale_expansion_zeroelim(abtlen, abt, cdxtail, cxtabt);
1759  temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, 2.0 * cdx, temp32a);
1760  temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
1761  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
1762  finswap = finnow;
1763  finnow = finother;
1764  finother = finswap;
1765  if (adytail != 0.0) {
1766  temp8len = scale_expansion_zeroelim(4, bb, cdxtail, temp8);
1767  temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail, temp16a);
1768  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
1769  finswap = finnow;
1770  finnow = finother;
1771  finother = finswap;
1772  }
1773  if (bdytail != 0.0) {
1774  temp8len = scale_expansion_zeroelim(4, aa, -cdxtail, temp8);
1775  temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail, temp16a);
1776  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
1777  finswap = finnow;
1778  finnow = finother;
1779  finother = finswap;
1780  }
1781 
1782  temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, cdxtail, temp32a);
1783  cxtabttlen = scale_expansion_zeroelim(abttlen, abtt, cdxtail, cxtabtt);
1784  temp16alen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, 2.0 * cdx, temp16a);
1785  temp16blen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, cdxtail, temp16b);
1786  temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
1787  temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
1788  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
1789  finswap = finnow;
1790  finnow = finother;
1791  finother = finswap;
1792  }
1793  if (cdytail != 0.0) {
1794  temp16alen = scale_expansion_zeroelim(cytablen, cytab, cdytail, temp16a);
1795  cytabtlen = scale_expansion_zeroelim(abtlen, abt, cdytail, cytabt);
1796  temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, 2.0 * cdy, temp32a);
1797  temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
1798  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
1799  finswap = finnow;
1800  finnow = finother;
1801  finother = finswap;
1802 
1803  temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, cdytail, temp32a);
1804  cytabttlen = scale_expansion_zeroelim(abttlen, abtt, cdytail, cytabtt);
1805  temp16alen = scale_expansion_zeroelim(cytabttlen, cytabtt, 2.0 * cdy, temp16a);
1806  temp16blen = scale_expansion_zeroelim(cytabttlen, cytabtt, cdytail, temp16b);
1807  temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
1808  temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
1809  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
1810  finswap = finnow;
1811  finnow = finother;
1812  finother = finswap;
1813  }
1814  }
1815 
1816  return finnow[finlength - 1];
1817 }
1818 
1819 double incircle(const double *pa, const double *pb, const double *pc, const double *pd)
1820 {
1821  double adx, bdx, cdx, ady, bdy, cdy;
1822  double bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
1823  double alift, blift, clift;
1824  double det;
1825  double permanent, errbound;
1826 
1827  adx = pa[0] - pd[0];
1828  bdx = pb[0] - pd[0];
1829  cdx = pc[0] - pd[0];
1830  ady = pa[1] - pd[1];
1831  bdy = pb[1] - pd[1];
1832  cdy = pc[1] - pd[1];
1833 
1834  bdxcdy = bdx * cdy;
1835  cdxbdy = cdx * bdy;
1836  alift = adx * adx + ady * ady;
1837 
1838  cdxady = cdx * ady;
1839  adxcdy = adx * cdy;
1840  blift = bdx * bdx + bdy * bdy;
1841 
1842  adxbdy = adx * bdy;
1843  bdxady = bdx * ady;
1844  clift = cdx * cdx + cdy * cdy;
1845 
1846  det = alift * (bdxcdy - cdxbdy) + blift * (cdxady - adxcdy) + clift * (adxbdy - bdxady);
1847 
1848  permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift +
1849  (Absolute(cdxady) + Absolute(adxcdy)) * blift +
1850  (Absolute(adxbdy) + Absolute(bdxady)) * clift;
1851  errbound = iccerrboundA * permanent;
1852  if ((det > errbound) || (-det > errbound)) {
1853  return det;
1854  }
1855 
1856  return incircleadapt(pa, pb, pc, pd, permanent);
1857 }
1858 
1880  const double *pa, const double *pb, const double *pc, const double *pd, const double *pe)
1881 {
1882  double aex, bex, cex, dex;
1883  double aey, bey, cey, dey;
1884  double aez, bez, cez, dez;
1885  double alift, blift, clift, dlift;
1886  double ab, bc, cd, da, ac, bd;
1887  double abc, bcd, cda, dab;
1888 
1889  aex = pa[0] - pe[0];
1890  bex = pb[0] - pe[0];
1891  cex = pc[0] - pe[0];
1892  dex = pd[0] - pe[0];
1893  aey = pa[1] - pe[1];
1894  bey = pb[1] - pe[1];
1895  cey = pc[1] - pe[1];
1896  dey = pd[1] - pe[1];
1897  aez = pa[2] - pe[2];
1898  bez = pb[2] - pe[2];
1899  cez = pc[2] - pe[2];
1900  dez = pd[2] - pe[2];
1901 
1902  ab = aex * bey - bex * aey;
1903  bc = bex * cey - cex * bey;
1904  cd = cex * dey - dex * cey;
1905  da = dex * aey - aex * dey;
1906 
1907  ac = aex * cey - cex * aey;
1908  bd = bex * dey - dex * bey;
1909 
1910  abc = aez * bc - bez * ac + cez * ab;
1911  bcd = bez * cd - cez * bd + dez * bc;
1912  cda = cez * da + dez * ac + aez * cd;
1913  dab = dez * ab + aez * bd + bez * da;
1914 
1915  alift = aex * aex + aey * aey + aez * aez;
1916  blift = bex * bex + bey * bey + bez * bez;
1917  clift = cex * cex + cey * cey + cez * cez;
1918  dlift = dex * dex + dey * dey + dez * dez;
1919 
1920  return (dlift * abc - clift * dab) + (blift * cda - alift * bcd);
1921 }
1922 
1923 static double insphereexact(
1924  const double *pa, const double *pb, const double *pc, const double *pd, const double *pe)
1925 {
1926  INEXACT double axby1, bxcy1, cxdy1, dxey1, exay1;
1927  INEXACT double bxay1, cxby1, dxcy1, exdy1, axey1;
1928  INEXACT double axcy1, bxdy1, cxey1, dxay1, exby1;
1929  INEXACT double cxay1, dxby1, excy1, axdy1, bxey1;
1930  double axby0, bxcy0, cxdy0, dxey0, exay0;
1931  double bxay0, cxby0, dxcy0, exdy0, axey0;
1932  double axcy0, bxdy0, cxey0, dxay0, exby0;
1933  double cxay0, dxby0, excy0, axdy0, bxey0;
1934  double ab[4], bc[4], cd[4], de[4], ea[4];
1935  double ac[4], bd[4], ce[4], da[4], eb[4];
1936  double temp8a[8], temp8b[8], temp16[16];
1937  int temp8alen, temp8blen, temp16len;
1938  double abc[24], bcd[24], cde[24], dea[24], eab[24];
1939  double abd[24], bce[24], cda[24], deb[24], eac[24];
1940  int abclen, bcdlen, cdelen, dealen, eablen;
1941  int abdlen, bcelen, cdalen, deblen, eaclen;
1942  double temp48a[48], temp48b[48];
1943  int temp48alen, temp48blen;
1944  double abcd[96], bcde[96], cdea[96], deab[96], eabc[96];
1945  int abcdlen, bcdelen, cdealen, deablen, eabclen;
1946  double temp192[192];
1947  double det384x[384], det384y[384], det384z[384];
1948  int xlen, ylen, zlen;
1949  double detxy[768];
1950  int xylen;
1951  double adet[1152], bdet[1152], cdet[1152], ddet[1152], edet[1152];
1952  int alen, blen, clen, dlen, elen;
1953  double abdet[2304], cddet[2304], cdedet[3456];
1954  int ablen, cdlen;
1955  double deter[5760];
1956  int deterlen;
1957  int i;
1958 
1959  INEXACT double bvirt;
1960  double avirt, bround, around;
1961  INEXACT double c;
1962  INEXACT double abig;
1963  double ahi, alo, bhi, blo;
1964  double err1, err2, err3;
1965  INEXACT double _i, _j;
1966  double _0;
1967 
1968  Two_Product(pa[0], pb[1], axby1, axby0);
1969  Two_Product(pb[0], pa[1], bxay1, bxay0);
1970  Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]);
1971 
1972  Two_Product(pb[0], pc[1], bxcy1, bxcy0);
1973  Two_Product(pc[0], pb[1], cxby1, cxby0);
1974  Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]);
1975 
1976  Two_Product(pc[0], pd[1], cxdy1, cxdy0);
1977  Two_Product(pd[0], pc[1], dxcy1, dxcy0);
1978  Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]);
1979 
1980  Two_Product(pd[0], pe[1], dxey1, dxey0);
1981  Two_Product(pe[0], pd[1], exdy1, exdy0);
1982  Two_Two_Diff(dxey1, dxey0, exdy1, exdy0, de[3], de[2], de[1], de[0]);
1983 
1984  Two_Product(pe[0], pa[1], exay1, exay0);
1985  Two_Product(pa[0], pe[1], axey1, axey0);
1986  Two_Two_Diff(exay1, exay0, axey1, axey0, ea[3], ea[2], ea[1], ea[0]);
1987 
1988  Two_Product(pa[0], pc[1], axcy1, axcy0);
1989  Two_Product(pc[0], pa[1], cxay1, cxay0);
1990  Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]);
1991 
1992  Two_Product(pb[0], pd[1], bxdy1, bxdy0);
1993  Two_Product(pd[0], pb[1], dxby1, dxby0);
1994  Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]);
1995 
1996  Two_Product(pc[0], pe[1], cxey1, cxey0);
1997  Two_Product(pe[0], pc[1], excy1, excy0);
1998  Two_Two_Diff(cxey1, cxey0, excy1, excy0, ce[3], ce[2], ce[1], ce[0]);
1999 
2000  Two_Product(pd[0], pa[1], dxay1, dxay0);
2001  Two_Product(pa[0], pd[1], axdy1, axdy0);
2002  Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]);
2003 
2004  Two_Product(pe[0], pb[1], exby1, exby0);
2005  Two_Product(pb[0], pe[1], bxey1, bxey0);
2006  Two_Two_Diff(exby1, exby0, bxey1, bxey0, eb[3], eb[2], eb[1], eb[0]);
2007 
2008  temp8alen = scale_expansion_zeroelim(4, bc, pa[2], temp8a);
2009  temp8blen = scale_expansion_zeroelim(4, ac, -pb[2], temp8b);
2010  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, temp16);
2011  temp8alen = scale_expansion_zeroelim(4, ab, pc[2], temp8a);
2012  abclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, abc);
2013 
2014  temp8alen = scale_expansion_zeroelim(4, cd, pb[2], temp8a);
2015  temp8blen = scale_expansion_zeroelim(4, bd, -pc[2], temp8b);
2016  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, temp16);
2017  temp8alen = scale_expansion_zeroelim(4, bc, pd[2], temp8a);
2018  bcdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, bcd);
2019 
2020  temp8alen = scale_expansion_zeroelim(4, de, pc[2], temp8a);
2021  temp8blen = scale_expansion_zeroelim(4, ce, -pd[2], temp8b);
2022  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, temp16);
2023  temp8alen = scale_expansion_zeroelim(4, cd, pe[2], temp8a);
2024  cdelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, cde);
2025 
2026  temp8alen = scale_expansion_zeroelim(4, ea, pd[2], temp8a);
2027  temp8blen = scale_expansion_zeroelim(4, da, -pe[2], temp8b);
2028  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, temp16);
2029  temp8alen = scale_expansion_zeroelim(4, de, pa[2], temp8a);
2030  dealen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, dea);
2031 
2032  temp8alen = scale_expansion_zeroelim(4, ab, pe[2], temp8a);
2033  temp8blen = scale_expansion_zeroelim(4, eb, -pa[2], temp8b);
2034  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, temp16);
2035  temp8alen = scale_expansion_zeroelim(4, ea, pb[2], temp8a);
2036  eablen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, eab);
2037 
2038  temp8alen = scale_expansion_zeroelim(4, bd, pa[2], temp8a);
2039  temp8blen = scale_expansion_zeroelim(4, da, pb[2], temp8b);
2040  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, temp16);
2041  temp8alen = scale_expansion_zeroelim(4, ab, pd[2], temp8a);
2042  abdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, abd);
2043 
2044  temp8alen = scale_expansion_zeroelim(4, ce, pb[2], temp8a);
2045  temp8blen = scale_expansion_zeroelim(4, eb, pc[2], temp8b);
2046  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, temp16);
2047  temp8alen = scale_expansion_zeroelim(4, bc, pe[2], temp8a);
2048  bcelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, bce);
2049 
2050  temp8alen = scale_expansion_zeroelim(4, da, pc[2], temp8a);
2051  temp8blen = scale_expansion_zeroelim(4, ac, pd[2], temp8b);
2052  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, temp16);
2053  temp8alen = scale_expansion_zeroelim(4, cd, pa[2], temp8a);
2054  cdalen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, cda);
2055 
2056  temp8alen = scale_expansion_zeroelim(4, eb, pd[2], temp8a);
2057  temp8blen = scale_expansion_zeroelim(4, bd, pe[2], temp8b);
2058  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, temp16);
2059  temp8alen = scale_expansion_zeroelim(4, de, pb[2], temp8a);
2060  deblen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, deb);
2061 
2062  temp8alen = scale_expansion_zeroelim(4, ac, pe[2], temp8a);
2063  temp8blen = scale_expansion_zeroelim(4, ce, pa[2], temp8b);
2064  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, temp16);
2065  temp8alen = scale_expansion_zeroelim(4, ea, pc[2], temp8a);
2066  eaclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16, eac);
2067 
2068  temp48alen = fast_expansion_sum_zeroelim(cdelen, cde, bcelen, bce, temp48a);
2069  temp48blen = fast_expansion_sum_zeroelim(deblen, deb, bcdlen, bcd, temp48b);
2070  for (i = 0; i < temp48blen; i++) {
2071  temp48b[i] = -temp48b[i];
2072  }
2073  bcdelen = fast_expansion_sum_zeroelim(temp48alen, temp48a, temp48blen, temp48b, bcde);
2074  xlen = scale_expansion_zeroelim(bcdelen, bcde, pa[0], temp192);
2075  xlen = scale_expansion_zeroelim(xlen, temp192, pa[0], det384x);
2076  ylen = scale_expansion_zeroelim(bcdelen, bcde, pa[1], temp192);
2077  ylen = scale_expansion_zeroelim(ylen, temp192, pa[1], det384y);
2078  zlen = scale_expansion_zeroelim(bcdelen, bcde, pa[2], temp192);
2079  zlen = scale_expansion_zeroelim(zlen, temp192, pa[2], det384z);
2080  xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
2081  alen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, adet);
2082 
2083  temp48alen = fast_expansion_sum_zeroelim(dealen, dea, cdalen, cda, temp48a);
2084  temp48blen = fast_expansion_sum_zeroelim(eaclen, eac, cdelen, cde, temp48b);
2085  for (i = 0; i < temp48blen; i++) {
2086  temp48b[i] = -temp48b[i];
2087  }
2088  cdealen = fast_expansion_sum_zeroelim(temp48alen, temp48a, temp48blen, temp48b, cdea);
2089  xlen = scale_expansion_zeroelim(cdealen, cdea, pb[0], temp192);
2090  xlen = scale_expansion_zeroelim(xlen, temp192, pb[0], det384x);
2091  ylen = scale_expansion_zeroelim(cdealen, cdea, pb[1], temp192);
2092  ylen = scale_expansion_zeroelim(ylen, temp192, pb[1], det384y);
2093  zlen = scale_expansion_zeroelim(cdealen, cdea, pb[2], temp192);
2094  zlen = scale_expansion_zeroelim(zlen, temp192, pb[2], det384z);
2095  xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
2096  blen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, bdet);
2097 
2098  temp48alen = fast_expansion_sum_zeroelim(eablen, eab, deblen, deb, temp48a);
2099  temp48blen = fast_expansion_sum_zeroelim(abdlen, abd, dealen, dea, temp48b);
2100  for (i = 0; i < temp48blen; i++) {
2101  temp48b[i] = -temp48b[i];
2102  }
2103  deablen = fast_expansion_sum_zeroelim(temp48alen, temp48a, temp48blen, temp48b, deab);
2104  xlen = scale_expansion_zeroelim(deablen, deab, pc[0], temp192);
2105  xlen = scale_expansion_zeroelim(xlen, temp192, pc[0], det384x);
2106  ylen = scale_expansion_zeroelim(deablen, deab, pc[1], temp192);
2107  ylen = scale_expansion_zeroelim(ylen, temp192, pc[1], det384y);
2108  zlen = scale_expansion_zeroelim(deablen, deab, pc[2], temp192);
2109  zlen = scale_expansion_zeroelim(zlen, temp192, pc[2], det384z);
2110  xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
2111  clen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, cdet);
2112 
2113  temp48alen = fast_expansion_sum_zeroelim(abclen, abc, eaclen, eac, temp48a);
2114  temp48blen = fast_expansion_sum_zeroelim(bcelen, bce, eablen, eab, temp48b);
2115  for (i = 0; i < temp48blen; i++) {
2116  temp48b[i] = -temp48b[i];
2117  }
2118  eabclen = fast_expansion_sum_zeroelim(temp48alen, temp48a, temp48blen, temp48b, eabc);
2119  xlen = scale_expansion_zeroelim(eabclen, eabc, pd[0], temp192);
2120  xlen = scale_expansion_zeroelim(xlen, temp192, pd[0], det384x);
2121  ylen = scale_expansion_zeroelim(eabclen, eabc, pd[1], temp192);
2122  ylen = scale_expansion_zeroelim(ylen, temp192, pd[1], det384y);
2123  zlen = scale_expansion_zeroelim(eabclen, eabc, pd[2], temp192);
2124  zlen = scale_expansion_zeroelim(zlen, temp192, pd[2], det384z);
2125  xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
2126  dlen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, ddet);
2127 
2128  temp48alen = fast_expansion_sum_zeroelim(bcdlen, bcd, abdlen, abd, temp48a);
2129  temp48blen = fast_expansion_sum_zeroelim(cdalen, cda, abclen, abc, temp48b);
2130  for (i = 0; i < temp48blen; i++) {
2131  temp48b[i] = -temp48b[i];
2132  }
2133  abcdlen = fast_expansion_sum_zeroelim(temp48alen, temp48a, temp48blen, temp48b, abcd);
2134  xlen = scale_expansion_zeroelim(abcdlen, abcd, pe[0], temp192);
2135  xlen = scale_expansion_zeroelim(xlen, temp192, pe[0], det384x);
2136  ylen = scale_expansion_zeroelim(abcdlen, abcd, pe[1], temp192);
2137  ylen = scale_expansion_zeroelim(ylen, temp192, pe[1], det384y);
2138  zlen = scale_expansion_zeroelim(abcdlen, abcd, pe[2], temp192);
2139  zlen = scale_expansion_zeroelim(zlen, temp192, pe[2], det384z);
2140  xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
2141  elen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, edet);
2142 
2143  ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
2144  cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
2145  cdelen = fast_expansion_sum_zeroelim(cdlen, cddet, elen, edet, cdedet);
2146  deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdelen, cdedet, deter);
2147 
2148  return deter[deterlen - 1];
2149 }
2150 
2151 static double insphereadapt(const double *pa,
2152  const double *pb,
2153  const double *pc,
2154  const double *pd,
2155  const double *pe,
2156  double permanent)
2157 {
2158  INEXACT double aex, bex, cex, dex, aey, bey, cey, dey, aez, bez, cez, dez;
2159  double det, errbound;
2160 
2161  INEXACT double aexbey1, bexaey1, bexcey1, cexbey1;
2162  INEXACT double cexdey1, dexcey1, dexaey1, aexdey1;
2163  INEXACT double aexcey1, cexaey1, bexdey1, dexbey1;
2164  double aexbey0, bexaey0, bexcey0, cexbey0;
2165  double cexdey0, dexcey0, dexaey0, aexdey0;
2166  double aexcey0, cexaey0, bexdey0, dexbey0;
2167  double ab[4], bc[4], cd[4], da[4], ac[4], bd[4];
2168  INEXACT double ab3, bc3, cd3, da3, ac3, bd3;
2169  double abeps, bceps, cdeps, daeps, aceps, bdeps;
2170  double temp8a[8], temp8b[8], temp8c[8], temp16[16], temp24[24], temp48[48];
2171  int temp8alen, temp8blen, temp8clen, temp16len, temp24len, temp48len;
2172  double xdet[96], ydet[96], zdet[96], xydet[192];
2173  int xlen, ylen, zlen, xylen;
2174  double adet[288], bdet[288], cdet[288], ddet[288];
2175  int alen, blen, clen, dlen;
2176  double abdet[576], cddet[576];
2177  int ablen, cdlen;
2178  double fin1[1152];
2179  int finlength;
2180 
2181  double aextail, bextail, cextail, dextail;
2182  double aeytail, beytail, ceytail, deytail;
2183  double aeztail, beztail, ceztail, deztail;
2184 
2185  INEXACT double bvirt;
2186  double avirt, bround, around;
2187  INEXACT double c;
2188  INEXACT double abig;
2189  double ahi, alo, bhi, blo;
2190  double err1, err2, err3;
2191  INEXACT double _i, _j;
2192  double _0;
2193 
2194  aex = (double)(pa[0] - pe[0]);
2195  bex = (double)(pb[0] - pe[0]);
2196  cex = (double)(pc[0] - pe[0]);
2197  dex = (double)(pd[0] - pe[0]);
2198  aey = (double)(pa[1] - pe[1]);
2199  bey = (double)(pb[1] - pe[1]);
2200  cey = (double)(pc[1] - pe[1]);
2201  dey = (double)(pd[1] - pe[1]);
2202  aez = (double)(pa[2] - pe[2]);
2203  bez = (double)(pb[2] - pe[2]);
2204  cez = (double)(pc[2] - pe[2]);
2205  dez = (double)(pd[2] - pe[2]);
2206 
2207  Two_Product(aex, bey, aexbey1, aexbey0);
2208  Two_Product(bex, aey, bexaey1, bexaey0);
2209  Two_Two_Diff(aexbey1, aexbey0, bexaey1, bexaey0, ab3, ab[2], ab[1], ab[0]);
2210  ab[3] = ab3;
2211 
2212  Two_Product(bex, cey, bexcey1, bexcey0);
2213  Two_Product(cex, bey, cexbey1, cexbey0);
2214  Two_Two_Diff(bexcey1, bexcey0, cexbey1, cexbey0, bc3, bc[2], bc[1], bc[0]);
2215  bc[3] = bc3;
2216 
2217  Two_Product(cex, dey, cexdey1, cexdey0);
2218  Two_Product(dex, cey, dexcey1, dexcey0);
2219  Two_Two_Diff(cexdey1, cexdey0, dexcey1, dexcey0, cd3, cd[2], cd[1], cd[0]);
2220  cd[3] = cd3;
2221 
2222  Two_Product(dex, aey, dexaey1, dexaey0);
2223  Two_Product(aex, dey, aexdey1, aexdey0);
2224  Two_Two_Diff(dexaey1, dexaey0, aexdey1, aexdey0, da3, da[2], da[1], da[0]);
2225  da[3] = da3;
2226 
2227  Two_Product(aex, cey, aexcey1, aexcey0);
2228  Two_Product(cex, aey, cexaey1, cexaey0);
2229  Two_Two_Diff(aexcey1, aexcey0, cexaey1, cexaey0, ac3, ac[2], ac[1], ac[0]);
2230  ac[3] = ac3;
2231 
2232  Two_Product(bex, dey, bexdey1, bexdey0);
2233  Two_Product(dex, bey, dexbey1, dexbey0);
2234  Two_Two_Diff(bexdey1, bexdey0, dexbey1, dexbey0, bd3, bd[2], bd[1], bd[0]);
2235  bd[3] = bd3;
2236 
2237  temp8alen = scale_expansion_zeroelim(4, cd, bez, temp8a);
2238  temp8blen = scale_expansion_zeroelim(4, bd, -cez, temp8b);
2239  temp8clen = scale_expansion_zeroelim(4, bc, dez, temp8c);
2240  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, temp16);
2241  temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, temp16len, temp16, temp24);
2242  temp48len = scale_expansion_zeroelim(temp24len, temp24, aex, temp48);
2243  xlen = scale_expansion_zeroelim(temp48len, temp48, -aex, xdet);
2244  temp48len = scale_expansion_zeroelim(temp24len, temp24, aey, temp48);
2245  ylen = scale_expansion_zeroelim(temp48len, temp48, -aey, ydet);
2246  temp48len = scale_expansion_zeroelim(temp24len, temp24, aez, temp48);
2247  zlen = scale_expansion_zeroelim(temp48len, temp48, -aez, zdet);
2248  xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
2249  alen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, adet);
2250 
2251  temp8alen = scale_expansion_zeroelim(4, da, cez, temp8a);
2252  temp8blen = scale_expansion_zeroelim(4, ac, dez, temp8b);
2253  temp8clen = scale_expansion_zeroelim(4, cd, aez, temp8c);
2254  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, temp16);
2255  temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, temp16len, temp16, temp24);
2256  temp48len = scale_expansion_zeroelim(temp24len, temp24, bex, temp48);
2257  xlen = scale_expansion_zeroelim(temp48len, temp48, bex, xdet);
2258  temp48len = scale_expansion_zeroelim(temp24len, temp24, bey, temp48);
2259  ylen = scale_expansion_zeroelim(temp48len, temp48, bey, ydet);
2260  temp48len = scale_expansion_zeroelim(temp24len, temp24, bez, temp48);
2261  zlen = scale_expansion_zeroelim(temp48len, temp48, bez, zdet);
2262  xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
2263  blen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, bdet);
2264 
2265  temp8alen = scale_expansion_zeroelim(4, ab, dez, temp8a);
2266  temp8blen = scale_expansion_zeroelim(4, bd, aez, temp8b);
2267  temp8clen = scale_expansion_zeroelim(4, da, bez, temp8c);
2268  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, temp16);
2269  temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, temp16len, temp16, temp24);
2270  temp48len = scale_expansion_zeroelim(temp24len, temp24, cex, temp48);
2271  xlen = scale_expansion_zeroelim(temp48len, temp48, -cex, xdet);
2272  temp48len = scale_expansion_zeroelim(temp24len, temp24, cey, temp48);
2273  ylen = scale_expansion_zeroelim(temp48len, temp48, -cey, ydet);
2274  temp48len = scale_expansion_zeroelim(temp24len, temp24, cez, temp48);
2275  zlen = scale_expansion_zeroelim(temp48len, temp48, -cez, zdet);
2276  xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
2277  clen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, cdet);
2278 
2279  temp8alen = scale_expansion_zeroelim(4, bc, aez, temp8a);
2280  temp8blen = scale_expansion_zeroelim(4, ac, -bez, temp8b);
2281  temp8clen = scale_expansion_zeroelim(4, ab, cez, temp8c);
2282  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b, temp16);
2283  temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c, temp16len, temp16, temp24);
2284  temp48len = scale_expansion_zeroelim(temp24len, temp24, dex, temp48);
2285  xlen = scale_expansion_zeroelim(temp48len, temp48, dex, xdet);
2286  temp48len = scale_expansion_zeroelim(temp24len, temp24, dey, temp48);
2287  ylen = scale_expansion_zeroelim(temp48len, temp48, dey, ydet);
2288  temp48len = scale_expansion_zeroelim(temp24len, temp24, dez, temp48);
2289  zlen = scale_expansion_zeroelim(temp48len, temp48, dez, zdet);
2290  xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
2291  dlen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, ddet);
2292 
2293  ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
2294  cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
2295  finlength = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, fin1);
2296 
2297  det = estimate(finlength, fin1);
2298  errbound = isperrboundB * permanent;
2299  if ((det >= errbound) || (-det >= errbound)) {
2300  return det;
2301  }
2302 
2303  Two_Diff_Tail(pa[0], pe[0], aex, aextail);
2304  Two_Diff_Tail(pa[1], pe[1], aey, aeytail);
2305  Two_Diff_Tail(pa[2], pe[2], aez, aeztail);
2306  Two_Diff_Tail(pb[0], pe[0], bex, bextail);
2307  Two_Diff_Tail(pb[1], pe[1], bey, beytail);
2308  Two_Diff_Tail(pb[2], pe[2], bez, beztail);
2309  Two_Diff_Tail(pc[0], pe[0], cex, cextail);
2310  Two_Diff_Tail(pc[1], pe[1], cey, ceytail);
2311  Two_Diff_Tail(pc[2], pe[2], cez, ceztail);
2312  Two_Diff_Tail(pd[0], pe[0], dex, dextail);
2313  Two_Diff_Tail(pd[1], pe[1], dey, deytail);
2314  Two_Diff_Tail(pd[2], pe[2], dez, deztail);
2315  if ((aextail == 0.0) && (aeytail == 0.0) && (aeztail == 0.0) && (bextail == 0.0) &&
2316  (beytail == 0.0) && (beztail == 0.0) && (cextail == 0.0) && (ceytail == 0.0) &&
2317  (ceztail == 0.0) && (dextail == 0.0) && (deytail == 0.0) && (deztail == 0.0)) {
2318  return det;
2319  }
2320 
2321  errbound = isperrboundC * permanent + resulterrbound * Absolute(det);
2322  abeps = (aex * beytail + bey * aextail) - (aey * bextail + bex * aeytail);
2323  bceps = (bex * ceytail + cey * bextail) - (bey * cextail + cex * beytail);
2324  cdeps = (cex * deytail + dey * cextail) - (cey * dextail + dex * ceytail);
2325  daeps = (dex * aeytail + aey * dextail) - (dey * aextail + aex * deytail);
2326  aceps = (aex * ceytail + cey * aextail) - (aey * cextail + cex * aeytail);
2327  bdeps = (bex * deytail + dey * bextail) - (bey * dextail + dex * beytail);
2328  det +=
2329  (((bex * bex + bey * bey + bez * bez) * ((cez * daeps + dez * aceps + aez * cdeps) +
2330  (ceztail * da3 + deztail * ac3 + aeztail * cd3)) +
2331  (dex * dex + dey * dey + dez * dez) * ((aez * bceps - bez * aceps + cez * abeps) +
2332  (aeztail * bc3 - beztail * ac3 + ceztail * ab3))) -
2333  ((aex * aex + aey * aey + aez * aez) * ((bez * cdeps - cez * bdeps + dez * bceps) +
2334  (beztail * cd3 - ceztail * bd3 + deztail * bc3)) +
2335  (cex * cex + cey * cey + cez * cez) * ((dez * abeps + aez * bdeps + bez * daeps) +
2336  (deztail * ab3 + aeztail * bd3 + beztail * da3)))) +
2337  2.0 *
2338  (((bex * bextail + bey * beytail + bez * beztail) * (cez * da3 + dez * ac3 + aez * cd3) +
2339  (dex * dextail + dey * deytail + dez * deztail) *
2340  (aez * bc3 - bez * ac3 + cez * ab3)) -
2341  ((aex * aextail + aey * aeytail + aez * aeztail) * (bez * cd3 - cez * bd3 + dez * bc3) +
2342  (cex * cextail + cey * ceytail + cez * ceztail) *
2343  (dez * ab3 + aez * bd3 + bez * da3)));
2344  if ((det >= errbound) || (-det >= errbound)) {
2345  return det;
2346  }
2347 
2348  return insphereexact(pa, pb, pc, pd, pe);
2349 }
2350 
2351 double insphere(
2352  const double *pa, const double *pb, const double *pc, const double *pd, const double *pe)
2353 {
2354  double aex, bex, cex, dex;
2355  double aey, bey, cey, dey;
2356  double aez, bez, cez, dez;
2357  double aexbey, bexaey, bexcey, cexbey, cexdey, dexcey, dexaey, aexdey;
2358  double aexcey, cexaey, bexdey, dexbey;
2359  double alift, blift, clift, dlift;
2360  double ab, bc, cd, da, ac, bd;
2361  double abc, bcd, cda, dab;
2362  double aezplus, bezplus, cezplus, dezplus;
2363  double aexbeyplus, bexaeyplus, bexceyplus, cexbeyplus;
2364  double cexdeyplus, dexceyplus, dexaeyplus, aexdeyplus;
2365  double aexceyplus, cexaeyplus, bexdeyplus, dexbeyplus;
2366  double det;
2367  double permanent, errbound;
2368 
2369  aex = pa[0] - pe[0];
2370  bex = pb[0] - pe[0];
2371  cex = pc[0] - pe[0];
2372  dex = pd[0] - pe[0];
2373  aey = pa[1] - pe[1];
2374  bey = pb[1] - pe[1];
2375  cey = pc[1] - pe[1];
2376  dey = pd[1] - pe[1];
2377  aez = pa[2] - pe[2];
2378  bez = pb[2] - pe[2];
2379  cez = pc[2] - pe[2];
2380  dez = pd[2] - pe[2];
2381 
2382  aexbey = aex * bey;
2383  bexaey = bex * aey;
2384  ab = aexbey - bexaey;
2385  bexcey = bex * cey;
2386  cexbey = cex * bey;
2387  bc = bexcey - cexbey;
2388  cexdey = cex * dey;
2389  dexcey = dex * cey;
2390  cd = cexdey - dexcey;
2391  dexaey = dex * aey;
2392  aexdey = aex * dey;
2393  da = dexaey - aexdey;
2394 
2395  aexcey = aex * cey;
2396  cexaey = cex * aey;
2397  ac = aexcey - cexaey;
2398  bexdey = bex * dey;
2399  dexbey = dex * bey;
2400  bd = bexdey - dexbey;
2401 
2402  abc = aez * bc - bez * ac + cez * ab;
2403  bcd = bez * cd - cez * bd + dez * bc;
2404  cda = cez * da + dez * ac + aez * cd;
2405  dab = dez * ab + aez * bd + bez * da;
2406 
2407  alift = aex * aex + aey * aey + aez * aez;
2408  blift = bex * bex + bey * bey + bez * bez;
2409  clift = cex * cex + cey * cey + cez * cez;
2410  dlift = dex * dex + dey * dey + dez * dez;
2411 
2412  det = (dlift * abc - clift * dab) + (blift * cda - alift * bcd);
2413 
2414  aezplus = Absolute(aez);
2415  bezplus = Absolute(bez);
2416  cezplus = Absolute(cez);
2417  dezplus = Absolute(dez);
2418  aexbeyplus = Absolute(aexbey);
2419  bexaeyplus = Absolute(bexaey);
2420  bexceyplus = Absolute(bexcey);
2421  cexbeyplus = Absolute(cexbey);
2422  cexdeyplus = Absolute(cexdey);
2423  dexceyplus = Absolute(dexcey);
2424  dexaeyplus = Absolute(dexaey);
2425  aexdeyplus = Absolute(aexdey);
2426  aexceyplus = Absolute(aexcey);
2427  cexaeyplus = Absolute(cexaey);
2428  bexdeyplus = Absolute(bexdey);
2429  dexbeyplus = Absolute(dexbey);
2430  permanent = ((cexdeyplus + dexceyplus) * bezplus + (dexbeyplus + bexdeyplus) * cezplus +
2431  (bexceyplus + cexbeyplus) * dezplus) *
2432  alift +
2433  ((dexaeyplus + aexdeyplus) * cezplus + (aexceyplus + cexaeyplus) * dezplus +
2434  (cexdeyplus + dexceyplus) * aezplus) *
2435  blift +
2436  ((aexbeyplus + bexaeyplus) * dezplus + (bexdeyplus + dexbeyplus) * aezplus +
2437  (dexaeyplus + aexdeyplus) * bezplus) *
2438  clift +
2439  ((bexceyplus + cexbeyplus) * aezplus + (cexaeyplus + aexceyplus) * bezplus +
2440  (aexbeyplus + bexaeyplus) * cezplus) *
2441  dlift;
2442  errbound = isperrboundA * permanent;
2443  if ((det > errbound) || (-det > errbound)) {
2444  return det;
2445  }
2446 
2447  return insphereadapt(pa, pb, pc, pd, pe, permanent);
2448 }
2449 
2450 } /* namespace robust_pred */
2451 
2452 static int sgn(double x)
2453 {
2454  return (x > 0) ? 1 : ((x < 0) ? -1 : 0);
2455 }
2456 
2457 int orient2d(const double2 &a, const double2 &b, const double2 &c)
2458 {
2459  return sgn(blender::robust_pred::orient2d(a, b, c));
2460 }
2461 
2462 int orient2d_fast(const double2 &a, const double2 &b, const double2 &c)
2463 {
2465 }
2466 
2467 int incircle(const double2 &a, const double2 &b, const double2 &c, const double2 &d)
2468 {
2469  return sgn(robust_pred::incircle(a, b, c, d));
2470 }
2471 
2472 int incircle_fast(const double2 &a, const double2 &b, const double2 &c, const double2 &d)
2473 {
2474  return sgn(robust_pred::incirclefast(a, b, c, d));
2475 }
2476 
2477 int orient3d(const double3 &a, const double3 &b, const double3 &c, const double3 &d)
2478 {
2479  return sgn(robust_pred::orient3d(a, b, c, d));
2480 }
2481 
2482 int orient3d_fast(const double3 &a, const double3 &b, const double3 &c, const double3 &d)
2483 {
2484  return sgn(robust_pred::orient3dfast(a, b, c, d));
2485 }
2486 
2488  const double3 &a, const double3 &b, const double3 &c, const double3 &d, const double3 &e)
2489 {
2490  return sgn(robust_pred::insphere(a, b, c, d, e));
2491 }
2492 
2494  const double3 &a, const double3 &b, const double3 &c, const double3 &d, const double3 &e)
2495 {
2496  return sgn(robust_pred::inspherefast(a, b, c, d, e));
2497 }
2498 
2499 } // namespace blender
#define D
Math vector functions needed specifically for mesh intersect and boolean.
#define ELEM(...)
typedef double(DMatrix)[4][4]
ATTR_WARN_UNUSED_RESULT const BMVert const BMEdge * e
ATTR_WARN_UNUSED_RESULT const BMVert * v
SIMD_FORCE_INLINE const btScalar & w() const
Return the w value.
Definition: btQuadWord.h:119
static T sum(const btAlignedObjectArray< T > &items)
Definition: half.h:41
unsigned short half
Definition: cuda/compat.h:110
#define Square(a, x, y)
#define Two_Product(a, b, x, y)
#define Two_Product_Presplit(a, b, bhi, blo, x, y)
#define Two_Diff_Tail(a, b, x, y)
#define Two_Sum(a, b, x, y)
#define Two_One_Product(a1, a0, b, x3, x2, x1, x0)
#define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0)
#define INEXACT
#define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0)
#define Fast_Two_Sum(a, b, x, y)
#define Split(a, ahi, alo)
#define Absolute(a)
#define B
static unsigned c
Definition: RandGen.cpp:83
static double B3(double u)
Definition: FitCurve.cpp:316
static unsigned a[3]
Definition: RandGen.cpp:78
static int fast_expansion_sum_zeroelim(int elen, const double *e, int flen, const double *f, double *h)
static double iccerrboundA
double orient2dfast(const double *pa, const double *pb, const double *pc)
static double resulterrbound
static double ccwerrboundB
static double orient2dadapt(const double *pa, const double *pb, const double *pc, double detsum)
static double o3derrboundC
double incirclefast(const double *pa, const double *pb, const double *pc, const double *pd)
static double isperrboundB
static double o3derrboundB
double insphere(const double *pa, const double *pb, const double *pc, const double *pd, const double *pe)
static double isperrboundC
static double epsilon
static double splitter
static double ccwerrboundC
static double iccerrboundB
static double estimate(int elen, const double *e)
double orient3dfast(const double *pa, const double *pb, const double *pc, const double *pd)
static double isperrboundA
static RobustInitCaller init_caller
static double incircleadapt(const double *pa, const double *pb, const double *pc, const double *pd, double permanent)
static double insphereadapt(const double *pa, const double *pb, const double *pc, const double *pd, const double *pe, double permanent)
static double o3derrboundA
static double orient3dadapt(const double *pa, const double *pb, const double *pc, const double *pd, double permanent)
double orient3d(const double *pa, const double *pb, const double *pc, const double *pd)
static double iccerrboundC
double incircle(const double *pa, const double *pb, const double *pc, const double *pd)
static int scale_expansion_zeroelim(int elen, const double *e, double b, double *h)
double inspherefast(const double *pa, const double *pb, const double *pc, const double *pd, const double *pe)
double orient2d(const double *pa, const double *pb, const double *pc)
static double ccwerrboundA
static double insphereexact(const double *pa, const double *pb, const double *pc, const double *pd, const double *pe)
int orient3d_fast(const double3 &a, const double3 &b, const double3 &c, const double3 &d)
int orient2d(const double2 &a, const double2 &b, const double2 &c)
static int sgn(double x)
int incircle(const double2 &a, const double2 &b, const double2 &c, const double2 &d)
int orient3d(const double3 &a, const double3 &b, const double3 &c, const double3 &d)
int incircle_fast(const double2 &a, const double2 &b, const double2 &c, const double2 &d)
int orient2d_fast(const double2 &a, const double2 &b, const double2 &c)
int insphere_fast(const double3 &a, const double3 &b, const double3 &c, const double3 &d, const double3 &e)
int insphere(const double3 &a, const double3 &b, const double3 &c, const double3 &d, const double3 &e)
static const pxr::TfToken b("b", pxr::TfToken::Immortal)