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