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