diff --git a/gpu/compute.go b/gpu/compute.go
index ba10f824..0ab2f23b 100644
--- a/gpu/compute.go
+++ b/gpu/compute.go
@@ -706,6 +706,13 @@ func encodePath(pathData []byte, stroke clip.StrokeStyle, dashes dashOp) encoder
 			}
 			hasPrev = true
 			prevTo = to
+		case scene.OpFillCubic:
+			from, _, _, to := scene.DecodeCubic(cmd)
+			if hasPrev && from != prevTo {
+				enc.scene[len(enc.scene)-1][0] |= (flagEndPath << 16)
+			}
+			hasPrev = true
+			prevTo = to
 		default:
 			panic("unsupported path scene command")
 		}
diff --git a/gpu/gpu.go b/gpu/gpu.go
index 0ecc2f92..65686ce1 100644
--- a/gpu/gpu.go
+++ b/gpu/gpu.go
@@ -1404,6 +1404,11 @@ func decodeToOutlineQuads(qs *quadSplitter, tr f32.Affine2D, pathData []byte) {
 			q.From, q.Ctrl, q.To = scene.DecodeQuad(cmd)
 			q = q.Transform(tr)
 			qs.splitAndEncode(q)
+		case scene.OpFillCubic:
+			for _, q := range splitCubic(scene.DecodeCubic(cmd)) {
+				q = q.Transform(tr)
+				qs.splitAndEncode(q)
+			}
 		default:
 			panic("unsupported scene command")
 		}
@@ -1436,6 +1441,14 @@ func decodeToStrokeQuads(pathData []byte) strokeQuads {
 				quad:    q,
 			}
 			quads = append(quads, quad)
+		case scene.OpFillCubic:
+			for _, q := range splitCubic(scene.DecodeCubic(cmd)) {
+				quad := strokeQuad{
+					contour: contour,
+					quad:    q,
+				}
+				quads = append(quads, quad)
+			}
 		default:
 			panic("unsupported scene command")
 		}
@@ -1536,3 +1549,78 @@ func decodeQuad(d []byte) (q quadSegment) {
 	q.From, q.Ctrl, q.To = scene.DecodeQuad(cmd)
 	return
 }
+
+func splitCubic(from, ctrl0, ctrl1, to f32.Point) []quadSegment {
+	quads := make([]quadSegment, 0, 10)
+	// Set the maximum distance proportionally to the longest side
+	// of the bounding rectangle.
+	hull := f32.Rectangle{
+		Min: from,
+		Max: ctrl0,
+	}.Canon().Add(ctrl1).Add(to)
+	l := hull.Dx()
+	if h := hull.Dy(); h > l {
+		l = h
+	}
+	approxCubeTo(&quads, 0, l*0.001, from, ctrl0, ctrl1, to)
+	return quads
+}
+
+// approxCube approximates a cubic Bézier by a series of quadratic
+// curves.
+func approxCubeTo(quads *[]quadSegment, splits int, maxDist float32, from, ctrl0, ctrl1, to f32.Point) int {
+	// The idea is from
+	// https://caffeineowl.com/graphics/2d/vectorial/cubic2quad01.html
+	// where a quadratic approximates a cubic by eliminating its t³ term
+	// from its polynomial expression anchored at the starting point:
+	//
+	// P(t) = pen + 3t(ctrl0 - pen) + 3t²(ctrl1 - 2ctrl0 + pen) + t³(to - 3ctrl1 + 3ctrl0 - pen)
+	//
+	// The control point for the new quadratic Q1 that shares starting point, pen, with P is
+	//
+	// C1 = (3ctrl0 - pen)/2
+	//
+	// The reverse cubic anchored at the end point has the polynomial
+	//
+	// P'(t) = to + 3t(ctrl1 - to) + 3t²(ctrl0 - 2ctrl1 + to) + t³(pen - 3ctrl0 + 3ctrl1 - to)
+	//
+	// The corresponding quadratic Q2 that shares the end point, to, with P has control
+	// point
+	//
+	// C2 = (3ctrl1 - to)/2
+	//
+	// The combined quadratic Bézier, Q, shares both start and end points with its cubic
+	// and use the midpoint between the two curves Q1 and Q2 as control point:
+	//
+	// C = (3ctrl0 - pen + 3ctrl1 - to)/4
+	c := ctrl0.Mul(3).Sub(from).Add(ctrl1.Mul(3)).Sub(to).Mul(1.0 / 4.0)
+	const maxSplits = 32
+	if splits >= maxSplits {
+		*quads = append(*quads, quadSegment{From: from, Ctrl: c, To: to})
+		return splits
+	}
+	// The maximum distance between the cubic P and its approximation Q given t
+	// can be shown to be
+	//
+	// d = sqrt(3)/36*|to - 3ctrl1 + 3ctrl0 - pen|
+	//
+	// To save a square root, compare d² with the squared tolerance.
+	v := to.Sub(ctrl1.Mul(3)).Add(ctrl0.Mul(3)).Sub(from)
+	d2 := (v.X*v.X + v.Y*v.Y) * 3 / (36 * 36)
+	if d2 <= maxDist*maxDist {
+		*quads = append(*quads, quadSegment{From: from, Ctrl: c, To: to})
+		return splits
+	}
+	// De Casteljau split the curve and approximate the halves.
+	t := float32(0.5)
+	c0 := from.Add(ctrl0.Sub(from).Mul(t))
+	c1 := ctrl0.Add(ctrl1.Sub(ctrl0).Mul(t))
+	c2 := ctrl1.Add(to.Sub(ctrl1).Mul(t))
+	c01 := c0.Add(c1.Sub(c0).Mul(t))
+	c12 := c1.Add(c2.Sub(c1).Mul(t))
+	c0112 := c01.Add(c12.Sub(c01).Mul(t))
+	splits++
+	splits = approxCubeTo(quads, splits, maxDist, from, c0, c01, c0112)
+	splits = approxCubeTo(quads, splits, maxDist, c0112, c12, c2, to)
+	return splits
+}
diff --git a/internal/scene/scene.go b/internal/scene/scene.go
index db663bf7..11717344 100644
--- a/internal/scene/scene.go
+++ b/internal/scene/scene.go
@@ -56,6 +56,24 @@ func Line(start, end f32.Point, stroke bool, flags uint32) Command {
 	}
 }
 
+func Cubic(start, ctrl0, ctrl1, end f32.Point, stroke bool) Command {
+	tag := uint32(OpFillCubic)
+	if stroke {
+		tag = uint32(OpStrokeCubic)
+	}
+	return Command{
+		0: tag,
+		1: math.Float32bits(start.X),
+		2: math.Float32bits(start.Y),
+		3: math.Float32bits(ctrl0.X),
+		4: math.Float32bits(ctrl0.Y),
+		5: math.Float32bits(ctrl1.X),
+		6: math.Float32bits(ctrl1.Y),
+		7: math.Float32bits(end.X),
+		8: math.Float32bits(end.Y),
+	}
+}
+
 func Quad(start, ctrl, end f32.Point, stroke bool) Command {
 	tag := uint32(OpFillQuad)
 	if stroke {
@@ -151,3 +169,14 @@ func DecodeQuad(cmd Command) (from, ctrl, to f32.Point) {
 	to = f32.Pt(math.Float32frombits(cmd[5]), math.Float32frombits(cmd[6]))
 	return
 }
+
+func DecodeCubic(cmd Command) (from, ctrl0, ctrl1, to f32.Point) {
+	if cmd[0] != uint32(OpFillCubic) {
+		panic("invalid command")
+	}
+	from = f32.Pt(math.Float32frombits(cmd[1]), math.Float32frombits(cmd[2]))
+	ctrl0 = f32.Pt(math.Float32frombits(cmd[3]), math.Float32frombits(cmd[4]))
+	ctrl1 = f32.Pt(math.Float32frombits(cmd[5]), math.Float32frombits(cmd[6]))
+	to = f32.Pt(math.Float32frombits(cmd[7]), math.Float32frombits(cmd[8]))
+	return
+}
diff --git a/op/clip/clip.go b/op/clip/clip.go
index 0d79975e..282142b5 100644
--- a/op/clip/clip.go
+++ b/op/clip/clip.go
@@ -300,76 +300,12 @@ func (p *Path) Cube(ctrl0, ctrl1, to f32.Point) {
 	ctrl0 = ctrl0.Add(p.pen)
 	ctrl1 = ctrl1.Add(p.pen)
 	to = to.Add(p.pen)
-	// Set the maximum distance proportionally to the longest side
-	// of the bounding rectangle.
-	hull := f32.Rectangle{
-		Min: p.pen,
-		Max: ctrl0,
-	}.Canon().Add(ctrl1).Add(to)
-	l := hull.Dx()
-	if h := hull.Dy(); h > l {
-		l = h
-	}
-	p.approxCubeTo(0, l*0.001, ctrl0, ctrl1, to)
-}
-
-// approxCube approximates a cubic Bézier by a series of quadratic
-// curves.
-func (p *Path) approxCubeTo(splits int, maxDist float32, ctrl0, ctrl1, to f32.Point) int {
-	// The idea is from
-	// https://caffeineowl.com/graphics/2d/vectorial/cubic2quad01.html
-	// where a quadratic approximates a cubic by eliminating its t³ term
-	// from its polynomial expression anchored at the starting point:
-	//
-	// P(t) = pen + 3t(ctrl0 - pen) + 3t²(ctrl1 - 2ctrl0 + pen) + t³(to - 3ctrl1 + 3ctrl0 - pen)
-	//
-	// The control point for the new quadratic Q1 that shares starting point, pen, with P is
-	//
-	// C1 = (3ctrl0 - pen)/2
-	//
-	// The reverse cubic anchored at the end point has the polynomial
-	//
-	// P'(t) = to + 3t(ctrl1 - to) + 3t²(ctrl0 - 2ctrl1 + to) + t³(pen - 3ctrl0 + 3ctrl1 - to)
-	//
-	// The corresponding quadratic Q2 that shares the end point, to, with P has control
-	// point
-	//
-	// C2 = (3ctrl1 - to)/2
-	//
-	// The combined quadratic Bézier, Q, shares both start and end points with its cubic
-	// and use the midpoint between the two curves Q1 and Q2 as control point:
-	//
-	// C = (3ctrl0 - pen + 3ctrl1 - to)/4
-	c := ctrl0.Mul(3).Sub(p.pen).Add(ctrl1.Mul(3)).Sub(to).Mul(1.0 / 4.0)
-	const maxSplits = 32
-	if splits >= maxSplits {
-		p.QuadTo(c, to)
-		return splits
-	}
-	// The maximum distance between the cubic P and its approximation Q given t
-	// can be shown to be
-	//
-	// d = sqrt(3)/36*|to - 3ctrl1 + 3ctrl0 - pen|
-	//
-	// To save a square root, compare d² with the squared tolerance.
-	v := to.Sub(ctrl1.Mul(3)).Add(ctrl0.Mul(3)).Sub(p.pen)
-	d2 := (v.X*v.X + v.Y*v.Y) * 3 / (36 * 36)
-	if d2 <= maxDist*maxDist {
-		p.QuadTo(c, to)
-		return splits
-	}
-	// De Casteljau split the curve and approximate the halves.
-	t := float32(0.5)
-	c0 := p.pen.Add(ctrl0.Sub(p.pen).Mul(t))
-	c1 := ctrl0.Add(ctrl1.Sub(ctrl0).Mul(t))
-	c2 := ctrl1.Add(to.Sub(ctrl1).Mul(t))
-	c01 := c0.Add(c1.Sub(c0).Mul(t))
-	c12 := c1.Add(c2.Sub(c1).Mul(t))
-	c0112 := c01.Add(c12.Sub(c01).Mul(t))
-	splits++
-	splits = p.approxCubeTo(splits, maxDist, c0, c01, c0112)
-	splits = p.approxCubeTo(splits, maxDist, c12, c2, to)
-	return splits
+	data := p.ops.Write(scene.CommandSize + 4)
+	bo := binary.LittleEndian
+	bo.PutUint32(data[0:], uint32(p.contour))
+	ops.EncodeCommand(data[4:], scene.Cubic(p.pen, ctrl0, ctrl1, to, false))
+	p.pen = to
+	p.hasSegments = true
 }
 
 // Close closes the path contour.