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