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