[base] Add SkBezierCubic::subdivide

This extracts SkDCubic::chopAt [1] into base and tests to
make sure the two implementations are identical. This also
adds documentation.

It's not clear to me if the "t == 0.5" conditional in the
original is a "fast path" or something for float precision
or something else. In any case, there is no explanation and
no unit tests that fail, so maybe it's ok.

[1] https://github.com/google/skia/blob/8106368d1f70d24e612afe7f6de8b6e6d1c6dd27/src/pathops/SkPathOpsCubic.cpp#L110

Change-Id: Iaa71b42b0b0b78259254711b54e8188c68fee63f
Bug: skia:13983
Reviewed-on: https://skia-review.googlesource.com/c/skia/+/632583
Reviewed-by: Tyler Denniston <tdenniston@google.com>
Reviewed-by: Michael Ludwig <michaelludwig@google.com>
diff --git a/gn/core.gni b/gn/core.gni
index a8c302b..e456219 100644
--- a/gn/core.gni
+++ b/gn/core.gni
@@ -199,6 +199,8 @@
   "$_src/base/SkArenaAlloc.h",
   "$_src/base/SkArenaAllocList.h",
   "$_src/base/SkAutoMalloc.h",
+  "$_src/base/SkBezierCurves.cpp",
+  "$_src/base/SkBezierCurves.h",
   "$_src/base/SkBlockAllocator.cpp",
   "$_src/base/SkBlockAllocator.h",
   "$_src/base/SkBuffer.cpp",
diff --git a/public.bzl b/public.bzl
index 1ed5150..5eaf6a6 100644
--- a/public.bzl
+++ b/public.bzl
@@ -309,6 +309,8 @@
     "src/base/SkArenaAlloc.h",
     "src/base/SkArenaAllocList.h",
     "src/base/SkAutoMalloc.h",
+    "src/base/SkBezierCurves.cpp",
+    "src/base/SkBezierCurves.h",
     "src/base/SkBlockAllocator.cpp",
     "src/base/SkBlockAllocator.h",
     "src/base/SkBuffer.cpp",
diff --git a/src/base/BUILD.bazel b/src/base/BUILD.bazel
index fec86a4..38ce287 100644
--- a/src/base/BUILD.bazel
+++ b/src/base/BUILD.bazel
@@ -29,6 +29,7 @@
     srcs = IWYU_HDRS + [
         "SkArenaAlloc.h",
         "SkAutoMalloc.h",
+        "SkBezierCurves.h",
         "SkBlockAllocator.h",
         "SkBuffer.h",
         "SkCubics.h",
@@ -68,6 +69,7 @@
 skia_filegroup(
     name = "srcs",
     srcs = [
+        "SkBezierCurves.cpp",
         "SkBuffer.cpp",
         "SkCubics.cpp",
         "SkDeque.cpp",
diff --git a/src/base/SkBezierCurves.cpp b/src/base/SkBezierCurves.cpp
new file mode 100644
index 0000000..95679b8
--- /dev/null
+++ b/src/base/SkBezierCurves.cpp
@@ -0,0 +1,58 @@
+/*
+ * Copyright 2012 Google LLC
+ *
+ * Use of this source code is governed by a BSD-style license that can be
+ * found in the LICENSE file.
+ */
+
+#include "src/base/SkBezierCurves.h"
+
+#include "include/private/base/SkAssert.h"
+
+#include <cstddef>
+
+static inline double interpolate(double A, double B, double t) {
+    return A + (B - A) * t;
+}
+
+// Perform subdivision using De Casteljau's algorithm, that is, repeated linear
+// interpolation between adjacent points.
+void SkBezierCubic::Subdivide(const double curve[8], double t,
+                              double twoCurves[14]) {
+    SkASSERT(0.0 <= t && t <= 1.0);
+    // We split the curve "in" into two curves "alpha" and "beta"
+    const auto in_X = [&curve](size_t n) { return curve[2*n]; };
+    const auto in_Y = [&curve](size_t n) { return curve[2*n + 1]; };
+    const auto alpha_X = [&twoCurves](size_t n) -> double& { return twoCurves[2*n]; };
+    const auto alpha_Y = [&twoCurves](size_t n) -> double& { return twoCurves[2*n + 1]; };
+    const auto beta_X = [&twoCurves](size_t n) -> double& { return twoCurves[2*n + 6]; };
+    const auto beta_Y = [&twoCurves](size_t n) -> double& { return twoCurves[2*n + 7]; };
+
+    alpha_X(0) = in_X(0);
+    alpha_Y(0) = in_Y(0);
+
+    beta_X(3) = in_X(3);
+    beta_Y(3) = in_Y(3);
+
+    double x01 = interpolate(in_X(0), in_X(1), t);
+    double y01 = interpolate(in_Y(0), in_Y(1), t);
+    double x12 = interpolate(in_X(1), in_X(2), t);
+    double y12 = interpolate(in_Y(1), in_Y(2), t);
+    double x23 = interpolate(in_X(2), in_X(3), t);
+    double y23 = interpolate(in_Y(2), in_Y(3), t);
+
+    alpha_X(1) = x01;
+    alpha_Y(1) = y01;
+
+    beta_X(2) = x23;
+    beta_Y(2) = y23;
+
+    alpha_X(2) = interpolate(x01, x12, t);
+    alpha_Y(2) = interpolate(y01, y12, t);
+
+    beta_X(1) = interpolate(x12, x23, t);
+    beta_Y(1) = interpolate(y12, y23, t);
+
+    alpha_X(3) /*= beta_X(0) */ = interpolate(alpha_X(2), beta_X(1), t);
+    alpha_Y(3) /*= beta_Y(0) */ = interpolate(alpha_Y(2), beta_Y(1), t);
+}
diff --git a/src/base/SkBezierCurves.h b/src/base/SkBezierCurves.h
new file mode 100644
index 0000000..c2e8b84
--- /dev/null
+++ b/src/base/SkBezierCurves.h
@@ -0,0 +1,38 @@
+/*
+ * Copyright 2023 Google LLC
+ *
+ * Use of this source code is governed by a BSD-style license that can be
+ * found in the LICENSE file.
+ */
+#ifndef SkBezierCurves_DEFINED
+#define SkBezierCurves_DEFINED
+
+/**
+ * Utilities for dealing with cubic Bézier curves. These have a start XY
+ * point, an end XY point, and two control XY points in between. They take
+ * a parameter t which is between 0 and 1 (inclusive) which is used to
+ * interpolate between the start and end points, via a route dictated by
+ * the control points, and return a new XY point.
+ *
+ * We store a Bézier curve as an array of 8 floats or doubles, where
+ * the even indices are the X coordinates, and the odd indices are the Y
+ * coordinates.
+ */
+class SkBezierCubic {
+public:
+    /**
+     * Splits the provided curve at the location t, resulting in two
+     * bezier curves that share a point (the end point from curve 1
+     * and the start point from curve 2 are the same).
+     *
+     * t must be in the interval [0, 1].
+     *
+     * The provided twoCurves array will be filled such that indices
+     * 0-7 are the first curve (representing the interval [0, t]), and
+     * indices 6-13 are the second curve (representing [t, 1]).
+     */
+    static void Subdivide(const double curve[8], double t,
+                          double twoCurves[14]);
+};
+
+#endif
diff --git a/tests/CubicChopTest.cpp b/tests/CubicChopTest.cpp
index 07f8d43..ed2fe77 100644
--- a/tests/CubicChopTest.cpp
+++ b/tests/CubicChopTest.cpp
@@ -8,6 +8,7 @@
 #include "include/core/SkSpan.h"
 #include "include/core/SkTypes.h"
 #include "include/private/base/SkFloatingPoint.h"
+#include "src/base/SkBezierCurves.h"
 #include "src/pathops/SkPathOpsCubic.h"
 #include "src/pathops/SkPathOpsPoint.h"
 #include "tests/Test.h"
@@ -15,7 +16,8 @@
 #include <cstring>
 #include <string>
 
-struct DoublePoint{
+// Grouping the test inputs into DoublePoints makes the test cases easier to read.
+struct DoublePoint {
     double x;
     double y;
 };
@@ -40,17 +42,36 @@
                     "Invalid test case. Should have 7 output points.");
 
 
-    SkDCubic input;
-    std::memcpy(&input.fPts[0], curveInputs.begin(), 8 * sizeof(double));
-    SkDCubicPair output = input.chopAt(t);
+    {
+        skiatest::ReporterContext subsubtest(reporter, "Pathops Implementation");
+        SkDCubic input;
+        std::memcpy(&input.fPts[0], curveInputs.begin(), 8 * sizeof(double));
+        SkDCubicPair output = input.chopAt(t);
 
-    for (int i = 0; i < 7; ++i) {
-        REPORTER_ASSERT(reporter,
-                        nearly_equal(expectedOutputs[i].x, output.pts[i].fX) &&
-                        nearly_equal(expectedOutputs[i].y, output.pts[i].fY),
-                        "(%.16f, %.16f) != (%.16f, %.16f) at index %d",
-                        expectedOutputs[i].x, expectedOutputs[i].y,
-                        output.pts[i].fX, output.pts[i].fY, i);
+        for (int i = 0; i < 7; ++i) {
+            REPORTER_ASSERT(reporter,
+                            nearly_equal(expectedOutputs[i].x, output.pts[i].fX) &&
+                            nearly_equal(expectedOutputs[i].y, output.pts[i].fY),
+                            "(%.16f, %.16f) != (%.16f, %.16f) at index %d",
+                            expectedOutputs[i].x, expectedOutputs[i].y,
+                            output.pts[i].fX, output.pts[i].fY, i);
+        }
+    }
+    {
+        skiatest::ReporterContext subsubtest(reporter, "SkBezier Implementation");
+        double input[8];
+        double output[14];
+        std::memcpy(input, curveInputs.begin(), 8 * sizeof(double));
+        SkBezierCubic::Subdivide(input, t, output);
+
+        for (int i = 0; i < 7; ++i) {
+            REPORTER_ASSERT(reporter,
+                            nearly_equal(expectedOutputs[i].x, output[i*2]) &&
+                            nearly_equal(expectedOutputs[i].y, output[i*2 + 1]),
+                            "(%.16f, %.16f) != (%.16f, %.16f) at index %d",
+                            expectedOutputs[i].x, expectedOutputs[i].y,
+                            output[i*2], output[i*2 + 1], i);
+        }
     }
 }