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