diff --git a/experimental/Intersection/BezierClip.cpp b/experimental/Intersection/BezierClip.cpp
new file mode 100644
index 0000000..1b5c60c
--- /dev/null
+++ b/experimental/Intersection/BezierClip.cpp
@@ -0,0 +1,83 @@
+#include "CubicIntersection.h"
+#include "LineParameters.h"
+#include <algorithm> // used for std::swap
+
+// return false if unable to clip (e.g., unable to create implicit line)
+// caller should subdivide, or create degenerate if the values are too small
+bool bezier_clip(const Cubic& cubic1, const Cubic& cubic2, double& minT, double& maxT) {
+    minT = 1;
+    maxT = 0;
+    // determine normalized implicit line equation for pt[0] to pt[3]
+    //   of the form ax + by + c = 0, where a*a + b*b == 1
+    
+    // find the implicit line equation parameters
+    LineParameters endLine;
+    endLine.cubicEndPoints(cubic1);
+    if (!endLine.normalize()) {
+        printf("line cannot be normalized: need more code here\n");
+        return false;
+    }
+
+    double distance[2];
+    endLine.controlPtDistance(cubic1, distance);
+    
+    // find fat line
+    double top = distance[0];
+    double bottom = distance[1];
+    if (top > bottom) {
+        std::swap(top, bottom);
+    }
+    if (top * bottom >= 0) {
+        const double scale = 3/4.0; // http://cagd.cs.byu.edu/~tom/papers/bezclip.pdf (13)
+        if (top < 0) {
+            top *= scale;
+            bottom = 0;
+        } else {
+            top = 0;
+            bottom *= scale;
+        }
+    } else {
+        const double scale = 4/9.0; // http://cagd.cs.byu.edu/~tom/papers/bezclip.pdf (15)
+        top *= scale;
+        bottom *= scale;
+    }
+    
+    // compute intersecting candidate distance
+    Cubic distance2y; // points with X of (0, 1/3, 2/3, 1)
+    endLine.cubicDistanceY(cubic2, distance2y);
+    
+    int flags = 0;
+    if (approximately_lesser(distance2y[0].y, top)) {
+        flags |= kFindTopMin;
+    } else if (approximately_greater(distance2y[0].y, bottom)) {
+        flags |= kFindBottomMin;
+    } else {
+        minT = 0;
+    }
+
+    if (approximately_lesser(distance2y[3].y, top)) {
+        flags |= kFindTopMax;
+    } else if (approximately_greater(distance2y[3].y, bottom)) {
+        flags |= kFindBottomMax;
+    } else {
+        maxT = 1;
+    }
+    // Find the intersection of distance convex hull and fat line.
+    char to_0[2];
+    char to_3[2];
+    bool do_1_2_edge = convex_x_hull(distance2y, to_0, to_3);
+    x_at(distance2y[0], distance2y[to_0[0]], top, bottom, flags, minT, maxT);
+    if (to_0[0] != to_0[1]) {
+        x_at(distance2y[0], distance2y[to_0[1]], top, bottom, flags, minT, maxT);
+    }
+    x_at(distance2y[to_3[0]], distance2y[3], top, bottom, flags, minT, maxT);
+    if (to_3[0] != to_3[1]) {
+        x_at(distance2y[to_3[1]], distance2y[3], top, bottom, flags, minT, maxT);
+    }
+    if (do_1_2_edge) {
+        x_at(distance2y[1], distance2y[2], top, bottom, flags, minT, maxT);
+    }
+    
+    return minT < maxT; // returns false if distance shows no intersection
+}
+
diff --git a/experimental/Intersection/BezierClip_Test.cpp b/experimental/Intersection/BezierClip_Test.cpp
new file mode 100644
index 0000000..3228d89
--- /dev/null
+++ b/experimental/Intersection/BezierClip_Test.cpp
@@ -0,0 +1,24 @@
+#include "CubicIntersection.h"
+#include "CubicIntersection_TestData.h"
+#include "CubicIntersection_Tests.h"
+
+void BezierClip_Test() {
+    for (size_t index = 0; index < tests_count; ++index) {
+        const Cubic& cubic1 = tests[index][0];
+        const Cubic& cubic2 = tests[index][1];
+        Cubic reduce1, reduce2;
+        int order1 = reduceOrder(cubic1, reduce1, kReduceOrder_NoQuadraticsAllowed);
+        int order2 = reduceOrder(cubic2, reduce2, kReduceOrder_NoQuadraticsAllowed);
+        if (order1 < 4) {
+            printf("%s [%d] cubic1 order=%d\n", __FUNCTION__, (int) index, order1);
+        }
+        if (order2 < 4) {
+            printf("%s [%d] cubic2 order=%d\n", __FUNCTION__, (int) index, order2);
+        }
+        if (order1 == 4 && order2 == 4) {
+            double minT = 0;
+            double maxT = 1;
+            bezier_clip(reduce1, reduce2, minT, maxT);
+        }
+    }
+}
diff --git a/experimental/Intersection/ConvexHull.cpp b/experimental/Intersection/ConvexHull.cpp
new file mode 100644
index 0000000..b03eb46
--- /dev/null
+++ b/experimental/Intersection/ConvexHull.cpp
@@ -0,0 +1,137 @@
+#include "CubicIntersection.h"
+#include "IntersectionUtilities.h"
+
+/* Given a cubic, find the convex hull described by the end and control points.
+   The hull may have 3 or 4 points. Cubics that degenerate into a point or line
+   are not considered.
+  
+   The hull is computed by assuming that three points, if unique and non-linear,
+   form a triangle. The fourth point may replace one of the first three, may be
+   discarded if in the triangle or on an edge, or may be inserted between any of
+   the three to form a convex quadralateral.
+   
+   The indices returned in order describe the convex hull.
+*/
+int convex_hull(const Cubic& cubic, char order[4]) {
+    size_t index;
+    // find top point
+    size_t yMin = 0;
+    for (index = 1; index < 4; ++index) {
+        if (cubic[yMin].y > cubic[index].y || (cubic[yMin].y == cubic[index].y
+                && cubic[yMin].x > cubic[index].x)) {
+            yMin = index;
+        }
+    }
+    order[0] = yMin;
+    int midX = -1;
+    int backupYMin = -1;
+    for (int pass = 0; pass < 2; ++pass) {
+        for (index = 0; index < 4; ++index) {
+            if (index == yMin) {
+                continue;
+            }
+            // rotate line from (yMin, index) to axis
+            // see if remaining two points are both above or below
+            // use this to find mid
+            int mask = other_two(yMin, index);
+            int side1 = yMin ^ mask;
+            int side2 = index ^ mask;
+            Cubic rotPath;
+            if (!rotate(cubic, yMin, index, rotPath)) { // ! if cbc[yMin]==cbc[idx]
+                order[1] = side1;
+                order[2] = side2;
+                return 3;
+            }
+            int sides = side(rotPath[side1].y - rotPath[yMin].y);
+            sides ^= side(rotPath[side2].y - rotPath[yMin].y);
+            if (sides == 2) { // '2' means one remaining point <0, one >0
+                if (midX >= 0) {
+                    printf("%s unexpected mid\n", __FUNCTION__); // there can be only one mid
+                }
+                midX = index;
+            } else if (sides == 0) { // '0' means both to one side or the other
+                backupYMin = index;
+            }
+        }
+        if (midX >= 0) {
+            break;
+        }
+        if (backupYMin < 0) {
+            break;
+        }
+        yMin = backupYMin;
+        backupYMin = -1;
+    }
+    if (midX < 0) {
+        midX = yMin ^ 3; // choose any other point
+    }
+    int mask = other_two(yMin, midX);
+    int least = yMin ^ mask;
+    int most = midX ^ mask;
+    order[0] = yMin;
+    order[1] = least;
+    
+    // see if mid value is on same side of line (least, most) as yMin
+    Cubic midPath;
+    if (!rotate(cubic, least, most, midPath)) { // ! if cbc[least]==cbc[most]
+        order[2] = midX;
+        return 3;
+    }
+    int midSides = side(midPath[yMin].y - midPath[least].y);
+    midSides ^= side(midPath[midX].y - midPath[least].y);
+    if (midSides != 2) {  // if mid point is not between 
+        order[2] = most;
+        return 3; // result is a triangle
+    }
+    order[2] = midX;
+    order[3] = most;
+    return 4; // result is a quadralateral
+}
+
+/* Find the convex hull for cubics with the x-axis interval regularly spaced.
+   Cubics computed as distance functions are formed this way.
+   
+   connectTo0[0], connectTo0[1] are the point indices that cubic[0] connects to.
+   connectTo3[0], connectTo3[1] are the point indices that cubic[3] connects to.
+   
+   Returns true if cubic[1] to cubic[2] also forms part of the hull.
+*/
+bool convex_x_hull(const Cubic& cubic, char connectTo0[2], char connectTo3[2]) {
+    double projectedY[4];
+    projectedY[0] = 0;
+    int index;
+    for (index = 1; index < 4; ++index) {
+        projectedY[index] = (cubic[index].y - cubic[0].y) * (3.0 / index);
+    }
+    int lower0Index = 1;
+    int upper0Index = 1;
+    for (index = 2; index < 4; ++index) {
+        if (approximately_greater(projectedY[lower0Index], projectedY[index])) {
+            lower0Index = index;
+        }
+        if (approximately_lesser(projectedY[upper0Index], projectedY[index])) {
+            upper0Index = index;
+        }
+    }
+    connectTo0[0] = lower0Index;
+    connectTo0[1] = upper0Index;
+    for (index = 0; index < 3; ++index) {
+        projectedY[index] = (cubic[3].y - cubic[index].y) * (3.0 / (3 - index));
+    }
+    projectedY[3] = 0;
+    int lower3Index = 2;
+    int upper3Index = 2;
+    for (index = 1; index > -1; --index) {
+        if (approximately_greater(projectedY[lower3Index], projectedY[index])) {
+            lower3Index = index;
+        }
+        if (approximately_lesser(projectedY[upper3Index], projectedY[index])) {
+            upper3Index = index;
+        }
+    }
+    connectTo3[0] = lower3Index;
+    connectTo3[1] = upper3Index;
+    return (1 << lower0Index | 1 << upper0Index
+            | 1 << lower3Index | 1 << upper3Index) == 0x0F;
+}
+
diff --git a/experimental/Intersection/ConvexHull_Test.cpp b/experimental/Intersection/ConvexHull_Test.cpp
new file mode 100644
index 0000000..dffca6d
--- /dev/null
+++ b/experimental/Intersection/ConvexHull_Test.cpp
@@ -0,0 +1,465 @@
+#include "CubicIntersection.h"
+#include "CubicIntersection_Tests.h"
+#include "IntersectionUtilities.h"
+
+const Cubic convex[] = {
+    {{0, 0}, {2, 0}, {2, 1}, {0, 1}},
+    {{1, 0}, {1, 1}, {0, 1}, {0, 0}},
+    {{1, 1}, {0, 1}, {0, 0}, {1, 0}},
+    {{0, 1}, {0, 0}, {1, 0}, {1, 1}},
+    {{0, 0}, {10, 0}, {10, 10}, {5, 6}},
+};
+
+size_t convex_count = sizeof(convex) / sizeof(convex[0]);
+
+const Cubic bowtie[] = {
+    {{0, 0}, {1, 1}, {1, 0}, {0, 1}},
+    {{1, 0}, {0, 1}, {1, 1}, {0, 0}},
+    {{1, 1}, {0, 0}, {0, 1}, {1, 0}},
+    {{0, 1}, {1, 0}, {0, 0}, {1, 1}},
+};
+
+size_t bowtie_count = sizeof(bowtie) / sizeof(bowtie[0]);
+
+const Cubic arrow[] = {
+    {{0, 0}, {10, 0}, {10, 10}, {5, 4}},
+    {{10, 0}, {10, 10}, {5, 4}, {0, 0}},
+    {{10, 10}, {5, 4}, {0, 0}, {10, 0}},
+    {{5, 4}, {0, 0}, {10, 0}, {10, 10}},
+};
+
+size_t arrow_count = sizeof(arrow) / sizeof(arrow[0]);
+
+const Cubic three[] = {
+    {{1, 0}, {1, 0}, {1, 1}, {0, 1}}, // 0 == 1
+    {{0, 0}, {1, 1}, {1, 1}, {0, 1}}, // 1 == 2
+    {{0, 0}, {1, 0}, {0, 1}, {0, 1}}, // 2 == 3
+    {{1, 0}, {1, 1}, {1, 0}, {0, 1}}, // 0 == 2
+    {{1, 0}, {1, 1}, {0, 1}, {1, 0}}, // 0 == 3
+    {{0, 0}, {1, 0}, {1, 1}, {1, 0}}, // 1 == 3
+};
+
+size_t three_count = sizeof(three) / sizeof(three[0]);
+
+const Cubic triangle[] = {
+    {{0, 0}, {1, 0}, {2, 0}, {0, 1}}, // extra point on horz
+    {{1, 0}, {2, 0}, {0, 1}, {0, 0}},
+    {{2, 0}, {0, 1}, {0, 0}, {1, 0}},
+    {{0, 1}, {0, 0}, {1, 0}, {2, 0}},
+    
+    {{0, 0}, {0, 1}, {0, 2}, {1, 1}}, // extra point on vert
+    {{0, 1}, {0, 2}, {1, 1}, {0, 0}},
+    {{0, 2}, {1, 1}, {0, 0}, {0, 1}},
+    {{1, 1}, {0, 0}, {0, 1}, {0, 2}},
+    
+    {{0, 0}, {1, 1}, {2, 2}, {2, 0}}, // extra point on diag
+    {{1, 1}, {2, 2}, {2, 0}, {0, 0}},
+    {{2, 2}, {2, 0}, {0, 0}, {1, 1}},
+    {{2, 0}, {0, 0}, {1, 1}, {2, 2}},
+    
+    {{0, 0}, {2, 0}, {2, 2}, {1, 1}}, // extra point on diag
+    {{2, 0}, {2, 2}, {1, 1}, {0, 0}},
+    {{2, 2}, {1, 1}, {0, 0}, {2, 0}},
+    {{1, 1}, {0, 0}, {2, 0}, {2, 2}},
+};
+
+size_t triangle_count = sizeof(triangle) / sizeof(triangle[0]);
+
+const struct CubicDataSet {
+    const Cubic* data;
+    size_t size;
+} cubicDataSet[] = {
+    { three, three_count },
+    { convex, convex_count },
+    { bowtie, bowtie_count },
+    { arrow, arrow_count },
+    { triangle, triangle_count },
+};
+
+size_t cubicDataSet_count = sizeof(cubicDataSet) / sizeof(cubicDataSet[0]);
+
+typedef double Matrix3x2[3][2];
+
+static bool rotateToAxis(const _Point& a, const _Point& b, Matrix3x2& matrix) {
+    double dx = b.x - a.x;
+    double dy = b.y - a.y;
+    double length = sqrt(dx * dx + dy * dy);
+    if (length == 0) {
+        return false;
+    }
+    double invLength = 1 / length;
+    matrix[0][0] = dx * invLength;
+    matrix[1][0] = dy * invLength;
+    matrix[2][0] = 0;
+    matrix[0][1] = -dy * invLength;
+    matrix[1][1] = dx * invLength;
+    matrix[2][1] = 0;
+    return true;
+}
+
+static void transform(const Cubic& cubic, const Matrix3x2& matrix, Cubic& rotPath) {
+    for (int index = 0; index < 4; ++index) {
+        rotPath[index].x = cubic[index].x * matrix[0][0] 
+                + cubic[index].y * matrix[1][0] + matrix[2][0];
+        rotPath[index].y = cubic[index].x * matrix[0][1] 
+                + cubic[index].y * matrix[1][1] + matrix[2][1];
+    }
+}
+
+// brute force way to find convex hull:
+// pick two points
+// rotate all four until the two points are horizontal
+// are the remaining two points both above or below the horizontal line?
+// if so, the two points must be an edge of the convex hull
+static int rotate_to_hull(const Cubic& cubic, char order[4], size_t idx, size_t inr) {
+    bool debug_rotate_to_hull = false;
+    int outsidePtSet[4];
+    memset(outsidePtSet, -1, sizeof(outsidePtSet));
+    for (int outer = 0; outer < 3; ++outer) {
+        for (int priorOuter = 0; priorOuter < outer; ++priorOuter) {
+            if (cubic[outer].approximatelyEqual(cubic[priorOuter])) {
+                goto skip;
+            }
+        }
+        for (int inner = outer + 1; inner < 4; ++inner) {
+            for (int priorInner = outer + 1; priorInner < inner; ++priorInner) {
+                if (cubic[inner].approximatelyEqual(cubic[priorInner])) {
+                    goto skipInner;
+                }
+            }
+            if (cubic[outer].approximatelyEqual(cubic[inner])) {
+                continue;
+            }
+            Matrix3x2 matrix;
+            if (!rotateToAxis(cubic[outer], cubic[inner], matrix)) {
+                continue;
+            }
+            Cubic rotPath;
+            transform(cubic, matrix, rotPath);
+            int sides[3];
+            int zeroes;
+            zeroes = -1;
+            bzero(sides, sizeof(sides));
+            if (debug_rotate_to_hull) printf("%s [%d,%d] [o=%d,i=%d] src=(%g,%g) rot=", __FUNCTION__,
+                    (int)idx, (int)inr, (int)outer, (int)inner,
+                    cubic[inner].x, cubic[inner].y);
+            for (int index = 0; index < 4; ++index) {
+                if (debug_rotate_to_hull) printf("(%g,%g) ", rotPath[index].x, rotPath[index].y);
+                sides[side(rotPath[index].y - rotPath[inner].y)]++;
+                if (index != outer && index != inner 
+                        && side(rotPath[index].y - rotPath[inner].y) == 1)
+                    zeroes = index;
+            }
+            if (debug_rotate_to_hull) printf("sides=(%d,%d,%d)\n", sides[0], sides[1], sides[2]);
+            if (sides[0] && sides[2]) {
+                continue;
+            }
+            if (sides[1] == 3 && zeroes >= 0) {
+                // verify that third point is between outer, inner
+                // if either of remaining two equals outer or equal, pick lower
+                if (rotPath[zeroes].approximatelyEqual(rotPath[inner])
+                        && zeroes < inner) {
+                    if (debug_rotate_to_hull) printf("%s [%d,%d] [o=%d,i=%d] zeroes < inner\n",
+                        __FUNCTION__, (int)idx, (int)inr, (int)outer, (int)inner);
+                    continue;
+                }
+                 if (rotPath[zeroes].approximatelyEqual(rotPath[outer])
+                        && zeroes < outer) {
+                    if (debug_rotate_to_hull) printf("%s [%d,%d] [o=%d,i=%d] zeroes < outer\n",
+                        __FUNCTION__, (int)idx, (int)inr, (int)outer, (int)inner);
+                    continue;
+                }
+                if (rotPath[zeroes].x < rotPath[inner].x 
+                        && rotPath[zeroes].x < rotPath[outer].x) {
+                    if (debug_rotate_to_hull) printf("%s [%d,%d] [o=%d,i=%d] zeroes < inner && outer\n",
+                        __FUNCTION__, (int)idx, (int)inr, (int)outer, (int)inner);
+                    continue;
+                }
+                if (rotPath[zeroes].x > rotPath[inner].x 
+                        && rotPath[zeroes].x > rotPath[outer].x) {
+                    if (debug_rotate_to_hull) printf("%s [%d,%d] [o=%d,i=%d] zeroes > inner && outer\n",
+                        __FUNCTION__, (int)idx, (int)inr, (int)outer, (int)inner);
+                    continue;
+                }
+            }
+            if (outsidePtSet[outer] < 0) {
+                outsidePtSet[outer] = inner;
+            } else {
+                if (outsidePtSet[inner] > 0) {
+                    if (debug_rotate_to_hull) printf("%s [%d,%d] [o=%d,i=%d] too many rays from one point\n",
+                        __FUNCTION__, (int)idx, (int)inr, (int)outer, (int)inner);
+                }
+                outsidePtSet[inner] = outer;
+            }
+skipInner:
+            ;
+        }
+skip:
+        ;
+    }
+    int totalSides = 0;
+    int first = 0;
+    for (; first < 4; ++first) {
+        if (outsidePtSet[first] >= 0) {
+            break;
+        }
+    }
+    if (first > 3) {
+        order[0] = 0;
+        return 1;
+    }
+    int next = first;
+    do {
+        order[totalSides++] = next;
+        next = outsidePtSet[next];
+    } while (next != -1 && next != first);
+    return totalSides;
+}
+
+int firstIndex = 0;
+int firstInner = 0;
+
+void ConvexHull_Test() {
+    for (size_t index = firstIndex; index < cubicDataSet_count; ++index) {
+        const CubicDataSet& set = cubicDataSet[index];
+        for (size_t inner = firstInner; inner < set.size; ++inner) {
+            const Cubic& cubic = set.data[inner];
+            char order[4], cmpOrder[4];
+            int cmp = rotate_to_hull(cubic, cmpOrder, index, inner);
+            if (cmp < 3) {
+                continue;
+            }
+            int result = convex_hull(cubic, order);
+            if (cmp != result) {
+                printf("%s [%d,%d] result=%d cmp=%d\n", __FUNCTION__,
+                    (int)index, (int)inner, result, cmp);
+                continue;
+            }
+            // check for same indices
+            char pts = 0;
+            char cmpPts = 0;
+            int pt, bit;
+            for (pt = 0; pt < cmp; ++pt) {
+                if (pts & 1 << order[pt]) {
+                    printf("%s [%d,%d] duplicate index in order: %d,%d,%d",
+                            __FUNCTION__, (int)index, (int)inner, 
+                            order[0], order[1], order[2]);
+                    if (cmp == 4) {
+                        printf(",%d", order[3]);
+                    }
+                    printf("\n");
+                    goto next;
+                }
+                if (cmpPts & 1 << cmpOrder[pt]) {
+                    printf("%s [%d,%d] duplicate index in order: %d,%d,%d",
+                            __FUNCTION__, (int)index, (int)inner, 
+                            cmpOrder[0], cmpOrder[1], cmpOrder[2]);
+                    if (cmp == 4) {
+                        printf(",%d", cmpOrder[3]);
+                    }
+                    printf("\n");
+                    goto next;
+                }
+                pts |= 1 << order[pt];
+                cmpPts |= 1 << cmpOrder[pt];
+            }
+            for (bit = 0; bit < 4; ++bit) {
+                if (pts & 1 << bit) {
+                    continue;
+                }
+                for (pt = 0; pt < cmp; ++pt) {
+                    if (order[pt] == bit) {
+                        continue;
+                    }
+                    if (cubic[order[pt]] == cubic[bit]) {
+                        pts |= 1 << bit;
+                    }
+                }
+            }
+            for (bit = 0; bit < 4; ++bit) {
+                if (cmpPts & 1 << bit) {
+                    continue;
+                }
+                for (pt = 0; pt < cmp; ++pt) {
+                    if (cmpOrder[pt] == bit) {
+                        continue;
+                    }
+                    if (cubic[cmpOrder[pt]] == cubic[bit]) {
+                        cmpPts |= 1 << bit;
+                    }
+                }
+            }
+            if (pts != cmpPts) {
+                printf("%s [%d,%d] mismatch indices: order=%d,%d,%d",
+                        __FUNCTION__, (int)index, (int)inner, 
+                        order[0], order[1], order[2]);
+                if (cmp == 4) {
+                    printf(",%d", order[3]);
+                }
+                printf(" cmpOrder=%d,%d,%d", cmpOrder[0], cmpOrder[1], cmpOrder[2]);
+                if (cmp == 4) {
+                    printf(",%d", cmpOrder[3]);
+                }
+                printf("\n");
+                continue;
+            }
+            if (cmp == 4) { // check for bow ties
+                int match = 0;
+                while (cmpOrder[match] != order[0]) {
+                    ++match;
+                }
+                if (cmpOrder[match ^ 2] != order[2]) {
+                    printf("%s [%d,%d] bowtie mismatch: order=%d,%d,%d,%d"
+                            " cmpOrder=%d,%d,%d,%d\n",
+                            __FUNCTION__, (int)index, (int)inner, 
+                            order[0], order[1], order[2], order[3],
+                            cmpOrder[0], cmpOrder[1], cmpOrder[2], cmpOrder[3]);
+                }
+            }
+    next:
+            ;
+        }
+    }
+}
+
+const double a = 1.0/3;
+const double b = 2.0/3;
+
+const Cubic x_cubic[] = {
+    {{0, 0}, {a, 0}, {b, 0}, {1, 0}}, // 0
+    {{0, 0}, {a, 0}, {b, 0}, {1, 1}}, // 1
+    {{0, 0}, {a, 0}, {b, 1}, {1, 0}}, // 2
+    {{0, 0}, {a, 0}, {b, 1}, {1, 1}}, // 3
+    {{0, 0}, {a, 1}, {b, 0}, {1, 0}}, // 4
+    {{0, 0}, {a, 1}, {b, 0}, {1, 1}}, // 5
+    {{0, 0}, {a, 1}, {b, 1}, {1, 0}}, // 6
+    {{0, 0}, {a, 1}, {b, 1}, {1, 1}}, // 7
+    {{0, 1}, {a, 0}, {b, 0}, {1, 0}}, // 8
+    {{0, 1}, {a, 0}, {b, 0}, {1, 1}}, // 9
+    {{0, 1}, {a, 0}, {b, 1}, {1, 0}}, // 10
+    {{0, 1}, {a, 0}, {b, 1}, {1, 1}}, // 11
+    {{0, 1}, {a, 1}, {b, 0}, {1, 0}}, // 12
+    {{0, 1}, {a, 1}, {b, 0}, {1, 1}}, // 13
+    {{0, 1}, {a, 1}, {b, 1}, {1, 0}}, // 14
+    {{0, 1}, {a, 1}, {b, 1}, {1, 1}}, // 15
+};
+
+size_t x_cubic_count = sizeof(x_cubic) / sizeof(x_cubic[0]);
+
+static int first_x_test = 0;
+
+void ConvexHull_X_Test() {
+    for (size_t index = first_x_test; index < x_cubic_count; ++index) {
+        const Cubic& cubic = x_cubic[index];
+        char connectTo0[2] = {-1, -1};
+        char connectTo3[2] = {-1, -1};
+        convex_x_hull(cubic, connectTo0, connectTo3);
+        int idx, cmp;
+        for (idx = 0; idx < 2; ++idx) {
+            if (connectTo0[idx] >= 1 && connectTo0[idx] < 4) {
+                continue;
+            } else {
+                printf("%s connectTo0[idx]=%d", __FUNCTION__, connectTo0[idx]);
+            }
+            if (connectTo3[idx] >= 0 && connectTo3[idx] < 3) {
+                continue;
+            } else {
+                printf("%s connectTo3[idx]=%d", __FUNCTION__, connectTo3[idx]);
+            }
+            goto nextTest;
+        }
+        char rOrder[4];
+        char cmpOrder[4];
+        cmp = rotate_to_hull(cubic, cmpOrder, index, 0);
+        if (index == 0 || index == 15) {
+            // FIXME: make rotate_to_hull work for degenerate 2 edge hull cases
+            cmpOrder[0] = 0;
+            cmpOrder[1] = 3;
+            cmp = 2;
+        }
+        if (cmp < 3) {
+            // FIXME: make rotate_to_hull work for index == 3 etc
+            continue;
+        }
+        for (idx = 0; idx < cmp; ++idx) {
+            if (cmpOrder[idx] == 0) {
+                rOrder[0] = cmpOrder[(idx + 1) % cmp]; 
+                rOrder[1] = cmpOrder[(idx + cmp - 1) % cmp];
+            } else if (cmpOrder[idx] == 3) {
+                rOrder[2] = cmpOrder[(idx + 1) % cmp]; 
+                rOrder[3] = cmpOrder[(idx + cmp - 1) % cmp];
+            }
+        }
+        if (connectTo0[0] != connectTo0[1]) {
+            if (rOrder[0] == rOrder[1]) {
+                printf("%s [%d] (1) order=(%d,%d,%d,%d) r_order=(%d,%d,%d,%d)\n",
+                    __FUNCTION__, (int)index, connectTo0[0], connectTo0[1],
+                    connectTo3[0], connectTo3[1],
+                    rOrder[0], rOrder[1], rOrder[2], rOrder[3]);
+                continue;
+            }
+            int unused = 6 - connectTo0[0] - connectTo0[1];
+            int rUnused = 6 - rOrder[0] - rOrder[1];
+            if (unused != rUnused) {
+                printf("%s [%d] (2) order=(%d,%d,%d,%d) r_order=(%d,%d,%d,%d)\n",
+                    __FUNCTION__, (int)index, connectTo0[0], connectTo0[1],
+                    connectTo3[0], connectTo3[1],
+                    rOrder[0], rOrder[1], rOrder[2], rOrder[3]);
+                continue;
+            }
+        } else {
+            if (rOrder[0] != rOrder[1]) {
+                printf("%s [%d] (3) order=(%d,%d,%d,%d) r_order=(%d,%d,%d,%d)\n",
+                    __FUNCTION__, (int)index, connectTo0[0], connectTo0[1],
+                    connectTo3[0], connectTo3[1],
+                    rOrder[0], rOrder[1], rOrder[2], rOrder[3]);
+                continue;
+            }
+            if (connectTo0[0] != rOrder[0]) {
+                printf("%s [%d] (4) order=(%d,%d,%d,%d) r_order=(%d,%d,%d,%d)\n",
+                    __FUNCTION__, (int)index, connectTo0[0], connectTo0[1],
+                    connectTo3[0], connectTo3[1],
+                    rOrder[0], rOrder[1], rOrder[2], rOrder[3]);
+                continue;
+            }
+        }
+        if (connectTo3[0] != connectTo3[1]) {
+             if (rOrder[2] == rOrder[3]) {
+                printf("%s [%d] (5) order=(%d,%d,%d,%d) r_order=(%d,%d,%d,%d)\n",
+                    __FUNCTION__, (int)index, connectTo0[0], connectTo0[1],
+                    connectTo3[0], connectTo3[1],
+                    rOrder[0], rOrder[1], rOrder[2], rOrder[3]);
+                continue;
+            }
+           int unused = 6 - connectTo3[0] - connectTo3[1];
+           int rUnused = 6 - rOrder[2] - rOrder[3];
+            if (unused != rUnused) {
+                printf("%s [%d] (6) order=(%d,%d,%d,%d) r_order=(%d,%d,%d,%d)\n",
+                    __FUNCTION__, (int)index, connectTo0[0], connectTo0[1],
+                    connectTo3[0], connectTo3[1],
+                    rOrder[0], rOrder[1], rOrder[2], rOrder[3]);
+                continue;
+            }
+        } else {
+            if (rOrder[2] != rOrder[3]) {
+                printf("%s [%d] (7) order=(%d,%d,%d,%d) r_order=(%d,%d,%d,%d)\n",
+                    __FUNCTION__, (int)index, connectTo0[0], connectTo0[1],
+                    connectTo3[0], connectTo3[1],
+                    rOrder[0], rOrder[1], rOrder[2], rOrder[3]);
+                continue;
+            }
+            if (connectTo3[1] != rOrder[3]) {
+                printf("%s [%d] (8) order=(%d,%d,%d,%d) r_order=(%d,%d,%d,%d)\n",
+                    __FUNCTION__, (int)index, connectTo0[0], connectTo0[1],
+                    connectTo3[0], connectTo3[1],
+                    rOrder[0], rOrder[1], rOrder[2], rOrder[3]);
+                continue;
+            }
+        }
+nextTest:
+        ;
+    }
+}
+
+
+
diff --git a/experimental/Intersection/CubeRoot.cpp b/experimental/Intersection/CubeRoot.cpp
new file mode 100644
index 0000000..e188b34
--- /dev/null
+++ b/experimental/Intersection/CubeRoot.cpp
@@ -0,0 +1,387 @@
+// http://metamerist.com/cbrt/CubeRoot.cpp 
+//
+
+#include <math.h>
+#include "CubicIntersection.h"
+
+#define TEST_ALTERNATIVES 0
+#if TEST_ALTERNATIVES
+typedef float  (*cuberootfnf) (float);
+typedef double (*cuberootfnd) (double);
+
+// estimate bits of precision (32-bit float case)
+inline int bits_of_precision(float a, float b)
+{
+	const double kd = 1.0 / log(2.0);
+
+	if (a==b)
+		return 23;
+
+	const double kdmin = pow(2.0, -23.0);
+
+	double d = fabs(a-b);
+	if (d < kdmin)
+		return 23;
+
+	return int(-log(d)*kd);
+}
+
+// estiamte bits of precision (64-bit double case)
+inline int bits_of_precision(double a, double b)
+{
+	const double kd = 1.0 / log(2.0);
+
+	if (a==b)
+		return 52;
+
+	const double kdmin = pow(2.0, -52.0);
+
+	double d = fabs(a-b);
+	if (d < kdmin)
+		return 52;
+
+	return int(-log(d)*kd);
+}
+
+// cube root via x^(1/3)
+static float pow_cbrtf(float x)
+{
+	return (float) pow(x, 1.0f/3.0f);
+}
+
+// cube root via x^(1/3)
+static double pow_cbrtd(double x)
+{
+	return pow(x, 1.0/3.0);
+}
+
+// cube root approximation using bit hack for 32-bit float
+static  float cbrt_5f(float f)
+{
+	unsigned int* p = (unsigned int *) &f;
+	*p = *p/3 + 709921077;
+	return f;
+}
+#endif
+
+// cube root approximation using bit hack for 64-bit float 
+// adapted from Kahan's cbrt
+static  double cbrt_5d(double d)
+{
+	const unsigned int B1 = 715094163;
+	double t = 0.0;
+	unsigned int* pt = (unsigned int*) &t;
+	unsigned int* px = (unsigned int*) &d;
+	pt[1]=px[1]/3+B1;
+	return t;
+}
+
+#if TEST_ALTERNATIVES
+// cube root approximation using bit hack for 64-bit float 
+// adapted from Kahan's cbrt
+#if 0
+static  double quint_5d(double d)
+{
+	return sqrt(sqrt(d));
+
+	const unsigned int B1 = 71509416*5/3;
+	double t = 0.0;
+	unsigned int* pt = (unsigned int*) &t;
+	unsigned int* px = (unsigned int*) &d;
+	pt[1]=px[1]/5+B1;
+	return t;
+}
+#endif
+
+// iterative cube root approximation using Halley's method (float)
+static  float cbrta_halleyf(const float a, const float R)
+{
+	const float a3 = a*a*a;
+    const float b= a * (a3 + R + R) / (a3 + a3 + R);
+	return b;
+}
+#endif
+
+// iterative cube root approximation using Halley's method (double)
+static  double cbrta_halleyd(const double a, const double R)
+{
+	const double a3 = a*a*a;
+    const double b= a * (a3 + R + R) / (a3 + a3 + R);
+	return b;
+}
+
+#if TEST_ALTERNATIVES
+// iterative cube root approximation using Newton's method (float)
+static  float cbrta_newtonf(const float a, const float x)
+{
+//    return (1.0 / 3.0) * ((a + a) + x / (a * a));
+	return a - (1.0f / 3.0f) * (a - x / (a*a));
+}
+
+// iterative cube root approximation using Newton's method (double)
+static  double cbrta_newtond(const double a, const double x)
+{
+	return (1.0/3.0) * (x / (a*a) + 2*a);
+}
+
+// cube root approximation using 1 iteration of Halley's method (double)
+static double halley_cbrt1d(double d)
+{
+	double a = cbrt_5d(d);
+	return cbrta_halleyd(a, d);
+}
+
+// cube root approximation using 1 iteration of Halley's method (float)
+static float halley_cbrt1f(float d)
+{
+	float a = cbrt_5f(d);
+	return cbrta_halleyf(a, d);
+}
+
+// cube root approximation using 2 iterations of Halley's method (double)
+static double halley_cbrt2d(double d)
+{
+	double a = cbrt_5d(d);
+	a = cbrta_halleyd(a, d);
+	return cbrta_halleyd(a, d);
+}
+#endif
+
+// cube root approximation using 3 iterations of Halley's method (double)
+static double halley_cbrt3d(double d)
+{
+	double a = cbrt_5d(d);
+	a = cbrta_halleyd(a, d);
+	a = cbrta_halleyd(a, d);
+	return cbrta_halleyd(a, d);
+}
+
+#if TEST_ALTERNATIVES
+// cube root approximation using 2 iterations of Halley's method (float)
+static float halley_cbrt2f(float d)
+{
+	float a = cbrt_5f(d);
+	a = cbrta_halleyf(a, d);
+	return cbrta_halleyf(a, d);
+}
+
+// cube root approximation using 1 iteration of Newton's method (double)
+static double newton_cbrt1d(double d)
+{
+	double a = cbrt_5d(d);
+	return cbrta_newtond(a, d);
+}
+
+// cube root approximation using 2 iterations of Newton's method (double)
+static double newton_cbrt2d(double d)
+{
+	double a = cbrt_5d(d);
+	a = cbrta_newtond(a, d);
+	return cbrta_newtond(a, d);
+}
+
+// cube root approximation using 3 iterations of Newton's method (double)
+static double newton_cbrt3d(double d)
+{
+	double a = cbrt_5d(d);
+	a = cbrta_newtond(a, d);
+	a = cbrta_newtond(a, d);
+	return cbrta_newtond(a, d);
+}
+
+// cube root approximation using 4 iterations of Newton's method (double)
+static double newton_cbrt4d(double d)
+{
+	double a = cbrt_5d(d);
+	a = cbrta_newtond(a, d);
+	a = cbrta_newtond(a, d);
+	a = cbrta_newtond(a, d);
+	return cbrta_newtond(a, d);
+}
+
+// cube root approximation using 2 iterations of Newton's method (float)
+static float newton_cbrt1f(float d)
+{
+	float a = cbrt_5f(d);
+	return cbrta_newtonf(a, d);
+}
+
+// cube root approximation using 2 iterations of Newton's method (float)
+static float newton_cbrt2f(float d)
+{
+	float a = cbrt_5f(d);
+	a = cbrta_newtonf(a, d);
+	return cbrta_newtonf(a, d);
+}
+
+// cube root approximation using 3 iterations of Newton's method (float)
+static float newton_cbrt3f(float d)
+{
+	float a = cbrt_5f(d);
+	a = cbrta_newtonf(a, d);
+	a = cbrta_newtonf(a, d);
+	return cbrta_newtonf(a, d);
+}
+
+// cube root approximation using 4 iterations of Newton's method (float)
+static float newton_cbrt4f(float d)
+{
+	float a = cbrt_5f(d);
+	a = cbrta_newtonf(a, d);
+	a = cbrta_newtonf(a, d);
+	a = cbrta_newtonf(a, d);
+	return cbrta_newtonf(a, d);
+}
+
+static double TestCubeRootf(const char* szName, cuberootfnf cbrt, double rA, double rB, int rN)
+{
+	const int N = rN;
+ 	
+	float dd = float((rB-rA) / N);
+
+	// calculate 1M numbers
+	int i=0;
+	float d = (float) rA;
+
+	double s = 0.0;
+
+	for(d=(float) rA, i=0; i<N; i++, d += dd)
+	{
+		s += cbrt(d);
+	}
+
+	double bits = 0.0;
+	double worstx=0.0;
+	double worsty=0.0;
+	int minbits=64;
+
+	for(d=(float) rA, i=0; i<N; i++, d += dd)
+	{
+		float a = cbrt((float) d);	
+		float b = (float) pow((double) d, 1.0/3.0);
+
+		int bc = bits_of_precision(a, b);
+		bits += bc;
+
+		if (b > 1.0e-6)
+		{
+			if (bc < minbits)
+			{
+				minbits = bc;
+				worstx = d;
+				worsty = a;
+			}
+		}
+	}
+
+	bits /= N;
+
+    printf(" %3d mbp  %6.3f abp\n", minbits, bits);
+
+	return s;
+}
+
+
+static double TestCubeRootd(const char* szName, cuberootfnd cbrt, double rA, double rB, int rN)
+{
+	const int N = rN;
+	
+	double dd = (rB-rA) / N;
+
+	int i=0;
+	
+	double s = 0.0;
+	double d = 0.0;
+
+	for(d=rA, i=0; i<N; i++, d += dd)
+	{
+		s += cbrt(d);
+	}
+
+
+	double bits = 0.0;
+	double worstx = 0.0;
+	double worsty = 0.0;
+	int minbits = 64;
+	for(d=rA, i=0; i<N; i++, d += dd)
+	{
+		double a = cbrt(d);	
+		double b = pow(d, 1.0/3.0);
+
+		int bc = bits_of_precision(a, b); // min(53, count_matching_bitsd(a, b) - 12);
+		bits += bc;
+
+		if (b > 1.0e-6)
+		{
+			if (bc < minbits)
+			{
+				bits_of_precision(a, b);
+				minbits = bc; 
+				worstx = d;
+				worsty = a;
+			}
+		}
+	}
+
+	bits /= N;
+
+    printf(" %3d mbp  %6.3f abp\n", minbits, bits);
+
+	return s;
+}
+
+static int _tmain()
+{
+	// a million uniform steps through the range from 0.0 to 1.0
+	// (doing uniform steps in the log scale would be better)
+	double a = 0.0;
+	double b = 1.0;
+	int n = 1000000;
+
+	printf("32-bit float tests\n");
+	printf("----------------------------------------\n");
+	TestCubeRootf("cbrt_5f", cbrt_5f, a, b, n);
+	TestCubeRootf("pow", pow_cbrtf, a, b, n);
+	TestCubeRootf("halley x 1", halley_cbrt1f, a, b, n);
+	TestCubeRootf("halley x 2", halley_cbrt2f, a, b, n);
+	TestCubeRootf("newton x 1", newton_cbrt1f, a, b, n);
+	TestCubeRootf("newton x 2", newton_cbrt2f, a, b, n);
+	TestCubeRootf("newton x 3", newton_cbrt3f, a, b, n);
+	TestCubeRootf("newton x 4", newton_cbrt4f, a, b, n);
+	printf("\n\n");
+
+	printf("64-bit double tests\n");
+	printf("----------------------------------------\n");
+	TestCubeRootd("cbrt_5d", cbrt_5d, a, b, n);
+	TestCubeRootd("pow", pow_cbrtd, a, b, n);
+	TestCubeRootd("halley x 1", halley_cbrt1d, a, b, n);
+	TestCubeRootd("halley x 2", halley_cbrt2d, a, b, n);
+	TestCubeRootd("halley x 3", halley_cbrt3d, a, b, n);
+	TestCubeRootd("newton x 1", newton_cbrt1d, a, b, n);
+	TestCubeRootd("newton x 2", newton_cbrt2d, a, b, n);
+	TestCubeRootd("newton x 3", newton_cbrt3d, a, b, n);
+	TestCubeRootd("newton x 4", newton_cbrt4d, a, b, n);
+	printf("\n\n");
+
+	return 0;
+}
+#endif
+
+double cube_root(double x) {
+    return halley_cbrt3d(x);
+}
+
+#if TEST_ALTERNATIVES
+// http://bytes.com/topic/c/answers/754588-tips-find-cube-root-program-using-c
+/* cube root */
+int icbrt(int n) {
+    int t=0, x=(n+2)/3; /* works for n=0 and n>=1 */
+    for(; t!=x;) {
+        int x3=x*x*x;
+        t=x;
+        x*=(2*n + x3);
+        x/=(2*x3 + n);
+    }
+    return x ; /* always(?) equal to floor(n^(1/3)) */
+}
+#endif
diff --git a/experimental/Intersection/CubicBezierClip.cpp b/experimental/Intersection/CubicBezierClip.cpp
new file mode 100644
index 0000000..1b5c60c
--- /dev/null
+++ b/experimental/Intersection/CubicBezierClip.cpp
@@ -0,0 +1,83 @@
+#include "CubicIntersection.h"
+#include "LineParameters.h"
+#include <algorithm> // used for std::swap
+
+// return false if unable to clip (e.g., unable to create implicit line)
+// caller should subdivide, or create degenerate if the values are too small
+bool bezier_clip(const Cubic& cubic1, const Cubic& cubic2, double& minT, double& maxT) {
+    minT = 1;
+    maxT = 0;
+    // determine normalized implicit line equation for pt[0] to pt[3]
+    //   of the form ax + by + c = 0, where a*a + b*b == 1
+    
+    // find the implicit line equation parameters
+    LineParameters endLine;
+    endLine.cubicEndPoints(cubic1);
+    if (!endLine.normalize()) {
+        printf("line cannot be normalized: need more code here\n");
+        return false;
+    }
+
+    double distance[2];
+    endLine.controlPtDistance(cubic1, distance);
+    
+    // find fat line
+    double top = distance[0];
+    double bottom = distance[1];
+    if (top > bottom) {
+        std::swap(top, bottom);
+    }
+    if (top * bottom >= 0) {
+        const double scale = 3/4.0; // http://cagd.cs.byu.edu/~tom/papers/bezclip.pdf (13)
+        if (top < 0) {
+            top *= scale;
+            bottom = 0;
+        } else {
+            top = 0;
+            bottom *= scale;
+        }
+    } else {
+        const double scale = 4/9.0; // http://cagd.cs.byu.edu/~tom/papers/bezclip.pdf (15)
+        top *= scale;
+        bottom *= scale;
+    }
+    
+    // compute intersecting candidate distance
+    Cubic distance2y; // points with X of (0, 1/3, 2/3, 1)
+    endLine.cubicDistanceY(cubic2, distance2y);
+    
+    int flags = 0;
+    if (approximately_lesser(distance2y[0].y, top)) {
+        flags |= kFindTopMin;
+    } else if (approximately_greater(distance2y[0].y, bottom)) {
+        flags |= kFindBottomMin;
+    } else {
+        minT = 0;
+    }
+
+    if (approximately_lesser(distance2y[3].y, top)) {
+        flags |= kFindTopMax;
+    } else if (approximately_greater(distance2y[3].y, bottom)) {
+        flags |= kFindBottomMax;
+    } else {
+        maxT = 1;
+    }
+    // Find the intersection of distance convex hull and fat line.
+    char to_0[2];
+    char to_3[2];
+    bool do_1_2_edge = convex_x_hull(distance2y, to_0, to_3);
+    x_at(distance2y[0], distance2y[to_0[0]], top, bottom, flags, minT, maxT);
+    if (to_0[0] != to_0[1]) {
+        x_at(distance2y[0], distance2y[to_0[1]], top, bottom, flags, minT, maxT);
+    }
+    x_at(distance2y[to_3[0]], distance2y[3], top, bottom, flags, minT, maxT);
+    if (to_3[0] != to_3[1]) {
+        x_at(distance2y[to_3[1]], distance2y[3], top, bottom, flags, minT, maxT);
+    }
+    if (do_1_2_edge) {
+        x_at(distance2y[1], distance2y[2], top, bottom, flags, minT, maxT);
+    }
+    
+    return minT < maxT; // returns false if distance shows no intersection
+}
+
diff --git a/experimental/Intersection/CubicBezierClip_Test.cpp b/experimental/Intersection/CubicBezierClip_Test.cpp
new file mode 100644
index 0000000..3228d89
--- /dev/null
+++ b/experimental/Intersection/CubicBezierClip_Test.cpp
@@ -0,0 +1,24 @@
+#include "CubicIntersection.h"
+#include "CubicIntersection_TestData.h"
+#include "CubicIntersection_Tests.h"
+
+void BezierClip_Test() {
+    for (size_t index = 0; index < tests_count; ++index) {
+        const Cubic& cubic1 = tests[index][0];
+        const Cubic& cubic2 = tests[index][1];
+        Cubic reduce1, reduce2;
+        int order1 = reduceOrder(cubic1, reduce1, kReduceOrder_NoQuadraticsAllowed);
+        int order2 = reduceOrder(cubic2, reduce2, kReduceOrder_NoQuadraticsAllowed);
+        if (order1 < 4) {
+            printf("%s [%d] cubic1 order=%d\n", __FUNCTION__, (int) index, order1);
+        }
+        if (order2 < 4) {
+            printf("%s [%d] cubic2 order=%d\n", __FUNCTION__, (int) index, order2);
+        }
+        if (order1 == 4 && order2 == 4) {
+            double minT = 0;
+            double maxT = 1;
+            bezier_clip(reduce1, reduce2, minT, maxT);
+        }
+    }
+}
diff --git a/experimental/Intersection/CubicCoincidence.cpp b/experimental/Intersection/CubicCoincidence.cpp
new file mode 100644
index 0000000..e2ec641
--- /dev/null
+++ b/experimental/Intersection/CubicCoincidence.cpp
@@ -0,0 +1,45 @@
+/*
+ Suppose two cubics are coincident. Then a third cubic exists described by two
+ of the four endpoints. The coincident endpoints must be on or inside the convex
+ hulls of both cubics. 
+ 
+ The coincident control points, while unknown, must lie on the line segment from
+ the coincident end point to its original control point. 
+ 
+ Given a cubic c1, c2, A, and D:
+ A = c1[0]*(1 - t1)*(1 - t1)*(1 - t1) + 3*c1[1]*t1*(1 - t1)*(1 - t1) + 3*c1[2]*t1*t1*(1 - t1) + c1[3]*t1*t1*t1
+ D = c2[0]*(1 - t2)*(1 - t2)*(1 - t2) + 3*c2[1]*t2*(1 - t2)*(1 - t2) + 3*c2[2]*t2*t2*(1 - t2) + c213]*t2*t2*t2
+ 
+ Assuming that A was originally from c2:
+ 
+ B = c2[0]*(1 - t2) + c2[1]*t2
+ C = c1[0]*(1 - t1) + c1[0]*t1
+ 
+ 
+ 
+ If both ends of the same cubic is contained in the convex hull of the other,
+ then, are there t values of the larger cubic that describe the end points; is
+ the mid-value of those ts on the smaller cubic, and, are the tangents at the
+ smaller ends the same as on both.
+ 
+ This all requires knowing the t values.
+ 
+ 
+ 
+ Maybe solving the cubic is possible. Given x, find t. Given t, find y.
+  
+see en.wikipedia.org/wiki/Cubic_polynomial
+ 
+ Another way of writing the solution may be obtained by noting that the proof of above formula shows that the product of the two cube roots is rational. This gives the following formula in which  or  stands for any choice of the square or cube root, if 
+
+If  and , the sign of  has to be chosen to have .
+If  and , the three roots are equal:
+
+If Q = 0 and , the above expression for the roots is correct but misleading, hiding the fact that no radical is needed to represent the roots. In fact, in this case, there is a double root,
+
+and a simple root
+
+ 
+
+ */
+
diff --git a/experimental/Intersection/CubicIntersection.cpp b/experimental/Intersection/CubicIntersection.cpp
new file mode 100644
index 0000000..c41e0d7
--- /dev/null
+++ b/experimental/Intersection/CubicIntersection.cpp
@@ -0,0 +1,80 @@
+#include "CubicIntersection.h"
+#include "Intersections.h"
+#include "LineIntersection.h"
+
+static bool chop(const Cubic& smaller, const Cubic& larger,
+        Intersections& intersections, double minT, double maxT);
+        
+static bool intersect(const Cubic& smaller, const Cubic& larger,
+        Intersections& intersections) {
+    // FIXME: carry order with cubic so we don't call it repeatedly
+    Cubic smallResult;
+    if (reduceOrder(smaller, smallResult, kReduceOrder_NoQuadraticsAllowed) <= 2) {
+        Cubic largeResult;
+        if (reduceOrder(larger, largeResult, kReduceOrder_NoQuadraticsAllowed) <= 2) {
+            _Point pt;
+            const _Line& smallLine = (const _Line&) smallResult;
+            const _Line& largeLine = (const _Line&) largeResult;
+            if (!lineIntersect(smallLine, largeLine, &pt)) {
+                return false;
+            }
+            double smallT = t_at(smallLine, pt);
+            double largeT = t_at(largeLine, pt);
+            intersections.add(smallT, largeT);
+            return true;
+        }
+    }
+    double minT, maxT;
+    if (!bezier_clip(smaller, larger, minT, maxT)) {
+        if (minT == maxT) {
+            intersections.add(minT, 0.5);
+            return true;
+        }
+        return false;
+    }
+    return chop(larger, smaller, intersections, minT, maxT);
+}
+
+bool chop(const Cubic& smaller, const Cubic& larger,
+        Intersections& intersections, double minT, double maxT) {
+    intersections.swap();
+    if (maxT - minT > 0.8) { // http://cagd.cs.byu.edu/~tom/papers/bezclip.pdf Multiple intersections
+        CubicPair cubicPair;
+        chop_at(larger, cubicPair, 0.5);
+        int used = intersections.used();
+        if (intersect(cubicPair.first(), smaller, intersections)) {
+            intersections.offset(used, 0, 0.5);
+        }
+        used = intersections.used();
+        if (intersect(cubicPair.second(), smaller, intersections)) {
+            intersections.offset(used, 0.5, 1);
+        }
+        intersections.swap();
+        return intersections.intersected();
+    }
+    Cubic cut;
+    sub_divide(smaller, minT, maxT, cut);
+    int used = intersections.used();
+    bool result = intersect(cut, larger, intersections);
+    if (result) {
+        intersections.offset(used, minT, maxT);
+    }
+    intersections.swap();
+    return result;
+}
+
+bool intersectStart(const Cubic& cubic1, const Cubic& cubic2,
+        Intersections& intersections) {
+    double minT1, minT2, maxT1, maxT2;
+    if (!bezier_clip(cubic2, cubic1, minT1, maxT1)) {
+        return false;
+    }
+    if (!bezier_clip(cubic1, cubic2, minT2, maxT2)) {
+        return false;
+    }
+    if (maxT1 - minT1 < maxT2 - minT2) {
+        intersections.swap();
+        return chop(cubic1, cubic2, intersections, minT1, maxT1);
+    } 
+    return chop(cubic2, cubic1, intersections, minT2, maxT2);
+}
diff --git a/experimental/Intersection/CubicIntersection.h b/experimental/Intersection/CubicIntersection.h
new file mode 100644
index 0000000..dd6c52e
--- /dev/null
+++ b/experimental/Intersection/CubicIntersection.h
@@ -0,0 +1,31 @@
+#include "DataTypes.h"
+
+class Intersections;
+
+// unit-testable utilities
+bool bezier_clip(const Cubic& cubic1, const Cubic& cubic2, double& minT, double& maxT);
+bool bezier_clip(const Quadratic& q1, const Quadratic& q2, double& minT, double& maxT);
+void chop_at(const Cubic& src, CubicPair& dst, double t);
+void chop_at(const Quadratic& src, QuadraticPair& dst, double t);
+int convex_hull(const Cubic& cubic, char order[4]);
+bool convex_x_hull(const Cubic& cubic, char connectTo0[2], char connectTo3[2]);
+double cube_root(double x);
+int cubic_roots(const double coeff[4], double tValues[3]);
+bool implicit_matches(const Cubic& cubic1, const Cubic& cubic2);
+bool implicit_matches(const Quadratic& quad1, const Quadratic& quad2);
+int quadratic_roots(const double coeff[3], double tValues[2]);
+void sub_divide(const Cubic& src, double t1, double t2, Cubic& dst);
+void sub_divide(const Quadratic& src, double t1, double t2, Quadratic& dst);
+void tangent(const Cubic& cubic, double t, _Point& result);
+void tangent(const Quadratic& cubic, double t, _Point& result);
+
+// main functions
+enum ReduceOrder_Flags {
+    kReduceOrder_NoQuadraticsAllowed,
+    kReduceOrder_QuadraticsAllowed
+};
+int reduceOrder(const Cubic& cubic, Cubic& reduction, ReduceOrder_Flags );
+int reduceOrder(const Quadratic& quad, Quadratic& reduction);
+bool intersectStart(const Cubic& cubic1, const Cubic& cubic2, Intersections& );
+bool intersectStartT(const Cubic& cubic1, const Cubic& cubic2, Intersections& );
+bool intersectStart(const Quadratic& q1, const Quadratic& q2, Intersections& );
diff --git a/experimental/Intersection/CubicIntersectionT.cpp b/experimental/Intersection/CubicIntersectionT.cpp
new file mode 100644
index 0000000..d58bcb9
--- /dev/null
+++ b/experimental/Intersection/CubicIntersectionT.cpp
@@ -0,0 +1,145 @@
+#include "CubicIntersection.h"
+#include "Intersections.h"
+#include "IntersectionUtilities.h"
+#include "LineIntersection.h"
+
+class CubicIntersections : public Intersections {
+public:
+
+CubicIntersections(const Cubic& c1, const Cubic& c2, Intersections& i) 
+    : cubic1(c1)
+    , cubic2(c2)
+    , intersections(i)
+    , depth(0) 
+    , splits(0) {
+}
+
+bool intersect() {
+    double minT1, minT2, maxT1, maxT2;
+    if (!bezier_clip(cubic2, cubic1, minT1, maxT1)) {
+        return false;
+    }
+    if (!bezier_clip(cubic1, cubic2, minT2, maxT2)) {
+        return false;
+    }
+    int split;
+    if (maxT1 - minT1 < maxT2 - minT2) {
+        intersections.swap();
+        minT2 = 0;
+        maxT2 = 1;
+        split = maxT1 - minT1 > tClipLimit;
+    } else {
+        minT1 = 0;
+        maxT1 = 1;
+        split = (maxT2 - minT2 > tClipLimit) << 1;
+    }
+    return chop(minT1, maxT1, minT2, maxT2, split);
+}
+
+protected:
+        
+bool intersect(double minT1, double maxT1, double minT2, double maxT2) {
+    Cubic smaller, larger;
+    // FIXME: carry last subdivide and reduceOrder result with cubic 
+    sub_divide(cubic1, minT1, maxT1, intersections.swapped() ? larger : smaller);
+    sub_divide(cubic2, minT2, maxT2, intersections.swapped() ? smaller : larger);
+    Cubic smallResult;
+    if (reduceOrder(smaller, smallResult, kReduceOrder_NoQuadraticsAllowed) <= 2) {
+        Cubic largeResult;
+        if (reduceOrder(larger, largeResult, kReduceOrder_NoQuadraticsAllowed) <= 2) {
+            _Point pt;
+            const _Line& smallLine = (const _Line&) smallResult;
+            const _Line& largeLine = (const _Line&) largeResult;
+            if (!lineIntersect(smallLine, largeLine, &pt)) {
+                return false;
+            }
+            double smallT = t_at(smallLine, pt);
+            double largeT = t_at(largeLine, pt);
+            if (intersections.swapped()) {
+                smallT = interp(minT2, maxT2, smallT); 
+                largeT = interp(minT1, maxT1, largeT); 
+            } else {
+                smallT = interp(minT1, maxT1, smallT); 
+                largeT = interp(minT2, maxT2, largeT); 
+            }
+            intersections.add(smallT, largeT);
+            return true;
+        }
+    }
+    double minT, maxT;
+    if (!bezier_clip(smaller, larger, minT, maxT)) {
+        if (minT == maxT) {
+            if (intersections.swapped()) {
+                minT1 = (minT1 + maxT1) / 2;
+                minT2 = interp(minT2, maxT2, minT);
+            } else {
+                minT1 = interp(minT1, maxT1, minT);
+                minT2 = (minT2 + maxT2) / 2;
+            }
+            intersections.add(minT1, minT2);
+            return true;
+        }
+        return false;
+    }
+    
+    int split;
+    if (intersections.swapped()) {
+        double newMinT1 = interp(minT1, maxT1, minT);
+        double newMaxT1 = interp(minT1, maxT1, maxT);
+        split = (newMaxT1 - newMinT1 > (maxT1 - minT1) * tClipLimit) << 1;
+        printf("%s d=%d s=%d new1=(%g,%g) old1=(%g,%g) split=%d\n", __FUNCTION__, depth,
+            splits, newMinT1, newMaxT1, minT1, maxT1, split);
+        minT1 = newMinT1;
+        maxT1 = newMaxT1;
+    } else {
+        double newMinT2 = interp(minT2, maxT2, minT);
+        double newMaxT2 = interp(minT2, maxT2, maxT);
+        split = newMaxT2 - newMinT2 > (maxT2 - minT2) * tClipLimit;
+        printf("%s d=%d s=%d new2=(%g,%g) old2=(%g,%g) split=%d\n", __FUNCTION__, depth,
+            splits, newMinT2, newMaxT2, minT2, maxT2, split);
+        minT2 = newMinT2;
+        maxT2 = newMaxT2;
+    }
+    return chop(minT1, maxT1, minT2, maxT2, split);
+}
+
+bool chop(double minT1, double maxT1, double minT2, double maxT2, int split) {
+    ++depth;
+    intersections.swap();
+    if (split) {
+        ++splits;
+        if (split & 2) {
+            double middle1 = (maxT1 + minT1) / 2;
+            intersect(minT1, middle1, minT2, maxT2);
+            intersect(middle1, maxT1, minT2, maxT2);
+        } else {
+            double middle2 = (maxT2 + minT2) / 2;
+            intersect(minT1, maxT1, minT2, middle2);
+            intersect(minT1, maxT1, middle2, maxT2);
+        }
+        --splits;
+        intersections.swap();
+        --depth;
+        return intersections.intersected();
+    }
+    bool result = intersect(minT1, maxT1, minT2, maxT2);
+    intersections.swap();
+    --depth;
+    return result;
+}
+
+private:
+
+static const double tClipLimit = 0.8; // http://cagd.cs.byu.edu/~tom/papers/bezclip.pdf see Multiple intersections
+const Cubic& cubic1;
+const Cubic& cubic2;
+Intersections& intersections;
+int depth;
+int splits;
+};
+
+bool intersectStartT(const Cubic& c1, const Cubic& c2, Intersections& i) {
+    CubicIntersections c(c1, c2, i);
+    return c.intersect();
+}
+
diff --git a/experimental/Intersection/CubicIntersection_Test.cpp b/experimental/Intersection/CubicIntersection_Test.cpp
new file mode 100644
index 0000000..0ae1749
--- /dev/null
+++ b/experimental/Intersection/CubicIntersection_Test.cpp
@@ -0,0 +1,65 @@
+#include "CubicIntersection.h"
+#include "CubicIntersection_TestData.h"
+#include "CubicIntersection_Tests.h"
+#include "Intersections.h"
+#include "TestUtilities.h"
+
+const int firstCubicIntersectionTest = 9;
+
+void CubicIntersection_Test() {
+    for (size_t index = firstCubicIntersectionTest; index < tests_count; ++index) {
+        const Cubic& cubic1 = tests[index][0];
+        const Cubic& cubic2 = tests[index][1];
+        Cubic reduce1, reduce2;
+        int order1 = reduceOrder(cubic1, reduce1, kReduceOrder_NoQuadraticsAllowed);
+        int order2 = reduceOrder(cubic2, reduce2, kReduceOrder_NoQuadraticsAllowed);
+        if (order1 < 4) {
+            printf("[%d] cubic1 order=%d\n", (int) index, order1);
+        }
+        if (order2 < 4) {
+            printf("[%d] cubic2 order=%d\n", (int) index, order2);
+        }
+        if (order1 == 4 && order2 == 4) {
+            Intersections tIntersections;
+            intersectStartT(reduce1, reduce2, tIntersections);
+#ifdef COMPUTE_BY_CUBIC_SUBDIVISION
+            Intersections intersections;
+            intersectStart(reduce1, reduce2, intersections);
+#endif
+            if (tIntersections.intersected()) {
+                for (int pt = 0; pt < tIntersections.used(); ++pt) {
+#ifdef COMPUTE_BY_CUBIC_SUBDIVISION
+                    double t1 = intersections.fT[0][pt];
+                    double x1, y1;
+                    xy_at_t(cubic1, t1, x1, y1);
+                    double t2 = intersections.fT[1][pt];
+                    double x2, y2;
+                    xy_at_t(cubic2, t2, x2, y2);
+                    if (!approximately_equal(x1, x2)) {
+                        printf("%s [%d] (1) t1=%g (%g,%g) t2=%g (%g,%g)\n",
+                            __FUNCTION__, pt, t1, x1, y1, t2, x2, y2);
+                    }
+                    if (!approximately_equal(y1, y2)) {
+                        printf("%s [%d] (2) t1=%g (%g,%g) t2=%g (%g,%g)\n",
+                            __FUNCTION__, pt, t1, x1, y1, t2, x2, y2);
+                    }
+#endif
+                    double tt1 = tIntersections.fT[0][pt];
+                    double tx1, ty1;
+                    xy_at_t(cubic1, tt1, tx1, ty1);
+                    double tt2 = tIntersections.fT[1][pt];
+                    double tx2, ty2;
+                    xy_at_t(cubic2, tt2, tx2, ty2);
+                    if (!approximately_equal(tx1, tx2)) {
+                        printf("%s [%d,%d] x!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
+                            __FUNCTION__, (int)index, pt, tt1, tx1, ty1, tt2, tx2, ty2);
+                    }
+                    if (!approximately_equal(ty1, ty2)) {
+                        printf("%s [%d,%d] y!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
+                            __FUNCTION__, (int)index, pt, tt1, tx1, ty1, tt2, tx2, ty2);
+                    }
+                }
+            }
+        }
+    }
+}
diff --git a/experimental/Intersection/CubicIntersection_TestData.cpp b/experimental/Intersection/CubicIntersection_TestData.cpp
new file mode 100644
index 0000000..95d44d7
--- /dev/null
+++ b/experimental/Intersection/CubicIntersection_TestData.cpp
@@ -0,0 +1,373 @@
+#include "CubicIntersection_TestData.h"
+#include <limits>
+
+const Cubic pointDegenerates[] = {
+    {{0, 0}, {0, 0}, {0, 0}, {0, 0}},
+    {{1, 1}, {1, 1}, {1, 1}, {1, 1}},
+    {{1 + PointEpsilon - std::numeric_limits<double>::epsilon(), 1},
+     {1, 1 + PointEpsilon - std::numeric_limits<double>::epsilon()}, {1, 1}, {1, 1}},
+    {{1 + PointEpsilon/2 - std::numeric_limits<double>::epsilon(), 1},
+     {1 - (PointEpsilon/2 - std::numeric_limits<double>::epsilon()), 1}, {1, 1}, {1, 1}}
+};
+
+const size_t pointDegenerates_count = sizeof(pointDegenerates) / sizeof(pointDegenerates[0]);
+    
+const Cubic notPointDegenerates[] = {
+    {{1 + PointEpsilon + std::numeric_limits<double>::epsilon(), 1}, {1, 1 + PointEpsilon}, {1, 1}, {1, 1}},
+    {{1 + PointEpsilon/2 + std::numeric_limits<double>::epsilon(), 1}, {1 - PointEpsilon/2, 1}, {1, 1}, {1, 1}}
+};
+
+const size_t notPointDegenerates_count = sizeof(notPointDegenerates) / sizeof(notPointDegenerates[0]);
+
+// from http://www.truetex.com/bezint.htm
+const Cubic tests[][2] = {
+    { // intersects in one place (data gives bezier clip fits
+     {{0, 45},
+      {6.0094158284751593, 51.610357411322688},
+      {12.741093228940867, 55.981703949474607},
+      {20.021417396476362, 58.652245509710262}},
+     {{2.2070737699246674, 52.703494107327209},
+      {31.591482272629477, 23.811002295222025},
+      {76.824588616426425, 44.049473790502674},
+      {119.25488947221436, 55.599248272955073}}
+    },
+    { // intersects in three places
+        {{0, 45}, {50, 100}, {150,   0}, {200, 55}},
+        {{0, 55}, {50,   0}, {150, 100}, {200, 45}}
+    },
+    { // intersects in one place, cross over is nearly parallel
+        {{0,   0}, {0, 100}, {200,   0}, {200, 100}},
+        {{0, 100}, {0,   0}, {200, 100}, {200,   0}}
+    },
+    { // intersects in two places
+        {{0,   0}, {0, 100}, {200, 100}, {200,   0}},
+        {{0, 100}, {0,   0}, {200,   0}, {200, 100}}
+    },
+    {
+        {{150, 100}, {150 + 0.1, 150}, {150, 200}, {150, 250}},
+        {{250, 150}, {200, 150 + 0.1}, {150, 150}, {100, 150}}
+    },
+    { // single intersection around 168,185
+        {{200, 100}, {150, 100}, {150, 150}, {200, 150}},
+        {{250, 150}, {250, 100}, {100, 100}, {100, 150}}
+    },
+    {
+        {{1.0, 1.5}, {15.5, 0.5}, {-8.0, 3.5}, {5.0, 1.5}},
+        {{4.0, 0.5}, {5.0, 15.0}, {2.0, -8.5}, {4.0, 4.5}}
+    },
+    {
+        {{664.00168, 0},       {726.11545, 124.22757}, {736.89069, 267.89743}, {694.0017, 400.0002}},
+        {{850.66843, 115.55563}, {728.515, 115.55563}, {725.21347, 275.15309}, {694.0017, 400.0002}}
+    },
+    {
+        {{1, 1},   {12.5, 6.5}, {-4, 6.5}, {7.5, 1}},
+        {{1, 6.5}, {12.5, 1},   {-4, 1},   {.5, 6}}
+    },
+    {
+        {{315.748, 312.84}, {312.644, 318.134}, {305.836, 319.909}, {300.542, 316.804}},
+        {{317.122, 309.05}, {316.112, 315.102}, {310.385, 319.19},  {304.332, 318.179}}
+    },
+    {
+        {{1046.604051, 172.937967},  {1046.604051, 178.9763059}, {1041.76745,  183.9279165}, {1035.703842, 184.0432409}},        
+        {{1046.452235, 174.7640504}, {1045.544872, 180.1973817}, {1040.837966, 184.0469882}, {1035.505925, 184.0469882}}
+    },
+    {
+        {{125.79356, 199.57382}, {51.16556, 128.93575}, {87.494,  16.67848}, {167.29361, 16.67848}},
+        {{167.29361, 55.81876}, {100.36128, 55.81876}, {68.64099, 145.4755}, {125.7942, 199.57309}}
+    }
+};
+
+const size_t tests_count = sizeof(tests) / sizeof(tests[0]);
+
+Cubic hexTests[][2] = {
+    {   
+        {{0}} // placeholder for hex converted below
+    }
+};
+
+const size_t hexTests_count = sizeof(hexTests) / sizeof(hexTests[0]);
+
+static const uint64_t testx[2][8] = {
+    {
+        0xf0d0d1ca63075a40LLU, 0x9408ce996a237740LLU, 0x6d5675460fbe5e40LLU, 0x6ef501e1b7487940LLU,
+        0x9a71d2f8143d6540LLU, 0x6bc18bbe02907a40LLU, 0x5b94d92093aa6b40LLU, 0x6ac18bbe02907a40LLU
+    },
+    {
+        0x92c56ed7b6145d40LLU, 0xede4f1255edb7740LLU, 0x1138c1101af75940LLU, 0x42e4f1255edb7740LLU,
+        0x408e51603ad95640LLU, 0x1e2e8fe9dd927740LLU, 0x1cb4777cd3a75440LLU, 0x212e1390de017740LLU 
+    }
+};
+
+void convert_testx() {
+    const uint64_t* inPtr = testx[0];
+    double* outPtr = &hexTests[sizeof(tests) / sizeof(tests[0]) - 1][0][0].x;
+    for (unsigned index = 0; index < sizeof(testx) / sizeof(testx[0][0]); ++index) {
+        uint64_t input = *inPtr++;
+        unsigned char* output = (unsigned char*) outPtr++;
+        for (unsigned byte = 0; byte < sizeof(input); ++byte) {
+            output[byte] = input >> (7 - byte) * 8;
+        }
+    }
+}
+
+const Cubic lines[] = {
+    {{0, 0}, {0, 0}, {0, 0}, {1, 0}}, // 0: horizontal
+    {{0, 0}, {0, 0}, {1, 0}, {0, 0}},
+    {{0, 0}, {1, 0}, {0, 0}, {0, 0}},
+    {{1, 0}, {0, 0}, {0, 0}, {0, 0}},
+    {{1, 0}, {2, 0}, {3, 0}, {4, 0}},
+    {{0, 0}, {0, 0}, {0, 0}, {0, 1}}, // 5: vertical
+    {{0, 0}, {0, 0}, {0, 1}, {0, 0}},
+    {{0, 0}, {0, 1}, {0, 0}, {0, 0}},
+    {{0, 1}, {0, 0}, {0, 0}, {0, 0}},
+    {{0, 1}, {0, 2}, {0, 3}, {0, 4}},
+    {{0, 0}, {0, 0}, {0, 0}, {1, 1}}, // 10: 3 coincident
+    {{0, 0}, {0, 0}, {1, 1}, {0, 0}},
+    {{0, 0}, {1, 1}, {0, 0}, {0, 0}},
+    {{1, 1}, {0, 0}, {0, 0}, {0, 0}},
+    {{0, 0}, {0, 0}, {1, 1}, {2, 2}}, // 14: 2 coincident
+    {{0, 0}, {1, 1}, {0, 0}, {2, 2}},
+    {{0, 0}, {1, 1}, {2, 2}, {0, 0}},
+    {{1, 1}, {0, 0}, {0, 0}, {2, 2}}, // 17:
+    {{1, 1}, {0, 0}, {2, 2}, {0, 0}},
+    {{1, 1}, {2, 2}, {0, 0}, {0, 0}},
+    {{1, 1}, {2, 2}, {3, 3}, {2, 2}}, // middle-last coincident
+    {{1, 1}, {2, 2}, {3, 3}, {3, 3}}, // middle-last coincident
+    {{1, 1}, {1, 1}, {2, 2}, {2, 2}}, // 2 pairs coincident
+    {{1, 1}, {2, 2}, {1, 1}, {2, 2}},
+    {{1, 1}, {2, 2}, {2, 2}, {1, 1}},
+    {{1, 1}, {1, 1}, {3, 3}, {3, 3}}, // first-middle middle-last coincident
+    {{1, 1}, {2, 2}, {3, 3}, {4, 4}}, // no coincident
+    {{1, 1}, {3, 3}, {2, 2}, {4, 4}},
+    {{1, 1}, {2, 2}, {4, 4}, {3, 3}},
+    {{1, 1}, {3, 3}, {4, 4}, {2, 2}},
+    {{1, 1}, {4, 4}, {2, 2}, {3, 3}},
+    {{1, 1}, {4, 4}, {3, 3}, {2, 2}},
+    {{2, 2}, {1, 1}, {3, 3}, {4, 4}},
+    {{2, 2}, {1, 1}, {4, 4}, {3, 3}},
+    {{2, 2}, {3, 3}, {1, 1}, {4, 4}},
+    {{2, 2}, {3, 3}, {4, 4}, {1, 1}},
+    {{2, 2}, {4, 4}, {1, 1}, {3, 3}},
+    {{2, 2}, {4, 4}, {3, 3}, {1, 1}},
+};
+
+const size_t lines_count = sizeof(lines) / sizeof(lines[0]);
+
+// 'not a line' tries to fool the line detection code
+const Cubic notLines[] = {
+    {{0, 0}, {0, 0}, {0, 1}, {1, 0}},
+    {{0, 0}, {0, 1}, {0, 0}, {1, 0}},
+    {{0, 0}, {0, 1}, {1, 0}, {0, 0}},
+    {{0, 1}, {0, 0}, {0, 0}, {1, 0}},
+    {{0, 1}, {0, 0}, {1, 0}, {0, 0}},
+    {{0, 1}, {1, 0}, {0, 0}, {0, 0}},
+};
+
+const size_t notLines_count = sizeof(notLines) / sizeof(notLines[0]);
+
+static const double E = PointEpsilon * 2;
+static const double F = PointEpsilon * 3;
+static const double H = PointEpsilon * 4;
+static const double J = PointEpsilon * 5;
+static const double K = PointEpsilon * 8; // INVESTIGATE: why are larger multiples necessary?
+
+const Cubic modEpsilonLines[] = {
+    {{0, E}, {0, 0}, {0, 0}, {1, 0}}, // horizontal
+    {{0, 0}, {0, E}, {1, 0}, {0, 0}},
+    {{0, 0}, {1, 0}, {0, E}, {0, 0}},
+    {{1, 0}, {0, 0}, {0, 0}, {0, E}},
+    {{1, E}, {2, 0}, {3, 0}, {4, 0}},
+    {{E, 0}, {0, 0}, {0, 0}, {0, 1}}, // vertical
+    {{0, 0}, {E, 0}, {0, 1}, {0, 0}},
+    {{0, 0}, {0, 1}, {E, 0}, {0, 0}},
+    {{0, 1}, {0, 0}, {0, 0}, {E, 0}},
+    {{E, 1}, {0, 2}, {0, 3}, {0, 4}},
+    {{E, 0}, {0, 0}, {0, 0}, {1, 1}}, // 3 coincident
+    {{0, 0}, {E, 0}, {1, 1}, {0, 0}},
+    {{0, 0}, {1, 1}, {E, 0}, {0, 0}},
+    {{1, 1}, {0, 0}, {0, 0}, {E, 0}},
+    {{0, E}, {0, 0}, {1, 1}, {2, 2}}, // 2 coincident
+    {{0, 0}, {1, 1}, {0, E}, {2, 2}},
+    {{0, 0}, {1, 1}, {2, 2}, {0, E}},
+    {{1, 1}, {0, E}, {0, 0}, {2, 2}},
+    {{1, 1}, {0, E}, {2, 2}, {0, 0}},
+    {{1, 1}, {2, 2}, {E, 0}, {0, 0}},
+    {{1, 1}, {2, 2+E}, {3, 3}, {2, 2}}, // middle-last coincident
+    {{1, 1}, {2+E, 2}, {3, 3}, {3, 3}}, // middle-last coincident
+    {{1, 1}, {1, 1}, {2, 2}, {2+E, 2}}, // 2 pairs coincident
+    {{1, 1}, {2, 2}, {1, 1}, {2+E, 2}},
+    {{1, 1}, {2, 2}, {2, 2+E}, {1, 1}},
+    {{1, 1}, {1, 1+E}, {3, 3}, {3, 3}}, // first-middle middle-last coincident
+    {{1, 1}, {2+E, 2}, {3, 3}, {4, 4}}, // no coincident
+    {{1, 1}, {3, 3}, {2, 2}, {4, 4+F}}, // INVESTIGATE: why the epsilon is bigger
+    {{1, 1+F}, {2, 2}, {4, 4}, {3, 3}}, // INVESTIGATE: why the epsilon is bigger
+    {{1, 1}, {3, 3}, {4, 4+E}, {2, 2}},
+    {{1, 1}, {4, 4}, {2, 2}, {3, 3+E}},
+    {{1, 1}, {4, 4}, {3, 3}, {2+E, 2}},
+    {{2, 2}, {1, 1}, {3+E, 3}, {4, 4}},
+    {{2, 2}, {1+E, 1}, {4, 4}, {3, 3}},
+    {{2, 2+E}, {3, 3}, {1, 1}, {4, 4}},
+    {{2+E, 2}, {3, 3}, {4, 4}, {1, 1}},
+    {{2, 2}, {4+E, 4}, {1, 1}, {3, 3}},
+    {{2, 2}, {4, 4}, {3, 3}, {1, 1+E}},
+};
+
+const size_t modEpsilonLines_count = sizeof(modEpsilonLines) / sizeof(modEpsilonLines[0]);
+
+static const double D = PointEpsilon / 2;
+static const double G = PointEpsilon / 3;
+
+const Cubic lessEpsilonLines[] = {
+    {{0, D}, {0, 0}, {0, 0}, {1, 0}}, // horizontal
+    {{0, 0}, {0, D}, {1, 0}, {0, 0}},
+    {{0, 0}, {1, 0}, {0, D}, {0, 0}},
+    {{1, 0}, {0, 0}, {0, 0}, {0, D}},
+    {{1, D}, {2, 0}, {3, 0}, {4, 0}},
+    {{D, 0}, {0, 0}, {0, 0}, {0, 1}}, // vertical
+    {{0, 0}, {D, 0}, {0, 1}, {0, 0}},
+    {{0, 0}, {0, 1}, {D, 0}, {0, 0}},
+    {{0, 1}, {0, 0}, {0, 0}, {D, 0}},
+    {{D, 1}, {0, 2}, {0, 3}, {0, 4}},
+    {{D, 0}, {0, 0}, {0, 0}, {1, 1}}, // 3 coincident
+    {{0, 0}, {D, 0}, {1, 1}, {0, 0}},
+    {{0, 0}, {1, 1}, {D, 0}, {0, 0}},
+    {{1, 1}, {0, 0}, {0, 0}, {D, 0}},
+    {{0, D}, {0, 0}, {1, 1}, {2, 2}}, // 2 coincident
+    {{0, 0}, {1, 1}, {0, D}, {2, 2}},
+    {{0, 0}, {1, 1}, {2, 2}, {0, D}},
+    {{1, 1}, {0, D}, {0, 0}, {2, 2}},
+    {{1, 1}, {0, D}, {2, 2}, {0, 0}},
+    {{1, 1}, {2, 2}, {D, 0}, {0, 0}},
+    {{1, 1}, {2, 2+D}, {3, 3}, {2, 2}}, // middle-last coincident
+    {{1, 1}, {2+D, 2}, {3, 3}, {3, 3}}, // middle-last coincident
+    {{1, 1}, {1, 1}, {2, 2}, {2+D, 2}}, // 2 pairs coincident
+    {{1, 1}, {2, 2}, {1, 1}, {2+D, 2}},
+    {{1, 1}, {2, 2}, {2, 2+D}, {1, 1}},
+    {{1, 1}, {1, 1+D}, {3, 3}, {3, 3}}, // first-middle middle-last coincident
+    {{1, 1}, {2+D/2, 2}, {3, 3}, {4, 4}}, // no coincident (FIXME: N as opposed to N/2 failed)
+    {{1, 1}, {3, 3}, {2, 2}, {4, 4+D}},
+    {{1, 1+D}, {2, 2}, {4, 4}, {3, 3}},
+    {{1, 1}, {3, 3}, {4, 4+D}, {2, 2}},
+    {{1, 1}, {4, 4}, {2, 2}, {3, 3+D}},
+    {{1, 1}, {4, 4}, {3, 3}, {2+G, 2}}, // INVESTIGATE: why the epsilon is smaller
+    {{2, 2}, {1, 1}, {3+D, 3}, {4, 4}},
+    {{2, 2}, {1+D, 1}, {4, 4}, {3, 3}},
+    {{2, 2+D}, {3, 3}, {1, 1}, {4, 4}},
+    {{2+G, 2}, {3, 3}, {4, 4}, {1, 1}}, // INVESTIGATE: why the epsilon is smaller
+    {{2, 2}, {4+D, 4}, {1, 1}, {3, 3}},
+    {{2, 2}, {4, 4}, {3, 3}, {1, 1+D}},
+};
+
+const size_t lessEpsilonLines_count = sizeof(lessEpsilonLines) / sizeof(lessEpsilonLines[0]);
+
+static const double N = -PointEpsilon / 2;
+static const double M = -PointEpsilon / 3;
+
+const Cubic negEpsilonLines[] = {
+    {{0, N}, {0, 0}, {0, 0}, {1, 0}}, // horizontal
+    {{0, 0}, {0, N}, {1, 0}, {0, 0}},
+    {{0, 0}, {1, 0}, {0, N}, {0, 0}},
+    {{1, 0}, {0, 0}, {0, 0}, {0, N}},
+    {{1, N}, {2, 0}, {3, 0}, {4, 0}},
+    {{N, 0}, {0, 0}, {0, 0}, {0, 1}}, // vertical
+    {{0, 0}, {N, 0}, {0, 1}, {0, 0}},
+    {{0, 0}, {0, 1}, {N, 0}, {0, 0}},
+    {{0, 1}, {0, 0}, {0, 0}, {N, 0}},
+    {{N, 1}, {0, 2}, {0, 3}, {0, 4}},
+    {{N, 0}, {0, 0}, {0, 0}, {1, 1}}, // 3 coincident
+    {{0, 0}, {N, 0}, {1, 1}, {0, 0}},
+    {{0, 0}, {1, 1}, {N, 0}, {0, 0}},
+    {{1, 1}, {0, 0}, {0, 0}, {N, 0}},
+    {{0, N}, {0, 0}, {1, 1}, {2, 2}}, // 2 coincident
+    {{0, 0}, {1, 1}, {0, N}, {2, 2}},
+    {{0, 0}, {1, 1}, {2, 2}, {0, N}},
+    {{1, 1}, {0, N}, {0, 0}, {2, 2}},
+    {{1, 1}, {0, N}, {2, 2}, {0, 0}},
+    {{1, 1}, {2, 2}, {N, 0}, {0, 0}},
+    {{1, 1}, {2, 2+N}, {3, 3}, {2, 2}}, // middle-last coincident
+    {{1, 1}, {2+N, 2}, {3, 3}, {3, 3}}, // middle-last coincident
+    {{1, 1}, {1, 1}, {2, 2}, {2+N, 2}}, // 2 pairs coincident
+    {{1, 1}, {2, 2}, {1, 1}, {2+N, 2}},
+    {{1, 1}, {2, 2}, {2, 2+N}, {1, 1}},
+    {{1, 1}, {1, 1+N}, {3, 3}, {3, 3}}, // first-middle middle-last coincident
+    {{1, 1}, {2+N/2, 2}, {3, 3}, {4, 4}}, // no coincident (FIXME: N as opposed to N/2 failed)
+    {{1, 1}, {3, 3}, {2, 2}, {4, 4+N}},
+    {{1, 1+N}, {2, 2}, {4, 4}, {3, 3}},
+    {{1, 1}, {3, 3}, {4, 4+N}, {2, 2}},
+    {{1, 1}, {4, 4}, {2, 2}, {3, 3+N}},
+    {{1, 1}, {4, 4}, {3, 3}, {2+M, 2}}, // INVESTIGATE: why the epsilon is smaller
+    {{2, 2}, {1, 1}, {3+N, 3}, {4, 4}},
+    {{2, 2}, {1+N, 1}, {4, 4}, {3, 3}},
+    {{2, 2+N}, {3, 3}, {1, 1}, {4, 4}},
+    {{2+M, 2}, {3, 3}, {4, 4}, {1, 1}}, // INVESTIGATE: why the epsilon is smaller
+    {{2, 2}, {4+N, 4}, {1, 1}, {3, 3}},
+    {{2, 2}, {4, 4}, {3, 3}, {1, 1+N}},
+};
+
+const size_t negEpsilonLines_count = sizeof(negEpsilonLines) / sizeof(negEpsilonLines[0]);
+
+const Quadratic quadraticLines[] = {
+    {{0, 0}, {0, 0}, {1, 0}},
+    {{0, 0}, {1, 0}, {0, 0}},
+    {{1, 0}, {0, 0}, {0, 0}},
+    {{1, 0}, {2, 0}, {3, 0}},
+    {{0, 0}, {0, 0}, {0, 1}},
+    {{0, 0}, {0, 1}, {0, 0}},
+    {{0, 1}, {0, 0}, {0, 0}},
+    {{0, 1}, {0, 2}, {0, 3}},
+    {{0, 0}, {0, 0}, {1, 1}},
+    {{0, 0}, {1, 1}, {0, 0}},
+    {{1, 1}, {0, 0}, {0, 0}},
+    {{1, 1}, {2, 2}, {3, 3}},
+    {{1, 1}, {3, 3}, {3, 3}},
+    {{1, 1}, {1, 1}, {2, 2}},
+    {{1, 1}, {2, 2}, {1, 1}},
+    {{1, 1}, {1, 1}, {3, 3}},
+    {{1, 1}, {2, 2}, {4, 4}}, // no coincident
+    {{1, 1}, {3, 3}, {4, 4}},
+    {{1, 1}, {3, 3}, {2, 2}},
+    {{1, 1}, {4, 4}, {2, 2}},
+    {{1, 1}, {4, 4}, {3, 3}},
+    {{2, 2}, {1, 1}, {3, 3}},
+    {{2, 2}, {1, 1}, {4, 4}},
+    {{2, 2}, {3, 3}, {1, 1}},
+    {{2, 2}, {3, 3}, {4, 4}},
+    {{2, 2}, {4, 4}, {1, 1}},
+    {{2, 2}, {4, 4}, {3, 3}},
+};
+
+const size_t quadraticLines_count = sizeof(quadraticLines) / sizeof(quadraticLines[0]);
+
+const Quadratic quadraticModEpsilonLines[] = {
+    {{0, F}, {0, 0}, {1, 0}},
+    {{0, 0}, {1, 0}, {0, F}},
+    {{1, 0}, {0, F}, {0, 0}},
+    {{1, H}, {2, 0}, {3, 0}},
+    {{F, 0}, {0, 0}, {0, 1}},
+    {{0, 0}, {0, 1}, {F, 0}},
+    {{0, 1}, {F, 0}, {0, 0}},
+    {{H, 1}, {0, 2}, {0, 3}},
+    {{0, F}, {0, 0}, {1, 1}},
+    {{0, 0}, {1, 1}, {F, 0}},
+    {{1, 1}, {F, 0}, {0, 0}},
+    {{1, 1+J}, {2, 2}, {3, 3}},
+    {{1, 1}, {3, 3}, {3+F, 3}},
+    {{1, 1}, {1+F, 1}, {2, 2}},
+    {{1, 1}, {2, 2}, {1, 1+F}},
+    {{1, 1}, {1, 1+F}, {3, 3}},
+    {{1+H, 1}, {2, 2}, {4, 4}}, // no coincident
+    {{1, 1+K}, {3, 3}, {4, 4}},
+    {{1, 1}, {3+F, 3}, {2, 2}},
+    {{1, 1}, {4, 4+F}, {2, 2}},
+    {{1, 1}, {4, 4}, {3+F, 3}},
+    {{2, 2}, {1, 1}, {3, 3+F}},
+    {{2+F, 2}, {1, 1}, {4, 4}},
+    {{2, 2+F}, {3, 3}, {1, 1}},
+    {{2, 2}, {3+F, 3}, {4, 4}},
+    {{2, 2}, {4, 4+F}, {1, 1}},
+    {{2, 2}, {4, 4}, {3+F, 3}},
+};
+
+const size_t quadraticModEpsilonLines_count = sizeof(quadraticModEpsilonLines) / sizeof(quadraticModEpsilonLines[0]);
+
+
diff --git a/experimental/Intersection/CubicIntersection_TestData.h b/experimental/Intersection/CubicIntersection_TestData.h
new file mode 100644
index 0000000..931a0c7
--- /dev/null
+++ b/experimental/Intersection/CubicIntersection_TestData.h
@@ -0,0 +1,28 @@
+#include "DataTypes.h"
+
+extern const Cubic pointDegenerates[];
+extern const Cubic notPointDegenerates[];
+extern const Cubic tests[][2];
+extern Cubic hexTests[][2];
+
+extern void convert_testx();
+
+extern const Cubic lines[];
+extern const Cubic notLines[];
+extern const Cubic modEpsilonLines[];
+extern const Cubic lessEpsilonLines[];
+extern const Cubic negEpsilonLines[];
+extern const Quadratic quadraticLines[];
+extern const Quadratic quadraticModEpsilonLines[];
+
+extern const size_t pointDegenerates_count;
+extern const size_t notPointDegenerates_count;
+extern const size_t tests_count;
+extern const size_t hexTests_count;
+extern const size_t lines_count;
+extern const size_t notLines_count;
+extern const size_t modEpsilonLines_count;
+extern const size_t lessEpsilonLines_count;
+extern const size_t negEpsilonLines_count;
+extern const size_t quadraticLines_count;
+extern const size_t quadraticModEpsilonLines_count;
diff --git a/experimental/Intersection/CubicIntersection_Tests.cpp b/experimental/Intersection/CubicIntersection_Tests.cpp
new file mode 100644
index 0000000..14e5821
--- /dev/null
+++ b/experimental/Intersection/CubicIntersection_Tests.cpp
@@ -0,0 +1,15 @@
+#include "CubicIntersection_Tests.h"
+
+void CubicIntersection_Tests() {
+    // tests are in dependency order
+    Inline_Tests();
+    ConvexHull_Test();
+    ConvexHull_X_Test();
+    LineParameter_Test();
+    LineIntersection_Test();
+    QuadraticCoincidence_Test();
+    CubicCoincidence_Test();
+    ReduceOrder_Test();
+ //   BezierClip_Test();
+    CubicIntersection_Test();
+}
diff --git a/experimental/Intersection/CubicIntersection_Tests.h b/experimental/Intersection/CubicIntersection_Tests.h
new file mode 100644
index 0000000..7c85284
--- /dev/null
+++ b/experimental/Intersection/CubicIntersection_Tests.h
@@ -0,0 +1,11 @@
+void BezierClip_Test();
+void ConvexHull_Test();
+void ConvexHull_X_Test();
+void CubicCoincidence_Test();
+void CubicIntersection_Test();
+void CubicIntersection_Tests();
+void Inline_Tests();
+void LineIntersection_Test();
+void LineParameter_Test();
+void QuadraticCoincidence_Test();
+void ReduceOrder_Test();
diff --git a/experimental/Intersection/CubicParameterization.cpp b/experimental/Intersection/CubicParameterization.cpp
new file mode 100644
index 0000000..9f45ba2
--- /dev/null
+++ b/experimental/Intersection/CubicParameterization.cpp
@@ -0,0 +1,306 @@
+#include "CubicIntersection.h"
+
+/* from http://tom.cs.byu.edu/~tom/papers/cvgip84.pdf 4.1
+ *
+ * This paper proves that Syvester's method can compute the implicit form of 
+ * the quadratic from the parameterzied form.
+ *
+ * Given x = a*t*t*t + b*t*t + c*t + d  (the parameterized form)
+ *       y = e*t*t*t + f*t*t + g*t + h
+ *
+ * we want to find an equation of the implicit form:
+ *
+ * A*x^3 + B*x*x*y + C*x*y*y + D*y^3 + E*x*x + F*x*y + G*y*y + H*x + I*y + J = 0
+ *
+ * The implicit form can be expressed as a 6x6 determinant, as shown.
+ *
+ * The resultant obtained by Syvester's method is
+ *
+ * |   a   b   c  (d - x)     0        0     |
+ * |   0   a   b     c     (d - x)     0     |
+ * |   0   0   a     b        c     (d - x)  |
+ * |   e   f   g  (h - y)     0        0     |
+ * |   0   e   f     g     (h - y)     0     |
+ * |   0   0   e     f        g     (h - y)  |
+ *
+ * which, according to Mathematica, expands as shown below.
+ *
+ * Resultant[a*t^3 + b*t^2 + c*t + d - x, e*t^3 + f*t^2 + g*t + h - y, t]
+ *
+ *  -d^3 e^3 + c d^2 e^2 f - b d^2 e f^2 + a d^2 f^3 - c^2 d e^2 g + 
+ *  2 b d^2 e^2 g + b c d e f g - 3 a d^2 e f g - a c d f^2 g - 
+ *  b^2 d e g^2 + 2 a c d e g^2 + a b d f g^2 - a^2 d g^3 + c^3 e^2 h - 
+ *  3 b c d e^2 h + 3 a d^2 e^2 h - b c^2 e f h + 2 b^2 d e f h + 
+ *  a c d e f h + a c^2 f^2 h - 2 a b d f^2 h + b^2 c e g h - 
+ *  2 a c^2 e g h - a b d e g h - a b c f g h + 3 a^2 d f g h + 
+ *  a^2 c g^2 h - b^3 e h^2 + 3 a b c e h^2 - 3 a^2 d e h^2 + 
+ *  a b^2 f h^2 - 2 a^2 c f h^2 - a^2 b g h^2 + a^3 h^3 + 3 d^2 e^3 x - 
+ *  2 c d e^2 f x + 2 b d e f^2 x - 2 a d f^3 x + c^2 e^2 g x - 
+ *  4 b d e^2 g x - b c e f g x + 6 a d e f g x + a c f^2 g x + 
+ *  b^2 e g^2 x - 2 a c e g^2 x - a b f g^2 x + a^2 g^3 x + 
+ *  3 b c e^2 h x - 6 a d e^2 h x - 2 b^2 e f h x - a c e f h x + 
+ *  2 a b f^2 h x + a b e g h x - 3 a^2 f g h x + 3 a^2 e h^2 x - 
+ *  3 d e^3 x^2 + c e^2 f x^2 - b e f^2 x^2 + a f^3 x^2 + 
+ *  2 b e^2 g x^2 - 3 a e f g x^2 + 3 a e^2 h x^2 + e^3 x^3 - 
+ *  c^3 e^2 y + 3 b c d e^2 y - 3 a d^2 e^2 y + b c^2 e f y - 
+ *  2 b^2 d e f y - a c d e f y - a c^2 f^2 y + 2 a b d f^2 y - 
+ *  b^2 c e g y + 2 a c^2 e g y + a b d e g y + a b c f g y - 
+ *  3 a^2 d f g y - a^2 c g^2 y + 2 b^3 e h y - 6 a b c e h y + 
+ *  6 a^2 d e h y - 2 a b^2 f h y + 4 a^2 c f h y + 2 a^2 b g h y - 
+ *  3 a^3 h^2 y - 3 b c e^2 x y + 6 a d e^2 x y + 2 b^2 e f x y + 
+ *  a c e f x y - 2 a b f^2 x y - a b e g x y + 3 a^2 f g x y - 
+ *  6 a^2 e h x y - 3 a e^2 x^2 y - b^3 e y^2 + 3 a b c e y^2 - 
+ *  3 a^2 d e y^2 + a b^2 f y^2 - 2 a^2 c f y^2 - a^2 b g y^2 + 
+ *  3 a^3 h y^2 + 3 a^2 e x y^2 - a^3 y^3
+ */
+
+enum {
+    xxx_coeff,
+    xxy_coeff,
+    xyy_coeff,
+    yyy_coeff,
+    xx_coeff,
+    xy_coeff,
+    yy_coeff,
+    x_coeff,
+    y_coeff,
+    c_coeff,
+    coeff_count
+};
+
+// FIXME: factoring version unwritten
+// static bool straight_forward = true;
+
+/* from CubicParameterizationCode.cpp output: 
+ *  double A =      e * e * e;
+ *  double B = -3 * a * e * e;
+ *  double C =  3 * a * a * e;
+ *  double D =     -a * a * a;
+ */
+static void calc_ABCD(double a, double e, double p[coeff_count]) {
+    double ee = e * e;
+    p[xxx_coeff] = e * ee;
+    p[xxy_coeff] = -3 * a * ee;
+    double aa = a * a;
+    p[xyy_coeff] = 3 * aa * e;
+    p[yyy_coeff] = -aa * a;
+}
+
+/* CubicParameterizationCode.cpp turns Mathematica output into C.
+ * Rather than edit the lines below, please edit the code there instead.
+ */
+// start of generated code
+static double calc_E(double a, double b, double c, double d,
+                     double e, double f, double g, double h) {
+    return 
+         -3 * d * e * e * e
+        +     c * e * e * f
+        -     b * e * f * f
+        +     a * f * f * f
+        + 2 * b * e * e * g
+        - 3 * a * e * f * g
+        + 3 * a * e * e * h;
+}
+
+static double calc_F(double a, double b, double c, double d,
+                     double e, double f, double g, double h) {
+    return 
+         -3 * b * c * e * e
+        + 6 * a * d * e * e
+        + 2 * b * b * e * f
+        +     a * c * e * f
+        - 2 * a * b * f * f
+        -     a * b * e * g
+        + 3 * a * a * f * g
+        - 6 * a * a * e * h;
+}
+
+static double calc_G(double a, double b, double c, double d,
+                     double e, double f, double g, double h) {
+    return 
+             -b * b * b * e
+        + 3 * a * b * c * e
+        - 3 * a * a * d * e
+        +     a * b * b * f
+        - 2 * a * a * c * f
+        -     a * a * b * g
+        + 3 * a * a * a * h;
+}
+
+static double calc_H(double a, double b, double c, double d,
+                     double e, double f, double g, double h) {
+    return 
+          3 * d * d * e * e * e
+        - 2 * c * d * e * e * f
+        + 2 * b * d * e * f * f
+        - 2 * a * d * f * f * f
+        +     c * c * e * e * g
+        - 4 * b * d * e * e * g
+        -     b * c * e * f * g
+        + 6 * a * d * e * f * g
+        +     a * c * f * f * g
+        +     b * b * e * g * g
+        - 2 * a * c * e * g * g
+        -     a * b * f * g * g
+        +     a * a * g * g * g
+        + 3 * b * c * e * e * h
+        - 6 * a * d * e * e * h
+        - 2 * b * b * e * f * h
+        -     a * c * e * f * h
+        + 2 * a * b * f * f * h
+        +     a * b * e * g * h
+        - 3 * a * a * f * g * h
+        + 3 * a * a * e * h * h;
+}
+
+static double calc_I(double a, double b, double c, double d,
+                     double e, double f, double g, double h) {
+    return 
+             -c * c * c * e * e
+        + 3 * b * c * d * e * e
+        - 3 * a * d * d * e * e
+        +     b * c * c * e * f
+        - 2 * b * b * d * e * f
+        -     a * c * d * e * f
+        -     a * c * c * f * f
+        + 2 * a * b * d * f * f
+        -     b * b * c * e * g
+        + 2 * a * c * c * e * g
+        +     a * b * d * e * g
+        +     a * b * c * f * g
+        - 3 * a * a * d * f * g
+        -     a * a * c * g * g
+        + 2 * b * b * b * e * h
+        - 6 * a * b * c * e * h
+        + 6 * a * a * d * e * h
+        - 2 * a * b * b * f * h
+        + 4 * a * a * c * f * h
+        + 2 * a * a * b * g * h
+        - 3 * a * a * a * h * h;
+}
+
+static double calc_J(double a, double b, double c, double d,
+                     double e, double f, double g, double h) {
+    return 
+             -d * d * d * e * e * e
+        +     c * d * d * e * e * f
+        -     b * d * d * e * f * f
+        +     a * d * d * f * f * f
+        -     c * c * d * e * e * g
+        + 2 * b * d * d * e * e * g
+        +     b * c * d * e * f * g
+        - 3 * a * d * d * e * f * g
+        -     a * c * d * f * f * g
+        -     b * b * d * e * g * g
+        + 2 * a * c * d * e * g * g
+        +     a * b * d * f * g * g
+        -     a * a * d * g * g * g
+        +     c * c * c * e * e * h
+        - 3 * b * c * d * e * e * h
+        + 3 * a * d * d * e * e * h
+        -     b * c * c * e * f * h
+        + 2 * b * b * d * e * f * h
+        +     a * c * d * e * f * h
+        +     a * c * c * f * f * h
+        - 2 * a * b * d * f * f * h
+        +     b * b * c * e * g * h
+        - 2 * a * c * c * e * g * h
+        -     a * b * d * e * g * h
+        -     a * b * c * f * g * h
+        + 3 * a * a * d * f * g * h
+        +     a * a * c * g * g * h
+        -     b * b * b * e * h * h
+        + 3 * a * b * c * e * h * h
+        - 3 * a * a * d * e * h * h
+        +     a * b * b * f * h * h
+        - 2 * a * a * c * f * h * h
+        -     a * a * b * g * h * h
+        +     a * a * a * h * h * h;
+}
+// end of generated code
+
+static double (*calc_proc[])(double a, double b, double c, double d,
+                             double e, double f, double g, double h) = {
+    calc_E, calc_F, calc_G, calc_H, calc_I, calc_J
+};
+
+/* Control points to parametric coefficients
+    s = 1 - t
+    Attt + 3Btt2 + 3Ctss + Dsss ==
+    Attt + 3B(1 - t)tt + 3C(1 - t)(t - tt) + D(1 - t)(1 - 2t + tt) ==
+    Attt + 3B(tt - ttt) + 3C(t - tt - tt + ttt) + D(1-2t+tt-t+2tt-ttt) ==
+    Attt + 3Btt - 3Bttt + 3Ct - 6Ctt + 3Cttt + D - 3Dt + 3Dtt - Dttt ==
+    D + (3C - 3D)t + (3B - 6C + 3D)tt + (A - 3B + 3C - D)ttt
+    a = A - 3*B + 3*C -   D
+    b =     3*B - 6*C + 3*D
+    c =           3*C - 3*D
+    d =                   D
+ */
+static void set_abcd(const double* cubic, double& a, double& b, double& c,
+                     double& d) {
+    a = cubic[0];     // a = A
+    b = 3 * cubic[2]; // b = 3*B (compute rest of b lazily)
+    c = 3 * cubic[4]; // c = 3*C (compute rest of c lazily)
+    d = cubic[6];     // d = D
+    a += -b + c - d;  // a = A - 3*B + 3*C - D
+}
+
+static void calc_bc(const double d, double& b, double& c) {
+    b -= 3 * c; // b = 3*B - 3*C
+    c -= 3 * d; // c = 3*C - 3*D
+    b -= c;     // b = 3*B - 6*C + 3*D
+}
+
+bool implicit_matches(const Cubic& one, const Cubic& two) {
+    double p1[coeff_count]; // a'xxx , b'xxy , c'xyy , d'xx , e'xy , f'yy, etc.
+    double p2[coeff_count];
+    double a1, b1, c1, d1;
+    set_abcd(&one[0].x, a1, b1, c1, d1);
+    double e1, f1, g1, h1;
+    set_abcd(&one[0].y, e1, f1, g1, h1);
+    calc_ABCD(a1, e1, p1);
+    double a2, b2, c2, d2;
+    set_abcd(&two[0].x, a2, b2, c2, d2);
+    double e2, f2, g2, h2;
+    set_abcd(&two[0].y, e2, f2, g2, h2);
+    calc_ABCD(a2, e2, p2);
+    int first = 0;
+    for (int index = 0; index < coeff_count; ++index) {
+        if (index == xx_coeff) {
+            calc_bc(d1, b1, c1);
+            calc_bc(h1, f1, g1);
+            calc_bc(d2, b2, c2);
+            calc_bc(h2, f2, g2);
+        }
+        if (index >= xx_coeff) {
+            int procIndex = index - xx_coeff;
+            p1[index] = (*calc_proc[procIndex])(a1, b1, c1, d1, e1, f1, g1, h1);
+            p2[index] = (*calc_proc[procIndex])(a2, b2, c2, d2, e2, f2, g2, h2);
+        }
+        if (approximately_zero(p1[index]) || approximately_zero(p2[index])) {
+            first += first == index;
+            continue;
+        }
+        if (first == index) {
+            continue;
+        }
+        if (!approximately_equal(p1[index] * p2[first],
+                p1[first] * p2[index])) {
+            return false;
+        }
+    }
+    return true;
+}
+
+static double tangent(const double* cubic, double t) {
+    double a, b, c, d;
+    set_abcd(cubic, a, b, c, d);
+    calc_bc(d, b, c);
+    return 3 * a * t * t + 2 * b * t + c;
+}
+
+void tangent(const Cubic& cubic, double t, _Point& result) {
+    result.x = tangent(&cubic[0].x, t);
+    result.y = tangent(&cubic[0].y, t);
+}
+
diff --git a/experimental/Intersection/CubicParameterizationCode.cpp b/experimental/Intersection/CubicParameterizationCode.cpp
new file mode 100644
index 0000000..9e7a3f6
--- /dev/null
+++ b/experimental/Intersection/CubicParameterizationCode.cpp
@@ -0,0 +1,282 @@
+#include <vector>
+
+/* Given:
+ * Resultant[a*t^3 + b*t^2 + c*t + d - x, e*t^3 + f*t^2 + g*t + h - y, t]
+ */
+
+const char result[] =
+"-d^3 e^3 + c d^2 e^2 f - b d^2 e f^2 + a d^2 f^3 - c^2 d e^2 g + "
+" 2 b d^2 e^2 g + b c d e f g - 3 a d^2 e f g - a c d f^2 g - "
+" b^2 d e g^2 + 2 a c d e g^2 + a b d f g^2 - a^2 d g^3 + c^3 e^2 h - "
+" 3 b c d e^2 h + 3 a d^2 e^2 h - b c^2 e f h + 2 b^2 d e f h + "
+" a c d e f h + a c^2 f^2 h - 2 a b d f^2 h + b^2 c e g h - "
+" 2 a c^2 e g h - a b d e g h - a b c f g h + 3 a^2 d f g h + "
+" a^2 c g^2 h - b^3 e h^2 + 3 a b c e h^2 - 3 a^2 d e h^2 + "
+" a b^2 f h^2 - 2 a^2 c f h^2 - a^2 b g h^2 + a^3 h^3 + 3 d^2 e^3 x - "
+" 2 c d e^2 f x + 2 b d e f^2 x - 2 a d f^3 x + c^2 e^2 g x - "
+" 4 b d e^2 g x - b c e f g x + 6 a d e f g x + a c f^2 g x + "
+" b^2 e g^2 x - 2 a c e g^2 x - a b f g^2 x + a^2 g^3 x + "
+" 3 b c e^2 h x - 6 a d e^2 h x - 2 b^2 e f h x - a c e f h x + "
+" 2 a b f^2 h x + a b e g h x - 3 a^2 f g h x + 3 a^2 e h^2 x - "
+" 3 d e^3 x^2 + c e^2 f x^2 - b e f^2 x^2 + a f^3 x^2 + "
+" 2 b e^2 g x^2 - 3 a e f g x^2 + 3 a e^2 h x^2 + e^3 x^3 - "
+" c^3 e^2 y + 3 b c d e^2 y - 3 a d^2 e^2 y + b c^2 e f y - "
+" 2 b^2 d e f y - a c d e f y - a c^2 f^2 y + 2 a b d f^2 y - "
+" b^2 c e g y + 2 a c^2 e g y + a b d e g y + a b c f g y - "
+" 3 a^2 d f g y - a^2 c g^2 y + 2 b^3 e h y - 6 a b c e h y + "
+" 6 a^2 d e h y - 2 a b^2 f h y + 4 a^2 c f h y + 2 a^2 b g h y - "
+" 3 a^3 h^2 y - 3 b c e^2 x y + 6 a d e^2 x y + 2 b^2 e f x y + "
+" a c e f x y - 2 a b f^2 x y - a b e g x y + 3 a^2 f g x y - "
+" 6 a^2 e h x y - 3 a e^2 x^2 y - b^3 e y^2 + 3 a b c e y^2 - "
+" 3 a^2 d e y^2 + a b^2 f y^2 - 2 a^2 c f y^2 - a^2 b g y^2 + "
+" 3 a^3 h y^2 + 3 a^2 e x y^2 - a^3 y^3";
+
+const size_t len = sizeof(result) - 1;
+const int factors = 8;
+
+struct coeff {
+    int s; // constant and coefficient sign
+    int n[factors]; // 0 or power of a (1, 2, or 3) for a through h
+};
+
+enum {
+    xxx_coeff,
+    xxy_coeff,
+    xyy_coeff,
+    yyy_coeff,
+    xx_coeff,
+    xy_coeff,
+    yy_coeff,
+    x_coeff,
+    y_coeff,
+    c_coeff,
+    coeff_count
+};
+
+typedef std::vector<coeff> coeffs;
+typedef std::vector<coeffs> n_coeffs;
+
+static char skipSpace(size_t& index) {
+    do {
+        ++index;
+    } while (result[index] == ' ');
+    return result[index];
+}
+
+static char backSkipSpace(size_t& end) {
+    while (result[end - 1] == ' ') {
+        --end;
+    }
+    return result[end - 1];
+}
+
+static void match(coeffs& co, const char pattern[]) {
+    size_t patternLen = strlen(pattern);
+    size_t index = 0;
+    while (index < len) {
+        char ch = result[index];
+        if (ch != '-' && ch != '+') {
+            printf("missing sign\n");
+        }
+        size_t end = index + 1;
+        while (result[end] != '+' && result[end] != '-' && ++end < len) {
+            ;
+        }
+        backSkipSpace(end);
+        size_t idx = index;
+        index = end;
+        skipSpace(index);
+        if (!strncmp(&result[end - patternLen], pattern, patternLen) == 0) {
+            continue;
+        }
+        size_t endCoeff = end - patternLen;
+        char last = backSkipSpace(endCoeff);
+        if (last == '2' || last == '3') {
+            last = result[endCoeff - 3]; // skip ^2
+        }
+        if (last == 'x' || last == 'y') {
+            continue;
+        }
+        coeff c;
+        c.s = result[idx] == '-' ? -1 : 1;
+        bzero(c.n, sizeof(c.n));
+        ch = skipSpace(idx);
+        if (ch >= '2' && ch <= '6') {
+            c.s *= ch - '0';
+            ch = skipSpace(idx);
+        }
+        while (idx < endCoeff) {
+            char x = result[idx];
+            if (x < 'a' || x > 'a' + factors) {
+                printf("expected factor\n");
+            }
+            idx++;
+            int pow = 1;
+            if (result[idx] == '^') {
+                idx++;
+                char exp = result[idx];
+                if (exp < '2' || exp > '3') {
+                    printf("expected exponent\n");
+                }
+                pow = exp - '0';
+            }
+            skipSpace(idx);
+            c.n[x - 'a'] = pow;
+        }
+        co.push_back(c);
+    }
+}
+
+void cubecode_test();
+
+void cubecode_test() {
+    n_coeffs c(coeff_count);
+    match(c[xxx_coeff], "x^3");   // 1 factor
+    match(c[xxy_coeff], "x^2 y"); // 1 factor
+    match(c[xyy_coeff], "x y^2"); // 1 factor
+    match(c[yyy_coeff], "y^3");   // 1 factor
+    match(c[xx_coeff], "x^2");    // 7 factors
+    match(c[xy_coeff], "x y");    // 8 factors
+    match(c[yy_coeff], "y^2");    // 7 factors
+    match(c[x_coeff], "x");       // 21 factors
+    match(c[y_coeff], "y");       // 21 factors
+    match(c[c_coeff], "");        // 34 factors
+#define COMPUTE_MOST_FREQUENT_EXPRESSION_TRIPLETS 0
+#define WRITE_AS_NONOPTIMIZED_C_CODE 0
+#if COMPUTE_MOST_FREQUENT_EXPRESSION_TRIPLETS
+    int count[factors][factors][factors];
+    bzero(count, sizeof(count));
+#endif
+#if WRITE_AS_NONOPTIMIZED_C_CODE
+    printf("// start of generated code");
+#endif
+    for (n_coeffs::iterator it = c.begin(); it < c.end(); ++it) {
+        coeffs& co = *it;
+#if WRITE_AS_NONOPTIMIZED_C_CODE
+        printf("\nstatic double calc_%c(double a, double b, double c, double d,"
+               "\n                     double e, double f, double g, double h) {"
+               "\n    return"
+               "\n ", 'A' + (it - c.begin()));
+        if (co[0].s > 0) {
+            printf(" ");
+        } 
+        if (abs(co[0].s) == 1) {
+            printf("    ");
+        }
+#endif
+        for (coeffs::iterator ct = co.begin(); ct < co.end(); ++ct) {
+            const coeff& cf = *ct;
+#if WRITE_AS_NONOPTIMIZED_C_CODE
+            printf("        ");
+            bool firstFactor = false;
+            if (ct - co.begin() > 0 || cf.s < 0) {
+                printf("%c", cf.s < 0 ? '-' : '+');
+            }
+            if (ct - co.begin() > 0) {
+                printf(" ");
+            }
+            if (abs(cf.s) > 1) {
+                printf("%d * ", abs(cf.s));
+            } else {
+                if (ct - co.begin() > 0) {
+                    printf("    ");
+                }
+            }
+#endif
+            for (int x = 0; x < factors; ++x) {
+                if (cf.n[x] == 0) {
+                    continue;
+                }
+#if WRITE_AS_NONOPTIMIZED_C_CODE
+                for (int y = 0 ; y < cf.n[x]; ++y) {
+                    if (y > 0 || firstFactor) {
+                        printf(" * ");
+                    }
+                    printf("%c", 'a' + x);
+                }
+                firstFactor = true;
+#endif
+#if COMPUTE_MOST_FREQUENT_EXPRESSION_TRIPLETS
+                for (int y = x; y < factors; ++y) {
+                    if (cf.n[y] == 0) {
+                        continue;
+                    }
+                    if (x == y && cf.n[y] == 1) {
+                        continue;
+                    }
+                    for (int z = y; z < factors; ++z) {
+                        if (cf.n[z] == 0) {
+                            continue;
+                        }
+                        if ((x == z || y == z) && cf.n[z] == 1) {
+                            continue;
+                        }
+                        if (x == y && y == z && cf.n[z] == 2) {
+                            continue;
+                        }
+                        count[x][y][z]++;
+                    }
+                }
+#endif
+            }
+#if WRITE_AS_NONOPTIMIZED_C_CODE
+            if (ct + 1 < co.end()) {
+                printf("\n");
+            }
+#endif
+        }
+#if WRITE_AS_NONOPTIMIZED_C_CODE
+            printf(";\n}\n");
+#endif
+    }
+#if WRITE_AS_NONOPTIMIZED_C_CODE
+    printf("// end of generated code\n");
+#endif
+#if COMPUTE_MOST_FREQUENT_EXPRESSION_TRIPLETS
+    const int bestCount = 20;
+    int best[bestCount][4];
+    bzero(best, sizeof(best));
+    for (int x = 0; x < factors; ++x) {
+        for (int y = x; y < factors; ++y) {
+            for (int z = y; z < factors; ++z) {
+                if (!count[x][y][z]) {
+                    continue;
+                }
+                for (int w = 0; w < bestCount; ++w) {
+                    if (best[w][0] < count[x][y][z]) {
+                        best[w][0] = count[x][y][z];
+                        best[w][1] = x;
+                        best[w][2] = y;
+                        best[w][3] = z;
+                        break;
+                    }
+                }
+            }
+        }
+    }
+    for (int w = 0; w < bestCount; ++w) {
+        printf("%c%c%c=%d\n", 'a' + best[w][1], 'a'  + best[w][2],
+            'a' + best[w][3], best[w][0]);
+    }
+#endif
+#if WRITE_AS_NONOPTIMIZED_C_CODE
+    printf("\n");
+#endif
+}
+
+/* results: variable triplets used 10 or more times:
+aah=14
+ade=14
+aeh=14
+dee=14
+bce=13
+beg=13
+beh=12
+bbe=11
+bef=11
+cee=11
+cef=11
+def=11
+ceh=10
+deg=10
+*/
diff --git a/experimental/Intersection/CubicParameterization_Test.cpp b/experimental/Intersection/CubicParameterization_Test.cpp
new file mode 100644
index 0000000..9bb2253
--- /dev/null
+++ b/experimental/Intersection/CubicParameterization_Test.cpp
@@ -0,0 +1,44 @@
+#include "CubicIntersection.h"
+#include "CubicIntersection_Tests.h"
+#include "TestUtilities.h"
+
+const Quadratic quadratics[] = {
+    {{0, 0}, {1, 0}, {1, 1}},
+};
+
+const size_t quadratics_count = sizeof(quadratics) / sizeof(quadratics[0]);
+
+int firstCubicCoincidenceTest = 0;
+
+void CubicCoincidence_Test() {
+    // split large quadratic
+    // upscale quadratics to cubics
+    // compare original, parts, to see if the are coincident
+    for (size_t index = firstCubicCoincidenceTest; index < quadratics_count; ++index) {
+        const Quadratic& test = quadratics[index];
+        QuadraticPair split;
+        chop_at(test, split, 0.5);
+        Quadratic midThird;
+        sub_divide(test, 1.0/3, 2.0/3, midThird);
+        Cubic whole, first, second, mid;
+        quad_to_cubic(test, whole);
+        quad_to_cubic(split.first(), first);
+        quad_to_cubic(split.second(), second);
+        quad_to_cubic(midThird, mid);
+        if (!implicit_matches(whole, first)) {
+            printf("%s-1 %d\n", __FUNCTION__, (int)index);
+        }
+        if (!implicit_matches(whole, second)) {
+            printf("%s-2 %d\n", __FUNCTION__, (int)index);
+        }
+        if (!implicit_matches(mid, first)) {
+            printf("%s-3 %d\n", __FUNCTION__, (int)index);
+        }
+        if (!implicit_matches(mid, second)) {
+            printf("%s-4 %d\n", __FUNCTION__, (int)index);
+        }
+        if (!implicit_matches(first, second)) {
+            printf("%s-5 %d\n", __FUNCTION__, (int)index);
+        }
+    }
+}
diff --git a/experimental/Intersection/CubicReduceOrder.cpp b/experimental/Intersection/CubicReduceOrder.cpp
new file mode 100644
index 0000000..b7e0ea9
--- /dev/null
+++ b/experimental/Intersection/CubicReduceOrder.cpp
@@ -0,0 +1,229 @@
+#include "CubicIntersection.h"
+#include "Extrema.h"
+#include "IntersectionUtilities.h"
+#include "LineParameters.h"
+
+static double interp_cubic_coords(const double* src, double t)
+{
+    double ab = interp(src[0], src[2], t);
+    double bc = interp(src[2], src[4], t);
+    double cd = interp(src[4], src[6], t);
+    double abc = interp(ab, bc, t);
+    double bcd = interp(bc, cd, t);
+    return interp(abc, bcd, t);
+}
+
+static int coincident_line(const Cubic& cubic, Cubic& reduction) {
+    reduction[0] = reduction[1] = cubic[0];
+    return 1;
+}
+
+static int vertical_line(const Cubic& cubic, Cubic& reduction) {
+    double tValues[2];
+    reduction[0] = cubic[0];
+    reduction[1] = cubic[3];
+    int smaller = reduction[1].y > reduction[0].y;
+    int larger = smaller ^ 1;
+    int roots = SkFindCubicExtrema(cubic[0].y, cubic[1].y, cubic[2].y, cubic[3].y, tValues);
+    for (int index = 0; index < roots; ++index) {
+        double yExtrema = interp_cubic_coords(&cubic[0].y, tValues[index]);
+        if (reduction[smaller].y > yExtrema) {
+            reduction[smaller].y = yExtrema;
+            continue;
+        } 
+        if (reduction[larger].y < yExtrema) {
+            reduction[larger].y = yExtrema;
+        }
+    }
+    return 2;
+}
+
+static int horizontal_line(const Cubic& cubic, Cubic& reduction) {
+    double tValues[2];
+    reduction[0] = cubic[0];
+    reduction[1] = cubic[3];
+    int smaller = reduction[1].x > reduction[0].x;
+    int larger = smaller ^ 1;
+    int roots = SkFindCubicExtrema(cubic[0].x, cubic[1].x, cubic[2].x, cubic[3].x, tValues);
+    for (int index = 0; index < roots; ++index) {
+        double xExtrema = interp_cubic_coords(&cubic[0].x, tValues[index]);
+        if (reduction[smaller].x > xExtrema) {
+            reduction[smaller].x = xExtrema;
+            continue;
+        } 
+        if (reduction[larger].x < xExtrema) {
+            reduction[larger].x = xExtrema;
+        }
+    }
+    return 2;
+}
+
+// check to see if it is a quadratic or a line
+static int check_quadratic(const Cubic& cubic, Cubic& reduction,
+        int minX, int maxX, int minY, int maxY) {
+    double dx10 = cubic[1].x - cubic[0].x;
+    double dx23 = cubic[2].x - cubic[3].x;
+    double midX = cubic[0].x + dx10 * 3 / 2;
+    if (!approximately_equal(midX - cubic[3].x, dx23 * 3 / 2)) {
+        return 0;
+    }
+    double dy10 = cubic[1].y - cubic[0].y;
+    double dy23 = cubic[2].y - cubic[3].y;
+    double midY = cubic[0].y + dy10 * 3 / 2;
+    if (!approximately_equal(midY - cubic[3].y, dy23 * 3 / 2)) {
+        return 0;
+    }
+    reduction[0] = cubic[0];
+    reduction[1].x = midX;
+    reduction[1].y = midY;
+    reduction[2] = cubic[3];
+    return 3;
+}
+
+static int check_linear(const Cubic& cubic, Cubic& reduction,
+        int minX, int maxX, int minY, int maxY) {
+    int startIndex = 0;
+    int endIndex = 3;
+    while (cubic[startIndex].approximatelyEqual(cubic[endIndex])) {
+        --endIndex;
+        if (endIndex == 0) {
+            printf("%s shouldn't get here if all four points are about equal", __FUNCTION__);
+            assert(0);
+        }
+    }
+    LineParameters lineParameters;
+    lineParameters.cubicEndPoints(cubic, startIndex, endIndex);
+    double normalSquared = lineParameters.normalSquared();
+    double distance[2]; // distance is not normalized
+    int mask = other_two(startIndex, endIndex);
+    int inner1 = startIndex ^ mask;
+    int inner2 = endIndex ^ mask;
+    lineParameters.controlPtDistance(cubic, inner1, inner2, distance);
+    double limit = normalSquared * SquaredEpsilon;
+    int index;
+    for (index = 0; index < 2; ++index) {
+        double distSq = distance[index];
+        distSq *= distSq;
+        if (distSq > limit) {
+            return 0;
+        }
+    }
+    // four are colinear: return line formed by outside
+    reduction[0] = cubic[0];
+    reduction[1] = cubic[3];
+    int sameSide1;
+    int sameSide2;
+    bool useX = cubic[maxX].x - cubic[minX].x >= cubic[maxY].y - cubic[minY].y;
+    if (useX) {
+        sameSide1 = sign(cubic[0].x - cubic[1].x) + sign(cubic[3].x - cubic[1].x);
+        sameSide2 = sign(cubic[0].x - cubic[2].x) + sign(cubic[3].x - cubic[2].x);
+    } else {
+        sameSide1 = sign(cubic[0].y - cubic[1].y) + sign(cubic[3].y - cubic[1].y);
+        sameSide2 = sign(cubic[0].y - cubic[2].y) + sign(cubic[3].y - cubic[2].y);
+    }
+    if (sameSide1 == sameSide2 && (sameSide1 & 3) != 2) {
+        return 2;
+    }
+    double tValues[2];
+    int roots;
+    if (useX) {
+        roots = SkFindCubicExtrema(cubic[0].x, cubic[1].x, cubic[2].x, cubic[3].x, tValues);
+    } else {
+        roots = SkFindCubicExtrema(cubic[0].y, cubic[1].y, cubic[2].y, cubic[3].y, tValues);
+    }
+    for (index = 0; index < roots; ++index) {
+        _Point extrema;
+        extrema.x = interp_cubic_coords(&cubic[0].x, tValues[index]);
+        extrema.y = interp_cubic_coords(&cubic[0].y, tValues[index]);
+        // sameSide > 0 means mid is smaller than either [0] or [3], so replace smaller
+        int replace;
+        if (useX) {
+            if (extrema.x < cubic[0].x ^ extrema.x < cubic[3].x) {
+                continue;
+            }
+            replace = (extrema.x < cubic[0].x | extrema.x < cubic[3].x)
+                    ^ cubic[0].x < cubic[3].x;
+        } else {
+            if (extrema.y < cubic[0].y ^ extrema.y < cubic[3].y) {
+                continue;
+            }
+            replace = (extrema.y < cubic[0].y | extrema.y < cubic[3].y)
+                    ^ cubic[0].y < cubic[3].y;
+        }
+        reduction[replace] = extrema;
+    }
+    return 2;
+}
+
+/* food for thought:
+http://objectmix.com/graphics/132906-fast-precision-driven-cubic-quadratic-piecewise-degree-reduction-algos-2-a.html
+
+Given points c1, c2, c3 and c4 of a cubic Bezier, the points of the
+corresponding quadratic Bezier are (given in convex combinations of
+points):
+
+q1 = (11/13)c1 + (3/13)c2 -(3/13)c3 + (2/13)c4
+q2 = -c1 + (3/2)c2 + (3/2)c3 - c4
+q3 = (2/13)c1 - (3/13)c2 + (3/13)c3 + (11/13)c4
+
+Of course, this curve does not interpolate the end-points, but it would
+be interesting to see the behaviour of such a curve in an applet.
+
+--
+Kalle Rutanen
+http://kaba.hilvi.org
+
+*/
+
+// reduce to a quadratic or smaller
+// look for identical points
+// look for all four points in a line 
+    // note that three points in a line doesn't simplify a cubic
+// look for approximation with single quadratic
+    // save approximation with multiple quadratics for later
+int reduceOrder(const Cubic& cubic, Cubic& reduction, ReduceOrder_Flags allowQuadratics) {
+    int index, minX, maxX, minY, maxY;
+    int minXSet, minYSet;
+    minX = maxX = minY = maxY = 0;
+    minXSet = minYSet = 0;
+    for (index = 1; index < 4; ++index) {
+        if (cubic[minX].x > cubic[index].x) {
+            minX = index;
+        }
+        if (cubic[minY].y > cubic[index].y) {
+            minY = index;
+        }
+        if (cubic[maxX].x < cubic[index].x) {
+            maxX = index;
+        }
+        if (cubic[maxY].y < cubic[index].y) {
+            maxY = index;
+        }
+    }
+    for (index = 0; index < 4; ++index) {
+        if (approximately_equal(cubic[index].x, cubic[minX].x)) {
+            minXSet |= 1 << index;
+        }
+        if (approximately_equal(cubic[index].y, cubic[minY].y)) {
+            minYSet |= 1 << index;
+        }
+    }
+    if (minXSet == 0xF) { // test for vertical line
+        if (minYSet == 0xF) { // return 1 if all four are coincident
+            return coincident_line(cubic, reduction);
+        }
+        return vertical_line(cubic, reduction);
+    }
+    if (minYSet == 0xF) { // test for horizontal line
+        return horizontal_line(cubic, reduction);
+    }
+    int result = check_linear(cubic, reduction, minX, maxX, minY, maxY);
+    if (result) {
+        return result;
+    }
+    if (allowQuadratics && (result = check_quadratic(cubic, reduction, minX, maxX, minY, maxY))) {
+        return result;
+    }
+    memcpy(reduction, cubic, sizeof(Cubic));
+    return 4;
+}
diff --git a/experimental/Intersection/CubicRoots.cpp b/experimental/Intersection/CubicRoots.cpp
new file mode 100644
index 0000000..a3da700
--- /dev/null
+++ b/experimental/Intersection/CubicRoots.cpp
@@ -0,0 +1,72 @@
+#include "CubicIntersection.h"
+
+//http://planetmath.org/encyclopedia/CubicEquation.html
+/* the roots of x^3 + ax^2 + bx + c are
+j = -2a^3 + 9ab - 27c
+k = sqrt((2a^3 - 9ab + 27c)^2 + 4(-a^2 + 3b)^3)
+t1 = -a/3 + cuberoot((j + k) / 54) + cuberoot((j - k) / 54)
+t2 = -a/3 - ( 1 + i*cuberoot(3))/2 * cuberoot((j + k) / 54)
+          + (-1 + i*cuberoot(3))/2 * cuberoot((j - k) / 54)
+t3 = -a/3 + (-1 + i*cuberoot(3))/2 * cuberoot((j + k) / 54)
+          - ( 1 + i*cuberoot(3))/2 * cuberoot((j - k) / 54)
+*/
+
+
+static bool is_unit_interval(double x) {
+    return x > 0 && x < 1;
+}
+
+const double PI = 4 * atan(1);
+
+// from SkGeometry.cpp
+int cubic_roots(const double coeff[4], double tValues[3]) {
+    if (approximately_zero(coeff[0]))   // we're just a quadratic
+    {
+        return quadratic_roots(&coeff[1], tValues);
+    }
+    double inva = 1 / coeff[0];
+    double a = coeff[1] * inva;
+    double b = coeff[2] * inva;
+    double c = coeff[3] * 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* roots = tValues;
+    double r;
+
+    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;
+        if (is_unit_interval(r))
+            *roots++ = r;
+
+        r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3;
+        if (is_unit_interval(r))
+            *roots++ = r;
+
+        r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3;
+        if (is_unit_interval(r))
+            *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;
+        if (is_unit_interval(r))
+            *roots++ = r;
+    }
+    return (int)(roots - tValues);
+}
diff --git a/experimental/Intersection/CubicSubDivide.cpp b/experimental/Intersection/CubicSubDivide.cpp
new file mode 100644
index 0000000..b8202fe
--- /dev/null
+++ b/experimental/Intersection/CubicSubDivide.cpp
@@ -0,0 +1,103 @@
+#include "CubicIntersection.h"
+#include "IntersectionUtilities.h"
+
+/*
+ Given a cubic c, t1, and t2, find a small cubic segment.
+ 
+ The new cubic is defined as points A, B, C, and D, where
+ s1 = 1 - t1
+ s2 = 1 - t2
+ A = c[0]*s1*s1*s1 + 3*c[1]*s1*s1*t1 + 3*c[2]*s1*t1*t1 + c[3]*t1*t1*t1
+ D = c[0]*s2*s2*s2 + 3*c[1]*s2*s2*t2 + 3*c[2]*s2*t2*t2 + c[3]*t2*t2*t2
+ 
+ We don't have B or C. So We define two equations to isolate them.
+ First, compute two reference T values 1/3 and 2/3 from t1 to t2:
+ 
+ c(at (2*t1 + t2)/3) == E
+ c(at (t1 + 2*t2)/3) == F
+ 
+ Next, compute where those values must be if we know the values of B and C:
+ 
+ _12   =  A*2/3 + B*1/3
+ 12_   =  A*1/3 + B*2/3
+ _23   =  B*2/3 + C*1/3
+ 23_   =  B*1/3 + C*2/3
+ _34   =  C*2/3 + D*1/3
+ 34_   =  C*1/3 + D*2/3
+ _123  = (A*2/3 + B*1/3)*2/3 + (B*2/3 + C*1/3)*1/3 = A*4/9 + B*4/9 + C*1/9
+ 123_  = (A*1/3 + B*2/3)*1/3 + (B*1/3 + C*2/3)*2/3 = A*1/9 + B*4/9 + C*4/9
+ _234  = (B*2/3 + C*1/3)*2/3 + (C*2/3 + D*1/3)*1/3 = B*4/9 + C*4/9 + D*1/9
+ 234_  = (B*1/3 + C*2/3)*1/3 + (C*1/3 + D*2/3)*2/3 = B*1/9 + C*4/9 + D*4/9
+ _1234 = (A*4/9 + B*4/9 + C*1/9)*2/3 + (B*4/9 + C*4/9 + D*1/9)*1/3
+       =  A*8/27 + B*12/27 + C*6/27 + D*1/27
+       =  E
+ 1234_ = (A*1/9 + B*4/9 + C*4/9)*1/3 + (B*1/9 + C*4/9 + D*4/9)*2/3
+       =  A*1/27 + B*6/27 + C*12/27 + D*8/27
+       =  F
+ E*27  =  A*8    + B*12   + C*6     + D
+ F*27  =  A      + B*6    + C*12    + D*8
+ 
+Group the known values on one side:
+       
+ M       = E*27 - A*8 - D     = B*12 + C* 6
+ N       = F*27 - A   - D*8   = B* 6 + C*12
+ M*2 - N = B*18
+ N*2 - M = C*18
+ B       = (M*2 - N)/18
+ C       = (N*2 - M)/18
+ */
+ 
+static double interp_cubic_coords(const double* src, double t)
+{
+    double ab = interp(src[0], src[2], t);
+    double bc = interp(src[2], src[4], t);
+    double cd = interp(src[4], src[6], t);
+    double abc = interp(ab, bc, t);
+    double bcd = interp(bc, cd, t);
+    double abcd = interp(abc, bcd, t);
+    return abcd;
+}
+ 
+void sub_divide(const Cubic& src, double t1, double t2, Cubic& dst) {
+    double ax = dst[0].x = interp_cubic_coords(&src[0].x, t1);
+    double ay = dst[0].y = interp_cubic_coords(&src[0].y, t1);
+    double ex = interp_cubic_coords(&src[0].x, (t1*2+t2)/3);
+    double ey = interp_cubic_coords(&src[0].y, (t1*2+t2)/3);
+    double fx = interp_cubic_coords(&src[0].x, (t1+t2*2)/3);
+    double fy = interp_cubic_coords(&src[0].y, (t1+t2*2)/3);
+    double dx = dst[3].x = interp_cubic_coords(&src[0].x, t2);
+    double dy = dst[3].y = interp_cubic_coords(&src[0].y, t2);
+    double mx = ex * 27 - ax * 8 - dx;
+    double my = ey * 27 - ay * 8 - dy;
+    double nx = fx * 27 - ax - dx * 8;
+    double ny = fy * 27 - ay - dy * 8;
+    /* bx = */ dst[1].x = (mx * 2 - nx) / 18;
+    /* by = */ dst[1].y = (my * 2 - ny) / 18;
+    /* cx = */ dst[2].x = (nx * 2 - mx) / 18;
+    /* cy = */ dst[2].y = (ny * 2 - my) / 18;
+}
+
+/* classic one t subdivision */
+static void interp_cubic_coords(const double* src, double* dst, double t)
+{
+    double ab = interp(src[0], src[2], t);
+    double bc = interp(src[2], src[4], t);
+    double cd = interp(src[4], src[6], t);
+    double abc = interp(ab, bc, t);
+    double bcd = interp(bc, cd, t);
+    double abcd = interp(abc, bcd, t);
+
+    dst[0] = src[0];
+    dst[2] = ab;
+    dst[4] = abc;
+    dst[6] = abcd;
+    dst[8] = bcd;
+    dst[10] = cd;
+    dst[12] = src[6];
+}
+
+void chop_at(const Cubic& src, CubicPair& dst, double t)
+{
+    interp_cubic_coords(&src[0].x, &dst.pts[0].x, t);
+    interp_cubic_coords(&src[0].y, &dst.pts[0].y, t);
+}
diff --git a/experimental/Intersection/DataTypes.cpp b/experimental/Intersection/DataTypes.cpp
new file mode 100644
index 0000000..4cb2275
--- /dev/null
+++ b/experimental/Intersection/DataTypes.cpp
@@ -0,0 +1,74 @@
+#include "DataTypes.h"
+
+const double PointEpsilon = 0.000001;
+const double SquaredEpsilon = PointEpsilon * PointEpsilon;
+
+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)) {
+            return false;
+        }
+        memcpy(rotPath, cubic, sizeof(Cubic));
+        return true;
+    }
+    for (int index = 0; index < 4; ++index) {
+        rotPath[index].x = cubic[index].x * dx + cubic[index].y * dy;
+        rotPath[index].y = cubic[index].y * dx - cubic[index].x * dy;
+    }
+    return true;
+}
+
+double t_at(const _Line& line, const _Point& pt) {
+    double dx = line[1].x - line[0].x;
+    double dy = line[1].y - line[0].y;
+    if (fabs(dx) > fabs(dy)) {
+        if (approximately_zero(dx)) {
+            return 0;
+        }
+        return (pt.x - line[0].x) / dx;
+    }
+    if (approximately_zero(dy)) {
+        return 0;
+    }
+    return (pt.y - line[0].y) / dy;
+}
+
+static void setMinMax(double x, double y1, double y2, int flags,
+        double& minX, double& maxX) {
+    if (minX > x && (flags & (kFindTopMin | kFindBottomMin))) {
+        minX = x;
+    }
+    if (maxX < x && (flags & (kFindTopMax | kFindBottomMax))) {
+        maxX = x;
+    }
+}
+
+void x_at(const _Point& p1, const _Point& p2, double top, double bottom,
+        int flags, double& minX, double& maxX) {
+    if (approximately_equal(p1.y, p2.y)) {
+        // It should be OK to bail early in this case. There's another edge
+        // which shares this end point which can intersect without failing to 
+        // have a slope ... maybe
+        return;
+    }
+    
+    // p2.x is always greater than p1.x -- the part of points (p1, p2) are
+    // moving from the start of the cubic towards its end.
+    // if p1.y < p2.y, minX can be affected
+    // if p1.y > p2.y, maxX can be affected 
+    double slope = (p2.x - p1.x) / (p2.y - p1.y);
+    int topFlags = flags & (kFindTopMin | kFindTopMax);
+    if (topFlags && (top <= p1.y && top >= p2.y
+            || top >= p1.y && top <= p2.y)) {
+        double x = p1.x + (top - p1.y) * slope;
+        setMinMax(x, p1.y, p2.y, topFlags, minX, maxX);
+    }
+    int bottomFlags = flags & (kFindBottomMin | kFindBottomMax);
+    if (bottomFlags && (bottom <= p1.y && bottom >= p2.y
+            || bottom >= p1.y && bottom <= p2.y)) {
+        double x = p1.x + (bottom - p1.y) * slope;
+        setMinMax(x, p1.y, p2.y, bottomFlags, minX, maxX);
+    }
+}
diff --git a/experimental/Intersection/DataTypes.h b/experimental/Intersection/DataTypes.h
new file mode 100644
index 0000000..b1cbdf7
--- /dev/null
+++ b/experimental/Intersection/DataTypes.h
@@ -0,0 +1,114 @@
+#ifndef __DataTypes_h__
+#define __DataTypes_h__
+
+extern const double PointEpsilon;
+extern const double SquaredEpsilon;
+
+inline bool approximately_equal(double x, double y) {
+    return fabs(x - y) < PointEpsilon;
+}
+
+inline bool approximately_equal_squared(double x, double y) {
+    return fabs(x - y) < SquaredEpsilon;
+}
+
+inline bool approximately_greater(double x, double y) {
+    return x > y - PointEpsilon;
+}
+
+inline bool approximately_lesser(double x, double y) {
+    return x < y + PointEpsilon;
+}
+
+inline bool approximately_zero(double x) {
+    return fabs(x) < PointEpsilon;
+}
+
+inline bool approximately_zero_squared(double x) {
+    return fabs(x) < SquaredEpsilon;
+}
+
+inline bool approximately_negative(double x) {
+    return x < PointEpsilon;
+}
+
+struct _Point {
+    double x;
+    double y;
+
+    void operator-=(const _Point& v) {
+        x -= v.x;
+        y -= v.y;
+    }
+
+    friend bool operator==(const _Point& a, const _Point& b) {
+        return a.x == b.x && a.y == b.y;
+    }
+
+    friend bool operator!=(const _Point& a, const _Point& b) {
+        return a.x!= b.x || a.y != b.y;
+    }
+    
+    bool approximatelyEqual(const _Point& a) const {
+        return approximately_equal(a.y, y) && approximately_equal(a.x, x);
+    }
+
+};
+
+typedef _Point _Line[2];
+typedef _Point Quadratic[3];
+typedef _Point Cubic[4];
+
+struct _Rect {
+    double left;
+    double top;
+    double right;
+    double bottom;
+    
+    void add(const _Point& pt) {
+        if (left > pt.x) {
+            left = pt.x;
+        }
+        if (top > pt.y) {
+            top = pt.y;
+        }
+        if (right < pt.x) {
+            right = pt.x;
+        }
+        if (bottom < pt.y) {
+            bottom = pt.y;
+        }
+    }
+    
+    void setBounds(const _Line& line) {
+        left = right = line[0].x;
+        top = bottom = line[0].y;
+        add(line[1]);
+    }
+};
+
+struct CubicPair {
+    const Cubic& first() const { return (const Cubic&) pts[0]; }
+    const Cubic& second() const { return (const Cubic&) pts[3]; }
+    _Point pts[7];
+};
+
+struct QuadraticPair {
+    const Quadratic& first() const { return (const Quadratic&) pts[0]; }
+    const Quadratic& second() const { return (const Quadratic&) pts[2]; }
+    _Point pts[5];
+};
+
+enum x_at_flags {
+    kFindTopMin = 1,
+    kFindTopMax = 2,
+    kFindBottomMin = 4,
+    kFindBottomMax = 8
+};
+
+bool rotate(const Cubic& cubic, int zero, int index, Cubic& rotPath);
+double t_at(const _Line&, const _Point& );
+void x_at(const _Point& p1, const _Point& p2, double minY, double maxY,
+        int flags, double& tMin, double& tMax);
+
+#endif // __DataTypes_h__
diff --git a/experimental/Intersection/EdgeApp.cpp b/experimental/Intersection/EdgeApp.cpp
new file mode 100644
index 0000000..58ff153
--- /dev/null
+++ b/experimental/Intersection/EdgeApp.cpp
@@ -0,0 +1,182 @@
+/*
+ *  EdgeApp.cpp
+ *  edge
+ *
+ *  Created by Cary Clark on 7/6/11.
+ *  Copyright 2011 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "SkCanvas.h"
+#include "SkDevice.h"
+#include "SkGraphics.h"
+#include "SkImageEncoder.h"
+#include "SkPaint.h"
+#include "SkPicture.h"
+#include "SkStream.h"
+#include "SkTime.h"
+#include "SkWindow.h"
+
+#include "SkTouchGesture.h"
+#include "SkTypeface.h"
+
+#include "CubicIntersection_Tests.h"
+#include "CubicIntersection_TestData.h"
+
+extern void CreateSweep(SkBitmap* , float width);
+extern void CreateHorz(SkBitmap* );
+extern void CreateVert(SkBitmap* );
+extern void CreateAngle(SkBitmap* sweep, float angle);
+extern void SkAntiEdge_Test();
+
+
+static const char gCharEvtName[] = "Char_Event";
+static const char gKeyEvtName[] = "Key_Event";
+
+class SkEdgeView : public SkView {
+public:
+    SkEdgeView() {
+        CreateSweep(&fSweep_1_0, 1);
+        CreateSweep(&fSweep_1_2, 1.2f);
+        CreateSweep(&fSweep_1_4, 1.4f);
+        CreateSweep(&fSweep_1_6, 1.6f);
+        CreateHorz(&fBitmapH);
+        CreateVert(&fBitmapV);
+        CreateAngle(&fAngle_12, 12);
+        CreateAngle(&fAngle_45, 45);
+    }
+    virtual ~SkEdgeView() {}
+
+protected:
+    virtual void onDraw(SkCanvas* canvas) {
+        canvas->drawColor(SK_ColorWHITE);
+        canvas->drawBitmap(fSweep_1_0, 0, 10);
+        canvas->drawBitmap(fBitmapH, 110, 10);
+        canvas->drawBitmap(fBitmapV, 220, 10);
+        canvas->drawBitmap(fSweep_1_2, 0, 110);
+        canvas->drawBitmap(fSweep_1_4, 100, 110);
+        canvas->drawBitmap(fSweep_1_6, 200, 110);
+        canvas->drawBitmap(fAngle_12, 0, 200);
+        canvas->drawBitmap(fAngle_45, 124, 220);
+    }
+
+private:
+    SkBitmap fSweep_1_0;
+    SkBitmap fSweep_1_2;
+    SkBitmap fSweep_1_4;
+    SkBitmap fSweep_1_6;
+    SkBitmap fBitmapH;
+    SkBitmap fBitmapV;
+    SkBitmap fAngle_12;
+    SkBitmap fAngle_45;
+    typedef SkView INHERITED;
+};
+
+class EdgeWindow : public SkOSWindow {
+public:
+    EdgeWindow(void* hwnd) : INHERITED(hwnd) {
+        this->setConfig(SkBitmap::kARGB_8888_Config);
+        this->setVisibleP(true);
+        fView.setVisibleP(true);
+        fView.setClipToBounds(false);
+        this->attachChildToFront(&fView)->unref();
+    }
+    virtual ~EdgeWindow() {}
+
+    virtual void draw(SkCanvas* canvas){
+        this->INHERITED::draw(canvas);
+    }
+
+
+protected:
+    SkEdgeView fView;
+    
+    virtual void onDraw(SkCanvas* canvas) {
+    }
+
+    virtual bool onHandleKey(SkKey key) {
+        SkEvent evt(gKeyEvtName);
+        evt.setFast32(key);
+        if (fView.doQuery(&evt)) {
+            return true;
+        }
+        return this->INHERITED::onHandleKey(key);
+    }
+
+    virtual bool onHandleChar(SkUnichar uni) {
+        SkEvent evt(gCharEvtName);
+        evt.setFast32(uni);
+        if (fView.doQuery(&evt)) {
+                return true;
+        }
+        return this->INHERITED::onHandleChar(uni);
+    }
+
+    virtual void onSizeChange() {
+        fView.setSize(this->width(), this->height());
+        this->INHERITED::onSizeChange();
+    }
+
+    virtual SkCanvas* beforeChildren(SkCanvas* canvas) {
+        return this->INHERITED::beforeChildren(canvas);
+    }
+    
+    virtual void afterChildren(SkCanvas*) {}
+    virtual void beforeChild(SkView* child, SkCanvas* canvas) {}
+    virtual void afterChild(SkView* child, SkCanvas* canvas) {}
+
+    virtual bool onEvent(const SkEvent& evt) {
+        return this->INHERITED::onEvent(evt);
+    }
+    
+    virtual bool onQuery(SkEvent* evt) {
+        return this->INHERITED::onQuery(evt);
+    }
+
+    virtual bool onDispatchClick(int x, int y, Click::State state, void* owner) {
+        int w = SkScalarRound(this->width());
+        int h = SkScalarRound(this->height());
+        // check for the resize-box
+        if (w - x < 16 && h - y < 16) {
+            return false;   // let the OS handle the click
+        } else {
+            return this->INHERITED::onDispatchClick(x, y, state, owner);
+        }
+    }
+
+    virtual bool onClick(Click* click) {
+        return false;
+    }
+
+    virtual Click* onFindClickHandler(SkScalar x, SkScalar y) {
+        return 0;
+    }
+
+    typedef SkOSWindow INHERITED;
+};
+
+
+#include "SkApplication.h"
+
+SkOSWindow* create_sk_window(void* hwnd, int argc, char** argv) {
+    return new EdgeWindow(hwnd);
+}
+
+void cubecode_test();
+
+void application_init() {
+    unsigned foo = 4;
+    SkGraphics::Init();
+    SkEvent::Init();
+    cubecode_test();
+
+    convert_testx();
+    CubicIntersection_Tests();
+    SkAntiEdge_Test();
+}
+
+void application_term() {
+    SkEvent::Term();
+    SkGraphics::Term();
+}
+
diff --git a/experimental/Intersection/Extrema.cpp b/experimental/Intersection/Extrema.cpp
new file mode 100644
index 0000000..0b85060
--- /dev/null
+++ b/experimental/Intersection/Extrema.cpp
@@ -0,0 +1,77 @@
+#include "DataTypes.h"
+#include "Extrema.h"
+
+static int valid_unit_divide(double numer, double denom, double* ratio)
+{
+    if (numer < 0)
+    {
+        numer = -numer;
+        denom = -denom;
+    }
+
+    if (denom == 0 || numer == 0 || numer >= denom)
+        return 0;
+
+    double r = numer / denom;
+    if (r == 0) // catch underflow if numer <<<< denom
+        return 0;
+    *ratio = r;
+    return 1;
+}
+
+/** From Numerical Recipes in C.
+
+    Q = -1/2 (B + sign(B) sqrt[B*B - 4*A*C])
+    x1 = Q / A
+    x2 = C / Q
+*/
+static int SkFindUnitQuadRoots(double A, double B, double C, double roots[2])
+{
+    if (A == 0)
+        return valid_unit_divide(-C, B, roots);
+
+    double* r = roots;
+
+    double R = B*B - 4*A*C;
+    if (R < 0) {  // complex roots
+        return 0;
+    }
+    R = sqrt(R);
+
+    double Q = (B < 0) ? -(B-R)/2 : -(B+R)/2;
+    r += valid_unit_divide(Q, A, r);
+    r += valid_unit_divide(C, Q, r);
+    if (r - roots == 2 && approximately_equal(roots[0], roots[1])) { // nearly-equal?
+        r -= 1; // skip the double root
+    }
+    return (int)(r - roots);
+}
+
+/** Cubic'(t) = At^2 + Bt + C, where
+    A = 3(-a + 3(b - c) + d)
+    B = 6(a - 2b + c)
+    C = 3(b - a)
+    Solve for t, keeping only those that fit betwee 0 < t < 1
+*/
+int SkFindCubicExtrema(double a, double b, double c, double d, double tValues[2])
+{
+    // we divide A,B,C by 3 to simplify
+    double A = d - a + 3*(b - c);
+    double B = 2*(a - b - b + c);
+    double C = b - a;
+
+    return SkFindUnitQuadRoots(A, B, C, tValues);
+}
+
+/** Quad'(t) = At + B, where
+    A = 2(a - 2b + c)
+    B = 2(b - a)
+    Solve for t, only if it fits between 0 < t < 1
+*/
+int SkFindQuadExtrema(double a, double b, double c, double tValue[1])
+{
+    /*  At + B == 0
+        t = -B / A
+    */
+    return valid_unit_divide(a - b, a - b - b + c, tValue);
+}
diff --git a/experimental/Intersection/Extrema.h b/experimental/Intersection/Extrema.h
new file mode 100644
index 0000000..aae422e
--- /dev/null
+++ b/experimental/Intersection/Extrema.h
@@ -0,0 +1,2 @@
+int SkFindCubicExtrema(double a, double b, double c, double d, double tValues[2]);
+int SkFindQuadExtrema(double a, double b, double c, double tValue[1]);
diff --git a/experimental/Intersection/Inline_Tests.cpp b/experimental/Intersection/Inline_Tests.cpp
new file mode 100644
index 0000000..59caf30
--- /dev/null
+++ b/experimental/Intersection/Inline_Tests.cpp
@@ -0,0 +1,48 @@
+#include "CubicIntersection.h"
+#include "CubicIntersection_Tests.h"
+#include "IntersectionUtilities.h"
+
+static void assert_that(int x, int y, const char* s) {
+    if (x == y) {
+        return;
+    }
+    printf("result=%d expected=%d %s\n", x, y, s);
+}
+
+static void side_test() {
+    assert_that(side(-1), 0, "side(-1) != 0");
+    assert_that(side(0), 1, "side(0) != 1");
+    assert_that(side(1), 2, "side(1) != 2");
+}
+
+static void sideBit_test() {
+    assert_that(sideBit(-1), 1, "sideBit(-1) != 1");
+    assert_that(sideBit(0), 2, "sideBit(0) != 2");
+    assert_that(sideBit(1), 4, "sideBit(1) != 4");
+}
+
+static void other_two_test() {
+    for (int x = 0; x < 4; ++x) {
+        for (int y = 0; y < 4; ++y) {
+            if (x == y) {
+                continue;
+            }
+            int mask = other_two(x, y);
+            int all = 1 << x;
+            all |= 1 << y;
+            all |= 1 << (x ^ mask);
+            all |= 1 << (y ^ mask);
+            if (all == 0x0F) {
+                continue;
+            }
+            printf("[%d,%d] other_two failed mask=%d [%d,%d]\n",
+                x, y, mask, x ^ mask, y ^ mask);
+        }
+    }
+}
+
+void Inline_Tests() {
+    side_test();
+    sideBit_test();
+    other_two_test();
+}
diff --git a/experimental/Intersection/IntersectionUtilities.cpp b/experimental/Intersection/IntersectionUtilities.cpp
new file mode 100644
index 0000000..d05afdf
--- /dev/null
+++ b/experimental/Intersection/IntersectionUtilities.cpp
@@ -0,0 +1,40 @@
+
+// snippets that one day may be useful, unused for now...
+
+// get sign, exponent, mantissa from double
+// Translate the double into sign, exponent and mantissa.
+    long bits = BitConverter.DoubleToInt64Bits(d);
+    // Note that the shift is sign-extended, hence the test against -1 not 1
+    bool negative = (bits < 0);
+    int exponent = (int) ((bits >> 52) & 0x7ffL);
+    long mantissa = bits & 0xfffffffffffffL;
+
+    // Subnormal numbers; exponent is effectively one higher,
+    // but there's no extra normalisation bit in the mantissa
+    if (exponent==0)
+    {
+        exponent++;
+    }
+    // Normal numbers; leave exponent as it is but add extra
+    // bit to the front of the mantissa
+    else
+    {
+        mantissa = mantissa | (1L<<52);
+    }
+
+    // Bias the exponent. It's actually biased by 1023, but we're
+    // treating the mantissa as m.0 rather than 0.m, so we need
+    // to subtract another 52 from it.
+    exponent -= 1075;
+
+    if (mantissa == 0) 
+    {
+        return "0";
+    }
+
+    /* Normalize */
+    while((mantissa & 1) == 0) 
+    {    /*  i.e., Mantissa is even */
+        mantissa >>= 1;
+        exponent++;
+    }
\ No newline at end of file
diff --git a/experimental/Intersection/IntersectionUtilities.h b/experimental/Intersection/IntersectionUtilities.h
new file mode 100644
index 0000000..3dfd311
--- /dev/null
+++ b/experimental/Intersection/IntersectionUtilities.h
@@ -0,0 +1,41 @@
+
+// inline utilities
+/* Returns 0 if negative, 1 if zero, 2 if positive
+*/
+inline int side(double x) {
+    return (x > 0) + (x >= 0);
+}
+
+/* Returns 1 if negative, 2 if zero, 4 if positive
+*/
+inline int sideBit(double x) {
+    return 1 << side(x);
+}
+
+/* Given the set [0, 1, 2, 3], and two of the four members, compute an XOR mask
+   that computes the other two. Note that:
+   
+   one ^ two == 3 for (0, 3), (1, 2)
+   one ^ two <  3 for (0, 1), (0, 2), (1, 3), (2, 3)
+   3 - (one ^ two) is either 0, 1, or 2
+   1 >> 3 - (one ^ two) is either 0 or 1
+thus:
+   returned == 2 for (0, 3), (1, 2)
+   returned == 3 for (0, 1), (0, 2), (1, 3), (2, 3)
+given that:
+   (0, 3) ^ 2 -> (2, 1)  (1, 2) ^ 2 -> (3, 0)
+   (0, 1) ^ 3 -> (3, 2)  (0, 2) ^ 3 -> (3, 1)  (1, 3) ^ 3 -> (2, 0)  (2, 3) ^ 3 -> (1, 0)
+*/
+inline int other_two(int one, int two) {
+    return 1 >> 3 - (one ^ two) ^ 3;
+}
+
+/* Returns -1 if negative, 0 if zero, 1 if positive
+*/
+inline int sign(double x) {
+    return (x > 0) - (x < 0);
+}
+
+inline double interp(double A, double B, double t) {
+    return A + (B - A) * t;
+}
diff --git a/experimental/Intersection/Intersections.h b/experimental/Intersection/Intersections.h
new file mode 100644
index 0000000..928926e
--- /dev/null
+++ b/experimental/Intersection/Intersections.h
@@ -0,0 +1,49 @@
+class Intersections {
+public:
+    Intersections()
+        : fUsed(0)
+        , fSwap(0)
+    {
+        bzero(fT, sizeof(fT));
+    }
+
+    void add(double one, double two) {
+        if (fUsed > 0 && approximately_equal(fT[fSwap][fUsed - 1], one)
+                && approximately_equal(fT[fSwap ^ 1][fUsed - 1], two)) {
+            return;
+        }
+        fT[fSwap][fUsed] = one;
+        fT[fSwap ^ 1][fUsed] = two;
+        ++fUsed;
+    }
+
+    void offset(int base, double start, double end) {
+        for (int index = base; index < fUsed; ++index) {
+            double val = fT[fSwap][index];
+            val *= end - start;
+            val += start;
+            fT[fSwap][index] = val;
+        }
+    }
+
+    bool intersected() {
+        return fUsed > 0;
+    }
+
+    void swap() {
+        fSwap ^= 1;
+    }
+    
+    bool swapped() {
+        return fSwap;
+    }
+
+    int used() {
+        return fUsed;
+    }
+
+    double fT[2][9];
+private:
+    int fUsed;
+    int fSwap;
+};
diff --git a/experimental/Intersection/LineIntersection.cpp b/experimental/Intersection/LineIntersection.cpp
new file mode 100644
index 0000000..d23401e
--- /dev/null
+++ b/experimental/Intersection/LineIntersection.cpp
@@ -0,0 +1,120 @@
+#include "DataTypes.h"
+#include "LineIntersection.h"
+#include "LineParameters.h"
+
+
+static bool no_intersection(_Point* result) {
+    result->x = 0;
+    result->y = 0;
+    return false;
+}
+
+/*
+   Determine the intersection point of two line segments
+   Return FALSE if the lines don't intersect
+   from: http://paulbourke.net/geometry/lineline2d/
+   min: 8 subs, 4 muls, 4 cmps
+*/
+
+bool lineIntersect(const _Line& a, const _Line& b, _Point* result) {
+    LineParameters paramsA, paramsB;
+    double axLen = a[1].x - a[0].x;
+    double ayLen = a[1].y - a[0].y;
+    double bxLen = b[1].x - b[0].x;
+    double byLen = b[1].y - b[0].y;
+    double denom  = byLen * axLen - ayLen * bxLen;
+    if (approximately_zero_squared(denom)) { // is either 'a' or 'b' a point?
+        bool aIsPoint = approximately_zero(axLen) && approximately_zero(ayLen);
+        bool bIsPoint = approximately_zero(bxLen) && approximately_zero(byLen);
+        if (aIsPoint & bIsPoint) {
+            if (!approximately_equal(a[0].x, b[0].x)
+                    || !approximately_equal(a[0].y, b[0].y)) {
+                return no_intersection(result);
+            }
+        } else {
+            double ptToLineDistance;
+            if (aIsPoint) {
+                paramsB.lineEndPoints(b);
+                ptToLineDistance = paramsB.pointDistance(a[0]);
+            } else {
+                paramsA.lineEndPoints(a);
+                ptToLineDistance = paramsA.pointDistance(b[0]);
+            }
+            if (!approximately_zero(ptToLineDistance)) {
+                return no_intersection(result);
+            }
+        }
+        if (aIsPoint) {
+            result->x = a[0].x;
+            result->y = a[0].y;
+        } else {
+            result->x = b[0].x;
+            result->y = b[0].y;
+        }
+        return true;
+    }
+    double ab0y = a[0].y - b[0].y;
+    double ab0x = a[0].x - b[0].x;
+    double numerA = ab0y * bxLen - byLen * ab0x;
+    double numerB;
+    if (numerA < 0 && denom > numerA || numerA > 0 && denom < numerA) {
+        return no_intersection(result);
+    }
+    numerB = ab0y * axLen - ayLen * ab0x;
+    if (numerB < 0 && denom > numerB || numerB > 0 && denom < numerB) {
+        return no_intersection(result);
+    }
+    /* Are the line coincident? See if they overlap */
+    paramsA.lineEndPoints(a);
+    double b0dist = paramsA.pointDistance(b[0]);
+    bool b0on = approximately_zero_squared(b0dist);
+    double b1dist = paramsA.pointDistance(b[1]);
+    bool b1on = approximately_zero_squared(b1dist);
+    if (b0on | b1on) {
+        if (b0on & b1on) {
+            result->x = (b[0].x + b[1].x) / 2;
+            result->y = (b[0].y + b[1].y) / 2;
+            return true;
+        }
+        paramsB.lineEndPoints(b);
+        double a0dist = paramsB.pointDistance(a[0]);
+        bool a0on = approximately_zero_squared(a0dist);
+        double a1dist = paramsB.pointDistance(a[1]);
+        bool a1on = approximately_zero_squared(a1dist);
+        if (a0on | a1on) {
+            if (a0on & a1on) {
+                result->x = (a[0].x + a[1].x) / 2;
+                result->y = (a[0].y + a[1].y) / 2;
+                return true;
+            }
+            const _Point& aPt = a0on ? a[0] : a[1];
+            const _Point& bPt = b0on ? b[0] : b[1];
+            if (aPt.approximatelyEqual(bPt)) {
+                *result = aPt;
+                return true;
+            }
+            result->x = (aPt.x + bPt.x) / 2;
+            result->y = (aPt.y + bPt.y) / 2;
+            return true;
+        }
+    }
+
+    /* Is the intersection along the the segments */
+    double mua = numerA / denom;
+    result->x = a[0].x + mua * (a[1].x - a[0].x);
+    result->y = a[0].y + mua * (a[1].y - a[0].y);
+    return true;
+}
+
+
+// from http://www.bryceboe.com/wordpress/wp-content/uploads/2006/10/intersect.py
+// 4 subs, 2 muls, 1 cmp
+static bool ccw(const _Point& A, const _Point& B, const _Point& C) {
+	return (C.y - A.y) * (B.x - A.x) > (B.y - A.y) * (C.x - A.x);
+}
+
+// 16 subs, 8 muls, 6 cmps
+bool testIntersect(const _Line& a, const _Line& b) {
+	return ccw(a[0], b[0], b[1]) != ccw(a[1], b[0], b[1])
+            && ccw(a[0], a[1], b[0]) != ccw(a[0], a[1], b[1]);
+}
diff --git a/experimental/Intersection/LineIntersection.h b/experimental/Intersection/LineIntersection.h
new file mode 100644
index 0000000..c1246a5
--- /dev/null
+++ b/experimental/Intersection/LineIntersection.h
@@ -0,0 +1,5 @@
+
+#include "DataTypes.h"
+
+bool lineIntersect(const _Line& a, const _Line& b, _Point* result);
+bool testIntersect(const _Line& a, const _Line& b);
diff --git a/experimental/Intersection/LineIntersection_Test.cpp b/experimental/Intersection/LineIntersection_Test.cpp
new file mode 100644
index 0000000..f7d00b5
--- /dev/null
+++ b/experimental/Intersection/LineIntersection_Test.cpp
@@ -0,0 +1,21 @@
+#include "CubicIntersection_Tests.h"
+#include "LineIntersection.h"
+
+const _Line tests[][2] = {
+    {{{166.86950047022856, 112.69654129527828}, {166.86948801592692, 112.69655741235339}}, 
+     {{166.86960700313026, 112.6965477747386},  {166.86925794355412, 112.69656471103423}}},
+};
+
+const size_t tests_count = sizeof(tests) / sizeof(tests[0]);
+
+static size_t firstLineIntersectionTest = 0;
+
+void LineIntersection_Test() {
+    for (size_t index = firstLineIntersectionTest; index < tests_count; ++index) {
+        const _Line& line1 = tests[index][0];
+        const _Line& line2 = tests[index][1];
+        _Point result;
+        lineIntersect(line1, line2, &result);
+        printf("%s (%g,%g)\n", __FUNCTION__, result.x, result.y);
+    }
+}
diff --git a/experimental/Intersection/LineParameters.h b/experimental/Intersection/LineParameters.h
new file mode 100644
index 0000000..1889622
--- /dev/null
+++ b/experimental/Intersection/LineParameters.h
@@ -0,0 +1,94 @@
+#include "DataTypes.h"
+
+// Sources
+// computer-aided design - volume 22 number 9 november 1990 pp 538 - 549
+// online at http://cagd.cs.byu.edu/~tom/papers/bezclip.pdf
+
+class LineParameters {
+public:
+    void cubicEndPoints(const Cubic& pts) {
+        a = pts[0].y - pts[3].y;
+        b = pts[3].x - pts[0].x;
+        c = pts[0].x * pts[3].y - pts[3].x * pts[0].y;
+    }
+    
+    void cubicEndPoints(const Cubic& pts, int s, int e) {
+        a = pts[s].y - pts[e].y;
+        b = pts[e].x - pts[s].x;
+        c = pts[s].x * pts[e].y - pts[e].x * pts[s].y;
+    }
+    
+    void lineEndPoints(const _Line& pts) {
+        a = pts[0].y - pts[1].y;
+        b = pts[1].x - pts[0].x;
+        c = pts[0].x * pts[1].y - pts[1].x * pts[0].y;
+    }
+    
+    void quadEndPoints(const Quadratic& pts) {
+        a = pts[0].y - pts[2].y;
+        b = pts[2].x - pts[0].x;
+        c = pts[0].x * pts[2].y - pts[2].x * pts[0].y;
+    }
+
+    void quadEndPoints(const Quadratic& pts, int s, int e) {
+        a = pts[s].y - pts[e].y;
+        b = pts[e].x - pts[s].x;
+        c = pts[s].x * pts[e].y - pts[e].x * pts[s].y;
+    }
+
+    double normalSquared() {
+        return a * a + b * b;
+    }
+
+    bool normalize() {
+        double normal = sqrt(normalSquared());
+        if (normal < SquaredEpsilon) {
+            a = b = c = 0;
+            return false;
+        }
+        double reciprocal = 1 / normal;
+        a *= reciprocal;
+        b *= reciprocal;
+        c *= reciprocal;
+        return true;
+    }
+    
+    void cubicDistanceY(const Cubic& pts, Cubic& distance) {
+        double oneThird = 1 / 3.0;
+        for (int index = 0; index < 4; ++index) {
+            distance[index].x = index * oneThird;
+            distance[index].y = a * pts[index].x + b * pts[index].y + c;
+        }
+    }
+
+    void quadDistanceY(const Quadratic& pts, Quadratic& distance) {
+        double oneHalf = 1 / 2.0;
+        for (int index = 0; index < 3; ++index) {
+            distance[index].x = index * oneHalf;
+            distance[index].y = a * pts[index].x + b * pts[index].y + c;
+        }
+    }
+
+    void controlPtDistance(const Cubic& pts, double distance[2]) {
+        for (int index = 0; index < 2; ++index) {
+            distance[index] = a * pts[index + 1].x + b * pts[index + 1].y + c;
+        }
+    }
+    
+    void controlPtDistance(const Cubic& pts, int i, int j, double distance[2]) {
+        distance[0] = a * pts[i].x + b * pts[i].y + c;
+        distance[1] = a * pts[j].x + b * pts[j].y + c;
+    }
+    
+    double controlPtDistance(const Quadratic& pts) {
+        return a * pts[1].x + b * pts[1].y + c;
+    }
+    
+    double pointDistance(const _Point& pt) {
+        return a * pt.x + b * pt.y + c;
+    }
+private:
+    double a;
+    double b;
+    double c;
+};
diff --git a/experimental/Intersection/LineParameteters_Test.cpp b/experimental/Intersection/LineParameteters_Test.cpp
new file mode 100644
index 0000000..5e0b4ac
--- /dev/null
+++ b/experimental/Intersection/LineParameteters_Test.cpp
@@ -0,0 +1,72 @@
+#include "CubicIntersection_Tests.h"
+#include "LineParameters.h"
+
+
+// tests to verify that distance calculations are coded correctly
+const Cubic tests[] = {
+    {{0, 0}, {1, 1}, {2, 2}, {0, 3}},
+    {{0, 0}, {1, 1}, {2, 2}, {3, 0}},
+    {{0, 0}, {5, 0}, {-2,4}, {3, 4}},
+    {{0, 2}, {1, 0}, {2, 0}, {3, 0}},
+    {{0, .2}, {1, 0}, {2, 0}, {3, 0}},
+    {{0, .02}, {1, 0}, {2, 0}, {3, 0}},
+    {{0, .002}, {1, 0}, {2, 0}, {3, 0}},
+    {{0, .0002}, {1, 0}, {2, 0}, {3, 0}},
+    {{0, .00002}, {1, 0}, {2, 0}, {3, 0}},
+    {{0, PointEpsilon * 2}, {1, 0}, {2, 0}, {3, 0}}, 
+};
+
+const double answers[][2] = {
+    {1, 2},
+    {1, 2},
+    {4, 4},
+    {1.1094003924, 0.5547001962},
+    {0.133038021, 0.06651901052},
+    {0.0133330370, 0.006666518523},
+    {0.001333333037, 0.0006666665185},
+    {0.000133333333, 6.666666652e-05},
+    {1.333333333e-05, 6.666666667e-06},
+    {1.333333333e-06, 6.666666667e-07},
+};
+
+const size_t tests_count = sizeof(tests) / sizeof(tests[0]);
+
+static size_t firstLineParameterTest = 0;
+
+void LineParameter_Test() {
+    for (size_t index = firstLineParameterTest; index < tests_count; ++index) {
+        LineParameters lineParameters;
+        const Cubic& cubic = tests[index];
+        lineParameters.cubicEndPoints(cubic);
+        double denormalizedDistance[2];
+        lineParameters.controlPtDistance(cubic, denormalizedDistance);
+        double normalSquared = lineParameters.normalSquared();
+        size_t inner;
+        for (inner = 0; inner < 2; ++inner) {
+            double distSq = denormalizedDistance[inner];
+            distSq *= distSq;
+            double answersSq = answers[index][inner];
+            answersSq *= answersSq;
+            if (approximately_equal(distSq, normalSquared * answersSq)) {
+                continue;
+            }
+            printf("%s [%d,%d] denormalizedDistance:%g != answer:%g"
+                    " distSq:%g answerSq:%g normalSquared:%g\n",
+                    __FUNCTION__, (int)index, (int)inner,
+                    denormalizedDistance[inner], answers[index][inner],
+                    distSq, answersSq, normalSquared);
+        }
+        lineParameters.normalize();
+        double normalizedDistance[2];
+        lineParameters.controlPtDistance(cubic, normalizedDistance);
+        for (inner = 0; inner < 2; ++inner) {
+            if (approximately_equal(fabs(normalizedDistance[inner]), 
+                    answers[index][inner])) {
+                continue;
+            }
+            printf("%s [%d,%d] normalizedDistance:%1.10g != answer:%g\n",
+                    __FUNCTION__, (int)index, (int)inner,
+                    normalizedDistance[inner], answers[index][inner]);
+        }
+    }
+}
diff --git a/experimental/Intersection/QuadraticBezierClip.cpp b/experimental/Intersection/QuadraticBezierClip.cpp
new file mode 100644
index 0000000..1114d74
--- /dev/null
+++ b/experimental/Intersection/QuadraticBezierClip.cpp
@@ -0,0 +1,57 @@
+#include "CubicIntersection.h"
+#include "LineParameters.h"
+#include <algorithm> // used for std::swap
+
+// return false if unable to clip (e.g., unable to create implicit line)
+// caller should subdivide, or create degenerate if the values are too small
+bool bezier_clip(const Quadratic& q1, const Quadratic& q2, double& minT, double& maxT) {
+    minT = 1;
+    maxT = 0;
+    // determine normalized implicit line equation for pt[0] to pt[3]
+    //   of the form ax + by + c = 0, where a*a + b*b == 1
+    
+    // find the implicit line equation parameters
+    LineParameters endLine;
+    endLine.quadEndPoints(q1);
+    if (!endLine.normalize()) {
+        printf("line cannot be normalized: need more code here\n");
+        return false;
+    }
+
+    double distance = endLine.controlPtDistance(q1);
+    
+    // find fat line
+    double top = 0;
+    double bottom = distance / 2; // http://students.cs.byu.edu/~tom/557/text/cic.pdf (7.6)
+    if (top > bottom) {
+        std::swap(top, bottom);
+    }
+    
+    // compute intersecting candidate distance
+    Quadratic distance2y; // points with X of (0, 1/2, 1)
+    endLine.quadDistanceY(q2, distance2y);
+    
+    int flags = 0;
+    if (approximately_lesser(distance2y[0].y, top)) {
+        flags |= kFindTopMin;
+    } else if (approximately_greater(distance2y[0].y, bottom)) {
+        flags |= kFindBottomMin;
+    } else {
+        minT = 0;
+    }
+
+    if (approximately_lesser(distance2y[2].y, top)) {
+        flags |= kFindTopMax;
+    } else if (approximately_greater(distance2y[2].y, bottom)) {
+        flags |= kFindBottomMax;
+    } else {
+        maxT = 1;
+    }
+    // Find the intersection of distance convex hull and fat line.
+    x_at(distance2y[0], distance2y[1], top, bottom, flags, minT, maxT);
+    x_at(distance2y[1], distance2y[2], top, bottom, flags, minT, maxT);
+    x_at(distance2y[2], distance2y[0], top, bottom, flags, minT, maxT);
+    
+    return minT < maxT; // returns false if distance shows no intersection
+}
+
diff --git a/experimental/Intersection/QuadraticBezierClip_Test.cpp b/experimental/Intersection/QuadraticBezierClip_Test.cpp
new file mode 100644
index 0000000..4fe35ba
--- /dev/null
+++ b/experimental/Intersection/QuadraticBezierClip_Test.cpp
@@ -0,0 +1,10 @@
+/*
+ *  QuadraticBezierClip_Test.cpp
+ *  edge
+ *
+ *  Created by Cary Clark on 1/10/12.
+ *  Copyright 2012 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+
diff --git a/experimental/Intersection/QuadraticIntersection.cpp b/experimental/Intersection/QuadraticIntersection.cpp
new file mode 100644
index 0000000..9a92c69
--- /dev/null
+++ b/experimental/Intersection/QuadraticIntersection.cpp
@@ -0,0 +1,188 @@
+#include "CubicIntersection.h"
+#include "Intersections.h"
+#include "IntersectionUtilities.h"
+#include "LineIntersection.h"
+
+class QuadraticIntersections : public Intersections {
+public:
+
+QuadraticIntersections(const Quadratic& q1, const Quadratic& q2, Intersections& i) 
+    : quad1(q1)
+    , quad2(q2)
+    , intersections(i)
+    , depth(0) 
+    , splits(0) {
+}
+
+bool intersect() {
+    double minT1, minT2, maxT1, maxT2;
+    if (!bezier_clip(quad2, quad1, minT1, maxT1)) {
+        return false;
+    }
+    if (!bezier_clip(quad1, quad2, minT2, maxT2)) {
+        return false;
+    }
+    int split;
+    if (maxT1 - minT1 < maxT2 - minT2) {
+        intersections.swap();
+        minT2 = 0;
+        maxT2 = 1;
+        split = maxT1 - minT1 > tClipLimit;
+    } else {
+        minT1 = 0;
+        maxT1 = 1;
+        split = (maxT2 - minT2 > tClipLimit) << 1;
+    }
+    return chop(minT1, maxT1, minT2, maxT2, split);
+}
+
+protected:
+        
+bool intersect(double minT1, double maxT1, double minT2, double maxT2) {
+    Quadratic smaller, larger;
+    // FIXME: carry last subdivide and reduceOrder result with quad 
+    sub_divide(quad1, minT1, maxT1, intersections.swapped() ? larger : smaller);
+    sub_divide(quad2, minT2, maxT2, intersections.swapped() ? smaller : larger);
+    Quadratic smallResult;
+    if (reduceOrder(smaller, smallResult) <= 2) {
+        Quadratic largeResult;
+        if (reduceOrder(larger, largeResult) <= 2) {
+            _Point pt;
+            const _Line& smallLine = (const _Line&) smallResult;
+            const _Line& largeLine = (const _Line&) largeResult;
+            if (!lineIntersect(smallLine, largeLine, &pt)) {
+                return false;
+            }
+            double smallT = t_at(smallLine, pt);
+            double largeT = t_at(largeLine, pt);
+            if (intersections.swapped()) {
+                smallT = interp(minT2, maxT2, smallT); 
+                largeT = interp(minT1, maxT1, largeT); 
+            } else {
+                smallT = interp(minT1, maxT1, smallT); 
+                largeT = interp(minT2, maxT2, largeT); 
+            }
+            intersections.add(smallT, largeT);
+            return true;
+        }
+    }
+    double minT, maxT;
+    if (!bezier_clip(smaller, larger, minT, maxT)) {
+        if (minT == maxT) {
+            if (intersections.swapped()) {
+                minT1 = (minT1 + maxT1) / 2;
+                minT2 = interp(minT2, maxT2, minT);
+            } else {
+                minT1 = interp(minT1, maxT1, minT);
+                minT2 = (minT2 + maxT2) / 2;
+            }
+            intersections.add(minT1, minT2);
+            return true;
+        }
+        return false;
+    }
+    
+    int split;
+    if (intersections.swapped()) {
+        double newMinT1 = interp(minT1, maxT1, minT);
+        double newMaxT1 = interp(minT1, maxT1, maxT);
+        split = (newMaxT1 - newMinT1 > (maxT1 - minT1) * tClipLimit) << 1;
+        printf("%s d=%d s=%d new1=(%g,%g) old1=(%g,%g) split=%d\n", __FUNCTION__, depth,
+            splits, newMinT1, newMaxT1, minT1, maxT1, split);
+        minT1 = newMinT1;
+        maxT1 = newMaxT1;
+    } else {
+        double newMinT2 = interp(minT2, maxT2, minT);
+        double newMaxT2 = interp(minT2, maxT2, maxT);
+        split = newMaxT2 - newMinT2 > (maxT2 - minT2) * tClipLimit;
+        printf("%s d=%d s=%d new2=(%g,%g) old2=(%g,%g) split=%d\n", __FUNCTION__, depth,
+            splits, newMinT2, newMaxT2, minT2, maxT2, split);
+        minT2 = newMinT2;
+        maxT2 = newMaxT2;
+    }
+    return chop(minT1, maxT1, minT2, maxT2, split);
+}
+
+bool chop(double minT1, double maxT1, double minT2, double maxT2, int split) {
+    ++depth;
+    intersections.swap();
+    if (split) {
+        ++splits;
+        if (split & 2) {
+            double middle1 = (maxT1 + minT1) / 2;
+            intersect(minT1, middle1, minT2, maxT2);
+            intersect(middle1, maxT1, minT2, maxT2);
+        } else {
+            double middle2 = (maxT2 + minT2) / 2;
+            intersect(minT1, maxT1, minT2, middle2);
+            intersect(minT1, maxT1, middle2, maxT2);
+        }
+        --splits;
+        intersections.swap();
+        --depth;
+        return intersections.intersected();
+    }
+    bool result = intersect(minT1, maxT1, minT2, maxT2);
+    intersections.swap();
+    --depth;
+    return result;
+}
+
+private:
+
+static const double tClipLimit = 0.8; // http://cagd.cs.byu.edu/~tom/papers/bezclip.pdf see Multiple intersections
+const Quadratic& quad1;
+const Quadratic& quad2;
+Intersections& intersections;
+int depth;
+int splits;
+};
+
+bool intersectStart(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
+    QuadraticIntersections q(q1, q2, i);
+    return q.intersect();
+}
+
+
+// Another approach is to start with the implicit form of one curve and solve
+// by substituting in the parametric form of the other.
+// The downside of this approach is that early rejects are difficult to come by.
+// http://planetmath.org/encyclopedia/GaloisTheoreticDerivationOfTheQuarticFormula.html#step
+/*
+given x^4 + ax^3 + bx^2 + cx + d
+the resolvent cubic is x^3 - 2bx^2 + (b^2 + ac - 4d)x + (c^2 + a^2d - abc)
+use the cubic formula (CubicRoots.cpp) to find the radical expressions t1, t2, and t3.
+
+(x - r1 r2) (x - r3 r4) = x^2 - (t2 + t3 - t1) / 2 x + d
+s = r1*r2 = ((t2 + t3 - t1) + sqrt((t2 + t3 - t1)^2 - 16*d)) / 4
+t = r3*r4 = ((t2 + t3 - t1) - sqrt((t2 + t3 - t1)^2 - 16*d)) / 4
+
+u = r1+r2 = (-a + sqrt(a^2 - 4*t1)) / 2
+v = r3+r4 = (-a - sqrt(a^2 - 4*t1)) / 2
+
+r1 = (u + sqrt(u^2 - 4*s)) / 2
+r2 = (u - sqrt(u^2 - 4*s)) / 2
+r3 = (v + sqrt(v^2 - 4*t)) / 2
+r4 = (v - sqrt(v^2 - 4*t)) / 2
+*/
+
+
+/* square root of complex number
+http://en.wikipedia.org/wiki/Square_root#Square_roots_of_negative_and_complex_numbers
+Algebraic formula
+When the number is expressed using Cartesian coordinates the following formula
+ can be used for the principal square root:[5][6]
+
+sqrt(x + iy) = sqrt((r + x) / 2) +/- i*sqrt((r - x) / 2)
+
+where the sign of the imaginary part of the root is taken to be same as the sign
+ of the imaginary part of the original number, and
+ 
+r = abs(x + iy) = sqrt(x^2 + y^2)
+
+is the absolute value or modulus of the original number. The real part of the 
+principal value is always non-negative.
+The other square root is simply –1 times the principal square root; in other 
+words, the two square roots of a number sum to 0.
+ */
+ 
\ No newline at end of file
diff --git a/experimental/Intersection/QuadraticIntersection_TestData.cpp b/experimental/Intersection/QuadraticIntersection_TestData.cpp
new file mode 100644
index 0000000..0d876bf
--- /dev/null
+++ b/experimental/Intersection/QuadraticIntersection_TestData.cpp
@@ -0,0 +1,11 @@
+/*
+ *  QuadraticIntersection_TestData.cpp
+ *  edge
+ *
+ *  Created by Cary Clark on 1/10/12.
+ *  Copyright 2012 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "QuadraticIntersection_TestData.h"
+
diff --git a/experimental/Intersection/QuadraticIntersection_TestData.h b/experimental/Intersection/QuadraticIntersection_TestData.h
new file mode 100644
index 0000000..003c376
--- /dev/null
+++ b/experimental/Intersection/QuadraticIntersection_TestData.h
@@ -0,0 +1,9 @@
+/*
+ *  QuadraticIntersection_TestData.h
+ *  edge
+ *
+ *  Created by Cary Clark on 1/10/12.
+ *  Copyright 2012 __MyCompanyName__. All rights reserved.
+ *
+ */
+
diff --git a/experimental/Intersection/QuadraticParameterization.cpp b/experimental/Intersection/QuadraticParameterization.cpp
new file mode 100644
index 0000000..a78aa9f
--- /dev/null
+++ b/experimental/Intersection/QuadraticParameterization.cpp
@@ -0,0 +1,132 @@
+#include "CubicIntersection.h"
+#include "QuadraticUtilities.h"
+
+/* from http://tom.cs.byu.edu/~tom/papers/cvgip84.pdf 4.1
+ *
+ * This paper proves that Syvester's method can compute the implicit form of 
+ * the quadratic from the parameterzied form.
+ *
+ * Given x = a*t*t + b*t + c  (the parameterized form)
+ *       y = d*t*t + e*t + f
+ *
+ * we want to find an equation of the implicit form:
+ *
+ * A*x*x + B*x*y + C*y*y + D*x + E*y + F = 0
+ *
+ * The implicit form can be expressed as a 4x4 determinant, as shown.
+ *
+ * The resultant obtained by Syvester's method is
+ *
+ * |   a   b   (c - x)     0     |
+ * |   0   a      b     (c - x)  |
+ * |   d   e   (f - y)     0     |
+ * |   0   d      e     (f - y)  |
+ *
+ * which expands to
+ *
+ * d*d*x*x + -2*a*d*x*y + a*a*y*y
+ *         + (-2*c*d*d + b*e*d - a*e*e + 2*a*f*d)*x
+ *         + (-2*f*a*a + e*b*a - d*b*b + 2*d*c*a)*y
+ *         +
+ * |   a   b   c   0   |
+ * |   0   a   b   c   | == 0.
+ * |   d   e   f   0   |
+ * |   0   d   e   f   |
+ *
+ * Expanding the constant determinant results in
+ *
+ *   | a b c |     | b c 0 |
+ * a*| e f 0 | + d*| a b c | ==
+ *   | d e f |     | d e f |
+ *
+ * a*(a*f*f + c*e*e - c*f*d - b*e*f) + d*(b*b*f + c*c*d - c*a*f - c*e*b)
+ *
+ */
+
+enum {
+    xx_coeff,
+    xy_coeff,
+    yy_coeff,
+    x_coeff,
+    y_coeff,
+    c_coeff,
+    coeff_count
+};
+
+static bool straight_forward = true;
+
+static void implicit_coefficients(const Quadratic& q, double p[coeff_count]) {
+    double a, b, c;
+    set_abc(&q[0].x, a, b, c);
+    double d, e, f;
+    set_abc(&q[0].y, d, e, f);
+    // compute the implicit coefficients
+    if (straight_forward) { // 42 muls, 13 adds
+        p[xx_coeff] = d * d;
+        p[xy_coeff] = -2 * a * d;
+        p[yy_coeff] = a * a;
+        p[x_coeff] = -2*c*d*d + b*e*d - a*e*e + 2*a*f*d;
+        p[y_coeff] = -2*f*a*a + e*b*a - d*b*b + 2*d*c*a;
+        p[c_coeff] = a*(a*f*f + c*e*e - c*f*d - b*e*f)
+                   + d*(b*b*f + c*c*d - c*a*f - c*e*b);
+    } else { // 26 muls, 11 adds
+        double aa = a * a;
+        double ad = a * d;
+        double dd = d * d;
+        p[xx_coeff] = dd;
+        p[xy_coeff] = -2 * ad;
+        p[yy_coeff] = aa;
+        double be = b * e;
+        double bde = be * d;
+        double cdd = c * dd;
+        double ee = e * e;
+        p[x_coeff] =  -2*cdd + bde - a*ee + 2*ad*f;
+        double aaf = aa * f;
+        double abe = a * be;
+        double ac = a * c;
+        double bb_2ac = b*b - 2*ac;
+        p[y_coeff] = -2*aaf + abe - d*bb_2ac;
+        p[c_coeff] = aaf*f + ac*ee + d*f*bb_2ac - abe*f + c*cdd - c*bde;
+    }
+}
+
+ /* Given a pair of quadratics, determine their parametric coefficients.
+  * If the scaled coefficients are nearly equal, then the part of the quadratics
+  * may be coincident.
+  * FIXME: optimization -- since comparison short-circuits on no match,
+  * lazily compute the coefficients, comparing the easiest to compute first.
+  * xx and yy first; then xy; and so on.
+  */
+bool implicit_matches(const Quadratic& one, const Quadratic& two) {
+    double p1[coeff_count]; // a'xx , b'xy , c'yy , d'x , e'y , f
+    double p2[coeff_count];
+    implicit_coefficients(one, p1);
+    implicit_coefficients(two, p2);
+    int first = 0;
+    for (int index = 0; index < coeff_count; ++index) {
+        if (approximately_zero(p1[index]) || approximately_zero(p2[index])) {
+            first += first == index;
+            continue;
+        }
+        if (first == index) {
+            continue;
+        }
+        if (!approximately_equal(p1[index] * p2[first],
+                p1[first] * p2[index])) {
+            return false;
+        }
+    }
+    return true;
+}
+
+static double tangent(const double* quadratic, double t) {
+    double a, b, c;
+    set_abc(quadratic, a, b, c);
+    return 2 * a * t + b;
+}
+
+void tangent(const Quadratic& quadratic, double t, _Point& result) {
+    result.x = tangent(&quadratic[0].x, t);
+    result.y = tangent(&quadratic[0].y, t);
+}
+
diff --git a/experimental/Intersection/QuadraticParameterization_Test.cpp b/experimental/Intersection/QuadraticParameterization_Test.cpp
new file mode 100644
index 0000000..c796685
--- /dev/null
+++ b/experimental/Intersection/QuadraticParameterization_Test.cpp
@@ -0,0 +1,37 @@
+#include "CubicIntersection.h"
+#include "CubicIntersection_Tests.h"
+
+const Quadratic quadratics[] = {
+    {{0, 0}, {1, 0}, {1, 1}},
+};
+
+const size_t quadratics_count = sizeof(quadratics) / sizeof(quadratics[0]);
+
+int firstQuadraticCoincidenceTest = 0;
+
+void QuadraticCoincidence_Test() {
+    // split large quadratic
+    // compare original, parts, to see if the are coincident
+    for (size_t index = firstQuadraticCoincidenceTest; index < quadratics_count; ++index) {
+        const Quadratic& test = quadratics[index];
+        QuadraticPair split;
+        chop_at(test, split, 0.5);
+        Quadratic midThird;
+        sub_divide(test, 1.0/3, 2.0/3, midThird);
+        if (!implicit_matches(test, split.first())) {
+            printf("%s-1 %d", __FUNCTION__, (int)index);
+        }
+        if (!implicit_matches(test, split.second())) {
+            printf("%s-2 %d", __FUNCTION__, (int)index);
+        }
+        if (!implicit_matches(midThird, split.first())) {
+            printf("%s-3 %d", __FUNCTION__, (int)index);
+        }
+        if (!implicit_matches(midThird, split.second())) {
+            printf("%s-4 %d", __FUNCTION__, (int)index);
+        }
+        if (!implicit_matches(split.first(), split.second())) {
+            printf("%s-5 %d", __FUNCTION__, (int)index);
+        }
+    }
+}
diff --git a/experimental/Intersection/QuadraticReduceOrder.cpp b/experimental/Intersection/QuadraticReduceOrder.cpp
new file mode 100644
index 0000000..02d1956
--- /dev/null
+++ b/experimental/Intersection/QuadraticReduceOrder.cpp
@@ -0,0 +1,164 @@
+#include "CubicIntersection.h"
+#include "Extrema.h"
+#include "IntersectionUtilities.h"
+#include "LineParameters.h"
+
+static double interp_quad_coords(double a, double b, double c, double t)
+{
+    double ab = interp(a, b, t);
+    double bc = interp(b, c, t);
+    return interp(ab, bc, t);
+}
+
+static int coincident_line(const Quadratic& quad, Quadratic& reduction) {
+    reduction[0] = reduction[1] = quad[0];
+    return 1;
+}
+
+static int vertical_line(const Quadratic& quad, Quadratic& reduction) {
+    double tValue;
+    reduction[0] = quad[0];
+    reduction[1] = quad[2];
+    int smaller = reduction[1].y > reduction[0].y;
+    int larger = smaller ^ 1;
+    if (SkFindQuadExtrema(quad[0].y, quad[1].y, quad[2].y, &tValue)) {
+        double yExtrema = interp_quad_coords(quad[0].y, quad[1].y, quad[2].y, tValue);
+        if (reduction[smaller].y > yExtrema) {
+            reduction[smaller].y = yExtrema;
+        } else if (reduction[larger].y < yExtrema) {
+            reduction[larger].y = yExtrema;
+        }
+    }
+    return 2;
+}
+
+static int horizontal_line(const Quadratic& quad, Quadratic& reduction) {
+    double tValue;
+    reduction[0] = quad[0];
+    reduction[1] = quad[2];
+    int smaller = reduction[1].x > reduction[0].x;
+    int larger = smaller ^ 1;
+    if (SkFindQuadExtrema(quad[0].x, quad[1].x, quad[2].x, &tValue)) {
+        double xExtrema = interp_quad_coords(quad[0].x, quad[1].x, quad[2].x, tValue);
+        if (reduction[smaller].x > xExtrema) {
+            reduction[smaller].x = xExtrema;
+        }  else if (reduction[larger].x < xExtrema) {
+            reduction[larger].x = xExtrema;
+        }
+    }
+    return 2;
+}
+
+static int check_linear(const Quadratic& quad, Quadratic& reduction,
+        int minX, int maxX, int minY, int maxY) {
+    int startIndex = 0;
+    int endIndex = 2;
+    while (quad[startIndex].approximatelyEqual(quad[endIndex])) {
+        --endIndex;
+        if (endIndex == 0) {
+            printf("%s shouldn't get here if all four points are about equal", __FUNCTION__);
+            assert(0);
+        }
+    }
+    LineParameters lineParameters;
+    lineParameters.quadEndPoints(quad, startIndex, endIndex);
+    double normalSquared = lineParameters.normalSquared();
+    double distance = lineParameters.controlPtDistance(quad); // not normalized
+    double limit = normalSquared * SquaredEpsilon;
+    double distSq = distance * distance;
+    if (distSq > limit) {
+        return 0;
+    }
+    // four are colinear: return line formed by outside
+    reduction[0] = quad[0];
+    reduction[1] = quad[2];
+    int sameSide;
+    bool useX = quad[maxX].x - quad[minX].x >= quad[maxY].y - quad[minY].y;
+    if (useX) {
+        sameSide = sign(quad[0].x - quad[1].x) + sign(quad[2].x - quad[1].x);
+    } else {
+        sameSide = sign(quad[0].y - quad[1].y) + sign(quad[2].y - quad[1].y);
+    }
+    if ((sameSide & 3) != 2) {
+        return 2;
+    }
+    double tValue;
+    int root;
+    if (useX) {
+        root = SkFindQuadExtrema(quad[0].x, quad[1].x, quad[2].x, &tValue);
+    } else {
+        root = SkFindQuadExtrema(quad[0].y, quad[1].y, quad[2].y, &tValue);
+    }
+    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);
+        // sameSide > 0 means mid is smaller than either [0] or [2], so replace smaller
+        int replace;
+        if (useX) {
+            if (extrema.x < quad[0].x ^ extrema.x < quad[2].x) {
+                return 2;
+            }
+            replace = (extrema.x < quad[0].x | extrema.x < quad[2].x)
+                    ^ quad[0].x < quad[2].x;
+        } else {
+            if (extrema.y < quad[0].y ^ extrema.y < quad[2].y) {
+                return 2;
+            }
+            replace = (extrema.y < quad[0].y | extrema.y < quad[2].y)
+                    ^ quad[0].y < quad[2].y;
+        }
+        reduction[replace] = extrema;
+    }
+    return 2;
+}
+
+// reduce to a quadratic or smaller
+// look for identical points
+// look for all four points in a line 
+    // note that three points in a line doesn't simplify a cubic
+// look for approximation with single quadratic
+    // save approximation with multiple quadratics for later
+int reduceOrder(const Quadratic& quad, Quadratic& reduction) {
+    int index, minX, maxX, minY, maxY;
+    int minXSet, minYSet;
+    minX = maxX = minY = maxY = 0;
+    minXSet = minYSet = 0;
+    for (index = 1; index < 3; ++index) {
+        if (quad[minX].x > quad[index].x) {
+            minX = index;
+        }
+        if (quad[minY].y > quad[index].y) {
+            minY = index;
+        }
+        if (quad[maxX].x < quad[index].x) {
+            maxX = index;
+        }
+        if (quad[maxY].y < quad[index].y) {
+            maxY = index;
+        }
+    }
+    for (index = 0; index < 3; ++index) {
+        if (approximately_equal(quad[index].x, quad[minX].x)) {
+            minXSet |= 1 << index;
+        }
+        if (approximately_equal(quad[index].y, quad[minY].y)) {
+            minYSet |= 1 << index;
+        }
+    }
+    if (minXSet == 0xF) { // test for vertical line
+        if (minYSet == 0xF) { // return 1 if all four are coincident
+            return coincident_line(quad, reduction);
+        }
+        return vertical_line(quad, reduction);
+    }
+    if (minYSet == 0xF) { // test for horizontal line
+        return horizontal_line(quad, reduction);
+    }
+    int result = check_linear(quad, reduction, minX, maxX, minY, maxY);
+    if (result) {
+        return result;
+    }
+    memcpy(reduction, quad, sizeof(Quadratic));
+    return 3;
+}
diff --git a/experimental/Intersection/QuadraticSubDivide.cpp b/experimental/Intersection/QuadraticSubDivide.cpp
new file mode 100644
index 0000000..947f7d8
--- /dev/null
+++ b/experimental/Intersection/QuadraticSubDivide.cpp
@@ -0,0 +1,63 @@
+#include "CubicIntersection.h"
+#include "IntersectionUtilities.h"
+
+/*
+Given a quadratic q, t1, and t2, find a small quadratic segment.
+
+The new quadratic is defined by A, B, and C, where
+ A = c[0]*(1 - t1)*(1 - t1) + 2*c[1]*t1*(1 - t1) + c[2]*t1*t1
+ C = c[3]*(1 - t1)*(1 - t1) + 2*c[2]*t1*(1 - t1) + c[1]*t1*t1
+
+To find B, compute the point halfway between t1 and t2:
+
+q(at (t1 + t2)/2) == D
+
+Next, compute where D must be if we know the value of B:
+
+_12 = A/2 + B/2
+12_ = B/2 + C/2
+123 = A/4 + B/2 + C/4
+    = D
+    
+Group the known values on one side:
+
+B   = D*2 - A/2 - C/2
+*/
+
+static double interp_quad_coords(const double* src, double t)
+{
+    double ab = interp(src[0], src[2], t);
+    double bc = interp(src[2], src[4], t);
+    double abc = interp(ab, bc, t);
+    return abc;
+}
+
+void sub_divide(const Quadratic& src, double t1, double t2, Quadratic& dst) {
+    double ax = dst[0].x = interp_quad_coords(&src[0].x, t1);
+    double ay = dst[0].y = interp_quad_coords(&src[0].y, t1);
+    double dx = interp_quad_coords(&src[0].x, (t1 + t2) / 2);
+    double dy = interp_quad_coords(&src[0].y, (t1 + t2) / 2);
+    double cx = dst[2].x = interp_quad_coords(&src[0].x, t2);
+    double cy = dst[2].y = interp_quad_coords(&src[0].y, t2);
+    /* bx = */ dst[1].x = 2*dx - (ax + cx)/2;
+    /* by = */ dst[1].y = 2*dy - (ay + cy)/2;
+}
+
+/* classic one t subdivision */
+static void interp_quad_coords(const double* src, double* dst, double t)
+{
+    double ab = interp(src[0], src[2], t);
+    double bc = interp(src[2], src[4], t);
+
+    dst[0] = src[0];
+    dst[2] = ab;
+    dst[4] = interp(ab, bc, t);
+    dst[6] = bc;
+    dst[8] = src[4];
+}
+
+void chop_at(const Quadratic& src, QuadraticPair& dst, double t)
+{
+    interp_quad_coords(&src[0].x, &dst.pts[0].x, t);
+    interp_quad_coords(&src[0].y, &dst.pts[0].y, t);
+}
diff --git a/experimental/Intersection/QuadraticUtilities.h b/experimental/Intersection/QuadraticUtilities.h
new file mode 100644
index 0000000..46f6745
--- /dev/null
+++ b/experimental/Intersection/QuadraticUtilities.h
@@ -0,0 +1,14 @@
+/* Parameterization form, given A*t*t + 2*B*t*(1-t) + C*(1-t)*(1-t)
+ *
+ * a = A - 2*B +   C
+ * b =     2*B - 2*C
+ * c =             C
+ */
+static void set_abc(const double* quad, double& a, double& b, double& c) {
+    a = quad[0];     // a = A
+    b = 2 * quad[2]; // b =     2*B
+    c = quad[4];     // c =             C
+    b -= c;          // b =     2*B -   C
+    a -= b;          // a = A - 2*B +   C
+    b -= c;          // b =     2*B - 2*C
+}
diff --git a/experimental/Intersection/ReduceOrder.cpp b/experimental/Intersection/ReduceOrder.cpp
new file mode 100644
index 0000000..e742a26
--- /dev/null
+++ b/experimental/Intersection/ReduceOrder.cpp
@@ -0,0 +1,238 @@
+#include "CubicIntersection.h"
+#include "Extrema.h"
+#include "IntersectionUtilities.h"
+#include "LineParameters.h"
+
+#ifdef MAYBE_USEFUL_IN_THE_FUTURE
+static double interp_quad_coords(double a, double b, double c, double t)
+{
+    double ab = interp(a, b, t);
+    double bc = interp(b, c, t);
+    return interp(ab, bc, t);
+}
+#endif
+
+static double interp_cubic_coords(const double* src, double t)
+{
+    double ab = interp(src[0], src[2], t);
+    double bc = interp(src[2], src[4], t);
+    double cd = interp(src[4], src[6], t);
+    double abc = interp(ab, bc, t);
+    double bcd = interp(bc, cd, t);
+    return interp(abc, bcd, t);
+}
+
+static int coincident_line(const Cubic& cubic, Cubic& reduction) {
+    reduction[0] = reduction[1] = cubic[0];
+    return 1;
+}
+
+static int vertical_line(const Cubic& cubic, Cubic& reduction) {
+    double tValues[2];
+    reduction[0] = cubic[0];
+    reduction[1] = cubic[3];
+    int smaller = reduction[1].y > reduction[0].y;
+    int larger = smaller ^ 1;
+    int roots = SkFindCubicExtrema(cubic[0].y, cubic[1].y, cubic[2].y, cubic[3].y, tValues);
+    for (int index = 0; index < roots; ++index) {
+        double yExtrema = interp_cubic_coords(&cubic[0].y, tValues[index]);
+        if (reduction[smaller].y > yExtrema) {
+            reduction[smaller].y = yExtrema;
+            continue;
+        } 
+        if (reduction[larger].y < yExtrema) {
+            reduction[larger].y = yExtrema;
+        }
+    }
+    return 2;
+}
+
+static int horizontal_line(const Cubic& cubic, Cubic& reduction) {
+    double tValues[2];
+    reduction[0] = cubic[0];
+    reduction[1] = cubic[3];
+    int smaller = reduction[1].x > reduction[0].x;
+    int larger = smaller ^ 1;
+    int roots = SkFindCubicExtrema(cubic[0].x, cubic[1].x, cubic[2].x, cubic[3].x, tValues);
+    for (int index = 0; index < roots; ++index) {
+        double xExtrema = interp_cubic_coords(&cubic[0].x, tValues[index]);
+        if (reduction[smaller].x > xExtrema) {
+            reduction[smaller].x = xExtrema;
+            continue;
+        } 
+        if (reduction[larger].x < xExtrema) {
+            reduction[larger].x = xExtrema;
+        }
+    }
+    return 2;
+}
+
+// check to see if it is a quadratic or a line
+static int check_quadratic(const Cubic& cubic, Cubic& reduction,
+        int minX, int maxX, int minY, int maxY) {
+    double dx10 = cubic[1].x - cubic[0].x;
+    double dx23 = cubic[2].x - cubic[3].x;
+    double midX = cubic[0].x + dx10 * 3 / 2;
+    if (!approximately_equal(midX - cubic[3].x, dx23 * 3 / 2)) {
+        return 0;
+    }
+    double dy10 = cubic[1].y - cubic[0].y;
+    double dy23 = cubic[2].y - cubic[3].y;
+    double midY = cubic[0].y + dy10 * 3 / 2;
+    if (!approximately_equal(midY - cubic[3].y, dy23 * 3 / 2)) {
+        return 0;
+    }
+    reduction[0] = cubic[0];
+    reduction[1].x = midX;
+    reduction[1].y = midY;
+    reduction[2] = cubic[3];
+    return 3;
+}
+
+static int check_linear(const Cubic& cubic, Cubic& reduction,
+        int minX, int maxX, int minY, int maxY) {
+    int startIndex = 0;
+    int endIndex = 3;
+    while (cubic[startIndex].approximatelyEqual(cubic[endIndex])) {
+        --endIndex;
+        if (endIndex == 0) {
+            printf("%s shouldn't get here if all four points are about equal", __FUNCTION__);
+            assert(0);
+        }
+    }
+    LineParameters lineParameters;
+    lineParameters.cubicEndPoints(cubic, startIndex, endIndex);
+    double normalSquared = lineParameters.normalSquared();
+    double distance[2]; // distance is not normalized
+    int mask = other_two(startIndex, endIndex);
+    int inner1 = startIndex ^ mask;
+    int inner2 = endIndex ^ mask;
+    lineParameters.controlPtDistance(cubic, inner1, inner2, distance);
+    double limit = normalSquared * SquaredEpsilon;
+    int index;
+    for (index = 0; index < 2; ++index) {
+        double distSq = distance[index];
+        distSq *= distSq;
+        if (distSq > limit) {
+            return 0;
+        }
+    }
+    // four are colinear: return line formed by outside
+    reduction[0] = cubic[0];
+    reduction[1] = cubic[3];
+    int sameSide1;
+    int sameSide2;
+    bool useX = cubic[maxX].x - cubic[minX].x >= cubic[maxY].y - cubic[minY].y;
+    if (useX) {
+        sameSide1 = sign(cubic[0].x - cubic[1].x) + sign(cubic[3].x - cubic[1].x);
+        sameSide2 = sign(cubic[0].x - cubic[2].x) + sign(cubic[3].x - cubic[2].x);
+    } else {
+        sameSide1 = sign(cubic[0].y - cubic[1].y) + sign(cubic[3].y - cubic[1].y);
+        sameSide2 = sign(cubic[0].y - cubic[2].y) + sign(cubic[3].y - cubic[2].y);
+    }
+    if (sameSide1 == sameSide2 && (sameSide1 & 3) != 2) {
+        return 2;
+    }
+    double tValues[2];
+    int roots;
+    if (useX) {
+        roots = SkFindCubicExtrema(cubic[0].x, cubic[1].x, cubic[2].x, cubic[3].x, tValues);
+    } else {
+        roots = SkFindCubicExtrema(cubic[0].y, cubic[1].y, cubic[2].y, cubic[3].y, tValues);
+    }
+    for (index = 0; index < roots; ++index) {
+        _Point extrema;
+        extrema.x = interp_cubic_coords(&cubic[0].x, tValues[index]);
+        extrema.y = interp_cubic_coords(&cubic[0].y, tValues[index]);
+        // sameSide > 0 means mid is smaller than either [0] or [3], so replace smaller
+        int replace;
+        if (useX) {
+            if (extrema.x < cubic[0].x ^ extrema.x < cubic[3].x) {
+                continue;
+            }
+            replace = (extrema.x < cubic[0].x | extrema.x < cubic[3].x)
+                    ^ cubic[0].x < cubic[3].x;
+        } else {
+            if (extrema.y < cubic[0].y ^ extrema.y < cubic[3].y) {
+                continue;
+            }
+            replace = (extrema.y < cubic[0].y | extrema.y < cubic[3].y)
+                    ^ cubic[0].y < cubic[3].y;
+        }
+        reduction[replace] = extrema;
+    }
+    return 2;
+}
+
+/* food for thought:
+http://objectmix.com/graphics/132906-fast-precision-driven-cubic-quadratic-piecewise-degree-reduction-algos-2-a.html
+
+Given points c1, c2, c3 and c4 of a cubic Bezier, the points of the
+corresponding quadratic Bezier are (given in convex combinations of
+points):
+
+q1 = (11/13)c1 + (3/13)c2 -(3/13)c3 + (2/13)c4
+q2 = -c1 + (3/2)c2 + (3/2)c3 - c4
+q3 = (2/13)c1 - (3/13)c2 + (3/13)c3 + (11/13)c4
+
+Of course, this curve does not interpolate the end-points, but it would
+be interesting to see the behaviour of such a curve in an applet.
+
+--
+Kalle Rutanen
+http://kaba.hilvi.org
+
+*/
+
+// reduce to a quadratic or smaller
+// look for identical points
+// look for all four points in a line 
+    // note that three points in a line doesn't simplify a cubic
+// look for approximation with single quadratic
+    // save approximation with multiple quadratics for later
+int reduceOrder(const Cubic& cubic, Cubic& reduction, ReduceOrder_Flags allowQuadratics) {
+    int index, minX, maxX, minY, maxY;
+    int minXSet, minYSet;
+    minX = maxX = minY = maxY = 0;
+    minXSet = minYSet = 0;
+    for (index = 1; index < 4; ++index) {
+        if (cubic[minX].x > cubic[index].x) {
+            minX = index;
+        }
+        if (cubic[minY].y > cubic[index].y) {
+            minY = index;
+        }
+        if (cubic[maxX].x < cubic[index].x) {
+            maxX = index;
+        }
+        if (cubic[maxY].y < cubic[index].y) {
+            maxY = index;
+        }
+    }
+    for (index = 0; index < 4; ++index) {
+        if (approximately_equal(cubic[index].x, cubic[minX].x)) {
+            minXSet |= 1 << index;
+        }
+        if (approximately_equal(cubic[index].y, cubic[minY].y)) {
+            minYSet |= 1 << index;
+        }
+    }
+    if (minXSet == 0xF) { // test for vertical line
+        if (minYSet == 0xF) { // return 1 if all four are coincident
+            return coincident_line(cubic, reduction);
+        }
+        return vertical_line(cubic, reduction);
+    }
+    if (minYSet == 0xF) { // test for horizontal line
+        return horizontal_line(cubic, reduction);
+    }
+    int result = check_linear(cubic, reduction, minX, maxX, minY, maxY);
+    if (result) {
+        return result;
+    }
+    if (allowQuadratics && (result = check_quadratic(cubic, reduction, minX, maxX, minY, maxY))) {
+        return result;
+    }
+    memcpy(reduction, cubic, sizeof(Cubic));
+    return 4;
+}
diff --git a/experimental/Intersection/ReduceOrder_Test.cpp b/experimental/Intersection/ReduceOrder_Test.cpp
new file mode 100644
index 0000000..04aea1a
--- /dev/null
+++ b/experimental/Intersection/ReduceOrder_Test.cpp
@@ -0,0 +1,137 @@
+#include "CubicIntersection.h"
+#include "CubicIntersection_TestData.h"
+#include "CubicIntersection_Tests.h"
+#include "TestUtilities.h"
+
+void ReduceOrder_Test() {
+    size_t index;
+    Cubic reduce;
+    int order;
+    enum {
+        RunAll,
+        RunPointDegenerates,
+        RunNotPointDegenerates,
+        RunLines,
+        RunNotLines,
+        RunModEpsilonLines,
+        RunLessEpsilonLines,
+        RunNegEpsilonLines,
+        RunQuadraticLines,
+        RunQuadraticModLines,
+        RunComputedLines,
+        RunNone
+    } run = RunAll;
+    int firstTestIndex = 0;
+#if 0
+    run = RunComputedLines;
+    firstTestIndex = 18;
+#endif
+    int firstPointDegeneratesTest = run == RunAll ? 0 : run == RunPointDegenerates ? firstTestIndex : INT_MAX;
+    int firstNotPointDegeneratesTest = run == RunAll ? 0 : run == RunNotPointDegenerates ? firstTestIndex : INT_MAX;
+    int firstLinesTest = run == RunAll ? 0 : run == RunLines ? firstTestIndex : INT_MAX;
+    int firstNotLinesTest = run == RunAll ? 0 : run == RunNotLines ? firstTestIndex : INT_MAX;
+    int firstModEpsilonTest = run == RunAll ? 0 : run == RunModEpsilonLines ? firstTestIndex : INT_MAX;
+    int firstLessEpsilonTest = run == RunAll ? 0 : run == RunLessEpsilonLines ? firstTestIndex : INT_MAX;
+    int firstNegEpsilonTest = run == RunAll ? 0 : run == RunNegEpsilonLines ? firstTestIndex : INT_MAX;
+    int firstQuadraticLineTest = run == RunAll ? 0 : run == RunQuadraticLines ? firstTestIndex : INT_MAX;
+    int firstQuadraticModLineTest = run == RunAll ? 0 : run == RunQuadraticModLines ? firstTestIndex : INT_MAX;
+    int firstComputedLinesTest = run == RunAll ? 0 : run == RunComputedLines ? firstTestIndex : INT_MAX;
+    
+    for (index = firstPointDegeneratesTest; index < pointDegenerates_count; ++index) {
+        const Cubic& cubic = pointDegenerates[index];
+        order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
+        if (order != 1) {
+            printf("[%d] pointDegenerates order=%d\n", (int) index, order);
+        }
+    }
+    for (index = firstNotPointDegeneratesTest; index < notPointDegenerates_count; ++index) {
+        const Cubic& cubic = notPointDegenerates[index];
+        order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
+        if (order == 1) {
+            printf("[%d] notPointDegenerates order=%d\n", (int) index, order);
+        }
+    }
+    for (index = firstLinesTest; index < lines_count; ++index) {
+        const Cubic& cubic = lines[index];
+        order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
+        if (order != 2) {
+            printf("[%d] lines order=%d\n", (int) index, order);
+        }
+    }
+    for (index = firstNotLinesTest; index < notLines_count; ++index) {
+        const Cubic& cubic = notLines[index];
+        order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
+        if (order == 2) {
+            printf("[%d] notLines order=%d\n", (int) index, order);
+        }
+    }
+    for (index = firstModEpsilonTest; index < modEpsilonLines_count; ++index) {
+        const Cubic& cubic = modEpsilonLines[index];
+        order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
+        if (order == 2) {
+            printf("[%d] line mod by epsilon order=%d\n", (int) index, order);
+        }
+    }
+    for (index = firstLessEpsilonTest; index < lessEpsilonLines_count; ++index) {
+        const Cubic& cubic = lessEpsilonLines[index];
+        order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
+        if (order != 2) {
+            printf("[%d] line less by epsilon/2 order=%d\n", (int) index, order);
+        }
+    }
+    for (index = firstNegEpsilonTest; index < negEpsilonLines_count; ++index) {
+        const Cubic& cubic = negEpsilonLines[index];
+        order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
+        if (order != 2) {
+            printf("[%d] line neg by epsilon/2 order=%d\n", (int) index, order);
+        }
+    }
+    for (index = firstQuadraticLineTest; index < quadraticLines_count; ++index) {
+        const Quadratic& quad = quadraticLines[index];
+        Cubic cubic;
+        quad_to_cubic(quad, cubic);
+        order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
+        if (order != 2) {
+            printf("[%d] line quad order=%d\n", (int) index, order);
+        }
+    }
+    for (index = firstQuadraticModLineTest; index < quadraticModEpsilonLines_count; ++index) {
+        const Quadratic& quad = quadraticModEpsilonLines[index];
+        Cubic cubic;
+        quad_to_cubic(quad, cubic);
+        order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
+        if (order != 3) {
+            printf("[%d] line mod quad order=%d\n", (int) index, order);
+        }
+    }
+    
+    // test if computed line end points are valid
+    for (index = firstComputedLinesTest; index < lines_count; ++index) {
+        const Cubic& cubic = lines[index];
+        bool controlsInside = controls_inside(cubic);        
+        order = reduceOrder(cubic, reduce, kReduceOrder_QuadraticsAllowed);
+        if (reduce[0].x == reduce[1].x && reduce[0].y == reduce[1].y) {
+            printf("[%d] line computed ends match order=%d\n", (int) index, order);
+        }
+        if (controlsInside) {
+            if (       reduce[0].x != cubic[0].x && reduce[0].x != cubic[3].x
+                    || reduce[0].y != cubic[0].y && reduce[0].y != cubic[3].y
+                    || reduce[1].x != cubic[0].x && reduce[1].x != cubic[3].x
+                    || reduce[1].y != cubic[0].y && reduce[1].y != cubic[3].y) {
+                printf("[%d] line computed ends order=%d\n", (int) index, order);
+            }
+        } else {
+            // binary search for extrema, compare against actual results
+                // while a control point is outside of bounding box formed by end points, split
+            _Rect bounds = {DBL_MAX, DBL_MAX, -DBL_MAX, -DBL_MAX};
+            find_tight_bounds(cubic, bounds);
+            if (       !approximately_equal(reduce[0].x, bounds.left) && !approximately_equal(reduce[0].x, bounds.right)
+                    || !approximately_equal(reduce[0].y, bounds.top) && !approximately_equal(reduce[0].y, bounds.bottom)
+                    || !approximately_equal(reduce[1].x, bounds.left) && !approximately_equal(reduce[1].x, bounds.right)
+                    || !approximately_equal(reduce[1].y, bounds.top) && !approximately_equal(reduce[1].y, bounds.bottom)) {
+                printf("[%d] line computed tight bounds order=%d\n", (int) index, order);
+            }
+            
+        }
+    }
+}
diff --git a/experimental/Intersection/SkAntiEdge.cpp b/experimental/Intersection/SkAntiEdge.cpp
new file mode 100644
index 0000000..7e22588
--- /dev/null
+++ b/experimental/Intersection/SkAntiEdge.cpp
@@ -0,0 +1,1081 @@
+/*
+ *  SkAntiEdge.cpp
+ *  core
+ *
+ *  Created by Cary Clark on 5/6/11.
+ *  Copyright 2011 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "SkAntiEdge.h"
+#include "SkPoint.h"
+
+void SkAntiEdge::pointOnLine(SkFixed x, SkFixed y) {
+    float x0 = SkFixedToFloat(x);
+    float y0 = SkFixedToFloat(y);
+    float x1 = SkFixedToFloat(fFirstX);
+    float y1 = SkFixedToFloat(fFirstY);
+    float x2 = SkFixedToFloat(fLastX);
+    float y2 = SkFixedToFloat(fLastY);
+    float numer = (x2 - x1) * (y1 - y0) - (x1 - x0) * (y2 - y1);
+    float denom = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1);
+    double dist = fabs(numer) / sqrt(denom);
+    SkAssertResult(dist < 0.01);
+}
+
+void SkAntiEdge::pointInLine(SkFixed x, SkFixed y) {
+    if (y == SK_MaxS32) {
+        return;
+    }
+    pointOnLine(x, y);
+    SkAssertResult(y >= fFirstY && y <= fLastY);
+}
+
+void SkAntiEdge::validate() {
+    pointOnLine(fWalkX, fY);
+    pointOnLine(fX, fWalkY);
+}
+
+bool SkAntiEdge::setLine(const SkPoint& p0, const SkPoint& p1) {
+    fFirstY = SkScalarToFixed(p0.fY);
+    fLastY = SkScalarToFixed(p1.fY);
+    if (fFirstY == fLastY) {
+        return false;
+    }   
+    fFirstX = SkScalarToFixed(p0.fX);
+    fLastX = SkScalarToFixed(p1.fX);
+    if (fFirstY > fLastY) {
+        SkTSwap(fFirstX, fLastX);
+        SkTSwap(fFirstY, fLastY);
+        fWinding = -1;
+    } else {
+        fWinding = 1;
+    }
+    SkFixed dx = fLastX - fFirstX;
+    fDXFlipped = dx < 0;
+    SkFixed dy = fLastY - fFirstY;
+    fDX = SkFixedDiv(dx, dy);
+    fDY = dx == 0 ? SK_MaxS32 : SkFixedDiv(dy, SkFixedAbs(dx));
+    fLink = NULL;
+    fLinkSet = false;
+    return true;
+}
+
+void SkAntiEdge::calcLine() {
+    SkFixed yStartFrac = SkFixedFraction(fFirstY);
+    if (fDXFlipped) {
+        SkFixed vert = SK_Fixed1 - yStartFrac; // distance from y start to x-axis
+        fX0 = fFirstX + SkFixedMul(fDX, vert);
+        SkFixed backupX = fFirstX + SkFixedMul(vert, fDX); // x cell to back up to
+        SkFixed cellX = SkIntToFixed(SkFixedFloor(backupX));
+        SkFixed endX = SkIntToFixed(SkFixedFloor(fLastX));
+        if (cellX < endX) {
+            cellX = endX;
+        }
+        SkFixed distX = fFirstX - cellX; // to y-axis
+        fY0 = fFirstY + SkFixedMul(fDY, distX);
+        SkFixed rowBottom = SkIntToFixed(SkFixedCeil(fFirstY + 1));
+        if (fLastY > rowBottom) {
+            fPartialY = 0;
+            fX = fX0;
+            fY = rowBottom;
+        } else {
+            fPartialY = SkFixedFraction(fLastY);
+            fX = fLastX;
+            fY = fLastY;
+        }
+    } else {
+        fPartialY = yStartFrac;
+        fX0 = fFirstX - SkFixedMul(fDX, yStartFrac);
+        fY0 = fFirstY;
+        if (fDY != SK_MaxS32) {
+            SkFixed xStartFrac = SkFixedFraction(fFirstX);
+            fY0 -= SkFixedMul(fDY, xStartFrac);
+        }
+        fX = fFirstX;
+        fY = fFirstY;
+    }
+    fWalkX = fX;
+    fWalkY = fY;
+    fFinished = false;
+}
+
+static SkFixed SkFixedAddPin(SkFixed a, SkFixed b) {
+    SkFixed result = a + b;
+    if (((a ^ ~b) & (a ^ result)) >= 0) { // one positive, one negative
+        return result;                    //  or all three same sign
+    }
+    return a < 0 ? -SK_FixedMax : SK_FixedMax;
+}
+
+// edge is increasing in x and y
+uint16_t SkAntiEdge::advanceX(SkFixed left) {
+    validate();
+    SkFixed x = SkFixedAddPin(fX0, fDX);
+    SkFixed wy = SkIntToFixed(SkFixedFloor(fWalkY + SK_Fixed1));
+    pointOnLine(x, wy);
+    SkFixed partial = SK_Fixed1 - fPartialY;
+    SkFixed bottomPartial = wy - fLastY;
+    if (bottomPartial > 0) {
+        partial -= bottomPartial;
+    }
+    if (x > fLastX) {
+        x = fLastX;
+        wy = fLastY;
+    }
+    uint16_t coverage;
+    if (left >= x) {
+        fFinished = true;
+        coverage = partial - 1; // walker is to the right of edge
+    } else {
+        SkFixed y = SkFixedAddPin(fY0, fDY);
+        SkFixed wx = SkIntToFixed(SkFixedFloor(fWalkX + SK_Fixed1));
+        if (fDY != SK_MaxS32) {
+            pointOnLine(wx, y);
+        }
+        if (y > fLastY) {
+            y = fLastY;
+            wx = fLastX;
+        }
+        bool topCorner = fWalkX <= fX;
+        bool bottomCorner = x <= wx;
+        bool halfPlane = !(topCorner ^ bottomCorner);
+        if (halfPlane) {
+            if (x - SkIntToFixed(SkFixedFloor(fX)) <= SK_Fixed1) {
+                coverage = ~((fX + x) >> 1); // avg of fx, fx+dx
+                fFinished = true;
+                if (x >= left + SK_Fixed1) {
+                    fWalkX = wx;
+                    fY = fY0 = y;
+                }
+            } else {
+                SkAssertResult(y - SkIntToFixed(SkFixedFloor(fY)) <= SK_Fixed1);
+                coverage = ((fY + y) >> 1);
+                fFinished = y == fLastY;
+                fWalkX = wx;
+                fY = fY0 = y;
+            }
+            coverage = coverage * partial >> 16;
+        } else if (topCorner) {
+            SkFixed xDiff = wx - fX;
+            SkAssertResult(xDiff >= 0);
+            SkAssertResult(xDiff <= SK_Fixed1);
+            SkFixed yDiff = y - fWalkY;
+            // This may be a very small negative number if error accumulates
+            // FIXME: for now, try setting it to zero in that case.
+            if (yDiff < 0) {
+                fX = fX0 = SkIntToFixed(SkFixedCeil(fX));
+                yDiff = 0;
+            }
+            SkAssertResult(yDiff >= 0);
+            SkAssertResult(yDiff <= SK_Fixed1);
+            int xCoverage = xDiff >> 1; // throw away 1 bit so multiply
+            int yCoverage = yDiff >> 1; //  stays in range
+            int triangle = xCoverage * yCoverage; // 30 bits
+            SkFixed bottomPartial = y - fLastY;
+            fFinished = bottomPartial >= 0;
+            if (fFinished) {
+                yCoverage = bottomPartial >> 1;
+                xCoverage = (wx - fLastX) >> 1;
+                triangle -= xCoverage * yCoverage;
+            }
+            coverage = triangle >> 15;
+            fWalkX = wx;
+            fY = fY0 = y;
+        } else {
+            SkAssertResult(bottomCorner);
+            SkFixed xDiff = x - fWalkX;
+            SkAssertResult(xDiff >= 0);
+            SkAssertResult(xDiff <= SK_Fixed1);
+            SkFixed yDiff = wy - fY;
+            SkAssertResult(yDiff >= 0);
+            SkAssertResult(yDiff <= SK_Fixed1);
+            int xCoverage = xDiff >> 1; // throw away 1 bit so multiply 
+            int yCoverage = yDiff >> 1; //  stays in range
+            int triangle = xCoverage * yCoverage >> 15;
+            coverage = partial - 1 - triangle;
+            fFinished = true;
+        }
+    }
+    validate();
+    return coverage;
+}
+
+// edge is increasing in x, but decreasing in y
+uint16_t SkAntiEdge::advanceFlippedX(SkFixed left) {
+    validate();
+    SkFixed x = SkFixedAddPin(fX0, -fDX);
+    SkFixed wy = SkIntToFixed(SkFixedFloor(fWalkY - 1));
+    pointOnLine(x, wy);
+    SkFixed partial = fPartialY ? fPartialY : SK_Fixed1;
+    SkFixed topPartial = fFirstY - wy;
+    if (topPartial > 0) {
+        partial -= topPartial;
+    }
+    if (x > fFirstX) {
+        x = fFirstX;
+        wy = fFirstY;
+    }
+    uint16_t coverage;
+    if (left >= x) {
+        fFinished = true;
+        coverage = partial - 1; // walker is to the right of edge
+    } else {
+        SkFixed y = SkFixedAddPin(fY0, -fDY);
+        SkFixed wx = SkIntToFixed(SkFixedFloor(fWalkX + SK_Fixed1));
+        pointOnLine(wx, y);
+        if (y < fFirstY) {
+            y = fFirstY;
+            wx = fFirstX;
+        }
+        bool bottomCorner = fWalkX <= fX;
+        bool topCorner = x <= wx;
+        bool halfPlane = !(topCorner ^ bottomCorner);
+        if (halfPlane) {
+            if (x - SkIntToFixed(SkFixedFloor(fX)) <= SK_Fixed1) {
+                coverage = ~((fX + x) >> 1); // avg of fx, fx+dx
+                fFinished = true;
+            } else {
+                SkAssertResult(y - SkIntToFixed(SkFixedFloor(fY)) <= SK_Fixed1);
+                coverage = ~((fY + y) >> 1);
+                fFinished = y == fY;
+                fWalkX = wx;
+                fY = fY0 = y;
+            }
+            coverage = coverage * partial >> 16;
+        } else if (bottomCorner) {
+            SkFixed xDiff = wx - fX;
+            SkAssertResult(xDiff >= 0);
+            SkAssertResult(xDiff <= SK_Fixed1);
+            SkFixed yDiff = fWalkY - y;
+            SkAssertResult(yDiff >= 0);
+            SkAssertResult(yDiff <= SK_Fixed1);
+            int xCoverage = xDiff >> 1; // throw away 1 bit so multiply
+            int yCoverage = yDiff >> 1; //  stays in range
+            int triangle = xCoverage * yCoverage; // 30 bits
+            SkFixed bottomPartial = fFirstY - y;
+            fFinished = bottomPartial >= 0;
+            if (fFinished) {
+                yCoverage = bottomPartial >> 1;
+                xCoverage = (wx - fFirstX) >> 1;
+                triangle -= xCoverage * yCoverage;
+            }
+            coverage = triangle >> 15;
+            fWalkX = wx;
+            fY = fY0 = y;
+        } else {
+            SkAssertResult(topCorner);
+            SkFixed xDiff = x - fWalkX;
+            SkAssertResult(xDiff >= 0);
+            SkAssertResult(xDiff <= SK_Fixed1);
+            SkFixed yDiff = fY - wy;
+            SkAssertResult(yDiff >= 0);
+            SkAssertResult(yDiff <= SK_Fixed1);
+            int xCoverage = xDiff >> 1; // throw away 1 bit so multiply 
+            int yCoverage = yDiff >> 1; //  stays in range
+            int triangle = xCoverage * yCoverage >> 15;
+            coverage = partial - 1 - triangle;
+            fFinished = true;
+        }
+    }
+    validate();
+    return coverage;
+}
+
+void SkAntiEdge::advanceY(SkFixed top) {
+    validate();
+    fX0 = SkFixedAddPin(fX0, fDX);
+    fPartialY = 0;
+    if (fDXFlipped) {
+        if (fX0 < fLastX) {
+            fWalkX = fX = fLastX;
+        } else {
+            fWalkX = fX = fX0;
+        }
+        SkFixed bottom = top + SK_Fixed1;
+        if (bottom > fLastY) {
+            bottom = fLastY;
+        }
+        SkFixed vert = bottom - fFirstY; // distance from y start to x-axis
+        SkFixed backupX = fFirstX + SkFixedMul(vert, fDX); // x cell to back up to
+        SkFixed distX = fFirstX - SkIntToFixed(SkFixedFloor(backupX)); // to y-axis
+        fY0 = fFirstY + SkFixedMul(fDY, distX);
+
+        fY = top + SK_Fixed1;
+        if (fY > fLastY) {
+            fY = fLastY;
+        }
+        if (fLastY < top + SK_Fixed1) {
+            fPartialY = SkFixedFraction(fLastY);
+        }
+    } else {
+        if (fX0 > fLastX) {
+            fX0 = fLastX;
+        }
+        fX = fX0;
+    }
+    fWalkY = SkIntToFixed(SkFixedFloor(fWalkY + SK_Fixed1));
+    if (fWalkY > fLastY) {
+        fWalkY = fLastY;
+    }
+    validate();
+    fFinished = false;
+}
+
+int SkAntiEdgeBuilder::build(const SkPoint pts[], int count) {
+    SkAntiEdge* edge = fEdges.append();
+    for (int index = 0; index < count; ++index) {
+        if (edge->setLine(pts[index], pts[(index + 1) % count])) {
+            edge = fEdges.append();
+        }
+    }
+    int result = fEdges.count();
+    fEdges.setCount(--result);
+    if (result > 0) {
+        sk_bzero(&fHeadEdge, sizeof(fHeadEdge));
+        sk_bzero(&fTailEdge, sizeof(fTailEdge));
+        for (int index = 0; index < result; ++index) {
+            *fList.append() = &fEdges[index];
+        }
+    }
+    return result;
+}
+
+void SkAntiEdgeBuilder::calc() {
+    for (SkAntiEdge* active = fEdges.begin(); active != fEdges.end(); ++active) {
+        active->calcLine();
+    }
+    // compute winding sum for edges
+    SkAntiEdge* first = fHeadEdge.fNext;
+    SkAntiEdge* active;
+    SkAntiEdge* listTop = first;
+    for (active = first; active != &fTailEdge; active = active->fNext) {
+        active->fWindingSum = active->fWinding;
+        while (listTop->fLastY < active->fFirstY) {
+            listTop = listTop->fNext;
+        }
+        for (SkAntiEdge* check = listTop; check->fFirstY <= active->fFirstY; check = check->fNext) {
+            if (check == active) {
+                continue;
+            }
+            if (check->fLastY <= active->fFirstY) {
+                continue;
+            }
+            if (check->fFirstX > active->fFirstX) {
+                continue;
+            }
+            if (check->fFirstX == active->fFirstX && check->fDX > active->fDX) {
+                continue;
+            }
+            active->fWindingSum += check->fWinding;
+        }
+    }
+}
+
+extern "C" {
+    static int edge_compare(const void* a, const void* b) {
+        const SkAntiEdge* edgea = *(const SkAntiEdge**)a;
+        const SkAntiEdge* edgeb = *(const SkAntiEdge**)b;
+
+        int valuea = edgea->fFirstY;
+        int valueb = edgeb->fFirstY;
+
+        if (valuea == valueb) {
+            valuea = edgea->fFirstX;
+            valueb = edgeb->fFirstX;
+        }
+        
+        if (valuea == valueb) {
+            valuea = edgea->fDX;
+            valueb = edgeb->fDX;
+        }
+
+        return valuea - valueb;
+    }
+}
+
+void SkAntiEdgeBuilder::sort(SkTDArray<SkAntiEdge*>& listOfEdges) {
+    SkAntiEdge** list = listOfEdges.begin();
+    int count = listOfEdges.count();
+    qsort(list, count, sizeof(SkAntiEdge*), edge_compare);
+
+    // link the edges in sorted order
+    for (int i = 1; i < count; i++) {
+        list[i - 1]->fNext = list[i];
+        list[i]->fPrev = list[i - 1];
+    }
+}
+
+#define kEDGE_HEAD_XY    SK_MinS32
+#define kEDGE_TAIL_XY    SK_MaxS32
+
+void SkAntiEdgeBuilder::sort() {
+    sort(fList);
+    SkAntiEdge* last = fList.end()[-1];
+    fHeadEdge.fNext = fList[0];
+    fHeadEdge.fFirstX = fHeadEdge.fFirstY = fHeadEdge.fWalkY = fHeadEdge.fLastY = kEDGE_HEAD_XY;
+    fList[0]->fPrev = &fHeadEdge;
+
+    fTailEdge.fPrev = last;
+    fTailEdge.fFirstX = fTailEdge.fFirstY = fTailEdge.fWalkY = fTailEdge.fLastY = kEDGE_TAIL_XY;
+    last->fNext = &fTailEdge;
+}
+
+static inline void remove_edge(SkAntiEdge* edge) {
+    edge->fPrev->fNext = edge->fNext;
+    edge->fNext->fPrev = edge->fPrev;
+}
+
+static inline void swap_edges(SkAntiEdge* prev, SkAntiEdge* next) {
+    SkASSERT(prev->fNext == next && next->fPrev == prev);
+
+    // remove prev from the list
+    prev->fPrev->fNext = next;
+    next->fPrev = prev->fPrev;
+
+    // insert prev after next
+    prev->fNext = next->fNext;
+    next->fNext->fPrev = prev;
+    next->fNext = prev;
+    prev->fPrev = next;
+}
+
+static void backward_insert_edge_based_on_x(SkAntiEdge* edge SkDECLAREPARAM(int, y)) {
+    SkFixed x = edge->fFirstX;
+
+    for (;;) {
+        SkAntiEdge* prev = edge->fPrev;
+
+        // add 1 to curr_y since we may have added new edges (built from curves)
+        // that start on the next scanline
+        SkASSERT(prev && SkFixedFloor(prev->fWalkY - prev->fDXFlipped) <= y + 1);
+
+        if (prev->fFirstX <= x) {
+            break;
+        }
+        swap_edges(prev, edge);
+    }
+}
+
+static void insert_new_edges(SkAntiEdge* newEdge, SkFixed curr_y) {
+    int y = SkFixedFloor(curr_y);
+    if (SkFixedFloor(newEdge->fWalkY - newEdge->fDXFlipped) < y) {
+        return;
+    }
+    while (SkFixedFloor(newEdge->fWalkY - newEdge->fDXFlipped) == y) {
+        SkAntiEdge* next = newEdge->fNext;
+        backward_insert_edge_based_on_x(newEdge  SkPARAM(y));
+        newEdge = next;
+    }
+}
+
+static int find_active_edges(int y, SkAntiEdge** activeLeft,
+                             SkAntiEdge** activeLast) {
+    SkAntiEdge* first = *activeLeft;
+    SkFixed bottom = first->fLastY;
+    SkAntiEdge* active = first->fNext;
+    first->fLinkSet = false;
+    SkFixed yLimit = SkIntToFixed(y + 1); // limiting pixel edge
+    for ( ; active->fWalkY != kEDGE_TAIL_XY; active = active->fNext) {
+        active->fLinkSet = false;
+        if (yLimit <= active->fWalkY - active->fDXFlipped) {
+            break;
+        }
+        if ((*activeLeft)->fWalkX > active->fWalkX) {
+            *activeLeft = active;
+        }
+        if (bottom > active->fLastY) {
+            bottom = active->fLastY;
+        }
+    }
+    *activeLast = active;
+    return SkFixedCeil(bottom);
+}
+
+// All edges are oriented to increase in y. Link edges with common tops and
+// bottoms so the links can share their winding sum.
+void SkAntiEdgeBuilder::link() {
+    SkAntiEdge* tail = fEdges.end();
+    // look for links forwards and backwards
+    SkAntiEdge* prev = fEdges.begin();
+    SkAntiEdge* active;
+    for (active = prev + 1; active != tail; ++active) {
+        if (prev->fWinding == active->fWinding) {
+            if (prev->fLastX == active->fFirstX && prev->fLastY == active->fFirstY) {
+                prev->fLink = active;
+                active->fLinkSet = true;
+            } else if (active->fLastX == prev->fFirstX && active->fLastY == prev->fFirstY) {
+                active->fLink = prev;
+                prev->fLinkSet = true;
+            }
+        }
+        prev = active;
+    }
+    // look for stragglers
+    prev = fEdges.begin() - 1;
+    do {
+        do {
+            if (++prev == tail) {
+                return;
+            }
+        } while (prev->fLinkSet || NULL != prev->fLink);
+        for (active = prev + 1; active != tail; ++active) {
+            if (active->fLinkSet || NULL != active->fLink) {
+                continue;
+            }
+            if (prev->fWinding != active->fWinding) {
+                continue;
+            }
+            if (prev->fLastX == active->fFirstX && prev->fLastY == active->fFirstY) {
+                prev->fLink = active;
+                active->fLinkSet = true;
+                break;
+            }
+            if (active->fLastX == prev->fFirstX && active->fLastY == prev->fFirstY) {
+                active->fLink = prev;
+                prev->fLinkSet = true;
+                break;
+            }
+        }
+    } while (true);
+}
+
+void SkAntiEdgeBuilder::split(SkAntiEdge* edge, SkFixed y) {
+    SkPoint upperPoint = {edge->fFirstX, edge->fFirstY};
+    SkPoint midPoint = {edge->fFirstX + SkMulDiv(y - edge->fFirstY,
+            edge->fLastX - edge->fFirstX, edge->fLastY - edge->fFirstY), y};
+    SkPoint lowerPoint = {edge->fLastX, edge->fLastY};
+    int8_t winding = edge->fWinding;
+    edge->setLine(upperPoint, midPoint);
+    edge->fWinding = winding;
+    SkAntiEdge* lower = fEdges.append();
+    lower->setLine(midPoint, lowerPoint);
+    lower->fWinding = winding;
+    insert_new_edges(lower, y);
+}
+
+// An edge computes pixel coverage by considering the integral winding value
+// to its left. If an edge is enclosed by fractional winding, split it.
+// FIXME: This is also a good time to find crossing edges and split them, too.
+void SkAntiEdgeBuilder::split() {
+    // create a new set of edges that describe the whole link
+    SkTDArray<SkAntiEdge> links;
+    SkAntiEdge* first = fHeadEdge.fNext;
+    SkAntiEdge* active;
+    for (active = first; active != &fTailEdge; active = active->fNext) {
+        if (active->fLinkSet || NULL == active->fLink) {
+            continue;
+        }
+        SkAntiEdge* link = links.append();
+        link->fFirstX = active->fFirstX;
+        link->fFirstY = active->fFirstY;
+        SkAntiEdge* linkEnd;
+        SkAntiEdge* next = active;
+        do {
+            linkEnd = next;
+            next = next->fLink;
+        } while (NULL != next);
+        link->fLastX = linkEnd->fLastX;
+        link->fLastY = linkEnd->fLastY;
+    }
+    // create a list of all edges, links and singletons
+    SkTDArray<SkAntiEdge*> list;
+    for (active = links.begin(); active != links.end(); ++active) {
+        *list.append() = active;
+    }
+    for (active = first; active != &fTailEdge; active = active->fNext) {
+        if (!active->fLinkSet && NULL == active->fLink) {
+            SkAntiEdge* link = links.append();
+            link->fFirstX = active->fFirstX;
+            link->fFirstY = active->fFirstY;
+            link->fLastX = active->fLastX;
+            link->fLastY = active->fLastY;
+            *list.append() = link;
+        }
+    }
+    SkAntiEdge tail;
+    tail.fFirstY = tail.fLastY = kEDGE_TAIL_XY;
+    *list.append() = &tail;
+    sort(list);
+    // walk the list, splitting edges partially occluded on the left
+    SkAntiEdge* listTop = list[0];
+    for (active = first; active != &fTailEdge; active = active->fNext) {
+        while (listTop->fLastY < active->fFirstY) {
+            listTop = listTop->fNext;
+        }
+        for (SkAntiEdge* check = listTop; check->fFirstY < active->fLastY; check = check->fNext) {
+            if (check->fFirstX > active->fFirstX) {
+                continue;
+            }
+            if (check->fFirstX == active->fFirstX && check->fDX > active->fDX) {
+                continue;
+            }
+            if (check->fFirstY > active->fFirstY) {
+                split(active, check->fFirstY);
+            }
+            if (check->fLastY < active->fLastY) {
+                split(active, check->fLastY);
+            }
+        }
+    }
+}
+
+static inline uint8_t coverage_to_8(int coverage) {
+    uint16_t x = coverage < 0 ? 0 : coverage > 0xFFFF ? 0xFFFF : coverage;
+    // for values 0x7FFF and smaller, add (0x7F - high byte) and trunc
+    // for values 0x8000 and larger, subtract (high byte - 0x80) and trunc
+    return (x + 0x7f + (x >> 15) - (x >> 8)) >> 8;
+}
+
+void SkAntiEdgeBuilder::walk(uint8_t* result, int rowBytes, int height) {
+    SkAntiEdge* first = fHeadEdge.fNext;
+    SkFixed top = first->fWalkY - first->fDXFlipped;
+    int y = SkFixedFloor(top);
+    do {
+        SkAntiEdge* activeLeft = first;
+        SkAntiEdge* activeLast, * active;
+        int yLast = find_active_edges(y, &activeLeft, &activeLast);
+        while (y < yLast) {
+            SkAssertResult(y >= 0);
+            SkAssertResult(y < height);
+            SkFixed left = activeLeft->fWalkX;
+            int x = SkFixedFloor(left);
+            uint8_t* resultPtr = &result[y * rowBytes + x];
+            bool finished;
+            do {
+                left = SkIntToFixed(x);
+                SkAssertResult(x >= 0);
+              //  SkAssertResult(x < pixelCol);
+                if (x >= rowBytes) { // FIXME: cumulative error in fX += fDX
+                    break;           // fails to set fFinished early enough
+                }                    // see test 6 (dy<dx)
+                finished = true;
+                int coverage = 0;
+                for (active = first; active != activeLast; active = active->fNext) {
+                    if (left + SK_Fixed1 <= active->fX) {
+                        finished = false;
+                        continue; // walker is to the left of edge
+                    }
+                    int cover = active->fDXFlipped ?
+                        active->advanceFlippedX(left) : active->advanceX(left);
+                    if (0 == active->fWindingSum) {
+                        cover = -cover;
+                    }
+                    coverage += cover;
+                    finished &= active->fFinished;
+                }
+                uint8_t old = *resultPtr;
+                uint8_t pix = coverage_to_8(coverage);
+                uint8_t blend = old > pix ? old : pix; 
+                *resultPtr++ = blend;
+                ++x;
+            } while (!finished);
+            ++y;
+            top = SkIntToFixed(y);
+            SkFixed topLimit = top + SK_Fixed1;
+            SkFixed xSort = -SK_FixedMax;
+            for (active = first; active != activeLast; active = active->fNext) {
+                if (xSort > active->fX || topLimit > active->fLastY) {
+                    yLast = y; // recompute bottom after all Ys are advanced
+                }
+                xSort = active->fX;
+                if (active->fWalkY < active->fLastY) {
+                    active->advanceY(top);
+                }
+            }
+            for (active = first; active != activeLast; ) {
+                SkAntiEdge* next = active->fNext;
+                if (top >= active->fLastY) {
+                    remove_edge(active);
+                }
+                active = next;
+            }
+            first = fHeadEdge.fNext;
+        }
+        SkAntiEdge* prev = activeLast->fPrev;
+        if (prev != &fHeadEdge) {
+            insert_new_edges(prev, top);
+            first = fHeadEdge.fNext;
+        }
+    } while (first->fWalkY < kEDGE_TAIL_XY);
+}
+
+void SkAntiEdgeBuilder::process(const SkPoint* points, int ptCount,
+        uint8_t* result, int pixelCol, int pixelRow) {
+    if (ptCount < 3) {
+        return;
+    }
+    int count = build(points, ptCount);
+    if (count == 0) {
+        return;
+    }
+    SkAssertResult(count > 1);
+    link();
+    sort();
+    split();
+    calc();
+    walk(result, pixelCol, pixelRow);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+int test3by3_test;
+
+// input is a rectangle
+static void test_3_by_3() {
+    const int pixelRow = 3;
+    const int pixelCol = 3;
+    const int ptCount = 4;
+    const int pixelCount = pixelRow * pixelCol;
+    const SkPoint tests[][ptCount] = {
+        {{2.0f, 1.0f}, {1.0f, 1.0f}, {1.0f, 2.0f}, {2.0f, 2.0f}}, // 0: full rect
+        {{2.5f, 1.0f}, {1.5f, 1.0f}, {1.5f, 2.0f}, {2.5f, 2.0f}}, // 1: y edge
+        {{2.0f, 1.5f}, {1.0f, 1.5f}, {1.0f, 2.5f}, {2.0f, 2.5f}}, // 2: x edge
+        {{2.5f, 1.5f}, {1.5f, 1.5f}, {1.5f, 2.5f}, {2.5f, 2.5f}}, // 3: x/y edge
+        {{2.8f, 0.2f}, {0.2f, 0.2f}, {0.2f, 2.8f}, {2.8f, 2.8f}}, // 4: large
+        {{1.8f, 1.2f}, {1.2f, 1.2f}, {1.2f, 1.8f}, {1.8f, 1.8f}}, // 5: small
+        {{0.0f, 0.0f}, {0.0f, 1.0f}, {3.0f, 2.0f}, {3.0f, 1.0f}}, // 6: dy<dx
+        {{3.0f, 0.0f}, {0.0f, 1.0f}, {0.0f, 2.0f}, {3.0f, 1.0f}}, // 7: dy<-dx
+        {{1.0f, 0.0f}, {0.0f, 0.0f}, {1.0f, 3.0f}, {2.0f, 3.0f}}, // 8: dy>dx
+        {{2.0f, 0.0f}, {1.0f, 0.0f}, {0.0f, 3.0f}, {1.0f, 3.0f}}, // 9: dy>-dx
+        {{0.5f, 0.5f}, {0.5f, 1.5f}, {2.5f, 2.5f}, {2.5f, 1.5f}}, // 10: dy<dx 2
+        {{2.5f, 0.5f}, {0.5f, 1.5f}, {0.5f, 2.5f}, {2.5f, 1.5f}}, // 11: dy<-dx 2
+        {{0.0f, 0.0f}, {2.0f, 0.0f}, {2.0f, 2.0f}, {0.0f, 2.0f}}, // 12: 2x2
+        {{0.0f, 0.0f}, {3.0f, 0.0f}, {3.0f, 3.0f}, {0.0f, 3.0f}}, // 13: 3x3
+        {{1.75f, 0.25f}, {2.75f, 1.25f}, {1.25f, 2.75f}, {0.25f, 1.75f}}, // 14
+        {{2.25f, 0.25f}, {2.75f, 0.75f}, {0.75f, 2.75f}, {0.25f, 2.25f}}, // 15
+        {{0.25f, 0.75f}, {0.75f, 0.25f}, {2.75f, 2.25f}, {2.25f, 2.75f}}, // 16
+        {{1.25f, 0.50f}, {1.75f, 0.25f}, {2.75f, 2.25f}, {2.25f, 2.50f}}, // 17
+        {{1.00f, 0.75f}, {2.00f, 0.50f}, {2.00f, 1.50f}, {1.00f, 1.75f}}, // 18
+        {{1.00f, 0.50f}, {2.00f, 0.75f}, {2.00f, 1.75f}, {1.00f, 1.50f}}, // 19
+        {{1.00f, 0.75f}, {1.00f, 1.75f}, {2.00f, 1.50f}, {2.00f, 0.50f}}, // 20
+        {{1.00f, 0.50f}, {1.00f, 1.50f}, {2.00f, 1.75f}, {2.00f, 0.75f}}, // 21
+    };
+    const uint8_t results[][pixelCount] = {
+        {0x00, 0x00, 0x00, // 0: 1 pixel rect
+         0x00, 0xFF, 0x00,
+         0x00, 0x00, 0x00},
+        {0x00, 0x00, 0x00, // 1: y edge
+         0x00, 0x7F, 0x80,
+         0x00, 0x00, 0x00},
+        {0x00, 0x00, 0x00, // 2: x edge
+         0x00, 0x7F, 0x00,
+         0x00, 0x7F, 0x00},
+        {0x00, 0x00, 0x00, // 3: x/y edge
+         0x00, 0x40, 0x40,
+         0x00, 0x40, 0x40},
+        {0xA3, 0xCC, 0xA3, // 4: large
+         0xCC, 0xFF, 0xCC,
+         0xA3, 0xCC, 0xA3},
+        {0x00, 0x00, 0x00, // 5: small
+         0x00, 0x5C, 0x00,
+         0x00, 0x00, 0x00},
+        {0xD5, 0x80, 0x2B, // 6: dy<dx
+         0x2A, 0x7F, 0xD4,
+         0x00, 0x00, 0x00},
+        {0x2B, 0x80, 0xD5, // 7: dy<-dx
+         0xD4, 0x7F, 0x2A,
+         0x00, 0x00, 0x00},
+        {0xD5, 0x2A, 0x00, // 8: dy>dx
+         0x80, 0x7F, 0x00,
+         0x2B, 0xD4, 0x00},
+        {0x2A, 0xD5, 0x00, // 9: dy>-dx
+         0x7F, 0x80, 0x00,
+         0xD4, 0x2B, 0x00},
+        {0x30, 0x10, 0x00, // 10: dy<dx 2
+         0x50, 0xDF, 0x50,
+         0x00, 0x10, 0x30},
+        {0x00, 0x10, 0x30, // 11: dy<-dx 2
+         0x50, 0xDF, 0x50,
+         0x30, 0x10, 0x00},
+        {0xFF, 0xFF, 0x00, // 12: 2x2
+         0xFF, 0xFF, 0x00,
+         0x00, 0x00, 0x00},
+        {0xFF, 0xFF, 0xFF, // 13: 3x3
+         0xFF, 0xFF, 0xFF,
+         0xFF, 0xFF, 0xFF},
+        {0x00, 0x70, 0x20, // 14
+         0x70, 0xFF, 0x70,
+         0x20, 0x70, 0x00},
+        {0x00, 0x20, 0x60, // 15
+         0x20, 0xBF, 0x20,
+         0x60, 0x20, 0x00},
+        {0x60, 0x20, 0x00, // 16
+         0x20, 0xBF, 0x20,
+         0x00, 0x20, 0x60},
+        {0x00, 0x60, 0x04, // 17
+         0x00, 0x40, 0x60,
+         0x00, 0x00, 0x3C},
+        {0x00, 0x60, 0x00, // 18
+         0x00, 0x9F, 0x00,
+         0x00, 0x00, 0x00},
+        {0x00, 0x60, 0x00, // 19
+         0x00, 0x9F, 0x00,
+         0x00, 0x00, 0x00},
+        {0x00, 0x60, 0x00, // 20
+         0x00, 0x9F, 0x00,
+         0x00, 0x00, 0x00},
+        {0x00, 0x60, 0x00, // 21
+         0x00, 0x9F, 0x00,
+         0x00, 0x00, 0x00},
+    };
+    const int testCount = sizeof(tests) / sizeof(tests[0]);
+    SkAssertResult(testCount == sizeof(results) / sizeof(results[0]));
+    int testFirst = test3by3_test < 0 ? 0 : test3by3_test;
+    int testLast = test3by3_test < 0 ? testCount : test3by3_test + 1;
+    for (int testIndex = testFirst; testIndex < testLast; ++testIndex) {
+        uint8_t result[pixelRow][pixelCol];
+        sk_bzero(result, sizeof(result));
+        const SkPoint* rect = tests[testIndex];
+        SkAntiEdgeBuilder builder;
+        builder.process(rect, ptCount, result[0], pixelCol, pixelRow);
+        SkAssertResult(memcmp(results[testIndex], result[0], pixelCount) == 0);
+    }
+}
+
+// input has arbitrary number of points
+static void test_arbitrary_3_by_3() {
+    const int pixelRow = 3;
+    const int pixelCol = 3;
+    const int pixelCount = pixelRow * pixelCol;
+    const SkPoint t1[] = { {1,1}, {2,1}, {2,1.5f}, {1,1.5f}, {1,2}, {2,2},
+        {2,1.5f}, {1,1.5f}, {1,1} };
+    const SkPoint* tests[] = { t1 };
+    size_t testPts[] = { sizeof(t1) / sizeof(t1[0]) };
+    const uint8_t results[][pixelCount] = {
+        {0x00, 0x00, 0x00, // 0: 1 pixel rect
+         0x00, 0xFF, 0x00,
+         0x00, 0x00, 0x00},
+    };
+    const int testCount = sizeof(tests) / sizeof(tests[0]);
+    SkAssertResult(testCount == sizeof(results) / sizeof(results[0]));
+    int testFirst = test3by3_test < 0 ? 0 : test3by3_test;
+    int testLast = test3by3_test < 0 ? testCount : test3by3_test + 1;
+    for (int testIndex = testFirst; testIndex < testLast; ++testIndex) {
+        uint8_t result[pixelRow][pixelCol];
+        sk_bzero(result, sizeof(result));
+        const SkPoint* pts = tests[testIndex];
+        size_t ptCount = testPts[testIndex];
+        SkAntiEdgeBuilder builder;
+        builder.process(pts, ptCount, result[0], pixelCol, pixelRow);
+        SkAssertResult(memcmp(results[testIndex], result[0], pixelCount) == 0);
+    }
+}
+
+#include "SkRect.h"
+#include "SkPath.h"
+
+int testsweep_test;
+
+static void create_sweep(uint8_t* result, int pixelRow, int pixelCol, SkScalar rectWidth) {
+    const int ptCount = 4;
+    SkRect refRect = {pixelCol / 2 - rectWidth / 2, 5,
+                      pixelCol / 2 + rectWidth / 2, pixelRow / 2 - 5};
+    SkPath refPath;
+    refPath.addRect(refRect);
+    SkScalar angleFirst = testsweep_test < 0 ? 0 : testsweep_test;
+    SkScalar angleLast = testsweep_test < 0 ? 360 : testsweep_test + 1;
+    for (SkScalar angle = angleFirst; angle < angleLast; angle += 12) {
+        SkPath rotPath;
+        SkMatrix matrix;
+        matrix.setRotate(angle, SkIntToScalar(pixelCol) / 2,
+            SkIntToScalar(pixelRow) / 2);
+        refPath.transform(matrix, &rotPath);
+        SkPoint rect[ptCount], temp[2];
+        SkPath::Iter iter(rotPath, false);
+        int index = 0;
+        for (;;) {
+            SkPath::Verb verb = iter.next(temp);
+            if (verb == SkPath::kMove_Verb) {
+                continue;
+            }
+            if (verb == SkPath::kClose_Verb) {
+                break;
+            }
+            SkAssertResult(SkPath::kLine_Verb == verb);
+            rect[index++] = temp[0];
+        }
+        SkAntiEdgeBuilder builder;
+        builder.process(rect, ptCount, result, pixelCol, pixelRow);
+    }
+}
+
+static void create_horz(uint8_t* result, int pixelRow, int pixelCol) {
+    const int ptCount = 4;
+    for (SkScalar x = 0; x < 100; x += 5) {
+        SkPoint rect[ptCount];
+        rect[0].fX = 0;     rect[0].fY = x;
+        rect[1].fX = 100;   rect[1].fY = x;
+        rect[2].fX = 100;   rect[2].fY = x + x / 50;
+        rect[3].fX = 0;     rect[3].fY = x + x / 50;
+        SkAntiEdgeBuilder builder;
+        builder.process(rect, ptCount, result, pixelCol, pixelRow);
+    }
+}
+
+static void create_vert(uint8_t* result, int pixelRow, int pixelCol) {
+    const int ptCount = 4;
+    for (SkScalar x = 0; x < 100; x += 5) {
+        SkPoint rect[ptCount];
+        rect[0].fY = 0;     rect[0].fX = x;
+        rect[1].fY = 100;   rect[1].fX = x;
+        rect[2].fY = 100;   rect[2].fX = x + x / 50;
+        rect[3].fY = 0;     rect[3].fX = x + x / 50;
+        SkAntiEdgeBuilder builder;
+        builder.process(rect, ptCount, result, pixelCol, pixelRow);
+    }
+}
+
+static void create_angle(uint8_t* result, int pixelRow, int pixelCol, SkScalar angle) {
+    const int ptCount = 4;
+    SkRect refRect = {25, 25, 125, 125};
+    SkPath refPath;
+    for (SkScalar x = 30; x < 125; x += 5) {
+        refRect.fTop = x;
+        refRect.fBottom = x + (x - 25) / 50;
+        refPath.addRect(refRect);
+    }
+    SkPath rotPath;
+    SkMatrix matrix;
+    matrix.setRotate(angle, 75, 75);
+    refPath.transform(matrix, &rotPath);
+    SkPath::Iter iter(rotPath, false);
+    for (SkScalar x = 30; x < 125; x += 5) {
+        SkPoint rect[ptCount], temp[2];
+        int index = 0;
+        for (;;) {
+            SkPath::Verb verb = iter.next(temp);
+            if (verb == SkPath::kMove_Verb) {
+                continue;
+            }
+            if (verb == SkPath::kClose_Verb) {
+                break;
+            }
+            SkAssertResult(SkPath::kLine_Verb == verb);
+            rect[index++] = temp[0];
+        }
+    //    if ((x == 30 || x == 75) && angle == 12) continue;
+        SkAntiEdgeBuilder builder;
+        builder.process(rect, ptCount, result, pixelCol, pixelRow);
+    }
+}
+
+static void test_sweep() {
+    const int pixelRow = 100;
+    const int pixelCol = 100;
+    uint8_t result[pixelRow][pixelCol];
+    sk_bzero(result, sizeof(result));
+    create_sweep(result[0], pixelRow, pixelCol, 1);
+}
+
+static void test_horz() {
+    const int pixelRow = 100;
+    const int pixelCol = 100;
+    uint8_t result[pixelRow][pixelCol];
+    sk_bzero(result, sizeof(result));
+    create_horz(result[0], pixelRow, pixelCol);
+}
+
+static void test_vert() {
+    const int pixelRow = 100;
+    const int pixelCol = 100;
+    uint8_t result[pixelRow][pixelCol];
+    sk_bzero(result, sizeof(result));
+    create_vert(result[0], pixelRow, pixelCol);
+}
+
+static void test_angle(SkScalar angle) {
+    const int pixelRow = 150;
+    const int pixelCol = 150;
+    uint8_t result[pixelRow][pixelCol];
+    sk_bzero(result, sizeof(result));
+    create_angle(result[0], pixelRow, pixelCol, angle);
+}
+
+#include "SkBitmap.h"
+
+void CreateSweep(SkBitmap* sweep, SkScalar rectWidth) {
+    const int pixelRow = 100;
+    const int pixelCol = 100;
+    sweep->setConfig(SkBitmap::kA8_Config, pixelCol, pixelRow);
+    sweep->allocPixels();
+    sweep->eraseColor(0);
+    sweep->lockPixels();
+    void* pixels = sweep->getPixels();
+    create_sweep((uint8_t*) pixels, pixelRow, pixelCol, rectWidth);
+    sweep->unlockPixels();
+}
+
+void CreateHorz(SkBitmap* sweep) {
+    const int pixelRow = 100;
+    const int pixelCol = 100;
+    sweep->setConfig(SkBitmap::kA8_Config, pixelCol, pixelRow);
+    sweep->allocPixels();
+    sweep->eraseColor(0);
+    sweep->lockPixels();
+    void* pixels = sweep->getPixels();
+    create_horz((uint8_t*) pixels, pixelRow, pixelCol);
+    sweep->unlockPixels();
+}
+
+void CreateVert(SkBitmap* sweep) {
+    const int pixelRow = 100;
+    const int pixelCol = 100;
+    sweep->setConfig(SkBitmap::kA8_Config, pixelCol, pixelRow);
+    sweep->allocPixels();
+    sweep->eraseColor(0);
+    sweep->lockPixels();
+    void* pixels = sweep->getPixels();
+    create_vert((uint8_t*) pixels, pixelRow, pixelCol);
+    sweep->unlockPixels();
+}
+
+void CreateAngle(SkBitmap* sweep, SkScalar angle) {
+    const int pixelRow = 150;
+    const int pixelCol = 150;
+    sweep->setConfig(SkBitmap::kA8_Config, pixelCol, pixelRow);
+    sweep->allocPixels();
+    sweep->eraseColor(0);
+    sweep->lockPixels();
+    void* pixels = sweep->getPixels();
+    create_angle((uint8_t*) pixels, pixelRow, pixelCol, angle);
+    sweep->unlockPixels();
+}
+
+#include "SkCanvas.h"
+
+static void testPng() {
+    SkCanvas canvas;
+    SkBitmap device;
+    device.setConfig(SkBitmap::kARGB_8888_Config, 4, 4);
+    device.allocPixels();
+    device.eraseColor(0xFFFFFFFF);
+    canvas.setBitmapDevice(device);
+    canvas.drawARGB(167, 0, 0, 0);
+    device.lockPixels();
+    unsigned char* pixels = (unsigned char*) device.getPixels();
+    SkDebugf("%02x%02x%02x%02x", pixels[3], pixels[2], pixels[1], pixels[0]);
+}
+
+void SkAntiEdge_Test() {
+    testPng();
+    test_arbitrary_3_by_3();
+    test_angle(12);
+#if 0
+    test3by3_test = 18;
+#else
+    test3by3_test = -1;
+#endif
+#if 0
+    testsweep_test = 7 * 12;
+#else
+    testsweep_test = -1;
+#endif
+    if (testsweep_test == -1) {
+        test_3_by_3();
+    }
+    test_sweep();
+    test_horz();
+    test_vert();
+}
+
diff --git a/experimental/Intersection/SkAntiEdge.h b/experimental/Intersection/SkAntiEdge.h
new file mode 100644
index 0000000..a85c518
--- /dev/null
+++ b/experimental/Intersection/SkAntiEdge.h
@@ -0,0 +1,79 @@
+/*
+ *  SkAntiEdge.h
+ *  core
+ *
+ *  Created by Cary Clark on 5/6/11.
+ *  Copyright 2011 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#ifndef SkAntiEdge_DEFINED
+#define SkAntiEdge_DEFINED
+
+#include "SkFixed.h"
+#include "SkTDArray.h"
+
+struct SkBitmap;
+struct SkPoint;
+
+struct SkAntiEdge {
+    SkAntiEdge* fNext; // list in walking order (y, then x, then diag)
+    SkAntiEdge* fPrev; // reverse in walking order
+    SkAntiEdge* fLink; // list in connected order, top to bottom
+
+    SkFixed     fFirstX; // starting X
+    SkFixed     fFirstY; // starting Y
+    SkFixed     fLastX; // ending X
+    SkFixed     fLastY; // ending Y
+    SkFixed     fX0; // computed edge current value (may be off end)
+    SkFixed     fY0;
+    SkFixed     fX; // edge current value (always on edge)
+    SkFixed     fY;
+    SkFixed     fDX; // change in X per unit step in Y
+    SkFixed     fDY; // change in Y per unit step in X
+    SkFixed     fWalkX; // unit step position (integer after initial step)
+    SkFixed     fWalkY;
+    uint16_t    fPartialY; // initial partial coverage in Y (0 .. SkFixed1]
+    int16_t     fWindingSum; // winding including contributions to the left
+    int8_t      fWinding; // 1 or -1 (could be 2 bits)
+    bool        fFinished : 1;
+    unsigned    fDXFlipped : 1; // used as bool and to adjust calculations (0/1)
+    bool        fLinkSet : 1; // set if edge has been attached to another edge
+
+    void calcLine();
+    bool setLine(const SkPoint& p0, const SkPoint& p1);
+    uint16_t advanceX(SkFixed left);
+    uint16_t advanceFlippedX(SkFixed left);
+    void advanceY(SkFixed top);
+// FIXME: mark DEBUG
+    void pointInLine(SkFixed x, SkFixed y);
+    void pointOnLine(SkFixed x, SkFixed y);
+    void validate();
+};
+
+class SkAntiEdgeBuilder {
+public:
+void process(const SkPoint* points, int ptCount,
+        uint8_t* result, int pixelCol, int pixelRow);
+private:
+    int build(const SkPoint pts[], int count);
+    void calc();
+    void link();
+    void sort();
+    void sort(SkTDArray<SkAntiEdge*>&);
+    void split();
+    void split(SkAntiEdge* edge, SkFixed y);
+    void walk(uint8_t* result, int rowBytes, int height);
+    SkAntiEdge fHeadEdge;
+    SkAntiEdge fTailEdge;
+    SkTDArray<SkAntiEdge> fEdges;
+    SkTDArray<SkAntiEdge*> fList;
+};
+
+void SkAntiEdge_Test();
+void CreateSweep(SkBitmap* , float width);
+void CreateHorz(SkBitmap* );
+void CreateVert(SkBitmap* );
+void CreateAngle(SkBitmap* sweep, float angle);
+
+#endif
diff --git a/experimental/Intersection/TestUtilities.cpp b/experimental/Intersection/TestUtilities.cpp
new file mode 100644
index 0000000..5e3c18e
--- /dev/null
+++ b/experimental/Intersection/TestUtilities.cpp
@@ -0,0 +1,69 @@
+#include "CubicIntersection.h"
+#include "TestUtilities.h"
+
+void quad_to_cubic(const Quadratic& quad, Cubic& cubic) {
+    cubic[0] = quad[0];
+    cubic[1].x = quad[0].x / 3 + quad[1].x * 2 / 3;
+    cubic[1].y = quad[0].y / 3 + quad[1].y * 2 / 3;
+    cubic[2].x = quad[2].x / 3 + quad[1].x * 2 / 3;
+    cubic[2].y = quad[2].y / 3 + quad[1].y * 2 / 3;
+    cubic[3] = quad[2];
+}
+
+static bool tiny(const Cubic& cubic) {
+    int index, minX, maxX, minY, maxY;
+    minX = maxX = minY = maxY = 0;
+    for (index = 1; index < 4; ++index) {
+        if (cubic[minX].x > cubic[index].x) {
+            minX = index;
+        }
+        if (cubic[minY].y > cubic[index].y) {
+            minY = index;
+        }
+        if (cubic[maxX].x < cubic[index].x) {
+            maxX = index;
+        }
+        if (cubic[maxY].y < cubic[index].y) {
+            maxY = index;
+        }
+    }
+    return     approximately_equal(cubic[maxX].x, cubic[minX].x)
+            && approximately_equal(cubic[maxY].y, cubic[minY].y);
+}
+
+void find_tight_bounds(const Cubic& cubic, _Rect& bounds) {
+    CubicPair cubicPair;
+    chop_at(cubic, cubicPair, 0.5);
+    if (!tiny(cubicPair.first()) && !controls_inside(cubicPair.first())) {
+        find_tight_bounds(cubicPair.first(), bounds);
+    } else {
+        bounds.add(cubicPair.first()[0]);
+        bounds.add(cubicPair.first()[3]);
+    }
+    if (!tiny(cubicPair.second()) && !controls_inside(cubicPair.second())) {
+        find_tight_bounds(cubicPair.second(), bounds);
+    } else {
+        bounds.add(cubicPair.second()[0]);
+        bounds.add(cubicPair.second()[3]);
+    }
+}
+
+bool controls_inside(const Cubic& cubic) {
+    return
+        ((cubic[0].x <= cubic[1].x && cubic[0].x <= cubic[2].x && cubic[1].x <= cubic[3].x && cubic[2].x <= cubic[3].x)
+     ||  (cubic[0].x >= cubic[1].x && cubic[0].x >= cubic[2].x && cubic[1].x >= cubic[3].x && cubic[2].x >= cubic[3].x))
+     && ((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));
+}
+
+void xy_at_t(const Cubic& cubic, double t, double& x, double& y) {
+    double one_t = 1 - t;
+    double one_t2 = one_t * one_t;
+    double a = one_t2 * one_t;
+    double b = 3 * one_t2 * t;
+    double t2 = t * t;
+    double c = 3 * one_t * t2;
+    double d = t2 * t;
+    x = a * cubic[0].x + b * cubic[1].x + c * cubic[2].x + d * cubic[3].x;
+    y = a * cubic[0].y + b * cubic[1].y + c * cubic[2].y + d * cubic[3].y;
+}
diff --git a/experimental/Intersection/TestUtilities.h b/experimental/Intersection/TestUtilities.h
new file mode 100644
index 0000000..c67d702
--- /dev/null
+++ b/experimental/Intersection/TestUtilities.h
@@ -0,0 +1,6 @@
+#include "DataTypes.h"
+
+bool controls_inside(const Cubic& );
+void find_tight_bounds(const Cubic& , _Rect& );
+void quad_to_cubic(const Quadratic& , Cubic& );
+void xy_at_t(const Cubic& , double t, double& x, double& y);