Chaste Release::3.1
predicates.cpp
00001 /*****************************************************************************/
00002 /*                                                                           */
00003 /*  Routines for Arbitrary Precision Floating-point Arithmetic               */
00004 /*  and Fast Robust Geometric Predicates                                     */
00005 /*  (predicates.c)                                                           */
00006 /*                                                                           */
00007 /*  May 18, 1996                                                             */
00008 /*                                                                           */
00009 /*  Placed in the public domain by                                           */
00010 /*  Jonathan Richard Shewchuk                                                */
00011 /*  School of Computer Science                                               */
00012 /*  Carnegie Mellon University                                               */
00013 /*  5000 Forbes Avenue                                                       */
00014 /*  Pittsburgh, Pennsylvania  15213-3891                                     */
00015 /*  jrs@cs.cmu.edu                                                           */
00016 /*                                                                           */
00017 /*  This file contains C implementation of algorithms for exact addition     */
00018 /*    and multiplication of floating-point numbers, and predicates for       */
00019 /*    robustly performing the orientation and incircle tests used in         */
00020 /*    computational geometry.  The algorithms and underlying theory are      */
00021 /*    described in Jonathan Richard Shewchuk.  "Adaptive Precision Floating- */
00022 /*    Point Arithmetic and Fast Robust Geometric Predicates."  Technical     */
00023 /*    Report CMU-CS-96-140, School of Computer Science, Carnegie Mellon      */
00024 /*    University, Pittsburgh, Pennsylvania, May 1996.  (Submitted to         */
00025 /*    Discrete & Computational Geometry.)                                    */
00026 /*                                                                           */
00027 /*  This file, the paper listed above, and other information are available   */
00028 /*    from the Web page http://www.cs.cmu.edu/~quake/robust.html .           */
00029 /*                                                                           */
00030 /*****************************************************************************/
00031 
00032 /*****************************************************************************/
00033 /*                                                                           */
00034 /*  Using this code:                                                         */
00035 /*                                                                           */
00036 /*  First, read the short or long version of the paper (from the Web page    */
00037 /*    above).                                                                */
00038 /*                                                                           */
00039 /*  Be sure to call exactinit() once, before calling any of the arithmetic   */
00040 /*    functions or geometric predicates.  Also be sure to turn on the        */
00041 /*    optimizer when compiling this file.                                    */
00042 /*                                                                           */
00043 /*                                                                           */
00044 /*  Several geometric predicates are defined.  Their parameters are all      */
00045 /*    points.  Each point is an array of two or three floating-point         */
00046 /*    numbers.  The geometric predicates, described in the papers, are       */
00047 /*                                                                           */
00048 /*    orient2d(pa, pb, pc)                                                   */
00049 /*    orient2dfast(pa, pb, pc)                                               */
00050 /*    orient3d(pa, pb, pc, pd)                                               */
00051 /*    orient3dfast(pa, pb, pc, pd)                                           */
00052 /*    incircle(pa, pb, pc, pd)                                               */
00053 /*    incirclefast(pa, pb, pc, pd)                                           */
00054 /*    insphere(pa, pb, pc, pd, pe)                                           */
00055 /*    inspherefast(pa, pb, pc, pd, pe)                                       */
00056 /*                                                                           */
00057 /*  Those with suffix "fast" are approximate, non-robust versions.  Those    */
00058 /*    without the suffix are adaptive precision, robust versions.  There     */
00059 /*    are also versions with the suffices "exact" and "slow", which are      */
00060 /*    non-adaptive, exact arithmetic versions, which I use only for timings  */
00061 /*    in my arithmetic papers.                                               */
00062 /*                                                                           */
00063 /*                                                                           */
00064 /*  An expansion is represented by an array of floating-point numbers,       */
00065 /*    sorted from smallest to largest magnitude (possibly with interspersed  */
00066 /*    zeros).  The length of each expansion is stored as a separate integer, */
00067 /*    and each arithmetic function returns an integer which is the length    */
00068 /*    of the expansion it created.                                           */
00069 /*                                                                           */
00070 /*  Several arithmetic functions are defined.  Their parameters are          */
00071 /*                                                                           */
00072 /*    e, f           Input expansions                                        */
00073 /*    elen, flen     Lengths of input expansions (must be >= 1)              */
00074 /*    h              Output expansion                                        */
00075 /*    b              Input scalar                                            */
00076 /*                                                                           */
00077 /*  The arithmetic functions are                                             */
00078 /*                                                                           */
00079 /*    grow_expansion(elen, e, b, h)                                          */
00080 /*    grow_expansion_zeroelim(elen, e, b, h)                                 */
00081 /*    expansion_sum(elen, e, flen, f, h)                                     */
00082 /*    expansion_sum_zeroelim1(elen, e, flen, f, h)                           */
00083 /*    expansion_sum_zeroelim2(elen, e, flen, f, h)                           */
00084 /*    fast_expansion_sum(elen, e, flen, f, h)                                */
00085 /*    fast_expansion_sum_zeroelim(elen, e, flen, f, h)                       */
00086 /*    linear_expansion_sum(elen, e, flen, f, h)                              */
00087 /*    linear_expansion_sum_zeroelim(elen, e, flen, f, h)                     */
00088 /*    scale_expansion(elen, e, b, h)                                         */
00089 /*    scale_expansion_zeroelim(elen, e, b, h)                                */
00090 /*    compress(elen, e, h)                                                   */
00091 /*                                                                           */
00092 /*  All of these are described in the long version of the paper; some are    */
00093 /*    described in the short version.  All return an integer that is the     */
00094 /*    length of h.  Those with suffix _zeroelim perform zero elimination,    */
00095 /*    and are recommended over their counterparts.  The procedure            */
00096 /*    fast_expansion_sum_zeroelim() (or linear_expansion_sum_zeroelim() on   */
00097 /*    processors that do not use the round-to-even tiebreaking rule) is      */
00098 /*    recommended over expansion_sum_zeroelim().  Each procedure has a       */
00099 /*    little note next to it (in the code below) that tells you whether or   */
00100 /*    not the output expansion may be the same array as one of the input     */
00101 /*    expansions.                                                            */
00102 /*                                                                           */
00103 /*                                                                           */
00104 /*  If you look around below, you'll also find macros for a bunch of         */
00105 /*    simple unrolled arithmetic operations, and procedures for printing     */
00106 /*    expansions (commented out because they don't work with all C           */
00107 /*    compilers) and for generating random floating-point numbers whose      */
00108 /*    significand bits are all random.  Most of the macros have undocumented */
00109 /*    requirements that certain of their parameters should not be the same   */
00110 /*    variable; for safety, better to make sure all the parameters are       */
00111 /*    distinct variables.  Feel free to send email to jrs@cs.cmu.edu if you  */
00112 /*    have questions.                                                        */
00113 /*                                                                           */
00114 /*****************************************************************************/
00115 
00116 #include <stdio.h>
00117 #include <stdlib.h>
00118 #include <math.h>
00119 #ifdef CPU86
00120 #include <float.h>
00121 #endif /* CPU86 */
00122 #ifdef LINUX
00123 #include <fpu_control.h>
00124 #endif /* LINUX */
00125 
00126 #include "tetgen.h"            // Defines the symbol REAL (float or double).
00127 
00128 /* On some machines, the exact arithmetic routines might be defeated by the  */
00129 /*   use of internal extended precision floating-point registers.  Sometimes */
00130 /*   this problem can be fixed by defining certain values to be volatile,    */
00131 /*   thus forcing them to be stored to memory and rounded off.  This isn't   */
00132 /*   a great solution, though, as it slows the arithmetic down.              */
00133 /*                                                                           */
00134 /* To try this out, write "#define INEXACT volatile" below.  Normally,       */
00135 /*   however, INEXACT should be defined to be nothing.  ("#define INEXACT".) */
00136 
00137 #define INEXACT                          /* Nothing */
00138 /* #define INEXACT volatile */
00139 
00140 /* #define REAL double */                      /* float or double */
00141 #define REALPRINT doubleprint
00142 #define REALRAND doublerand
00143 #define NARROWRAND narrowdoublerand
00144 #define UNIFORMRAND uniformdoublerand
00145 
00146 /* Which of the following two methods of finding the absolute values is      */
00147 /*   fastest is compiler-dependent.  A few compilers can inline and optimize */
00148 /*   the fabs() call; but most will incur the overhead of a function call,   */
00149 /*   which is disastrously slow.  A faster way on IEEE machines might be to  */
00150 /*   mask the appropriate bit, but that's difficult to do in C.              */
00151 
00152 #define Absolute(a)  ((a) >= 0.0 ? (a) : -(a))
00153 /* #define Absolute(a)  fabs(a) */
00154 
00155 /* Many of the operations are broken up into two pieces, a main part that    */
00156 /*   performs an approximate operation, and a "tail" that computes the       */
00157 /*   roundoff error of that operation.                                       */
00158 /*                                                                           */
00159 /* The operations Fast_Two_Sum(), Fast_Two_Diff(), Two_Sum(), Two_Diff(),    */
00160 /*   Split(), and Two_Product() are all implemented as described in the      */
00161 /*   reference.  Each of these macros requires certain variables to be       */
00162 /*   defined in the calling routine.  The variables `bvirt', `c', `abig',    */
00163 /*   `_i', `_j', `_k', `_l', `_m', and `_n' are declared `INEXACT' because   */
00164 /*   they store the result of an operation that may incur roundoff error.    */
00165 /*   The input parameter `x' (or the highest numbered `x_' parameter) must   */
00166 /*   also be declared `INEXACT'.                                             */
00167 
00168 #define Fast_Two_Sum_Tail(a, b, x, y) \
00169   bvirt = x - a; \
00170   y = b - bvirt
00171 
00172 #define Fast_Two_Sum(a, b, x, y) \
00173   x = (REAL) (a + b); \
00174   Fast_Two_Sum_Tail(a, b, x, y)
00175 
00176 #define Fast_Two_Diff_Tail(a, b, x, y) \
00177   bvirt = a - x; \
00178   y = bvirt - b
00179 
00180 #define Fast_Two_Diff(a, b, x, y) \
00181   x = (REAL) (a - b); \
00182   Fast_Two_Diff_Tail(a, b, x, y)
00183 
00184 #define Two_Sum_Tail(a, b, x, y) \
00185   bvirt = (REAL) (x - a); \
00186   avirt = x - bvirt; \
00187   bround = b - bvirt; \
00188   around = a - avirt; \
00189   y = around + bround
00190 
00191 #define Two_Sum(a, b, x, y) \
00192   x = (REAL) (a + b); \
00193   Two_Sum_Tail(a, b, x, y)
00194 
00195 #define Two_Diff_Tail(a, b, x, y) \
00196   bvirt = (REAL) (a - x); \
00197   avirt = x + bvirt; \
00198   bround = bvirt - b; \
00199   around = a - avirt; \
00200   y = around + bround
00201 
00202 #define Two_Diff(a, b, x, y) \
00203   x = (REAL) (a - b); \
00204   Two_Diff_Tail(a, b, x, y)
00205 
00206 #define Split(a, ahi, alo) \
00207   c = (REAL) (splitter * a); \
00208   abig = (REAL) (c - a); \
00209   ahi = c - abig; \
00210   alo = a - ahi
00211 
00212 #define Two_Product_Tail(a, b, x, y) \
00213   Split(a, ahi, alo); \
00214   Split(b, bhi, blo); \
00215   err1 = x - (ahi * bhi); \
00216   err2 = err1 - (alo * bhi); \
00217   err3 = err2 - (ahi * blo); \
00218   y = (alo * blo) - err3
00219 
00220 #define Two_Product(a, b, x, y) \
00221   x = (REAL) (a * b); \
00222   Two_Product_Tail(a, b, x, y)
00223 
00224 /* Two_Product_Presplit() is Two_Product() where one of the inputs has       */
00225 /*   already been split.  Avoids redundant splitting.                        */
00226 
00227 #define Two_Product_Presplit(a, b, bhi, blo, x, y) \
00228   x = (REAL) (a * b); \
00229   Split(a, ahi, alo); \
00230   err1 = x - (ahi * bhi); \
00231   err2 = err1 - (alo * bhi); \
00232   err3 = err2 - (ahi * blo); \
00233   y = (alo * blo) - err3
00234 
00235 /* Two_Product_2Presplit() is Two_Product() where both of the inputs have    */
00236 /*   already been split.  Avoids redundant splitting.                        */
00237 
00238 #define Two_Product_2Presplit(a, ahi, alo, b, bhi, blo, x, y) \
00239   x = (REAL) (a * b); \
00240   err1 = x - (ahi * bhi); \
00241   err2 = err1 - (alo * bhi); \
00242   err3 = err2 - (ahi * blo); \
00243   y = (alo * blo) - err3
00244 
00245 /* Square() can be done more quickly than Two_Product().                     */
00246 
00247 #define Square_Tail(a, x, y) \
00248   Split(a, ahi, alo); \
00249   err1 = x - (ahi * ahi); \
00250   err3 = err1 - ((ahi + ahi) * alo); \
00251   y = (alo * alo) - err3
00252 
00253 #define Square(a, x, y) \
00254   x = (REAL) (a * a); \
00255   Square_Tail(a, x, y)
00256 
00257 /* Macros for summing expansions of various fixed lengths.  These are all    */
00258 /*   unrolled versions of Expansion_Sum().                                   */
00259 
00260 #define Two_One_Sum(a1, a0, b, x2, x1, x0) \
00261   Two_Sum(a0, b , _i, x0); \
00262   Two_Sum(a1, _i, x2, x1)
00263 
00264 #define Two_One_Diff(a1, a0, b, x2, x1, x0) \
00265   Two_Diff(a0, b , _i, x0); \
00266   Two_Sum( a1, _i, x2, x1)
00267 
00268 #define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \
00269   Two_One_Sum(a1, a0, b0, _j, _0, x0); \
00270   Two_One_Sum(_j, _0, b1, x3, x2, x1)
00271 
00272 #define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \
00273   Two_One_Diff(a1, a0, b0, _j, _0, x0); \
00274   Two_One_Diff(_j, _0, b1, x3, x2, x1)
00275 
00276 #define Four_One_Sum(a3, a2, a1, a0, b, x4, x3, x2, x1, x0) \
00277   Two_One_Sum(a1, a0, b , _j, x1, x0); \
00278   Two_One_Sum(a3, a2, _j, x4, x3, x2)
00279 
00280 #define Four_Two_Sum(a3, a2, a1, a0, b1, b0, x5, x4, x3, x2, x1, x0) \
00281   Four_One_Sum(a3, a2, a1, a0, b0, _k, _2, _1, _0, x0); \
00282   Four_One_Sum(_k, _2, _1, _0, b1, x5, x4, x3, x2, x1)
00283 
00284 #define Four_Four_Sum(a3, a2, a1, a0, b4, b3, b1, b0, x7, x6, x5, x4, x3, x2, \
00285                       x1, x0) \
00286   Four_Two_Sum(a3, a2, a1, a0, b1, b0, _l, _2, _1, _0, x1, x0); \
00287   Four_Two_Sum(_l, _2, _1, _0, b4, b3, x7, x6, x5, x4, x3, x2)
00288 
00289 #define Eight_One_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b, x8, x7, x6, x5, x4, \
00290                       x3, x2, x1, x0) \
00291   Four_One_Sum(a3, a2, a1, a0, b , _j, x3, x2, x1, x0); \
00292   Four_One_Sum(a7, a6, a5, a4, _j, x8, x7, x6, x5, x4)
00293 
00294 #define Eight_Two_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b1, b0, x9, x8, x7, \
00295                       x6, x5, x4, x3, x2, x1, x0) \
00296   Eight_One_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b0, _k, _6, _5, _4, _3, _2, \
00297                 _1, _0, x0); \
00298   Eight_One_Sum(_k, _6, _5, _4, _3, _2, _1, _0, b1, x9, x8, x7, x6, x5, x4, \
00299                 x3, x2, x1)
00300 
00301 #define Eight_Four_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b4, b3, b1, b0, x11, \
00302                        x10, x9, x8, x7, x6, x5, x4, x3, x2, x1, x0) \
00303   Eight_Two_Sum(a7, a6, a5, a4, a3, a2, a1, a0, b1, b0, _l, _6, _5, _4, _3, \
00304                 _2, _1, _0, x1, x0); \
00305   Eight_Two_Sum(_l, _6, _5, _4, _3, _2, _1, _0, b4, b3, x11, x10, x9, x8, \
00306                 x7, x6, x5, x4, x3, x2)
00307 
00308 /* Macros for multiplying expansions of various fixed lengths.               */
00309 
00310 #define Two_One_Product(a1, a0, b, x3, x2, x1, x0) \
00311   Split(b, bhi, blo); \
00312   Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \
00313   Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
00314   Two_Sum(_i, _0, _k, x1); \
00315   Fast_Two_Sum(_j, _k, x3, x2)
00316 
00317 #define Four_One_Product(a3, a2, a1, a0, b, x7, x6, x5, x4, x3, x2, x1, x0) \
00318   Split(b, bhi, blo); \
00319   Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \
00320   Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
00321   Two_Sum(_i, _0, _k, x1); \
00322   Fast_Two_Sum(_j, _k, _i, x2); \
00323   Two_Product_Presplit(a2, b, bhi, blo, _j, _0); \
00324   Two_Sum(_i, _0, _k, x3); \
00325   Fast_Two_Sum(_j, _k, _i, x4); \
00326   Two_Product_Presplit(a3, b, bhi, blo, _j, _0); \
00327   Two_Sum(_i, _0, _k, x5); \
00328   Fast_Two_Sum(_j, _k, x7, x6)
00329 
00330 #define Two_Two_Product(a1, a0, b1, b0, x7, x6, x5, x4, x3, x2, x1, x0) \
00331   Split(a0, a0hi, a0lo); \
00332   Split(b0, bhi, blo); \
00333   Two_Product_2Presplit(a0, a0hi, a0lo, b0, bhi, blo, _i, x0); \
00334   Split(a1, a1hi, a1lo); \
00335   Two_Product_2Presplit(a1, a1hi, a1lo, b0, bhi, blo, _j, _0); \
00336   Two_Sum(_i, _0, _k, _1); \
00337   Fast_Two_Sum(_j, _k, _l, _2); \
00338   Split(b1, bhi, blo); \
00339   Two_Product_2Presplit(a0, a0hi, a0lo, b1, bhi, blo, _i, _0); \
00340   Two_Sum(_1, _0, _k, x1); \
00341   Two_Sum(_2, _k, _j, _1); \
00342   Two_Sum(_l, _j, _m, _2); \
00343   Two_Product_2Presplit(a1, a1hi, a1lo, b1, bhi, blo, _j, _0); \
00344   Two_Sum(_i, _0, _n, _0); \
00345   Two_Sum(_1, _0, _i, x2); \
00346   Two_Sum(_2, _i, _k, _1); \
00347   Two_Sum(_m, _k, _l, _2); \
00348   Two_Sum(_j, _n, _k, _0); \
00349   Two_Sum(_1, _0, _j, x3); \
00350   Two_Sum(_2, _j, _i, _1); \
00351   Two_Sum(_l, _i, _m, _2); \
00352   Two_Sum(_1, _k, _i, x4); \
00353   Two_Sum(_2, _i, _k, x5); \
00354   Two_Sum(_m, _k, x7, x6)
00355 
00356 /* An expansion of length two can be squared more quickly than finding the   */
00357 /*   product of two different expansions of length two, and the result is    */
00358 /*   guaranteed to have no more than six (rather than eight) components.     */
00359 
00360 #define Two_Square(a1, a0, x5, x4, x3, x2, x1, x0) \
00361   Square(a0, _j, x0); \
00362   _0 = a0 + a0; \
00363   Two_Product(a1, _0, _k, _1); \
00364   Two_One_Sum(_k, _1, _j, _l, _2, x1); \
00365   Square(a1, _j, _1); \
00366   Two_Two_Sum(_j, _1, _l, _2, x5, x4, x3, x2)
00367 
00368 /* splitter = 2^ceiling(p / 2) + 1.  Used to split floats in half.           */
00369 static REAL splitter;
00370 static REAL epsilon;         /* = 2^(-p).  Used to estimate roundoff errors. */
00371 /* A set of coefficients used to calculate maximum roundoff errors.          */
00372 static REAL resulterrbound;
00373 static REAL ccwerrboundA, ccwerrboundB, ccwerrboundC;
00374 static REAL o3derrboundA, o3derrboundB, o3derrboundC;
00375 static REAL iccerrboundA, iccerrboundB, iccerrboundC;
00376 static REAL isperrboundA, isperrboundB, isperrboundC;
00377 
00378 /*****************************************************************************/
00379 /*                                                                           */
00380 /*  doubleprint()   Print the bit representation of a double.                */
00381 /*                                                                           */
00382 /*  Useful for debugging exact arithmetic routines.                          */
00383 /*                                                                           */
00384 /*****************************************************************************/
00385 
00386 /*
00387 void doubleprint(number)
00388 double number;
00389 {
00390   unsigned long long no;
00391   unsigned long long sign, expo;
00392   int exponent;
00393   int i, bottomi;
00394 
00395   no = *(unsigned long long *) &number;
00396   sign = no & 0x8000000000000000ll;
00397   expo = (no >> 52) & 0x7ffll;
00398   exponent = (int) expo;
00399   exponent = exponent - 1023;
00400   if (sign) {
00401     printf("-");
00402   } else {
00403     printf(" ");
00404   }
00405   if (exponent == -1023) {
00406     printf(
00407       "0.0000000000000000000000000000000000000000000000000000_     (   )");
00408   } else {
00409     printf("1.");
00410     bottomi = -1;
00411     for (i = 0; i < 52; i++) {
00412       if (no & 0x0008000000000000ll) {
00413         printf("1");
00414         bottomi = i;
00415       } else {
00416         printf("0");
00417       }
00418       no <<= 1;
00419     }
00420     printf("_%d  (%d)", exponent, exponent - 1 - bottomi);
00421   }
00422 }
00423 */
00424 
00425 /*****************************************************************************/
00426 /*                                                                           */
00427 /*  floatprint()   Print the bit representation of a float.                  */
00428 /*                                                                           */
00429 /*  Useful for debugging exact arithmetic routines.                          */
00430 /*                                                                           */
00431 /*****************************************************************************/
00432 
00433 /*
00434 void floatprint(number)
00435 float number;
00436 {
00437   unsigned no;
00438   unsigned sign, expo;
00439   int exponent;
00440   int i, bottomi;
00441 
00442   no = *(unsigned *) &number;
00443   sign = no & 0x80000000;
00444   expo = (no >> 23) & 0xff;
00445   exponent = (int) expo;
00446   exponent = exponent - 127;
00447   if (sign) {
00448     printf("-");
00449   } else {
00450     printf(" ");
00451   }
00452   if (exponent == -127) {
00453     printf("0.00000000000000000000000_     (   )");
00454   } else {
00455     printf("1.");
00456     bottomi = -1;
00457     for (i = 0; i < 23; i++) {
00458       if (no & 0x00400000) {
00459         printf("1");
00460         bottomi = i;
00461       } else {
00462         printf("0");
00463       }
00464       no <<= 1;
00465     }
00466     printf("_%3d  (%3d)", exponent, exponent - 1 - bottomi);
00467   }
00468 }
00469 */
00470 
00471 /*****************************************************************************/
00472 /*                                                                           */
00473 /*  expansion_print()   Print the bit representation of an expansion.        */
00474 /*                                                                           */
00475 /*  Useful for debugging exact arithmetic routines.                          */
00476 /*                                                                           */
00477 /*****************************************************************************/
00478 
00479 /*
00480 void expansion_print(elen, e)
00481 int elen;
00482 REAL *e;
00483 {
00484   int i;
00485 
00486   for (i = elen - 1; i >= 0; i--) {
00487     REALPRINT(e[i]);
00488     if (i > 0) {
00489       printf(" +\n");
00490     } else {
00491       printf("\n");
00492     }
00493   }
00494 }
00495 */
00496 
00497 /*****************************************************************************/
00498 /*                                                                           */
00499 /*  doublerand()   Generate a double with random 53-bit significand and a    */
00500 /*                 random exponent in [0, 511].                              */
00501 /*                                                                           */
00502 /*****************************************************************************/
00503 
00504 /*
00505 double doublerand()
00506 {
00507   double result;
00508   double expo;
00509   long a, b, c;
00510   long i;
00511 
00512   a = random();
00513   b = random();
00514   c = random();
00515   result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8);
00516   for (i = 512, expo = 2; i <= 131072; i *= 2, expo = expo * expo) {
00517     if (c & i) {
00518       result *= expo;
00519     }
00520   }
00521   return result;
00522 }
00523 */
00524 
00525 /*****************************************************************************/
00526 /*                                                                           */
00527 /*  narrowdoublerand()   Generate a double with random 53-bit significand    */
00528 /*                       and a random exponent in [0, 7].                    */
00529 /*                                                                           */
00530 /*****************************************************************************/
00531 
00532 /*
00533 double narrowdoublerand()
00534 {
00535   double result;
00536   double expo;
00537   long a, b, c;
00538   long i;
00539 
00540   a = random();
00541   b = random();
00542   c = random();
00543   result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8);
00544   for (i = 512, expo = 2; i <= 2048; i *= 2, expo = expo * expo) {
00545     if (c & i) {
00546       result *= expo;
00547     }
00548   }
00549   return result;
00550 }
00551 */
00552 
00553 /*****************************************************************************/
00554 /*                                                                           */
00555 /*  uniformdoublerand()   Generate a double with random 53-bit significand.  */
00556 /*                                                                           */
00557 /*****************************************************************************/
00558 
00559 /*
00560 double uniformdoublerand()
00561 {
00562   double result;
00563   long a, b;
00564 
00565   a = random();
00566   b = random();
00567   result = (double) (a - 1073741824) * 8388608.0 + (double) (b >> 8);
00568   return result;
00569 }
00570 */
00571 
00572 /*****************************************************************************/
00573 /*                                                                           */
00574 /*  floatrand()   Generate a float with random 24-bit significand and a      */
00575 /*                random exponent in [0, 63].                                */
00576 /*                                                                           */
00577 /*****************************************************************************/
00578 
00579 /*
00580 float floatrand()
00581 {
00582   float result;
00583   float expo;
00584   long a, c;
00585   long i;
00586 
00587   a = random();
00588   c = random();
00589   result = (float) ((a - 1073741824) >> 6);
00590   for (i = 512, expo = 2; i <= 16384; i *= 2, expo = expo * expo) {
00591     if (c & i) {
00592       result *= expo;
00593     }
00594   }
00595   return result;
00596 }
00597 */
00598 
00599 /*****************************************************************************/
00600 /*                                                                           */
00601 /*  narrowfloatrand()   Generate a float with random 24-bit significand and  */
00602 /*                      a random exponent in [0, 7].                         */
00603 /*                                                                           */
00604 /*****************************************************************************/
00605 
00606 /*
00607 float narrowfloatrand()
00608 {
00609   float result;
00610   float expo;
00611   long a, c;
00612   long i;
00613 
00614   a = random();
00615   c = random();
00616   result = (float) ((a - 1073741824) >> 6);
00617   for (i = 512, expo = 2; i <= 2048; i *= 2, expo = expo * expo) {
00618     if (c & i) {
00619       result *= expo;
00620     }
00621   }
00622   return result;
00623 }
00624 */
00625 
00626 /*****************************************************************************/
00627 /*                                                                           */
00628 /*  uniformfloatrand()   Generate a float with random 24-bit significand.    */
00629 /*                                                                           */
00630 /*****************************************************************************/
00631 
00632 /*
00633 float uniformfloatrand()
00634 {
00635   float result;
00636   long a;
00637 
00638   a = random();
00639   result = (float) ((a - 1073741824) >> 6);
00640   return result;
00641 }
00642 */
00643 namespace tetgen
00644 { //Added namespace to avoid clash with triangle
00645 
00646 /*****************************************************************************/
00647 /*                                                                           */
00648 /*  exactinit()   Initialize the variables used for exact arithmetic.        */
00649 /*                                                                           */
00650 /*  `epsilon' is the largest power of two such that 1.0 + epsilon = 1.0 in   */
00651 /*  floating-point arithmetic.  `epsilon' bounds the relative roundoff       */
00652 /*  error.  It is used for floating-point error analysis.                    */
00653 /*                                                                           */
00654 /*  `splitter' is used to split floating-point numbers into two half-        */
00655 /*  length significands for exact multiplication.                            */
00656 /*                                                                           */
00657 /*  I imagine that a highly optimizing compiler might be too smart for its   */
00658 /*  own good, and somehow cause this routine to fail, if it pretends that    */
00659 /*  floating-point arithmetic is too much like real arithmetic.              */
00660 /*                                                                           */
00661 /*  Don't change this routine unless you fully understand it.                */
00662 /*                                                                           */
00663 /*****************************************************************************/
00664 
00665 REAL exactinit()
00666 {
00667   REAL half;
00668   REAL check, lastcheck;
00669   int every_other;
00670 #ifdef LINUX
00671   int cword;
00672 #endif /* LINUX */
00673 
00674 #ifdef CPU86
00675 #ifdef SINGLE
00676   _control87(_PC_24, _MCW_PC); /* Set FPU control word for single precision. */
00677 #else /* not SINGLE */
00678   _control87(_PC_53, _MCW_PC); /* Set FPU control word for double precision. */
00679 #endif /* not SINGLE */
00680 #endif /* CPU86 */
00681 #ifdef LINUX
00682 #ifdef SINGLE
00683   /*  cword = 4223; */
00684   cword = 4210;                 /* set FPU control word for single precision */
00685 #else /* not SINGLE */
00686   /*  cword = 4735; */
00687   cword = 4722;                 /* set FPU control word for double precision */
00688 #endif /* not SINGLE */
00689   _FPU_SETCW(cword);
00690 #endif /* LINUX */
00691 
00692   every_other = 1;
00693   half = 0.5;
00694   epsilon = 1.0;
00695   splitter = 1.0;
00696   check = 1.0;
00697   /* Repeatedly divide `epsilon' by two until it is too small to add to    */
00698   /*   one without causing roundoff.  (Also check if the sum is equal to   */
00699   /*   the previous sum, for machines that round up instead of using exact */
00700   /*   rounding.  Not that this library will work on such machines anyway. */
00701   do {
00702     lastcheck = check;
00703     epsilon *= half;
00704     if (every_other) {
00705       splitter *= 2.0;
00706     }
00707     every_other = !every_other;
00708     check = 1.0 + epsilon;
00709   } while ((check != 1.0) && (check != lastcheck));
00710   splitter += 1.0;
00711 
00712   /* Error bounds for orientation and incircle tests. */
00713   resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
00714   ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon;
00715   ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon;
00716   ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon;
00717   o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon;
00718   o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon;
00719   o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;
00720   iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon;
00721   iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon;
00722   iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon;
00723   isperrboundA = (16.0 + 224.0 * epsilon) * epsilon;
00724   isperrboundB = (5.0 + 72.0 * epsilon) * epsilon;
00725   isperrboundC = (71.0 + 1408.0 * epsilon) * epsilon * epsilon;
00726 
00727   return epsilon; /* Added by H. Si 30 Juli, 2004. */
00728 }
00729 
00730 /*****************************************************************************/
00731 /*                                                                           */
00732 /*  grow_expansion()   Add a scalar to an expansion.                         */
00733 /*                                                                           */
00734 /*  Sets h = e + b.  See the long version of my paper for details.           */
00735 /*                                                                           */
00736 /*  Maintains the nonoverlapping property.  If round-to-even is used (as     */
00737 /*  with IEEE 754), maintains the strongly nonoverlapping and nonadjacent    */
00738 /*  properties as well.  (That is, if e has one of these properties, so      */
00739 /*  will h.)                                                                 */
00740 /*                                                                           */
00741 /*****************************************************************************/
00742 
00743 int grow_expansion(int elen, REAL *e, REAL b, REAL *h)
00744 /* e and h can be the same. */
00745 {
00746   REAL Q;
00747   INEXACT REAL Qnew;
00748   int eindex;
00749   REAL enow;
00750   INEXACT REAL bvirt;
00751   REAL avirt, bround, around;
00752 
00753   Q = b;
00754   for (eindex = 0; eindex < elen; eindex++) {
00755     enow = e[eindex];
00756     Two_Sum(Q, enow, Qnew, h[eindex]);
00757     Q = Qnew;
00758   }
00759   h[eindex] = Q;
00760   return eindex + 1;
00761 }
00762 
00763 /*****************************************************************************/
00764 /*                                                                           */
00765 /*  grow_expansion_zeroelim()   Add a scalar to an expansion, eliminating    */
00766 /*                              zero components from the output expansion.   */
00767 /*                                                                           */
00768 /*  Sets h = e + b.  See the long version of my paper for details.           */
00769 /*                                                                           */
00770 /*  Maintains the nonoverlapping property.  If round-to-even is used (as     */
00771 /*  with IEEE 754), maintains the strongly nonoverlapping and nonadjacent    */
00772 /*  properties as well.  (That is, if e has one of these properties, so      */
00773 /*  will h.)                                                                 */
00774 /*                                                                           */
00775 /*****************************************************************************/
00776 
00777 int grow_expansion_zeroelim(int elen, REAL *e, REAL b, REAL *h)
00778 /* e and h can be the same. */
00779 {
00780   REAL Q, hh;
00781   INEXACT REAL Qnew;
00782   int eindex, hindex;
00783   REAL enow;
00784   INEXACT REAL bvirt;
00785   REAL avirt, bround, around;
00786 
00787   hindex = 0;
00788   Q = b;
00789   for (eindex = 0; eindex < elen; eindex++) {
00790     enow = e[eindex];
00791     Two_Sum(Q, enow, Qnew, hh);
00792     Q = Qnew;
00793     if (hh != 0.0) {
00794       h[hindex++] = hh;
00795     }
00796   }
00797   if ((Q != 0.0) || (hindex == 0)) {
00798     h[hindex++] = Q;
00799   }
00800   return hindex;
00801 }
00802 
00803 /*****************************************************************************/
00804 /*                                                                           */
00805 /*  expansion_sum()   Sum two expansions.                                    */
00806 /*                                                                           */
00807 /*  Sets h = e + f.  See the long version of my paper for details.           */
00808 /*                                                                           */
00809 /*  Maintains the nonoverlapping property.  If round-to-even is used (as     */
00810 /*  with IEEE 754), maintains the nonadjacent property as well.  (That is,   */
00811 /*  if e has one of these properties, so will h.)  Does NOT maintain the     */
00812 /*  strongly nonoverlapping property.                                        */
00813 /*                                                                           */
00814 /*****************************************************************************/
00815 
00816 int expansion_sum(int elen, REAL *e, int flen, REAL *f, REAL *h)
00817 /* e and h can be the same, but f and h cannot. */
00818 {
00819   REAL Q;
00820   INEXACT REAL Qnew;
00821   int findex, hindex, hlast;
00822   REAL hnow;
00823   INEXACT REAL bvirt;
00824   REAL avirt, bround, around;
00825 
00826   Q = f[0];
00827   for (hindex = 0; hindex < elen; hindex++) {
00828     hnow = e[hindex];
00829     Two_Sum(Q, hnow, Qnew, h[hindex]);
00830     Q = Qnew;
00831   }
00832   h[hindex] = Q;
00833   hlast = hindex;
00834   for (findex = 1; findex < flen; findex++) {
00835     Q = f[findex];
00836     for (hindex = findex; hindex <= hlast; hindex++) {
00837       hnow = h[hindex];
00838       Two_Sum(Q, hnow, Qnew, h[hindex]);
00839       Q = Qnew;
00840     }
00841     h[++hlast] = Q;
00842   }
00843   return hlast + 1;
00844 }
00845 
00846 /*****************************************************************************/
00847 /*                                                                           */
00848 /*  expansion_sum_zeroelim1()   Sum two expansions, eliminating zero         */
00849 /*                              components from the output expansion.        */
00850 /*                                                                           */
00851 /*  Sets h = e + f.  See the long version of my paper for details.           */
00852 /*                                                                           */
00853 /*  Maintains the nonoverlapping property.  If round-to-even is used (as     */
00854 /*  with IEEE 754), maintains the nonadjacent property as well.  (That is,   */
00855 /*  if e has one of these properties, so will h.)  Does NOT maintain the     */
00856 /*  strongly nonoverlapping property.                                        */
00857 /*                                                                           */
00858 /*****************************************************************************/
00859 
00860 int expansion_sum_zeroelim1(int elen, REAL *e, int flen, REAL *f, REAL *h)
00861 /* e and h can be the same, but f and h cannot. */
00862 {
00863   REAL Q;
00864   INEXACT REAL Qnew;
00865   int index, findex, hindex, hlast;
00866   REAL hnow;
00867   INEXACT REAL bvirt;
00868   REAL avirt, bround, around;
00869 
00870   Q = f[0];
00871   for (hindex = 0; hindex < elen; hindex++) {
00872     hnow = e[hindex];
00873     Two_Sum(Q, hnow, Qnew, h[hindex]);
00874     Q = Qnew;
00875   }
00876   h[hindex] = Q;
00877   hlast = hindex;
00878   for (findex = 1; findex < flen; findex++) {
00879     Q = f[findex];
00880     for (hindex = findex; hindex <= hlast; hindex++) {
00881       hnow = h[hindex];
00882       Two_Sum(Q, hnow, Qnew, h[hindex]);
00883       Q = Qnew;
00884     }
00885     h[++hlast] = Q;
00886   }
00887   hindex = -1;
00888   for (index = 0; index <= hlast; index++) {
00889     hnow = h[index];
00890     if (hnow != 0.0) {
00891       h[++hindex] = hnow;
00892     }
00893   }
00894   if (hindex == -1) {
00895     return 1;
00896   } else {
00897     return hindex + 1;
00898   }
00899 }
00900 
00901 /*****************************************************************************/
00902 /*                                                                           */
00903 /*  expansion_sum_zeroelim2()   Sum two expansions, eliminating zero         */
00904 /*                              components from the output expansion.        */
00905 /*                                                                           */
00906 /*  Sets h = e + f.  See the long version of my paper for details.           */
00907 /*                                                                           */
00908 /*  Maintains the nonoverlapping property.  If round-to-even is used (as     */
00909 /*  with IEEE 754), maintains the nonadjacent property as well.  (That is,   */
00910 /*  if e has one of these properties, so will h.)  Does NOT maintain the     */
00911 /*  strongly nonoverlapping property.                                        */
00912 /*                                                                           */
00913 /*****************************************************************************/
00914 
00915 int expansion_sum_zeroelim2(int elen, REAL *e, int flen, REAL *f, REAL *h)
00916 /* e and h can be the same, but f and h cannot. */
00917 {
00918   REAL Q, hh;
00919   INEXACT REAL Qnew;
00920   int eindex, findex, hindex, hlast;
00921   REAL enow;
00922   INEXACT REAL bvirt;
00923   REAL avirt, bround, around;
00924 
00925   hindex = 0;
00926   Q = f[0];
00927   for (eindex = 0; eindex < elen; eindex++) {
00928     enow = e[eindex];
00929     Two_Sum(Q, enow, Qnew, hh);
00930     Q = Qnew;
00931     if (hh != 0.0) {
00932       h[hindex++] = hh;
00933     }
00934   }
00935   h[hindex] = Q;
00936   hlast = hindex;
00937   for (findex = 1; findex < flen; findex++) {
00938     hindex = 0;
00939     Q = f[findex];
00940     for (eindex = 0; eindex <= hlast; eindex++) {
00941       enow = h[eindex];
00942       Two_Sum(Q, enow, Qnew, hh);
00943       Q = Qnew;
00944       if (hh != 0) {
00945         h[hindex++] = hh;
00946       }
00947     }
00948     h[hindex] = Q;
00949     hlast = hindex;
00950   }
00951   return hlast + 1;
00952 }
00953 
00954 /*****************************************************************************/
00955 /*                                                                           */
00956 /*  fast_expansion_sum()   Sum two expansions.                               */
00957 /*                                                                           */
00958 /*  Sets h = e + f.  See the long version of my paper for details.           */
00959 /*                                                                           */
00960 /*  If round-to-even is used (as with IEEE 754), maintains the strongly      */
00961 /*  nonoverlapping property.  (That is, if e is strongly nonoverlapping, h   */
00962 /*  will be also.)  Does NOT maintain the nonoverlapping or nonadjacent      */
00963 /*  properties.                                                              */
00964 /*                                                                           */
00965 /*****************************************************************************/
00966 
00967 int fast_expansion_sum(int elen, REAL *e, int flen, REAL *f, REAL *h)
00968 /* h cannot be e or f. */
00969 {
00970   REAL Q;
00971   INEXACT REAL Qnew;
00972   INEXACT REAL bvirt;
00973   REAL avirt, bround, around;
00974   int eindex, findex, hindex;
00975   REAL enow, fnow;
00976 
00977   enow = e[0];
00978   fnow = f[0];
00979   eindex = findex = 0;
00980   if ((fnow > enow) == (fnow > -enow)) {
00981     Q = enow;
00982     enow = e[++eindex];
00983   } else {
00984     Q = fnow;
00985     fnow = f[++findex];
00986   }
00987   hindex = 0;
00988   if ((eindex < elen) && (findex < flen)) {
00989     if ((fnow > enow) == (fnow > -enow)) {
00990       Fast_Two_Sum(enow, Q, Qnew, h[0]);
00991       enow = e[++eindex];
00992     } else {
00993       Fast_Two_Sum(fnow, Q, Qnew, h[0]);
00994       fnow = f[++findex];
00995     }
00996     Q = Qnew;
00997     hindex = 1;
00998     while ((eindex < elen) && (findex < flen)) {
00999       if ((fnow > enow) == (fnow > -enow)) {
01000         Two_Sum(Q, enow, Qnew, h[hindex]);
01001         enow = e[++eindex];
01002       } else {
01003         Two_Sum(Q, fnow, Qnew, h[hindex]);
01004         fnow = f[++findex];
01005       }
01006       Q = Qnew;
01007       hindex++;
01008     }
01009   }
01010   while (eindex < elen) {
01011     Two_Sum(Q, enow, Qnew, h[hindex]);
01012     enow = e[++eindex];
01013     Q = Qnew;
01014     hindex++;
01015   }
01016   while (findex < flen) {
01017     Two_Sum(Q, fnow, Qnew, h[hindex]);
01018     fnow = f[++findex];
01019     Q = Qnew;
01020     hindex++;
01021   }
01022   h[hindex] = Q;
01023   return hindex + 1;
01024 }
01025 
01026 /*****************************************************************************/
01027 /*                                                                           */
01028 /*  fast_expansion_sum_zeroelim()   Sum two expansions, eliminating zero     */
01029 /*                                  components from the output expansion.    */
01030 /*                                                                           */
01031 /*  Sets h = e + f.  See the long version of my paper for details.           */
01032 /*                                                                           */
01033 /*  If round-to-even is used (as with IEEE 754), maintains the strongly      */
01034 /*  nonoverlapping property.  (That is, if e is strongly nonoverlapping, h   */
01035 /*  will be also.)  Does NOT maintain the nonoverlapping or nonadjacent      */
01036 /*  properties.                                                              */
01037 /*                                                                           */
01038 /*****************************************************************************/
01039 
01040 int fast_expansion_sum_zeroelim(int elen, REAL *e, int flen, REAL *f, REAL *h)
01041 /* h cannot be e or f. */
01042 {
01043   REAL Q;
01044   INEXACT REAL Qnew;
01045   INEXACT REAL hh;
01046   INEXACT REAL bvirt;
01047   REAL avirt, bround, around;
01048   int eindex, findex, hindex;
01049   REAL enow, fnow;
01050 
01051   enow = e[0];
01052   fnow = f[0];
01053   eindex = findex = 0;
01054   if ((fnow > enow) == (fnow > -enow)) {
01055     Q = enow;
01056     enow = e[++eindex];
01057   } else {
01058     Q = fnow;
01059     fnow = f[++findex];
01060   }
01061   hindex = 0;
01062   if ((eindex < elen) && (findex < flen)) {
01063     if ((fnow > enow) == (fnow > -enow)) {
01064       Fast_Two_Sum(enow, Q, Qnew, hh);
01065       enow = e[++eindex];
01066     } else {
01067       Fast_Two_Sum(fnow, Q, Qnew, hh);
01068       fnow = f[++findex];
01069     }
01070     Q = Qnew;
01071     if (hh != 0.0) {
01072       h[hindex++] = hh;
01073     }
01074     while ((eindex < elen) && (findex < flen)) {
01075       if ((fnow > enow) == (fnow > -enow)) {
01076         Two_Sum(Q, enow, Qnew, hh);
01077         enow = e[++eindex];
01078       } else {
01079         Two_Sum(Q, fnow, Qnew, hh);
01080         fnow = f[++findex];
01081       }
01082       Q = Qnew;
01083       if (hh != 0.0) {
01084         h[hindex++] = hh;
01085       }
01086     }
01087   }
01088   while (eindex < elen) {
01089     Two_Sum(Q, enow, Qnew, hh);
01090     enow = e[++eindex];
01091     Q = Qnew;
01092     if (hh != 0.0) {
01093       h[hindex++] = hh;
01094     }
01095   }
01096   while (findex < flen) {
01097     Two_Sum(Q, fnow, Qnew, hh);
01098     fnow = f[++findex];
01099     Q = Qnew;
01100     if (hh != 0.0) {
01101       h[hindex++] = hh;
01102     }
01103   }
01104   if ((Q != 0.0) || (hindex == 0)) {
01105     h[hindex++] = Q;
01106   }
01107   return hindex;
01108 }
01109 
01110 /*****************************************************************************/
01111 /*                                                                           */
01112 /*  linear_expansion_sum()   Sum two expansions.                             */
01113 /*                                                                           */
01114 /*  Sets h = e + f.  See either version of my paper for details.             */
01115 /*                                                                           */
01116 /*  Maintains the nonoverlapping property.  (That is, if e is                */
01117 /*  nonoverlapping, h will be also.)                                         */
01118 /*                                                                           */
01119 /*****************************************************************************/
01120 
01121 int linear_expansion_sum(int elen, REAL *e, int flen, REAL *f, REAL *h)
01122 /* h cannot be e or f. */
01123 {
01124   REAL Q, q;
01125   INEXACT REAL Qnew;
01126   INEXACT REAL R;
01127   INEXACT REAL bvirt;
01128   REAL avirt, bround, around;
01129   int eindex, findex, hindex;
01130   REAL enow, fnow;
01131   REAL g0;
01132 
01133   enow = e[0];
01134   fnow = f[0];
01135   eindex = findex = 0;
01136   if ((fnow > enow) == (fnow > -enow)) {
01137     g0 = enow;
01138     enow = e[++eindex];
01139   } else {
01140     g0 = fnow;
01141     fnow = f[++findex];
01142   }
01143   if ((eindex < elen) && ((findex >= flen)
01144                           || ((fnow > enow) == (fnow > -enow)))) {
01145     Fast_Two_Sum(enow, g0, Qnew, q);
01146     enow = e[++eindex];
01147   } else {
01148     Fast_Two_Sum(fnow, g0, Qnew, q);
01149     fnow = f[++findex];
01150   }
01151   Q = Qnew;
01152   for (hindex = 0; hindex < elen + flen - 2; hindex++) {
01153     if ((eindex < elen) && ((findex >= flen)
01154                             || ((fnow > enow) == (fnow > -enow)))) {
01155       Fast_Two_Sum(enow, q, R, h[hindex]);
01156       enow = e[++eindex];
01157     } else {
01158       Fast_Two_Sum(fnow, q, R, h[hindex]);
01159       fnow = f[++findex];
01160     }
01161     Two_Sum(Q, R, Qnew, q);
01162     Q = Qnew;
01163   }
01164   h[hindex] = q;
01165   h[hindex + 1] = Q;
01166   return hindex + 2;
01167 }
01168 
01169 /*****************************************************************************/
01170 /*                                                                           */
01171 /*  linear_expansion_sum_zeroelim()   Sum two expansions, eliminating zero   */
01172 /*                                    components from the output expansion.  */
01173 /*                                                                           */
01174 /*  Sets h = e + f.  See either version of my paper for details.             */
01175 /*                                                                           */
01176 /*  Maintains the nonoverlapping property.  (That is, if e is                */
01177 /*  nonoverlapping, h will be also.)                                         */
01178 /*                                                                           */
01179 /*****************************************************************************/
01180 
01181 int linear_expansion_sum_zeroelim(int elen, REAL *e, int flen, REAL *f,
01182                                   REAL *h)
01183 /* h cannot be e or f. */
01184 {
01185   REAL Q, q, hh;
01186   INEXACT REAL Qnew;
01187   INEXACT REAL R;
01188   INEXACT REAL bvirt;
01189   REAL avirt, bround, around;
01190   int eindex, findex, hindex;
01191   int count;
01192   REAL enow, fnow;
01193   REAL g0;
01194 
01195   enow = e[0];
01196   fnow = f[0];
01197   eindex = findex = 0;
01198   hindex = 0;
01199   if ((fnow > enow) == (fnow > -enow)) {
01200     g0 = enow;
01201     enow = e[++eindex];
01202   } else {
01203     g0 = fnow;
01204     fnow = f[++findex];
01205   }
01206   if ((eindex < elen) && ((findex >= flen)
01207                           || ((fnow > enow) == (fnow > -enow)))) {
01208     Fast_Two_Sum(enow, g0, Qnew, q);
01209     enow = e[++eindex];
01210   } else {
01211     Fast_Two_Sum(fnow, g0, Qnew, q);
01212     fnow = f[++findex];
01213   }
01214   Q = Qnew;
01215   for (count = 2; count < elen + flen; count++) {
01216     if ((eindex < elen) && ((findex >= flen)
01217                             || ((fnow > enow) == (fnow > -enow)))) {
01218       Fast_Two_Sum(enow, q, R, hh);
01219       enow = e[++eindex];
01220     } else {
01221       Fast_Two_Sum(fnow, q, R, hh);
01222       fnow = f[++findex];
01223     }
01224     Two_Sum(Q, R, Qnew, q);
01225     Q = Qnew;
01226     if (hh != 0) {
01227       h[hindex++] = hh;
01228     }
01229   }
01230   if (q != 0) {
01231     h[hindex++] = q;
01232   }
01233   if ((Q != 0.0) || (hindex == 0)) {
01234     h[hindex++] = Q;
01235   }
01236   return hindex;
01237 }
01238 
01239 /*****************************************************************************/
01240 /*                                                                           */
01241 /*  scale_expansion()   Multiply an expansion by a scalar.                   */
01242 /*                                                                           */
01243 /*  Sets h = be.  See either version of my paper for details.                */
01244 /*                                                                           */
01245 /*  Maintains the nonoverlapping property.  If round-to-even is used (as     */
01246 /*  with IEEE 754), maintains the strongly nonoverlapping and nonadjacent    */
01247 /*  properties as well.  (That is, if e has one of these properties, so      */
01248 /*  will h.)                                                                 */
01249 /*                                                                           */
01250 /*****************************************************************************/
01251 
01252 int scale_expansion(int elen, REAL *e, REAL b, REAL *h)
01253 /* e and h cannot be the same. */
01254 {
01255   INEXACT REAL Q;
01256   INEXACT REAL sum;
01257   INEXACT REAL product1;
01258   REAL product0;
01259   int eindex, hindex;
01260   REAL enow;
01261   INEXACT REAL bvirt;
01262   REAL avirt, bround, around;
01263   INEXACT REAL c;
01264   INEXACT REAL abig;
01265   REAL ahi, alo, bhi, blo;
01266   REAL err1, err2, err3;
01267 
01268   Split(b, bhi, blo);
01269   Two_Product_Presplit(e[0], b, bhi, blo, Q, h[0]);
01270   hindex = 1;
01271   for (eindex = 1; eindex < elen; eindex++) {
01272     enow = e[eindex];
01273     Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
01274     Two_Sum(Q, product0, sum, h[hindex]);
01275     hindex++;
01276     Two_Sum(product1, sum, Q, h[hindex]);
01277     hindex++;
01278   }
01279   h[hindex] = Q;
01280   return elen + elen;
01281 }
01282 
01283 /*****************************************************************************/
01284 /*                                                                           */
01285 /*  scale_expansion_zeroelim()   Multiply an expansion by a scalar,          */
01286 /*                               eliminating zero components from the        */
01287 /*                               output expansion.                           */
01288 /*                                                                           */
01289 /*  Sets h = be.  See either version of my paper for details.                */
01290 /*                                                                           */
01291 /*  Maintains the nonoverlapping property.  If round-to-even is used (as     */
01292 /*  with IEEE 754), maintains the strongly nonoverlapping and nonadjacent    */
01293 /*  properties as well.  (That is, if e has one of these properties, so      */
01294 /*  will h.)                                                                 */
01295 /*                                                                           */
01296 /*****************************************************************************/
01297 
01298 int scale_expansion_zeroelim(int elen, REAL *e, REAL b, REAL *h)
01299 /* e and h cannot be the same. */
01300 {
01301   INEXACT REAL Q, sum;
01302   REAL hh;
01303   INEXACT REAL product1;
01304   REAL product0;
01305   int eindex, hindex;
01306   REAL enow;
01307   INEXACT REAL bvirt;
01308   REAL avirt, bround, around;
01309   INEXACT REAL c;
01310   INEXACT REAL abig;
01311   REAL ahi, alo, bhi, blo;
01312   REAL err1, err2, err3;
01313 
01314   Split(b, bhi, blo);
01315   Two_Product_Presplit(e[0], b, bhi, blo, Q, hh);
01316   hindex = 0;
01317   if (hh != 0) {
01318     h[hindex++] = hh;
01319   }
01320   for (eindex = 1; eindex < elen; eindex++) {
01321     enow = e[eindex];
01322     Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
01323     Two_Sum(Q, product0, sum, hh);
01324     if (hh != 0) {
01325       h[hindex++] = hh;
01326     }
01327     Fast_Two_Sum(product1, sum, Q, hh);
01328     if (hh != 0) {
01329       h[hindex++] = hh;
01330     }
01331   }
01332   if ((Q != 0.0) || (hindex == 0)) {
01333     h[hindex++] = Q;
01334   }
01335   return hindex;
01336 }
01337 
01338 /*****************************************************************************/
01339 /*                                                                           */
01340 /*  compress()   Compress an expansion.                                      */
01341 /*                                                                           */
01342 /*  See the long version of my paper for details.                            */
01343 /*                                                                           */
01344 /*  Maintains the nonoverlapping property.  If round-to-even is used (as     */
01345 /*  with IEEE 754), then any nonoverlapping expansion is converted to a      */
01346 /*  nonadjacent expansion.                                                   */
01347 /*                                                                           */
01348 /*****************************************************************************/
01349 
01350 int compress(int elen, REAL *e, REAL *h)
01351 /* e and h may be the same. */
01352 {
01353   REAL Q, q;
01354   INEXACT REAL Qnew;
01355   int eindex, hindex;
01356   INEXACT REAL bvirt;
01357   REAL enow, hnow;
01358   int top, bottom;
01359 
01360   bottom = elen - 1;
01361   Q = e[bottom];
01362   for (eindex = elen - 2; eindex >= 0; eindex--) {
01363     enow = e[eindex];
01364     Fast_Two_Sum(Q, enow, Qnew, q);
01365     if (q != 0) {
01366       h[bottom--] = Qnew;
01367       Q = q;
01368     } else {
01369       Q = Qnew;
01370     }
01371   }
01372   top = 0;
01373   for (hindex = bottom + 1; hindex < elen; hindex++) {
01374     hnow = h[hindex];
01375     Fast_Two_Sum(hnow, Q, Qnew, q);
01376     if (q != 0) {
01377       h[top++] = q;
01378     }
01379     Q = Qnew;
01380   }
01381   h[top] = Q;
01382   return top + 1;
01383 }
01384 
01385 /*****************************************************************************/
01386 /*                                                                           */
01387 /*  estimate()   Produce a one-word estimate of an expansion's value.        */
01388 /*                                                                           */
01389 /*  See either version of my paper for details.                              */
01390 /*                                                                           */
01391 /*****************************************************************************/
01392 
01393 REAL estimate(int elen, REAL *e)
01394 {
01395   REAL Q;
01396   int eindex;
01397 
01398   Q = e[0];
01399   for (eindex = 1; eindex < elen; eindex++) {
01400     Q += e[eindex];
01401   }
01402   return Q;
01403 }
01404 
01405 /*****************************************************************************/
01406 /*                                                                           */
01407 /*  orient2dfast()   Approximate 2D orientation test.  Nonrobust.            */
01408 /*  orient2dexact()   Exact 2D orientation test.  Robust.                    */
01409 /*  orient2dslow()   Another exact 2D orientation test.  Robust.             */
01410 /*  orient2d()   Adaptive exact 2D orientation test.  Robust.                */
01411 /*                                                                           */
01412 /*               Return a positive value if the points pa, pb, and pc occur  */
01413 /*               in counterclockwise order; a negative value if they occur   */
01414 /*               in clockwise order; and zero if they are collinear.  The    */
01415 /*               result is also a rough approximation of twice the signed    */
01416 /*               area of the triangle defined by the three points.           */
01417 /*                                                                           */
01418 /*  Only the first and last routine should be used; the middle two are for   */
01419 /*  timings.                                                                 */
01420 /*                                                                           */
01421 /*  The last three use exact arithmetic to ensure a correct answer.  The     */
01422 /*  result returned is the determinant of a matrix.  In orient2d() only,     */
01423 /*  this determinant is computed adaptively, in the sense that exact         */
01424 /*  arithmetic is used only to the degree it is needed to ensure that the    */
01425 /*  returned value has the correct sign.  Hence, orient2d() is usually quite */
01426 /*  fast, but will run more slowly when the input points are collinear or    */
01427 /*  nearly so.                                                               */
01428 /*                                                                           */
01429 /*****************************************************************************/
01430 
01431 REAL orient2dfast(REAL *pa, REAL *pb, REAL *pc)
01432 {
01433   REAL acx, bcx, acy, bcy;
01434 
01435   acx = pa[0] - pc[0];
01436   bcx = pb[0] - pc[0];
01437   acy = pa[1] - pc[1];
01438   bcy = pb[1] - pc[1];
01439   return acx * bcy - acy * bcx;
01440 }
01441 
01442 REAL orient2dexact(REAL *pa, REAL *pb, REAL *pc)
01443 {
01444   INEXACT REAL axby1, axcy1, bxcy1, bxay1, cxay1, cxby1;
01445   REAL axby0, axcy0, bxcy0, bxay0, cxay0, cxby0;
01446   REAL aterms[4], bterms[4], cterms[4];
01447   INEXACT REAL aterms3, bterms3, cterms3;
01448   REAL v[8], w[12];
01449   int vlength, wlength;
01450 
01451   INEXACT REAL bvirt;
01452   REAL avirt, bround, around;
01453   INEXACT REAL c;
01454   INEXACT REAL abig;
01455   REAL ahi, alo, bhi, blo;
01456   REAL err1, err2, err3;
01457   INEXACT REAL _i, _j;
01458   REAL _0;
01459 
01460   Two_Product(pa[0], pb[1], axby1, axby0);
01461   Two_Product(pa[0], pc[1], axcy1, axcy0);
01462   Two_Two_Diff(axby1, axby0, axcy1, axcy0,
01463                aterms3, aterms[2], aterms[1], aterms[0]);
01464   aterms[3] = aterms3;
01465 
01466   Two_Product(pb[0], pc[1], bxcy1, bxcy0);
01467   Two_Product(pb[0], pa[1], bxay1, bxay0);
01468   Two_Two_Diff(bxcy1, bxcy0, bxay1, bxay0,
01469                bterms3, bterms[2], bterms[1], bterms[0]);
01470   bterms[3] = bterms3;
01471 
01472   Two_Product(pc[0], pa[1], cxay1, cxay0);
01473   Two_Product(pc[0], pb[1], cxby1, cxby0);
01474   Two_Two_Diff(cxay1, cxay0, cxby1, cxby0,
01475                cterms3, cterms[2], cterms[1], cterms[0]);
01476   cterms[3] = cterms3;
01477 
01478   vlength = fast_expansion_sum_zeroelim(4, aterms, 4, bterms, v);
01479   wlength = fast_expansion_sum_zeroelim(vlength, v, 4, cterms, w);
01480 
01481   return w[wlength - 1];
01482 }
01483 
01484 REAL orient2dslow(REAL *pa, REAL *pb, REAL *pc)
01485 {
01486   INEXACT REAL acx, acy, bcx, bcy;
01487   REAL acxtail, acytail;
01488   REAL bcxtail, bcytail;
01489   REAL negate, negatetail;
01490   REAL axby[8], bxay[8];
01491   INEXACT REAL axby7, bxay7;
01492   REAL deter[16];
01493   int deterlen;
01494 
01495   INEXACT REAL bvirt;
01496   REAL avirt, bround, around;
01497   INEXACT REAL c;
01498   INEXACT REAL abig;
01499   REAL a0hi, a0lo, a1hi, a1lo, bhi, blo;
01500   REAL err1, err2, err3;
01501   INEXACT REAL _i, _j, _k, _l, _m, _n;
01502   REAL _0, _1, _2;
01503 
01504   Two_Diff(pa[0], pc[0], acx, acxtail);
01505   Two_Diff(pa[1], pc[1], acy, acytail);
01506   Two_Diff(pb[0], pc[0], bcx, bcxtail);
01507   Two_Diff(pb[1], pc[1], bcy, bcytail);
01508 
01509   Two_Two_Product(acx, acxtail, bcy, bcytail,
01510                   axby7, axby[6], axby[5], axby[4],
01511                   axby[3], axby[2], axby[1], axby[0]);
01512   axby[7] = axby7;
01513   negate = -acy;
01514   negatetail = -acytail;
01515   Two_Two_Product(bcx, bcxtail, negate, negatetail,
01516                   bxay7, bxay[6], bxay[5], bxay[4],
01517                   bxay[3], bxay[2], bxay[1], bxay[0]);
01518   bxay[7] = bxay7;
01519 
01520   deterlen = fast_expansion_sum_zeroelim(8, axby, 8, bxay, deter);
01521 
01522   return deter[deterlen - 1];
01523 }
01524 
01525 REAL orient2dadapt(REAL *pa, REAL *pb, REAL *pc, REAL detsum)
01526 {
01527   INEXACT REAL acx, acy, bcx, bcy;
01528   REAL acxtail, acytail, bcxtail, bcytail;
01529   INEXACT REAL detleft, detright;
01530   REAL detlefttail, detrighttail;
01531   REAL det, errbound;
01532   REAL B[4], C1[8], C2[12], D[16];
01533   INEXACT REAL B3;
01534   int C1length, C2length, Dlength;
01535   REAL u[4];
01536   INEXACT REAL u3;
01537   INEXACT REAL s1, t1;
01538   REAL s0, t0;
01539 
01540   INEXACT REAL bvirt;
01541   REAL avirt, bround, around;
01542   INEXACT REAL c;
01543   INEXACT REAL abig;
01544   REAL ahi, alo, bhi, blo;
01545   REAL err1, err2, err3;
01546   INEXACT REAL _i, _j;
01547   REAL _0;
01548 
01549   acx = (REAL) (pa[0] - pc[0]);
01550   bcx = (REAL) (pb[0] - pc[0]);
01551   acy = (REAL) (pa[1] - pc[1]);
01552   bcy = (REAL) (pb[1] - pc[1]);
01553 
01554   Two_Product(acx, bcy, detleft, detlefttail);
01555   Two_Product(acy, bcx, detright, detrighttail);
01556 
01557   Two_Two_Diff(detleft, detlefttail, detright, detrighttail,
01558                B3, B[2], B[1], B[0]);
01559   B[3] = B3;
01560 
01561   det = estimate(4, B);
01562   errbound = ccwerrboundB * detsum;
01563   if ((det >= errbound) || (-det >= errbound)) {
01564     return det;
01565   }
01566 
01567   Two_Diff_Tail(pa[0], pc[0], acx, acxtail);
01568   Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail);
01569   Two_Diff_Tail(pa[1], pc[1], acy, acytail);
01570   Two_Diff_Tail(pb[1], pc[1], bcy, bcytail);
01571 
01572   if ((acxtail == 0.0) && (acytail == 0.0)
01573       && (bcxtail == 0.0) && (bcytail == 0.0)) {
01574     return det;
01575   }
01576 
01577   errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det);
01578   det += (acx * bcytail + bcy * acxtail)
01579        - (acy * bcxtail + bcx * acytail);
01580   if ((det >= errbound) || (-det >= errbound)) {
01581     return det;
01582   }
01583 
01584   Two_Product(acxtail, bcy, s1, s0);
01585   Two_Product(acytail, bcx, t1, t0);
01586   Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
01587   u[3] = u3;
01588   C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1);
01589 
01590   Two_Product(acx, bcytail, s1, s0);
01591   Two_Product(acy, bcxtail, t1, t0);
01592   Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
01593   u[3] = u3;
01594   C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2);
01595 
01596   Two_Product(acxtail, bcytail, s1, s0);
01597   Two_Product(acytail, bcxtail, t1, t0);
01598   Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
01599   u[3] = u3;
01600   Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D);
01601 
01602   return(D[Dlength - 1]);
01603 }
01604 
01605 REAL orient2d(REAL *pa, REAL *pb, REAL *pc)
01606 {
01607   REAL detleft, detright, det;
01608   REAL detsum, errbound;
01609 
01610   detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]);
01611   detright = (pa[1] - pc[1]) * (pb[0] - pc[0]);
01612   det = detleft - detright;
01613 
01614   if (detleft > 0.0) {
01615     if (detright <= 0.0) {
01616       return det;
01617     } else {
01618       detsum = detleft + detright;
01619     }
01620   } else if (detleft < 0.0) {
01621     if (detright >= 0.0) {
01622       return det;
01623     } else {
01624       detsum = -detleft - detright;
01625     }
01626   } else {
01627     return det;
01628   }
01629 
01630   errbound = ccwerrboundA * detsum;
01631   if ((det >= errbound) || (-det >= errbound)) {
01632     return det;
01633   }
01634 
01635   return orient2dadapt(pa, pb, pc, detsum);
01636 }
01637 
01638 /*****************************************************************************/
01639 /*                                                                           */
01640 /*  orient3dfast()   Approximate 3D orientation test.  Nonrobust.            */
01641 /*  orient3dexact()   Exact 3D orientation test.  Robust.                    */
01642 /*  orient3dslow()   Another exact 3D orientation test.  Robust.             */
01643 /*  orient3d()   Adaptive exact 3D orientation test.  Robust.                */
01644 /*                                                                           */
01645 /*               Return a positive value if the point pd lies below the      */
01646 /*               plane passing through pa, pb, and pc; "below" is defined so */
01647 /*               that pa, pb, and pc appear in counterclockwise order when   */
01648 /*               viewed from above the plane.  Returns a negative value if   */
01649 /*               pd lies above the plane.  Returns zero if the points are    */
01650 /*               coplanar.  The result is also a rough approximation of six  */
01651 /*               times the signed volume of the tetrahedron defined by the   */
01652 /*               four points.                                                */
01653 /*                                                                           */
01654 /*  Only the first and last routine should be used; the middle two are for   */
01655 /*  timings.                                                                 */
01656 /*                                                                           */
01657 /*  The last three use exact arithmetic to ensure a correct answer.  The     */
01658 /*  result returned is the determinant of a matrix.  In orient3d() only,     */
01659 /*  this determinant is computed adaptively, in the sense that exact         */
01660 /*  arithmetic is used only to the degree it is needed to ensure that the    */
01661 /*  returned value has the correct sign.  Hence, orient3d() is usually quite */
01662 /*  fast, but will run more slowly when the input points are coplanar or     */
01663 /*  nearly so.                                                               */
01664 /*                                                                           */
01665 /*****************************************************************************/
01666 
01667 REAL orient3dfast(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
01668 {
01669   REAL adx, bdx, cdx;
01670   REAL ady, bdy, cdy;
01671   REAL adz, bdz, cdz;
01672 
01673   adx = pa[0] - pd[0];
01674   bdx = pb[0] - pd[0];
01675   cdx = pc[0] - pd[0];
01676   ady = pa[1] - pd[1];
01677   bdy = pb[1] - pd[1];
01678   cdy = pc[1] - pd[1];
01679   adz = pa[2] - pd[2];
01680   bdz = pb[2] - pd[2];
01681   cdz = pc[2] - pd[2];
01682 
01683   return adx * (bdy * cdz - bdz * cdy)
01684        + bdx * (cdy * adz - cdz * ady)
01685        + cdx * (ady * bdz - adz * bdy);
01686 }
01687 
01688 REAL orient3dexact(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
01689 {
01690   INEXACT REAL axby1, bxcy1, cxdy1, dxay1, axcy1, bxdy1;
01691   INEXACT REAL bxay1, cxby1, dxcy1, axdy1, cxay1, dxby1;
01692   REAL axby0, bxcy0, cxdy0, dxay0, axcy0, bxdy0;
01693   REAL bxay0, cxby0, dxcy0, axdy0, cxay0, dxby0;
01694   REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4];
01695   REAL temp8[8];
01696   int templen;
01697   REAL abc[12], bcd[12], cda[12], dab[12];
01698   int abclen, bcdlen, cdalen, dablen;
01699   REAL adet[24], bdet[24], cdet[24], ddet[24];
01700   int alen, blen, clen, dlen;
01701   REAL abdet[48], cddet[48];
01702   int ablen, cdlen;
01703   REAL deter[96];
01704   int deterlen;
01705   int i;
01706 
01707   INEXACT REAL bvirt;
01708   REAL avirt, bround, around;
01709   INEXACT REAL c;
01710   INEXACT REAL abig;
01711   REAL ahi, alo, bhi, blo;
01712   REAL err1, err2, err3;
01713   INEXACT REAL _i, _j;
01714   REAL _0;
01715 
01716   Two_Product(pa[0], pb[1], axby1, axby0);
01717   Two_Product(pb[0], pa[1], bxay1, bxay0);
01718   Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]);
01719 
01720   Two_Product(pb[0], pc[1], bxcy1, bxcy0);
01721   Two_Product(pc[0], pb[1], cxby1, cxby0);
01722   Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]);
01723 
01724   Two_Product(pc[0], pd[1], cxdy1, cxdy0);
01725   Two_Product(pd[0], pc[1], dxcy1, dxcy0);
01726   Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]);
01727 
01728   Two_Product(pd[0], pa[1], dxay1, dxay0);
01729   Two_Product(pa[0], pd[1], axdy1, axdy0);
01730   Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]);
01731 
01732   Two_Product(pa[0], pc[1], axcy1, axcy0);
01733   Two_Product(pc[0], pa[1], cxay1, cxay0);
01734   Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]);
01735 
01736   Two_Product(pb[0], pd[1], bxdy1, bxdy0);
01737   Two_Product(pd[0], pb[1], dxby1, dxby0);
01738   Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]);
01739 
01740   templen = fast_expansion_sum_zeroelim(4, cd, 4, da, temp8);
01741   cdalen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, cda);
01742   templen = fast_expansion_sum_zeroelim(4, da, 4, ab, temp8);
01743   dablen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, dab);
01744   for (i = 0; i < 4; i++) {
01745     bd[i] = -bd[i];
01746     ac[i] = -ac[i];
01747   }
01748   templen = fast_expansion_sum_zeroelim(4, ab, 4, bc, temp8);
01749   abclen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, abc);
01750   templen = fast_expansion_sum_zeroelim(4, bc, 4, cd, temp8);
01751   bcdlen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, bcd);
01752 
01753   alen = scale_expansion_zeroelim(bcdlen, bcd, pa[2], adet);
01754   blen = scale_expansion_zeroelim(cdalen, cda, -pb[2], bdet);
01755   clen = scale_expansion_zeroelim(dablen, dab, pc[2], cdet);
01756   dlen = scale_expansion_zeroelim(abclen, abc, -pd[2], ddet);
01757 
01758   ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
01759   cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
01760   deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, deter);
01761 
01762   return deter[deterlen - 1];
01763 }
01764 
01765 REAL orient3dslow(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
01766 {
01767   INEXACT REAL adx, ady, adz, bdx, bdy, bdz, cdx, cdy, cdz;
01768   REAL adxtail, adytail, adztail;
01769   REAL bdxtail, bdytail, bdztail;
01770   REAL cdxtail, cdytail, cdztail;
01771   REAL negate, negatetail;
01772   INEXACT REAL axby7, bxcy7, axcy7, bxay7, cxby7, cxay7;
01773   REAL axby[8], bxcy[8], axcy[8], bxay[8], cxby[8], cxay[8];
01774   REAL temp16[16], temp32[32], temp32t[32];
01775   int temp16len, temp32len, temp32tlen;
01776   REAL adet[64], bdet[64], cdet[64];
01777   int alen, blen, clen;
01778   REAL abdet[128];
01779   int ablen;
01780   REAL deter[192];
01781   int deterlen;
01782 
01783   INEXACT REAL bvirt;
01784   REAL avirt, bround, around;
01785   INEXACT REAL c;
01786   INEXACT REAL abig;
01787   REAL a0hi, a0lo, a1hi, a1lo, bhi, blo;
01788   REAL err1, err2, err3;
01789   INEXACT REAL _i, _j, _k, _l, _m, _n;
01790   REAL _0, _1, _2;
01791 
01792   Two_Diff(pa[0], pd[0], adx, adxtail);
01793   Two_Diff(pa[1], pd[1], ady, adytail);
01794   Two_Diff(pa[2], pd[2], adz, adztail);
01795   Two_Diff(pb[0], pd[0], bdx, bdxtail);
01796   Two_Diff(pb[1], pd[1], bdy, bdytail);
01797   Two_Diff(pb[2], pd[2], bdz, bdztail);
01798   Two_Diff(pc[0], pd[0], cdx, cdxtail);
01799   Two_Diff(pc[1], pd[1], cdy, cdytail);
01800   Two_Diff(pc[2], pd[2], cdz, cdztail);
01801 
01802   Two_Two_Product(adx, adxtail, bdy, bdytail,
01803                   axby7, axby[6], axby[5], axby[4],
01804                   axby[3], axby[2], axby[1], axby[0]);
01805   axby[7] = axby7;
01806   negate = -ady;
01807   negatetail = -adytail;
01808   Two_Two_Product(bdx, bdxtail, negate, negatetail,
01809                   bxay7, bxay[6], bxay[5], bxay[4],
01810                   bxay[3], bxay[2], bxay[1], bxay[0]);
01811   bxay[7] = bxay7;
01812   Two_Two_Product(bdx, bdxtail, cdy, cdytail,
01813                   bxcy7, bxcy[6], bxcy[5], bxcy[4],
01814                   bxcy[3], bxcy[2], bxcy[1], bxcy[0]);
01815   bxcy[7] = bxcy7;
01816   negate = -bdy;
01817   negatetail = -bdytail;
01818   Two_Two_Product(cdx, cdxtail, negate, negatetail,
01819                   cxby7, cxby[6], cxby[5], cxby[4],
01820                   cxby[3], cxby[2], cxby[1], cxby[0]);
01821   cxby[7] = cxby7;
01822   Two_Two_Product(cdx, cdxtail, ady, adytail,
01823                   cxay7, cxay[6], cxay[5], cxay[4],
01824                   cxay[3], cxay[2], cxay[1], cxay[0]);
01825   cxay[7] = cxay7;
01826   negate = -cdy;
01827   negatetail = -cdytail;
01828   Two_Two_Product(adx, adxtail, negate, negatetail,
01829                   axcy7, axcy[6], axcy[5], axcy[4],
01830                   axcy[3], axcy[2], axcy[1], axcy[0]);
01831   axcy[7] = axcy7;
01832 
01833   temp16len = fast_expansion_sum_zeroelim(8, bxcy, 8, cxby, temp16);
01834   temp32len = scale_expansion_zeroelim(temp16len, temp16, adz, temp32);
01835   temp32tlen = scale_expansion_zeroelim(temp16len, temp16, adztail, temp32t);
01836   alen = fast_expansion_sum_zeroelim(temp32len, temp32, temp32tlen, temp32t,
01837                                      adet);
01838 
01839   temp16len = fast_expansion_sum_zeroelim(8, cxay, 8, axcy, temp16);
01840   temp32len = scale_expansion_zeroelim(temp16len, temp16, bdz, temp32);
01841   temp32tlen = scale_expansion_zeroelim(temp16len, temp16, bdztail, temp32t);
01842   blen = fast_expansion_sum_zeroelim(temp32len, temp32, temp32tlen, temp32t,
01843                                      bdet);
01844 
01845   temp16len = fast_expansion_sum_zeroelim(8, axby, 8, bxay, temp16);
01846   temp32len = scale_expansion_zeroelim(temp16len, temp16, cdz, temp32);
01847   temp32tlen = scale_expansion_zeroelim(temp16len, temp16, cdztail, temp32t);
01848   clen = fast_expansion_sum_zeroelim(temp32len, temp32, temp32tlen, temp32t,
01849                                      cdet);
01850 
01851   ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
01852   deterlen = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, deter);
01853 
01854   return deter[deterlen - 1];
01855 }
01856 
01857 REAL orient3dadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL permanent)
01858 {
01859   INEXACT REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
01860   REAL det, errbound;
01861 
01862   INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
01863   REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
01864   REAL bc[4], ca[4], ab[4];
01865   INEXACT REAL bc3, ca3, ab3;
01866   REAL adet[8], bdet[8], cdet[8];
01867   int alen, blen, clen;
01868   REAL abdet[16];
01869   int ablen;
01870   REAL *finnow, *finother, *finswap;
01871   REAL fin1[192], fin2[192];
01872   int finlength;
01873 
01875   // To avoid uninitialized warnings reported by valgrind.
01876   int i;
01877   for (i = 0; i < 8; i++) {
01878     adet[i] = bdet[i] = cdet[i] = 0.0;
01879   }
01880   for (i = 0; i < 16; i++) {
01881     abdet[i] = 0.0;
01882   }
01884 
01885   REAL adxtail, bdxtail, cdxtail;
01886   REAL adytail, bdytail, cdytail;
01887   REAL adztail, bdztail, cdztail;
01888   INEXACT REAL at_blarge, at_clarge;
01889   INEXACT REAL bt_clarge, bt_alarge;
01890   INEXACT REAL ct_alarge, ct_blarge;
01891   REAL at_b[4], at_c[4], bt_c[4], bt_a[4], ct_a[4], ct_b[4];
01892   int at_blen, at_clen, bt_clen, bt_alen, ct_alen, ct_blen;
01893   INEXACT REAL bdxt_cdy1, cdxt_bdy1, cdxt_ady1;
01894   INEXACT REAL adxt_cdy1, adxt_bdy1, bdxt_ady1;
01895   REAL bdxt_cdy0, cdxt_bdy0, cdxt_ady0;
01896   REAL adxt_cdy0, adxt_bdy0, bdxt_ady0;
01897   INEXACT REAL bdyt_cdx1, cdyt_bdx1, cdyt_adx1;
01898   INEXACT REAL adyt_cdx1, adyt_bdx1, bdyt_adx1;
01899   REAL bdyt_cdx0, cdyt_bdx0, cdyt_adx0;
01900   REAL adyt_cdx0, adyt_bdx0, bdyt_adx0;
01901   REAL bct[8], cat[8], abt[8];
01902   int bctlen, catlen, abtlen;
01903   INEXACT REAL bdxt_cdyt1, cdxt_bdyt1, cdxt_adyt1;
01904   INEXACT REAL adxt_cdyt1, adxt_bdyt1, bdxt_adyt1;
01905   REAL bdxt_cdyt0, cdxt_bdyt0, cdxt_adyt0;
01906   REAL adxt_cdyt0, adxt_bdyt0, bdxt_adyt0;
01907   REAL u[4], v[12], w[16];
01908   INEXACT REAL u3;
01909   int vlength, wlength;
01910   REAL negate;
01911 
01912   INEXACT REAL bvirt;
01913   REAL avirt, bround, around;
01914   INEXACT REAL c;
01915   INEXACT REAL abig;
01916   REAL ahi, alo, bhi, blo;
01917   REAL err1, err2, err3;
01918   INEXACT REAL _i, _j, _k;
01919   REAL _0;
01920 
01921   adx = (REAL) (pa[0] - pd[0]);
01922   bdx = (REAL) (pb[0] - pd[0]);
01923   cdx = (REAL) (pc[0] - pd[0]);
01924   ady = (REAL) (pa[1] - pd[1]);
01925   bdy = (REAL) (pb[1] - pd[1]);
01926   cdy = (REAL) (pc[1] - pd[1]);
01927   adz = (REAL) (pa[2] - pd[2]);
01928   bdz = (REAL) (pb[2] - pd[2]);
01929   cdz = (REAL) (pc[2] - pd[2]);
01930 
01931   Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
01932   Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
01933   Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
01934   bc[3] = bc3;
01935   alen = scale_expansion_zeroelim(4, bc, adz, adet);
01936 
01937   Two_Product(cdx, ady, cdxady1, cdxady0);
01938   Two_Product(adx, cdy, adxcdy1, adxcdy0);
01939   Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
01940   ca[3] = ca3;
01941   blen = scale_expansion_zeroelim(4, ca, bdz, bdet);
01942 
01943   Two_Product(adx, bdy, adxbdy1, adxbdy0);
01944   Two_Product(bdx, ady, bdxady1, bdxady0);
01945   Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
01946   ab[3] = ab3;
01947   clen = scale_expansion_zeroelim(4, ab, cdz, cdet);
01948 
01949   ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
01950   finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
01951 
01952   det = estimate(finlength, fin1);
01953   errbound = o3derrboundB * permanent;
01954   if ((det >= errbound) || (-det >= errbound)) {
01955     return det;
01956   }
01957 
01958   Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
01959   Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
01960   Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
01961   Two_Diff_Tail(pa[1], pd[1], ady, adytail);
01962   Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
01963   Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
01964   Two_Diff_Tail(pa[2], pd[2], adz, adztail);
01965   Two_Diff_Tail(pb[2], pd[2], bdz, bdztail);
01966   Two_Diff_Tail(pc[2], pd[2], cdz, cdztail);
01967 
01968   if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
01969       && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)
01970       && (adztail == 0.0) && (bdztail == 0.0) && (cdztail == 0.0)) {
01971     return det;
01972   }
01973 
01974   errbound = o3derrboundC * permanent + resulterrbound * Absolute(det);
01975   det += (adz * ((bdx * cdytail + cdy * bdxtail)
01976                  - (bdy * cdxtail + cdx * bdytail))
01977           + adztail * (bdx * cdy - bdy * cdx))
01978        + (bdz * ((cdx * adytail + ady * cdxtail)
01979                  - (cdy * adxtail + adx * cdytail))
01980           + bdztail * (cdx * ady - cdy * adx))
01981        + (cdz * ((adx * bdytail + bdy * adxtail)
01982                  - (ady * bdxtail + bdx * adytail))
01983           + cdztail * (adx * bdy - ady * bdx));
01984   if ((det >= errbound) || (-det >= errbound)) {
01985     return det;
01986   }
01987 
01988   finnow = fin1;
01989   finother = fin2;
01990 
01991   if (adxtail == 0.0) {
01992     if (adytail == 0.0) {
01993       at_b[0] = 0.0;
01994       at_blen = 1;
01995       at_c[0] = 0.0;
01996       at_clen = 1;
01997     } else {
01998       negate = -adytail;
01999       Two_Product(negate, bdx, at_blarge, at_b[0]);
02000       at_b[1] = at_blarge;
02001       at_blen = 2;
02002       Two_Product(adytail, cdx, at_clarge, at_c[0]);
02003       at_c[1] = at_clarge;
02004       at_clen = 2;
02005     }
02006   } else {
02007     if (adytail == 0.0) {
02008       Two_Product(adxtail, bdy, at_blarge, at_b[0]);
02009       at_b[1] = at_blarge;
02010       at_blen = 2;
02011       negate = -adxtail;
02012       Two_Product(negate, cdy, at_clarge, at_c[0]);
02013       at_c[1] = at_clarge;
02014       at_clen = 2;
02015     } else {
02016       Two_Product(adxtail, bdy, adxt_bdy1, adxt_bdy0);
02017       Two_Product(adytail, bdx, adyt_bdx1, adyt_bdx0);
02018       Two_Two_Diff(adxt_bdy1, adxt_bdy0, adyt_bdx1, adyt_bdx0,
02019                    at_blarge, at_b[2], at_b[1], at_b[0]);
02020       at_b[3] = at_blarge;
02021       at_blen = 4;
02022       Two_Product(adytail, cdx, adyt_cdx1, adyt_cdx0);
02023       Two_Product(adxtail, cdy, adxt_cdy1, adxt_cdy0);
02024       Two_Two_Diff(adyt_cdx1, adyt_cdx0, adxt_cdy1, adxt_cdy0,
02025                    at_clarge, at_c[2], at_c[1], at_c[0]);
02026       at_c[3] = at_clarge;
02027       at_clen = 4;
02028     }
02029   }
02030   if (bdxtail == 0.0) {
02031     if (bdytail == 0.0) {
02032       bt_c[0] = 0.0;
02033       bt_clen = 1;
02034       bt_a[0] = 0.0;
02035       bt_alen = 1;
02036     } else {
02037       negate = -bdytail;
02038       Two_Product(negate, cdx, bt_clarge, bt_c[0]);
02039       bt_c[1] = bt_clarge;
02040       bt_clen = 2;
02041       Two_Product(bdytail, adx, bt_alarge, bt_a[0]);
02042       bt_a[1] = bt_alarge;
02043       bt_alen = 2;
02044     }
02045   } else {
02046     if (bdytail == 0.0) {
02047       Two_Product(bdxtail, cdy, bt_clarge, bt_c[0]);
02048       bt_c[1] = bt_clarge;
02049       bt_clen = 2;
02050       negate = -bdxtail;
02051       Two_Product(negate, ady, bt_alarge, bt_a[0]);
02052       bt_a[1] = bt_alarge;
02053       bt_alen = 2;
02054     } else {
02055       Two_Product(bdxtail, cdy, bdxt_cdy1, bdxt_cdy0);
02056       Two_Product(bdytail, cdx, bdyt_cdx1, bdyt_cdx0);
02057       Two_Two_Diff(bdxt_cdy1, bdxt_cdy0, bdyt_cdx1, bdyt_cdx0,
02058                    bt_clarge, bt_c[2], bt_c[1], bt_c[0]);
02059       bt_c[3] = bt_clarge;
02060       bt_clen = 4;
02061       Two_Product(bdytail, adx, bdyt_adx1, bdyt_adx0);
02062       Two_Product(bdxtail, ady, bdxt_ady1, bdxt_ady0);
02063       Two_Two_Diff(bdyt_adx1, bdyt_adx0, bdxt_ady1, bdxt_ady0,
02064                   bt_alarge, bt_a[2], bt_a[1], bt_a[0]);
02065       bt_a[3] = bt_alarge;
02066       bt_alen = 4;
02067     }
02068   }
02069   if (cdxtail == 0.0) {
02070     if (cdytail == 0.0) {
02071       ct_a[0] = 0.0;
02072       ct_alen = 1;
02073       ct_b[0] = 0.0;
02074       ct_blen = 1;
02075     } else {
02076       negate = -cdytail;
02077       Two_Product(negate, adx, ct_alarge, ct_a[0]);
02078       ct_a[1] = ct_alarge;
02079       ct_alen = 2;
02080       Two_Product(cdytail, bdx, ct_blarge, ct_b[0]);
02081       ct_b[1] = ct_blarge;
02082       ct_blen = 2;
02083     }
02084   } else {
02085     if (cdytail == 0.0) {
02086       Two_Product(cdxtail, ady, ct_alarge, ct_a[0]);
02087       ct_a[1] = ct_alarge;
02088       ct_alen = 2;
02089       negate = -cdxtail;
02090       Two_Product(negate, bdy, ct_blarge, ct_b[0]);
02091       ct_b[1] = ct_blarge;
02092       ct_blen = 2;
02093     } else {
02094       Two_Product(cdxtail, ady, cdxt_ady1, cdxt_ady0);
02095       Two_Product(cdytail, adx, cdyt_adx1, cdyt_adx0);
02096       Two_Two_Diff(cdxt_ady1, cdxt_ady0, cdyt_adx1, cdyt_adx0,
02097                    ct_alarge, ct_a[2], ct_a[1], ct_a[0]);
02098       ct_a[3] = ct_alarge;
02099       ct_alen = 4;
02100       Two_Product(cdytail, bdx, cdyt_bdx1, cdyt_bdx0);
02101       Two_Product(cdxtail, bdy, cdxt_bdy1, cdxt_bdy0);
02102       Two_Two_Diff(cdyt_bdx1, cdyt_bdx0, cdxt_bdy1, cdxt_bdy0,
02103                    ct_blarge, ct_b[2], ct_b[1], ct_b[0]);
02104       ct_b[3] = ct_blarge;
02105       ct_blen = 4;
02106     }
02107   }
02108 
02109   bctlen = fast_expansion_sum_zeroelim(bt_clen, bt_c, ct_blen, ct_b, bct);
02110   wlength = scale_expansion_zeroelim(bctlen, bct, adz, w);
02111   finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
02112                                           finother);
02113   finswap = finnow; finnow = finother; finother = finswap;
02114 
02115   catlen = fast_expansion_sum_zeroelim(ct_alen, ct_a, at_clen, at_c, cat);
02116   wlength = scale_expansion_zeroelim(catlen, cat, bdz, w);
02117   finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
02118                                           finother);
02119   finswap = finnow; finnow = finother; finother = finswap;
02120 
02121   abtlen = fast_expansion_sum_zeroelim(at_blen, at_b, bt_alen, bt_a, abt);
02122   wlength = scale_expansion_zeroelim(abtlen, abt, cdz, w);
02123   finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
02124                                           finother);
02125   finswap = finnow; finnow = finother; finother = finswap;
02126 
02127   if (adztail != 0.0) {
02128     vlength = scale_expansion_zeroelim(4, bc, adztail, v);
02129     finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
02130                                             finother);
02131     finswap = finnow; finnow = finother; finother = finswap;
02132   }
02133   if (bdztail != 0.0) {
02134     vlength = scale_expansion_zeroelim(4, ca, bdztail, v);
02135     finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
02136                                             finother);
02137     finswap = finnow; finnow = finother; finother = finswap;
02138   }
02139   if (cdztail != 0.0) {
02140     vlength = scale_expansion_zeroelim(4, ab, cdztail, v);
02141     finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
02142                                             finother);
02143     finswap = finnow; finnow = finother; finother = finswap;
02144   }
02145 
02146   if (adxtail != 0.0) {
02147     if (bdytail != 0.0) {
02148       Two_Product(adxtail, bdytail, adxt_bdyt1, adxt_bdyt0);
02149       Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdz, u3, u[2], u[1], u[0]);
02150       u[3] = u3;
02151       finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
02152                                               finother);
02153       finswap = finnow; finnow = finother; finother = finswap;
02154       if (cdztail != 0.0) {
02155         Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdztail, u3, u[2], u[1], u[0]);
02156         u[3] = u3;
02157         finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
02158                                                 finother);
02159         finswap = finnow; finnow = finother; finother = finswap;
02160       }
02161     }
02162     if (cdytail != 0.0) {
02163       negate = -adxtail;
02164       Two_Product(negate, cdytail, adxt_cdyt1, adxt_cdyt0);
02165       Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdz, u3, u[2], u[1], u[0]);
02166       u[3] = u3;
02167       finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
02168                                               finother);
02169       finswap = finnow; finnow = finother; finother = finswap;
02170       if (bdztail != 0.0) {
02171         Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdztail, u3, u[2], u[1], u[0]);
02172         u[3] = u3;
02173         finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
02174                                                 finother);
02175         finswap = finnow; finnow = finother; finother = finswap;
02176       }
02177     }
02178   }
02179   if (bdxtail != 0.0) {
02180     if (cdytail != 0.0) {
02181       Two_Product(bdxtail, cdytail, bdxt_cdyt1, bdxt_cdyt0);
02182       Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adz, u3, u[2], u[1], u[0]);
02183       u[3] = u3;
02184       finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
02185                                               finother);
02186       finswap = finnow; finnow = finother; finother = finswap;
02187       if (adztail != 0.0) {
02188         Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adztail, u3, u[2], u[1], u[0]);
02189         u[3] = u3;
02190         finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
02191                                                 finother);
02192         finswap = finnow; finnow = finother; finother = finswap;
02193       }
02194     }
02195     if (adytail != 0.0) {
02196       negate = -bdxtail;
02197       Two_Product(negate, adytail, bdxt_adyt1, bdxt_adyt0);
02198       Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdz, u3, u[2], u[1], u[0]);
02199       u[3] = u3;
02200       finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
02201                                               finother);
02202       finswap = finnow; finnow = finother; finother = finswap;
02203       if (cdztail != 0.0) {
02204         Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdztail, u3, u[2], u[1], u[0]);
02205         u[3] = u3;
02206         finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
02207                                                 finother);
02208         finswap = finnow; finnow = finother; finother = finswap;
02209       }
02210     }
02211   }
02212   if (cdxtail != 0.0) {
02213     if (adytail != 0.0) {
02214       Two_Product(cdxtail, adytail, cdxt_adyt1, cdxt_adyt0);
02215       Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdz, u3, u[2], u[1], u[0]);
02216       u[3] = u3;
02217       finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
02218                                               finother);
02219       finswap = finnow; finnow = finother; finother = finswap;
02220       if (bdztail != 0.0) {
02221         Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdztail, u3, u[2], u[1], u[0]);
02222         u[3] = u3;
02223         finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
02224                                                 finother);
02225         finswap = finnow; finnow = finother; finother = finswap;
02226       }
02227     }
02228     if (bdytail != 0.0) {
02229       negate = -cdxtail;
02230       Two_Product(negate, bdytail, cdxt_bdyt1, cdxt_bdyt0);
02231       Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adz, u3, u[2], u[1], u[0]);
02232       u[3] = u3;
02233       finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
02234                                               finother);
02235       finswap = finnow; finnow = finother; finother = finswap;
02236       if (adztail != 0.0) {
02237         Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adztail, u3, u[2], u[1], u[0]);
02238         u[3] = u3;
02239         finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
02240                                                 finother);
02241         finswap = finnow; finnow = finother; finother = finswap;
02242       }
02243     }
02244   }
02245 
02246   if (adztail != 0.0) {
02247     wlength = scale_expansion_zeroelim(bctlen, bct, adztail, w);
02248     finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
02249                                             finother);
02250     finswap = finnow; finnow = finother; finother = finswap;
02251   }
02252   if (bdztail != 0.0) {
02253     wlength = scale_expansion_zeroelim(catlen, cat, bdztail, w);
02254     finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
02255                                             finother);
02256     finswap = finnow; finnow = finother; finother = finswap;
02257   }
02258   if (cdztail != 0.0) {
02259     wlength = scale_expansion_zeroelim(abtlen, abt, cdztail, w);
02260     finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
02261                                             finother);
02262     finswap = finnow; finnow = finother; finother = finswap;
02263   }
02264 
02265   return finnow[finlength - 1];
02266 }
02267 
02268 REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
02269 {
02270   REAL adx, bdx, cdx, ady, bdy, cdy, adz, bdz, cdz;
02271   REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
02272   REAL det;
02273   REAL permanent, errbound;
02274 
02275   adx = pa[0] - pd[0];
02276   bdx = pb[0] - pd[0];
02277   cdx = pc[0] - pd[0];
02278   ady = pa[1] - pd[1];
02279   bdy = pb[1] - pd[1];
02280   cdy = pc[1] - pd[1];
02281   adz = pa[2] - pd[2];
02282   bdz = pb[2] - pd[2];
02283   cdz = pc[2] - pd[2];
02284 
02285   bdxcdy = bdx * cdy;
02286   cdxbdy = cdx * bdy;
02287 
02288   cdxady = cdx * ady;
02289   adxcdy = adx * cdy;
02290 
02291   adxbdy = adx * bdy;
02292   bdxady = bdx * ady;
02293 
02294   det = adz * (bdxcdy - cdxbdy)
02295       + bdz * (cdxady - adxcdy)
02296       + cdz * (adxbdy - bdxady);
02297 
02298   permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * Absolute(adz)
02299             + (Absolute(cdxady) + Absolute(adxcdy)) * Absolute(bdz)
02300             + (Absolute(adxbdy) + Absolute(bdxady)) * Absolute(cdz);
02301   errbound = o3derrboundA * permanent;
02302   if ((det > errbound) || (-det > errbound)) {
02303     return det;
02304   }
02305 
02306   return orient3dadapt(pa, pb, pc, pd, permanent);
02307 }
02308 
02309 /*****************************************************************************/
02310 /*                                                                           */
02311 /*  incirclefast()   Approximate 2D incircle test.  Nonrobust.               */
02312 /*  incircleexact()   Exact 2D incircle test.  Robust.                       */
02313 /*  incircleslow()   Another exact 2D incircle test.  Robust.                */
02314 /*  incircle()   Adaptive exact 2D incircle test.  Robust.                   */
02315 /*                                                                           */
02316 /*               Return a positive value if the point pd lies inside the     */
02317 /*               circle passing through pa, pb, and pc; a negative value if  */
02318 /*               it lies outside; and zero if the four points are cocircular.*/
02319 /*               The points pa, pb, and pc must be in counterclockwise       */
02320 /*               order, or the sign of the result will be reversed.          */
02321 /*                                                                           */
02322 /*  Only the first and last routine should be used; the middle two are for   */
02323 /*  timings.                                                                 */
02324 /*                                                                           */
02325 /*  The last three use exact arithmetic to ensure a correct answer.  The     */
02326 /*  result returned is the determinant of a matrix.  In incircle() only,     */
02327 /*  this determinant is computed adaptively, in the sense that exact         */
02328 /*  arithmetic is used only to the degree it is needed to ensure that the    */
02329 /*  returned value has the correct sign.  Hence, incircle() is usually quite */
02330 /*  fast, but will run more slowly when the input points are cocircular or   */
02331 /*  nearly so.                                                               */
02332 /*                                                                           */
02333 /*****************************************************************************/
02334 
02335 REAL incirclefast(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
02336 {
02337   REAL adx, ady, bdx, bdy, cdx, cdy;
02338   REAL abdet, bcdet, cadet;
02339   REAL alift, blift, clift;
02340 
02341   adx = pa[0] - pd[0];
02342   ady = pa[1] - pd[1];
02343   bdx = pb[0] - pd[0];
02344   bdy = pb[1] - pd[1];
02345   cdx = pc[0] - pd[0];
02346   cdy = pc[1] - pd[1];
02347 
02348   abdet = adx * bdy - bdx * ady;
02349   bcdet = bdx * cdy - cdx * bdy;
02350   cadet = cdx * ady - adx * cdy;
02351   alift = adx * adx + ady * ady;
02352   blift = bdx * bdx + bdy * bdy;
02353   clift = cdx * cdx + cdy * cdy;
02354 
02355   return alift * bcdet + blift * cadet + clift * abdet;
02356 }
02357 
02358 REAL incircleexact(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
02359 {
02360   INEXACT REAL axby1, bxcy1, cxdy1, dxay1, axcy1, bxdy1;
02361   INEXACT REAL bxay1, cxby1, dxcy1, axdy1, cxay1, dxby1;
02362   REAL axby0, bxcy0, cxdy0, dxay0, axcy0, bxdy0;
02363   REAL bxay0, cxby0, dxcy0, axdy0, cxay0, dxby0;
02364   REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4];
02365   REAL temp8[8];
02366   int templen;
02367   REAL abc[12], bcd[12], cda[12], dab[12];
02368   int abclen, bcdlen, cdalen, dablen;
02369   REAL det24x[24], det24y[24], det48x[48], det48y[48];
02370   int xlen, ylen;
02371   REAL adet[96], bdet[96], cdet[96], ddet[96];
02372   int alen, blen, clen, dlen;
02373   REAL abdet[192], cddet[192];
02374   int ablen, cdlen;
02375   REAL deter[384];
02376   int deterlen;
02377   int i;
02378 
02379   INEXACT REAL bvirt;
02380   REAL avirt, bround, around;
02381   INEXACT REAL c;
02382   INEXACT REAL abig;
02383   REAL ahi, alo, bhi, blo;
02384   REAL err1, err2, err3;
02385   INEXACT REAL _i, _j;
02386   REAL _0;
02387 
02388   Two_Product(pa[0], pb[1], axby1, axby0);
02389   Two_Product(pb[0], pa[1], bxay1, bxay0);
02390   Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]);
02391 
02392   Two_Product(pb[0], pc[1], bxcy1, bxcy0);
02393   Two_Product(pc[0], pb[1], cxby1, cxby0);
02394   Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]);
02395 
02396   Two_Product(pc[0], pd[1], cxdy1, cxdy0);
02397   Two_Product(pd[0], pc[1], dxcy1, dxcy0);
02398   Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]);
02399 
02400   Two_Product(pd[0], pa[1], dxay1, dxay0);
02401   Two_Product(pa[0], pd[1], axdy1, axdy0);
02402   Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]);
02403 
02404   Two_Product(pa[0], pc[1], axcy1, axcy0);
02405   Two_Product(pc[0], pa[1], cxay1, cxay0);
02406   Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]);
02407 
02408   Two_Product(pb[0], pd[1], bxdy1, bxdy0);
02409   Two_Product(pd[0], pb[1], dxby1, dxby0);
02410   Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]);
02411 
02412   templen = fast_expansion_sum_zeroelim(4, cd, 4, da, temp8);
02413   cdalen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, cda);
02414   templen = fast_expansion_sum_zeroelim(4, da, 4, ab, temp8);
02415   dablen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, dab);
02416   for (i = 0; i < 4; i++) {
02417     bd[i] = -bd[i];
02418     ac[i] = -ac[i];
02419   }
02420   templen = fast_expansion_sum_zeroelim(4, ab, 4, bc, temp8);
02421   abclen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, abc);
02422   templen = fast_expansion_sum_zeroelim(4, bc, 4, cd, temp8);
02423   bcdlen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, bcd);
02424 
02425   xlen = scale_expansion_zeroelim(bcdlen, bcd, pa[0], det24x);
02426   xlen = scale_expansion_zeroelim(xlen, det24x, pa[0], det48x);
02427   ylen = scale_expansion_zeroelim(bcdlen, bcd, pa[1], det24y);
02428   ylen = scale_expansion_zeroelim(ylen, det24y, pa[1], det48y);
02429   alen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, adet);
02430 
02431   xlen = scale_expansion_zeroelim(cdalen, cda, pb[0], det24x);
02432   xlen = scale_expansion_zeroelim(xlen, det24x, -pb[0], det48x);
02433   ylen = scale_expansion_zeroelim(cdalen, cda, pb[1], det24y);
02434   ylen = scale_expansion_zeroelim(ylen, det24y, -pb[1], det48y);
02435   blen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, bdet);
02436 
02437   xlen = scale_expansion_zeroelim(dablen, dab, pc[0], det24x);
02438   xlen = scale_expansion_zeroelim(xlen, det24x, pc[0], det48x);
02439   ylen = scale_expansion_zeroelim(dablen, dab, pc[1], det24y);
02440   ylen = scale_expansion_zeroelim(ylen, det24y, pc[1], det48y);
02441   clen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, cdet);
02442 
02443   xlen = scale_expansion_zeroelim(abclen, abc, pd[0], det24x);
02444   xlen = scale_expansion_zeroelim(xlen, det24x, -pd[0], det48x);
02445   ylen = scale_expansion_zeroelim(abclen, abc, pd[1], det24y);
02446   ylen = scale_expansion_zeroelim(ylen, det24y, -pd[1], det48y);
02447   dlen = fast_expansion_sum_zeroelim(xlen, det48x, ylen, det48y, ddet);
02448 
02449   ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
02450   cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
02451   deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, deter);
02452 
02453   return deter[deterlen - 1];
02454 }
02455 
02456 REAL incircleslow(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
02457 {
02458   INEXACT REAL adx, bdx, cdx, ady, bdy, cdy;
02459   REAL adxtail, bdxtail, cdxtail;
02460   REAL adytail, bdytail, cdytail;
02461   REAL negate, negatetail;
02462   INEXACT REAL axby7, bxcy7, axcy7, bxay7, cxby7, cxay7;
02463   REAL axby[8], bxcy[8], axcy[8], bxay[8], cxby[8], cxay[8];
02464   REAL temp16[16];
02465   int temp16len;
02466   REAL detx[32], detxx[64], detxt[32], detxxt[64], detxtxt[64];
02467   int xlen, xxlen, xtlen, xxtlen, xtxtlen;
02468   REAL x1[128], x2[192];
02469   int x1len, x2len;
02470   REAL dety[32], detyy[64], detyt[32], detyyt[64], detytyt[64];
02471   int ylen, yylen, ytlen, yytlen, ytytlen;
02472   REAL y1[128], y2[192];
02473   int y1len, y2len;
02474   REAL adet[384], bdet[384], cdet[384], abdet[768], deter[1152];
02475   int alen, blen, clen, ablen, deterlen;
02476   int i;
02477 
02478   INEXACT REAL bvirt;
02479   REAL avirt, bround, around;
02480   INEXACT REAL c;
02481   INEXACT REAL abig;
02482   REAL a0hi, a0lo, a1hi, a1lo, bhi, blo;
02483   REAL err1, err2, err3;
02484   INEXACT REAL _i, _j, _k, _l, _m, _n;
02485   REAL _0, _1, _2;
02486 
02487   Two_Diff(pa[0], pd[0], adx, adxtail);
02488   Two_Diff(pa[1], pd[1], ady, adytail);
02489   Two_Diff(pb[0], pd[0], bdx, bdxtail);
02490   Two_Diff(pb[1], pd[1], bdy, bdytail);
02491   Two_Diff(pc[0], pd[0], cdx, cdxtail);
02492   Two_Diff(pc[1], pd[1], cdy, cdytail);
02493 
02494   Two_Two_Product(adx, adxtail, bdy, bdytail,
02495                   axby7, axby[6], axby[5], axby[4],
02496                   axby[3], axby[2], axby[1], axby[0]);
02497   axby[7] = axby7;
02498   negate = -ady;
02499   negatetail = -adytail;
02500   Two_Two_Product(bdx, bdxtail, negate, negatetail,
02501                   bxay7, bxay[6], bxay[5], bxay[4],
02502                   bxay[3], bxay[2], bxay[1], bxay[0]);
02503   bxay[7] = bxay7;
02504   Two_Two_Product(bdx, bdxtail, cdy, cdytail,
02505                   bxcy7, bxcy[6], bxcy[5], bxcy[4],
02506                   bxcy[3], bxcy[2], bxcy[1], bxcy[0]);
02507   bxcy[7] = bxcy7;
02508   negate = -bdy;
02509   negatetail = -bdytail;
02510   Two_Two_Product(cdx, cdxtail, negate, negatetail,
02511                   cxby7, cxby[6], cxby[5], cxby[4],
02512                   cxby[3], cxby[2], cxby[1], cxby[0]);
02513   cxby[7] = cxby7;
02514   Two_Two_Product(cdx, cdxtail, ady, adytail,
02515                   cxay7, cxay[6], cxay[5], cxay[4],
02516                   cxay[3], cxay[2], cxay[1], cxay[0]);
02517   cxay[7] = cxay7;
02518   negate = -cdy;
02519   negatetail = -cdytail;
02520   Two_Two_Product(adx, adxtail, negate, negatetail,
02521                   axcy7, axcy[6], axcy[5], axcy[4],
02522                   axcy[3], axcy[2], axcy[1], axcy[0]);
02523   axcy[7] = axcy7;
02524 
02525 
02526   temp16len = fast_expansion_sum_zeroelim(8, bxcy, 8, cxby, temp16);
02527 
02528   xlen = scale_expansion_zeroelim(temp16len, temp16, adx, detx);
02529   xxlen = scale_expansion_zeroelim(xlen, detx, adx, detxx);
02530   xtlen = scale_expansion_zeroelim(temp16len, temp16, adxtail, detxt);
02531   xxtlen = scale_expansion_zeroelim(xtlen, detxt, adx, detxxt);
02532   for (i = 0; i < xxtlen; i++) {
02533     detxxt[i] *= 2.0;
02534   }
02535   xtxtlen = scale_expansion_zeroelim(xtlen, detxt, adxtail, detxtxt);
02536   x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
02537   x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
02538 
02539   ylen = scale_expansion_zeroelim(temp16len, temp16, ady, dety);
02540   yylen = scale_expansion_zeroelim(ylen, dety, ady, detyy);
02541   ytlen = scale_expansion_zeroelim(temp16len, temp16, adytail, detyt);
02542   yytlen = scale_expansion_zeroelim(ytlen, detyt, ady, detyyt);
02543   for (i = 0; i < yytlen; i++) {
02544     detyyt[i] *= 2.0;
02545   }
02546   ytytlen = scale_expansion_zeroelim(ytlen, detyt, adytail, detytyt);
02547   y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
02548   y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
02549 
02550   alen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, adet);
02551 
02552 
02553   temp16len = fast_expansion_sum_zeroelim(8, cxay, 8, axcy, temp16);
02554 
02555   xlen = scale_expansion_zeroelim(temp16len, temp16, bdx, detx);
02556   xxlen = scale_expansion_zeroelim(xlen, detx, bdx, detxx);
02557   xtlen = scale_expansion_zeroelim(temp16len, temp16, bdxtail, detxt);
02558   xxtlen = scale_expansion_zeroelim(xtlen, detxt, bdx, detxxt);
02559   for (i = 0; i < xxtlen; i++) {
02560     detxxt[i] *= 2.0;
02561   }
02562   xtxtlen = scale_expansion_zeroelim(xtlen, detxt, bdxtail, detxtxt);
02563   x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
02564   x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
02565 
02566   ylen = scale_expansion_zeroelim(temp16len, temp16, bdy, dety);
02567   yylen = scale_expansion_zeroelim(ylen, dety, bdy, detyy);
02568   ytlen = scale_expansion_zeroelim(temp16len, temp16, bdytail, detyt);
02569   yytlen = scale_expansion_zeroelim(ytlen, detyt, bdy, detyyt);
02570   for (i = 0; i < yytlen; i++) {
02571     detyyt[i] *= 2.0;
02572   }
02573   ytytlen = scale_expansion_zeroelim(ytlen, detyt, bdytail, detytyt);
02574   y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
02575   y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
02576 
02577   blen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, bdet);
02578 
02579 
02580   temp16len = fast_expansion_sum_zeroelim(8, axby, 8, bxay, temp16);
02581 
02582   xlen = scale_expansion_zeroelim(temp16len, temp16, cdx, detx);
02583   xxlen = scale_expansion_zeroelim(xlen, detx, cdx, detxx);
02584   xtlen = scale_expansion_zeroelim(temp16len, temp16, cdxtail, detxt);
02585   xxtlen = scale_expansion_zeroelim(xtlen, detxt, cdx, detxxt);
02586   for (i = 0; i < xxtlen; i++) {
02587     detxxt[i] *= 2.0;
02588   }
02589   xtxtlen = scale_expansion_zeroelim(xtlen, detxt, cdxtail, detxtxt);
02590   x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
02591   x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
02592 
02593   ylen = scale_expansion_zeroelim(temp16len, temp16, cdy, dety);
02594   yylen = scale_expansion_zeroelim(ylen, dety, cdy, detyy);
02595   ytlen = scale_expansion_zeroelim(temp16len, temp16, cdytail, detyt);
02596   yytlen = scale_expansion_zeroelim(ytlen, detyt, cdy, detyyt);
02597   for (i = 0; i < yytlen; i++) {
02598     detyyt[i] *= 2.0;
02599   }
02600   ytytlen = scale_expansion_zeroelim(ytlen, detyt, cdytail, detytyt);
02601   y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
02602   y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
02603 
02604   clen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, cdet);
02605 
02606   ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
02607   deterlen = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, deter);
02608 
02609   return deter[deterlen - 1];
02610 }
02611 
02612 REAL incircleadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL permanent)
02613 {
02614   INEXACT REAL adx, bdx, cdx, ady, bdy, cdy;
02615   REAL det, errbound;
02616 
02617   INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
02618   REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
02619   REAL bc[4], ca[4], ab[4];
02620   INEXACT REAL bc3, ca3, ab3;
02621   REAL axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32];
02622   int axbclen, axxbclen, aybclen, ayybclen, alen;
02623   REAL bxca[8], bxxca[16], byca[8], byyca[16], bdet[32];
02624   int bxcalen, bxxcalen, bycalen, byycalen, blen;
02625   REAL cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32];
02626   int cxablen, cxxablen, cyablen, cyyablen, clen;
02627   REAL abdet[64];
02628   int ablen;
02629   REAL fin1[1152], fin2[1152];
02630   REAL *finnow, *finother, *finswap;
02631   int finlength;
02632 
02633   REAL adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
02634   INEXACT REAL adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1;
02635   REAL adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0;
02636   REAL aa[4], bb[4], cc[4];
02637   INEXACT REAL aa3, bb3, cc3;
02638   INEXACT REAL ti1, tj1;
02639   REAL ti0, tj0;
02640   REAL u[4], v[4];
02641   INEXACT REAL u3, v3;
02642   REAL temp8[8], temp16a[16], temp16b[16], temp16c[16];
02643   REAL temp32a[32], temp32b[32], temp48[48], temp64[64];
02644   int temp8len, temp16alen, temp16blen, temp16clen;
02645   int temp32alen, temp32blen, temp48len, temp64len;
02646   REAL axtbb[8], axtcc[8], aytbb[8], aytcc[8];
02647   int axtbblen, axtcclen, aytbblen, aytcclen;
02648   REAL bxtaa[8], bxtcc[8], bytaa[8], bytcc[8];
02649   int bxtaalen, bxtcclen, bytaalen, bytcclen;
02650   REAL cxtaa[8], cxtbb[8], cytaa[8], cytbb[8];
02651   int cxtaalen, cxtbblen, cytaalen, cytbblen;
02652   REAL axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8];
02653   int axtbclen, aytbclen, bxtcalen, bytcalen, cxtablen, cytablen;
02654   REAL axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16];
02655   int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen;
02656   REAL axtbctt[8], aytbctt[8], bxtcatt[8];
02657   REAL bytcatt[8], cxtabtt[8], cytabtt[8];
02658   int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen;
02659   REAL abt[8], bct[8], cat[8];
02660   int abtlen, bctlen, catlen;
02661   REAL abtt[4], bctt[4], catt[4];
02662   int abttlen, bcttlen, cattlen;
02663   INEXACT REAL abtt3, bctt3, catt3;
02664   REAL negate;
02665 
02666   INEXACT REAL bvirt;
02667   REAL avirt, bround, around;
02668   INEXACT REAL c;
02669   INEXACT REAL abig;
02670   REAL ahi, alo, bhi, blo;
02671   REAL err1, err2, err3;
02672   INEXACT REAL _i, _j;
02673   REAL _0;
02674 
02675   axtbclen=0;
02676   aytbclen=0;
02677   bxtcalen=0;
02678   bytcalen=0;
02679   cxtablen=0;
02680   cytablen=0;
02681   adx = (REAL) (pa[0] - pd[0]);
02682   bdx = (REAL) (pb[0] - pd[0]);
02683   cdx = (REAL) (pc[0] - pd[0]);
02684   ady = (REAL) (pa[1] - pd[1]);
02685   bdy = (REAL) (pb[1] - pd[1]);
02686   cdy = (REAL) (pc[1] - pd[1]);
02687 
02688   Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
02689   Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
02690   Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
02691   bc[3] = bc3;
02692   axbclen = scale_expansion_zeroelim(4, bc, adx, axbc);
02693   axxbclen = scale_expansion_zeroelim(axbclen, axbc, adx, axxbc);
02694   aybclen = scale_expansion_zeroelim(4, bc, ady, aybc);
02695   ayybclen = scale_expansion_zeroelim(aybclen, aybc, ady, ayybc);
02696   alen = fast_expansion_sum_zeroelim(axxbclen, axxbc, ayybclen, ayybc, adet);
02697 
02698   Two_Product(cdx, ady, cdxady1, cdxady0);
02699   Two_Product(adx, cdy, adxcdy1, adxcdy0);
02700   Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
02701   ca[3] = ca3;
02702   bxcalen = scale_expansion_zeroelim(4, ca, bdx, bxca);
02703   bxxcalen = scale_expansion_zeroelim(bxcalen, bxca, bdx, bxxca);
02704   bycalen = scale_expansion_zeroelim(4, ca, bdy, byca);
02705   byycalen = scale_expansion_zeroelim(bycalen, byca, bdy, byyca);
02706   blen = fast_expansion_sum_zeroelim(bxxcalen, bxxca, byycalen, byyca, bdet);
02707 
02708   Two_Product(adx, bdy, adxbdy1, adxbdy0);
02709   Two_Product(bdx, ady, bdxady1, bdxady0);
02710   Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
02711   ab[3] = ab3;
02712   cxablen = scale_expansion_zeroelim(4, ab, cdx, cxab);
02713   cxxablen = scale_expansion_zeroelim(cxablen, cxab, cdx, cxxab);
02714   cyablen = scale_expansion_zeroelim(4, ab, cdy, cyab);
02715   cyyablen = scale_expansion_zeroelim(cyablen, cyab, cdy, cyyab);
02716   clen = fast_expansion_sum_zeroelim(cxxablen, cxxab, cyyablen, cyyab, cdet);
02717 
02718   ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
02719   finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
02720 
02721   det = estimate(finlength, fin1);
02722   errbound = iccerrboundB * permanent;
02723   if ((det >= errbound) || (-det >= errbound)) {
02724     return det;
02725   }
02726 
02727   Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
02728   Two_Diff_Tail(pa[1], pd[1], ady, adytail);
02729   Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
02730   Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
02731   Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
02732   Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
02733   if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
02734       && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)) {
02735     return det;
02736   }
02737 
02738   errbound = iccerrboundC * permanent + resulterrbound * Absolute(det);
02739   det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail)
02740                                      - (bdy * cdxtail + cdx * bdytail))
02741           + 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx))
02742        + ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail)
02743                                      - (cdy * adxtail + adx * cdytail))
02744           + 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx))
02745        + ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail)
02746                                      - (ady * bdxtail + bdx * adytail))
02747           + 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));
02748   if ((det >= errbound) || (-det >= errbound)) {
02749     return det;
02750   }
02751 
02752   finnow = fin1;
02753   finother = fin2;
02754 
02755   if ((bdxtail != 0.0) || (bdytail != 0.0)
02756       || (cdxtail != 0.0) || (cdytail != 0.0)) {
02757     Square(adx, adxadx1, adxadx0);
02758     Square(ady, adyady1, adyady0);
02759     Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]);
02760     aa[3] = aa3;
02761   }
02762   if ((cdxtail != 0.0) || (cdytail != 0.0)
02763       || (adxtail != 0.0) || (adytail != 0.0)) {
02764     Square(bdx, bdxbdx1, bdxbdx0);
02765     Square(bdy, bdybdy1, bdybdy0);
02766     Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]);
02767     bb[3] = bb3;
02768   }
02769   if ((adxtail != 0.0) || (adytail != 0.0)
02770       || (bdxtail != 0.0) || (bdytail != 0.0)) {
02771     Square(cdx, cdxcdx1, cdxcdx0);
02772     Square(cdy, cdycdy1, cdycdy0);
02773     Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]);
02774     cc[3] = cc3;
02775   }
02776 
02777   if (adxtail != 0.0) {
02778     axtbclen = scale_expansion_zeroelim(4, bc, adxtail, axtbc);
02779     temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, 2.0 * adx,
02780                                           temp16a);
02781 
02782     axtcclen = scale_expansion_zeroelim(4, cc, adxtail, axtcc);
02783     temp16blen = scale_expansion_zeroelim(axtcclen, axtcc, bdy, temp16b);
02784 
02785     axtbblen = scale_expansion_zeroelim(4, bb, adxtail, axtbb);
02786     temp16clen = scale_expansion_zeroelim(axtbblen, axtbb, -cdy, temp16c);
02787 
02788     temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
02789                                             temp16blen, temp16b, temp32a);
02790     temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
02791                                             temp32alen, temp32a, temp48);
02792     finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
02793                                             temp48, finother);
02794     finswap = finnow; finnow = finother; finother = finswap;
02795   }
02796   if (adytail != 0.0) {
02797     aytbclen = scale_expansion_zeroelim(4, bc, adytail, aytbc);
02798     temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, 2.0 * ady,
02799                                           temp16a);
02800 
02801     aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb);
02802     temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b);
02803 
02804     aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc);
02805     temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c);
02806 
02807     temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
02808                                             temp16blen, temp16b, temp32a);
02809     temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
02810                                             temp32alen, temp32a, temp48);
02811     finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
02812                                             temp48, finother);
02813     finswap = finnow; finnow = finother; finother = finswap;
02814   }
02815   if (bdxtail != 0.0) {
02816     bxtcalen = scale_expansion_zeroelim(4, ca, bdxtail, bxtca);
02817     temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, 2.0 * bdx,
02818                                           temp16a);
02819 
02820     bxtaalen = scale_expansion_zeroelim(4, aa, bdxtail, bxtaa);
02821     temp16blen = scale_expansion_zeroelim(bxtaalen, bxtaa, cdy, temp16b);
02822 
02823     bxtcclen = scale_expansion_zeroelim(4, cc, bdxtail, bxtcc);
02824     temp16clen = scale_expansion_zeroelim(bxtcclen, bxtcc, -ady, temp16c);
02825 
02826     temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
02827                                             temp16blen, temp16b, temp32a);
02828     temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
02829                                             temp32alen, temp32a, temp48);
02830     finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
02831                                             temp48, finother);
02832     finswap = finnow; finnow = finother; finother = finswap;
02833   }
02834   if (bdytail != 0.0) {
02835     bytcalen = scale_expansion_zeroelim(4, ca, bdytail, bytca);
02836     temp16alen = scale_expansion_zeroelim(bytcalen, bytca, 2.0 * bdy,
02837                                           temp16a);
02838 
02839     bytcclen = scale_expansion_zeroelim(4, cc, bdytail, bytcc);
02840     temp16blen = scale_expansion_zeroelim(bytcclen, bytcc, adx, temp16b);
02841 
02842     bytaalen = scale_expansion_zeroelim(4, aa, bdytail, bytaa);
02843     temp16clen = scale_expansion_zeroelim(bytaalen, bytaa, -cdx, temp16c);
02844 
02845     temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
02846                                             temp16blen, temp16b, temp32a);
02847     temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
02848                                             temp32alen, temp32a, temp48);
02849     finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
02850                                             temp48, finother);
02851     finswap = finnow; finnow = finother; finother = finswap;
02852   }
02853   if (cdxtail != 0.0) {
02854     cxtablen = scale_expansion_zeroelim(4, ab, cdxtail, cxtab);
02855     temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, 2.0 * cdx,
02856                                           temp16a);
02857 
02858     cxtbblen = scale_expansion_zeroelim(4, bb, cdxtail, cxtbb);
02859     temp16blen = scale_expansion_zeroelim(cxtbblen, cxtbb, ady, temp16b);
02860 
02861     cxtaalen = scale_expansion_zeroelim(4, aa, cdxtail, cxtaa);
02862     temp16clen = scale_expansion_zeroelim(cxtaalen, cxtaa, -bdy, temp16c);
02863 
02864     temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
02865                                             temp16blen, temp16b, temp32a);
02866     temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
02867                                             temp32alen, temp32a, temp48);
02868     finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
02869                                             temp48, finother);
02870     finswap = finnow; finnow = finother; finother = finswap;
02871   }
02872   if (cdytail != 0.0) {
02873     cytablen = scale_expansion_zeroelim(4, ab, cdytail, cytab);
02874     temp16alen = scale_expansion_zeroelim(cytablen, cytab, 2.0 * cdy,
02875                                           temp16a);
02876 
02877     cytaalen = scale_expansion_zeroelim(4, aa, cdytail, cytaa);
02878     temp16blen = scale_expansion_zeroelim(cytaalen, cytaa, bdx, temp16b);
02879 
02880     cytbblen = scale_expansion_zeroelim(4, bb, cdytail, cytbb);
02881     temp16clen = scale_expansion_zeroelim(cytbblen, cytbb, -adx, temp16c);
02882 
02883     temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
02884                                             temp16blen, temp16b, temp32a);
02885     temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
02886                                             temp32alen, temp32a, temp48);
02887     finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
02888                                             temp48, finother);
02889     finswap = finnow; finnow = finother; finother = finswap;
02890   }
02891 
02892   if ((adxtail != 0.0) || (adytail != 0.0)) {
02893     if ((bdxtail != 0.0) || (bdytail != 0.0)
02894         || (cdxtail != 0.0) || (cdytail != 0.0)) {
02895       Two_Product(bdxtail, cdy, ti1, ti0);
02896       Two_Product(bdx, cdytail, tj1, tj0);
02897       Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
02898       u[3] = u3;
02899       negate = -bdy;
02900       Two_Product(cdxtail, negate, ti1, ti0);
02901       negate = -bdytail;
02902       Two_Product(cdx, negate, tj1, tj0);
02903       Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
02904       v[3] = v3;
02905       bctlen = fast_expansion_sum_zeroelim(4, u, 4, v, bct);
02906 
02907       Two_Product(bdxtail, cdytail, ti1, ti0);
02908       Two_Product(cdxtail, bdytail, tj1, tj0);
02909       Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]);
02910       bctt[3] = bctt3;
02911       bcttlen = 4;
02912     } else {
02913       bct[0] = 0.0;
02914       bctlen = 1;
02915       bctt[0] = 0.0;
02916       bcttlen = 1;
02917     }
02918 
02919     if (adxtail != 0.0) {
02920       temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, adxtail, temp16a);
02921       axtbctlen = scale_expansion_zeroelim(bctlen, bct, adxtail, axtbct);
02922       temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, 2.0 * adx,
02923                                             temp32a);
02924       temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
02925                                               temp32alen, temp32a, temp48);
02926       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
02927                                               temp48, finother);
02928       finswap = finnow; finnow = finother; finother = finswap;
02929       if (bdytail != 0.0) {
02930         temp8len = scale_expansion_zeroelim(4, cc, adxtail, temp8);
02931         temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
02932                                               temp16a);
02933         finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
02934                                                 temp16a, finother);
02935         finswap = finnow; finnow = finother; finother = finswap;
02936       }
02937       if (cdytail != 0.0) {
02938         temp8len = scale_expansion_zeroelim(4, bb, -adxtail, temp8);
02939         temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
02940                                               temp16a);
02941         finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
02942                                                 temp16a, finother);
02943         finswap = finnow; finnow = finother; finother = finswap;
02944       }
02945 
02946       temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, adxtail,
02947                                             temp32a);
02948       axtbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adxtail, axtbctt);
02949       temp16alen = scale_expansion_zeroelim(axtbcttlen, axtbctt, 2.0 * adx,
02950                                             temp16a);
02951       temp16blen = scale_expansion_zeroelim(axtbcttlen, axtbctt, adxtail,
02952                                             temp16b);
02953       temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
02954                                               temp16blen, temp16b, temp32b);
02955       temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
02956                                               temp32blen, temp32b, temp64);
02957       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
02958                                               temp64, finother);
02959       finswap = finnow; finnow = finother; finother = finswap;
02960     }
02961     if (adytail != 0.0) {
02962       temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, adytail, temp16a);
02963       aytbctlen = scale_expansion_zeroelim(bctlen, bct, adytail, aytbct);
02964       temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, 2.0 * ady,
02965                                             temp32a);
02966       temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
02967                                               temp32alen, temp32a, temp48);
02968       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
02969                                               temp48, finother);
02970       finswap = finnow; finnow = finother; finother = finswap;
02971 
02972 
02973       temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, adytail,
02974                                             temp32a);
02975       aytbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adytail, aytbctt);
02976       temp16alen = scale_expansion_zeroelim(aytbcttlen, aytbctt, 2.0 * ady,
02977                                             temp16a);
02978       temp16blen = scale_expansion_zeroelim(aytbcttlen, aytbctt, adytail,
02979                                             temp16b);
02980       temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
02981                                               temp16blen, temp16b, temp32b);
02982       temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
02983                                               temp32blen, temp32b, temp64);
02984       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
02985                                               temp64, finother);
02986       finswap = finnow; finnow = finother; finother = finswap;
02987     }
02988   }
02989   if ((bdxtail != 0.0) || (bdytail != 0.0)) {
02990     if ((cdxtail != 0.0) || (cdytail != 0.0)
02991         || (adxtail != 0.0) || (adytail != 0.0)) {
02992       Two_Product(cdxtail, ady, ti1, ti0);
02993       Two_Product(cdx, adytail, tj1, tj0);
02994       Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
02995       u[3] = u3;
02996       negate = -cdy;
02997       Two_Product(adxtail, negate, ti1, ti0);
02998       negate = -cdytail;
02999       Two_Product(adx, negate, tj1, tj0);
03000       Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
03001       v[3] = v3;
03002       catlen = fast_expansion_sum_zeroelim(4, u, 4, v, cat);
03003 
03004       Two_Product(cdxtail, adytail, ti1, ti0);
03005       Two_Product(adxtail, cdytail, tj1, tj0);
03006       Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]);
03007       catt[3] = catt3;
03008       cattlen = 4;
03009     } else {
03010       cat[0] = 0.0;
03011       catlen = 1;
03012       catt[0] = 0.0;
03013       cattlen = 1;
03014     }
03015 
03016     if (bdxtail != 0.0) {
03017       temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, bdxtail, temp16a);
03018       bxtcatlen = scale_expansion_zeroelim(catlen, cat, bdxtail, bxtcat);
03019       temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, 2.0 * bdx,
03020                                             temp32a);
03021       temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03022                                               temp32alen, temp32a, temp48);
03023       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
03024                                               temp48, finother);
03025       finswap = finnow; finnow = finother; finother = finswap;
03026       if (cdytail != 0.0) {
03027         temp8len = scale_expansion_zeroelim(4, aa, bdxtail, temp8);
03028         temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
03029                                               temp16a);
03030         finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
03031                                                 temp16a, finother);
03032         finswap = finnow; finnow = finother; finother = finswap;
03033       }
03034       if (adytail != 0.0) {
03035         temp8len = scale_expansion_zeroelim(4, cc, -bdxtail, temp8);
03036         temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
03037                                               temp16a);
03038         finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
03039                                                 temp16a, finother);
03040         finswap = finnow; finnow = finother; finother = finswap;
03041       }
03042 
03043       temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, bdxtail,
03044                                             temp32a);
03045       bxtcattlen = scale_expansion_zeroelim(cattlen, catt, bdxtail, bxtcatt);
03046       temp16alen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, 2.0 * bdx,
03047                                             temp16a);
03048       temp16blen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, bdxtail,
03049                                             temp16b);
03050       temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03051                                               temp16blen, temp16b, temp32b);
03052       temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03053                                               temp32blen, temp32b, temp64);
03054       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
03055                                               temp64, finother);
03056       finswap = finnow; finnow = finother; finother = finswap;
03057     }
03058     if (bdytail != 0.0) {
03059       temp16alen = scale_expansion_zeroelim(bytcalen, bytca, bdytail, temp16a);
03060       bytcatlen = scale_expansion_zeroelim(catlen, cat, bdytail, bytcat);
03061       temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, 2.0 * bdy,
03062                                             temp32a);
03063       temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03064                                               temp32alen, temp32a, temp48);
03065       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
03066                                               temp48, finother);
03067       finswap = finnow; finnow = finother; finother = finswap;
03068 
03069 
03070       temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, bdytail,
03071                                             temp32a);
03072       bytcattlen = scale_expansion_zeroelim(cattlen, catt, bdytail, bytcatt);
03073       temp16alen = scale_expansion_zeroelim(bytcattlen, bytcatt, 2.0 * bdy,
03074                                             temp16a);
03075       temp16blen = scale_expansion_zeroelim(bytcattlen, bytcatt, bdytail,
03076                                             temp16b);
03077       temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03078                                               temp16blen, temp16b, temp32b);
03079       temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03080                                               temp32blen, temp32b, temp64);
03081       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
03082                                               temp64, finother);
03083       finswap = finnow; finnow = finother; finother = finswap;
03084     }
03085   }
03086   if ((cdxtail != 0.0) || (cdytail != 0.0)) {
03087     if ((adxtail != 0.0) || (adytail != 0.0)
03088         || (bdxtail != 0.0) || (bdytail != 0.0)) {
03089       Two_Product(adxtail, bdy, ti1, ti0);
03090       Two_Product(adx, bdytail, tj1, tj0);
03091       Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
03092       u[3] = u3;
03093       negate = -ady;
03094       Two_Product(bdxtail, negate, ti1, ti0);
03095       negate = -adytail;
03096       Two_Product(bdx, negate, tj1, tj0);
03097       Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
03098       v[3] = v3;
03099       abtlen = fast_expansion_sum_zeroelim(4, u, 4, v, abt);
03100 
03101       Two_Product(adxtail, bdytail, ti1, ti0);
03102       Two_Product(bdxtail, adytail, tj1, tj0);
03103       Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]);
03104       abtt[3] = abtt3;
03105       abttlen = 4;
03106     } else {
03107       abt[0] = 0.0;
03108       abtlen = 1;
03109       abtt[0] = 0.0;
03110       abttlen = 1;
03111     }
03112 
03113     if (cdxtail != 0.0) {
03114       temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, cdxtail, temp16a);
03115       cxtabtlen = scale_expansion_zeroelim(abtlen, abt, cdxtail, cxtabt);
03116       temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, 2.0 * cdx,
03117                                             temp32a);
03118       temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03119                                               temp32alen, temp32a, temp48);
03120       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
03121                                               temp48, finother);
03122       finswap = finnow; finnow = finother; finother = finswap;
03123       if (adytail != 0.0) {
03124         temp8len = scale_expansion_zeroelim(4, bb, cdxtail, temp8);
03125         temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
03126                                               temp16a);
03127         finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
03128                                                 temp16a, finother);
03129         finswap = finnow; finnow = finother; finother = finswap;
03130       }
03131       if (bdytail != 0.0) {
03132         temp8len = scale_expansion_zeroelim(4, aa, -cdxtail, temp8);
03133         temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
03134                                               temp16a);
03135         finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
03136                                                 temp16a, finother);
03137         finswap = finnow; finnow = finother; finother = finswap;
03138       }
03139 
03140       temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, cdxtail,
03141                                             temp32a);
03142       cxtabttlen = scale_expansion_zeroelim(abttlen, abtt, cdxtail, cxtabtt);
03143       temp16alen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, 2.0 * cdx,
03144                                             temp16a);
03145       temp16blen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, cdxtail,
03146                                             temp16b);
03147       temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03148                                               temp16blen, temp16b, temp32b);
03149       temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03150                                               temp32blen, temp32b, temp64);
03151       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
03152                                               temp64, finother);
03153       finswap = finnow; finnow = finother; finother = finswap;
03154     }
03155     if (cdytail != 0.0) {
03156       temp16alen = scale_expansion_zeroelim(cytablen, cytab, cdytail, temp16a);
03157       cytabtlen = scale_expansion_zeroelim(abtlen, abt, cdytail, cytabt);
03158       temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, 2.0 * cdy,
03159                                             temp32a);
03160       temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03161                                               temp32alen, temp32a, temp48);
03162       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
03163                                               temp48, finother);
03164       finswap = finnow; finnow = finother; finother = finswap;
03165 
03166 
03167       temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, cdytail,
03168                                             temp32a);
03169       cytabttlen = scale_expansion_zeroelim(abttlen, abtt, cdytail, cytabtt);
03170       temp16alen = scale_expansion_zeroelim(cytabttlen, cytabtt, 2.0 * cdy,
03171                                             temp16a);
03172       temp16blen = scale_expansion_zeroelim(cytabttlen, cytabtt, cdytail,
03173                                             temp16b);
03174       temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
03175                                               temp16blen, temp16b, temp32b);
03176       temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03177                                               temp32blen, temp32b, temp64);
03178       finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
03179                                               temp64, finother);
03180       finswap = finnow; finnow = finother; finother = finswap;
03181     }
03182   }
03183 
03184   return finnow[finlength - 1];
03185 }
03186 
03187 REAL incircle(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
03188 {
03189   REAL adx, bdx, cdx, ady, bdy, cdy;
03190   REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
03191   REAL alift, blift, clift;
03192   REAL det;
03193   REAL permanent, errbound;
03194 
03195   adx = pa[0] - pd[0];
03196   bdx = pb[0] - pd[0];
03197   cdx = pc[0] - pd[0];
03198   ady = pa[1] - pd[1];
03199   bdy = pb[1] - pd[1];
03200   cdy = pc[1] - pd[1];
03201 
03202   bdxcdy = bdx * cdy;
03203   cdxbdy = cdx * bdy;
03204   alift = adx * adx + ady * ady;
03205 
03206   cdxady = cdx * ady;
03207   adxcdy = adx * cdy;
03208   blift = bdx * bdx + bdy * bdy;
03209 
03210   adxbdy = adx * bdy;
03211   bdxady = bdx * ady;
03212   clift = cdx * cdx + cdy * cdy;
03213 
03214   det = alift * (bdxcdy - cdxbdy)
03215       + blift * (cdxady - adxcdy)
03216       + clift * (adxbdy - bdxady);
03217 
03218   permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift
03219             + (Absolute(cdxady) + Absolute(adxcdy)) * blift
03220             + (Absolute(adxbdy) + Absolute(bdxady)) * clift;
03221   errbound = iccerrboundA * permanent;
03222   if ((det > errbound) || (-det > errbound)) {
03223     return det;
03224   }
03225 
03226   return incircleadapt(pa, pb, pc, pd, permanent);
03227 }
03228 
03229 /*****************************************************************************/
03230 /*                                                                           */
03231 /*  inspherefast()   Approximate 3D insphere test.  Nonrobust.               */
03232 /*  insphereexact()   Exact 3D insphere test.  Robust.                       */
03233 /*  insphereslow()   Another exact 3D insphere test.  Robust.                */
03234 /*  insphere()   Adaptive exact 3D insphere test.  Robust.                   */
03235 /*                                                                           */
03236 /*               Return a positive value if the point pe lies inside the     */
03237 /*               sphere passing through pa, pb, pc, and pd; a negative value */
03238 /*               if it lies outside; and zero if the five points are         */
03239 /*               cospherical.  The points pa, pb, pc, and pd must be ordered */
03240 /*               so that they have a positive orientation (as defined by     */
03241 /*               orient3d()), or the sign of the result will be reversed.    */
03242 /*                                                                           */
03243 /*  Only the first and last routine should be used; the middle two are for   */
03244 /*  timings.                                                                 */
03245 /*                                                                           */
03246 /*  The last three use exact arithmetic to ensure a correct answer.  The     */
03247 /*  result returned is the determinant of a matrix.  In insphere() only,     */
03248 /*  this determinant is computed adaptively, in the sense that exact         */
03249 /*  arithmetic is used only to the degree it is needed to ensure that the    */
03250 /*  returned value has the correct sign.  Hence, insphere() is usually quite */
03251 /*  fast, but will run more slowly when the input points are cospherical or  */
03252 /*  nearly so.                                                               */
03253 /*                                                                           */
03254 /*****************************************************************************/
03255 
03256 REAL inspherefast(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
03257 {
03258   REAL aex, bex, cex, dex;
03259   REAL aey, bey, cey, dey;
03260   REAL aez, bez, cez, dez;
03261   REAL alift, blift, clift, dlift;
03262   REAL ab, bc, cd, da, ac, bd;
03263   REAL abc, bcd, cda, dab;
03264 
03265   aex = pa[0] - pe[0];
03266   bex = pb[0] - pe[0];
03267   cex = pc[0] - pe[0];
03268   dex = pd[0] - pe[0];
03269   aey = pa[1] - pe[1];
03270   bey = pb[1] - pe[1];
03271   cey = pc[1] - pe[1];
03272   dey = pd[1] - pe[1];
03273   aez = pa[2] - pe[2];
03274   bez = pb[2] - pe[2];
03275   cez = pc[2] - pe[2];
03276   dez = pd[2] - pe[2];
03277 
03278   ab = aex * bey - bex * aey;
03279   bc = bex * cey - cex * bey;
03280   cd = cex * dey - dex * cey;
03281   da = dex * aey - aex * dey;
03282 
03283   ac = aex * cey - cex * aey;
03284   bd = bex * dey - dex * bey;
03285 
03286   abc = aez * bc - bez * ac + cez * ab;
03287   bcd = bez * cd - cez * bd + dez * bc;
03288   cda = cez * da + dez * ac + aez * cd;
03289   dab = dez * ab + aez * bd + bez * da;
03290 
03291   alift = aex * aex + aey * aey + aez * aez;
03292   blift = bex * bex + bey * bey + bez * bez;
03293   clift = cex * cex + cey * cey + cez * cez;
03294   dlift = dex * dex + dey * dey + dez * dez;
03295 
03296   return (dlift * abc - clift * dab) + (blift * cda - alift * bcd);
03297 }
03298 
03299 REAL insphereexact(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
03300 {
03301   INEXACT REAL axby1, bxcy1, cxdy1, dxey1, exay1;
03302   INEXACT REAL bxay1, cxby1, dxcy1, exdy1, axey1;
03303   INEXACT REAL axcy1, bxdy1, cxey1, dxay1, exby1;
03304   INEXACT REAL cxay1, dxby1, excy1, axdy1, bxey1;
03305   REAL axby0, bxcy0, cxdy0, dxey0, exay0;
03306   REAL bxay0, cxby0, dxcy0, exdy0, axey0;
03307   REAL axcy0, bxdy0, cxey0, dxay0, exby0;
03308   REAL cxay0, dxby0, excy0, axdy0, bxey0;
03309   REAL ab[4], bc[4], cd[4], de[4], ea[4];
03310   REAL ac[4], bd[4], ce[4], da[4], eb[4];
03311   REAL temp8a[8], temp8b[8], temp16[16];
03312   int temp8alen, temp8blen, temp16len;
03313   REAL abc[24], bcd[24], cde[24], dea[24], eab[24];
03314   REAL abd[24], bce[24], cda[24], deb[24], eac[24];
03315   int abclen, bcdlen, cdelen, dealen, eablen;
03316   int abdlen, bcelen, cdalen, deblen, eaclen;
03317   REAL temp48a[48], temp48b[48];
03318   int temp48alen, temp48blen;
03319   REAL abcd[96], bcde[96], cdea[96], deab[96], eabc[96];
03320   int abcdlen, bcdelen, cdealen, deablen, eabclen;
03321   REAL temp192[192];
03322   REAL det384x[384], det384y[384], det384z[384];
03323   int xlen, ylen, zlen;
03324   REAL detxy[768];
03325   int xylen;
03326   REAL adet[1152], bdet[1152], cdet[1152], ddet[1152], edet[1152];
03327   int alen, blen, clen, dlen, elen;
03328   REAL abdet[2304], cddet[2304], cdedet[3456];
03329   int ablen, cdlen;
03330   REAL deter[5760];
03331   int deterlen;
03332   int i;
03333 
03334   INEXACT REAL bvirt;
03335   REAL avirt, bround, around;
03336   INEXACT REAL c;
03337   INEXACT REAL abig;
03338   REAL ahi, alo, bhi, blo;
03339   REAL err1, err2, err3;
03340   INEXACT REAL _i, _j;
03341   REAL _0;
03342 
03343   Two_Product(pa[0], pb[1], axby1, axby0);
03344   Two_Product(pb[0], pa[1], bxay1, bxay0);
03345   Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]);
03346 
03347   Two_Product(pb[0], pc[1], bxcy1, bxcy0);
03348   Two_Product(pc[0], pb[1], cxby1, cxby0);
03349   Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]);
03350 
03351   Two_Product(pc[0], pd[1], cxdy1, cxdy0);
03352   Two_Product(pd[0], pc[1], dxcy1, dxcy0);
03353   Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]);
03354 
03355   Two_Product(pd[0], pe[1], dxey1, dxey0);
03356   Two_Product(pe[0], pd[1], exdy1, exdy0);
03357   Two_Two_Diff(dxey1, dxey0, exdy1, exdy0, de[3], de[2], de[1], de[0]);
03358 
03359   Two_Product(pe[0], pa[1], exay1, exay0);
03360   Two_Product(pa[0], pe[1], axey1, axey0);
03361   Two_Two_Diff(exay1, exay0, axey1, axey0, ea[3], ea[2], ea[1], ea[0]);
03362 
03363   Two_Product(pa[0], pc[1], axcy1, axcy0);
03364   Two_Product(pc[0], pa[1], cxay1, cxay0);
03365   Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]);
03366 
03367   Two_Product(pb[0], pd[1], bxdy1, bxdy0);
03368   Two_Product(pd[0], pb[1], dxby1, dxby0);
03369   Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]);
03370 
03371   Two_Product(pc[0], pe[1], cxey1, cxey0);
03372   Two_Product(pe[0], pc[1], excy1, excy0);
03373   Two_Two_Diff(cxey1, cxey0, excy1, excy0, ce[3], ce[2], ce[1], ce[0]);
03374 
03375   Two_Product(pd[0], pa[1], dxay1, dxay0);
03376   Two_Product(pa[0], pd[1], axdy1, axdy0);
03377   Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]);
03378 
03379   Two_Product(pe[0], pb[1], exby1, exby0);
03380   Two_Product(pb[0], pe[1], bxey1, bxey0);
03381   Two_Two_Diff(exby1, exby0, bxey1, bxey0, eb[3], eb[2], eb[1], eb[0]);
03382 
03383   temp8alen = scale_expansion_zeroelim(4, bc, pa[2], temp8a);
03384   temp8blen = scale_expansion_zeroelim(4, ac, -pb[2], temp8b);
03385   temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
03386                                           temp16);
03387   temp8alen = scale_expansion_zeroelim(4, ab, pc[2], temp8a);
03388   abclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
03389                                        abc);
03390 
03391   temp8alen = scale_expansion_zeroelim(4, cd, pb[2], temp8a);
03392   temp8blen = scale_expansion_zeroelim(4, bd, -pc[2], temp8b);
03393   temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
03394                                           temp16);
03395   temp8alen = scale_expansion_zeroelim(4, bc, pd[2], temp8a);
03396   bcdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
03397                                        bcd);
03398 
03399   temp8alen = scale_expansion_zeroelim(4, de, pc[2], temp8a);
03400   temp8blen = scale_expansion_zeroelim(4, ce, -pd[2], temp8b);
03401   temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
03402                                           temp16);
03403   temp8alen = scale_expansion_zeroelim(4, cd, pe[2], temp8a);
03404   cdelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
03405                                        cde);
03406 
03407   temp8alen = scale_expansion_zeroelim(4, ea, pd[2], temp8a);
03408   temp8blen = scale_expansion_zeroelim(4, da, -pe[2], temp8b);
03409   temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
03410                                           temp16);
03411   temp8alen = scale_expansion_zeroelim(4, de, pa[2], temp8a);
03412   dealen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
03413                                        dea);
03414 
03415   temp8alen = scale_expansion_zeroelim(4, ab, pe[2], temp8a);
03416   temp8blen = scale_expansion_zeroelim(4, eb, -pa[2], temp8b);
03417   temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
03418                                           temp16);
03419   temp8alen = scale_expansion_zeroelim(4, ea, pb[2], temp8a);
03420   eablen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
03421                                        eab);
03422 
03423   temp8alen = scale_expansion_zeroelim(4, bd, pa[2], temp8a);
03424   temp8blen = scale_expansion_zeroelim(4, da, pb[2], temp8b);
03425   temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
03426                                           temp16);
03427   temp8alen = scale_expansion_zeroelim(4, ab, pd[2], temp8a);
03428   abdlen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
03429                                        abd);
03430 
03431   temp8alen = scale_expansion_zeroelim(4, ce, pb[2], temp8a);
03432   temp8blen = scale_expansion_zeroelim(4, eb, pc[2], temp8b);
03433   temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
03434                                           temp16);
03435   temp8alen = scale_expansion_zeroelim(4, bc, pe[2], temp8a);
03436   bcelen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
03437                                        bce);
03438 
03439   temp8alen = scale_expansion_zeroelim(4, da, pc[2], temp8a);
03440   temp8blen = scale_expansion_zeroelim(4, ac, pd[2], temp8b);
03441   temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
03442                                           temp16);
03443   temp8alen = scale_expansion_zeroelim(4, cd, pa[2], temp8a);
03444   cdalen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
03445                                        cda);
03446 
03447   temp8alen = scale_expansion_zeroelim(4, eb, pd[2], temp8a);
03448   temp8blen = scale_expansion_zeroelim(4, bd, pe[2], temp8b);
03449   temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
03450                                           temp16);
03451   temp8alen = scale_expansion_zeroelim(4, de, pb[2], temp8a);
03452   deblen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
03453                                        deb);
03454 
03455   temp8alen = scale_expansion_zeroelim(4, ac, pe[2], temp8a);
03456   temp8blen = scale_expansion_zeroelim(4, ce, pa[2], temp8b);
03457   temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp8blen, temp8b,
03458                                           temp16);
03459   temp8alen = scale_expansion_zeroelim(4, ea, pc[2], temp8a);
03460   eaclen = fast_expansion_sum_zeroelim(temp8alen, temp8a, temp16len, temp16,
03461                                        eac);
03462 
03463   temp48alen = fast_expansion_sum_zeroelim(cdelen, cde, bcelen, bce, temp48a);
03464   temp48blen = fast_expansion_sum_zeroelim(deblen, deb, bcdlen, bcd, temp48b);
03465   for (i = 0; i < temp48blen; i++) {
03466     temp48b[i] = -temp48b[i];
03467   }
03468   bcdelen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
03469                                         temp48blen, temp48b, bcde);
03470   xlen = scale_expansion_zeroelim(bcdelen, bcde, pa[0], temp192);
03471   xlen = scale_expansion_zeroelim(xlen, temp192, pa[0], det384x);
03472   ylen = scale_expansion_zeroelim(bcdelen, bcde, pa[1], temp192);
03473   ylen = scale_expansion_zeroelim(ylen, temp192, pa[1], det384y);
03474   zlen = scale_expansion_zeroelim(bcdelen, bcde, pa[2], temp192);
03475   zlen = scale_expansion_zeroelim(zlen, temp192, pa[2], det384z);
03476   xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
03477   alen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, adet);
03478 
03479   temp48alen = fast_expansion_sum_zeroelim(dealen, dea, cdalen, cda, temp48a);
03480   temp48blen = fast_expansion_sum_zeroelim(eaclen, eac, cdelen, cde, temp48b);
03481   for (i = 0; i < temp48blen; i++) {
03482     temp48b[i] = -temp48b[i];
03483   }
03484   cdealen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
03485                                         temp48blen, temp48b, cdea);
03486   xlen = scale_expansion_zeroelim(cdealen, cdea, pb[0], temp192);
03487   xlen = scale_expansion_zeroelim(xlen, temp192, pb[0], det384x);
03488   ylen = scale_expansion_zeroelim(cdealen, cdea, pb[1], temp192);
03489   ylen = scale_expansion_zeroelim(ylen, temp192, pb[1], det384y);
03490   zlen = scale_expansion_zeroelim(cdealen, cdea, pb[2], temp192);
03491   zlen = scale_expansion_zeroelim(zlen, temp192, pb[2], det384z);
03492   xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
03493   blen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, bdet);
03494 
03495   temp48alen = fast_expansion_sum_zeroelim(eablen, eab, deblen, deb, temp48a);
03496   temp48blen = fast_expansion_sum_zeroelim(abdlen, abd, dealen, dea, temp48b);
03497   for (i = 0; i < temp48blen; i++) {
03498     temp48b[i] = -temp48b[i];
03499   }
03500   deablen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
03501                                         temp48blen, temp48b, deab);
03502   xlen = scale_expansion_zeroelim(deablen, deab, pc[0], temp192);
03503   xlen = scale_expansion_zeroelim(xlen, temp192, pc[0], det384x);
03504   ylen = scale_expansion_zeroelim(deablen, deab, pc[1], temp192);
03505   ylen = scale_expansion_zeroelim(ylen, temp192, pc[1], det384y);
03506   zlen = scale_expansion_zeroelim(deablen, deab, pc[2], temp192);
03507   zlen = scale_expansion_zeroelim(zlen, temp192, pc[2], det384z);
03508   xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
03509   clen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, cdet);
03510 
03511   temp48alen = fast_expansion_sum_zeroelim(abclen, abc, eaclen, eac, temp48a);
03512   temp48blen = fast_expansion_sum_zeroelim(bcelen, bce, eablen, eab, temp48b);
03513   for (i = 0; i < temp48blen; i++) {
03514     temp48b[i] = -temp48b[i];
03515   }
03516   eabclen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
03517                                         temp48blen, temp48b, eabc);
03518   xlen = scale_expansion_zeroelim(eabclen, eabc, pd[0], temp192);
03519   xlen = scale_expansion_zeroelim(xlen, temp192, pd[0], det384x);
03520   ylen = scale_expansion_zeroelim(eabclen, eabc, pd[1], temp192);
03521   ylen = scale_expansion_zeroelim(ylen, temp192, pd[1], det384y);
03522   zlen = scale_expansion_zeroelim(eabclen, eabc, pd[2], temp192);
03523   zlen = scale_expansion_zeroelim(zlen, temp192, pd[2], det384z);
03524   xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
03525   dlen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, ddet);
03526 
03527   temp48alen = fast_expansion_sum_zeroelim(bcdlen, bcd, abdlen, abd, temp48a);
03528   temp48blen = fast_expansion_sum_zeroelim(cdalen, cda, abclen, abc, temp48b);
03529   for (i = 0; i < temp48blen; i++) {
03530     temp48b[i] = -temp48b[i];
03531   }
03532   abcdlen = fast_expansion_sum_zeroelim(temp48alen, temp48a,
03533                                         temp48blen, temp48b, abcd);
03534   xlen = scale_expansion_zeroelim(abcdlen, abcd, pe[0], temp192);
03535   xlen = scale_expansion_zeroelim(xlen, temp192, pe[0], det384x);
03536   ylen = scale_expansion_zeroelim(abcdlen, abcd, pe[1], temp192);
03537   ylen = scale_expansion_zeroelim(ylen, temp192, pe[1], det384y);
03538   zlen = scale_expansion_zeroelim(abcdlen, abcd, pe[2], temp192);
03539   zlen = scale_expansion_zeroelim(zlen, temp192, pe[2], det384z);
03540   xylen = fast_expansion_sum_zeroelim(xlen, det384x, ylen, det384y, detxy);
03541   elen = fast_expansion_sum_zeroelim(xylen, detxy, zlen, det384z, edet);
03542 
03543   ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
03544   cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
03545   cdelen = fast_expansion_sum_zeroelim(cdlen, cddet, elen, edet, cdedet);
03546   deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdelen, cdedet, deter);
03547 
03548   return deter[deterlen - 1];
03549 }
03550 
03551 REAL insphereslow(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
03552 {
03553   INEXACT REAL aex, bex, cex, dex, aey, bey, cey, dey, aez, bez, cez, dez;
03554   REAL aextail, bextail, cextail, dextail;
03555   REAL aeytail, beytail, ceytail, deytail;
03556   REAL aeztail, beztail, ceztail, deztail;
03557   REAL negate, negatetail;
03558   INEXACT REAL axby7, bxcy7, cxdy7, dxay7, axcy7, bxdy7;
03559   INEXACT REAL bxay7, cxby7, dxcy7, axdy7, cxay7, dxby7;
03560   REAL axby[8], bxcy[8], cxdy[8], dxay[8], axcy[8], bxdy[8];
03561   REAL bxay[8], cxby[8], dxcy[8], axdy[8], cxay[8], dxby[8];
03562   REAL ab[16], bc[16], cd[16], da[16], ac[16], bd[16];
03563   int ablen, bclen, cdlen, dalen, aclen, bdlen;
03564   REAL temp32a[32], temp32b[32], temp64a[64], temp64b[64], temp64c[64];
03565   int temp32alen, temp32blen, temp64alen, temp64blen, temp64clen;
03566   REAL temp128[128], temp192[192];
03567   int temp128len, temp192len;
03568   REAL detx[384], detxx[768], detxt[384], detxxt[768], detxtxt[768];
03569   int xlen, xxlen, xtlen, xxtlen, xtxtlen;
03570   REAL x1[1536], x2[2304];
03571   int x1len, x2len;
03572   REAL dety[384], detyy[768], detyt[384], detyyt[768], detytyt[768];
03573   int ylen, yylen, ytlen, yytlen, ytytlen;
03574   REAL y1[1536], y2[2304];
03575   int y1len, y2len;
03576   REAL detz[384], detzz[768], detzt[384], detzzt[768], detztzt[768];
03577   int zlen, zzlen, ztlen, zztlen, ztztlen;
03578   REAL z1[1536], z2[2304];
03579   int z1len, z2len;
03580   REAL detxy[4608];
03581   int xylen;
03582   REAL adet[6912], bdet[6912], cdet[6912], ddet[6912];
03583   int alen, blen, clen, dlen;
03584   REAL abdet[13824], cddet[13824], deter[27648];
03585   int deterlen;
03586   int i;
03587 
03588   INEXACT REAL bvirt;
03589   REAL avirt, bround, around;
03590   INEXACT REAL c;
03591   INEXACT REAL abig;
03592   REAL a0hi, a0lo, a1hi, a1lo, bhi, blo;
03593   REAL err1, err2, err3;
03594   INEXACT REAL _i, _j, _k, _l, _m, _n;
03595   REAL _0, _1, _2;
03596 
03597   Two_Diff(pa[0], pe[0], aex, aextail);
03598   Two_Diff(pa[1], pe[1], aey, aeytail);
03599   Two_Diff(pa[2], pe[2], aez, aeztail);
03600   Two_Diff(pb[0], pe[0], bex, bextail);
03601   Two_Diff(pb[1], pe[1], bey, beytail);
03602   Two_Diff(pb[2], pe[2], bez, beztail);
03603   Two_Diff(pc[0], pe[0], cex, cextail);
03604   Two_Diff(pc[1], pe[1], cey, ceytail);
03605   Two_Diff(pc[2], pe[2], cez, ceztail);
03606   Two_Diff(pd[0], pe[0], dex, dextail);
03607   Two_Diff(pd[1], pe[1], dey, deytail);
03608   Two_Diff(pd[2], pe[2], dez, deztail);
03609 
03610   Two_Two_Product(aex, aextail, bey, beytail,
03611                   axby7, axby[6], axby[5], axby[4],
03612                   axby[3], axby[2], axby[1], axby[0]);
03613   axby[7] = axby7;
03614   negate = -aey;
03615   negatetail = -aeytail;
03616   Two_Two_Product(bex, bextail, negate, negatetail,
03617                   bxay7, bxay[6], bxay[5], bxay[4],
03618                   bxay[3], bxay[2], bxay[1], bxay[0]);
03619   bxay[7] = bxay7;
03620   ablen = fast_expansion_sum_zeroelim(8, axby, 8, bxay, ab);
03621   Two_Two_Product(bex, bextail, cey, ceytail,
03622                   bxcy7, bxcy[6], bxcy[5], bxcy[4],
03623                   bxcy[3], bxcy[2], bxcy[1], bxcy[0]);
03624   bxcy[7] = bxcy7;
03625   negate = -bey;
03626   negatetail = -beytail;
03627   Two_Two_Product(cex, cextail, negate, negatetail,
03628                   cxby7, cxby[6], cxby[5], cxby[4],
03629                   cxby[3], cxby[2], cxby[1], cxby[0]);
03630   cxby[7] = cxby7;
03631   bclen = fast_expansion_sum_zeroelim(8, bxcy, 8, cxby, bc);
03632   Two_Two_Product(cex, cextail, dey, deytail,
03633                   cxdy7, cxdy[6], cxdy[5], cxdy[4],
03634                   cxdy[3], cxdy[2], cxdy[1], cxdy[0]);
03635   cxdy[7] = cxdy7;
03636   negate = -cey;
03637   negatetail = -ceytail;
03638   Two_Two_Product(dex, dextail, negate, negatetail,
03639                   dxcy7, dxcy[6], dxcy[5], dxcy[4],
03640                   dxcy[3], dxcy[2], dxcy[1], dxcy[0]);
03641   dxcy[7] = dxcy7;
03642   cdlen = fast_expansion_sum_zeroelim(8, cxdy, 8, dxcy, cd);
03643   Two_Two_Product(dex, dextail, aey, aeytail,
03644                   dxay7, dxay[6], dxay[5], dxay[4],
03645                   dxay[3], dxay[2], dxay[1], dxay[0]);
03646   dxay[7] = dxay7;
03647   negate = -dey;
03648   negatetail = -deytail;
03649   Two_Two_Product(aex, aextail, negate, negatetail,
03650                   axdy7, axdy[6], axdy[5], axdy[4],
03651                   axdy[3], axdy[2], axdy[1], axdy[0]);
03652   axdy[7] = axdy7;
03653   dalen = fast_expansion_sum_zeroelim(8, dxay, 8, axdy, da);
03654   Two_Two_Product(aex, aextail, cey, ceytail,
03655                   axcy7, axcy[6], axcy[5], axcy[4],
03656                   axcy[3], axcy[2], axcy[1], axcy[0]);
03657   axcy[7] = axcy7;
03658   negate = -aey;
03659   negatetail = -aeytail;
03660   Two_Two_Product(cex, cextail, negate, negatetail,
03661                   cxay7, cxay[6], cxay[5], cxay[4],
03662                   cxay[3], cxay[2], cxay[1], cxay[0]);
03663   cxay[7] = cxay7;
03664   aclen = fast_expansion_sum_zeroelim(8, axcy, 8, cxay, ac);
03665   Two_Two_Product(bex, bextail, dey, deytail,
03666                   bxdy7, bxdy[6], bxdy[5], bxdy[4],
03667                   bxdy[3], bxdy[2], bxdy[1], bxdy[0]);
03668   bxdy[7] = bxdy7;
03669   negate = -bey;
03670   negatetail = -beytail;
03671   Two_Two_Product(dex, dextail, negate, negatetail,
03672                   dxby7, dxby[6], dxby[5], dxby[4],
03673                   dxby[3], dxby[2], dxby[1], dxby[0]);
03674   dxby[7] = dxby7;
03675   bdlen = fast_expansion_sum_zeroelim(8, bxdy, 8, dxby, bd);
03676 
03677   temp32alen = scale_expansion_zeroelim(cdlen, cd, -bez, temp32a);
03678   temp32blen = scale_expansion_zeroelim(cdlen, cd, -beztail, temp32b);
03679   temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03680                                            temp32blen, temp32b, temp64a);
03681   temp32alen = scale_expansion_zeroelim(bdlen, bd, cez, temp32a);
03682   temp32blen = scale_expansion_zeroelim(bdlen, bd, ceztail, temp32b);
03683   temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03684                                            temp32blen, temp32b, temp64b);
03685   temp32alen = scale_expansion_zeroelim(bclen, bc, -dez, temp32a);
03686   temp32blen = scale_expansion_zeroelim(bclen, bc, -deztail, temp32b);
03687   temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03688                                            temp32blen, temp32b, temp64c);
03689   temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a,
03690                                            temp64blen, temp64b, temp128);
03691   temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c,
03692                                            temp128len, temp128, temp192);
03693   xlen = scale_expansion_zeroelim(temp192len, temp192, aex, detx);
03694   xxlen = scale_expansion_zeroelim(xlen, detx, aex, detxx);
03695   xtlen = scale_expansion_zeroelim(temp192len, temp192, aextail, detxt);
03696   xxtlen = scale_expansion_zeroelim(xtlen, detxt, aex, detxxt);
03697   for (i = 0; i < xxtlen; i++) {
03698     detxxt[i] *= 2.0;
03699   }
03700   xtxtlen = scale_expansion_zeroelim(xtlen, detxt, aextail, detxtxt);
03701   x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
03702   x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
03703   ylen = scale_expansion_zeroelim(temp192len, temp192, aey, dety);
03704   yylen = scale_expansion_zeroelim(ylen, dety, aey, detyy);
03705   ytlen = scale_expansion_zeroelim(temp192len, temp192, aeytail, detyt);
03706   yytlen = scale_expansion_zeroelim(ytlen, detyt, aey, detyyt);
03707   for (i = 0; i < yytlen; i++) {
03708     detyyt[i] *= 2.0;
03709   }
03710   ytytlen = scale_expansion_zeroelim(ytlen, detyt, aeytail, detytyt);
03711   y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
03712   y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
03713   zlen = scale_expansion_zeroelim(temp192len, temp192, aez, detz);
03714   zzlen = scale_expansion_zeroelim(zlen, detz, aez, detzz);
03715   ztlen = scale_expansion_zeroelim(temp192len, temp192, aeztail, detzt);
03716   zztlen = scale_expansion_zeroelim(ztlen, detzt, aez, detzzt);
03717   for (i = 0; i < zztlen; i++) {
03718     detzzt[i] *= 2.0;
03719   }
03720   ztztlen = scale_expansion_zeroelim(ztlen, detzt, aeztail, detztzt);
03721   z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1);
03722   z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2);
03723   xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy);
03724   alen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, adet);
03725 
03726   temp32alen = scale_expansion_zeroelim(dalen, da, cez, temp32a);
03727   temp32blen = scale_expansion_zeroelim(dalen, da, ceztail, temp32b);
03728   temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03729                                            temp32blen, temp32b, temp64a);
03730   temp32alen = scale_expansion_zeroelim(aclen, ac, dez, temp32a);
03731   temp32blen = scale_expansion_zeroelim(aclen, ac, deztail, temp32b);
03732   temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03733                                            temp32blen, temp32b, temp64b);
03734   temp32alen = scale_expansion_zeroelim(cdlen, cd, aez, temp32a);
03735   temp32blen = scale_expansion_zeroelim(cdlen, cd, aeztail, temp32b);
03736   temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03737                                            temp32blen, temp32b, temp64c);
03738   temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a,
03739                                            temp64blen, temp64b, temp128);
03740   temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c,
03741                                            temp128len, temp128, temp192);
03742   xlen = scale_expansion_zeroelim(temp192len, temp192, bex, detx);
03743   xxlen = scale_expansion_zeroelim(xlen, detx, bex, detxx);
03744   xtlen = scale_expansion_zeroelim(temp192len, temp192, bextail, detxt);
03745   xxtlen = scale_expansion_zeroelim(xtlen, detxt, bex, detxxt);
03746   for (i = 0; i < xxtlen; i++) {
03747     detxxt[i] *= 2.0;
03748   }
03749   xtxtlen = scale_expansion_zeroelim(xtlen, detxt, bextail, detxtxt);
03750   x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
03751   x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
03752   ylen = scale_expansion_zeroelim(temp192len, temp192, bey, dety);
03753   yylen = scale_expansion_zeroelim(ylen, dety, bey, detyy);
03754   ytlen = scale_expansion_zeroelim(temp192len, temp192, beytail, detyt);
03755   yytlen = scale_expansion_zeroelim(ytlen, detyt, bey, detyyt);
03756   for (i = 0; i < yytlen; i++) {
03757     detyyt[i] *= 2.0;
03758   }
03759   ytytlen = scale_expansion_zeroelim(ytlen, detyt, beytail, detytyt);
03760   y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
03761   y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
03762   zlen = scale_expansion_zeroelim(temp192len, temp192, bez, detz);
03763   zzlen = scale_expansion_zeroelim(zlen, detz, bez, detzz);
03764   ztlen = scale_expansion_zeroelim(temp192len, temp192, beztail, detzt);
03765   zztlen = scale_expansion_zeroelim(ztlen, detzt, bez, detzzt);
03766   for (i = 0; i < zztlen; i++) {
03767     detzzt[i] *= 2.0;
03768   }
03769   ztztlen = scale_expansion_zeroelim(ztlen, detzt, beztail, detztzt);
03770   z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1);
03771   z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2);
03772   xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy);
03773   blen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, bdet);
03774 
03775   temp32alen = scale_expansion_zeroelim(ablen, ab, -dez, temp32a);
03776   temp32blen = scale_expansion_zeroelim(ablen, ab, -deztail, temp32b);
03777   temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03778                                            temp32blen, temp32b, temp64a);
03779   temp32alen = scale_expansion_zeroelim(bdlen, bd, -aez, temp32a);
03780   temp32blen = scale_expansion_zeroelim(bdlen, bd, -aeztail, temp32b);
03781   temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03782                                            temp32blen, temp32b, temp64b);
03783   temp32alen = scale_expansion_zeroelim(dalen, da, -bez, temp32a);
03784   temp32blen = scale_expansion_zeroelim(dalen, da, -beztail, temp32b);
03785   temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03786                                            temp32blen, temp32b, temp64c);
03787   temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a,
03788                                            temp64blen, temp64b, temp128);
03789   temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c,
03790                                            temp128len, temp128, temp192);
03791   xlen = scale_expansion_zeroelim(temp192len, temp192, cex, detx);
03792   xxlen = scale_expansion_zeroelim(xlen, detx, cex, detxx);
03793   xtlen = scale_expansion_zeroelim(temp192len, temp192, cextail, detxt);
03794   xxtlen = scale_expansion_zeroelim(xtlen, detxt, cex, detxxt);
03795   for (i = 0; i < xxtlen; i++) {
03796     detxxt[i] *= 2.0;
03797   }
03798   xtxtlen = scale_expansion_zeroelim(xtlen, detxt, cextail, detxtxt);
03799   x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
03800   x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
03801   ylen = scale_expansion_zeroelim(temp192len, temp192, cey, dety);
03802   yylen = scale_expansion_zeroelim(ylen, dety, cey, detyy);
03803   ytlen = scale_expansion_zeroelim(temp192len, temp192, ceytail, detyt);
03804   yytlen = scale_expansion_zeroelim(ytlen, detyt, cey, detyyt);
03805   for (i = 0; i < yytlen; i++) {
03806     detyyt[i] *= 2.0;
03807   }
03808   ytytlen = scale_expansion_zeroelim(ytlen, detyt, ceytail, detytyt);
03809   y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
03810   y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
03811   zlen = scale_expansion_zeroelim(temp192len, temp192, cez, detz);
03812   zzlen = scale_expansion_zeroelim(zlen, detz, cez, detzz);
03813   ztlen = scale_expansion_zeroelim(temp192len, temp192, ceztail, detzt);
03814   zztlen = scale_expansion_zeroelim(ztlen, detzt, cez, detzzt);
03815   for (i = 0; i < zztlen; i++) {
03816     detzzt[i] *= 2.0;
03817   }
03818   ztztlen = scale_expansion_zeroelim(ztlen, detzt, ceztail, detztzt);
03819   z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1);
03820   z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2);
03821   xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy);
03822   clen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, cdet);
03823 
03824   temp32alen = scale_expansion_zeroelim(bclen, bc, aez, temp32a);
03825   temp32blen = scale_expansion_zeroelim(bclen, bc, aeztail, temp32b);
03826   temp64alen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03827                                            temp32blen, temp32b, temp64a);
03828   temp32alen = scale_expansion_zeroelim(aclen, ac, -bez, temp32a);
03829   temp32blen = scale_expansion_zeroelim(aclen, ac, -beztail, temp32b);
03830   temp64blen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03831                                            temp32blen, temp32b, temp64b);
03832   temp32alen = scale_expansion_zeroelim(ablen, ab, cez, temp32a);
03833   temp32blen = scale_expansion_zeroelim(ablen, ab, ceztail, temp32b);
03834   temp64clen = fast_expansion_sum_zeroelim(temp32alen, temp32a,
03835                                            temp32blen, temp32b, temp64c);
03836   temp128len = fast_expansion_sum_zeroelim(temp64alen, temp64a,
03837                                            temp64blen, temp64b, temp128);
03838   temp192len = fast_expansion_sum_zeroelim(temp64clen, temp64c,
03839                                            temp128len, temp128, temp192);
03840   xlen = scale_expansion_zeroelim(temp192len, temp192, dex, detx);
03841   xxlen = scale_expansion_zeroelim(xlen, detx, dex, detxx);
03842   xtlen = scale_expansion_zeroelim(temp192len, temp192, dextail, detxt);
03843   xxtlen = scale_expansion_zeroelim(xtlen, detxt, dex, detxxt);
03844   for (i = 0; i < xxtlen; i++) {
03845     detxxt[i] *= 2.0;
03846   }
03847   xtxtlen = scale_expansion_zeroelim(xtlen, detxt, dextail, detxtxt);
03848   x1len = fast_expansion_sum_zeroelim(xxlen, detxx, xxtlen, detxxt, x1);
03849   x2len = fast_expansion_sum_zeroelim(x1len, x1, xtxtlen, detxtxt, x2);
03850   ylen = scale_expansion_zeroelim(temp192len, temp192, dey, dety);
03851   yylen = scale_expansion_zeroelim(ylen, dety, dey, detyy);
03852   ytlen = scale_expansion_zeroelim(temp192len, temp192, deytail, detyt);
03853   yytlen = scale_expansion_zeroelim(ytlen, detyt, dey, detyyt);
03854   for (i = 0; i < yytlen; i++) {
03855     detyyt[i] *= 2.0;
03856   }
03857   ytytlen = scale_expansion_zeroelim(ytlen, detyt, deytail, detytyt);
03858   y1len = fast_expansion_sum_zeroelim(yylen, detyy, yytlen, detyyt, y1);
03859   y2len = fast_expansion_sum_zeroelim(y1len, y1, ytytlen, detytyt, y2);
03860   zlen = scale_expansion_zeroelim(temp192len, temp192, dez, detz);
03861   zzlen = scale_expansion_zeroelim(zlen, detz, dez, detzz);
03862   ztlen = scale_expansion_zeroelim(temp192len, temp192, deztail, detzt);
03863   zztlen = scale_expansion_zeroelim(ztlen, detzt, dez, detzzt);
03864   for (i = 0; i < zztlen; i++) {
03865     detzzt[i] *= 2.0;
03866   }
03867   ztztlen = scale_expansion_zeroelim(ztlen, detzt, deztail, detztzt);
03868   z1len = fast_expansion_sum_zeroelim(zzlen, detzz, zztlen, detzzt, z1);
03869   z2len = fast_expansion_sum_zeroelim(z1len, z1, ztztlen, detztzt, z2);
03870   xylen = fast_expansion_sum_zeroelim(x2len, x2, y2len, y2, detxy);
03871   dlen = fast_expansion_sum_zeroelim(z2len, z2, xylen, detxy, ddet);
03872 
03873   ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
03874   cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
03875   deterlen = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, deter);
03876 
03877   return deter[deterlen - 1];
03878 }
03879 
03880 REAL insphereadapt(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe,
03881                    REAL permanent)
03882 {
03883   INEXACT REAL aex, bex, cex, dex, aey, bey, cey, dey, aez, bez, cez, dez;
03884   REAL det, errbound;
03885 
03886   INEXACT REAL aexbey1, bexaey1, bexcey1, cexbey1;
03887   INEXACT REAL cexdey1, dexcey1, dexaey1, aexdey1;
03888   INEXACT REAL aexcey1, cexaey1, bexdey1, dexbey1;
03889   REAL aexbey0, bexaey0, bexcey0, cexbey0;
03890   REAL cexdey0, dexcey0, dexaey0, aexdey0;
03891   REAL aexcey0, cexaey0, bexdey0, dexbey0;
03892   REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4];
03893   INEXACT REAL ab3, bc3, cd3, da3, ac3, bd3;
03894   REAL abeps, bceps, cdeps, daeps, aceps, bdeps;
03895   REAL temp8a[8], temp8b[8], temp8c[8], temp16[16], temp24[24], temp48[48];
03896   int temp8alen, temp8blen, temp8clen, temp16len, temp24len, temp48len;
03897   REAL xdet[96], ydet[96], zdet[96], xydet[192];
03898   int xlen, ylen, zlen, xylen;
03899   REAL adet[288], bdet[288], cdet[288], ddet[288];
03900   int alen, blen, clen, dlen;
03901   REAL abdet[576], cddet[576];
03902   int ablen, cdlen;
03903   REAL fin1[1152];
03904   int finlength;
03905 
03906   REAL aextail, bextail, cextail, dextail;
03907   REAL aeytail, beytail, ceytail, deytail;
03908   REAL aeztail, beztail, ceztail, deztail;
03909 
03910   INEXACT REAL bvirt;
03911   REAL avirt, bround, around;
03912   INEXACT REAL c;
03913   INEXACT REAL abig;
03914   REAL ahi, alo, bhi, blo;
03915   REAL err1, err2, err3;
03916   INEXACT REAL _i, _j;
03917   REAL _0;
03918 
03919   aex = (REAL) (pa[0] - pe[0]);
03920   bex = (REAL) (pb[0] - pe[0]);
03921   cex = (REAL) (pc[0] - pe[0]);
03922   dex = (REAL) (pd[0] - pe[0]);
03923   aey = (REAL) (pa[1] - pe[1]);
03924   bey = (REAL) (pb[1] - pe[1]);
03925   cey = (REAL) (pc[1] - pe[1]);
03926   dey = (REAL) (pd[1] - pe[1]);
03927   aez = (REAL) (pa[2] - pe[2]);
03928   bez = (REAL) (pb[2] - pe[2]);
03929   cez = (REAL) (pc[2] - pe[2]);
03930   dez = (REAL) (pd[2] - pe[2]);
03931 
03932   Two_Product(aex, bey, aexbey1, aexbey0);
03933   Two_Product(bex, aey, bexaey1, bexaey0);
03934   Two_Two_Diff(aexbey1, aexbey0, bexaey1, bexaey0, ab3, ab[2], ab[1], ab[0]);
03935   ab[3] = ab3;
03936 
03937   Two_Product(bex, cey, bexcey1, bexcey0);
03938   Two_Product(cex, bey, cexbey1, cexbey0);
03939   Two_Two_Diff(bexcey1, bexcey0, cexbey1, cexbey0, bc3, bc[2], bc[1], bc[0]);
03940   bc[3] = bc3;
03941 
03942   Two_Product(cex, dey, cexdey1, cexdey0);
03943   Two_Product(dex, cey, dexcey1, dexcey0);
03944   Two_Two_Diff(cexdey1, cexdey0, dexcey1, dexcey0, cd3, cd[2], cd[1], cd[0]);
03945   cd[3] = cd3;
03946 
03947   Two_Product(dex, aey, dexaey1, dexaey0);
03948   Two_Product(aex, dey, aexdey1, aexdey0);
03949   Two_Two_Diff(dexaey1, dexaey0, aexdey1, aexdey0, da3, da[2], da[1], da[0]);
03950   da[3] = da3;
03951 
03952   Two_Product(aex, cey, aexcey1, aexcey0);
03953   Two_Product(cex, aey, cexaey1, cexaey0);
03954   Two_Two_Diff(aexcey1, aexcey0, cexaey1, cexaey0, ac3, ac[2], ac[1], ac[0]);
03955   ac[3] = ac3;
03956 
03957   Two_Product(bex, dey, bexdey1, bexdey0);
03958   Two_Product(dex, bey, dexbey1, dexbey0);
03959   Two_Two_Diff(bexdey1, bexdey0, dexbey1, dexbey0, bd3, bd[2], bd[1], bd[0]);
03960   bd[3] = bd3;
03961 
03962   temp8alen = scale_expansion_zeroelim(4, cd, bez, temp8a);
03963   temp8blen = scale_expansion_zeroelim(4, bd, -cez, temp8b);
03964   temp8clen = scale_expansion_zeroelim(4, bc, dez, temp8c);
03965   temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
03966                                           temp8blen, temp8b, temp16);
03967   temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
03968                                           temp16len, temp16, temp24);
03969   temp48len = scale_expansion_zeroelim(temp24len, temp24, aex, temp48);
03970   xlen = scale_expansion_zeroelim(temp48len, temp48, -aex, xdet);
03971   temp48len = scale_expansion_zeroelim(temp24len, temp24, aey, temp48);
03972   ylen = scale_expansion_zeroelim(temp48len, temp48, -aey, ydet);
03973   temp48len = scale_expansion_zeroelim(temp24len, temp24, aez, temp48);
03974   zlen = scale_expansion_zeroelim(temp48len, temp48, -aez, zdet);
03975   xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
03976   alen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, adet);
03977 
03978   temp8alen = scale_expansion_zeroelim(4, da, cez, temp8a);
03979   temp8blen = scale_expansion_zeroelim(4, ac, dez, temp8b);
03980   temp8clen = scale_expansion_zeroelim(4, cd, aez, temp8c);
03981   temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
03982                                           temp8blen, temp8b, temp16);
03983   temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
03984                                           temp16len, temp16, temp24);
03985   temp48len = scale_expansion_zeroelim(temp24len, temp24, bex, temp48);
03986   xlen = scale_expansion_zeroelim(temp48len, temp48, bex, xdet);
03987   temp48len = scale_expansion_zeroelim(temp24len, temp24, bey, temp48);
03988   ylen = scale_expansion_zeroelim(temp48len, temp48, bey, ydet);
03989   temp48len = scale_expansion_zeroelim(temp24len, temp24, bez, temp48);
03990   zlen = scale_expansion_zeroelim(temp48len, temp48, bez, zdet);
03991   xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
03992   blen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, bdet);
03993 
03994   temp8alen = scale_expansion_zeroelim(4, ab, dez, temp8a);
03995   temp8blen = scale_expansion_zeroelim(4, bd, aez, temp8b);
03996   temp8clen = scale_expansion_zeroelim(4, da, bez, temp8c);
03997   temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
03998                                           temp8blen, temp8b, temp16);
03999   temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
04000                                           temp16len, temp16, temp24);
04001   temp48len = scale_expansion_zeroelim(temp24len, temp24, cex, temp48);
04002   xlen = scale_expansion_zeroelim(temp48len, temp48, -cex, xdet);
04003   temp48len = scale_expansion_zeroelim(temp24len, temp24, cey, temp48);
04004   ylen = scale_expansion_zeroelim(temp48len, temp48, -cey, ydet);
04005   temp48len = scale_expansion_zeroelim(temp24len, temp24, cez, temp48);
04006   zlen = scale_expansion_zeroelim(temp48len, temp48, -cez, zdet);
04007   xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
04008   clen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, cdet);
04009 
04010   temp8alen = scale_expansion_zeroelim(4, bc, aez, temp8a);
04011   temp8blen = scale_expansion_zeroelim(4, ac, -bez, temp8b);
04012   temp8clen = scale_expansion_zeroelim(4, ab, cez, temp8c);
04013   temp16len = fast_expansion_sum_zeroelim(temp8alen, temp8a,
04014                                           temp8blen, temp8b, temp16);
04015   temp24len = fast_expansion_sum_zeroelim(temp8clen, temp8c,
04016                                           temp16len, temp16, temp24);
04017   temp48len = scale_expansion_zeroelim(temp24len, temp24, dex, temp48);
04018   xlen = scale_expansion_zeroelim(temp48len, temp48, dex, xdet);
04019   temp48len = scale_expansion_zeroelim(temp24len, temp24, dey, temp48);
04020   ylen = scale_expansion_zeroelim(temp48len, temp48, dey, ydet);
04021   temp48len = scale_expansion_zeroelim(temp24len, temp24, dez, temp48);
04022   zlen = scale_expansion_zeroelim(temp48len, temp48, dez, zdet);
04023   xylen = fast_expansion_sum_zeroelim(xlen, xdet, ylen, ydet, xydet);
04024   dlen = fast_expansion_sum_zeroelim(xylen, xydet, zlen, zdet, ddet);
04025 
04026   ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
04027   cdlen = fast_expansion_sum_zeroelim(clen, cdet, dlen, ddet, cddet);
04028   finlength = fast_expansion_sum_zeroelim(ablen, abdet, cdlen, cddet, fin1);
04029 
04030   det = estimate(finlength, fin1);
04031   errbound = isperrboundB * permanent;
04032   if ((det >= errbound) || (-det >= errbound)) {
04033     return det;
04034   }
04035 
04036   Two_Diff_Tail(pa[0], pe[0], aex, aextail);
04037   Two_Diff_Tail(pa[1], pe[1], aey, aeytail);
04038   Two_Diff_Tail(pa[2], pe[2], aez, aeztail);
04039   Two_Diff_Tail(pb[0], pe[0], bex, bextail);
04040   Two_Diff_Tail(pb[1], pe[1], bey, beytail);
04041   Two_Diff_Tail(pb[2], pe[2], bez, beztail);
04042   Two_Diff_Tail(pc[0], pe[0], cex, cextail);
04043   Two_Diff_Tail(pc[1], pe[1], cey, ceytail);
04044   Two_Diff_Tail(pc[2], pe[2], cez, ceztail);
04045   Two_Diff_Tail(pd[0], pe[0], dex, dextail);
04046   Two_Diff_Tail(pd[1], pe[1], dey, deytail);
04047   Two_Diff_Tail(pd[2], pe[2], dez, deztail);
04048   if ((aextail == 0.0) && (aeytail == 0.0) && (aeztail == 0.0)
04049       && (bextail == 0.0) && (beytail == 0.0) && (beztail == 0.0)
04050       && (cextail == 0.0) && (ceytail == 0.0) && (ceztail == 0.0)
04051       && (dextail == 0.0) && (deytail == 0.0) && (deztail == 0.0)) {
04052     return det;
04053   }
04054 
04055   errbound = isperrboundC * permanent + resulterrbound * Absolute(det);
04056   abeps = (aex * beytail + bey * aextail)
04057         - (aey * bextail + bex * aeytail);
04058   bceps = (bex * ceytail + cey * bextail)
04059         - (bey * cextail + cex * beytail);
04060   cdeps = (cex * deytail + dey * cextail)
04061         - (cey * dextail + dex * ceytail);
04062   daeps = (dex * aeytail + aey * dextail)
04063         - (dey * aextail + aex * deytail);
04064   aceps = (aex * ceytail + cey * aextail)
04065         - (aey * cextail + cex * aeytail);
04066   bdeps = (bex * deytail + dey * bextail)
04067         - (bey * dextail + dex * beytail);
04068   det += (((bex * bex + bey * bey + bez * bez)
04069            * ((cez * daeps + dez * aceps + aez * cdeps)
04070               + (ceztail * da3 + deztail * ac3 + aeztail * cd3))
04071            + (dex * dex + dey * dey + dez * dez)
04072            * ((aez * bceps - bez * aceps + cez * abeps)
04073               + (aeztail * bc3 - beztail * ac3 + ceztail * ab3)))
04074           - ((aex * aex + aey * aey + aez * aez)
04075            * ((bez * cdeps - cez * bdeps + dez * bceps)
04076               + (beztail * cd3 - ceztail * bd3 + deztail * bc3))
04077            + (cex * cex + cey * cey + cez * cez)
04078            * ((dez * abeps + aez * bdeps + bez * daeps)
04079               + (deztail * ab3 + aeztail * bd3 + beztail * da3))))
04080        + 2.0 * (((bex * bextail + bey * beytail + bez * beztail)
04081                  * (cez * da3 + dez * ac3 + aez * cd3)
04082                  + (dex * dextail + dey * deytail + dez * deztail)
04083                  * (aez * bc3 - bez * ac3 + cez * ab3))
04084                 - ((aex * aextail + aey * aeytail + aez * aeztail)
04085                  * (bez * cd3 - cez * bd3 + dez * bc3)
04086                  + (cex * cextail + cey * ceytail + cez * ceztail)
04087                  * (dez * ab3 + aez * bd3 + bez * da3)));
04088   if ((det >= errbound) || (-det >= errbound)) {
04089     return det;
04090   }
04091 
04092   return insphereexact(pa, pb, pc, pd, pe);
04093 }
04094 
04095 REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe)
04096 {
04097   REAL aex, bex, cex, dex;
04098   REAL aey, bey, cey, dey;
04099   REAL aez, bez, cez, dez;
04100   REAL aexbey, bexaey, bexcey, cexbey, cexdey, dexcey, dexaey, aexdey;
04101   REAL aexcey, cexaey, bexdey, dexbey;
04102   REAL alift, blift, clift, dlift;
04103   REAL ab, bc, cd, da, ac, bd;
04104   REAL abc, bcd, cda, dab;
04105   REAL aezplus, bezplus, cezplus, dezplus;
04106   REAL aexbeyplus, bexaeyplus, bexceyplus, cexbeyplus;
04107   REAL cexdeyplus, dexceyplus, dexaeyplus, aexdeyplus;
04108   REAL aexceyplus, cexaeyplus, bexdeyplus, dexbeyplus;
04109   REAL det;
04110   REAL permanent, errbound;
04111 
04112   aex = pa[0] - pe[0];
04113   bex = pb[0] - pe[0];
04114   cex = pc[0] - pe[0];
04115   dex = pd[0] - pe[0];
04116   aey = pa[1] - pe[1];
04117   bey = pb[1] - pe[1];
04118   cey = pc[1] - pe[1];
04119   dey = pd[1] - pe[1];
04120   aez = pa[2] - pe[2];
04121   bez = pb[2] - pe[2];
04122   cez = pc[2] - pe[2];
04123   dez = pd[2] - pe[2];
04124 
04125   aexbey = aex * bey;
04126   bexaey = bex * aey;
04127   ab = aexbey - bexaey;
04128   bexcey = bex * cey;
04129   cexbey = cex * bey;
04130   bc = bexcey - cexbey;
04131   cexdey = cex * dey;
04132   dexcey = dex * cey;
04133   cd = cexdey - dexcey;
04134   dexaey = dex * aey;
04135   aexdey = aex * dey;
04136   da = dexaey - aexdey;
04137 
04138   aexcey = aex * cey;
04139   cexaey = cex * aey;
04140   ac = aexcey - cexaey;
04141   bexdey = bex * dey;
04142   dexbey = dex * bey;
04143   bd = bexdey - dexbey;
04144 
04145   abc = aez * bc - bez * ac + cez * ab;
04146   bcd = bez * cd - cez * bd + dez * bc;
04147   cda = cez * da + dez * ac + aez * cd;
04148   dab = dez * ab + aez * bd + bez * da;
04149 
04150   alift = aex * aex + aey * aey + aez * aez;
04151   blift = bex * bex + bey * bey + bez * bez;
04152   clift = cex * cex + cey * cey + cez * cez;
04153   dlift = dex * dex + dey * dey + dez * dez;
04154 
04155   det = (dlift * abc - clift * dab) + (blift * cda - alift * bcd);
04156 
04157   aezplus = Absolute(aez);
04158   bezplus = Absolute(bez);
04159   cezplus = Absolute(cez);
04160   dezplus = Absolute(dez);
04161   aexbeyplus = Absolute(aexbey);
04162   bexaeyplus = Absolute(bexaey);
04163   bexceyplus = Absolute(bexcey);
04164   cexbeyplus = Absolute(cexbey);
04165   cexdeyplus = Absolute(cexdey);
04166   dexceyplus = Absolute(dexcey);
04167   dexaeyplus = Absolute(dexaey);
04168   aexdeyplus = Absolute(aexdey);
04169   aexceyplus = Absolute(aexcey);
04170   cexaeyplus = Absolute(cexaey);
04171   bexdeyplus = Absolute(bexdey);
04172   dexbeyplus = Absolute(dexbey);
04173   permanent = ((cexdeyplus + dexceyplus) * bezplus
04174                + (dexbeyplus + bexdeyplus) * cezplus
04175                + (bexceyplus + cexbeyplus) * dezplus)
04176             * alift
04177             + ((dexaeyplus + aexdeyplus) * cezplus
04178                + (aexceyplus + cexaeyplus) * dezplus
04179                + (cexdeyplus + dexceyplus) * aezplus)
04180             * blift
04181             + ((aexbeyplus + bexaeyplus) * dezplus
04182                + (bexdeyplus + dexbeyplus) * aezplus
04183                + (dexaeyplus + aexdeyplus) * bezplus)
04184             * clift
04185             + ((bexceyplus + cexbeyplus) * aezplus
04186                + (cexaeyplus + aexceyplus) * bezplus
04187                + (aexbeyplus + bexaeyplus) * cezplus)
04188             * dlift;
04189   errbound = isperrboundA * permanent;
04190   if ((det > errbound) || (-det > errbound)) {
04191     return det;
04192   }
04193 
04194   return insphereadapt(pa, pb, pc, pd, pe, permanent);
04195 }
04196 } //Added namespace to avoid clash with triangle