shape ops work in progress

first 100,000 random cubic/cubic intersections working

git-svn-id: http://skia.googlecode.com/svn/trunk@7380 2bbb7eff-a529-9590-31e7-b0007b416f81

20 files changed, 817 insertions, 410 deletions

diff --git a/experimental/Intersection/CubeRoot.cpp b/experimental/Intersection/CubeRoot.cpp
index 82d2732f6d..5f785a0358 100644
--- a/experimental/Intersection/CubeRoot.cpp
+++ b/experimental/Intersection/CubeRoot.cpp
@@ -374,7 +374,7 @@ static int _tmain()
 #endif
 
 double cube_root(double x) {
-    if (approximately_zero(x)) {
+    if (approximately_zero_cubed(x)) {
         return 0;
     }
     double result = halley_cbrt3d(fabs(x));
diff --git a/experimental/Intersection/CubicIntersection.cpp b/experimental/Intersection/CubicIntersection.cpp
index 4ae0d84e5f..a5b4261dfe 100644
--- a/experimental/Intersection/CubicIntersection.cpp
+++ b/experimental/Intersection/CubicIntersection.cpp
@@ -10,6 +10,7 @@
 #include "Intersections.h"
 #include "IntersectionUtilities.h"
 #include "LineIntersection.h"
+#include "LineUtilities.h"
 
 static const double tClipLimit = 0.8; // http://cagd.cs.byu.edu/~tom/papers/bezclip.pdf see Multiple intersections
 
@@ -163,27 +164,49 @@ bool intersect(const Cubic& c1, const Cubic& c2, Intersections& i) {
 
 #include "CubicUtilities.h"
 
-// FIXME: ? if needed, compute the error term from the tangent length
-start here;
-// need better delta computation -- assert fails
-static double computeDelta(const Cubic& cubic, double t, double scale) {
-    double attempt = scale / precisionUnit * 2;
-#if SK_DEBUG
-    double precision = calcPrecision(cubic, t, scale);
-    _Point dxy;
+static void cubicTangent(const Cubic& cubic, double t, _Line& tangent, _Point& pt, _Point& dxy) {
+    xy_at_t(cubic, t, tangent[0].x, tangent[0].y);
+    pt = tangent[1] = tangent[0];
     dxdy_at_t(cubic, t, dxy);
-    _Point p1, p2;
-    xy_at_t(cubic, std::max(t - attempt, 0.), p1.x, p1.y);
-    xy_at_t(cubic, std::min(t + attempt, 1.), p2.x, p2.y);
-    double dx = p1.x - p2.x;
-    double dy = p1.y - p2.y;
-    double distSq = dx * dx + dy * dy;
-    double dist = sqrt(distSq);
-    assert(dist > precision);
+    tangent[0] -= dxy;
+    tangent[1] += dxy;
+}
+
+static double cubicDelta(const _Point& dxy, _Line& tangent, double scale)  {
+    double tangentLen = dxy.length();
+    tangent[0] -= tangent[1];
+    double intersectLen = tangent[0].length();
+    double result = intersectLen / tangentLen + scale;
+    return result;
+}
+
+// FIXME: after testing, make this static
+void computeDelta(const Cubic& c1, double t1, double scale1, const Cubic& c2, double t2,
+        double scale2, double& delta1, double& delta2) {
+    _Line tangent1, tangent2, line1, line2;
+    _Point dxy1, dxy2;
+    cubicTangent(c1, t1, line1, tangent1[0], dxy1);
+    cubicTangent(c2, t2, line2, tangent2[0], dxy2);
+    double range1[2], range2[2];
+    int found = intersect(line1, line2, range1, range2);
+    if (found == 0) {
+        range1[0] = 0.5;
+    } else {
+        SkASSERT(found == 1);
+    }
+    xy_at_t(line1, range1[0], tangent1[1].x, tangent1[1].y);
+#if SK_DEBUG
+    if (found == 1) {
+        xy_at_t(line2, range2[0], tangent2[1].x, tangent2[1].y);
+        SkASSERT(tangent2[1].approximatelyEqual(tangent1[1]));
+    }
 #endif
-    return attempt;
+    tangent2[1] = tangent1[1];
+    delta1 = cubicDelta(dxy1, tangent1, scale1 / precisionUnit);
+    delta2 = cubicDelta(dxy2, tangent2, scale2 / precisionUnit);
 }
 
+int debugDepth;
 // this flavor approximates the cubics with quads to find the intersecting ts
 // OPTIMIZE: if this strategy proves successful, the quad approximations, or the ts used
 // to create the approximations, could be stored in the cubic segment
@@ -252,12 +275,16 @@ static bool intersect2(const Cubic& cubic1, double t1s, double t1e, const Cubic&
                 if (p1.approximatelyEqual(p2)) {
                     i.insert(i.swapped() ? to2 : to1, i.swapped() ? to1 : to2);
                 } else {
-                    double dt1 = computeDelta(cubic1, to1, t1e - t1s);
-                    double dt2 = computeDelta(cubic2, to2, t2e - t2s);
+                    double dt1, dt2;
+                    computeDelta(cubic1, to1, (t1e - t1s), cubic2, to2, (t2e - t2s), dt1, dt2);
+                    ++debugDepth;
+                    assert(debugDepth < 10);
                     i.swap();
-                    intersect2(cubic2, std::max(to2 - dt2, 0.), std::min(to2 + dt2, 1.),
-                            cubic1, std::max(to1 - dt1, 0.), std::min(to1 + dt1, 1.), i);
+                    intersect2(cubic2, SkTMax(to2 - dt2, 0.), SkTMin(to2 + dt2, 1.),
+                            cubic1, SkTMax(to1 - dt1, 0.), SkTMin(to1 + dt1, 1.), i);
                     i.swap();
+                    --debugDepth;
+                    
                 }
             }
             t2Start = t2;
@@ -309,6 +336,7 @@ static bool intersectEnd(const Cubic& cubic1, bool start, const Cubic& cubic2, c
         tMin = std::min(tMin, local2.fT[0][index]);
         tMax = std::max(tMax, local2.fT[0][index]);
     }
+    debugDepth = 0;
     return intersect2(cubic1, start ? 0 : 1, start ? 1.0 / precisionUnit : 1 - 1.0 / precisionUnit,
             cubic2, tMin, tMax, i);
 }
@@ -318,6 +346,7 @@ static bool intersectEnd(const Cubic& cubic1, bool start, const Cubic& cubic2, c
 // line segments intersect the cubic, then use the intersections to construct a subdivision for
 // quadratic curve fitting.
 bool intersect2(const Cubic& c1, const Cubic& c2, Intersections& i) {
+    debugDepth = 0;
     bool result = intersect2(c1, 0, 1, c2, 0, 1, i);
     // FIXME: pass in cached bounds from caller
     _Rect c1Bounds, c2Bounds;
diff --git a/experimental/Intersection/CubicIntersection_Test.cpp b/experimental/Intersection/CubicIntersection_Test.cpp
index 8c2263e3a2..2e738fc0e4 100644
--- a/experimental/Intersection/CubicIntersection_Test.cpp
+++ b/experimental/Intersection/CubicIntersection_Test.cpp
@@ -57,23 +57,29 @@ void CubicIntersection_Test() {
     }
 }
 
+#define ONE_OFF_DEBUG 1
+
 static void oneOff(const Cubic& cubic1, const Cubic& cubic2) {
     SkTDArray<Quadratic> quads1;
     cubic_to_quadratics(cubic1, calcPrecision(cubic1), quads1);
+#if ONE_OFF_DEBUG
     for (int index = 0; index < quads1.count(); ++index) {
         const Quadratic& q = quads1[index];
-        SkDebugf("{{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n", q[0].x, q[0].y,
+        SkDebugf("  {{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n", q[0].x, q[0].y,
                  q[1].x, q[1].y,  q[2].x, q[2].y);
     }
     SkDebugf("\n");
+#endif
     SkTDArray<Quadratic> quads2;
     cubic_to_quadratics(cubic2, calcPrecision(cubic2), quads2);
+#if ONE_OFF_DEBUG
     for (int index = 0; index < quads2.count(); ++index) {
         const Quadratic& q = quads2[index];
-        SkDebugf("{{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n", q[0].x, q[0].y,
+        SkDebugf("  {{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n", q[0].x, q[0].y,
                  q[1].x, q[1].y,  q[2].x, q[2].y);
     }
     SkDebugf("\n");
+#endif
     Intersections intersections2;
     intersect2(cubic1, cubic2, intersections2);
     for (int pt = 0; pt < intersections2.used(); ++pt) {
@@ -83,13 +89,30 @@ static void oneOff(const Cubic& cubic1, const Cubic& cubic2) {
         int pt2 = intersections2.fFlip ? intersections2.used() - pt - 1 : pt;
         double tt2 = intersections2.fT[1][pt2];
         xy_at_t(cubic2, tt2, xy2.x, xy2.y);
+#if ONE_OFF_DEBUG
         SkDebugf("%s t1=%1.9g (%1.9g, %1.9g) (%1.9g, %1.9g) t2=%1.9g\n", __FUNCTION__,
             tt1, xy1.x, xy1.y, xy2.x, xy2.y, tt2);
+#endif
         assert(xy1.approximatelyEqual(xy2));
     }
 }
 
 static const Cubic testSet[] = {
+{{65.454505973241524, 93.881892270353575}, {45.867360264932437, 92.723972719499827}, {2.1464054482739447, 74.636369140183717}, {33.774068594804994, 40.770872887582925}},
+{{72.963387832494163, 95.659300729473728}, {11.809496633619768, 82.209921247423594}, {13.456139067865974, 57.329313623406605}, {36.060621606214262, 70.867335643091849}},
+
+{{32.484981432782945, 75.082940782924624}, {42.467313093350882, 48.131159948246157}, {3.5963115764764657, 43.208665839959245}, {79.442476890721579, 89.709102357602262}},
+{{18.98573861410177, 93.308887208490106}, {40.405250173250792, 91.039661826118675}, {8.0467721950480584, 42.100282172719147}, {40.883324221187891, 26.030185504830527}},
+
+{{7.5374809128872498, 82.441702896003477}, {22.444346930107265, 22.138854312775123}, {66.76091829629658, 50.753805856571446}, {78.193478508942519, 97.7932997968948}},
+{{97.700573130371311, 53.53260215070685}, {87.72443481149358, 84.575876772671876}, {19.215031396232092, 47.032676472809484}, {11.989686410869325, 10.659507480757082}},
+
+{{26.192053931854691, 9.8504326817814416}, {10.174241480498686, 98.476562741434464}, {21.177712558385782, 33.814968789841501}, {75.329030899018534, 55.02231980442177}},
+{{56.222082700683771, 24.54395039218662}, {95.589995289030483, 81.050822735322086}, {28.180450866082897, 28.837706255185282}, {60.128952916771617, 87.311672180570511}},
+
+{{42.449716172390481, 52.379709366885805}, {27.896043159019225, 48.797373636065686}, {92.770268299044233, 89.899302036454571}, {12.102066544863426, 99.43241951960718}},
+{{45.77532924980639, 45.958701495993274}, {37.458701356062065, 68.393691335056758}, {37.569326692060258, 27.673713456687381}, {60.674866037757539, 62.47349659096146}},
+
 {{67.426548091427676, 37.993772624988935}, {23.483695892376684, 90.476863174921306}, {35.597065061143162, 79.872482633158796}, {75.38634169631932, 18.244890038969412}},
 {{61.336508189019057, 82.693132843213675}, {44.639380902349664, 54.074825790745592}, {16.815615499771951, 20.049704667203923}, {41.866884958868326, 56.735503699973002}},
 
@@ -104,10 +127,14 @@ const size_t testSetCount = sizeof(testSet) / sizeof(testSet[0]);
 
 void CubicIntersection_OneOffTest() {
     for (size_t outer = 0; outer < testSetCount - 1; ++outer) {
+#if ONE_OFF_DEBUG
         SkDebugf("%s quads1[%d]\n", __FUNCTION__, outer);
+#endif
         const Cubic& cubic1 = testSet[outer];
         for (size_t inner = outer + 1; inner < testSetCount; ++inner) {
+#if ONE_OFF_DEBUG
         SkDebugf("%s quads2[%d]\n", __FUNCTION__, inner);
+#endif
             const Cubic& cubic2 = testSet[inner];
             oneOff(cubic1, cubic2);
         }
@@ -309,9 +336,56 @@ void CubicIntersection_RandTest() {
             int pt2 = intersections2.fFlip ? intersections2.used() - pt - 1 : pt;
             double tt2 = intersections2.fT[1][pt2];
             xy_at_t(cubic2, tt2, xy2.x, xy2.y);
+        #if 0
             SkDebugf("%s t1=%1.9g (%1.9g, %1.9g) (%1.9g, %1.9g) t2=%1.9g\n", __FUNCTION__,
                 tt1, xy1.x, xy1.y, xy2.x, xy2.y, tt2);
+        #endif
             assert(xy1.approximatelyEqual(xy2));
         }
     }
 }
+
+static Cubic deltaTestSet[] = {
+  {{1, 4}, {1, 4.*2/3}, {1, 4.*1/3}, {1, 0}},
+  {{0, 3}, {1, 2}, {2, 1}, {3, 0}},
+  {{1, 4}, {1, 4.*2/3}, {1, 4.*1/3}, {1, 0}},
+  {{3.5, 1}, {2.5, 2}, {1.5, 3}, {0.5, 4}}
+};
+
+size_t deltaTestSetLen = sizeof(deltaTestSet) / sizeof(deltaTestSet[0]);
+
+static double deltaTestSetT[] = {
+    3./8,
+    5./12,
+    6./8,
+    9./12
+};
+
+size_t deltaTestSetTLen = sizeof(deltaTestSetT) / sizeof(deltaTestSetT[0]);
+
+static double expectedT[] = {
+    0.5,
+    1./3,
+    1./8,
+    5./6
+};
+
+size_t expectedTLen = sizeof(expectedT) / sizeof(expectedT[0]);
+
+// FIXME: this test no longer valid -- does not take minimum scale contribution into account
+void CubicIntersection_ComputeDeltaTest() {
+    SkASSERT(deltaTestSetLen == deltaTestSetTLen);
+    SkASSERT(expectedTLen == deltaTestSetTLen);
+    for (size_t index = 0; index < deltaTestSetLen; index += 2) {
+        const Cubic& c1 = deltaTestSet[index];
+        const Cubic& c2 = deltaTestSet[index + 1];
+        double t1 = deltaTestSetT[index];
+        double t2 = deltaTestSetT[index + 1];
+        double d1, d2;
+        computeDelta(c1, t1, 1, c2, t2, 1, d1, d2);
+        SkASSERT(approximately_equal(t1 + d1, expectedT[index])
+            || approximately_equal(t1 - d1, expectedT[index]));
+        SkASSERT(approximately_equal(t2 + d2, expectedT[index + 1])
+            || approximately_equal(t2 - d2, expectedT[index + 1]));
+    }
+}
diff --git a/experimental/Intersection/CubicUtilities.cpp b/experimental/Intersection/CubicUtilities.cpp
index 19f16c6aaa..3e2f474d7c 100644
--- a/experimental/Intersection/CubicUtilities.cpp
+++ b/experimental/Intersection/CubicUtilities.cpp
@@ -21,7 +21,7 @@ double calcPrecision(const Cubic& cubic) {
 #if SK_DEBUG
 double calcPrecision(const Cubic& cubic, double t, double scale) {
     Cubic part;
-    sub_divide(cubic, SkMax32(0, t - scale), SkMin32(1, t + scale), part);
+    sub_divide(cubic, SkTMax(0., t - scale), SkTMin(1., t + scale), part);
     return calcPrecision(part);
 }
 #endif
@@ -41,14 +41,11 @@ void coefficients(const double* cubic, double& A, double& B, double& C, double&
 
 const double PI = 4 * atan(1);
 
-static bool is_unit_interval(double x) {
-    return x > 0 && x < 1;
-}
-
 // from SkGeometry.cpp (and Numeric Solutions, 5.6)
-int cubicRoots(double A, double B, double C, double D, double t[3]) {
+int cubicRootsValidT(double A, double B, double C, double D, double t[3]) {
+#if 0
     if (approximately_zero(A)) {  // we're just a quadratic
-        return quadraticRoots(B, C, D, t);
+        return quadraticRootsValidT(B, C, D, t);
     }
     double a, b, c;
     {
@@ -98,6 +95,113 @@ int cubicRoots(double A, double B, double C, double D, double t[3]) {
             *roots++ = r;
     }
     return (int)(roots - t);
+#else
+    double s[3];
+    int realRoots = cubicRootsReal(A, B, C, D, s);
+    int foundRoots = add_valid_ts(s, realRoots, t);
+    return foundRoots;
+#endif
+}
+
+int cubicRootsReal(double A, double B, double C, double D, double s[3]) {
+#if SK_DEBUG
+    // create a string mathematica understands
+    // GDB set print repe 15 # if repeated digits is a bother
+    //     set print elements 400 # if line doesn't fit
+    char str[1024];
+    bzero(str, sizeof(str));
+    sprintf(str, "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", A, B, C, D);
+#endif
+    if (approximately_zero(A)) {  // we're just a quadratic
+        return quadraticRootsReal(B, C, D, s);
+    }
+    if (approximately_zero(D)) { // 0 is one root
+        int num = quadraticRootsReal(A, B, C, s);
+        for (int i = 0; i < num; ++i) {
+            if (approximately_zero(s[i])) {
+                return num;
+            }
+        }
+        s[num++] = 0;
+        return num;
+    }
+    if (approximately_zero(A + B + C + D)) { // 1 is one root
+        int num = quadraticRootsReal(A, A + B, -D, s);
+        for (int i = 0; i < num; ++i) {
+            if (AlmostEqualUlps(s[i], 1)) {
+                return num;
+            }
+        }
+        s[num++] = 1;
+        return num;
+    }
+    double a, b, c;
+    {
+        double invA = 1 / A;
+        a = B * invA;
+        b = C * invA;
+        c = D * invA;
+    }
+    double a2 = a * a;
+    double Q = (a2 - b * 3) / 9;
+    double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
+    double R2 = R * R;
+    double Q3 = Q * Q * Q;
+    double R2MinusQ3 = R2 - Q3;
+    double adiv3 = a / 3;
+    double r;
+    double* roots = s;
+#if 0
+    if (approximately_zero_squared(R2MinusQ3) && AlmostEqualUlps(R2, Q3)) {
+        if (approximately_zero_squared(R)) {/* one triple solution */
+            *roots++ = -adiv3;
+        } else { /* one single and one double solution */
+
+            double u = cube_root(-R);
+            *roots++ = 2 * u - adiv3;
+            *roots++ = -u - adiv3;
+        }
+    }
+    else 
+#endif
+    if (R2MinusQ3 < 0)   // we have 3 real roots
+    {
+        double theta = acos(R / sqrt(Q3));
+        double neg2RootQ = -2 * sqrt(Q);
+
+        r = neg2RootQ * cos(theta / 3) - adiv3;
+        *roots++ = r;
+
+        r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3;
+        if (!AlmostEqualUlps(s[0], r)) {
+            *roots++ = r;
+        }
+        r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3;
+        if (!AlmostEqualUlps(s[0], r) && (roots - s == 1 || !AlmostEqualUlps(s[1], r))) {
+            *roots++ = r;
+        }
+    }
+    else                // we have 1 real root
+    {
+        double sqrtR2MinusQ3 = sqrt(R2MinusQ3);
+        double A = fabs(R) + sqrtR2MinusQ3;
+        A = cube_root(A);
+        if (R > 0) {
+            A = -A;
+        }
+        if (A != 0) {
+            A += Q / A;
+        }
+        r = A - adiv3;
+        *roots++ = r;
+        if (AlmostEqualUlps(R2, Q3)) {
+            r = -A / 2 - adiv3;
+            if (!AlmostEqualUlps(s[0], r)) {
+                *roots++ = r;
+            }
+        }
+    }
+    return (int)(roots - s);
 }
 
 // from http://www.cs.sunysb.edu/~qin/courses/geometry/4.pdf
@@ -136,14 +240,14 @@ int find_cubic_inflections(const Cubic& src, double tValues[])
     double By = src[2].y - 2 * src[1].y + src[0].y;
     double Cx = src[3].x + 3 * (src[1].x - src[2].x) - src[0].x;
     double Cy = src[3].y + 3 * (src[1].y - src[2].y) - src[0].y;
-    return quadraticRoots(Bx * Cy - By * Cx, (Ax * Cy - Ay * Cx) / 2, Ax * By - Ay * Bx, tValues);
+    return quadraticRootsValidT(Bx * Cy - By * Cx, (Ax * Cy - Ay * Cx), Ax * By - Ay * Bx, tValues);
 }
 
 bool rotate(const Cubic& cubic, int zero, int index, Cubic& rotPath) {
     double dy = cubic[index].y - cubic[zero].y;
     double dx = cubic[index].x - cubic[zero].x;
-    if (approximately_equal(dy, 0)) {
-        if (approximately_equal(dx, 0)) {
+    if (approximately_zero(dy)) {
+        if (approximately_zero(dx)) {
             return false;
         }
         memcpy(rotPath, cubic, sizeof(Cubic));
diff --git a/experimental/Intersection/CubicUtilities.h b/experimental/Intersection/CubicUtilities.h
index 5205574d29..186ed43b6a 100644
--- a/experimental/Intersection/CubicUtilities.h
+++ b/experimental/Intersection/CubicUtilities.h
@@ -15,13 +15,16 @@ double calcPrecision(const Cubic& cubic);
 double calcPrecision(const Cubic& cubic, double t, double scale);
 #endif
 void chop_at(const Cubic& src, CubicPair& dst, double t);
+// FIXME: should be private static but public here for testing
+void computeDelta(const Cubic& c1, double t1, double scale1, const Cubic& c2, double t2,
+    double scale2, double& delta1, double& delta2);
 double cube_root(double x);
 int cubic_to_quadratics(const Cubic& cubic, double precision,
         SkTDArray<Quadratic>& quadratics);
 void cubic_to_quadratics(const Cubic& cubic, double precision, SkTDArray<double>& ts);
 void coefficients(const double* cubic, double& A, double& B, double& C, double& D);
-int cubicRoots(double A, double B, double C, double D, double t[3]);
-int cubicRootsX(double A, double B, double C, double D, double s[3]);
+int cubicRootsValidT(double A, double B, double C, double D, double t[3]);
+int cubicRootsReal(double A, double B, double C, double D, double s[3]);
 void demote_cubic_to_quad(const Cubic& cubic, Quadratic& quad);
 double dx_at_t(const Cubic& , double t);
 double dy_at_t(const Cubic& , double t);
diff --git a/experimental/Intersection/DataTypes.h b/experimental/Intersection/DataTypes.h
index afdb8e9182..de763edcc2 100644
--- a/experimental/Intersection/DataTypes.h
+++ b/experimental/Intersection/DataTypes.h
@@ -24,6 +24,8 @@ inline bool AlmostEqualUlps(double A, double B) { return AlmostEqualUlps((float)
 // FIXME: delete
 int UlpsDiff(float A, float B);
 
+// FLT_EPSILON == 1.19209290E-07 == 1 / (2 ^ 23)
+const double FLT_EPSILON_CUBED = FLT_EPSILON * FLT_EPSILON * FLT_EPSILON;
 const double FLT_EPSILON_SQUARED = FLT_EPSILON * FLT_EPSILON;
 const double FLT_EPSILON_SQRT = sqrt(FLT_EPSILON);
 
@@ -42,6 +44,10 @@ inline bool approximately_zero(float x) {
     return fabs(x) < FLT_EPSILON;
 }
 
+inline bool approximately_zero_cubed(double x) {
+    return fabs(x) < FLT_EPSILON_CUBED;
+}
+
 inline bool approximately_zero_squared(double x) {
     return fabs(x) < FLT_EPSILON_SQUARED;
 }
@@ -50,8 +56,28 @@ inline bool approximately_zero_sqrt(double x) {
     return fabs(x) < FLT_EPSILON_SQRT;
 }
 
+// Use this for comparing Ts in the range of 0 to 1. For general numbers (larger and smaller) use
+// AlmostEqualUlps instead.
 inline bool approximately_equal(double x, double y) {
+#if 1
     return approximately_zero(x - y);
+#else
+// see http://visualstudiomagazine.com/blogs/tool-tracker/2011/11/compare-floating-point-numbers.aspx
+// this allows very small (e.g. degenerate) values to compare unequally, but in this case, 
+// AlmostEqualUlps should be used instead.
+    if (x == y) {
+        return true;
+    }
+    double absY = fabs(y);
+    if (x == 0) {
+        return absY < FLT_EPSILON;
+    }
+    double absX = fabs(x);
+    if (y == 0) {
+        return absX < FLT_EPSILON;
+    }
+    return fabs(x - y) < (absX > absY ? absX : absY) * FLT_EPSILON;
+#endif
 }
 
 inline bool approximately_equal_squared(double x, double y) {
@@ -160,7 +186,7 @@ struct _Point {
     }
 
     friend bool operator!=(const _Point& a, const _Point& b) {
-        return a.x!= b.x || a.y != b.y;
+        return a.x != b.x || a.y != b.y;
     }
 
     // note: this can not be implemented with
@@ -171,9 +197,22 @@ struct _Point {
                 && AlmostEqualUlps((float) y, (float) a.y);
     }
 
-    double dot(const _Point& a) {
+    double cross(const _Point& a) const {
+        return x * a.y - y * a.x;
+    }
+
+    double dot(const _Point& a) const {
         return x * a.x + y * a.y;
     }
+
+    double length() const {
+        return sqrt(lengthSquared());
+    }
+
+    double lengthSquared() const {
+        return x * x + y * y;
+    }
+
 };
 
 typedef _Point _Line[2];
@@ -243,4 +282,17 @@ struct QuadraticPair {
     _Point pts[5];
 };
 
+// FIXME: move these into SkTypes.h
+template <typename T> inline T SkTMax(T a, T b) {
+    if (a < b)
+        a = b;
+    return a;
+}
+
+template <typename T> inline T SkTMin(T a, T b) {
+    if (a > b)
+        a = b;
+    return a;
+}
+
 #endif // __DataTypes_h__
diff --git a/experimental/Intersection/Intersection_Tests.cpp b/experimental/Intersection/Intersection_Tests.cpp
index cfa340f596..bf86c07c79 100644
--- a/experimental/Intersection/Intersection_Tests.cpp
+++ b/experimental/Intersection/Intersection_Tests.cpp
@@ -15,9 +15,10 @@ void cubecode_test(int test);
 void Intersection_Tests() {
     int testsRun = 0;
 
-    CubicIntersection_RandTest();
     QuadraticIntersection_Test();
     CubicIntersection_OneOffTest();
+    QuarticRoot_Test();
+    CubicIntersection_RandTest();
     SimplifyNew_Test();
     CubicsToQuadratics_RandTest();
     CubicToQuadratics_Test();
@@ -32,7 +33,6 @@ void Intersection_Tests() {
     LineQuadraticIntersection_Test();
     MiniSimplify_Test();
     SimplifyAngle_Test();
-    QuarticRoot_Test();
     QuadraticBezierClip_Test();
     SimplifyFindNext_Test();
     SimplifyFindTop_Test();
diff --git a/experimental/Intersection/Intersection_Tests.h b/experimental/Intersection/Intersection_Tests.h
index 6fc0c754ae..4bf3068ef0 100644
--- a/experimental/Intersection/Intersection_Tests.h
+++ b/experimental/Intersection/Intersection_Tests.h
@@ -11,6 +11,7 @@ void ConvexHull_Test();
 void ConvexHull_X_Test();
 void CubicBezierClip_Test();
 void CubicCoincidence_Test();
+void CubicIntersection_ComputeDeltaTest();
 void CubicIntersection_OneOffTest();
 void CubicIntersection_Test();
 void CubicIntersection_RandTest();
diff --git a/experimental/Intersection/LineCubicIntersection.cpp b/experimental/Intersection/LineCubicIntersection.cpp
index 333e30f0e3..cc63099fd1 100644
--- a/experimental/Intersection/LineCubicIntersection.cpp
+++ b/experimental/Intersection/LineCubicIntersection.cpp
@@ -96,7 +96,7 @@ int intersectRay(double roots[3]) {
     }
     double A, B, C, D;
     coefficients(&r[0].x, A, B, C, D);
-    return cubicRoots(A, B, C, D, roots);
+    return cubicRootsValidT(A, B, C, D, roots);
 }
 
 int intersect() {
@@ -117,7 +117,7 @@ int horizontalIntersect(double axisIntercept, double roots[3]) {
     double A, B, C, D;
     coefficients(&cubic[0].y, A, B, C, D);
     D -= axisIntercept;
-    return cubicRoots(A, B, C, D, roots);
+    return cubicRootsValidT(A, B, C, D, roots);
 }
 
 int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
@@ -143,7 +143,7 @@ int verticalIntersect(double axisIntercept, double roots[3]) {
     double A, B, C, D;
     coefficients(&cubic[0].x, A, B, C, D);
     D -= axisIntercept;
-    return cubicRoots(A, B, C, D, roots);
+    return cubicRootsValidT(A, B, C, D, roots);
 }
 
 int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
diff --git a/experimental/Intersection/LineQuadraticIntersection.cpp b/experimental/Intersection/LineQuadraticIntersection.cpp
index 00c3554327..9ec3fb75e0 100644
--- a/experimental/Intersection/LineQuadraticIntersection.cpp
+++ b/experimental/Intersection/LineQuadraticIntersection.cpp
@@ -124,7 +124,7 @@ int intersectRay(double roots[2]) {
     double C = r[0];
     A += C - 2 * B; // A = a - 2*b + c
     B -= C; // B = -(b - c)
-    return quadraticRoots(A, B, C, roots);
+    return quadraticRootsValidT(A, 2 * B, C, roots);
 }
 
 int intersect() {
@@ -148,7 +148,7 @@ int horizontalIntersect(double axisIntercept, double roots[2]) {
     D += F - 2 * E; // D = d - 2*e + f
     E -= F; // E = -(d - e)
     F -= axisIntercept;
-    return quadraticRoots(D, E, F, roots);
+    return quadraticRootsValidT(D, 2 * E, F, roots);
 }
 
 int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
@@ -177,7 +177,7 @@ int verticalIntersect(double axisIntercept, double roots[2]) {
     D += F - 2 * E; // D = d - 2*e + f
     E -= F; // E = -(d - e)
     F -= axisIntercept;
-    return quadraticRoots(D, E, F, roots);
+    return quadraticRootsValidT(D, 2 * E, F, roots);
 }
 
 int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
diff --git a/experimental/Intersection/QuadraticImplicit.cpp b/experimental/Intersection/QuadraticImplicit.cpp
index 2d9d9b9a88..9cd24e9f53 100644
--- a/experimental/Intersection/QuadraticImplicit.cpp
+++ b/experimental/Intersection/QuadraticImplicit.cpp
@@ -13,8 +13,6 @@
 #include "QuadraticUtilities.h"
 #include "TSearch.h"
 
-#include <algorithm> // for std::min, max
-
 /* given the implicit form 0 = Ax^2 + Bxy + Cy^2 + Dx + Ey + F
  * and given x = at^2 + bt + c  (the parameterized form)
  *           y = dt^2 + et + f
@@ -22,14 +20,8 @@
  * 0 = A(at^2+bt+c)(at^2+bt+c)+B(at^2+bt+c)(dt^2+et+f)+C(dt^2+et+f)(dt^2+et+f)+D(at^2+bt+c)+E(dt^2+et+f)+F
  */
 
-#if SK_DEBUG
-#define QUARTIC_DEBUG 1
-#else
-#define QUARTIC_DEBUG 0
-#endif
-
 static int findRoots(const QuadImplicitForm& i, const Quadratic& q2, double roots[4],
-        bool useCubic, bool& disregardCount) {
+        bool oneHint) {
     double a, b, c;
     set_abc(&q2[0].x, a, b, c);
     double d, e, f;
@@ -56,43 +48,11 @@ static int findRoots(const QuadImplicitForm& i, const Quadratic& q2, double root
                     +     i.x()  *  c
                     +     i.y()  *  f
                     +     i.c();
-#if QUARTIC_DEBUG
-    // create a string mathematica understands
-    char str[1024];
-    bzero(str, sizeof(str));
-    sprintf(str, "Solve[%1.19g x^4 + %1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]",
-        t4, t3, t2, t1, t0);
-#endif
-    if (approximately_zero(t4)) {
-        disregardCount = true;
-        if (approximately_zero(t3)) {
-            return quadraticRootsX(t2, t1, t0, roots);
-        }
-        return cubicRootsX(t3, t2, t1, t0, roots);
-    }
-    if (approximately_zero(t0)) { // 0 is one root
-        disregardCount = true;
-        int num = cubicRootsX(t4, t3, t2, t1, roots);
-        for (int i = 0; i < num; ++i) {
-            if (approximately_zero(roots[i])) {
-                return num;
-            }
-        }
-        roots[num++] = 0;
-        return num;
-    }
-    if (useCubic) {
-        assert(approximately_zero(t4 + t3 + t2 + t1 + t0)); // 1 is one root
-        int num = cubicRootsX(t4, t4 + t3, -(t1 + t0), -t0, roots); // note that -C==A+B+D+E
-        for (int i = 0; i < num; ++i) {
-            if (approximately_equal(roots[i], 1)) {
-                return num;
-            }
-        }
-        roots[num++] = 1;
-        return num;
+    int rootCount = reducedQuarticRoots(t4, t3, t2, t1, t0, oneHint, roots);
+    if (rootCount >= 0) {
+        return rootCount;
     }
-    return quarticRoots(t4, t3, t2, t1, t0, roots);
+    return quarticRootsReal(t4, t3, t2, t1, t0, roots);
 }
 
 static void addValidRoots(const double roots[4], const int count, const int side, Intersections& i) {
@@ -196,7 +156,10 @@ static bool addIntercept(const Quadratic& q1, const Quadratic& q2, double tMin,
     line[1].y += dxdy.y;
     Intersections rootTs;
     int roots = intersect(q1, line, rootTs);
-    assert(roots == 1);
+    if (roots == 2) {
+        return false;
+    }
+    SkASSERT(roots == 1);
     _Point pt2;
     xy_at_t(q1, rootTs.fT[0][0], pt2.x, pt2.y);
     if (!pt2.approximatelyEqual(mid)) {
@@ -288,12 +251,12 @@ static bool isLinearInner(const Quadratic& q1, double t1s, double t1e, const Qua
 }
 
 static double flatMeasure(const Quadratic& q) {
-    _Point mid;
-    xy_at_t(q, 0.5, mid.x, mid.y);
-    double dx = q[2].x - q[0].x;
-    double dy = q[2].y - q[0].y;
-    double length = sqrt(dx * dx + dy * dy); // OPTIMIZE: get rid of sqrt
-    return ((mid.x - q[0].x) * dy - (mid.y - q[0].y) * dx) / length;
+    _Point mid = q[1];
+    mid -= q[0];
+    _Point dxy = q[2];
+    dxy -= q[0];
+    double length = dxy.length(); // OPTIMIZE: get rid of sqrt
+    return fabs(mid.cross(dxy) / length);
 }
 
 // FIXME ? should this measure both and then use the quad that is the flattest as the line?
@@ -306,11 +269,18 @@ static bool isLinear(const Quadratic& q1, const Quadratic& q2, Intersections& i)
     return isLinearInner(q1, 0, 1, q2, 0, 1, i);
 }
 
+// FIXME: if flat measure is sufficiently large, then probably the quartic solution failed
 static bool relaxedIsLinear(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
     double m1 = flatMeasure(q1);
     double m2 = flatMeasure(q2);
+#if SK_DEBUG
+    double min = SkTMin(m1, m2);
+    if (min > 5) {
+        SkDebugf("%s maybe not flat enough.. %1.9g\n", __FUNCTION__, min);
+    }
+#endif
     i.reset();
-    if (fabs(m1) < fabs(m2)) {
+    if (m1 < m2) {
         isLinearInner(q1, 0, 1, q2, 0, 1, i);
         return false;
     } else {
@@ -400,18 +370,10 @@ bool intersect2(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
         return i.fCoincidentUsed > 0;
     }
     double roots1[4], roots2[4];
-    bool disregardCount1 = false;
-    bool disregardCount2 = false;
     bool useCubic = q1[0] == q2[0] || q1[0] == q2[2] || q1[2] == q2[0];
-    int rootCount = findRoots(i2, q1, roots1, useCubic, disregardCount1);
+    int rootCount = findRoots(i2, q1, roots1, useCubic);
     // OPTIMIZATION: could short circuit here if all roots are < 0 or > 1
-    int rootCount2 = findRoots(i1, q2, roots2, useCubic, disregardCount2);
- #if 0
-    if (rootCount != rootCount2 && !disregardCount1 && !disregardCount2) {
-        unsortableExpanse(q1, q2, i);
-        return false;
-    }
- #endif
+    int rootCount2 = findRoots(i1, q2, roots2, useCubic);
     addValidRoots(roots1, rootCount, 0, i);
     addValidRoots(roots2, rootCount2, 1, i);
     if (i.insertBalanced() && i.fUsed <= 1) {
diff --git a/experimental/Intersection/QuadraticIntersection_Test.cpp b/experimental/Intersection/QuadraticIntersection_Test.cpp
index 432f614282..055a18e825 100644
--- a/experimental/Intersection/QuadraticIntersection_Test.cpp
+++ b/experimental/Intersection/QuadraticIntersection_Test.cpp
@@ -59,6 +59,21 @@ static void standardTestCases() {
 }
 
 static const Quadratic testSet[] = {
+  {{41.5072916,87.1234036}, {28.2747836,80.9545395}, {23.5780771,69.3344126}},
+  {{72.9633878,95.6593007}, {42.7738746,88.4730382}, {31.1932785,80.2458029}},
+
+  {{31.1663962,54.7302484}, {31.1662882,54.7301074}, {31.1663969,54.7302485}},
+  {{26.0404936,45.4260361}, {27.7887523,33.1863051}, {40.8833242,26.0301855}},
+
+  {{29.9404074,49.1672596}, {44.3131071,45.3915253}, {58.1067559,59.5061814}},
+  {{72.6510251,64.2972928}, {53.6989659,60.1862397}, {35.2053722,44.8391126}},
+
+{{52.14807018377202, 65.012420045148644}, {44.778669050208237, 66.315562705604378}, {51.619118408823567, 63.787827046262684}},
+{{30.004993234763383, 93.921296668202288}, {53.384822003076991, 60.732180341802753}, {58.652998934338584, 43.111073088306185}},
+
+{{80.897794748143198, 49.236332042718459}, {81.082078218891212, 64.066749904488631}, {69.972305057149981, 72.968595519850993}},
+{{72.503745601281395, 32.952320736577882}, {88.030880716061645, 38.137194847810164}, {73.193774825517906, 67.773492479591397}},
+
 {{67.426548091427676, 37.993772624988935}, {51.129513170665035, 57.542281234563646}, {44.594748190899189, 65.644267382683879}},
 {{61.336508189019057, 82.693132843213675}, {54.825078921449354, 71.663932799212432}, {47.727444217558926, 61.4049645128392}},
 
@@ -119,37 +134,47 @@ static const Quadratic testSet[] = {
 
 const size_t testSetCount = sizeof(testSet) / sizeof(testSet[0]);
 
+#define ONE_OFF_DEBUG 0
+
+static void oneOffTest1(size_t outer, size_t inner) {
+    const Quadratic& quad1 = testSet[outer];
+    const Quadratic& quad2 = testSet[inner];
+    Intersections intersections2;
+    intersect2(quad1, quad2, intersections2);
+    if (intersections2.fUnsortable) {
+        SkASSERT(0);
+        return;
+    }
+    for (int pt = 0; pt < intersections2.used(); ++pt) {
+        double tt1 = intersections2.fT[0][pt];
+        double tx1, ty1;
+        xy_at_t(quad1, tt1, tx1, ty1);
+        int pt2 = intersections2.fFlip ? intersections2.used() - pt - 1 : pt;
+        double tt2 = intersections2.fT[1][pt2];
+        double tx2, ty2;
+        xy_at_t(quad2, tt2, tx2, ty2);
+        if (!AlmostEqualUlps(tx1, tx2)) {
+            SkDebugf("%s [%d,%d] x!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
+                __FUNCTION__, (int)outer, (int)inner, tt1, tx1, ty1, tt2, tx2, ty2);
+            SkASSERT(0);
+        }
+        if (!AlmostEqualUlps(ty1, ty2)) {
+            SkDebugf("%s [%d,%d] y!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
+                __FUNCTION__, (int)outer, (int)inner, tt1, tx1, ty1, tt2, tx2, ty2);
+            SkASSERT(0);
+        }
+#if ONE_OFF_DEBUG
+        SkDebugf("%s [%d][%d] t1=%1.9g (%1.9g, %1.9g) t2=%1.9g\n", __FUNCTION__,
+            outer, inner, tt1, tx1, tx2, tt2);
+#endif
+    }
+}
+
 static void oneOffTest() {
+//    oneOffTest1(0, 1);
     for (size_t outer = 0; outer < testSetCount - 1; ++outer) {
         for (size_t inner = outer + 1; inner < testSetCount; ++inner) {
-            const Quadratic& quad1 = testSet[outer];
-            const Quadratic& quad2 = testSet[inner];
-            Intersections intersections2;
-            intersect2(quad1, quad2, intersections2);
-            if (intersections2.fUnsortable) {
-                continue;
-            }
-            for (int pt = 0; pt < intersections2.used(); ++pt) {
-                double tt1 = intersections2.fT[0][pt];
-                double tx1, ty1;
-                xy_at_t(quad1, tt1, tx1, ty1);
-                int pt2 = intersections2.fFlip ? intersections2.used() - pt - 1 : pt;
-                double tt2 = intersections2.fT[1][pt2];
-                double tx2, ty2;
-                xy_at_t(quad2, tt2, tx2, ty2);
-                if (!AlmostEqualUlps(tx1, tx2)) {
-                    SkDebugf("%s [%d,%d] x!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
-                        __FUNCTION__, (int)outer, (int)inner, tt1, tx1, ty1, tt2, tx2, ty2);
-                    SkASSERT(0);
-                }
-                if (!AlmostEqualUlps(ty1, ty2)) {
-                    SkDebugf("%s [%d,%d] y!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
-                        __FUNCTION__, (int)outer, (int)inner, tt1, tx1, ty1, tt2, tx2, ty2);
-                    SkASSERT(0);
-                }
-                SkDebugf("%s [%d][%d] t1=%1.9g (%1.9g, %1.9g) t2=%1.9g\n", __FUNCTION__,
-                    outer, inner, tt1, tx1, tx2, tt2);
-            }
+            oneOffTest1(outer, inner);
         }
     }
 }
diff --git a/experimental/Intersection/QuadraticReduceOrder.cpp b/experimental/Intersection/QuadraticReduceOrder.cpp
index aa1d0585dd..096f6b6643 100644
--- a/experimental/Intersection/QuadraticReduceOrder.cpp
+++ b/experimental/Intersection/QuadraticReduceOrder.cpp
@@ -92,7 +92,7 @@ static int check_linear(const Quadratic& quad, Quadratic& reduction,
     if (root) {
         _Point extrema;
         extrema.x = interp_quad_coords(quad[0].x, quad[1].x, quad[2].x, tValue);
-        extrema.y = interp_quad_coords(quad[0].x, quad[1].x, quad[2].x, tValue);
+        extrema.y = interp_quad_coords(quad[0].y, quad[1].y, quad[2].y, tValue);
         // sameSide > 0 means mid is smaller than either [0] or [2], so replace smaller
         int replace;
         if (useX) {
diff --git a/experimental/Intersection/QuadraticUtilities.cpp b/experimental/Intersection/QuadraticUtilities.cpp
index 8b02b4ce6e..475ec13b7d 100644
--- a/experimental/Intersection/QuadraticUtilities.cpp
+++ b/experimental/Intersection/QuadraticUtilities.cpp
@@ -5,6 +5,7 @@
  * found in the LICENSE file.
  */
 #include "QuadraticUtilities.h"
+#include "SkTypes.h"
 #include <math.h>
 
 /*
@@ -20,12 +21,37 @@ and using the roots
 
 */
 
+    
+int add_valid_ts(double s[], int realRoots, double* t) {
+    int foundRoots = 0;
+    for (int index = 0; index < realRoots; ++index) {
+        double tValue = s[index];
+        if (approximately_zero_or_more(tValue) && approximately_one_or_less(tValue)) {
+            if (approximately_less_than_zero(tValue)) {
+                tValue = 0;
+            } else if (approximately_greater_than_one(tValue)) {
+                tValue = 1;
+            }
+            for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
+                if (approximately_equal(t[idx2], tValue)) {
+                    goto nextRoot;
+                }
+            }
+            t[foundRoots++] = tValue;
+        }
+nextRoot:
+        ;
+    }
+    return foundRoots;
+}
+
 // note: caller expects multiple results to be sorted smaller first
 // note: http://en.wikipedia.org/wiki/Loss_of_significance has an interesting
 //  analysis of the quadratic equation, suggesting why the following looks at
 //  the sign of B -- and further suggesting that the greatest loss of precision
 //  is in b squared less two a c
-int quadraticRoots(double A, double B, double C, double t[2]) {
+int quadraticRootsValidT(double A, double B, double C, double t[2]) {
+#if 0
     B *= 2;
     double square = B * B - 4 * A * C;
     if (approximately_negative(square)) {
@@ -61,9 +87,64 @@ int quadraticRoots(double A, double B, double C, double t[2]) {
             t[0] = ratio;
         }
     }
+#else
+    double s[2];
+    int realRoots = quadraticRootsReal(A, B, C, s);
+    int foundRoots = add_valid_ts(s, realRoots, t);
+#endif
     return foundRoots;
 }
 
+// unlike quadratic roots, this does not discard real roots <= 0 or >= 1
+int quadraticRootsReal(const double A, const double B, const double C, double s[2]) {
+    if (approximately_zero(A)) {
+        if (approximately_zero(B)) {
+            s[0] = 0;
+            return C == 0;
+        }
+        s[0] = -C / B;
+        return 1;
+    }
+    /* normal form: x^2 + px + q = 0 */
+    const double p = B / (2 * A);
+    const double q = C / A;
+    const double p2 = p * p;
+#if 0
+    double D = AlmostEqualUlps(p2, q) ? 0 : p2 - q;
+    if (D <= 0) {
+        if (D < 0) {
+            return 0;
+        }
+        s[0] = -p;
+        SkDebugf("[%d] %1.9g\n", 1, s[0]);
+        return 1;
+    }
+    double sqrt_D = sqrt(D);
+    s[0] = sqrt_D - p;
+    s[1] = -sqrt_D - p;
+    SkDebugf("[%d] %1.9g %1.9g\n", 2, s[0], s[1]);
+    return 2;
+#else
+    if (!AlmostEqualUlps(p2, q) && p2 < q) {
+        return 0;
+    }
+    double sqrt_D = 0;
+    if (p2 > q) {
+        sqrt_D = sqrt(p2 - q);
+    }
+    s[0] = sqrt_D - p;
+    s[1] = -sqrt_D - p;
+#if 0
+    if (AlmostEqualUlps(s[0], s[1])) {
+        SkDebugf("[%d] %1.9g\n", 1, s[0]);
+    } else {
+        SkDebugf("[%d] %1.9g %1.9g\n", 2, s[0], s[1]);
+    }
+#endif
+    return 1 + !AlmostEqualUlps(s[0], s[1]);
+#endif
+}
+
 static double derivativeAtT(const double* quad, double t) {
     double a = t - 1;
     double b = 1 - 2 * t;
diff --git a/experimental/Intersection/QuadraticUtilities.h b/experimental/Intersection/QuadraticUtilities.h
index 0fe8be8371..f423dc6ecc 100644
--- a/experimental/Intersection/QuadraticUtilities.h
+++ b/experimental/Intersection/QuadraticUtilities.h
@@ -10,6 +10,7 @@
 
 #include "DataTypes.h"
 
+int add_valid_ts(double s[], int realRoots, double* t);
 void chop_at(const Quadratic& src, QuadraticPair& dst, double t);
 double dx_at_t(const Quadratic& , double t);
 double dy_at_t(const Quadratic& , double t);
@@ -30,8 +31,8 @@ inline void set_abc(const double* quad, double& a, double& b, double& c) {
     b -= c;          // b =     2*B - 2*C
 }
 
-int quadraticRoots(double A, double B, double C, double t[2]);
-int quadraticRootsX(const double A, const double B, const double C, double s[2]);
+int quadraticRootsReal(double A, double B, double C, double t[2]);
+int quadraticRootsValidT(const double A, const double B, const double C, double s[2]);
 void sub_divide(const Quadratic& src, double t1, double t2, Quadratic& dst);
 void xy_at_t(const Quadratic& , double t, double& x, double& y);
 
diff --git a/experimental/Intersection/QuarticRoot.cpp b/experimental/Intersection/QuarticRoot.cpp
index 86ea7a63fd..6941935f43 100644
--- a/experimental/Intersection/QuarticRoot.cpp
+++ b/experimental/Intersection/QuarticRoot.cpp
@@ -30,190 +30,48 @@
 #include "QuadraticUtilities.h"
 #include "QuarticRoot.h"
 
+int reducedQuarticRoots(const double t4, const double t3, const double t2, const double t1,
+        const double t0, const bool oneHint, double roots[4]) {
 #if SK_DEBUG
-#define QUARTIC_DEBUG 1
-#else
-#define QUARTIC_DEBUG 0
-#endif
-
-const double PI = 4 * atan(1);
-
-// unlike quadraticRoots in QuadraticUtilities.cpp, this does not discard
-// real roots <= 0 or >= 1
-int quadraticRootsX(const double A, const double B, const double C,
-        double s[2]) {
-    if (approximately_zero(A)) {
-        if (approximately_zero(B)) {
-            s[0] = 0;
-            return C == 0;
-        }
-        s[0] = -C / B;
-        return 1;
-    }
-    /* normal form: x^2 + px + q = 0 */
-    const double p = B / (2 * A);
-    const double q = C / A;
-    double D = p * p - q;
-    if (D < 0) {
-        if (approximately_positive_squared(D)) {
-            D = 0;
-        } else {
-            return 0;
-        }
-    }
-    double sqrt_D = sqrt(D);
-    if (approximately_less_than_zero(sqrt_D)) {
-        s[0] = -p;
-        return 1;
-    }
-    s[0] = sqrt_D - p;
-    s[1] = -sqrt_D - p;
-    return 2;
-}
-
-#define USE_GEMS 0
-#if USE_GEMS
-// unlike cubicRoots in CubicUtilities.cpp, this does not discard
-// real roots <= 0 or >= 1
-int cubicRootsX(const double A, const double B, const double C,
-        const double D, double s[3]) {
-    int num;
-    /* normal form: x^3 + Ax^2 + Bx + C = 0 */
-    const double invA = 1 / A;
-    const double a = B * invA;
-    const double b = C * invA;
-    const double c = D * invA;
-    /*  substitute x = y - a/3 to eliminate quadric term:
-    x^3 +px + q = 0 */
-    const double a2 = a * a;
-    const double Q = (-a2 + b * 3) / 9;
-    const double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
-    /* use Cardano's formula */
-    const double Q3 = Q * Q * Q;
-    const double R2plusQ3 = R * R + Q3;
-    if (approximately_zero(R2plusQ3)) {
-        if (approximately_zero(R)) {/* one triple solution */
-            s[0] = 0;
-            num = 1;
-        } else { /* one single and one double solution */
-
-            double u = cube_root(-R);
-            s[0] = 2 * u;
-            s[1] = -u;
-            num = 2;
-        }
-    }
-    else if (R2plusQ3 < 0) { /* Casus irreducibilis: three real solutions */
-        const double theta = acos(-R / sqrt(-Q3)) / 3;
-        const double _2RootQ = 2 * sqrt(-Q);
-        s[0] = _2RootQ * cos(theta);
-        s[1] = -_2RootQ * cos(theta + PI / 3);
-        s[2] = -_2RootQ * cos(theta - PI / 3);
-        num = 3;
-    } else { /* one real solution */
-        const double sqrt_D = sqrt(R2plusQ3);
-        const double u = cube_root(sqrt_D - R);
-        const double v = -cube_root(sqrt_D + R);
-        s[0] = u + v;
-        num = 1;
-    }
-    /* resubstitute */
-    const double sub = a / 3;
-    for (int i = 0; i < num; ++i) {
-        s[i] -= sub;
-    }
-    return num;
-}
-#else
-
-int cubicRootsX(double A, double B, double C, double D, double s[3]) {
-#if QUARTIC_DEBUG
     // create a string mathematica understands
+    // GDB set print repe 15 # if repeated digits is a bother
+    //     set print elements 400 # if line doesn't fit
     char str[1024];
     bzero(str, sizeof(str));
-    sprintf(str, "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", A, B, C, D);
+    sprintf(str, "Solve[%1.19g x^4 + %1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]",
+        t4, t3, t2, t1, t0);
 #endif
-    if (approximately_zero(A)) {  // we're just a quadratic
-        return quadraticRootsX(B, C, D, s);
+    if (approximately_zero(t4)) {
+        if (approximately_zero(t3)) {
+            return quadraticRootsReal(t2, t1, t0, roots);
+        }
+        return cubicRootsReal(t3, t2, t1, t0, roots);
     }
-    if (approximately_zero(D)) { // 0 is one root
-        int num = quadraticRootsX(A, B, C, s);
+    if (approximately_zero(t0)) { // 0 is one root
+        int num = cubicRootsReal(t4, t3, t2, t1, roots);
         for (int i = 0; i < num; ++i) {
-            if (approximately_zero(s[i])) {
+            if (approximately_zero(roots[i])) {
                 return num;
             }
         }
-        s[num++] = 0;
+        roots[num++] = 0;
         return num;
     }
-    if (approximately_zero(A + B + C + D)) { // 1 is one root
-        int num = quadraticRootsX(A, A + B, -D, s);
+    if (oneHint) {
+        assert(approximately_zero(t4 + t3 + t2 + t1 + t0)); // 1 is one root
+        int num = cubicRootsReal(t4, t4 + t3, -(t1 + t0), -t0, roots); // note that -C==A+B+D+E
         for (int i = 0; i < num; ++i) {
-            if (approximately_equal(s[i], 1)) {
+            if (approximately_equal(roots[i], 1)) {
                 return num;
             }
         }
-        s[num++] = 1;
+        roots[num++] = 1;
         return num;
     }
-    double a, b, c;
-    {
-        double invA = 1 / A;
-        a = B * invA;
-        b = C * invA;
-        c = D * invA;
-    }
-    double a2 = a * a;
-    double Q = (a2 - b * 3) / 9;
-    double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
-    double Q3 = Q * Q * Q;
-    double R2MinusQ3 = R * R - Q3;
-    double adiv3 = a / 3;
-    double r;
-    double* roots = s;
-
-    if (approximately_zero_squared(R2MinusQ3)) {
-        if (approximately_zero(R)) {/* one triple solution */
-            *roots++ = -adiv3;
-        } else { /* one single and one double solution */
-
-            double u = cube_root(-R);
-            *roots++ = 2 * u - adiv3;
-            *roots++ = -u - adiv3;
-        }
-    }
-    else if (R2MinusQ3 < 0)   // we have 3 real roots
-    {
-        double theta = acos(R / sqrt(Q3));
-        double neg2RootQ = -2 * sqrt(Q);
-
-        r = neg2RootQ * cos(theta / 3) - adiv3;
-        *roots++ = r;
-
-        r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3;
-        *roots++ = r;
-
-        r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3;
-        *roots++ = r;
-    }
-    else                // we have 1 real root
-    {
-        double A = fabs(R) + sqrt(R2MinusQ3);
-        A = cube_root(A);
-        if (R > 0) {
-            A = -A;
-        }
-        if (A != 0) {
-            A += Q / A;
-        }
-        r = A - adiv3;
-        *roots++ = r;
-    }
-    return (int)(roots - s);
+    return -1;
 }
-#endif
 
-int quarticRoots(const double A, const double B, const double C, const double D,
+int quarticRootsReal(const double A, const double B, const double C, const double D,
         const double E, double s[4]) {
     double  u, v;
     /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */
@@ -231,37 +89,97 @@ int quarticRoots(const double A, const double B, const double C, const double D,
     int num;
     if (approximately_zero(r)) {
     /* no absolute term: y(y^3 + py + q) = 0 */
-        num = cubicRootsX(1, 0, p, q, s);
+        num = cubicRootsReal(1, 0, p, q, s);
         s[num++] = 0;
     } else {
         /* solve the resolvent cubic ... */
-        (void) cubicRootsX(1, -p / 2, -r, r * p / 2 - q * q / 8, s);
+        double cubicRoots[3];
+        int roots = cubicRootsReal(1, -p / 2, -r, r * p / 2 - q * q / 8, cubicRoots);
+        int index;
+    #if 0 && SK_DEBUG // enable to verify that any cubic root is as good as any other
+        double tries[3][4];
+        int nums[3];
+        for (index = 0; index < roots; ++index) {
+            /* ... and take one real solution ... */
+            const double z = cubicRoots[index];
+            /* ... to build two quadric equations */
+            u = z * z - r;
+            v = 2 * z - p;
+            if (approximately_zero_squared(u)) {
+                u = 0;
+            } else if (u > 0) {
+                u = sqrt(u);
+            } else {
+                SkDebugf("%s u=%1.9g <0\n", __FUNCTION__, u);
+                continue;
+            }
+            if (approximately_zero_squared(v)) {
+                v = 0;
+            } else if (v > 0) {
+                v = sqrt(v);
+            } else {
+                SkDebugf("%s v=%1.9g <0\n", __FUNCTION__, v);
+                continue;
+            }
+            nums[index] = quadraticRootsReal(1, q < 0 ? -v : v, z - u, tries[index]);
+            nums[index] += quadraticRootsReal(1, q < 0 ? v : -v, z + u, tries[index] + nums[index]);
+            /* resubstitute */
+            const double sub = a / 4;
+            for (int i = 0; i < nums[index]; ++i) {
+                tries[index][i] -= sub;
+            }
+        }
+        for (index = 0; index < roots; ++index) {
+            SkDebugf("%s", __FUNCTION__);
+            for (int idx2 = 0; idx2 < nums[index]; ++idx2) {
+                SkDebugf(" %1.9g", tries[index][idx2]);
+            }
+            SkDebugf("\n");
+        }
+    #endif
         /* ... and take one real solution ... */
-        const double z = s[0];
-        /* ... to build two quadric equations */
-        u = z * z - r;
-        v = 2 * z - p;
-        if (approximately_zero_squared(u)) {
-            u = 0;
-        } else if (u > 0) {
-            u = sqrt(u);
-        } else {
-            return 0;
+        double z;
+        num = 0;
+        int num2 = 0;
+        for (index = 0; index < roots; ++index) {
+            z = cubicRoots[index];
+            /* ... to build two quadric equations */
+            u = z * z - r;
+            v = 2 * z - p;
+            if (approximately_zero_squared(u)) {
+                u = 0;
+            } else if (u > 0) {
+                u = sqrt(u);
+            } else {
+                continue;
+            }
+            if (approximately_zero_squared(v)) {
+                v = 0;
+            } else if (v > 0) {
+                v = sqrt(v);
+            } else {
+                continue;
+            }
+            num = quadraticRootsReal(1, q < 0 ? -v : v, z - u, s);
+            num2 = quadraticRootsReal(1, q < 0 ? v : -v, z + u, s + num);
+            if (!((num | num2) & 1)) {
+                break; // prefer solutions without single quad roots
+            }
         }
-        if (approximately_zero_squared(v)) {
-            v = 0;
-        } else if (v > 0) {
-            v = sqrt(v);
-        } else {
-            return 0;
+        num += num2;
+        if (!num) {
+            return 0; // no valid cubic root
         }
-        num = quadraticRootsX(1, q < 0 ? -v : v, z - u, s);
-        num += quadraticRootsX(1, q < 0 ? v : -v, z + u, s + num);
+    }
+    /* resubstitute */
+    const double sub = a / 4;
+    for (int i = 0; i < num; ++i) {
+        s[i] -= sub;
     }
     // eliminate duplicates
     for (int i = 0; i < num - 1; ++i) {
         for (int j = i + 1; j < num; ) {
-            if (approximately_equal(s[i], s[j])) {
+            if (AlmostEqualUlps(s[i], s[j])) {
                 if (j < --num) {
                     s[j] = s[num];
                 }
@@ -270,10 +188,5 @@ int quarticRoots(const double A, const double B, const double C, const double D,
             }
         }
     }
-    /* resubstitute */
-    const double sub = a / 4;
-    for (int i = 0; i < num; ++i) {
-        s[i] -= sub;
-    }
     return num;
 }
diff --git a/experimental/Intersection/QuarticRoot.h b/experimental/Intersection/QuarticRoot.h
index 8baa842f3c..f76509e2f4 100644
--- a/experimental/Intersection/QuarticRoot.h
+++ b/experimental/Intersection/QuarticRoot.h
@@ -1,2 +1,5 @@
-int quarticRoots(const double A, const double B, const double C, const double D,
+int reducedQuarticRoots(const double t4, const double t3, const double t2, const double t1,
+        const double t0, const bool oneHint, double s[4]);
+        
+int quarticRootsReal(const double A, const double B, const double C, const double D,
         const double E, double s[4]);
diff --git a/experimental/Intersection/QuarticRoot_Test.cpp b/experimental/Intersection/QuarticRoot_Test.cpp
index 0f310bb08e..027d19bd75 100644
--- a/experimental/Intersection/QuarticRoot_Test.cpp
+++ b/experimental/Intersection/QuarticRoot_Test.cpp
@@ -17,23 +17,37 @@ double rootE[] = {-5, -1, 0, 1, 7};
 size_t rootECount = sizeof(rootE) / sizeof(rootE[0]);
 
 
-static void quadraticTest() {
+static void quadraticTest(bool limit) {
     // (x - a)(x - b) == x^2 - (a + b)x + ab
     for (size_t aIndex = 0; aIndex < mulACount; ++aIndex) {
         for (size_t bIndex = 0; bIndex < rootBCount; ++bIndex) {
             for (size_t cIndex = 0; cIndex < rootCCount; ++cIndex) {
                 const double A = mulA[aIndex];
-                const double B = rootB[bIndex];
-                const double C = rootC[cIndex];
+                double B = rootB[bIndex];
+                double C = rootC[cIndex];
+                if (limit) {
+                    B = (B - 6) / 12;
+                    C = (C - 6) / 12;
+                }
                 const double b = A * (B + C);
                 const double c = A * B * C;
                 double roots[2];
-                const int rootCount = quadraticRootsX(A, b, c, roots);
-                const int expected = 1 + (B != C);
+                const int rootCount = limit ? quadraticRootsValidT(A, b, c, roots)
+                    : quadraticRootsReal(A, b, c, roots);
+                int expected;
+                if (limit) {
+                    expected = B <= 0 && B >= -1;
+                    expected += B != C && C <= 0 && C >= -1;
+                } else {
+                    expected = 1 + (B != C);
+                }
                 assert(rootCount == expected);
+                if (!rootCount) {
+                    continue;
+                }
                 assert(approximately_equal(roots[0], -B)
                         || approximately_equal(roots[0], -C));
-                if (B != C) {
+                if (expected > 1) {
                     assert(!approximately_equal(roots[0], roots[1]));
                     assert(approximately_equal(roots[1], -B)
                             || approximately_equal(roots[1], -C));
@@ -43,45 +57,119 @@ static void quadraticTest() {
     }
 }
 
-static void cubicTest() {
+static void testOneCubic(bool limit, size_t aIndex, size_t bIndex, size_t cIndex, size_t dIndex) {
+    const double A = mulA[aIndex];
+    double B = rootB[bIndex];
+    double C = rootC[cIndex];
+    double D = rootD[dIndex];
+    if (limit) {
+        B = (B - 6) / 12;
+        C = (C - 6) / 12;
+        D = (C - 2) / 6;
+    }
+    const double b = A * (B + C + D);
+    const double c = A * (B * C + C * D + B * D);
+    const double d = A * B * C * D;
+    double roots[3];
+    const int rootCount = limit ? cubicRootsValidT(A, b, c, d, roots)
+            : cubicRootsReal(A, b, c, d, roots);
+    int expected;
+    if (limit) {
+        expected = B <= 0 && B >= -1;
+        expected += B != C && C <= 0 && C >= -1;
+        expected += B != D && C != D && D <= 0 && D >= -1;
+    } else {
+        expected = 1 + (B != C) + (B != D && C != D);
+    }
+    assert(rootCount == expected);
+    if (!rootCount) {
+        return;
+    }
+    assert(approximately_equal(roots[0], -B)
+            || approximately_equal(roots[0], -C)
+            || approximately_equal(roots[0], -D));
+    if (expected <= 1) {
+        return;
+    }
+    assert(!approximately_equal(roots[0], roots[1]));
+    assert(approximately_equal(roots[1], -B)
+            || approximately_equal(roots[1], -C)
+            || approximately_equal(roots[1], -D));
+    if (expected <= 2) {
+        return;
+    }
+    assert(!approximately_equal(roots[0], roots[2])
+            && !approximately_equal(roots[1], roots[2]));
+    assert(approximately_equal(roots[2], -B)
+            || approximately_equal(roots[2], -C)
+            || approximately_equal(roots[2], -D));
+}
+
+static void cubicTest(bool limit) {
     // (x - a)(x - b)(x - c) == x^3 - (a + b + c)x^2 + (ab + bc + ac)x - abc
     for (size_t aIndex = 0; aIndex < mulACount; ++aIndex) {
         for (size_t bIndex = 0; bIndex < rootBCount; ++bIndex) {
             for (size_t cIndex = 0; cIndex < rootCCount; ++cIndex) {
                 for (size_t dIndex = 0; dIndex < rootDCount; ++dIndex) {
-                    const double A = mulA[aIndex];
-                    const double B = rootB[bIndex];
-                    const double C = rootC[cIndex];
-                    const double D = rootD[dIndex];
-                    const double b = A * (B + C + D);
-                    const double c = A * (B * C + C * D + B * D);
-                    const double d = A * B * C * D;
-                    double roots[3];
-                    const int rootCount = cubicRootsX(A, b, c, d, roots);
-                    const int expected = 1 + (B != C) + (B != D && C != D);
-                    assert(rootCount == expected);
-                    assert(approximately_equal(roots[0], -B)
-                            || approximately_equal(roots[0], -C)
-                            || approximately_equal(roots[0], -D));
-                    if (expected > 1) {
-                        assert(!approximately_equal(roots[0], roots[1]));
-                        assert(approximately_equal(roots[1], -B)
-                                || approximately_equal(roots[1], -C)
-                                || approximately_equal(roots[1], -D));
-                        if (expected > 2) {
-                            assert(!approximately_equal(roots[0], roots[2])
-                                    && !approximately_equal(roots[1], roots[2]));
-                            assert(approximately_equal(roots[2], -B)
-                                    || approximately_equal(roots[2], -C)
-                                    || approximately_equal(roots[2], -D));
-                        }
-                    }
+                    testOneCubic(limit, aIndex, bIndex, cIndex, dIndex);
                 }
             }
         }
     }
 }
 
+static void testOneQuartic(size_t aIndex, size_t bIndex, size_t cIndex, size_t dIndex,
+        size_t eIndex) {
+    const double A = mulA[aIndex];
+    const double B = rootB[bIndex];
+    const double C = rootC[cIndex];
+    const double D = rootD[dIndex];
+    const double E = rootE[eIndex];
+    const double b = A * (B + C + D + E);
+    const double c = A * (B * C + C * D + B * D + B * E + C * E + D * E);
+    const double d = A * (B * C * D + B * C * E + B * D * E + C * D * E);
+    const double e = A * B * C * D * E;
+    double roots[4];
+    bool oneHint = approximately_zero(A + b + c + d + e);
+    int rootCount = reducedQuarticRoots(A, b, c, d, e, oneHint, roots);
+    if (rootCount < 0) {
+        rootCount = quarticRootsReal(A, b, c, d, e, roots);
+    }
+    const int expected = 1 + (B != C) + (B != D && C != D) + (B != E && C != E && D != E);
+    assert(rootCount == expected);
+    assert(AlmostEqualUlps(roots[0], -B)
+            || AlmostEqualUlps(roots[0], -C)
+            || AlmostEqualUlps(roots[0], -D)
+            || AlmostEqualUlps(roots[0], -E));
+    if (expected <= 1) {
+        return;
+    }
+    assert(!AlmostEqualUlps(roots[0], roots[1]));
+    assert(AlmostEqualUlps(roots[1], -B)
+            || AlmostEqualUlps(roots[1], -C)
+            || AlmostEqualUlps(roots[1], -D)
+            || AlmostEqualUlps(roots[1], -E));
+    if (expected <= 2) {
+        return;
+    }
+    assert(!AlmostEqualUlps(roots[0], roots[2])
+            && !AlmostEqualUlps(roots[1], roots[2]));
+    assert(AlmostEqualUlps(roots[2], -B)
+            || AlmostEqualUlps(roots[2], -C)
+            || AlmostEqualUlps(roots[2], -D)
+            || AlmostEqualUlps(roots[2], -E));
+    if (expected <= 3) {
+        return;
+    }
+    assert(!AlmostEqualUlps(roots[0], roots[3])
+            && !AlmostEqualUlps(roots[1], roots[3])
+            && !AlmostEqualUlps(roots[2], roots[3]));
+    assert(AlmostEqualUlps(roots[3], -B)
+            || AlmostEqualUlps(roots[3], -C)
+            || AlmostEqualUlps(roots[3], -D)
+            || AlmostEqualUlps(roots[3], -E));
+}
+
 static void quarticTest() {
     // (x - a)(x - b)(x - c)(x - d) == x^4 - (a + b + c + d)x^3
     //   + (ab + bc + cd + ac + bd + cd)x^2 - (abc + bcd + abd + acd) * x + abcd
@@ -90,47 +178,7 @@ static void quarticTest() {
             for (size_t cIndex = 0; cIndex < rootCCount; ++cIndex) {
                 for (size_t dIndex = 0; dIndex < rootDCount; ++dIndex) {
                     for (size_t eIndex = 0; eIndex < rootECount; ++eIndex) {
-                        const double A = mulA[aIndex];
-                        const double B = rootB[bIndex];
-                        const double C = rootC[cIndex];
-                        const double D = rootD[dIndex];
-                        const double E = rootE[eIndex];
-                        const double b = A * (B + C + D + E);
-                        const double c = A * (B * C + C * D + B * D + B * E + C * E + D * E);
-                        const double d = A * (B * C * D + B * C * E + B * D * E + C * D * E);
-                        const double e = A * B * C * D * E;
-                        double roots[4];
-                        const int rootCount = quarticRoots(A, b, c, d, e, roots);
-                        const int expected = 1 + (B != C) + (B != D && C != D) + (B != E && C != E && D != E);
-                        assert(rootCount == expected);
-                        assert(approximately_equal(roots[0], -B)
-                                || approximately_equal(roots[0], -C)
-                                || approximately_equal(roots[0], -D)
-                                || approximately_equal(roots[0], -E));
-                        if (expected > 1) {
-                            assert(!approximately_equal(roots[0], roots[1]));
-                            assert(approximately_equal(roots[1], -B)
-                                    || approximately_equal(roots[1], -C)
-                                    || approximately_equal(roots[1], -D)
-                                    || approximately_equal(roots[1], -E));
-                            if (expected > 2) {
-                                assert(!approximately_equal(roots[0], roots[2])
-                                        && !approximately_equal(roots[1], roots[2]));
-                                assert(approximately_equal(roots[2], -B)
-                                        || approximately_equal(roots[2], -C)
-                                        || approximately_equal(roots[2], -D)
-                                        || approximately_equal(roots[2], -E));
-                                if (expected > 3) {
-                                    assert(!approximately_equal(roots[0], roots[3])
-                                            && !approximately_equal(roots[1], roots[3])
-                                            && !approximately_equal(roots[2], roots[3]));
-                                    assert(approximately_equal(roots[3], -B)
-                                            || approximately_equal(roots[3], -C)
-                                            || approximately_equal(roots[3], -D)
-                                            || approximately_equal(roots[3], -E));
-                                }
-                            }
-                        }
+                        testOneQuartic(aIndex, bIndex, cIndex, dIndex, eIndex);
                     }
                 }
             }
@@ -139,7 +187,11 @@ static void quarticTest() {
 }
 
 void QuarticRoot_Test() {
-    quadraticTest();
-    cubicTest();
+    testOneCubic(false, 0, 5, 5, 4);
+    testOneQuartic(0, 0, 2, 4, 3);
+    quadraticTest(true);
+    quadraticTest(false);
+    cubicTest(true);
+    cubicTest(false);
     quarticTest();
 }
diff --git a/experimental/Intersection/TestUtilities.cpp b/experimental/Intersection/TestUtilities.cpp
index 0477ed48e1..3ae10d216a 100644
--- a/experimental/Intersection/TestUtilities.cpp
+++ b/experimental/Intersection/TestUtilities.cpp
@@ -63,4 +63,3 @@ bool controls_inside(const Cubic& cubic) {
      && ((cubic[0].y <= cubic[1].y && cubic[0].y <= cubic[2].y && cubic[1].y <= cubic[3].y && cubic[2].y <= cubic[3].y)
      ||  (cubic[0].y >= cubic[1].y && cubic[0].y >= cubic[2].y && cubic[1].y >= cubic[3].y && cubic[2].x >= cubic[3].y));
 }
-
diff --git a/experimental/Intersection/qc.htm b/experimental/Intersection/qc.htm
index 3fe3998b0c..2cf3592e56 100644
--- a/experimental/Intersection/qc.htm
+++ b/experimental/Intersection/qc.htm
@@ -1545,11 +1545,119 @@ $7 = {{x = 24.006224853920855, y = 72.621119847810419}, {x = 29.758671200376888,
 {{61.3365082,82.6931328}, {54.8250789,71.6639328}, {47.7274442,61.4049645}},
 </div>
 
+<div id="quad15">
+{{x = 80.897794748143198, y = 49.236332042718459}, {x = 81.082078218891212, y = 64.066749904488631}, {x = 69.972305057149981, y = 72.968595519850993}}
+{{x = 72.503745601281395, y = 32.952320736577882}, {x = 88.030880716061645, y = 38.137194847810164}, {x = 73.193774825517906, y = 67.773492479591397}}
+</div>
+
+<div id="cubic19">
+{{x = 34.560092601254624, y = 51.476349286491221}, {x = 27.498466254909744, y = 66.722346267999313}, {x = 42.500359724508769, y = 3.5458898188294325}, {x = 73.37353619438295, y = 89.022818994253328}}
+{{x = 63.002458057833124, y = 82.312578001205154}, {x = 2.4737262644217006, y = 75.917326135522373}, {x = 95.77018506628005, y = 9.5004089686555826}, {x = 6.5188364156143912, y = 62.083637231068508}}
+</div>
+
+<div id="cubic20">
+ {{x = 42.449716172390481, y = 52.379709366885805}, {x = 27.896043159019225, y = 48.797373636065686}, {x = 92.770268299044233, y = 89.899302036454571}, {x = 12.102066544863426, y = 99.43241951960718}}
+{{x = 45.77532924980639, y = 45.958701495993274}, {x = 37.458701356062065, y = 68.393691335056758}, {x = 37.569326692060258, y = 27.673713456687381}, {x = 60.674866037757539, y = 62.47349659096146}}
+</div>
+
+<div id="cubic21">
+{{x = 26.192053931854691, y = 9.8504326817814416}, {x = 10.174241480498686, y = 98.476562741434464}, {x = 21.177712558385782, y = 33.814968789841501}, {x = 75.329030899018534, y = 55.02231980442177}}
+{{x = 56.222082700683771, y = 24.54395039218662}, {x = 95.589995289030483, y = 81.050822735322086}, {x = 28.180450866082897, y = 28.837706255185282}, {x = 60.128952916771617, y = 87.311672180570511}}
+</div>
+
+<div id="quad16">
+{{x = 67.965974918365831, y = 52.573040929556633}, {x = 67.973015821010591, y = 52.57495862082331}, {x = 67.980057838863502, y = 52.576878275262274}}
+{{x = 67.975025709349239, y = 52.572750461020817}, {x = 67.973101328974863, y = 52.57506284863603}, {x = 67.971173663444745, y = 52.577372136133093}}
+</div>
+
+<div id="quad17">
+{{x = 52.14807018377202, y = 65.012420045148644}, {x = 44.778669050208237, y = 66.315562705604378}, {x = 51.619118408823567, y = 63.787827046262684}}
+{{x = 30.004993234763383, y = 93.921296668202288}, {x = 53.384822003076991, y = 60.732180341802753}, {x = 58.652998934338584, y = 43.111073088306185}}
+</div>
+
+<div id="quad18">
+{{x = 369.850525, y = 145.67596399999999}, {x = 382.36291499999999, y = 121.29286999999999}, {x = 406.21127300000001, y = 121.29286999999999}}
+{{x = 369.962311, y = 137.976044}, {x = 383.97189300000002, y = 121.29286999999999}, {x = 406.21612499999998, y = 121.29286999999999}}
+</div>
+
+<div id="quad19">
+{{x = 406.23635899999999, y = 121.254936}, {x = 409.44567899999998, y = 121.254936}, {x = 412.97595200000001, y = 121.789818}}
+{{x = 406.23599200000001, y = 121.254936}, {x = 425.70590199999998, y = 121.254936}, {x = 439.71994000000001, y = 137.087616}}
+</div>
+
+<div id="cubic22">
+{{x = 7.5374809128872498, y = 82.441702896003477}, {x = 22.444346930107265, y = 22.138854312775123}, {x = 66.76091829629658, y = 50.753805856571446}, {x = 78.193478508942519, y = 97.7932997968948}}
+{{x = 97.700573130371311, y = 53.53260215070685}, {x = 87.72443481149358, y = 84.575876772671876}, {x = 19.215031396232092, y = 47.032676472809484}, {x = 11.989686410869325, y = 10.659507480757082}}
+  {{7.53748091,82.4417029}, {15.5677076,52.942994}, {29.9404074,49.1672596}},
+  {{29.9404074,49.1672596}, {44.3131071,45.3915253}, {58.1067559,59.5061814}},
+  {{58.1067559,59.5061814}, {71.9004047,73.6208375}, {78.1934785,97.7932998}},
+
+  {{97.7005731,53.5326022}, {91.6030843,68.4083459}, {72.6510251,64.2972928}},
+  {{72.6510251,64.2972928}, {53.6989659,60.1862397}, {35.2053722,44.8391126}},
+  {{35.2053722,44.8391126}, {16.7117786,29.4919856}, {11.9896864,10.6595075}},
+</div>
+
+<div id="quad20">
+  {{29.9404074,49.1672596}, {44.3131071,45.3915253}, {58.1067559,59.5061814}},
+  {{72.6510251,64.2972928}, {53.6989659,60.1862397}, {35.2053722,44.8391126}},
+</div>
+
+<div id="cubic23">
+{{x = 32.484981432782945, y = 75.082940782924624}, {x = 42.467313093350882, y = 48.131159948246157}, {x = 3.5963115764764657, y = 43.208665839959245}, {x = 79.442476890721579, y = 89.709102357602262}}
+{{x = 18.98573861410177, y = 93.308887208490106}, {x = 40.405250173250792, y = 91.039661826118675}, {x = 8.0467721950480584, y = 42.100282172719147}, {x = 40.883324221187891, y = 26.030185504830527}}
+  {{32.4849814,75.0829408}, {35.4553509,65.5763004}, {33.5767697,60.2097835}},
+  {{33.5767697,60.2097835}, {31.6981886,54.8432666}, {31.1663962,54.7302484}},
+  {{31.1663962,54.7302484}, {31.1662882,54.7301074}, {31.1663969,54.7302485}},
+  {{31.1663969,54.7302485}, {30.4117445,54.6146017}, {40.1631726,62.9428436}},
+  {{40.1631726,62.9428436}, {49.9146008,71.2710854}, {79.4424769,89.7091024}},
+
+  {{18.9857386,93.3088872}, {25.7662938,92.3417699}, {26.5917262,85.8225583}},
+  {{26.5917262,85.8225583}, {27.4171586,79.3033467}, {26.141946,69.8089528}},
+  {{26.141946,69.8089528}, {24.2922348,57.665767}, {26.0404936,45.4260361}},
+  {{26.0404936,45.4260361}, {27.7887523,33.1863051}, {40.8833242,26.0301855}},
+</div>
+
+<div id="quad21">
+  {{31.1663962,54.7302484}, {31.1662882,54.7301074}, {31.1663969,54.7302485}},
+  {{26.0404936,45.4260361}, {27.7887523,33.1863051}, {40.8833242,26.0301855}},
+</div>
+
+<div id="cubic24">
+{{x = 65.454505973241524, y = 93.881892270353575}, {x = 45.867360264932437, y = 92.723972719499827}, {x = 2.1464054482739447, y = 74.636369140183717}, {x = 33.774068594804994, y = 40.770872887582925}}
+{{x = 72.963387832494163, y = 95.659300729473728}, {x = 11.809496633619768, y = 82.209921247423594}, {x = 13.456139067865974, y = 57.329313623406605}, {x = 36.060621606214262, y = 70.867335643091849}}
+  {{65.454506,93.8818923}, {54.7397995,93.2922678}, {41.5072916,87.1234036}},
+  {{41.5072916,87.1234036}, {28.2747836,80.9545395}, {23.5780771,69.3344126}},
+  {{23.5780771,69.3344126}, {18.8813706,57.7142857}, {33.7740686,40.7708729}},
+
+  {{72.9633878,95.6593007}, {42.7738746,88.4730382}, {31.1932785,80.2458029}},
+  {{31.1932785,80.2458029}, {19.6126823,72.0185676}, {21.9918152,68.2892325}},
+  {{21.9918152,68.2892325}, {24.370948,64.5598974}, {36.0606216,70.8673356}},
+</div>
+
+<div id="quad22">
+  {{41.5072916,87.1234036}, {28.2747836,80.9545395}, {23.5780771,69.3344126}},
+  {{72.9633878,95.6593007}, {42.7738746,88.4730382}, {31.1932785,80.2458029}},
+</div>
+
 </div>
 
 <script type="text/javascript">
 
 var testDivs = [
+    quad22,
+    cubic24,
+    quad21,
+    cubic23,
+    quad20,
+    cubic22,
+    quad19,
+    quad18,
+    quad17,
+    quad16,
+    cubic21,
+    cubic20,
+    cubic19,
+    quad15,
     quad14,
     cubic18,
     cubic17,