Chaste  Release::2017.1
predicates.cpp
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 #include <stdio.h>
117 #include <stdlib.h>
118 #include <math.h>
119 #ifdef CPU86
120 #include <float.h>
121 #endif /* CPU86 */
122 #ifdef LINUX
123 #include <fpu_control.h>
124 #endif /* LINUX */
125 
126 #include "tetgen.h" // Defines the symbol REAL (float or double).
127 
128 /* On some machines, the exact arithmetic routines might be defeated by the */
129 /* use of internal extended precision floating-point registers. Sometimes */
130 /* this problem can be fixed by defining certain values to be volatile, */
131 /* thus forcing them to be stored to memory and rounded off. This isn't */
132 /* a great solution, though, as it slows the arithmetic down. */
133 /* */
134 /* To try this out, write "#define INEXACT volatile" below. Normally, */
135 /* however, INEXACT should be defined to be nothing. ("#define INEXACT".) */
136 
137 #define INEXACT /* Nothing */
138 /* #define INEXACT volatile */
139 
140 /* #define REAL double */ /* float or double */
141 #define REALPRINT doubleprint
142 #define REALRAND doublerand
143 #define NARROWRAND narrowdoublerand
144 #define UNIFORMRAND uniformdoublerand
145 
146 /* Which of the following two methods of finding the absolute values is */
147 /* fastest is compiler-dependent. A few compilers can inline and optimize */
148 /* the fabs() call; but most will incur the overhead of a function call, */
149 /* which is disastrously slow. A faster way on IEEE machines might be to */
150 /* mask the appropriate bit, but that's difficult to do in C. */
151 
152 #define Absolute(a) ((a) >= 0.0 ? (a) : -(a))
153 /* #define Absolute(a) fabs(a) */
154 
155 /* Many of the operations are broken up into two pieces, a main part that */
156 /* performs an approximate operation, and a "tail" that computes the */
157 /* roundoff error of that operation. */
158 /* */
159 /* The operations Fast_Two_Sum(), Fast_Two_Diff(), Two_Sum(), Two_Diff(), */
160 /* Split(), and Two_Product() are all implemented as described in the */
161 /* reference. Each of these macros requires certain variables to be */
162 /* defined in the calling routine. The variables `bvirt', `c', `abig', */
163 /* `_i', `_j', `_k', `_l', `_m', and `_n' are declared `INEXACT' because */
164 /* they store the result of an operation that may incur roundoff error. */
165 /* The input parameter `x' (or the highest numbered `x_' parameter) must */
166 /* also be declared `INEXACT'. */
167 
168 #define Fast_Two_Sum_Tail(a, b, x, y) \
169  bvirt = x - a; \
170  y = b - bvirt
171 
172 #define Fast_Two_Sum(a, b, x, y) \
173  x = (REAL) (a + b); \
174  Fast_Two_Sum_Tail(a, b, x, y)
175 
176 #define Fast_Two_Diff_Tail(a, b, x, y) \
177  bvirt = a - x; \
178  y = bvirt - b
179 
180 #define Fast_Two_Diff(a, b, x, y) \
181  x = (REAL) (a - b); \
182  Fast_Two_Diff_Tail(a, b, x, y)
183 
184 #define Two_Sum_Tail(a, b, x, y) \
185  bvirt = (REAL) (x - a); \
186  avirt = x - bvirt; \
187  bround = b - bvirt; \
188  around = a - avirt; \
189  y = around + bround
190 
191 #define Two_Sum(a, b, x, y) \
192  x = (REAL) (a + b); \
193  Two_Sum_Tail(a, b, x, y)
194 
195 #define Two_Diff_Tail(a, b, x, y) \
196  bvirt = (REAL) (a - x); \
197  avirt = x + bvirt; \
198  bround = bvirt - b; \
199  around = a - avirt; \
200  y = around + bround
201 
202 #define Two_Diff(a, b, x, y) \
203  x = (REAL) (a - b); \
204  Two_Diff_Tail(a, b, x, y)
205 
206 #define Split(a, ahi, alo) \
207  c = (REAL) (splitter * a); \
208  abig = (REAL) (c - a); \
209  ahi = c - abig; \
210  alo = a - ahi
211 
212 #define Two_Product_Tail(a, b, x, y) \
213  Split(a, ahi, alo); \
214  Split(b, bhi, blo); \
215  err1 = x - (ahi * bhi); \
216  err2 = err1 - (alo * bhi); \
217  err3 = err2 - (ahi * blo); \
218  y = (alo * blo) - err3
219 
220 #define Two_Product(a, b, x, y) \
221  x = (REAL) (a * b); \
222  Two_Product_Tail(a, b, x, y)
223 
224 /* Two_Product_Presplit() is Two_Product() where one of the inputs has */
225 /* already been split. Avoids redundant splitting. */
226 
227 #define Two_Product_Presplit(a, b, bhi, blo, x, y) \
228  x = (REAL) (a * b); \
229  Split(a, ahi, alo); \
230  err1 = x - (ahi * bhi); \
231  err2 = err1 - (alo * bhi); \
232  err3 = err2 - (ahi * blo); \
233  y = (alo * blo) - err3
234 
235 /* Two_Product_2Presplit() is Two_Product() where both of the inputs have */
236 /* already been split. Avoids redundant splitting. */
237 
238 #define Two_Product_2Presplit(a, ahi, alo, b, bhi, blo, x, y) \
239  x = (REAL) (a * b); \
240  err1 = x - (ahi * bhi); \
241  err2 = err1 - (alo * bhi); \
242  err3 = err2 - (ahi * blo); \
243  y = (alo * blo) - err3
244 
245 /* Square() can be done more quickly than Two_Product(). */
246 
247 #define Square_Tail(a, x, y) \
248  Split(a, ahi, alo); \
249  err1 = x - (ahi * ahi); \
250  err3 = err1 - ((ahi + ahi) * alo); \
251  y = (alo * alo) - err3
252 
253 #define Square(a, x, y) \
254  x = (REAL) (a * a); \
255  Square_Tail(a, x, y)
256 
257 /* Macros for summing expansions of various fixed lengths. These are all */
258 /* unrolled versions of Expansion_Sum(). */
259 
260 #define Two_One_Sum(a1, a0, b, x2, x1, x0) \
261  Two_Sum(a0, b , _i, x0); \
262  Two_Sum(a1, _i, x2, x1)
263 
264 #define Two_One_Diff(a1, a0, b, x2, x1, x0) \
265  Two_Diff(a0, b , _i, x0); \
266  Two_Sum( a1, _i, x2, x1)
267 
268 #define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \
269  Two_One_Sum(a1, a0, b0, _j, _0, x0); \
270  Two_One_Sum(_j, _0, b1, x3, x2, x1)
271 
272 #define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \
273  Two_One_Diff(a1, a0, b0, _j, _0, x0); \
274  Two_One_Diff(_j, _0, b1, x3, x2, x1)
275 
276 #define Four_One_Sum(a3, a2, a1, a0, b, x4, x3, x2, x1, x0) \
277  Two_One_Sum(a1, a0, b , _j, x1, x0); \
278  Two_One_Sum(a3, a2, _j, x4, x3, x2)
279 
280 #define Four_Two_Sum(a3, a2, a1, a0, b1, b0, x5, x4, x3, x2, x1, x0) \
281  Four_One_Sum(a3, a2, a1, a0, b0, _k, _2, _1, _0, x0); \
282  Four_One_Sum(_k, _2, _1, _0, b1, x5, x4, x3, x2, x1)
283 
284 #define Four_Four_Sum(a3, a2, a1, a0, b4, b3, b1, b0, x7, x6, x5, x4, x3, x2, \
285  x1, x0) \
286  Four_Two_Sum(a3, a2, a1, a0, b1, b0, _l, _2, _1, _0, x1, x0); \
287  Four_Two_Sum(_l, _2, _1, _0, b4, b3, x7, x6, x5, x4, x3, x2)
288 
289 #define Eight_One_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b, x8, x7, x6, x5, x4, \
290  x3, x2, x1, x0) \
291  Four_One_Sum(a3, a2, a1, a0, b , _j, x3, x2, x1, x0); \
292  Four_One_Sum(a7, a6, a5, a4, _j, x8, x7, x6, x5, x4)
293 
294 #define Eight_Two_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b1, b0, x9, x8, x7, \
295  x6, x5, x4, x3, x2, x1, x0) \
296  Eight_One_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b0, _k, _6, _5, _4, _3, _2, \
297  _1, _0, x0); \
298  Eight_One_Sum(_k, _6, _5, _4, _3, _2, _1, _0, b1, x9, x8, x7, x6, x5, x4, \
299  x3, x2, x1)
300 
301 #define Eight_Four_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b4, b3, b1, b0, x11, \
302  x10, x9, x8, x7, x6, x5, x4, x3, x2, x1, x0) \
303  Eight_Two_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b1, b0, _l, _6, _5, _4, _3, \
304  _2, _1, _0, x1, x0); \
305  Eight_Two_Sum(_l, _6, _5, _4, _3, _2, _1, _0, b4, b3, x11, x10, x9, x8, \
306  x7, x6, x5, x4, x3, x2)
307 
308 /* Macros for multiplying expansions of various fixed lengths. */
309 
310 #define Two_One_Product(a1, a0, b, x3, x2, x1, x0) \
311  Split(b, bhi, blo); \
312  Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \
313  Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
314  Two_Sum(_i, _0, _k, x1); \
315  Fast_Two_Sum(_j, _k, x3, x2)
316 
317 #define Four_One_Product(a3, a2, a1, a0, b, x7, x6, x5, x4, x3, x2, x1, x0) \
318  Split(b, bhi, blo); \
319  Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \
320  Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
321  Two_Sum(_i, _0, _k, x1); \
322  Fast_Two_Sum(_j, _k, _i, x2); \
323  Two_Product_Presplit(a2, b, bhi, blo, _j, _0); \
324  Two_Sum(_i, _0, _k, x3); \
325  Fast_Two_Sum(_j, _k, _i, x4); \
326  Two_Product_Presplit(a3, b, bhi, blo, _j, _0); \
327  Two_Sum(_i, _0, _k, x5); \
328  Fast_Two_Sum(_j, _k, x7, x6)
329 
330 #define Two_Two_Product(a1, a0, b1, b0, x7, x6, x5, x4, x3, x2, x1, x0) \
331  Split(a0, a0hi, a0lo); \
332  Split(b0, bhi, blo); \
333  Two_Product_2Presplit(a0, a0hi, a0lo, b0, bhi, blo, _i, x0); \
334  Split(a1, a1hi, a1lo); \
335  Two_Product_2Presplit(a1, a1hi, a1lo, b0, bhi, blo, _j, _0); \
336  Two_Sum(_i, _0, _k, _1); \
337  Fast_Two_Sum(_j, _k, _l, _2); \
338  Split(b1, bhi, blo); \
339  Two_Product_2Presplit(a0, a0hi, a0lo, b1, bhi, blo, _i, _0); \
340  Two_Sum(_1, _0, _k, x1); \
341  Two_Sum(_2, _k, _j, _1); \
342  Two_Sum(_l, _j, _m, _2); \
343  Two_Product_2Presplit(a1, a1hi, a1lo, b1, bhi, blo, _j, _0); \
344  Two_Sum(_i, _0, _n, _0); \
345  Two_Sum(_1, _0, _i, x2); \
346  Two_Sum(_2, _i, _k, _1); \
347  Two_Sum(_m, _k, _l, _2); \
348  Two_Sum(_j, _n, _k, _0); \
349  Two_Sum(_1, _0, _j, x3); \
350  Two_Sum(_2, _j, _i, _1); \
351  Two_Sum(_l, _i, _m, _2); \
352  Two_Sum(_1, _k, _i, x4); \
353  Two_Sum(_2, _i, _k, x5); \
354  Two_Sum(_m, _k, x7, x6)
355 
356 /* An expansion of length two can be squared more quickly than finding the */
357 /* product of two different expansions of length two, and the result is */
358 /* guaranteed to have no more than six (rather than eight) components. */
359 
360 #define Two_Square(a1, a0, x5, x4, x3, x2, x1, x0) \
361  Square(a0, _j, x0); \
362  _0 = a0 + a0; \
363  Two_Product(a1, _0, _k, _1); \
364  Two_One_Sum(_k, _1, _j, _l, _2, x1); \
365  Square(a1, _j, _1); \
366  Two_Two_Sum(_j, _1, _l, _2, x5, x4, x3, x2)
367 
368 /* splitter = 2^ceiling(p / 2) + 1. Used to split floats in half. */
369 static REAL splitter;
370 static REAL epsilon; /* = 2^(-p). Used to estimate roundoff errors. */
371 /* A set of coefficients used to calculate maximum roundoff errors. */
372 static REAL resulterrbound;
373 static REAL ccwerrboundA, ccwerrboundB, ccwerrboundC;
374 static REAL o3derrboundA, o3derrboundB, o3derrboundC;
375 static REAL iccerrboundA, iccerrboundB, iccerrboundC;
376 static REAL isperrboundA, isperrboundB, isperrboundC;
377 
378 /*****************************************************************************/
379 /* */
380 /* doubleprint() Print the bit representation of a double. */
381 /* */
382 /* Useful for debugging exact arithmetic routines. */
383 /* */
384 /*****************************************************************************/
385 
386 /*
387 void doubleprint(number)
388 double number;
389 {
390  unsigned long long no;
391  unsigned long long sign, expo;
392  int exponent;
393  int i, bottomi;
394 
395  no = *(unsigned long long *) &number;
396  sign = no & 0x8000000000000000ll;
397  expo = (no >> 52) & 0x7ffll;
398  exponent = (int) expo;
399  exponent = exponent - 1023;
400  if (sign) {
401  printf("-");
402  } else {
403  printf(" ");
404  }
405  if (exponent == -1023) {
406  printf(
407  "0.0000000000000000000000000000000000000000000000000000_ ( )");
408  } else {
409  printf("1.");
410  bottomi = -1;
411  for (i = 0; i < 52; i++) {
412  if (no & 0x0008000000000000ll) {
413  printf("1");
414  bottomi = i;
415  } else {
416  printf("0");
417  }
418  no <<= 1;
419  }
420  printf("_%d (%d)", exponent, exponent - 1 - bottomi);
421  }
422 }
423 */
424 
425 /*****************************************************************************/
426 /* */
427 /* floatprint() Print the bit representation of a float. */
428 /* */
429 /* Useful for debugging exact arithmetic routines. */
430 /* */
431 /*****************************************************************************/
432 
433 /*
434 void floatprint(number)
435 float number;
436 {
437  unsigned no;
438  unsigned sign, expo;
439  int exponent;
440  int i, bottomi;
441 
442  no = *(unsigned *) &number;
443  sign = no & 0x80000000;
444  expo = (no >> 23) & 0xff;
445  exponent = (int) expo;
446  exponent = exponent - 127;
447  if (sign) {
448  printf("-");
449  } else {
450  printf(" ");
451  }
452  if (exponent == -127) {
453  printf("0.00000000000000000000000_ ( )");
454  } else {
455  printf("1.");
456  bottomi = -1;
457  for (i = 0; i < 23; i++) {
458  if (no & 0x00400000) {
459  printf("1");
460  bottomi = i;
461  } else {
462  printf("0");
463  }
464  no <<= 1;
465  }
466  printf("_%3d (%3d)", exponent, exponent - 1 - bottomi);
467  }
468 }
469 */
470 
471 /*****************************************************************************/
472 /* */
473 /* expansion_print() Print the bit representation of an expansion. */
474 /* */
475 /* Useful for debugging exact arithmetic routines. */
476 /* */
477 /*****************************************************************************/
478 
479 /*
480 void expansion_print(elen, e)
481 int elen;
482 REAL *e;
483 {
484  int i;
485 
486  for (i = elen - 1; i >= 0; i--) {
487  REALPRINT(e[i]);
488  if (i > 0) {
489  printf(" +\n");
490  } else {
491  printf("\n");
492  }
493  }
494 }
495 */
496 
497 /*****************************************************************************/
498 /* */
499 /* doublerand() Generate a double with random 53-bit significand and a */
500 /* random exponent in [0, 511]. */
501 /* */
502 /*****************************************************************************/
503 
504 /*
505 double doublerand()
506 {
507  double result;
508  double expo;
509  long a, b, c;
510  long i;
511 
512  a = random();
513  b = random();
514  c = random();
515  result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8);
516  for (i = 512, expo = 2; i <= 131072; i *= 2, expo = expo * expo) {
517  if (c & i) {
518  result *= expo;
519  }
520  }
521  return result;
522 }
523 */
524 
525 /*****************************************************************************/
526 /* */
527 /* narrowdoublerand() Generate a double with random 53-bit significand */
528 /* and a random exponent in [0, 7]. */
529 /* */
530 /*****************************************************************************/
531 
532 /*
533 double narrowdoublerand()
534 {
535  double result;
536  double expo;
537  long a, b, c;
538  long i;
539 
540  a = random();
541  b = random();
542  c = random();
543  result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8);
544  for (i = 512, expo = 2; i <= 2048; i *= 2, expo = expo * expo) {
545  if (c & i) {
546  result *= expo;
547  }
548  }
549  return result;
550 }
551 */
552 
553 /*****************************************************************************/
554 /* */
555 /* uniformdoublerand() Generate a double with random 53-bit significand. */
556 /* */
557 /*****************************************************************************/
558 
559 /*
560 double uniformdoublerand()
561 {
562  double result;
563  long a, b;
564 
565  a = random();
566  b = random();
567  result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8);
568  return result;
569 }
570 */
571 
572 /*****************************************************************************/
573 /* */
574 /* floatrand() Generate a float with random 24-bit significand and a */
575 /* random exponent in [0, 63]. */
576 /* */
577 /*****************************************************************************/
578 
579 /*
580 float floatrand()
581 {
582  float result;
583  float expo;
584  long a, c;
585  long i;
586 
587  a = random();
588  c = random();
589  result = (float) ((a - 1073741824) >> 6);
590  for (i = 512, expo = 2; i <= 16384; i *= 2, expo = expo * expo) {
591  if (c & i) {
592  result *= expo;
593  }
594  }
595  return result;
596 }
597 */
598 
599 /*****************************************************************************/
600 /* */
601 /* narrowfloatrand() Generate a float with random 24-bit significand and */
602 /* a random exponent in [0, 7]. */
603 /* */
604 /*****************************************************************************/
605 
606 /*
607 float narrowfloatrand()
608 {
609  float result;
610  float expo;
611  long a, c;
612  long i;
613 
614  a = random();
615  c = random();
616  result = (float) ((a - 1073741824) >> 6);
617  for (i = 512, expo = 2; i <= 2048; i *= 2, expo = expo * expo) {
618  if (c & i) {
619  result *= expo;
620  }
621  }
622  return result;
623 }
624 */
625 
626 /*****************************************************************************/
627 /* */
628 /* uniformfloatrand() Generate a float with random 24-bit significand. */
629 /* */
630 /*****************************************************************************/
631 
632 /*
633 float uniformfloatrand()
634 {
635  float result;
636  long a;
637 
638  a = random();
639  result = (float) ((a - 1073741824) >> 6);
640  return result;
641 }
642 */
643 namespace tetgen
644 { //Added namespace to avoid clash with triangle
645 
646 /*****************************************************************************/
647 /* */
648 /* exactinit() Initialize the variables used for exact arithmetic. */
649 /* */
650 /* `epsilon' is the largest power of two such that 1.0 + epsilon = 1.0 in */
651 /* floating-point arithmetic. `epsilon' bounds the relative roundoff */
652 /* error. It is used for floating-point error analysis. */
653 /* */
654 /* `splitter' is used to split floating-point numbers into two half- */
655 /* length significands for exact multiplication. */
656 /* */
657 /* I imagine that a highly optimizing compiler might be too smart for its */
658 /* own good, and somehow cause this routine to fail, if it pretends that */
659 /* floating-point arithmetic is too much like real arithmetic. */
660 /* */
661 /* Don't change this routine unless you fully understand it. */
662 /* */
663 /*****************************************************************************/
664 
665 REAL exactinit()
666 {
667  REAL half;
668  REAL check, lastcheck;
669  int every_other;
670 #ifdef LINUX
671  int cword;
672 #endif /* LINUX */
673 
674 #ifdef CPU86
675 #ifdef SINGLE
676  _control87(_PC_24, _MCW_PC); /* Set FPU control word for single precision. */
677 #else /* not SINGLE */
678  _control87(_PC_53, _MCW_PC); /* Set FPU control word for double precision. */
679 #endif /* not SINGLE */
680 #endif /* CPU86 */
681 #ifdef LINUX
682 #ifdef SINGLE
683  /* cword = 4223; */
684  cword = 4210; /* set FPU control word for single precision */
685 #else /* not SINGLE */
686  /* cword = 4735; */
687  cword = 4722; /* set FPU control word for double precision */
688 #endif /* not SINGLE */
689  _FPU_SETCW(cword);
690 #endif /* LINUX */
691 
692  every_other = 1;
693  half = 0.5;
694  epsilon = 1.0;
695  splitter = 1.0;
696  check = 1.0;
697  /* Repeatedly divide `epsilon' by two until it is too small to add to */
698  /* one without causing roundoff. (Also check if the sum is equal to */
699  /* the previous sum, for machines that round up instead of using exact */
700  /* rounding. Not that this library will work on such machines anyway. */
701  do {
702  lastcheck = check;
703  epsilon *= half;
704  if (every_other) {
705  splitter *= 2.0;
706  }
707  every_other = !every_other;
708  check = 1.0 + epsilon;
709  } while ((check != 1.0) && (check != lastcheck));
710  splitter += 1.0;
711 
712  /* Error bounds for orientation and incircle tests. */
713  resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
714  ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon;
715  ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon;
716  ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon;
717  o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon;
718  o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon;
719  o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;
720  iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon;
721  iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon;
722  iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon;
723  isperrboundA = (16.0 + 224.0 * epsilon) * epsilon;
724  isperrboundB = (5.0 + 72.0 * epsilon) * epsilon;
725  isperrboundC = (71.0 + 1408.0 * epsilon) * epsilon * epsilon;
726 
727  return epsilon; /* Added by H. Si 30 Juli, 2004. */
728 }
729 
730 /*****************************************************************************/
731 /* */
732 /* grow_expansion() Add a scalar to an expansion. */
733 /* */
734 /* Sets h = e + b. See the long version of my paper for details. */
735 /* */
736 /* Maintains the nonoverlapping property. If round-to-even is used (as */
737 /* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */
738 /* properties as well. (That is, if e has one of these properties, so */
739 /* will h.) */
740 /* */
741 /*****************************************************************************/
742 
743 int grow_expansion(int elen, REAL *e, REAL b, REAL *h)
744 /* e and h can be the same. */
745 {
746  REAL Q;
747  INEXACT REAL Qnew;
748  int eindex;
749  REAL enow;
750  INEXACT REAL bvirt;
751  REAL avirt, bround, around;
752 
753  Q = b;
754  for (eindex = 0; eindex < elen; eindex++) {
755  enow = e[eindex];
756  Two_Sum(Q, enow, Qnew, h[eindex]);
757  Q = Qnew;
758  }
759  h[eindex] = Q;
760  return eindex + 1;
761 }
762 
763 /*****************************************************************************/
764 /* */
765 /* grow_expansion_zeroelim() Add a scalar to an expansion, eliminating */
766 /* zero components from the output expansion. */
767 /* */
768 /* Sets h = e + b. See the long version of my paper for details. */
769 /* */
770 /* Maintains the nonoverlapping property. If round-to-even is used (as */
771 /* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */
772 /* properties as well. (That is, if e has one of these properties, so */
773 /* will h.) */
774 /* */
775 /*****************************************************************************/
776 
777 int grow_expansion_zeroelim(int elen, REAL *e, REAL b, REAL *h)
778 /* e and h can be the same. */
779 {
780  REAL Q, hh;
781  INEXACT REAL Qnew;
782  int eindex, hindex;
783  REAL enow;
784  INEXACT REAL bvirt;
785  REAL avirt, bround, around;
786 
787  hindex = 0;
788  Q = b;
789  for (eindex = 0; eindex < elen; eindex++) {
790  enow = e[eindex];
791  Two_Sum(Q, enow, Qnew, hh);
792  Q = Qnew;
793  if (hh != 0.0) {
794  h[hindex++] = hh;
795  }
796  }
797  if ((Q != 0.0) || (hindex == 0)) {
798  h[hindex++] = Q;
799  }
800  return hindex;
801 }
802 
803 /*****************************************************************************/
804 /* */
805 /* expansion_sum() Sum two expansions. */
806 /* */
807 /* Sets h = e + f. See the long version of my paper for details. */
808 /* */
809 /* Maintains the nonoverlapping property. If round-to-even is used (as */
810 /* with IEEE 754), maintains the nonadjacent property as well. (That is, */
811 /* if e has one of these properties, so will h.) Does NOT maintain the */
812 /* strongly nonoverlapping property. */
813 /* */
814 /*****************************************************************************/
815 
816 int expansion_sum(int elen, REAL *e, int flen, REAL *f, REAL *h)
817 /* e and h can be the same, but f and h cannot. */
818 {
819  REAL Q;
820  INEXACT REAL Qnew;
821  int findex, hindex, hlast;
822  REAL hnow;
823  INEXACT REAL bvirt;
824  REAL avirt, bround, around;
825 
826  Q = f[0];
827  for (hindex = 0; hindex < elen; hindex++) {
828  hnow = e[hindex];
829  Two_Sum(Q, hnow, Qnew, h[hindex]);
830  Q = Qnew;
831  }
832  h[hindex] = Q;
833  hlast = hindex;
834  for (findex = 1; findex < flen; findex++) {
835  Q = f[findex];
836  for (hindex = findex; hindex <= hlast; hindex++) {
837  hnow = h[hindex];
838  Two_Sum(Q, hnow, Qnew, h[hindex]);
839  Q = Qnew;
840  }
841  h[++hlast] = Q;
842  }
843  return hlast + 1;
844 }
845 
846 /*****************************************************************************/
847 /* */
848 /* expansion_sum_zeroelim1() Sum two expansions, eliminating zero */
849 /* components from the output expansion. */
850 /* */
851 /* Sets h = e + f. See the long version of my paper for details. */
852 /* */
853 /* Maintains the nonoverlapping property. If round-to-even is used (as */
854 /* with IEEE 754), maintains the nonadjacent property as well. (That is, */
855 /* if e has one of these properties, so will h.) Does NOT maintain the */
856 /* strongly nonoverlapping property. */
857 /* */
858 /*****************************************************************************/
859 
860 int expansion_sum_zeroelim1(int elen, REAL *e, int flen, REAL *f, REAL *h)
861 /* e and h can be the same, but f and h cannot. */
862 {
863  REAL Q;
864  INEXACT REAL Qnew;
865  int index, findex, hindex, hlast;
866  REAL hnow;
867  INEXACT REAL bvirt;
868  REAL avirt, bround, around;
869 
870  Q = f[0];
871  for (hindex = 0; hindex < elen; hindex++) {
872  hnow = e[hindex];
873  Two_Sum(Q, hnow, Qnew, h[hindex]);
874  Q = Qnew;
875  }
876  h[hindex] = Q;
877  hlast = hindex;
878  for (findex = 1; findex < flen; findex++) {
879  Q = f[findex];
880  for (hindex = findex; hindex <= hlast; hindex++) {
881  hnow = h[hindex];
882  Two_Sum(Q, hnow, Qnew, h[hindex]);
883  Q = Qnew;
884  }
885  h[++hlast] = Q;
886  }
887  hindex = -1;
888  for (index = 0; index <= hlast; index++) {
889  hnow = h[index];
890  if (hnow != 0.0) {
891  h[++hindex] = hnow;
892  }
893  }
894  if (hindex == -1) {
895  return 1;
896  } else {
897  return hindex + 1;
898  }
899 }
900 
901 /*****************************************************************************/
902 /* */
903 /* expansion_sum_zeroelim2() Sum two expansions, eliminating zero */
904 /* components from the output expansion. */
905 /* */
906 /* Sets h = e + f. See the long version of my paper for details. */
907 /* */
908 /* Maintains the nonoverlapping property. If round-to-even is used (as */
909 /* with IEEE 754), maintains the nonadjacent property as well. (That is, */
910 /* if e has one of these properties, so will h.) Does NOT maintain the */
911 /* strongly nonoverlapping property. */
912 /* */
913 /*****************************************************************************/
914 
915 int expansion_sum_zeroelim2(int elen, REAL *e, int flen, REAL *f, REAL *h)
916 /* e and h can be the same, but f and h cannot. */
917 {
918  REAL Q, hh;
919  INEXACT REAL Qnew;
920  int eindex, findex, hindex, hlast;
921  REAL enow;
922  INEXACT REAL bvirt;
923  REAL avirt, bround, around;
924 
925  hindex = 0;
926  Q = f[0];
927  for (eindex = 0; eindex < elen; eindex++) {
928  enow = e[eindex];
929  Two_Sum(Q, enow, Qnew, hh);
930  Q = Qnew;
931  if (hh != 0.0) {
932  h[hindex++] = hh;
933  }
934  }
935  h[hindex] = Q;
936  hlast = hindex;
937  for (findex = 1; findex < flen; findex++) {
938  hindex = 0;
939  Q = f[findex];
940  for (eindex = 0; eindex <= hlast; eindex++) {
941  enow = h[eindex];
942  Two_Sum(Q, enow, Qnew, hh);
943  Q = Qnew;
944  if (hh != 0) {
945  h[hindex++] = hh;
946  }
947  }
948  h[hindex] = Q;
949  hlast = hindex;
950  }
951  return hlast + 1;
952 }
953 
954 /*****************************************************************************/
955 /* */
956 /* fast_expansion_sum() Sum two expansions. */
957 /* */
958 /* Sets h = e + f. See the long version of my paper for details. */
959 /* */
960 /* If round-to-even is used (as with IEEE 754), maintains the strongly */
961 /* nonoverlapping property. (That is, if e is strongly nonoverlapping, h */
962 /* will be also.) Does NOT maintain the nonoverlapping or nonadjacent */
963 /* properties. */
964 /* */
965 /*****************************************************************************/
966 
967 int fast_expansion_sum(int elen, REAL *e, int flen, REAL *f, REAL *h)
968 /* h cannot be e or f. */
969 {
970  REAL Q;
971  INEXACT REAL Qnew;
972  INEXACT REAL bvirt;
973  REAL avirt, bround, around;
974  int eindex, findex, hindex;
975  REAL enow, fnow;
976 
977  enow = e[0];
978  fnow = f[0];
979  eindex = findex = 0;
980  if ((fnow > enow) == (fnow > -enow)) {
981  Q = enow;
982  enow = e[++eindex];
983  } else {
984  Q = fnow;
985  fnow = f[++findex];
986  }
987  hindex = 0;
988  if ((eindex < elen) && (findex < flen)) {
989  if ((fnow > enow) == (fnow > -enow)) {
990  Fast_Two_Sum(enow, Q, Qnew, h[0]);
991  enow = e[++eindex];
992  } else {
993  Fast_Two_Sum(fnow, Q, Qnew, h[0]);
994  fnow = f[++findex];
995  }
996  Q = Qnew;
997  hindex = 1;
998  while ((eindex < elen) && (findex < flen)) {
999  if ((fnow > enow) == (fnow > -enow)) {
1000  Two_Sum(Q, enow, Qnew, h[hindex]);
1001  enow = e[++eindex];
1002  } else {
1003  Two_Sum(Q, fnow, Qnew, h[hindex]);
1004  fnow = f[++findex];
1005  }
1006  Q = Qnew;
1007  hindex++;
1008  }
1009  }
1010  while (eindex < elen) {
1011  Two_Sum(Q, enow, Qnew, h[hindex]);
1012  enow = e[++eindex];
1013  Q = Qnew;
1014  hindex++;
1015  }
1016  while (findex < flen) {
1017  Two_Sum(Q, fnow, Qnew, h[hindex]);
1018  fnow = f[++findex];
1019  Q = Qnew;
1020  hindex++;
1021  }
1022  h[hindex] = Q;
1023  return hindex + 1;
1024 }
1025 
1026 /*****************************************************************************/
1027 /* */
1028 /* fast_expansion_sum_zeroelim() Sum two expansions, eliminating zero */
1029 /* components from the output expansion. */
1030 /* */
1031 /* Sets h = e + f. See the long version of my paper for details. */
1032 /* */
1033 /* If round-to-even is used (as with IEEE 754), maintains the strongly */
1034 /* nonoverlapping property. (That is, if e is strongly nonoverlapping, h */
1035 /* will be also.) Does NOT maintain the nonoverlapping or nonadjacent */
1036 /* properties. */
1037 /* */
1038 /*****************************************************************************/
1039 
1040 int fast_expansion_sum_zeroelim(int elen, REAL *e, int flen, REAL *f, REAL *h)
1041 /* h cannot be e or f. */
1042 {
1043  REAL Q;
1044  INEXACT REAL Qnew;
1045  INEXACT REAL hh;
1046  INEXACT REAL bvirt;
1047  REAL avirt, bround, around;
1048  int eindex, findex, hindex;
1049  REAL enow, fnow;
1050 
1051  enow = e[0];
1052  fnow = f[0];
1053  eindex = findex = 0;
1054  if ((fnow > enow) == (fnow > -enow)) {
1055  Q = enow;
1056  enow = e[++eindex];
1057  } else {
1058  Q = fnow;
1059  fnow = f[++findex];
1060  }
1061  hindex = 0;
1062  if ((eindex < elen) && (findex < flen)) {
1063  if ((fnow > enow) == (fnow > -enow)) {
1064  Fast_Two_Sum(enow, Q, Qnew, hh);
1065  enow = e[++eindex];
1066  } else {
1067  Fast_Two_Sum(fnow, Q, Qnew, hh);
1068  fnow = f[++findex];
1069  }
1070  Q = Qnew;
1071  if (hh != 0.0) {
1072  h[hindex++] = hh;
1073  }
1074  while ((eindex < elen) && (findex < flen)) {
1075  if ((fnow > enow) == (fnow > -enow)) {
1076  Two_Sum(Q, enow, Qnew, hh);
1077  enow = e[++eindex];
1078  } else {
1079  Two_Sum(Q, fnow, Qnew, hh);
1080  fnow = f[++findex];
1081  }
1082  Q = Qnew;
1083  if (hh != 0.0) {
1084  h[hindex++] = hh;
1085  }
1086  }
1087  }
1088  while (eindex < elen) {
1089  Two_Sum(Q, enow, Qnew, hh);
1090  enow = e[++eindex];
1091  Q = Qnew;
1092  if (hh != 0.0) {
1093  h[hindex++] = hh;
1094  }
1095  }
1096  while (findex < flen) {
1097  Two_Sum(Q, fnow, Qnew, hh);
1098  fnow = f[++findex];
1099  Q = Qnew;
1100  if (hh != 0.0) {
1101  h[hindex++] = hh;
1102  }
1103  }
1104  if ((Q != 0.0) || (hindex == 0)) {
1105  h[hindex++] = Q;
1106  }
1107  return hindex;
1108 }
1109 
1110 /*****************************************************************************/
1111 /* */
1112 /* linear_expansion_sum() Sum two expansions. */
1113 /* */
1114 /* Sets h = e + f. See either version of my paper for details. */
1115 /* */
1116 /* Maintains the nonoverlapping property. (That is, if e is */
1117 /* nonoverlapping, h will be also.) */
1118 /* */
1119 /*****************************************************************************/
1120 
1121 int linear_expansion_sum(int elen, REAL *e, int flen, REAL *f, REAL *h)
1122 /* h cannot be e or f. */
1123 {
1124  REAL Q, q;
1125  INEXACT REAL Qnew;
1126  INEXACT REAL R;
1127  INEXACT REAL bvirt;
1128  REAL avirt, bround, around;
1129  int eindex, findex, hindex;
1130  REAL enow, fnow;
1131  REAL g0;
1132 
1133  enow = e[0];
1134  fnow = f[0];
1135  eindex = findex = 0;
1136  if ((fnow > enow) == (fnow > -enow)) {
1137  g0 = enow;
1138  enow = e[++eindex];
1139  } else {
1140  g0 = fnow;
1141  fnow = f[++findex];
1142  }
1143  if ((eindex < elen) && ((findex >= flen)
1144  || ((fnow > enow) == (fnow > -enow)))) {
1145  Fast_Two_Sum(enow, g0, Qnew, q);
1146  enow = e[++eindex];
1147  } else {
1148  Fast_Two_Sum(fnow, g0, Qnew, q);
1149  fnow = f[++findex];
1150  }
1151  Q = Qnew;
1152  for (hindex = 0; hindex < elen + flen - 2; hindex++) {
1153  if ((eindex < elen) && ((findex >= flen)
1154  || ((fnow > enow) == (fnow > -enow)))) {
1155  Fast_Two_Sum(enow, q, R, h[hindex]);
1156  enow = e[++eindex];
1157  } else {
1158  Fast_Two_Sum(fnow, q, R, h[hindex]);
1159  fnow = f[++findex];
1160  }
1161  Two_Sum(Q, R, Qnew, q);
1162  Q = Qnew;
1163  }
1164  h[hindex] = q;
1165  h[hindex + 1] = Q;
1166  return hindex + 2;
1167 }
1168 
1169 /*****************************************************************************/
1170 /* */
1171 /* linear_expansion_sum_zeroelim() Sum two expansions, eliminating zero */
1172 /* components from the output expansion. */
1173 /* */
1174 /* Sets h = e + f. See either version of my paper for details. */
1175 /* */
1176 /* Maintains the nonoverlapping property. (That is, if e is */
1177 /* nonoverlapping, h will be also.) */
1178 /* */
1179 /*****************************************************************************/
1180 
1181 int linear_expansion_sum_zeroelim(int elen, REAL *e, int flen, REAL *f,
1182  REAL *h)
1183 /* h cannot be e or f. */
1184 {
1185  REAL Q, q, hh;
1186  INEXACT REAL Qnew;
1187  INEXACT REAL R;
1188  INEXACT REAL bvirt;
1189  REAL avirt, bround, around;
1190  int eindex, findex, hindex;
1191  int count;
1192  REAL enow, fnow;
1193  REAL g0;
1194 
1195  enow = e[0];
1196  fnow = f[0];
1197  eindex = findex = 0;
1198  hindex = 0;
1199  if ((fnow > enow) == (fnow > -enow)) {
1200  g0 = enow;
1201  enow = e[++eindex];
1202  } else {
1203  g0 = fnow;
1204  fnow = f[++findex];
1205  }
1206  if ((eindex < elen) && ((findex >= flen)
1207  || ((fnow > enow) == (fnow > -enow)))) {
1208  Fast_Two_Sum(enow, g0, Qnew, q);
1209  enow = e[++eindex];
1210  } else {
1211  Fast_Two_Sum(fnow, g0, Qnew, q);
1212  fnow = f[++findex];
1213  }
1214  Q = Qnew;
1215  for (count = 2; count < elen + flen; count++) {
1216  if ((eindex < elen) && ((findex >= flen)
1217  || ((fnow > enow) == (fnow > -enow)))) {
1218  Fast_Two_Sum(enow, q, R, hh);
1219  enow = e[++eindex];
1220  } else {
1221  Fast_Two_Sum(fnow, q, R, hh);
1222  fnow = f[++findex];
1223  }
1224  Two_Sum(Q, R, Qnew, q);
1225  Q = Qnew;
1226  if (hh != 0) {
1227  h[hindex++] = hh;
1228  }
1229  }
1230  if (q != 0) {
1231  h[hindex++] = q;
1232  }
1233  if ((Q != 0.0) || (hindex == 0)) {
1234  h[hindex++] = Q;
1235  }
1236  return hindex;
1237 }
1238 
1239 /*****************************************************************************/
1240 /* */
1241 /* scale_expansion() Multiply an expansion by a scalar. */
1242 /* */
1243 /* Sets h = be. See either version of my paper for details. */
1244 /* */
1245 /* Maintains the nonoverlapping property. If round-to-even is used (as */
1246 /* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */
1247 /* properties as well. (That is, if e has one of these properties, so */
1248 /* will h.) */
1249 /* */
1250 /*****************************************************************************/
1251 
1252 int scale_expansion(int elen, REAL *e, REAL b, REAL *h)
1253 /* e and h cannot be the same. */
1254 {
1255  INEXACT REAL Q;
1256  INEXACT REAL sum;
1257  INEXACT REAL product1;
1258  REAL product0;
1259  int eindex, hindex;
1260  REAL enow;
1261  INEXACT REAL bvirt;
1262  REAL avirt, bround, around;
1263  INEXACT REAL c;
1264  INEXACT REAL abig;
1265  REAL ahi, alo, bhi, blo;
1266  REAL err1, err2, err3;
1267 
1268  Split(b, bhi, blo);
1269  Two_Product_Presplit(e[0], b, bhi, blo, Q, h[0]);
1270  hindex = 1;
1271  for (eindex = 1; eindex < elen; eindex++) {
1272  enow = e[eindex];
1273  Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
1274  Two_Sum(Q, product0, sum, h[hindex]);
1275  hindex++;
1276  Two_Sum(product1, sum, Q, h[hindex]);
1277  hindex++;
1278  }
1279  h[hindex] = Q;
1280  return elen + elen;
1281 }
1282 
1283 /*****************************************************************************/
1284 /* */
1285 /* scale_expansion_zeroelim() Multiply an expansion by a scalar, */
1286 /* eliminating zero components from the */
1287 /* output expansion. */
1288 /* */
1289 /* Sets h = be. See either version of my paper for details. */
1290 /* */
1291 /* Maintains the nonoverlapping property. If round-to-even is used (as */
1292 /* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */
1293 /* properties as well. (That is, if e has one of these properties, so */
1294 /* will h.) */
1295 /* */
1296 /*****************************************************************************/
1297 
1298 int scale_expansion_zeroelim(int elen, REAL *e, REAL b, REAL *h)
1299 /* e and h cannot be the same. */
1300 {
1301  INEXACT REAL Q, sum;
1302  REAL hh;
1303  INEXACT REAL product1;
1304  REAL product0;
1305  int eindex, hindex;
1306  REAL enow;
1307  INEXACT REAL bvirt;
1308  REAL avirt, bround, around;
1309  INEXACT REAL c;
1310  INEXACT REAL abig;
1311  REAL ahi, alo, bhi, blo;
1312  REAL err1, err2, err3;
1313 
1314  Split(b, bhi, blo);
1315  Two_Product_Presplit(e[0], b, bhi, blo, Q, hh);
1316  hindex = 0;
1317  if (hh != 0) {
1318  h[hindex++] = hh;
1319  }
1320  for (eindex = 1; eindex < elen; eindex++) {
1321  enow = e[eindex];
1322  Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
1323  Two_Sum(Q, product0, sum, hh);
1324  if (hh != 0) {
1325  h[hindex++] = hh;
1326  }
1327  Fast_Two_Sum(product1, sum, Q, hh);
1328  if (hh != 0) {
1329  h[hindex++] = hh;
1330  }
1331  }
1332  if ((Q != 0.0) || (hindex == 0)) {
1333  h[hindex++] = Q;
1334  }
1335  return hindex;
1336 }
1337 
1338 /*****************************************************************************/
1339 /* */
1340 /* compress() Compress an expansion. */
1341 /* */
1342 /* See the long version of my paper for details. */
1343 /* */
1344 /* Maintains the nonoverlapping property. If round-to-even is used (as */
1345 /* with IEEE 754), then any nonoverlapping expansion is converted to a */
1346 /* nonadjacent expansion. */
1347 /* */
1348 /*****************************************************************************/
1349 
1350 int compress(int elen, REAL *e, REAL *h)
1351 /* e and h may be the same. */
1352 {
1353  REAL Q, q;
1354  INEXACT REAL Qnew;
1355  int eindex, hindex;
1356  INEXACT REAL bvirt;
1357  REAL enow, hnow;
1358  int top, bottom;
1359 
1360  bottom = elen - 1;
1361  Q = e[bottom];
1362  for (eindex = elen - 2; eindex >= 0; eindex--) {
1363  enow = e[eindex];
1364  Fast_Two_Sum(Q, enow, Qnew, q);
1365  if (q != 0) {
1366  h[bottom--] = Qnew;
1367  Q = q;
1368  } else {
1369  Q = Qnew;
1370  }
1371  }
1372  top = 0;
1373  for (hindex = bottom + 1; hindex < elen; hindex++) {
1374  hnow = h[hindex];
1375  Fast_Two_Sum(hnow, Q, Qnew, q);
1376  if (q != 0) {
1377  h[top++] = q;
1378  }
1379  Q = Qnew;
1380  }
1381  h[top] = Q;
1382  return top + 1;
1383 }
1384 
1385 /*****************************************************************************/
1386 /* */
1387 /* estimate() Produce a one-word estimate of an expansion's value. */
1388 /* */
1389 /* See either version of my paper for details. */
1390 /* */
1391 /*****************************************************************************/
1392 
1393 REAL estimate(int elen, REAL *e)
1394 {
1395  REAL Q;
1396  int eindex;
1397 
1398  Q = e[0];
1399  for (eindex = 1; eindex < elen; eindex++) {
1400  Q += e[eindex];
1401  }
1402  return Q;
1403 }
1404 
1405 /*****************************************************************************/
1406 /* */
1407 /* orient2dfast() Approximate 2D orientation test. Nonrobust. */
1408 /* orient2dexact() Exact 2D orientation test. Robust. */
1409 /* orient2dslow() Another exact 2D orientation test. Robust. */
1410 /* orient2d() Adaptive exact 2D orientation test. Robust. */
1411 /* */
1412 /* Return a positive value if the points pa, pb, and pc occur */
1413 /* in counterclockwise order; a negative value if they occur */
1414 /* in clockwise order; and zero if they are collinear. The */
1415 /* result is also a rough approximation of twice the signed */
1416 /* area of the triangle defined by the three points. */
1417 /* */
1418 /* Only the first and last routine should be used; the middle two are for */
1419 /* timings. */
1420 /* */
1421 /* The last three use exact arithmetic to ensure a correct answer. The */
1422 /* result returned is the determinant of a matrix. In orient2d() only, */
1423 /* this determinant is computed adaptively, in the sense that exact */
1424 /* arithmetic is used only to the degree it is needed to ensure that the */
1425 /* returned value has the correct sign. Hence, orient2d() is usually quite */
1426 /* fast, but will run more slowly when the input points are collinear or */
1427 /* nearly so. */
1428 /* */
1429 /*****************************************************************************/
1430 
1431 REAL orient2dfast(REAL *pa, REAL *pb, REAL *pc)
1432 {
1433  REAL acx, bcx, acy, bcy;
1434 
1435  acx = pa[0] - pc[0];
1436  bcx = pb[0] - pc[0];
1437  acy = pa[1] - pc[1];
1438  bcy = pb[1] - pc[1];
1439  return acx * bcy - acy * bcx;
1440 }
1441 
1442 REAL orient2dexact(REAL *pa, REAL *pb, REAL *pc)
1443 {
1444  INEXACT REAL axby1, axcy1, bxcy1, bxay1, cxay1, cxby1;
1445  REAL axby0, axcy0, bxcy0, bxay0, cxay0, cxby0;
1446  REAL aterms[4], bterms[4], cterms[4];
1447  INEXACT REAL aterms3, bterms3, cterms3;
1448  REAL v[8], w[12];
1449  int vlength, wlength;
1450 
1451  INEXACT REAL bvirt;
1452  REAL avirt, bround, around;
1453  INEXACT REAL c;
1454  INEXACT REAL abig;
1455  REAL ahi, alo, bhi, blo;
1456  REAL err1, err2, err3;
1457  INEXACT REAL _i, _j;
1458  REAL _0;
1459 
1460  Two_Product(pa[0], pb[1], axby1, axby0);
1461  Two_Product(pa[0], pc[1], axcy1, axcy0);
1462  Two_Two_Diff(axby1, axby0, axcy1, axcy0,
1463  aterms3, aterms[2], aterms[1], aterms[0]);
1464  aterms[3] = aterms3;
1465 
1466  Two_Product(pb[0], pc[1], bxcy1, bxcy0);
1467  Two_Product(pb[0], pa[1], bxay1, bxay0);
1468  Two_Two_Diff(bxcy1, bxcy0, bxay1, bxay0,
1469  bterms3, bterms[2], bterms[1], bterms[0]);
1470  bterms[3] = bterms3;
1471 
1472  Two_Product(pc[0], pa[1], cxay1, cxay0);
1473  Two_Product(pc[0], pb[1], cxby1, cxby0);
1474  Two_Two_Diff(cxay1, cxay0, cxby1, cxby0,
1475  cterms3, cterms[2], cterms[1], cterms[0]);
1476  cterms[3] = cterms3;
1477 
1478  vlength = fast_expansion_sum_zeroelim(4, aterms, 4, bterms, v);
1479  wlength = fast_expansion_sum_zeroelim(vlength, v, 4, cterms, w);
1480 
1481  return w[wlength - 1];
1482 }
1483 
1484 REAL orient2dslow(REAL *pa, REAL *pb, REAL *pc)
1485 {
1486  INEXACT REAL acx, acy, bcx, bcy;
1487  REAL acxtail, acytail;
1488  REAL bcxtail, bcytail;
1489  REAL negate, negatetail;
1490  REAL axby[8], bxay[8];
1491  INEXACT REAL axby7, bxay7;
1492  REAL deter[16];
1493  int deterlen;
1494 
1495  INEXACT REAL bvirt;
1496  REAL avirt, bround, around;
1497  INEXACT REAL c;
1498  INEXACT REAL abig;
1499  REAL a0hi, a0lo, a1hi, a1lo, bhi, blo;
1500  REAL err1, err2, err3;
1501  INEXACT REAL _i, _j, _k, _l, _m, _n;
1502  REAL _0, _1, _2;
1503 
1504  Two_Diff(pa[0], pc[0], acx, acxtail);
1505  Two_Diff(pa[1], pc[1], acy, acytail);
1506  Two_Diff(pb[0], pc[0], bcx, bcxtail);
1507  Two_Diff(pb[1], pc[1], bcy, bcytail);
1508 
1509  Two_Two_Product(acx, acxtail, bcy, bcytail,
1510  axby7, axby[6], axby[5], axby[4],
1511  axby[3], axby[2], axby[1], axby[0]);
1512  axby[7] = axby7;
1513  negate = -acy;
1514  negatetail = -acytail;
1515  Two_Two_Product(bcx, bcxtail, negate, negatetail,
1516  bxay7, bxay[6], bxay[5], bxay[4],
1517  bxay[3], bxay[2], bxay[1], bxay[0]);
1518  bxay[7] = bxay7;
1519 
1520  deterlen = fast_expansion_sum_zeroelim(8, axby, 8, bxay, deter);
1521 
1522  return deter[deterlen - 1];
1523 }
1524 
1525 REAL orient2dadapt(REAL *pa, REAL *pb, REAL *pc, REAL detsum)
1526 {
1527  INEXACT REAL acx, acy, bcx, bcy;
1528  REAL acxtail, acytail, bcxtail, bcytail;
1529  INEXACT REAL detleft, detright;
1530  REAL detlefttail, detrighttail;
1531  REAL det, errbound;
1532  REAL B[4], C1[8], C2[12], D[16];
1533  INEXACT REAL B3;
1534  int C1length, C2length, Dlength;
1535  REAL u[4];
1536  INEXACT REAL u3;
1537  INEXACT REAL s1, t1;
1538  REAL s0, t0;
1539 
1540  INEXACT REAL bvirt;
1541  REAL avirt, bround, around;
1542  INEXACT REAL c;
1543  INEXACT REAL abig;
1544  REAL ahi, alo, bhi, blo;
1545  REAL err1, err2, err3;
1546  INEXACT REAL _i, _j;
1547  REAL _0;
1548 
1549  acx = (REAL) (pa[0] - pc[0]);
1550  bcx = (REAL) (pb[0] - pc[0]);
1551  acy = (REAL) (pa[1] - pc[1]);
1552  bcy = (REAL) (pb[1] - pc[1]);
1553 
1554  Two_Product(acx, bcy, detleft, detlefttail);
1555  Two_Product(acy, bcx, detright, detrighttail);
1556 
1557  Two_Two_Diff(detleft, detlefttail, detright, detrighttail,
1558  B3, B[2], B[1], B[0]);
1559  B[3] = B3;
1560 
1561  det = estimate(4, B);
1562  errbound = ccwerrboundB * detsum;
1563  if ((det >= errbound) || (-det >= errbound)) {
1564  return det;
1565  }
1566 
1567  Two_Diff_Tail(pa[0], pc[0], acx, acxtail);
1568  Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail);
1569  Two_Diff_Tail(pa[1], pc[1], acy, acytail);
1570  Two_Diff_Tail(pb[1], pc[1], bcy, bcytail);
1571 
1572  if ((acxtail == 0.0) && (acytail == 0.0)
1573  && (bcxtail == 0.0) && (bcytail == 0.0)) {
1574  return det;
1575  }
1576 
1577  errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det);
1578  det += (acx * bcytail + bcy * acxtail)
1579  - (acy * bcxtail + bcx * acytail);
1580  if ((det >= errbound) || (-det >= errbound)) {
1581  return det;
1582  }
1583 
1584  Two_Product(acxtail, bcy, s1, s0);
1585  Two_Product(acytail, bcx, t1, t0);
1586  Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
1587  u[3] = u3;
1588  C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1);
1589 
1590  Two_Product(acx, bcytail, s1, s0);
1591  Two_Product(acy, bcxtail, t1, t0);
1592  Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
1593  u[3] = u3;
1594  C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2);
1595 
1596  Two_Product(acxtail, bcytail, s1, s0);
1597  Two_Product(acytail, bcxtail, t1, t0);
1598  Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
1599  u[3] = u3;
1600  Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D);
1601 
1602  return(D[Dlength - 1]);
1603 }
1604 
1605 REAL orient2d(REAL *pa, REAL *pb, REAL *pc)
1606 {
1607  REAL detleft, detright, det;
1608  REAL detsum, errbound;
1609 
1610  detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]);
1611  detright = (pa[1] - pc[1]) * (pb[0] - pc[0]);
1612  det = detleft - detright;
1613 
1614  if (detleft > 0.0) {
1615  if (detright <= 0.0) {
1616  return det;
1617  } else {
1618  detsum = detleft + detright;
1619  }
1620  } else if (detleft < 0.0) {
1621  if (detright >= 0.0) {
1622  return det;
1623  } else {
1624  detsum = -detleft - detright;
1625  }
1626  } else {
1627  return det;
1628  }
1629 
1630  errbound = ccwerrboundA * detsum;
1631  if ((det >= errbound) || (-det >= errbound)) {
1632  return det;
1633  }
1634 
1635  return orient2dadapt(pa, pb, pc, detsum);
1636 }
1637 
1638 /*****************************************************************************/
1639 /* */
1640 /* orient3dfast() Approximate 3D orientation test. Nonrobust. */
1641 /* orient3dexact() Exact 3D orientation test. Robust. */
1642 /* orient3dslow() Another exact 3D orientation test. Robust. */
1643 /* orient3d() Adaptive exact 3D orientation test. Robust. */
1644 /* */
1645 /* Return a positive value if the point pd lies below the */
1646 /* plane passing through pa, pb, and pc; "below" is defined so */
1647 /* that pa, pb, and pc appear in counterclockwise order when */
1648 /* viewed from above the plane. Returns a negative value if */
1649 /* pd lies above the plane. Returns zero if the points are */
1650 /* coplanar. The result is also a rough approximation of six */
1651 /* times the signed volume of the tetrahedron defined by the */
1652 /* four points. */
1653 /* */
1654 /* Only the first and last routine should be used; the middle two are for */
1655 /* timings. */
1656 /* */
1657 /* The last three use exact arithmetic to ensure a correct answer. The */
1658 /* result returned is the determinant of a matrix. In orient3d() only, */
1659 /* this determinant is computed adaptively, in the sense that exact */
1660 /* arithmetic is used only to the degree it is needed to ensure that the */
1661 /* returned value has the correct sign. Hence, orient3d() is usually quite */
1662 /* fast, but will run more slowly when the input points are coplanar or */
1663 /* nearly so. */
1664 /* */
1665 /*****************************************************************************/
1666 
1667 REAL orient3dfast(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
1668 {
1669  REAL adx, bdx, cdx;
1670  REAL ady, bdy, cdy;
1671  REAL adz, bdz, cdz;
1672 
1673  adx = pa[0] - pd[0];
1674  bdx = pb[0] - pd[0];
1675  cdx = pc[0] - pd[0];
1676  ady = pa[1] - pd[1];
1677  bdy = pb[1] - pd[1];
1678  cdy = pc[1] - pd[1];
1679  adz = pa[2] - pd[2];
1680  bdz = pb[2] - pd[2];
1681  cdz = pc[2] - pd[2];
1682 
1683  return adx * (bdy * cdz - bdz * cdy)
1684  + bdx * (cdy * adz - cdz * ady)
1685  + cdx * (ady * bdz - adz * bdy);
1686 }
1687 
1688 REAL orient3dexact(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
1689 {
1690  INEXACT REAL axby1, bxcy1, cxdy1, dxay1, axcy1, bxdy1;
1691  INEXACT REAL bxay1, cxby1, dxcy1, axdy1, cxay1, dxby1;
1692  REAL axby0, bxcy0, cxdy0, dxay0, axcy0, bxdy0;
1693  REAL bxay0, cxby0, dxcy0, axdy0, cxay0, dxby0;
1694  REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4];
1695  REAL temp8[8];
1696  int templen;
1697  REAL abc[12], bcd[12], cda[12], dab[12];
1698  int abclen, bcdlen, cdalen, dablen;
1699  REAL adet[24], bdet[24], cdet[24], ddet[24];
1700  int alen, blen, clen, dlen;
1701  REAL abdet[48], cddet[48];
1702  int ablen, cdlen;
1703  REAL deter[96];
1704  int deterlen;
1705  int i;
1706 
1707  INEXACT REAL bvirt;
1708  REAL avirt, bround, around;
1709  INEXACT REAL c;
1710  INEXACT REAL abig;
1711  REAL ahi, alo, bhi, blo;
1712  REAL err1, err2, err3;
1713  INEXACT REAL _i, _j;
1714  REAL _0;
1715 
1716  Two_Product(pa[0], pb[1], axby1, axby0);
1717  Two_Product(pb[0], pa[1], bxay1, bxay0);
1718  Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]);
1719 
1720  Two_Product(pb[0], pc[1], bxcy1, bxcy0);
1721  Two_Product(pc[0], pb[1], cxby1, cxby0);
1722  Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]);
1723 
1724  Two_Product(pc[0], pd[1], cxdy1, cxdy0);
1725  Two_Product(pd[0], pc[1], dxcy1, dxcy0);
1726  Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]);
1727 
1728  Two_Product(pd[0], pa[1], dxay1, dxay0);
1729  Two_Product(pa[0], pd[1], axdy1, axdy0);
1730  Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]);
1731 
1732  Two_Product(pa[0], pc[1], axcy1, axcy0);
1733  Two_Product(pc[0], pa[1], cxay1, cxay0);
1734  Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]);
1735 
1736  Two_Product(pb[0], pd[1], bxdy1, bxdy0);
1737  Two_Product(pd[0], pb[1], dxby1, dxby0);
1738  Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]);
1739 
1740  templen = fast_expansion_sum_zeroelim(4, cd, 4, da, temp8);
1741  cdalen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, cda);
1742  templen = fast_expansion_sum_zeroelim(4, da, 4, ab, temp8);
1743  dablen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, dab);
1744  for (i = 0; i < 4; i++) {
1745  bd[i] = -bd[i];
1746  ac[i] = -ac[i];
1747  }
1748  templen = fast_expansion_sum_zeroelim(4, ab, 4, bc, temp8);
1749  abclen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, abc);
1750  templen = fast_expansion_sum_zeroelim(4, bc, 4, cd, temp8);
1751  bcdlen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, bcd);
1752 
1753  alen = scale_expansion_zeroelim(bcdlen, bcd, pa[2], adet);
1754  blen = scale_expansion_zeroelim(cdalen, cda, -pb[2], bdet);
1755  clen = scale_expansion_zeroelim(dablen, dab, pc[2], cdet);
1756  dlen = scale_expansion_zeroelim(abclen, abc, -pd[2], ddet);
1757 
1758  ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
1759  cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
1760  deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, deter);
1761 
1762  return deter[deterlen - 1];
1763 }
1764 
1765 REAL orient3dslow(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
1766 {
1767  INEXACT REAL adx, ady, adz, bdx, bdy, bdz, cdx, cdy, cdz;
1768  REAL adxtail, adytail, adztail;
1769  REAL bdxtail, bdytail, bdztail;
1770  REAL cdxtail, cdytail, cdztail;
1771  REAL negate, negatetail;
1772  INEXACT REAL axby7, bxcy7, axcy7, bxay7, cxby7, cxay7;
1773  REAL axby[8], bxcy[8], axcy[8], bxay[8], cxby[8], cxay[8];
1774  REAL temp16[16], temp32[32], temp32t[32];
1775  int temp16len, temp32len, temp32tlen;
1776  REAL adet[64], bdet[64], cdet[64];
1777  int alen, blen, clen;
1778  REAL abdet[128];
1779  int ablen;
1780  REAL deter[192];
1781  int deterlen;
1782 
1783  INEXACT REAL bvirt;
1784  REAL avirt, bround, around;
1785  INEXACT REAL c;
1786  INEXACT REAL abig;
1787  REAL a0hi, a0lo, a1hi, a1lo, bhi, blo;
1788  REAL err1, err2, err3;
1789  INEXACT REAL _i, _j, _k, _l, _m, _n;
1790  REAL _0, _1, _2;
1791 
1792  Two_Diff(pa[0], pd[0], adx, adxtail);
1793  Two_Diff(pa[1], pd[1], ady, adytail);
1794  Two_Diff(pa[2], pd[2], adz, adztail);
1795  Two_Diff(pb[0], pd[0], bdx, bdxtail);
1796  Two_Diff(pb[1], pd[1], bdy, bdytail);
1797  Two_Diff(pb[2], pd[2], bdz, bdztail);
1798  Two_Diff(pc[0], pd[0], cdx, cdxtail);
1799  Two_Diff(pc[1], pd[1], cdy, cdytail);
1800  Two_Diff(pc[2], pd[2], cdz, cdztail);
1801 
1802  Two_Two_Product(adx, adxtail, bdy, bdytail,
1803  axby7, axby[6], axby[5], axby[4],
1804  axby[3], axby[2], axby[1], axby[0]);
1805  axby[7] = axby7;
1806  negate = -ady;
1807  negatetail = -adytail;
1808  Two_Two_Product(bdx, bdxtail, negate, negatetail,
1809  bxay7, bxay[6], bxay[5], bxay[4],
1810  bxay[3], bxay[2], bxay[1], bxay[0]);
1811  bxay[7] = bxay7;
1812  Two_Two_Product(bdx, bdxtail, cdy, cdytail,
1813  bxcy7, bxcy[6], bxcy[5], bxcy[4],
1814  bxcy[3], bxcy[2], bxcy[1], bxcy[0]);
1815  bxcy[7] = bxcy7;
1816  negate = -bdy;
1817  negatetail = -bdytail;
1818  Two_Two_Product(cdx, cdxtail, negate, negatetail,
1819  cxby7, cxby[6], cxby[5], cxby[4],
1820  cxby[3], cxby[2], cxby[1], cxby[0]);
1821  cxby[7] = cxby7;
1822  Two_Two_Product(cdx, cdxtail, ady, adytail,
1823  cxay7, cxay[6], cxay[5], cxay[4],
1824  cxay[3], cxay[2], cxay[1], cxay[0]);
1825  cxay[7] = cxay7;
1826  negate = -cdy;
1827  negatetail = -cdytail;
1828  Two_Two_Product(adx, adxtail, negate, negatetail,
1829  axcy7, axcy[6], axcy[5], axcy[4],
1830  axcy[3], axcy[2], axcy[1], axcy[0]);
1831  axcy[7] = axcy7;
1832 
1833  temp16len = fast_expansion_sum_zeroelim(8, bxcy, 8, cxby, temp16);
1834  temp32len = scale_expansion_zeroelim(temp16len, temp16, adz, temp32);
1835  temp32tlen = scale_expansion_zeroelim(temp16len, temp16, adztail, temp32t);
1836  alen = fast_expansion_sum_zeroelim(temp32len, temp32, temp32tlen, temp32t,
1837  adet);
1838 
1839  temp16len = fast_expansion_sum_zeroelim(8, cxay, 8, axcy, temp16);
1840  temp32len = scale_expansion_zeroelim(temp16len, temp16, bdz, temp32);
1841  temp32tlen = scale_expansion_zeroelim(temp16len, temp16, bdztail, temp32t);
1842  blen = fast_expansion_sum_zeroelim(temp32len, temp32, temp32tlen, temp32t,
1843  bdet);
1844 
1845  temp16len = fast_expansion_sum_zeroelim(8, axby, 8, bxay, temp16);
1846  temp32len = scale_expansion_zeroelim(temp16len, temp16, cdz, temp32);
1847  temp32tlen = scale_expansion_zeroelim(temp16len, temp16, cdztail, temp32t);
1848  clen = fast_expansion_sum_zeroelim(temp32len, temp32, temp32tlen, temp32t,
1849  cdet);
1850 
1851  ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
1852  deterlen = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, deter);
1853 
1854  return deter[deterlen - 1];
1855 }
1856 
1857 REAL orient3dadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL permanent)
1858 {
1859  INEXACT REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
1860  REAL det, errbound;
1861 
1862  INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
1863  REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
1864  REAL bc[4], ca[4], ab[4];
1865  INEXACT REAL bc3, ca3, ab3;
1866  REAL adet[8], bdet[8], cdet[8];
1867  int alen, blen, clen;
1868  REAL abdet[16];
1869  int ablen;
1870  REAL *finnow, *finother, *finswap;
1871  REAL fin1[192], fin2[192];
1872  int finlength;
1873 
1875  // To avoid uninitialized warnings reported by valgrind.
1876  int i;
1877  for (i = 0; i < 8; i++) {
1878  adet[i] = bdet[i] = cdet[i] = 0.0;
1879  }
1880  for (i = 0; i < 16; i++) {
1881  abdet[i] = 0.0;
1882  }
1884 
1885  REAL adxtail, bdxtail, cdxtail;
1886  REAL adytail, bdytail, cdytail;
1887  REAL adztail, bdztail, cdztail;
1888  INEXACT REAL at_blarge, at_clarge;
1889  INEXACT REAL bt_clarge, bt_alarge;
1890  INEXACT REAL ct_alarge, ct_blarge;
1891  REAL at_b[4], at_c[4], bt_c[4], bt_a[4], ct_a[4], ct_b[4];
1892  int at_blen, at_clen, bt_clen, bt_alen, ct_alen, ct_blen;
1893  INEXACT REAL bdxt_cdy1, cdxt_bdy1, cdxt_ady1;
1894  INEXACT REAL adxt_cdy1, adxt_bdy1, bdxt_ady1;
1895  REAL bdxt_cdy0, cdxt_bdy0, cdxt_ady0;
1896  REAL adxt_cdy0, adxt_bdy0, bdxt_ady0;
1897  INEXACT REAL bdyt_cdx1, cdyt_bdx1, cdyt_adx1;
1898  INEXACT REAL adyt_cdx1, adyt_bdx1, bdyt_adx1;
1899  REAL bdyt_cdx0, cdyt_bdx0, cdyt_adx0;
1900  REAL adyt_cdx0, adyt_bdx0, bdyt_adx0;
1901  REAL bct[8], cat[8], abt[8];
1902  int bctlen, catlen, abtlen;
1903  INEXACT REAL bdxt_cdyt1, cdxt_bdyt1, cdxt_adyt1;
1904  INEXACT REAL adxt_cdyt1, adxt_bdyt1, bdxt_adyt1;
1905  REAL bdxt_cdyt0, cdxt_bdyt0, cdxt_adyt0;
1906  REAL adxt_cdyt0, adxt_bdyt0, bdxt_adyt0;
1907  REAL u[4], v[12], w[16];
1908  INEXACT REAL u3;
1909  int vlength, wlength;
1910  REAL negate;
1911 
1912  INEXACT REAL bvirt;
1913  REAL avirt, bround, around;
1914  INEXACT REAL c;
1915  INEXACT REAL abig;
1916  REAL ahi, alo, bhi, blo;
1917  REAL err1, err2, err3;
1918  INEXACT REAL _i, _j, _k;
1919  REAL _0;
1920 
1921  adx = (REAL) (pa[0] - pd[0]);
1922  bdx = (REAL) (pb[0] - pd[0]);
1923  cdx = (REAL) (pc[0] - pd[0]);
1924  ady = (REAL) (pa[1] - pd[1]);
1925  bdy = (REAL) (pb[1] - pd[1]);
1926  cdy = (REAL) (pc[1] - pd[1]);
1927  adz = (REAL) (pa[2] - pd[2]);
1928  bdz = (REAL) (pb[2] - pd[2]);
1929  cdz = (REAL) (pc[2] - pd[2]);
1930 
1931  Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
1932  Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
1933  Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
1934  bc[3] = bc3;
1935  alen = scale_expansion_zeroelim(4, bc, adz, adet);
1936 
1937  Two_Product(cdx, ady, cdxady1, cdxady0);
1938  Two_Product(adx, cdy, adxcdy1, adxcdy0);
1939  Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
1940  ca[3] = ca3;
1941  blen = scale_expansion_zeroelim(4, ca, bdz, bdet);
1942 
1943  Two_Product(adx, bdy, adxbdy1, adxbdy0);
1944  Two_Product(bdx, ady, bdxady1, bdxady0);
1945  Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
1946  ab[3] = ab3;
1947  clen = scale_expansion_zeroelim(4, ab, cdz, cdet);
1948 
1949  ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
1950  finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
1951 
1952  det = estimate(finlength, fin1);
1953  errbound = o3derrboundB * permanent;
1954  if ((det >= errbound) || (-det >= errbound)) {
1955  return det;
1956  }
1957 
1958  Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
1959  Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
1960  Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
1961  Two_Diff_Tail(pa[1], pd[1], ady, adytail);
1962  Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
1963  Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
1964  Two_Diff_Tail(pa[2], pd[2], adz, adztail);
1965  Two_Diff_Tail(pb[2], pd[2], bdz, bdztail);
1966  Two_Diff_Tail(pc[2], pd[2], cdz, cdztail);
1967 
1968  if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
1969  && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)
1970  && (adztail == 0.0) && (bdztail == 0.0) && (cdztail == 0.0)) {
1971  return det;
1972  }
1973 
1974  errbound = o3derrboundC * permanent + resulterrbound * Absolute(det);
1975  det += (adz * ((bdx * cdytail + cdy * bdxtail)
1976  - (bdy * cdxtail + cdx * bdytail))
1977  + adztail * (bdx * cdy - bdy * cdx))
1978  + (bdz * ((cdx * adytail + ady * cdxtail)
1979  - (cdy * adxtail + adx * cdytail))
1980  + bdztail * (cdx * ady - cdy * adx))
1981  + (cdz * ((adx * bdytail + bdy * adxtail)
1982  - (ady * bdxtail + bdx * adytail))
1983  + cdztail * (adx * bdy - ady * bdx));
1984  if ((det >= errbound) || (-det >= errbound)) {
1985  return det;
1986  }
1987 
1988  finnow = fin1;
1989  finother = fin2;
1990 
1991  if (adxtail == 0.0) {
1992  if (adytail == 0.0) {
1993  at_b[0] = 0.0;
1994  at_blen = 1;
1995  at_c[0] = 0.0;
1996  at_clen = 1;
1997  } else {
1998  negate = -adytail;
1999  Two_Product(negate, bdx, at_blarge, at_b[0]);
2000  at_b[1] = at_blarge;
2001  at_blen = 2;
2002  Two_Product(adytail, cdx, at_clarge, at_c[0]);
2003  at_c[1] = at_clarge;
2004  at_clen = 2;
2005  }
2006  } else {
2007  if (adytail == 0.0) {
2008  Two_Product(adxtail, bdy, at_blarge, at_b[0]);
2009  at_b[1] = at_blarge;
2010  at_blen = 2;
2011  negate = -adxtail;
2012  Two_Product(negate, cdy, at_clarge, at_c[0]);
2013  at_c[1] = at_clarge;
2014  at_clen = 2;
2015  } else {
2016  Two_Product(adxtail, bdy, adxt_bdy1, adxt_bdy0);
2017  Two_Product(adytail, bdx, adyt_bdx1, adyt_bdx0);
2018  Two_Two_Diff(adxt_bdy1, adxt_bdy0, adyt_bdx1, adyt_bdx0,
2019  at_blarge, at_b[2], at_b[1], at_b[0]);
2020  at_b[3] = at_blarge;
2021  at_blen = 4;
2022  Two_Product(adytail, cdx, adyt_cdx1, adyt_cdx0);
2023  Two_Product(adxtail, cdy, adxt_cdy1, adxt_cdy0);
2024  Two_Two_Diff(adyt_cdx1, adyt_cdx0, adxt_cdy1, adxt_cdy0,
2025  at_clarge, at_c[2], at_c[1], at_c[0]);
2026  at_c[3] = at_clarge;
2027  at_clen = 4;
2028  }
2029  }
2030  if (bdxtail == 0.0) {
2031  if (bdytail == 0.0) {
2032  bt_c[0] = 0.0;
2033  bt_clen = 1;
2034  bt_a[0] = 0.0;
2035  bt_alen = 1;
2036  } else {
2037  negate = -bdytail;
2038  Two_Product(negate, cdx, bt_clarge, bt_c[0]);
2039  bt_c[1] = bt_clarge;
2040  bt_clen = 2;
2041  Two_Product(bdytail, adx, bt_alarge, bt_a[0]);
2042  bt_a[1] = bt_alarge;
2043  bt_alen = 2;
2044  }
2045  } else {
2046  if (bdytail == 0.0) {
2047  Two_Product(bdxtail, cdy, bt_clarge, bt_c[0]);
2048  bt_c[1] = bt_clarge;
2049  bt_clen = 2;
2050  negate = -bdxtail;
2051  Two_Product(negate, ady, bt_alarge, bt_a[0]);
2052  bt_a[1] = bt_alarge;
2053  bt_alen = 2;
2054  } else {
2055  Two_Product(bdxtail, cdy, bdxt_cdy1, bdxt_cdy0);
2056  Two_Product(bdytail, cdx, bdyt_cdx1, bdyt_cdx0);
2057  Two_Two_Diff(bdxt_cdy1, bdxt_cdy0, bdyt_cdx1, bdyt_cdx0,
2058  bt_clarge, bt_c[2], bt_c[1], bt_c[0]);
2059  bt_c[3] = bt_clarge;
2060  bt_clen = 4;
2061  Two_Product(bdytail, adx, bdyt_adx1, bdyt_adx0);
2062  Two_Product(bdxtail, ady, bdxt_ady1, bdxt_ady0);
2063  Two_Two_Diff(bdyt_adx1, bdyt_adx0, bdxt_ady1, bdxt_ady0,
2064  bt_alarge, bt_a[2], bt_a[1], bt_a[0]);
2065  bt_a[3] = bt_alarge;
2066  bt_alen = 4;
2067  }
2068  }
2069  if (cdxtail == 0.0) {
2070  if (cdytail == 0.0) {
2071  ct_a[0] = 0.0;
2072  ct_alen = 1;
2073  ct_b[0] = 0.0;
2074  ct_blen = 1;
2075  } else {
2076  negate = -cdytail;
2077  Two_Product(negate, adx, ct_alarge, ct_a[0]);
2078  ct_a[1] = ct_alarge;
2079  ct_alen = 2;
2080  Two_Product(cdytail, bdx, ct_blarge, ct_b[0]);
2081  ct_b[1] = ct_blarge;
2082  ct_blen = 2;
2083  }
2084  } else {
2085  if (cdytail == 0.0) {
2086  Two_Product(cdxtail, ady, ct_alarge, ct_a[0]);
2087  ct_a[1] = ct_alarge;
2088  ct_alen = 2;
2089  negate = -cdxtail;
2090  Two_Product(negate, bdy, ct_blarge, ct_b[0]);
2091  ct_b[1] = ct_blarge;
2092  ct_blen = 2;
2093  } else {
2094  Two_Product(cdxtail, ady, cdxt_ady1, cdxt_ady0);
2095  Two_Product(cdytail, adx, cdyt_adx1, cdyt_adx0);
2096  Two_Two_Diff(cdxt_ady1, cdxt_ady0, cdyt_adx1, cdyt_adx0,
2097  ct_alarge, ct_a[2], ct_a[1], ct_a[0]);
2098  ct_a[3] = ct_alarge;
2099  ct_alen = 4;
2100  Two_Product(cdytail, bdx, cdyt_bdx1, cdyt_bdx0);
2101  Two_Product(cdxtail, bdy, cdxt_bdy1, cdxt_bdy0);
2102  Two_Two_Diff(cdyt_bdx1, cdyt_bdx0, cdxt_bdy1, cdxt_bdy0,
2103  ct_blarge, ct_b[2], ct_b[1], ct_b[0]);
2104  ct_b[3] = ct_blarge;
2105  ct_blen = 4;
2106  }
2107  }
2108 
2109  bctlen = fast_expansion_sum_zeroelim(bt_clen, bt_c, ct_blen, ct_b, bct);
2110  wlength = scale_expansion_zeroelim(bctlen, bct, adz, w);
2111  finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
2112  finother);
2113  finswap = finnow; finnow = finother; finother = finswap;
2114 
2115  catlen = fast_expansion_sum_zeroelim(ct_alen, ct_a, at_clen, at_c, cat);
2116  wlength = scale_expansion_zeroelim(catlen, cat, bdz, w);
2117  finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
2118  finother);
2119  finswap = finnow; finnow = finother; finother = finswap;
2120 
2121  abtlen = fast_expansion_sum_zeroelim(at_blen, at_b, bt_alen, bt_a, abt);
2122  wlength = scale_expansion_zeroelim(abtlen, abt, cdz, w);
2123  finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
2124  finother);
2125  finswap = finnow; finnow = finother; finother = finswap;
2126 
2127  if (adztail != 0.0) {
2128  vlength = scale_expansion_zeroelim(4, bc, adztail, v);
2129  finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
2130  finother);
2131  finswap = finnow; finnow = finother; finother = finswap;
2132  }
2133  if (bdztail != 0.0) {
2134  vlength = scale_expansion_zeroelim(4, ca, bdztail, v);
2135  finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
2136  finother);
2137  finswap = finnow; finnow = finother; finother = finswap;
2138  }
2139  if (cdztail != 0.0) {
2140  vlength = scale_expansion_zeroelim(4, ab, cdztail, v);
2141  finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
2142  finother);
2143  finswap = finnow; finnow = finother; finother = finswap;
2144  }
2145 
2146  if (adxtail != 0.0) {
2147  if (bdytail != 0.0) {
2148  Two_Product(adxtail, bdytail, adxt_bdyt1, adxt_bdyt0);
2149  Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdz, u3, u[2], u[1], u[0]);
2150  u[3] = u3;
2151  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2152  finother);
2153  finswap = finnow; finnow = finother; finother = finswap;
2154  if (cdztail != 0.0) {
2155  Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdztail, u3, u[2], u[1], u[0]);
2156  u[3] = u3;
2157  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2158  finother);
2159  finswap = finnow; finnow = finother; finother = finswap;
2160  }
2161  }
2162  if (cdytail != 0.0) {
2163  negate = -adxtail;
2164  Two_Product(negate, cdytail, adxt_cdyt1, adxt_cdyt0);
2165  Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdz, u3, u[2], u[1], u[0]);
2166  u[3] = u3;
2167  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2168  finother);
2169  finswap = finnow; finnow = finother; finother = finswap;
2170  if (bdztail != 0.0) {
2171  Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdztail, u3, u[2], u[1], u[0]);
2172  u[3] = u3;
2173  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2174  finother);
2175  finswap = finnow; finnow = finother; finother = finswap;
2176  }
2177  }
2178  }
2179  if (bdxtail != 0.0) {
2180  if (cdytail != 0.0) {
2181  Two_Product(bdxtail, cdytail, bdxt_cdyt1, bdxt_cdyt0);
2182  Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adz, u3, u[2], u[1], u[0]);
2183  u[3] = u3;
2184  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2185  finother);
2186  finswap = finnow; finnow = finother; finother = finswap;
2187  if (adztail != 0.0) {
2188  Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adztail, u3, u[2], u[1], u[0]);
2189  u[3] = u3;
2190  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2191  finother);
2192  finswap = finnow; finnow = finother; finother = finswap;
2193  }
2194  }
2195  if (adytail != 0.0) {
2196  negate = -bdxtail;
2197  Two_Product(negate, adytail, bdxt_adyt1, bdxt_adyt0);
2198  Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdz, u3, u[2], u[1], u[0]);
2199  u[3] = u3;
2200  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2201  finother);
2202  finswap = finnow; finnow = finother; finother = finswap;
2203  if (cdztail != 0.0) {
2204  Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdztail, u3, u[2], u[1], u[0]);
2205  u[3] = u3;
2206  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2207  finother);
2208  finswap = finnow; finnow = finother; finother = finswap;
2209  }
2210  }
2211  }
2212  if (cdxtail != 0.0) {
2213  if (adytail != 0.0) {
2214  Two_Product(cdxtail, adytail, cdxt_adyt1, cdxt_adyt0);
2215  Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdz, u3, u[2], u[1], u[0]);
2216  u[3] = u3;
2217  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2218  finother);
2219  finswap = finnow; finnow = finother; finother = finswap;
2220  if (bdztail != 0.0) {
2221  Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdztail, u3, u[2], u[1], u[0]);
2222  u[3] = u3;
2223  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2224  finother);
2225  finswap = finnow; finnow = finother; finother = finswap;
2226  }
2227  }
2228  if (bdytail != 0.0) {
2229  negate = -cdxtail;
2230  Two_Product(negate, bdytail, cdxt_bdyt1, cdxt_bdyt0);
2231  Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adz, u3, u[2], u[1], u[0]);
2232  u[3] = u3;
2233  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2234  finother);
2235  finswap = finnow; finnow = finother; finother = finswap;
2236  if (adztail != 0.0) {
2237  Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adztail, u3, u[2], u[1], u[0]);
2238  u[3] = u3;
2239  finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2240  finother);
2241  finswap = finnow; finnow = finother; finother = finswap;
2242  }
2243  }
2244  }
2245 
2246  if (adztail != 0.0) {
2247  wlength = scale_expansion_zeroelim(bctlen, bct, adztail, w);
2248  finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
2249  finother);
2250  finswap = finnow; finnow = finother; finother = finswap;
2251  }
2252  if (bdztail != 0.0) {
2253  wlength = scale_expansion_zeroelim(catlen, cat, bdztail, w);
2254  finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
2255  finother);
2256  finswap = finnow; finnow = finother; finother = finswap;
2257  }
2258  if (cdztail != 0.0) {
2259  wlength = scale_expansion_zeroelim(abtlen, abt, cdztail, w);
2260  finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
2261  finother);
2262  finswap = finnow; finnow = finother; finother = finswap;
2263  }
2264 
2265  return finnow[finlength - 1];
2266 }
2267 
2268 REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
2269 {
2270  REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
2271  REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
2272  REAL det;
2273  REAL permanent, errbound;
2274 
2275  adx = pa[0] - pd[0];
2276  bdx = pb[0] - pd[0];
2277  cdx = pc[0] - pd[0];
2278  ady = pa[1] - pd[1];
2279  bdy = pb[1] - pd[1];
2280  cdy = pc[1] - pd[1];
2281  adz = pa[2] - pd[2];
2282  bdz = pb[2] - pd[2];
2283  cdz = pc[2] - pd[2];
2284 
2285  bdxcdy = bdx * cdy;
2286  cdxbdy = cdx * bdy;
2287 
2288  cdxady = cdx * ady;
2289  adxcdy = adx * cdy;
2290 
2291  adxbdy = adx * bdy;
2292  bdxady = bdx * ady;
2293 
2294  det = adz * (bdxcdy - cdxbdy)
2295  + bdz * (cdxady - adxcdy)
2296  + cdz * (adxbdy - bdxady);
2297 
2298  permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * Absolute(adz)
2299  + (Absolute(cdxady) + Absolute(adxcdy)) * Absolute(bdz)
2300  + (Absolute(adxbdy) + Absolute(bdxady)) * Absolute(cdz);
2301  errbound = o3derrboundA * permanent;
2302  if ((det > errbound) || (-det > errbound)) {
2303  return det;
2304  }
2305 
2306  return orient3dadapt(pa, pb, pc, pd, permanent);
2307 }
2308 
2309 /*****************************************************************************/
2310 /* */
2311 /* incirclefast() Approximate 2D incircle test. Nonrobust. */
2312 /* incircleexact() Exact 2D incircle test. Robust. */
2313 /* incircleslow() Another exact 2D incircle test. Robust. */
2314 /* incircle() Adaptive exact 2D incircle test. Robust. */
2315 /* */
2316 /* Return a positive value if the point pd lies inside the */
2317 /* circle passing through pa, pb, and pc; a negative value if */
2318 /* it lies outside; and zero if the four points are cocircular.*/
2319 /* The points pa, pb, and pc must be in counterclockwise */
2320 /* order, or the sign of the result will be reversed. */
2321 /* */
2322 /* Only the first and last routine should be used; the middle two are for */
2323 /* timings. */
2324 /* */
2325 /* The last three use exact arithmetic to ensure a correct answer. The */
2326 /* result returned is the determinant of a matrix. In incircle() only, */
2327 /* this determinant is computed adaptively, in the sense that exact */
2328 /* arithmetic is used only to the degree it is needed to ensure that the */
2329 /* returned value has the correct sign. Hence, incircle() is usually quite */
2330 /* fast, but will run more slowly when the input points are cocircular or */
2331 /* nearly so. */
2332 /* */
2333 /*****************************************************************************/
2334 
2335 REAL incirclefast(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
2336 {
2337  REAL adx, ady, bdx, bdy, cdx, cdy;
2338  REAL abdet, bcdet, cadet;
2339  REAL alift, blift, clift;
2340 
2341  adx = pa[0] - pd[0];
2342  ady = pa[1] - pd[1];
2343  bdx = pb[0] - pd[0];
2344  bdy = pb[1] - pd[1];
2345  cdx = pc[0] - pd[0];
2346  cdy = pc[1] - pd[1];
2347 
2348  abdet = adx * bdy - bdx * ady;
2349  bcdet = bdx * cdy - cdx * bdy;
2350  cadet = cdx * ady - adx * cdy;
2351  alift = adx * adx + ady * ady;
2352  blift = bdx * bdx + bdy * bdy;
2353  clift = cdx * cdx + cdy * cdy;
2354 
2355  return alift * bcdet + blift * cadet + clift * abdet;
2356 }
2357 
2358 REAL incircleexact(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
2359 {
2360  INEXACT REAL axby1, bxcy1, cxdy1, dxay1, axcy1, bxdy1;
2361  INEXACT REAL bxay1, cxby1, dxcy1, axdy1, cxay1, dxby1;
2362  REAL axby0, bxcy0, cxdy0, dxay0, axcy0, bxdy0;
2363  REAL bxay0, cxby0, dxcy0, axdy0, cxay0, dxby0;
2364  REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4];
2365  REAL temp8[8];
2366  int templen;
2367  REAL abc[12], bcd[12], cda[12], dab[12];
2368  int abclen, bcdlen, cdalen, dablen;
2369  REAL det24x[24], det24y[24], det48x[48], det48y[48];
2370  int xlen, ylen;
2371  REAL adet[96], bdet[96], cdet[96], ddet[96];
2372  int alen, blen, clen, dlen;
2373  REAL abdet[192], cddet[192];
2374  int ablen, cdlen;
2375  REAL deter[384];
2376  int deterlen;
2377  int i;
2378 
2379  INEXACT REAL bvirt;
2380  REAL avirt, bround, around;
2381  INEXACT REAL c;
2382  INEXACT REAL abig;
2383  REAL ahi, alo, bhi, blo;
2384  REAL err1, err2, err3;
2385  INEXACT REAL _i, _j;
2386  REAL _0;
2387 
2388  Two_Product(pa[0], pb[1], axby1, axby0);
2389  Two_Product(pb[0], pa[1], bxay1, bxay0);
2390  Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]);
2391 
2392  Two_Product(pb[0], pc[1], bxcy1, bxcy0);
2393  Two_Product(pc[0], pb[1], cxby1, cxby0);
2394  Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]);
2395 
2396  Two_Product(pc[0], pd[1], cxdy1, cxdy0);
2397  Two_Product(pd[0], pc[1], dxcy1, dxcy0);
2398  Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]);
2399 
2400  Two_Product(pd[0], pa[1], dxay1, dxay0);
2401  Two_Product(pa[0], pd[1], axdy1, axdy0);
2402  Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]);
2403 
2404  Two_Product(pa[0], pc[1], axcy1, axcy0);
2405  Two_Product(pc[0], pa[1], cxay1, cxay0);
2406  Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]);
2407 
2408  Two_Product(pb[0], pd[1], bxdy1, bxdy0);
2409  Two_Product(pd[0], pb[1], dxby1, dxby0);
2410  Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]);
2411 
2412  templen = fast_expansion_sum_zeroelim(4, cd, 4, da, temp8);
2413  cdalen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, cda);
2414  templen = fast_expansion_sum_zeroelim(4, da, 4, ab, temp8);
2415  dablen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, dab);
2416  for (i = 0; i < 4; i++) {
2417  bd[i] = -bd[i];
2418  ac[i] = -ac[i];
2419  }
2420  templen = fast_expansion_sum_zeroelim(4, ab, 4, bc, temp8);
2421  abclen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, abc);
2422  templen = fast_expansion_sum_zeroelim(4, bc, 4, cd, temp8);
2423  bcdlen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, bcd);
2424 
2425  xlen = scale_expansion_zeroelim(bcdlen, bcd, pa[0], det24x);
2426  xlen = scale_expansion_zeroelim(xlen, det24x, pa[0], det48x);
2427  ylen = scale_expansion_zeroelim(bcdlen, bcd, pa[1], det24y);
2428  ylen = scale_expansion_zeroelim(ylen, det24y, pa[1], det48y);
2429  alen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, adet);
2430 
2431  xlen = scale_expansion_zeroelim(cdalen, cda, pb[0], det24x);
2432  xlen = scale_expansion_zeroelim(xlen, det24x, -pb[0], det48x);
2433  ylen = scale_expansion_zeroelim(cdalen, cda, pb[1], det24y);
2434  ylen = scale_expansion_zeroelim(ylen, det24y, -pb[1], det48y);
2435  blen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, bdet);
2436 
2437  xlen = scale_expansion_zeroelim(dablen, dab, pc[0], det24x);
2438  xlen = scale_expansion_zeroelim(xlen, det24x, pc[0], det48x);
2439  ylen = scale_expansion_zeroelim(dablen, dab, pc[1], det24y);
2440  ylen = scale_expansion_zeroelim(ylen, det24y, pc[1], det48y);
2441  clen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, cdet);
2442 
2443  xlen = scale_expansion_zeroelim(abclen, abc, pd[0], det24x);
2444  xlen = scale_expansion_zeroelim(xlen, det24x, -pd[0], det48x);
2445  ylen = scale_expansion_zeroelim(abclen, abc, pd[1], det24y);
2446  ylen = scale_expansion_zeroelim(ylen, det24y, -pd[1], det48y);
2447  dlen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, ddet);
2448 
2449  ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
2450  cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
2451  deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, deter);
2452 
2453  return deter[deterlen - 1];
2454 }
2455 
2456 REAL incircleslow(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
2457 {
2458  INEXACT REAL adx, bdx, cdx, ady, bdy, cdy;
2459  REAL adxtail, bdxtail, cdxtail;
2460  REAL adytail, bdytail, cdytail;
2461  REAL negate, negatetail;
2462  INEXACT REAL axby7, bxcy7, axcy7, bxay7, cxby7, cxay7;
2463  REAL axby[8], bxcy[8], axcy[8], bxay[8], cxby[8], cxay[8];
2464  REAL temp16[16];
2465  int temp16len;
2466  REAL detx[32], detxx[64], detxt[32], detxxt[64], detxtxt[64];
2467  int xlen, xxlen, xtlen, xxtlen, xtxtlen;
2468  REAL x1[128], x2[192];
2469  int x1len, x2len;
2470  REAL dety[32], detyy[64], detyt[32], detyyt[64], detytyt[64];
2471  int ylen, yylen, ytlen, yytlen, ytytlen;
2472  REAL y1[128], y2[192];
2473  int y1len, y2len;
2474  REAL adet[384], bdet[384], cdet[384], abdet[768], deter[1152];
2475  int alen, blen, clen, ablen, deterlen;
2476  int i;
2477 
2478  INEXACT REAL bvirt;
2479  REAL avirt, bround, around;
2480  INEXACT REAL c;
2481  INEXACT REAL abig;
2482  REAL a0hi, a0lo, a1hi, a1lo, bhi, blo;
2483  REAL err1, err2, err3;
2484  INEXACT REAL _i, _j, _k, _l, _m, _n;
2485  REAL _0, _1, _2;
2486 
2487  Two_Diff(pa[0], pd[0], adx, adxtail);
2488  Two_Diff(pa[1], pd[1], ady, adytail);
2489  Two_Diff(pb[0], pd[0], bdx, bdxtail);
2490  Two_Diff(pb[1], pd[1], bdy, bdytail);
2491  Two_Diff(pc[0], pd[0], cdx, cdxtail);
2492  Two_Diff(pc[1], pd[1], cdy, cdytail);
2493 
2494  Two_Two_Product(adx, adxtail, bdy, bdytail,
2495  axby7, axby[6], axby[5], axby[4],
2496  axby[3], axby[2], axby[1], axby[0]);
2497  axby[7] = axby7;
2498  negate = -ady;
2499  negatetail = -adytail;
2500  Two_Two_Product(bdx, bdxtail, negate, negatetail,
2501  bxay7, bxay[6], bxay[5], bxay[4],
2502  bxay[3], bxay[2], bxay[1], bxay[0]);
2503  bxay[7] = bxay7;
2504  Two_Two_Product(bdx, bdxtail, cdy, cdytail,
2505  bxcy7, bxcy[6], bxcy[5], bxcy[4],
2506  bxcy[3], bxcy[2], bxcy[1], bxcy[0]);
2507  bxcy[7] = bxcy7;
2508  negate = -bdy;
2509  negatetail = -bdytail;
2510  Two_Two_Product(cdx, cdxtail, negate, negatetail,
2511  cxby7, cxby[6], cxby[5], cxby[4],
2512  cxby[3], cxby[2], cxby[1], cxby[0]);
2513  cxby[7] = cxby7;
2514  Two_Two_Product(cdx, cdxtail, ady, adytail,
2515  cxay7, cxay[6], cxay[5], cxay[4],
2516  cxay[3], cxay[2], cxay[1], cxay[0]);
2517  cxay[7] = cxay7;
2518  negate = -cdy;
2519  negatetail = -cdytail;
2520  Two_Two_Product(adx, adxtail, negate, negatetail,
2521  axcy7, axcy[6], axcy[5], axcy[4],
2522  axcy[3], axcy[2], axcy[1], axcy[0]);
2523  axcy[7] = axcy7;
2524 
2525 
2526  temp16len = fast_expansion_sum_zeroelim(8, bxcy, 8, cxby, temp16);
2527 
2528  xlen = scale_expansion_zeroelim(temp16len, temp16, adx, detx);
2529  xxlen = scale_expansion_zeroelim(xlen, detx, adx, detxx);
2530  xtlen = scale_expansion_zeroelim(temp16len, temp16, adxtail, detxt);
2531  xxtlen = scale_expansion_zeroelim(xtlen, detxt, adx, detxxt);
2532  for (i = 0; i < xxtlen; i++) {
2533  detxxt[i] *= 2.0;
2534  }
2535  xtxtlen = scale_expansion_zeroelim(xtlen, detxt, adxtail, detxtxt);
2536  x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
2537  x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
2538 
2539  ylen = scale_expansion_zeroelim(temp16len, temp16, ady, dety);
2540  yylen = scale_expansion_zeroelim(ylen, dety, ady, detyy);
2541  ytlen = scale_expansion_zeroelim(temp16len, temp16, adytail, detyt);
2542  yytlen = scale_expansion_zeroelim(ytlen, detyt, ady, detyyt);
2543  for (i = 0; i < yytlen; i++) {
2544  detyyt[i] *= 2.0;
2545  }
2546  ytytlen = scale_expansion_zeroelim(ytlen, detyt, adytail, detytyt);
2547  y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
2548  y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
2549 
2550  alen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, adet);
2551 
2552 
2553  temp16len = fast_expansion_sum_zeroelim(8, cxay, 8, axcy, temp16);
2554 
2555  xlen = scale_expansion_zeroelim(temp16len, temp16, bdx, detx);
2556  xxlen = scale_expansion_zeroelim(xlen, detx, bdx, detxx);
2557  xtlen = scale_expansion_zeroelim(temp16len, temp16, bdxtail, detxt);
2558  xxtlen = scale_expansion_zeroelim(xtlen, detxt, bdx, detxxt);
2559  for (i = 0; i < xxtlen; i++) {
2560  detxxt[i] *= 2.0;
2561  }
2562  xtxtlen = scale_expansion_zeroelim(xtlen, detxt, bdxtail, detxtxt);
2563  x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
2564  x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
2565 
2566  ylen = scale_expansion_zeroelim(temp16len, temp16, bdy, dety);
2567  yylen = scale_expansion_zeroelim(ylen, dety, bdy, detyy);
2568  ytlen = scale_expansion_zeroelim(temp16len, temp16, bdytail, detyt);
2569  yytlen = scale_expansion_zeroelim(ytlen, detyt, bdy, detyyt);
2570  for (i = 0; i < yytlen; i++) {
2571  detyyt[i] *= 2.0;
2572  }
2573  ytytlen = scale_expansion_zeroelim(ytlen, detyt, bdytail, detytyt);
2574  y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
2575  y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
2576 
2577  blen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, bdet);
2578 
2579 
2580  temp16len = fast_expansion_sum_zeroelim(8, axby, 8, bxay, temp16);
2581 
2582  xlen = scale_expansion_zeroelim(temp16len, temp16, cdx, detx);
2583  xxlen = scale_expansion_zeroelim(xlen, detx, cdx, detxx);
2584  xtlen = scale_expansion_zeroelim(temp16len, temp16, cdxtail, detxt);
2585  xxtlen = scale_expansion_zeroelim(xtlen, detxt, cdx, detxxt);
2586  for (i = 0; i < xxtlen; i++) {
2587  detxxt[i] *= 2.0;
2588  }
2589  xtxtlen = scale_expansion_zeroelim(xtlen, detxt, cdxtail, detxtxt);
2590  x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
2591  x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
2592 
2593  ylen = scale_expansion_zeroelim(temp16len, temp16, cdy, dety);
2594  yylen = scale_expansion_zeroelim(ylen, dety, cdy, detyy);
2595  ytlen = scale_expansion_zeroelim(temp16len, temp16, cdytail, detyt);
2596  yytlen = scale_expansion_zeroelim(ytlen, detyt, cdy, detyyt);
2597  for (i = 0; i < yytlen; i++) {
2598  detyyt[i] *= 2.0;
2599  }
2600  ytytlen = scale_expansion_zeroelim(ytlen, detyt, cdytail, detytyt);
2601  y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
2602  y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
2603 
2604  clen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, cdet);
2605 
2606  ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
2607  deterlen = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, deter);
2608 
2609  return deter[deterlen - 1];
2610 }
2611 
2612 REAL incircleadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL permanent)
2613 {
2614  INEXACT REAL adx, bdx, cdx, ady, bdy, cdy;
2615  REAL det, errbound;
2616 
2617  INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
2618  REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
2619  REAL bc[4], ca[4], ab[4];
2620  INEXACT REAL bc3, ca3, ab3;
2621  REAL axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32];
2622  int axbclen, axxbclen, aybclen, ayybclen, alen;
2623  REAL bxca[8], bxxca[16], byca[8], byyca[16], bdet[32];
2624  int bxcalen, bxxcalen, bycalen, byycalen, blen;
2625  REAL cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32];
2626  int cxablen, cxxablen, cyablen, cyyablen, clen;
2627  REAL abdet[64];
2628  int ablen;
2629  REAL fin1[1152], fin2[1152];
2630  REAL *finnow, *finother, *finswap;
2631  int finlength;
2632 
2633  REAL adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
2634  INEXACT REAL adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1;
2635  REAL adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0;
2636  REAL aa[4], bb[4], cc[4];
2637  INEXACT REAL aa3, bb3, cc3;
2638  INEXACT REAL ti1, tj1;
2639  REAL ti0, tj0;
2640  REAL u[4], v[4];
2641  INEXACT REAL u3, v3;
2642  REAL temp8[8], temp16a[16], temp16b[16], temp16c[16];
2643  REAL temp32a[32], temp32b[32], temp48[48], temp64[64];
2644  int temp8len, temp16alen, temp16blen, temp16clen;
2645  int temp32alen, temp32blen, temp48len, temp64len;
2646  REAL axtbb[8], axtcc[8], aytbb[8], aytcc[8];
2647  int axtbblen, axtcclen, aytbblen, aytcclen;
2648  REAL bxtaa[8], bxtcc[8], bytaa[8], bytcc[8];
2649  int bxtaalen, bxtcclen, bytaalen, bytcclen;
2650  REAL cxtaa[8], cxtbb[8], cytaa[8], cytbb[8];
2651  int cxtaalen, cxtbblen, cytaalen, cytbblen;
2652  REAL axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8];
2653  int axtbclen, aytbclen, bxtcalen, bytcalen, cxtablen, cytablen;
2654  REAL axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16];
2655  int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen;
2656  REAL axtbctt[8], aytbctt[8], bxtcatt[8];
2657  REAL bytcatt[8], cxtabtt[8], cytabtt[8];
2658  int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen;
2659  REAL abt[8], bct[8], cat[8];
2660  int abtlen, bctlen, catlen;
2661  REAL abtt[4], bctt[4], catt[4];
2662  int abttlen, bcttlen, cattlen;
2663  INEXACT REAL abtt3, bctt3, catt3;
2664  REAL negate;
2665 
2666  INEXACT REAL bvirt;
2667  REAL avirt, bround, around;
2668  INEXACT REAL c;
2669  INEXACT REAL abig;
2670  REAL ahi, alo, bhi, blo;
2671  REAL err1, err2, err3;
2672  INEXACT REAL _i, _j;
2673  REAL _0;
2674 
2675  axtbclen=0;
2676  aytbclen=0;
2677  bxtcalen=0;
2678  bytcalen=0;
2679  cxtablen=0;
2680  cytablen=0;
2681  adx = (REAL) (pa[0] - pd[0]);
2682  bdx = (REAL) (pb[0] - pd[0]);
2683  cdx = (REAL) (pc[0] - pd[0]);
2684  ady = (REAL) (pa[1] - pd[1]);
2685  bdy = (REAL) (pb[1] - pd[1]);
2686  cdy = (REAL) (pc[1] - pd[1]);
2687 
2688  Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
2689  Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
2690  Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
2691  bc[3] = bc3;
2692  axbclen = scale_expansion_zeroelim(4, bc, adx, axbc);
2693  axxbclen = scale_expansion_zeroelim(axbclen, axbc, adx, axxbc);
2694  aybclen = scale_expansion_zeroelim(4, bc, ady, aybc);
2695  ayybclen = scale_expansion_zeroelim(aybclen, aybc, ady, ayybc);
2696  alen = fast_expansion_sum_zeroelim(axxbclen, axxbc, ayybclen, ayybc, adet);
2697 
2698  Two_Product(cdx, ady, cdxady1, cdxady0);
2699  Two_Product(adx, cdy, adxcdy1, adxcdy0);
2700  Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
2701  ca[3] = ca3;
2702  bxcalen = scale_expansion_zeroelim(4, ca, bdx, bxca);
2703  bxxcalen = scale_expansion_zeroelim(bxcalen, bxca, bdx, bxxca);
2704  bycalen = scale_expansion_zeroelim(4, ca, bdy, byca);
2705  byycalen = scale_expansion_zeroelim(bycalen, byca, bdy, byyca);
2706  blen = fast_expansion_sum_zeroelim(bxxcalen, bxxca, byycalen, byyca, bdet);
2707 
2708  Two_Product(adx, bdy, adxbdy1, adxbdy0);
2709  Two_Product(bdx, ady, bdxady1, bdxady0);
2710  Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
2711  ab[3] = ab3;
2712  cxablen = scale_expansion_zeroelim(4, ab, cdx, cxab);
2713  cxxablen = scale_expansion_zeroelim(cxablen, cxab, cdx, cxxab);
2714  cyablen = scale_expansion_zeroelim(4, ab, cdy, cyab);
2715  cyyablen = scale_expansion_zeroelim(cyablen, cyab, cdy, cyyab);
2716  clen = fast_expansion_sum_zeroelim(cxxablen, cxxab, cyyablen, cyyab, cdet);
2717 
2718  ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
2719  finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
2720 
2721  det = estimate(finlength, fin1);
2722  errbound = iccerrboundB * permanent;
2723  if ((det >= errbound) || (-det >= errbound)) {
2724  return det;
2725  }
2726 
2727  Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
2728  Two_Diff_Tail(pa[1], pd[1], ady, adytail);
2729  Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
2730  Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
2731  Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
2732  Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
2733  if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
2734  && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)) {
2735  return det;
2736  }
2737 
2738  errbound = iccerrboundC * permanent + resulterrbound * Absolute(det);
2739  det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail)
2740  - (bdy * cdxtail + cdx * bdytail))
2741  + 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx))
2742  + ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail)
2743  - (cdy * adxtail + adx * cdytail))
2744  + 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx))
2745  + ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail)
2746  - (ady * bdxtail + bdx * adytail))
2747  + 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));
2748  if ((det >= errbound) || (-det >= errbound)) {
2749  return det;
2750  }
2751 
2752  finnow = fin1;
2753  finother = fin2;
2754 
2755  if ((bdxtail != 0.0) || (bdytail != 0.0)
2756  || (cdxtail != 0.0) || (cdytail != 0.0)) {
2757  Square(adx, adxadx1, adxadx0);
2758  Square(ady, adyady1, adyady0);
2759  Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]);
2760  aa[3] = aa3;
2761  }
2762  if ((cdxtail != 0.0) || (cdytail != 0.0)
2763  || (adxtail != 0.0) || (adytail != 0.0)) {
2764  Square(bdx, bdxbdx1, bdxbdx0);
2765  Square(bdy, bdybdy1, bdybdy0);
2766  Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]);
2767  bb[3] = bb3;
2768  }
2769  if ((adxtail != 0.0) || (adytail != 0.0)
2770  || (bdxtail != 0.0) || (bdytail != 0.0)) {
2771  Square(cdx, cdxcdx1, cdxcdx0);
2772  Square(cdy, cdycdy1, cdycdy0);
2773  Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]);
2774  cc[3] = cc3;
2775  }
2776 
2777  if (adxtail != 0.0) {
2778  axtbclen = scale_expansion_zeroelim(4, bc, adxtail, axtbc);
2779  temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, 2.0 * adx,
2780  temp16a);
2781 
2782  axtcclen = scale_expansion_zeroelim(4, cc, adxtail, axtcc);
2783  temp16blen = scale_expansion_zeroelim(axtcclen, axtcc, bdy, temp16b);
2784 
2785  axtbblen = scale_expansion_zeroelim(4, bb, adxtail, axtbb);
2786  temp16clen = scale_expansion_zeroelim(axtbblen, axtbb, -cdy, temp16c);
2787 
2788  temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2789  temp16blen, temp16b, temp32a);
2790  temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
2791  temp32alen, temp32a, temp48);
2792  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2793  temp48, finother);
2794  finswap = finnow; finnow = finother; finother = finswap;
2795  }
2796  if (adytail != 0.0) {
2797  aytbclen = scale_expansion_zeroelim(4, bc, adytail, aytbc);
2798  temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, 2.0 * ady,
2799  temp16a);
2800 
2801  aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb);
2802  temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b);
2803 
2804  aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc);
2805  temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c);
2806 
2807  temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2808  temp16blen, temp16b, temp32a);
2809  temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
2810  temp32alen, temp32a, temp48);
2811  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2812  temp48, finother);
2813  finswap = finnow; finnow = finother; finother = finswap;
2814  }
2815  if (bdxtail != 0.0) {
2816  bxtcalen = scale_expansion_zeroelim(4, ca, bdxtail, bxtca);
2817  temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, 2.0 * bdx,
2818  temp16a);
2819 
2820  bxtaalen = scale_expansion_zeroelim(4, aa, bdxtail, bxtaa);
2821  temp16blen = scale_expansion_zeroelim(bxtaalen, bxtaa, cdy, temp16b);
2822 
2823  bxtcclen = scale_expansion_zeroelim(4, cc, bdxtail, bxtcc);
2824  temp16clen = scale_expansion_zeroelim(bxtcclen, bxtcc, -ady, temp16c);
2825 
2826  temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2827  temp16blen, temp16b, temp32a);
2828  temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
2829  temp32alen, temp32a, temp48);
2830  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2831  temp48, finother);
2832  finswap = finnow; finnow = finother; finother = finswap;
2833  }
2834  if (bdytail != 0.0) {
2835  bytcalen = scale_expansion_zeroelim(4, ca, bdytail, bytca);
2836  temp16alen = scale_expansion_zeroelim(bytcalen, bytca, 2.0 * bdy,
2837  temp16a);
2838 
2839  bytcclen = scale_expansion_zeroelim(4, cc, bdytail, bytcc);
2840  temp16blen = scale_expansion_zeroelim(bytcclen, bytcc, adx, temp16b);
2841 
2842  bytaalen = scale_expansion_zeroelim(4, aa, bdytail, bytaa);
2843  temp16clen = scale_expansion_zeroelim(bytaalen, bytaa, -cdx, temp16c);
2844 
2845  temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2846  temp16blen, temp16b, temp32a);
2847  temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
2848  temp32alen, temp32a, temp48);
2849  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2850  temp48, finother);
2851  finswap = finnow; finnow = finother; finother = finswap;
2852  }
2853  if (cdxtail != 0.0) {
2854  cxtablen = scale_expansion_zeroelim(4, ab, cdxtail, cxtab);
2855  temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, 2.0 * cdx,
2856  temp16a);
2857 
2858  cxtbblen = scale_expansion_zeroelim(4, bb, cdxtail, cxtbb);
2859  temp16blen = scale_expansion_zeroelim(cxtbblen, cxtbb, ady, temp16b);
2860 
2861  cxtaalen = scale_expansion_zeroelim(4, aa, cdxtail, cxtaa);
2862  temp16clen = scale_expansion_zeroelim(cxtaalen, cxtaa, -bdy, temp16c);
2863 
2864  temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2865  temp16blen, temp16b, temp32a);
2866  temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
2867  temp32alen, temp32a, temp48);
2868  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2869  temp48, finother);
2870  finswap = finnow; finnow = finother; finother = finswap;
2871  }
2872  if (cdytail != 0.0) {
2873  cytablen = scale_expansion_zeroelim(4, ab, cdytail, cytab);
2874  temp16alen = scale_expansion_zeroelim(cytablen, cytab, 2.0 * cdy,
2875  temp16a);
2876 
2877  cytaalen = scale_expansion_zeroelim(4, aa, cdytail, cytaa);
2878  temp16blen = scale_expansion_zeroelim(cytaalen, cytaa, bdx, temp16b);
2879 
2880  cytbblen = scale_expansion_zeroelim(4, bb, cdytail, cytbb);
2881  temp16clen = scale_expansion_zeroelim(cytbblen, cytbb, -adx, temp16c);
2882 
2883  temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2884  temp16blen, temp16b, temp32a);
2885  temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
2886  temp32alen, temp32a, temp48);
2887  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2888  temp48, finother);
2889  finswap = finnow; finnow = finother; finother = finswap;
2890  }
2891 
2892  if ((adxtail != 0.0) || (adytail != 0.0)) {
2893  if ((bdxtail != 0.0) || (bdytail != 0.0)
2894  || (cdxtail != 0.0) || (cdytail != 0.0)) {
2895  Two_Product(bdxtail, cdy, ti1, ti0);
2896  Two_Product(bdx, cdytail, tj1, tj0);
2897  Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
2898  u[3] = u3;
2899  negate = -bdy;
2900  Two_Product(cdxtail, negate, ti1, ti0);
2901  negate = -bdytail;
2902  Two_Product(cdx, negate, tj1, tj0);
2903  Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
2904  v[3] = v3;
2905  bctlen = fast_expansion_sum_zeroelim(4, u, 4, v, bct);
2906 
2907  Two_Product(bdxtail, cdytail, ti1, ti0);
2908  Two_Product(cdxtail, bdytail, tj1, tj0);
2909  Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]);
2910  bctt[3] = bctt3;
2911  bcttlen = 4;
2912  } else {
2913  bct[0] = 0.0;
2914  bctlen = 1;
2915  bctt[0] = 0.0;
2916  bcttlen = 1;
2917  }
2918 
2919  if (adxtail != 0.0) {
2920  temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, adxtail, temp16a);
2921  axtbctlen = scale_expansion_zeroelim(bctlen, bct, adxtail, axtbct);
2922  temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, 2.0 * adx,
2923  temp32a);
2924  temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2925  temp32alen, temp32a, temp48);
2926  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2927  temp48, finother);
2928  finswap = finnow; finnow = finother; finother = finswap;
2929  if (bdytail != 0.0) {
2930  temp8len = scale_expansion_zeroelim(4, cc, adxtail, temp8);
2931  temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
2932  temp16a);
2933  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
2934  temp16a, finother);
2935  finswap = finnow; finnow = finother; finother = finswap;
2936  }
2937  if (cdytail != 0.0) {
2938  temp8len = scale_expansion_zeroelim(4, bb, -adxtail, temp8);
2939  temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
2940  temp16a);
2941  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
2942  temp16a, finother);
2943  finswap = finnow; finnow = finother; finother = finswap;
2944  }
2945 
2946  temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, adxtail,
2947  temp32a);
2948  axtbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adxtail, axtbctt);
2949  temp16alen = scale_expansion_zeroelim(axtbcttlen, axtbctt, 2.0 * adx,
2950  temp16a);
2951  temp16blen = scale_expansion_zeroelim(axtbcttlen, axtbctt, adxtail,
2952  temp16b);
2953  temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2954  temp16blen, temp16b, temp32b);
2955  temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
2956  temp32blen, temp32b, temp64);
2957  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
2958  temp64, finother);
2959  finswap = finnow; finnow = finother; finother = finswap;
2960  }
2961  if (adytail != 0.0) {
2962  temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, adytail, temp16a);
2963  aytbctlen = scale_expansion_zeroelim(bctlen, bct, adytail, aytbct);
2964  temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, 2.0 * ady,
2965  temp32a);
2966  temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2967  temp32alen, temp32a, temp48);
2968  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2969  temp48, finother);
2970  finswap = finnow; finnow = finother; finother = finswap;
2971 
2972 
2973  temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, adytail,
2974  temp32a);
2975  aytbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adytail, aytbctt);
2976  temp16alen = scale_expansion_zeroelim(aytbcttlen, aytbctt, 2.0 * ady,
2977  temp16a);
2978  temp16blen = scale_expansion_zeroelim(aytbcttlen, aytbctt, adytail,
2979  temp16b);
2980  temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2981  temp16blen, temp16b, temp32b);
2982  temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
2983  temp32blen, temp32b, temp64);
2984  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
2985  temp64, finother);
2986  finswap = finnow; finnow = finother; finother = finswap;
2987  }
2988  }
2989  if ((bdxtail != 0.0) || (bdytail != 0.0)) {
2990  if ((cdxtail != 0.0) || (cdytail != 0.0)
2991  || (adxtail != 0.0) || (adytail != 0.0)) {
2992  Two_Product(cdxtail, ady, ti1, ti0);
2993  Two_Product(cdx, adytail, tj1, tj0);
2994  Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
2995  u[3] = u3;
2996  negate = -cdy;
2997  Two_Product(adxtail, negate, ti1, ti0);
2998  negate = -cdytail;
2999  Two_Product(adx, negate, tj1, tj0);
3000  Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
3001  v[3] = v3;
3002  catlen = fast_expansion_sum_zeroelim(4, u, 4, v, cat);
3003 
3004  Two_Product(cdxtail, adytail, ti1, ti0);
3005  Two_Product(adxtail, cdytail, tj1, tj0);
3006  Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]);
3007  catt[3] = catt3;
3008  cattlen = 4;
3009  } else {
3010  cat[0] = 0.0;
3011  catlen = 1;
3012  catt[0] = 0.0;
3013  cattlen = 1;
3014  }
3015 
3016  if (bdxtail != 0.0) {
3017  temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, bdxtail, temp16a);
3018  bxtcatlen = scale_expansion_zeroelim(catlen, cat, bdxtail, bxtcat);
3019  temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, 2.0 * bdx,
3020  temp32a);
3021  temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3022  temp32alen, temp32a, temp48);
3023  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
3024  temp48, finother);
3025  finswap = finnow; finnow = finother; finother = finswap;
3026  if (cdytail != 0.0) {
3027  temp8len = scale_expansion_zeroelim(4, aa, bdxtail, temp8);
3028  temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
3029  temp16a);
3030  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
3031  temp16a, finother);
3032  finswap = finnow; finnow = finother; finother = finswap;
3033  }
3034  if (adytail != 0.0) {
3035  temp8len = scale_expansion_zeroelim(4, cc, -bdxtail, temp8);
3036  temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
3037  temp16a);
3038  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
3039  temp16a, finother);
3040  finswap = finnow; finnow = finother; finother = finswap;
3041  }
3042 
3043  temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, bdxtail,
3044  temp32a);
3045  bxtcattlen = scale_expansion_zeroelim(cattlen, catt, bdxtail, bxtcatt);
3046  temp16alen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, 2.0 * bdx,
3047  temp16a);
3048  temp16blen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, bdxtail,
3049  temp16b);
3050  temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3051  temp16blen, temp16b, temp32b);
3052  temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3053  temp32blen, temp32b, temp64);
3054  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
3055  temp64, finother);
3056  finswap = finnow; finnow = finother; finother = finswap;
3057  }
3058  if (bdytail != 0.0) {
3059  temp16alen = scale_expansion_zeroelim(bytcalen, bytca, bdytail, temp16a);
3060  bytcatlen = scale_expansion_zeroelim(catlen, cat, bdytail, bytcat);
3061  temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, 2.0 * bdy,
3062  temp32a);
3063  temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3064  temp32alen, temp32a, temp48);
3065  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
3066  temp48, finother);
3067  finswap = finnow; finnow = finother; finother = finswap;
3068 
3069 
3070  temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, bdytail,
3071  temp32a);
3072  bytcattlen = scale_expansion_zeroelim(cattlen, catt, bdytail, bytcatt);
3073  temp16alen = scale_expansion_zeroelim(bytcattlen, bytcatt, 2.0 * bdy,
3074  temp16a);
3075  temp16blen = scale_expansion_zeroelim(bytcattlen, bytcatt, bdytail,
3076  temp16b);
3077  temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3078  temp16blen, temp16b, temp32b);
3079  temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3080  temp32blen, temp32b, temp64);
3081  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
3082  temp64, finother);
3083  finswap = finnow; finnow = finother; finother = finswap;
3084  }
3085  }
3086  if ((cdxtail != 0.0) || (cdytail != 0.0)) {
3087  if ((adxtail != 0.0) || (adytail != 0.0)
3088  || (bdxtail != 0.0) || (bdytail != 0.0)) {
3089  Two_Product(adxtail, bdy, ti1, ti0);
3090  Two_Product(adx, bdytail, tj1, tj0);
3091  Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
3092  u[3] = u3;
3093  negate = -ady;
3094  Two_Product(bdxtail, negate, ti1, ti0);
3095  negate = -adytail;
3096  Two_Product(bdx, negate, tj1, tj0);
3097  Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
3098  v[3] = v3;
3099  abtlen = fast_expansion_sum_zeroelim(4, u, 4, v, abt);
3100 
3101  Two_Product(adxtail, bdytail, ti1, ti0);
3102  Two_Product(bdxtail, adytail, tj1, tj0);
3103  Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]);
3104  abtt[3] = abtt3;
3105  abttlen = 4;
3106  } else {
3107  abt[0] = 0.0;
3108  abtlen = 1;
3109  abtt[0] = 0.0;
3110  abttlen = 1;
3111  }
3112 
3113  if (cdxtail != 0.0) {
3114  temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, cdxtail, temp16a);
3115  cxtabtlen = scale_expansion_zeroelim(abtlen, abt, cdxtail, cxtabt);
3116  temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, 2.0 * cdx,
3117  temp32a);
3118  temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3119  temp32alen, temp32a, temp48);
3120  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
3121  temp48, finother);
3122  finswap = finnow; finnow = finother; finother = finswap;
3123  if (adytail != 0.0) {
3124  temp8len = scale_expansion_zeroelim(4, bb, cdxtail, temp8);
3125  temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
3126  temp16a);
3127  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
3128  temp16a, finother);
3129  finswap = finnow; finnow = finother; finother = finswap;
3130  }
3131  if (bdytail != 0.0) {
3132  temp8len = scale_expansion_zeroelim(4, aa, -cdxtail, temp8);
3133  temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
3134  temp16a);
3135  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
3136  temp16a, finother);
3137  finswap = finnow; finnow = finother; finother = finswap;
3138  }
3139 
3140  temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, cdxtail,
3141  temp32a);
3142  cxtabttlen = scale_expansion_zeroelim(abttlen, abtt, cdxtail, cxtabtt);
3143  temp16alen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, 2.0 * cdx,
3144  temp16a);
3145  temp16blen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, cdxtail,
3146  temp16b);
3147  temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3148  temp16blen, temp16b, temp32b);
3149  temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3150  temp32blen, temp32b, temp64);
3151  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
3152  temp64, finother);
3153  finswap = finnow; finnow = finother; finother = finswap;
3154  }
3155  if (cdytail != 0.0) {
3156  temp16alen = scale_expansion_zeroelim(cytablen, cytab, cdytail, temp16a);
3157  cytabtlen = scale_expansion_zeroelim(abtlen, abt, cdytail, cytabt);
3158  temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, 2.0 * cdy,
3159  temp32a);
3160  temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3161  temp32alen, temp32a, temp48);
3162  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
3163  temp48, finother);
3164  finswap = finnow; finnow = finother; finother = finswap;
3165 
3166 
3167  temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, cdytail,
3168  temp32a);
3169  cytabttlen = scale_expansion_zeroelim(abttlen, abtt, cdytail, cytabtt);
3170  temp16alen = scale_expansion_zeroelim(cytabttlen, cytabtt, 2.0 * cdy,
3171  temp16a);
3172  temp16blen = scale_expansion_zeroelim(cytabttlen, cytabtt, cdytail,
3173  temp16b);
3174  temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3175  temp16blen, temp16b, temp32b);
3176  temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3177  temp32blen, temp32b, temp64);
3178  finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
3179  temp64, finother);
3180  finswap = finnow; finnow = finother; finother = finswap;
3181  }
3182  }
3183 
3184  return finnow[finlength - 1];
3185 }
3186 
3187 REAL incircle(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
3188 {
3189  REAL adx, bdx, cdx, ady, bdy, cdy;
3190  REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
3191  REAL alift, blift, clift;
3192  REAL det;
3193  REAL permanent, errbound;
3194 
3195  adx = pa[0] - pd[0];
3196  bdx = pb[0] - pd[0];
3197  cdx = pc[0] - pd[0];
3198  ady = pa[1] - pd[1];
3199  bdy = pb[1] - pd[1];
3200  cdy = pc[1] - pd[1];
3201 
3202  bdxcdy = bdx * cdy;
3203  cdxbdy = cdx * bdy;
3204  alift = adx * adx + ady * ady;
3205 
3206  cdxady = cdx * ady;
3207  adxcdy = adx * cdy;
3208  blift = bdx * bdx + bdy * bdy;
3209 
3210  adxbdy = adx * bdy;
3211  bdxady = bdx * ady;
3212  clift = cdx * cdx + cdy * cdy;
3213 
3214  det = alift * (bdxcdy - cdxbdy)
3215  + blift * (cdxady - adxcdy)
3216  + clift * (adxbdy - bdxady);
3217 
3218  permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift
3219  + (Absolute(cdxady) + Absolute(adxcdy)) * blift
3220  + (Absolute(adxbdy) + Absolute(bdxady)) * clift;
3221  errbound = iccerrboundA * permanent;
3222  if ((det > errbound) || (-det > errbound)) {
3223  return det;
3224  }
3225 
3226  return incircleadapt(pa, pb, pc, pd, permanent);
3227 }
3228 
3229 /*****************************************************************************/
3230 /* */
3231 /* inspherefast() Approximate 3D insphere test. Nonrobust. */
3232 /* insphereexact() Exact 3D insphere test. Robust. */
3233 /* insphereslow() Another exact 3D insphere test. Robust. */
3234 /* insphere() Adaptive exact 3D insphere test. Robust. */
3235 /* */
3236 /* Return a positive value if the point pe lies inside the */
3237 /* sphere passing through pa, pb, pc, and pd; a negative value */
3238 /* if it lies outside; and zero if the five points are */
3239 /* cospherical. The points pa, pb, pc, and pd must be ordered */
3240 /* so that they have a positive orientation (as defined by */
3241 /* orient3d()), or the sign of the result will be reversed. */
3242 /* */
3243 /* Only the first and last routine should be used; the middle two are for */
3244 /* timings. */
3245 /* */
3246 /* The last three use exact arithmetic to ensure a correct answer. The */
3247 /* result returned is the determinant of a matrix. In insphere() only, */
3248 /* this determinant is computed adaptively, in the sense that exact */
3249 /* arithmetic is used only to the degree it is needed to ensure that the */
3250 /* returned value has the correct sign. Hence, insphere() is usually quite */
3251 /* fast, but will run more slowly when the input points are cospherical or */
3252 /* nearly so. */
3253 /* */
3254 /*****************************************************************************/
3255 
3256 REAL inspherefast(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
3257 {
3258  REAL aex, bex, cex, dex;
3259  REAL aey, bey, cey, dey;
3260  REAL aez, bez, cez, dez;
3261  REAL alift, blift, clift, dlift;
3262  REAL ab, bc, cd, da, ac, bd;
3263  REAL abc, bcd, cda, dab;
3264 
3265  aex = pa[0] - pe[0];
3266  bex = pb[0] - pe[0];
3267  cex = pc[0] - pe[0];
3268  dex = pd[0] - pe[0];
3269  aey = pa[1] - pe[1];
3270  bey = pb[1] - pe[1];
3271  cey = pc[1] - pe[1];
3272  dey = pd[1] - pe[1];
3273  aez = pa[2] - pe[2];
3274  bez = pb[2] - pe[2];
3275  cez = pc[2] - pe[2];
3276  dez = pd[2] - pe[2];
3277 
3278  ab = aex * bey - bex * aey;
3279  bc = bex * cey - cex * bey;
3280  cd = cex * dey - dex * cey;
3281  da = dex * aey - aex * dey;
3282 
3283  ac = aex * cey - cex * aey;
3284  bd = bex * dey - dex * bey;
3285 
3286  abc = aez * bc - bez * ac + cez * ab;
3287  bcd = bez * cd - cez * bd + dez * bc;
3288  cda = cez * da + dez * ac + aez * cd;
3289  dab = dez * ab + aez * bd + bez * da;
3290 
3291  alift = aex * aex + aey * aey + aez * aez;
3292  blift = bex * bex + bey * bey + bez * bez;
3293  clift = cex * cex + cey * cey + cez * cez;
3294  dlift = dex * dex + dey * dey + dez * dez;
3295 
3296  return (dlift * abc - clift * dab) + (blift * cda - alift * bcd);
3297 }
3298 
3299 REAL insphereexact(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
3300 {
3301  INEXACT REAL axby1, bxcy1, cxdy1, dxey1, exay1;
3302  INEXACT REAL bxay1, cxby1, dxcy1, exdy1, axey1;
3303  INEXACT REAL axcy1, bxdy1, cxey1, dxay1, exby1;
3304  INEXACT REAL cxay1, dxby1, excy1, axdy1, bxey1;
3305  REAL axby0, bxcy0, cxdy0, dxey0, exay0;
3306  REAL bxay0, cxby0, dxcy0, exdy0, axey0;
3307  REAL axcy0, bxdy0, cxey0, dxay0, exby0;
3308  REAL cxay0, dxby0, excy0, axdy0, bxey0;
3309  REAL ab[4], bc[4], cd[4], de[4], ea[4];
3310  REAL ac[4], bd[4], ce[4], da[4], eb[4];
3311  REAL temp8a[8], temp8b[8], temp16[16];
3312  int temp8alen, temp8blen, temp16len;
3313  REAL abc[24], bcd[24], cde[24], dea[24], eab[24];
3314  REAL abd[24], bce[24], cda[24], deb[24], eac[24];
3315  int abclen, bcdlen, cdelen, dealen, eablen;
3316  int abdlen, bcelen, cdalen, deblen, eaclen;
3317  REAL temp48a[48], temp48b[48];
3318  int temp48alen, temp48blen;
3319  REAL abcd[96], bcde[96], cdea[96], deab[96], eabc[96];
3320  int abcdlen, bcdelen, cdealen, deablen, eabclen;
3321  REAL temp192[192];
3322  REAL det384x[384], det384y[384], det384z[384];
3323  int xlen, ylen, zlen;
3324  REAL detxy[768];
3325  int xylen;
3326  REAL adet[1152], bdet[1152], cdet[1152], ddet[1152], edet[1152];
3327  int alen, blen, clen, dlen, elen;
3328  REAL abdet[2304], cddet[2304], cdedet[3456];
3329  int ablen, cdlen;
3330  REAL deter[5760];
3331  int deterlen;
3332  int i;
3333 
3334  INEXACT REAL bvirt;
3335  REAL avirt, bround, around;
3336  INEXACT REAL c;
3337  INEXACT REAL abig;
3338  REAL ahi, alo, bhi, blo;
3339  REAL err1, err2, err3;
3340  INEXACT REAL _i, _j;
3341  REAL _0;
3342 
3343  Two_Product(pa[0], pb[1], axby1, axby0);
3344  Two_Product(pb[0], pa[1], bxay1, bxay0);
3345  Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]);
3346 
3347  Two_Product(pb[0], pc[1], bxcy1, bxcy0);
3348  Two_Product(pc[0], pb[1], cxby1, cxby0);
3349  Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]);
3350 
3351  Two_Product(pc[0], pd[1], cxdy1, cxdy0);
3352  Two_Product(pd[0], pc[1], dxcy1, dxcy0);
3353  Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]);
3354 
3355  Two_Product(pd[0], pe[1], dxey1, dxey0);
3356  Two_Product(pe[0], pd[1], exdy1, exdy0);
3357  Two_Two_Diff(dxey1, dxey0, exdy1, exdy0, de[3], de[2], de[1], de[0]);
3358 
3359  Two_Product(pe[0], pa[1], exay1, exay0);
3360  Two_Product(pa[0], pe[1], axey1, axey0);
3361  Two_Two_Diff(exay1, exay0, axey1, axey0, ea[3], ea[2], ea[1], ea[0]);
3362 
3363  Two_Product(pa[0], pc[1], axcy1, axcy0);
3364  Two_Product(pc[0], pa[1], cxay1, cxay0);
3365  Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]);
3366 
3367  Two_Product(pb[0], pd[1], bxdy1, bxdy0);
3368  Two_Product(pd[0], pb[1], dxby1, dxby0);
3369  Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]);
3370 
3371  Two_Product(pc[0], pe[1], cxey1, cxey0);
3372  Two_Product(pe[0], pc[1], excy1, excy0);
3373  Two_Two_Diff(cxey1, cxey0, excy1, excy0, ce[3], ce[2], ce[1], ce[0]);
3374 
3375  Two_Product(pd[0], pa[1], dxay1, dxay0);
3376  Two_Product(pa[0], pd[1], axdy1, axdy0);
3377  Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]);
3378 
3379  Two_Product(pe[0], pb[1], exby1, exby0);
3380  Two_Product(pb[0], pe[1], bxey1, bxey0);
3381  Two_Two_Diff(exby1, exby0, bxey1, bxey0, eb[3], eb[2], eb[1], eb[0]);
3382 
3383  temp8alen = scale_expansion_zeroelim(4, bc, pa[2], temp8a);
3384  temp8blen = scale_expansion_zeroelim(4, ac, -pb[2], temp8b);
3385  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3386  temp16);
3387  temp8alen = scale_expansion_zeroelim(4, ab, pc[2], temp8a);
3388  abclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3389  abc);
3390 
3391  temp8alen = scale_expansion_zeroelim(4, cd, pb[2], temp8a);
3392  temp8blen = scale_expansion_zeroelim(4, bd, -pc[2], temp8b);
3393  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3394  temp16);
3395  temp8alen = scale_expansion_zeroelim(4, bc, pd[2], temp8a);
3396  bcdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3397  bcd);
3398 
3399  temp8alen = scale_expansion_zeroelim(4, de, pc[2], temp8a);
3400  temp8blen = scale_expansion_zeroelim(4, ce, -pd[2], temp8b);
3401  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3402  temp16);
3403  temp8alen = scale_expansion_zeroelim(4, cd, pe[2], temp8a);
3404  cdelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3405  cde);
3406 
3407  temp8alen = scale_expansion_zeroelim(4, ea, pd[2], temp8a);
3408  temp8blen = scale_expansion_zeroelim(4, da, -pe[2], temp8b);
3409  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3410  temp16);
3411  temp8alen = scale_expansion_zeroelim(4, de, pa[2], temp8a);
3412  dealen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3413  dea);
3414 
3415  temp8alen = scale_expansion_zeroelim(4, ab, pe[2], temp8a);
3416  temp8blen = scale_expansion_zeroelim(4, eb, -pa[2], temp8b);
3417  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3418  temp16);
3419  temp8alen = scale_expansion_zeroelim(4, ea, pb[2], temp8a);
3420  eablen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3421  eab);
3422 
3423  temp8alen = scale_expansion_zeroelim(4, bd, pa[2], temp8a);
3424  temp8blen = scale_expansion_zeroelim(4, da, pb[2], temp8b);
3425  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3426  temp16);
3427  temp8alen = scale_expansion_zeroelim(4, ab, pd[2], temp8a);
3428  abdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3429  abd);
3430 
3431  temp8alen = scale_expansion_zeroelim(4, ce, pb[2], temp8a);
3432  temp8blen = scale_expansion_zeroelim(4, eb, pc[2], temp8b);
3433  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3434  temp16);
3435  temp8alen = scale_expansion_zeroelim(4, bc, pe[2], temp8a);
3436  bcelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3437  bce);
3438 
3439  temp8alen = scale_expansion_zeroelim(4, da, pc[2], temp8a);
3440  temp8blen = scale_expansion_zeroelim(4, ac, pd[2], temp8b);
3441  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3442  temp16);
3443  temp8alen = scale_expansion_zeroelim(4, cd, pa[2], temp8a);
3444  cdalen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3445  cda);
3446 
3447  temp8alen = scale_expansion_zeroelim(4, eb, pd[2], temp8a);
3448  temp8blen = scale_expansion_zeroelim(4, bd, pe[2], temp8b);
3449  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3450  temp16);
3451  temp8alen = scale_expansion_zeroelim(4, de, pb[2], temp8a);
3452  deblen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3453  deb);
3454 
3455  temp8alen = scale_expansion_zeroelim(4, ac, pe[2], temp8a);
3456  temp8blen = scale_expansion_zeroelim(4, ce, pa[2], temp8b);
3457  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3458  temp16);
3459  temp8alen = scale_expansion_zeroelim(4, ea, pc[2], temp8a);
3460  eaclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3461  eac);
3462 
3463  temp48alen = fast_expansion_sum_zeroelim(cdelen, cde, bcelen, bce, temp48a);
3464  temp48blen = fast_expansion_sum_zeroelim(deblen, deb, bcdlen, bcd, temp48b);
3465  for (i = 0; i < temp48blen; i++) {
3466  temp48b[i] = -temp48b[i];
3467  }
3468  bcdelen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
3469  temp48blen, temp48b, bcde);
3470  xlen = scale_expansion_zeroelim(bcdelen, bcde, pa[0], temp192);
3471  xlen = scale_expansion_zeroelim(xlen, temp192, pa[0], det384x);
3472  ylen = scale_expansion_zeroelim(bcdelen, bcde, pa[1], temp192);
3473  ylen = scale_expansion_zeroelim(ylen, temp192, pa[1], det384y);
3474  zlen = scale_expansion_zeroelim(bcdelen, bcde, pa[2], temp192);
3475  zlen = scale_expansion_zeroelim(zlen, temp192, pa[2], det384z);
3476  xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
3477  alen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, adet);
3478 
3479  temp48alen = fast_expansion_sum_zeroelim(dealen, dea, cdalen, cda, temp48a);
3480  temp48blen = fast_expansion_sum_zeroelim(eaclen, eac, cdelen, cde, temp48b);
3481  for (i = 0; i < temp48blen; i++) {
3482  temp48b[i] = -temp48b[i];
3483  }
3484  cdealen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
3485  temp48blen, temp48b, cdea);
3486  xlen = scale_expansion_zeroelim(cdealen, cdea, pb[0], temp192);
3487  xlen = scale_expansion_zeroelim(xlen, temp192, pb[0], det384x);
3488  ylen = scale_expansion_zeroelim(cdealen, cdea, pb[1], temp192);
3489  ylen = scale_expansion_zeroelim(ylen, temp192, pb[1], det384y);
3490  zlen = scale_expansion_zeroelim(cdealen, cdea, pb[2], temp192);
3491  zlen = scale_expansion_zeroelim(zlen, temp192, pb[2], det384z);
3492  xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
3493  blen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, bdet);
3494 
3495  temp48alen = fast_expansion_sum_zeroelim(eablen, eab, deblen, deb, temp48a);
3496  temp48blen = fast_expansion_sum_zeroelim(abdlen, abd, dealen, dea, temp48b);
3497  for (i = 0; i < temp48blen; i++) {
3498  temp48b[i] = -temp48b[i];
3499  }
3500  deablen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
3501  temp48blen, temp48b, deab);
3502  xlen = scale_expansion_zeroelim(deablen, deab, pc[0], temp192);
3503  xlen = scale_expansion_zeroelim(xlen, temp192, pc[0], det384x);
3504  ylen = scale_expansion_zeroelim(deablen, deab, pc[1], temp192);
3505  ylen = scale_expansion_zeroelim(ylen, temp192, pc[1], det384y);
3506  zlen = scale_expansion_zeroelim(deablen, deab, pc[2], temp192);
3507  zlen = scale_expansion_zeroelim(zlen, temp192, pc[2], det384z);
3508  xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
3509  clen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, cdet);
3510 
3511  temp48alen = fast_expansion_sum_zeroelim(abclen, abc, eaclen, eac, temp48a);
3512  temp48blen = fast_expansion_sum_zeroelim(bcelen, bce, eablen, eab, temp48b);
3513  for (i = 0; i < temp48blen; i++) {
3514  temp48b[i] = -temp48b[i];
3515  }
3516  eabclen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
3517  temp48blen, temp48b, eabc);
3518  xlen = scale_expansion_zeroelim(eabclen, eabc, pd[0], temp192);
3519  xlen = scale_expansion_zeroelim(xlen, temp192, pd[0], det384x);
3520  ylen = scale_expansion_zeroelim(eabclen, eabc, pd[1], temp192);
3521  ylen = scale_expansion_zeroelim(ylen, temp192, pd[1], det384y);
3522  zlen = scale_expansion_zeroelim(eabclen, eabc, pd[2], temp192);
3523  zlen = scale_expansion_zeroelim(zlen, temp192, pd[2], det384z);
3524  xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
3525  dlen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, ddet);
3526 
3527  temp48alen = fast_expansion_sum_zeroelim(bcdlen, bcd, abdlen, abd, temp48a);
3528  temp48blen = fast_expansion_sum_zeroelim(cdalen, cda, abclen, abc, temp48b);
3529  for (i = 0; i < temp48blen; i++) {
3530  temp48b[i] = -temp48b[i];
3531  }
3532  abcdlen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
3533  temp48blen, temp48b, abcd);
3534  xlen = scale_expansion_zeroelim(abcdlen, abcd, pe[0], temp192);
3535  xlen = scale_expansion_zeroelim(xlen, temp192, pe[0], det384x);
3536  ylen = scale_expansion_zeroelim(abcdlen, abcd, pe[1], temp192);
3537  ylen = scale_expansion_zeroelim(ylen, temp192, pe[1], det384y);
3538  zlen = scale_expansion_zeroelim(abcdlen, abcd, pe[2], temp192);
3539  zlen = scale_expansion_zeroelim(zlen, temp192, pe[2], det384z);
3540  xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
3541  elen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, edet);
3542 
3543  ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
3544  cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
3545  cdelen = fast_expansion_sum_zeroelim(cdlen, cddet, elen, edet, cdedet);
3546  deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdelen, cdedet, deter);
3547 
3548  return deter[deterlen - 1];
3549 }
3550 
3551 REAL insphereslow(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
3552 {
3553  INEXACT REAL aex, bex, cex, dex, aey, bey, cey, dey, aez, bez, cez, dez;
3554  REAL aextail, bextail, cextail, dextail;
3555  REAL aeytail, beytail, ceytail, deytail;
3556  REAL aeztail, beztail, ceztail, deztail;
3557  REAL negate, negatetail;
3558  INEXACT REAL axby7, bxcy7, cxdy7, dxay7, axcy7, bxdy7;
3559  INEXACT REAL bxay7, cxby7, dxcy7, axdy7, cxay7, dxby7;
3560  REAL axby[8], bxcy[8], cxdy[8], dxay[8], axcy[8], bxdy[8];
3561  REAL bxay[8], cxby[8], dxcy[8], axdy[8], cxay[8], dxby[8];
3562  REAL ab[16], bc[16], cd[16], da[16], ac[16], bd[16];
3563  int ablen, bclen, cdlen, dalen, aclen, bdlen;
3564  REAL temp32a[32], temp32b[32], temp64a[64], temp64b[64], temp64c[64];
3565  int temp32alen, temp32blen, temp64alen, temp64blen, temp64clen;
3566  REAL temp128[128], temp192[192];
3567  int temp128len, temp192len;
3568  REAL detx[384], detxx[768], detxt[384], detxxt[768], detxtxt[768];
3569  int xlen, xxlen, xtlen, xxtlen, xtxtlen;
3570  REAL x1[1536], x2[2304];
3571  int x1len, x2len;
3572  REAL dety[384], detyy[768], detyt[384], detyyt[768], detytyt[768];
3573  int ylen, yylen, ytlen, yytlen, ytytlen;
3574  REAL y1[1536], y2[2304];
3575  int y1len, y2len;
3576  REAL detz[384], detzz[768], detzt[384], detzzt[768], detztzt[768];
3577  int zlen, zzlen, ztlen, zztlen, ztztlen;
3578  REAL z1[1536], z2[2304];
3579  int z1len, z2len;
3580  REAL detxy[4608];
3581  int xylen;
3582  REAL adet[6912], bdet[6912], cdet[6912], ddet[6912];
3583  int alen, blen, clen, dlen;
3584  REAL abdet[13824], cddet[13824], deter[27648];
3585  int deterlen;
3586  int i;
3587 
3588  INEXACT REAL bvirt;
3589  REAL avirt, bround, around;
3590  INEXACT REAL c;
3591  INEXACT REAL abig;
3592  REAL a0hi, a0lo, a1hi, a1lo, bhi, blo;
3593  REAL err1, err2, err3;
3594  INEXACT REAL _i, _j, _k, _l, _m, _n;
3595  REAL _0, _1, _2;
3596 
3597  Two_Diff(pa[0], pe[0], aex, aextail);
3598  Two_Diff(pa[1], pe[1], aey, aeytail);
3599  Two_Diff(pa[2], pe[2], aez, aeztail);
3600  Two_Diff(pb[0], pe[0], bex, bextail);
3601  Two_Diff(pb[1], pe[1], bey, beytail);
3602  Two_Diff(pb[2], pe[2], bez, beztail);
3603  Two_Diff(pc[0], pe[0], cex, cextail);
3604  Two_Diff(pc[1], pe[1], cey, ceytail);
3605  Two_Diff(pc[2], pe[2], cez, ceztail);
3606  Two_Diff(pd[0], pe[0], dex, dextail);
3607  Two_Diff(pd[1], pe[1], dey, deytail);
3608  Two_Diff(pd[2], pe[2], dez, deztail);
3609 
3610  Two_Two_Product(aex, aextail, bey, beytail,
3611  axby7, axby[6], axby[5], axby[4],
3612  axby[3], axby[2], axby[1], axby[0]);
3613  axby[7] = axby7;
3614  negate = -aey;
3615  negatetail = -aeytail;
3616  Two_Two_Product(bex, bextail, negate, negatetail,
3617  bxay7, bxay[6], bxay[5], bxay[4],
3618  bxay[3], bxay[2], bxay[1], bxay[0]);
3619  bxay[7] = bxay7;
3620  ablen = fast_expansion_sum_zeroelim(8, axby, 8, bxay, ab);
3621  Two_Two_Product(bex, bextail, cey, ceytail,
3622  bxcy7, bxcy[6], bxcy[5], bxcy[4],
3623  bxcy[3], bxcy[2], bxcy[1], bxcy[0]);
3624  bxcy[7] = bxcy7;
3625  negate = -bey;
3626  negatetail = -beytail;
3627  Two_Two_Product(cex, cextail, negate, negatetail,
3628  cxby7, cxby[6], cxby[5], cxby[4],
3629  cxby[3], cxby[2], cxby[1], cxby[0]);
3630  cxby[7] = cxby7;
3631  bclen = fast_expansion_sum_zeroelim(8, bxcy, 8, cxby, bc);
3632  Two_Two_Product(cex, cextail, dey, deytail,
3633  cxdy7, cxdy[6], cxdy[5], cxdy[4],
3634  cxdy[3], cxdy[2], cxdy[1], cxdy[0]);
3635  cxdy[7] = cxdy7;
3636  negate = -cey;
3637  negatetail = -ceytail;
3638  Two_Two_Product(dex, dextail, negate, negatetail,
3639  dxcy7, dxcy[6], dxcy[5], dxcy[4],
3640  dxcy[3], dxcy[2], dxcy[1], dxcy[0]);
3641  dxcy[7] = dxcy7;
3642  cdlen = fast_expansion_sum_zeroelim(8, cxdy, 8, dxcy, cd);
3643  Two_Two_Product(dex, dextail, aey, aeytail,
3644  dxay7, dxay[6], dxay[5], dxay[4],
3645  dxay[3], dxay[2], dxay[1], dxay[0]);
3646  dxay[7] = dxay7;
3647  negate = -dey;
3648  negatetail = -deytail;
3649  Two_Two_Product(aex, aextail, negate, negatetail,
3650  axdy7, axdy[6], axdy[5], axdy[4],
3651  axdy[3], axdy[2], axdy[1], axdy[0]);
3652  axdy[7] = axdy7;
3653  dalen = fast_expansion_sum_zeroelim(8, dxay, 8, axdy, da);
3654  Two_Two_Product(aex, aextail, cey, ceytail,
3655  axcy7, axcy[6], axcy[5], axcy[4],
3656  axcy[3], axcy[2], axcy[1], axcy[0]);
3657  axcy[7] = axcy7;
3658  negate = -aey;
3659  negatetail = -aeytail;
3660  Two_Two_Product(cex, cextail, negate, negatetail,
3661  cxay7, cxay[6], cxay[5], cxay[4],
3662  cxay[3], cxay[2], cxay[1], cxay[0]);
3663  cxay[7] = cxay7;
3664  aclen = fast_expansion_sum_zeroelim(8, axcy, 8, cxay, ac);
3665  Two_Two_Product(bex, bextail, dey, deytail,
3666  bxdy7, bxdy[6], bxdy[5], bxdy[4],
3667  bxdy[3], bxdy[2], bxdy[1], bxdy[0]);
3668  bxdy[7] = bxdy7;
3669  negate = -bey;
3670  negatetail = -beytail;
3671  Two_Two_Product(dex, dextail, negate, negatetail,
3672  dxby7, dxby[6], dxby[5], dxby[4],
3673  dxby[3], dxby[2], dxby[1], dxby[0]);
3674  dxby[7] = dxby7;
3675  bdlen = fast_expansion_sum_zeroelim(8, bxdy, 8, dxby, bd);
3676 
3677  temp32alen = scale_expansion_zeroelim(cdlen, cd, -bez, temp32a);
3678  temp32blen = scale_expansion_zeroelim(cdlen, cd, -beztail, temp32b);
3679  temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3680  temp32blen, temp32b, temp64a);
3681  temp32alen = scale_expansion_zeroelim(bdlen, bd, cez, temp32a);
3682  temp32blen = scale_expansion_zeroelim(bdlen, bd, ceztail, temp32b);
3683  temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3684  temp32blen, temp32b, temp64b);
3685  temp32alen = scale_expansion_zeroelim(bclen, bc, -dez, temp32a);
3686  temp32blen = scale_expansion_zeroelim(bclen, bc, -deztail, temp32b);
3687  temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3688  temp32blen, temp32b, temp64c);
3689  temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a,
3690  temp64blen, temp64b, temp128);
3691  temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c,
3692  temp128len, temp128, temp192);
3693  xlen = scale_expansion_zeroelim(temp192len, temp192, aex, detx);
3694  xxlen = scale_expansion_zeroelim(xlen, detx, aex, detxx);
3695  xtlen = scale_expansion_zeroelim(temp192len, temp192, aextail, detxt);
3696  xxtlen = scale_expansion_zeroelim(xtlen, detxt, aex, detxxt);
3697  for (i = 0; i < xxtlen; i++) {
3698  detxxt[i] *= 2.0;
3699  }
3700  xtxtlen = scale_expansion_zeroelim(xtlen, detxt, aextail, detxtxt);
3701  x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
3702  x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
3703  ylen = scale_expansion_zeroelim(temp192len, temp192, aey, dety);
3704  yylen = scale_expansion_zeroelim(ylen, dety, aey, detyy);
3705  ytlen = scale_expansion_zeroelim(temp192len, temp192, aeytail, detyt);
3706  yytlen = scale_expansion_zeroelim(ytlen, detyt, aey, detyyt);
3707  for (i = 0; i < yytlen; i++) {
3708  detyyt[i] *= 2.0;
3709  }
3710  ytytlen = scale_expansion_zeroelim(ytlen, detyt, aeytail, detytyt);
3711  y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
3712  y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
3713  zlen = scale_expansion_zeroelim(temp192len, temp192, aez, detz);
3714  zzlen = scale_expansion_zeroelim(zlen, detz, aez, detzz);
3715  ztlen = scale_expansion_zeroelim(temp192len, temp192, aeztail, detzt);
3716  zztlen = scale_expansion_zeroelim(ztlen, detzt, aez, detzzt);
3717  for (i = 0; i < zztlen; i++) {
3718  detzzt[i] *= 2.0;
3719  }
3720  ztztlen = scale_expansion_zeroelim(ztlen, detzt, aeztail, detztzt);
3721  z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1);
3722  z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2);
3723  xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy);
3724  alen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, adet);
3725 
3726  temp32alen = scale_expansion_zeroelim(dalen, da, cez, temp32a);
3727  temp32blen = scale_expansion_zeroelim(dalen, da, ceztail, temp32b);
3728  temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3729  temp32blen, temp32b, temp64a);
3730  temp32alen = scale_expansion_zeroelim(aclen, ac, dez, temp32a);
3731  temp32blen = scale_expansion_zeroelim(aclen, ac, deztail, temp32b);
3732  temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3733  temp32blen, temp32b, temp64b);
3734  temp32alen = scale_expansion_zeroelim(cdlen, cd, aez, temp32a);
3735  temp32blen = scale_expansion_zeroelim(cdlen, cd, aeztail, temp32b);
3736  temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3737  temp32blen, temp32b, temp64c);
3738  temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a,
3739  temp64blen, temp64b, temp128);
3740  temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c,
3741  temp128len, temp128, temp192);
3742  xlen = scale_expansion_zeroelim(temp192len, temp192, bex, detx);
3743  xxlen = scale_expansion_zeroelim(xlen, detx, bex, detxx);
3744  xtlen = scale_expansion_zeroelim(temp192len, temp192, bextail, detxt);
3745  xxtlen = scale_expansion_zeroelim(xtlen, detxt, bex, detxxt);
3746  for (i = 0; i < xxtlen; i++) {
3747  detxxt[i] *= 2.0;
3748  }
3749  xtxtlen = scale_expansion_zeroelim(xtlen, detxt, bextail, detxtxt);
3750  x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
3751  x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
3752  ylen = scale_expansion_zeroelim(temp192len, temp192, bey, dety);
3753  yylen = scale_expansion_zeroelim(ylen, dety, bey, detyy);
3754  ytlen = scale_expansion_zeroelim(temp192len, temp192, beytail, detyt);
3755  yytlen = scale_expansion_zeroelim(ytlen, detyt, bey, detyyt);
3756  for (i = 0; i < yytlen; i++) {
3757  detyyt[i] *= 2.0;
3758  }
3759  ytytlen = scale_expansion_zeroelim(ytlen, detyt, beytail, detytyt);
3760  y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
3761  y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
3762  zlen = scale_expansion_zeroelim(temp192len, temp192, bez, detz);
3763  zzlen = scale_expansion_zeroelim(zlen, detz, bez, detzz);
3764  ztlen = scale_expansion_zeroelim(temp192len, temp192, beztail, detzt);
3765  zztlen = scale_expansion_zeroelim(ztlen, detzt, bez, detzzt);
3766  for (i = 0; i < zztlen; i++) {
3767  detzzt[i] *= 2.0;
3768  }
3769  ztztlen = scale_expansion_zeroelim(ztlen, detzt, beztail, detztzt);
3770  z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1);
3771  z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2);
3772  xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy);
3773  blen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, bdet);
3774 
3775  temp32alen = scale_expansion_zeroelim(ablen, ab, -dez, temp32a);
3776  temp32blen = scale_expansion_zeroelim(ablen, ab, -deztail, temp32b);
3777  temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3778  temp32blen, temp32b, temp64a);
3779  temp32alen = scale_expansion_zeroelim(bdlen, bd, -aez, temp32a);
3780  temp32blen = scale_expansion_zeroelim(bdlen, bd, -aeztail, temp32b);
3781  temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3782  temp32blen, temp32b, temp64b);
3783  temp32alen = scale_expansion_zeroelim(dalen, da, -bez, temp32a);
3784  temp32blen = scale_expansion_zeroelim(dalen, da, -beztail, temp32b);
3785  temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3786  temp32blen, temp32b, temp64c);
3787  temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a,
3788  temp64blen, temp64b, temp128);
3789  temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c,
3790  temp128len, temp128, temp192);
3791  xlen = scale_expansion_zeroelim(temp192len, temp192, cex, detx);
3792  xxlen = scale_expansion_zeroelim(xlen, detx, cex, detxx);
3793  xtlen = scale_expansion_zeroelim(temp192len, temp192, cextail, detxt);
3794  xxtlen = scale_expansion_zeroelim(xtlen, detxt, cex, detxxt);
3795  for (i = 0; i < xxtlen; i++) {
3796  detxxt[i] *= 2.0;
3797  }
3798  xtxtlen = scale_expansion_zeroelim(xtlen, detxt, cextail, detxtxt);
3799  x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
3800  x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
3801  ylen = scale_expansion_zeroelim(temp192len, temp192, cey, dety);
3802  yylen = scale_expansion_zeroelim(ylen, dety, cey, detyy);
3803  ytlen = scale_expansion_zeroelim(temp192len, temp192, ceytail, detyt);
3804  yytlen = scale_expansion_zeroelim(ytlen, detyt, cey, detyyt);
3805  for (i = 0; i < yytlen; i++) {
3806  detyyt[i] *= 2.0;
3807  }
3808  ytytlen = scale_expansion_zeroelim(ytlen, detyt, ceytail, detytyt);
3809  y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
3810  y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
3811  zlen = scale_expansion_zeroelim(temp192len, temp192, cez, detz);
3812  zzlen = scale_expansion_zeroelim(zlen, detz, cez, detzz);
3813  ztlen = scale_expansion_zeroelim(temp192len, temp192, ceztail, detzt);
3814  zztlen = scale_expansion_zeroelim(ztlen, detzt, cez, detzzt);
3815  for (i = 0; i < zztlen; i++) {
3816  detzzt[i] *= 2.0;
3817  }
3818  ztztlen = scale_expansion_zeroelim(ztlen, detzt, ceztail, detztzt);
3819  z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1);
3820  z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2);
3821  xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy);
3822  clen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, cdet);
3823 
3824  temp32alen = scale_expansion_zeroelim(bclen, bc, aez, temp32a);
3825  temp32blen = scale_expansion_zeroelim(bclen, bc, aeztail, temp32b);
3826  temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3827  temp32blen, temp32b, temp64a);
3828  temp32alen = scale_expansion_zeroelim(aclen, ac, -bez, temp32a);
3829  temp32blen = scale_expansion_zeroelim(aclen, ac, -beztail, temp32b);
3830  temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3831  temp32blen, temp32b, temp64b);
3832  temp32alen = scale_expansion_zeroelim(ablen, ab, cez, temp32a);
3833  temp32blen = scale_expansion_zeroelim(ablen, ab, ceztail, temp32b);
3834  temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3835  temp32blen, temp32b, temp64c);
3836  temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a,
3837  temp64blen, temp64b, temp128);
3838  temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c,
3839  temp128len, temp128, temp192);
3840  xlen = scale_expansion_zeroelim(temp192len, temp192, dex, detx);
3841  xxlen = scale_expansion_zeroelim(xlen, detx, dex, detxx);
3842  xtlen = scale_expansion_zeroelim(temp192len, temp192, dextail, detxt);
3843  xxtlen = scale_expansion_zeroelim(xtlen, detxt, dex, detxxt);
3844  for (i = 0; i < xxtlen; i++) {
3845  detxxt[i] *= 2.0;
3846  }
3847  xtxtlen = scale_expansion_zeroelim(xtlen, detxt, dextail, detxtxt);
3848  x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
3849  x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
3850  ylen = scale_expansion_zeroelim(temp192len, temp192, dey, dety);
3851  yylen = scale_expansion_zeroelim(ylen, dety, dey, detyy);
3852  ytlen = scale_expansion_zeroelim(temp192len, temp192, deytail, detyt);
3853  yytlen = scale_expansion_zeroelim(ytlen, detyt, dey, detyyt);
3854  for (i = 0; i < yytlen; i++) {
3855  detyyt[i] *= 2.0;
3856  }
3857  ytytlen = scale_expansion_zeroelim(ytlen, detyt, deytail, detytyt);
3858  y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
3859  y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
3860  zlen = scale_expansion_zeroelim(temp192len, temp192, dez, detz);
3861  zzlen = scale_expansion_zeroelim(zlen, detz, dez, detzz);
3862  ztlen = scale_expansion_zeroelim(temp192len, temp192, deztail, detzt);
3863  zztlen = scale_expansion_zeroelim(ztlen, detzt, dez, detzzt);
3864  for (i = 0; i < zztlen; i++) {
3865  detzzt[i] *= 2.0;
3866  }
3867  ztztlen = scale_expansion_zeroelim(ztlen, detzt, deztail, detztzt);
3868  z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1);
3869  z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2);
3870  xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy);
3871  dlen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, ddet);
3872 
3873  ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
3874  cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
3875  deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, deter);
3876 
3877  return deter[deterlen - 1];
3878 }
3879 
3880 REAL insphereadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe,
3881  REAL permanent)
3882 {
3883  INEXACT REAL aex, bex, cex, dex, aey, bey, cey, dey, aez, bez, cez, dez;
3884  REAL det, errbound;
3885 
3886  INEXACT REAL aexbey1, bexaey1, bexcey1, cexbey1;
3887  INEXACT REAL cexdey1, dexcey1, dexaey1, aexdey1;
3888  INEXACT REAL aexcey1, cexaey1, bexdey1, dexbey1;
3889  REAL aexbey0, bexaey0, bexcey0, cexbey0;
3890  REAL cexdey0, dexcey0, dexaey0, aexdey0;
3891  REAL aexcey0, cexaey0, bexdey0, dexbey0;
3892  REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4];
3893  INEXACT REAL ab3, bc3, cd3, da3, ac3, bd3;
3894  REAL abeps, bceps, cdeps, daeps, aceps, bdeps;
3895  REAL temp8a[8], temp8b[8], temp8c[8], temp16[16], temp24[24], temp48[48];
3896  int temp8alen, temp8blen, temp8clen, temp16len, temp24len, temp48len;
3897  REAL xdet[96], ydet[96], zdet[96], xydet[192];
3898  int xlen, ylen, zlen, xylen;
3899  REAL adet[288], bdet[288], cdet[288], ddet[288];
3900  int alen, blen, clen, dlen;
3901  REAL abdet[576], cddet[576];
3902  int ablen, cdlen;
3903  REAL fin1[1152];
3904  int finlength;
3905 
3906  REAL aextail, bextail, cextail, dextail;
3907  REAL aeytail, beytail, ceytail, deytail;
3908  REAL aeztail, beztail, ceztail, deztail;
3909 
3910  INEXACT REAL bvirt;
3911  REAL avirt, bround, around;
3912  INEXACT REAL c;
3913  INEXACT REAL abig;
3914  REAL ahi, alo, bhi, blo;
3915  REAL err1, err2, err3;
3916  INEXACT REAL _i, _j;
3917  REAL _0;
3918 
3919  aex = (REAL) (pa[0] - pe[0]);
3920  bex = (REAL) (pb[0] - pe[0]);
3921  cex = (REAL) (pc[0] - pe[0]);
3922  dex = (REAL) (pd[0] - pe[0]);
3923  aey = (REAL) (pa[1] - pe[1]);
3924  bey = (REAL) (pb[1] - pe[1]);
3925  cey = (REAL) (pc[1] - pe[1]);
3926  dey = (REAL) (pd[1] - pe[1]);
3927  aez = (REAL) (pa[2] - pe[2]);
3928  bez = (REAL) (pb[2] - pe[2]);
3929  cez = (REAL) (pc[2] - pe[2]);
3930  dez = (REAL) (pd[2] - pe[2]);
3931 
3932  Two_Product(aex, bey, aexbey1, aexbey0);
3933  Two_Product(bex, aey, bexaey1, bexaey0);
3934  Two_Two_Diff(aexbey1, aexbey0, bexaey1, bexaey0, ab3, ab[2], ab[1], ab[0]);
3935  ab[3] = ab3;
3936 
3937  Two_Product(bex, cey, bexcey1, bexcey0);
3938  Two_Product(cex, bey, cexbey1, cexbey0);
3939  Two_Two_Diff(bexcey1, bexcey0, cexbey1, cexbey0, bc3, bc[2], bc[1], bc[0]);
3940  bc[3] = bc3;
3941 
3942  Two_Product(cex, dey, cexdey1, cexdey0);
3943  Two_Product(dex, cey, dexcey1, dexcey0);
3944  Two_Two_Diff(cexdey1, cexdey0, dexcey1, dexcey0, cd3, cd[2], cd[1], cd[0]);
3945  cd[3] = cd3;
3946 
3947  Two_Product(dex, aey, dexaey1, dexaey0);
3948  Two_Product(aex, dey, aexdey1, aexdey0);
3949  Two_Two_Diff(dexaey1, dexaey0, aexdey1, aexdey0, da3, da[2], da[1], da[0]);
3950  da[3] = da3;
3951 
3952  Two_Product(aex, cey, aexcey1, aexcey0);
3953  Two_Product(cex, aey, cexaey1, cexaey0);
3954  Two_Two_Diff(aexcey1, aexcey0, cexaey1, cexaey0, ac3, ac[2], ac[1], ac[0]);
3955  ac[3] = ac3;
3956 
3957  Two_Product(bex, dey, bexdey1, bexdey0);
3958  Two_Product(dex, bey, dexbey1, dexbey0);
3959  Two_Two_Diff(bexdey1, bexdey0, dexbey1, dexbey0, bd3, bd[2], bd[1], bd[0]);
3960  bd[3] = bd3;
3961 
3962  temp8alen = scale_expansion_zeroelim(4, cd, bez, temp8a);
3963  temp8blen = scale_expansion_zeroelim(4, bd, -cez, temp8b);
3964  temp8clen = scale_expansion_zeroelim(4, bc, dez, temp8c);
3965  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
3966  temp8blen, temp8b, temp16);
3967  temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
3968  temp16len, temp16, temp24);
3969  temp48len = scale_expansion_zeroelim(temp24len, temp24, aex, temp48);
3970  xlen = scale_expansion_zeroelim(temp48len, temp48, -aex, xdet);
3971  temp48len = scale_expansion_zeroelim(temp24len, temp24, aey, temp48);
3972  ylen = scale_expansion_zeroelim(temp48len, temp48, -aey, ydet);
3973  temp48len = scale_expansion_zeroelim(temp24len, temp24, aez, temp48);
3974  zlen = scale_expansion_zeroelim(temp48len, temp48, -aez, zdet);
3975  xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
3976  alen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, adet);
3977 
3978  temp8alen = scale_expansion_zeroelim(4, da, cez, temp8a);
3979  temp8blen = scale_expansion_zeroelim(4, ac, dez, temp8b);
3980  temp8clen = scale_expansion_zeroelim(4, cd, aez, temp8c);
3981  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
3982  temp8blen, temp8b, temp16);
3983  temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
3984  temp16len, temp16, temp24);
3985  temp48len = scale_expansion_zeroelim(temp24len, temp24, bex, temp48);
3986  xlen = scale_expansion_zeroelim(temp48len, temp48, bex, xdet);
3987  temp48len = scale_expansion_zeroelim(temp24len, temp24, bey, temp48);
3988  ylen = scale_expansion_zeroelim(temp48len, temp48, bey, ydet);
3989  temp48len = scale_expansion_zeroelim(temp24len, temp24, bez, temp48);
3990  zlen = scale_expansion_zeroelim(temp48len, temp48, bez, zdet);
3991  xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
3992  blen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, bdet);
3993 
3994  temp8alen = scale_expansion_zeroelim(4, ab, dez, temp8a);
3995  temp8blen = scale_expansion_zeroelim(4, bd, aez, temp8b);
3996  temp8clen = scale_expansion_zeroelim(4, da, bez, temp8c);
3997  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
3998  temp8blen, temp8b, temp16);
3999  temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
4000  temp16len, temp16, temp24);
4001  temp48len = scale_expansion_zeroelim(temp24len, temp24, cex, temp48);
4002  xlen = scale_expansion_zeroelim(temp48len, temp48, -cex, xdet);
4003  temp48len = scale_expansion_zeroelim(temp24len, temp24, cey, temp48);
4004  ylen = scale_expansion_zeroelim(temp48len, temp48, -cey, ydet);
4005  temp48len = scale_expansion_zeroelim(temp24len, temp24, cez, temp48);
4006  zlen = scale_expansion_zeroelim(temp48len, temp48, -cez, zdet);
4007  xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
4008  clen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, cdet);
4009 
4010  temp8alen = scale_expansion_zeroelim(4, bc, aez, temp8a);
4011  temp8blen = scale_expansion_zeroelim(4, ac, -bez, temp8b);
4012  temp8clen = scale_expansion_zeroelim(4, ab, cez, temp8c);
4013  temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
4014  temp8blen, temp8b, temp16);
4015  temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
4016  temp16len, temp16, temp24);
4017  temp48len = scale_expansion_zeroelim(temp24len, temp24, dex, temp48);
4018  xlen = scale_expansion_zeroelim(temp48len, temp48, dex, xdet);
4019  temp48len = scale_expansion_zeroelim(temp24len, temp24, dey, temp48);
4020  ylen = scale_expansion_zeroelim(temp48len, temp48, dey, ydet);
4021  temp48len = scale_expansion_zeroelim(temp24len, temp24, dez, temp48);
4022  zlen = scale_expansion_zeroelim(temp48len, temp48, dez, zdet);
4023  xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
4024  dlen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, ddet);
4025 
4026  ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
4027  cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
4028  finlength = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, fin1);
4029 
4030  det = estimate(finlength, fin1);
4031  errbound = isperrboundB * permanent;
4032  if ((det >= errbound) || (-det >= errbound)) {
4033  return det;
4034  }
4035 
4036  Two_Diff_Tail(pa[0], pe[0], aex, aextail);
4037  Two_Diff_Tail(pa[1], pe[1], aey, aeytail);
4038  Two_Diff_Tail(pa[2], pe[2], aez, aeztail);
4039  Two_Diff_Tail(pb[0], pe[0], bex, bextail);
4040  Two_Diff_Tail(pb[1], pe[1], bey, beytail);
4041  Two_Diff_Tail(pb[2], pe[2], bez, beztail);
4042  Two_Diff_Tail(pc[0], pe[0], cex, cextail);
4043  Two_Diff_Tail(pc[1], pe[1], cey, ceytail);
4044  Two_Diff_Tail(pc[2], pe[2], cez, ceztail);
4045  Two_Diff_Tail(pd[0], pe[0], dex, dextail);
4046  Two_Diff_Tail(pd[1], pe[1], dey, deytail);
4047  Two_Diff_Tail(pd[2], pe[2], dez, deztail);
4048  if ((aextail == 0.0) && (aeytail == 0.0) && (aeztail == 0.0)
4049  && (bextail == 0.0) && (beytail == 0.0) && (beztail == 0.0)
4050  && (cextail == 0.0) && (ceytail == 0.0) && (ceztail == 0.0)
4051  && (dextail == 0.0) && (deytail == 0.0) && (deztail == 0.0)) {
4052  return det;
4053  }
4054 
4055  errbound = isperrboundC * permanent + resulterrbound * Absolute(det);
4056  abeps = (aex * beytail + bey * aextail)
4057  - (aey * bextail + bex * aeytail);
4058  bceps = (bex * ceytail + cey * bextail)
4059  - (bey * cextail + cex * beytail);
4060  cdeps = (cex * deytail + dey * cextail)
4061  - (cey * dextail + dex * ceytail);
4062  daeps = (dex * aeytail + aey * dextail)
4063  - (dey * aextail + aex * deytail);
4064  aceps = (aex * ceytail + cey * aextail)
4065  - (aey * cextail + cex * aeytail);
4066  bdeps = (bex * deytail + dey * bextail)
4067  - (bey * dextail + dex * beytail);
4068  det += (((bex * bex + bey * bey + bez * bez)
4069  * ((cez * daeps + dez * aceps + aez * cdeps)
4070  + (ceztail * da3 + deztail * ac3 + aeztail * cd3))
4071  + (dex * dex + dey * dey + dez * dez)
4072  * ((aez * bceps - bez * aceps + cez * abeps)
4073  + (aeztail * bc3 - beztail * ac3 + ceztail * ab3)))
4074  - ((aex * aex + aey * aey + aez * aez)
4075  * ((bez * cdeps - cez * bdeps + dez * bceps)
4076  + (beztail * cd3 - ceztail * bd3 + deztail * bc3))
4077  + (cex * cex + cey * cey + cez * cez)
4078  * ((dez * abeps + aez * bdeps + bez * daeps)
4079  + (deztail * ab3 + aeztail * bd3 + beztail * da3))))
4080  + 2.0 * (((bex * bextail + bey * beytail + bez * beztail)
4081  * (cez * da3 + dez * ac3 + aez * cd3)
4082  + (dex * dextail + dey * deytail + dez * deztail)
4083  * (aez * bc3 - bez * ac3 + cez * ab3))
4084  - ((aex * aextail + aey * aeytail + aez * aeztail)
4085  * (bez * cd3 - cez * bd3 + dez * bc3)
4086  + (cex * cextail + cey * ceytail + cez * ceztail)
4087  * (dez * ab3 + aez * bd3 + bez * da3)));
4088  if ((det >= errbound) || (-det >= errbound)) {
4089  return det;
4090  }
4091 
4092  return insphereexact(pa, pb, pc, pd, pe);
4093 }
4094 
4095 REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
4096 {
4097  REAL aex, bex, cex, dex;
4098  REAL aey, bey, cey, dey;
4099  REAL aez, bez, cez, dez;
4100  REAL aexbey, bexaey, bexcey, cexbey, cexdey, dexcey, dexaey, aexdey;
4101  REAL aexcey, cexaey, bexdey, dexbey;
4102  REAL alift, blift, clift, dlift;
4103  REAL ab, bc, cd, da, ac, bd;
4104  REAL abc, bcd, cda, dab;
4105  REAL aezplus, bezplus, cezplus, dezplus;
4106  REAL aexbeyplus, bexaeyplus, bexceyplus, cexbeyplus;
4107  REAL cexdeyplus, dexceyplus, dexaeyplus, aexdeyplus;
4108  REAL aexceyplus, cexaeyplus, bexdeyplus, dexbeyplus;
4109  REAL det;
4110  REAL permanent, errbound;
4111 
4112  aex = pa[0] - pe[0];
4113  bex = pb[0] - pe[0];
4114  cex = pc[0] - pe[0];
4115  dex = pd[0] - pe[0];
4116  aey = pa[1] - pe[1];
4117  bey = pb[1] - pe[1];
4118  cey = pc[1] - pe[1];
4119  dey = pd[1] - pe[1];
4120  aez = pa[2] - pe[2];
4121  bez = pb[2] - pe[2];
4122  cez = pc[2] - pe[2];
4123  dez = pd[2] - pe[2];
4124 
4125  aexbey = aex * bey;
4126  bexaey = bex * aey;
4127  ab = aexbey - bexaey;
4128  bexcey = bex * cey;
4129  cexbey = cex * bey;
4130  bc = bexcey - cexbey;
4131  cexdey = cex * dey;
4132  dexcey = dex * cey;
4133  cd = cexdey - dexcey;
4134  dexaey = dex * aey;
4135  aexdey = aex * dey;
4136  da = dexaey - aexdey;
4137 
4138  aexcey = aex * cey;
4139  cexaey = cex * aey;
4140  ac = aexcey - cexaey;
4141  bexdey = bex * dey;
4142  dexbey = dex * bey;
4143  bd = bexdey - dexbey;
4144 
4145  abc = aez * bc - bez * ac + cez * ab;
4146  bcd = bez * cd - cez * bd + dez * bc;
4147  cda = cez * da + dez * ac + aez * cd;
4148  dab = dez * ab + aez * bd + bez * da;
4149 
4150  alift = aex * aex + aey * aey + aez * aez;
4151  blift = bex * bex + bey * bey + bez * bez;
4152  clift = cex * cex + cey * cey + cez * cez;
4153  dlift = dex * dex + dey * dey + dez * dez;
4154 
4155  det = (dlift * abc - clift * dab) + (blift * cda - alift * bcd);
4156 
4157  aezplus = Absolute(aez);
4158  bezplus = Absolute(bez);
4159  cezplus = Absolute(cez);
4160  dezplus = Absolute(dez);
4161  aexbeyplus = Absolute(aexbey);
4162  bexaeyplus = Absolute(bexaey);
4163  bexceyplus = Absolute(bexcey);
4164  cexbeyplus = Absolute(cexbey);
4165  cexdeyplus = Absolute(cexdey);
4166  dexceyplus = Absolute(dexcey);
4167  dexaeyplus = Absolute(dexaey);
4168  aexdeyplus = Absolute(aexdey);
4169  aexceyplus = Absolute(aexcey);
4170  cexaeyplus = Absolute(cexaey);
4171  bexdeyplus = Absolute(bexdey);
4172  dexbeyplus = Absolute(dexbey);
4173  permanent = ((cexdeyplus + dexceyplus) * bezplus
4174  + (dexbeyplus + bexdeyplus) * cezplus
4175  + (bexceyplus + cexbeyplus) * dezplus)
4176  * alift
4177  + ((dexaeyplus + aexdeyplus) * cezplus
4178  + (aexceyplus + cexaeyplus) * dezplus
4179  + (cexdeyplus + dexceyplus) * aezplus)
4180  * blift
4181  + ((aexbeyplus + bexaeyplus) * dezplus
4182  + (bexdeyplus + dexbeyplus) * aezplus
4183  + (dexaeyplus + aexdeyplus) * bezplus)
4184  * clift
4185  + ((bexceyplus + cexbeyplus) * aezplus
4186  + (cexaeyplus + aexceyplus) * bezplus
4187  + (aexbeyplus + bexaeyplus) * cezplus)
4188  * dlift;
4189  errbound = isperrboundA * permanent;
4190  if ((det > errbound) || (-det > errbound)) {
4191  return det;
4192  }
4193 
4194  return insphereadapt(pa, pb, pc, pd, pe, permanent);
4195 }
4196 } //Added namespace to avoid clash with triangle