}/*****************************************************************************//*                                                                           *//*  orient2dfast()   Approximate 2D orientation test.  Nonrobust.            *//*  orient2dexact()   Exact 2D orientation test.  Robust.                    *//*  orient2dslow()   Another exact 2D orientation test.  Robust.             *//*  orient2d()   Adaptive exact 2D orientation test.  Robust.                *//*                                                                           *//*               Return a positive value if the points pa, pb, and pc occur  *//*               in counterclockwise order; a negative value if they occur   *//*               in clockwise order; and zero if they are collinear.  The    *//*               result is also a rough approximation of twice the signed    *//*               area of the triangle defined by the three points.           *//*                                                                           *//*  Only the first and last routine should be used; the middle two are for   *//*  timings.                                                                 *//*                                                                           *//*  The last three use exact arithmetic to ensure a correct answer.  The     *//*  result returned is the determinant of a matrix.  In orient2d() only,     *//*  this determinant is computed adaptively, in the sense that exact         *//*  arithmetic is used only to the degree it is needed to ensure that the    *//*  returned value has the correct sign.  Hence, orient2d() is usually quite *//*  fast, but will run more slowly when the input points are collinear or    *//*  nearly so.                                                               *//*                                                                           *//*****************************************************************************/REAL orient2dfast(REAL *pa, REAL *pb, REAL *pc){  REAL acx, bcx, acy, bcy;  acx = pa[0] - pc[0];  bcx = pb[0] - pc[0];  acy = pa[1] - pc[1];  bcy = pb[1] - pc[1];  return acx * bcy - acy * bcx;}REAL orient2dexact(REAL *pa, REAL *pb, REAL *pc){  INEXACT REAL axby1, axcy1, bxcy1, bxay1, cxay1, cxby1;  REAL axby0, axcy0, bxcy0, bxay0, cxay0, cxby0;  REAL aterms[4], bterms[4], cterms[4];  INEXACT REAL aterms3, bterms3, cterms3;  REAL v[8], w[12];  int vlength, wlength;  INEXACT REAL bvirt;  REAL avirt, bround, around;  INEXACT REAL c;  INEXACT REAL abig;  REAL ahi, alo, bhi, blo;  REAL err1, err2, err3;  INEXACT REAL _i, _j;  REAL _0;  Two_Product(pa[0], pb[1], axby1, axby0);  Two_Product(pa[0], pc[1], axcy1, axcy0);  Two_Two_Diff(axby1, axby0, axcy1, axcy0,               aterms3, aterms[2], aterms[1], aterms[0]);  aterms[3] = aterms3;  Two_Product(pb[0], pc[1], bxcy1, bxcy0);  Two_Product(pb[0], pa[1], bxay1, bxay0);  Two_Two_Diff(bxcy1, bxcy0, bxay1, bxay0,               bterms3, bterms[2], bterms[1], bterms[0]);  bterms[3] = bterms3;  Two_Product(pc[0], pa[1], cxay1, cxay0);  Two_Product(pc[0], pb[1], cxby1, cxby0);  Two_Two_Diff(cxay1, cxay0, cxby1, cxby0,               cterms3, cterms[2], cterms[1], cterms[0]);  cterms[3] = cterms3;  vlength = fast_expansion_sum_zeroelim(4, aterms, 4, bterms, v);  wlength = fast_expansion_sum_zeroelim(vlength, v, 4, cterms, w);  return w[wlength - 1];}REAL orient2dslow(REAL *pa, REAL *pb, REAL *pc){  INEXACT REAL acx, acy, bcx, bcy;  REAL acxtail, acytail;  REAL bcxtail, bcytail;  REAL negate, negatetail;  REAL axby[8], bxay[8];  INEXACT REAL axby7, bxay7;  REAL deter[16];  int deterlen;  INEXACT REAL bvirt;  REAL avirt, bround, around;  INEXACT REAL c;  INEXACT REAL abig;  REAL a0hi, a0lo, a1hi, a1lo, bhi, blo;  REAL err1, err2, err3;  INEXACT REAL _i, _j, _k, _l, _m, _n;  REAL _0, _1, _2;  Two_Diff(pa[0], pc[0], acx, acxtail);  Two_Diff(pa[1], pc[1], acy, acytail);  Two_Diff(pb[0], pc[0], bcx, bcxtail);  Two_Diff(pb[1], pc[1], bcy, bcytail);  Two_Two_Product(acx, acxtail, bcy, bcytail,                  axby7, axby[6], axby[5], axby[4],                  axby[3], axby[2], axby[1], axby[0]);  axby[7] = axby7;  negate = -acy;  negatetail = -acytail;  Two_Two_Product(bcx, bcxtail, negate, negatetail,                  bxay7, bxay[6], bxay[5], bxay[4],                  bxay[3], bxay[2], bxay[1], bxay[0]);  bxay[7] = bxay7;  deterlen = fast_expansion_sum_zeroelim(8, axby, 8, bxay, deter);  return deter[deterlen - 1];}REAL orient2dadapt(REAL *pa, REAL *pb, REAL *pc, REAL detsum){  INEXACT REAL acx, acy, bcx, bcy;  REAL acxtail, acytail, bcxtail, bcytail;  INEXACT REAL detleft, detright;  REAL detlefttail, detrighttail;  REAL det, errbound;  REAL B[4], C1[8], C2[12], D[16];  INEXACT REAL B3;  int C1length, C2length, Dlength;  REAL u[4];  INEXACT REAL u3;  INEXACT REAL s1, t1;  REAL s0, t0;  INEXACT REAL bvirt;  REAL avirt, bround, around;  INEXACT REAL c;  INEXACT REAL abig;  REAL ahi, alo, bhi, blo;  REAL err1, err2, err3;  INEXACT REAL _i, _j;  REAL _0;  acx = (REAL) (pa[0] - pc[0]);  bcx = (REAL) (pb[0] - pc[0]);  acy = (REAL) (pa[1] - pc[1]);  bcy = (REAL) (pb[1] - pc[1]);  Two_Product(acx, bcy, detleft, detlefttail);  Two_Product(acy, bcx, detright, detrighttail);  Two_Two_Diff(detleft, detlefttail, detright, detrighttail,               B3, B[2], B[1], B[0]);  B[3] = B3;  det = estimate(4, B);  errbound = ccwerrboundB * detsum;  if ((det >= errbound) || (-det >= errbound)) {    return det;  }  Two_Diff_Tail(pa[0], pc[0], acx, acxtail);  Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail);  Two_Diff_Tail(pa[1], pc[1], acy, acytail);  Two_Diff_Tail(pb[1], pc[1], bcy, bcytail);  if ((acxtail == 0.0) && (acytail == 0.0)      && (bcxtail == 0.0) && (bcytail == 0.0)) {    return det;  }  errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det);  det += (acx * bcytail + bcy * acxtail)       - (acy * bcxtail + bcx * acytail);  if ((det >= errbound) || (-det >= errbound)) {    return det;  }  Two_Product(acxtail, bcy, s1, s0);  Two_Product(acytail, bcx, t1, t0);  Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);  u[3] = u3;  C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1);  Two_Product(acx, bcytail, s1, s0);  Two_Product(acy, bcxtail, t1, t0);  Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);  u[3] = u3;  C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2);  Two_Product(acxtail, bcytail, s1, s0);  Two_Product(acytail, bcxtail, t1, t0);  Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);  u[3] = u3;  Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D);  return(D[Dlength - 1]);}REAL orient2d(REAL *pa, REAL *pb, REAL *pc){  REAL detleft, detright, det;  REAL detsum, errbound;  detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]);  detright = (pa[1] - pc[1]) * (pb[0] - pc[0]);  det = detleft - detright;  if (detleft > 0.0) {    if (detright <= 0.0) {      return det;    } else {      detsum = detleft + detright;    }  } else if (detleft < 0.0) {    if (detright >= 0.0) {      return det;    } else {      detsum = -detleft - detright;    }  } else {    return det;  }  errbound = ccwerrboundA * detsum;  if ((det >= errbound) || (-det >= errbound)) {    return det;  }  return orient2dadapt(pa, pb, pc, detsum);}/*****************************************************************************//*                                                                           *//*  orient3dfast()   Approximate 3D orientation test.  Nonrobust.            *//*  orient3dexact()   Exact 3D orientation test.  Robust.                    *//*  orient3dslow()   Another exact 3D orientation test.  Robust.             *//*  orient3d()   Adaptive exact 3D orientation test.  Robust.                *//*                                                                           *//*               Return a positive value if the point pd lies below the      *//*               plane passing through pa, pb, and pc; "below" is defined so *//*               that pa, pb, and pc appear in counterclockwise order when   *//*               viewed from above the plane.  Returns a negative value if   *//*               pd lies above the plane.  Returns zero if the points are    *//*               coplanar.  The result is also a rough approximation of six  *//*               times the signed volume of the tetrahedron defined by the   *//*               four points.                                                *//*                                                                           *//*  Only the first and last routine should be used; the middle two are for   *//*  timings.                                                                 *//*                                                                           *//*  The last three use exact arithmetic to ensure a correct answer.  The     *//*  result returned is the determinant of a matrix.  In orient3d() only,     *//*  this determinant is computed adaptively, in the sense that exact         *//*  arithmetic is used only to the degree it is needed to ensure that the    *//*  returned value has the correct sign.  Hence, orient3d() is usually quite *//*  fast, but will run more slowly when the input points are coplanar or     *//*  nearly so.                                                               *//*                                                                           *//*****************************************************************************/REAL orient3dfast(REAL *pa, REAL *pb, REAL *pc, REAL *pd){  REAL adx, bdx, cdx;  REAL ady, bdy, cdy;  REAL adz, bdz, cdz;  adx = pa[0] - pd[0];  bdx = pb[0] - pd[0];  cdx = pc[0] - pd[0];  ady = pa[1] - pd[1];  bdy = pb[1] - pd[1];  cdy = pc[1] - pd[1];  adz = pa[2] - pd[2];  bdz = pb[2] - pd[2];  cdz = pc[2] - pd[2];  return adx * (bdy * cdz - bdz * cdy)       + bdx * (cdy * adz - cdz * ady)       + cdx * (ady * bdz - adz * bdy);}REAL orient3dexact(REAL *pa, REAL *pb, REAL *pc, REAL *pd){  INEXACT REAL axby1, bxcy1, cxdy1, dxay1, axcy1, bxdy1;  INEXACT REAL bxay1, cxby1, dxcy1, axdy1, cxay1, dxby1;  REAL axby0, bxcy0, cxdy0, dxay0, axcy0, bxdy0;  REAL bxay0, cxby0, dxcy0, axdy0, cxay0, dxby0;  REAL ab[4], bc[4], cd[4], da[4], ac[4], bd[4];  REAL temp8[8];  int templen;  REAL abc[12], bcd[12], cda[12], dab[12];  int abclen, bcdlen, cdalen, dablen;  REAL adet[24], bdet[24], cdet[24], ddet[24];  int alen, blen, clen, dlen;  REAL abdet[48], cddet[48];  int ablen, cdlen;  REAL deter[96];  int deterlen;  int i;  INEXACT REAL bvirt;  REAL avirt, bround, around;  INEXACT REAL c;  INEXACT REAL abig;  REAL ahi, alo, bhi, blo;  REAL err1, err2, err3;  INEXACT REAL _i, _j;  REAL _0;  Two_Product(pa[0], pb[1], axby1, axby0);  Two_Product(pb[0], pa[1], bxay1, bxay0);  Two_Two_Diff(axby1, axby0, bxay1, bxay0, ab[3], ab[2], ab[1], ab[0]);  Two_Product(pb[0], pc[1], bxcy1, bxcy0);  Two_Product(pc[0], pb[1], cxby1, cxby0);  Two_Two_Diff(bxcy1, bxcy0, cxby1, cxby0, bc[3], bc[2], bc[1], bc[0]);  Two_Product(pc[0], pd[1], cxdy1, cxdy0);  Two_Product(pd[0], pc[1], dxcy1, dxcy0);  Two_Two_Diff(cxdy1, cxdy0, dxcy1, dxcy0, cd[3], cd[2], cd[1], cd[0]);  Two_Product(pd[0], pa[1], dxay1, dxay0);  Two_Product(pa[0], pd[1], axdy1, axdy0);  Two_Two_Diff(dxay1, dxay0, axdy1, axdy0, da[3], da[2], da[1], da[0]);  Two_Product(pa[0], pc[1], axcy1, axcy0);  Two_Product(pc[0], pa[1], cxay1, cxay0);  Two_Two_Diff(axcy1, axcy0, cxay1, cxay0, ac[3], ac[2], ac[1], ac[0]);  Two_Product(pb[0], pd[1], bxdy1, bxdy0);  Two_Product(pd[0], pb[1], dxby1, dxby0);  Two_Two_Diff(bxdy1, bxdy0, dxby1, dxby0, bd[3], bd[2], bd[1], bd[0]);  templen = fast_expansion_sum_zeroelim(4, cd, 4, da, temp8);  cdalen = fast_expansion_sum_zeroelim(templen, temp8, 4, ac, cda);  templen = fast_expansion_sum_zeroelim(4, da, 4, ab, temp8);  dablen = fast_expansion_sum_zeroelim(templen, temp8, 4, bd, dab);  for (i = 0; i < 4; i++) {    bd[i] = -bd[i];    ac[i] = -ac[i];  }  templen = fast_expansion_sum_zeroelim(4, ab, 4, bc, temp8);  abclen = fast_expansion_sum_