Find cubic KLM functionals directly

- Updates GrPathUtils to computes the KLM functionals directly instead
  of deriving them from their explicit values at the control points.
- Updates the utility to return these functionals as a matrix
  rather than an array of scalar values.
- Adds a benchmark for chopCubicAtLoopIntersection.

BUG=skia:

Change-Id: I97a9b5cf610d33e15c9af96b9d9a8eb4a94b1ca7
Reviewed-on: https://skia-review.googlesource.com/9951
Commit-Queue: Chris Dalton <csmartdalton@google.com>
Reviewed-by: Greg Daniel <egdaniel@google.com>

6 files changed, 404 insertions, 251 deletions

diff --git a/bench/CubicKLMBench.cpp b/bench/CubicKLMBench.cpp
new file mode 100644
index 0000000000..3c8f740bc7
--- /dev/null
+++ b/bench/CubicKLMBench.cpp
@@ -0,0 +1,68 @@
+/*
+ * Copyright 2017 Google Inc.
+ *
+ * Use of this source code is governed by a BSD-style license that can be
+ * found in the LICENSE file.
+ */
+
+#include "Benchmark.h"
+
+#if SK_SUPPORT_GPU
+
+#include "GrPathUtils.h"
+#include "SkGeometry.h"
+
+class CubicKLMBench : public Benchmark {
+public:
+    CubicKLMBench(SkScalar x0, SkScalar y0, SkScalar x1, SkScalar y1,
+                  SkScalar x2, SkScalar y2, SkScalar x3, SkScalar y3)  {
+        fPoints[0].set(x0, y0);
+        fPoints[1].set(x1, y1);
+        fPoints[2].set(x2, y2);
+        fPoints[3].set(x3, y3);
+
+        fName = "cubic_klm_";
+        SkScalar d[3];
+        switch (SkClassifyCubic(fPoints, d)) {
+            case kSerpentine_SkCubicType:
+                fName.append("serp");
+                break;
+            case kLoop_SkCubicType:
+                fName.append("loop");
+                break;
+            default:
+                SkFAIL("Unexpected cubic type");
+                break;
+        }
+    }
+
+    bool isSuitableFor(Backend backend) override {
+        return backend == kNonRendering_Backend;
+    }
+
+    const char* onGetName() override {
+        return fName.c_str();
+    }
+
+    void onDraw(int loops, SkCanvas*) override {
+        SkPoint dst[10];
+        SkMatrix klm;
+        int loopIdx;
+        for (int i = 0; i < loops * 50000; ++i) {
+            GrPathUtils::chopCubicAtLoopIntersection(fPoints, dst, &klm, &loopIdx);
+        }
+    }
+
+private:
+    SkPoint     fPoints[4];
+    SkString    fName;
+
+    typedef Benchmark INHERITED;
+};
+
+DEF_BENCH( return new CubicKLMBench(285.625f, 499.687f, 411.625f, 808.188f,
+                                    1064.62f, 135.688f, 1042.63f, 585.187f); )
+DEF_BENCH( return new CubicKLMBench(635.625f, 614.687f, 171.625f, 236.188f,
+                                    1064.62f, 135.688f, 516.625f, 570.187f); )
+
+#endif
diff --git a/gm/beziereffects.cpp b/gm/beziereffects.cpp
index 9167410e9f..523ccd31f8 100644
--- a/gm/beziereffects.cpp
+++ b/gm/beziereffects.cpp
@@ -22,10 +22,6 @@
 
 #include "effects/GrBezierEffect.h"
 
-static inline SkScalar eval_line(const SkPoint& p, const SkScalar lineEq[3], SkScalar sign) {
-    return sign * (lineEq[0] * p.fX + lineEq[1] * p.fY + lineEq[2]);
-}
-
 namespace skiagm {
 
 class BezierCubicOrConicTestOp : public GrTestMeshDrawOp {
@@ -35,20 +31,21 @@ public:
     const char* name() const override { return "BezierCubicOrConicTestOp"; }
 
     static std::unique_ptr<GrMeshDrawOp> Make(sk_sp<GrGeometryProcessor> gp, const SkRect& rect,
-                                              GrColor color, const SkScalar klmEqs[9],
-                                              SkScalar sign) {
+                                              GrColor color, const SkMatrix& klm, SkScalar sign) {
         return std::unique_ptr<GrMeshDrawOp>(
-                new BezierCubicOrConicTestOp(gp, rect, color, klmEqs, sign));
+                new BezierCubicOrConicTestOp(gp, rect, color, klm, sign));
     }
 
 private:
     BezierCubicOrConicTestOp(sk_sp<GrGeometryProcessor> gp, const SkRect& rect, GrColor color,
-                             const SkScalar klmEqs[9], SkScalar sign)
-            : INHERITED(ClassID(), rect, color), fRect(rect), fGeometryProcessor(std::move(gp)) {
-        for (int i = 0; i < 9; i++) {
-            fKlmEqs[i] = klmEqs[i];
+                             const SkMatrix& klm, SkScalar sign)
+        : INHERITED(ClassID(), rect, color)
+        , fKLM(klm)
+        , fRect(rect)
+        , fGeometryProcessor(std::move(gp)) {
+        if (1 != sign) {
+            fKLM.postScale(sign, sign);
         }
-        fSign = sign;
     }
     struct Vertex {
         SkPoint fPosition;
@@ -66,15 +63,13 @@ private:
         verts[0].fPosition.setRectFan(fRect.fLeft, fRect.fTop, fRect.fRight, fRect.fBottom,
                                       sizeof(Vertex));
         for (int v = 0; v < 4; ++v) {
-            verts[v].fKLM[0] = eval_line(verts[v].fPosition, fKlmEqs + 0, fSign);
-            verts[v].fKLM[1] = eval_line(verts[v].fPosition, fKlmEqs + 3, fSign);
-            verts[v].fKLM[2] = eval_line(verts[v].fPosition, fKlmEqs + 6, 1.f);
+            SkScalar pt3[3] = {verts[v].fPosition.x(), verts[v].fPosition.y(), 1.f};
+            fKLM.mapHomogeneousPoints(verts[v].fKLM, pt3, 1);
         }
         helper.recordDraw(target, fGeometryProcessor.get());
     }
 
-    SkScalar fKlmEqs[9];
-    SkScalar fSign;
+    SkMatrix fKLM;
     SkRect fRect;
     sk_sp<GrGeometryProcessor> fGeometryProcessor;
 
@@ -155,11 +150,11 @@ protected:
                     {x + baseControlPts[3].fX, y + baseControlPts[3].fY}
                 };
                 SkPoint chopped[10];
-                SkScalar klmEqs[9];
+                SkMatrix klm;
                 int loopIndex;
                 int cnt = GrPathUtils::chopCubicAtLoopIntersection(controlPts,
                                                                    chopped,
-                                                                   klmEqs,
+                                                                   &klm,
                                                                    &loopIndex);
 
                 SkPaint ctrlPtPaint;
@@ -203,7 +198,7 @@ protected:
                     }
 
                     std::unique_ptr<GrMeshDrawOp> op =
-                            BezierCubicOrConicTestOp::Make(gp, bounds, color, klmEqs, sign);
+                            BezierCubicOrConicTestOp::Make(gp, bounds, color, klm, sign);
 
                     renderTargetContext->priv().testingOnly_addMeshDrawOp(
                             std::move(grPaint), GrAAType::kNone, std::move(op));
@@ -296,9 +291,9 @@ protected:
                     {x + baseControlPts[2].fX, y + baseControlPts[2].fY}
                 };
                 SkConic dst[4];
-                SkScalar klmEqs[9];
+                SkMatrix klm;
                 int cnt = chop_conic(controlPts, dst, weight);
-                GrPathUtils::getConicKLM(controlPts, weight, klmEqs);
+                GrPathUtils::getConicKLM(controlPts, weight, &klm);
 
                 SkPaint ctrlPtPaint;
                 ctrlPtPaint.setColor(rand.nextU() | 0xFF000000);
@@ -336,7 +331,7 @@ protected:
                     grPaint.setXPFactory(GrPorterDuffXPFactory::Get(SkBlendMode::kSrc));
 
                     std::unique_ptr<GrMeshDrawOp> op =
-                            BezierCubicOrConicTestOp::Make(gp, bounds, color, klmEqs, 1.f);
+                            BezierCubicOrConicTestOp::Make(gp, bounds, color, klm, 1.f);
 
                     renderTargetContext->priv().testingOnly_addMeshDrawOp(
                             std::move(grPaint), GrAAType::kNone, std::move(op));
diff --git a/gn/bench.gni b/gn/bench.gni
index cf6ce881b3..d2d233fe55 100644
--- a/gn/bench.gni
+++ b/gn/bench.gni
@@ -34,6 +34,7 @@ bench_sources = [
   "$_bench/ColorPrivBench.cpp",
   "$_bench/ControlBench.cpp",
   "$_bench/CoverageBench.cpp",
+  "$_bench/CubicKLMBench.cpp",
   "$_bench/DashBench.cpp",
   "$_bench/DisplacementBench.cpp",
   "$_bench/DrawBitmapAABench.cpp",
diff --git a/src/gpu/GrPathUtils.cpp b/src/gpu/GrPathUtils.cpp
index 0bc84795b3..6dbac62fc1 100644
--- a/src/gpu/GrPathUtils.cpp
+++ b/src/gpu/GrPathUtils.cpp
@@ -323,7 +323,8 @@ void GrPathUtils::QuadUVMatrix::set(const SkPoint qPts[3]) {
 // k = (y2 - y0, x0 - x2, x2*y0 - x0*y2)
 // l = (y1 - y0, x0 - x1, x1*y0 - x0*y1) * 2*w
 // m = (y2 - y1, x1 - x2, x2*y1 - x1*y2) * 2*w
-void GrPathUtils::getConicKLM(const SkPoint p[3], const SkScalar weight, SkScalar klm[9]) {
+void GrPathUtils::getConicKLM(const SkPoint p[3], const SkScalar weight, SkMatrix* out) {
+    SkMatrix& klm = *out;
     const SkScalar w2 = 2.f * weight;
     klm[0] = p[2].fY - p[0].fY;
     klm[1] = p[0].fX - p[2].fX;
@@ -575,157 +576,272 @@ void GrPathUtils::convertCubicToQuadsConstrainToTangents(const SkPoint p[4],
 
 ////////////////////////////////////////////////////////////////////////////////
 
-// Solves linear system to extract klm
-// P.K = k (similarly for l, m)
-// Where P is matrix of control points
-// K is coefficients for the line K
-// k is vector of values of K evaluated at the control points
-// Solving for K, thus K = P^(-1) . k
-static void calc_cubic_klm(const SkPoint p[4], const SkScalar controlK[4],
-                           const SkScalar controlL[4], const SkScalar controlM[4],
-                           SkScalar k[3], SkScalar l[3], SkScalar m[3]) {
-    SkMatrix matrix;
-    matrix.setAll(p[0].fX, p[0].fY, 1.f,
-                  p[1].fX, p[1].fY, 1.f,
-                  p[2].fX, p[2].fY, 1.f);
-    SkMatrix inverse;
-    if (matrix.invert(&inverse)) {
-       inverse.mapHomogeneousPoints(k, controlK, 1);
-       inverse.mapHomogeneousPoints(l, controlL, 1);
-       inverse.mapHomogeneousPoints(m, controlM, 1);
+/**
+ * Computes an SkMatrix that can find the cubic KLM functionals as follows:
+ *
+ *     | ..K.. |   | ..kcoeffs.. |
+ *     | ..L.. | = | ..lcoeffs.. | * inverse_transpose_power_basis_matrix
+ *     | ..M.. |   | ..mcoeffs.. |
+ *
+ * 'kcoeffs' are the power basis coefficients to a scalar valued cubic function that returns the
+ * signed distance to line K from a given point on the curve:
+ *
+ *     k(t,s) = C(t,s) * K   [C(t,s) is defined in the following comment]
+ *
+ * The same applies for lcoeffs and mcoeffs. These are found separately, depending on the type of
+ * curve. There are 4 coefficients but 3 rows in the matrix, so in order to do this calculation the
+ * caller must first remove a specific column of coefficients.
+ *
+ * @return which column of klm coefficients to exclude from the calculation.
+ */
+static int calc_inverse_transpose_power_basis_matrix(const SkPoint pts[4], SkMatrix* out) {
+    using SkScalar4 = SkNx<4, SkScalar>;
+
+    // First we convert the bezier coordinates 'pts' to power basis coefficients X,Y,W=[0 0 0 1].
+    // M3 is the matrix that does this conversion. The homogeneous equation for the cubic becomes:
+    //
+    //                                     | X   Y   0 |
+    // C(t,s) = [t^3  t^2*s  t*s^2  s^3] * | .   .   0 |
+    //                                     | .   .   0 |
+    //                                     | .   .   1 |
+    //
+    const SkScalar4 M3[3] = {SkScalar4(-1, 3, -3, 1),
+                             SkScalar4(3, -6, 3, 0),
+                             SkScalar4(-3, 3, 0, 0)};
+    // 4th column of M3   =  SkScalar4(1, 0, 0, 0)};
+    SkScalar4 X(pts[3].x(), 0, 0, 0);
+    SkScalar4 Y(pts[3].y(), 0, 0, 0);
+    for (int i = 2; i >= 0; --i) {
+        X += M3[i] * pts[i].x();
+        Y += M3[i] * pts[i].y();
     }
 
+    // The matrix is 3x4. In order to invert it, we first need to make it square by throwing out one
+    // of the top three rows. We toss the row that leaves us with the largest determinant. Since the
+    // right column will be [0 0 1], the determinant reduces to x0*y1 - y0*x1.
+    SkScalar det[4];
+    SkScalar4 DETX1 = SkNx_shuffle<1,0,0,3>(X), DETY1 = SkNx_shuffle<1,0,0,3>(Y);
+    SkScalar4 DETX2 = SkNx_shuffle<2,2,1,3>(X), DETY2 = SkNx_shuffle<2,2,1,3>(Y);
+    (DETX1 * DETY2 - DETY1 * DETX2).store(det);
+    const int skipRow = det[0] > det[2] ? (det[0] > det[1] ? 0 : 1)
+                                        : (det[1] > det[2] ? 1 : 2);
+    const SkScalar rdet = 1 / det[skipRow];
+    const int row0 = (0 != skipRow) ? 0 : 1;
+    const int row1 = (2 == skipRow) ? 1 : 2;
+
+    // Compute the inverse-transpose of the power basis matrix with the 'skipRow'th row removed.
+    // Since W=[0 0 0 1], it follows that our corresponding solution will be equal to:
+    //
+    //             |  y1  -x1   x1*y2 - y1*x2 |
+    //     1/det * | -y0   x0  -x0*y2 + y0*x2 |
+    //             |   0    0             det |
+    //
+    const SkScalar4 R(rdet, rdet, rdet, 1);
+    X *= R;
+    Y *= R;
+
+    SkScalar x[4], y[4], z[4];
+    X.store(x);
+    Y.store(y);
+    (X * SkNx_shuffle<3,3,3,3>(Y) - Y * SkNx_shuffle<3,3,3,3>(X)).store(z);
+
+    out->setAll( y[row1], -x[row1],  z[row1],
+                -y[row0],  x[row0], -z[row0],
+                       0,        0,        1);
+
+    return skipRow;
 }
 
-static void set_serp_klm(const SkScalar d[3], SkScalar k[4], SkScalar l[4], SkScalar m[4]) {
-    SkScalar tempSqrt = SkScalarSqrt(9.f * d[1] * d[1] - 12.f * d[0] * d[2]);
-    SkScalar ls = 3.f * d[1] - tempSqrt;
-    SkScalar lt = 6.f * d[0];
-    SkScalar ms = 3.f * d[1] + tempSqrt;
-    SkScalar mt = 6.f * d[0];
-
-    k[0] = ls * ms;
-    k[1] = (3.f * ls * ms - ls * mt - lt * ms) / 3.f;
-    k[2] = (lt * (mt - 2.f * ms) + ls * (3.f * ms - 2.f * mt)) / 3.f;
-    k[3] = (lt - ls) * (mt - ms);
-
-    l[0] = ls * ls * ls;
-    const SkScalar lt_ls = lt - ls;
-    l[1] = ls * ls * lt_ls * -1.f;
-    l[2] = lt_ls * lt_ls * ls;
-    l[3] = -1.f * lt_ls * lt_ls * lt_ls;
-
-    m[0] = ms * ms * ms;
-    const SkScalar mt_ms = mt - ms;
-    m[1] = ms * ms * mt_ms * -1.f;
-    m[2] = mt_ms * mt_ms * ms;
-    m[3] = -1.f * mt_ms * mt_ms * mt_ms;
+static void negate_kl(SkMatrix* klm) {
+    // We could use klm->postScale(-1, -1), but it ends up doing a full matrix multiply.
+    for (int i = 0; i < 6; ++i) {
+        (*klm)[i] = -(*klm)[i];
+    }
+}
+
+static void calc_serp_klm(const SkPoint pts[4], const SkScalar d[3], SkMatrix* klm) {
+    SkMatrix CIT;
+    int skipCol = calc_inverse_transpose_power_basis_matrix(pts, &CIT);
+
+    const SkScalar root = SkScalarSqrt(9 * d[1] * d[1] - 12 * d[0] * d[2]);
+
+    const SkScalar tl = 3 * d[1] + root;
+    const SkScalar sl = 6 * d[0];
+    const SkScalar tm = 3 * d[1] - root;
+    const SkScalar sm = 6 * d[0];
+
+    SkMatrix klmCoeffs;
+    int col = 0;
+    if (0 != skipCol) {
+        klmCoeffs[0] = 0;
+        klmCoeffs[3] = -sl * sl * sl;
+        klmCoeffs[6] = -sm * sm * sm;
+        ++col;
+    }
+    if (1 != skipCol) {
+        klmCoeffs[col + 0] = sl * sm;
+        klmCoeffs[col + 3] = 3 * sl * sl * tl;
+        klmCoeffs[col + 6] = 3 * sm * sm * tm;
+        ++col;
+    }
+    if (2 != skipCol) {
+        klmCoeffs[col + 0] = -tl * sm - tm * sl;
+        klmCoeffs[col + 3] = -3 * sl * tl * tl;
+        klmCoeffs[col + 6] = -3 * sm * tm * tm;
+        ++col;
+    }
+
+    SkASSERT(2 == col);
+    klmCoeffs[2] = tl * tm;
+    klmCoeffs[5] = tl * tl * tl;
+    klmCoeffs[8] = tm * tm * tm;
+
+    klm->setConcat(klmCoeffs, CIT);
 
     // If d0 > 0 we need to flip the orientation of our curve
     // This is done by negating the k and l values
     // We want negative distance values to be on the inside
-    if ( d[0] > 0) {
-        for (int i = 0; i < 4; ++i) {
-            k[i] = -k[i];
-            l[i] = -l[i];
-        }
+    if (d[0] > 0) {
+        negate_kl(klm);
     }
 }
 
-static void set_loop_klm(const SkScalar d[3], SkScalar k[4], SkScalar l[4], SkScalar m[4]) {
-    SkScalar tempSqrt = SkScalarSqrt(4.f * d[0] * d[2] - 3.f * d[1] * d[1]);
-    SkScalar ls = d[1] - tempSqrt;
-    SkScalar lt = 2.f * d[0];
-    SkScalar ms = d[1] + tempSqrt;
-    SkScalar mt = 2.f * d[0];
-
-    k[0] = ls * ms;
-    k[1] = (3.f * ls*ms - ls * mt - lt * ms) / 3.f;
-    k[2] = (lt * (mt - 2.f * ms) + ls * (3.f * ms - 2.f * mt)) / 3.f;
-    k[3] = (lt - ls) * (mt - ms);
-
-    l[0] = ls * ls * ms;
-    l[1] = (ls * (ls * (mt - 3.f * ms) + 2.f * lt * ms))/-3.f;
-    l[2] = ((lt - ls) * (ls * (2.f * mt - 3.f * ms) + lt * ms))/3.f;
-    l[3] = -1.f * (lt - ls) * (lt - ls) * (mt - ms);
+static void calc_loop_klm(const SkPoint pts[4], SkScalar d1, SkScalar td, SkScalar sd,
+                          SkScalar te, SkScalar se, SkMatrix* klm) {
+    SkMatrix CIT;
+    int skipCol = calc_inverse_transpose_power_basis_matrix(pts, &CIT);
+
+    const SkScalar tesd = te * sd;
+    const SkScalar tdse = td * se;
+
+    SkMatrix klmCoeffs;
+    int col = 0;
+    if (0 != skipCol) {
+        klmCoeffs[0] = 0;
+        klmCoeffs[3] = -sd * sd * se;
+        klmCoeffs[6] = -se * se * sd;
+        ++col;
+    }
+    if (1 != skipCol) {
+        klmCoeffs[col + 0] = sd * se;
+        klmCoeffs[col + 3] = sd * (2 * tdse + tesd);
+        klmCoeffs[col + 6] = se * (2 * tesd + tdse);
+        ++col;
+    }
+    if (2 != skipCol) {
+        klmCoeffs[col + 0] = -tdse - tesd;
+        klmCoeffs[col + 3] = -td * (tdse + 2 * tesd);
+        klmCoeffs[col + 6] = -te * (tesd + 2 * tdse);
+        ++col;
+    }
 
-    m[0] = ls * ms * ms;
-    m[1] = (ms * (ls * (2.f * mt - 3.f * ms) + lt * ms))/-3.f;
-    m[2] = ((mt - ms) * (ls * (mt - 3.f * ms) + 2.f * lt * ms))/3.f;
-    m[3] = -1.f * (lt - ls) * (mt - ms) * (mt - ms);
+    SkASSERT(2 == col);
+    klmCoeffs[2] = td * te;
+    klmCoeffs[5] = td * td * te;
+    klmCoeffs[8] = te * te * td;
 
+    klm->setConcat(klmCoeffs, CIT);
 
     // For the general loop curve, we flip the orientation in the same pattern as the serp case
-    // above. Thus we only check d[0]. Technically we should check the value of the hessian as well
-    // cause we care about the sign of d[0]*Hessian. However, the Hessian is always negative outside
+    // above. Thus we only check d1. Technically we should check the value of the hessian as well
+    // cause we care about the sign of d1*Hessian. However, the Hessian is always negative outside
     // the loop section and positive inside. We take care of the flipping for the loop sections
     // later on.
-    if (d[0] > 0) {
-        for (int i = 0; i < 4; ++i) {
-            k[i] = -k[i];
-            l[i] = -l[i];
-        }
+    if (d1 > 0) {
+        negate_kl(klm);
     }
 }
 
-static void set_cusp_klm(const SkScalar d[3], SkScalar k[4], SkScalar l[4], SkScalar m[4]) {
-    const SkScalar ls = d[2];
-    const SkScalar lt = 3.f * d[1];
-
-    k[0] = ls;
-    k[1] = ls - lt / 3.f;
-    k[2] = ls - 2.f * lt / 3.f;
-    k[3] = ls - lt;
-
-    l[0] = ls * ls * ls;
-    const SkScalar ls_lt = ls - lt;
-    l[1] = ls * ls * ls_lt;
-    l[2] = ls_lt * ls_lt * ls;
-    l[3] = ls_lt * ls_lt * ls_lt;
-
-    m[0] = 1.f;
-    m[1] = 1.f;
-    m[2] = 1.f;
-    m[3] = 1.f;
+// For the case when we have a cusp at a parameter value of infinity (discr == 0, d1 == 0).
+static void calc_inf_cusp_klm(const SkPoint pts[4], SkScalar d2, SkScalar d3, SkMatrix* klm) {
+    SkMatrix CIT;
+    int skipCol = calc_inverse_transpose_power_basis_matrix(pts, &CIT);
+
+    const SkScalar tn = d3;
+    const SkScalar sn = 3 * d2;
+
+    SkMatrix klmCoeffs;
+    int col = 0;
+    if (0 != skipCol) {
+        klmCoeffs[0] = 0;
+        klmCoeffs[3] = -sn * sn * sn;
+        ++col;
+    }
+    if (1 != skipCol) {
+        klmCoeffs[col + 0] = 0;
+        klmCoeffs[col + 3] = 3 * sn * sn * tn;
+        ++col;
+    }
+    if (2 != skipCol) {
+        klmCoeffs[col + 0] = -sn;
+        klmCoeffs[col + 3] = -3 * sn * tn * tn;
+        ++col;
+    }
+
+    SkASSERT(2 == col);
+    klmCoeffs[2] = tn;
+    klmCoeffs[5] = tn * tn * tn;
+
+    klmCoeffs[6] = 0;
+    klmCoeffs[7] = 0;
+    klmCoeffs[8] = 1;
+
+    klm->setConcat(klmCoeffs, CIT);
 }
 
-// For the case when a cubic is actually a quadratic
-// M =
-// 0     0     0
-// 1/3   0     1/3
-// 2/3   1/3   2/3
-// 1     1     1
-static void set_quadratic_klm(const SkScalar d[3], SkScalar k[4], SkScalar l[4], SkScalar m[4]) {
-    k[0] = 0.f;
-    k[1] = 1.f/3.f;
-    k[2] = 2.f/3.f;
-    k[3] = 1.f;
-
-    l[0] = 0.f;
-    l[1] = 0.f;
-    l[2] = 1.f/3.f;
-    l[3] = 1.f;
-
-    m[0] = 0.f;
-    m[1] = 1.f/3.f;
-    m[2] = 2.f/3.f;
-    m[3] = 1.f;
-
-    // If d2 < 0 we need to flip the orientation of our curve since we want negative values to be on
-    // the "inside" of the curve. This is done by negating the k and l values
-    if ( d[2] > 0) {
-        for (int i = 0; i < 4; ++i) {
-            k[i] = -k[i];
-            l[i] = -l[i];
-        }
+// For the case when a cubic bezier is actually a quadratic. We duplicate k in l so that the
+// implicit becomes:
+//
+//     k^3 - l*m == k^3 - l*k == k * (k^2 - l)
+//
+// In the quadratic case we can simply assign fixed values at each control point:
+//
+//     | ..K.. |     | pts[0]  pts[1]  pts[2]  pts[3] |      | 0   1/3  2/3  1 |
+//     | ..L.. |  *  |   .       .       .       .    |  ==  | 0     0  1/3  1 |
+//     | ..K.. |     |   1       1       1       1    |      | 0   1/3  2/3  1 |
+//
+static void calc_quadratic_klm(const SkPoint pts[4], SkScalar d3, SkMatrix* klm) {
+    SkMatrix klmAtPts;
+    klmAtPts.setAll(0,  1.f/3,  1,
+                    0,      0,  1,
+                    0,  1.f/3,  1);
+
+    SkMatrix inversePts;
+    inversePts.setAll(pts[0].x(),  pts[1].x(),  pts[3].x(),
+                      pts[0].y(),  pts[1].y(),  pts[3].y(),
+                               1,           1,           1);
+    SkAssertResult(inversePts.invert(&inversePts));
+
+    klm->setConcat(klmAtPts, inversePts);
+
+    // If d3 > 0 we need to flip the orientation of our curve
+    // This is done by negating the k and l values
+    if (d3 > 0) {
+        negate_kl(klm);
     }
 }
 
-int GrPathUtils::chopCubicAtLoopIntersection(const SkPoint src[4], SkPoint dst[10], SkScalar klm[9],
+// For the case when a cubic bezier is actually a line. We set K=0, L=1, M=-line, which results in
+// the following implicit:
+//
+//     k^3 - l*m == 0^3 - 1*(-line) == -(-line) == line
+//
+static void calc_line_klm(const SkPoint pts[4], SkMatrix* klm) {
+    SkScalar ny = pts[0].x() - pts[3].x();
+    SkScalar nx = pts[3].y() - pts[0].y();
+    SkScalar k = nx * pts[0].x() + ny * pts[0].y();
+    klm->setAll(  0,   0, 0,
+                  0,   0, 1,
+                -nx, -ny, k);
+}
+
+int GrPathUtils::chopCubicAtLoopIntersection(const SkPoint src[4], SkPoint dst[10], SkMatrix* klm,
                                              int* loopIndex) {
-    // Variable to store the two parametric values at the loop double point
-    SkScalar smallS = 0.f;
-    SkScalar largeS = 0.f;
+    // Variables to store the two parametric values at the loop double point.
+    SkScalar t1 = 0, t2 = 0;
+
+    // Homogeneous parametric values at the loop double point.
+    SkScalar td, sd, te, se;
 
     SkScalar d[3];
     SkCubicType cType = SkClassifyCubic(src, d);
@@ -733,27 +849,24 @@ int GrPathUtils::chopCubicAtLoopIntersection(const SkPoint src[4], SkPoint dst[1
     int chop_count = 0;
     if (kLoop_SkCubicType == cType) {
         SkScalar tempSqrt = SkScalarSqrt(4.f * d[0] * d[2] - 3.f * d[1] * d[1]);
-        SkScalar ls = d[1] - tempSqrt;
-        SkScalar lt = 2.f * d[0];
-        SkScalar ms = d[1] + tempSqrt;
-        SkScalar mt = 2.f * d[0];
-        ls = ls / lt;
-        ms = ms / mt;
+        td = d[1] + tempSqrt;
+        sd = 2.f * d[0];
+        te = d[1] - tempSqrt;
+        se = 2.f * d[0];
+
+        t1 = td / sd;
+        t2 = te / se;
         // need to have t values sorted since this is what is expected by SkChopCubicAt
-        if (ls <= ms) {
-            smallS = ls;
-            largeS = ms;
-        } else {
-            smallS = ms;
-            largeS = ls;
+        if (t1 > t2) {
+            SkTSwap(t1, t2);
         }
 
         SkScalar chop_ts[2];
-        if (smallS > 0.f && smallS < 1.f) {
-            chop_ts[chop_count++] = smallS;
+        if (t1 > 0.f && t1 < 1.f) {
+            chop_ts[chop_count++] = t1;
         }
-        if (largeS > 0.f && largeS < 1.f) {
-            chop_ts[chop_count++] = largeS;
+        if (t2 > 0.f && t2 < 1.f) {
+            chop_ts[chop_count++] = t2;
         }
         if(dst) {
             SkChopCubicAt(src, dst, chop_ts, chop_count);
@@ -768,58 +881,45 @@ int GrPathUtils::chopCubicAtLoopIntersection(const SkPoint src[4], SkPoint dst[1
         if (2 == chop_count) {
             *loopIndex = 1;
         } else if (1 == chop_count) {
-            if (smallS < 0.f) {
+            if (t1 < 0.f) {
                 *loopIndex = 0;
             } else {
                 *loopIndex = 1;
             }
         } else {
-            if (smallS < 0.f && largeS > 1.f) {
+            if (t1 < 0.f && t2 > 1.f) {
                 *loopIndex = 0;
             } else {
                 *loopIndex = -1;
             }
         }
     }
-    if (klm) {
-        SkScalar controlK[4];
-        SkScalar controlL[4];
-        SkScalar controlM[4];
-
-        if (kSerpentine_SkCubicType == cType || (kCusp_SkCubicType == cType && 0.f != d[0])) {
-            set_serp_klm(d, controlK, controlL, controlM);
-        } else if (kLoop_SkCubicType == cType) {
-            set_loop_klm(d, controlK, controlL, controlM);
-        } else if (kCusp_SkCubicType == cType) {
-            SkASSERT(0.f == d[0]);
-            set_cusp_klm(d, controlK, controlL, controlM);
-        } else if (kQuadratic_SkCubicType == cType) {
-            set_quadratic_klm(d, controlK, controlL, controlM);
-        }
 
-        calc_cubic_klm(src, controlK, controlL, controlM, klm, &klm[3], &klm[6]);
+    if (klm) {
+        switch (cType) {
+            case kSerpentine_SkCubicType:
+                calc_serp_klm(src, d, klm);
+                break;
+            case kLoop_SkCubicType:
+                calc_loop_klm(src, d[0], td, sd, te, se, klm);
+                break;
+            case kCusp_SkCubicType:
+                if (0 != d[0]) {
+                    // FIXME: SkClassifyCubic has a tolerance, but we need an exact classification
+                    // here to be sure we won't get a negative in the square root.
+                    calc_serp_klm(src, d, klm);
+                } else {
+                    calc_inf_cusp_klm(src, d[1], d[2], klm);
+                }
+                break;
+            case kQuadratic_SkCubicType:
+                calc_quadratic_klm(src, d[2], klm);
+                break;
+            case kLine_SkCubicType:
+            case kPoint_SkCubicType:
+                calc_line_klm(src, klm);
+                break;
+        };
     }
     return chop_count + 1;
 }
-
-void GrPathUtils::getCubicKLM(const SkPoint p[4], SkScalar klm[9]) {
-    SkScalar d[3];
-    SkCubicType cType = SkClassifyCubic(p, d);
-
-    SkScalar controlK[4];
-    SkScalar controlL[4];
-    SkScalar controlM[4];
-
-    if (kSerpentine_SkCubicType == cType || (kCusp_SkCubicType == cType && 0.f != d[0])) {
-        set_serp_klm(d, controlK, controlL, controlM);
-    } else if (kLoop_SkCubicType == cType) {
-        set_loop_klm(d, controlK, controlL, controlM);
-    } else if (kCusp_SkCubicType == cType) {
-        SkASSERT(0.f == d[0]);
-        set_cusp_klm(d, controlK, controlL, controlM);
-    } else if (kQuadratic_SkCubicType == cType) {
-        set_quadratic_klm(d, controlK, controlL, controlM);
-    }
-
-    calc_cubic_klm(p, controlK, controlL, controlM, klm, &klm[3], &klm[6]);
-}
diff --git a/src/gpu/GrPathUtils.h b/src/gpu/GrPathUtils.h
index 8b844180b7..7dea3a1353 100644
--- a/src/gpu/GrPathUtils.h
+++ b/src/gpu/GrPathUtils.h
@@ -98,13 +98,15 @@ namespace GrPathUtils {
 
     // Input is 3 control points and a weight for a bezier conic. Calculates the
     // three linear functionals (K,L,M) that represent the implicit equation of the
-    // conic, K^2 - LM.
+    // conic, k^2 - lm.
     //
-    // Output:
-    //  K = (klm[0], klm[1], klm[2])
-    //  L = (klm[3], klm[4], klm[5])
-    //  M = (klm[6], klm[7], klm[8])
-    void getConicKLM(const SkPoint p[3], const SkScalar weight, SkScalar klm[9]);
+    // Output: klm holds the linear functionals K,L,M as row vectors:
+    //
+    //     | ..K.. |   | x |      | k |
+    //     | ..L.. | * | y |  ==  | l |
+    //     | ..M.. |   | 1 |      | m |
+    //
+    void getConicKLM(const SkPoint p[3], const SkScalar weight, SkMatrix* klm);
 
     // Converts a cubic into a sequence of quads. If working in device space
     // use tolScale = 1, otherwise set based on stretchiness of the matrix. The
@@ -127,48 +129,39 @@ namespace GrPathUtils {
 
     // Chops the cubic bezier passed in by src, at the double point (intersection point)
     // if the curve is a cubic loop. If it is a loop, there will be two parametric values for
-    // the double point: ls and ms. We chop the cubic at these values if they are between 0 and 1.
+    // the double point: t1 and t2. We chop the cubic at these values if they are between 0 and 1.
     // Return value:
-    // Value of 3: ls and ms are both between (0,1), and dst will contain the three cubics,
+    // Value of 3: t1 and t2 are both between (0,1), and dst will contain the three cubics,
     //             dst[0..3], dst[3..6], and dst[6..9] if dst is not nullptr
-    // Value of 2: Only one of ls and ms are between (0,1), and dst will contain the two cubics,
+    // Value of 2: Only one of t1 and t2 are between (0,1), and dst will contain the two cubics,
     //             dst[0..3] and dst[3..6] if dst is not nullptr
-    // Value of 1: Neither ls or ms are between (0,1), and dst will contain the one original cubic,
+    // Value of 1: Neither t1 nor t2 are between (0,1), and dst will contain the one original cubic,
     //             dst[0..3] if dst is not nullptr
     //
     // Optional KLM Calculation:
-    // The function can also return the KLM linear functionals for the chopped cubic implicit form
-    // of K^3 - LM.
-    // It will calculate a single set of KLM values that can be shared by all sub cubics, except
-    // for the subsection that is "the loop" the K and L values need to be negated.
+    // The function can also return the KLM linear functionals for the cubic implicit form of
+    // k^3 - lm. This can be shared by all chopped cubics.
+    //
     // Output:
-    // klm:     Holds the values for the linear functionals as:
-    //          K = (klm[0], klm[1], klm[2])
-    //          L = (klm[3], klm[4], klm[5])
-    //          M = (klm[6], klm[7], klm[8])
-    // loopIndex: This value will tell the caller which of the chopped sections are the the actual
-    //            loop once we've chopped. A value of -1 means there is no loop section. The caller
-    //            can then use this value to decide how/if they want to flip the orientation of this
-    //            section. This flip should be done by negating the K and L values.
     //
-    // Notice that the klm lines are calculated in the same space as the input control points.
-    // If you transform the points the lines will also need to be transformed. This can be done
-    // by mapping the lines with the inverse-transpose of the matrix used to map the points.
-    int chopCubicAtLoopIntersection(const SkPoint src[4], SkPoint dst[10] = nullptr,
-                                    SkScalar klm[9] = nullptr, int* loopIndex = nullptr);
-
-    // Input is p which holds the 4 control points of a non-rational cubic Bezier curve.
-    // Output is the coefficients of the three linear functionals K, L, & M which
-    // represent the implicit form of the cubic as f(x,y,w) = K^3 - LM. The w term
-    // will always be 1. The output is stored in the array klm, where the values are:
-    // K = (klm[0], klm[1], klm[2])
-    // L = (klm[3], klm[4], klm[5])
-    // M = (klm[6], klm[7], klm[8])
+    // klm: Holds the linear functionals K,L,M as row vectors:
+    //
+    //          | ..K.. |   | x |      | k |
+    //          | ..L.. | * | y |  ==  | l |
+    //          | ..M.. |   | 1 |      | m |
+    //
+    // loopIndex: This value will tell the caller which of the chopped sections (if any) are the
+    //            actual loop. A value of -1 means there is no loop section. The caller can then use
+    //            this value to decide how/if they want to flip the orientation of this section.
+    //            The flip should be done by negating the k and l values as follows:
     //
-    // Notice that the klm lines are calculated in the same space as the input control points.
+    //                KLM.postScale(-1, -1)
+    //
+    // Notice that the KLM lines are calculated in the same space as the input control points.
     // If you transform the points the lines will also need to be transformed. This can be done
     // by mapping the lines with the inverse-transpose of the matrix used to map the points.
-    void getCubicKLM(const SkPoint p[4], SkScalar klm[9]);
+    int chopCubicAtLoopIntersection(const SkPoint src[4], SkPoint dst[10] = nullptr,
+                                    SkMatrix* klm = nullptr, int* loopIndex = nullptr);
 
     // When tessellating curved paths into linear segments, this defines the maximum distance
     // in screen space which a segment may deviate from the mathmatically correct value.
diff --git a/src/gpu/ops/GrAAHairLinePathRenderer.cpp b/src/gpu/ops/GrAAHairLinePathRenderer.cpp
index e373b97d05..ec2610400f 100644
--- a/src/gpu/ops/GrAAHairLinePathRenderer.cpp
+++ b/src/gpu/ops/GrAAHairLinePathRenderer.cpp
@@ -409,9 +409,7 @@ struct BezierVertex {
     SkPoint fPos;
     union {
         struct {
-            SkScalar fK;
-            SkScalar fL;
-            SkScalar fM;
+            SkScalar fKLM[3];
         } fConic;
         SkVector   fQuadCoord;
         struct {
@@ -526,15 +524,13 @@ static void bloat_quad(const SkPoint qpts[3], const SkMatrix* toDevice,
 // k, l, m are calculated in function GrPathUtils::getConicKLM
 static void set_conic_coeffs(const SkPoint p[3], BezierVertex verts[kQuadNumVertices],
                              const SkScalar weight) {
-    SkScalar klm[9];
+    SkMatrix klm;
 
-    GrPathUtils::getConicKLM(p, weight, klm);
+    GrPathUtils::getConicKLM(p, weight, &klm);
 
     for (int i = 0; i < kQuadNumVertices; ++i) {
-        const SkPoint pnt = verts[i].fPos;
-        verts[i].fConic.fK = pnt.fX * klm[0] + pnt.fY * klm[1] + klm[2];
-        verts[i].fConic.fL = pnt.fX * klm[3] + pnt.fY * klm[4] + klm[5];
-        verts[i].fConic.fM = pnt.fX * klm[6] + pnt.fY * klm[7] + klm[8];
+        const SkScalar pt3[3] = {verts[i].fPos.x(), verts[i].fPos.y(), 1.f};
+        klm.mapHomogeneousPoints(verts[i].fConic.fKLM, pt3, 1);
     }
 }