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 /**
995 * [JWP] Increment 'eindex,' but only update 'enow' if it is
996 * safe. Otherwise, 'enow' remains unchanged from the previous
997 * iteration...
998 */
999 ++eindex;
1000 if (eindex < elen)
1001 enow = e[eindex];
1002 } else {
1003 Two_Sum(Q, fnow, Qnew, hh);
1004 /**
1005 * [JWP] Increment 'findex,' but only update 'fnow' if it is
1006 * safe. Otherwise, 'fnow' remains unchanged from the previous
1007 * iteration...
1008 */
1009 ++findex;
1010 if (findex < flen)
1011 fnow = f[findex];
1012 }
1013 Q = Qnew;
1014 if (hh != 0.0) {
1015 h[hindex++] = hh;
1016 }
1017 }
1018 }
1019 while (eindex < elen) {
1020 Two_Sum(Q, enow, Qnew, hh);
1021 /**
1022 * [JWP] Increment 'eindex,' but only update 'enow' if it is
1023 * safe. Otherwise, 'enow' remains unchanged from the previous
1024 * iteration...
1025 */
1026 ++eindex;
1027 if (eindex < elen)
1028 enow = e[eindex];
1029 Q = Qnew;
1030 if (hh != 0.0) {
1031 h[hindex++] = hh;
1032 }
1033 }
1034 while (findex < flen) {
1035 Two_Sum(Q, fnow, Qnew, hh);
1036 /**
1037 * [JWP] Increment 'findex,' but only update 'fnow' if it is
1038 * safe. Otherwise, 'fnow' remains unchanged from the previous
1039 * iteration...
1040 */
1041 ++findex;
1042 if (findex < flen)
1043 fnow = f[findex];
1044 Q = Qnew;
1045 if (hh != 0.0) {
1046 h[hindex++] = hh;
1047 }
1048 }
1049 if ((Q != 0.0) || (hindex == 0)) {
1050 h[hindex++] = Q;
1051 }
1052 return hindex;
1053 }
1054
1055 /*****************************************************************************/
1056 /* */
1057 /* linear_expansion_sum() Sum two expansions. */
1058 /* */
1059 /* Sets h = e + f. See either version of my paper for details. */
1060 /* */
1061 /* Maintains the nonoverlapping property. (That is, if e is */
1062 /* nonoverlapping, h will be also.) */
1063 /* */
1064 /*****************************************************************************/
1065
linear_expansion_sum(int elen,REAL * e,int flen,REAL * f,REAL * h)1066 int linear_expansion_sum(int elen, REAL *e, int flen, REAL *f, REAL *h)
1067 /* h cannot be e or f. */
1068 {
1069 REAL Q, q;
1070 INEXACT REAL Qnew;
1071 INEXACT REAL R;
1072 INEXACT REAL bvirt;
1073 REAL avirt, bround, around;
1074 int eindex, findex, hindex;
1075 REAL enow, fnow;
1076 REAL g0;
1077
1078 enow = e[0];
1079 fnow = f[0];
1080 eindex = findex = 0;
1081 if ((fnow > enow) == (fnow > -enow)) {
1082 g0 = enow;
1083 enow = e[++eindex];
1084 } else {
1085 g0 = fnow;
1086 fnow = f[++findex];
1087 }
1088 if ((eindex < elen) && ((findex >= flen)
1089 || ((fnow > enow) == (fnow > -enow)))) {
1090 Fast_Two_Sum(enow, g0, Qnew, q);
1091 enow = e[++eindex];
1092 } else {
1093 Fast_Two_Sum(fnow, g0, Qnew, q);
1094 fnow = f[++findex];
1095 }
1096 Q = Qnew;
1097 for (hindex = 0; hindex < elen + flen - 2; hindex++) {
1098 if ((eindex < elen) && ((findex >= flen)
1099 || ((fnow > enow) == (fnow > -enow)))) {
1100 Fast_Two_Sum(enow, q, R, h[hindex]);
1101 enow = e[++eindex];
1102 } else {
1103 Fast_Two_Sum(fnow, q, R, h[hindex]);
1104 fnow = f[++findex];
1105 }
1106 Two_Sum(Q, R, Qnew, q);
1107 Q = Qnew;
1108 }
1109 h[hindex] = q;
1110 h[hindex + 1] = Q;
1111 return hindex + 2;
1112 }
1113
1114 /*****************************************************************************/
1115 /* */
1116 /* linear_expansion_sum_zeroelim() Sum two expansions, eliminating zero */
1117 /* components from the output expansion. */
1118 /* */
1119 /* Sets h = e + f. See either version of my paper for details. */
1120 /* */
1121 /* Maintains the nonoverlapping property. (That is, if e is */
1122 /* nonoverlapping, h will be also.) */
1123 /* */
1124 /*****************************************************************************/
1125
linear_expansion_sum_zeroelim(int elen,REAL * e,int flen,REAL * f,REAL * h)1126 int linear_expansion_sum_zeroelim(int elen, REAL *e, int flen, REAL *f,
1127 REAL *h)
1128 /* h cannot be e or f. */
1129 {
1130 REAL Q, q, hh;
1131 INEXACT REAL Qnew;
1132 INEXACT REAL R;
1133 INEXACT REAL bvirt;
1134 REAL avirt, bround, around;
1135 int eindex, findex, hindex;
1136 int count;
1137 REAL enow, fnow;
1138 REAL g0;
1139
1140 enow = e[0];
1141 fnow = f[0];
1142 eindex = findex = 0;
1143 hindex = 0;
1144 if ((fnow > enow) == (fnow > -enow)) {
1145 g0 = enow;
1146 enow = e[++eindex];
1147 } else {
1148 g0 = fnow;
1149 fnow = f[++findex];
1150 }
1151 if ((eindex < elen) && ((findex >= flen)
1152 || ((fnow > enow) == (fnow > -enow)))) {
1153 Fast_Two_Sum(enow, g0, Qnew, q);
1154 enow = e[++eindex];
1155 } else {
1156 Fast_Two_Sum(fnow, g0, Qnew, q);
1157 fnow = f[++findex];
1158 }
1159 Q = Qnew;
1160 for (count = 2; count < elen + flen; count++) {
1161 if ((eindex < elen) && ((findex >= flen)
1162 || ((fnow > enow) == (fnow > -enow)))) {
1163 Fast_Two_Sum(enow, q, R, hh);
1164 enow = e[++eindex];
1165 } else {
1166 Fast_Two_Sum(fnow, q, R, hh);
1167 fnow = f[++findex];
1168 }
1169 Two_Sum(Q, R, Qnew, q);
1170 Q = Qnew;
1171 if (hh != 0) {
1172 h[hindex++] = hh;
1173 }
1174 }
1175 if (q != 0) {
1176 h[hindex++] = q;
1177 }
1178 if ((Q != 0.0) || (hindex == 0)) {
1179 h[hindex++] = Q;
1180 }
1181 return hindex;
1182 }
1183
1184 /*****************************************************************************/
1185 /* */
1186 /* scale_expansion() Multiply an expansion by a scalar. */
1187 /* */
1188 /* Sets h = be. See either version of my paper for details. */
1189 /* */
1190 /* Maintains the nonoverlapping property. If round-to-even is used (as */
1191 /* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */
1192 /* properties as well. (That is, if e has one of these properties, so */
1193 /* will h.) */
1194 /* */
1195 /*****************************************************************************/
1196
scale_expansion(int elen,REAL * e,REAL b,REAL * h)1197 int scale_expansion(int elen, REAL *e, REAL b, REAL *h)
1198 /* e and h cannot be the same. */
1199 {
1200 INEXACT REAL Q;
1201 INEXACT REAL sum;
1202 INEXACT REAL product1;
1203 REAL product0;
1204 int eindex, hindex;
1205 REAL enow;
1206 INEXACT REAL bvirt;
1207 REAL avirt, bround, around;
1208 INEXACT REAL c;
1209 INEXACT REAL abig;
1210 REAL ahi, alo, bhi, blo;
1211 REAL err1, err2, err3;
1212
1213 Split(b, bhi, blo);
1214 Two_Product_Presplit(e[0], b, bhi, blo, Q, h[0]);
1215 hindex = 1;
1216 for (eindex = 1; eindex < elen; eindex++) {
1217 enow = e[eindex];
1218 Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
1219 Two_Sum(Q, product0, sum, h[hindex]);
1220 hindex++;
1221 Two_Sum(product1, sum, Q, h[hindex]);
1222 hindex++;
1223 }
1224 h[hindex] = Q;
1225 return elen + elen;
1226 }
1227
1228 /*****************************************************************************/
1229 /* */
1230 /* scale_expansion_zeroelim() Multiply an expansion by a scalar, */
1231 /* eliminating zero components from the */
1232 /* output expansion. */
1233 /* */
1234 /* Sets h = be. See either version of my paper for details. */
1235 /* */
1236 /* Maintains the nonoverlapping property. If round-to-even is used (as */
1237 /* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */
1238 /* properties as well. (That is, if e has one of these properties, so */
1239 /* will h.) */
1240 /* */
1241 /*****************************************************************************/
1242
scale_expansion_zeroelim(int elen,REAL * e,REAL b,REAL * h)1243 int scale_expansion_zeroelim(int elen, REAL *e, REAL b, REAL *h)
1244 /* e and h cannot be the same. */
1245 {
1246 INEXACT REAL Q, sum;
1247 REAL hh;
1248 INEXACT REAL product1;
1249 REAL product0;
1250 int eindex, hindex;
1251 REAL enow;
1252 INEXACT REAL bvirt;
1253 REAL avirt, bround, around;
1254 INEXACT REAL c;
1255 INEXACT REAL abig;
1256 REAL ahi, alo, bhi, blo;
1257 REAL err1, err2, err3;
1258
1259 Split(b, bhi, blo);
1260 Two_Product_Presplit(e[0], b, bhi, blo, Q, hh);
1261 hindex = 0;
1262 if (hh != 0) {
1263 h[hindex++] = hh;
1264 }
1265 for (eindex = 1; eindex < elen; eindex++) {
1266 enow = e[eindex];
1267 Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
1268 Two_Sum(Q, product0, sum, hh);
1269 if (hh != 0) {
1270 h[hindex++] = hh;
1271 }
1272 Fast_Two_Sum(product1, sum, Q, hh);
1273 if (hh != 0) {
1274 h[hindex++] = hh;
1275 }
1276 }
1277 if ((Q != 0.0) || (hindex == 0)) {
1278 h[hindex++] = Q;
1279 }
1280 return hindex;
1281 }
1282
1283 /*****************************************************************************/
1284 /* */
1285 /* compress() Compress an expansion. */
1286 /* */
1287 /* See the long version of my paper for details. */
1288 /* */
1289 /* Maintains the nonoverlapping property. If round-to-even is used (as */
1290 /* with IEEE 754), then any nonoverlapping expansion is converted to a */
1291 /* nonadjacent expansion. */
1292 /* */
1293 /*****************************************************************************/
1294
compress(int elen,REAL * e,REAL * h)1295 int compress(int elen, REAL *e, REAL *h)
1296 /* e and h may be the same. */
1297 {
1298 REAL Q, q;
1299 INEXACT REAL Qnew;
1300 int eindex, hindex;
1301 INEXACT REAL bvirt;
1302 REAL enow, hnow;
1303 int top, bottom;
1304
1305 bottom = elen - 1;
1306 Q = e[bottom];
1307 for (eindex = elen - 2; eindex >= 0; eindex--) {
1308 enow = e[eindex];
1309 Fast_Two_Sum(Q, enow, Qnew, q);
1310 if (q != 0) {
1311 h[bottom--] = Qnew;
1312 Q = q;
1313 } else {
1314 Q = Qnew;
1315 }
1316 }
1317 top = 0;
1318 for (hindex = bottom + 1; hindex < elen; hindex++) {
1319 hnow = h[hindex];
1320 Fast_Two_Sum(hnow, Q, Qnew, q);
1321 if (q != 0) {
1322 h[top++] = q;
1323 }
1324 Q = Qnew;
1325 }
1326 h[top] = Q;
1327 return top + 1;
1328 }
1329
1330 /*****************************************************************************/
1331 /* */
1332 /* estimate() Produce a one-word estimate of an expansion's value. */
1333 /* */
1334 /* See either version of my paper for details. */
1335 /* */
1336 /*****************************************************************************/
1337
estimate(int elen,REAL * e)1338 REAL estimate(int elen, REAL *e)
1339 {
1340 REAL Q;
1341 int eindex;
1342
1343 Q = e[0];
1344 for (eindex = 1; eindex < elen; eindex++) {
1345 Q += e[eindex];
1346 }
1347 return Q;
1348 }
1349
1350 /*****************************************************************************/
1351 /* */
1352 /* orient2dfast() Approximate 2D orientation test. Nonrobust. */
1353 /* orient2dexact() Exact 2D orientation test. Robust. */
1354 /* orient2dslow() Another exact 2D orientation test. Robust. */
1355 /* orient2d() Adaptive exact 2D orientation test. Robust. */
1356 /* */
1357 /* Return a positive value if the points pa, pb, and pc occur */
1358 /* in counterclockwise order; a negative value if they occur */
1359 /* in clockwise order; and zero if they are collinear. The */
1360 /* result is also a rough approximation of twice the signed */
1361 /* area of the triangle defined by the three points. */
1362 /* */
1363 /* Only the first and last routine should be used; the middle two are for */
1364 /* timings. */
1365 /* */
1366 /* The last three use exact arithmetic to ensure a correct answer. The */
1367 /* result returned is the determinant of a matrix. In orient2d() only, */
1368 /* this determinant is computed adaptively, in the sense that exact */
1369 /* arithmetic is used only to the degree it is needed to ensure that the */
1370 /* returned value has the correct sign. Hence, orient2d() is usually quite */
1371 /* fast, but will run more slowly when the input points are collinear or */
1372 /* nearly so. */
1373 /* */
1374 /*****************************************************************************/
1375
orient2dfast(REAL * pa,REAL * pb,REAL * pc)1376 REAL orient2dfast(REAL *pa, REAL *pb, REAL *pc)
1377 {
1378 REAL acx, bcx, acy, bcy;
1379
1380 acx = pa[0] - pc[0];
1381 bcx = pb[0] - pc[0];
1382 acy = pa[1] - pc[1];
1383 bcy = pb[1] - pc[1];
1384 return acx * bcy - acy * bcx;
1385 }
1386
orient2dexact(REAL * pa,REAL * pb,REAL * pc)1387 REAL orient2dexact(REAL *pa, REAL *pb, REAL *pc)
1388 {
1389 INEXACT REAL axby1, axcy1, bxcy1, bxay1, cxay1, cxby1;
1390 REAL axby0, axcy0, bxcy0, bxay0, cxay0, cxby0;
1391 REAL aterms[4], bterms[4], cterms[4];
1392 INEXACT REAL aterms3, bterms3, cterms3;
1393 REAL v[8], w[12];
1394 int vlength, wlength;
1395
1396 INEXACT REAL bvirt;
1397 REAL avirt, bround, around;
1398 INEXACT REAL c;
1399 INEXACT REAL abig;
1400 REAL ahi, alo, bhi, blo;
1401 REAL err1, err2, err3;
1402 INEXACT REAL _i, _j;
1403 REAL _0;
1404
1405 Two_Product(pa[0], pb[1], axby1, axby0);
1406 Two_Product(pa[0], pc[1], axcy1, axcy0);
1407 Two_Two_Diff(axby1, axby0, axcy1, axcy0,
1408 aterms3, aterms[2], aterms[1], aterms[0]);
1409 aterms[3] = aterms3;
1410
1411 Two_Product(pb[0], pc[1], bxcy1, bxcy0);
1412 Two_Product(pb[0], pa[1], bxay1, bxay0);
1413 Two_Two_Diff(bxcy1, bxcy0, bxay1, bxay0,
1414 bterms3, bterms[2], bterms[1], bterms[0]);
1415 bterms[3] = bterms3;
1416
1417 Two_Product(pc[0], pa[1], cxay1, cxay0);
1418 Two_Product(pc[0], pb[1], cxby1, cxby0);
1419 Two_Two_Diff(cxay1, cxay0, cxby1, cxby0,
1420 cterms3, cterms[2], cterms[1], cterms[0]);
1421 cterms[3] = cterms3;
1422
1423 vlength = fast_expansion_sum_zeroelim(4, aterms, 4, bterms, v);
1424 wlength = fast_expansion_sum_zeroelim(vlength, v, 4, cterms, w);
1425
1426 return w[wlength - 1];
1427 }
1428
orient2dslow(REAL * pa,REAL * pb,REAL * pc)1429 REAL orient2dslow(REAL *pa, REAL *pb, REAL *pc)
1430 {
1431 INEXACT REAL acx, acy, bcx, bcy;
1432 REAL acxtail, acytail;
1433 REAL bcxtail, bcytail;
1434 REAL negate, negatetail;
1435 REAL axby[8], bxay[8];
1436 INEXACT REAL axby7, bxay7;
1437 REAL deter[16];
1438 int deterlen;
1439
1440 INEXACT REAL bvirt;
1441 REAL avirt, bround, around;
1442 INEXACT REAL c;
1443 INEXACT REAL abig;
1444 REAL a0hi, a0lo, a1hi, a1lo, bhi, blo;
1445 REAL err1, err2, err3;
1446 INEXACT REAL _i, _j, _k, _l, _m, _n;
1447 REAL _0, _1, _2;
1448
1449 Two_Diff(pa[0], pc[0], acx, acxtail);
1450 Two_Diff(pa[1], pc[1], acy, acytail);
1451 Two_Diff(pb[0], pc[0], bcx, bcxtail);
1452 Two_Diff(pb[1], pc[1], bcy, bcytail);
1453
1454 Two_Two_Product(acx, acxtail, bcy, bcytail,
1455 axby7, axby[6], axby[5], axby[4],
1456 axby[3], axby[2], axby[1], axby[0]);
1457 axby[7] = axby7;
1458 negate = -acy;
1459 negatetail = -acytail;
1460 Two_Two_Product(bcx, bcxtail, negate, negatetail,
1461 bxay7, bxay[6], bxay[5], bxay[4],
1462 bxay[3], bxay[2], bxay[1], bxay[0]);
1463 bxay[7] = bxay7;
1464
1465 deterlen = fast_expansion_sum_zeroelim(8, axby, 8, bxay, deter);
1466
1467 return deter[deterlen - 1];
1468 }
1469
orient2dadapt(REAL * pa,REAL * pb,REAL * pc,REAL detsum)1470 REAL orient2dadapt(REAL *pa, REAL *pb, REAL *pc, REAL detsum)
1471 {
1472 INEXACT REAL acx, acy, bcx, bcy;
1473 REAL acxtail, acytail, bcxtail, bcytail;
1474 INEXACT REAL detleft, detright;
1475 REAL detlefttail, detrighttail;
1476 REAL det, errbound;
1477 REAL B[4], C1[8], C2[12], D[16];
1478 INEXACT REAL B3;
1479 int C1length, C2length, Dlength;
1480 REAL u[4];
1481 INEXACT REAL u3;
1482 INEXACT REAL s1, t1;
1483 REAL s0, t0;
1484
1485 INEXACT REAL bvirt;
1486 REAL avirt, bround, around;
1487 INEXACT REAL c;
1488 INEXACT REAL abig;
1489 REAL ahi, alo, bhi, blo;
1490 REAL err1, err2, err3;
1491 INEXACT REAL _i, _j;
1492 REAL _0;
1493
1494 acx = (REAL) (pa[0] - pc[0]);
1495 bcx = (REAL) (pb[0] - pc[0]);
1496 acy = (REAL) (pa[1] - pc[1]);
1497 bcy = (REAL) (pb[1] - pc[1]);
1498
1499 Two_Product(acx, bcy, detleft, detlefttail);
1500 Two_Product(acy, bcx, detright, detrighttail);
1501
1502 Two_Two_Diff(detleft, detlefttail, detright, detrighttail,
1503 B3, B[2], B[1], B[0]);
1504 B[3] = B3;
1505
1506 det = estimate(4, B);
1507 errbound = ccwerrboundB * detsum;
1508 if ((det >= errbound) || (-det >= errbound)) {
1509 return det;
1510 }
1511
1512 Two_Diff_Tail(pa[0], pc[0], acx, acxtail);
1513 Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail);
1514 Two_Diff_Tail(pa[1], pc[1], acy, acytail);
1515 Two_Diff_Tail(pb[1], pc[1], bcy, bcytail);
1516
1517 if ((acxtail == 0.0) && (acytail == 0.0)
1518 && (bcxtail == 0.0) && (bcytail == 0.0)) {
1519 return det;
1520 }
1521
1522 errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det);
1523 det += (acx * bcytail + bcy * acxtail)
1524 - (acy * bcxtail + bcx * acytail);
1525 if ((det >= errbound) || (-det >= errbound)) {
1526 return det;
1527 }
1528
1529 Two_Product(acxtail, bcy, s1, s0);
1530 Two_Product(acytail, bcx, t1, t0);
1531 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
1532 u[3] = u3;
1533 C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1);
1534
1535 Two_Product(acx, bcytail, s1, s0);
1536 Two_Product(acy, bcxtail, t1, t0);
1537 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
1538 u[3] = u3;
1539 C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2);
1540
1541 Two_Product(acxtail, bcytail, s1, s0);
1542 Two_Product(acytail, bcxtail, t1, t0);
1543 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
1544 u[3] = u3;
1545 Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D);
1546
1547 return(D[Dlength - 1]);
1548 }
1549
orient2d(REAL * pa,REAL * pb,REAL * pc)1550 REAL orient2d(REAL *pa, REAL *pb, REAL *pc)
1551 {
1552 REAL detleft, detright, det;
1553 REAL detsum, errbound;
1554
1555 detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]);
1556 detright = (pa[1] - pc[1]) * (pb[0] - pc[0]);
1557 det = detleft - detright;
1558
1559 if (detleft > 0.0) {
1560 if (detright <= 0.0) {
1561 return det;
1562 } else {
1563 detsum = detleft + detright;
1564 }
1565 } else if (detleft < 0.0) {
1566 if (detright >= 0.0) {
1567 return det;
1568 } else {
1569 detsum = -detleft - detright;
1570 }
1571 } else {
1572 return det;
1573 }
1574
1575 errbound = ccwerrboundA * detsum;
1576 if ((det >= errbound) || (-det >= errbound)) {
1577 return det;
1578 }
1579
1580 return orient2dadapt(pa, pb, pc, detsum);
1581 }
1582
1583 /*****************************************************************************/
1584 /* */
1585 /* orient3dfast() Approximate 3D orientation test. Nonrobust. */
1586 /* orient3dexact() Exact 3D orientation test. Robust. */
1587 /* orient3dslow() Another exact 3D orientation test. Robust. */
1588 /* orient3d() Adaptive exact 3D orientation test. Robust. */
1589 /* */
1590 /* Return a positive value if the point pd lies below the */
1591 /* plane passing through pa, pb, and pc; "below" is defined so */
1592 /* that pa, pb, and pc appear in counterclockwise order when */
1593 /* viewed from above the plane. Returns a negative value if */
1594 /* pd lies above the plane. Returns zero if the points are */
1595 /* coplanar. The result is also a rough approximation of six */
1596 /* times the signed volume of the tetrahedron defined by the */
1597 /* four points. */
1598 /* */
1599 /* Only the first and last routine should be used; the middle two are for */
1600 /* timings. */
1601 /* */
1602 /* The last three use exact arithmetic to ensure a correct answer. The */
1603 /* result returned is the determinant of a matrix. In orient3d() only, */
1604 /* this determinant is computed adaptively, in the sense that exact */
1605 /* arithmetic is used only to the degree it is needed to ensure that the */
1606 /* returned value has the correct sign. Hence, orient3d() is usually quite */
1607 /* fast, but will run more slowly when the input points are coplanar or */
1608 /* nearly so. */
1609 /* */
1610 /*****************************************************************************/
1611
orient3dfast(REAL * pa,REAL * pb,REAL * pc,REAL * pd)1612 REAL orient3dfast(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
1613 {
1614 REAL adx, bdx, cdx;
1615 REAL ady, bdy, cdy;
1616 REAL adz, bdz, cdz;
1617
1618 adx = pa[0] - pd[0];
1619 bdx = pb[0] - pd[0];
1620 cdx = pc[0] - pd[0];
1621 ady = pa[1] - pd[1];
1622 bdy = pb[1] - pd[1];
1623 cdy = pc[1] - pd[1];
1624 adz = pa[2] - pd[2];
1625 bdz = pb[2] - pd[2];
1626 cdz = pc[2] - pd[2];
1627
1628 return adx * (bdy * cdz - bdz * cdy)
1629 + bdx * (cdy * adz - cdz * ady)
1630 + cdx * (ady * bdz - adz * bdy);
1631 }
1632
orient3dexact(REAL * pa,REAL * pb,REAL * pc,REAL * pd)1633 REAL orient3dexact(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
1634 {
1635 INEXACT REAL axby1, bxcy1, cxdy1, dxay1, axcy1, bxdy1;
1636 INEXACT REAL bxay1, cxby1, dxcy1, axdy1, cxay1, dxby1;
1637 REAL axby0, bxcy0, cxdy0, dxay0, axcy0, bxdy0;
1638 REAL bxay0, cxby0, dxcy0, axdy0, cxay0, dxby0;
1639 REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4];
1640 REAL temp8[8];
1641 int templen;
1642 REAL abc[12], bcd[12], cda[12], dab[12];
1643 int abclen, bcdlen, cdalen, dablen;
1644 REAL adet[24], bdet[24], cdet[24], ddet[24];
1645 int alen, blen, clen, dlen;
1646 REAL abdet[48], cddet[48];
1647 int ablen, cdlen;
1648 REAL deter[96];
1649 int deterlen;
1650 int i;
1651
1652 INEXACT REAL bvirt;
1653 REAL avirt, bround, around;
1654 INEXACT REAL c;
1655 INEXACT REAL abig;
1656 REAL ahi, alo, bhi, blo;
1657 REAL err1, err2, err3;
1658 INEXACT REAL _i, _j;
1659 REAL _0;
1660
1661 Two_Product(pa[0], pb[1], axby1, axby0);
1662 Two_Product(pb[0], pa[1], bxay1, bxay0);
1663 Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]);
1664
1665 Two_Product(pb[0], pc[1], bxcy1, bxcy0);
1666 Two_Product(pc[0], pb[1], cxby1, cxby0);
1667 Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]);
1668
1669 Two_Product(pc[0], pd[1], cxdy1, cxdy0);
1670 Two_Product(pd[0], pc[1], dxcy1, dxcy0);
1671 Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]);
1672
1673 Two_Product(pd[0], pa[1], dxay1, dxay0);
1674 Two_Product(pa[0], pd[1], axdy1, axdy0);
1675 Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]);
1676
1677 Two_Product(pa[0], pc[1], axcy1, axcy0);
1678 Two_Product(pc[0], pa[1], cxay1, cxay0);
1679 Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]);
1680
1681 Two_Product(pb[0], pd[1], bxdy1, bxdy0);
1682 Two_Product(pd[0], pb[1], dxby1, dxby0);
1683 Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]);
1684
1685 templen = fast_expansion_sum_zeroelim(4, cd, 4, da, temp8);
1686 cdalen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, cda);
1687 templen = fast_expansion_sum_zeroelim(4, da, 4, ab, temp8);
1688 dablen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, dab);
1689 for (i = 0; i < 4; i++) {
1690 bd[i] = -bd[i];
1691 ac[i] = -ac[i];
1692 }
1693 templen = fast_expansion_sum_zeroelim(4, ab, 4, bc, temp8);
1694 abclen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, abc);
1695 templen = fast_expansion_sum_zeroelim(4, bc, 4, cd, temp8);
1696 bcdlen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, bcd);
1697
1698 alen = scale_expansion_zeroelim(bcdlen, bcd, pa[2], adet);
1699 blen = scale_expansion_zeroelim(cdalen, cda, -pb[2], bdet);
1700 clen = scale_expansion_zeroelim(dablen, dab, pc[2], cdet);
1701 dlen = scale_expansion_zeroelim(abclen, abc, -pd[2], ddet);
1702
1703 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
1704 cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
1705 deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, deter);
1706
1707 return deter[deterlen - 1];
1708 }
1709
orient3dslow(REAL * pa,REAL * pb,REAL * pc,REAL * pd)1710 REAL orient3dslow(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
1711 {
1712 INEXACT REAL adx, ady, adz, bdx, bdy, bdz, cdx, cdy, cdz;
1713 REAL adxtail, adytail, adztail;
1714 REAL bdxtail, bdytail, bdztail;
1715 REAL cdxtail, cdytail, cdztail;
1716 REAL negate, negatetail;
1717 INEXACT REAL axby7, bxcy7, axcy7, bxay7, cxby7, cxay7;
1718 REAL axby[8], bxcy[8], axcy[8], bxay[8], cxby[8], cxay[8];
1719 REAL temp16[16], temp32[32], temp32t[32];
1720 int temp16len, temp32len, temp32tlen;
1721 REAL adet[64], bdet[64], cdet[64];
1722 int alen, blen, clen;
1723 REAL abdet[128];
1724 int ablen;
1725 REAL deter[192];
1726 int deterlen;
1727
1728 INEXACT REAL bvirt;
1729 REAL avirt, bround, around;
1730 INEXACT REAL c;
1731 INEXACT REAL abig;
1732 REAL a0hi, a0lo, a1hi, a1lo, bhi, blo;
1733 REAL err1, err2, err3;
1734 INEXACT REAL _i, _j, _k, _l, _m, _n;
1735 REAL _0, _1, _2;
1736
1737 Two_Diff(pa[0], pd[0], adx, adxtail);
1738 Two_Diff(pa[1], pd[1], ady, adytail);
1739 Two_Diff(pa[2], pd[2], adz, adztail);
1740 Two_Diff(pb[0], pd[0], bdx, bdxtail);
1741 Two_Diff(pb[1], pd[1], bdy, bdytail);
1742 Two_Diff(pb[2], pd[2], bdz, bdztail);
1743 Two_Diff(pc[0], pd[0], cdx, cdxtail);
1744 Two_Diff(pc[1], pd[1], cdy, cdytail);
1745 Two_Diff(pc[2], pd[2], cdz, cdztail);
1746
1747 Two_Two_Product(adx, adxtail, bdy, bdytail,
1748 axby7, axby[6], axby[5], axby[4],
1749 axby[3], axby[2], axby[1], axby[0]);
1750 axby[7] = axby7;
1751 negate = -ady;
1752 negatetail = -adytail;
1753 Two_Two_Product(bdx, bdxtail, negate, negatetail,
1754 bxay7, bxay[6], bxay[5], bxay[4],
1755 bxay[3], bxay[2], bxay[1], bxay[0]);
1756 bxay[7] = bxay7;
1757 Two_Two_Product(bdx, bdxtail, cdy, cdytail,
1758 bxcy7, bxcy[6], bxcy[5], bxcy[4],
1759 bxcy[3], bxcy[2], bxcy[1], bxcy[0]);
1760 bxcy[7] = bxcy7;
1761 negate = -bdy;
1762 negatetail = -bdytail;
1763 Two_Two_Product(cdx, cdxtail, negate, negatetail,
1764 cxby7, cxby[6], cxby[5], cxby[4],
1765 cxby[3], cxby[2], cxby[1], cxby[0]);
1766 cxby[7] = cxby7;
1767 Two_Two_Product(cdx, cdxtail, ady, adytail,
1768 cxay7, cxay[6], cxay[5], cxay[4],
1769 cxay[3], cxay[2], cxay[1], cxay[0]);
1770 cxay[7] = cxay7;
1771 negate = -cdy;
1772 negatetail = -cdytail;
1773 Two_Two_Product(adx, adxtail, negate, negatetail,
1774 axcy7, axcy[6], axcy[5], axcy[4],
1775 axcy[3], axcy[2], axcy[1], axcy[0]);
1776 axcy[7] = axcy7;
1777
1778 temp16len = fast_expansion_sum_zeroelim(8, bxcy, 8, cxby, temp16);
1779 temp32len = scale_expansion_zeroelim(temp16len, temp16, adz, temp32);
1780 temp32tlen = scale_expansion_zeroelim(temp16len, temp16, adztail, temp32t);
1781 alen = fast_expansion_sum_zeroelim(temp32len, temp32, temp32tlen, temp32t,
1782 adet);
1783
1784 temp16len = fast_expansion_sum_zeroelim(8, cxay, 8, axcy, temp16);
1785 temp32len = scale_expansion_zeroelim(temp16len, temp16, bdz, temp32);
1786 temp32tlen = scale_expansion_zeroelim(temp16len, temp16, bdztail, temp32t);
1787 blen = fast_expansion_sum_zeroelim(temp32len, temp32, temp32tlen, temp32t,
1788 bdet);
1789
1790 temp16len = fast_expansion_sum_zeroelim(8, axby, 8, bxay, temp16);
1791 temp32len = scale_expansion_zeroelim(temp16len, temp16, cdz, temp32);
1792 temp32tlen = scale_expansion_zeroelim(temp16len, temp16, cdztail, temp32t);
1793 clen = fast_expansion_sum_zeroelim(temp32len, temp32, temp32tlen, temp32t,
1794 cdet);
1795
1796 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
1797 deterlen = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, deter);
1798
1799 return deter[deterlen - 1];
1800 }
1801
orient3dadapt(REAL * pa,REAL * pb,REAL * pc,REAL * pd,REAL permanent)1802 REAL orient3dadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL permanent)
1803 {
1804 INEXACT REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
1805 REAL det, errbound;
1806
1807 INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
1808 REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
1809 REAL bc[4], ca[4], ab[4];
1810 INEXACT REAL bc3, ca3, ab3;
1811 REAL adet[8], bdet[8], cdet[8];
1812 int alen, blen, clen;
1813 REAL abdet[16];
1814 int ablen;
1815 REAL *finnow, *finother, *finswap;
1816 REAL fin1[192], fin2[192];
1817 int finlength;
1818
1819
1820 REAL adxtail, bdxtail, cdxtail;
1821 REAL adytail, bdytail, cdytail;
1822 REAL adztail, bdztail, cdztail;
1823 INEXACT REAL at_blarge, at_clarge;
1824 INEXACT REAL bt_clarge, bt_alarge;
1825 INEXACT REAL ct_alarge, ct_blarge;
1826 REAL at_b[4], at_c[4], bt_c[4], bt_a[4], ct_a[4], ct_b[4];
1827 int at_blen, at_clen, bt_clen, bt_alen, ct_alen, ct_blen;
1828 INEXACT REAL bdxt_cdy1, cdxt_bdy1, cdxt_ady1;
1829 INEXACT REAL adxt_cdy1, adxt_bdy1, bdxt_ady1;
1830 REAL bdxt_cdy0, cdxt_bdy0, cdxt_ady0;
1831 REAL adxt_cdy0, adxt_bdy0, bdxt_ady0;
1832 INEXACT REAL bdyt_cdx1, cdyt_bdx1, cdyt_adx1;
1833 INEXACT REAL adyt_cdx1, adyt_bdx1, bdyt_adx1;
1834 REAL bdyt_cdx0, cdyt_bdx0, cdyt_adx0;
1835 REAL adyt_cdx0, adyt_bdx0, bdyt_adx0;
1836 REAL bct[8], cat[8], abt[8];
1837 int bctlen, catlen, abtlen;
1838 INEXACT REAL bdxt_cdyt1, cdxt_bdyt1, cdxt_adyt1;
1839 INEXACT REAL adxt_cdyt1, adxt_bdyt1, bdxt_adyt1;
1840 REAL bdxt_cdyt0, cdxt_bdyt0, cdxt_adyt0;
1841 REAL adxt_cdyt0, adxt_bdyt0, bdxt_adyt0;
1842 REAL u[4], v[12], w[16];
1843 INEXACT REAL u3;
1844 int vlength, wlength;
1845 REAL negate;
1846
1847 INEXACT REAL bvirt;
1848 REAL avirt, bround, around;
1849 INEXACT REAL c;
1850 INEXACT REAL abig;
1851 REAL ahi, alo, bhi, blo;
1852 REAL err1, err2, err3;
1853 INEXACT REAL _i, _j, _k;
1854 REAL _0;
1855
1856
1857 adx = (REAL) (pa[0] - pd[0]);
1858 bdx = (REAL) (pb[0] - pd[0]);
1859 cdx = (REAL) (pc[0] - pd[0]);
1860 ady = (REAL) (pa[1] - pd[1]);
1861 bdy = (REAL) (pb[1] - pd[1]);
1862 cdy = (REAL) (pc[1] - pd[1]);
1863 adz = (REAL) (pa[2] - pd[2]);
1864 bdz = (REAL) (pb[2] - pd[2]);
1865 cdz = (REAL) (pc[2] - pd[2]);
1866
1867 Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
1868 Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
1869 Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
1870 bc[3] = bc3;
1871 alen = scale_expansion_zeroelim(4, bc, adz, adet);
1872
1873 Two_Product(cdx, ady, cdxady1, cdxady0);
1874 Two_Product(adx, cdy, adxcdy1, adxcdy0);
1875 Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
1876 ca[3] = ca3;
1877 blen = scale_expansion_zeroelim(4, ca, bdz, bdet);
1878
1879 Two_Product(adx, bdy, adxbdy1, adxbdy0);
1880 Two_Product(bdx, ady, bdxady1, bdxady0);
1881 Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
1882 ab[3] = ab3;
1883 clen = scale_expansion_zeroelim(4, ab, cdz, cdet);
1884
1885 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
1886 finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
1887
1888 det = estimate(finlength, fin1);
1889 errbound = o3derrboundB * permanent;
1890 if ((det >= errbound) || (-det >= errbound)) {
1891 return det;
1892 }
1893
1894 Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
1895 Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
1896 Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
1897 Two_Diff_Tail(pa[1], pd[1], ady, adytail);
1898 Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
1899 Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
1900 Two_Diff_Tail(pa[2], pd[2], adz, adztail);
1901 Two_Diff_Tail(pb[2], pd[2], bdz, bdztail);
1902 Two_Diff_Tail(pc[2], pd[2], cdz, cdztail);
1903
1904 if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
1905 && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)
1906 && (adztail == 0.0) && (bdztail == 0.0) && (cdztail == 0.0)) {
1907 return det;
1908 }
1909
1910 errbound = o3derrboundC * permanent + resulterrbound * Absolute(det);
1911 det += (adz * ((bdx * cdytail + cdy * bdxtail)
1912 - (bdy * cdxtail + cdx * bdytail))
1913 + adztail * (bdx * cdy - bdy * cdx))
1914 + (bdz * ((cdx * adytail + ady * cdxtail)
1915 - (cdy * adxtail + adx * cdytail))
1916 + bdztail * (cdx * ady - cdy * adx))
1917 + (cdz * ((adx * bdytail + bdy * adxtail)
1918 - (ady * bdxtail + bdx * adytail))
1919 + cdztail * (adx * bdy - ady * bdx));
1920 if ((det >= errbound) || (-det >= errbound)) {
1921 return det;
1922 }
1923
1924 finnow = fin1;
1925 finother = fin2;
1926
1927 if (adxtail == 0.0) {
1928 if (adytail == 0.0) {
1929 at_b[0] = 0.0;
1930 at_blen = 1;
1931 at_c[0] = 0.0;
1932 at_clen = 1;
1933 } else {
1934 negate = -adytail;
1935 Two_Product(negate, bdx, at_blarge, at_b[0]);
1936 at_b[1] = at_blarge;
1937 at_blen = 2;
1938 Two_Product(adytail, cdx, at_clarge, at_c[0]);
1939 at_c[1] = at_clarge;
1940 at_clen = 2;
1941 }
1942 } else {
1943 if (adytail == 0.0) {
1944 Two_Product(adxtail, bdy, at_blarge, at_b[0]);
1945 at_b[1] = at_blarge;
1946 at_blen = 2;
1947 negate = -adxtail;
1948 Two_Product(negate, cdy, at_clarge, at_c[0]);
1949 at_c[1] = at_clarge;
1950 at_clen = 2;
1951 } else {
1952 Two_Product(adxtail, bdy, adxt_bdy1, adxt_bdy0);
1953 Two_Product(adytail, bdx, adyt_bdx1, adyt_bdx0);
1954 Two_Two_Diff(adxt_bdy1, adxt_bdy0, adyt_bdx1, adyt_bdx0,
1955 at_blarge, at_b[2], at_b[1], at_b[0]);
1956 at_b[3] = at_blarge;
1957 at_blen = 4;
1958 Two_Product(adytail, cdx, adyt_cdx1, adyt_cdx0);
1959 Two_Product(adxtail, cdy, adxt_cdy1, adxt_cdy0);
1960 Two_Two_Diff(adyt_cdx1, adyt_cdx0, adxt_cdy1, adxt_cdy0,
1961 at_clarge, at_c[2], at_c[1], at_c[0]);
1962 at_c[3] = at_clarge;
1963 at_clen = 4;
1964 }
1965 }
1966 if (bdxtail == 0.0) {
1967 if (bdytail == 0.0) {
1968 bt_c[0] = 0.0;
1969 bt_clen = 1;
1970 bt_a[0] = 0.0;
1971 bt_alen = 1;
1972 } else {
1973 negate = -bdytail;
1974 Two_Product(negate, cdx, bt_clarge, bt_c[0]);
1975 bt_c[1] = bt_clarge;
1976 bt_clen = 2;
1977 Two_Product(bdytail, adx, bt_alarge, bt_a[0]);
1978 bt_a[1] = bt_alarge;
1979 bt_alen = 2;
1980 }
1981 } else {
1982 if (bdytail == 0.0) {
1983 Two_Product(bdxtail, cdy, bt_clarge, bt_c[0]);
1984 bt_c[1] = bt_clarge;
1985 bt_clen = 2;
1986 negate = -bdxtail;
1987 Two_Product(negate, ady, bt_alarge, bt_a[0]);
1988 bt_a[1] = bt_alarge;
1989 bt_alen = 2;
1990 } else {
1991 Two_Product(bdxtail, cdy, bdxt_cdy1, bdxt_cdy0);
1992 Two_Product(bdytail, cdx, bdyt_cdx1, bdyt_cdx0);
1993 Two_Two_Diff(bdxt_cdy1, bdxt_cdy0, bdyt_cdx1, bdyt_cdx0,
1994 bt_clarge, bt_c[2], bt_c[1], bt_c[0]);
1995 bt_c[3] = bt_clarge;
1996 bt_clen = 4;
1997 Two_Product(bdytail, adx, bdyt_adx1, bdyt_adx0);
1998 Two_Product(bdxtail, ady, bdxt_ady1, bdxt_ady0);
1999 Two_Two_Diff(bdyt_adx1, bdyt_adx0, bdxt_ady1, bdxt_ady0,
2000 bt_alarge, bt_a[2], bt_a[1], bt_a[0]);
2001 bt_a[3] = bt_alarge;
2002 bt_alen = 4;
2003 }
2004 }
2005 if (cdxtail == 0.0) {
2006 if (cdytail == 0.0) {
2007 ct_a[0] = 0.0;
2008 ct_alen = 1;
2009 ct_b[0] = 0.0;
2010 ct_blen = 1;
2011 } else {
2012 negate = -cdytail;
2013 Two_Product(negate, adx, ct_alarge, ct_a[0]);
2014 ct_a[1] = ct_alarge;
2015 ct_alen = 2;
2016 Two_Product(cdytail, bdx, ct_blarge, ct_b[0]);
2017 ct_b[1] = ct_blarge;
2018 ct_blen = 2;
2019 }
2020 } else {
2021 if (cdytail == 0.0) {
2022 Two_Product(cdxtail, ady, ct_alarge, ct_a[0]);
2023 ct_a[1] = ct_alarge;
2024 ct_alen = 2;
2025 negate = -cdxtail;
2026 Two_Product(negate, bdy, ct_blarge, ct_b[0]);
2027 ct_b[1] = ct_blarge;
2028 ct_blen = 2;
2029 } else {
2030 Two_Product(cdxtail, ady, cdxt_ady1, cdxt_ady0);
2031 Two_Product(cdytail, adx, cdyt_adx1, cdyt_adx0);
2032 Two_Two_Diff(cdxt_ady1, cdxt_ady0, cdyt_adx1, cdyt_adx0,
2033 ct_alarge, ct_a[2], ct_a[1], ct_a[0]);
2034 ct_a[3] = ct_alarge;
2035 ct_alen = 4;
2036 Two_Product(cdytail, bdx, cdyt_bdx1, cdyt_bdx0);
2037 Two_Product(cdxtail, bdy, cdxt_bdy1, cdxt_bdy0);
2038 Two_Two_Diff(cdyt_bdx1, cdyt_bdx0, cdxt_bdy1, cdxt_bdy0,
2039 ct_blarge, ct_b[2], ct_b[1], ct_b[0]);
2040 ct_b[3] = ct_blarge;
2041 ct_blen = 4;
2042 }
2043 }
2044
2045 bctlen = fast_expansion_sum_zeroelim(bt_clen, bt_c, ct_blen, ct_b, bct);
2046 wlength = scale_expansion_zeroelim(bctlen, bct, adz, w);
2047 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
2048 finother);
2049 finswap = finnow; finnow = finother; finother = finswap;
2050
2051 catlen = fast_expansion_sum_zeroelim(ct_alen, ct_a, at_clen, at_c, cat);
2052 wlength = scale_expansion_zeroelim(catlen, cat, bdz, w);
2053 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
2054 finother);
2055 finswap = finnow; finnow = finother; finother = finswap;
2056
2057 abtlen = fast_expansion_sum_zeroelim(at_blen, at_b, bt_alen, bt_a, abt);
2058 wlength = scale_expansion_zeroelim(abtlen, abt, cdz, w);
2059 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
2060 finother);
2061 finswap = finnow; finnow = finother; finother = finswap;
2062
2063 if (adztail != 0.0) {
2064 vlength = scale_expansion_zeroelim(4, bc, adztail, v);
2065 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
2066 finother);
2067 finswap = finnow; finnow = finother; finother = finswap;
2068 }
2069 if (bdztail != 0.0) {
2070 vlength = scale_expansion_zeroelim(4, ca, bdztail, v);
2071 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
2072 finother);
2073 finswap = finnow; finnow = finother; finother = finswap;
2074 }
2075 if (cdztail != 0.0) {
2076 vlength = scale_expansion_zeroelim(4, ab, cdztail, v);
2077 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
2078 finother);
2079 finswap = finnow; finnow = finother; finother = finswap;
2080 }
2081
2082 if (adxtail != 0.0) {
2083 if (bdytail != 0.0) {
2084 Two_Product(adxtail, bdytail, adxt_bdyt1, adxt_bdyt0);
2085 Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdz, u3, u[2], u[1], u[0]);
2086 u[3] = u3;
2087 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2088 finother);
2089 finswap = finnow; finnow = finother; finother = finswap;
2090 if (cdztail != 0.0) {
2091 Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdztail, u3, u[2], u[1], u[0]);
2092 u[3] = u3;
2093 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2094 finother);
2095 finswap = finnow; finnow = finother; finother = finswap;
2096 }
2097 }
2098 if (cdytail != 0.0) {
2099 negate = -adxtail;
2100 Two_Product(negate, cdytail, adxt_cdyt1, adxt_cdyt0);
2101 Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdz, u3, u[2], u[1], u[0]);
2102 u[3] = u3;
2103 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2104 finother);
2105 finswap = finnow; finnow = finother; finother = finswap;
2106 if (bdztail != 0.0) {
2107 Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdztail, u3, u[2], u[1], u[0]);
2108 u[3] = u3;
2109 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2110 finother);
2111 finswap = finnow; finnow = finother; finother = finswap;
2112 }
2113 }
2114 }
2115 if (bdxtail != 0.0) {
2116 if (cdytail != 0.0) {
2117 Two_Product(bdxtail, cdytail, bdxt_cdyt1, bdxt_cdyt0);
2118 Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adz, u3, u[2], u[1], u[0]);
2119 u[3] = u3;
2120 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2121 finother);
2122 finswap = finnow; finnow = finother; finother = finswap;
2123 if (adztail != 0.0) {
2124 Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adztail, u3, u[2], u[1], u[0]);
2125 u[3] = u3;
2126 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2127 finother);
2128 finswap = finnow; finnow = finother; finother = finswap;
2129 }
2130 }
2131 if (adytail != 0.0) {
2132 negate = -bdxtail;
2133 Two_Product(negate, adytail, bdxt_adyt1, bdxt_adyt0);
2134 Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdz, u3, u[2], u[1], u[0]);
2135 u[3] = u3;
2136 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2137 finother);
2138 finswap = finnow; finnow = finother; finother = finswap;
2139 if (cdztail != 0.0) {
2140 Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdztail, u3, u[2], u[1], u[0]);
2141 u[3] = u3;
2142 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2143 finother);
2144 finswap = finnow; finnow = finother; finother = finswap;
2145 }
2146 }
2147 }
2148 if (cdxtail != 0.0) {
2149 if (adytail != 0.0) {
2150 Two_Product(cdxtail, adytail, cdxt_adyt1, cdxt_adyt0);
2151 Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdz, u3, u[2], u[1], u[0]);
2152 u[3] = u3;
2153 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2154 finother);
2155 finswap = finnow; finnow = finother; finother = finswap;
2156 if (bdztail != 0.0) {
2157 Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdztail, u3, u[2], u[1], u[0]);
2158 u[3] = u3;
2159 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2160 finother);
2161 finswap = finnow; finnow = finother; finother = finswap;
2162 }
2163 }
2164 if (bdytail != 0.0) {
2165 negate = -cdxtail;
2166 Two_Product(negate, bdytail, cdxt_bdyt1, cdxt_bdyt0);
2167 Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adz, u3, u[2], u[1], u[0]);
2168 u[3] = u3;
2169 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2170 finother);
2171 finswap = finnow; finnow = finother; finother = finswap;
2172 if (adztail != 0.0) {
2173 Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adztail, u3, u[2], u[1], u[0]);
2174 u[3] = u3;
2175 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
2176 finother);
2177 finswap = finnow; finnow = finother; finother = finswap;
2178 }
2179 }
2180 }
2181
2182 if (adztail != 0.0) {
2183 wlength = scale_expansion_zeroelim(bctlen, bct, adztail, w);
2184 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
2185 finother);
2186 finswap = finnow; finnow = finother; finother = finswap;
2187 }
2188 if (bdztail != 0.0) {
2189 wlength = scale_expansion_zeroelim(catlen, cat, bdztail, w);
2190 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
2191 finother);
2192 finswap = finnow; finnow = finother; finother = finswap;
2193 }
2194 if (cdztail != 0.0) {
2195 wlength = scale_expansion_zeroelim(abtlen, abt, cdztail, w);
2196 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
2197 finother);
2198 finswap = finnow; finnow = finother; finother = finswap;
2199 }
2200
2201 return finnow[finlength - 1];
2202 }
2203
2204 #ifdef USE_CGAL_PREDICATES
2205
orient3d(REAL * pa,REAL * pb,REAL * pc,REAL * pd)2206 REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
2207 {
2208 return (REAL)
2209 - cgal_pred_obj.orientation_3_object()
2210 (Point(pa[0], pa[1], pa[2]),
2211 Point(pb[0], pb[1], pb[2]),
2212 Point(pc[0], pc[1], pc[2]),
2213 Point(pd[0], pd[1], pd[2]));
2214 }
2215
2216 #else
2217
orient3d(REAL * pa,REAL * pb,REAL * pc,REAL * pd)2218 REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
2219 {
2220 REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
2221 REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
2222 REAL det;
2223
2224
2225 adx = pa[0] - pd[0];
2226 ady = pa[1] - pd[1];
2227 adz = pa[2] - pd[2];
2228 bdx = pb[0] - pd[0];
2229 bdy = pb[1] - pd[1];
2230 bdz = pb[2] - pd[2];
2231 cdx = pc[0] - pd[0];
2232 cdy = pc[1] - pd[1];
2233 cdz = pc[2] - pd[2];
2234
2235 bdxcdy = bdx * cdy;
2236 cdxbdy = cdx * bdy;
2237
2238 cdxady = cdx * ady;
2239 adxcdy = adx * cdy;
2240
2241 adxbdy = adx * bdy;
2242 bdxady = bdx * ady;
2243
2244 det = adz * (bdxcdy - cdxbdy)
2245 + bdz * (cdxady - adxcdy)
2246 + cdz * (adxbdy - bdxady);
2247
2248 if (_use_inexact_arith) {
2249 return det;
2250 }
2251
2252 if (_use_static_filter) {
2253 //if (fabs(det) > o3dstaticfilter) return det;
2254 if (det > o3dstaticfilter) return det;
2255 if (det < -o3dstaticfilter) return det;
2256 }
2257
2258
2259 REAL permanent, errbound;
2260
2261 permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * Absolute(adz)
2262 + (Absolute(cdxady) + Absolute(adxcdy)) * Absolute(bdz)
2263 + (Absolute(adxbdy) + Absolute(bdxady)) * Absolute(cdz);
2264 errbound = o3derrboundA * permanent;
2265 if ((det > errbound) || (-det > errbound)) {
2266 return det;
2267 }
2268
2269 return orient3dadapt(pa, pb, pc, pd, permanent);
2270 }
2271
2272 #endif // #ifdef USE_CGAL_PREDICATES
2273
2274 /*****************************************************************************/
2275 /* */
2276 /* incirclefast() Approximate 2D incircle test. Nonrobust. */
2277 /* incircleexact() Exact 2D incircle test. Robust. */
2278 /* incircleslow() Another exact 2D incircle test. Robust. */
2279 /* incircle() Adaptive exact 2D incircle test. Robust. */
2280 /* */
2281 /* Return a positive value if the point pd lies inside the */
2282 /* circle passing through pa, pb, and pc; a negative value if */
2283 /* it lies outside; and zero if the four points are cocircular.*/
2284 /* The points pa, pb, and pc must be in counterclockwise */
2285 /* order, or the sign of the result will be reversed. */
2286 /* */
2287 /* Only the first and last routine should be used; the middle two are for */
2288 /* timings. */
2289 /* */
2290 /* The last three use exact arithmetic to ensure a correct answer. The */
2291 /* result returned is the determinant of a matrix. In incircle() only, */
2292 /* this determinant is computed adaptively, in the sense that exact */
2293 /* arithmetic is used only to the degree it is needed to ensure that the */
2294 /* returned value has the correct sign. Hence, incircle() is usually quite */
2295 /* fast, but will run more slowly when the input points are cocircular or */
2296 /* nearly so. */
2297 /* */
2298 /*****************************************************************************/
2299
incirclefast(REAL * pa,REAL * pb,REAL * pc,REAL * pd)2300 REAL incirclefast(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
2301 {
2302 REAL adx, ady, bdx, bdy, cdx, cdy;
2303 REAL abdet, bcdet, cadet;
2304 REAL alift, blift, clift;
2305
2306 adx = pa[0] - pd[0];
2307 ady = pa[1] - pd[1];
2308 bdx = pb[0] - pd[0];
2309 bdy = pb[1] - pd[1];
2310 cdx = pc[0] - pd[0];
2311 cdy = pc[1] - pd[1];
2312
2313 abdet = adx * bdy - bdx * ady;
2314 bcdet = bdx * cdy - cdx * bdy;
2315 cadet = cdx * ady - adx * cdy;
2316 alift = adx * adx + ady * ady;
2317 blift = bdx * bdx + bdy * bdy;
2318 clift = cdx * cdx + cdy * cdy;
2319
2320 return alift * bcdet + blift * cadet + clift * abdet;
2321 }
2322
incircleexact(REAL * pa,REAL * pb,REAL * pc,REAL * pd)2323 REAL incircleexact(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
2324 {
2325 INEXACT REAL axby1, bxcy1, cxdy1, dxay1, axcy1, bxdy1;
2326 INEXACT REAL bxay1, cxby1, dxcy1, axdy1, cxay1, dxby1;
2327 REAL axby0, bxcy0, cxdy0, dxay0, axcy0, bxdy0;
2328 REAL bxay0, cxby0, dxcy0, axdy0, cxay0, dxby0;
2329 REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4];
2330 REAL temp8[8];
2331 int templen;
2332 REAL abc[12], bcd[12], cda[12], dab[12];
2333 int abclen, bcdlen, cdalen, dablen;
2334 REAL det24x[24], det24y[24], det48x[48], det48y[48];
2335 int xlen, ylen;
2336 REAL adet[96], bdet[96], cdet[96], ddet[96];
2337 int alen, blen, clen, dlen;
2338 REAL abdet[192], cddet[192];
2339 int ablen, cdlen;
2340 REAL deter[384];
2341 int deterlen;
2342 int i;
2343
2344 INEXACT REAL bvirt;
2345 REAL avirt, bround, around;
2346 INEXACT REAL c;
2347 INEXACT REAL abig;
2348 REAL ahi, alo, bhi, blo;
2349 REAL err1, err2, err3;
2350 INEXACT REAL _i, _j;
2351 REAL _0;
2352
2353 Two_Product(pa[0], pb[1], axby1, axby0);
2354 Two_Product(pb[0], pa[1], bxay1, bxay0);
2355 Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]);
2356
2357 Two_Product(pb[0], pc[1], bxcy1, bxcy0);
2358 Two_Product(pc[0], pb[1], cxby1, cxby0);
2359 Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]);
2360
2361 Two_Product(pc[0], pd[1], cxdy1, cxdy0);
2362 Two_Product(pd[0], pc[1], dxcy1, dxcy0);
2363 Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]);
2364
2365 Two_Product(pd[0], pa[1], dxay1, dxay0);
2366 Two_Product(pa[0], pd[1], axdy1, axdy0);
2367 Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]);
2368
2369 Two_Product(pa[0], pc[1], axcy1, axcy0);
2370 Two_Product(pc[0], pa[1], cxay1, cxay0);
2371 Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]);
2372
2373 Two_Product(pb[0], pd[1], bxdy1, bxdy0);
2374 Two_Product(pd[0], pb[1], dxby1, dxby0);
2375 Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]);
2376
2377 templen = fast_expansion_sum_zeroelim(4, cd, 4, da, temp8);
2378 cdalen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, cda);
2379 templen = fast_expansion_sum_zeroelim(4, da, 4, ab, temp8);
2380 dablen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, dab);
2381 for (i = 0; i < 4; i++) {
2382 bd[i] = -bd[i];
2383 ac[i] = -ac[i];
2384 }
2385 templen = fast_expansion_sum_zeroelim(4, ab, 4, bc, temp8);
2386 abclen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, abc);
2387 templen = fast_expansion_sum_zeroelim(4, bc, 4, cd, temp8);
2388 bcdlen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, bcd);
2389
2390 xlen = scale_expansion_zeroelim(bcdlen, bcd, pa[0], det24x);
2391 xlen = scale_expansion_zeroelim(xlen, det24x, pa[0], det48x);
2392 ylen = scale_expansion_zeroelim(bcdlen, bcd, pa[1], det24y);
2393 ylen = scale_expansion_zeroelim(ylen, det24y, pa[1], det48y);
2394 alen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, adet);
2395
2396 xlen = scale_expansion_zeroelim(cdalen, cda, pb[0], det24x);
2397 xlen = scale_expansion_zeroelim(xlen, det24x, -pb[0], det48x);
2398 ylen = scale_expansion_zeroelim(cdalen, cda, pb[1], det24y);
2399 ylen = scale_expansion_zeroelim(ylen, det24y, -pb[1], det48y);
2400 blen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, bdet);
2401
2402 xlen = scale_expansion_zeroelim(dablen, dab, pc[0], det24x);
2403 xlen = scale_expansion_zeroelim(xlen, det24x, pc[0], det48x);
2404 ylen = scale_expansion_zeroelim(dablen, dab, pc[1], det24y);
2405 ylen = scale_expansion_zeroelim(ylen, det24y, pc[1], det48y);
2406 clen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, cdet);
2407
2408 xlen = scale_expansion_zeroelim(abclen, abc, pd[0], det24x);
2409 xlen = scale_expansion_zeroelim(xlen, det24x, -pd[0], det48x);
2410 ylen = scale_expansion_zeroelim(abclen, abc, pd[1], det24y);
2411 ylen = scale_expansion_zeroelim(ylen, det24y, -pd[1], det48y);
2412 dlen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, ddet);
2413
2414 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
2415 cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
2416 deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, deter);
2417
2418 return deter[deterlen - 1];
2419 }
2420
incircleslow(REAL * pa,REAL * pb,REAL * pc,REAL * pd)2421 REAL incircleslow(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
2422 {
2423 INEXACT REAL adx, bdx, cdx, ady, bdy, cdy;
2424 REAL adxtail, bdxtail, cdxtail;
2425 REAL adytail, bdytail, cdytail;
2426 REAL negate, negatetail;
2427 INEXACT REAL axby7, bxcy7, axcy7, bxay7, cxby7, cxay7;
2428 REAL axby[8], bxcy[8], axcy[8], bxay[8], cxby[8], cxay[8];
2429 REAL temp16[16];
2430 int temp16len;
2431 REAL detx[32], detxx[64], detxt[32], detxxt[64], detxtxt[64];
2432 int xlen, xxlen, xtlen, xxtlen, xtxtlen;
2433 REAL x1[128], x2[192];
2434 int x1len, x2len;
2435 REAL dety[32], detyy[64], detyt[32], detyyt[64], detytyt[64];
2436 int ylen, yylen, ytlen, yytlen, ytytlen;
2437 REAL y1[128], y2[192];
2438 int y1len, y2len;
2439 REAL adet[384], bdet[384], cdet[384], abdet[768], deter[1152];
2440 int alen, blen, clen, ablen, deterlen;
2441 int i;
2442
2443 INEXACT REAL bvirt;
2444 REAL avirt, bround, around;
2445 INEXACT REAL c;
2446 INEXACT REAL abig;
2447 REAL a0hi, a0lo, a1hi, a1lo, bhi, blo;
2448 REAL err1, err2, err3;
2449 INEXACT REAL _i, _j, _k, _l, _m, _n;
2450 REAL _0, _1, _2;
2451
2452 Two_Diff(pa[0], pd[0], adx, adxtail);
2453 Two_Diff(pa[1], pd[1], ady, adytail);
2454 Two_Diff(pb[0], pd[0], bdx, bdxtail);
2455 Two_Diff(pb[1], pd[1], bdy, bdytail);
2456 Two_Diff(pc[0], pd[0], cdx, cdxtail);
2457 Two_Diff(pc[1], pd[1], cdy, cdytail);
2458
2459 Two_Two_Product(adx, adxtail, bdy, bdytail,
2460 axby7, axby[6], axby[5], axby[4],
2461 axby[3], axby[2], axby[1], axby[0]);
2462 axby[7] = axby7;
2463 negate = -ady;
2464 negatetail = -adytail;
2465 Two_Two_Product(bdx, bdxtail, negate, negatetail,
2466 bxay7, bxay[6], bxay[5], bxay[4],
2467 bxay[3], bxay[2], bxay[1], bxay[0]);
2468 bxay[7] = bxay7;
2469 Two_Two_Product(bdx, bdxtail, cdy, cdytail,
2470 bxcy7, bxcy[6], bxcy[5], bxcy[4],
2471 bxcy[3], bxcy[2], bxcy[1], bxcy[0]);
2472 bxcy[7] = bxcy7;
2473 negate = -bdy;
2474 negatetail = -bdytail;
2475 Two_Two_Product(cdx, cdxtail, negate, negatetail,
2476 cxby7, cxby[6], cxby[5], cxby[4],
2477 cxby[3], cxby[2], cxby[1], cxby[0]);
2478 cxby[7] = cxby7;
2479 Two_Two_Product(cdx, cdxtail, ady, adytail,
2480 cxay7, cxay[6], cxay[5], cxay[4],
2481 cxay[3], cxay[2], cxay[1], cxay[0]);
2482 cxay[7] = cxay7;
2483 negate = -cdy;
2484 negatetail = -cdytail;
2485 Two_Two_Product(adx, adxtail, negate, negatetail,
2486 axcy7, axcy[6], axcy[5], axcy[4],
2487 axcy[3], axcy[2], axcy[1], axcy[0]);
2488 axcy[7] = axcy7;
2489
2490
2491 temp16len = fast_expansion_sum_zeroelim(8, bxcy, 8, cxby, temp16);
2492
2493 xlen = scale_expansion_zeroelim(temp16len, temp16, adx, detx);
2494 xxlen = scale_expansion_zeroelim(xlen, detx, adx, detxx);
2495 xtlen = scale_expansion_zeroelim(temp16len, temp16, adxtail, detxt);
2496 xxtlen = scale_expansion_zeroelim(xtlen, detxt, adx, detxxt);
2497 for (i = 0; i < xxtlen; i++) {
2498 detxxt[i] *= 2.0;
2499 }
2500 xtxtlen = scale_expansion_zeroelim(xtlen, detxt, adxtail, detxtxt);
2501 x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
2502 x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
2503
2504 ylen = scale_expansion_zeroelim(temp16len, temp16, ady, dety);
2505 yylen = scale_expansion_zeroelim(ylen, dety, ady, detyy);
2506 ytlen = scale_expansion_zeroelim(temp16len, temp16, adytail, detyt);
2507 yytlen = scale_expansion_zeroelim(ytlen, detyt, ady, detyyt);
2508 for (i = 0; i < yytlen; i++) {
2509 detyyt[i] *= 2.0;
2510 }
2511 ytytlen = scale_expansion_zeroelim(ytlen, detyt, adytail, detytyt);
2512 y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
2513 y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
2514
2515 alen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, adet);
2516
2517
2518 temp16len = fast_expansion_sum_zeroelim(8, cxay, 8, axcy, temp16);
2519
2520 xlen = scale_expansion_zeroelim(temp16len, temp16, bdx, detx);
2521 xxlen = scale_expansion_zeroelim(xlen, detx, bdx, detxx);
2522 xtlen = scale_expansion_zeroelim(temp16len, temp16, bdxtail, detxt);
2523 xxtlen = scale_expansion_zeroelim(xtlen, detxt, bdx, detxxt);
2524 for (i = 0; i < xxtlen; i++) {
2525 detxxt[i] *= 2.0;
2526 }
2527 xtxtlen = scale_expansion_zeroelim(xtlen, detxt, bdxtail, detxtxt);
2528 x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
2529 x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
2530
2531 ylen = scale_expansion_zeroelim(temp16len, temp16, bdy, dety);
2532 yylen = scale_expansion_zeroelim(ylen, dety, bdy, detyy);
2533 ytlen = scale_expansion_zeroelim(temp16len, temp16, bdytail, detyt);
2534 yytlen = scale_expansion_zeroelim(ytlen, detyt, bdy, detyyt);
2535 for (i = 0; i < yytlen; i++) {
2536 detyyt[i] *= 2.0;
2537 }
2538 ytytlen = scale_expansion_zeroelim(ytlen, detyt, bdytail, detytyt);
2539 y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
2540 y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
2541
2542 blen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, bdet);
2543
2544
2545 temp16len = fast_expansion_sum_zeroelim(8, axby, 8, bxay, temp16);
2546
2547 xlen = scale_expansion_zeroelim(temp16len, temp16, cdx, detx);
2548 xxlen = scale_expansion_zeroelim(xlen, detx, cdx, detxx);
2549 xtlen = scale_expansion_zeroelim(temp16len, temp16, cdxtail, detxt);
2550 xxtlen = scale_expansion_zeroelim(xtlen, detxt, cdx, detxxt);
2551 for (i = 0; i < xxtlen; i++) {
2552 detxxt[i] *= 2.0;
2553 }
2554 xtxtlen = scale_expansion_zeroelim(xtlen, detxt, cdxtail, detxtxt);
2555 x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
2556 x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
2557
2558 ylen = scale_expansion_zeroelim(temp16len, temp16, cdy, dety);
2559 yylen = scale_expansion_zeroelim(ylen, dety, cdy, detyy);
2560 ytlen = scale_expansion_zeroelim(temp16len, temp16, cdytail, detyt);
2561 yytlen = scale_expansion_zeroelim(ytlen, detyt, cdy, detyyt);
2562 for (i = 0; i < yytlen; i++) {
2563 detyyt[i] *= 2.0;
2564 }
2565 ytytlen = scale_expansion_zeroelim(ytlen, detyt, cdytail, detytyt);
2566 y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
2567 y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
2568
2569 clen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, cdet);
2570
2571 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
2572 deterlen = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, deter);
2573
2574 return deter[deterlen - 1];
2575 }
2576
incircleadapt(REAL * pa,REAL * pb,REAL * pc,REAL * pd,REAL permanent)2577 REAL incircleadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL permanent)
2578 {
2579 INEXACT REAL adx, bdx, cdx, ady, bdy, cdy;
2580 REAL det, errbound;
2581
2582 INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
2583 REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
2584 REAL bc[4], ca[4], ab[4];
2585 INEXACT REAL bc3, ca3, ab3;
2586 REAL axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32];
2587 int axbclen, axxbclen, aybclen, ayybclen, alen;
2588 REAL bxca[8], bxxca[16], byca[8], byyca[16], bdet[32];
2589 int bxcalen, bxxcalen, bycalen, byycalen, blen;
2590 REAL cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32];
2591 int cxablen, cxxablen, cyablen, cyyablen, clen;
2592 REAL abdet[64];
2593 int ablen;
2594 REAL fin1[1152], fin2[1152];
2595 REAL *finnow, *finother, *finswap;
2596 int finlength;
2597
2598 REAL adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
2599 INEXACT REAL adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1;
2600 REAL adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0;
2601 REAL aa[4], bb[4], cc[4];
2602 INEXACT REAL aa3, bb3, cc3;
2603 INEXACT REAL ti1, tj1;
2604 REAL ti0, tj0;
2605 REAL u[4], v[4];
2606 INEXACT REAL u3, v3;
2607 REAL temp8[8], temp16a[16], temp16b[16], temp16c[16];
2608 REAL temp32a[32], temp32b[32], temp48[48], temp64[64];
2609 int temp8len, temp16alen, temp16blen, temp16clen;
2610 int temp32alen, temp32blen, temp48len, temp64len;
2611 REAL axtbb[8], axtcc[8], aytbb[8], aytcc[8];
2612 int axtbblen, axtcclen, aytbblen, aytcclen;
2613 REAL bxtaa[8], bxtcc[8], bytaa[8], bytcc[8];
2614 int bxtaalen, bxtcclen, bytaalen, bytcclen;
2615 REAL cxtaa[8], cxtbb[8], cytaa[8], cytbb[8];
2616 int cxtaalen, cxtbblen, cytaalen, cytbblen;
2617 REAL axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8];
2618 int axtbclen, aytbclen, bxtcalen, bytcalen, cxtablen, cytablen;
2619 REAL axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16];
2620 int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen;
2621 REAL axtbctt[8], aytbctt[8], bxtcatt[8];
2622 REAL bytcatt[8], cxtabtt[8], cytabtt[8];
2623 int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen;
2624 REAL abt[8], bct[8], cat[8];
2625 int abtlen, bctlen, catlen;
2626 REAL abtt[4], bctt[4], catt[4];
2627 int abttlen, bcttlen, cattlen;
2628 INEXACT REAL abtt3, bctt3, catt3;
2629 REAL negate;
2630
2631 INEXACT REAL bvirt;
2632 REAL avirt, bround, around;
2633 INEXACT REAL c;
2634 INEXACT REAL abig;
2635 REAL ahi, alo, bhi, blo;
2636 REAL err1, err2, err3;
2637 INEXACT REAL _i, _j;
2638 REAL _0;
2639
2640 // Avoid compiler warnings. H. Si, 2012-02-16.
2641 axtbclen = aytbclen = bxtcalen = bytcalen = cxtablen = cytablen = 0;
2642
2643 adx = (REAL) (pa[0] - pd[0]);
2644 bdx = (REAL) (pb[0] - pd[0]);
2645 cdx = (REAL) (pc[0] - pd[0]);
2646 ady = (REAL) (pa[1] - pd[1]);
2647 bdy = (REAL) (pb[1] - pd[1]);
2648 cdy = (REAL) (pc[1] - pd[1]);
2649
2650 Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
2651 Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
2652 Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
2653 bc[3] = bc3;
2654 axbclen = scale_expansion_zeroelim(4, bc, adx, axbc);
2655 axxbclen = scale_expansion_zeroelim(axbclen, axbc, adx, axxbc);
2656 aybclen = scale_expansion_zeroelim(4, bc, ady, aybc);
2657 ayybclen = scale_expansion_zeroelim(aybclen, aybc, ady, ayybc);
2658 alen = fast_expansion_sum_zeroelim(axxbclen, axxbc, ayybclen, ayybc, adet);
2659
2660 Two_Product(cdx, ady, cdxady1, cdxady0);
2661 Two_Product(adx, cdy, adxcdy1, adxcdy0);
2662 Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
2663 ca[3] = ca3;
2664 bxcalen = scale_expansion_zeroelim(4, ca, bdx, bxca);
2665 bxxcalen = scale_expansion_zeroelim(bxcalen, bxca, bdx, bxxca);
2666 bycalen = scale_expansion_zeroelim(4, ca, bdy, byca);
2667 byycalen = scale_expansion_zeroelim(bycalen, byca, bdy, byyca);
2668 blen = fast_expansion_sum_zeroelim(bxxcalen, bxxca, byycalen, byyca, bdet);
2669
2670 Two_Product(adx, bdy, adxbdy1, adxbdy0);
2671 Two_Product(bdx, ady, bdxady1, bdxady0);
2672 Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
2673 ab[3] = ab3;
2674 cxablen = scale_expansion_zeroelim(4, ab, cdx, cxab);
2675 cxxablen = scale_expansion_zeroelim(cxablen, cxab, cdx, cxxab);
2676 cyablen = scale_expansion_zeroelim(4, ab, cdy, cyab);
2677 cyyablen = scale_expansion_zeroelim(cyablen, cyab, cdy, cyyab);
2678 clen = fast_expansion_sum_zeroelim(cxxablen, cxxab, cyyablen, cyyab, cdet);
2679
2680 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
2681 finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
2682
2683 det = estimate(finlength, fin1);
2684 errbound = iccerrboundB * permanent;
2685 if ((det >= errbound) || (-det >= errbound)) {
2686 return det;
2687 }
2688
2689 Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
2690 Two_Diff_Tail(pa[1], pd[1], ady, adytail);
2691 Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
2692 Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
2693 Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
2694 Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
2695 if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
2696 && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)) {
2697 return det;
2698 }
2699
2700 errbound = iccerrboundC * permanent + resulterrbound * Absolute(det);
2701 det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail)
2702 - (bdy * cdxtail + cdx * bdytail))
2703 + 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx))
2704 + ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail)
2705 - (cdy * adxtail + adx * cdytail))
2706 + 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx))
2707 + ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail)
2708 - (ady * bdxtail + bdx * adytail))
2709 + 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));
2710 if ((det >= errbound) || (-det >= errbound)) {
2711 return det;
2712 }
2713
2714 finnow = fin1;
2715 finother = fin2;
2716
2717 if ((bdxtail != 0.0) || (bdytail != 0.0)
2718 || (cdxtail != 0.0) || (cdytail != 0.0)) {
2719 Square(adx, adxadx1, adxadx0);
2720 Square(ady, adyady1, adyady0);
2721 Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]);
2722 aa[3] = aa3;
2723 }
2724 if ((cdxtail != 0.0) || (cdytail != 0.0)
2725 || (adxtail != 0.0) || (adytail != 0.0)) {
2726 Square(bdx, bdxbdx1, bdxbdx0);
2727 Square(bdy, bdybdy1, bdybdy0);
2728 Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]);
2729 bb[3] = bb3;
2730 }
2731 if ((adxtail != 0.0) || (adytail != 0.0)
2732 || (bdxtail != 0.0) || (bdytail != 0.0)) {
2733 Square(cdx, cdxcdx1, cdxcdx0);
2734 Square(cdy, cdycdy1, cdycdy0);
2735 Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]);
2736 cc[3] = cc3;
2737 }
2738
2739 if (adxtail != 0.0) {
2740 axtbclen = scale_expansion_zeroelim(4, bc, adxtail, axtbc);
2741 temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, 2.0 * adx,
2742 temp16a);
2743
2744 axtcclen = scale_expansion_zeroelim(4, cc, adxtail, axtcc);
2745 temp16blen = scale_expansion_zeroelim(axtcclen, axtcc, bdy, temp16b);
2746
2747 axtbblen = scale_expansion_zeroelim(4, bb, adxtail, axtbb);
2748 temp16clen = scale_expansion_zeroelim(axtbblen, axtbb, -cdy, temp16c);
2749
2750 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2751 temp16blen, temp16b, temp32a);
2752 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
2753 temp32alen, temp32a, temp48);
2754 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2755 temp48, finother);
2756 finswap = finnow; finnow = finother; finother = finswap;
2757 }
2758 if (adytail != 0.0) {
2759 aytbclen = scale_expansion_zeroelim(4, bc, adytail, aytbc);
2760 temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, 2.0 * ady,
2761 temp16a);
2762
2763 aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb);
2764 temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b);
2765
2766 aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc);
2767 temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c);
2768
2769 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2770 temp16blen, temp16b, temp32a);
2771 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
2772 temp32alen, temp32a, temp48);
2773 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2774 temp48, finother);
2775 finswap = finnow; finnow = finother; finother = finswap;
2776 }
2777 if (bdxtail != 0.0) {
2778 bxtcalen = scale_expansion_zeroelim(4, ca, bdxtail, bxtca);
2779 temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, 2.0 * bdx,
2780 temp16a);
2781
2782 bxtaalen = scale_expansion_zeroelim(4, aa, bdxtail, bxtaa);
2783 temp16blen = scale_expansion_zeroelim(bxtaalen, bxtaa, cdy, temp16b);
2784
2785 bxtcclen = scale_expansion_zeroelim(4, cc, bdxtail, bxtcc);
2786 temp16clen = scale_expansion_zeroelim(bxtcclen, bxtcc, -ady, temp16c);
2787
2788 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2789 temp16blen, temp16b, temp32a);
2790 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
2791 temp32alen, temp32a, temp48);
2792 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2793 temp48, finother);
2794 finswap = finnow; finnow = finother; finother = finswap;
2795 }
2796 if (bdytail != 0.0) {
2797 bytcalen = scale_expansion_zeroelim(4, ca, bdytail, bytca);
2798 temp16alen = scale_expansion_zeroelim(bytcalen, bytca, 2.0 * bdy,
2799 temp16a);
2800
2801 bytcclen = scale_expansion_zeroelim(4, cc, bdytail, bytcc);
2802 temp16blen = scale_expansion_zeroelim(bytcclen, bytcc, adx, temp16b);
2803
2804 bytaalen = scale_expansion_zeroelim(4, aa, bdytail, bytaa);
2805 temp16clen = scale_expansion_zeroelim(bytaalen, bytaa, -cdx, temp16c);
2806
2807 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2808 temp16blen, temp16b, temp32a);
2809 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
2810 temp32alen, temp32a, temp48);
2811 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2812 temp48, finother);
2813 finswap = finnow; finnow = finother; finother = finswap;
2814 }
2815 if (cdxtail != 0.0) {
2816 cxtablen = scale_expansion_zeroelim(4, ab, cdxtail, cxtab);
2817 temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, 2.0 * cdx,
2818 temp16a);
2819
2820 cxtbblen = scale_expansion_zeroelim(4, bb, cdxtail, cxtbb);
2821 temp16blen = scale_expansion_zeroelim(cxtbblen, cxtbb, ady, temp16b);
2822
2823 cxtaalen = scale_expansion_zeroelim(4, aa, cdxtail, cxtaa);
2824 temp16clen = scale_expansion_zeroelim(cxtaalen, cxtaa, -bdy, temp16c);
2825
2826 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2827 temp16blen, temp16b, temp32a);
2828 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
2829 temp32alen, temp32a, temp48);
2830 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2831 temp48, finother);
2832 finswap = finnow; finnow = finother; finother = finswap;
2833 }
2834 if (cdytail != 0.0) {
2835 cytablen = scale_expansion_zeroelim(4, ab, cdytail, cytab);
2836 temp16alen = scale_expansion_zeroelim(cytablen, cytab, 2.0 * cdy,
2837 temp16a);
2838
2839 cytaalen = scale_expansion_zeroelim(4, aa, cdytail, cytaa);
2840 temp16blen = scale_expansion_zeroelim(cytaalen, cytaa, bdx, temp16b);
2841
2842 cytbblen = scale_expansion_zeroelim(4, bb, cdytail, cytbb);
2843 temp16clen = scale_expansion_zeroelim(cytbblen, cytbb, -adx, temp16c);
2844
2845 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2846 temp16blen, temp16b, temp32a);
2847 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
2848 temp32alen, temp32a, temp48);
2849 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2850 temp48, finother);
2851 finswap = finnow; finnow = finother; finother = finswap;
2852 }
2853
2854 if ((adxtail != 0.0) || (adytail != 0.0)) {
2855 if ((bdxtail != 0.0) || (bdytail != 0.0)
2856 || (cdxtail != 0.0) || (cdytail != 0.0)) {
2857 Two_Product(bdxtail, cdy, ti1, ti0);
2858 Two_Product(bdx, cdytail, tj1, tj0);
2859 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
2860 u[3] = u3;
2861 negate = -bdy;
2862 Two_Product(cdxtail, negate, ti1, ti0);
2863 negate = -bdytail;
2864 Two_Product(cdx, negate, tj1, tj0);
2865 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
2866 v[3] = v3;
2867 bctlen = fast_expansion_sum_zeroelim(4, u, 4, v, bct);
2868
2869 Two_Product(bdxtail, cdytail, ti1, ti0);
2870 Two_Product(cdxtail, bdytail, tj1, tj0);
2871 Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]);
2872 bctt[3] = bctt3;
2873 bcttlen = 4;
2874 } else {
2875 bct[0] = 0.0;
2876 bctlen = 1;
2877 bctt[0] = 0.0;
2878 bcttlen = 1;
2879 }
2880
2881 if (adxtail != 0.0) {
2882 temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, adxtail, temp16a);
2883 axtbctlen = scale_expansion_zeroelim(bctlen, bct, adxtail, axtbct);
2884 temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, 2.0 * adx,
2885 temp32a);
2886 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2887 temp32alen, temp32a, temp48);
2888 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2889 temp48, finother);
2890 finswap = finnow; finnow = finother; finother = finswap;
2891 if (bdytail != 0.0) {
2892 temp8len = scale_expansion_zeroelim(4, cc, adxtail, temp8);
2893 temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
2894 temp16a);
2895 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
2896 temp16a, finother);
2897 finswap = finnow; finnow = finother; finother = finswap;
2898 }
2899 if (cdytail != 0.0) {
2900 temp8len = scale_expansion_zeroelim(4, bb, -adxtail, temp8);
2901 temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
2902 temp16a);
2903 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
2904 temp16a, finother);
2905 finswap = finnow; finnow = finother; finother = finswap;
2906 }
2907
2908 temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, adxtail,
2909 temp32a);
2910 axtbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adxtail, axtbctt);
2911 temp16alen = scale_expansion_zeroelim(axtbcttlen, axtbctt, 2.0 * adx,
2912 temp16a);
2913 temp16blen = scale_expansion_zeroelim(axtbcttlen, axtbctt, adxtail,
2914 temp16b);
2915 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2916 temp16blen, temp16b, temp32b);
2917 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
2918 temp32blen, temp32b, temp64);
2919 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
2920 temp64, finother);
2921 finswap = finnow; finnow = finother; finother = finswap;
2922 }
2923 if (adytail != 0.0) {
2924 temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, adytail, temp16a);
2925 aytbctlen = scale_expansion_zeroelim(bctlen, bct, adytail, aytbct);
2926 temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, 2.0 * ady,
2927 temp32a);
2928 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2929 temp32alen, temp32a, temp48);
2930 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2931 temp48, finother);
2932 finswap = finnow; finnow = finother; finother = finswap;
2933
2934
2935 temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, adytail,
2936 temp32a);
2937 aytbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adytail, aytbctt);
2938 temp16alen = scale_expansion_zeroelim(aytbcttlen, aytbctt, 2.0 * ady,
2939 temp16a);
2940 temp16blen = scale_expansion_zeroelim(aytbcttlen, aytbctt, adytail,
2941 temp16b);
2942 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2943 temp16blen, temp16b, temp32b);
2944 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
2945 temp32blen, temp32b, temp64);
2946 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
2947 temp64, finother);
2948 finswap = finnow; finnow = finother; finother = finswap;
2949 }
2950 }
2951 if ((bdxtail != 0.0) || (bdytail != 0.0)) {
2952 if ((cdxtail != 0.0) || (cdytail != 0.0)
2953 || (adxtail != 0.0) || (adytail != 0.0)) {
2954 Two_Product(cdxtail, ady, ti1, ti0);
2955 Two_Product(cdx, adytail, tj1, tj0);
2956 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
2957 u[3] = u3;
2958 negate = -cdy;
2959 Two_Product(adxtail, negate, ti1, ti0);
2960 negate = -cdytail;
2961 Two_Product(adx, negate, tj1, tj0);
2962 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
2963 v[3] = v3;
2964 catlen = fast_expansion_sum_zeroelim(4, u, 4, v, cat);
2965
2966 Two_Product(cdxtail, adytail, ti1, ti0);
2967 Two_Product(adxtail, cdytail, tj1, tj0);
2968 Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]);
2969 catt[3] = catt3;
2970 cattlen = 4;
2971 } else {
2972 cat[0] = 0.0;
2973 catlen = 1;
2974 catt[0] = 0.0;
2975 cattlen = 1;
2976 }
2977
2978 if (bdxtail != 0.0) {
2979 temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, bdxtail, temp16a);
2980 bxtcatlen = scale_expansion_zeroelim(catlen, cat, bdxtail, bxtcat);
2981 temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, 2.0 * bdx,
2982 temp32a);
2983 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2984 temp32alen, temp32a, temp48);
2985 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2986 temp48, finother);
2987 finswap = finnow; finnow = finother; finother = finswap;
2988 if (cdytail != 0.0) {
2989 temp8len = scale_expansion_zeroelim(4, aa, bdxtail, temp8);
2990 temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
2991 temp16a);
2992 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
2993 temp16a, finother);
2994 finswap = finnow; finnow = finother; finother = finswap;
2995 }
2996 if (adytail != 0.0) {
2997 temp8len = scale_expansion_zeroelim(4, cc, -bdxtail, temp8);
2998 temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
2999 temp16a);
3000 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
3001 temp16a, finother);
3002 finswap = finnow; finnow = finother; finother = finswap;
3003 }
3004
3005 temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, bdxtail,
3006 temp32a);
3007 bxtcattlen = scale_expansion_zeroelim(cattlen, catt, bdxtail, bxtcatt);
3008 temp16alen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, 2.0 * bdx,
3009 temp16a);
3010 temp16blen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, bdxtail,
3011 temp16b);
3012 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3013 temp16blen, temp16b, temp32b);
3014 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3015 temp32blen, temp32b, temp64);
3016 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
3017 temp64, finother);
3018 finswap = finnow; finnow = finother; finother = finswap;
3019 }
3020 if (bdytail != 0.0) {
3021 temp16alen = scale_expansion_zeroelim(bytcalen, bytca, bdytail, temp16a);
3022 bytcatlen = scale_expansion_zeroelim(catlen, cat, bdytail, bytcat);
3023 temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, 2.0 * bdy,
3024 temp32a);
3025 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3026 temp32alen, temp32a, temp48);
3027 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
3028 temp48, finother);
3029 finswap = finnow; finnow = finother; finother = finswap;
3030
3031
3032 temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, bdytail,
3033 temp32a);
3034 bytcattlen = scale_expansion_zeroelim(cattlen, catt, bdytail, bytcatt);
3035 temp16alen = scale_expansion_zeroelim(bytcattlen, bytcatt, 2.0 * bdy,
3036 temp16a);
3037 temp16blen = scale_expansion_zeroelim(bytcattlen, bytcatt, bdytail,
3038 temp16b);
3039 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3040 temp16blen, temp16b, temp32b);
3041 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3042 temp32blen, temp32b, temp64);
3043 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
3044 temp64, finother);
3045 finswap = finnow; finnow = finother; finother = finswap;
3046 }
3047 }
3048 if ((cdxtail != 0.0) || (cdytail != 0.0)) {
3049 if ((adxtail != 0.0) || (adytail != 0.0)
3050 || (bdxtail != 0.0) || (bdytail != 0.0)) {
3051 Two_Product(adxtail, bdy, ti1, ti0);
3052 Two_Product(adx, bdytail, tj1, tj0);
3053 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
3054 u[3] = u3;
3055 negate = -ady;
3056 Two_Product(bdxtail, negate, ti1, ti0);
3057 negate = -adytail;
3058 Two_Product(bdx, negate, tj1, tj0);
3059 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
3060 v[3] = v3;
3061 abtlen = fast_expansion_sum_zeroelim(4, u, 4, v, abt);
3062
3063 Two_Product(adxtail, bdytail, ti1, ti0);
3064 Two_Product(bdxtail, adytail, tj1, tj0);
3065 Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]);
3066 abtt[3] = abtt3;
3067 abttlen = 4;
3068 } else {
3069 abt[0] = 0.0;
3070 abtlen = 1;
3071 abtt[0] = 0.0;
3072 abttlen = 1;
3073 }
3074
3075 if (cdxtail != 0.0) {
3076 temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, cdxtail, temp16a);
3077 cxtabtlen = scale_expansion_zeroelim(abtlen, abt, cdxtail, cxtabt);
3078 temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, 2.0 * cdx,
3079 temp32a);
3080 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3081 temp32alen, temp32a, temp48);
3082 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
3083 temp48, finother);
3084 finswap = finnow; finnow = finother; finother = finswap;
3085 if (adytail != 0.0) {
3086 temp8len = scale_expansion_zeroelim(4, bb, cdxtail, temp8);
3087 temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
3088 temp16a);
3089 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
3090 temp16a, finother);
3091 finswap = finnow; finnow = finother; finother = finswap;
3092 }
3093 if (bdytail != 0.0) {
3094 temp8len = scale_expansion_zeroelim(4, aa, -cdxtail, temp8);
3095 temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
3096 temp16a);
3097 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
3098 temp16a, finother);
3099 finswap = finnow; finnow = finother; finother = finswap;
3100 }
3101
3102 temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, cdxtail,
3103 temp32a);
3104 cxtabttlen = scale_expansion_zeroelim(abttlen, abtt, cdxtail, cxtabtt);
3105 temp16alen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, 2.0 * cdx,
3106 temp16a);
3107 temp16blen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, cdxtail,
3108 temp16b);
3109 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3110 temp16blen, temp16b, temp32b);
3111 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3112 temp32blen, temp32b, temp64);
3113 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
3114 temp64, finother);
3115 finswap = finnow; finnow = finother; finother = finswap;
3116 }
3117 if (cdytail != 0.0) {
3118 temp16alen = scale_expansion_zeroelim(cytablen, cytab, cdytail, temp16a);
3119 cytabtlen = scale_expansion_zeroelim(abtlen, abt, cdytail, cytabt);
3120 temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, 2.0 * cdy,
3121 temp32a);
3122 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3123 temp32alen, temp32a, temp48);
3124 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
3125 temp48, finother);
3126 finswap = finnow; finnow = finother; finother = finswap;
3127
3128
3129 temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, cdytail,
3130 temp32a);
3131 cytabttlen = scale_expansion_zeroelim(abttlen, abtt, cdytail, cytabtt);
3132 temp16alen = scale_expansion_zeroelim(cytabttlen, cytabtt, 2.0 * cdy,
3133 temp16a);
3134 temp16blen = scale_expansion_zeroelim(cytabttlen, cytabtt, cdytail,
3135 temp16b);
3136 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
3137 temp16blen, temp16b, temp32b);
3138 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3139 temp32blen, temp32b, temp64);
3140 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
3141 temp64, finother);
3142 finswap = finnow; finnow = finother; finother = finswap;
3143 }
3144 }
3145
3146 return finnow[finlength - 1];
3147 }
3148
incircle(REAL * pa,REAL * pb,REAL * pc,REAL * pd)3149 REAL incircle(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
3150 {
3151 REAL adx, bdx, cdx, ady, bdy, cdy;
3152 REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
3153 REAL alift, blift, clift;
3154 REAL det;
3155 REAL permanent, errbound;
3156
3157 adx = pa[0] - pd[0];
3158 bdx = pb[0] - pd[0];
3159 cdx = pc[0] - pd[0];
3160 ady = pa[1] - pd[1];
3161 bdy = pb[1] - pd[1];
3162 cdy = pc[1] - pd[1];
3163
3164 bdxcdy = bdx * cdy;
3165 cdxbdy = cdx * bdy;
3166 alift = adx * adx + ady * ady;
3167
3168 cdxady = cdx * ady;
3169 adxcdy = adx * cdy;
3170 blift = bdx * bdx + bdy * bdy;
3171
3172 adxbdy = adx * bdy;
3173 bdxady = bdx * ady;
3174 clift = cdx * cdx + cdy * cdy;
3175
3176 det = alift * (bdxcdy - cdxbdy)
3177 + blift * (cdxady - adxcdy)
3178 + clift * (adxbdy - bdxady);
3179
3180 permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift
3181 + (Absolute(cdxady) + Absolute(adxcdy)) * blift
3182 + (Absolute(adxbdy) + Absolute(bdxady)) * clift;
3183 errbound = iccerrboundA * permanent;
3184 if ((det > errbound) || (-det > errbound)) {
3185 return det;
3186 }
3187
3188 return incircleadapt(pa, pb, pc, pd, permanent);
3189 }
3190
3191 /*****************************************************************************/
3192 /* */
3193 /* inspherefast() Approximate 3D insphere test. Nonrobust. */
3194 /* insphereexact() Exact 3D insphere test. Robust. */
3195 /* insphereslow() Another exact 3D insphere test. Robust. */
3196 /* insphere() Adaptive exact 3D insphere test. Robust. */
3197 /* */
3198 /* Return a positive value if the point pe lies inside the */
3199 /* sphere passing through pa, pb, pc, and pd; a negative value */
3200 /* if it lies outside; and zero if the five points are */
3201 /* cospherical. The points pa, pb, pc, and pd must be ordered */
3202 /* so that they have a positive orientation (as defined by */
3203 /* orient3d()), or the sign of the result will be reversed. */
3204 /* */
3205 /* Only the first and last routine should be used; the middle two are for */
3206 /* timings. */
3207 /* */
3208 /* The last three use exact arithmetic to ensure a correct answer. The */
3209 /* result returned is the determinant of a matrix. In insphere() only, */
3210 /* this determinant is computed adaptively, in the sense that exact */
3211 /* arithmetic is used only to the degree it is needed to ensure that the */
3212 /* returned value has the correct sign. Hence, insphere() is usually quite */
3213 /* fast, but will run more slowly when the input points are cospherical or */
3214 /* nearly so. */
3215 /* */
3216 /*****************************************************************************/
3217
inspherefast(REAL * pa,REAL * pb,REAL * pc,REAL * pd,REAL * pe)3218 REAL inspherefast(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
3219 {
3220 REAL aex, bex, cex, dex;
3221 REAL aey, bey, cey, dey;
3222 REAL aez, bez, cez, dez;
3223 REAL alift, blift, clift, dlift;
3224 REAL ab, bc, cd, da, ac, bd;
3225 REAL abc, bcd, cda, dab;
3226
3227 aex = pa[0] - pe[0];
3228 bex = pb[0] - pe[0];
3229 cex = pc[0] - pe[0];
3230 dex = pd[0] - pe[0];
3231 aey = pa[1] - pe[1];
3232 bey = pb[1] - pe[1];
3233 cey = pc[1] - pe[1];
3234 dey = pd[1] - pe[1];
3235 aez = pa[2] - pe[2];
3236 bez = pb[2] - pe[2];
3237 cez = pc[2] - pe[2];
3238 dez = pd[2] - pe[2];
3239
3240 ab = aex * bey - bex * aey;
3241 bc = bex * cey - cex * bey;
3242 cd = cex * dey - dex * cey;
3243 da = dex * aey - aex * dey;
3244
3245 ac = aex * cey - cex * aey;
3246 bd = bex * dey - dex * bey;
3247
3248 abc = aez * bc - bez * ac + cez * ab;
3249 bcd = bez * cd - cez * bd + dez * bc;
3250 cda = cez * da + dez * ac + aez * cd;
3251 dab = dez * ab + aez * bd + bez * da;
3252
3253 alift = aex * aex + aey * aey + aez * aez;
3254 blift = bex * bex + bey * bey + bez * bez;
3255 clift = cex * cex + cey * cey + cez * cez;
3256 dlift = dex * dex + dey * dey + dez * dez;
3257
3258 return (dlift * abc - clift * dab) + (blift * cda - alift * bcd);
3259 }
3260
insphereexact(REAL * pa,REAL * pb,REAL * pc,REAL * pd,REAL * pe)3261 REAL insphereexact(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
3262 {
3263 INEXACT REAL axby1, bxcy1, cxdy1, dxey1, exay1;
3264 INEXACT REAL bxay1, cxby1, dxcy1, exdy1, axey1;
3265 INEXACT REAL axcy1, bxdy1, cxey1, dxay1, exby1;
3266 INEXACT REAL cxay1, dxby1, excy1, axdy1, bxey1;
3267 REAL axby0, bxcy0, cxdy0, dxey0, exay0;
3268 REAL bxay0, cxby0, dxcy0, exdy0, axey0;
3269 REAL axcy0, bxdy0, cxey0, dxay0, exby0;
3270 REAL cxay0, dxby0, excy0, axdy0, bxey0;
3271 REAL ab[4], bc[4], cd[4], de[4], ea[4];
3272 REAL ac[4], bd[4], ce[4], da[4], eb[4];
3273 REAL temp8a[8], temp8b[8], temp16[16];
3274 int temp8alen, temp8blen, temp16len;
3275 REAL abc[24], bcd[24], cde[24], dea[24], eab[24];
3276 REAL abd[24], bce[24], cda[24], deb[24], eac[24];
3277 int abclen, bcdlen, cdelen, dealen, eablen;
3278 int abdlen, bcelen, cdalen, deblen, eaclen;
3279 REAL temp48a[48], temp48b[48];
3280 int temp48alen, temp48blen;
3281 REAL abcd[96], bcde[96], cdea[96], deab[96], eabc[96];
3282 int abcdlen, bcdelen, cdealen, deablen, eabclen;
3283 REAL temp192[192];
3284 REAL det384x[384], det384y[384], det384z[384];
3285 int xlen, ylen, zlen;
3286 REAL detxy[768];
3287 int xylen;
3288 REAL adet[1152], bdet[1152], cdet[1152], ddet[1152], edet[1152];
3289 int alen, blen, clen, dlen, elen;
3290 REAL abdet[2304], cddet[2304], cdedet[3456];
3291 int ablen, cdlen;
3292 REAL deter[5760];
3293 int deterlen;
3294 int i;
3295
3296 INEXACT REAL bvirt;
3297 REAL avirt, bround, around;
3298 INEXACT REAL c;
3299 INEXACT REAL abig;
3300 REAL ahi, alo, bhi, blo;
3301 REAL err1, err2, err3;
3302 INEXACT REAL _i, _j;
3303 REAL _0;
3304
3305
3306 Two_Product(pa[0], pb[1], axby1, axby0);
3307 Two_Product(pb[0], pa[1], bxay1, bxay0);
3308 Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]);
3309
3310 Two_Product(pb[0], pc[1], bxcy1, bxcy0);
3311 Two_Product(pc[0], pb[1], cxby1, cxby0);
3312 Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]);
3313
3314 Two_Product(pc[0], pd[1], cxdy1, cxdy0);
3315 Two_Product(pd[0], pc[1], dxcy1, dxcy0);
3316 Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]);
3317
3318 Two_Product(pd[0], pe[1], dxey1, dxey0);
3319 Two_Product(pe[0], pd[1], exdy1, exdy0);
3320 Two_Two_Diff(dxey1, dxey0, exdy1, exdy0, de[3], de[2], de[1], de[0]);
3321
3322 Two_Product(pe[0], pa[1], exay1, exay0);
3323 Two_Product(pa[0], pe[1], axey1, axey0);
3324 Two_Two_Diff(exay1, exay0, axey1, axey0, ea[3], ea[2], ea[1], ea[0]);
3325
3326 Two_Product(pa[0], pc[1], axcy1, axcy0);
3327 Two_Product(pc[0], pa[1], cxay1, cxay0);
3328 Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]);
3329
3330 Two_Product(pb[0], pd[1], bxdy1, bxdy0);
3331 Two_Product(pd[0], pb[1], dxby1, dxby0);
3332 Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]);
3333
3334 Two_Product(pc[0], pe[1], cxey1, cxey0);
3335 Two_Product(pe[0], pc[1], excy1, excy0);
3336 Two_Two_Diff(cxey1, cxey0, excy1, excy0, ce[3], ce[2], ce[1], ce[0]);
3337
3338 Two_Product(pd[0], pa[1], dxay1, dxay0);
3339 Two_Product(pa[0], pd[1], axdy1, axdy0);
3340 Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]);
3341
3342 Two_Product(pe[0], pb[1], exby1, exby0);
3343 Two_Product(pb[0], pe[1], bxey1, bxey0);
3344 Two_Two_Diff(exby1, exby0, bxey1, bxey0, eb[3], eb[2], eb[1], eb[0]);
3345
3346 temp8alen = scale_expansion_zeroelim(4, bc, pa[2], temp8a);
3347 temp8blen = scale_expansion_zeroelim(4, ac, -pb[2], temp8b);
3348 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3349 temp16);
3350 temp8alen = scale_expansion_zeroelim(4, ab, pc[2], temp8a);
3351 abclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3352 abc);
3353
3354 temp8alen = scale_expansion_zeroelim(4, cd, pb[2], temp8a);
3355 temp8blen = scale_expansion_zeroelim(4, bd, -pc[2], temp8b);
3356 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3357 temp16);
3358 temp8alen = scale_expansion_zeroelim(4, bc, pd[2], temp8a);
3359 bcdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3360 bcd);
3361
3362 temp8alen = scale_expansion_zeroelim(4, de, pc[2], temp8a);
3363 temp8blen = scale_expansion_zeroelim(4, ce, -pd[2], temp8b);
3364 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3365 temp16);
3366 temp8alen = scale_expansion_zeroelim(4, cd, pe[2], temp8a);
3367 cdelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3368 cde);
3369
3370 temp8alen = scale_expansion_zeroelim(4, ea, pd[2], temp8a);
3371 temp8blen = scale_expansion_zeroelim(4, da, -pe[2], temp8b);
3372 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3373 temp16);
3374 temp8alen = scale_expansion_zeroelim(4, de, pa[2], temp8a);
3375 dealen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3376 dea);
3377
3378 temp8alen = scale_expansion_zeroelim(4, ab, pe[2], temp8a);
3379 temp8blen = scale_expansion_zeroelim(4, eb, -pa[2], temp8b);
3380 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3381 temp16);
3382 temp8alen = scale_expansion_zeroelim(4, ea, pb[2], temp8a);
3383 eablen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3384 eab);
3385
3386 temp8alen = scale_expansion_zeroelim(4, bd, pa[2], temp8a);
3387 temp8blen = scale_expansion_zeroelim(4, da, pb[2], temp8b);
3388 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3389 temp16);
3390 temp8alen = scale_expansion_zeroelim(4, ab, pd[2], temp8a);
3391 abdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3392 abd);
3393
3394 temp8alen = scale_expansion_zeroelim(4, ce, pb[2], temp8a);
3395 temp8blen = scale_expansion_zeroelim(4, eb, pc[2], temp8b);
3396 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3397 temp16);
3398 temp8alen = scale_expansion_zeroelim(4, bc, pe[2], temp8a);
3399 bcelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3400 bce);
3401
3402 temp8alen = scale_expansion_zeroelim(4, da, pc[2], temp8a);
3403 temp8blen = scale_expansion_zeroelim(4, ac, pd[2], temp8b);
3404 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3405 temp16);
3406 temp8alen = scale_expansion_zeroelim(4, cd, pa[2], temp8a);
3407 cdalen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3408 cda);
3409
3410 temp8alen = scale_expansion_zeroelim(4, eb, pd[2], temp8a);
3411 temp8blen = scale_expansion_zeroelim(4, bd, pe[2], temp8b);
3412 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3413 temp16);
3414 temp8alen = scale_expansion_zeroelim(4, de, pb[2], temp8a);
3415 deblen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3416 deb);
3417
3418 temp8alen = scale_expansion_zeroelim(4, ac, pe[2], temp8a);
3419 temp8blen = scale_expansion_zeroelim(4, ce, pa[2], temp8b);
3420 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
3421 temp16);
3422 temp8alen = scale_expansion_zeroelim(4, ea, pc[2], temp8a);
3423 eaclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
3424 eac);
3425
3426 temp48alen = fast_expansion_sum_zeroelim(cdelen, cde, bcelen, bce, temp48a);
3427 temp48blen = fast_expansion_sum_zeroelim(deblen, deb, bcdlen, bcd, temp48b);
3428 for (i = 0; i < temp48blen; i++) {
3429 temp48b[i] = -temp48b[i];
3430 }
3431 bcdelen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
3432 temp48blen, temp48b, bcde);
3433 xlen = scale_expansion_zeroelim(bcdelen, bcde, pa[0], temp192);
3434 xlen = scale_expansion_zeroelim(xlen, temp192, pa[0], det384x);
3435 ylen = scale_expansion_zeroelim(bcdelen, bcde, pa[1], temp192);
3436 ylen = scale_expansion_zeroelim(ylen, temp192, pa[1], det384y);
3437 zlen = scale_expansion_zeroelim(bcdelen, bcde, pa[2], temp192);
3438 zlen = scale_expansion_zeroelim(zlen, temp192, pa[2], det384z);
3439 xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
3440 alen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, adet);
3441
3442 temp48alen = fast_expansion_sum_zeroelim(dealen, dea, cdalen, cda, temp48a);
3443 temp48blen = fast_expansion_sum_zeroelim(eaclen, eac, cdelen, cde, temp48b);
3444 for (i = 0; i < temp48blen; i++) {
3445 temp48b[i] = -temp48b[i];
3446 }
3447 cdealen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
3448 temp48blen, temp48b, cdea);
3449 xlen = scale_expansion_zeroelim(cdealen, cdea, pb[0], temp192);
3450 xlen = scale_expansion_zeroelim(xlen, temp192, pb[0], det384x);
3451 ylen = scale_expansion_zeroelim(cdealen, cdea, pb[1], temp192);
3452 ylen = scale_expansion_zeroelim(ylen, temp192, pb[1], det384y);
3453 zlen = scale_expansion_zeroelim(cdealen, cdea, pb[2], temp192);
3454 zlen = scale_expansion_zeroelim(zlen, temp192, pb[2], det384z);
3455 xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
3456 blen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, bdet);
3457
3458 temp48alen = fast_expansion_sum_zeroelim(eablen, eab, deblen, deb, temp48a);
3459 temp48blen = fast_expansion_sum_zeroelim(abdlen, abd, dealen, dea, temp48b);
3460 for (i = 0; i < temp48blen; i++) {
3461 temp48b[i] = -temp48b[i];
3462 }
3463 deablen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
3464 temp48blen, temp48b, deab);
3465 xlen = scale_expansion_zeroelim(deablen, deab, pc[0], temp192);
3466 xlen = scale_expansion_zeroelim(xlen, temp192, pc[0], det384x);
3467 ylen = scale_expansion_zeroelim(deablen, deab, pc[1], temp192);
3468 ylen = scale_expansion_zeroelim(ylen, temp192, pc[1], det384y);
3469 zlen = scale_expansion_zeroelim(deablen, deab, pc[2], temp192);
3470 zlen = scale_expansion_zeroelim(zlen, temp192, pc[2], det384z);
3471 xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
3472 clen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, cdet);
3473
3474 temp48alen = fast_expansion_sum_zeroelim(abclen, abc, eaclen, eac, temp48a);
3475 temp48blen = fast_expansion_sum_zeroelim(bcelen, bce, eablen, eab, temp48b);
3476 for (i = 0; i < temp48blen; i++) {
3477 temp48b[i] = -temp48b[i];
3478 }
3479 eabclen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
3480 temp48blen, temp48b, eabc);
3481 xlen = scale_expansion_zeroelim(eabclen, eabc, pd[0], temp192);
3482 xlen = scale_expansion_zeroelim(xlen, temp192, pd[0], det384x);
3483 ylen = scale_expansion_zeroelim(eabclen, eabc, pd[1], temp192);
3484 ylen = scale_expansion_zeroelim(ylen, temp192, pd[1], det384y);
3485 zlen = scale_expansion_zeroelim(eabclen, eabc, pd[2], temp192);
3486 zlen = scale_expansion_zeroelim(zlen, temp192, pd[2], det384z);
3487 xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
3488 dlen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, ddet);
3489
3490 temp48alen = fast_expansion_sum_zeroelim(bcdlen, bcd, abdlen, abd, temp48a);
3491 temp48blen = fast_expansion_sum_zeroelim(cdalen, cda, abclen, abc, temp48b);
3492 for (i = 0; i < temp48blen; i++) {
3493 temp48b[i] = -temp48b[i];
3494 }
3495 abcdlen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
3496 temp48blen, temp48b, abcd);
3497 xlen = scale_expansion_zeroelim(abcdlen, abcd, pe[0], temp192);
3498 xlen = scale_expansion_zeroelim(xlen, temp192, pe[0], det384x);
3499 ylen = scale_expansion_zeroelim(abcdlen, abcd, pe[1], temp192);
3500 ylen = scale_expansion_zeroelim(ylen, temp192, pe[1], det384y);
3501 zlen = scale_expansion_zeroelim(abcdlen, abcd, pe[2], temp192);
3502 zlen = scale_expansion_zeroelim(zlen, temp192, pe[2], det384z);
3503 xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
3504 elen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, edet);
3505
3506 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
3507 cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
3508 cdelen = fast_expansion_sum_zeroelim(cdlen, cddet, elen, edet, cdedet);
3509 deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdelen, cdedet, deter);
3510
3511 return deter[deterlen - 1];
3512 }
3513
insphereslow(REAL * pa,REAL * pb,REAL * pc,REAL * pd,REAL * pe)3514 REAL insphereslow(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
3515 {
3516 INEXACT REAL aex, bex, cex, dex, aey, bey, cey, dey, aez, bez, cez, dez;
3517 REAL aextail, bextail, cextail, dextail;
3518 REAL aeytail, beytail, ceytail, deytail;
3519 REAL aeztail, beztail, ceztail, deztail;
3520 REAL negate, negatetail;
3521 INEXACT REAL axby7, bxcy7, cxdy7, dxay7, axcy7, bxdy7;
3522 INEXACT REAL bxay7, cxby7, dxcy7, axdy7, cxay7, dxby7;
3523 REAL axby[8], bxcy[8], cxdy[8], dxay[8], axcy[8], bxdy[8];
3524 REAL bxay[8], cxby[8], dxcy[8], axdy[8], cxay[8], dxby[8];
3525 REAL ab[16], bc[16], cd[16], da[16], ac[16], bd[16];
3526 int ablen, bclen, cdlen, dalen, aclen, bdlen;
3527 REAL temp32a[32], temp32b[32], temp64a[64], temp64b[64], temp64c[64];
3528 int temp32alen, temp32blen, temp64alen, temp64blen, temp64clen;
3529 REAL temp128[128], temp192[192];
3530 int temp128len, temp192len;
3531 REAL detx[384], detxx[768], detxt[384], detxxt[768], detxtxt[768];
3532 int xlen, xxlen, xtlen, xxtlen, xtxtlen;
3533 REAL x1[1536], x2[2304];
3534 int x1len, x2len;
3535 REAL dety[384], detyy[768], detyt[384], detyyt[768], detytyt[768];
3536 int ylen, yylen, ytlen, yytlen, ytytlen;
3537 REAL y1[1536], y2[2304];
3538 int y1len, y2len;
3539 REAL detz[384], detzz[768], detzt[384], detzzt[768], detztzt[768];
3540 int zlen, zzlen, ztlen, zztlen, ztztlen;
3541 REAL z1[1536], z2[2304];
3542 int z1len, z2len;
3543 REAL detxy[4608];
3544 int xylen;
3545 REAL adet[6912], bdet[6912], cdet[6912], ddet[6912];
3546 int alen, blen, clen, dlen;
3547 REAL abdet[13824], cddet[13824], deter[27648];
3548 int deterlen;
3549 int i;
3550
3551 INEXACT REAL bvirt;
3552 REAL avirt, bround, around;
3553 INEXACT REAL c;
3554 INEXACT REAL abig;
3555 REAL a0hi, a0lo, a1hi, a1lo, bhi, blo;
3556 REAL err1, err2, err3;
3557 INEXACT REAL _i, _j, _k, _l, _m, _n;
3558 REAL _0, _1, _2;
3559
3560 Two_Diff(pa[0], pe[0], aex, aextail);
3561 Two_Diff(pa[1], pe[1], aey, aeytail);
3562 Two_Diff(pa[2], pe[2], aez, aeztail);
3563 Two_Diff(pb[0], pe[0], bex, bextail);
3564 Two_Diff(pb[1], pe[1], bey, beytail);
3565 Two_Diff(pb[2], pe[2], bez, beztail);
3566 Two_Diff(pc[0], pe[0], cex, cextail);
3567 Two_Diff(pc[1], pe[1], cey, ceytail);
3568 Two_Diff(pc[2], pe[2], cez, ceztail);
3569 Two_Diff(pd[0], pe[0], dex, dextail);
3570 Two_Diff(pd[1], pe[1], dey, deytail);
3571 Two_Diff(pd[2], pe[2], dez, deztail);
3572
3573 Two_Two_Product(aex, aextail, bey, beytail,
3574 axby7, axby[6], axby[5], axby[4],
3575 axby[3], axby[2], axby[1], axby[0]);
3576 axby[7] = axby7;
3577 negate = -aey;
3578 negatetail = -aeytail;
3579 Two_Two_Product(bex, bextail, negate, negatetail,
3580 bxay7, bxay[6], bxay[5], bxay[4],
3581 bxay[3], bxay[2], bxay[1], bxay[0]);
3582 bxay[7] = bxay7;
3583 ablen = fast_expansion_sum_zeroelim(8, axby, 8, bxay, ab);
3584 Two_Two_Product(bex, bextail, cey, ceytail,
3585 bxcy7, bxcy[6], bxcy[5], bxcy[4],
3586 bxcy[3], bxcy[2], bxcy[1], bxcy[0]);
3587 bxcy[7] = bxcy7;
3588 negate = -bey;
3589 negatetail = -beytail;
3590 Two_Two_Product(cex, cextail, negate, negatetail,
3591 cxby7, cxby[6], cxby[5], cxby[4],
3592 cxby[3], cxby[2], cxby[1], cxby[0]);
3593 cxby[7] = cxby7;
3594 bclen = fast_expansion_sum_zeroelim(8, bxcy, 8, cxby, bc);
3595 Two_Two_Product(cex, cextail, dey, deytail,
3596 cxdy7, cxdy[6], cxdy[5], cxdy[4],
3597 cxdy[3], cxdy[2], cxdy[1], cxdy[0]);
3598 cxdy[7] = cxdy7;
3599 negate = -cey;
3600 negatetail = -ceytail;
3601 Two_Two_Product(dex, dextail, negate, negatetail,
3602 dxcy7, dxcy[6], dxcy[5], dxcy[4],
3603 dxcy[3], dxcy[2], dxcy[1], dxcy[0]);
3604 dxcy[7] = dxcy7;
3605 cdlen = fast_expansion_sum_zeroelim(8, cxdy, 8, dxcy, cd);
3606 Two_Two_Product(dex, dextail, aey, aeytail,
3607 dxay7, dxay[6], dxay[5], dxay[4],
3608 dxay[3], dxay[2], dxay[1], dxay[0]);
3609 dxay[7] = dxay7;
3610 negate = -dey;
3611 negatetail = -deytail;
3612 Two_Two_Product(aex, aextail, negate, negatetail,
3613 axdy7, axdy[6], axdy[5], axdy[4],
3614 axdy[3], axdy[2], axdy[1], axdy[0]);
3615 axdy[7] = axdy7;
3616 dalen = fast_expansion_sum_zeroelim(8, dxay, 8, axdy, da);
3617 Two_Two_Product(aex, aextail, cey, ceytail,
3618 axcy7, axcy[6], axcy[5], axcy[4],
3619 axcy[3], axcy[2], axcy[1], axcy[0]);
3620 axcy[7] = axcy7;
3621 negate = -aey;
3622 negatetail = -aeytail;
3623 Two_Two_Product(cex, cextail, negate, negatetail,
3624 cxay7, cxay[6], cxay[5], cxay[4],
3625 cxay[3], cxay[2], cxay[1], cxay[0]);
3626 cxay[7] = cxay7;
3627 aclen = fast_expansion_sum_zeroelim(8, axcy, 8, cxay, ac);
3628 Two_Two_Product(bex, bextail, dey, deytail,
3629 bxdy7, bxdy[6], bxdy[5], bxdy[4],
3630 bxdy[3], bxdy[2], bxdy[1], bxdy[0]);
3631 bxdy[7] = bxdy7;
3632 negate = -bey;
3633 negatetail = -beytail;
3634 Two_Two_Product(dex, dextail, negate, negatetail,
3635 dxby7, dxby[6], dxby[5], dxby[4],
3636 dxby[3], dxby[2], dxby[1], dxby[0]);
3637 dxby[7] = dxby7;
3638 bdlen = fast_expansion_sum_zeroelim(8, bxdy, 8, dxby, bd);
3639
3640 temp32alen = scale_expansion_zeroelim(cdlen, cd, -bez, temp32a);
3641 temp32blen = scale_expansion_zeroelim(cdlen, cd, -beztail, temp32b);
3642 temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3643 temp32blen, temp32b, temp64a);
3644 temp32alen = scale_expansion_zeroelim(bdlen, bd, cez, temp32a);
3645 temp32blen = scale_expansion_zeroelim(bdlen, bd, ceztail, temp32b);
3646 temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3647 temp32blen, temp32b, temp64b);
3648 temp32alen = scale_expansion_zeroelim(bclen, bc, -dez, temp32a);
3649 temp32blen = scale_expansion_zeroelim(bclen, bc, -deztail, temp32b);
3650 temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3651 temp32blen, temp32b, temp64c);
3652 temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a,
3653 temp64blen, temp64b, temp128);
3654 temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c,
3655 temp128len, temp128, temp192);
3656 xlen = scale_expansion_zeroelim(temp192len, temp192, aex, detx);
3657 xxlen = scale_expansion_zeroelim(xlen, detx, aex, detxx);
3658 xtlen = scale_expansion_zeroelim(temp192len, temp192, aextail, detxt);
3659 xxtlen = scale_expansion_zeroelim(xtlen, detxt, aex, detxxt);
3660 for (i = 0; i < xxtlen; i++) {
3661 detxxt[i] *= 2.0;
3662 }
3663 xtxtlen = scale_expansion_zeroelim(xtlen, detxt, aextail, detxtxt);
3664 x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
3665 x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
3666 ylen = scale_expansion_zeroelim(temp192len, temp192, aey, dety);
3667 yylen = scale_expansion_zeroelim(ylen, dety, aey, detyy);
3668 ytlen = scale_expansion_zeroelim(temp192len, temp192, aeytail, detyt);
3669 yytlen = scale_expansion_zeroelim(ytlen, detyt, aey, detyyt);
3670 for (i = 0; i < yytlen; i++) {
3671 detyyt[i] *= 2.0;
3672 }
3673 ytytlen = scale_expansion_zeroelim(ytlen, detyt, aeytail, detytyt);
3674 y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
3675 y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
3676 zlen = scale_expansion_zeroelim(temp192len, temp192, aez, detz);
3677 zzlen = scale_expansion_zeroelim(zlen, detz, aez, detzz);
3678 ztlen = scale_expansion_zeroelim(temp192len, temp192, aeztail, detzt);
3679 zztlen = scale_expansion_zeroelim(ztlen, detzt, aez, detzzt);
3680 for (i = 0; i < zztlen; i++) {
3681 detzzt[i] *= 2.0;
3682 }
3683 ztztlen = scale_expansion_zeroelim(ztlen, detzt, aeztail, detztzt);
3684 z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1);
3685 z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2);
3686 xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy);
3687 alen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, adet);
3688
3689 temp32alen = scale_expansion_zeroelim(dalen, da, cez, temp32a);
3690 temp32blen = scale_expansion_zeroelim(dalen, da, ceztail, temp32b);
3691 temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3692 temp32blen, temp32b, temp64a);
3693 temp32alen = scale_expansion_zeroelim(aclen, ac, dez, temp32a);
3694 temp32blen = scale_expansion_zeroelim(aclen, ac, deztail, temp32b);
3695 temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3696 temp32blen, temp32b, temp64b);
3697 temp32alen = scale_expansion_zeroelim(cdlen, cd, aez, temp32a);
3698 temp32blen = scale_expansion_zeroelim(cdlen, cd, aeztail, temp32b);
3699 temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3700 temp32blen, temp32b, temp64c);
3701 temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a,
3702 temp64blen, temp64b, temp128);
3703 temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c,
3704 temp128len, temp128, temp192);
3705 xlen = scale_expansion_zeroelim(temp192len, temp192, bex, detx);
3706 xxlen = scale_expansion_zeroelim(xlen, detx, bex, detxx);
3707 xtlen = scale_expansion_zeroelim(temp192len, temp192, bextail, detxt);
3708 xxtlen = scale_expansion_zeroelim(xtlen, detxt, bex, detxxt);
3709 for (i = 0; i < xxtlen; i++) {
3710 detxxt[i] *= 2.0;
3711 }
3712 xtxtlen = scale_expansion_zeroelim(xtlen, detxt, bextail, detxtxt);
3713 x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
3714 x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
3715 ylen = scale_expansion_zeroelim(temp192len, temp192, bey, dety);
3716 yylen = scale_expansion_zeroelim(ylen, dety, bey, detyy);
3717 ytlen = scale_expansion_zeroelim(temp192len, temp192, beytail, detyt);
3718 yytlen = scale_expansion_zeroelim(ytlen, detyt, bey, detyyt);
3719 for (i = 0; i < yytlen; i++) {
3720 detyyt[i] *= 2.0;
3721 }
3722 ytytlen = scale_expansion_zeroelim(ytlen, detyt, beytail, detytyt);
3723 y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
3724 y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
3725 zlen = scale_expansion_zeroelim(temp192len, temp192, bez, detz);
3726 zzlen = scale_expansion_zeroelim(zlen, detz, bez, detzz);
3727 ztlen = scale_expansion_zeroelim(temp192len, temp192, beztail, detzt);
3728 zztlen = scale_expansion_zeroelim(ztlen, detzt, bez, detzzt);
3729 for (i = 0; i < zztlen; i++) {
3730 detzzt[i] *= 2.0;
3731 }
3732 ztztlen = scale_expansion_zeroelim(ztlen, detzt, beztail, detztzt);
3733 z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1);
3734 z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2);
3735 xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy);
3736 blen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, bdet);
3737
3738 temp32alen = scale_expansion_zeroelim(ablen, ab, -dez, temp32a);
3739 temp32blen = scale_expansion_zeroelim(ablen, ab, -deztail, temp32b);
3740 temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3741 temp32blen, temp32b, temp64a);
3742 temp32alen = scale_expansion_zeroelim(bdlen, bd, -aez, temp32a);
3743 temp32blen = scale_expansion_zeroelim(bdlen, bd, -aeztail, temp32b);
3744 temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3745 temp32blen, temp32b, temp64b);
3746 temp32alen = scale_expansion_zeroelim(dalen, da, -bez, temp32a);
3747 temp32blen = scale_expansion_zeroelim(dalen, da, -beztail, temp32b);
3748 temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3749 temp32blen, temp32b, temp64c);
3750 temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a,
3751 temp64blen, temp64b, temp128);
3752 temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c,
3753 temp128len, temp128, temp192);
3754 xlen = scale_expansion_zeroelim(temp192len, temp192, cex, detx);
3755 xxlen = scale_expansion_zeroelim(xlen, detx, cex, detxx);
3756 xtlen = scale_expansion_zeroelim(temp192len, temp192, cextail, detxt);
3757 xxtlen = scale_expansion_zeroelim(xtlen, detxt, cex, detxxt);
3758 for (i = 0; i < xxtlen; i++) {
3759 detxxt[i] *= 2.0;
3760 }
3761 xtxtlen = scale_expansion_zeroelim(xtlen, detxt, cextail, detxtxt);
3762 x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
3763 x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
3764 ylen = scale_expansion_zeroelim(temp192len, temp192, cey, dety);
3765 yylen = scale_expansion_zeroelim(ylen, dety, cey, detyy);
3766 ytlen = scale_expansion_zeroelim(temp192len, temp192, ceytail, detyt);
3767 yytlen = scale_expansion_zeroelim(ytlen, detyt, cey, detyyt);
3768 for (i = 0; i < yytlen; i++) {
3769 detyyt[i] *= 2.0;
3770 }
3771 ytytlen = scale_expansion_zeroelim(ytlen, detyt, ceytail, detytyt);
3772 y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
3773 y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
3774 zlen = scale_expansion_zeroelim(temp192len, temp192, cez, detz);
3775 zzlen = scale_expansion_zeroelim(zlen, detz, cez, detzz);
3776 ztlen = scale_expansion_zeroelim(temp192len, temp192, ceztail, detzt);
3777 zztlen = scale_expansion_zeroelim(ztlen, detzt, cez, detzzt);
3778 for (i = 0; i < zztlen; i++) {
3779 detzzt[i] *= 2.0;
3780 }
3781 ztztlen = scale_expansion_zeroelim(ztlen, detzt, ceztail, detztzt);
3782 z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1);
3783 z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2);
3784 xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy);
3785 clen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, cdet);
3786
3787 temp32alen = scale_expansion_zeroelim(bclen, bc, aez, temp32a);
3788 temp32blen = scale_expansion_zeroelim(bclen, bc, aeztail, temp32b);
3789 temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3790 temp32blen, temp32b, temp64a);
3791 temp32alen = scale_expansion_zeroelim(aclen, ac, -bez, temp32a);
3792 temp32blen = scale_expansion_zeroelim(aclen, ac, -beztail, temp32b);
3793 temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3794 temp32blen, temp32b, temp64b);
3795 temp32alen = scale_expansion_zeroelim(ablen, ab, cez, temp32a);
3796 temp32blen = scale_expansion_zeroelim(ablen, ab, ceztail, temp32b);
3797 temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
3798 temp32blen, temp32b, temp64c);
3799 temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a,
3800 temp64blen, temp64b, temp128);
3801 temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c,
3802 temp128len, temp128, temp192);
3803 xlen = scale_expansion_zeroelim(temp192len, temp192, dex, detx);
3804 xxlen = scale_expansion_zeroelim(xlen, detx, dex, detxx);
3805 xtlen = scale_expansion_zeroelim(temp192len, temp192, dextail, detxt);
3806 xxtlen = scale_expansion_zeroelim(xtlen, detxt, dex, detxxt);
3807 for (i = 0; i < xxtlen; i++) {
3808 detxxt[i] *= 2.0;
3809 }
3810 xtxtlen = scale_expansion_zeroelim(xtlen, detxt, dextail, detxtxt);
3811 x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
3812 x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
3813 ylen = scale_expansion_zeroelim(temp192len, temp192, dey, dety);
3814 yylen = scale_expansion_zeroelim(ylen, dety, dey, detyy);
3815 ytlen = scale_expansion_zeroelim(temp192len, temp192, deytail, detyt);
3816 yytlen = scale_expansion_zeroelim(ytlen, detyt, dey, detyyt);
3817 for (i = 0; i < yytlen; i++) {
3818 detyyt[i] *= 2.0;
3819 }
3820 ytytlen = scale_expansion_zeroelim(ytlen, detyt, deytail, detytyt);
3821 y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
3822 y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
3823 zlen = scale_expansion_zeroelim(temp192len, temp192, dez, detz);
3824 zzlen = scale_expansion_zeroelim(zlen, detz, dez, detzz);
3825 ztlen = scale_expansion_zeroelim(temp192len, temp192, deztail, detzt);
3826 zztlen = scale_expansion_zeroelim(ztlen, detzt, dez, detzzt);
3827 for (i = 0; i < zztlen; i++) {
3828 detzzt[i] *= 2.0;
3829 }
3830 ztztlen = scale_expansion_zeroelim(ztlen, detzt, deztail, detztzt);
3831 z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1);
3832 z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2);
3833 xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy);
3834 dlen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, ddet);
3835
3836 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
3837 cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
3838 deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, deter);
3839
3840 return deter[deterlen - 1];
3841 }
3842
insphereadapt(REAL * pa,REAL * pb,REAL * pc,REAL * pd,REAL * pe,REAL permanent)3843 REAL insphereadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe,
3844 REAL permanent)
3845 {
3846 INEXACT REAL aex, bex, cex, dex, aey, bey, cey, dey, aez, bez, cez, dez;
3847 REAL det, errbound;
3848
3849 INEXACT REAL aexbey1, bexaey1, bexcey1, cexbey1;
3850 INEXACT REAL cexdey1, dexcey1, dexaey1, aexdey1;
3851 INEXACT REAL aexcey1, cexaey1, bexdey1, dexbey1;
3852 REAL aexbey0, bexaey0, bexcey0, cexbey0;
3853 REAL cexdey0, dexcey0, dexaey0, aexdey0;
3854 REAL aexcey0, cexaey0, bexdey0, dexbey0;
3855 REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4];
3856 INEXACT REAL ab3, bc3, cd3, da3, ac3, bd3;
3857 REAL abeps, bceps, cdeps, daeps, aceps, bdeps;
3858 REAL temp8a[8], temp8b[8], temp8c[8], temp16[16], temp24[24], temp48[48];
3859 int temp8alen, temp8blen, temp8clen, temp16len, temp24len, temp48len;
3860 REAL xdet[96], ydet[96], zdet[96], xydet[192];
3861 int xlen, ylen, zlen, xylen;
3862 REAL adet[288], bdet[288], cdet[288], ddet[288];
3863 int alen, blen, clen, dlen;
3864 REAL abdet[576], cddet[576];
3865 int ablen, cdlen;
3866 REAL fin1[1152];
3867 int finlength;
3868
3869 REAL aextail, bextail, cextail, dextail;
3870 REAL aeytail, beytail, ceytail, deytail;
3871 REAL aeztail, beztail, ceztail, deztail;
3872
3873 INEXACT REAL bvirt;
3874 REAL avirt, bround, around;
3875 INEXACT REAL c;
3876 INEXACT REAL abig;
3877 REAL ahi, alo, bhi, blo;
3878 REAL err1, err2, err3;
3879 INEXACT REAL _i, _j;
3880 REAL _0;
3881
3882
3883 aex = (REAL) (pa[0] - pe[0]);
3884 bex = (REAL) (pb[0] - pe[0]);
3885 cex = (REAL) (pc[0] - pe[0]);
3886 dex = (REAL) (pd[0] - pe[0]);
3887 aey = (REAL) (pa[1] - pe[1]);
3888 bey = (REAL) (pb[1] - pe[1]);
3889 cey = (REAL) (pc[1] - pe[1]);
3890 dey = (REAL) (pd[1] - pe[1]);
3891 aez = (REAL) (pa[2] - pe[2]);
3892 bez = (REAL) (pb[2] - pe[2]);
3893 cez = (REAL) (pc[2] - pe[2]);
3894 dez = (REAL) (pd[2] - pe[2]);
3895
3896 Two_Product(aex, bey, aexbey1, aexbey0);
3897 Two_Product(bex, aey, bexaey1, bexaey0);
3898 Two_Two_Diff(aexbey1, aexbey0, bexaey1, bexaey0, ab3, ab[2], ab[1], ab[0]);
3899 ab[3] = ab3;
3900
3901 Two_Product(bex, cey, bexcey1, bexcey0);
3902 Two_Product(cex, bey, cexbey1, cexbey0);
3903 Two_Two_Diff(bexcey1, bexcey0, cexbey1, cexbey0, bc3, bc[2], bc[1], bc[0]);
3904 bc[3] = bc3;
3905
3906 Two_Product(cex, dey, cexdey1, cexdey0);
3907 Two_Product(dex, cey, dexcey1, dexcey0);
3908 Two_Two_Diff(cexdey1, cexdey0, dexcey1, dexcey0, cd3, cd[2], cd[1], cd[0]);
3909 cd[3] = cd3;
3910
3911 Two_Product(dex, aey, dexaey1, dexaey0);
3912 Two_Product(aex, dey, aexdey1, aexdey0);
3913 Two_Two_Diff(dexaey1, dexaey0, aexdey1, aexdey0, da3, da[2], da[1], da[0]);
3914 da[3] = da3;
3915
3916 Two_Product(aex, cey, aexcey1, aexcey0);
3917 Two_Product(cex, aey, cexaey1, cexaey0);
3918 Two_Two_Diff(aexcey1, aexcey0, cexaey1, cexaey0, ac3, ac[2], ac[1], ac[0]);
3919 ac[3] = ac3;
3920
3921 Two_Product(bex, dey, bexdey1, bexdey0);
3922 Two_Product(dex, bey, dexbey1, dexbey0);
3923 Two_Two_Diff(bexdey1, bexdey0, dexbey1, dexbey0, bd3, bd[2], bd[1], bd[0]);
3924 bd[3] = bd3;
3925
3926 temp8alen = scale_expansion_zeroelim(4, cd, bez, temp8a);
3927 temp8blen = scale_expansion_zeroelim(4, bd, -cez, temp8b);
3928 temp8clen = scale_expansion_zeroelim(4, bc, dez, temp8c);
3929 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
3930 temp8blen, temp8b, temp16);
3931 temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
3932 temp16len, temp16, temp24);
3933 temp48len = scale_expansion_zeroelim(temp24len, temp24, aex, temp48);
3934 xlen = scale_expansion_zeroelim(temp48len, temp48, -aex, xdet);
3935 temp48len = scale_expansion_zeroelim(temp24len, temp24, aey, temp48);
3936 ylen = scale_expansion_zeroelim(temp48len, temp48, -aey, ydet);
3937 temp48len = scale_expansion_zeroelim(temp24len, temp24, aez, temp48);
3938 zlen = scale_expansion_zeroelim(temp48len, temp48, -aez, zdet);
3939 xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
3940 alen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, adet);
3941
3942 temp8alen = scale_expansion_zeroelim(4, da, cez, temp8a);
3943 temp8blen = scale_expansion_zeroelim(4, ac, dez, temp8b);
3944 temp8clen = scale_expansion_zeroelim(4, cd, aez, temp8c);
3945 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
3946 temp8blen, temp8b, temp16);
3947 temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
3948 temp16len, temp16, temp24);
3949 temp48len = scale_expansion_zeroelim(temp24len, temp24, bex, temp48);
3950 xlen = scale_expansion_zeroelim(temp48len, temp48, bex, xdet);
3951 temp48len = scale_expansion_zeroelim(temp24len, temp24, bey, temp48);
3952 ylen = scale_expansion_zeroelim(temp48len, temp48, bey, ydet);
3953 temp48len = scale_expansion_zeroelim(temp24len, temp24, bez, temp48);
3954 zlen = scale_expansion_zeroelim(temp48len, temp48, bez, zdet);
3955 xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
3956 blen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, bdet);
3957
3958 temp8alen = scale_expansion_zeroelim(4, ab, dez, temp8a);
3959 temp8blen = scale_expansion_zeroelim(4, bd, aez, temp8b);
3960 temp8clen = scale_expansion_zeroelim(4, da, bez, temp8c);
3961 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
3962 temp8blen, temp8b, temp16);
3963 temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
3964 temp16len, temp16, temp24);
3965 temp48len = scale_expansion_zeroelim(temp24len, temp24, cex, temp48);
3966 xlen = scale_expansion_zeroelim(temp48len, temp48, -cex, xdet);
3967 temp48len = scale_expansion_zeroelim(temp24len, temp24, cey, temp48);
3968 ylen = scale_expansion_zeroelim(temp48len, temp48, -cey, ydet);
3969 temp48len = scale_expansion_zeroelim(temp24len, temp24, cez, temp48);
3970 zlen = scale_expansion_zeroelim(temp48len, temp48, -cez, zdet);
3971 xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
3972 clen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, cdet);
3973
3974 temp8alen = scale_expansion_zeroelim(4, bc, aez, temp8a);
3975 temp8blen = scale_expansion_zeroelim(4, ac, -bez, temp8b);
3976 temp8clen = scale_expansion_zeroelim(4, ab, cez, temp8c);
3977 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
3978 temp8blen, temp8b, temp16);
3979 temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
3980 temp16len, temp16, temp24);
3981 temp48len = scale_expansion_zeroelim(temp24len, temp24, dex, temp48);
3982 xlen = scale_expansion_zeroelim(temp48len, temp48, dex, xdet);
3983 temp48len = scale_expansion_zeroelim(temp24len, temp24, dey, temp48);
3984 ylen = scale_expansion_zeroelim(temp48len, temp48, dey, ydet);
3985 temp48len = scale_expansion_zeroelim(temp24len, temp24, dez, temp48);
3986 zlen = scale_expansion_zeroelim(temp48len, temp48, dez, zdet);
3987 xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
3988 dlen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, ddet);
3989
3990 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
3991 cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
3992 finlength = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, fin1);
3993
3994 det = estimate(finlength, fin1);
3995 errbound = isperrboundB * permanent;
3996 if ((det >= errbound) || (-det >= errbound)) {
3997 return det;
3998 }
3999
4000 Two_Diff_Tail(pa[0], pe[0], aex, aextail);
4001 Two_Diff_Tail(pa[1], pe[1], aey, aeytail);
4002 Two_Diff_Tail(pa[2], pe[2], aez, aeztail);
4003 Two_Diff_Tail(pb[0], pe[0], bex, bextail);
4004 Two_Diff_Tail(pb[1], pe[1], bey, beytail);
4005 Two_Diff_Tail(pb[2], pe[2], bez, beztail);
4006 Two_Diff_Tail(pc[0], pe[0], cex, cextail);
4007 Two_Diff_Tail(pc[1], pe[1], cey, ceytail);
4008 Two_Diff_Tail(pc[2], pe[2], cez, ceztail);
4009 Two_Diff_Tail(pd[0], pe[0], dex, dextail);
4010 Two_Diff_Tail(pd[1], pe[1], dey, deytail);
4011 Two_Diff_Tail(pd[2], pe[2], dez, deztail);
4012 if ((aextail == 0.0) && (aeytail == 0.0) && (aeztail == 0.0)
4013 && (bextail == 0.0) && (beytail == 0.0) && (beztail == 0.0)
4014 && (cextail == 0.0) && (ceytail == 0.0) && (ceztail == 0.0)
4015 && (dextail == 0.0) && (deytail == 0.0) && (deztail == 0.0)) {
4016 return det;
4017 }
4018
4019 errbound = isperrboundC * permanent + resulterrbound * Absolute(det);
4020 abeps = (aex * beytail + bey * aextail)
4021 - (aey * bextail + bex * aeytail);
4022 bceps = (bex * ceytail + cey * bextail)
4023 - (bey * cextail + cex * beytail);
4024 cdeps = (cex * deytail + dey * cextail)
4025 - (cey * dextail + dex * ceytail);
4026 daeps = (dex * aeytail + aey * dextail)
4027 - (dey * aextail + aex * deytail);
4028 aceps = (aex * ceytail + cey * aextail)
4029 - (aey * cextail + cex * aeytail);
4030 bdeps = (bex * deytail + dey * bextail)
4031 - (bey * dextail + dex * beytail);
4032 det += (((bex * bex + bey * bey + bez * bez)
4033 * ((cez * daeps + dez * aceps + aez * cdeps)
4034 + (ceztail * da3 + deztail * ac3 + aeztail * cd3))
4035 + (dex * dex + dey * dey + dez * dez)
4036 * ((aez * bceps - bez * aceps + cez * abeps)
4037 + (aeztail * bc3 - beztail * ac3 + ceztail * ab3)))
4038 - ((aex * aex + aey * aey + aez * aez)
4039 * ((bez * cdeps - cez * bdeps + dez * bceps)
4040 + (beztail * cd3 - ceztail * bd3 + deztail * bc3))
4041 + (cex * cex + cey * cey + cez * cez)
4042 * ((dez * abeps + aez * bdeps + bez * daeps)
4043 + (deztail * ab3 + aeztail * bd3 + beztail * da3))))
4044 + 2.0 * (((bex * bextail + bey * beytail + bez * beztail)
4045 * (cez * da3 + dez * ac3 + aez * cd3)
4046 + (dex * dextail + dey * deytail + dez * deztail)
4047 * (aez * bc3 - bez * ac3 + cez * ab3))
4048 - ((aex * aextail + aey * aeytail + aez * aeztail)
4049 * (bez * cd3 - cez * bd3 + dez * bc3)
4050 + (cex * cextail + cey * ceytail + cez * ceztail)
4051 * (dez * ab3 + aez * bd3 + bez * da3)));
4052 if ((det >= errbound) || (-det >= errbound)) {
4053 return det;
4054 }
4055
4056 return insphereexact(pa, pb, pc, pd, pe);
4057 }
4058
4059 #ifdef USE_CGAL_PREDICATES
4060
insphere(REAL * pa,REAL * pb,REAL * pc,REAL * pd,REAL * pe)4061 REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
4062 {
4063 return (REAL)
4064 - cgal_pred_obj.side_of_oriented_sphere_3_object()
4065 (Point(pa[0], pa[1], pa[2]),
4066 Point(pb[0], pb[1], pb[2]),
4067 Point(pc[0], pc[1], pc[2]),
4068 Point(pd[0], pd[1], pd[2]),
4069 Point(pe[0], pe[1], pe[2]));
4070 }
4071
4072 #else
4073
insphere(REAL * pa,REAL * pb,REAL * pc,REAL * pd,REAL * pe)4074 REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
4075 {
4076 REAL aex, bex, cex, dex;
4077 REAL aey, bey, cey, dey;
4078 REAL aez, bez, cez, dez;
4079 REAL aexbey, bexaey, bexcey, cexbey, cexdey, dexcey, dexaey, aexdey;
4080 REAL aexcey, cexaey, bexdey, dexbey;
4081 REAL alift, blift, clift, dlift;
4082 REAL ab, bc, cd, da, ac, bd;
4083 REAL abc, bcd, cda, dab;
4084 REAL det;
4085
4086
4087 aex = pa[0] - pe[0];
4088 bex = pb[0] - pe[0];
4089 cex = pc[0] - pe[0];
4090 dex = pd[0] - pe[0];
4091 aey = pa[1] - pe[1];
4092 bey = pb[1] - pe[1];
4093 cey = pc[1] - pe[1];
4094 dey = pd[1] - pe[1];
4095 aez = pa[2] - pe[2];
4096 bez = pb[2] - pe[2];
4097 cez = pc[2] - pe[2];
4098 dez = pd[2] - pe[2];
4099
4100 aexbey = aex * bey;
4101 bexaey = bex * aey;
4102 ab = aexbey - bexaey;
4103 bexcey = bex * cey;
4104 cexbey = cex * bey;
4105 bc = bexcey - cexbey;
4106 cexdey = cex * dey;
4107 dexcey = dex * cey;
4108 cd = cexdey - dexcey;
4109 dexaey = dex * aey;
4110 aexdey = aex * dey;
4111 da = dexaey - aexdey;
4112
4113 aexcey = aex * cey;
4114 cexaey = cex * aey;
4115 ac = aexcey - cexaey;
4116 bexdey = bex * dey;
4117 dexbey = dex * bey;
4118 bd = bexdey - dexbey;
4119
4120 abc = aez * bc - bez * ac + cez * ab;
4121 bcd = bez * cd - cez * bd + dez * bc;
4122 cda = cez * da + dez * ac + aez * cd;
4123 dab = dez * ab + aez * bd + bez * da;
4124
4125 alift = aex * aex + aey * aey + aez * aez;
4126 blift = bex * bex + bey * bey + bez * bez;
4127 clift = cex * cex + cey * cey + cez * cez;
4128 dlift = dex * dex + dey * dey + dez * dez;
4129
4130 det = (dlift * abc - clift * dab) + (blift * cda - alift * bcd);
4131
4132 if (_use_inexact_arith) {
4133 return det;
4134 }
4135
4136 if (_use_static_filter) {
4137 if (fabs(det) > ispstaticfilter) return det;
4138 //if (det > ispstaticfilter) return det;
4139 //if (det < minus_ispstaticfilter) return det;
4140
4141 }
4142
4143 REAL aezplus, bezplus, cezplus, dezplus;
4144 REAL aexbeyplus, bexaeyplus, bexceyplus, cexbeyplus;
4145 REAL cexdeyplus, dexceyplus, dexaeyplus, aexdeyplus;
4146 REAL aexceyplus, cexaeyplus, bexdeyplus, dexbeyplus;
4147 REAL permanent, errbound;
4148
4149 aezplus = Absolute(aez);
4150 bezplus = Absolute(bez);
4151 cezplus = Absolute(cez);
4152 dezplus = Absolute(dez);
4153 aexbeyplus = Absolute(aexbey);
4154 bexaeyplus = Absolute(bexaey);
4155 bexceyplus = Absolute(bexcey);
4156 cexbeyplus = Absolute(cexbey);
4157 cexdeyplus = Absolute(cexdey);
4158 dexceyplus = Absolute(dexcey);
4159 dexaeyplus = Absolute(dexaey);
4160 aexdeyplus = Absolute(aexdey);
4161 aexceyplus = Absolute(aexcey);
4162 cexaeyplus = Absolute(cexaey);
4163 bexdeyplus = Absolute(bexdey);
4164 dexbeyplus = Absolute(dexbey);
4165 permanent = ((cexdeyplus + dexceyplus) * bezplus
4166 + (dexbeyplus + bexdeyplus) * cezplus
4167 + (bexceyplus + cexbeyplus) * dezplus)
4168 * alift
4169 + ((dexaeyplus + aexdeyplus) * cezplus
4170 + (aexceyplus + cexaeyplus) * dezplus
4171 + (cexdeyplus + dexceyplus) * aezplus)
4172 * blift
4173 + ((aexbeyplus + bexaeyplus) * dezplus
4174 + (bexdeyplus + dexbeyplus) * aezplus
4175 + (dexaeyplus + aexdeyplus) * bezplus)
4176 * clift
4177 + ((bexceyplus + cexbeyplus) * aezplus
4178 + (cexaeyplus + aexceyplus) * bezplus
4179 + (aexbeyplus + bexaeyplus) * cezplus)
4180 * dlift;
4181 errbound = isperrboundA * permanent;
4182 if ((det > errbound) || (-det > errbound)) {
4183 return det;
4184 }
4185
4186 return insphereadapt(pa, pb, pc, pd, pe, permanent);
4187 }
4188
4189 #endif // #ifdef USE_CGAL_PREDICATES
4190
4191 /*****************************************************************************/
4192 /* */
4193 /* orient4d() Return a positive value if the point pe lies above the */
4194 /* hyperplane passing through pa, pb, pc, and pd; "above" is */
4195 /* defined in a manner best found by trial-and-error. Returns */
4196 /* a negative value if pe lies below the hyperplane. Returns */
4197 /* zero if the points are co-hyperplanar (not affinely */
4198 /* independent). The result is also a rough approximation of */
4199 /* 24 times the signed volume of the 4-simplex defined by the */
4200 /* five points. */
4201 /* */
4202 /* Uses exact arithmetic if necessary to ensure a correct answer. The */
4203 /* result returned is the determinant of a matrix. This determinant is */
4204 /* computed adaptively, in the sense that exact arithmetic is used only to */
4205 /* the degree it is needed to ensure that the returned value has the */
4206 /* correct sign. Hence, orient4d() is usually quite fast, but will run */
4207 /* more slowly when the input points are hyper-coplanar or nearly so. */
4208 /* */
4209 /* See my Robust Predicates paper for details. */
4210 /* */
4211 /*****************************************************************************/
4212
orient4dexact(REAL * pa,REAL * pb,REAL * pc,REAL * pd,REAL * pe,REAL aheight,REAL bheight,REAL cheight,REAL dheight,REAL eheight)4213 REAL orient4dexact(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL* pe,
4214 REAL aheight, REAL bheight, REAL cheight, REAL dheight,
4215 REAL eheight)
4216 {
4217 INEXACT REAL axby1, bxcy1, cxdy1, dxey1, exay1;
4218 INEXACT REAL bxay1, cxby1, dxcy1, exdy1, axey1;
4219 INEXACT REAL axcy1, bxdy1, cxey1, dxay1, exby1;
4220 INEXACT REAL cxay1, dxby1, excy1, axdy1, bxey1;
4221 REAL axby0, bxcy0, cxdy0, dxey0, exay0;
4222 REAL bxay0, cxby0, dxcy0, exdy0, axey0;
4223 REAL axcy0, bxdy0, cxey0, dxay0, exby0;
4224 REAL cxay0, dxby0, excy0, axdy0, bxey0;
4225 REAL ab[4], bc[4], cd[4], de[4], ea[4];
4226 REAL ac[4], bd[4], ce[4], da[4], eb[4];
4227 REAL temp8a[8], temp8b[8], temp16[16];
4228 int temp8alen, temp8blen, temp16len;
4229 REAL abc[24], bcd[24], cde[24], dea[24], eab[24];
4230 REAL abd[24], bce[24], cda[24], deb[24], eac[24];
4231 int abclen, bcdlen, cdelen, dealen, eablen;
4232 int abdlen, bcelen, cdalen, deblen, eaclen;
4233 REAL temp48a[48], temp48b[48];
4234 int temp48alen, temp48blen;
4235 REAL abcd[96], bcde[96], cdea[96], deab[96], eabc[96];
4236 int abcdlen, bcdelen, cdealen, deablen, eabclen;
4237 REAL adet[192], bdet[192], cdet[192], ddet[192], edet[192];
4238 int alen, blen, clen, dlen, elen;
4239 REAL abdet[384], cddet[384], cdedet[576];
4240 int ablen, cdlen;
4241 REAL deter[960];
4242 int deterlen;
4243 int i;
4244
4245 INEXACT REAL bvirt;
4246 REAL avirt, bround, around;
4247 INEXACT REAL c;
4248 INEXACT REAL abig;
4249 REAL ahi, alo, bhi, blo;
4250 REAL err1, err2, err3;
4251 INEXACT REAL _i, _j;
4252 REAL _0;
4253
4254
4255 Two_Product(pa[0], pb[1], axby1, axby0);
4256 Two_Product(pb[0], pa[1], bxay1, bxay0);
4257 Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]);
4258
4259 Two_Product(pb[0], pc[1], bxcy1, bxcy0);
4260 Two_Product(pc[0], pb[1], cxby1, cxby0);
4261 Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]);
4262
4263 Two_Product(pc[0], pd[1], cxdy1, cxdy0);
4264 Two_Product(pd[0], pc[1], dxcy1, dxcy0);
4265 Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]);
4266
4267 Two_Product(pd[0], pe[1], dxey1, dxey0);
4268 Two_Product(pe[0], pd[1], exdy1, exdy0);
4269 Two_Two_Diff(dxey1, dxey0, exdy1, exdy0, de[3], de[2], de[1], de[0]);
4270
4271 Two_Product(pe[0], pa[1], exay1, exay0);
4272 Two_Product(pa[0], pe[1], axey1, axey0);
4273 Two_Two_Diff(exay1, exay0, axey1, axey0, ea[3], ea[2], ea[1], ea[0]);
4274
4275 Two_Product(pa[0], pc[1], axcy1, axcy0);
4276 Two_Product(pc[0], pa[1], cxay1, cxay0);
4277 Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]);
4278
4279 Two_Product(pb[0], pd[1], bxdy1, bxdy0);
4280 Two_Product(pd[0], pb[1], dxby1, dxby0);
4281 Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]);
4282
4283 Two_Product(pc[0], pe[1], cxey1, cxey0);
4284 Two_Product(pe[0], pc[1], excy1, excy0);
4285 Two_Two_Diff(cxey1, cxey0, excy1, excy0, ce[3], ce[2], ce[1], ce[0]);
4286
4287 Two_Product(pd[0], pa[1], dxay1, dxay0);
4288 Two_Product(pa[0], pd[1], axdy1, axdy0);
4289 Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]);
4290
4291 Two_Product(pe[0], pb[1], exby1, exby0);
4292 Two_Product(pb[0], pe[1], bxey1, bxey0);
4293 Two_Two_Diff(exby1, exby0, bxey1, bxey0, eb[3], eb[2], eb[1], eb[0]);
4294
4295 temp8alen = scale_expansion_zeroelim(4, bc, pa[2], temp8a);
4296 temp8blen = scale_expansion_zeroelim(4, ac, -pb[2], temp8b);
4297 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
4298 temp16);
4299 temp8alen = scale_expansion_zeroelim(4, ab, pc[2], temp8a);
4300 abclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
4301 abc);
4302
4303 temp8alen = scale_expansion_zeroelim(4, cd, pb[2], temp8a);
4304 temp8blen = scale_expansion_zeroelim(4, bd, -pc[2], temp8b);
4305 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
4306 temp16);
4307 temp8alen = scale_expansion_zeroelim(4, bc, pd[2], temp8a);
4308 bcdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
4309 bcd);
4310
4311 temp8alen = scale_expansion_zeroelim(4, de, pc[2], temp8a);
4312 temp8blen = scale_expansion_zeroelim(4, ce, -pd[2], temp8b);
4313 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
4314 temp16);
4315 temp8alen = scale_expansion_zeroelim(4, cd, pe[2], temp8a);
4316 cdelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
4317 cde);
4318
4319 temp8alen = scale_expansion_zeroelim(4, ea, pd[2], temp8a);
4320 temp8blen = scale_expansion_zeroelim(4, da, -pe[2], temp8b);
4321 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
4322 temp16);
4323 temp8alen = scale_expansion_zeroelim(4, de, pa[2], temp8a);
4324 dealen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
4325 dea);
4326
4327 temp8alen = scale_expansion_zeroelim(4, ab, pe[2], temp8a);
4328 temp8blen = scale_expansion_zeroelim(4, eb, -pa[2], temp8b);
4329 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
4330 temp16);
4331 temp8alen = scale_expansion_zeroelim(4, ea, pb[2], temp8a);
4332 eablen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
4333 eab);
4334
4335 temp8alen = scale_expansion_zeroelim(4, bd, pa[2], temp8a);
4336 temp8blen = scale_expansion_zeroelim(4, da, pb[2], temp8b);
4337 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
4338 temp16);
4339 temp8alen = scale_expansion_zeroelim(4, ab, pd[2], temp8a);
4340 abdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
4341 abd);
4342
4343 temp8alen = scale_expansion_zeroelim(4, ce, pb[2], temp8a);
4344 temp8blen = scale_expansion_zeroelim(4, eb, pc[2], temp8b);
4345 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
4346 temp16);
4347 temp8alen = scale_expansion_zeroelim(4, bc, pe[2], temp8a);
4348 bcelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
4349 bce);
4350
4351 temp8alen = scale_expansion_zeroelim(4, da, pc[2], temp8a);
4352 temp8blen = scale_expansion_zeroelim(4, ac, pd[2], temp8b);
4353 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
4354 temp16);
4355 temp8alen = scale_expansion_zeroelim(4, cd, pa[2], temp8a);
4356 cdalen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
4357 cda);
4358
4359 temp8alen = scale_expansion_zeroelim(4, eb, pd[2], temp8a);
4360 temp8blen = scale_expansion_zeroelim(4, bd, pe[2], temp8b);
4361 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
4362 temp16);
4363 temp8alen = scale_expansion_zeroelim(4, de, pb[2], temp8a);
4364 deblen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
4365 deb);
4366
4367 temp8alen = scale_expansion_zeroelim(4, ac, pe[2], temp8a);
4368 temp8blen = scale_expansion_zeroelim(4, ce, pa[2], temp8b);
4369 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
4370 temp16);
4371 temp8alen = scale_expansion_zeroelim(4, ea, pc[2], temp8a);
4372 eaclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
4373 eac);
4374
4375 temp48alen = fast_expansion_sum_zeroelim(cdelen, cde, bcelen, bce, temp48a);
4376 temp48blen = fast_expansion_sum_zeroelim(deblen, deb, bcdlen, bcd, temp48b);
4377 for (i = 0; i < temp48blen; i++) {
4378 temp48b[i] = -temp48b[i];
4379 }
4380 bcdelen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
4381 temp48blen, temp48b, bcde);
4382 alen = scale_expansion_zeroelim(bcdelen, bcde, aheight, adet);
4383
4384 temp48alen = fast_expansion_sum_zeroelim(dealen, dea, cdalen, cda, temp48a);
4385 temp48blen = fast_expansion_sum_zeroelim(eaclen, eac, cdelen, cde, temp48b);
4386 for (i = 0; i < temp48blen; i++) {
4387 temp48b[i] = -temp48b[i];
4388 }
4389 cdealen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
4390 temp48blen, temp48b, cdea);
4391 blen = scale_expansion_zeroelim(cdealen, cdea, bheight, bdet);
4392
4393 temp48alen = fast_expansion_sum_zeroelim(eablen, eab, deblen, deb, temp48a);
4394 temp48blen = fast_expansion_sum_zeroelim(abdlen, abd, dealen, dea, temp48b);
4395 for (i = 0; i < temp48blen; i++) {
4396 temp48b[i] = -temp48b[i];
4397 }
4398 deablen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
4399 temp48blen, temp48b, deab);
4400 clen = scale_expansion_zeroelim(deablen, deab, cheight, cdet);
4401
4402 temp48alen = fast_expansion_sum_zeroelim(abclen, abc, eaclen, eac, temp48a);
4403 temp48blen = fast_expansion_sum_zeroelim(bcelen, bce, eablen, eab, temp48b);
4404 for (i = 0; i < temp48blen; i++) {
4405 temp48b[i] = -temp48b[i];
4406 }
4407 eabclen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
4408 temp48blen, temp48b, eabc);
4409 dlen = scale_expansion_zeroelim(eabclen, eabc, dheight, ddet);
4410
4411 temp48alen = fast_expansion_sum_zeroelim(bcdlen, bcd, abdlen, abd, temp48a);
4412 temp48blen = fast_expansion_sum_zeroelim(cdalen, cda, abclen, abc, temp48b);
4413 for (i = 0; i < temp48blen; i++) {
4414 temp48b[i] = -temp48b[i];
4415 }
4416 abcdlen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
4417 temp48blen, temp48b, abcd);
4418 elen = scale_expansion_zeroelim(abcdlen, abcd, eheight, edet);
4419
4420 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
4421 cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
4422 cdelen = fast_expansion_sum_zeroelim(cdlen, cddet, elen, edet, cdedet);
4423 deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdelen, cdedet, deter);
4424
4425 return deter[deterlen - 1];
4426 }
4427
orient4dadapt(REAL * pa,REAL * pb,REAL * pc,REAL * pd,REAL * pe,REAL aheight,REAL bheight,REAL cheight,REAL dheight,REAL eheight,REAL permanent)4428 REAL orient4dadapt(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL* pe,
4429 REAL aheight, REAL bheight, REAL cheight, REAL dheight,
4430 REAL eheight, REAL permanent)
4431 {
4432 INEXACT REAL aex, bex, cex, dex, aey, bey, cey, dey, aez, bez, cez, dez;
4433 INEXACT REAL aeheight, beheight, ceheight, deheight;
4434 REAL det, errbound;
4435
4436 INEXACT REAL aexbey1, bexaey1, bexcey1, cexbey1;
4437 INEXACT REAL cexdey1, dexcey1, dexaey1, aexdey1;
4438 INEXACT REAL aexcey1, cexaey1, bexdey1, dexbey1;
4439 REAL aexbey0, bexaey0, bexcey0, cexbey0;
4440 REAL cexdey0, dexcey0, dexaey0, aexdey0;
4441 REAL aexcey0, cexaey0, bexdey0, dexbey0;
4442 REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4];
4443 INEXACT REAL ab3, bc3, cd3, da3, ac3, bd3;
4444 REAL abeps, bceps, cdeps, daeps, aceps, bdeps;
4445 REAL temp8a[8], temp8b[8], temp8c[8], temp16[16], temp24[24];
4446 int temp8alen, temp8blen, temp8clen, temp16len, temp24len;
4447 REAL adet[48], bdet[48], cdet[48], ddet[48];
4448 int alen, blen, clen, dlen;
4449 REAL abdet[96], cddet[96];
4450 int ablen, cdlen;
4451 REAL fin1[192];
4452 int finlength;
4453
4454 REAL aextail, bextail, cextail, dextail;
4455 REAL aeytail, beytail, ceytail, deytail;
4456 REAL aeztail, beztail, ceztail, deztail;
4457 REAL aeheighttail, beheighttail, ceheighttail, deheighttail;
4458
4459 INEXACT REAL bvirt;
4460 REAL avirt, bround, around;
4461 INEXACT REAL c;
4462 INEXACT REAL abig;
4463 REAL ahi, alo, bhi, blo;
4464 REAL err1, err2, err3;
4465 INEXACT REAL _i, _j;
4466 REAL _0;
4467
4468
4469 aex = (REAL) (pa[0] - pe[0]);
4470 bex = (REAL) (pb[0] - pe[0]);
4471 cex = (REAL) (pc[0] - pe[0]);
4472 dex = (REAL) (pd[0] - pe[0]);
4473 aey = (REAL) (pa[1] - pe[1]);
4474 bey = (REAL) (pb[1] - pe[1]);
4475 cey = (REAL) (pc[1] - pe[1]);
4476 dey = (REAL) (pd[1] - pe[1]);
4477 aez = (REAL) (pa[2] - pe[2]);
4478 bez = (REAL) (pb[2] - pe[2]);
4479 cez = (REAL) (pc[2] - pe[2]);
4480 dez = (REAL) (pd[2] - pe[2]);
4481 aeheight = (REAL) (aheight - eheight);
4482 beheight = (REAL) (bheight - eheight);
4483 ceheight = (REAL) (cheight - eheight);
4484 deheight = (REAL) (dheight - eheight);
4485
4486 Two_Product(aex, bey, aexbey1, aexbey0);
4487 Two_Product(bex, aey, bexaey1, bexaey0);
4488 Two_Two_Diff(aexbey1, aexbey0, bexaey1, bexaey0, ab3, ab[2], ab[1], ab[0]);
4489 ab[3] = ab3;
4490
4491 Two_Product(bex, cey, bexcey1, bexcey0);
4492 Two_Product(cex, bey, cexbey1, cexbey0);
4493 Two_Two_Diff(bexcey1, bexcey0, cexbey1, cexbey0, bc3, bc[2], bc[1], bc[0]);
4494 bc[3] = bc3;
4495
4496 Two_Product(cex, dey, cexdey1, cexdey0);
4497 Two_Product(dex, cey, dexcey1, dexcey0);
4498 Two_Two_Diff(cexdey1, cexdey0, dexcey1, dexcey0, cd3, cd[2], cd[1], cd[0]);
4499 cd[3] = cd3;
4500
4501 Two_Product(dex, aey, dexaey1, dexaey0);
4502 Two_Product(aex, dey, aexdey1, aexdey0);
4503 Two_Two_Diff(dexaey1, dexaey0, aexdey1, aexdey0, da3, da[2], da[1], da[0]);
4504 da[3] = da3;
4505
4506 Two_Product(aex, cey, aexcey1, aexcey0);
4507 Two_Product(cex, aey, cexaey1, cexaey0);
4508 Two_Two_Diff(aexcey1, aexcey0, cexaey1, cexaey0, ac3, ac[2], ac[1], ac[0]);
4509 ac[3] = ac3;
4510
4511 Two_Product(bex, dey, bexdey1, bexdey0);
4512 Two_Product(dex, bey, dexbey1, dexbey0);
4513 Two_Two_Diff(bexdey1, bexdey0, dexbey1, dexbey0, bd3, bd[2], bd[1], bd[0]);
4514 bd[3] = bd3;
4515
4516 temp8alen = scale_expansion_zeroelim(4, cd, bez, temp8a);
4517 temp8blen = scale_expansion_zeroelim(4, bd, -cez, temp8b);
4518 temp8clen = scale_expansion_zeroelim(4, bc, dez, temp8c);
4519 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
4520 temp8blen, temp8b, temp16);
4521 temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
4522 temp16len, temp16, temp24);
4523 alen = scale_expansion_zeroelim(temp24len, temp24, -aeheight, adet);
4524
4525 temp8alen = scale_expansion_zeroelim(4, da, cez, temp8a);
4526 temp8blen = scale_expansion_zeroelim(4, ac, dez, temp8b);
4527 temp8clen = scale_expansion_zeroelim(4, cd, aez, temp8c);
4528 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
4529 temp8blen, temp8b, temp16);
4530 temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
4531 temp16len, temp16, temp24);
4532 blen = scale_expansion_zeroelim(temp24len, temp24, beheight, bdet);
4533
4534 temp8alen = scale_expansion_zeroelim(4, ab, dez, temp8a);
4535 temp8blen = scale_expansion_zeroelim(4, bd, aez, temp8b);
4536 temp8clen = scale_expansion_zeroelim(4, da, bez, temp8c);
4537 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
4538 temp8blen, temp8b, temp16);
4539 temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
4540 temp16len, temp16, temp24);
4541 clen = scale_expansion_zeroelim(temp24len, temp24, -ceheight, cdet);
4542
4543 temp8alen = scale_expansion_zeroelim(4, bc, aez, temp8a);
4544 temp8blen = scale_expansion_zeroelim(4, ac, -bez, temp8b);
4545 temp8clen = scale_expansion_zeroelim(4, ab, cez, temp8c);
4546 temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
4547 temp8blen, temp8b, temp16);
4548 temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
4549 temp16len, temp16, temp24);
4550 dlen = scale_expansion_zeroelim(temp24len, temp24, deheight, ddet);
4551
4552 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
4553 cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
4554 finlength = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, fin1);
4555
4556 det = estimate(finlength, fin1);
4557 errbound = isperrboundB * permanent;
4558 if ((det >= errbound) || (-det >= errbound)) {
4559 return det;
4560 }
4561
4562 Two_Diff_Tail(pa[0], pe[0], aex, aextail);
4563 Two_Diff_Tail(pa[1], pe[1], aey, aeytail);
4564 Two_Diff_Tail(pa[2], pe[2], aez, aeztail);
4565 Two_Diff_Tail(aheight, eheight, aeheight, aeheighttail);
4566 Two_Diff_Tail(pb[0], pe[0], bex, bextail);
4567 Two_Diff_Tail(pb[1], pe[1], bey, beytail);
4568 Two_Diff_Tail(pb[2], pe[2], bez, beztail);
4569 Two_Diff_Tail(bheight, eheight, beheight, beheighttail);
4570 Two_Diff_Tail(pc[0], pe[0], cex, cextail);
4571 Two_Diff_Tail(pc[1], pe[1], cey, ceytail);
4572 Two_Diff_Tail(pc[2], pe[2], cez, ceztail);
4573 Two_Diff_Tail(cheight, eheight, ceheight, ceheighttail);
4574 Two_Diff_Tail(pd[0], pe[0], dex, dextail);
4575 Two_Diff_Tail(pd[1], pe[1], dey, deytail);
4576 Two_Diff_Tail(pd[2], pe[2], dez, deztail);
4577 Two_Diff_Tail(dheight, eheight, deheight, deheighttail);
4578 if ((aextail == 0.0) && (aeytail == 0.0) && (aeztail == 0.0)
4579 && (bextail == 0.0) && (beytail == 0.0) && (beztail == 0.0)
4580 && (cextail == 0.0) && (ceytail == 0.0) && (ceztail == 0.0)
4581 && (dextail == 0.0) && (deytail == 0.0) && (deztail == 0.0)
4582 && (aeheighttail == 0.0) && (beheighttail == 0.0)
4583 && (ceheighttail == 0.0) && (deheighttail == 0.0)) {
4584 return det;
4585 }
4586
4587 errbound = isperrboundC * permanent + resulterrbound * Absolute(det);
4588 abeps = (aex * beytail + bey * aextail)
4589 - (aey * bextail + bex * aeytail);
4590 bceps = (bex * ceytail + cey * bextail)
4591 - (bey * cextail + cex * beytail);
4592 cdeps = (cex * deytail + dey * cextail)
4593 - (cey * dextail + dex * ceytail);
4594 daeps = (dex * aeytail + aey * dextail)
4595 - (dey * aextail + aex * deytail);
4596 aceps = (aex * ceytail + cey * aextail)
4597 - (aey * cextail + cex * aeytail);
4598 bdeps = (bex * deytail + dey * bextail)
4599 - (bey * dextail + dex * beytail);
4600 det += ((beheight
4601 * ((cez * daeps + dez * aceps + aez * cdeps)
4602 + (ceztail * da3 + deztail * ac3 + aeztail * cd3))
4603 + deheight
4604 * ((aez * bceps - bez * aceps + cez * abeps)
4605 + (aeztail * bc3 - beztail * ac3 + ceztail * ab3)))
4606 - (aeheight
4607 * ((bez * cdeps - cez * bdeps + dez * bceps)
4608 + (beztail * cd3 - ceztail * bd3 + deztail * bc3))
4609 + ceheight
4610 * ((dez * abeps + aez * bdeps + bez * daeps)
4611 + (deztail * ab3 + aeztail * bd3 + beztail * da3))))
4612 + ((beheighttail * (cez * da3 + dez * ac3 + aez * cd3)
4613 + deheighttail * (aez * bc3 - bez * ac3 + cez * ab3))
4614 - (aeheighttail * (bez * cd3 - cez * bd3 + dez * bc3)
4615 + ceheighttail * (dez * ab3 + aez * bd3 + bez * da3)));
4616 if ((det >= errbound) || (-det >= errbound)) {
4617 return det;
4618 }
4619
4620 return orient4dexact(pa, pb, pc, pd, pe,
4621 aheight, bheight, cheight, dheight, eheight);
4622 }
4623
orient4d(REAL * pa,REAL * pb,REAL * pc,REAL * pd,REAL * pe,REAL aheight,REAL bheight,REAL cheight,REAL dheight,REAL eheight)4624 REAL orient4d(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL* pe,
4625 REAL aheight, REAL bheight, REAL cheight, REAL dheight,
4626 REAL eheight)
4627 {
4628 REAL aex, bex, cex, dex;
4629 REAL aey, bey, cey, dey;
4630 REAL aez, bez, cez, dez;
4631 REAL aexbey, bexaey, bexcey, cexbey, cexdey, dexcey, dexaey, aexdey;
4632 REAL aexcey, cexaey, bexdey, dexbey;
4633 REAL aeheight, beheight, ceheight, deheight;
4634 REAL ab, bc, cd, da, ac, bd;
4635 REAL abc, bcd, cda, dab;
4636 REAL aezplus, bezplus, cezplus, dezplus;
4637 REAL aexbeyplus, bexaeyplus, bexceyplus, cexbeyplus;
4638 REAL cexdeyplus, dexceyplus, dexaeyplus, aexdeyplus;
4639 REAL aexceyplus, cexaeyplus, bexdeyplus, dexbeyplus;
4640 REAL det;
4641 REAL permanent, errbound;
4642
4643
4644 aex = pa[0] - pe[0];
4645 bex = pb[0] - pe[0];
4646 cex = pc[0] - pe[0];
4647 dex = pd[0] - pe[0];
4648 aey = pa[1] - pe[1];
4649 bey = pb[1] - pe[1];
4650 cey = pc[1] - pe[1];
4651 dey = pd[1] - pe[1];
4652 aez = pa[2] - pe[2];
4653 bez = pb[2] - pe[2];
4654 cez = pc[2] - pe[2];
4655 dez = pd[2] - pe[2];
4656 aeheight = aheight - eheight;
4657 beheight = bheight - eheight;
4658 ceheight = cheight - eheight;
4659 deheight = dheight - eheight;
4660
4661 aexbey = aex * bey;
4662 bexaey = bex * aey;
4663 ab = aexbey - bexaey;
4664 bexcey = bex * cey;
4665 cexbey = cex * bey;
4666 bc = bexcey - cexbey;
4667 cexdey = cex * dey;
4668 dexcey = dex * cey;
4669 cd = cexdey - dexcey;
4670 dexaey = dex * aey;
4671 aexdey = aex * dey;
4672 da = dexaey - aexdey;
4673
4674 aexcey = aex * cey;
4675 cexaey = cex * aey;
4676 ac = aexcey - cexaey;
4677 bexdey = bex * dey;
4678 dexbey = dex * bey;
4679 bd = bexdey - dexbey;
4680
4681 abc = aez * bc - bez * ac + cez * ab;
4682 bcd = bez * cd - cez * bd + dez * bc;
4683 cda = cez * da + dez * ac + aez * cd;
4684 dab = dez * ab + aez * bd + bez * da;
4685
4686 det = (deheight * abc - ceheight * dab) + (beheight * cda - aeheight * bcd);
4687
4688 aezplus = Absolute(aez);
4689 bezplus = Absolute(bez);
4690 cezplus = Absolute(cez);
4691 dezplus = Absolute(dez);
4692 aexbeyplus = Absolute(aexbey);
4693 bexaeyplus = Absolute(bexaey);
4694 bexceyplus = Absolute(bexcey);
4695 cexbeyplus = Absolute(cexbey);
4696 cexdeyplus = Absolute(cexdey);
4697 dexceyplus = Absolute(dexcey);
4698 dexaeyplus = Absolute(dexaey);
4699 aexdeyplus = Absolute(aexdey);
4700 aexceyplus = Absolute(aexcey);
4701 cexaeyplus = Absolute(cexaey);
4702 bexdeyplus = Absolute(bexdey);
4703 dexbeyplus = Absolute(dexbey);
4704 permanent = ((cexdeyplus + dexceyplus) * bezplus
4705 + (dexbeyplus + bexdeyplus) * cezplus
4706 + (bexceyplus + cexbeyplus) * dezplus)
4707 * Absolute(aeheight)
4708 + ((dexaeyplus + aexdeyplus) * cezplus
4709 + (aexceyplus + cexaeyplus) * dezplus
4710 + (cexdeyplus + dexceyplus) * aezplus)
4711 * Absolute(beheight)
4712 + ((aexbeyplus + bexaeyplus) * dezplus
4713 + (bexdeyplus + dexbeyplus) * aezplus
4714 + (dexaeyplus + aexdeyplus) * bezplus)
4715 * Absolute(ceheight)
4716 + ((bexceyplus + cexbeyplus) * aezplus
4717 + (cexaeyplus + aexceyplus) * bezplus
4718 + (aexbeyplus + bexaeyplus) * cezplus)
4719 * Absolute(deheight);
4720 errbound = isperrboundA * permanent;
4721 if ((det > errbound) || (-det > errbound)) {
4722 return det;
4723 }
4724
4725 return orient4dadapt(pa, pb, pc, pd, pe,
4726 aheight, bheight, cheight, dheight, eheight, permanent);
4727 }
4728
4729
4730
4731