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

asr_approx_mvbb
Author(s): Gassner Nikolai
autogenerated on Mon Jun 10 2019 12:38:08