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