From 8b3c7456f0b70257cacd8eb92e1616d379b70983 Mon Sep 17 00:00:00 2001
From: Tim <6703801+timzifer@users.noreply.github.com>
Date: Tue, 13 Oct 2026 12:22:00 +0200
Subject: [PATCH 1/2] A pixel gets the same value in whatever region it is
 drawn
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

A scene drawn in bands or tiles was not, byte for byte, the scene drawn
at once: one or two levels in a few antialiased pixels, depending on
where the clip cut. Of the first pages of 965 files of the pdf.js corpus
drawn by cera at 72 dpi, 85 differed between a render at once and one in
bands of 100 rows or tiles of 256 pixels; none does now. Four causes,
independent of one another:

Edges. AddLine cut a segment at the clip and rounded the ends of what
was left to 1/256 pixel; the slope came from the rounded ends, so the
same segment had another slope in every clip that cut it elsewhere. An
edge is now the line of its segment in device space (c + y·dxdy,
evaluated at device y and moved to the clip in integers) with the rows
of it that are drawn. Nothing is cut: an edge that leaves the clip at a
side is drawn by rasterWide, which sums what lies left of the clip into
its first cell and drops what lies right of it, cell by cell as a wider
clip would have it. Edges that are not nearly level (|dx/dy| <= 256)
are "plain": their x needs no bounding and they cross cells by the
arithmetic there was. Nearly level ones are bounded by the ends of
their segment and cross cells where their line does, since such an end
may have been moved by a caller that shortened the segment outside the
clip (the stroker's chords for culled curves).

Runs. SolidBlitter.BlitRun with a coverage below 255 scaled the color
and composited it, two roundings where BlitCoverage has one. How a row
is cut into runs and coverage depends on the clip. An opaque color is
now blended by runLerp, the arithmetic of covOpaque (scalar, AVX2 and
AVX-512 kernels).

Analytic strokes and culled curves. Which rows of a segment are filled
analytically is set by the vertices of the whole polyline; a curve
outside the clip that takes its chord has other vertices, and rows of
its neighbour inside the clip changed sides. Curves are not culled on
the analytic path any more (subpaths still are).

Analytic rows. They stepped n·c by nx from the first pixel drawn, which
is the clip's; every pixel is now measured from its own place, by an
exact count of x.

Tests: TestFillTiles (paths in rows, columns, tiles and single pixels
against the fill at once), TestCanvasTiles (scenes of fills, strokes,
clips, colors and shaders), TestRunLikeCoverage.

Speed, old against new, best of six alternating runs: the scenes of
BenchmarkScenes that are all edges cost 4 % more (uncached glyphs 5 to
9 %, short strokes and hatches 3 to 6 %), BenchmarkStrokeFixedCost 5 %,
wide fills nothing. The first pages of the pdf.js corpus drawn by cera
at 144 dpi take 1451 ms before and 1426 ms after: no difference that
can be measured.

Fixes #36
---
 README.md       |   7 +
 blit.go         |  20 +-
 doc.go          |   5 +
 raster.go       | 559 ++++++++++++++++++++++++++++++++++++++++--------
 simd_amd64.go   |  41 ++++
 simd_generic.go |   1 +
 strip.go        |  81 ++++---
 stroke.go       |   7 +-
 tiles_test.go   | 259 ++++++++++++++++++++++
 9 files changed, 853 insertions(+), 127 deletions(-)
 create mode 100644 tiles_test.go

diff --git a/README.md b/README.md
index 29097c7..49e764e 100644
--- a/README.md
+++ b/README.md
@@ -131,6 +131,13 @@ stroked text painted with a shader.
   reach it. They are made by the first band that draws the shape, so the
   work stays in the workers. Byte-identical to `Fill`/`Stroke`; 2 000
   hatch lines drawn as 64 bands on one core: 87 → 54 ms (one band: 46 ms).
+- **Bands and tiles**: a pixel gets the same value in whatever region it
+  is drawn (`Canvas.Reset`). An edge is the line of its segment in device
+  space, not the segment cut at the clip and rounded there; a run blends
+  as coverage does; an analytic row measures every pixel from its own
+  place. A scene drawn in tiles is byte-identical to the scene drawn at
+  once (`TestCanvasTiles`), so neither the number of workers nor the
+  tiling of a viewer shows in a page.
 - **Rectangle clips** (91 % of real PDF clips) cost nothing per pixel;
   other clips are rasterized once into a mask that is never cleared as a
   whole and keeps fully opaque runs as runs.
diff --git a/blit.go b/blit.go
index a8ef7a7..79e619f 100644
--- a/blit.go
+++ b/blit.go
@@ -157,8 +157,14 @@ func (b *SolidBlitter) SetColor(c color.RGBA) {
 
 func (b *SolidBlitter) BlitRun(y, x0, x1 int, alpha uint8) {
 	d := b.t.row(y, x0, x1)
-	if alpha == 255 && b.opaque {
-		fill32(d, b.c)
+	if b.opaque {
+		if alpha == 255 {
+			fill32(d, b.c)
+		} else if alpha != 0 {
+			// As BlitCoverage blends it: how a row is cut into runs and
+			// coverage depends on the clip, the pixels must not.
+			runLerp(d, b.c, uint32(alpha))
+		}
 		return
 	}
 	s := b.c
@@ -193,6 +199,16 @@ func runOverScalar(d []uint32, s uint32) {
 	}
 }
 
+// runLerpScalar blends the opaque color c into every pixel with coverage
+// a: the pixels of covOpaqueScalar for a constant coverage.
+func runLerpScalar(d []uint32, c, a uint32) {
+	ca := expand(c) * uint64(a)
+	inv := uint64(255 - a)
+	for i, v := range d {
+		d[i] = compact(div255x(expand(v)*inv + ca))
+	}
+}
+
 // covOpaqueScalar blends an opaque color c (cx expanded) with coverage.
 func covOpaqueScalar(d []uint32, cov []uint8, c uint32, cx uint64) {
 	d = d[:len(cov)]
diff --git a/doc.go b/doc.go
index 47b6550..804fac4 100644
--- a/doc.go
+++ b/doc.go
@@ -54,6 +54,11 @@
 // concurrent use; use one per worker and give each worker its own band or
 // tile of the destination (Canvas.Reset takes a region).
 //
+// A pixel gets the same value in whatever region it is drawn: a scene
+// drawn in bands or tiles is, byte for byte, the scene drawn at once. How
+// many workers draw a page, and how a viewer cuts it into tiles, does not
+// show in it.
+//
 // # Robustness
 //
 // NaN and infinite coordinates are rejected, huge ones are clipped
diff --git a/raster.go b/raster.go
index 5753753..a5e7416 100644
--- a/raster.go
+++ b/raster.go
@@ -52,23 +52,90 @@ const (
 	maxSegs    = 1024
 )
 
-// edge is a line segment in fixed point 24.8, x relative to the clip's left
-// edge and y relative to its top edge, oriented so that y0 < y1.
+// edge is a line segment of a path, oriented downward, as the line it
+// lies on and the rows of it that are drawn. Coordinates are fixed point
+// 24.8, x relative to the clip's left edge and y to its top edge.
+//
+// The line is that of the segment as the path gives it, not of its end
+// points rounded to fixed point, and it is measured in device space, not
+// from the clip: where it passes a row depends neither on where the
+// segment ends nor on the clip. A segment the clip cuts, or one a caller
+// shortened outside the clip, is thereby drawn like the whole one, and a
+// path covers a pixel the same whether it is drawn at once, in bands or
+// in tiles.
 type edge struct {
-	x0, y0, x1, y1 int32
-	dir            int32   // +1 downward in the original path, -1 upward
-	dxdy           float64 // (x1-x0)/(y1-y0)
-	dydx           float64 // |(y1-y0)/(x1-x0)|; 0 until first needed (vertical edges never need it)
+	// In device fixed point, x at y is c + y·dxdy - xAbs, rounded down
+	// (see xAt): c holds the half that rounds and the bias that makes
+	// the conversion a floor.
+	c    float64
+	dxdy float64
+	dydx float64 // 1/dxdy; 0 until first needed (vertical edges never need it)
+	// For a plain edge, lo and hi are its x at y0 and at y1. For any
+	// other they bound x: the x of the segment's ends, xSlack beyond
+	// them. They are equal for an edge whose x is the same in all its
+	// rows, a vertical one above all: it is drawn as that.
+	lo, hi int32
+	y0, y1 int32 // the rows drawn, y0 < y1
+	dir    int32 // +1 downward in the original path, -1 upward
+	// wide says that the edge leaves the clip at its left or right within
+	// its rows: its x is then not a cell of the accumulator everywhere
+	// (see rasterWide). The x of any other edge is in [0, w·256].
+	wide bool
+	// plain says that the edge is not nearly level (see plainSlope): its
+	// line stays near the segment in all its rows, whatever the clip, and
+	// its x needs no bounding (xPlain).
+	plain bool
 }
 
-// xAt returns the edge's x at fixed-point y (y0 <= y <= y1). The same formula
-// is used on both sides of every band and row boundary, so adjacent pieces of
-// an edge always meet exactly.
-func (e *edge) xAt(y int32) int32 {
-	if y >= e.y1 {
-		return e.x1
-	}
-	return e.x0 + int32(float64(y-e.y0)*e.dxdy+0.5+float64(1<<20)) - (1 << 20)
+const (
+	// xFar bounds the x of an edge, in fixed point from the clip's left.
+	xFar = 1 << 29
+	// xAbs bounds the device x of a line, in fixed point, where it is
+	// rounded; what lies beyond is at an end of its edge.
+	xAbs = 1 << 36
+	// xSlack is how far beyond the ends of its segment the x of an edge
+	// may lie, in fixed point. The rows of an edge start and end on
+	// 1/256 pixel, up to half of that beside the segment's ends, where
+	// the line of a nearly level segment is far beyond them.
+	xSlack = 256
+	// plainSlope is the dx/dy up to which an edge is plain: half a step
+	// of y moves its line by half of xSlack at most, so it never leaves
+	// the slack. Only edges within a quarter of a degree of level are
+	// steeper in x.
+	plainSlope = 256
+	// nearCoord bounds the device coordinates up to which a line is
+	// measured from the segment's own point, alike for every clip.
+	nearCoord = 1 << 26
+)
+
+// xAt returns the x at fixed-point y (y0 <= y <= y1) of an edge that is
+// not plain, both relative to the clip, whose top is oy in device fixed
+// point and whose left, with xAbs added, bo. The line is evaluated in
+// device space and the result moved to the clip in integers: it is the
+// same pixel position for every clip. The same formula is used on both
+// sides of every band and row boundary, so adjacent pieces of an edge
+// always meet exactly.
+func (e *edge) xAt(y int32, bo, oy int64) int32 {
+	v := e.c + float64(int64(y)+oy)*e.dxdy
+	if !(v > 0) { // also NaN
+		return e.lo
+	}
+	if v > 2*xAbs {
+		return e.hi
+	}
+	x := int64(v) - bo
+	if x < int64(e.lo) {
+		return e.lo
+	}
+	if x > int64(e.hi) {
+		return e.hi
+	}
+	return int32(x)
+}
+
+// xPlain is xAt for a plain edge: its line as it is.
+func (e *edge) xPlain(y int32, bo, oy int64) int32 {
+	return int32(int64(e.c+float64(int64(y)+oy)*e.dxdy) - bo)
 }
 
 // Rasterizer converts paths into coverage spans. Its cost per path is the
@@ -86,6 +153,9 @@ type Rasterizer struct {
 	cy0    float64
 	cx1    float64
 	cy1    float64
+	ox, oy int64 // the clip's corner in device fixed point
+	bo, yo int64 // ox and oy with xAbs added
+	near   bool  // the clip lies within nearCoord: its rows convert without checks
 	w, h   int
 	stride int // cells per accumulation row: w + 2
 	nw     int // dirty words per row
@@ -140,6 +210,9 @@ func (r *Rasterizer) SetClip(clip image.Rectangle) {
 	r.cx0, r.cy0 = float64(clip.Min.X), float64(clip.Min.Y)
 	r.cx1, r.cy1 = float64(clip.Max.X), float64(clip.Max.Y)
 	r.ry0, r.ry1 = r.cy0, r.cy1
+	r.ox, r.oy = int64(clip.Min.X)<<8, int64(clip.Min.Y)<<8
+	r.bo, r.yo = r.ox+xAbs, r.oy+xAbs
+	r.near = clip.Min.X > -nearCoord && clip.Min.Y > -nearCoord && clip.Max.X < nearCoord && clip.Max.Y < nearCoord
 	r.w, r.h = clip.Dx(), clip.Dy()
 	r.stride = r.w + 2
 	r.nw = (r.stride>>blkShift)/64 + 1
@@ -437,10 +510,10 @@ func segCount(v float64) int {
 	return int(math.Ceil(math.Sqrt(v)))
 }
 
-// AddLine adds a device-space line segment. Segments are clipped
-// analytically: parts above or below the clip are dropped, parts right of it
-// are dropped, and parts left of it become vertical edges on the clip's left
-// border, which preserves the winding of everything inside.
+// AddLine adds a device-space line segment. Of a segment the rows inside
+// the clip are kept; where it lies right of the clip it adds nothing, and
+// where it lies left of it it counts as an edge on the clip's left border,
+// which preserves the winding of everything inside.
 func (r *Rasterizer) AddLine(x0, y0, x1, y1 float64) {
 	r.begin()
 	if y0 == y1 || !(math.Abs(x0) <= math.MaxFloat64 && math.Abs(x1) <= math.MaxFloat64 &&
@@ -456,41 +529,79 @@ func (r *Rasterizer) AddLine(x0, y0, x1, y1 float64) {
 	if y1 <= cy0 || y0 >= cy1 {
 		return
 	}
-	if y0 < cy0 {
-		x0 = lineAt(cy0, y0, y1, x0, x1)
-		y0 = cy0
+	fy0 := r.fixedY(max(y0, cy0))
+	fy1 := r.fixedY(min(y1, cy1))
+	if fy0 >= fy1 {
+		return
 	}
-	if y1 > cy1 {
-		x1 = lineAt(cy1, y0, y1, x0, x1)
-		y1 = cy1
+	right := int32(r.w << 8)
+	var e edge
+	var xa, xb int32
+	if x0 == x1 {
+		x := r.fixedFar(x0)
+		if x >= right {
+			return
+		}
+		x = max(x, 0)
+		e = edge{lo: x, hi: x, y0: fy0, y1: fy1, dir: dir}
+		xa, xb = x, x
+	} else {
+		dy := y1 - y0
+		dxdy := (x1 - x0) / dy
+		if dy > math.MaxFloat64 {
+			dxdy = (x1*0.5 - x0*0.5) / (y1*0.5 - y0*0.5)
+		}
+		// The line is measured from the segment's upper end: from the
+		// same point for every clip. A segment that comes from very far
+		// away is measured from where it enters the clip, which keeps
+		// the numbers small.
+		ax, ay := x0, y0
+		if y0 < cy0 && !(math.Abs(x0) < nearCoord && math.Abs(y0) < nearCoord && math.Abs(x1) < nearCoord && math.Abs(y1) < nearCoord) {
+			ax, ay = lineAt(cy0, y0, y1, x0, x1), cy0
+		}
+		e = edge{c: (ax-ay*dxdy)*256 + (0.5 + xAbs), dxdy: dxdy, y0: fy0, y1: fy1, dir: dir}
+		// Plain by its slope, not by this clip: an edge is drawn by the
+		// same arithmetic in every clip (but one it is too far from to
+		// be counted in fixed point: a conversion out of range gives a
+		// number out of range here, whatever it is).
+		a := int64(e.c+float64(int64(fy0)+r.oy)*dxdy) - r.bo
+		b := int64(e.c+float64(int64(fy1)+r.oy)*dxdy) - r.bo
+		if math.Abs(dxdy) <= plainSlope && uint64(a+xFar) < 2*xFar && uint64(b+xFar) < 2*xFar {
+			xa, xb = int32(a), int32(b)
+			e.lo, e.hi, e.plain = xa, xb, true
+		} else {
+			e.lo = r.fixedFar(min(x0, x1)) - xSlack
+			e.hi = r.fixedFar(max(x0, x1)) + xSlack
+			xa, xb = e.xAt(fy0, r.bo, r.oy), e.xAt(fy1, r.bo, r.oy)
+		}
+		switch {
+		case xa >= right && xb >= right:
+			return
+		case xa <= 0 && xb <= 0:
+			e = edge{y0: fy0, y1: fy1, dir: dir}
+			xa, xb = 0, 0
+		case xa < 0 || xb < 0 || xa > right || xb > right:
+			e.wide = true
+			xa, xb = min(max(xa, 0), right), min(max(xb, 0), right)
+		}
 	}
-	cx0, cx1 := r.cx0, r.cx1
-	if x0 >= cx1 && x1 >= cx1 {
-		return
+	limit := r.MaxEdges
+	if limit == 0 {
+		limit = DefaultMaxEdges
 	}
-	if x0 <= cx0 && x1 <= cx0 {
-		r.emit(cx0, y0, cx0, y1, dir)
+	if len(r.edges) >= limit {
+		r.truncated = true
 		return
 	}
-	if x0 < cx0 || x1 < cx0 {
-		yc := lineAt(cx0, x0, x1, y0, y1)
-		if x0 < cx0 {
-			r.emit(cx0, y0, cx0, yc, dir)
-			x0, y0 = cx0, yc
-		} else {
-			r.emit(cx0, yc, cx0, y1, dir)
-			x1, y1 = cx0, yc
-		}
+	r.edges = append(r.edges, e)
+	r.minX = min(r.minX, xa, xb)
+	r.maxX = max(r.maxX, xa, xb)
+	if fy0 < r.minY {
+		r.minY = fy0
 	}
-	if x0 > cx1 || x1 > cx1 {
-		yc := lineAt(cx1, x0, x1, y0, y1)
-		if x0 > cx1 {
-			x0, y0 = cx1, yc
-		} else {
-			x1, y1 = cx1, yc
-		}
+	if fy1 > r.maxY {
+		r.maxY = fy1
 	}
-	r.emit(x0, y0, x1, y1, dir)
 }
 
 // lineAt returns the value at v of the line through (a, fa) and (b, fb),
@@ -518,37 +629,35 @@ func fixed(v, o float64, n int) int32 {
 	return int32(f)
 }
 
-// emit stores a clipped segment (y0 <= y1) as a fixed-point edge.
-func (r *Rasterizer) emit(x0, y0, x1, y1 float64, dir int32) {
-	fy0 := fixed(y0, r.cy0, r.h)
-	fy1 := fixed(y1, r.cy0, r.h)
-	if fy0 >= fy1 {
-		return
-	}
-	limit := r.MaxEdges
-	if limit == 0 {
-		limit = DefaultMaxEdges
-	}
-	if len(r.edges) >= limit {
-		r.truncated = true
-		return
-	}
-	fx0 := fixed(x0, r.cx0, r.w)
-	fx1 := fixed(x1, r.cx0, r.w)
-	// dydx is computed when a row of the edge first crosses a cell border
-	// (rasterEdge); most edges of small shapes never do.
-	r.edges = append(r.edges, edge{
-		x0: fx0, y0: fy0, x1: fx1, y1: fy1, dir: dir,
-		dxdy: float64(fx1-fx0) / float64(fy1-fy0),
-	})
-	r.minX = min(r.minX, fx0, fx1)
-	r.maxX = max(r.maxX, fx0, fx1)
-	if fy0 < r.minY {
-		r.minY = fy0
+// fixedY converts a device y inside the clip's rows into 24.8 fixed point
+// relative to the clip's top. It is rounded in device space and moved to
+// the clip then: the same y lands on the same place for every clip.
+func (r *Rasterizer) fixedY(v float64) int32 {
+	if r.near {
+		// v is a y of the clip: the sum is positive, and its floor lies
+		// in the clip's rows.
+		return int32(int64(v*256+(0.5+xAbs)) - r.yo)
+	}
+	var y int64
+	if f := v*256 + (0.5 + xAbs); f > 0 && f < 2*xAbs {
+		y = int64(f) - xAbs - r.oy
+	} else if g := math.Floor(v*256+0.5) - float64(r.oy); g > 0 { // not NaN
+		y = int64(min(g, float64(r.h<<8)))
+	}
+	return int32(min(max(y, 0), int64(r.h<<8)))
+}
+
+// fixedFar converts a device x into 24.8 fixed point relative to the
+// clip's left, clamped to ±xFar: an x that may lie outside the clip.
+func (r *Rasterizer) fixedFar(v float64) int32 {
+	f := v*256 + (0.5 + xAbs)
+	if !(f > 0) { // also NaN
+		return -xFar
 	}
-	if fy1 > r.maxY {
-		r.maxY = fy1
+	if f > 2*xAbs {
+		return xFar
 	}
+	return int32(min(max(int64(f)-r.bo, -xFar), xFar))
 }
 
 // reverse reverses the orientation of the edges added since there were n.
@@ -804,23 +913,27 @@ func (r *Rasterizer) rasterStrip(s *accStrip, bandRow int) uint64 {
 				a, b = ia, ib
 			}
 		}
+		// Each pixel from its own place (see rowCtx.centre).
+		x, g := rc.rowAt(i0, fy)
 		var prev int64
 		for i := 0; i < a; i++ {
+			d := nx*x + g
+			x++
 			v := int64((tr.area(k1-d) - tr.area(k2-d)) * fullArea)
 			acc[i] += v - prev
 			prev = v
-			d += nx
 		}
 		if b > a {
 			acc[a] += fullArea - prev
 			prev = fullArea
-			d += float64(b-a) * nx
 		}
+		x, _ = rc.rowAt(i0+b, fy)
 		for i := b; i < n; i++ {
+			d := nx*x + g
+			x++
 			v := int64((tr.area(k1-d) - tr.area(k2-d)) * fullArea)
 			acc[i] += v - prev
 			prev = v
-			d += nx
 		}
 		acc[n] -= prev
 		if !r.narrow {
@@ -847,6 +960,9 @@ func (r *Rasterizer) rasterEdge(e *edge, bandRow int) uint64 {
 	if yA >= yB {
 		return 0
 	}
+	if e.wide {
+		return r.rasterWide(e, bandRow, yA, yB)
+	}
 	stride := r.stride
 	acc := r.acc
 	dirty := r.dirty
@@ -854,10 +970,10 @@ func (r *Rasterizer) rasterEdge(e *edge, bandRow int) uint64 {
 	dir := e.dir
 	var mask uint64
 
-	if e.x0 == e.x1 {
-		// Vertical edge: one cell per row.
-		c := int(e.x0 >> 8)
-		fx2 := (e.x0 & 255) * 2
+	if e.lo == e.hi {
+		// Vertical edge, or one that stays at its x: one cell per row.
+		c := int(e.lo >> 8)
+		fx2 := (e.lo & 255) * 2
 		wA := dir * (512 - fx2)
 		wB := dir * fx2
 		bw, bb := c>>(blkShift+6), uint64(1)<<(uint(c>>blkShift)&63)
@@ -889,7 +1005,19 @@ func (r *Rasterizer) rasterEdge(e *edge, bandRow int) uint64 {
 		return mask
 	}
 
-	x := e.xAt(yA)
+	bo, oy := r.bo, r.oy
+	// The ends of a plain edge are known; most edges of small shapes
+	// have no rows between them.
+	plain := e.plain
+	var x int32
+	switch {
+	case !plain:
+		x = e.xAt(yA, bo, oy)
+	case yA == e.y0:
+		x = e.lo
+	default:
+		x = e.xPlain(yA, bo, oy)
+	}
 	// Indices below are in range by construction: x ∈ [0, w·256], so a cell
 	// and its right neighbour are < stride, and ri < bandH. The single-cell
 	// case, by far the most common, is written without bounds checks.
@@ -900,7 +1028,15 @@ func (r *Rasterizer) rasterEdge(e *edge, bandRow int) uint64 {
 		if rowEnd > yB {
 			rowEnd = yB
 		}
-		xn := e.xAt(rowEnd)
+		var xn int32
+		switch {
+		case !plain:
+			xn = e.xAt(rowEnd, bo, oy)
+		case rowEnd == e.y1:
+			xn = e.hi
+		default:
+			xn = e.xPlain(rowEnd, bo, oy)
+		}
 		ri := int(y>>8) - bandRow
 		c0, c1 := x>>8, xn>>8
 		var lo, hi int32
@@ -925,9 +1061,87 @@ func (r *Rasterizer) rasterEdge(e *edge, bandRow int) uint64 {
 				dr = dirty[ri*nw : (ri+1)*nw]
 			}
 			if e.dydx == 0 {
-				e.dydx = math.Abs(float64(e.y1-e.y0) / float64(e.x1-e.x0))
+				e.dydx = 1 / e.dxdy
 			}
-			lo, hi = cellLine(acc[ri*stride:(ri+1)*stride], dr, x, y-base, xn, rowEnd-base, dir, e.dydx)
+			if plain {
+				lo, hi = cellLine(acc[ri*stride:(ri+1)*stride], dr, x, y-base, xn, rowEnd-base, dir, math.Abs(e.dydx))
+			} else {
+				l := e.at(x, y-base, base, r.ox, oy)
+				lo, hi = cellLineFlat(acc[ri*stride:(ri+1)*stride], dr, x, y-base, xn, rowEnd-base, dir, &l)
+			}
+		}
+		if lo < r.rmin[ri] {
+			r.rmin[ri] = lo
+		}
+		if hi > r.rmax[ri] {
+			r.rmax[ri] = hi
+		}
+		mask |= 1 << uint(ri)
+		x, y = xn, rowEnd
+	}
+	return mask
+}
+
+// rasterWide is rasterEdge for an edge that leaves the clip at its left
+// or right, rows yA to yB of it. It accumulates what the edge would in a
+// clip wide enough for it and keeps the part of the clip's cells: what
+// lies left of the clip is summed into its first cell, as an edge on the
+// border would be, and what lies right of it is dropped. The cells
+// inside are thereby the same for every clip.
+func (r *Rasterizer) rasterWide(e *edge, bandRow int, yA, yB int32) uint64 {
+	stride := r.stride
+	nw := r.nw
+	dir := e.dir
+	last := int32(r.w) // the cell after the clip's; it and its neighbour exist
+	var mask uint64
+	if e.dydx == 0 && e.dxdy != 0 {
+		e.dydx = 1 / e.dxdy
+	}
+	// A plain edge has the x of a plain one here too: in a clip that
+	// holds it, it is drawn by rasterEdge.
+	bo, oy := r.bo, r.oy
+	plain := e.plain
+	at := func(y int32) int32 {
+		if plain {
+			return e.xPlain(y, bo, oy)
+		}
+		return e.xAt(y, bo, oy)
+	}
+	x := at(yA)
+	for y := yA; y < yB; {
+		rowEnd := (y | 255) + 1
+		if rowEnd > yB {
+			rowEnd = yB
+		}
+		xn := at(rowEnd)
+		ri := int(y>>8) - bandRow
+		c0, c1 := x>>8, xn>>8
+		if c0 > last && c1 > last {
+			x, y = xn, rowEnd
+			continue
+		}
+		acc := r.acc[ri*stride : (ri+1)*stride]
+		var dr []uint64
+		if !r.narrow {
+			dr = r.dirty[ri*nw : (ri+1)*nw]
+		}
+		var lo, hi int32
+		switch {
+		case c0 < 0 && c1 < 0:
+			acc[0] += int64((rowEnd - y) * dir * 512)
+			markRange(dr, 0, 1)
+			lo, hi = 0, 1
+		case c0 == c1:
+			d := (rowEnd - y) * dir
+			sx := (x & 255) + (xn & 255)
+			acc[c0] += int64(d * (512 - sx))
+			acc[c0+1] += int64(d * sx)
+			markRange(dr, int(c0), int(c0)+1)
+			lo, hi = c0, c0+1
+		default:
+			base := y &^ 255
+			l := e.at(x, y-base, base, r.ox, oy)
+			lo, hi = cellLineWide(acc, dr, x, y-base, xn, rowEnd-base, dir, &l, last)
 		}
 		if lo < r.rmin[ri] {
 			r.rmin[ri] = lo
@@ -941,9 +1155,57 @@ func (r *Rasterizer) rasterEdge(e *edge, bandRow int) uint64 {
 	return mask
 }
 
-// cellLine accumulates a segment that lies within one pixel row, from
-// (x0, fy0) to (x1, fy1), with 0 <= fy0 < fy1 <= 256, and returns the range
-// of touched cells [lo, hi].
+// rowLine says where an edge crosses the borders of the cells in a pixel
+// row.
+//
+// A plain edge crosses them as the straight piece from where it enters
+// the row, (x0, fy0), with its slope: both are the same for every clip.
+//
+// Any other edge may enter the row at an end of its bounds, off its line,
+// and that end depends on where a caller cut the segment outside the
+// clip. Its crossings are taken from the line itself: in the row that
+// starts at fixed-point y base of a clip whose corner is (ox, oy), it
+// reaches the clip's x at (x + ox - c)·dydx - base, counted from the
+// row's top, with c and base in device fixed point.
+type rowLine struct {
+	plain         bool
+	x0, fy0       int32
+	dydx          float64 // plain: |dy/dx|
+	c, base, offx float64 // base less the half that rounds
+}
+
+// at returns the crossings of e, whose dydx is set, in the row at base,
+// which it enters at (x0, fy0).
+func (e *edge) at(x0, fy0, base int32, ox, oy int64) rowLine {
+	if e.plain {
+		return rowLine{plain: true, x0: x0, fy0: fy0, dydx: math.Abs(e.dydx)}
+	}
+	return rowLine{dydx: e.dydx, c: e.c - (0.5 + xAbs), base: float64(int64(base)+oy) - 0.5, offx: float64(ox)}
+}
+
+// crossing returns the y in the row at which the edge reaches x, no less
+// than ya and no more than yb.
+func (l *rowLine) crossing(x, ya, yb int32) int32 {
+	if l.plain {
+		d := x - l.x0
+		if d < 0 {
+			d = -d
+		}
+		return max(min(l.fy0+int32(float64(d)*l.dydx+0.5), yb), ya)
+	}
+	t := (float64(x)+l.offx-l.c)*l.dydx - l.base
+	if !(t > float64(ya)) { // also NaN
+		return ya
+	}
+	if t >= float64(yb) {
+		return yb
+	}
+	return int32(t)
+}
+
+// cellLine accumulates a segment of a plain edge that lies within one
+// pixel row, from (x0, fy0) to (x1, fy1), with 0 <= fy0 < fy1 <= 256, and
+// returns the range of touched cells [lo, hi]. dydx is |dy/dx|.
 func cellLine(acc []int64, dirty []uint64, x0, fy0, x1, fy1, dir int32, dydx float64) (lo, hi int32) {
 	c0 := x0 >> 8
 	c1 := x1 >> 8
@@ -991,6 +1253,125 @@ func cellLine(acc []int64, dirty []uint64, x0, fy0, x1, fy1, dir int32, dydx flo
 	return c1, c0 + 1
 }
 
+// cellLineFlat is cellLine for an edge that is not plain: it crosses the
+// borders of the cells where l says.
+func cellLineFlat(acc []int64, dirty []uint64, x0, fy0, x1, fy1, dir int32, l *rowLine) (lo, hi int32) {
+	c0 := x0 >> 8
+	c1 := x1 >> 8
+	if c0 == c1 {
+		d := (fy1 - fy0) * dir
+		s := (x0 & 255) + (x1 & 255)
+		acc[c0] += int64(d * (512 - s))
+		acc[c0+1] += int64(d * s)
+		markRange(dirty, int(c0), int(c0)+1)
+		return c0, c0 + 1
+	}
+	if x1 > x0 {
+		xa, ya := x0, fy0
+		for c := c0; c < c1; c++ {
+			bx := (c + 1) << 8
+			yb := l.crossing(bx, ya, fy1)
+			d := (yb - ya) * dir
+			s := (xa - c<<8) + 256
+			acc[c] += int64(d * (512 - s))
+			acc[c+1] += int64(d * s)
+			xa, ya = bx, yb
+		}
+		d := (fy1 - ya) * dir
+		s := x1 - c1<<8
+		acc[c1] += int64(d * (512 - s))
+		acc[c1+1] += int64(d * s)
+		markRange(dirty, int(c0), int(c1)+1)
+		return c0, c1 + 1
+	}
+	xa, ya := x0, fy0
+	for c := c0; c > c1; c-- {
+		bx := c << 8
+		yb := l.crossing(bx, ya, fy1)
+		d := (yb - ya) * dir
+		s := xa - c<<8
+		acc[c] += int64(d * (512 - s))
+		acc[c+1] += int64(d * s)
+		xa, ya = bx, yb
+	}
+	d := (fy1 - ya) * dir
+	s := 256 + (x1 - c1<<8)
+	acc[c1] += int64(d * (512 - s))
+	acc[c1+1] += int64(d * s)
+	markRange(dirty, int(c1), int(c0)+1)
+	return c1, c0 + 1
+}
+
+// cellLineWide is cellLine for a segment whose cells are not all the
+// clip's, of which last is the one after them: cells left of the first
+// are summed into it, cells after last are left out. Each cell gets what
+// cellLine or cellLineFlat gives it.
+func cellLineWide(acc []int64, dirty []uint64, x0, fy0, x1, fy1, dir int32, l *rowLine, last int32) (lo, hi int32) {
+	c0 := x0 >> 8
+	c1 := x1 >> 8
+	put := func(c int32, a, b int32) {
+		switch {
+		case c < 0:
+			acc[0] += int64(a) + int64(b)
+		case c <= last:
+			acc[c] += int64(a)
+			acc[c+1] += int64(b)
+		}
+	}
+	if x1 > x0 {
+		xa, ya := x0, fy0
+		c := c0
+		if c < -1 {
+			// The cells before the one left of the clip, at once.
+			c = -1
+			xa, ya = -256, l.crossing(-256, fy0, fy1)
+			acc[0] += int64((ya - fy0) * dir * 512)
+		}
+		for ; c < c1 && c <= last; c++ {
+			bx := (c + 1) << 8
+			yb := l.crossing(bx, ya, fy1)
+			d := (yb - ya) * dir
+			s := (xa - c<<8) + 256
+			put(c, d*(512-s), d*s)
+			xa, ya = bx, yb
+		}
+		if c == c1 {
+			d := (fy1 - ya) * dir
+			s := x1 - c1<<8
+			put(c1, d*(512-s), d*s)
+		}
+		lo, hi = max(c0, 0), min(c1, last)+1
+		markRange(dirty, int(lo), int(hi))
+		return lo, hi
+	}
+	xa, ya := x0, fy0
+	c := c0
+	if c > last {
+		// The cells after last add nothing.
+		c = last
+		xa, ya = (last+1)<<8, l.crossing((last+1)<<8, fy0, fy1)
+	}
+	for ; c > c1 && c >= 0; c-- {
+		bx := c << 8
+		yb := l.crossing(bx, ya, fy1)
+		d := (yb - ya) * dir
+		s := xa - c<<8
+		put(c, d*(512-s), d*s)
+		xa, ya = bx, yb
+	}
+	if c == c1 {
+		d := (fy1 - ya) * dir
+		s := 256 + (x1 - c1<<8)
+		put(c1, d*(512-s), d*s)
+	} else {
+		// The rest lies left of the clip.
+		acc[0] += int64((fy1 - ya) * dir * 512)
+	}
+	lo, hi = max(c1, 0), min(c0, last)+1
+	markRange(dirty, int(lo), int(hi))
+	return lo, hi
+}
+
 // markRange sets the dirty bits covering cells [lo, hi].
 func markRange(dirty []uint64, lo, hi int) {
 	if dirty == nil {
diff --git a/simd_amd64.go b/simd_amd64.go
index 545c934..f6d1f6e 100644
--- a/simd_amd64.go
+++ b/simd_amd64.go
@@ -171,6 +171,31 @@ func runOver(d []uint32, s uint32) {
 	runOverScalar(d[i:], s)
 }
 
+func runLerp(d []uint32, c, a uint32) {
+	if !useSIMD || len(d) < 4 {
+		runLerpScalar(d, c, a)
+		return
+	}
+	if use512 && len(d) >= 8 {
+		runLerp512(d, c, a)
+		return
+	}
+	db := bytesOf(d)
+	k128 := archsimd.BroadcastUint16x16(128)
+	k257 := archsimd.BroadcastUint16x16(257)
+	pk := loadInt8x32(&packIdx)
+	ca := archsimd.BroadcastUint32x4(c).AsUint8x16().ExtendToUint16().Mul(archsimd.BroadcastUint16x16(uint16(a)))
+	inv := archsimd.BroadcastUint16x16(uint16(255 - a))
+	i := 0
+	for ; i+4 <= len(d); i += 4 {
+		px := db[i*4 : i*4+16]
+		d16 := loadUint8x16(px).ExtendToUint16()
+		storeUint8x16(packLanes(div255v(d16.Mul(inv).Add(ca), k128, k257), pk), px)
+	}
+	archsimd.ClearAVXUpperBits()
+	runLerpScalar(d[i:], c, a)
+}
+
 // AVX-512 kernels: eight pixels in 16-bit lanes of one ZMM register.
 
 func div255w(x, k128, k257 archsimd.Uint16x32) archsimd.Uint16x32 {
@@ -248,6 +273,22 @@ func runOver512(d []uint32, s uint32) {
 	runOverScalar(d[i:], s)
 }
 
+func runLerp512(d []uint32, c, a uint32) {
+	db := bytesOf(d)
+	k128 := archsimd.BroadcastUint16x32(128)
+	k257 := archsimd.BroadcastUint16x32(257)
+	ca := archsimd.BroadcastUint32x8(c).AsUint8x32().ExtendToUint16().Mul(archsimd.BroadcastUint16x32(uint16(a)))
+	inv := archsimd.BroadcastUint16x32(uint16(255 - a))
+	i := 0
+	for ; i+8 <= len(d); i += 8 {
+		px := db[i*4 : i*4+32]
+		d16 := loadUint8x32(px).ExtendToUint16()
+		storeUint8x32(truncToUint8(div255w(d16.Mul(inv).Add(ca), k128, k257)), px)
+	}
+	archsimd.ClearAVXUpperBits()
+	runLerpScalar(d[i:], c, a)
+}
+
 // spanOver and spanOverCov composite varying source pixels: a shader's
 // span, or a layer's or an image's own row.
 
diff --git a/simd_generic.go b/simd_generic.go
index 273860e..ae27275 100644
--- a/simd_generic.go
+++ b/simd_generic.go
@@ -7,6 +7,7 @@ package stilus
 func SIMD() bool { return false }
 
 func runOver(d []uint32, s uint32)                           { runOverScalar(d, s) }
+func runLerp(d []uint32, c, a uint32)                        { runLerpScalar(d, c, a) }
 func covOpaque(d []uint32, cov []uint8, c uint32, cx uint64) { covOpaqueScalar(d, cov, c, cx) }
 func covOver(d []uint32, cov []uint8, cx uint64)             { covOverScalar(d, cov, cx) }
 func spanOver(d, s []uint32, k, a uint32)                    { spanOverScalar(d, s, k, a) }
diff --git a/strip.go b/strip.go
index c478a12..203c65f 100644
--- a/strip.go
+++ b/strip.go
@@ -491,11 +491,22 @@ func (rc *rowCtx) span(j, cx0, cx1 int) (i0, i1 int, fy float64) {
 	return max(ffloor(rc.xl+fy*rc.sl), cx0), min(ffloor(rc.xr+fy*rc.sl)+1, cx1), fy
 }
 
-// centre returns n·c of the centre of pixel (i0, fy).
+// centre returns n·c of the centre of pixel (i0, fy). Every pixel of a
+// row is measured this way, from its own place, and not by steps of nx
+// from the first one drawn: where a row starts depends on the clip, what
+// a pixel gets must not. Loops over a row count the centres' x instead,
+// which is exact (see rowAt).
 func (rc *rowCtx) centre(i0 int, fy float64) float64 {
 	return rc.nx*(float64(i0)+0.5) + rc.ny*(fy+0.5)
 }
 
+// rowAt returns the x of the centre of pixel i0 and the part of n·c that
+// row fy adds: the centre of the pixel k further has n·c = nx·(x+k) + g,
+// which is centre(i0+k, fy) to the bit.
+func (rc *rowCtx) rowAt(i0 int, fy float64) (x, g float64) {
+	return float64(i0) + 0.5, rc.ny * (fy + 0.5)
+}
+
 // borders returns the clip's partially covered border columns (the first
 // and last ones, or math.MinInt when fully covered): rows touching them
 // are composited through the frac blitter or summed in the accumulator
@@ -522,7 +533,7 @@ func (f *segFast) solidRows(y0, y1 int, rc *rowCtx) {
 	// last; they are taken out of the loop.
 	if y0 == clip.Min.Y && f.frac[1] != 255 {
 		if i0, i1, fy := rc.span(y0, cx0, cx1); i0 < i1 {
-			f.row(y0, i0, i1, rc.centre(i0, fy), rc)
+			f.row(y0, i0, i1, fy, rc)
 		}
 		y0++
 	}
@@ -539,15 +550,16 @@ func (f *segFast) solidRows(y0, y1 int, rc *rowCtx) {
 		if i0 >= i1 {
 			continue
 		}
-		d := rc.centre(i0, fy)
 		if i0 == rc.col0 || i1 == rc.col1 || (SIMD() && i1-i0 >= 16) {
-			f.row(j, i0, i1, d, rc)
+			f.row(j, i0, i1, fy, rc)
 			continue
 		}
 		row := solid.t.row(j, i0, i1)
+		x, g := rc.rowAt(i0, fy)
 		for i, dst := range row {
+			d := nx*x + g
+			x++
 			a := quant(tr.area(k1-d) - tr.area(k2-d))
-			d += nx
 			switch a {
 			case 0:
 			case 255:
@@ -559,7 +571,7 @@ func (f *segFast) solidRows(y0, y1 int, rc *rowCtx) {
 	}
 	if last >= 0 {
 		if i0, i1, fy := rc.span(last, cx0, cx1); i0 < i1 {
-			f.row(last, i0, i1, rc.centre(i0, fy), rc)
+			f.row(last, i0, i1, fy, rc)
 		}
 	}
 }
@@ -576,50 +588,50 @@ func (f *segFast) blitRows(y0, y1 int, rc *rowCtx) {
 		if i0 >= i1 {
 			continue
 		}
-		d0 := rc.centre(i0, fy)
 		if j == f.by0 || j == f.by1 {
-			f.borderRow(j, i0, i1, d0, rc)
+			f.borderRow(j, i0, i1, fy, rc)
 			continue
 		}
 		if (f.bx0 >= i0 && f.bx0 < i1) || (f.bx1 >= i0 && f.bx1 < i1) {
-			f.borderCols(j, i0, i1, d0, rc)
+			f.borderCols(j, i0, i1, fy, rc)
 			continue
 		}
 		if rc.runs {
-			f.emitRow(j, i0, i1, d0, rc)
+			f.emitRow(j, i0, i1, fy, rc)
 			continue
 		}
 		// A thin strip's row: one coverage span (emitRow, inline).
 		cov := f.cov[:i1-i0]
-		d := d0
+		x, g := rc.rowAt(i0, fy)
 		for i := range cov {
+			d := nx*x + g
+			x++
 			cov[i] = quant(tr.area(k1-d) - tr.area(k2-d))
-			d += nx
 		}
 		emitCoverage(f.b, j, i0, cov)
 	}
 }
 
-// row handles pixels [i0, i1) of row j (d0 = n·centre of pixel i0) that
+// row handles pixels [i0, i1) of row j (fy as rowCtx.span gives it) that
 // solidRows cannot take: rows in a partially covered clip border row or
 // crossing a border column go to the accumulator, the rest through the
 // blitter chain.
-func (f *segFast) row(j, i0, i1 int, d0 float64, rc *rowCtx) {
+func (f *segFast) row(j, i0, i1 int, fy float64, rc *rowCtx) {
 	if j == f.by0 || j == f.by1 {
-		f.borderRow(j, i0, i1, d0, rc)
+		f.borderRow(j, i0, i1, fy, rc)
 		return
 	}
 	if (f.bx0 >= i0 && f.bx0 < i1) || (f.bx1 >= i0 && f.bx1 < i1) {
-		f.borderCols(j, i0, i1, d0, rc)
+		f.borderCols(j, i0, i1, fy, rc)
 		return
 	}
-	f.emitRow(j, i0, i1, d0, rc)
+	f.emitRow(j, i0, i1, fy, rc)
 }
 
-// emitRow blits pixels [i0, i1) of row j (d0 = n·centre of pixel i0)
-// through the blitter chain: edge pixels as coverage, a fully covered
-// interior as a run.
-func (f *segFast) emitRow(j, i0, i1 int, d0 float64, rc *rowCtx) {
+// emitRow blits pixels [i0, i1) of row j through the blitter chain: edge
+// pixels as coverage, a fully covered interior as a run.
+func (f *segFast) emitRow(j, i0, i1 int, fy float64, rc *rowCtx) {
+	d0 := rc.centre(i0, fy)
 	n := i1 - i0
 	a, b := n, n // interior [a, b) relative to i0; none by default
 	if rc.runs {
@@ -634,10 +646,10 @@ func (f *segFast) emitRow(j, i0, i1 int, d0 float64, rc *rowCtx) {
 		}
 	}
 	cov := f.cov[:n]
-	d := d0
+	x, g := rc.rowAt(i0, fy)
 	for i := 0; i < a; i++ {
-		cov[i] = quant(rc.cov(d))
-		d += rc.nx
+		cov[i] = quant(rc.cov(rc.nx*x + g))
+		x++
 	}
 	if a > 0 {
 		emitCoverage(f.b, j, i0, cov[:a])
@@ -645,10 +657,10 @@ func (f *segFast) emitRow(j, i0, i1 int, d0 float64, rc *rowCtx) {
 	if b > a {
 		f.b.BlitRun(j, i0+a, i0+b, 255)
 	}
-	d = d0 + float64(b)*rc.nx
+	x, _ = rc.rowAt(i0+b, fy)
 	for i := b; i < n; i++ {
-		cov[i] = quant(rc.cov(d))
-		d += rc.nx
+		cov[i] = quant(rc.cov(rc.nx*x + g))
+		x++
 	}
 	if n > b {
 		emitCoverage(f.b, j, i0+b, cov[b:n])
@@ -658,33 +670,34 @@ func (f *segFast) emitRow(j, i0, i1 int, d0 float64, rc *rowCtx) {
 // borderCols handles a row crossing a partially covered clip border
 // column: the border pixels go to the accumulator, the pieces between them
 // are emitted like any other row.
-func (f *segFast) borderCols(j, i0, i1 int, d0 float64, rc *rowCtx) {
+func (f *segFast) borderCols(j, i0, i1 int, fy float64, rc *rowCtx) {
 	x := i0
 	for _, bx := range [2]int{min(f.bx0, f.bx1), max(f.bx0, f.bx1)} {
 		if bx < x || bx >= i1 {
 			continue
 		}
 		if bx > x {
-			f.emitRow(j, x, bx, d0+float64(x-i0)*rc.nx, rc)
+			f.emitRow(j, x, bx, fy, rc)
 		}
-		f.inject(bx, j, rc.cov(d0+float64(bx-i0)*rc.nx))
+		f.inject(bx, j, rc.cov(rc.centre(bx, fy)))
 		x = bx + 1
 	}
 	if x < i1 {
-		f.emitRow(j, x, i1, d0+float64(x-i0)*rc.nx, rc)
+		f.emitRow(j, x, i1, fy, rc)
 	}
 }
 
 // borderRow sends a row lying in a partially covered clip border row to
 // the accumulator: fully covered stretches as one rectangle each, edge
 // pixels one by one.
-func (f *segFast) borderRow(j, i0, i1 int, d float64, rc *rowCtx) {
+func (f *segFast) borderRow(j, i0, i1 int, fy float64, rc *rowCtx) {
 	// The current fully covered stretch starts at runStart while inRun;
 	// x can be negative, so no coordinate doubles as the marker.
 	inRun, runStart := false, 0
+	xc, g := rc.rowAt(i0, fy)
 	for x := i0; x < i1; x++ {
-		c := rc.cov(d)
-		d += rc.nx
+		c := rc.cov(rc.nx*xc + g)
+		xc++
 		if c >= 1-1e-12 {
 			if !inRun {
 				inRun, runStart = true, x
diff --git a/stroke.go b/stroke.go
index 23e5bdc..4ce6686 100644
--- a/stroke.go
+++ b/stroke.go
@@ -254,8 +254,11 @@ func (s *Stroker) run(p *Path, m Matrix, st *StrokeStyle, pr strokePrep) {
 	}
 
 	// A chord changes the outline only outside the clip, but it shortens
-	// the path a dash pattern walks.
-	s.cullCurves = s.cull && !dashed
+	// the path a dash pattern walks. On the analytic path it also moves
+	// the rows that are filled analytically, which the vertices of the
+	// polyline set, those outside the clip too: the pixels inside would
+	// depend on the clip.
+	s.cullCurves = s.cull && !dashed && s.seg == nil
 	pts := p.Points
 	pi := 0
 	closed := false
diff --git a/tiles_test.go b/tiles_test.go
new file mode 100644
index 0000000..4fc20b2
--- /dev/null
+++ b/tiles_test.go
@@ -0,0 +1,259 @@
+package stilus
+
+import (
+	"bytes"
+	"image"
+	"image/color"
+	"math/rand/v2"
+	"testing"
+)
+
+// randomPath returns a path of lines and curves in and around a box of
+// the given size, some of them far outside it.
+func randomPath(rng *rand.Rand, w, h float64) *Path {
+	f := func(n float64) float32 {
+		v := rng.Float64()*1.6 - 0.3
+		if rng.IntN(12) == 0 {
+			v = rng.Float64()*40 - 20 // far away
+		}
+		return float32(v * n)
+	}
+	var p Path
+	for range 1 + rng.IntN(3) {
+		p.MoveTo(f(w), f(h))
+		for range 1 + rng.IntN(6) {
+			switch rng.IntN(5) {
+			case 0:
+				p.QuadTo(f(w), f(h), f(w), f(h))
+			case 1:
+				p.CubicTo(f(w), f(h), f(w), f(h), f(w), f(h))
+			case 2:
+				// Nearly level, nearly upright.
+				x, y := f(w), f(h)
+				p.LineTo(x, y)
+				if rng.IntN(2) == 0 {
+					p.LineTo(f(w), y+float32(rng.Float64()*0.01))
+				} else {
+					p.LineTo(x+float32(rng.Float64()*0.01), f(h))
+				}
+			default:
+				p.LineTo(f(w), f(h))
+			}
+		}
+		p.Close()
+	}
+	return &p
+}
+
+// tilings are the ways the tests cut a clip: rows, columns, tiles.
+func tilings(r image.Rectangle) [][]image.Rectangle {
+	var out [][]image.Rectangle
+	for _, size := range []image.Point{{r.Dx(), 7}, {r.Dx(), 32}, {r.Dx(), 33}, {5, r.Dy()}, {64, r.Dy()}, {17, 13}, {1, 1}} {
+		var tiles []image.Rectangle
+		for y := r.Min.Y; y < r.Max.Y; y += size.Y {
+			for x := r.Min.X; x < r.Max.X; x += size.X {
+				tiles = append(tiles, image.Rect(x, y, x+size.X, y+size.Y).Intersect(r))
+			}
+		}
+		out = append(out, tiles)
+	}
+	return out
+}
+
+// TestFillTiles holds a path filled in tiles against the path filled at
+// once: every pixel has the same coverage, whatever the clip.
+func TestFillTiles(t *testing.T) {
+	rng := rand.New(rand.NewPCG(11, 12))
+	box := image.Rect(-3, 5, 97, 85)
+	for it := range 300 {
+		p := randomPath(rng, 100, 80)
+		m := []Matrix{Identity, {0.7, 0.4, -0.3, 1.2, 5, -3}, {-1, 0, 0, 1, 100, 0}, {3, 0, 0, 3, -100, -80}}[rng.IntN(4)]
+		rule := FillRule(rng.IntN(2))
+		want := image.NewAlpha(box)
+		r := NewRasterizer(box)
+		r.Fill(p, m, rule, &MaskBlitter{want})
+		for _, tiles := range tilings(box) {
+			if it%10 != 0 && tiles[0].Dx() == 1 {
+				continue // pixel by pixel is slow
+			}
+			got := image.NewAlpha(box)
+			for _, tile := range tiles {
+				r.SetClip(tile)
+				r.Fill(p, m, rule, &MaskBlitter{got})
+			}
+			if mean, max := diff(got, want); max > 0 {
+				t.Fatalf("iteration %d, tiles of %v: mean %.4f max %d", it, tiles[0].Size(), mean, max)
+			}
+		}
+	}
+}
+
+// hashShader paints every pixel a color of its place alone, opaque or
+// not.
+type hashShader struct{ opaque bool }
+
+func (h hashShader) ShadeSpan(y, x int, out []uint32) {
+	for i := range out {
+		v := uint32(x+i)*0x9e3779b1 ^ uint32(y)*0x85ebca6b
+		v ^= v >> 15
+		v *= 0x2c1b3c6d
+		v ^= v >> 12
+		c := color.RGBA{uint8(v), uint8(v >> 8), uint8(v >> 16), 255}
+		if !h.opaque {
+			a := uint32(uint8(v >> 24))
+			c = color.RGBA{uint8(uint32(c.R) * a / 255), uint8(uint32(c.G) * a / 255), uint8(uint32(c.B) * a / 255), uint8(a)}
+		}
+		out[i] = PackRGBA(c)
+	}
+}
+
+// randomPaint returns an opaque or a translucent color, or a shader.
+func randomPaint(rng *rand.Rand) *Paint {
+	switch rng.IntN(5) {
+	case 0:
+		return &Paint{Shader: hashShader{opaque: true}}
+	case 1:
+		return &Paint{Shader: hashShader{}}
+	case 2:
+		a := uint32(1 + rng.IntN(254))
+		ch := func() uint8 { return uint8(uint32(rng.IntN(256)) * a / 255) }
+		return &Paint{Color: color.RGBA{ch(), ch(), ch(), uint8(a)}}
+	}
+	return &Paint{Color: color.RGBA{uint8(rng.IntN(256)), uint8(rng.IntN(256)), uint8(rng.IntN(256)), 255}}
+}
+
+// TestCanvasTiles holds scenes drawn in tiles against the scenes drawn at
+// once: fills and strokes of all kinds, in colors and shaders, under
+// clips, over a backdrop that is not plain. A viewer that draws a page in
+// tiles, or in bands by several workers, gets the pixels of the page.
+func TestCanvasTiles(t *testing.T) {
+	rng := rand.New(rand.NewPCG(21, 22))
+	box := image.Rect(-3, 5, 97, 85)
+	backdrop := image.NewRGBA(box)
+	hashShader{opaque: true}.ShadeSpan(0, 0, pixels(backdrop))
+	for it := range 400 {
+		var ops []func(c *Canvas)
+		m := []Matrix{Identity, {0.7, 0.4, -0.3, 1.2, 5, -3}, {-1, 0, 0, 1, 100, 0}, {3, 0, 0, 3, -100, -80}, {1.5, 0, 0, -1.5, 0, 90}}[rng.IntN(5)]
+		clips := 0
+		for range 1 + rng.IntN(4) {
+			p := randomPath(rng, 100, 80)
+			if rng.IntN(3) == 0 {
+				// A box, as most of what a page draws.
+				x, y := float32(rng.Float64()*60), float32(rng.Float64()*40)
+				w, h := float32(rng.Float64()*70), float32(rng.Float64()*50)
+				p = &Path{}
+				p.MoveTo(x, y)
+				p.LineTo(x+w, y)
+				p.LineTo(x+w, y+h)
+				p.LineTo(x, y+h)
+				p.Close()
+			}
+			st := &StrokeStyle{
+				Width:      []float64{0, 0.3, 0.75, 1, 2.83, 9}[rng.IntN(6)],
+				Cap:        Cap(rng.IntN(3)),
+				Join:       Join(rng.IntN(3)),
+				MiterLimit: []float64{1, 4, 20}[rng.IntN(3)],
+			}
+			if rng.IntN(4) == 0 {
+				st.Dash = [][]float64{{7, 5}, {1, 1}, {0.2, 0.3}, {12, 3, 2, 3}}[rng.IntN(4)]
+				st.DashPhase = rng.Float64() * 5
+			}
+			rule := FillRule(rng.IntN(2))
+			paint := randomPaint(rng)
+			switch k := rng.IntN(10); {
+			case k < 4:
+				ops = append(ops, func(c *Canvas) { c.Fill(p, m, rule, paint) })
+			case k < 8:
+				ops = append(ops, func(c *Canvas) { c.Stroke(p, m, st, paint) })
+			case k == 8:
+				clips++
+				ops = append(ops, func(c *Canvas) { c.ClipPath(p, m, rule) })
+			default:
+				clips++
+				r := Rect{rng.Float64() * 40, rng.Float64() * 30, 40 + rng.Float64()*60, 30 + rng.Float64()*50}
+				ops = append(ops, func(c *Canvas) { c.ClipRect(r, m) })
+			}
+		}
+		draw := func(dst *image.RGBA, region image.Rectangle) {
+			c := NewCanvas(dst)
+			c.Reset(dst, region)
+			for _, op := range ops {
+				op(c)
+			}
+			for range clips {
+				c.PopClip()
+			}
+			if err := c.Err(); err != nil {
+				t.Fatalf("iteration %d: %v", it, err)
+			}
+		}
+		want := image.NewRGBA(box)
+		copy(want.Pix, backdrop.Pix)
+		draw(want, box)
+		for _, tiles := range tilings(box) {
+			if it%10 != 0 && tiles[0].Dx() == 1 {
+				continue
+			}
+			got := image.NewRGBA(box)
+			copy(got.Pix, backdrop.Pix)
+			for _, tile := range tiles {
+				draw(got, tile)
+			}
+			if !bytes.Equal(got.Pix, want.Pix) {
+				n, at := 0, image.Point{}
+				for i := 0; i < len(got.Pix); i += 4 {
+					if !bytes.Equal(got.Pix[i:i+4], want.Pix[i:i+4]) {
+						if n == 0 {
+							at = image.Pt(box.Min.X+(i/4)%box.Dx(), box.Min.Y+(i/4)/box.Dx())
+						}
+						n++
+					}
+				}
+				t.Errorf("iteration %d, tiles of %v: %d pixels differ, the first at %v: %v, want %v", it, tiles[0].Size(), n, at, got.RGBAAt(at.X, at.Y), want.RGBAAt(at.X, at.Y))
+				break
+			}
+		}
+	}
+}
+
+// TestRunLikeCoverage holds a run of constant coverage against the same
+// coverage given pixel by pixel: the rasterizer cuts a row into runs and
+// coverage as the clip and the shape have it, so the two must blend
+// alike.
+func TestRunLikeCoverage(t *testing.T) {
+	rng := rand.New(rand.NewPCG(31, 32))
+	const n = 41
+	rect := image.Rect(0, 0, n, 1)
+	for _, tt := range []struct {
+		name string
+		new  func(dst *image.RGBA) Blitter
+	}{
+		{"opaque color", func(dst *image.RGBA) Blitter { return NewSolidBlitter(dst, color.RGBA{10, 128, 250, 255}) }},
+		{"translucent color", func(dst *image.RGBA) Blitter { return NewSolidBlitter(dst, color.RGBA{10, 100, 90, 130}) }},
+		{"opaque shader", func(dst *image.RGBA) Blitter { return NewShaderBlitter(dst, hashShader{opaque: true}) }},
+		{"translucent shader", func(dst *image.RGBA) Blitter { return NewShaderBlitter(dst, hashShader{}) }},
+	} {
+		for alpha := 0; alpha < 256; alpha++ {
+			run, cov := image.NewRGBA(rect), image.NewRGBA(rect)
+			for i := range run.Pix {
+				run.Pix[i] = uint8(rng.IntN(256))
+			}
+			for i := 3; i < len(run.Pix); i += 4 {
+				run.Pix[i] = 255 // an opaque backdrop, as a page is
+			}
+			copy(cov.Pix, run.Pix)
+			for _, span := range [][2]int{{0, n}, {3, 9}, {7, 8}} {
+				tt.new(run).BlitRun(0, span[0], span[1], uint8(alpha))
+				c := make([]uint8, span[1]-span[0])
+				for i := range c {
+					c[i] = uint8(alpha)
+				}
+				tt.new(cov).BlitCoverage(0, span[0], c)
+			}
+			if !bytes.Equal(run.Pix, cov.Pix) {
+				t.Errorf("%s, coverage %d: a run blends unlike coverage", tt.name, alpha)
+				break
+			}
+		}
+	}
+}

From 658ebcc8ca110243ce5d15b6487168a372c91b8b Mon Sep 17 00:00:00 2001
From: Tim <6703801+timzifer@users.noreply.github.com>
Date: Tue, 13 Oct 2026 12:30:31 +0200
Subject: [PATCH 2/2] Round the products of the line arithmetic before they are
 added

A compiler may fuse a multiplication and an addition into one operation
with one rounding (arm64, amd64 at GOAMD64=v3), and need not do so alike
wherever the same expression stands. The x of an edge is computed in
AddLine, which keeps it at the ends, and again row by row; the cells an
edge writes to are in range because its x lies between the ends. The
conversions make every one of these the same arithmetic on every
machine, and the centres of analytic rows with them.
---
 raster.go | 27 +++++++++++++++++----------
 strip.go  | 15 ++++++++-------
 2 files changed, 25 insertions(+), 17 deletions(-)

diff --git a/raster.go b/raster.go
index a5e7416..db10398 100644
--- a/raster.go
+++ b/raster.go
@@ -116,7 +116,7 @@ const (
 // sides of every band and row boundary, so adjacent pieces of an edge
 // always meet exactly.
 func (e *edge) xAt(y int32, bo, oy int64) int32 {
-	v := e.c + float64(int64(y)+oy)*e.dxdy
+	v := e.c + float64(float64(int64(y)+oy)*e.dxdy)
 	if !(v > 0) { // also NaN
 		return e.lo
 	}
@@ -134,8 +134,15 @@ func (e *edge) xAt(y int32, bo, oy int64) int32 {
 }
 
 // xPlain is xAt for a plain edge: its line as it is.
+//
+// The product is rounded before it is added, here and wherever a line is
+// evaluated (the conversion says so): a compiler may fuse a multiplication
+// and an addition into one operation with one rounding, and need not do
+// so alike everywhere. The x of an edge must be the same wherever it is
+// computed — AddLine keeps what it is at the ends, and the cells an edge
+// writes to are in range because it lies between them.
 func (e *edge) xPlain(y int32, bo, oy int64) int32 {
-	return int32(int64(e.c+float64(int64(y)+oy)*e.dxdy) - bo)
+	return int32(int64(e.c+float64(float64(int64(y)+oy)*e.dxdy)) - bo)
 }
 
 // Rasterizer converts paths into coverage spans. Its cost per path is the
@@ -564,8 +571,8 @@ func (r *Rasterizer) AddLine(x0, y0, x1, y1 float64) {
 		// same arithmetic in every clip (but one it is too far from to
 		// be counted in fixed point: a conversion out of range gives a
 		// number out of range here, whatever it is).
-		a := int64(e.c+float64(int64(fy0)+r.oy)*dxdy) - r.bo
-		b := int64(e.c+float64(int64(fy1)+r.oy)*dxdy) - r.bo
+		a := int64(e.c+float64(float64(int64(fy0)+r.oy)*dxdy)) - r.bo
+		b := int64(e.c+float64(float64(int64(fy1)+r.oy)*dxdy)) - r.bo
 		if math.Abs(dxdy) <= plainSlope && uint64(a+xFar) < 2*xFar && uint64(b+xFar) < 2*xFar {
 			xa, xb = int32(a), int32(b)
 			e.lo, e.hi, e.plain = xa, xb, true
@@ -917,7 +924,7 @@ func (r *Rasterizer) rasterStrip(s *accStrip, bandRow int) uint64 {
 		x, g := rc.rowAt(i0, fy)
 		var prev int64
 		for i := 0; i < a; i++ {
-			d := nx*x + g
+			d := float64(nx*x) + g
 			x++
 			v := int64((tr.area(k1-d) - tr.area(k2-d)) * fullArea)
 			acc[i] += v - prev
@@ -929,7 +936,7 @@ func (r *Rasterizer) rasterStrip(s *accStrip, bandRow int) uint64 {
 		}
 		x, _ = rc.rowAt(i0+b, fy)
 		for i := b; i < n; i++ {
-			d := nx*x + g
+			d := float64(nx*x) + g
 			x++
 			v := int64((tr.area(k1-d) - tr.area(k2-d)) * fullArea)
 			acc[i] += v - prev
@@ -1191,9 +1198,9 @@ func (l *rowLine) crossing(x, ya, yb int32) int32 {
 		if d < 0 {
 			d = -d
 		}
-		return max(min(l.fy0+int32(float64(d)*l.dydx+0.5), yb), ya)
+		return max(min(l.fy0+int32(float64(float64(d)*l.dydx)+0.5), yb), ya)
 	}
-	t := (float64(x)+l.offx-l.c)*l.dydx - l.base
+	t := float64((float64(x)+l.offx-l.c)*l.dydx) - l.base
 	if !(t > float64(ya)) { // also NaN
 		return ya
 	}
@@ -1221,7 +1228,7 @@ func cellLine(acc []int64, dirty []uint64, x0, fy0, x1, fy1, dir int32, dydx flo
 		xa, ya := x0, fy0
 		for c := c0; c < c1; c++ {
 			bx := (c + 1) << 8
-			yb := min(fy0+int32(float64(bx-x0)*dydx+0.5), fy1)
+			yb := min(fy0+int32(float64(float64(bx-x0)*dydx)+0.5), fy1)
 			d := (yb - ya) * dir
 			s := (xa - c<<8) + 256
 			acc[c] += int64(d * (512 - s))
@@ -1238,7 +1245,7 @@ func cellLine(acc []int64, dirty []uint64, x0, fy0, x1, fy1, dir int32, dydx flo
 	xa, ya := x0, fy0
 	for c := c0; c > c1; c-- {
 		bx := c << 8
-		yb := min(fy0+int32(float64(x0-bx)*dydx+0.5), fy1)
+		yb := min(fy0+int32(float64(float64(x0-bx)*dydx)+0.5), fy1)
 		d := (yb - ya) * dir
 		s := xa - c<<8
 		acc[c] += int64(d * (512 - s))
diff --git a/strip.go b/strip.go
index 203c65f..e541fc4 100644
--- a/strip.go
+++ b/strip.go
@@ -497,12 +497,13 @@ func (rc *rowCtx) span(j, cx0, cx1 int) (i0, i1 int, fy float64) {
 // a pixel gets must not. Loops over a row count the centres' x instead,
 // which is exact (see rowAt).
 func (rc *rowCtx) centre(i0 int, fy float64) float64 {
-	return rc.nx*(float64(i0)+0.5) + rc.ny*(fy+0.5)
+	return float64(rc.nx*(float64(i0)+0.5)) + float64(rc.ny*(fy+0.5))
 }
 
 // rowAt returns the x of the centre of pixel i0 and the part of n·c that
 // row fy adds: the centre of the pixel k further has n·c = nx·(x+k) + g,
-// which is centre(i0+k, fy) to the bit.
+// which is centre(i0+k, fy) to the bit (the products rounded before they
+// are added, on every machine).
 func (rc *rowCtx) rowAt(i0 int, fy float64) (x, g float64) {
 	return float64(i0) + 0.5, rc.ny * (fy + 0.5)
 }
@@ -557,7 +558,7 @@ func (f *segFast) solidRows(y0, y1 int, rc *rowCtx) {
 		row := solid.t.row(j, i0, i1)
 		x, g := rc.rowAt(i0, fy)
 		for i, dst := range row {
-			d := nx*x + g
+			d := float64(nx*x) + g
 			x++
 			a := quant(tr.area(k1-d) - tr.area(k2-d))
 			switch a {
@@ -604,7 +605,7 @@ func (f *segFast) blitRows(y0, y1 int, rc *rowCtx) {
 		cov := f.cov[:i1-i0]
 		x, g := rc.rowAt(i0, fy)
 		for i := range cov {
-			d := nx*x + g
+			d := float64(nx*x) + g
 			x++
 			cov[i] = quant(tr.area(k1-d) - tr.area(k2-d))
 		}
@@ -648,7 +649,7 @@ func (f *segFast) emitRow(j, i0, i1 int, fy float64, rc *rowCtx) {
 	cov := f.cov[:n]
 	x, g := rc.rowAt(i0, fy)
 	for i := 0; i < a; i++ {
-		cov[i] = quant(rc.cov(rc.nx*x + g))
+		cov[i] = quant(rc.cov(float64(rc.nx*x) + g))
 		x++
 	}
 	if a > 0 {
@@ -659,7 +660,7 @@ func (f *segFast) emitRow(j, i0, i1 int, fy float64, rc *rowCtx) {
 	}
 	x, _ = rc.rowAt(i0+b, fy)
 	for i := b; i < n; i++ {
-		cov[i] = quant(rc.cov(rc.nx*x + g))
+		cov[i] = quant(rc.cov(float64(rc.nx*x) + g))
 		x++
 	}
 	if n > b {
@@ -696,7 +697,7 @@ func (f *segFast) borderRow(j, i0, i1 int, fy float64, rc *rowCtx) {
 	inRun, runStart := false, 0
 	xc, g := rc.rowAt(i0, fy)
 	for x := i0; x < i1; x++ {
-		c := rc.cov(rc.nx*xc + g)
+		c := rc.cov(float64(rc.nx*xc) + g)
 		xc++
 		if c >= 1-1e-12 {
 			if !inRun {