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