1
// Copyright (C) 1999 Jonathan Richard Shewchuk
2
// Licensed under the GNU LGPL Version 2.1.
4
// Modified by Johan Jansson, 2006.
9
#include <dolfin/GeometricPredicates.h>
11
using namespace dolfin;
17
// "predicates_init.h" should be generated for the target architecture.
18
// Look into "rounding.h", why does it need "config.h"?
20
/*****************************************************************************/
22
/* Routines for Arbitrary Precision Floating-point Arithmetic */
23
/* and Fast Robust Geometric Predicates */
28
/* Placed in the public domain by */
29
/* Jonathan Richard Shewchuk */
30
/* School of Computer Science */
31
/* Carnegie Mellon University */
32
/* 5000 Forbes Avenue */
33
/* Pittsburgh, Pennsylvania 15213-3891 */
36
/* This file contains C implementation of algorithms for exact addition */
37
/* and multiplication of floating-point numbers, and predicates for */
38
/* robustly performing the orientation and incircle tests used in */
39
/* computational geometry. The algorithms and underlying theory are */
40
/* described in Jonathan Richard Shewchuk. "Adaptive Precision Floating- */
41
/* Point Arithmetic and Fast Robust Geometric Predicates." Technical */
42
/* Report CMU-CS-96-140, School of Computer Science, Carnegie Mellon */
43
/* University, Pittsburgh, Pennsylvania, May 1996. (Submitted to */
44
/* Discrete & Computational Geometry.) */
46
/* This file, the paper listed above, and other information are available */
47
/* from the Web page http://www.cs.cmu.edu/~quake/robust.html . */
49
/*****************************************************************************/
51
/*****************************************************************************/
53
/* Using this code: */
55
/* First, read the short or long version of the paper (from the Web page */
58
/* Be sure to call exactinit() once, before calling any of the arithmetic */
59
/* functions or geometric predicates. Also be sure to turn on the */
60
/* optimizer when compiling this file. */
63
/* Several geometric predicates are defined. Their parameters are all */
64
/* points. Each point is an array of two or three floating-point */
65
/* numbers. The geometric predicates, described in the papers, are */
67
/* orient2d(pa, pb, pc) */
68
/* orient2dfast(pa, pb, pc) */
69
/* orient3d(pa, pb, pc, pd) */
70
/* orient3dfast(pa, pb, pc, pd) */
71
/* incircle(pa, pb, pc, pd) */
72
/* incirclefast(pa, pb, pc, pd) */
73
/* insphere(pa, pb, pc, pd, pe) */
74
/* inspherefast(pa, pb, pc, pd, pe) */
76
/* Those with suffix "fast" are approximate, non-robust versions. Those */
77
/* without the suffix are adaptive precision, robust versions. There */
78
/* are also versions with the suffices "exact" and "slow", which are */
79
/* non-adaptive, exact arithmetic versions, which I use only for timings */
80
/* in my arithmetic papers. */
83
/* An expansion is represented by an array of floating-point numbers, */
84
/* sorted from smallest to largest magnitude (possibly with interspersed */
85
/* zeros). The length of each expansion is stored as a separate integer, */
86
/* and each arithmetic function returns an integer which is the length */
87
/* of the expansion it created. */
89
/* Several arithmetic functions are defined. Their parameters are */
91
/* e, f Input expansions */
92
/* elen, flen Lengths of input expansions (must be >= 1) */
93
/* h Output expansion */
96
/* The arithmetic functions are */
98
/* grow_expansion(elen, e, b, h) */
99
/* grow_expansion_zeroelim(elen, e, b, h) */
100
/* expansion_sum(elen, e, flen, f, h) */
101
/* expansion_sum_zeroelim1(elen, e, flen, f, h) */
102
/* expansion_sum_zeroelim2(elen, e, flen, f, h) */
103
/* fast_expansion_sum(elen, e, flen, f, h) */
104
/* fast_expansion_sum_zeroelim(elen, e, flen, f, h) */
105
/* linear_expansion_sum(elen, e, flen, f, h) */
106
/* linear_expansion_sum_zeroelim(elen, e, flen, f, h) */
107
/* scale_expansion(elen, e, b, h) */
108
/* scale_expansion_zeroelim(elen, e, b, h) */
109
/* compress(elen, e, h) */
111
/* All of these are described in the long version of the paper; some are */
112
/* described in the short version. All return an integer that is the */
113
/* length of h. Those with suffix _zeroelim perform zero elimination, */
114
/* and are recommended over their counterparts. The procedure */
115
/* fast_expansion_sum_zeroelim() (or linear_expansion_sum_zeroelim() on */
116
/* processors that do not use the round-to-even tiebreaking rule) is */
117
/* recommended over expansion_sum_zeroelim(). Each procedure has a */
118
/* little note next to it (in the code below) that tells you whether or */
119
/* not the output expansion may be the same array as one of the input */
123
/* If you look around below, you'll also find macros for a bunch of */
124
/* simple unrolled arithmetic operations, and procedures for printing */
125
/* expansions (commented out because they don't work with all C */
126
/* compilers) and for generating random floating-point numbers whose */
127
/* significand bits are all random. Most of the macros have undocumented */
128
/* requirements that certain of their parameters should not be the same */
129
/* variable; for safety, better to make sure all the parameters are */
130
/* distinct variables. Feel free to send email to jrs@cs.cmu.edu if you */
131
/* have questions. */
133
/*****************************************************************************/
136
/* Use header file generated automatically by predicates_init. */
137
#define USE_PREDICATES_INIT
139
#ifdef USE_PREDICATES_INIT
140
#include "predicates_init.h"
141
#endif /* USE_PREDICATES_INIT */
143
/* FPU control. We MUST have only double precision (not extended precision) */
144
#include "rounding.h"
146
/* On some machines, the exact arithmetic routines might be defeated by the */
147
/* use of internal extended precision floating-point registers. Sometimes */
148
/* this problem can be fixed by defining certain values to be volatile, */
149
/* thus forcing them to be stored to memory and rounded off. This isn't */
150
/* a great solution, though, as it slows the arithmetic down. */
152
/* To try this out, write "#define INEXACT volatile" below. Normally, */
153
/* however, INEXACT should be defined to be nothing. ("#define INEXACT".) */
155
#define INEXACT /* Nothing */
156
/* #define INEXACT volatile */
158
#define REAL double /* float or double */
159
#define REALPRINT doubleprint
160
#define REALRAND doublerand
161
#define NARROWRAND narrowdoublerand
162
#define UNIFORMRAND uniformdoublerand
164
/* Which of the following two methods of finding the absolute values is */
165
/* fastest is compiler-dependent. A few compilers can inline and optimize */
166
/* the fabs() call; but most will incur the overhead of a function call, */
167
/* which is disastrously slow. A faster way on IEEE machines might be to */
168
/* mask the appropriate bit, but that's difficult to do in C. */
170
#define Absolute(a) ((a) >= 0.0 ? (a) : -(a))
171
/* #define Absolute(a) fabs(a) */
173
/* Many of the operations are broken up into two pieces, a main part that */
174
/* performs an approximate operation, and a "tail" that computes the */
175
/* roundoff error of that operation. */
177
/* The operations Fast_Two_Sum(), Fast_Two_Diff(), Two_Sum(), Two_Diff(), */
178
/* Split(), and Two_Product() are all implemented as described in the */
179
/* reference. Each of these macros requires certain variables to be */
180
/* defined in the calling routine. The variables `bvirt', `c', `abig', */
181
/* `_i', `_j', `_k', `_l', `_m', and `_n' are declared `INEXACT' because */
182
/* they store the result of an operation that may incur roundoff error. */
183
/* The input parameter `x' (or the highest numbered `x_' parameter) must */
184
/* also be declared `INEXACT'. */
186
#define Fast_Two_Sum_Tail(a, b, x, y) \
190
#define Fast_Two_Sum(a, b, x, y) \
191
x = (REAL) (a + b); \
192
Fast_Two_Sum_Tail(a, b, x, y)
194
#define Fast_Two_Diff_Tail(a, b, x, y) \
198
#define Fast_Two_Diff(a, b, x, y) \
199
x = (REAL) (a - b); \
200
Fast_Two_Diff_Tail(a, b, x, y)
202
#define Two_Sum_Tail(a, b, x, y) \
203
bvirt = (REAL) (x - a); \
205
bround = b - bvirt; \
206
around = a - avirt; \
209
#define Two_Sum(a, b, x, y) \
210
x = (REAL) (a + b); \
211
Two_Sum_Tail(a, b, x, y)
213
#define Two_Diff_Tail(a, b, x, y) \
214
bvirt = (REAL) (a - x); \
216
bround = bvirt - b; \
217
around = a - avirt; \
220
#define Two_Diff(a, b, x, y) \
221
x = (REAL) (a - b); \
222
Two_Diff_Tail(a, b, x, y)
224
#define Split(a, ahi, alo) \
225
c = (REAL) (splitter * a); \
226
abig = (REAL) (c - a); \
230
#define Two_Product_Tail(a, b, x, y) \
231
Split(a, ahi, alo); \
232
Split(b, bhi, blo); \
233
err1 = x - (ahi * bhi); \
234
err2 = err1 - (alo * bhi); \
235
err3 = err2 - (ahi * blo); \
236
y = (alo * blo) - err3
238
#define Two_Product(a, b, x, y) \
239
x = (REAL) (a * b); \
240
Two_Product_Tail(a, b, x, y)
242
/* Two_Product_Presplit() is Two_Product() where one of the inputs has */
243
/* already been split. Avoids redundant splitting. */
245
#define Two_Product_Presplit(a, b, bhi, blo, x, y) \
246
x = (REAL) (a * b); \
247
Split(a, ahi, alo); \
248
err1 = x - (ahi * bhi); \
249
err2 = err1 - (alo * bhi); \
250
err3 = err2 - (ahi * blo); \
251
y = (alo * blo) - err3
253
/* Two_Product_2Presplit() is Two_Product() where both of the inputs have */
254
/* already been split. Avoids redundant splitting. */
256
#define Two_Product_2Presplit(a, ahi, alo, b, bhi, blo, x, y) \
257
x = (REAL) (a * b); \
258
err1 = x - (ahi * bhi); \
259
err2 = err1 - (alo * bhi); \
260
err3 = err2 - (ahi * blo); \
261
y = (alo * blo) - err3
263
/* Square() can be done more quickly than Two_Product(). */
265
#define Square_Tail(a, x, y) \
266
Split(a, ahi, alo); \
267
err1 = x - (ahi * ahi); \
268
err3 = err1 - ((ahi + ahi) * alo); \
269
y = (alo * alo) - err3
271
#define Square(a, x, y) \
272
x = (REAL) (a * a); \
275
/* Macros for summing expansions of various fixed lengths. These are all */
276
/* unrolled versions of Expansion_Sum(). */
278
#define Two_One_Sum(a1, a0, b, x2, x1, x0) \
279
Two_Sum(a0, b , _i, x0); \
280
Two_Sum(a1, _i, x2, x1)
282
#define Two_One_Diff(a1, a0, b, x2, x1, x0) \
283
Two_Diff(a0, b , _i, x0); \
284
Two_Sum( a1, _i, x2, x1)
286
#define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \
287
Two_One_Sum(a1, a0, b0, _j, _0, x0); \
288
Two_One_Sum(_j, _0, b1, x3, x2, x1)
290
#define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \
291
Two_One_Diff(a1, a0, b0, _j, _0, x0); \
292
Two_One_Diff(_j, _0, b1, x3, x2, x1)
294
#define Four_One_Sum(a3, a2, a1, a0, b, x4, x3, x2, x1, x0) \
295
Two_One_Sum(a1, a0, b , _j, x1, x0); \
296
Two_One_Sum(a3, a2, _j, x4, x3, x2)
298
#define Four_Two_Sum(a3, a2, a1, a0, b1, b0, x5, x4, x3, x2, x1, x0) \
299
Four_One_Sum(a3, a2, a1, a0, b0, _k, _2, _1, _0, x0); \
300
Four_One_Sum(_k, _2, _1, _0, b1, x5, x4, x3, x2, x1)
302
#define Four_Four_Sum(a3, a2, a1, a0, b4, b3, b1, b0, x7, x6, x5, x4, x3, x2, \
304
Four_Two_Sum(a3, a2, a1, a0, b1, b0, _l, _2, _1, _0, x1, x0); \
305
Four_Two_Sum(_l, _2, _1, _0, b4, b3, x7, x6, x5, x4, x3, x2)
307
#define Eight_One_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b, x8, x7, x6, x5, x4, \
309
Four_One_Sum(a3, a2, a1, a0, b , _j, x3, x2, x1, x0); \
310
Four_One_Sum(a7, a6, a5, a4, _j, x8, x7, x6, x5, x4)
312
#define Eight_Two_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b1, b0, x9, x8, x7, \
313
x6, x5, x4, x3, x2, x1, x0) \
314
Eight_One_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b0, _k, _6, _5, _4, _3, _2, \
316
Eight_One_Sum(_k, _6, _5, _4, _3, _2, _1, _0, b1, x9, x8, x7, x6, x5, x4, \
319
#define Eight_Four_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b4, b3, b1, b0, x11, \
320
x10, x9, x8, x7, x6, x5, x4, x3, x2, x1, x0) \
321
Eight_Two_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b1, b0, _l, _6, _5, _4, _3, \
322
_2, _1, _0, x1, x0); \
323
Eight_Two_Sum(_l, _6, _5, _4, _3, _2, _1, _0, b4, b3, x11, x10, x9, x8, \
324
x7, x6, x5, x4, x3, x2)
326
/* Macros for multiplying expansions of various fixed lengths. */
328
#define Two_One_Product(a1, a0, b, x3, x2, x1, x0) \
329
Split(b, bhi, blo); \
330
Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \
331
Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
332
Two_Sum(_i, _0, _k, x1); \
333
Fast_Two_Sum(_j, _k, x3, x2)
335
#define Four_One_Product(a3, a2, a1, a0, b, x7, x6, x5, x4, x3, x2, x1, x0) \
336
Split(b, bhi, blo); \
337
Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \
338
Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
339
Two_Sum(_i, _0, _k, x1); \
340
Fast_Two_Sum(_j, _k, _i, x2); \
341
Two_Product_Presplit(a2, b, bhi, blo, _j, _0); \
342
Two_Sum(_i, _0, _k, x3); \
343
Fast_Two_Sum(_j, _k, _i, x4); \
344
Two_Product_Presplit(a3, b, bhi, blo, _j, _0); \
345
Two_Sum(_i, _0, _k, x5); \
346
Fast_Two_Sum(_j, _k, x7, x6)
348
#define Two_Two_Product(a1, a0, b1, b0, x7, x6, x5, x4, x3, x2, x1, x0) \
349
Split(a0, a0hi, a0lo); \
350
Split(b0, bhi, blo); \
351
Two_Product_2Presplit(a0, a0hi, a0lo, b0, bhi, blo, _i, x0); \
352
Split(a1, a1hi, a1lo); \
353
Two_Product_2Presplit(a1, a1hi, a1lo, b0, bhi, blo, _j, _0); \
354
Two_Sum(_i, _0, _k, _1); \
355
Fast_Two_Sum(_j, _k, _l, _2); \
356
Split(b1, bhi, blo); \
357
Two_Product_2Presplit(a0, a0hi, a0lo, b1, bhi, blo, _i, _0); \
358
Two_Sum(_1, _0, _k, x1); \
359
Two_Sum(_2, _k, _j, _1); \
360
Two_Sum(_l, _j, _m, _2); \
361
Two_Product_2Presplit(a1, a1hi, a1lo, b1, bhi, blo, _j, _0); \
362
Two_Sum(_i, _0, _n, _0); \
363
Two_Sum(_1, _0, _i, x2); \
364
Two_Sum(_2, _i, _k, _1); \
365
Two_Sum(_m, _k, _l, _2); \
366
Two_Sum(_j, _n, _k, _0); \
367
Two_Sum(_1, _0, _j, x3); \
368
Two_Sum(_2, _j, _i, _1); \
369
Two_Sum(_l, _i, _m, _2); \
370
Two_Sum(_1, _k, _i, x4); \
371
Two_Sum(_2, _i, _k, x5); \
372
Two_Sum(_m, _k, x7, x6)
374
/* An expansion of length two can be squared more quickly than finding the */
375
/* product of two different expansions of length two, and the result is */
376
/* guaranteed to have no more than six (rather than eight) components. */
378
#define Two_Square(a1, a0, x5, x4, x3, x2, x1, x0) \
379
Square(a0, _j, x0); \
381
Two_Product(a1, _0, _k, _1); \
382
Two_One_Sum(_k, _1, _j, _l, _2, x1); \
383
Square(a1, _j, _1); \
384
Two_Two_Sum(_j, _1, _l, _2, x5, x4, x3, x2)
386
#ifndef USE_PREDICATES_INIT
388
static REAL splitter; /* = 2^ceiling(p / 2) + 1. Used to split floats in half. */
389
/* A set of coefficients used to calculate maximum roundoff errors. */
390
static REAL resulterrbound;
391
static REAL ccwerrboundA, ccwerrboundB, ccwerrboundC;
392
static REAL o3derrboundA, o3derrboundB, o3derrboundC;
393
static REAL iccerrboundA, iccerrboundB, iccerrboundC;
394
static REAL isperrboundA, isperrboundB, isperrboundC;
396
#endif /* USE_PREDICATES_INIT */
398
/*****************************************************************************/
400
/* doubleprint() Print the bit representation of a double. */
402
/* Useful for debugging exact arithmetic routines. */
404
/*****************************************************************************/
407
void doubleprint(number)
410
unsigned long long no;
411
unsigned long long sign, expo;
415
no = *(unsigned long long *) &number;
416
sign = no & 0x8000000000000000ll;
417
expo = (no >> 52) & 0x7ffll;
418
exponent = (int) expo;
419
exponent = exponent - 1023;
425
if (exponent == -1023) {
427
"0.0000000000000000000000000000000000000000000000000000_ ( )");
431
for (i = 0; i < 52; i++) {
432
if (no & 0x0008000000000000ll) {
440
printf("_%d (%d)", exponent, exponent - 1 - bottomi);
445
/*****************************************************************************/
447
/* floatprint() Print the bit representation of a float. */
449
/* Useful for debugging exact arithmetic routines. */
451
/*****************************************************************************/
454
void floatprint(number)
462
no = *(unsigned *) &number;
463
sign = no & 0x80000000;
464
expo = (no >> 23) & 0xff;
465
exponent = (int) expo;
466
exponent = exponent - 127;
472
if (exponent == -127) {
473
printf("0.00000000000000000000000_ ( )");
477
for (i = 0; i < 23; i++) {
478
if (no & 0x00400000) {
486
printf("_%3d (%3d)", exponent, exponent - 1 - bottomi);
491
/*****************************************************************************/
493
/* expansion_print() Print the bit representation of an expansion. */
495
/* Useful for debugging exact arithmetic routines. */
497
/*****************************************************************************/
500
void expansion_print(elen, e)
506
for (i = elen - 1; i >= 0; i--) {
517
/*****************************************************************************/
519
/* doublerand() Generate a double with random 53-bit significand and a */
520
/* random exponent in [0, 511]. */
522
/*****************************************************************************/
525
static double doublerand()
535
result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8);
536
for (i = 512, expo = 2; i <= 131072; i *= 2, expo = expo * expo) {
545
/*****************************************************************************/
547
/* narrowdoublerand() Generate a double with random 53-bit significand */
548
/* and a random exponent in [0, 7]. */
550
/*****************************************************************************/
553
static double narrowdoublerand()
563
result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8);
564
for (i = 512, expo = 2; i <= 2048; i *= 2, expo = expo * expo) {
573
/*****************************************************************************/
575
/* uniformdoublerand() Generate a double with random 53-bit significand. */
577
/*****************************************************************************/
580
static double uniformdoublerand()
587
result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8);
592
/*****************************************************************************/
594
/* floatrand() Generate a float with random 24-bit significand and a */
595
/* random exponent in [0, 63]. */
597
/*****************************************************************************/
600
static float floatrand()
609
result = (float) ((a - 1073741824) >> 6);
610
for (i = 512, expo = 2; i <= 16384; i *= 2, expo = expo * expo) {
619
/*****************************************************************************/
621
/* narrowfloatrand() Generate a float with random 24-bit significand and */
622
/* a random exponent in [0, 7]. */
624
/*****************************************************************************/
627
static float narrowfloatrand()
636
result = (float) ((a - 1073741824) >> 6);
637
for (i = 512, expo = 2; i <= 2048; i *= 2, expo = expo * expo) {
646
/*****************************************************************************/
648
/* uniformfloatrand() Generate a float with random 24-bit significand. */
650
/*****************************************************************************/
653
static float uniformfloatrand()
659
result = (float) ((a - 1073741824) >> 6);
664
/*****************************************************************************/
666
/* fast_expansion_sum_zeroelim() Sum two expansions, eliminating zero */
667
/* components from the output expansion. */
669
/* Sets h = e + f. See the long version of my paper for details. */
671
/* If round-to-even is used (as with IEEE 754), maintains the strongly */
672
/* nonoverlapping property. (That is, if e is strongly nonoverlapping, h */
673
/* will be also.) Does NOT maintain the nonoverlapping or nonadjacent */
676
/*****************************************************************************/
678
static int fast_expansion_sum_zeroelim(int elen, REAL *e,
679
int flen, REAL *f, REAL *h)
680
/* h cannot be e or f. */
686
REAL avirt, bround, around;
687
int eindex, findex, hindex;
693
if ((fnow > enow) == (fnow > -enow)) {
701
if ((eindex < elen) && (findex < flen)) {
702
if ((fnow > enow) == (fnow > -enow)) {
703
Fast_Two_Sum(enow, Q, Qnew, hh);
706
Fast_Two_Sum(fnow, Q, Qnew, hh);
713
while ((eindex < elen) && (findex < flen)) {
714
if ((fnow > enow) == (fnow > -enow)) {
715
Two_Sum(Q, enow, Qnew, hh);
718
Two_Sum(Q, fnow, Qnew, hh);
727
while (eindex < elen) {
728
Two_Sum(Q, enow, Qnew, hh);
735
while (findex < flen) {
736
Two_Sum(Q, fnow, Qnew, hh);
743
if ((Q != 0.0) || (hindex == 0)) {
749
/*****************************************************************************/
751
/* scale_expansion_zeroelim() Multiply an expansion by a scalar, */
752
/* eliminating zero components from the */
753
/* output expansion. */
755
/* Sets h = be. See either version of my paper for details. */
757
/* Maintains the nonoverlapping property. If round-to-even is used (as */
758
/* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */
759
/* properties as well. (That is, if e has one of these properties, so */
762
/*****************************************************************************/
764
static int scale_expansion_zeroelim(int elen, REAL *e, REAL b, REAL *h)
765
/* e and h cannot be the same. */
769
INEXACT REAL product1;
774
REAL avirt, bround, around;
777
REAL ahi, alo, bhi, blo;
778
REAL err1, err2, err3;
781
Two_Product_Presplit(e[0], b, bhi, blo, Q, hh);
786
for (eindex = 1; eindex < elen; eindex++) {
788
Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
789
Two_Sum(Q, product0, sum, hh);
793
Fast_Two_Sum(product1, sum, Q, hh);
798
if ((Q != 0.0) || (hindex == 0)) {
804
/*****************************************************************************/
806
/* estimate() Produce a one-word estimate of an expansion's value. */
808
/* See either version of my paper for details. */
810
/*****************************************************************************/
812
static REAL estimate(int elen, REAL *e)
818
for (eindex = 1; eindex < elen; eindex++) {
824
/*****************************************************************************/
826
/* orient2dfast() Approximate 2D orientation test. Nonrobust. */
827
/* orient2dexact() Exact 2D orientation test. Robust. */
828
/* orient2dslow() Another exact 2D orientation test. Robust. */
829
/* orient2d() Adaptive exact 2D orientation test. Robust. */
831
/* Return a positive value if the points pa, pb, and pc occur */
832
/* in counterclockwise order; a negative value if they occur */
833
/* in clockwise order; and zero if they are collinear. The */
834
/* result is also a rough approximation of twice the signed */
835
/* area of the triangle defined by the three points. */
837
/* Only the first and last routine should be used; the middle two are for */
840
/* The last three use exact arithmetic to ensure a correct answer. The */
841
/* result returned is the determinant of a matrix. In orient2d() only, */
842
/* this determinant is computed adaptively, in the sense that exact */
843
/* arithmetic is used only to the degree it is needed to ensure that the */
844
/* returned value has the correct sign. Hence, orient2d() is usually quite */
845
/* fast, but will run more slowly when the input points are collinear or */
848
/*****************************************************************************/
850
static REAL orient2dadapt(REAL *pa, REAL *pb, REAL *pc, REAL detsum)
852
INEXACT REAL acx, acy, bcx, bcy;
853
REAL acxtail, acytail, bcxtail, bcytail;
854
INEXACT REAL detleft, detright;
855
REAL detlefttail, detrighttail;
857
REAL B[4], C1[8], C2[12], D[16];
859
int C1length, C2length, Dlength;
866
REAL avirt, bround, around;
869
REAL ahi, alo, bhi, blo;
870
REAL err1, err2, err3;
874
acx = (REAL) (pa[0] - pc[0]);
875
bcx = (REAL) (pb[0] - pc[0]);
876
acy = (REAL) (pa[1] - pc[1]);
877
bcy = (REAL) (pb[1] - pc[1]);
879
Two_Product(acx, bcy, detleft, detlefttail);
880
Two_Product(acy, bcx, detright, detrighttail);
882
Two_Two_Diff(detleft, detlefttail, detright, detrighttail,
883
B3, B[2], B[1], B[0]);
886
det = estimate(4, B);
887
errbound = ccwerrboundB * detsum;
888
if ((det >= errbound) || (-det >= errbound)) {
892
Two_Diff_Tail(pa[0], pc[0], acx, acxtail);
893
Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail);
894
Two_Diff_Tail(pa[1], pc[1], acy, acytail);
895
Two_Diff_Tail(pb[1], pc[1], bcy, bcytail);
897
if ((acxtail == 0.0) && (acytail == 0.0)
898
&& (bcxtail == 0.0) && (bcytail == 0.0)) {
902
errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det);
903
det += (acx * bcytail + bcy * acxtail)
904
- (acy * bcxtail + bcx * acytail);
905
if ((det >= errbound) || (-det >= errbound)) {
909
Two_Product(acxtail, bcy, s1, s0);
910
Two_Product(acytail, bcx, t1, t0);
911
Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
913
C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1);
915
Two_Product(acx, bcytail, s1, s0);
916
Two_Product(acy, bcxtail, t1, t0);
917
Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
919
C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2);
921
Two_Product(acxtail, bcytail, s1, s0);
922
Two_Product(acytail, bcxtail, t1, t0);
923
Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
925
Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D);
927
return(D[Dlength - 1]);
930
REAL orient2d(REAL *pa, REAL *pb, REAL *pc)
932
REAL detleft, detright, det;
933
REAL detsum, errbound;
938
detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]);
939
detright = (pa[1] - pc[1]) * (pb[0] - pc[0]);
940
det = detleft - detright;
943
if (detright <= 0.0) {
947
detsum = detleft + detright;
949
} else if (detleft < 0.0) {
950
if (detright >= 0.0) {
954
detsum = -detleft - detright;
961
errbound = ccwerrboundA * detsum;
962
if ((det >= errbound) || (-det >= errbound)) {
967
orient = orient2dadapt(pa, pb, pc, detsum);
972
/*****************************************************************************/
974
/* orient3dfast() Approximate 3D orientation test. Nonrobust. */
975
/* orient3dexact() Exact 3D orientation test. Robust. */
976
/* orient3dslow() Another exact 3D orientation test. Robust. */
977
/* orient3d() Adaptive exact 3D orientation test. Robust. */
979
/* Return a positive value if the point pd lies below the */
980
/* plane passing through pa, pb, and pc; "below" is defined so */
981
/* that pa, pb, and pc appear in counterclockwise order when */
982
/* viewed from above the plane. Returns a negative value if */
983
/* pd lies above the plane. Returns zero if the points are */
984
/* coplanar. The result is also a rough approximation of six */
985
/* times the signed volume of the tetrahedron defined by the */
988
/* Only the first and last routine should be used; the middle two are for */
991
/* The last three use exact arithmetic to ensure a correct answer. The */
992
/* result returned is the determinant of a matrix. In orient3d() only, */
993
/* this determinant is computed adaptively, in the sense that exact */
994
/* arithmetic is used only to the degree it is needed to ensure that the */
995
/* returned value has the correct sign. Hence, orient3d() is usually quite */
996
/* fast, but will run more slowly when the input points are coplanar or */
999
/*****************************************************************************/
1001
static REAL orient3dadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd,
1004
INEXACT REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
1007
INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
1008
REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
1009
REAL bc[4], ca[4], ab[4];
1010
INEXACT REAL bc3, ca3, ab3;
1011
REAL adet[8], bdet[8], cdet[8];
1012
int alen, blen, clen;
1015
REAL *finnow, *finother, *finswap;
1016
REAL fin1[192], fin2[192];
1019
REAL adxtail, bdxtail, cdxtail;
1020
REAL adytail, bdytail, cdytail;
1021
REAL adztail, bdztail, cdztail;
1022
INEXACT REAL at_blarge, at_clarge;
1023
INEXACT REAL bt_clarge, bt_alarge;
1024
INEXACT REAL ct_alarge, ct_blarge;
1025
REAL at_b[4], at_c[4], bt_c[4], bt_a[4], ct_a[4], ct_b[4];
1026
int at_blen, at_clen, bt_clen, bt_alen, ct_alen, ct_blen;
1027
INEXACT REAL bdxt_cdy1, cdxt_bdy1, cdxt_ady1;
1028
INEXACT REAL adxt_cdy1, adxt_bdy1, bdxt_ady1;
1029
REAL bdxt_cdy0, cdxt_bdy0, cdxt_ady0;
1030
REAL adxt_cdy0, adxt_bdy0, bdxt_ady0;
1031
INEXACT REAL bdyt_cdx1, cdyt_bdx1, cdyt_adx1;
1032
INEXACT REAL adyt_cdx1, adyt_bdx1, bdyt_adx1;
1033
REAL bdyt_cdx0, cdyt_bdx0, cdyt_adx0;
1034
REAL adyt_cdx0, adyt_bdx0, bdyt_adx0;
1035
REAL bct[8], cat[8], abt[8];
1036
int bctlen, catlen, abtlen;
1037
INEXACT REAL bdxt_cdyt1, cdxt_bdyt1, cdxt_adyt1;
1038
INEXACT REAL adxt_cdyt1, adxt_bdyt1, bdxt_adyt1;
1039
REAL bdxt_cdyt0, cdxt_bdyt0, cdxt_adyt0;
1040
REAL adxt_cdyt0, adxt_bdyt0, bdxt_adyt0;
1041
REAL u[4], v[12], w[16];
1043
int vlength, wlength;
1047
REAL avirt, bround, around;
1050
REAL ahi, alo, bhi, blo;
1051
REAL err1, err2, err3;
1052
INEXACT REAL _i, _j, _k;
1055
adx = (REAL) (pa[0] - pd[0]);
1056
bdx = (REAL) (pb[0] - pd[0]);
1057
cdx = (REAL) (pc[0] - pd[0]);
1058
ady = (REAL) (pa[1] - pd[1]);
1059
bdy = (REAL) (pb[1] - pd[1]);
1060
cdy = (REAL) (pc[1] - pd[1]);
1061
adz = (REAL) (pa[2] - pd[2]);
1062
bdz = (REAL) (pb[2] - pd[2]);
1063
cdz = (REAL) (pc[2] - pd[2]);
1065
Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
1066
Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
1067
Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
1069
alen = scale_expansion_zeroelim(4, bc, adz, adet);
1071
Two_Product(cdx, ady, cdxady1, cdxady0);
1072
Two_Product(adx, cdy, adxcdy1, adxcdy0);
1073
Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
1075
blen = scale_expansion_zeroelim(4, ca, bdz, bdet);
1077
Two_Product(adx, bdy, adxbdy1, adxbdy0);
1078
Two_Product(bdx, ady, bdxady1, bdxady0);
1079
Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
1081
clen = scale_expansion_zeroelim(4, ab, cdz, cdet);
1083
ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
1084
finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
1086
det = estimate(finlength, fin1);
1087
errbound = o3derrboundB * permanent;
1088
if ((det >= errbound) || (-det >= errbound)) {
1092
Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
1093
Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
1094
Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
1095
Two_Diff_Tail(pa[1], pd[1], ady, adytail);
1096
Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
1097
Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
1098
Two_Diff_Tail(pa[2], pd[2], adz, adztail);
1099
Two_Diff_Tail(pb[2], pd[2], bdz, bdztail);
1100
Two_Diff_Tail(pc[2], pd[2], cdz, cdztail);
1102
if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
1103
&& (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)
1104
&& (adztail == 0.0) && (bdztail == 0.0) && (cdztail == 0.0)) {
1108
errbound = o3derrboundC * permanent + resulterrbound * Absolute(det);
1109
det += (adz * ((bdx * cdytail + cdy * bdxtail)
1110
- (bdy * cdxtail + cdx * bdytail))
1111
+ adztail * (bdx * cdy - bdy * cdx))
1112
+ (bdz * ((cdx * adytail + ady * cdxtail)
1113
- (cdy * adxtail + adx * cdytail))
1114
+ bdztail * (cdx * ady - cdy * adx))
1115
+ (cdz * ((adx * bdytail + bdy * adxtail)
1116
- (ady * bdxtail + bdx * adytail))
1117
+ cdztail * (adx * bdy - ady * bdx));
1118
if ((det >= errbound) || (-det >= errbound)) {
1125
if (adxtail == 0.0) {
1126
if (adytail == 0.0) {
1133
Two_Product(negate, bdx, at_blarge, at_b[0]);
1134
at_b[1] = at_blarge;
1136
Two_Product(adytail, cdx, at_clarge, at_c[0]);
1137
at_c[1] = at_clarge;
1141
if (adytail == 0.0) {
1142
Two_Product(adxtail, bdy, at_blarge, at_b[0]);
1143
at_b[1] = at_blarge;
1146
Two_Product(negate, cdy, at_clarge, at_c[0]);
1147
at_c[1] = at_clarge;
1150
Two_Product(adxtail, bdy, adxt_bdy1, adxt_bdy0);
1151
Two_Product(adytail, bdx, adyt_bdx1, adyt_bdx0);
1152
Two_Two_Diff(adxt_bdy1, adxt_bdy0, adyt_bdx1, adyt_bdx0,
1153
at_blarge, at_b[2], at_b[1], at_b[0]);
1154
at_b[3] = at_blarge;
1156
Two_Product(adytail, cdx, adyt_cdx1, adyt_cdx0);
1157
Two_Product(adxtail, cdy, adxt_cdy1, adxt_cdy0);
1158
Two_Two_Diff(adyt_cdx1, adyt_cdx0, adxt_cdy1, adxt_cdy0,
1159
at_clarge, at_c[2], at_c[1], at_c[0]);
1160
at_c[3] = at_clarge;
1164
if (bdxtail == 0.0) {
1165
if (bdytail == 0.0) {
1172
Two_Product(negate, cdx, bt_clarge, bt_c[0]);
1173
bt_c[1] = bt_clarge;
1175
Two_Product(bdytail, adx, bt_alarge, bt_a[0]);
1176
bt_a[1] = bt_alarge;
1180
if (bdytail == 0.0) {
1181
Two_Product(bdxtail, cdy, bt_clarge, bt_c[0]);
1182
bt_c[1] = bt_clarge;
1185
Two_Product(negate, ady, bt_alarge, bt_a[0]);
1186
bt_a[1] = bt_alarge;
1189
Two_Product(bdxtail, cdy, bdxt_cdy1, bdxt_cdy0);
1190
Two_Product(bdytail, cdx, bdyt_cdx1, bdyt_cdx0);
1191
Two_Two_Diff(bdxt_cdy1, bdxt_cdy0, bdyt_cdx1, bdyt_cdx0,
1192
bt_clarge, bt_c[2], bt_c[1], bt_c[0]);
1193
bt_c[3] = bt_clarge;
1195
Two_Product(bdytail, adx, bdyt_adx1, bdyt_adx0);
1196
Two_Product(bdxtail, ady, bdxt_ady1, bdxt_ady0);
1197
Two_Two_Diff(bdyt_adx1, bdyt_adx0, bdxt_ady1, bdxt_ady0,
1198
bt_alarge, bt_a[2], bt_a[1], bt_a[0]);
1199
bt_a[3] = bt_alarge;
1203
if (cdxtail == 0.0) {
1204
if (cdytail == 0.0) {
1211
Two_Product(negate, adx, ct_alarge, ct_a[0]);
1212
ct_a[1] = ct_alarge;
1214
Two_Product(cdytail, bdx, ct_blarge, ct_b[0]);
1215
ct_b[1] = ct_blarge;
1219
if (cdytail == 0.0) {
1220
Two_Product(cdxtail, ady, ct_alarge, ct_a[0]);
1221
ct_a[1] = ct_alarge;
1224
Two_Product(negate, bdy, ct_blarge, ct_b[0]);
1225
ct_b[1] = ct_blarge;
1228
Two_Product(cdxtail, ady, cdxt_ady1, cdxt_ady0);
1229
Two_Product(cdytail, adx, cdyt_adx1, cdyt_adx0);
1230
Two_Two_Diff(cdxt_ady1, cdxt_ady0, cdyt_adx1, cdyt_adx0,
1231
ct_alarge, ct_a[2], ct_a[1], ct_a[0]);
1232
ct_a[3] = ct_alarge;
1234
Two_Product(cdytail, bdx, cdyt_bdx1, cdyt_bdx0);
1235
Two_Product(cdxtail, bdy, cdxt_bdy1, cdxt_bdy0);
1236
Two_Two_Diff(cdyt_bdx1, cdyt_bdx0, cdxt_bdy1, cdxt_bdy0,
1237
ct_blarge, ct_b[2], ct_b[1], ct_b[0]);
1238
ct_b[3] = ct_blarge;
1243
bctlen = fast_expansion_sum_zeroelim(bt_clen, bt_c, ct_blen, ct_b, bct);
1244
wlength = scale_expansion_zeroelim(bctlen, bct, adz, w);
1245
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
1247
finswap = finnow; finnow = finother; finother = finswap;
1249
catlen = fast_expansion_sum_zeroelim(ct_alen, ct_a, at_clen, at_c, cat);
1250
wlength = scale_expansion_zeroelim(catlen, cat, bdz, w);
1251
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
1253
finswap = finnow; finnow = finother; finother = finswap;
1255
abtlen = fast_expansion_sum_zeroelim(at_blen, at_b, bt_alen, bt_a, abt);
1256
wlength = scale_expansion_zeroelim(abtlen, abt, cdz, w);
1257
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
1259
finswap = finnow; finnow = finother; finother = finswap;
1261
if (adztail != 0.0) {
1262
vlength = scale_expansion_zeroelim(4, bc, adztail, v);
1263
finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
1265
finswap = finnow; finnow = finother; finother = finswap;
1267
if (bdztail != 0.0) {
1268
vlength = scale_expansion_zeroelim(4, ca, bdztail, v);
1269
finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
1271
finswap = finnow; finnow = finother; finother = finswap;
1273
if (cdztail != 0.0) {
1274
vlength = scale_expansion_zeroelim(4, ab, cdztail, v);
1275
finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
1277
finswap = finnow; finnow = finother; finother = finswap;
1280
if (adxtail != 0.0) {
1281
if (bdytail != 0.0) {
1282
Two_Product(adxtail, bdytail, adxt_bdyt1, adxt_bdyt0);
1283
Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdz, u3, u[2], u[1], u[0]);
1285
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
1287
finswap = finnow; finnow = finother; finother = finswap;
1288
if (cdztail != 0.0) {
1289
Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdztail, u3, u[2], u[1], u[0]);
1291
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
1293
finswap = finnow; finnow = finother; finother = finswap;
1296
if (cdytail != 0.0) {
1298
Two_Product(negate, cdytail, adxt_cdyt1, adxt_cdyt0);
1299
Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdz, u3, u[2], u[1], u[0]);
1301
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
1303
finswap = finnow; finnow = finother; finother = finswap;
1304
if (bdztail != 0.0) {
1305
Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdztail, u3, u[2], u[1], u[0]);
1307
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
1309
finswap = finnow; finnow = finother; finother = finswap;
1313
if (bdxtail != 0.0) {
1314
if (cdytail != 0.0) {
1315
Two_Product(bdxtail, cdytail, bdxt_cdyt1, bdxt_cdyt0);
1316
Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adz, u3, u[2], u[1], u[0]);
1318
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
1320
finswap = finnow; finnow = finother; finother = finswap;
1321
if (adztail != 0.0) {
1322
Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adztail, u3, u[2], u[1], u[0]);
1324
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
1326
finswap = finnow; finnow = finother; finother = finswap;
1329
if (adytail != 0.0) {
1331
Two_Product(negate, adytail, bdxt_adyt1, bdxt_adyt0);
1332
Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdz, u3, u[2], u[1], u[0]);
1334
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
1336
finswap = finnow; finnow = finother; finother = finswap;
1337
if (cdztail != 0.0) {
1338
Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdztail, u3, u[2], u[1], u[0]);
1340
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
1342
finswap = finnow; finnow = finother; finother = finswap;
1346
if (cdxtail != 0.0) {
1347
if (adytail != 0.0) {
1348
Two_Product(cdxtail, adytail, cdxt_adyt1, cdxt_adyt0);
1349
Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdz, u3, u[2], u[1], u[0]);
1351
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
1353
finswap = finnow; finnow = finother; finother = finswap;
1354
if (bdztail != 0.0) {
1355
Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdztail, u3, u[2], u[1], u[0]);
1357
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
1359
finswap = finnow; finnow = finother; finother = finswap;
1362
if (bdytail != 0.0) {
1364
Two_Product(negate, bdytail, cdxt_bdyt1, cdxt_bdyt0);
1365
Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adz, u3, u[2], u[1], u[0]);
1367
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
1369
finswap = finnow; finnow = finother; finother = finswap;
1370
if (adztail != 0.0) {
1371
Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adztail, u3, u[2], u[1], u[0]);
1373
finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
1375
finswap = finnow; finnow = finother; finother = finswap;
1380
if (adztail != 0.0) {
1381
wlength = scale_expansion_zeroelim(bctlen, bct, adztail, w);
1382
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
1384
finswap = finnow; finnow = finother; finother = finswap;
1386
if (bdztail != 0.0) {
1387
wlength = scale_expansion_zeroelim(catlen, cat, bdztail, w);
1388
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
1390
finswap = finnow; finnow = finother; finother = finswap;
1392
if (cdztail != 0.0) {
1393
wlength = scale_expansion_zeroelim(abtlen, abt, cdztail, w);
1394
finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
1396
finswap = finnow; finnow = finother; finother = finswap;
1399
return finnow[finlength - 1];
1402
REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
1404
REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
1405
REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
1407
REAL permanent, errbound;
1412
adx = pa[0] - pd[0];
1413
bdx = pb[0] - pd[0];
1414
cdx = pc[0] - pd[0];
1415
ady = pa[1] - pd[1];
1416
bdy = pb[1] - pd[1];
1417
cdy = pc[1] - pd[1];
1418
adz = pa[2] - pd[2];
1419
bdz = pb[2] - pd[2];
1420
cdz = pc[2] - pd[2];
1431
det = adz * (bdxcdy - cdxbdy)
1432
+ bdz * (cdxady - adxcdy)
1433
+ cdz * (adxbdy - bdxady);
1435
permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * Absolute(adz)
1436
+ (Absolute(cdxady) + Absolute(adxcdy)) * Absolute(bdz)
1437
+ (Absolute(adxbdy) + Absolute(bdxady)) * Absolute(cdz);
1438
errbound = o3derrboundA * permanent;
1439
if ((det > errbound) || (-det > errbound)) {
1444
orient = orient3dadapt(pa, pb, pc, pd, permanent);
1449
/*****************************************************************************/
1451
/* incirclefast() Approximate 2D incircle test. Nonrobust. */
1452
/* incircleexact() Exact 2D incircle test. Robust. */
1453
/* incircleslow() Another exact 2D incircle test. Robust. */
1454
/* incircle() Adaptive exact 2D incircle test. Robust. */
1456
/* Return a positive value if the point pd lies inside the */
1457
/* circle passing through pa, pb, and pc; a negative value if */
1458
/* it lies outside; and zero if the four points are cocircular.*/
1459
/* The points pa, pb, and pc must be in counterclockwise */
1460
/* order, or the sign of the result will be reversed. */
1462
/* Only the first and last routine should be used; the middle two are for */
1465
/* The last three use exact arithmetic to ensure a correct answer. The */
1466
/* result returned is the determinant of a matrix. In incircle() only, */
1467
/* this determinant is computed adaptively, in the sense that exact */
1468
/* arithmetic is used only to the degree it is needed to ensure that the */
1469
/* returned value has the correct sign. Hence, incircle() is usually quite */
1470
/* fast, but will run more slowly when the input points are cocircular or */
1473
/*****************************************************************************/
1475
static REAL incircleadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd,
1478
INEXACT REAL adx, bdx, cdx, ady, bdy, cdy;
1481
INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
1482
REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
1483
REAL bc[4], ca[4], ab[4];
1484
INEXACT REAL bc3, ca3, ab3;
1485
REAL axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32];
1486
int axbclen, axxbclen, aybclen, ayybclen, alen;
1487
REAL bxca[8], bxxca[16], byca[8], byyca[16], bdet[32];
1488
int bxcalen, bxxcalen, bycalen, byycalen, blen;
1489
REAL cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32];
1490
int cxablen, cxxablen, cyablen, cyyablen, clen;
1493
REAL fin1[1152], fin2[1152];
1494
REAL *finnow, *finother, *finswap;
1497
REAL adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
1498
INEXACT REAL adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1;
1499
REAL adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0;
1500
REAL aa[4], bb[4], cc[4];
1501
INEXACT REAL aa3, bb3, cc3;
1502
INEXACT REAL ti1, tj1;
1505
INEXACT REAL u3, v3;
1506
REAL temp8[8], temp16a[16], temp16b[16], temp16c[16];
1507
REAL temp32a[32], temp32b[32], temp48[48], temp64[64];
1508
int temp8len, temp16alen, temp16blen, temp16clen;
1509
int temp32alen, temp32blen, temp48len, temp64len;
1510
REAL axtbb[8], axtcc[8], aytbb[8], aytcc[8];
1511
int axtbblen, axtcclen, aytbblen, aytcclen;
1512
REAL bxtaa[8], bxtcc[8], bytaa[8], bytcc[8];
1513
int bxtaalen, bxtcclen, bytaalen, bytcclen;
1514
REAL cxtaa[8], cxtbb[8], cytaa[8], cytbb[8];
1515
int cxtaalen, cxtbblen, cytaalen, cytbblen;
1516
REAL axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8];
1517
int axtbclen = 0, aytbclen = 0;
1518
int bxtcalen = 0, bytcalen = 0;
1519
int cxtablen = 0, cytablen = 0;
1520
REAL axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16];
1521
int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen;
1522
REAL axtbctt[8], aytbctt[8], bxtcatt[8];
1523
REAL bytcatt[8], cxtabtt[8], cytabtt[8];
1524
int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen;
1525
REAL abt[8], bct[8], cat[8];
1526
int abtlen, bctlen, catlen;
1527
REAL abtt[4], bctt[4], catt[4];
1528
int abttlen, bcttlen, cattlen;
1529
INEXACT REAL abtt3, bctt3, catt3;
1533
REAL avirt, bround, around;
1536
REAL ahi, alo, bhi, blo;
1537
REAL err1, err2, err3;
1538
INEXACT REAL _i, _j;
1541
adx = (REAL) (pa[0] - pd[0]);
1542
bdx = (REAL) (pb[0] - pd[0]);
1543
cdx = (REAL) (pc[0] - pd[0]);
1544
ady = (REAL) (pa[1] - pd[1]);
1545
bdy = (REAL) (pb[1] - pd[1]);
1546
cdy = (REAL) (pc[1] - pd[1]);
1548
Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
1549
Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
1550
Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
1552
axbclen = scale_expansion_zeroelim(4, bc, adx, axbc);
1553
axxbclen = scale_expansion_zeroelim(axbclen, axbc, adx, axxbc);
1554
aybclen = scale_expansion_zeroelim(4, bc, ady, aybc);
1555
ayybclen = scale_expansion_zeroelim(aybclen, aybc, ady, ayybc);
1556
alen = fast_expansion_sum_zeroelim(axxbclen, axxbc, ayybclen, ayybc, adet);
1558
Two_Product(cdx, ady, cdxady1, cdxady0);
1559
Two_Product(adx, cdy, adxcdy1, adxcdy0);
1560
Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
1562
bxcalen = scale_expansion_zeroelim(4, ca, bdx, bxca);
1563
bxxcalen = scale_expansion_zeroelim(bxcalen, bxca, bdx, bxxca);
1564
bycalen = scale_expansion_zeroelim(4, ca, bdy, byca);
1565
byycalen = scale_expansion_zeroelim(bycalen, byca, bdy, byyca);
1566
blen = fast_expansion_sum_zeroelim(bxxcalen, bxxca, byycalen, byyca, bdet);
1568
Two_Product(adx, bdy, adxbdy1, adxbdy0);
1569
Two_Product(bdx, ady, bdxady1, bdxady0);
1570
Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
1572
cxablen = scale_expansion_zeroelim(4, ab, cdx, cxab);
1573
cxxablen = scale_expansion_zeroelim(cxablen, cxab, cdx, cxxab);
1574
cyablen = scale_expansion_zeroelim(4, ab, cdy, cyab);
1575
cyyablen = scale_expansion_zeroelim(cyablen, cyab, cdy, cyyab);
1576
clen = fast_expansion_sum_zeroelim(cxxablen, cxxab, cyyablen, cyyab, cdet);
1578
ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
1579
finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
1581
det = estimate(finlength, fin1);
1582
errbound = iccerrboundB * permanent;
1583
if ((det >= errbound) || (-det >= errbound)) {
1587
Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
1588
Two_Diff_Tail(pa[1], pd[1], ady, adytail);
1589
Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
1590
Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
1591
Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
1592
Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
1593
if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
1594
&& (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)) {
1598
errbound = iccerrboundC * permanent + resulterrbound * Absolute(det);
1599
det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail)
1600
- (bdy * cdxtail + cdx * bdytail))
1601
+ 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx))
1602
+ ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail)
1603
- (cdy * adxtail + adx * cdytail))
1604
+ 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx))
1605
+ ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail)
1606
- (ady * bdxtail + bdx * adytail))
1607
+ 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));
1608
if ((det >= errbound) || (-det >= errbound)) {
1615
if ((bdxtail != 0.0) || (bdytail != 0.0)
1616
|| (cdxtail != 0.0) || (cdytail != 0.0)) {
1617
Square(adx, adxadx1, adxadx0);
1618
Square(ady, adyady1, adyady0);
1619
Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]);
1622
if ((cdxtail != 0.0) || (cdytail != 0.0)
1623
|| (adxtail != 0.0) || (adytail != 0.0)) {
1624
Square(bdx, bdxbdx1, bdxbdx0);
1625
Square(bdy, bdybdy1, bdybdy0);
1626
Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]);
1629
if ((adxtail != 0.0) || (adytail != 0.0)
1630
|| (bdxtail != 0.0) || (bdytail != 0.0)) {
1631
Square(cdx, cdxcdx1, cdxcdx0);
1632
Square(cdy, cdycdy1, cdycdy0);
1633
Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]);
1637
if (adxtail != 0.0) {
1638
axtbclen = scale_expansion_zeroelim(4, bc, adxtail, axtbc);
1639
temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, 2.0 * adx,
1642
axtcclen = scale_expansion_zeroelim(4, cc, adxtail, axtcc);
1643
temp16blen = scale_expansion_zeroelim(axtcclen, axtcc, bdy, temp16b);
1645
axtbblen = scale_expansion_zeroelim(4, bb, adxtail, axtbb);
1646
temp16clen = scale_expansion_zeroelim(axtbblen, axtbb, -cdy, temp16c);
1648
temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
1649
temp16blen, temp16b, temp32a);
1650
temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
1651
temp32alen, temp32a, temp48);
1652
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
1654
finswap = finnow; finnow = finother; finother = finswap;
1656
if (adytail != 0.0) {
1657
aytbclen = scale_expansion_zeroelim(4, bc, adytail, aytbc);
1658
temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, 2.0 * ady,
1661
aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb);
1662
temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b);
1664
aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc);
1665
temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c);
1667
temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
1668
temp16blen, temp16b, temp32a);
1669
temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
1670
temp32alen, temp32a, temp48);
1671
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
1673
finswap = finnow; finnow = finother; finother = finswap;
1675
if (bdxtail != 0.0) {
1676
bxtcalen = scale_expansion_zeroelim(4, ca, bdxtail, bxtca);
1677
temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, 2.0 * bdx,
1680
bxtaalen = scale_expansion_zeroelim(4, aa, bdxtail, bxtaa);
1681
temp16blen = scale_expansion_zeroelim(bxtaalen, bxtaa, cdy, temp16b);
1683
bxtcclen = scale_expansion_zeroelim(4, cc, bdxtail, bxtcc);
1684
temp16clen = scale_expansion_zeroelim(bxtcclen, bxtcc, -ady, temp16c);
1686
temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
1687
temp16blen, temp16b, temp32a);
1688
temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
1689
temp32alen, temp32a, temp48);
1690
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
1692
finswap = finnow; finnow = finother; finother = finswap;
1694
if (bdytail != 0.0) {
1695
bytcalen = scale_expansion_zeroelim(4, ca, bdytail, bytca);
1696
temp16alen = scale_expansion_zeroelim(bytcalen, bytca, 2.0 * bdy,
1699
bytcclen = scale_expansion_zeroelim(4, cc, bdytail, bytcc);
1700
temp16blen = scale_expansion_zeroelim(bytcclen, bytcc, adx, temp16b);
1702
bytaalen = scale_expansion_zeroelim(4, aa, bdytail, bytaa);
1703
temp16clen = scale_expansion_zeroelim(bytaalen, bytaa, -cdx, temp16c);
1705
temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
1706
temp16blen, temp16b, temp32a);
1707
temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
1708
temp32alen, temp32a, temp48);
1709
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
1711
finswap = finnow; finnow = finother; finother = finswap;
1713
if (cdxtail != 0.0) {
1714
cxtablen = scale_expansion_zeroelim(4, ab, cdxtail, cxtab);
1715
temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, 2.0 * cdx,
1718
cxtbblen = scale_expansion_zeroelim(4, bb, cdxtail, cxtbb);
1719
temp16blen = scale_expansion_zeroelim(cxtbblen, cxtbb, ady, temp16b);
1721
cxtaalen = scale_expansion_zeroelim(4, aa, cdxtail, cxtaa);
1722
temp16clen = scale_expansion_zeroelim(cxtaalen, cxtaa, -bdy, temp16c);
1724
temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
1725
temp16blen, temp16b, temp32a);
1726
temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
1727
temp32alen, temp32a, temp48);
1728
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
1730
finswap = finnow; finnow = finother; finother = finswap;
1732
if (cdytail != 0.0) {
1733
cytablen = scale_expansion_zeroelim(4, ab, cdytail, cytab);
1734
temp16alen = scale_expansion_zeroelim(cytablen, cytab, 2.0 * cdy,
1737
cytaalen = scale_expansion_zeroelim(4, aa, cdytail, cytaa);
1738
temp16blen = scale_expansion_zeroelim(cytaalen, cytaa, bdx, temp16b);
1740
cytbblen = scale_expansion_zeroelim(4, bb, cdytail, cytbb);
1741
temp16clen = scale_expansion_zeroelim(cytbblen, cytbb, -adx, temp16c);
1743
temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
1744
temp16blen, temp16b, temp32a);
1745
temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
1746
temp32alen, temp32a, temp48);
1747
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
1749
finswap = finnow; finnow = finother; finother = finswap;
1752
if ((adxtail != 0.0) || (adytail != 0.0)) {
1753
if ((bdxtail != 0.0) || (bdytail != 0.0)
1754
|| (cdxtail != 0.0) || (cdytail != 0.0)) {
1755
Two_Product(bdxtail, cdy, ti1, ti0);
1756
Two_Product(bdx, cdytail, tj1, tj0);
1757
Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
1760
Two_Product(cdxtail, negate, ti1, ti0);
1762
Two_Product(cdx, negate, tj1, tj0);
1763
Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
1765
bctlen = fast_expansion_sum_zeroelim(4, u, 4, v, bct);
1767
Two_Product(bdxtail, cdytail, ti1, ti0);
1768
Two_Product(cdxtail, bdytail, tj1, tj0);
1769
Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]);
1779
if (adxtail != 0.0) {
1780
temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, adxtail, temp16a);
1781
axtbctlen = scale_expansion_zeroelim(bctlen, bct, adxtail, axtbct);
1782
temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, 2.0 * adx,
1784
temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
1785
temp32alen, temp32a, temp48);
1786
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
1788
finswap = finnow; finnow = finother; finother = finswap;
1789
if (bdytail != 0.0) {
1790
temp8len = scale_expansion_zeroelim(4, cc, adxtail, temp8);
1791
temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
1793
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
1795
finswap = finnow; finnow = finother; finother = finswap;
1797
if (cdytail != 0.0) {
1798
temp8len = scale_expansion_zeroelim(4, bb, -adxtail, temp8);
1799
temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
1801
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
1803
finswap = finnow; finnow = finother; finother = finswap;
1806
temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, adxtail,
1808
axtbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adxtail, axtbctt);
1809
temp16alen = scale_expansion_zeroelim(axtbcttlen, axtbctt, 2.0 * adx,
1811
temp16blen = scale_expansion_zeroelim(axtbcttlen, axtbctt, adxtail,
1813
temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
1814
temp16blen, temp16b, temp32b);
1815
temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
1816
temp32blen, temp32b, temp64);
1817
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
1819
finswap = finnow; finnow = finother; finother = finswap;
1821
if (adytail != 0.0) {
1822
temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, adytail, temp16a);
1823
aytbctlen = scale_expansion_zeroelim(bctlen, bct, adytail, aytbct);
1824
temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, 2.0 * ady,
1826
temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
1827
temp32alen, temp32a, temp48);
1828
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
1830
finswap = finnow; finnow = finother; finother = finswap;
1833
temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, adytail,
1835
aytbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adytail, aytbctt);
1836
temp16alen = scale_expansion_zeroelim(aytbcttlen, aytbctt, 2.0 * ady,
1838
temp16blen = scale_expansion_zeroelim(aytbcttlen, aytbctt, adytail,
1840
temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
1841
temp16blen, temp16b, temp32b);
1842
temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
1843
temp32blen, temp32b, temp64);
1844
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
1846
finswap = finnow; finnow = finother; finother = finswap;
1849
if ((bdxtail != 0.0) || (bdytail != 0.0)) {
1850
if ((cdxtail != 0.0) || (cdytail != 0.0)
1851
|| (adxtail != 0.0) || (adytail != 0.0)) {
1852
Two_Product(cdxtail, ady, ti1, ti0);
1853
Two_Product(cdx, adytail, tj1, tj0);
1854
Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
1857
Two_Product(adxtail, negate, ti1, ti0);
1859
Two_Product(adx, negate, tj1, tj0);
1860
Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
1862
catlen = fast_expansion_sum_zeroelim(4, u, 4, v, cat);
1864
Two_Product(cdxtail, adytail, ti1, ti0);
1865
Two_Product(adxtail, cdytail, tj1, tj0);
1866
Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]);
1876
if (bdxtail != 0.0) {
1877
temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, bdxtail, temp16a);
1878
bxtcatlen = scale_expansion_zeroelim(catlen, cat, bdxtail, bxtcat);
1879
temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, 2.0 * bdx,
1881
temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
1882
temp32alen, temp32a, temp48);
1883
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
1885
finswap = finnow; finnow = finother; finother = finswap;
1886
if (cdytail != 0.0) {
1887
temp8len = scale_expansion_zeroelim(4, aa, bdxtail, temp8);
1888
temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
1890
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
1892
finswap = finnow; finnow = finother; finother = finswap;
1894
if (adytail != 0.0) {
1895
temp8len = scale_expansion_zeroelim(4, cc, -bdxtail, temp8);
1896
temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
1898
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
1900
finswap = finnow; finnow = finother; finother = finswap;
1903
temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, bdxtail,
1905
bxtcattlen = scale_expansion_zeroelim(cattlen, catt, bdxtail, bxtcatt);
1906
temp16alen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, 2.0 * bdx,
1908
temp16blen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, bdxtail,
1910
temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
1911
temp16blen, temp16b, temp32b);
1912
temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
1913
temp32blen, temp32b, temp64);
1914
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
1916
finswap = finnow; finnow = finother; finother = finswap;
1918
if (bdytail != 0.0) {
1919
temp16alen = scale_expansion_zeroelim(bytcalen, bytca, bdytail, temp16a);
1920
bytcatlen = scale_expansion_zeroelim(catlen, cat, bdytail, bytcat);
1921
temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, 2.0 * bdy,
1923
temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
1924
temp32alen, temp32a, temp48);
1925
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
1927
finswap = finnow; finnow = finother; finother = finswap;
1930
temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, bdytail,
1932
bytcattlen = scale_expansion_zeroelim(cattlen, catt, bdytail, bytcatt);
1933
temp16alen = scale_expansion_zeroelim(bytcattlen, bytcatt, 2.0 * bdy,
1935
temp16blen = scale_expansion_zeroelim(bytcattlen, bytcatt, bdytail,
1937
temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
1938
temp16blen, temp16b, temp32b);
1939
temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
1940
temp32blen, temp32b, temp64);
1941
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
1943
finswap = finnow; finnow = finother; finother = finswap;
1946
if ((cdxtail != 0.0) || (cdytail != 0.0)) {
1947
if ((adxtail != 0.0) || (adytail != 0.0)
1948
|| (bdxtail != 0.0) || (bdytail != 0.0)) {
1949
Two_Product(adxtail, bdy, ti1, ti0);
1950
Two_Product(adx, bdytail, tj1, tj0);
1951
Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
1954
Two_Product(bdxtail, negate, ti1, ti0);
1956
Two_Product(bdx, negate, tj1, tj0);
1957
Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
1959
abtlen = fast_expansion_sum_zeroelim(4, u, 4, v, abt);
1961
Two_Product(adxtail, bdytail, ti1, ti0);
1962
Two_Product(bdxtail, adytail, tj1, tj0);
1963
Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]);
1973
if (cdxtail != 0.0) {
1974
temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, cdxtail, temp16a);
1975
cxtabtlen = scale_expansion_zeroelim(abtlen, abt, cdxtail, cxtabt);
1976
temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, 2.0 * cdx,
1978
temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
1979
temp32alen, temp32a, temp48);
1980
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
1982
finswap = finnow; finnow = finother; finother = finswap;
1983
if (adytail != 0.0) {
1984
temp8len = scale_expansion_zeroelim(4, bb, cdxtail, temp8);
1985
temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
1987
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
1989
finswap = finnow; finnow = finother; finother = finswap;
1991
if (bdytail != 0.0) {
1992
temp8len = scale_expansion_zeroelim(4, aa, -cdxtail, temp8);
1993
temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
1995
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
1997
finswap = finnow; finnow = finother; finother = finswap;
2000
temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, cdxtail,
2002
cxtabttlen = scale_expansion_zeroelim(abttlen, abtt, cdxtail, cxtabtt);
2003
temp16alen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, 2.0 * cdx,
2005
temp16blen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, cdxtail,
2007
temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2008
temp16blen, temp16b, temp32b);
2009
temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
2010
temp32blen, temp32b, temp64);
2011
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
2013
finswap = finnow; finnow = finother; finother = finswap;
2015
if (cdytail != 0.0) {
2016
temp16alen = scale_expansion_zeroelim(cytablen, cytab, cdytail, temp16a);
2017
cytabtlen = scale_expansion_zeroelim(abtlen, abt, cdytail, cytabt);
2018
temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, 2.0 * cdy,
2020
temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2021
temp32alen, temp32a, temp48);
2022
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
2024
finswap = finnow; finnow = finother; finother = finswap;
2027
temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, cdytail,
2029
cytabttlen = scale_expansion_zeroelim(abttlen, abtt, cdytail, cytabtt);
2030
temp16alen = scale_expansion_zeroelim(cytabttlen, cytabtt, 2.0 * cdy,
2032
temp16blen = scale_expansion_zeroelim(cytabttlen, cytabtt, cdytail,
2034
temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
2035
temp16blen, temp16b, temp32b);
2036
temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
2037
temp32blen, temp32b, temp64);
2038
finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
2040
finswap = finnow; finnow = finother; finother = finswap;
2044
return finnow[finlength - 1];
2047
REAL incircle(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
2049
REAL adx, bdx, cdx, ady, bdy, cdy;
2050
REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
2051
REAL alift, blift, clift;
2053
REAL permanent, errbound;
2058
adx = pa[0] - pd[0];
2059
bdx = pb[0] - pd[0];
2060
cdx = pc[0] - pd[0];
2061
ady = pa[1] - pd[1];
2062
bdy = pb[1] - pd[1];
2063
cdy = pc[1] - pd[1];
2067
alift = adx * adx + ady * ady;
2071
blift = bdx * bdx + bdy * bdy;
2075
clift = cdx * cdx + cdy * cdy;
2077
det = alift * (bdxcdy - cdxbdy)
2078
+ blift * (cdxady - adxcdy)
2079
+ clift * (adxbdy - bdxady);
2081
permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift
2082
+ (Absolute(cdxady) + Absolute(adxcdy)) * blift
2083
+ (Absolute(adxbdy) + Absolute(bdxady)) * clift;
2084
errbound = iccerrboundA * permanent;
2085
if ((det > errbound) || (-det > errbound)) {
2090
inc = incircleadapt(pa, pb, pc, pd, permanent);
2095
/*****************************************************************************/
2097
/* inspherefast() Approximate 3D insphere test. Nonrobust. */
2098
/* insphereexact() Exact 3D insphere test. Robust. */
2099
/* insphereslow() Another exact 3D insphere test. Robust. */
2100
/* insphere() Adaptive exact 3D insphere test. Robust. */
2102
/* Return a positive value if the point pe lies inside the */
2103
/* sphere passing through pa, pb, pc, and pd; a negative value */
2104
/* if it lies outside; and zero if the five points are */
2105
/* cospherical. The points pa, pb, pc, and pd must be ordered */
2106
/* so that they have a positive orientation (as defined by */
2107
/* orient3d()), or the sign of the result will be reversed. */
2109
/* Only the first and last routine should be used; the middle two are for */
2112
/* The last three use exact arithmetic to ensure a correct answer. The */
2113
/* result returned is the determinant of a matrix. In insphere() only, */
2114
/* this determinant is computed adaptively, in the sense that exact */
2115
/* arithmetic is used only to the degree it is needed to ensure that the */
2116
/* returned value has the correct sign. Hence, insphere() is usually quite */
2117
/* fast, but will run more slowly when the input points are cospherical or */
2120
/*****************************************************************************/
2122
static REAL insphereexact(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
2124
INEXACT REAL axby1, bxcy1, cxdy1, dxey1, exay1;
2125
INEXACT REAL bxay1, cxby1, dxcy1, exdy1, axey1;
2126
INEXACT REAL axcy1, bxdy1, cxey1, dxay1, exby1;
2127
INEXACT REAL cxay1, dxby1, excy1, axdy1, bxey1;
2128
REAL axby0, bxcy0, cxdy0, dxey0, exay0;
2129
REAL bxay0, cxby0, dxcy0, exdy0, axey0;
2130
REAL axcy0, bxdy0, cxey0, dxay0, exby0;
2131
REAL cxay0, dxby0, excy0, axdy0, bxey0;
2132
REAL ab[4], bc[4], cd[4], de[4], ea[4];
2133
REAL ac[4], bd[4], ce[4], da[4], eb[4];
2134
REAL temp8a[8], temp8b[8], temp16[16];
2135
int temp8alen, temp8blen, temp16len;
2136
REAL abc[24], bcd[24], cde[24], dea[24], eab[24];
2137
REAL abd[24], bce[24], cda[24], deb[24], eac[24];
2138
int abclen, bcdlen, cdelen, dealen, eablen;
2139
int abdlen, bcelen, cdalen, deblen, eaclen;
2140
REAL temp48a[48], temp48b[48];
2141
int temp48alen, temp48blen;
2142
REAL abcd[96], bcde[96], cdea[96], deab[96], eabc[96];
2143
int abcdlen, bcdelen, cdealen, deablen, eabclen;
2145
REAL det384x[384], det384y[384], det384z[384];
2146
int xlen, ylen, zlen;
2149
REAL adet[1152], bdet[1152], cdet[1152], ddet[1152], edet[1152];
2150
int alen, blen, clen, dlen, elen;
2151
REAL abdet[2304], cddet[2304], cdedet[3456];
2158
REAL avirt, bround, around;
2161
REAL ahi, alo, bhi, blo;
2162
REAL err1, err2, err3;
2163
INEXACT REAL _i, _j;
2166
Two_Product(pa[0], pb[1], axby1, axby0);
2167
Two_Product(pb[0], pa[1], bxay1, bxay0);
2168
Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]);
2170
Two_Product(pb[0], pc[1], bxcy1, bxcy0);
2171
Two_Product(pc[0], pb[1], cxby1, cxby0);
2172
Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]);
2174
Two_Product(pc[0], pd[1], cxdy1, cxdy0);
2175
Two_Product(pd[0], pc[1], dxcy1, dxcy0);
2176
Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]);
2178
Two_Product(pd[0], pe[1], dxey1, dxey0);
2179
Two_Product(pe[0], pd[1], exdy1, exdy0);
2180
Two_Two_Diff(dxey1, dxey0, exdy1, exdy0, de[3], de[2], de[1], de[0]);
2182
Two_Product(pe[0], pa[1], exay1, exay0);
2183
Two_Product(pa[0], pe[1], axey1, axey0);
2184
Two_Two_Diff(exay1, exay0, axey1, axey0, ea[3], ea[2], ea[1], ea[0]);
2186
Two_Product(pa[0], pc[1], axcy1, axcy0);
2187
Two_Product(pc[0], pa[1], cxay1, cxay0);
2188
Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]);
2190
Two_Product(pb[0], pd[1], bxdy1, bxdy0);
2191
Two_Product(pd[0], pb[1], dxby1, dxby0);
2192
Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]);
2194
Two_Product(pc[0], pe[1], cxey1, cxey0);
2195
Two_Product(pe[0], pc[1], excy1, excy0);
2196
Two_Two_Diff(cxey1, cxey0, excy1, excy0, ce[3], ce[2], ce[1], ce[0]);
2198
Two_Product(pd[0], pa[1], dxay1, dxay0);
2199
Two_Product(pa[0], pd[1], axdy1, axdy0);
2200
Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]);
2202
Two_Product(pe[0], pb[1], exby1, exby0);
2203
Two_Product(pb[0], pe[1], bxey1, bxey0);
2204
Two_Two_Diff(exby1, exby0, bxey1, bxey0, eb[3], eb[2], eb[1], eb[0]);
2206
temp8alen = scale_expansion_zeroelim(4, bc, pa[2], temp8a);
2207
temp8blen = scale_expansion_zeroelim(4, ac, -pb[2], temp8b);
2208
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
2210
temp8alen = scale_expansion_zeroelim(4, ab, pc[2], temp8a);
2211
abclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
2214
temp8alen = scale_expansion_zeroelim(4, cd, pb[2], temp8a);
2215
temp8blen = scale_expansion_zeroelim(4, bd, -pc[2], temp8b);
2216
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
2218
temp8alen = scale_expansion_zeroelim(4, bc, pd[2], temp8a);
2219
bcdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
2222
temp8alen = scale_expansion_zeroelim(4, de, pc[2], temp8a);
2223
temp8blen = scale_expansion_zeroelim(4, ce, -pd[2], temp8b);
2224
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
2226
temp8alen = scale_expansion_zeroelim(4, cd, pe[2], temp8a);
2227
cdelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
2230
temp8alen = scale_expansion_zeroelim(4, ea, pd[2], temp8a);
2231
temp8blen = scale_expansion_zeroelim(4, da, -pe[2], temp8b);
2232
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
2234
temp8alen = scale_expansion_zeroelim(4, de, pa[2], temp8a);
2235
dealen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
2238
temp8alen = scale_expansion_zeroelim(4, ab, pe[2], temp8a);
2239
temp8blen = scale_expansion_zeroelim(4, eb, -pa[2], temp8b);
2240
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
2242
temp8alen = scale_expansion_zeroelim(4, ea, pb[2], temp8a);
2243
eablen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
2246
temp8alen = scale_expansion_zeroelim(4, bd, pa[2], temp8a);
2247
temp8blen = scale_expansion_zeroelim(4, da, pb[2], temp8b);
2248
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
2250
temp8alen = scale_expansion_zeroelim(4, ab, pd[2], temp8a);
2251
abdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
2254
temp8alen = scale_expansion_zeroelim(4, ce, pb[2], temp8a);
2255
temp8blen = scale_expansion_zeroelim(4, eb, pc[2], temp8b);
2256
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
2258
temp8alen = scale_expansion_zeroelim(4, bc, pe[2], temp8a);
2259
bcelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
2262
temp8alen = scale_expansion_zeroelim(4, da, pc[2], temp8a);
2263
temp8blen = scale_expansion_zeroelim(4, ac, pd[2], temp8b);
2264
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
2266
temp8alen = scale_expansion_zeroelim(4, cd, pa[2], temp8a);
2267
cdalen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
2270
temp8alen = scale_expansion_zeroelim(4, eb, pd[2], temp8a);
2271
temp8blen = scale_expansion_zeroelim(4, bd, pe[2], temp8b);
2272
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
2274
temp8alen = scale_expansion_zeroelim(4, de, pb[2], temp8a);
2275
deblen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
2278
temp8alen = scale_expansion_zeroelim(4, ac, pe[2], temp8a);
2279
temp8blen = scale_expansion_zeroelim(4, ce, pa[2], temp8b);
2280
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
2282
temp8alen = scale_expansion_zeroelim(4, ea, pc[2], temp8a);
2283
eaclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
2286
temp48alen = fast_expansion_sum_zeroelim(cdelen, cde, bcelen, bce, temp48a);
2287
temp48blen = fast_expansion_sum_zeroelim(deblen, deb, bcdlen, bcd, temp48b);
2288
for (i = 0; i < temp48blen; i++) {
2289
temp48b[i] = -temp48b[i];
2291
bcdelen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
2292
temp48blen, temp48b, bcde);
2293
xlen = scale_expansion_zeroelim(bcdelen, bcde, pa[0], temp192);
2294
xlen = scale_expansion_zeroelim(xlen, temp192, pa[0], det384x);
2295
ylen = scale_expansion_zeroelim(bcdelen, bcde, pa[1], temp192);
2296
ylen = scale_expansion_zeroelim(ylen, temp192, pa[1], det384y);
2297
zlen = scale_expansion_zeroelim(bcdelen, bcde, pa[2], temp192);
2298
zlen = scale_expansion_zeroelim(zlen, temp192, pa[2], det384z);
2299
xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
2300
alen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, adet);
2302
temp48alen = fast_expansion_sum_zeroelim(dealen, dea, cdalen, cda, temp48a);
2303
temp48blen = fast_expansion_sum_zeroelim(eaclen, eac, cdelen, cde, temp48b);
2304
for (i = 0; i < temp48blen; i++) {
2305
temp48b[i] = -temp48b[i];
2307
cdealen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
2308
temp48blen, temp48b, cdea);
2309
xlen = scale_expansion_zeroelim(cdealen, cdea, pb[0], temp192);
2310
xlen = scale_expansion_zeroelim(xlen, temp192, pb[0], det384x);
2311
ylen = scale_expansion_zeroelim(cdealen, cdea, pb[1], temp192);
2312
ylen = scale_expansion_zeroelim(ylen, temp192, pb[1], det384y);
2313
zlen = scale_expansion_zeroelim(cdealen, cdea, pb[2], temp192);
2314
zlen = scale_expansion_zeroelim(zlen, temp192, pb[2], det384z);
2315
xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
2316
blen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, bdet);
2318
temp48alen = fast_expansion_sum_zeroelim(eablen, eab, deblen, deb, temp48a);
2319
temp48blen = fast_expansion_sum_zeroelim(abdlen, abd, dealen, dea, temp48b);
2320
for (i = 0; i < temp48blen; i++) {
2321
temp48b[i] = -temp48b[i];
2323
deablen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
2324
temp48blen, temp48b, deab);
2325
xlen = scale_expansion_zeroelim(deablen, deab, pc[0], temp192);
2326
xlen = scale_expansion_zeroelim(xlen, temp192, pc[0], det384x);
2327
ylen = scale_expansion_zeroelim(deablen, deab, pc[1], temp192);
2328
ylen = scale_expansion_zeroelim(ylen, temp192, pc[1], det384y);
2329
zlen = scale_expansion_zeroelim(deablen, deab, pc[2], temp192);
2330
zlen = scale_expansion_zeroelim(zlen, temp192, pc[2], det384z);
2331
xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
2332
clen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, cdet);
2334
temp48alen = fast_expansion_sum_zeroelim(abclen, abc, eaclen, eac, temp48a);
2335
temp48blen = fast_expansion_sum_zeroelim(bcelen, bce, eablen, eab, temp48b);
2336
for (i = 0; i < temp48blen; i++) {
2337
temp48b[i] = -temp48b[i];
2339
eabclen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
2340
temp48blen, temp48b, eabc);
2341
xlen = scale_expansion_zeroelim(eabclen, eabc, pd[0], temp192);
2342
xlen = scale_expansion_zeroelim(xlen, temp192, pd[0], det384x);
2343
ylen = scale_expansion_zeroelim(eabclen, eabc, pd[1], temp192);
2344
ylen = scale_expansion_zeroelim(ylen, temp192, pd[1], det384y);
2345
zlen = scale_expansion_zeroelim(eabclen, eabc, pd[2], temp192);
2346
zlen = scale_expansion_zeroelim(zlen, temp192, pd[2], det384z);
2347
xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
2348
dlen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, ddet);
2350
temp48alen = fast_expansion_sum_zeroelim(bcdlen, bcd, abdlen, abd, temp48a);
2351
temp48blen = fast_expansion_sum_zeroelim(cdalen, cda, abclen, abc, temp48b);
2352
for (i = 0; i < temp48blen; i++) {
2353
temp48b[i] = -temp48b[i];
2355
abcdlen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
2356
temp48blen, temp48b, abcd);
2357
xlen = scale_expansion_zeroelim(abcdlen, abcd, pe[0], temp192);
2358
xlen = scale_expansion_zeroelim(xlen, temp192, pe[0], det384x);
2359
ylen = scale_expansion_zeroelim(abcdlen, abcd, pe[1], temp192);
2360
ylen = scale_expansion_zeroelim(ylen, temp192, pe[1], det384y);
2361
zlen = scale_expansion_zeroelim(abcdlen, abcd, pe[2], temp192);
2362
zlen = scale_expansion_zeroelim(zlen, temp192, pe[2], det384z);
2363
xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
2364
elen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, edet);
2366
ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
2367
cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
2368
cdelen = fast_expansion_sum_zeroelim(cdlen, cddet, elen, edet, cdedet);
2369
deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdelen, cdedet, deter);
2371
return deter[deterlen - 1];
2374
static REAL insphereadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe,
2377
INEXACT REAL aex, bex, cex, dex, aey, bey, cey, dey, aez, bez, cez, dez;
2380
INEXACT REAL aexbey1, bexaey1, bexcey1, cexbey1;
2381
INEXACT REAL cexdey1, dexcey1, dexaey1, aexdey1;
2382
INEXACT REAL aexcey1, cexaey1, bexdey1, dexbey1;
2383
REAL aexbey0, bexaey0, bexcey0, cexbey0;
2384
REAL cexdey0, dexcey0, dexaey0, aexdey0;
2385
REAL aexcey0, cexaey0, bexdey0, dexbey0;
2386
REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4];
2387
INEXACT REAL ab3, bc3, cd3, da3, ac3, bd3;
2388
REAL abeps, bceps, cdeps, daeps, aceps, bdeps;
2389
REAL temp8a[8], temp8b[8], temp8c[8], temp16[16], temp24[24], temp48[48];
2390
int temp8alen, temp8blen, temp8clen, temp16len, temp24len, temp48len;
2391
REAL xdet[96], ydet[96], zdet[96], xydet[192];
2392
int xlen, ylen, zlen, xylen;
2393
REAL adet[288], bdet[288], cdet[288], ddet[288];
2394
int alen, blen, clen, dlen;
2395
REAL abdet[576], cddet[576];
2400
REAL aextail, bextail, cextail, dextail;
2401
REAL aeytail, beytail, ceytail, deytail;
2402
REAL aeztail, beztail, ceztail, deztail;
2405
REAL avirt, bround, around;
2408
REAL ahi, alo, bhi, blo;
2409
REAL err1, err2, err3;
2410
INEXACT REAL _i, _j;
2413
aex = (REAL) (pa[0] - pe[0]);
2414
bex = (REAL) (pb[0] - pe[0]);
2415
cex = (REAL) (pc[0] - pe[0]);
2416
dex = (REAL) (pd[0] - pe[0]);
2417
aey = (REAL) (pa[1] - pe[1]);
2418
bey = (REAL) (pb[1] - pe[1]);
2419
cey = (REAL) (pc[1] - pe[1]);
2420
dey = (REAL) (pd[1] - pe[1]);
2421
aez = (REAL) (pa[2] - pe[2]);
2422
bez = (REAL) (pb[2] - pe[2]);
2423
cez = (REAL) (pc[2] - pe[2]);
2424
dez = (REAL) (pd[2] - pe[2]);
2426
Two_Product(aex, bey, aexbey1, aexbey0);
2427
Two_Product(bex, aey, bexaey1, bexaey0);
2428
Two_Two_Diff(aexbey1, aexbey0, bexaey1, bexaey0, ab3, ab[2], ab[1], ab[0]);
2431
Two_Product(bex, cey, bexcey1, bexcey0);
2432
Two_Product(cex, bey, cexbey1, cexbey0);
2433
Two_Two_Diff(bexcey1, bexcey0, cexbey1, cexbey0, bc3, bc[2], bc[1], bc[0]);
2436
Two_Product(cex, dey, cexdey1, cexdey0);
2437
Two_Product(dex, cey, dexcey1, dexcey0);
2438
Two_Two_Diff(cexdey1, cexdey0, dexcey1, dexcey0, cd3, cd[2], cd[1], cd[0]);
2441
Two_Product(dex, aey, dexaey1, dexaey0);
2442
Two_Product(aex, dey, aexdey1, aexdey0);
2443
Two_Two_Diff(dexaey1, dexaey0, aexdey1, aexdey0, da3, da[2], da[1], da[0]);
2446
Two_Product(aex, cey, aexcey1, aexcey0);
2447
Two_Product(cex, aey, cexaey1, cexaey0);
2448
Two_Two_Diff(aexcey1, aexcey0, cexaey1, cexaey0, ac3, ac[2], ac[1], ac[0]);
2451
Two_Product(bex, dey, bexdey1, bexdey0);
2452
Two_Product(dex, bey, dexbey1, dexbey0);
2453
Two_Two_Diff(bexdey1, bexdey0, dexbey1, dexbey0, bd3, bd[2], bd[1], bd[0]);
2456
temp8alen = scale_expansion_zeroelim(4, cd, bez, temp8a);
2457
temp8blen = scale_expansion_zeroelim(4, bd, -cez, temp8b);
2458
temp8clen = scale_expansion_zeroelim(4, bc, dez, temp8c);
2459
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
2460
temp8blen, temp8b, temp16);
2461
temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
2462
temp16len, temp16, temp24);
2463
temp48len = scale_expansion_zeroelim(temp24len, temp24, aex, temp48);
2464
xlen = scale_expansion_zeroelim(temp48len, temp48, -aex, xdet);
2465
temp48len = scale_expansion_zeroelim(temp24len, temp24, aey, temp48);
2466
ylen = scale_expansion_zeroelim(temp48len, temp48, -aey, ydet);
2467
temp48len = scale_expansion_zeroelim(temp24len, temp24, aez, temp48);
2468
zlen = scale_expansion_zeroelim(temp48len, temp48, -aez, zdet);
2469
xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
2470
alen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, adet);
2472
temp8alen = scale_expansion_zeroelim(4, da, cez, temp8a);
2473
temp8blen = scale_expansion_zeroelim(4, ac, dez, temp8b);
2474
temp8clen = scale_expansion_zeroelim(4, cd, aez, temp8c);
2475
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
2476
temp8blen, temp8b, temp16);
2477
temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
2478
temp16len, temp16, temp24);
2479
temp48len = scale_expansion_zeroelim(temp24len, temp24, bex, temp48);
2480
xlen = scale_expansion_zeroelim(temp48len, temp48, bex, xdet);
2481
temp48len = scale_expansion_zeroelim(temp24len, temp24, bey, temp48);
2482
ylen = scale_expansion_zeroelim(temp48len, temp48, bey, ydet);
2483
temp48len = scale_expansion_zeroelim(temp24len, temp24, bez, temp48);
2484
zlen = scale_expansion_zeroelim(temp48len, temp48, bez, zdet);
2485
xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
2486
blen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, bdet);
2488
temp8alen = scale_expansion_zeroelim(4, ab, dez, temp8a);
2489
temp8blen = scale_expansion_zeroelim(4, bd, aez, temp8b);
2490
temp8clen = scale_expansion_zeroelim(4, da, bez, temp8c);
2491
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
2492
temp8blen, temp8b, temp16);
2493
temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
2494
temp16len, temp16, temp24);
2495
temp48len = scale_expansion_zeroelim(temp24len, temp24, cex, temp48);
2496
xlen = scale_expansion_zeroelim(temp48len, temp48, -cex, xdet);
2497
temp48len = scale_expansion_zeroelim(temp24len, temp24, cey, temp48);
2498
ylen = scale_expansion_zeroelim(temp48len, temp48, -cey, ydet);
2499
temp48len = scale_expansion_zeroelim(temp24len, temp24, cez, temp48);
2500
zlen = scale_expansion_zeroelim(temp48len, temp48, -cez, zdet);
2501
xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
2502
clen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, cdet);
2504
temp8alen = scale_expansion_zeroelim(4, bc, aez, temp8a);
2505
temp8blen = scale_expansion_zeroelim(4, ac, -bez, temp8b);
2506
temp8clen = scale_expansion_zeroelim(4, ab, cez, temp8c);
2507
temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
2508
temp8blen, temp8b, temp16);
2509
temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
2510
temp16len, temp16, temp24);
2511
temp48len = scale_expansion_zeroelim(temp24len, temp24, dex, temp48);
2512
xlen = scale_expansion_zeroelim(temp48len, temp48, dex, xdet);
2513
temp48len = scale_expansion_zeroelim(temp24len, temp24, dey, temp48);
2514
ylen = scale_expansion_zeroelim(temp48len, temp48, dey, ydet);
2515
temp48len = scale_expansion_zeroelim(temp24len, temp24, dez, temp48);
2516
zlen = scale_expansion_zeroelim(temp48len, temp48, dez, zdet);
2517
xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
2518
dlen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, ddet);
2520
ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
2521
cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
2522
finlength = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, fin1);
2524
det = estimate(finlength, fin1);
2525
errbound = isperrboundB * permanent;
2526
if ((det >= errbound) || (-det >= errbound)) {
2530
Two_Diff_Tail(pa[0], pe[0], aex, aextail);
2531
Two_Diff_Tail(pa[1], pe[1], aey, aeytail);
2532
Two_Diff_Tail(pa[2], pe[2], aez, aeztail);
2533
Two_Diff_Tail(pb[0], pe[0], bex, bextail);
2534
Two_Diff_Tail(pb[1], pe[1], bey, beytail);
2535
Two_Diff_Tail(pb[2], pe[2], bez, beztail);
2536
Two_Diff_Tail(pc[0], pe[0], cex, cextail);
2537
Two_Diff_Tail(pc[1], pe[1], cey, ceytail);
2538
Two_Diff_Tail(pc[2], pe[2], cez, ceztail);
2539
Two_Diff_Tail(pd[0], pe[0], dex, dextail);
2540
Two_Diff_Tail(pd[1], pe[1], dey, deytail);
2541
Two_Diff_Tail(pd[2], pe[2], dez, deztail);
2542
if ((aextail == 0.0) && (aeytail == 0.0) && (aeztail == 0.0)
2543
&& (bextail == 0.0) && (beytail == 0.0) && (beztail == 0.0)
2544
&& (cextail == 0.0) && (ceytail == 0.0) && (ceztail == 0.0)
2545
&& (dextail == 0.0) && (deytail == 0.0) && (deztail == 0.0)) {
2549
errbound = isperrboundC * permanent + resulterrbound * Absolute(det);
2550
abeps = (aex * beytail + bey * aextail)
2551
- (aey * bextail + bex * aeytail);
2552
bceps = (bex * ceytail + cey * bextail)
2553
- (bey * cextail + cex * beytail);
2554
cdeps = (cex * deytail + dey * cextail)
2555
- (cey * dextail + dex * ceytail);
2556
daeps = (dex * aeytail + aey * dextail)
2557
- (dey * aextail + aex * deytail);
2558
aceps = (aex * ceytail + cey * aextail)
2559
- (aey * cextail + cex * aeytail);
2560
bdeps = (bex * deytail + dey * bextail)
2561
- (bey * dextail + dex * beytail);
2562
det += (((bex * bex + bey * bey + bez * bez)
2563
* ((cez * daeps + dez * aceps + aez * cdeps)
2564
+ (ceztail * da3 + deztail * ac3 + aeztail * cd3))
2565
+ (dex * dex + dey * dey + dez * dez)
2566
* ((aez * bceps - bez * aceps + cez * abeps)
2567
+ (aeztail * bc3 - beztail * ac3 + ceztail * ab3)))
2568
- ((aex * aex + aey * aey + aez * aez)
2569
* ((bez * cdeps - cez * bdeps + dez * bceps)
2570
+ (beztail * cd3 - ceztail * bd3 + deztail * bc3))
2571
+ (cex * cex + cey * cey + cez * cez)
2572
* ((dez * abeps + aez * bdeps + bez * daeps)
2573
+ (deztail * ab3 + aeztail * bd3 + beztail * da3))))
2574
+ 2.0 * (((bex * bextail + bey * beytail + bez * beztail)
2575
* (cez * da3 + dez * ac3 + aez * cd3)
2576
+ (dex * dextail + dey * deytail + dez * deztail)
2577
* (aez * bc3 - bez * ac3 + cez * ab3))
2578
- ((aex * aextail + aey * aeytail + aez * aeztail)
2579
* (bez * cd3 - cez * bd3 + dez * bc3)
2580
+ (cex * cextail + cey * ceytail + cez * ceztail)
2581
* (dez * ab3 + aez * bd3 + bez * da3)));
2582
if ((det >= errbound) || (-det >= errbound)) {
2586
return insphereexact(pa, pb, pc, pd, pe);
2589
REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
2591
REAL aex, bex, cex, dex;
2592
REAL aey, bey, cey, dey;
2593
REAL aez, bez, cez, dez;
2594
REAL aexbey, bexaey, bexcey, cexbey, cexdey, dexcey, dexaey, aexdey;
2595
REAL aexcey, cexaey, bexdey, dexbey;
2596
REAL alift, blift, clift, dlift;
2597
REAL ab, bc, cd, da, ac, bd;
2598
REAL abc, bcd, cda, dab;
2599
REAL aezplus, bezplus, cezplus, dezplus;
2600
REAL aexbeyplus, bexaeyplus, bexceyplus, cexbeyplus;
2601
REAL cexdeyplus, dexceyplus, dexaeyplus, aexdeyplus;
2602
REAL aexceyplus, cexaeyplus, bexdeyplus, dexbeyplus;
2604
REAL permanent, errbound;
2609
aex = pa[0] - pe[0];
2610
bex = pb[0] - pe[0];
2611
cex = pc[0] - pe[0];
2612
dex = pd[0] - pe[0];
2613
aey = pa[1] - pe[1];
2614
bey = pb[1] - pe[1];
2615
cey = pc[1] - pe[1];
2616
dey = pd[1] - pe[1];
2617
aez = pa[2] - pe[2];
2618
bez = pb[2] - pe[2];
2619
cez = pc[2] - pe[2];
2620
dez = pd[2] - pe[2];
2624
ab = aexbey - bexaey;
2627
bc = bexcey - cexbey;
2630
cd = cexdey - dexcey;
2633
da = dexaey - aexdey;
2637
ac = aexcey - cexaey;
2640
bd = bexdey - dexbey;
2642
abc = aez * bc - bez * ac + cez * ab;
2643
bcd = bez * cd - cez * bd + dez * bc;
2644
cda = cez * da + dez * ac + aez * cd;
2645
dab = dez * ab + aez * bd + bez * da;
2647
alift = aex * aex + aey * aey + aez * aez;
2648
blift = bex * bex + bey * bey + bez * bez;
2649
clift = cex * cex + cey * cey + cez * cez;
2650
dlift = dex * dex + dey * dey + dez * dez;
2652
det = (dlift * abc - clift * dab) + (blift * cda - alift * bcd);
2654
aezplus = Absolute(aez);
2655
bezplus = Absolute(bez);
2656
cezplus = Absolute(cez);
2657
dezplus = Absolute(dez);
2658
aexbeyplus = Absolute(aexbey);
2659
bexaeyplus = Absolute(bexaey);
2660
bexceyplus = Absolute(bexcey);
2661
cexbeyplus = Absolute(cexbey);
2662
cexdeyplus = Absolute(cexdey);
2663
dexceyplus = Absolute(dexcey);
2664
dexaeyplus = Absolute(dexaey);
2665
aexdeyplus = Absolute(aexdey);
2666
aexceyplus = Absolute(aexcey);
2667
cexaeyplus = Absolute(cexaey);
2668
bexdeyplus = Absolute(bexdey);
2669
dexbeyplus = Absolute(dexbey);
2670
permanent = ((cexdeyplus + dexceyplus) * bezplus
2671
+ (dexbeyplus + bexdeyplus) * cezplus
2672
+ (bexceyplus + cexbeyplus) * dezplus)
2674
+ ((dexaeyplus + aexdeyplus) * cezplus
2675
+ (aexceyplus + cexaeyplus) * dezplus
2676
+ (cexdeyplus + dexceyplus) * aezplus)
2678
+ ((aexbeyplus + bexaeyplus) * dezplus
2679
+ (bexdeyplus + dexbeyplus) * aezplus
2680
+ (dexaeyplus + aexdeyplus) * bezplus)
2682
+ ((bexceyplus + cexbeyplus) * aezplus
2683
+ (cexaeyplus + aexceyplus) * bezplus
2684
+ (aexbeyplus + bexaeyplus) * cezplus)
2686
errbound = isperrboundA * permanent;
2687
if ((det > errbound) || (-det > errbound)) {
2692
ins = insphereadapt(pa, pb, pc, pd, pe, permanent);
added
removed
src/kernel/mesh/GeometricPredicates.cpp
