1 files changed, 76 insertions, 0 deletions
diff --git a/experimental/Intersection/CubicUtilities.cpp b/experimental/Intersection/CubicUtilities.cpp
new file mode 100644
index 0000000000..3fab29ec80
--- /dev/null
+++ b/experimental/Intersection/CubicUtilities.cpp
@@ -0,0 +1,76 @@
+#include "CubicUtilities.h"
+#include "DataTypes.h"
+#include "QuadraticUtilities.h"
+
+void coefficients(const double* cubic, double& A, double& B, double& C, double& D) {
+    A = cubic[6]; // d
+    B = cubic[4] * 3; // 3*c
+    C = cubic[2] * 3; // 3*b
+    D = cubic[0]; // a
+    A -= D - C + B;     // A =   -a + 3*b - 3*c + d
+    B += 3 * D - 2 * C; // B =  3*a - 6*b + 3*c
+    C -= 3 * D;         // C = -3*a + 3*b
+}
+
+// cubic roots
+
+const double PI = 4 * atan(1);
+
+static bool is_unit_interval(double x) {
+    return x > 0 && x < 1;
+}
+
+// from SkGeometry.cpp (and Numeric Solutions, 5.6)
+int cubicRoots(double A, double B, double C, double D, double t[3]) {
+    if (approximately_zero(A)) {  // we're just a quadratic
+        return quadraticRoots(B, C, D, t);
+    }
+    double a, b, c;
+    {
+        double invA = 1 / A;
+        a = B * invA;
+        b = C * invA;
+        c = D * invA;
+    }
+    double a2 = a * a;
+    double Q = (a2 - b * 3) / 9;
+    double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
+    double Q3 = Q * Q * Q;
+    double R2MinusQ3 = R * R - Q3;
+    double adiv3 = a / 3;
+    double* roots = t;
+    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 - t);
+}