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