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