epsilon *= half;    if (every_other) {      splitter *= 2.0;    }    every_other = !every_other;    check = 1.0 + epsilon;  } while ((check != 1.0) && (check != lastcheck));  splitter += 1.0;  /* Error bounds for orientation and incircle tests. */  resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;  ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon;  ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon;  ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon;  o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon;  o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon;  o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;  iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon;  iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon;  iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon;  isperrboundA = (16.0 + 224.0 * epsilon) * epsilon;  isperrboundB = (5.0 + 72.0 * epsilon) * epsilon;  isperrboundC = (71.0 + 1408.0 * epsilon) * epsilon * epsilon;  return epsilon; /* Added by H. Si 30 Juli, 2004. */}/*****************************************************************************//*                                                                           *//*  grow_expansion()   Add a scalar to an expansion.                         *//*                                                                           *//*  Sets h = e + b.  See the long version of my paper for details.           *//*                                                                           *//*  Maintains the nonoverlapping property.  If round-to-even is used (as     *//*  with IEEE 754), maintains the strongly nonoverlapping and nonadjacent    *//*  properties as well.  (That is, if e has one of these properties, so      *//*  will h.)                                                                 *//*                                                                           *//*****************************************************************************/int grow_expansion(int elen, REAL *e, REAL b, REAL *h)/* e and h can be the same. */{  REAL Q;  INEXACT REAL Qnew;  int eindex;  REAL enow;  INEXACT REAL bvirt;  REAL avirt, bround, around;  Q = b;  for (eindex = 0; eindex < elen; eindex++) {    enow = e[eindex];    Two_Sum(Q, enow, Qnew, h[eindex]);    Q = Qnew;  }  h[eindex] = Q;  return eindex + 1;}/*****************************************************************************//*                                                                           *//*  grow_expansion_zeroelim()   Add a scalar to an expansion, eliminating    *//*                              zero components from the output expansion.   *//*                                                                           *//*  Sets h = e + b.  See the long version of my paper for details.           *//*                                                                           *//*  Maintains the nonoverlapping property.  If round-to-even is used (as     *//*  with IEEE 754), maintains the strongly nonoverlapping and nonadjacent    *//*  properties as well.  (That is, if e has one of these properties, so      *//*  will h.)                                                                 *//*                                                                           *//*****************************************************************************/int grow_expansion_zeroelim(int elen, REAL *e, REAL b, REAL *h)/* e and h can be the same. */{  REAL Q, hh;  INEXACT REAL Qnew;  int eindex, hindex;  REAL enow;  INEXACT REAL bvirt;  REAL avirt, bround, around;  hindex = 0;  Q = b;  for (eindex = 0; eindex < elen; eindex++) {    enow = e[eindex];    Two_Sum(Q, enow, Qnew, hh);    Q = Qnew;    if (hh != 0.0) {      h[hindex++] = hh;    }  }  if ((Q != 0.0) || (hindex == 0)) {    h[hindex++] = Q;  }  return hindex;}/*****************************************************************************//*                                                                           *//*  expansion_sum()   Sum two expansions.                                    *//*                                                                           *//*  Sets h = e + f.  See the long version of my paper for details.           *//*                                                                           *//*  Maintains the nonoverlapping property.  If round-to-even is used (as     *//*  with IEEE 754), maintains the nonadjacent property as well.  (That is,   *//*  if e has one of these properties, so will h.)  Does NOT maintain the     *//*  strongly nonoverlapping property.                                        *//*                                                                           *//*****************************************************************************/int expansion_sum(int elen, REAL *e, int flen, REAL *f, REAL *h)/* e and h can be the same, but f and h cannot. */{  REAL Q;  INEXACT REAL Qnew;  int findex, hindex, hlast;  REAL hnow;  INEXACT REAL bvirt;  REAL avirt, bround, around;  Q = f[0];  for (hindex = 0; hindex < elen; hindex++) {    hnow = e[hindex];    Two_Sum(Q, hnow, Qnew, h[hindex]);    Q = Qnew;  }  h[hindex] = Q;  hlast = hindex;  for (findex = 1; findex < flen; findex++) {    Q = f[findex];    for (hindex = findex; hindex <= hlast; hindex++) {      hnow = h[hindex];      Two_Sum(Q, hnow, Qnew, h[hindex]);      Q = Qnew;    }    h[++hlast] = Q;  }  return hlast + 1;}/*****************************************************************************//*                                                                           *//*  expansion_sum_zeroelim1()   Sum two expansions, eliminating zero         *//*                              components from the output expansion.        *//*                                                                           *//*  Sets h = e + f.  See the long version of my paper for details.           *//*                                                                           *//*  Maintains the nonoverlapping property.  If round-to-even is used (as     *//*  with IEEE 754), maintains the nonadjacent property as well.  (That is,   *//*  if e has one of these properties, so will h.)  Does NOT maintain the     *//*  strongly nonoverlapping property.                                        *//*                                                                           *//*****************************************************************************/int expansion_sum_zeroelim1(int elen, REAL *e, int flen, REAL *f, REAL *h)/* e and h can be the same, but f and h cannot. */{  REAL Q;  INEXACT REAL Qnew;  int index, findex, hindex, hlast;  REAL hnow;  INEXACT REAL bvirt;  REAL avirt, bround, around;  Q = f[0];  for (hindex = 0; hindex < elen; hindex++) {    hnow = e[hindex];    Two_Sum(Q, hnow, Qnew, h[hindex]);    Q = Qnew;  }  h[hindex] = Q;  hlast = hindex;  for (findex = 1; findex < flen; findex++) {    Q = f[findex];    for (hindex = findex; hindex <= hlast; hindex++) {      hnow = h[hindex];      Two_Sum(Q, hnow, Qnew, h[hindex]);      Q = Qnew;    }    h[++hlast] = Q;  }  hindex = -1;  for (index = 0; index <= hlast; index++) {    hnow = h[index];    if (hnow != 0.0) {      h[++hindex] = hnow;    }  }  if (hindex == -1) {    return 1;  } else {    return hindex + 1;  }}/*****************************************************************************//*                                                                           *//*  expansion_sum_zeroelim2()   Sum two expansions, eliminating zero         *//*                              components from the output expansion.        *//*                                                                           *//*  Sets h = e + f.  See the long version of my paper for details.           *//*                                                                           *//*  Maintains the nonoverlapping property.  If round-to-even is used (as     *//*  with IEEE 754), maintains the nonadjacent property as well.  (That is,   *//*  if e has one of these properties, so will h.)  Does NOT maintain the     *//*  strongly nonoverlapping property.                                        *//*                                                                           *//*****************************************************************************/int expansion_sum_zeroelim2(int elen, REAL *e, int flen, REAL *f, REAL *h)/* e and h can be the same, but f and h cannot. */{  REAL Q, hh;  INEXACT REAL Qnew;  int eindex, findex, hindex, hlast;  REAL enow;  INEXACT REAL bvirt;  REAL avirt, bround, around;  hindex = 0;  Q = f[0];  for (eindex = 0; eindex < elen; eindex++) {    enow = e[eindex];    Two_Sum(Q, enow, Qnew, hh);    Q = Qnew;    if (hh != 0.0) {      h[hindex++] = hh;    }  }  h[hindex] = Q;  hlast = hindex;  for (findex = 1; findex < flen; findex++) {    hindex = 0;    Q = f[findex];    for (eindex = 0; eindex <= hlast; eindex++) {      enow = h[eindex];      Two_Sum(Q, enow, Qnew, hh);      Q = Qnew;      if (hh != 0) {        h[hindex++] = hh;      }    }    h[hindex] = Q;    hlast = hindex;  }  return hlast + 1;}/*****************************************************************************//*                                                                           *//*  fast_expansion_sum()   Sum two expansions.                               *//*                                                                           *//*  Sets h = e + f.  See the long version of my paper for details.           *//*                                                                           *//*  If round-to-even is used (as with IEEE 754), maintains the strongly      *//*  nonoverlapping property.  (That is, if e is strongly nonoverlapping, h   *//*  will be also.)  Does NOT maintain the nonoverlapping or nonadjacent      *//*  properties.                                                              *//*                                                                           *//*****************************************************************************/int fast_expansion_sum(int elen, REAL *e, int flen, REAL *f, REAL *h)/* h cannot be e or f. */{  REAL Q;  INEXACT REAL Qnew;  INEXACT REAL bvirt;  REAL avirt, bround, around;  int eindex, findex, hindex;  REAL enow, fnow;  enow = e[0];  fnow = f[0];  eindex = findex = 0;  if ((fnow > enow) == (fnow > -enow)) {    Q = enow;    enow = e[++eindex];  } else {    Q = fnow;    fnow = f[++findex];  }  hindex = 0;  if ((eindex < elen) && (findex < flen)) {    if ((fnow > enow) == (fnow > -enow)) {      Fast_Two_Sum(enow, Q, Qnew, h[0]);      enow = e[++eindex];    } else {      Fast_Two_Sum(fnow, Q, Qnew, h[0]);      fnow = f[++findex];    }    Q = Qnew;    hindex = 1;    while ((eindex < elen) && (findex < flen)) {      if ((fnow > enow) == (fnow > -enow)) {        Two_Sum(Q, enow, Qnew, h[hindex]);        enow = e[++eindex];      } else {        Two_Sum(Q, fnow, Qnew, h[hindex]);        fnow = f[++findex];      }      Q = Qnew;      hindex++;    }  }  while (eindex < elen) {    Two_Sum(Q, enow, Qnew, h[hindex]);    enow = e[++eindex];    Q = Qnew;    hindex++;  }  while (findex < flen) {    Two_Sum(Q, fnow, Qnew, h[hindex]);    fnow = f[++findex];    Q = Qnew;    hindex++;  }  h[hindex] = Q;  return hindex + 1;}/*****************************************************************************//*                                                                           *//*  fast_expansion_sum_zeroelim()   Sum two expansions, eliminating zero     *//*                                  components from the output expansion.    *//*                                                                           *//*  Sets h = e + f.  See the long version of my paper for details.           *//*                                                                           *//*  If round-to-even is used (as with IEEE 754), maintains the strongly      *//*  nonoverlapping property.  (That is, if e is strongly nonoverlapping, h   *//*  will be also.)  Does NOT maintain the nonoverlapping or nonadjacent      *//*  properties.                                                              *//*                                                                           *//*****************************************************************************/int fast_expansion_sum_zeroelim(int elen, REAL *e, int flen, REAL *f, REAL *h)/* h cannot be e or f. */{  REAL Q;  INEXACT REAL Qnew;  INEXACT REAL hh;  INEXACT REAL bvirt;  REAL avirt, bround, around;  int eindex, findex, hindex;  REAL enow, fnow;  enow = e[0];  fnow = f[0];