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