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