From f07cb1ded76ed428e1296f347a3fb3bab3f7491c Mon Sep 17 00:00:00 2001 From: "SND\\wo80_cp" Date: Tue, 20 Aug 2013 14:52:38 +0000 Subject: [PATCH] Added Shewchuk's robust geometric predicates Fixed bug in sweepline algorithm git-svn-id: https://triangle.svn.codeplex.com/svn@73454 0e2699bc-83d4-4a8f-98e7-55e24ab8c7a5 --- Triangle.NET/Triangle/Algorithm/SweepLine.cs | 69 +- Triangle.NET/Triangle/Mesh.cs | 3 +- Triangle.NET/Triangle/NewLocation.cs | 2 +- Triangle.NET/Triangle/Primitives.cs | 983 +++++++++++++++++-- 4 files changed, 914 insertions(+), 143 deletions(-) diff --git a/Triangle.NET/Triangle/Algorithm/SweepLine.cs b/Triangle.NET/Triangle/Algorithm/SweepLine.cs index f299469..aa0fd11 100644 --- a/Triangle.NET/Triangle/Algorithm/SweepLine.cs +++ b/Triangle.NET/Triangle/Algorithm/SweepLine.cs @@ -31,7 +31,7 @@ namespace TriangleNet.Algorithm } Mesh mesh; - double xminextreme; // Nonexistent x value used as a flag in sweepline. + double xminextreme; // Nonexistent x value used as a flag in sweepline. List splaynodes; #region Heap @@ -177,7 +177,6 @@ namespace TriangleNet.Algorithm evt.xkey = thisvertex.x; evt.ykey = thisvertex.y; HeapInsert(eventheap, i++, evt); - } } @@ -357,7 +356,23 @@ namespace TriangleNet.Algorithm return newsplaynode; } - #endregion + SplayNode FrontLocate(SplayNode splayroot, Otri bottommost, Vertex searchvertex, + ref Otri searchtri, ref bool farright) + { + bool farrightflag; + + bottommost.Copy(ref searchtri); + splayroot = Splay(splayroot, searchvertex, ref searchtri); + + farrightflag = false; + while (!farrightflag && RightOfHyperbola(ref searchtri, searchvertex)) + { + searchtri.OnextSelf(); + farrightflag = searchtri.Equal(bottommost); + } + farright = farrightflag; + return splayroot; + } SplayNode CircleTopInsert(SplayNode splayroot, Otri newkey, Vertex pa, Vertex pb, Vertex pc, double topy) @@ -380,6 +395,8 @@ namespace TriangleNet.Algorithm return SplayInsert(Splay(splayroot, searchpoint, ref dummytri), newkey, searchpoint); } + #endregion + bool RightOfHyperbola(ref Otri fronttri, Point newsite) { Vertex leftvertex, rightvertex; @@ -449,24 +466,6 @@ namespace TriangleNet.Algorithm } } - SplayNode FrontLocate(SplayNode splayroot, Otri bottommost, Vertex searchvertex, - ref Otri searchtri, ref bool farright) - { - bool farrightflag; - - bottommost.Copy(ref searchtri); - splayroot = Splay(splayroot, searchvertex, ref searchtri); - - farrightflag = false; - while (!farrightflag && RightOfHyperbola(ref searchtri, searchvertex)) - { - searchtri.OnextSelf(); - farrightflag = searchtri.Equal(bottommost); - } - farright = farrightflag; - return splayroot; - } - /// /// Removes ghost triangles. /// @@ -574,7 +573,7 @@ namespace TriangleNet.Algorithm { if (heapsize == 0) { - SimpleLog.Instance.Error("Input vertices are all identical.", "SweepLine.SweepLineDelaunay()"); + SimpleLog.Instance.Error("Input vertices are all identical.", "SweepLine.Triangulate()"); throw new Exception("Input vertices are all identical."); } secondvertex = eventheap[0].vertexEvent; @@ -585,8 +584,8 @@ namespace TriangleNet.Algorithm { if (Behavior.Verbose) { - SimpleLog.Instance.Warning("A duplicate vertex appeared and was ignored.", - "SweepLine.SweepLineDelaunay().1"); + SimpleLog.Instance.Warning("A duplicate vertex appeared and was ignored (ID " + secondvertex.id + ").", + "SweepLine.Triangulate().1"); } secondvertex.type = VertexType.UndeadVertex; mesh.undeads++; @@ -641,8 +640,8 @@ namespace TriangleNet.Algorithm { if (Behavior.Verbose) { - SimpleLog.Instance.Warning("A duplicate vertex appeared and was ignored.", - "SweepLine.SweepLineDelaunay().2"); + SimpleLog.Instance.Warning("A duplicate vertex appeared and was ignored (ID " + nextvertex.id + ").", + "SweepLine.Triangulate().2"); } nextvertex.type = VertexType.UndeadVertex; mesh.undeads++; @@ -652,17 +651,15 @@ namespace TriangleNet.Algorithm { lastvertex = nextvertex; - splayroot = FrontLocate(splayroot, bottommost, nextvertex, - ref searchtri, ref farrightflag); - // - bottommost.Copy(ref searchtri); - farrightflag = false; - while (!farrightflag && RightOfHyperbola(ref searchtri, nextvertex)) - { - searchtri.OnextSelf(); - farrightflag = searchtri.Equal(bottommost); - } + splayroot = FrontLocate(splayroot, bottommost, nextvertex, ref searchtri, ref farrightflag); + //bottommost.Copy(ref searchtri); + //farrightflag = false; + //while (!farrightflag && RightOfHyperbola(ref searchtri, nextvertex)) + //{ + // searchtri.OnextSelf(); + // farrightflag = searchtri.Equal(bottommost); + //} Check4DeadEvent(ref searchtri, eventheap, ref heapsize); diff --git a/Triangle.NET/Triangle/Mesh.cs b/Triangle.NET/Triangle/Mesh.cs index aaa353f..daf00a6 100644 --- a/Triangle.NET/Triangle/Mesh.cs +++ b/Triangle.NET/Triangle/Mesh.cs @@ -2705,7 +2705,8 @@ namespace TriangleNet { if (Behavior.Verbose) { - logger.Warning("Endpoints of segments are coincident.", "Mesh.FormSkeleton()"); + logger.Warning("Endpoints of segment (IDs " + end1 + "/" + end2 + ") are coincident.", + "Mesh.FormSkeleton()"); } } else diff --git a/Triangle.NET/Triangle/NewLocation.cs b/Triangle.NET/Triangle/NewLocation.cs index cd4036e..193ea33 100644 --- a/Triangle.NET/Triangle/NewLocation.cs +++ b/Triangle.NET/Triangle/NewLocation.cs @@ -4100,7 +4100,7 @@ namespace TriangleNet if (ahead < 0.0) { // Turn around so that 'searchpoint' is to the left of the - // edge specified by 'searchtri'. + // edge specified by 'searchtri'. searchtri.SymSelf(); searchtri.Copy(ref horiz); intersect = mesh.locator.PreciseLocate(newvertex, ref horiz, false); diff --git a/Triangle.NET/Triangle/Primitives.cs b/Triangle.NET/Triangle/Primitives.cs index 6d7855a..79b8902 100644 --- a/Triangle.NET/Triangle/Primitives.cs +++ b/Triangle.NET/Triangle/Primitives.cs @@ -13,18 +13,25 @@ namespace TriangleNet using TriangleNet.Tools; /// - /// Provides some primitives regularly used in computational geometry. + /// The adaptive exact arithmetic geometric predicates implemented herein are described in + /// detail in the paper "Adaptive Precision Floating-Point Arithmetic and Fast Robust + /// Geometric Predicates." by Jonathan Richard Shewchuk, see + /// http://www.cs.cmu.edu/~quake/robust.html /// + /// + /// The macros of the original C code were automatically expanded using the Visual Studio + /// command prompt with the command "CL /P /C EXACT.C", see + /// http://msdn.microsoft.com/en-us/library/8z9z0bx6.aspx + /// public static class Primitives { - static double splitter; // Used to split double factors for exact multiplication. - static double epsilon; // Floating-point machine epsilon. - //static double resulterrbound; - static double ccwerrboundA; // ccwerrboundB, ccwerrboundC; - static double iccerrboundA; // iccerrboundB, iccerrboundC; + private static double epsilon, splitter, resulterrbound; + private static double ccwerrboundA, ccwerrboundB, ccwerrboundC; + private static double iccerrboundA, iccerrboundB, iccerrboundC; + private static double o3derrboundA, o3derrboundB, o3derrboundC; /// - /// Initialize the variables used for exact arithmetic. + /// Initialize the variables used for exact arithmetic. /// /// /// 'epsilon' is the largest power of two such that 1.0 + epsilon = 1.0 in @@ -36,7 +43,7 @@ namespace TriangleNet /// /// I imagine that a highly optimizing compiler might be too smart for its /// own good, and somehow cause this routine to fail, if it pretends that - /// floating-point arithmetic is too much like real arithmetic. + /// floating-point arithmetic is too much like double arithmetic. /// /// Don't change this routine unless you fully understand it. /// @@ -52,7 +59,7 @@ namespace TriangleNet splitter = 1.0; check = 1.0; // Repeatedly divide 'epsilon' by two until it is too small to add to - // one without causing roundoff. (Also check if the sum is equal to + // one without causing roundoff. (Also check if the sum is equal to // the previous sum, for machines that round up instead of using exact // rounding. Not that these routines will work on such machines.) do @@ -67,14 +74,17 @@ namespace TriangleNet 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; + // 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; + ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon; + ccwerrboundC = (9.0 + 64.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; + iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon; + iccerrboundC = (44.0 + 576.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; } /// @@ -88,16 +98,6 @@ namespace TriangleNet /// 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. - /// - /// Uses exact arithmetic if necessary to ensure a correct answer. The - /// result returned is the determinant of a matrix. 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, this function is usually quite fast, but will run - /// more slowly when the input points are collinear or nearly so. - /// - /// See Robust Predicates paper for details. - /// public static double CounterClockwise(Point pa, Point pb, Point pc) { double detleft, detright, det; @@ -147,43 +147,7 @@ namespace TriangleNet return det; } - return (double)CounterClockwiseDecimal(pa, pb, pc); - } - - private static decimal CounterClockwiseDecimal(Point pa, Point pb, Point pc) - { - Statistic.CounterClockwiseCountDecimal++; - - decimal detleft, detright, det, detsum; - - detleft = ((decimal)pa.x - (decimal)pc.x) * ((decimal)pb.y - (decimal)pc.y); - detright = ((decimal)pa.y - (decimal)pc.y) * ((decimal)pb.x - (decimal)pc.x); - det = detleft - detright; - - if (detleft > 0.0m) - { - if (detright <= 0.0m) - { - return det; - } - else - { - detsum = detleft + detright; - } - } - else if (detleft < 0.0m) - { - if (detright >= 0.0m) - { - return det; - } - else - { - detsum = -detleft - detright; - } - } - - return det; + return CounterClockwiseAdapt(pa, pb, pc, detsum); } /// @@ -198,16 +162,6 @@ namespace TriangleNet /// Return a positive value if the point pd lies inside the circle passing through /// pa, pb, and pc; a negative value if it lies outside; and zero if the four points /// are cocircular. - /// - /// Uses exact arithmetic if necessary to ensure a correct answer. The - /// result returned is the determinant of a matrix. 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, this function is usually quite fast, but will run - /// more slowly when the input points are cocircular or nearly so. - /// - /// See Robust Predicates paper for details. - /// public static double InCircle(Point pa, Point pb, Point pc, Point pd) { double adx, bdx, cdx, ady, bdy, cdy; @@ -255,39 +209,7 @@ namespace TriangleNet return det; } - return (double)InCircleDecimal(pa, pb, pc, pd); - } - - private static decimal InCircleDecimal(Point pa, Point pb, Point pc, Point pd) - { - Statistic.InCircleCountDecimal++; - - decimal adx, bdx, cdx, ady, bdy, cdy; - decimal bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady; - decimal alift, blift, clift; - - adx = (decimal)pa.x - (decimal)pd.x; - bdx = (decimal)pb.x - (decimal)pd.x; - cdx = (decimal)pc.x - (decimal)pd.x; - ady = (decimal)pa.y - (decimal)pd.y; - bdy = (decimal)pb.y - (decimal)pd.y; - cdy = (decimal)pc.y - (decimal)pd.y; - - bdxcdy = bdx * cdy; - cdxbdy = cdx * bdy; - alift = adx * adx + ady * ady; - - cdxady = cdx * ady; - adxcdy = adx * cdy; - blift = bdx * bdx + bdy * bdy; - - adxbdy = adx * bdy; - bdxady = bdx * ady; - clift = cdx * cdx + cdy * cdy; - - return alift * (bdxcdy - cdxbdy) - + blift * (cdxady - adxcdy) - + clift * (adxbdy - bdxady); + return InCircleAdapt(pa, pb, pc, pd, permanent); } /// @@ -484,5 +406,856 @@ namespace TriangleNet return new Point(torg.x + dx, torg.y + dy); } + + #region Exact arithmetics + + /// + /// Sum two expansions, eliminating zero components from the output expansion. + /// + /// + /// + /// + /// + /// + /// + /// + /// Sets h = e + f. See my Robust Predicates 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. + /// + static int FastExpansionSumZeroElim(int elen, double[] e, int flen, double[] f, double[] h) + { + double Q; + double Qnew; + double hh; + double bvirt; + double avirt, bround, around; + int eindex, findex, hindex; + double 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)) + { + Qnew = (double)(enow + Q); bvirt = Qnew - enow; hh = Q - bvirt; + enow = e[++eindex]; + } + else + { + Qnew = (double)(fnow + Q); bvirt = Qnew - fnow; hh = Q - bvirt; + fnow = f[++findex]; + } + Q = Qnew; + if (hh != 0.0) + { + h[hindex++] = hh; + } + while ((eindex < elen) && (findex < flen)) + { + if ((fnow > enow) == (fnow > -enow)) + { + Qnew = (double)(Q + enow); + bvirt = (double)(Qnew - Q); + avirt = Qnew - bvirt; + bround = enow - bvirt; + around = Q - avirt; + hh = around + bround; + + enow = e[++eindex]; + } + else + { + Qnew = (double)(Q + fnow); + bvirt = (double)(Qnew - Q); + avirt = Qnew - bvirt; + bround = fnow - bvirt; + around = Q - avirt; + hh = around + bround; + + fnow = f[++findex]; + } + Q = Qnew; + if (hh != 0.0) + { + h[hindex++] = hh; + } + } + } + while (eindex < elen) + { + Qnew = (double)(Q + enow); + bvirt = (double)(Qnew - Q); + avirt = Qnew - bvirt; + bround = enow - bvirt; + around = Q - avirt; + hh = around + bround; + + enow = e[++eindex]; + Q = Qnew; + if (hh != 0.0) + { + h[hindex++] = hh; + } + } + while (findex < flen) + { + Qnew = (double)(Q + fnow); + bvirt = (double)(Qnew - Q); + avirt = Qnew - bvirt; + bround = fnow - bvirt; + around = Q - avirt; + hh = around + bround; + + fnow = f[++findex]; + Q = Qnew; + if (hh != 0.0) + { + h[hindex++] = hh; + } + } + if ((Q != 0.0) || (hindex == 0)) + { + h[hindex++] = Q; + } + return hindex; + } + + /// + /// Multiply an expansion by a scalar, eliminating zero components from the output expansion. + /// + /// + /// + /// + /// + /// + /// + /// Sets h = be. See my Robust Predicates 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.) + /// + static int ScaleExpansionZeroElim(int elen, double[] e, double b, double[] h) + { + double Q, sum; + double hh; + double product1; + double product0; + int eindex, hindex; + double enow; + double bvirt; + double avirt, bround, around; + double c; + double abig; + double ahi, alo, bhi, blo; + double err1, err2, err3; + + c = (double)(splitter * b); abig = (double)(c - b); bhi = c - abig; blo = b - bhi; + Q = (double)(e[0] * b); c = (double)(splitter * e[0]); abig = (double)(c - e[0]); ahi = c - abig; alo = e[0] - ahi; err1 = Q - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); hh = (alo * blo) - err3; + hindex = 0; + if (hh != 0) + { + h[hindex++] = hh; + } + for (eindex = 1; eindex < elen; eindex++) + { + enow = e[eindex]; + product1 = (double)(enow * b); c = (double)(splitter * enow); abig = (double)(c - enow); ahi = c - abig; alo = enow - ahi; err1 = product1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); product0 = (alo * blo) - err3; + sum = (double)(Q + product0); bvirt = (double)(sum - Q); avirt = sum - bvirt; bround = product0 - bvirt; around = Q - avirt; hh = around + bround; + if (hh != 0) + { + h[hindex++] = hh; + } + Q = (double)(product1 + sum); bvirt = Q - product1; hh = sum - bvirt; + if (hh != 0) + { + h[hindex++] = hh; + } + } + if ((Q != 0.0) || (hindex == 0)) + { + h[hindex++] = Q; + } + return hindex; + } + + /// + /// Produce a one-word estimate of an expansion's value. + /// + /// + /// + /// + static double Estimate(int elen, double[] e) + { + double Q; + int eindex; + + Q = e[0]; + for (eindex = 1; eindex < elen; eindex++) + { + Q += e[eindex]; + } + return Q; + } + + /// + /// 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. + /// + /// + /// + /// + /// + /// + /// + /// Uses exact arithmetic if necessary to ensure a correct answer. The result returned + /// is the determinant of a matrix. 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, this function is usually quite fast, + /// but will run more slowly when the input points are collinear or nearly so. + /// + static double CounterClockwiseAdapt(Point pa, Point pb, Point pc, double detsum) + { + double acx, acy, bcx, bcy; + double acxtail, acytail, bcxtail, bcytail; + double detleft, detright; + double detlefttail, detrighttail; + double det, errbound; + double[] B = new double[4], C1 = new double[8], C2 = new double[12], D = new double[16]; + double B3; + int C1length, C2length, Dlength; + double[] u = new double[4]; + double u3; + double s1, t1; + double s0, t0; + + double bvirt; + double avirt, bround, around; + double c; + double abig; + double ahi, alo, bhi, blo; + double err1, err2, err3; + double _i, _j; + double _0; + + acx = (double)(pa.X - pc.X); + bcx = (double)(pb.X - pc.X); + acy = (double)(pa.Y - pc.Y); + bcy = (double)(pb.Y - pc.Y); + + detleft = (double)(acx * bcy); c = (double)(splitter * acx); abig = (double)(c - acx); ahi = c - abig; alo = acx - ahi; c = (double)(splitter * bcy); abig = (double)(c - bcy); bhi = c - abig; blo = bcy - bhi; err1 = detleft - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); detlefttail = (alo * blo) - err3; + detright = (double)(acy * bcx); c = (double)(splitter * acy); abig = (double)(c - acy); ahi = c - abig; alo = acy - ahi; c = (double)(splitter * bcx); abig = (double)(c - bcx); bhi = c - abig; blo = bcx - bhi; err1 = detright - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); detrighttail = (alo * blo) - err3; + + _i = (double)(detlefttail - detrighttail); bvirt = (double)(detlefttail - _i); avirt = _i + bvirt; bround = bvirt - detrighttail; around = detlefttail - avirt; B[0] = around + bround; _j = (double)(detleft + _i); bvirt = (double)(_j - detleft); avirt = _j - bvirt; bround = _i - bvirt; around = detleft - avirt; _0 = around + bround; _i = (double)(_0 - detright); bvirt = (double)(_0 - _i); avirt = _i + bvirt; bround = bvirt - detright; around = _0 - avirt; B[1] = around + bround; B3 = (double)(_j + _i); bvirt = (double)(B3 - _j); avirt = B3 - bvirt; bround = _i - bvirt; around = _j - avirt; B[2] = around + bround; + + B[3] = B3; + + det = Estimate(4, B); + errbound = ccwerrboundB * detsum; + if ((det >= errbound) || (-det >= errbound)) + { + return det; + } + + bvirt = (double)(pa.X - acx); avirt = acx + bvirt; bround = bvirt - pc.X; around = pa.X - avirt; acxtail = around + bround; + bvirt = (double)(pb.X - bcx); avirt = bcx + bvirt; bround = bvirt - pc.X; around = pb.X - avirt; bcxtail = around + bround; + bvirt = (double)(pa.Y - acy); avirt = acy + bvirt; bround = bvirt - pc.Y; around = pa.Y - avirt; acytail = around + bround; + bvirt = (double)(pb.Y - bcy); avirt = bcy + bvirt; bround = bvirt - pc.Y; around = pb.Y - avirt; bcytail = around + bround; + + if ((acxtail == 0.0) && (acytail == 0.0) + && (bcxtail == 0.0) && (bcytail == 0.0)) + { + return det; + } + + errbound = ccwerrboundC * detsum + resulterrbound * ((det) >= 0.0 ? (det) : -(det)); + det += (acx * bcytail + bcy * acxtail) + - (acy * bcxtail + bcx * acytail); + if ((det >= errbound) || (-det >= errbound)) + { + return det; + } + + s1 = (double)(acxtail * bcy); c = (double)(splitter * acxtail); abig = (double)(c - acxtail); ahi = c - abig; alo = acxtail - ahi; c = (double)(splitter * bcy); abig = (double)(c - bcy); bhi = c - abig; blo = bcy - bhi; err1 = s1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); s0 = (alo * blo) - err3; + t1 = (double)(acytail * bcx); c = (double)(splitter * acytail); abig = (double)(c - acytail); ahi = c - abig; alo = acytail - ahi; c = (double)(splitter * bcx); abig = (double)(c - bcx); bhi = c - abig; blo = bcx - bhi; err1 = t1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); t0 = (alo * blo) - err3; + _i = (double)(s0 - t0); bvirt = (double)(s0 - _i); avirt = _i + bvirt; bround = bvirt - t0; around = s0 - avirt; u[0] = around + bround; _j = (double)(s1 + _i); bvirt = (double)(_j - s1); avirt = _j - bvirt; bround = _i - bvirt; around = s1 - avirt; _0 = around + bround; _i = (double)(_0 - t1); bvirt = (double)(_0 - _i); avirt = _i + bvirt; bround = bvirt - t1; around = _0 - avirt; u[1] = around + bround; u3 = (double)(_j + _i); bvirt = (double)(u3 - _j); avirt = u3 - bvirt; bround = _i - bvirt; around = _j - avirt; u[2] = around + bround; + u[3] = u3; + C1length = FastExpansionSumZeroElim(4, B, 4, u, C1); + + s1 = (double)(acx * bcytail); c = (double)(splitter * acx); abig = (double)(c - acx); ahi = c - abig; alo = acx - ahi; c = (double)(splitter * bcytail); abig = (double)(c - bcytail); bhi = c - abig; blo = bcytail - bhi; err1 = s1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); s0 = (alo * blo) - err3; + t1 = (double)(acy * bcxtail); c = (double)(splitter * acy); abig = (double)(c - acy); ahi = c - abig; alo = acy - ahi; c = (double)(splitter * bcxtail); abig = (double)(c - bcxtail); bhi = c - abig; blo = bcxtail - bhi; err1 = t1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); t0 = (alo * blo) - err3; + _i = (double)(s0 - t0); bvirt = (double)(s0 - _i); avirt = _i + bvirt; bround = bvirt - t0; around = s0 - avirt; u[0] = around + bround; _j = (double)(s1 + _i); bvirt = (double)(_j - s1); avirt = _j - bvirt; bround = _i - bvirt; around = s1 - avirt; _0 = around + bround; _i = (double)(_0 - t1); bvirt = (double)(_0 - _i); avirt = _i + bvirt; bround = bvirt - t1; around = _0 - avirt; u[1] = around + bround; u3 = (double)(_j + _i); bvirt = (double)(u3 - _j); avirt = u3 - bvirt; bround = _i - bvirt; around = _j - avirt; u[2] = around + bround; + u[3] = u3; + C2length = FastExpansionSumZeroElim(C1length, C1, 4, u, C2); + + s1 = (double)(acxtail * bcytail); c = (double)(splitter * acxtail); abig = (double)(c - acxtail); ahi = c - abig; alo = acxtail - ahi; c = (double)(splitter * bcytail); abig = (double)(c - bcytail); bhi = c - abig; blo = bcytail - bhi; err1 = s1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); s0 = (alo * blo) - err3; + t1 = (double)(acytail * bcxtail); c = (double)(splitter * acytail); abig = (double)(c - acytail); ahi = c - abig; alo = acytail - ahi; c = (double)(splitter * bcxtail); abig = (double)(c - bcxtail); bhi = c - abig; blo = bcxtail - bhi; err1 = t1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); t0 = (alo * blo) - err3; + _i = (double)(s0 - t0); bvirt = (double)(s0 - _i); avirt = _i + bvirt; bround = bvirt - t0; around = s0 - avirt; u[0] = around + bround; _j = (double)(s1 + _i); bvirt = (double)(_j - s1); avirt = _j - bvirt; bround = _i - bvirt; around = s1 - avirt; _0 = around + bround; _i = (double)(_0 - t1); bvirt = (double)(_0 - _i); avirt = _i + bvirt; bround = bvirt - t1; around = _0 - avirt; u[1] = around + bround; u3 = (double)(_j + _i); bvirt = (double)(u3 - _j); avirt = u3 - bvirt; bround = _i - bvirt; around = _j - avirt; u[2] = around + bround; + u[3] = u3; + Dlength = FastExpansionSumZeroElim(C2length, C2, 4, u, D); + + return (D[Dlength - 1]); + } + + /// + /// Return a positive value if the point pd lies inside the circle passing through + /// pa, pb, and pc; a negative value if it lies outside; and zero if the four points + /// are cocircular. The points pa, pb, and pc must be in counterclockwise order, or + /// the sign of the result will be reversed. + /// + /// + /// + /// + /// + /// + /// + /// + /// Uses exact arithmetic if necessary to ensure a correct answer. The result returned + /// is the determinant of a matrix. 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, this function is usually quite fast, + /// but will run more slowly when the input points are cocircular or nearly so. + /// + static double InCircleAdapt(Point pa, Point pb, Point pc, Point pd, double permanent) + { + double adx, bdx, cdx, ady, bdy, cdy; + double det, errbound; + + double bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1; + double bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0; + double[] bc = new double[4], ca = new double[4], ab = new double[4]; + double bc3, ca3, ab3; + double[] axbc = new double[8], axxbc = new double[16], aybc = new double[8], ayybc = new double[16], adet = new double[32]; + int axbclen, axxbclen, aybclen, ayybclen, alen; + double[] bxca = new double[8], bxxca = new double[16], byca = new double[8], byyca = new double[16], bdet = new double[32]; + int bxcalen, bxxcalen, bycalen, byycalen, blen; + double[] cxab = new double[8], cxxab = new double[16], cyab = new double[8], cyyab = new double[16], cdet = new double[32]; + int cxablen, cxxablen, cyablen, cyyablen, clen; + double[] abdet = new double[64]; + int ablen; + double[] fin1 = new double[1152], fin2 = new double[1152]; + double[] finnow, finother, finswap; + int finlength; + + double adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail; + double adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1; + double adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0; + double[] aa = new double[4], bb = new double[4], cc = new double[4]; + double aa3, bb3, cc3; + double ti1, tj1; + double ti0, tj0; + // Edited to work around rare index out of range exceptions (changed array length from 4 to 5) + // Reason: see loops in FastExpansionSumZeroElim + // while (eindex < elen) { [...] enow = e[++eindex]; [...] } + // Definitely not safe, but Triangle C code doesn't crash!? + double[] u = new double[5], v = new double[5]; + double u3, v3; + double[] temp8 = new double[8], temp16a = new double[16], temp16b = new double[16], temp16c = new double[16]; + double[] temp32a = new double[32], temp32b = new double[32], temp48 = new double[48], temp64 = new double[64]; + int temp8len, temp16alen, temp16blen, temp16clen; + int temp32alen, temp32blen, temp48len, temp64len; + double[] axtbb = new double[8], axtcc = new double[8], aytbb = new double[8], aytcc = new double[8]; + int axtbblen, axtcclen, aytbblen, aytcclen; + double[] bxtaa = new double[8], bxtcc = new double[8], bytaa = new double[8], bytcc = new double[8]; + int bxtaalen, bxtcclen, bytaalen, bytcclen; + double[] cxtaa = new double[8], cxtbb = new double[8], cytaa = new double[8], cytbb = new double[8]; + int cxtaalen, cxtbblen, cytaalen, cytbblen; + double[] axtbc = new double[8], aytbc = new double[8], bxtca = new double[8], bytca = new double[8], cxtab = new double[8], cytab = new double[8]; + int axtbclen = 0, aytbclen = 0, bxtcalen = 0, bytcalen = 0, cxtablen = 0, cytablen = 0; + double[] axtbct = new double[16], aytbct = new double[16], bxtcat = new double[16], bytcat = new double[16], cxtabt = new double[16], cytabt = new double[16]; + int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen; + double[] axtbctt = new double[8], aytbctt = new double[8], bxtcatt = new double[8]; + double[] bytcatt = new double[8], cxtabtt = new double[8], cytabtt = new double[8]; + int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen; + double[] abt = new double[8], bct = new double[8], cat = new double[8]; + int abtlen, bctlen, catlen; + double[] abtt = new double[4], bctt = new double[4], catt = new double[4]; + int abttlen, bcttlen, cattlen; + double abtt3, bctt3, catt3; + double negate; + + double bvirt; + double avirt, bround, around; + double c; + double abig; + double ahi, alo, bhi, blo; + double err1, err2, err3; + double _i, _j; + double _0; + + adx = (double)(pa.X - pd.X); + bdx = (double)(pb.X - pd.X); + cdx = (double)(pc.X - pd.X); + ady = (double)(pa.Y - pd.Y); + bdy = (double)(pb.Y - pd.Y); + cdy = (double)(pc.Y - pd.Y); + + adx = (double)(pa.X - pd.X); + bdx = (double)(pb.X - pd.X); + cdx = (double)(pc.X - pd.X); + ady = (double)(pa.Y - pd.Y); + bdy = (double)(pb.Y - pd.Y); + cdy = (double)(pc.Y - pd.Y); + + bdxcdy1 = (double)(bdx * cdy); c = (double)(splitter * bdx); abig = (double)(c - bdx); ahi = c - abig; alo = bdx - ahi; c = (double)(splitter * cdy); abig = (double)(c - cdy); bhi = c - abig; blo = cdy - bhi; err1 = bdxcdy1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); bdxcdy0 = (alo * blo) - err3; + cdxbdy1 = (double)(cdx * bdy); c = (double)(splitter * cdx); abig = (double)(c - cdx); ahi = c - abig; alo = cdx - ahi; c = (double)(splitter * bdy); abig = (double)(c - bdy); bhi = c - abig; blo = bdy - bhi; err1 = cdxbdy1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); cdxbdy0 = (alo * blo) - err3; + _i = (double)(bdxcdy0 - cdxbdy0); bvirt = (double)(bdxcdy0 - _i); avirt = _i + bvirt; bround = bvirt - cdxbdy0; around = bdxcdy0 - avirt; bc[0] = around + bround; _j = (double)(bdxcdy1 + _i); bvirt = (double)(_j - bdxcdy1); avirt = _j - bvirt; bround = _i - bvirt; around = bdxcdy1 - avirt; _0 = around + bround; _i = (double)(_0 - cdxbdy1); bvirt = (double)(_0 - _i); avirt = _i + bvirt; bround = bvirt - cdxbdy1; around = _0 - avirt; bc[1] = around + bround; bc3 = (double)(_j + _i); bvirt = (double)(bc3 - _j); avirt = bc3 - bvirt; bround = _i - bvirt; around = _j - avirt; bc[2] = around + bround; + bc[3] = bc3; + axbclen = ScaleExpansionZeroElim(4, bc, adx, axbc); + axxbclen = ScaleExpansionZeroElim(axbclen, axbc, adx, axxbc); + aybclen = ScaleExpansionZeroElim(4, bc, ady, aybc); + ayybclen = ScaleExpansionZeroElim(aybclen, aybc, ady, ayybc); + alen = FastExpansionSumZeroElim(axxbclen, axxbc, ayybclen, ayybc, adet); + + cdxady1 = (double)(cdx * ady); c = (double)(splitter * cdx); abig = (double)(c - cdx); ahi = c - abig; alo = cdx - ahi; c = (double)(splitter * ady); abig = (double)(c - ady); bhi = c - abig; blo = ady - bhi; err1 = cdxady1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); cdxady0 = (alo * blo) - err3; + adxcdy1 = (double)(adx * cdy); c = (double)(splitter * adx); abig = (double)(c - adx); ahi = c - abig; alo = adx - ahi; c = (double)(splitter * cdy); abig = (double)(c - cdy); bhi = c - abig; blo = cdy - bhi; err1 = adxcdy1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); adxcdy0 = (alo * blo) - err3; + _i = (double)(cdxady0 - adxcdy0); bvirt = (double)(cdxady0 - _i); avirt = _i + bvirt; bround = bvirt - adxcdy0; around = cdxady0 - avirt; ca[0] = around + bround; _j = (double)(cdxady1 + _i); bvirt = (double)(_j - cdxady1); avirt = _j - bvirt; bround = _i - bvirt; around = cdxady1 - avirt; _0 = around + bround; _i = (double)(_0 - adxcdy1); bvirt = (double)(_0 - _i); avirt = _i + bvirt; bround = bvirt - adxcdy1; around = _0 - avirt; ca[1] = around + bround; ca3 = (double)(_j + _i); bvirt = (double)(ca3 - _j); avirt = ca3 - bvirt; bround = _i - bvirt; around = _j - avirt; ca[2] = around + bround; + ca[3] = ca3; + bxcalen = ScaleExpansionZeroElim(4, ca, bdx, bxca); + bxxcalen = ScaleExpansionZeroElim(bxcalen, bxca, bdx, bxxca); + bycalen = ScaleExpansionZeroElim(4, ca, bdy, byca); + byycalen = ScaleExpansionZeroElim(bycalen, byca, bdy, byyca); + blen = FastExpansionSumZeroElim(bxxcalen, bxxca, byycalen, byyca, bdet); + + adxbdy1 = (double)(adx * bdy); c = (double)(splitter * adx); abig = (double)(c - adx); ahi = c - abig; alo = adx - ahi; c = (double)(splitter * bdy); abig = (double)(c - bdy); bhi = c - abig; blo = bdy - bhi; err1 = adxbdy1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); adxbdy0 = (alo * blo) - err3; + bdxady1 = (double)(bdx * ady); c = (double)(splitter * bdx); abig = (double)(c - bdx); ahi = c - abig; alo = bdx - ahi; c = (double)(splitter * ady); abig = (double)(c - ady); bhi = c - abig; blo = ady - bhi; err1 = bdxady1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); bdxady0 = (alo * blo) - err3; + _i = (double)(adxbdy0 - bdxady0); bvirt = (double)(adxbdy0 - _i); avirt = _i + bvirt; bround = bvirt - bdxady0; around = adxbdy0 - avirt; ab[0] = around + bround; _j = (double)(adxbdy1 + _i); bvirt = (double)(_j - adxbdy1); avirt = _j - bvirt; bround = _i - bvirt; around = adxbdy1 - avirt; _0 = around + bround; _i = (double)(_0 - bdxady1); bvirt = (double)(_0 - _i); avirt = _i + bvirt; bround = bvirt - bdxady1; around = _0 - avirt; ab[1] = around + bround; ab3 = (double)(_j + _i); bvirt = (double)(ab3 - _j); avirt = ab3 - bvirt; bround = _i - bvirt; around = _j - avirt; ab[2] = around + bround; + ab[3] = ab3; + cxablen = ScaleExpansionZeroElim(4, ab, cdx, cxab); + cxxablen = ScaleExpansionZeroElim(cxablen, cxab, cdx, cxxab); + cyablen = ScaleExpansionZeroElim(4, ab, cdy, cyab); + cyyablen = ScaleExpansionZeroElim(cyablen, cyab, cdy, cyyab); + clen = FastExpansionSumZeroElim(cxxablen, cxxab, cyyablen, cyyab, cdet); + + ablen = FastExpansionSumZeroElim(alen, adet, blen, bdet, abdet); + finlength = FastExpansionSumZeroElim(ablen, abdet, clen, cdet, fin1); + + det = Estimate(finlength, fin1); + errbound = iccerrboundB * permanent; + if ((det >= errbound) || (-det >= errbound)) + { + return det; + } + + bvirt = (double)(pa.X - adx); avirt = adx + bvirt; bround = bvirt - pd.X; around = pa.X - avirt; adxtail = around + bround; + bvirt = (double)(pa.Y - ady); avirt = ady + bvirt; bround = bvirt - pd.Y; around = pa.Y - avirt; adytail = around + bround; + bvirt = (double)(pb.X - bdx); avirt = bdx + bvirt; bround = bvirt - pd.X; around = pb.X - avirt; bdxtail = around + bround; + bvirt = (double)(pb.Y - bdy); avirt = bdy + bvirt; bround = bvirt - pd.Y; around = pb.Y - avirt; bdytail = around + bround; + bvirt = (double)(pc.X - cdx); avirt = cdx + bvirt; bround = bvirt - pd.X; around = pc.X - avirt; cdxtail = around + bround; + bvirt = (double)(pc.Y - cdy); avirt = cdy + bvirt; bround = bvirt - pd.Y; around = pc.Y - avirt; cdytail = around + bround; + if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) + && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)) + { + return det; + } + + errbound = iccerrboundC * permanent + resulterrbound * ((det) >= 0.0 ? (det) : -(det)); + det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail) - (bdy * cdxtail + cdx * bdytail)) + + 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx)) + + ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail) - (cdy * adxtail + adx * cdytail)) + + 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx)) + + ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail) - (ady * bdxtail + bdx * adytail)) + + 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx)); + if ((det >= errbound) || (-det >= errbound)) + { + return det; + } + + finnow = fin1; + finother = fin2; + + if ((bdxtail != 0.0) || (bdytail != 0.0) + || (cdxtail != 0.0) || (cdytail != 0.0)) + { + adxadx1 = (double)(adx * adx); c = (double)(splitter * adx); abig = (double)(c - adx); ahi = c - abig; alo = adx - ahi; err1 = adxadx1 - (ahi * ahi); err3 = err1 - ((ahi + ahi) * alo); adxadx0 = (alo * alo) - err3; + adyady1 = (double)(ady * ady); c = (double)(splitter * ady); abig = (double)(c - ady); ahi = c - abig; alo = ady - ahi; err1 = adyady1 - (ahi * ahi); err3 = err1 - ((ahi + ahi) * alo); adyady0 = (alo * alo) - err3; + _i = (double)(adxadx0 + adyady0); bvirt = (double)(_i - adxadx0); avirt = _i - bvirt; bround = adyady0 - bvirt; around = adxadx0 - avirt; aa[0] = around + bround; _j = (double)(adxadx1 + _i); bvirt = (double)(_j - adxadx1); avirt = _j - bvirt; bround = _i - bvirt; around = adxadx1 - avirt; _0 = around + bround; _i = (double)(_0 + adyady1); bvirt = (double)(_i - _0); avirt = _i - bvirt; bround = adyady1 - bvirt; around = _0 - avirt; aa[1] = around + bround; aa3 = (double)(_j + _i); bvirt = (double)(aa3 - _j); avirt = aa3 - bvirt; bround = _i - bvirt; around = _j - avirt; aa[2] = around + bround; + aa[3] = aa3; + } + if ((cdxtail != 0.0) || (cdytail != 0.0) + || (adxtail != 0.0) || (adytail != 0.0)) + { + bdxbdx1 = (double)(bdx * bdx); c = (double)(splitter * bdx); abig = (double)(c - bdx); ahi = c - abig; alo = bdx - ahi; err1 = bdxbdx1 - (ahi * ahi); err3 = err1 - ((ahi + ahi) * alo); bdxbdx0 = (alo * alo) - err3; + bdybdy1 = (double)(bdy * bdy); c = (double)(splitter * bdy); abig = (double)(c - bdy); ahi = c - abig; alo = bdy - ahi; err1 = bdybdy1 - (ahi * ahi); err3 = err1 - ((ahi + ahi) * alo); bdybdy0 = (alo * alo) - err3; + _i = (double)(bdxbdx0 + bdybdy0); bvirt = (double)(_i - bdxbdx0); avirt = _i - bvirt; bround = bdybdy0 - bvirt; around = bdxbdx0 - avirt; bb[0] = around + bround; _j = (double)(bdxbdx1 + _i); bvirt = (double)(_j - bdxbdx1); avirt = _j - bvirt; bround = _i - bvirt; around = bdxbdx1 - avirt; _0 = around + bround; _i = (double)(_0 + bdybdy1); bvirt = (double)(_i - _0); avirt = _i - bvirt; bround = bdybdy1 - bvirt; around = _0 - avirt; bb[1] = around + bround; bb3 = (double)(_j + _i); bvirt = (double)(bb3 - _j); avirt = bb3 - bvirt; bround = _i - bvirt; around = _j - avirt; bb[2] = around + bround; + bb[3] = bb3; + } + if ((adxtail != 0.0) || (adytail != 0.0) + || (bdxtail != 0.0) || (bdytail != 0.0)) + { + cdxcdx1 = (double)(cdx * cdx); c = (double)(splitter * cdx); abig = (double)(c - cdx); ahi = c - abig; alo = cdx - ahi; err1 = cdxcdx1 - (ahi * ahi); err3 = err1 - ((ahi + ahi) * alo); cdxcdx0 = (alo * alo) - err3; + cdycdy1 = (double)(cdy * cdy); c = (double)(splitter * cdy); abig = (double)(c - cdy); ahi = c - abig; alo = cdy - ahi; err1 = cdycdy1 - (ahi * ahi); err3 = err1 - ((ahi + ahi) * alo); cdycdy0 = (alo * alo) - err3; + _i = (double)(cdxcdx0 + cdycdy0); bvirt = (double)(_i - cdxcdx0); avirt = _i - bvirt; bround = cdycdy0 - bvirt; around = cdxcdx0 - avirt; cc[0] = around + bround; _j = (double)(cdxcdx1 + _i); bvirt = (double)(_j - cdxcdx1); avirt = _j - bvirt; bround = _i - bvirt; around = cdxcdx1 - avirt; _0 = around + bround; _i = (double)(_0 + cdycdy1); bvirt = (double)(_i - _0); avirt = _i - bvirt; bround = cdycdy1 - bvirt; around = _0 - avirt; cc[1] = around + bround; cc3 = (double)(_j + _i); bvirt = (double)(cc3 - _j); avirt = cc3 - bvirt; bround = _i - bvirt; around = _j - avirt; cc[2] = around + bround; + cc[3] = cc3; + } + + if (adxtail != 0.0) + { + axtbclen = ScaleExpansionZeroElim(4, bc, adxtail, axtbc); + temp16alen = ScaleExpansionZeroElim(axtbclen, axtbc, 2.0 * adx, temp16a); + + axtcclen = ScaleExpansionZeroElim(4, cc, adxtail, axtcc); + temp16blen = ScaleExpansionZeroElim(axtcclen, axtcc, bdy, temp16b); + + axtbblen = ScaleExpansionZeroElim(4, bb, adxtail, axtbb); + temp16clen = ScaleExpansionZeroElim(axtbblen, axtbb, -cdy, temp16c); + + temp32alen = FastExpansionSumZeroElim(temp16alen, temp16a, temp16blen, temp16b, temp32a); + temp48len = FastExpansionSumZeroElim(temp16clen, temp16c, temp32alen, temp32a, temp48); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (adytail != 0.0) + { + aytbclen = ScaleExpansionZeroElim(4, bc, adytail, aytbc); + temp16alen = ScaleExpansionZeroElim(aytbclen, aytbc, 2.0 * ady, temp16a); + + aytbblen = ScaleExpansionZeroElim(4, bb, adytail, aytbb); + temp16blen = ScaleExpansionZeroElim(aytbblen, aytbb, cdx, temp16b); + + aytcclen = ScaleExpansionZeroElim(4, cc, adytail, aytcc); + temp16clen = ScaleExpansionZeroElim(aytcclen, aytcc, -bdx, temp16c); + + temp32alen = FastExpansionSumZeroElim(temp16alen, temp16a, temp16blen, temp16b, temp32a); + temp48len = FastExpansionSumZeroElim(temp16clen, temp16c, temp32alen, temp32a, temp48); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (bdxtail != 0.0) + { + bxtcalen = ScaleExpansionZeroElim(4, ca, bdxtail, bxtca); + temp16alen = ScaleExpansionZeroElim(bxtcalen, bxtca, 2.0 * bdx, temp16a); + + bxtaalen = ScaleExpansionZeroElim(4, aa, bdxtail, bxtaa); + temp16blen = ScaleExpansionZeroElim(bxtaalen, bxtaa, cdy, temp16b); + + bxtcclen = ScaleExpansionZeroElim(4, cc, bdxtail, bxtcc); + temp16clen = ScaleExpansionZeroElim(bxtcclen, bxtcc, -ady, temp16c); + + temp32alen = FastExpansionSumZeroElim(temp16alen, temp16a, temp16blen, temp16b, temp32a); + temp48len = FastExpansionSumZeroElim(temp16clen, temp16c, temp32alen, temp32a, temp48); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (bdytail != 0.0) + { + bytcalen = ScaleExpansionZeroElim(4, ca, bdytail, bytca); + temp16alen = ScaleExpansionZeroElim(bytcalen, bytca, 2.0 * bdy, temp16a); + + bytcclen = ScaleExpansionZeroElim(4, cc, bdytail, bytcc); + temp16blen = ScaleExpansionZeroElim(bytcclen, bytcc, adx, temp16b); + + bytaalen = ScaleExpansionZeroElim(4, aa, bdytail, bytaa); + temp16clen = ScaleExpansionZeroElim(bytaalen, bytaa, -cdx, temp16c); + + temp32alen = FastExpansionSumZeroElim(temp16alen, temp16a, temp16blen, temp16b, temp32a); + temp48len = FastExpansionSumZeroElim(temp16clen, temp16c, temp32alen, temp32a, temp48); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (cdxtail != 0.0) + { + cxtablen = ScaleExpansionZeroElim(4, ab, cdxtail, cxtab); + temp16alen = ScaleExpansionZeroElim(cxtablen, cxtab, 2.0 * cdx, temp16a); + + cxtbblen = ScaleExpansionZeroElim(4, bb, cdxtail, cxtbb); + temp16blen = ScaleExpansionZeroElim(cxtbblen, cxtbb, ady, temp16b); + + cxtaalen = ScaleExpansionZeroElim(4, aa, cdxtail, cxtaa); + temp16clen = ScaleExpansionZeroElim(cxtaalen, cxtaa, -bdy, temp16c); + + temp32alen = FastExpansionSumZeroElim(temp16alen, temp16a, temp16blen, temp16b, temp32a); + temp48len = FastExpansionSumZeroElim(temp16clen, temp16c, temp32alen, temp32a, temp48); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (cdytail != 0.0) + { + cytablen = ScaleExpansionZeroElim(4, ab, cdytail, cytab); + temp16alen = ScaleExpansionZeroElim(cytablen, cytab, 2.0 * cdy, temp16a); + + cytaalen = ScaleExpansionZeroElim(4, aa, cdytail, cytaa); + temp16blen = ScaleExpansionZeroElim(cytaalen, cytaa, bdx, temp16b); + + cytbblen = ScaleExpansionZeroElim(4, bb, cdytail, cytbb); + temp16clen = ScaleExpansionZeroElim(cytbblen, cytbb, -adx, temp16c); + + temp32alen = FastExpansionSumZeroElim(temp16alen, temp16a, temp16blen, temp16b, temp32a); + temp48len = FastExpansionSumZeroElim(temp16clen, temp16c, temp32alen, temp32a, temp48); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + + if ((adxtail != 0.0) || (adytail != 0.0)) + { + if ((bdxtail != 0.0) || (bdytail != 0.0) + || (cdxtail != 0.0) || (cdytail != 0.0)) + { + ti1 = (double)(bdxtail * cdy); c = (double)(splitter * bdxtail); abig = (double)(c - bdxtail); ahi = c - abig; alo = bdxtail - ahi; c = (double)(splitter * cdy); abig = (double)(c - cdy); bhi = c - abig; blo = cdy - bhi; err1 = ti1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); ti0 = (alo * blo) - err3; + tj1 = (double)(bdx * cdytail); c = (double)(splitter * bdx); abig = (double)(c - bdx); ahi = c - abig; alo = bdx - ahi; c = (double)(splitter * cdytail); abig = (double)(c - cdytail); bhi = c - abig; blo = cdytail - bhi; err1 = tj1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); tj0 = (alo * blo) - err3; + _i = (double)(ti0 + tj0); bvirt = (double)(_i - ti0); avirt = _i - bvirt; bround = tj0 - bvirt; around = ti0 - avirt; u[0] = around + bround; _j = (double)(ti1 + _i); bvirt = (double)(_j - ti1); avirt = _j - bvirt; bround = _i - bvirt; around = ti1 - avirt; _0 = around + bround; _i = (double)(_0 + tj1); bvirt = (double)(_i - _0); avirt = _i - bvirt; bround = tj1 - bvirt; around = _0 - avirt; u[1] = around + bround; u3 = (double)(_j + _i); bvirt = (double)(u3 - _j); avirt = u3 - bvirt; bround = _i - bvirt; around = _j - avirt; u[2] = around + bround; + u[3] = u3; + negate = -bdy; + ti1 = (double)(cdxtail * negate); c = (double)(splitter * cdxtail); abig = (double)(c - cdxtail); ahi = c - abig; alo = cdxtail - ahi; c = (double)(splitter * negate); abig = (double)(c - negate); bhi = c - abig; blo = negate - bhi; err1 = ti1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); ti0 = (alo * blo) - err3; + negate = -bdytail; + tj1 = (double)(cdx * negate); c = (double)(splitter * cdx); abig = (double)(c - cdx); ahi = c - abig; alo = cdx - ahi; c = (double)(splitter * negate); abig = (double)(c - negate); bhi = c - abig; blo = negate - bhi; err1 = tj1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); tj0 = (alo * blo) - err3; + _i = (double)(ti0 + tj0); bvirt = (double)(_i - ti0); avirt = _i - bvirt; bround = tj0 - bvirt; around = ti0 - avirt; v[0] = around + bround; _j = (double)(ti1 + _i); bvirt = (double)(_j - ti1); avirt = _j - bvirt; bround = _i - bvirt; around = ti1 - avirt; _0 = around + bround; _i = (double)(_0 + tj1); bvirt = (double)(_i - _0); avirt = _i - bvirt; bround = tj1 - bvirt; around = _0 - avirt; v[1] = around + bround; v3 = (double)(_j + _i); bvirt = (double)(v3 - _j); avirt = v3 - bvirt; bround = _i - bvirt; around = _j - avirt; v[2] = around + bround; + v[3] = v3; + bctlen = FastExpansionSumZeroElim(4, u, 4, v, bct); + + ti1 = (double)(bdxtail * cdytail); c = (double)(splitter * bdxtail); abig = (double)(c - bdxtail); ahi = c - abig; alo = bdxtail - ahi; c = (double)(splitter * cdytail); abig = (double)(c - cdytail); bhi = c - abig; blo = cdytail - bhi; err1 = ti1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); ti0 = (alo * blo) - err3; + tj1 = (double)(cdxtail * bdytail); c = (double)(splitter * cdxtail); abig = (double)(c - cdxtail); ahi = c - abig; alo = cdxtail - ahi; c = (double)(splitter * bdytail); abig = (double)(c - bdytail); bhi = c - abig; blo = bdytail - bhi; err1 = tj1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); tj0 = (alo * blo) - err3; + _i = (double)(ti0 - tj0); bvirt = (double)(ti0 - _i); avirt = _i + bvirt; bround = bvirt - tj0; around = ti0 - avirt; bctt[0] = around + bround; _j = (double)(ti1 + _i); bvirt = (double)(_j - ti1); avirt = _j - bvirt; bround = _i - bvirt; around = ti1 - avirt; _0 = around + bround; _i = (double)(_0 - tj1); bvirt = (double)(_0 - _i); avirt = _i + bvirt; bround = bvirt - tj1; around = _0 - avirt; bctt[1] = around + bround; bctt3 = (double)(_j + _i); bvirt = (double)(bctt3 - _j); avirt = bctt3 - bvirt; bround = _i - bvirt; around = _j - avirt; bctt[2] = around + bround; + bctt[3] = bctt3; + bcttlen = 4; + } + else + { + bct[0] = 0.0; + bctlen = 1; + bctt[0] = 0.0; + bcttlen = 1; + } + + if (adxtail != 0.0) + { + temp16alen = ScaleExpansionZeroElim(axtbclen, axtbc, adxtail, temp16a); + axtbctlen = ScaleExpansionZeroElim(bctlen, bct, adxtail, axtbct); + temp32alen = ScaleExpansionZeroElim(axtbctlen, axtbct, 2.0 * adx, temp32a); + temp48len = FastExpansionSumZeroElim(temp16alen, temp16a, temp32alen, temp32a, temp48); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + if (bdytail != 0.0) + { + temp8len = ScaleExpansionZeroElim(4, cc, adxtail, temp8); + temp16alen = ScaleExpansionZeroElim(temp8len, temp8, bdytail, temp16a); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp16alen, temp16a, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (cdytail != 0.0) + { + temp8len = ScaleExpansionZeroElim(4, bb, -adxtail, temp8); + temp16alen = ScaleExpansionZeroElim(temp8len, temp8, cdytail, temp16a); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp16alen, temp16a, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + + temp32alen = ScaleExpansionZeroElim(axtbctlen, axtbct, adxtail, temp32a); + axtbcttlen = ScaleExpansionZeroElim(bcttlen, bctt, adxtail, axtbctt); + temp16alen = ScaleExpansionZeroElim(axtbcttlen, axtbctt, 2.0 * adx, temp16a); + temp16blen = ScaleExpansionZeroElim(axtbcttlen, axtbctt, adxtail, temp16b); + temp32blen = FastExpansionSumZeroElim(temp16alen, temp16a, temp16blen, temp16b, temp32b); + temp64len = FastExpansionSumZeroElim(temp32alen, temp32a, temp32blen, temp32b, temp64); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp64len, temp64, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (adytail != 0.0) + { + temp16alen = ScaleExpansionZeroElim(aytbclen, aytbc, adytail, temp16a); + aytbctlen = ScaleExpansionZeroElim(bctlen, bct, adytail, aytbct); + temp32alen = ScaleExpansionZeroElim(aytbctlen, aytbct, 2.0 * ady, temp32a); + temp48len = FastExpansionSumZeroElim(temp16alen, temp16a, temp32alen, temp32a, temp48); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + + + temp32alen = ScaleExpansionZeroElim(aytbctlen, aytbct, adytail, temp32a); + aytbcttlen = ScaleExpansionZeroElim(bcttlen, bctt, adytail, aytbctt); + temp16alen = ScaleExpansionZeroElim(aytbcttlen, aytbctt, 2.0 * ady, temp16a); + temp16blen = ScaleExpansionZeroElim(aytbcttlen, aytbctt, adytail, temp16b); + temp32blen = FastExpansionSumZeroElim(temp16alen, temp16a, temp16blen, temp16b, temp32b); + temp64len = FastExpansionSumZeroElim(temp32alen, temp32a, temp32blen, temp32b, temp64); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp64len, temp64, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + if ((bdxtail != 0.0) || (bdytail != 0.0)) + { + if ((cdxtail != 0.0) || (cdytail != 0.0) + || (adxtail != 0.0) || (adytail != 0.0)) + { + ti1 = (double)(cdxtail * ady); c = (double)(splitter * cdxtail); abig = (double)(c - cdxtail); ahi = c - abig; alo = cdxtail - ahi; c = (double)(splitter * ady); abig = (double)(c - ady); bhi = c - abig; blo = ady - bhi; err1 = ti1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); ti0 = (alo * blo) - err3; + tj1 = (double)(cdx * adytail); c = (double)(splitter * cdx); abig = (double)(c - cdx); ahi = c - abig; alo = cdx - ahi; c = (double)(splitter * adytail); abig = (double)(c - adytail); bhi = c - abig; blo = adytail - bhi; err1 = tj1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); tj0 = (alo * blo) - err3; + _i = (double)(ti0 + tj0); bvirt = (double)(_i - ti0); avirt = _i - bvirt; bround = tj0 - bvirt; around = ti0 - avirt; u[0] = around + bround; _j = (double)(ti1 + _i); bvirt = (double)(_j - ti1); avirt = _j - bvirt; bround = _i - bvirt; around = ti1 - avirt; _0 = around + bround; _i = (double)(_0 + tj1); bvirt = (double)(_i - _0); avirt = _i - bvirt; bround = tj1 - bvirt; around = _0 - avirt; u[1] = around + bround; u3 = (double)(_j + _i); bvirt = (double)(u3 - _j); avirt = u3 - bvirt; bround = _i - bvirt; around = _j - avirt; u[2] = around + bround; + u[3] = u3; + negate = -cdy; + ti1 = (double)(adxtail * negate); c = (double)(splitter * adxtail); abig = (double)(c - adxtail); ahi = c - abig; alo = adxtail - ahi; c = (double)(splitter * negate); abig = (double)(c - negate); bhi = c - abig; blo = negate - bhi; err1 = ti1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); ti0 = (alo * blo) - err3; + negate = -cdytail; + tj1 = (double)(adx * negate); c = (double)(splitter * adx); abig = (double)(c - adx); ahi = c - abig; alo = adx - ahi; c = (double)(splitter * negate); abig = (double)(c - negate); bhi = c - abig; blo = negate - bhi; err1 = tj1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); tj0 = (alo * blo) - err3; + _i = (double)(ti0 + tj0); bvirt = (double)(_i - ti0); avirt = _i - bvirt; bround = tj0 - bvirt; around = ti0 - avirt; v[0] = around + bround; _j = (double)(ti1 + _i); bvirt = (double)(_j - ti1); avirt = _j - bvirt; bround = _i - bvirt; around = ti1 - avirt; _0 = around + bround; _i = (double)(_0 + tj1); bvirt = (double)(_i - _0); avirt = _i - bvirt; bround = tj1 - bvirt; around = _0 - avirt; v[1] = around + bround; v3 = (double)(_j + _i); bvirt = (double)(v3 - _j); avirt = v3 - bvirt; bround = _i - bvirt; around = _j - avirt; v[2] = around + bround; + v[3] = v3; + catlen = FastExpansionSumZeroElim(4, u, 4, v, cat); + + ti1 = (double)(cdxtail * adytail); c = (double)(splitter * cdxtail); abig = (double)(c - cdxtail); ahi = c - abig; alo = cdxtail - ahi; c = (double)(splitter * adytail); abig = (double)(c - adytail); bhi = c - abig; blo = adytail - bhi; err1 = ti1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); ti0 = (alo * blo) - err3; + tj1 = (double)(adxtail * cdytail); c = (double)(splitter * adxtail); abig = (double)(c - adxtail); ahi = c - abig; alo = adxtail - ahi; c = (double)(splitter * cdytail); abig = (double)(c - cdytail); bhi = c - abig; blo = cdytail - bhi; err1 = tj1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); tj0 = (alo * blo) - err3; + _i = (double)(ti0 - tj0); bvirt = (double)(ti0 - _i); avirt = _i + bvirt; bround = bvirt - tj0; around = ti0 - avirt; catt[0] = around + bround; _j = (double)(ti1 + _i); bvirt = (double)(_j - ti1); avirt = _j - bvirt; bround = _i - bvirt; around = ti1 - avirt; _0 = around + bround; _i = (double)(_0 - tj1); bvirt = (double)(_0 - _i); avirt = _i + bvirt; bround = bvirt - tj1; around = _0 - avirt; catt[1] = around + bround; catt3 = (double)(_j + _i); bvirt = (double)(catt3 - _j); avirt = catt3 - bvirt; bround = _i - bvirt; around = _j - avirt; catt[2] = around + bround; + catt[3] = catt3; + cattlen = 4; + } + else + { + cat[0] = 0.0; + catlen = 1; + catt[0] = 0.0; + cattlen = 1; + } + + if (bdxtail != 0.0) + { + temp16alen = ScaleExpansionZeroElim(bxtcalen, bxtca, bdxtail, temp16a); + bxtcatlen = ScaleExpansionZeroElim(catlen, cat, bdxtail, bxtcat); + temp32alen = ScaleExpansionZeroElim(bxtcatlen, bxtcat, 2.0 * bdx, temp32a); + temp48len = FastExpansionSumZeroElim(temp16alen, temp16a, temp32alen, temp32a, temp48); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + if (cdytail != 0.0) + { + temp8len = ScaleExpansionZeroElim(4, aa, bdxtail, temp8); + temp16alen = ScaleExpansionZeroElim(temp8len, temp8, cdytail, temp16a); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp16alen, temp16a, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (adytail != 0.0) + { + temp8len = ScaleExpansionZeroElim(4, cc, -bdxtail, temp8); + temp16alen = ScaleExpansionZeroElim(temp8len, temp8, adytail, temp16a); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp16alen, temp16a, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + + temp32alen = ScaleExpansionZeroElim(bxtcatlen, bxtcat, bdxtail, temp32a); + bxtcattlen = ScaleExpansionZeroElim(cattlen, catt, bdxtail, bxtcatt); + temp16alen = ScaleExpansionZeroElim(bxtcattlen, bxtcatt, 2.0 * bdx, temp16a); + temp16blen = ScaleExpansionZeroElim(bxtcattlen, bxtcatt, bdxtail, temp16b); + temp32blen = FastExpansionSumZeroElim(temp16alen, temp16a, temp16blen, temp16b, temp32b); + temp64len = FastExpansionSumZeroElim(temp32alen, temp32a, temp32blen, temp32b, temp64); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp64len, temp64, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (bdytail != 0.0) + { + temp16alen = ScaleExpansionZeroElim(bytcalen, bytca, bdytail, temp16a); + bytcatlen = ScaleExpansionZeroElim(catlen, cat, bdytail, bytcat); + temp32alen = ScaleExpansionZeroElim(bytcatlen, bytcat, 2.0 * bdy, temp32a); + temp48len = FastExpansionSumZeroElim(temp16alen, temp16a, temp32alen, temp32a, temp48); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + + temp32alen = ScaleExpansionZeroElim(bytcatlen, bytcat, bdytail, temp32a); + bytcattlen = ScaleExpansionZeroElim(cattlen, catt, bdytail, bytcatt); + temp16alen = ScaleExpansionZeroElim(bytcattlen, bytcatt, 2.0 * bdy, temp16a); + temp16blen = ScaleExpansionZeroElim(bytcattlen, bytcatt, bdytail, temp16b); + temp32blen = FastExpansionSumZeroElim(temp16alen, temp16a, temp16blen, temp16b, temp32b); + temp64len = FastExpansionSumZeroElim(temp32alen, temp32a, temp32blen, temp32b, temp64); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp64len, temp64, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + if ((cdxtail != 0.0) || (cdytail != 0.0)) + { + if ((adxtail != 0.0) || (adytail != 0.0) + || (bdxtail != 0.0) || (bdytail != 0.0)) + { + ti1 = (double)(adxtail * bdy); c = (double)(splitter * adxtail); abig = (double)(c - adxtail); ahi = c - abig; alo = adxtail - ahi; c = (double)(splitter * bdy); abig = (double)(c - bdy); bhi = c - abig; blo = bdy - bhi; err1 = ti1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); ti0 = (alo * blo) - err3; + tj1 = (double)(adx * bdytail); c = (double)(splitter * adx); abig = (double)(c - adx); ahi = c - abig; alo = adx - ahi; c = (double)(splitter * bdytail); abig = (double)(c - bdytail); bhi = c - abig; blo = bdytail - bhi; err1 = tj1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); tj0 = (alo * blo) - err3; + _i = (double)(ti0 + tj0); bvirt = (double)(_i - ti0); avirt = _i - bvirt; bround = tj0 - bvirt; around = ti0 - avirt; u[0] = around + bround; _j = (double)(ti1 + _i); bvirt = (double)(_j - ti1); avirt = _j - bvirt; bround = _i - bvirt; around = ti1 - avirt; _0 = around + bround; _i = (double)(_0 + tj1); bvirt = (double)(_i - _0); avirt = _i - bvirt; bround = tj1 - bvirt; around = _0 - avirt; u[1] = around + bround; u3 = (double)(_j + _i); bvirt = (double)(u3 - _j); avirt = u3 - bvirt; bround = _i - bvirt; around = _j - avirt; u[2] = around + bround; + u[3] = u3; + negate = -ady; + ti1 = (double)(bdxtail * negate); c = (double)(splitter * bdxtail); abig = (double)(c - bdxtail); ahi = c - abig; alo = bdxtail - ahi; c = (double)(splitter * negate); abig = (double)(c - negate); bhi = c - abig; blo = negate - bhi; err1 = ti1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); ti0 = (alo * blo) - err3; + negate = -adytail; + tj1 = (double)(bdx * negate); c = (double)(splitter * bdx); abig = (double)(c - bdx); ahi = c - abig; alo = bdx - ahi; c = (double)(splitter * negate); abig = (double)(c - negate); bhi = c - abig; blo = negate - bhi; err1 = tj1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); tj0 = (alo * blo) - err3; + _i = (double)(ti0 + tj0); bvirt = (double)(_i - ti0); avirt = _i - bvirt; bround = tj0 - bvirt; around = ti0 - avirt; v[0] = around + bround; _j = (double)(ti1 + _i); bvirt = (double)(_j - ti1); avirt = _j - bvirt; bround = _i - bvirt; around = ti1 - avirt; _0 = around + bround; _i = (double)(_0 + tj1); bvirt = (double)(_i - _0); avirt = _i - bvirt; bround = tj1 - bvirt; around = _0 - avirt; v[1] = around + bround; v3 = (double)(_j + _i); bvirt = (double)(v3 - _j); avirt = v3 - bvirt; bround = _i - bvirt; around = _j - avirt; v[2] = around + bround; + v[3] = v3; + abtlen = FastExpansionSumZeroElim(4, u, 4, v, abt); + + ti1 = (double)(adxtail * bdytail); c = (double)(splitter * adxtail); abig = (double)(c - adxtail); ahi = c - abig; alo = adxtail - ahi; c = (double)(splitter * bdytail); abig = (double)(c - bdytail); bhi = c - abig; blo = bdytail - bhi; err1 = ti1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); ti0 = (alo * blo) - err3; + tj1 = (double)(bdxtail * adytail); c = (double)(splitter * bdxtail); abig = (double)(c - bdxtail); ahi = c - abig; alo = bdxtail - ahi; c = (double)(splitter * adytail); abig = (double)(c - adytail); bhi = c - abig; blo = adytail - bhi; err1 = tj1 - (ahi * bhi); err2 = err1 - (alo * bhi); err3 = err2 - (ahi * blo); tj0 = (alo * blo) - err3; + _i = (double)(ti0 - tj0); bvirt = (double)(ti0 - _i); avirt = _i + bvirt; bround = bvirt - tj0; around = ti0 - avirt; abtt[0] = around + bround; _j = (double)(ti1 + _i); bvirt = (double)(_j - ti1); avirt = _j - bvirt; bround = _i - bvirt; around = ti1 - avirt; _0 = around + bround; _i = (double)(_0 - tj1); bvirt = (double)(_0 - _i); avirt = _i + bvirt; bround = bvirt - tj1; around = _0 - avirt; abtt[1] = around + bround; abtt3 = (double)(_j + _i); bvirt = (double)(abtt3 - _j); avirt = abtt3 - bvirt; bround = _i - bvirt; around = _j - avirt; abtt[2] = around + bround; + abtt[3] = abtt3; + abttlen = 4; + } + else + { + abt[0] = 0.0; + abtlen = 1; + abtt[0] = 0.0; + abttlen = 1; + } + + if (cdxtail != 0.0) + { + temp16alen = ScaleExpansionZeroElim(cxtablen, cxtab, cdxtail, temp16a); + cxtabtlen = ScaleExpansionZeroElim(abtlen, abt, cdxtail, cxtabt); + temp32alen = ScaleExpansionZeroElim(cxtabtlen, cxtabt, 2.0 * cdx, temp32a); + temp48len = FastExpansionSumZeroElim(temp16alen, temp16a, temp32alen, temp32a, temp48); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + if (adytail != 0.0) + { + temp8len = ScaleExpansionZeroElim(4, bb, cdxtail, temp8); + temp16alen = ScaleExpansionZeroElim(temp8len, temp8, adytail, temp16a); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp16alen, temp16a, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (bdytail != 0.0) + { + temp8len = ScaleExpansionZeroElim(4, aa, -cdxtail, temp8); + temp16alen = ScaleExpansionZeroElim(temp8len, temp8, bdytail, temp16a); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp16alen, temp16a, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + + temp32alen = ScaleExpansionZeroElim(cxtabtlen, cxtabt, cdxtail, temp32a); + cxtabttlen = ScaleExpansionZeroElim(abttlen, abtt, cdxtail, cxtabtt); + temp16alen = ScaleExpansionZeroElim(cxtabttlen, cxtabtt, 2.0 * cdx, temp16a); + temp16blen = ScaleExpansionZeroElim(cxtabttlen, cxtabtt, cdxtail, temp16b); + temp32blen = FastExpansionSumZeroElim(temp16alen, temp16a, temp16blen, temp16b, temp32b); + temp64len = FastExpansionSumZeroElim(temp32alen, temp32a, temp32blen, temp32b, temp64); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp64len, temp64, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + if (cdytail != 0.0) + { + temp16alen = ScaleExpansionZeroElim(cytablen, cytab, cdytail, temp16a); + cytabtlen = ScaleExpansionZeroElim(abtlen, abt, cdytail, cytabt); + temp32alen = ScaleExpansionZeroElim(cytabtlen, cytabt, 2.0 * cdy, temp32a); + temp48len = FastExpansionSumZeroElim(temp16alen, temp16a, temp32alen, temp32a, temp48); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp48len, temp48, finother); + finswap = finnow; finnow = finother; finother = finswap; + + + temp32alen = ScaleExpansionZeroElim(cytabtlen, cytabt, cdytail, temp32a); + cytabttlen = ScaleExpansionZeroElim(abttlen, abtt, cdytail, cytabtt); + temp16alen = ScaleExpansionZeroElim(cytabttlen, cytabtt, 2.0 * cdy, temp16a); + temp16blen = ScaleExpansionZeroElim(cytabttlen, cytabtt, cdytail, temp16b); + temp32blen = FastExpansionSumZeroElim(temp16alen, temp16a, temp16blen, temp16b, temp32b); + temp64len = FastExpansionSumZeroElim(temp32alen, temp32a, temp32blen, temp32b, temp64); + finlength = FastExpansionSumZeroElim(finlength, finnow, temp64len, temp64, finother); + finswap = finnow; finnow = finother; finother = finswap; + } + } + + return finnow[finlength - 1]; + } + + #endregion } }