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