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