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