misc/predicates.c

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