diff --git a/x/h3go/area.go b/x/h3go/area.go
index faaf50f..6568ac7 100644
--- a/x/h3go/area.go
+++ b/x/h3go/area.go
@@ -72,13 +72,13 @@ func (b CellBoundary) areaRads2() float64 {
 // arc from ll to other (lat/lng in degrees), following the d3-geo spherical-area
 // formulation.
 func (ll LatLng) cagnoli(other LatLng) float64 {
-	lat := ll.Lat*degsToRads/2 + math.Pi/4
-	otherLat := other.Lat*degsToRads/2 + math.Pi/4
+	lat := ll.Lat*DegsToRads/2 + math.Pi/4
+	otherLat := other.Lat*DegsToRads/2 + math.Pi/4
 
 	sa := math.Sin(lat) * math.Sin(otherLat)
 	ca := math.Cos(lat) * math.Cos(otherLat)
 
-	delta := (other.Lng - ll.Lng) * degsToRads
+	delta := (other.Lng - ll.Lng) * DegsToRads
 	sinDelta := math.Sin(delta)
 	cosDelta := math.Cos(delta)
 
diff --git a/x/h3go/area_test.go b/x/h3go/area_test.go
index 096d364..cf7e0c8 100644
--- a/x/h3go/area_test.go
+++ b/x/h3go/area_test.go
@@ -98,12 +98,45 @@ func TestCellAreaRes0SumsToSphere(t *testing.T) {
 	assertRelClose(t, sum, 4*math.Pi, "res0 area sum")
 }
 
+// TestCellAreaNullIslandTable ports the testH3CellArea.c specific_cell_area
+// regression: the exact area in km² of the cell containing (0, 0) at each
+// resolution, to a tight absolute tolerance.
+func TestCellAreaNullIslandTable(t *testing.T) {
+	t.Parallel()
+
+	areasKm2 := []float64{
+		2.562182162955496e+06, 4.476842017201860e+05, 6.596162242711056e+04,
+		9.228872919002590e+03, 1.318694490797110e+03, 1.879593512281298e+02,
+		2.687164354763186e+01, 3.840848847060638e+00, 5.486939641329893e-01,
+		7.838600808637444e-02, 1.119834221989390e-02, 1.599777169186614e-03,
+		2.285390931423380e-04, 3.264850232091780e-05, 4.664070326136774e-06,
+		6.662957615868888e-07,
+	}
+
+	// The C fixture checks res 0..MAX_H3_RES-1.
+	for res := 0; res < MaxResolution; res++ {
+		cell, err := LatLngToCell(LatLng{Lat: 0, Lng: 0}, res)
+		if err != nil {
+			t.Fatalf("LatLngToCell(res %d): %v", res, err)
+		}
+
+		area, err := CellAreaKm2(cell)
+		if err != nil {
+			t.Fatalf("CellAreaKm2(res %d): %v", res, err)
+		}
+
+		if math.Abs(area-areasKm2[res]) >= 1e-8 {
+			t.Fatalf("CellAreaKm2(res %d): got %.15e, want %.15e", res, area, areasKm2[res])
+		}
+	}
+}
+
 // TestCellAreaInvalidBaseCell covers the error branch of CellAreaRads2/Km2/M2,
 // reached when the boundary cannot be computed for an out-of-range base cell.
 func TestCellAreaInvalidBaseCell(t *testing.T) {
 	t.Parallel()
 
-	bad := Cell(h3Init) | Cell(cellMode)<<modeOffset | Cell(numBaseCells)<<baseCellOffset
+	bad := Cell(h3Init) | Cell(cellMode)<<modeOffset | Cell(NumBaseCells)<<baseCellOffset
 
 	if _, err := CellAreaRads2(bad); err == nil {
 		t.Fatal("CellAreaRads2: got nil error, want failure")
diff --git a/x/h3go/bbox.go b/x/h3go/bbox.go
new file mode 100644
index 0000000..bb72b64
--- /dev/null
+++ b/x/h3go/bbox.go
@@ -0,0 +1,351 @@
+/*
+ * Copyright 2026 Uber Technologies, Inc.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package h3go
+
+import "math"
+
+// bbox is a geographic bounding box with degree coordinates, matching the units
+// of the public GeoPolygon. east < west indicates a box crossing the
+// antimeridian.
+type bbox struct {
+	north, south, east, west float64
+}
+
+// Longitude span constants in degrees, the units of the public API. They stand
+// in for the H3 C library's radian M_PI / M_2PI when detecting and normalizing
+// antimeridian-crossing geometry.
+const (
+	piDeg     = 180.0
+	twoPiDeg  = 360.0
+	halfPiDeg = 90.0
+)
+
+// twoPiRad is a full longitude turn in radians, used to normalize
+// antimeridian-crossing geometry where the math must run at the radian scale of
+// the H3 C library.
+const twoPiRad = 2 * math.Pi
+
+// longitudeNormalization selects how a longitude is shifted so two bounding
+// boxes, either of which may cross the antimeridian, can be compared in one
+// frame of reference.
+type longitudeNormalization int
+
+const (
+	normalizeNone longitudeNormalization = iota
+	normalizeEast
+	normalizeWest
+)
+
+// isTransmeridian reports whether the bounding box crosses the antimeridian.
+func (b bbox) isTransmeridian() bool {
+	return b.east < b.west
+}
+
+// contains reports whether the bounding box contains the point.
+func (b bbox) contains(point LatLng) bool {
+	if point.Lat < b.south || point.Lat > b.north {
+		return false
+	}
+
+	if b.isTransmeridian() {
+		return point.Lng >= b.west || point.Lng <= b.east
+	}
+
+	return point.Lng >= b.west && point.Lng <= b.east
+}
+
+// bboxFromGeoLoop computes the bounding box of a loop of coordinates. It does not
+// support loops with adjacent points more than 180° of longitude apart (treated
+// as antimeridian crossings) or loops containing a pole.
+func bboxFromGeoLoop(loop GeoLoop) bbox {
+	if len(loop) == 0 {
+		return bbox{}
+	}
+
+	out := bbox{north: -math.MaxFloat64, south: math.MaxFloat64, east: -math.MaxFloat64, west: math.MaxFloat64}
+	minPosLng := math.MaxFloat64
+	maxNegLng := -math.MaxFloat64
+	isTransmeridian := false
+
+	for i := range loop {
+		coord := loop[i]
+		next := loop[(i+1)%len(loop)]
+
+		out.south = min(out.south, coord.Lat)
+		out.west = min(out.west, coord.Lng)
+		out.north = max(out.north, coord.Lat)
+		out.east = max(out.east, coord.Lng)
+
+		if coord.Lng > 0 && coord.Lng < minPosLng {
+			minPosLng = coord.Lng
+		}
+
+		if coord.Lng < 0 && coord.Lng > maxNegLng {
+			maxNegLng = coord.Lng
+		}
+
+		if math.Abs(coord.Lng-next.Lng) > piDeg {
+			isTransmeridian = true
+		}
+	}
+
+	if isTransmeridian {
+		out.east = maxNegLng
+		out.west = minPosLng
+	}
+
+	return out
+}
+
+// bboxesFromGeoPolygon returns the bounding box for the outer loop followed by
+// one for each hole, in order.
+func bboxesFromGeoPolygon(polygon GeoPolygon) []bbox {
+	bboxes := make([]bbox, len(polygon.Holes)+1)
+	bboxes[0] = bboxFromGeoLoop(polygon.GeoLoop)
+
+	for i := range polygon.Holes {
+		bboxes[i+1] = bboxFromGeoLoop(polygon.Holes[i])
+	}
+
+	return bboxes
+}
+
+// hexRadiusKm returns the radius of a cell in kilometers, the distance from its
+// center to its first boundary vertex. It is only called with known-valid cells
+// (the pentagons of a resolution), so any projection error is ignored, matching
+// the H3 C library's error-free _hexRadiusKm.
+func (c Cell) hexRadiusKm() float64 {
+	fijk, _ := c.toFaceIjk()
+
+	res := c.Resolution()
+	center := fijk.toVec3(res).toLatLng()
+
+	var boundary CellBoundary
+	if c.IsPentagon() {
+		boundary = fijk.pentToCellBoundary(res, 0, numPentVerts)
+	} else {
+		boundary = fijk.toCellBoundary(res, 0, numHexVerts)
+	}
+
+	return GreatCircleDistanceKm(center, boundary[0])
+}
+
+// bboxHexEstimate estimates the number of cells of the given resolution that fit
+// within the Cartesian-projected bounding box.
+func bboxHexEstimate(box bbox, res int) (int, error) {
+	pentagons, err := Pentagons(res)
+	if err != nil {
+		return 0, err
+	}
+
+	pentagonRadiusKm := pentagons[0].hexRadiusKm()
+
+	// Area of a regular hexagon is 3/2*sqrt(3) * r * r. The pentagon has the most
+	// distortion (smallest edges), shrunk by 20% in case the box perfectly bounds
+	// a pentagon.
+	pentagonAreaKm2 := 0.8 * (2.59807621135 * pentagonRadiusKm * pentagonRadiusKm)
+
+	corner1 := LatLng{Lat: box.north, Lng: box.east}
+	corner2 := LatLng{Lat: box.south, Lng: box.west}
+	diagonalKm := GreatCircleDistanceKm(corner1, corner2)
+
+	lngDiff := math.Abs(corner1.Lng - corner2.Lng)
+	latDiff := math.Abs(corner1.Lat - corner2.Lat)
+
+	if lngDiff == 0 || latDiff == 0 {
+		return 0, ErrFailed
+	}
+
+	length := max(lngDiff, latDiff)
+	width := min(lngDiff, latDiff)
+	ratio := length / width
+
+	// Derived constant, clamped to 3 as higher values drag the estimate to zero.
+	area := diagonalKm * diagonalKm / min(3.0, ratio)
+
+	estimate := math.Ceil(area / pentagonAreaKm2)
+	if math.IsInf(estimate, 0) || math.IsNaN(estimate) {
+		return 0, ErrFailed
+	}
+
+	return max(int(estimate), 1), nil
+}
+
+// lineHexEstimate estimates the number of cells of the given resolution needed
+// to trace the Cartesian-projected line between two points.
+func lineHexEstimate(origin, destination LatLng, res int) (int, error) {
+	pentagons, err := Pentagons(res)
+	if err != nil {
+		return 0, err
+	}
+
+	pentagonRadiusKm := pentagons[0].hexRadiusKm()
+
+	distKm := GreatCircleDistanceKm(origin, destination)
+
+	distCeil := math.Ceil(distKm / (2 * pentagonRadiusKm))
+	if math.IsInf(distCeil, 0) || math.IsNaN(distCeil) {
+		return 0, ErrFailed
+	}
+
+	return max(int(distCeil), 1), nil
+}
+
+// width returns the longitude span of the bounding box in degrees, accounting
+// for an antimeridian crossing.
+func (b bbox) width() float64 {
+	if b.isTransmeridian() {
+		return b.east - b.west + twoPiDeg
+	}
+
+	return b.east - b.west
+}
+
+// height returns the latitude span of the bounding box in degrees.
+func (b bbox) height() float64 {
+	return b.north - b.south
+}
+
+// applyNormalization shifts a longitude east or west by a full turn so that two
+// boxes can be compared in a common frame; normalizeNone leaves it unchanged.
+func applyNormalization(lng float64, normalization longitudeNormalization) float64 {
+	switch normalization {
+	case normalizeEast:
+		if lng < 0 {
+			return lng + twoPiDeg
+		}
+	case normalizeWest:
+		if lng > 0 {
+			return lng - twoPiDeg
+		}
+	case normalizeNone:
+	}
+
+	return lng
+}
+
+// normalizationFor determines the longitude normalization scheme for two
+// bounding boxes, either or both of which may cross the antimeridian, so they
+// can be operated on with standard Cartesian comparisons.
+func normalizationFor(first, second bbox) (longitudeNormalization, longitudeNormalization) {
+	firstTrans := first.isTransmeridian()
+	secondTrans := second.isTransmeridian()
+	firstTrendsEast := first.west-second.east < second.west-first.east
+
+	var firstNorm, secondNorm longitudeNormalization
+
+	switch {
+	case !firstTrans:
+		firstNorm = normalizeNone
+	case secondTrans, firstTrendsEast:
+		firstNorm = normalizeEast
+	default:
+		firstNorm = normalizeWest
+	}
+
+	switch {
+	case !secondTrans:
+		secondNorm = normalizeNone
+	case firstTrans:
+		secondNorm = normalizeEast
+	case firstTrendsEast:
+		secondNorm = normalizeWest
+	default:
+		secondNorm = normalizeEast
+	}
+
+	return firstNorm, secondNorm
+}
+
+// overlaps reports whether two bounding boxes overlap, accounting for
+// antimeridian crossings.
+func (b bbox) overlaps(other bbox) bool {
+	if b.north < other.south || b.south > other.north {
+		return false
+	}
+
+	bNorm, otherNorm := normalizationFor(b, other)
+
+	if applyNormalization(b.east, bNorm) < applyNormalization(other.west, otherNorm) ||
+		applyNormalization(b.west, bNorm) > applyNormalization(other.east, otherNorm) {
+		return false
+	}
+
+	return true
+}
+
+// containsBBox reports whether the bounding box fully contains another box,
+// accounting for antimeridian crossings.
+func (b bbox) containsBBox(other bbox) bool {
+	if b.north < other.north || b.south > other.south {
+		return false
+	}
+
+	bNorm, otherNorm := normalizationFor(b, other)
+
+	return applyNormalization(b.west, bNorm) <= applyNormalization(other.west, otherNorm) &&
+		applyNormalization(b.east, bNorm) >= applyNormalization(other.east, otherNorm)
+}
+
+// toCellBoundary converts the bounding box to a four-vertex cell boundary in
+// counter-clockwise order.
+func (b bbox) toCellBoundary() CellBoundary {
+	return CellBoundary{
+		{Lat: b.north, Lng: b.east},
+		{Lat: b.north, Lng: b.west},
+		{Lat: b.south, Lng: b.west},
+		{Lat: b.south, Lng: b.east},
+	}
+}
+
+// scaled returns the bounding box scaled about its center by the given factor,
+// normalized to the latitude and longitude domains. Both width and height are
+// scaled by the factor (so the area scales by factor squared).
+func (b bbox) scaled(scale float64) bbox {
+	widthBuffer := (b.width()*scale - b.width()) * 0.5
+	heightBuffer := (b.height()*scale - b.height()) * 0.5
+
+	b.north += heightBuffer
+	if b.north > halfPiDeg {
+		b.north = halfPiDeg
+	}
+
+	b.south -= heightBuffer
+	if b.south < -halfPiDeg {
+		b.south = -halfPiDeg
+	}
+
+	b.east += widthBuffer
+	if b.east > piDeg {
+		b.east -= twoPiDeg
+	}
+
+	if b.east < -piDeg {
+		b.east += twoPiDeg
+	}
+
+	b.west -= widthBuffer
+	if b.west > piDeg {
+		b.west -= twoPiDeg
+	}
+
+	if b.west < -piDeg {
+		b.west += twoPiDeg
+	}
+
+	return b
+}
diff --git a/x/h3go/boundary_test.go b/x/h3go/boundary_test.go
index 548e25b..4e7c340 100644
--- a/x/h3go/boundary_test.go
+++ b/x/h3go/boundary_test.go
@@ -58,7 +58,7 @@ func reverseCorpus(t *testing.T) []Cell {
 		latOffset  = 90.0
 		lngSpan    = 360.0
 		lngOffset  = 180.0
-		resCount   = maxResolution + 1
+		resCount   = MaxResolution + 1
 	)
 
 	rng := rand.New(rand.NewSource(seed))
@@ -79,8 +79,8 @@ func reverseCorpus(t *testing.T) []Cell {
 
 // TestCellToBoundarySweep exercises the boundary builders over the reverse
 // corpus and checks each boundary is well-formed: a sane vertex count and every
-// vertex in geographic range. Exact-value correctness is covered by the cgo
-// parity tests.
+// vertex in geographic range. Exact-value correctness is covered by the parity
+// tests.
 func TestCellToBoundarySweep(t *testing.T) {
 	t.Parallel()
 
@@ -114,7 +114,7 @@ func TestCellToBoundarySweep(t *testing.T) {
 func TestCellToBoundaryInvalidBaseCell(t *testing.T) {
 	t.Parallel()
 
-	bad := Cell(h3Init) | Cell(cellMode)<<modeOffset | Cell(numBaseCells)<<baseCellOffset
+	bad := Cell(h3Init) | Cell(cellMode)<<modeOffset | Cell(NumBaseCells)<<baseCellOffset
 
 	if _, err := CellToBoundary(bad); err == nil {
 		t.Fatalf("CellToBoundary(%015x): got nil error, want ErrCellInvalid", uint64(bad))
@@ -210,11 +210,10 @@ func TestCellToBoundaryKnownValues(t *testing.T) {
 }
 
 // TestCellToBoundaryDoublePrecisionVertex is the regression for the double-
-// precision edge-intersection case in the cgo reference: a res-1 pentagon whose
-// distortion vertices shift between float and double precision. The boundary
-// must be geometrically consistent with indexing — a point that indexes to the
-// cell must fall inside the boundary, and a point that does not must fall
-// outside.
+// precision edge-intersection case: a res-1 pentagon whose distortion vertices
+// shift between float and double precision. The boundary must be geometrically
+// consistent with indexing — a point that indexes to the cell must fall inside
+// the boundary, and a point that does not must fall outside.
 func TestCellToBoundaryDoublePrecisionVertex(t *testing.T) {
 	t.Parallel()
 
diff --git a/x/h3go/convert_test.go b/x/h3go/convert_test.go
index c72a4a4..55cf084 100644
--- a/x/h3go/convert_test.go
+++ b/x/h3go/convert_test.go
@@ -42,8 +42,7 @@ func TestIndexToString(t *testing.T) {
 }
 
 // TestIndexFromString covers IndexFromString, including junk/empty input (which
-// yields 0, matching the cgo binding that discards the parse error) and the
-// "0x" prefix.
+// yields 0, discarding the parse error) and the "0x" prefix.
 func TestIndexFromString(t *testing.T) {
 	t.Parallel()
 
diff --git a/x/h3go/corpus_test.go b/x/h3go/corpus_test.go
index 9dd2ebd..39212f5 100644
--- a/x/h3go/corpus_test.go
+++ b/x/h3go/corpus_test.go
@@ -20,7 +20,7 @@ import "testing"
 
 // corpus builds a deterministic set of cells spanning base cells, hexagons, and
 // pentagons across every resolution, for the white-box property and round-trip
-// tests that exercise the public API without the cgo reference.
+// tests that exercise the public API directly.
 func corpus(t *testing.T) []Cell {
 	t.Helper()
 
@@ -36,7 +36,7 @@ func corpus(t *testing.T) []Cell {
 	lats := []float64{-89.9, -45, -23.43, 0, 11.7, 37.7749, 45, 67.1509, 89.9}
 	lngs := []float64{-179.9, -122.4194, -45, 0, 13.4, 100.5, 151.2093, 179.9}
 
-	for res := 0; res <= maxResolution; res++ {
+	for res := 0; res <= MaxResolution; res++ {
 		pents, err := Pentagons(res)
 		if err != nil {
 			t.Fatalf("Pentagons(%d): %v", res, err)
diff --git a/x/h3go/directededge.go b/x/h3go/directededge.go
index 0205b98..59c5d3d 100644
--- a/x/h3go/directededge.go
+++ b/x/h3go/directededge.go
@@ -164,7 +164,7 @@ func (e DirectedEdge) Resolution() int {
 }
 
 // IndexDigit returns the indexing digit of the directed edge at res, for res in
-// [1, maxResolution].
+// [1, MaxResolution].
 func (e DirectedEdge) IndexDigit(res int) (int, error) {
 	return Cell(e).IndexDigit(res)
 }
diff --git a/x/h3go/directededge_test.go b/x/h3go/directededge_test.go
index 68dc9a3..f1ceaf8 100644
--- a/x/h3go/directededge_test.go
+++ b/x/h3go/directededge_test.go
@@ -99,6 +99,17 @@ func TestDirectedEdgeNotNeighbors(t *testing.T) {
 	}
 }
 
+// TestReverseDirectedEdgeFuzzFail ports the testDirectedEdge.c fuzz_fail
+// regression: a malformed index whose decoded origin and destination are not
+// neighbors must report ErrNotNeighbors rather than producing an edge.
+func TestReverseDirectedEdgeFuzzFail(t *testing.T) {
+	t.Parallel()
+
+	if _, err := DirectedEdge(0x1001fff7ff2fbfff).Reverse(); !errors.Is(err, ErrNotNeighbors) {
+		t.Fatalf("Reverse(fuzz): got %v, want ErrNotNeighbors", err)
+	}
+}
+
 // TestDirectedEdgeBoundary checks the edge boundary against the cell boundary it
 // is derived from over the corpus.
 func TestDirectedEdgeBoundary(t *testing.T) {
@@ -223,7 +234,7 @@ func TestDirectedEdgeInvalid(t *testing.T) {
 func TestDirectedEdgeDestinationError(t *testing.T) {
 	t.Parallel()
 
-	corrupt := CellFromString("8928308280fffff").setBaseCell(numBaseCells)
+	corrupt := CellFromString("8928308280fffff").setBaseCell(NumBaseCells)
 	edge := DirectedEdge(corrupt.setMode(directedEdgeMode).setReservedBits(jAxesDigit))
 
 	if _, err := edge.Destination(); err == nil {
diff --git a/x/h3go/faceijk.go b/x/h3go/faceijk.go
index c8e2871..2108ce2 100644
--- a/x/h3go/faceijk.go
+++ b/x/h3go/faceijk.go
@@ -30,7 +30,7 @@ import (
 func (v vec3d) closestFace() (face int, sqd float64) {
 	sqd = 5.0
 
-	for f := range numIcosaFaces {
+	for f := range NumIcosaFaces {
 		s := faceCenterPoint[f].distSq(v)
 		if s < sqd {
 			face = f
@@ -721,7 +721,7 @@ func (c Cell) toFaceIjkWithInitializedFijk(fijk faceIJK) (faceIJK, bool) {
 // of range.
 func (c Cell) toFaceIjk() (faceIJK, error) {
 	baseCell := c.BaseCellNumber()
-	if baseCell >= numBaseCells {
+	if baseCell >= NumBaseCells {
 		return faceIJK{}, ErrCellInvalid
 	}
 
@@ -930,3 +930,90 @@ func (fijk faceIJK) pentToVerts(res int) (int, [numPentVerts]faceIJK) {
 
 	return res, out
 }
+
+// invalidFace marks an unused slot while collecting a cell's icosahedron faces.
+const invalidFace = -1
+
+// IcosahedronFaces returns the icosahedron faces (0-19) that the cell intersects,
+// in no particular order. A hexagon touches one or two faces; a pentagon touches
+// five.
+func (c Cell) IcosahedronFaces() ([]int, error) {
+	res := c.Resolution()
+	isPent := c.IsPentagon()
+
+	// Class II pentagons have every vertex on an icosahedron edge, so the
+	// vertex-based check is ambiguous. Their direct child pentagons cross the
+	// same faces, so use those instead. A Class II pentagon is at an even
+	// resolution below the maximum, so the center child always exists.
+	if isPent && !isResClassIII(res) {
+		childPentagon, _ := c.CenterChild(res + 1)
+
+		return childPentagon.IcosahedronFaces()
+	}
+
+	fijk, err := c.toFaceIjk()
+	if err != nil {
+		return nil, err
+	}
+
+	var (
+		vertexCount int
+		fijkVerts   [numHexVerts]faceIJK
+		adjRes      int
+	)
+
+	if isPent {
+		vertexCount = numPentVerts
+
+		var pentVerts [numPentVerts]faceIJK
+
+		adjRes, pentVerts = fijk.pentToVerts(res)
+		copy(fijkVerts[:], pentVerts[:])
+	} else {
+		vertexCount = numHexVerts
+		adjRes, fijkVerts = fijk.toVerts(res)
+	}
+
+	// A pentagon touches five faces, a hexagon at most two.
+	faceCount := numEdgeCells
+	if isPent {
+		faceCount = numPentVerts
+	}
+
+	faces := make([]int, faceCount)
+	for i := range faces {
+		faces[i] = invalidFace
+	}
+
+	for i := range vertexCount {
+		vert := fijkVerts[i]
+		if isPent {
+			vert, _ = vert.adjustPentVertOverage(adjRes)
+		} else {
+			vert, _ = vert.adjustOverageClassII(adjRes, false, true)
+		}
+
+		// Use the output array as a small hash set: find the first empty slot or
+		// the slot already holding this face.
+		pos := 0
+		for faces[pos] != invalidFace && faces[pos] != vert.face {
+			pos++
+
+			if pos >= faceCount {
+				return nil, ErrFailed
+			}
+		}
+
+		faces[pos] = vert.face
+	}
+
+	out := make([]int, 0, faceCount)
+
+	for _, face := range faces {
+		if face != invalidFace {
+			out = append(out, face)
+		}
+	}
+
+	return out, nil
+}
diff --git a/x/h3go/faceijk_test.go b/x/h3go/faceijk_test.go
index 6c2bd91..ba85398 100644
--- a/x/h3go/faceijk_test.go
+++ b/x/h3go/faceijk_test.go
@@ -19,6 +19,8 @@ package h3go
 import (
 	"errors"
 	"testing"
+
+	"github.com/uber/h3-go/v4/internal/h3core"
 )
 
 // farCoord is far beyond maxFaceCoord and stays out of range even after the
@@ -77,3 +79,117 @@ func TestFaceIjkToH3OutOfRange(t *testing.T) {
 		})
 	}
 }
+
+// TestIcosahedronFacesKnown ports the testGetIcosahedronFaces.c regression cases:
+// single-face and multi-face hexagons, and pentagons at several resolutions.
+func TestIcosahedronFacesKnown(t *testing.T) {
+	t.Parallel()
+
+	validFaces := func(t *testing.T, c Cell) []int {
+		t.Helper()
+
+		faces, err := c.IcosahedronFaces()
+		if err != nil {
+			t.Fatalf("IcosahedronFaces(%015x): %v", uint64(c), err)
+		}
+
+		for _, face := range faces {
+			if face < 0 || face > 19 {
+				t.Fatalf("face %d out of range for %015x", face, uint64(c))
+			}
+		}
+
+		return faces
+	}
+
+	t.Run("single_face_hexes", func(t *testing.T) {
+		t.Parallel()
+
+		// Base cell 16 sits at the center of an icosahedron face, so all of its
+		// children share that single face.
+		baseCell16 := setH3Index(0, 16, centerDigit)
+		for _, childRes := range []int{2, 3} {
+			children, err := baseCell16.Children(childRes)
+			if err != nil {
+				t.Fatalf("Children(%d): %v", childRes, err)
+			}
+
+			for _, child := range children {
+				if got := validFaces(t, child); len(got) != 1 {
+					t.Fatalf("child %015x: got %d faces, want 1", uint64(child), len(got))
+				}
+			}
+		}
+	})
+
+	t.Run("hexagon_with_edge_vertices", func(t *testing.T) {
+		t.Parallel()
+		// Class II pentagon neighbor: one face, two adjacent vertices on an edge.
+		if got := validFaces(t, CellFromString("821c37fffffffff")); len(got) != 1 {
+			t.Fatalf("got %d faces, want 1", len(got))
+		}
+	})
+
+	t.Run("hexagon_with_distortion", func(t *testing.T) {
+		t.Parallel()
+		// Class III pentagon neighbor: distortion spans two faces.
+		if got := validFaces(t, CellFromString("831c06fffffffff")); len(got) != 2 {
+			t.Fatalf("got %d faces, want 2", len(got))
+		}
+	})
+
+	t.Run("hexagon_crossing_faces", func(t *testing.T) {
+		t.Parallel()
+		// Class II hexagon with two vertices on an edge.
+		if got := validFaces(t, CellFromString("821ce7fffffffff")); len(got) != 2 {
+			t.Fatalf("got %d faces, want 2", len(got))
+		}
+	})
+
+	t.Run("pentagons", func(t *testing.T) {
+		t.Parallel()
+		// Class III (res 1), Class II (res 2), and res 15 pentagons on base cell 4.
+		for _, res := range []int{1, 2, 15} {
+			pentagon := setH3Index(res, 4, centerDigit)
+			if !pentagon.IsPentagon() {
+				t.Fatalf("setH3Index(%d,4,0) is not a pentagon", res)
+			}
+
+			if got := validFaces(t, pentagon); len(got) != 5 {
+				t.Fatalf("res %d pentagon: got %d faces, want 5", res, len(got))
+			}
+		}
+	})
+
+	t.Run("base_cell_hexagons", func(t *testing.T) {
+		t.Parallel()
+
+		for bc := range NumBaseCells {
+			if h3core.IsBaseCellPentagon[bc] {
+				continue
+			}
+
+			baseCell := setH3Index(0, bc, centerDigit)
+			if got := validFaces(t, baseCell); len(got) < 1 {
+				t.Fatalf("base cell %d: got no faces", bc)
+			}
+		}
+	})
+}
+
+// TestIcosahedronFacesInvalid covers the defensive error paths: an out-of-range
+// base cell (projection failure) and a malformed index whose vertices span more
+// faces than the maximum (the hash-set overflow guard).
+func TestIcosahedronFacesInvalid(t *testing.T) {
+	t.Parallel()
+
+	corrupt := CellFromString("8928308280fffff").setBaseCell(NumBaseCells)
+	if _, err := corrupt.IcosahedronFaces(); !errors.Is(err, ErrCellInvalid) {
+		t.Fatalf("corrupt base cell: got %v, want ErrCellInvalid", err)
+	}
+
+	overflow := Cell(0x08191d58a34080d2)
+	if _, err := overflow.IcosahedronFaces(); !errors.Is(err, ErrFailed) {
+		t.Fatalf("face-count overflow: got %v, want ErrFailed", err)
+	}
+}
diff --git a/x/h3go/grid.go b/x/h3go/grid.go
index 0537362..65b9213 100644
--- a/x/h3go/grid.go
+++ b/x/h3go/grid.go
@@ -67,7 +67,7 @@ func (c Cell) neighborRotations(dir, rotations int) (Cell, int, error) {
 	newRotations := 0
 
 	oldBaseCell := current.BaseCellNumber()
-	if oldBaseCell < 0 || oldBaseCell >= numBaseCells {
+	if oldBaseCell < 0 || oldBaseCell >= NumBaseCells {
 		return 0, 0, ErrCellInvalid
 	}
 
@@ -571,7 +571,7 @@ func (c Cell) directionForNeighbor(destination Cell) int {
 // base cell's coordinate system on the given face, or invalidRotations if the
 // base cell does not appear on that face.
 func baseCellToCCWrot60(baseCell, face int) int {
-	if face < 0 || face >= numIcosaFaces {
+	if face < 0 || face >= NumIcosaFaces {
 		return invalidRotations
 	}
 
diff --git a/x/h3go/grid_test.go b/x/h3go/grid_test.go
index 61fc24e..eac9a9d 100644
--- a/x/h3go/grid_test.go
+++ b/x/h3go/grid_test.go
@@ -18,6 +18,7 @@ package h3go
 
 import (
 	"errors"
+	"math"
 	"testing"
 )
 
@@ -305,7 +306,7 @@ func TestNeighborRotationsErrors(t *testing.T) {
 		t.Fatalf("neighborRotations(-1): got %v, want ErrFailed", err)
 	}
 
-	badBaseCell := Cell(h3Init) | Cell(cellMode)<<modeOffset | Cell(numBaseCells)<<baseCellOffset
+	badBaseCell := Cell(h3Init) | Cell(cellMode)<<modeOffset | Cell(NumBaseCells)<<baseCellOffset
 	if _, _, err := badBaseCell.neighborRotations(kAxesDigit, 0); !errors.Is(err, ErrCellInvalid) {
 		t.Fatalf("neighborRotations(bad base cell): got %v, want ErrCellInvalid", err)
 	}
@@ -568,3 +569,350 @@ func assertSameSet(t *testing.T, got, want []Cell, msg string) {
 		}
 	}
 }
+
+// TestGridRingUnsafeRegressions ports the testGridRingUnsafe.c cases: the exact
+// rings 1 and 2 around an SF cell, the identity and negative-k contracts, and
+// the pentagon failure modes (near a pentagon and starting on one).
+func TestGridRingUnsafeRegressions(t *testing.T) {
+	t.Parallel()
+
+	// The C fixture builds this cell from a truncated-pi SF coordinate; rebuild
+	// it the same way so the expected rings line up bit for bit.
+	sfHex, err := LatLngToCell(LatLng{
+		Lat: 0.659966917655 * RadsToDegs,
+		Lng: (2*3.14159 - 2.1364398519396) * RadsToDegs,
+	}, 9)
+	if err != nil {
+		t.Fatalf("LatLngToCell(sf): %v", err)
+	}
+
+	t.Run("negative_k", func(t *testing.T) {
+		t.Parallel()
+
+		if _, err := sfHex.GridRingUnsafe(-1); !errors.Is(err, ErrDomain) {
+			t.Fatalf("GridRingUnsafe(-1): got %v, want ErrDomain", err)
+		}
+	})
+
+	t.Run("identity", func(t *testing.T) {
+		t.Parallel()
+
+		ring, err := sfHex.GridRingUnsafe(0)
+		if err != nil {
+			t.Fatalf("GridRingUnsafe(0): %v", err)
+		}
+
+		if len(ring) != 1 || ring[0] != sfHex {
+			t.Fatalf("GridRingUnsafe(0): got %v, want [%015x]", ring, uint64(sfHex))
+		}
+	})
+
+	t.Run("ring1", func(t *testing.T) {
+		t.Parallel()
+
+		want := []Cell{
+			0x89283080ddbffff, 0x89283080c37ffff, 0x89283080c27ffff,
+			0x89283080d53ffff, 0x89283080dcfffff, 0x89283080dc3ffff,
+		}
+
+		got, err := sfHex.GridRingUnsafe(1)
+		if err != nil {
+			t.Fatalf("GridRingUnsafe(1): %v", err)
+		}
+
+		assertSameSet(t, got, want, "ring1")
+	})
+
+	t.Run("ring2", func(t *testing.T) {
+		t.Parallel()
+
+		want := []Cell{
+			0x89283080ca7ffff, 0x89283080cafffff, 0x89283080c33ffff,
+			0x89283080c23ffff, 0x89283080c2fffff, 0x89283080d5bffff,
+			0x89283080d43ffff, 0x89283080d57ffff, 0x89283080d1bffff,
+			0x89283080dc7ffff, 0x89283080dd7ffff, 0x89283080dd3ffff,
+		}
+
+		got, err := sfHex.GridRingUnsafe(2)
+		if err != nil {
+			t.Fatalf("GridRingUnsafe(2): %v", err)
+		}
+
+		assertSameSet(t, got, want, "ring2")
+	})
+
+	t.Run("near_pentagon", func(t *testing.T) {
+		t.Parallel()
+
+		nearPentagon := Cell(0x837405fffffffff)
+		for _, k := range []int{1, 2} {
+			if _, err := nearPentagon.GridRingUnsafe(k); !errors.Is(err, ErrPentagon) {
+				t.Fatalf("GridRingUnsafe(near, %d): got %v, want ErrPentagon", k, err)
+			}
+		}
+	})
+
+	t.Run("on_pentagon", func(t *testing.T) {
+		t.Parallel()
+
+		onPentagon := setH3Index(0, 4, centerDigit)
+		if _, err := onPentagon.GridRingUnsafe(2); !errors.Is(err, ErrPentagon) {
+			t.Fatalf("GridRingUnsafe(on, 2): got %v, want ErrPentagon", err)
+		}
+	})
+}
+
+// TestGridDiskRegressions ports the testGridDisk.c known-array cases: the disk
+// around an SF base cell, and the polar pentagon at res 0 and 1 (which have only
+// 5 neighbors, so the disk holds 6 cells). It checks the exact membership and
+// that GridDiskDistances places the origin at distance 0 and the rest at 1.
+func TestGridDiskRegressions(t *testing.T) {
+	t.Parallel()
+
+	sfHex0, err := LatLngToCell(LatLng{
+		Lat: 0.659966917655 * RadsToDegs,
+		Lng: (2*3.14159 - 2.1364398519396) * RadsToDegs,
+	}, 0)
+	if err != nil {
+		t.Fatalf("LatLngToCell(sf, 0): %v", err)
+	}
+
+	tests := map[string]struct {
+		giveOrigin Cell
+		want       []Cell
+	}{
+		"sf_res0": {
+			giveOrigin: sfHex0,
+			want: []Cell{
+				0x8029fffffffffff, 0x801dfffffffffff, 0x8013fffffffffff,
+				0x8027fffffffffff, 0x8049fffffffffff, 0x8051fffffffffff,
+				0x8037fffffffffff,
+			},
+		},
+		"polar_pentagon_res0": {
+			giveOrigin: setH3Index(0, 4, centerDigit),
+			want: []Cell{
+				0x8009fffffffffff, 0x8007fffffffffff, 0x8001fffffffffff,
+				0x8011fffffffffff, 0x801ffffffffffff, 0x8019fffffffffff,
+			},
+		},
+		"polar_pentagon_res1": {
+			giveOrigin: setH3Index(1, 4, centerDigit),
+			want: []Cell{
+				0x81083ffffffffff, 0x81093ffffffffff, 0x81097ffffffffff,
+				0x8108fffffffffff, 0x8108bffffffffff, 0x8109bffffffffff,
+			},
+		},
+	}
+
+	for name, tt := range tests {
+		t.Run(name, func(t *testing.T) {
+			t.Parallel()
+
+			disk, err := tt.giveOrigin.GridDisk(1)
+			if err != nil {
+				t.Fatalf("GridDisk(1): %v", err)
+			}
+
+			assertSameSet(t, disk, tt.want, "disk")
+
+			rings, err := tt.giveOrigin.GridDiskDistances(1)
+			if err != nil {
+				t.Fatalf("GridDiskDistances(1): %v", err)
+			}
+
+			if len(rings[0]) != 1 || rings[0][0] != tt.giveOrigin {
+				t.Fatalf("ring 0: got %v, want [origin]", rings[0])
+			}
+		})
+	}
+}
+
+// TestNeighborRotationsRegressions ports the testGridDiskInternal.c
+// h3NeighborRotations cases: large rotation counts must be reduced modulo 6
+// without overflow, and out-of-range directions must fail.
+func TestNeighborRotationsRegressions(t *testing.T) {
+	t.Parallel()
+
+	t.Run("rotations_overflow", func(t *testing.T) {
+		t.Parallel()
+
+		origin := setH3Index(0, 0, centerDigit)
+
+		out, rotations, err := origin.neighborRotations(kAxesDigit, 2147483646)
+		if err != nil {
+			t.Fatalf("neighborRotations: %v", err)
+		}
+
+		if want := setH3Index(0, 1, centerDigit); out != want {
+			t.Fatalf("out: got %015x, want %015x", uint64(out), uint64(want))
+		}
+
+		if rotations != 5 {
+			t.Fatalf("rotations: got %d, want 5", rotations)
+		}
+	})
+
+	t.Run("rotations_overflow2", func(t *testing.T) {
+		t.Parallel()
+
+		origin := setH3Index(0, 4, centerDigit)
+
+		out, rotations, err := origin.neighborRotations(jkAxesDigit, math.MaxInt32)
+		if err != nil {
+			t.Fatalf("neighborRotations: %v", err)
+		}
+
+		if want := setH3Index(0, 0, centerDigit); out != want {
+			t.Fatalf("out: got %015x, want %015x", uint64(out), uint64(want))
+		}
+
+		if rotations != 0 {
+			t.Fatalf("rotations: got %d, want 0", rotations)
+		}
+	})
+
+	t.Run("invalid_direction", func(t *testing.T) {
+		t.Parallel()
+
+		origin := Cell(0x811d7ffffffffff)
+		for _, dir := range []int{-1, 7, 100} {
+			if _, _, err := origin.neighborRotations(dir, 0); !errors.Is(err, ErrFailed) {
+				t.Fatalf("neighborRotations(dir %d): got %v, want ErrFailed", dir, err)
+			}
+		}
+	})
+}
+
+// TestGridDiskPentagonArrays ports the testGridDisk.c k3 polar-pentagon and k4
+// pentagon known-array cases. The polar k3 case also checks the per-cell grid
+// distance; the pentagon k4 case checks set membership.
+func TestGridDiskPentagonArrays(t *testing.T) {
+	t.Parallel()
+
+	t.Run("polar_pentagon_k3", func(t *testing.T) {
+		t.Parallel()
+
+		// cell -> expected grid distance from the origin setH3Index(1,4,0).
+		wantDist := map[Cell]int{
+			0x81013ffffffffff: 2, 0x811fbffffffffff: 3, 0x81193ffffffffff: 2,
+			0x81097ffffffffff: 1, 0x81003ffffffffff: 3, 0x81183ffffffffff: 3,
+			0x8111bffffffffff: 3, 0x81077ffffffffff: 2, 0x811f7ffffffffff: 2,
+			0x81067ffffffffff: 3, 0x81093ffffffffff: 1, 0x811e7ffffffffff: 3,
+			0x81083ffffffffff: 0, 0x81117ffffffffff: 2, 0x8101bffffffffff: 3,
+			0x81107ffffffffff: 3, 0x81073ffffffffff: 2, 0x811f3ffffffffff: 2,
+			0x81063ffffffffff: 3, 0x8108fffffffffff: 1, 0x811e3ffffffffff: 3,
+			0x8119bffffffffff: 3, 0x81113ffffffffff: 2, 0x81017ffffffffff: 2,
+			0x81103ffffffffff: 3, 0x8109bffffffffff: 1, 0x81197ffffffffff: 2,
+			0x81007ffffffffff: 3, 0x8108bffffffffff: 1, 0x81187ffffffffff: 3,
+			0x8107bffffffffff: 3,
+		}
+
+		rings, err := setH3Index(1, 4, centerDigit).GridDiskDistances(3)
+		if err != nil {
+			t.Fatalf("GridDiskDistances(3): %v", err)
+		}
+
+		gotDist := make(map[Cell]int)
+
+		for distance, ring := range rings {
+			for _, cell := range ring {
+				gotDist[cell] = distance
+			}
+		}
+
+		if len(gotDist) != len(wantDist) {
+			t.Fatalf("cell count: got %d, want %d", len(gotDist), len(wantDist))
+		}
+
+		for cell, want := range wantDist {
+			if got, ok := gotDist[cell]; !ok || got != want {
+				t.Fatalf("cell %015x: got distance %d (present=%v), want %d", uint64(cell), got, ok, want)
+			}
+		}
+	})
+
+	t.Run("pentagon_k4", func(t *testing.T) {
+		t.Parallel()
+
+		want := []Cell{
+			0x811d7ffffffffff, 0x810c7ffffffffff, 0x81227ffffffffff, 0x81293ffffffffff,
+			0x81133ffffffffff, 0x8136bffffffffff, 0x81167ffffffffff, 0x811d3ffffffffff,
+			0x810c3ffffffffff, 0x81223ffffffffff, 0x81477ffffffffff, 0x8128fffffffffff,
+			0x81367ffffffffff, 0x8112fffffffffff, 0x811cfffffffffff, 0x8123bffffffffff,
+			0x810dbffffffffff, 0x8112bffffffffff, 0x81473ffffffffff, 0x8128bffffffffff,
+			0x81363ffffffffff, 0x811cbffffffffff, 0x81237ffffffffff, 0x810d7ffffffffff,
+			0x81127ffffffffff, 0x8137bffffffffff, 0x81287ffffffffff, 0x8126bffffffffff,
+			0x81177ffffffffff, 0x810d3ffffffffff, 0x81233ffffffffff, 0x8150fffffffffff,
+			0x81123ffffffffff, 0x81377ffffffffff, 0x81283ffffffffff, 0x8102fffffffffff,
+			0x811c3ffffffffff, 0x810cfffffffffff, 0x8122fffffffffff, 0x8113bffffffffff,
+			0x81373ffffffffff, 0x8129bffffffffff, 0x8102bffffffffff, 0x811dbffffffffff,
+			0x810cbffffffffff, 0x8122bffffffffff, 0x81297ffffffffff, 0x81507ffffffffff,
+			0x8136fffffffffff, 0x8127bffffffffff, 0x81137ffffffffff,
+		}
+
+		got, err := setH3Index(1, 14, centerDigit).GridDisk(4)
+		if err != nil {
+			t.Fatalf("GridDisk(4): %v", err)
+		}
+
+		assertSameSet(t, got, want, "pentagon k4")
+	})
+}
+
+// TestGridRingPentagonArrays ports the testGridRing.c k3 polar-pentagon and k4
+// pentagon hollow-ring known-array cases.
+func TestGridRingPentagonArrays(t *testing.T) {
+	t.Parallel()
+
+	tests := map[string]struct {
+		giveOrigin Cell
+		giveK      int
+		want       []Cell
+	}{
+		"polar_pentagon_k3": {
+			giveOrigin: setH3Index(1, 4, centerDigit),
+			giveK:      3,
+			want: []Cell{
+				0x811fbffffffffff, 0x81003ffffffffff, 0x81183ffffffffff, 0x8111bffffffffff,
+				0x81067ffffffffff, 0x811e7ffffffffff, 0x8101bffffffffff, 0x81107ffffffffff,
+				0x81063ffffffffff, 0x811e3ffffffffff, 0x8119bffffffffff, 0x81103ffffffffff,
+				0x81007ffffffffff, 0x81187ffffffffff, 0x8107bffffffffff,
+			},
+		},
+		"pentagon_k4": {
+			giveOrigin: setH3Index(1, 14, centerDigit),
+			giveK:      4,
+			want: []Cell{
+				0x81227ffffffffff, 0x81293ffffffffff, 0x8136bffffffffff, 0x81167ffffffffff,
+				0x81477ffffffffff, 0x810dbffffffffff, 0x81473ffffffffff, 0x81237ffffffffff,
+				0x81127ffffffffff, 0x8126bffffffffff, 0x81177ffffffffff, 0x810d3ffffffffff,
+				0x8150fffffffffff, 0x8102fffffffffff, 0x8129bffffffffff, 0x8102bffffffffff,
+				0x81507ffffffffff, 0x8136fffffffffff, 0x8127bffffffffff, 0x81137ffffffffff,
+			},
+		},
+	}
+
+	for name, tt := range tests {
+		t.Run(name, func(t *testing.T) {
+			t.Parallel()
+
+			got, err := tt.giveOrigin.GridRing(tt.giveK)
+			if err != nil {
+				t.Fatalf("GridRing(%d): %v", tt.giveK, err)
+			}
+
+			assertSameSet(t, got, tt.want, name)
+		})
+	}
+}
+
+// TestGridDiskInvalidDigit ports the testGridDisk.c gridDiskInvalidDigit
+// regression: a malformed index must fail with ErrCellInvalid.
+func TestGridDiskInvalidDigit(t *testing.T) {
+	t.Parallel()
+
+	if _, err := GridDisk(Cell(0x4d4b00fe5c5c3030), 2); !errors.Is(err, ErrCellInvalid) {
+		t.Fatalf("GridDisk(invalid digit): got %v, want ErrCellInvalid", err)
+	}
+}
diff --git a/x/h3go/h3go.go b/x/h3go/h3go.go
index 05a32cf..e29151d 100644
--- a/x/h3go/h3go.go
+++ b/x/h3go/h3go.go
@@ -37,19 +37,47 @@ type LatLng struct {
 // distortion vertices.
 type CellBoundary []LatLng
 
-// Error codes. Messages mirror the cgo-backed h3 package so the two
-// implementations report equivalent failures.
+// Error codes. Messages mirror the H3 C library's error strings.
 var (
-	ErrFailed              = errors.New("the operation failed")
-	ErrDomain              = errors.New("argument was outside of acceptable range")
-	ErrLatLngDomain        = errors.New("latitude or longitude arguments were outside of acceptable range")
-	ErrResolutionDomain    = errors.New("resolution argument was outside of acceptable range")
-	ErrResolutionMismatch  = errors.New("H3Index cell arguments had incompatible resolutions")
-	ErrCellInvalid         = errors.New("H3Index cell argument was not valid")
-	ErrDirectedEdgeInvalid = errors.New("H3Index directed edge argument was not valid")
-	ErrNotNeighbors        = errors.New("H3Index cell arguments were not neighbors")
-	ErrDuplicateInput      = errors.New("duplicate input was encountered in the arguments")
-	ErrPentagon            = errors.New("pentagon distortion was encountered")
+	ErrFailed                = errors.New("the operation failed")
+	ErrDomain                = errors.New("argument was outside of acceptable range")
+	ErrLatLngDomain          = errors.New("latitude or longitude arguments were outside of acceptable range")
+	ErrResolutionDomain      = errors.New("resolution argument was outside of acceptable range")
+	ErrResolutionMismatch    = errors.New("H3Index cell arguments had incompatible resolutions")
+	ErrCellInvalid           = errors.New("H3Index cell argument was not valid")
+	ErrDirectedEdgeInvalid   = errors.New("H3Index directed edge argument was not valid")
+	ErrNotNeighbors          = errors.New("H3Index cell arguments were not neighbors")
+	ErrDuplicateInput        = errors.New("duplicate input was encountered in the arguments")
+	ErrPentagon              = errors.New("pentagon distortion was encountered")
+	ErrMemoryAlloc           = errors.New("necessary memory allocation failed")
+	ErrMemoryBounds          = errors.New("bounds of provided memory were not large enough")
+	ErrOptionInvalid         = errors.New("mode or flags argument was not valid")
+	ErrUndirectedEdgeInvalid = errors.New("H3Index undirected edge argument was not valid")
+	ErrVertexInvalid         = errors.New("H3Index vertex argument was not valid")
+	ErrIndexInvalid          = errors.New("index argument was not valid")
+	ErrBaseCellDomain        = errors.New("base cell number was outside of acceptable range")
+	ErrDigitDomain           = errors.New("child digits invalid")
+	ErrDeletedDigit          = errors.New("deleted subsequence indicates invalid index")
+)
+
+// Exported limits and conversion constants, matching the H3 C library.
+const (
+	// MaxResolution is the finest H3 resolution.
+	MaxResolution = h3core.MaxResolution
+	// MaxCellBndryVerts is the maximum number of vertices in a CellBoundary.
+	MaxCellBndryVerts = 10
+	// NumBaseCells is the number of resolution-0 base cells.
+	NumBaseCells = h3core.NumBaseCells
+	// NumIcosaFaces is the number of faces on the icosahedron.
+	NumIcosaFaces = 20
+	// NumPentagons is the number of pentagons at each resolution.
+	NumPentagons = h3core.NumPentagons
+	// InvalidH3Index is the zero value returned for an invalid index.
+	InvalidH3Index = 0
+	// DegsToRads converts degrees to radians when multiplied.
+	DegsToRads = math.Pi / 180.0
+	// RadsToDegs converts radians to degrees when multiplied.
+	RadsToDegs = 180.0 / math.Pi
 )
 
 // Internal types for the pure Go projection pipeline.
@@ -98,7 +126,6 @@ const (
 	invRes0UGnomonic = 2.61803398874989588842
 	res0UGnomonic    = 0.38196601125010500003
 	maxFaceCoord     = 2
-	numIcosaFaces    = 20
 
 	// numHexVerts and numPentVerts are the topological vertex counts of a
 	// hexagon and a pentagon cell, respectively.
@@ -106,8 +133,8 @@ const (
 	numPentVerts = 5
 
 	// fltEpsilon is the 32-bit float epsilon used to detect when a cell-boundary
-	// edge intersection coincides with an existing vertex (matching the cgo
-	// reference, which compares with FLT_EPSILON).
+	// edge intersection coincides with an existing vertex, matching the H3 C
+	// library's use of FLT_EPSILON.
 	fltEpsilon = 1.1920928955078125e-07
 
 	// faceNeighbors quadrant indices: the direction from a face to the adjacent
@@ -116,11 +143,6 @@ const (
 	dirKI = 2
 	dirJK = 3
 
-	// degsToRads converts degrees to radians by multiplying degrees by this constant.
-	degsToRads = math.Pi / 180.0
-	// radsToDegs converts radians to degrees by multiplying radians by this constant.
-	radsToDegs = 180.0 / math.Pi
-
 	// earthRadiusKm is the authalic (equal-area) radius of the Earth in
 	// kilometers, used to convert spherical measures to physical units.
 	earthRadiusKm = 6371.007180918475
@@ -138,7 +160,7 @@ const (
 	invalidDigit = 7
 	numDigits    = 7
 
-	// H3 index bit-layout offsets and masks, shared with the cgo h3 package.
+	// H3 index bit-layout offsets and masks.
 	cellMode         = h3core.CellMode
 	directedEdgeMode = h3core.DirectedEdgeMode
 	vertexMode       = h3core.VertexMode
@@ -148,7 +170,6 @@ const (
 	perDigitOffset   = h3core.PerDigitOffset
 	digitMask        = h3core.DigitMask
 	resolutionMask   = h3core.ResolutionMask
-	maxResolution    = h3core.MaxResolution
 
 	// numCellEdges is the number of directed edges originating at a cell.
 	numCellEdges = 6
@@ -176,7 +197,7 @@ var unitIjkToDigitLUT = [2][2][2]int{
 	},
 }
 
-var faceCenterPoint = [numIcosaFaces]vec3d{
+var faceCenterPoint = [NumIcosaFaces]vec3d{
 	{0.2199307791404606, 0.6583691780274996, 0.7198475378926182},
 	{-0.2139234834501421, 0.1478171829550703, 0.9656017935214205},
 	{0.1092625278784797, -0.4811951572873210, 0.8697775121287253},
@@ -199,7 +220,7 @@ var faceCenterPoint = [numIcosaFaces]vec3d{
 	{-0.1092625278784796, 0.4811951572873210, -0.8697775121287253},
 }
 
-var faceAxesAzRadsCII = [numIcosaFaces][3]float64{
+var faceAxesAzRadsCII = [NumIcosaFaces][3]float64{
 	{5.619958268523939882, 3.525563166130744542, 1.431168063737548730},
 	{5.760339081714187279, 3.665943979320991689, 1.571548876927796127},
 	{0.780213654393430055, 4.969003859179821079, 2.874608756786625655},
@@ -356,7 +377,7 @@ var unitVecs = [7]coordIJK{
 // coordinate into the adjacent face across each of the three axis-pair edges
 // (plus the identity transform for the face itself at index 0). It is indexed by
 // [face][dir] where dir is one of dirIJ, dirKI, dirJK.
-var faceNeighbors = [numIcosaFaces][4]faceOrientIJK{
+var faceNeighbors = [NumIcosaFaces][4]faceOrientIJK{
 	{ // face 0
 		{0, coordIJK{0, 0, 0}, 0}, // central face
 		{4, coordIJK{2, 0, 2}, 1}, // ij quadrant
@@ -482,7 +503,7 @@ var faceNeighbors = [numIcosaFaces][4]faceOrientIJK{
 // adjacentFaceDir gives the direction (dirIJ/dirKI/dirJK) from the origin face
 // to the destination face, in the origin face's coordinate system, 0 if they
 // are the same face, or -1 if the faces are not adjacent.
-var adjacentFaceDir = [numIcosaFaces][numIcosaFaces]int{
+var adjacentFaceDir = [NumIcosaFaces][NumIcosaFaces]int{
 	{0, dirKI, -1, -1, dirIJ, dirJK, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // face 0
 	{dirIJ, 0, dirKI, -1, -1, -1, dirJK, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // face 1
 	{-1, dirIJ, 0, dirKI, -1, -1, -1, dirJK, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}, // face 2
@@ -508,21 +529,21 @@ var adjacentFaceDir = [numIcosaFaces][numIcosaFaces]int{
 // maxDimByCIIres is the maximum IJK dimension value, by Class II resolution,
 // used to detect overage past a face edge. Odd (Class III) entries are -1
 // because the overage check operates on Class II grids only.
-var maxDimByCIIres = [maxResolution + 2]int{
+var maxDimByCIIres = [MaxResolution + 2]int{
 	2, -1, 14, -1, 98, -1, 686, -1, 4802, -1, 33614, -1, 235298, -1, 1647086, -1, 11529602,
 }
 
 // unitScaleByCIIres is the unit-scale distance, by Class II resolution, used to
 // translate IJK coordinates onto an adjacent face. Odd entries are -1 for the
 // same reason as maxDimByCIIres.
-var unitScaleByCIIres = [maxResolution + 2]int{
+var unitScaleByCIIres = [MaxResolution + 2]int{
 	1, -1, 7, -1, 49, -1, 343, -1, 2401, -1, 16807, -1, 117649, -1, 823543, -1, 5764801,
 }
 
 // baseCellHomeFijk maps each base cell to its "home" face and the normalized IJK
 // coordinates of its center on that face — the starting point for decoding a
 // cell back into a face-centered coordinate.
-var baseCellHomeFijk = [numBaseCells]faceIJK{
+var baseCellHomeFijk = [NumBaseCells]faceIJK{
 	{1, coordIJK{1, 0, 0}},  // base cell 0
 	{2, coordIJK{1, 1, 0}},  // base cell 1
 	{1, coordIJK{0, 0, 0}},  // base cell 2
@@ -653,7 +674,7 @@ const invalidBaseCell = 127
 
 // baseCellNeighbors[baseCell][dir] is the base cell reached by stepping from
 // baseCell in direction dir, or invalidBaseCell across a pentagon's deleted edge.
-var baseCellNeighbors = [numBaseCells][7]int{
+var baseCellNeighbors = [NumBaseCells][7]int{
 	{0, 1, 5, 2, 4, 3, 8},                       // base cell 0
 	{1, 7, 6, 9, 0, 3, 2},                       // base cell 1
 	{2, 6, 10, 11, 0, 1, 5},                     // base cell 2
@@ -782,7 +803,7 @@ var baseCellNeighbors = [numBaseCells][7]int{
 
 // baseCellNeighbor60CCWRots[baseCell][dir] is the number of 60° ccw rotations
 // to apply when stepping from baseCell in direction dir.
-var baseCellNeighbor60CCWRots = [numBaseCells][7]int{
+var baseCellNeighbor60CCWRots = [NumBaseCells][7]int{
 	{0, 5, 0, 0, 1, 5, 1},  // base cell 0
 	{0, 0, 1, 0, 1, 0, 1},  // base cell 1
 	{0, 0, 0, 0, 0, 5, 0},  // base cell 2
diff --git a/x/h3go/hierarchy.go b/x/h3go/hierarchy.go
index 7482b42..66e0cb7 100644
--- a/x/h3go/hierarchy.go
+++ b/x/h3go/hierarchy.go
@@ -39,7 +39,7 @@ func (c Cell) zeroIndexDigits(start, end int) Cell {
 	}
 	// Mask with 0s in the [start, end] digit slots and 1s everywhere else.
 	digits := Cell(1)<<(perDigitOffset*(end-start+1)) - 1
-	mask := ^(digits << (perDigitOffset * (maxResolution - end)))
+	mask := ^(digits << (perDigitOffset * (MaxResolution - end)))
 
 	return c & mask
 }
@@ -47,7 +47,7 @@ func (c Cell) zeroIndexDigits(start, end int) Cell {
 // hasChildAtRes reports whether childRes is a valid child resolution for c.
 func (c Cell) hasChildAtRes(childRes int) bool {
 	parentRes := c.Resolution()
-	return childRes >= parentRes && childRes <= maxResolution
+	return childRes >= parentRes && childRes <= MaxResolution
 }
 
 // childrenSize returns the exact number of children of c at childRes, handling
@@ -69,7 +69,7 @@ func (c Cell) childrenSize(childRes int) (int64, error) {
 func (c Cell) Parent(parentRes int) (Cell, error) {
 	childRes := c.Resolution()
 	switch {
-	case parentRes < 0 || parentRes > maxResolution:
+	case parentRes < 0 || parentRes > MaxResolution:
 		return 0, ErrResolutionDomain
 	case parentRes > childRes:
 		return 0, ErrResolutionMismatch
@@ -123,11 +123,11 @@ func (c Cell) ImmediateChildren() ([]Cell, error) {
 
 // childCells returns an iterator over the children of c at childRes. It yields
 // nothing for invalid input (c == 0, or childRes outside [resolution(c),
-// maxResolution]).
+// MaxResolution]).
 func (c Cell) childCells(childRes int) iter.Seq[Cell] {
 	return func(yield func(Cell) bool) {
 		parentRes := c.Resolution()
-		if c == 0 || childRes < parentRes || childRes > maxResolution {
+		if c == 0 || childRes < parentRes || childRes > MaxResolution {
 			return
 		}
 
@@ -182,7 +182,7 @@ func (c Cell) childCells(childRes int) iter.Seq[Cell] {
 // and carries into the next coarser digit.
 func (c Cell) incrementResDigit(res int) Cell {
 	var val Cell = 1
-	val <<= perDigitOffset * (maxResolution - res)
+	val <<= perDigitOffset * (MaxResolution - res)
 
 	return c + val
 }
@@ -270,7 +270,7 @@ func ChildPosToCell(position int, parent Cell, childRes int) (Cell, error) {
 // ChildPosToCell returns the child cell at the given position within an ordered
 // list of all children at childRes.
 func (c Cell) ChildPosToCell(position int, childRes int) (Cell, error) {
-	if childRes < 0 || childRes > maxResolution {
+	if childRes < 0 || childRes > MaxResolution {
 		return 0, ErrResolutionDomain
 	}
 
diff --git a/x/h3go/hierarchy_test.go b/x/h3go/hierarchy_test.go
index 79d6f37..e51b94d 100644
--- a/x/h3go/hierarchy_test.go
+++ b/x/h3go/hierarchy_test.go
@@ -81,7 +81,7 @@ func TestChildrenErrors(t *testing.T) {
 	t.Run("beyond_finest_resolution", func(t *testing.T) {
 		t.Parallel()
 
-		if _, err := c.Children(maxResolution + 1); err == nil {
+		if _, err := c.Children(MaxResolution + 1); err == nil {
 			t.Fatal("Children beyond finest res should fail")
 		}
 	})
@@ -268,7 +268,7 @@ func TestUncompactErrors(t *testing.T) {
 	}{
 		"coarser_than_input":  {giveRes: 4},
 		"negative_resolution": {giveRes: -1},
-		"too_high_resolution": {giveRes: maxResolution + 1},
+		"too_high_resolution": {giveRes: MaxResolution + 1},
 	}
 
 	for name, tt := range tests {
@@ -313,7 +313,7 @@ func TestChildCellsInvalidInput(t *testing.T) {
 	}{
 		"zero_cell":           {giveCell: 0, giveRes: 5},
 		"coarser_than_parent": {giveCell: setIndexCell(5, 0, 0), giveRes: 4},
-		"finer_than_max":      {giveCell: setIndexCell(5, 0, 0), giveRes: maxResolution + 1},
+		"finer_than_max":      {giveCell: setIndexCell(5, 0, 0), giveRes: MaxResolution + 1},
 	}
 
 	for name, tt := range tests {
@@ -447,7 +447,7 @@ func TestHierarchyCorpus(t *testing.T) {
 
 		assertParentRoundTrip(t, c, res)
 
-		if res >= maxResolution {
+		if res >= MaxResolution {
 			continue
 		}
 
@@ -557,7 +557,7 @@ func TestCenterChildErrors(t *testing.T) {
 		giveRes int
 	}{
 		"coarser_than_cell":   {giveRes: 4},
-		"too_high_resolution": {giveRes: maxResolution + 1},
+		"too_high_resolution": {giveRes: MaxResolution + 1},
 	}
 
 	for name, tt := range tests {
@@ -660,3 +660,89 @@ func assertSameCellSet[A ~int64, B ~int64](t *testing.T, got []A, want []B, msg
 		}
 	}
 }
+
+// TestCellToChildrenKnown ports the testCellToChildren.c oneResStep regression:
+// the exact seven res-9 children of a specific res-8 hexagon.
+func TestCellToChildrenKnown(t *testing.T) {
+	t.Parallel()
+
+	parent := Cell(0x88283080ddfffff)
+	want := []Cell{
+		0x89283080dc3ffff, 0x89283080dc7ffff, 0x89283080dcbffff,
+		0x89283080dcfffff, 0x89283080dd3ffff, 0x89283080dd7ffff,
+		0x89283080ddbffff,
+	}
+
+	got, err := parent.Children(9)
+	if err != nil {
+		t.Fatalf("Children(9): %v", err)
+	}
+
+	assertSameSet(t, got, want, "children")
+}
+
+// TestCellToChildrenMultipleResSteps ports the testCellToChildren.c
+// multipleResSteps regression: the exact 49 res-10 children of a res-8 hexagon.
+func TestCellToChildrenMultipleResSteps(t *testing.T) {
+	t.Parallel()
+
+	want := []Cell{
+		0x8a283080dd27fff, 0x8a283080dd37fff, 0x8a283080dc47fff, 0x8a283080dcdffff,
+		0x8a283080dc5ffff, 0x8a283080dc27fff, 0x8a283080ddb7fff, 0x8a283080dc07fff,
+		0x8a283080dd8ffff, 0x8a283080dd5ffff, 0x8a283080dc4ffff, 0x8a283080dd47fff,
+		0x8a283080dce7fff, 0x8a283080dd1ffff, 0x8a283080dceffff, 0x8a283080dc6ffff,
+		0x8a283080dc87fff, 0x8a283080dcaffff, 0x8a283080dd2ffff, 0x8a283080dcd7fff,
+		0x8a283080dd9ffff, 0x8a283080dd6ffff, 0x8a283080dcc7fff, 0x8a283080dca7fff,
+		0x8a283080dccffff, 0x8a283080dd77fff, 0x8a283080dc97fff, 0x8a283080dd4ffff,
+		0x8a283080dd97fff, 0x8a283080dc37fff, 0x8a283080dc8ffff, 0x8a283080dcb7fff,
+		0x8a283080dcf7fff, 0x8a283080dd87fff, 0x8a283080dda7fff, 0x8a283080dc9ffff,
+		0x8a283080dc77fff, 0x8a283080dc67fff, 0x8a283080dc57fff, 0x8a283080ddaffff,
+		0x8a283080dd17fff, 0x8a283080dc17fff, 0x8a283080dd57fff, 0x8a283080dc0ffff,
+		0x8a283080dd07fff, 0x8a283080dc1ffff, 0x8a283080dd0ffff, 0x8a283080dc2ffff,
+		0x8a283080dd67fff,
+	}
+
+	got, err := Cell(0x88283080ddfffff).Children(10)
+	if err != nil {
+		t.Fatalf("Children(10): %v", err)
+	}
+
+	assertSameSet(t, got, want, "multipleResSteps")
+}
+
+// TestCellToChildrenPentagon ports the testCellToChildren.c pentagonChildren
+// regression: the exact 41 res-3 children of a res-1 pentagon.
+func TestCellToChildrenPentagon(t *testing.T) {
+	t.Parallel()
+
+	want := []Cell{
+		0x830800fffffffff, 0x830802fffffffff, 0x830803fffffffff, 0x830804fffffffff,
+		0x830805fffffffff, 0x830806fffffffff, 0x830810fffffffff, 0x830811fffffffff,
+		0x830812fffffffff, 0x830813fffffffff, 0x830814fffffffff, 0x830815fffffffff,
+		0x830816fffffffff, 0x830818fffffffff, 0x830819fffffffff, 0x83081afffffffff,
+		0x83081bfffffffff, 0x83081cfffffffff, 0x83081dfffffffff, 0x83081efffffffff,
+		0x830820fffffffff, 0x830821fffffffff, 0x830822fffffffff, 0x830823fffffffff,
+		0x830824fffffffff, 0x830825fffffffff, 0x830826fffffffff, 0x830828fffffffff,
+		0x830829fffffffff, 0x83082afffffffff, 0x83082bfffffffff, 0x83082cfffffffff,
+		0x83082dfffffffff, 0x83082efffffffff, 0x830830fffffffff, 0x830831fffffffff,
+		0x830832fffffffff, 0x830833fffffffff, 0x830834fffffffff, 0x830835fffffffff,
+		0x830836fffffffff,
+	}
+
+	got, err := Cell(0x81083ffffffffff).Children(3)
+	if err != nil {
+		t.Fatalf("Children(3): %v", err)
+	}
+
+	assertSameSet(t, got, want, "pentagonChildren")
+}
+
+// TestCellToChildrenResTooFine ports the testCellToChildren.c childResTooFine
+// regression: requesting children beyond the maximum resolution fails.
+func TestCellToChildrenResTooFine(t *testing.T) {
+	t.Parallel()
+
+	if _, err := Cell(0x8f283080dcb0ae2).Children(MaxResolution + 1); !errors.Is(err, ErrResolutionDomain) {
+		t.Fatalf("Children(MaxResolution+1): got %v, want ErrResolutionDomain", err)
+	}
+}
diff --git a/x/h3go/index.go b/x/h3go/index.go
index 460bc8e..c920833 100644
--- a/x/h3go/index.go
+++ b/x/h3go/index.go
@@ -27,11 +27,6 @@ const (
 	base16  = 16
 	bitSize = 64
 
-	// numBaseCells is the number of H3 base cells (NUM_BASE_CELLS).
-	numBaseCells = h3core.NumBaseCells
-	// numPentagons is the number of H3 pentagons, the same at every resolution.
-	numPentagons = h3core.NumPentagons
-
 	// reservedOffset is the bit offset of the reserved field (H3_RESERVED_OFFSET).
 	reservedOffset = 56
 	// baseCellMask masks the 7-bit base cell field after shifting.
@@ -43,15 +38,15 @@ const (
 
 	// digitRegionOffset is the number of non-digit bits above the 15×3-bit digit
 	// region (high + mode + reserved + resolution + base cell = 19).
-	digitRegionOffset = bitSize - maxResolution*perDigitOffset
+	digitRegionOffset = bitSize - MaxResolution*perDigitOffset
 
 	// validCellTopBits is the expected value of the top 8 bits of a valid cell:
 	// high bit=0, mode=1 (cell), reserved=000.
 	validCellTopBits = 0b00001000
 )
 
-// pow7 holds precomputed powers of 7: pow7[i] == 7^i for i in [0, maxResolution].
-var pow7 = [maxResolution + 1]int64{
+// pow7 holds precomputed powers of 7: pow7[i] == 7^i for i in [0, MaxResolution].
+var pow7 = [MaxResolution + 1]int64{
 	1,
 	7,
 	49,
@@ -77,12 +72,12 @@ func (c Cell) mode() int {
 
 // indexDigit returns the indexing digit of the index at res.
 func (c Cell) indexDigit(res int) int {
-	return int((c >> ((maxResolution - res) * perDigitOffset)) & digitMask)
+	return int((c >> ((MaxResolution - res) * perDigitOffset)) & digitMask)
 }
 
 // setIndexDigit returns the index with the digit at res set to digit.
 func (c Cell) setIndexDigit(res, digit int) Cell {
-	shift := (maxResolution - res) * perDigitOffset
+	shift := (MaxResolution - res) * perDigitOffset
 	c &= ^(Cell(digitMask) << shift)
 	c |= Cell(digit) << shift
 
@@ -128,9 +123,9 @@ func (c Cell) BaseCellNumber() int {
 }
 
 // NumCells returns the number of cells at the given resolution. Resolutions
-// outside [0, maxResolution] return 0.
+// outside [0, MaxResolution] return 0.
 func NumCells(res int) int {
-	if res < 0 || res > maxResolution {
+	if res < 0 || res > MaxResolution {
 		return 0
 	}
 	// See h3api.h for the formula derivation.
@@ -173,7 +168,7 @@ func (c Cell) IsPentagon() bool {
 // IndexDigit returns the indexing digit of the cell at res, which starts at 1
 // for resolution 1 up to and including resolution 15.
 func (c Cell) IndexDigit(res int) (int, error) {
-	if res < 1 || res > maxResolution {
+	if res < 1 || res > MaxResolution {
 		return 0, ErrResolutionDomain
 	}
 
@@ -196,7 +191,7 @@ func (c Cell) IsValid() bool {
 	}
 
 	bc := c.BaseCellNumber()
-	if bc >= numBaseCells {
+	if bc >= NumBaseCells {
 		return false
 	}
 
@@ -220,7 +215,7 @@ func (c Cell) IsValid() bool {
 func hasAny7UptoRes(h uint64, res int) bool {
 	const mhi uint64 = 0b100100100100100100100100100100100100100100100
 	const mlo = mhi >> 2
-	shift := perDigitOffset * (maxResolution - res)
+	shift := perDigitOffset * (MaxResolution - res)
 	h >>= shift
 	h <<= shift
 	h = h & mhi & (^h - mlo)
@@ -230,7 +225,7 @@ func hasAny7UptoRes(h uint64, res int) bool {
 
 // hasAll7AfterRes reports whether all unused digits after res are set to 7.
 func hasAll7AfterRes(h uint64, res int) bool {
-	if res >= maxResolution {
+	if res >= MaxResolution {
 		return true
 	}
 	shift := digitRegionOffset + perDigitOffset*res
diff --git a/x/h3go/index_test.go b/x/h3go/index_test.go
index f21da91..c3a4f7d 100644
--- a/x/h3go/index_test.go
+++ b/x/h3go/index_test.go
@@ -61,7 +61,7 @@ func TestIsValidCellBitPatterns(t *testing.T) {
 	t.Run("base_cell_out_of_range", func(t *testing.T) {
 		t.Parallel()
 
-		h := setH3Index(0, numBaseCells, centerDigit)
+		h := setH3Index(0, NumBaseCells, centerDigit)
 		if h.IsValid() {
 			t.Fatal("isValidCell should fail for an out-of-range base cell")
 		}
@@ -90,7 +90,7 @@ func TestIsValidCellBitPatterns(t *testing.T) {
 	t.Run("more_deleted_subsequence", func(t *testing.T) {
 		t.Parallel()
 
-		for res := 1; res <= maxResolution; res++ {
+		for res := 1; res <= MaxResolution; res++ {
 			pent := setH3Index(res, 4, centerDigit) // pentagon center child
 			if !pent.IsValid() {
 				t.Fatalf("res %d: pentagon center child should be valid", res)
@@ -206,7 +206,7 @@ func TestIntrospectionCorpus(t *testing.T) {
 		}
 
 		res := c.Resolution()
-		if res < 0 || res > maxResolution {
+		if res < 0 || res > MaxResolution {
 			t.Fatalf("cell %015x: resolution %d out of range", uint64(c), res)
 		}
 
@@ -218,7 +218,7 @@ func TestIntrospectionCorpus(t *testing.T) {
 			t.Fatalf("cell %015x: BaseCellNumber method and func disagree", uint64(c))
 		}
 
-		if bc := c.BaseCellNumber(); bc < 0 || bc >= numBaseCells {
+		if bc := c.BaseCellNumber(); bc < 0 || bc >= NumBaseCells {
 			t.Fatalf("cell %015x: base cell %d out of range", uint64(c), bc)
 		}
 
diff --git a/x/h3go/latlng.go b/x/h3go/latlng.go
index f693662..1de34e4 100644
--- a/x/h3go/latlng.go
+++ b/x/h3go/latlng.go
@@ -16,11 +16,54 @@
 
 package h3go
 
-import "math"
+import (
+	"math"
+	"strconv"
+)
+
+// NewLatLng creates a LatLng from a latitude and longitude in degrees.
+func NewLatLng(lat, lng float64) LatLng {
+	return LatLng{Lat: lat, Lng: lng}
+}
+
+// Cell returns the cell at the given resolution that contains the coordinate.
+func (g LatLng) Cell(res int) (Cell, error) {
+	return LatLngToCell(g, res)
+}
+
+// String returns the coordinate formatted as "(lat, lng)" in degrees.
+func (g LatLng) String() string {
+	buf := make([]byte, 0, latLngStringSize)
+	buf = append(buf, '(')
+	buf = strconv.AppendFloat(buf, g.Lat, 'f', latLngFloatPrecision, 64)
+	buf = append(buf, ',', ' ')
+	buf = strconv.AppendFloat(buf, g.Lng, 'f', latLngFloatPrecision, 64)
+	buf = append(buf, ')')
+
+	return string(buf)
+}
+
+// LatLng string formatting parameters, matching the H3 C library.
+const (
+	latLngFloatPrecision = 5
+	latLngStringSize     = 32
+)
+
+// LatLngToCellString returns the string form of the cell at resolution that
+// contains the coordinate. It is a convenience wrapper for LatLngToCell followed
+// by Cell.String.
+func LatLngToCellString(latitude, longitude float64, res int) (string, error) {
+	cell, err := NewLatLng(latitude, longitude).Cell(res)
+	if err != nil {
+		return "", err
+	}
+
+	return cell.String(), nil
+}
 
 // LatLngToCell returns the Cell at resolution for a geographic coordinate.
 func LatLngToCell(latLng LatLng, res int) (Cell, error) {
-	if res < 0 || res > maxResolution {
+	if res < 0 || res > MaxResolution {
 		return 0, ErrResolutionDomain
 	}
 
@@ -29,8 +72,8 @@ func LatLngToCell(latLng LatLng, res int) (Cell, error) {
 		return 0, ErrLatLngDomain
 	}
 
-	lat := latLng.Lat * degsToRads
-	lng := latLng.Lng * degsToRads
+	lat := latLng.Lat * DegsToRads
+	lng := latLng.Lng * DegsToRads
 	v := latLngToVec3(lat, lng)
 	fijk := v.toFaceIjk(res)
 
@@ -57,10 +100,10 @@ func (c Cell) LatLng() (LatLng, error) {
 // GreatCircleDistanceRads returns the great-circle distance between two points
 // in radians, using the haversine formula. The points are given in degrees.
 func GreatCircleDistanceRads(a, b LatLng) float64 {
-	aLat := a.Lat * degsToRads
-	aLng := a.Lng * degsToRads
-	bLat := b.Lat * degsToRads
-	bLng := b.Lng * degsToRads
+	aLat := a.Lat * DegsToRads
+	aLng := a.Lng * DegsToRads
+	bLat := b.Lat * DegsToRads
+	bLng := b.Lng * DegsToRads
 
 	sinLat := math.Sin((bLat - aLat) / 2)
 	sinLng := math.Sin((bLng - aLng) / 2)
@@ -118,7 +161,7 @@ func EdgeLengthM(e DirectedEdge) (float64, error) {
 
 // hexAreaAvgKm2 holds the average hexagon area at each resolution in square
 // kilometers, indexed by resolution.
-var hexAreaAvgKm2 = [maxResolution + 1]float64{
+var hexAreaAvgKm2 = [MaxResolution + 1]float64{
 	4.357449416078383e+06, 6.097884417941332e+05, 8.680178039899720e+04,
 	1.239343465508816e+04, 1.770347654491307e+03, 2.529038581819449e+02,
 	3.612906216441245e+01, 5.161293359717191e+00, 7.373275975944177e-01,
@@ -129,7 +172,7 @@ var hexAreaAvgKm2 = [maxResolution + 1]float64{
 
 // hexAreaAvgM2 holds the average hexagon area at each resolution in square
 // meters, indexed by resolution.
-var hexAreaAvgM2 = [maxResolution + 1]float64{
+var hexAreaAvgM2 = [MaxResolution + 1]float64{
 	4.357449416078390e+12, 6.097884417941339e+11, 8.680178039899731e+10,
 	1.239343465508818e+10, 1.770347654491309e+09, 2.529038581819452e+08,
 	3.612906216441250e+07, 5.161293359717198e+06, 7.373275975944188e+05,
@@ -140,7 +183,7 @@ var hexAreaAvgM2 = [maxResolution + 1]float64{
 
 // hexEdgeLenAvgKm holds the average hexagon edge length at each resolution in
 // kilometers, indexed by resolution.
-var hexEdgeLenAvgKm = [maxResolution + 1]float64{
+var hexEdgeLenAvgKm = [MaxResolution + 1]float64{
 	1281.256011, 483.0568391, 182.5129565, 68.97922179,
 	26.07175968, 9.854090990, 3.724532667, 1.406475763,
 	0.531414010, 0.200786148, 0.075863783, 0.028663897,
@@ -149,7 +192,7 @@ var hexEdgeLenAvgKm = [maxResolution + 1]float64{
 
 // hexEdgeLenAvgM holds the average hexagon edge length at each resolution in
 // meters, indexed by resolution.
-var hexEdgeLenAvgM = [maxResolution + 1]float64{
+var hexEdgeLenAvgM = [MaxResolution + 1]float64{
 	1281256.011, 483056.8391, 182512.9565, 68979.22179,
 	26071.75968, 9854.090990, 3724.532667, 1406.475763,
 	531.4140101, 200.7861476, 75.86378287, 28.66389748,
@@ -159,7 +202,7 @@ var hexEdgeLenAvgM = [maxResolution + 1]float64{
 // HexagonAreaAvgKm2 returns the average area of a hexagon at the given
 // resolution in square kilometers.
 func HexagonAreaAvgKm2(res int) (float64, error) {
-	if res < 0 || res > maxResolution {
+	if res < 0 || res > MaxResolution {
 		return 0, ErrResolutionDomain
 	}
 
@@ -169,7 +212,7 @@ func HexagonAreaAvgKm2(res int) (float64, error) {
 // HexagonAreaAvgM2 returns the average area of a hexagon at the given resolution
 // in square meters.
 func HexagonAreaAvgM2(res int) (float64, error) {
-	if res < 0 || res > maxResolution {
+	if res < 0 || res > MaxResolution {
 		return 0, ErrResolutionDomain
 	}
 
@@ -179,7 +222,7 @@ func HexagonAreaAvgM2(res int) (float64, error) {
 // HexagonEdgeLengthAvgKm returns the average edge length of a hexagon at the
 // given resolution in kilometers.
 func HexagonEdgeLengthAvgKm(res int) (float64, error) {
-	if res < 0 || res > maxResolution {
+	if res < 0 || res > MaxResolution {
 		return 0, ErrResolutionDomain
 	}
 
@@ -189,7 +232,7 @@ func HexagonEdgeLengthAvgKm(res int) (float64, error) {
 // HexagonEdgeLengthAvgM returns the average edge length of a hexagon at the given
 // resolution in meters.
 func HexagonEdgeLengthAvgM(res int) (float64, error) {
-	if res < 0 || res > maxResolution {
+	if res < 0 || res > MaxResolution {
 		return 0, ErrResolutionDomain
 	}
 
@@ -212,8 +255,8 @@ func latLngToVec3(lat, lng float64) vec3d {
 // in degrees.
 func (v vec3d) toLatLng() LatLng {
 	return LatLng{
-		Lat: radsToDegs * math.Asin(v.z),
-		Lng: radsToDegs * math.Atan2(v.y, v.x),
+		Lat: RadsToDegs * math.Asin(v.z),
+		Lng: RadsToDegs * math.Atan2(v.y, v.x),
 	}
 }
 
diff --git a/x/h3go/latlng_test.go b/x/h3go/latlng_test.go
index d4f6b77..ccb06ca 100644
--- a/x/h3go/latlng_test.go
+++ b/x/h3go/latlng_test.go
@@ -65,7 +65,7 @@ func TestLatLngToCellErrors(t *testing.T) {
 		wantErr error
 	}{
 		"resolution_below_zero": {giveLat: 0, giveLng: 0, giveRes: -1, wantErr: ErrResolutionDomain},
-		"resolution_above_max":  {giveLat: 0, giveLng: 0, giveRes: maxResolution + 1, wantErr: ErrResolutionDomain},
+		"resolution_above_max":  {giveLat: 0, giveLng: 0, giveRes: MaxResolution + 1, wantErr: ErrResolutionDomain},
 		"lat_is_nan":            {giveLat: math.NaN(), giveLng: 0, giveRes: 5, wantErr: ErrLatLngDomain},
 		"lng_is_nan":            {giveLat: 0, giveLng: math.NaN(), giveRes: 5, wantErr: ErrLatLngDomain},
 		"lat_is_pos_inf":        {giveLat: math.Inf(1), giveLng: 0, giveRes: 5, wantErr: ErrLatLngDomain},
@@ -94,7 +94,7 @@ func TestLatLngToCellProjectionSweep(t *testing.T) {
 	lats := []float64{-90, -89.9, -67.5, -45, -23.43, 0, 11.7, 37.7749, 45, 67.1509, 89.9, 90}
 	lngs := []float64{-180, -179.9, -122.4194, -73.9857, -45, 0, 13.4, 100.5, 151.2093, 179.9, 180}
 
-	for res := 0; res <= maxResolution; res++ {
+	for res := 0; res <= MaxResolution; res++ {
 		for _, lat := range lats {
 			for _, lng := range lngs {
 				assertValidAtRes(t, lat, lng, res)
@@ -109,7 +109,7 @@ func TestLatLngToCellProjectionSweep(t *testing.T) {
 		latOffset  = 90.0
 		lngSpan    = 360.0
 		lngOffset  = 180.0
-		resCount   = maxResolution + 1
+		resCount   = MaxResolution + 1
 	)
 
 	rng := rand.New(rand.NewSource(seed))
@@ -203,15 +203,15 @@ func TestCellToLatLngRoundTrip(t *testing.T) {
 func TestCellToLatLngInvalidBaseCell(t *testing.T) {
 	t.Parallel()
 
-	bad := Cell(h3Init) | Cell(cellMode)<<modeOffset | Cell(numBaseCells)<<baseCellOffset
+	bad := Cell(h3Init) | Cell(cellMode)<<modeOffset | Cell(NumBaseCells)<<baseCellOffset
 
 	if _, err := CellToLatLng(bad); !errors.Is(err, ErrCellInvalid) {
 		t.Fatalf("CellToLatLng(%015x): got %v, want %v", uint64(bad), err, ErrCellInvalid)
 	}
 }
 
-// TestCellToLatLngInvalidIndex is the regression matching the cgo reference: an
-// all-ones index (whose base cell field is out of range) is rejected.
+// TestCellToLatLngInvalidIndex is the regression for an all-ones index (whose
+// base cell field is out of range), which must be rejected.
 func TestCellToLatLngInvalidIndex(t *testing.T) {
 	t.Parallel()
 
@@ -271,7 +271,7 @@ func TestHexagonAreaAvg(t *testing.T) {
 
 		prevKm, prevM := math.Inf(1), math.Inf(1)
 
-		for res := 0; res <= maxResolution; res++ {
+		for res := 0; res <= MaxResolution; res++ {
 			km, err := HexagonAreaAvgKm2(res)
 			if err != nil {
 				t.Fatalf("HexagonAreaAvgKm2(%d): %v", res, err)
@@ -294,7 +294,7 @@ func TestHexagonAreaAvg(t *testing.T) {
 	t.Run("out_of_range", func(t *testing.T) {
 		t.Parallel()
 
-		for _, res := range []int{-1, maxResolution + 1} {
+		for _, res := range []int{-1, MaxResolution + 1} {
 			if _, err := HexagonAreaAvgKm2(res); !errors.Is(err, ErrResolutionDomain) {
 				t.Fatalf("HexagonAreaAvgKm2(%d): got %v, want %v", res, err, ErrResolutionDomain)
 			}
@@ -324,7 +324,7 @@ func TestHexagonEdgeLengthAvg(t *testing.T) {
 
 		prevKm, prevM := math.Inf(1), math.Inf(1)
 
-		for res := 0; res <= maxResolution; res++ {
+		for res := 0; res <= MaxResolution; res++ {
 			km, err := HexagonEdgeLengthAvgKm(res)
 			if err != nil {
 				t.Fatalf("HexagonEdgeLengthAvgKm(%d): %v", res, err)
@@ -346,7 +346,7 @@ func TestHexagonEdgeLengthAvg(t *testing.T) {
 	t.Run("out_of_range", func(t *testing.T) {
 		t.Parallel()
 
-		for _, res := range []int{-1, maxResolution + 1} {
+		for _, res := range []int{-1, MaxResolution + 1} {
 			if _, err := HexagonEdgeLengthAvgKm(res); !errors.Is(err, ErrResolutionDomain) {
 				t.Fatalf("HexagonEdgeLengthAvgKm(%d): got %v, want %v", res, err, ErrResolutionDomain)
 			}
@@ -357,3 +357,114 @@ func TestHexagonEdgeLengthAvg(t *testing.T) {
 		}
 	})
 }
+
+// TestNewLatLngAndCell covers the NewLatLng constructor and the LatLng.Cell
+// method form of LatLngToCell.
+func TestNewLatLngAndCell(t *testing.T) {
+	t.Parallel()
+
+	ll := NewLatLng(37.775, -122.418)
+	if ll.Lat != 37.775 || ll.Lng != -122.418 {
+		t.Fatalf("NewLatLng: got %+v", ll)
+	}
+
+	viaMethod, err := ll.Cell(9)
+	if err != nil {
+		t.Fatalf("LatLng.Cell: %v", err)
+	}
+
+	viaFunc, err := LatLngToCell(ll, 9)
+	if err != nil {
+		t.Fatalf("LatLngToCell: %v", err)
+	}
+
+	if viaMethod != viaFunc {
+		t.Fatalf("LatLng.Cell %015x != LatLngToCell %015x", uint64(viaMethod), uint64(viaFunc))
+	}
+
+	if _, err := ll.Cell(-1); !errors.Is(err, ErrResolutionDomain) {
+		t.Fatalf("LatLng.Cell(-1): got %v, want ErrResolutionDomain", err)
+	}
+}
+
+// TestLatLngString covers the LatLng stringer formatting.
+func TestLatLngString(t *testing.T) {
+	t.Parallel()
+
+	if got := NewLatLng(37.775, -122.418).String(); got != "(37.77500, -122.41800)" {
+		t.Fatalf("LatLng.String: got %q", got)
+	}
+}
+
+// TestLatLngToCellString covers the coordinate-to-cell-string helper and its
+// error path.
+func TestLatLngToCellString(t *testing.T) {
+	t.Parallel()
+
+	got, err := LatLngToCellString(37.775, -122.418, 9)
+	if err != nil {
+		t.Fatalf("LatLngToCellString: %v", err)
+	}
+
+	want, err := LatLngToCell(NewLatLng(37.775, -122.418), 9)
+	if err != nil {
+		t.Fatalf("LatLngToCell: %v", err)
+	}
+
+	if got != want.String() {
+		t.Fatalf("LatLngToCellString: got %q, want %q", got, want.String())
+	}
+
+	if _, err := LatLngToCellString(0, 0, 16); !errors.Is(err, ErrResolutionDomain) {
+		t.Fatalf("LatLngToCellString(res 16): got %v, want ErrResolutionDomain", err)
+	}
+}
+
+// TestLatLngToCellExtremeCoordinates ports the testH3Index.c
+// latLngToCellExtremeCoordinates regression: absurd but finite coordinates must
+// not crash and must yield a cell without error.
+func TestLatLngToCellExtremeCoordinates(t *testing.T) {
+	t.Parallel()
+
+	tests := map[string]struct {
+		giveLatLng LatLng
+		giveRes    int
+	}{
+		"huge_lng":      {giveLatLng: LatLng{Lat: 0, Lng: 1e45}, giveRes: 14},
+		"huge_lat_lng":  {giveLatLng: LatLng{Lat: 1e46, Lng: 1e45}, giveRes: 15},
+		"huge_negative": {giveLatLng: LatLng{Lat: 2, Lng: -3e39}, giveRes: 0},
+	}
+
+	for name, tt := range tests {
+		t.Run(name, func(t *testing.T) {
+			t.Parallel()
+
+			if _, err := LatLngToCell(tt.giveLatLng, tt.giveRes); err != nil {
+				t.Fatalf("LatLngToCell: got %v, want success", err)
+			}
+		})
+	}
+}
+
+// TestGreatCircleDistanceWrappedLongitude ports the testLatLng.c
+// distanceRads_wrappedLongitude regression: a longitude difference greater than
+// 180° (here -270°) measures the short way around, π/2.
+func TestGreatCircleDistanceWrappedLongitude(t *testing.T) {
+	t.Parallel()
+
+	negativeLongitude := LatLng{Lat: 0, Lng: -270}
+	zero := LatLng{Lat: 0, Lng: 0}
+
+	const wantRads = math.Pi / 2
+
+	// Matches the C reference's EPSILON_RAD (EPSILON_DEG 1e-9, in radians).
+	const epsilonRad = 1e-9 * DegsToRads
+
+	if got := GreatCircleDistanceRads(negativeLongitude, zero); math.Abs(got-wantRads) >= epsilonRad {
+		t.Fatalf("GreatCircleDistanceRads: got %v, want %v", got, wantRads)
+	}
+
+	if got := GreatCircleDistanceRads(zero, negativeLongitude); math.Abs(got-wantRads) >= epsilonRad {
+		t.Fatalf("GreatCircleDistanceRads (swapped): got %v, want %v", got, wantRads)
+	}
+}
diff --git a/x/h3go/localij.go b/x/h3go/localij.go
index 028814d..09e4ed7 100644
--- a/x/h3go/localij.go
+++ b/x/h3go/localij.go
@@ -78,7 +78,7 @@ func (c Cell) cellToLocalIjk(target Cell) (coordIJK, error) {
 	originBaseCell := c.BaseCellNumber()
 	baseCell := target.BaseCellNumber()
 
-	if originBaseCell >= numBaseCells || baseCell >= numBaseCells {
+	if originBaseCell >= NumBaseCells || baseCell >= NumBaseCells {
 		return coordIJK{}, ErrCellInvalid
 	}
 
@@ -208,7 +208,7 @@ func (c Cell) localIjkToCell(ijk coordIJK) (Cell, error) {
 	res := c.Resolution()
 
 	originBaseCell := c.BaseCellNumber()
-	if originBaseCell >= numBaseCells {
+	if originBaseCell >= NumBaseCells {
 		return 0, ErrCellInvalid
 	}
 
diff --git a/x/h3go/localij_test.go b/x/h3go/localij_test.go
index 427c766..63139ef 100644
--- a/x/h3go/localij_test.go
+++ b/x/h3go/localij_test.go
@@ -18,6 +18,7 @@ package h3go
 
 import (
 	"errors"
+	"math"
 	"testing"
 )
 
@@ -91,12 +92,12 @@ func TestCellToLocalIJErrors(t *testing.T) {
 		t.Fatalf("CellToLocalIJ(non-neighbor base cells): got %v, want ErrFailed", err)
 	}
 
-	bad := Cell(h3Init) | Cell(cellMode)<<modeOffset | Cell(numBaseCells)<<baseCellOffset
+	bad := Cell(h3Init) | Cell(cellMode)<<modeOffset | Cell(NumBaseCells)<<baseCellOffset
 	if _, err := CellToLocalIJ(bad, bad); !errors.Is(err, ErrCellInvalid) {
 		t.Fatalf("CellToLocalIJ(bad base cell): got %v, want ErrCellInvalid", err)
 	}
 
-	if _, err := CellToLocalIJ(origin, origin.setBaseCell(numBaseCells)); !errors.Is(err, ErrCellInvalid) {
+	if _, err := CellToLocalIJ(origin, origin.setBaseCell(NumBaseCells)); !errors.Is(err, ErrCellInvalid) {
 		t.Fatalf("CellToLocalIJ(bad target base cell): got %v, want ErrCellInvalid", err)
 	}
 }
@@ -112,7 +113,7 @@ func TestLocalIJToCellErrors(t *testing.T) {
 		t.Fatal("LocalIJToCell(out of range): got nil error, want failure")
 	}
 
-	bad := Cell(h3Init) | Cell(cellMode)<<modeOffset | Cell(numBaseCells)<<baseCellOffset
+	bad := Cell(h3Init) | Cell(cellMode)<<modeOffset | Cell(NumBaseCells)<<baseCellOffset
 	if _, err := LocalIJToCell(bad, CoordIJ{I: 0, J: 0}); !errors.Is(err, ErrCellInvalid) {
 		t.Fatalf("LocalIJToCell(bad origin): got %v, want ErrCellInvalid", err)
 	}
@@ -250,12 +251,165 @@ func TestGridMetricsErrorsPropagate(t *testing.T) {
 		t.Fatal("GridPath(far): got nil error, want failure")
 	}
 
-	bad := Cell(h3Init) | Cell(cellMode)<<modeOffset | Cell(numBaseCells)<<baseCellOffset
+	bad := Cell(h3Init) | Cell(cellMode)<<modeOffset | Cell(NumBaseCells)<<baseCellOffset
 	if _, err := GridDistance(bad, bad); !errors.Is(err, ErrCellInvalid) {
 		t.Fatalf("GridDistance(bad): got %v, want ErrCellInvalid", err)
 	}
 }
 
+// TestGridDistanceRegressions ports the testGridDistance.c regression cases:
+// distances onto and across pentagons, base-cell neighbors, and the
+// resolution-mismatch and invalid-cell error contracts.
+func TestGridDistanceRegressions(t *testing.T) {
+	t.Parallel()
+
+	t.Run("onto_pentagon", func(t *testing.T) {
+		t.Parallel()
+
+		origin := setH3Index(1, 17, centerDigit)
+		tests := map[string]struct {
+			giveDigit int
+			want      int
+		}{
+			"digit_0": {giveDigit: 0, want: 3},
+			"digit_2": {giveDigit: 2, want: 2},
+			"digit_3": {giveDigit: 3, want: 3},
+			"digit_6": {giveDigit: 6, want: 2},
+		}
+
+		for name, tt := range tests {
+			t.Run(name, func(t *testing.T) {
+				t.Parallel()
+
+				target := setH3Index(1, 14, tt.giveDigit)
+
+				got, err := origin.GridDistance(target)
+				if err != nil {
+					t.Fatalf("GridDistance: %v", err)
+				}
+
+				if got != tt.want {
+					t.Fatalf("GridDistance: got %d, want %d", got, tt.want)
+				}
+			})
+		}
+	})
+
+	t.Run("across_pentagon_fails", func(t *testing.T) {
+		t.Parallel()
+		// Both directions are rejected because of pentagon distortion.
+		origin := Cell(0x820c4ffffffffff)
+		destination := Cell(0x821ce7fffffffff)
+
+		if _, err := GridDistance(destination, origin); err == nil {
+			t.Fatal("GridDistance across pentagon: got nil, want error")
+		}
+
+		if _, err := GridDistance(origin, destination); err == nil {
+			t.Fatal("GridDistance across pentagon (reversed): got nil, want error")
+		}
+	})
+
+	t.Run("base_cell_neighbors", func(t *testing.T) {
+		t.Parallel()
+
+		bc1 := setH3Index(0, 15, centerDigit)
+		bc2 := setH3Index(0, 8, centerDigit)
+		bc3 := setH3Index(0, 31, centerDigit)
+		pent1 := setH3Index(0, 4, centerDigit)
+
+		for _, target := range []Cell{pent1, bc2, bc3} {
+			got, err := bc1.GridDistance(target)
+			if err != nil {
+				t.Fatalf("GridDistance(15, %015x): %v", uint64(target), err)
+			}
+
+			if got != 1 {
+				t.Fatalf("GridDistance(15, %015x): got %d, want 1", uint64(target), got)
+			}
+		}
+
+		if _, err := pent1.GridDistance(bc3); err == nil {
+			t.Fatal("GridDistance(pent1, 31): got nil, want error")
+		}
+	})
+
+	t.Run("resolution_mismatch", func(t *testing.T) {
+		t.Parallel()
+
+		_, err := GridDistance(Cell(0x832830fffffffff), Cell(0x822837fffffffff))
+		if !errors.Is(err, ErrResolutionMismatch) {
+			t.Fatalf("got %v, want ErrResolutionMismatch", err)
+		}
+	})
+
+	t.Run("invalid", func(t *testing.T) {
+		t.Parallel()
+
+		invalid := ^Cell(0)
+		if _, err := GridDistance(invalid, invalid); !errors.Is(err, ErrCellInvalid) {
+			t.Fatalf("GridDistance(invalid, invalid): got %v, want ErrCellInvalid", err)
+		}
+
+		bc1 := setH3Index(0, 15, centerDigit)
+		if _, err := GridDistance(bc1, invalid); !errors.Is(err, ErrResolutionMismatch) {
+			t.Fatalf("GridDistance(bc1, invalid): got %v, want ErrResolutionMismatch", err)
+		}
+	})
+}
+
+// TestGridPathRegressions ports the testGridPathCells.c named cases: lines that
+// must fail (crossing multiple faces / pentagon distortion) and one that must
+// succeed where a naive forward interpolation would fail.
+func TestGridPathRegressions(t *testing.T) {
+	t.Parallel()
+
+	tests := map[string]struct {
+		giveOrigin Cell
+		giveDest   Cell
+		wantErr    bool
+	}{
+		"across_multiple_faces_fails": {
+			giveOrigin: Cell(0x85285aa7fffffff),
+			giveDest:   Cell(0x851d9b1bfffffff),
+			wantErr:    true,
+		},
+		"pentagon_reverse_interpolation": {
+			giveOrigin: Cell(0x820807fffffffff),
+			giveDest:   Cell(0x8208e7fffffffff),
+			wantErr:    false,
+		},
+		"known_failure": {
+			giveOrigin: Cell(0x8411b61ffffffff),
+			giveDest:   Cell(0x84016d3ffffffff),
+			wantErr:    true,
+		},
+	}
+
+	for name, tt := range tests {
+		t.Run(name, func(t *testing.T) {
+			t.Parallel()
+
+			path, err := GridPath(tt.giveOrigin, tt.giveDest)
+			if tt.wantErr {
+				if err == nil {
+					t.Fatalf("GridPath: got %v, want error", path)
+				}
+
+				return
+			}
+
+			if err != nil {
+				t.Fatalf("GridPath: %v", err)
+			}
+
+			if path[0] != tt.giveOrigin || path[len(path)-1] != tt.giveDest {
+				t.Fatalf("GridPath: endpoints %015x..%015x", uint64(path[0]), uint64(path[len(path)-1]))
+			}
+		})
+	}
+}
+
 // TestLocalIJPentagonUnfold exercises the pentagon unfolding paths in both
 // directions: each pentagon paired with the cells around it as both origin and
 // target. Reachable pairs round trip; unreachable ones (across pentagon
@@ -529,3 +683,125 @@ func TestGridPathInterpolateError(t *testing.T) {
 		t.Fatalf("gridPathInterpolate(corrupt anchor): got %v, want ErrCellInvalid", err)
 	}
 }
+
+// TestLocalIJToCellOverflow ports the testCellToLocalIj.c overflow regressions:
+// high-magnitude IJ coordinates (and a few hand-picked particular cases) must
+// fail rather than overflow the internal cube-coordinate conversion, at every
+// resolution.
+func TestLocalIJToCellOverflow(t *testing.T) {
+	t.Parallel()
+
+	coords := []CoordIJ{
+		{I: math.MinInt32, J: math.MaxInt32},
+		{I: math.MaxInt32, J: math.MinInt32},
+		{I: math.MinInt32, J: math.MinInt32},
+		{I: 553648127, J: -2145378272},
+		{I: math.MaxInt32 - 10, J: -11},
+		{I: math.MaxInt32 - 10, J: -10},
+		{I: math.MaxInt32 - 10, J: -9},
+	}
+
+	for res := 0; res <= MaxResolution; res++ {
+		origin := setH3Index(res, 2, centerDigit)
+		for _, ij := range coords {
+			if _, err := LocalIJToCell(origin, ij); err == nil {
+				t.Fatalf("LocalIJToCell(res %d, %+v): got nil error, want failure", res, ij)
+			}
+		}
+	}
+}
+
+// TestLocalIJToCellNegative ports the testCellToLocalIj.c invalid_negativeIj
+// regression: a specific index with large-negative IJ components must fail.
+func TestLocalIJToCellNegative(t *testing.T) {
+	t.Parallel()
+
+	index := Cell(0x200f202020202020)
+	ij := CoordIJ{I: -14671840, J: math.MinInt32}
+
+	if _, err := LocalIJToCell(index, ij); err == nil {
+		t.Fatal("LocalIJToCell(negative): got nil error, want failure")
+	}
+}
+
+// TestLocalIJToCellBaseCells ports the testCellToLocalIj.c ijBaseCells
+// regression: exact IJ->cell mappings around a res-0 base cell, including the
+// out-of-range failures.
+func TestLocalIJToCellBaseCells(t *testing.T) {
+	t.Parallel()
+
+	origin := Cell(0x8029fffffffffff)
+
+	tests := map[string]struct {
+		giveIJ  CoordIJ
+		want    Cell
+		wantErr bool
+	}{
+		"self":         {giveIJ: CoordIJ{I: 0, J: 0}, want: 0x8029fffffffffff},
+		"i1":           {giveIJ: CoordIJ{I: 1, J: 0}, want: 0x8051fffffffffff},
+		"i2_out":       {giveIJ: CoordIJ{I: 2, J: 0}, wantErr: true},
+		"j2_out":       {giveIJ: CoordIJ{I: 0, J: 2}, wantErr: true},
+		"negative_out": {giveIJ: CoordIJ{I: -2, J: -2}, wantErr: true},
+	}
+
+	for name, tt := range tests {
+		t.Run(name, func(t *testing.T) {
+			t.Parallel()
+
+			got, err := LocalIJToCell(origin, tt.giveIJ)
+			if tt.wantErr {
+				if err == nil {
+					t.Fatalf("got %015x, want error", uint64(got))
+				}
+
+				return
+			}
+
+			if err != nil || got != tt.want {
+				t.Fatalf("got %015x (%v), want %015x", uint64(got), err, uint64(tt.want))
+			}
+		})
+	}
+}
+
+// TestLocalIJToCellOutOfRange ports the testCellToLocalIj.c ijOutOfRange
+// regression: exact IJ->cell mappings along the i axis, with the far
+// coordinates failing.
+func TestLocalIJToCellOutOfRange(t *testing.T) {
+	t.Parallel()
+
+	origin := Cell(0x81283ffffffffff)
+
+	tests := map[string]struct {
+		giveIJ  CoordIJ
+		want    Cell
+		wantErr bool
+	}{
+		"i0":      {giveIJ: CoordIJ{I: 0, J: 0}, want: 0x81283ffffffffff},
+		"i1":      {giveIJ: CoordIJ{I: 1, J: 0}, want: 0x81293ffffffffff},
+		"i2":      {giveIJ: CoordIJ{I: 2, J: 0}, want: 0x8150bffffffffff},
+		"i3":      {giveIJ: CoordIJ{I: 3, J: 0}, want: 0x8151bffffffffff},
+		"i4_out":  {giveIJ: CoordIJ{I: 4, J: 0}, wantErr: true},
+		"in4_out": {giveIJ: CoordIJ{I: -4, J: 0}, wantErr: true},
+		"j4_out":  {giveIJ: CoordIJ{I: 0, J: 4}, wantErr: true},
+	}
+
+	for name, tt := range tests {
+		t.Run(name, func(t *testing.T) {
+			t.Parallel()
+
+			got, err := LocalIJToCell(origin, tt.giveIJ)
+			if tt.wantErr {
+				if err == nil {
+					t.Fatalf("got %015x, want error", uint64(got))
+				}
+
+				return
+			}
+
+			if err != nil || got != tt.want {
+				t.Fatalf("got %015x (%v), want %015x", uint64(got), err, uint64(tt.want))
+			}
+		})
+	}
+}
diff --git a/x/h3go/multipoly.go b/x/h3go/multipoly.go
new file mode 100644
index 0000000..fee87de
--- /dev/null
+++ b/x/h3go/multipoly.go
@@ -0,0 +1,312 @@
+/*
+ * Copyright 2026 Uber Technologies, Inc.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package h3go
+
+import (
+	"slices"
+	"sort"
+)
+
+// arc is one directed edge of a cell during multipolygon assembly. Arcs form
+// doubly-linked loops (counter-clockwise) and a union-find forest whose roots
+// identify connected components (each component becomes one polygon).
+type arc struct {
+	parent  *arc
+	prev    *arc
+	next    *arc
+	id      DirectedEdge
+	rank    int
+	removed bool
+	visited bool
+}
+
+// Counter-clockwise orderings of a cell's directed edges into its linked loop.
+// idxHex is for hexagons (six edges); idxPent is for pentagons (five edges).
+var (
+	idxHex  = [numCellEdges]int{0, 4, 3, 5, 1, 2}
+	idxPent = [numCellEdges - 1]int{0, 1, 3, 2, 4}
+)
+
+// root returns the representative arc of the connected component, compressing
+// the path along the way.
+func (a *arc) root() *arc {
+	if a.parent == a {
+		return a
+	}
+
+	a.parent = a.parent.root()
+
+	return a.parent
+}
+
+// union merges the connected components of two arcs, attaching the lower-rank
+// root under the higher-rank one.
+func union(first, second *arc) {
+	first = first.root()
+	second = second.root()
+
+	if first.rank < second.rank {
+		first, second = second, first
+	}
+
+	if first != second {
+		first.rank += second.rank
+		second.parent = first
+	}
+}
+
+// validateCellSet checks that the cells are all valid, share one resolution, and
+// contain no duplicates, matching the H3 C library's contract.
+func validateCellSet(cells []Cell) error {
+	res := cells[0].Resolution()
+	for _, cell := range cells {
+		if !cell.IsValid() {
+			return ErrCellInvalid
+		}
+
+		if cell.Resolution() != res {
+			return ErrResolutionMismatch
+		}
+	}
+
+	if len(cells) >= 2 {
+		sorted := slices.Clone(cells)
+		slices.Sort(sorted)
+
+		for i := 1; i < len(sorted); i++ {
+			if sorted[i] == sorted[i-1] {
+				return ErrDuplicateInput
+			}
+		}
+	}
+
+	return nil
+}
+
+// buildArcs creates the arcs for every cell, linking each cell's edges into a
+// counter-clockwise loop and indexing them by edge id for reverse lookup. Every
+// arc in a cell starts in the cell's own connected component.
+func buildArcs(cells []Cell) ([]*arc, map[DirectedEdge]*arc) {
+	var arcs []*arc
+
+	index := make(map[DirectedEdge]*arc)
+
+	for _, cell := range cells {
+		// cell is valid here, so enumerating its edges cannot fail.
+		edges, _ := cell.DirectedEdges()
+		count := len(edges)
+
+		block := make([]*arc, count)
+		for i := range block {
+			block[i] = &arc{id: edges[i], rank: 1}
+		}
+
+		for i := range block {
+			block[i].parent = block[0]
+		}
+
+		order := idxHex[:]
+		if count == numCellEdges-1 {
+			order = idxPent[:]
+		}
+
+		for i := range block {
+			cur := order[i]
+			prev := order[(i-1+count)%count]
+			next := order[(i+1)%count]
+			block[cur].prev = block[prev]
+			block[cur].next = block[next]
+		}
+
+		for i := range block {
+			arcs = append(arcs, block[i])
+			index[block[i].id] = block[i]
+		}
+	}
+
+	return arcs, index
+}
+
+// cancelArcPairs removes each pair of opposite edges shared by two adjacent
+// cells, stitching the linked loops back together and merging the two arcs'
+// connected components. What remains are the outline loops of the cell set.
+func cancelArcPairs(arcs []*arc, index map[DirectedEdge]*arc) {
+	for _, current := range arcs {
+		if current.removed {
+			continue
+		}
+
+		// current.id is a valid edge, so reversing it cannot fail.
+		reversed, _ := current.id.Reverse()
+
+		opposite, ok := index[reversed]
+		if !ok {
+			continue
+		}
+
+		current.removed = true
+		opposite.removed = true
+
+		current.next.prev = opposite.prev
+		current.prev.next = opposite.next
+		opposite.next.prev = current.prev
+		opposite.prev.next = current.next
+
+		union(current, opposite)
+	}
+}
+
+// outlineLoop holds one assembled boundary loop with the connected component it
+// belongs to and its area, used to order loops within and across polygons.
+type outlineLoop struct {
+	loop GeoLoop
+	root Cell
+	area float64
+}
+
+// buildOutlineLoops walks the remaining arcs into boundary loops, recording each
+// loop's connected component and enclosed area. Loops are sorted by component,
+// then by area so that each polygon's loops are contiguous with its outer loop
+// (smallest enclosed area) first.
+func buildOutlineLoops(arcs []*arc) []outlineLoop {
+	for _, current := range arcs {
+		current.visited = false
+	}
+
+	var loops []outlineLoop
+
+	for _, start := range arcs {
+		if start.visited || start.removed {
+			continue
+		}
+
+		var verts GeoLoop
+
+		current := start
+		for {
+			// current.id is valid, so its boundary cannot fail.
+			boundary, _ := current.id.Boundary()
+			verts = append(verts, boundary[:len(boundary)-1]...)
+			current.visited = true
+			current = current.next
+
+			if current.id == start.id {
+				break
+			}
+		}
+
+		loops = append(loops, outlineLoop{
+			root: Cell(start.root().id),
+			loop: verts,
+			area: CellBoundary(verts).areaRads2(),
+		})
+	}
+
+	sort.SliceStable(loops, func(i, j int) bool {
+		if loops[i].root != loops[j].root {
+			return loops[i].root < loops[j].root
+		}
+
+		return loops[i].area < loops[j].area
+	})
+
+	return loops
+}
+
+// assembleMultiPolygon groups contiguous same-component loops into polygons
+// (outer loop first, the rest holes) and orders the polygons by decreasing outer
+// loop area.
+func assembleMultiPolygon(loops []outlineLoop) []GeoPolygon {
+	type sortablePoly struct {
+		poly      GeoPolygon
+		outerArea float64
+	}
+
+	var polys []sortablePoly
+
+	for i := 0; i < len(loops); {
+		j := i
+		for j < len(loops) && loops[j].root == loops[i].root {
+			j++
+		}
+
+		poly := GeoPolygon{GeoLoop: loops[i].loop}
+		for k := i + 1; k < j; k++ {
+			poly.Holes = append(poly.Holes, loops[k].loop)
+		}
+
+		polys = append(polys, sortablePoly{poly: poly, outerArea: loops[i].area})
+		i = j
+	}
+
+	sort.SliceStable(polys, func(i, j int) bool {
+		return polys[i].outerArea > polys[j].outerArea
+	})
+
+	out := make([]GeoPolygon, len(polys))
+	for i := range polys {
+		out[i] = polys[i].poly
+	}
+
+	return out
+}
+
+// globeMultiPolygon returns the eight-triangle representation of the entire
+// globe, used when the cell set covers the whole sphere and leaves no outline.
+func globeMultiPolygon() []GeoPolygon {
+	verts := [8][3]LatLng{
+		{{Lat: halfPiDeg, Lng: 0}, {Lat: 0, Lng: 0}, {Lat: 0, Lng: halfPiDeg}},
+		{{Lat: halfPiDeg, Lng: 0}, {Lat: 0, Lng: halfPiDeg}, {Lat: 0, Lng: piDeg}},
+		{{Lat: halfPiDeg, Lng: 0}, {Lat: 0, Lng: piDeg}, {Lat: 0, Lng: -halfPiDeg}},
+		{{Lat: halfPiDeg, Lng: 0}, {Lat: 0, Lng: -halfPiDeg}, {Lat: 0, Lng: 0}},
+		{{Lat: -halfPiDeg, Lng: 0}, {Lat: 0, Lng: 0}, {Lat: 0, Lng: -halfPiDeg}},
+		{{Lat: -halfPiDeg, Lng: 0}, {Lat: 0, Lng: -halfPiDeg}, {Lat: 0, Lng: -piDeg}},
+		{{Lat: -halfPiDeg, Lng: 0}, {Lat: 0, Lng: -piDeg}, {Lat: 0, Lng: halfPiDeg}},
+		{{Lat: -halfPiDeg, Lng: 0}, {Lat: 0, Lng: halfPiDeg}, {Lat: 0, Lng: 0}},
+	}
+
+	out := make([]GeoPolygon, len(verts))
+	for i := range verts {
+		out[i] = GeoPolygon{GeoLoop: GeoLoop{verts[i][0], verts[i][1], verts[i][2]}}
+	}
+
+	return out
+}
+
+// CellsToMultiPolygon outlines a set of cells as GeoPolygons: each polygon has
+// one counter-clockwise outer loop followed by any clockwise holes, and polygons
+// are ordered by decreasing outer loop area. The cells must all be valid, share
+// one resolution, and contain no duplicates.
+func CellsToMultiPolygon(cells []Cell) ([]GeoPolygon, error) {
+	if len(cells) == 0 {
+		return nil, nil
+	}
+
+	if err := validateCellSet(cells); err != nil {
+		return nil, err
+	}
+
+	arcs, index := buildArcs(cells)
+	cancelArcPairs(arcs, index)
+
+	loops := buildOutlineLoops(arcs)
+	if len(loops) == 0 {
+		return globeMultiPolygon(), nil
+	}
+
+	return assembleMultiPolygon(loops), nil
+}
diff --git a/x/h3go/multipoly_test.go b/x/h3go/multipoly_test.go
new file mode 100644
index 0000000..e4bacbe
--- /dev/null
+++ b/x/h3go/multipoly_test.go
@@ -0,0 +1,448 @@
+/*
+ * Copyright 2026 Uber Technologies, Inc.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package h3go
+
+import (
+	"errors"
+	"testing"
+)
+
+// TestUnionFind covers the union-find rank-swap and already-merged branches.
+func TestUnionFind(t *testing.T) {
+	t.Parallel()
+
+	first := &arc{rank: 1}
+	first.parent = first
+	second := &arc{rank: 1}
+	second.parent = second
+	third := &arc{rank: 1}
+	third.parent = third
+
+	// Equal ranks: no swap; first absorbs second.
+	union(first, second)
+
+	// Lower-rank first absorbs into higher-rank component: triggers the swap.
+	union(third, first)
+
+	// Already in the same component: the merge body is skipped.
+	union(second, third)
+
+	if first.root() != second.root() || second.root() != third.root() {
+		t.Fatal("all arcs should share one root after the unions")
+	}
+}
+
+// TestCellsToMultiPolygonSingle covers a single hexagon: one polygon, one outer
+// loop, no holes.
+func TestCellsToMultiPolygonSingle(t *testing.T) {
+	t.Parallel()
+
+	origin, err := LatLngToCell(LatLng{Lat: 37.78, Lng: -122.42}, 7)
+	if err != nil {
+		t.Fatalf("LatLngToCell: %v", err)
+	}
+
+	polys, err := CellsToMultiPolygon([]Cell{origin})
+	if err != nil {
+		t.Fatalf("CellsToMultiPolygon: %v", err)
+	}
+
+	if len(polys) != 1 || len(polys[0].Holes) != 0 {
+		t.Fatalf("single cell: got %d polygons, %d holes", len(polys), len(polys[0].Holes))
+	}
+
+	if len(polys[0].GeoLoop) != 6 {
+		t.Fatalf("single hexagon outline: got %d verts, want 6", len(polys[0].GeoLoop))
+	}
+}
+
+// TestCellsToMultiPolygonDisk covers a contiguous blob, exercising edge
+// cancellation and union-find merging into a single outer loop.
+func TestCellsToMultiPolygonDisk(t *testing.T) {
+	t.Parallel()
+
+	origin, err := LatLngToCell(LatLng{Lat: 0, Lng: 0}, 6)
+	if err != nil {
+		t.Fatalf("LatLngToCell: %v", err)
+	}
+
+	disk, err := origin.GridDisk(2)
+	if err != nil {
+		t.Fatalf("GridDisk: %v", err)
+	}
+
+	polys, err := CellsToMultiPolygon(disk)
+	if err != nil {
+		t.Fatalf("CellsToMultiPolygon: %v", err)
+	}
+
+	if len(polys) != 1 || len(polys[0].Holes) != 0 {
+		t.Fatalf("disk: got %d polygons, %d holes", len(polys), len(polys[0].Holes))
+	}
+}
+
+// TestCellsToMultiPolygonRing covers a hollow ring, which yields one polygon with
+// a single hole.
+func TestCellsToMultiPolygonRing(t *testing.T) {
+	t.Parallel()
+
+	origin, err := LatLngToCell(LatLng{Lat: 10, Lng: 10}, 6)
+	if err != nil {
+		t.Fatalf("LatLngToCell: %v", err)
+	}
+
+	ring, err := origin.GridRing(2)
+	if err != nil {
+		t.Fatalf("GridRing: %v", err)
+	}
+
+	polys, err := CellsToMultiPolygon(ring)
+	if err != nil {
+		t.Fatalf("CellsToMultiPolygon: %v", err)
+	}
+
+	if len(polys) != 1 || len(polys[0].Holes) != 1 {
+		t.Fatalf("ring: got %d polygons, %d holes (want 1, 1)", len(polys), len(polys[0].Holes))
+	}
+}
+
+// TestCellsToMultiPolygonDisjoint covers two far-apart blobs, which yield two
+// separate polygons.
+func TestCellsToMultiPolygonDisjoint(t *testing.T) {
+	t.Parallel()
+
+	first, err := LatLngToCell(LatLng{Lat: 0, Lng: 0}, 6)
+	if err != nil {
+		t.Fatalf("LatLngToCell first: %v", err)
+	}
+
+	second, err := LatLngToCell(LatLng{Lat: -30, Lng: 60}, 6)
+	if err != nil {
+		t.Fatalf("LatLngToCell second: %v", err)
+	}
+
+	firstDisk, err := first.GridDisk(1)
+	if err != nil {
+		t.Fatalf("GridDisk first: %v", err)
+	}
+
+	secondDisk, err := second.GridDisk(1)
+	if err != nil {
+		t.Fatalf("GridDisk second: %v", err)
+	}
+
+	polys, err := CellsToMultiPolygon(append(firstDisk, secondDisk...))
+	if err != nil {
+		t.Fatalf("CellsToMultiPolygon: %v", err)
+	}
+
+	if len(polys) != 2 {
+		t.Fatalf("disjoint: got %d polygons, want 2", len(polys))
+	}
+}
+
+// TestCellsToMultiPolygonPentagon covers a pentagon-centered set, exercising the
+// five-edge arc ordering.
+func TestCellsToMultiPolygonPentagon(t *testing.T) {
+	t.Parallel()
+
+	pentagon := setH3Index(5, 4, centerDigit)
+	if !pentagon.IsPentagon() {
+		t.Fatalf("fixture %015x is not a pentagon", uint64(pentagon))
+	}
+
+	disk, err := pentagon.GridDisk(1)
+	if err != nil {
+		t.Fatalf("GridDisk: %v", err)
+	}
+
+	polys, err := CellsToMultiPolygon(disk)
+	if err != nil {
+		t.Fatalf("CellsToMultiPolygon: %v", err)
+	}
+
+	if len(polys) != 1 {
+		t.Fatalf("pentagon disk: got %d polygons, want 1", len(polys))
+	}
+}
+
+// TestCellsToMultiPolygonGlobe covers the whole-globe case: all base cells leave
+// no outline, producing the eight-triangle globe representation.
+func TestCellsToMultiPolygonGlobe(t *testing.T) {
+	t.Parallel()
+
+	cells, err := Res0Cells()
+	if err != nil {
+		t.Fatalf("Res0Cells: %v", err)
+	}
+
+	polys, err := CellsToMultiPolygon(cells)
+	if err != nil {
+		t.Fatalf("CellsToMultiPolygon: %v", err)
+	}
+
+	if len(polys) != 8 {
+		t.Fatalf("globe: got %d polygons, want 8", len(polys))
+	}
+
+	for _, poly := range polys {
+		if len(poly.GeoLoop) != 3 {
+			t.Fatalf("globe triangle: got %d verts, want 3", len(poly.GeoLoop))
+		}
+	}
+}
+
+// TestCellsToMultiPolygonEmpty covers the empty-input short circuit.
+func TestCellsToMultiPolygonEmpty(t *testing.T) {
+	t.Parallel()
+
+	polys, err := CellsToMultiPolygon(nil)
+	if err != nil || polys != nil {
+		t.Fatalf("empty: got %v (%v), want nil, nil", polys, err)
+	}
+}
+
+// TestCellsToMultiPolygonErrors covers the validation paths: invalid cell,
+// mismatched resolution, and duplicate input.
+func TestCellsToMultiPolygonErrors(t *testing.T) {
+	t.Parallel()
+
+	origin, err := LatLngToCell(LatLng{Lat: 0, Lng: 0}, 7)
+	if err != nil {
+		t.Fatalf("LatLngToCell: %v", err)
+	}
+
+	parent, err := origin.Parent(6)
+	if err != nil {
+		t.Fatalf("Parent: %v", err)
+	}
+
+	tests := map[string]struct {
+		cells   []Cell
+		wantErr error
+	}{
+		"invalid":   {[]Cell{0}, ErrCellInvalid},
+		"mixed_res": {[]Cell{origin, parent}, ErrResolutionMismatch},
+		"duplicate": {[]Cell{origin, origin}, ErrDuplicateInput},
+	}
+
+	for name, tt := range tests {
+		t.Run(name, func(t *testing.T) {
+			t.Parallel()
+
+			if _, err := CellsToMultiPolygon(tt.cells); !errors.Is(err, tt.wantErr) {
+				t.Fatalf("%s: got %v, want %v", name, err, tt.wantErr)
+			}
+		})
+	}
+}
+
+// TestCellsToMultiPolygonReported ports the testCellsToLinkedMultiPolygon.c
+// cases: exact polygon, loop, and coordinate counts for single cells, contiguous
+// and non-contiguous sets, a ring with a hole, and a pentagon.
+func TestCellsToMultiPolygonReported(t *testing.T) {
+	t.Parallel()
+
+	t.Run("invalid", func(t *testing.T) {
+		t.Parallel()
+
+		if _, err := CellsToMultiPolygon([]Cell{0xfffffffffffffff}); !errors.Is(err, ErrCellInvalid) {
+			t.Fatalf("got %v, want ErrCellInvalid", err)
+		}
+	})
+
+	tests := map[string]struct {
+		giveCells    []Cell
+		wantPolygons int
+		wantOuter    int // coords on the first polygon's outer loop
+		wantHole     int // coords on the first polygon's single hole, 0 if none
+	}{
+		"single_hex": {
+			giveCells:    []Cell{0x890dab6220bffff},
+			wantPolygons: 1,
+			wantOuter:    6,
+		},
+		"contiguous2": {
+			giveCells:    []Cell{0x8928308291bffff, 0x89283082957ffff},
+			wantPolygons: 1,
+			wantOuter:    10,
+		},
+		"non_contiguous2": {
+			giveCells:    []Cell{0x8928308291bffff, 0x89283082943ffff},
+			wantPolygons: 2,
+			wantOuter:    6,
+		},
+		"contiguous3": {
+			giveCells:    []Cell{0x8928308288bffff, 0x892830828d7ffff, 0x8928308289bffff},
+			wantPolygons: 1,
+			wantOuter:    12,
+		},
+		"hole": {
+			giveCells: []Cell{
+				0x892830828c7ffff, 0x892830828d7ffff, 0x8928308289bffff,
+				0x89283082813ffff, 0x8928308288fffff, 0x89283082883ffff,
+			},
+			wantPolygons: 1,
+			wantOuter:    18,
+			wantHole:     6,
+		},
+		"pentagon": {
+			giveCells:    []Cell{0x851c0003fffffff},
+			wantPolygons: 1,
+			wantOuter:    10,
+		},
+	}
+
+	for name, tt := range tests {
+		t.Run(name, func(t *testing.T) {
+			t.Parallel()
+
+			polygons, err := CellsToMultiPolygon(tt.giveCells)
+			if err != nil {
+				t.Fatalf("CellsToMultiPolygon: %v", err)
+			}
+
+			if len(polygons) != tt.wantPolygons {
+				t.Fatalf("polygons: got %d, want %d", len(polygons), tt.wantPolygons)
+			}
+
+			if got := len(polygons[0].GeoLoop); got != tt.wantOuter {
+				t.Fatalf("outer coords: got %d, want %d", got, tt.wantOuter)
+			}
+
+			if tt.wantHole == 0 {
+				if len(polygons[0].Holes) != 0 {
+					t.Fatalf("holes: got %d, want 0", len(polygons[0].Holes))
+				}
+
+				return
+			}
+
+			if len(polygons[0].Holes) != 1 {
+				t.Fatalf("holes: got %d, want 1", len(polygons[0].Holes))
+			}
+
+			if got := len(polygons[0].Holes[0]); got != tt.wantHole {
+				t.Fatalf("hole coords: got %d, want %d", got, tt.wantHole)
+			}
+		})
+	}
+}
+
+// TestCellsToMultiPolygonIssue1049 ports the testCellsToMultiPoly.c issue_1049
+// regression: a 168-cell res-2 set must assemble into exactly 12 polygons, each
+// with no holes.
+func TestCellsToMultiPolygonIssue1049(t *testing.T) {
+	t.Parallel()
+
+	cells := []Cell{
+		0x827487fffffffff, 0x82748ffffffffff, 0x827497fffffffff, 0x82749ffffffffff,
+		0x8274affffffffff, 0x8274c7fffffffff, 0x8274cffffffffff, 0x8274d7fffffffff,
+		0x8274e7fffffffff, 0x8274effffffffff, 0x8274f7fffffffff, 0x82754ffffffffff,
+		0x827c07fffffffff, 0x827c27fffffffff, 0x827c2ffffffffff, 0x827c37fffffffff,
+		0x827d87fffffffff, 0x827d8ffffffffff, 0x827d97fffffffff, 0x827d9ffffffffff,
+		0x827da7fffffffff, 0x827daffffffffff, 0x82801ffffffffff, 0x8280a7fffffffff,
+		0x8280affffffffff, 0x8280b7fffffffff, 0x828197fffffffff, 0x82819ffffffffff,
+		0x8281a7fffffffff, 0x8281b7fffffffff, 0x828207fffffffff, 0x82820ffffffffff,
+		0x828227fffffffff, 0x82822ffffffffff, 0x8282e7fffffffff, 0x828307fffffffff,
+		0x82830ffffffffff, 0x82831ffffffffff, 0x82832ffffffffff, 0x828347fffffffff,
+		0x82834ffffffffff, 0x828357fffffffff, 0x82835ffffffffff, 0x828367fffffffff,
+		0x828377fffffffff, 0x82a447fffffffff, 0x82a457fffffffff, 0x82a45ffffffffff,
+		0x82a467fffffffff, 0x82a46ffffffffff, 0x82a477fffffffff, 0x82a4c7fffffffff,
+		0x82a4cffffffffff, 0x82a4d7fffffffff, 0x82a4e7fffffffff, 0x82a4effffffffff,
+		0x82a4f7fffffffff, 0x82a547fffffffff, 0x82a54ffffffffff, 0x82a557fffffffff,
+		0x82a55ffffffffff, 0x82a567fffffffff, 0x82a577fffffffff, 0x82a837fffffffff,
+		0x82a897fffffffff, 0x82a8a7fffffffff, 0x82a8b7fffffffff, 0x82a917fffffffff,
+		0x82a927fffffffff, 0x82a937fffffffff, 0x82a987fffffffff, 0x82a98ffffffffff,
+		0x82a997fffffffff, 0x82a99ffffffffff, 0x82a9a7fffffffff, 0x82a9affffffffff,
+		0x82ac47fffffffff, 0x82ac57fffffffff, 0x82ac5ffffffffff, 0x82ac67fffffffff,
+		0x82ac6ffffffffff, 0x82ac77fffffffff, 0x82ad47fffffffff, 0x82ad4ffffffffff,
+		0x82ad57fffffffff, 0x82ad5ffffffffff, 0x82ad67fffffffff, 0x82ad77fffffffff,
+		0x82c207fffffffff, 0x82c217fffffffff, 0x82c227fffffffff, 0x82c237fffffffff,
+		0x82c287fffffffff, 0x82c28ffffffffff, 0x82c29ffffffffff, 0x82c2a7fffffffff,
+		0x82c2affffffffff, 0x82c2b7fffffffff, 0x82c307fffffffff, 0x82c317fffffffff,
+		0x82c31ffffffffff, 0x82c337fffffffff, 0x82cfb7fffffffff, 0x82d0c7fffffffff,
+		0x82d0d7fffffffff, 0x82d0dffffffffff, 0x82d0e7fffffffff, 0x82d0f7fffffffff,
+		0x82d147fffffffff, 0x82d157fffffffff, 0x82d15ffffffffff, 0x82d167fffffffff,
+		0x82d177fffffffff, 0x82d187fffffffff, 0x82d18ffffffffff, 0x82d197fffffffff,
+		0x82d19ffffffffff, 0x82d1a7fffffffff, 0x82d1affffffffff, 0x82dc47fffffffff,
+		0x82dc57fffffffff, 0x82dc5ffffffffff, 0x82dc67fffffffff, 0x82dc6ffffffffff,
+		0x82dc77fffffffff, 0x82dcc7fffffffff, 0x82dccffffffffff, 0x82dcd7fffffffff,
+		0x82dce7fffffffff, 0x82dceffffffffff, 0x82dcf7fffffffff, 0x82dd1ffffffffff,
+		0x82dd47fffffffff, 0x82dd4ffffffffff, 0x82dd57fffffffff, 0x82dd5ffffffffff,
+		0x82dd6ffffffffff, 0x82dd87fffffffff, 0x82dd8ffffffffff, 0x82dd97fffffffff,
+		0x82dd9ffffffffff, 0x82ddaffffffffff, 0x82ddb7fffffffff, 0x82dec7fffffffff,
+		0x82decffffffffff, 0x82ded7fffffffff, 0x82dee7fffffffff, 0x82deeffffffffff,
+		0x82def7fffffffff, 0x82df0ffffffffff, 0x82df1ffffffffff, 0x82df47fffffffff,
+		0x82df4ffffffffff, 0x82df57fffffffff, 0x82df5ffffffffff, 0x82df77fffffffff,
+		0x82df8ffffffffff, 0x82df9ffffffffff, 0x82e6c7fffffffff, 0x82e6cffffffffff,
+		0x82e6d7fffffffff, 0x82e6dffffffffff, 0x82e6effffffffff, 0x82e6f7fffffffff,
+	}
+
+	polygons, err := CellsToMultiPolygon(cells)
+	if err != nil {
+		t.Fatalf("CellsToMultiPolygon: %v", err)
+	}
+
+	if len(polygons) != 12 {
+		t.Fatalf("polygons: got %d, want 12", len(polygons))
+	}
+
+	for i := range polygons {
+		if len(polygons[i].Holes) != 0 {
+			t.Fatalf("polygon %d: got %d holes, want 0", i, len(polygons[i].Holes))
+		}
+	}
+}
+
+// TestCellsToMultiPolygonEquator ports the testCellsToMultiPoly.c equator_cells
+// regression: a global band of cells assembles into a single polygon with one
+// hole.
+func TestCellsToMultiPolygonEquator(t *testing.T) {
+	t.Parallel()
+
+	cells := []Cell{
+		0x81807ffffffffff, 0x817efffffffffff, 0x81723ffffffffff, 0x817ebffffffffff,
+		0x817c3ffffffffff, 0x817e3ffffffffff, 0x817a3ffffffffff, 0x8166fffffffffff,
+		0x8172bffffffffff, 0x816afffffffffff, 0x81933ffffffffff, 0x8168fffffffffff,
+		0x8188fffffffffff, 0x81853ffffffffff, 0x817f7ffffffffff, 0x8180bffffffffff,
+		0x81783ffffffffff, 0x81743ffffffffff, 0x8170bffffffffff, 0x8173bffffffffff,
+		0x8179bffffffffff, 0x817cbffffffffff, 0x8188bffffffffff, 0x81857ffffffffff,
+		0x816f7ffffffffff, 0x8177bffffffffff, 0x81617ffffffffff, 0x816f3ffffffffff,
+		0x8174bffffffffff, 0x8180fffffffffff, 0x817a7ffffffffff, 0x81767ffffffffff,
+		0x81757ffffffffff, 0x81957ffffffffff, 0x81787ffffffffff, 0x81847ffffffffff,
+		0x81653ffffffffff, 0x817bbffffffffff, 0x816cfffffffffff, 0x816abffffffffff,
+		0x815f3ffffffffff, 0x817c7ffffffffff, 0x8168bffffffffff, 0x818cbffffffffff,
+		0x818cfffffffffff, 0x818afffffffffff, 0x8174fffffffffff, 0x8172fffffffffff,
+		0x8170fffffffffff, 0x816fbffffffffff, 0x81657ffffffffff, 0x816c7ffffffffff,
+		0x8186bffffffffff, 0x81763ffffffffff, 0x818a7ffffffffff, 0x8186fffffffffff,
+		0x81707ffffffffff, 0x8182bffffffffff, 0x818f3ffffffffff, 0x8182fffffffffff,
+	}
+
+	polygons, err := CellsToMultiPolygon(cells)
+	if err != nil {
+		t.Fatalf("CellsToMultiPolygon: %v", err)
+	}
+
+	if len(polygons) != 1 {
+		t.Fatalf("polygons: got %d, want 1", len(polygons))
+	}
+
+	if len(polygons[0].Holes) != 1 {
+		t.Fatalf("holes: got %d, want 1", len(polygons[0].Holes))
+	}
+}
diff --git a/x/h3go/paritytest/faces_test.go b/x/h3go/paritytest/faces_test.go
new file mode 100644
index 0000000..a73a942
--- /dev/null
+++ b/x/h3go/paritytest/faces_test.go
@@ -0,0 +1,56 @@
+/*
+ * Copyright 2026 Uber Technologies, Inc.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package paritytest
+
+import (
+	"sort"
+	"testing"
+
+	"github.com/uber/h3-go/v4/x/h3go"
+)
+
+// TestIcosahedronFacesMatchesCgo asserts Cell.IcosahedronFaces matches the cgo
+// reference across the corpus, as a sorted set.
+func TestIcosahedronFacesMatchesCgo(t *testing.T) {
+	t.Parallel()
+
+	for _, cell := range referenceCorpus(t) {
+		want, wantErr := cell.IcosahedronFaces()
+		got, gotErr := h3go.Cell(cell).IcosahedronFaces()
+
+		if !bothErr(wantErr, gotErr) {
+			t.Fatalf("IcosahedronFaces(%015x) error: cgo=%v h3go=%v", uint64(cell), wantErr, gotErr)
+		}
+
+		if wantErr != nil {
+			continue
+		}
+
+		sort.Ints(want)
+		sort.Ints(got)
+
+		if len(want) != len(got) {
+			t.Fatalf("IcosahedronFaces(%015x): cgo=%v h3go=%v", uint64(cell), want, got)
+		}
+
+		for i := range want {
+			if want[i] != got[i] {
+				t.Fatalf("IcosahedronFaces(%015x): cgo=%v h3go=%v", uint64(cell), want, got)
+			}
+		}
+	}
+}
diff --git a/x/h3go/paritytest/latlng_test.go b/x/h3go/paritytest/latlng_test.go
index 021ea74..5a6787f 100644
--- a/x/h3go/paritytest/latlng_test.go
+++ b/x/h3go/paritytest/latlng_test.go
@@ -147,3 +147,31 @@ func TestLatLngToCellErrorParity(t *testing.T) {
 		})
 	}
 }
+
+// TestLatLngHelpersMatchCgo checks NewLatLng, LatLng.Cell, LatLng.String, and
+// LatLngToCellString against the cgo reference.
+func TestLatLngHelpersMatchCgo(t *testing.T) {
+	t.Parallel()
+
+	for _, p := range corpusPoints {
+		for res := 0; res <= h3.MaxResolution; res++ {
+			wantCell, wantErr := h3.NewLatLng(p.Lat, p.Lng).Cell(res)
+			gotCell, gotErr := h3go.NewLatLng(p.Lat, p.Lng).Cell(res)
+
+			if !bothErr(wantErr, gotErr) || h3go.Cell(wantCell) != gotCell {
+				t.Fatalf("LatLng.Cell(%v,%d): cgo=%015x(%v) h3go=%015x(%v)", p, res, uint64(wantCell), wantErr, uint64(gotCell), gotErr)
+			}
+
+			wantStr, wantErr := h3.LatLngToCellString(p.Lat, p.Lng, res)
+			gotStr, gotErr := h3go.LatLngToCellString(p.Lat, p.Lng, res)
+
+			if !bothErr(wantErr, gotErr) || wantStr != gotStr {
+				t.Fatalf("LatLngToCellString(%v,%d): cgo=%q h3go=%q", p, res, wantStr, gotStr)
+			}
+		}
+
+		if want, got := h3.NewLatLng(p.Lat, p.Lng).String(), h3go.NewLatLng(p.Lat, p.Lng).String(); want != got {
+			t.Fatalf("LatLng.String(%v): cgo=%q h3go=%q", p, want, got)
+		}
+	}
+}
diff --git a/x/h3go/paritytest/multipoly_test.go b/x/h3go/paritytest/multipoly_test.go
new file mode 100644
index 0000000..45ef7e4
--- /dev/null
+++ b/x/h3go/paritytest/multipoly_test.go
@@ -0,0 +1,245 @@
+/*
+ * Copyright 2026 Uber Technologies, Inc.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package paritytest
+
+import (
+	"fmt"
+	"sort"
+	"testing"
+
+	"github.com/uber/h3-go/v4"
+	"github.com/uber/h3-go/v4/x/h3go"
+)
+
+// cgoLoopKey renders a cgo loop as a rounded, ordered vertex string for set
+// comparison.
+func cgoLoopKey(loop h3.GeoLoop) string {
+	parts := make([]string, len(loop))
+	for i, point := range loop {
+		parts[i] = fmt.Sprintf("%.7f,%.7f", point.Lat, point.Lng)
+	}
+
+	return fmt.Sprint(parts)
+}
+
+// goLoopKey renders a pure-Go loop the same way as cgoLoopKey.
+func goLoopKey(loop h3go.GeoLoop) string {
+	parts := make([]string, len(loop))
+	for i, point := range loop {
+		parts[i] = fmt.Sprintf("%.7f,%.7f", point.Lat, point.Lng)
+	}
+
+	return fmt.Sprint(parts)
+}
+
+// cgoLoopSet returns every loop (outer and holes) across all polygons, sorted.
+func cgoLoopSet(polygons []h3.GeoPolygon) []string {
+	var keys []string
+	for _, poly := range polygons {
+		keys = append(keys, cgoLoopKey(poly.GeoLoop))
+		for _, hole := range poly.Holes {
+			keys = append(keys, cgoLoopKey(hole))
+		}
+	}
+
+	sort.Strings(keys)
+
+	return keys
+}
+
+// goLoopSet returns every loop (outer and holes) across all polygons, sorted.
+func goLoopSet(polygons []h3go.GeoPolygon) []string {
+	var keys []string
+	for _, poly := range polygons {
+		keys = append(keys, goLoopKey(poly.GeoLoop))
+		for _, hole := range poly.Holes {
+			keys = append(keys, goLoopKey(hole))
+		}
+	}
+
+	sort.Strings(keys)
+
+	return keys
+}
+
+// holeCounts returns the per-polygon hole counts, sorted, for either result.
+func holeCounts[P any](polygons []P, count func(P) int) []int {
+	out := make([]int, len(polygons))
+	for i, poly := range polygons {
+		out[i] = count(poly)
+	}
+
+	sort.Ints(out)
+
+	return out
+}
+
+// assertSameMultiPolygon compares two multipolygon results: same polygon count,
+// same per-polygon hole counts, and the same set of loops.
+func assertSameMultiPolygon(t *testing.T, got []h3go.GeoPolygon, want []h3.GeoPolygon, msg string) {
+	t.Helper()
+
+	if len(got) != len(want) {
+		t.Fatalf("%s: polygon count cgo=%d h3go=%d", msg, len(want), len(got))
+	}
+
+	gotHoles := holeCounts(got, func(p h3go.GeoPolygon) int { return len(p.Holes) })
+	wantHoles := holeCounts(want, func(p h3.GeoPolygon) int { return len(p.Holes) })
+
+	if fmt.Sprint(gotHoles) != fmt.Sprint(wantHoles) {
+		t.Fatalf("%s: hole counts cgo=%v h3go=%v", msg, wantHoles, gotHoles)
+	}
+
+	gotLoops := goLoopSet(got)
+	wantLoops := cgoLoopSet(want)
+
+	if len(gotLoops) != len(wantLoops) {
+		t.Fatalf("%s: loop count cgo=%d h3go=%d", msg, len(wantLoops), len(gotLoops))
+	}
+
+	for i := range wantLoops {
+		if gotLoops[i] != wantLoops[i] {
+			t.Fatalf("%s: loop mismatch:\ncgo  %s\nh3go %s", msg, wantLoops[i], gotLoops[i])
+		}
+	}
+}
+
+// multiPolygonCellSets returns named cell sets covering single cells, contiguous
+// blobs, rings (which produce a hole), disjoint regions, and a pentagon area.
+func multiPolygonCellSets(t *testing.T) map[string][]h3.Cell {
+	t.Helper()
+
+	origin, err := h3.LatLngToCell(h3.LatLng{Lat: 37.78, Lng: -122.42}, 7)
+	if err != nil {
+		t.Fatalf("LatLngToCell: %v", err)
+	}
+
+	disk, err := origin.GridDisk(2)
+	if err != nil {
+		t.Fatalf("GridDisk: %v", err)
+	}
+
+	ring, err := origin.GridRing(2)
+	if err != nil {
+		t.Fatalf("GridRing: %v", err)
+	}
+
+	far, err := h3.LatLngToCell(h3.LatLng{Lat: -20, Lng: 30}, 7)
+	if err != nil {
+		t.Fatalf("LatLngToCell far: %v", err)
+	}
+
+	farDisk, err := far.GridDisk(1)
+	if err != nil {
+		t.Fatalf("GridDisk far: %v", err)
+	}
+
+	disjoint := append(append([]h3.Cell{}, disk...), farDisk...)
+
+	// A pentagon base cell and its neighborhood.
+	pentagon := h3.Cell(0x85080003fffffff)
+
+	pentDisk, err := pentagon.GridDisk(1)
+	if err != nil {
+		t.Fatalf("GridDisk pentagon: %v", err)
+	}
+
+	return map[string][]h3.Cell{
+		"single":   {origin},
+		"disk":     disk,
+		"ring":     ring,
+		"disjoint": disjoint,
+		"pentagon": pentDisk,
+	}
+}
+
+// TestCellsToMultiPolygonMatchesCgo asserts CellsToMultiPolygon matches the cgo
+// reference across several cell sets.
+func TestCellsToMultiPolygonMatchesCgo(t *testing.T) {
+	t.Parallel()
+
+	for name, cells := range multiPolygonCellSets(t) {
+		want, wantErr := h3.CellsToMultiPolygon(cells)
+		got, gotErr := h3go.CellsToMultiPolygon(toGoCells(cells))
+
+		if !bothErr(wantErr, gotErr) {
+			t.Fatalf("%s: error mismatch cgo=%v h3go=%v", name, wantErr, gotErr)
+		}
+
+		if wantErr != nil {
+			continue
+		}
+
+		assertSameMultiPolygon(t, got, want, name)
+	}
+}
+
+// TestCellsToMultiPolygonGlobe asserts the whole-globe case (all base cells)
+// matches the cgo reference.
+func TestCellsToMultiPolygonGlobe(t *testing.T) {
+	t.Parallel()
+
+	cells, err := h3.Res0Cells()
+	if err != nil {
+		t.Fatalf("Res0Cells: %v", err)
+	}
+
+	want, wantErr := h3.CellsToMultiPolygon(cells)
+	got, gotErr := h3go.CellsToMultiPolygon(toGoCells(cells))
+
+	if !bothErr(wantErr, gotErr) {
+		t.Fatalf("globe error mismatch cgo=%v h3go=%v", wantErr, gotErr)
+	}
+
+	assertSameMultiPolygon(t, got, want, "globe")
+}
+
+// TestCellsToMultiPolygonErrors asserts validation errors match the cgo
+// reference for empty, mismatched-resolution, duplicate, and invalid inputs.
+func TestCellsToMultiPolygonErrors(t *testing.T) {
+	t.Parallel()
+
+	origin, err := h3.LatLngToCell(h3.LatLng{Lat: 0, Lng: 0}, 7)
+	if err != nil {
+		t.Fatalf("LatLngToCell: %v", err)
+	}
+
+	parent, err := origin.Parent(6)
+	if err != nil {
+		t.Fatalf("Parent: %v", err)
+	}
+
+	cases := map[string][]h3.Cell{
+		"empty":     {},
+		"mixed_res": {origin, parent},
+		"duplicate": {origin, origin},
+		"invalid":   {h3.Cell(0)},
+	}
+
+	for name, cells := range cases {
+		want, wantErr := h3.CellsToMultiPolygon(cells)
+		got, gotErr := h3go.CellsToMultiPolygon(toGoCells(cells))
+
+		if !bothErr(wantErr, gotErr) {
+			t.Fatalf("%s: error mismatch cgo=%v h3go=%v", name, wantErr, gotErr)
+		}
+
+		if wantErr == nil && len(got) != len(want) {
+			t.Fatalf("%s: polygon count cgo=%d h3go=%d", name, len(want), len(got))
+		}
+	}
+}
diff --git a/x/h3go/paritytest/region_test.go b/x/h3go/paritytest/region_test.go
new file mode 100644
index 0000000..4c40b22
--- /dev/null
+++ b/x/h3go/paritytest/region_test.go
@@ -0,0 +1,213 @@
+/*
+ * Copyright 2026 Uber Technologies, Inc.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package paritytest
+
+import (
+	"testing"
+
+	"github.com/uber/h3-go/v4"
+	"github.com/uber/h3-go/v4/x/h3go"
+)
+
+// regionPolygons returns a set of named test polygons spanning simple squares,
+// a polygon with a hole, a transmeridian polygon, and a pentagon-spanning area.
+func regionPolygons() map[string]h3.GeoPolygon {
+	sfSquare := h3.GeoLoop{
+		{Lat: 37.813, Lng: -122.513},
+		{Lat: 37.813, Lng: -122.345},
+		{Lat: 37.700, Lng: -122.345},
+		{Lat: 37.700, Lng: -122.513},
+	}
+
+	return map[string]h3.GeoPolygon{
+		"sf_square": {GeoLoop: sfSquare},
+		"with_hole": {
+			GeoLoop: sfSquare,
+			Holes: []h3.GeoLoop{{
+				{Lat: 37.78, Lng: -122.45},
+				{Lat: 37.78, Lng: -122.42},
+				{Lat: 37.76, Lng: -122.42},
+				{Lat: 37.76, Lng: -122.45},
+			}},
+		},
+		"equator": {GeoLoop: h3.GeoLoop{
+			{Lat: 1, Lng: -1},
+			{Lat: 1, Lng: 1},
+			{Lat: -1, Lng: 1},
+			{Lat: -1, Lng: -1},
+		}},
+		"transmeridian": {GeoLoop: h3.GeoLoop{
+			{Lat: 10, Lng: 178},
+			{Lat: 10, Lng: -178},
+			{Lat: -10, Lng: -178},
+			{Lat: -10, Lng: 178},
+		}},
+		"near_pentagon": {GeoLoop: h3.GeoLoop{
+			{Lat: 64.7, Lng: 10.5},
+			{Lat: 64.7, Lng: 11.5},
+			{Lat: 63.7, Lng: 11.5},
+			{Lat: 63.7, Lng: 10.5},
+		}},
+	}
+}
+
+// toGoPolygon converts a cgo GeoPolygon to the pure-Go type.
+func toGoPolygon(polygon h3.GeoPolygon) h3go.GeoPolygon {
+	out := h3go.GeoPolygon{GeoLoop: toGoLoop(polygon.GeoLoop)}
+	for _, hole := range polygon.Holes {
+		out.Holes = append(out.Holes, toGoLoop(hole))
+	}
+
+	return out
+}
+
+// toGoLoop converts a cgo GeoLoop to the pure-Go type.
+func toGoLoop(loop h3.GeoLoop) h3go.GeoLoop {
+	out := make(h3go.GeoLoop, len(loop))
+	for i, point := range loop {
+		out[i] = h3go.LatLng{Lat: point.Lat, Lng: point.Lng}
+	}
+
+	return out
+}
+
+// TestPolygonToCellsMatchesCgo asserts PolygonToCells matches the cgo reference
+// as a set across several polygons and resolutions, including error parity.
+func TestPolygonToCellsMatchesCgo(t *testing.T) {
+	t.Parallel()
+
+	for name, polygon := range regionPolygons() {
+		for res := 4; res <= 8; res++ {
+			want, wantErr := h3.PolygonToCells(polygon, res)
+			got, gotErr := h3go.PolygonToCells(toGoPolygon(polygon), res)
+
+			if !bothErr(wantErr, gotErr) {
+				t.Fatalf("PolygonToCells(%s, %d) error mismatch: cgo=%v h3go=%v", name, res, wantErr, gotErr)
+			}
+
+			if wantErr != nil {
+				continue
+			}
+
+			assertSameCellSet(t, got, dropZeroCells(want), name)
+		}
+	}
+}
+
+// TestPolygonToCellsResolutionError asserts an out-of-range resolution fails the
+// same way in both implementations.
+func TestPolygonToCellsResolutionError(t *testing.T) {
+	t.Parallel()
+
+	polygon := regionPolygons()["sf_square"]
+
+	for _, res := range []int{-1, 16} {
+		_, wantErr := h3.PolygonToCells(polygon, res)
+		_, gotErr := h3go.PolygonToCells(toGoPolygon(polygon), res)
+
+		if !bothErr(wantErr, gotErr) {
+			t.Fatalf("PolygonToCells(res %d) error mismatch: cgo=%v h3go=%v", res, wantErr, gotErr)
+		}
+	}
+}
+
+// dropZeroCells removes the zero padding the cgo reference leaves in its output.
+func dropZeroCells(cells []h3.Cell) []h3.Cell {
+	out := cells[:0:0]
+	for _, cell := range cells {
+		if cell != 0 {
+			out = append(out, cell)
+		}
+	}
+
+	return out
+}
+
+// experimentalModes pairs each cgo containment mode with the pure-Go equivalent.
+func experimentalModes() []struct {
+	name string
+	cgo  h3.ContainmentMode
+	h3go h3go.ContainmentMode
+} {
+	return []struct {
+		name string
+		cgo  h3.ContainmentMode
+		h3go h3go.ContainmentMode
+	}{
+		{"center", h3.ContainmentCenter, h3go.ContainmentCenter},
+		{"full", h3.ContainmentFull, h3go.ContainmentFull},
+		{"overlapping", h3.ContainmentOverlapping, h3go.ContainmentOverlapping},
+		{"overlapping_bbox", h3.ContainmentOverlappingBbox, h3go.ContainmentOverlappingBbox},
+	}
+}
+
+// TestPolygonToCellsExperimentalMatchesCgo asserts PolygonToCellsExperimental
+// matches the cgo reference as a set across polygons, resolutions, and every
+// containment mode.
+func TestPolygonToCellsExperimentalMatchesCgo(t *testing.T) {
+	t.Parallel()
+
+	for name, polygon := range regionPolygons() {
+		for _, mode := range experimentalModes() {
+			for res := 4; res <= 7; res++ {
+				want, wantErr := h3.PolygonToCellsExperimental(polygon, res, mode.cgo)
+				got, gotErr := h3go.PolygonToCellsExperimental(toGoPolygon(polygon), res, mode.h3go)
+
+				if !bothErr(wantErr, gotErr) {
+					t.Fatalf("PolygonToCellsExperimental(%s, %s, %d) error mismatch: cgo=%v h3go=%v", name, mode.name, res, wantErr, gotErr)
+				}
+
+				if wantErr != nil {
+					continue
+				}
+
+				assertSameCellSet(t, got, dropZeroCells(want), name+"/"+mode.name)
+			}
+		}
+	}
+}
+
+// TestPolygonToCellsExperimentalErrors asserts resolution and mode errors match
+// the cgo reference, and that a tight cell cap surfaces the same bounds error.
+func TestPolygonToCellsExperimentalErrors(t *testing.T) {
+	t.Parallel()
+
+	polygon := regionPolygons()["sf_square"]
+
+	for _, res := range []int{-1, 16} {
+		_, wantErr := h3.PolygonToCellsExperimental(polygon, res, h3.ContainmentCenter)
+		_, gotErr := h3go.PolygonToCellsExperimental(toGoPolygon(polygon), res, h3go.ContainmentCenter)
+
+		if !bothErr(wantErr, gotErr) {
+			t.Fatalf("res %d error mismatch: cgo=%v h3go=%v", res, wantErr, gotErr)
+		}
+	}
+
+	_, wantErr := h3.PolygonToCellsExperimental(polygon, 6, h3.ContainmentInvalid)
+	_, gotErr := h3go.PolygonToCellsExperimental(toGoPolygon(polygon), 6, h3go.ContainmentInvalid)
+
+	if !bothErr(wantErr, gotErr) {
+		t.Fatalf("invalid mode error mismatch: cgo=%v h3go=%v", wantErr, gotErr)
+	}
+
+	_, wantErr = h3.PolygonToCellsExperimental(polygon, 7, h3.ContainmentCenter, 1)
+	_, gotErr = h3go.PolygonToCellsExperimental(toGoPolygon(polygon), 7, h3go.ContainmentCenter, 1)
+
+	if !bothErr(wantErr, gotErr) {
+		t.Fatalf("bounds error mismatch: cgo=%v h3go=%v", wantErr, gotErr)
+	}
+}
diff --git a/x/h3go/polyfill.go b/x/h3go/polyfill.go
new file mode 100644
index 0000000..79d25a9
--- /dev/null
+++ b/x/h3go/polyfill.go
@@ -0,0 +1,492 @@
+/*
+ * Copyright 2026 Uber Technologies, Inc.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package h3go
+
+import (
+	"iter"
+	"math"
+)
+
+// ContainmentMode selects which cells PolygonToCellsExperimental includes
+// relative to the polygon.
+type ContainmentMode uint32
+
+const (
+	// ContainmentCenter includes cells whose center is inside the polygon.
+	ContainmentCenter ContainmentMode = iota
+	// ContainmentFull includes cells that are fully contained by the polygon.
+	ContainmentFull
+	// ContainmentOverlapping includes cells that overlap the polygon at any point.
+	ContainmentOverlapping
+	// ContainmentOverlappingBbox includes cells whose bounding box overlaps the
+	// polygon.
+	ContainmentOverlappingBbox
+	// ContainmentInvalid is the first invalid mode value and must not be used.
+	ContainmentInvalid
+)
+
+// flagContainmentModeMask masks the containment-mode bits out of the flags word.
+const flagContainmentModeMask = 15
+
+// cellScaleFactor scales a cell's bounding box to be sure it covers the cell.
+// childScaleFactor scales it to cover all of the cell's finer-resolution
+// children. Both were chosen empirically.
+const (
+	cellScaleFactor  = 1.1
+	childScaleFactor = 1.4
+)
+
+// maxEdgeLengthRads is the maximum cell edge length, in radians, for each
+// resolution, taken at the center of each base cell at that resolution.
+var maxEdgeLengthRads = [MaxResolution + 1]float64{
+	0.21577206265130, 0.08308767068495, 0.03148970436439, 0.01190662871439,
+	0.00450053330908, 0.00170105523619, 0.00064293917678, 0.00024300820659,
+	0.00009184847087, 0.00003471545901, 0.00001312121017, 0.00000495935129,
+	0.00000187445860, 0.00000070847876, 0.00000026777980, 0.00000010121125,
+}
+
+// northPoleCells and southPoleCells list the cell containing each pole at every
+// resolution, used to expand a pole-covering bounding box to a full circle.
+var (
+	northPoleCells = [MaxResolution + 1]Cell{
+		0x8001fffffffffff, 0x81033ffffffffff, 0x820327fffffffff, 0x830326fffffffff,
+		0x8403263ffffffff, 0x85032623fffffff, 0x860326237ffffff, 0x870326233ffffff,
+		0x880326233bfffff, 0x890326233abffff, 0x8a0326233ab7fff, 0x8b0326233ab0fff,
+		0x8c0326233ab03ff, 0x8d0326233ab03bf, 0x8e0326233ab039f, 0x8f0326233ab0399,
+	}
+	southPoleCells = [MaxResolution + 1]Cell{
+		0x80f3fffffffffff, 0x81f2bffffffffff, 0x82f297fffffffff, 0x83f293fffffffff,
+		0x84f2939ffffffff, 0x85f29383fffffff, 0x86f29380fffffff, 0x87f29380effffff,
+		0x88f29380e1fffff, 0x89f29380e0fffff, 0x8af29380e0d7fff, 0x8bf29380e0d0fff,
+		0x8cf29380e0d0dff, 0x8df29380e0d0cff, 0x8ef29380e0d0cc7, 0x8ff29380e0d0cc4,
+	}
+)
+
+// res0BBoxesRads holds a precomputed bounding box, in radians, for every
+// resolution-0 base cell, indexed by base cell number.
+var res0BBoxesRads = [NumBaseCells]bbox{
+	{north: 1.52480158339146, south: 1.20305471830087, east: -0.60664883654036, west: 0.00568297271999},
+	{north: 1.52480158339146, south: 1.17872424267511, east: -0.60664883654036, west: 2.54046980298264},
+	{north: 1.52480158339146, south: 1.09069387298096, east: -2.85286053297673, west: 1.64310689027893},
+	{north: 1.41845302535151, south: 1.01285145697208, east: 0.00568297271999, west: -1.16770379632602},
+	{north: 1.27950477868453, south: 0.97226652536306, east: 0.55556064983494, west: -0.18229924845326},
+	{north: 1.32929586572429, south: 0.91898920750071, east: 2.05622344943192, west: 1.08813154278274},
+	{north: 1.32899086063916, south: 0.94271815376360, east: -2.29875289606378, west: 3.01700008041993},
+	{north: 1.26020983864103, south: 0.84291228415618, east: -0.89971867664861, west: -1.75967359310997},
+	{north: 1.21114673854945, south: 0.86170600921069, east: 1.19129757609455, west: 0.43777608996454},
+	{north: 1.21075831414294, south: 0.83795331049498, east: -1.72022875779891, west: -2.43793861727138},
+	{north: 1.15546530929588, south: 0.78982455384253, east: 2.53659412229266, west: 1.85709133451243},
+	{north: 1.15528445067052, south: 0.76641428724335, east: -3.06738507202411, west: 2.53646110244042},
+	{north: 1.10121643537669, south: 0.71330093663066, east: 0.09640581900154, west: -0.52154514518248},
+	{north: 1.07042472765165, south: 0.67603948819406, east: -0.47984202840088, west: -1.10306159603090},
+	{north: 1.03270228748960, south: 0.72356358827215, east: -2.24990138725146, west: -2.74510220919157},
+	{north: 1.01929924623886, south: 0.65491232835426, east: 0.63035574240731, west: 0.03537030096470},
+	{north: 1.01786037568858, south: 0.58827636737638, east: 1.53192721817065, west: 0.93672682511233},
+	{north: 0.98081434136020, south: 0.61076063532947, east: -2.67100636598529, west: 3.06516463008733},
+	{north: 0.98106023192774, south: 0.58679836571570, east: 2.02829766214461, west: 1.51334374970280},
+	{north: 0.96374551790056, south: 0.55186491737474, east: -1.42976721313659, west: -1.96852202530104},
+	{north: 0.87536136210723, south: 0.50008952762292, east: -1.92435613571430, west: -2.41641343219793},
+	{north: 0.88611243445554, south: 0.52742963716774, east: -0.95781946324194, west: -1.47628966305930},
+	{north: 0.86881343251986, south: 0.50770567021439, east: 1.03236795495839, west: 0.50347284027426},
+	{north: 0.89235638181782, south: 0.48781264892508, east: 2.76430302119150, west: 2.29989716697031},
+	{north: 0.82570569254601, south: 0.52173101741059, east: 2.30921681461428, west: 1.93198541828980},
+	{north: 0.80599330438546, south: 0.40150819579319, east: -3.06417559403240, west: 2.70079300784409},
+	{north: 0.81612079704781, south: 0.38396800633226, east: -0.21614378891839, west: -0.70420149722178},
+	{north: 0.75822779851431, south: 0.39943555383751, east: -2.34059978084699, west: -2.82127373822444},
+	{north: 0.78861390967531, south: 0.38742018303868, east: 0.23115687731652, west: -0.22599491086066},
+	{north: 0.71515840341957, south: 0.33012478438475, east: -0.64847976163163, west: -1.08249728121219},
+	{north: 0.70359051048414, south: 0.29148673180722, east: 1.71441081857246, west: 1.28443348381696},
+	{north: 0.69190629544818, south: 0.28808313184381, east: 0.64863909244647, west: 0.16372369282557},
+	{north: 0.64863235654749, south: 0.26290420067147, east: 2.10318098268379, west: 1.69556122548344},
+	{north: 0.65722892279906, south: 0.28222653310929, east: 1.30918693285466, west: 0.87594416271685},
+	{north: 0.64750997738584, south: 0.24149865709850, east: -1.30272192474556, west: -1.68708570163242},
+	{north: 0.62380174028378, south: 0.25522080363509, east: -2.72428423026826, west: 3.10401473237630},
+	{north: 0.64228460410023, south: 0.21206753429148, east: -1.67639240992071, west: -2.11772366767341},
+	{north: 0.59919175361146, south: 0.21620460836570, east: 2.48592868387690, west: 2.07350353893591},
+	{north: 0.55637406851384, south: 0.25276557437230, east: -0.99885388505694, west: -1.32642489358939},
+	{north: 0.55648013300665, south: 0.15187401321019, east: 2.87032088421324, west: 2.44642320475367},
+	{north: 0.54603687970450, south: 0.15589091511369, east: -2.06789866067060, west: -2.49091419631961},
+	{north: 0.51206347752697, south: 0.15522020377124, east: 0.95446767315996, west: 0.54443262110414},
+	{north: 0.49767951537101, south: 0.10944898890579, east: -0.04335162263358, west: -0.42900268178569},
+	{north: 0.46538045483671, south: 0.06029968637720, east: -0.41240613713421, west: -0.80603623808166},
+	{north: 0.44686891066946, south: 0.06926857458503, east: 0.32053284794952, west: -0.07005748900849},
+	{north: 0.43208958202064, south: 0.07796440938140, east: -3.06232453079660, west: 2.80602499990282},
+	{north: 0.43103892586713, south: 0.02927431919853, east: -2.41589238618422, west: -2.85735809951951},
+	{north: 0.38073727558986, south: -0.00297016159959, east: -0.77039553861218, west: -1.14788248745028},
+	{north: 0.39113816687141, south: -0.01518764903038, east: 1.49130246958290, west: 1.14714731736311},
+	{north: 0.33421063142418, south: 0.02526613430348, east: 1.15141032578749, west: 0.85000706261644},
+	{north: 0.38915669778582, south: -0.04371359825454, east: 1.88046353933242, west: 1.48230231380717},
+	{north: 0.33787520825987, south: -0.04835090128296, east: -1.12274014380603, west: -1.49454408844749},
+	{north: 0.33601418932337, south: -0.06675068178541, east: 2.23792354204464, west: 1.85723423013211},
+	{north: 0.31838318078049, south: -0.05821955623722, east: 0.66058854060373, west: 0.25452572938783},
+	{north: 0.33630761471457, south: -0.07589541031521, east: -1.47957331741818, west: -1.85981735718264},
+	{north: 0.28924817322870, south: -0.09150638064667, east: -1.83561930288569, west: -2.21855897384292},
+	{north: 0.26678632252475, south: -0.10058088990867, east: -2.76808651991421, west: 3.12792953247061},
+	{north: 0.29285254112587, south: -0.13483165093783, east: 2.61406468380434, west: 2.20466422911705},
+	{north: 0.20150342788824, south: -0.10279852729762, east: 0.06881896344365, west: -0.23925229432978},
+	{north: 0.21283813275258, south: -0.18626835417891, east: 2.93800440256577, west: 2.57470747655623},
+	{north: 0.19587614179884, south: -0.17237030304155, east: -2.16941795427335, west: -2.55405165906601},
+	{north: 0.17237030304155, south: -0.19587614179884, east: 0.97217469931645, west: 0.58754099452378},
+	{north: 0.18626835417891, south: -0.21283813275258, east: -0.20358825102402, west: -0.56688517703356},
+	{north: 0.10279852729762, south: -0.20150342788824, east: -3.07277369014614, west: 2.90234035926002},
+	{north: 0.13483165093783, south: -0.29285254112587, east: -0.52752796978545, west: -0.93692842447275},
+	{north: 0.10058088990867, south: -0.26678632252475, east: 0.37350613367558, west: -0.01366312111919},
+	{north: 0.09150638064667, south: -0.28924817322870, east: 1.30597335070410, west: 0.92303367974687},
+	{north: 0.07589541031521, south: -0.33630761471457, east: 1.66201933617161, west: 1.28177529640715},
+	{north: 0.05821955623722, south: -0.31838318078049, east: -2.48100411298606, west: -2.88706692420196},
+	{north: 0.06675068178541, south: -0.33601418932337, east: -0.90366911154516, west: -1.28435842345769},
+	{north: 0.04835090128296, south: -0.33787520825987, east: 2.01885250978376, west: 1.64704856514230},
+	{north: 0.04371359825454, south: -0.38915669778582, east: -1.26112911425737, west: -1.65929033978262},
+	{north: -0.02526613430348, south: -0.33421063142418, east: -1.99018232780231, west: -2.29158559097336},
+	{north: 0.01518764903038, south: -0.39113816687140, east: -1.65029018400690, west: -1.99444533622668},
+	{north: 0.00297016159959, south: -0.38073727558986, east: 2.37119711497761, west: 1.99371016613951},
+	{north: -0.02927431919853, south: -0.43103892586713, east: 0.72570026740558, west: 0.28423455407029},
+	{north: -0.07796440938140, south: -0.43208958202064, east: 0.07926812279319, west: -0.33556765368697},
+	{north: -0.06926857458503, south: -0.44686891066946, east: -2.82105980564027, west: 3.07153516458131},
+	{north: -0.06029968637720, south: -0.46538045483671, east: 2.72918651645558, west: 2.33555641550814},
+	{north: -0.10944898890579, south: -0.49767951537101, east: 3.09824103095621, west: 2.71258997180410},
+	{north: -0.15522020377124, south: -0.51206347752697, east: -2.18712498042983, west: -2.59716003248565},
+	{north: -0.15589091511369, south: -0.54603687970450, east: 1.07369399291919, west: 0.65067845727018},
+	{north: -0.15187401321019, south: -0.55648013300665, east: -0.27127176937655, west: -0.69516944883612},
+	{north: -0.25276557437230, south: -0.55637406851385, east: 2.14273876853285, west: 1.81516776000041},
+	{north: -0.21620460836570, south: -0.59919175361146, east: -0.65566396971290, west: -1.06808911465388},
+	{north: -0.21206753429148, south: -0.64228460410023, east: 1.46520024366909, west: 1.02386898591638},
+	{north: -0.25522080363509, south: -0.62380174028378, east: 0.41730842332153, west: -0.03757792121350},
+	{north: -0.24149865709850, south: -0.64750997738584, east: 1.83887072884423, west: 1.45450695195737},
+	{north: -0.28222653310929, south: -0.65722892279906, east: -1.83240572073513, west: -2.26564849087294},
+	{north: -0.26290420067147, south: -0.64863235654749, east: -1.03841167090601, west: -1.44603142810635},
+	{north: -0.28808313184381, south: -0.69190629544818, east: -2.49295356114332, west: -2.97786896076422},
+	{north: -0.29148673180722, south: -0.70359051048414, east: -1.42718183501734, west: -1.85715916977284},
+	{north: -0.33012478438475, south: -0.71515840341957, east: 2.49311289195816, west: 2.05909537237761},
+	{north: -0.38742018303868, south: -0.78861390967531, east: -2.91043577627328, west: 2.91559774272914},
+	{north: -0.39943555383751, south: -0.75822779851431, east: 0.80099287274280, west: 0.32031891536535},
+	{north: -0.38396800633226, south: -0.81612079704781, east: 2.92544886467140, west: 2.43739115636801},
+	{north: -0.40150819579319, south: -0.80599330438546, east: 0.07741705955739, west: -0.44079964574570},
+	{north: -0.52173101741059, south: -0.82570569254601, east: -0.83237583897551, west: -1.20960723529999},
+	{north: -0.48781264892508, south: -0.89235638181782, east: -0.37728963239830, west: -0.84169548661948},
+	{north: -0.50770567021439, south: -0.86881343251986, east: -2.10922469863141, west: -2.63811981331554},
+	{north: -0.52742963716774, south: -0.88611243445554, east: 2.18377319034785, west: 1.66530299053050},
+	{north: -0.50008952762292, south: -0.87536136210723, east: 1.21723651787549, west: 0.72517922139186},
+	{north: -0.55186491737474, south: -0.96374551790056, east: 1.71182544045320, west: 1.17307062828876},
+	{north: -0.58679836571570, south: -0.98106023192774, east: -1.11329499144518, west: -1.62824890388699},
+	{north: -0.61076063532947, south: -0.98081434136020, east: 0.47058628760450, west: -0.07642802350246},
+	{north: -0.58827636737638, south: -1.01786037568858, east: -1.60966543541914, west: -2.20486582847747},
+	{north: -0.65491232835426, south: -1.01929924623886, east: -2.51123691118248, west: -3.10622235262510},
+	{north: -0.72356358827215, south: -1.03270228748960, east: 0.89169126633833, west: 0.39649044439822},
+	{north: -0.67603948819406, south: -1.07042472765165, east: 2.66175062518892, west: 2.03853105755889},
+	{north: -0.71330093663066, south: -1.10121643537669, east: -3.04518683458825, west: 2.62004750840731},
+	{north: -0.76641428724335, south: -1.15528445067052, east: 0.07420758156568, west: -0.60513155114938},
+	{north: -0.78982455384253, south: -1.15546530929588, east: -0.60499853129713, west: -1.28450131907736},
+	{north: -0.83795331049498, south: -1.21075831414294, east: 1.42136389579088, west: 0.70365403631841},
+	{north: -0.86170600921069, south: -1.21114673854945, east: -1.95029507749525, west: -2.70381656362525},
+	{north: -0.84291228415618, south: -1.26020983864103, east: 2.24187397694118, west: 1.38191906047983},
+	{north: -0.94271815376360, south: -1.32899086063916, east: 0.84283975752601, west: -0.12459257316986},
+	{north: -0.91898920750071, south: -1.32929586572429, east: -1.08536920415787, west: -2.05346111080706},
+	{north: -0.97226652536306, south: -1.27950477868453, east: -2.58603200375485, west: 2.95929340513654},
+	{north: -1.01285145697208, south: -1.41845302535151, east: -3.13590968086981, west: 1.97388885726377},
+	{north: -1.09069387298096, south: -1.52480158339146, east: 0.28873212061306, west: -1.49848576331087},
+	{north: -1.17872424267511, south: -1.52480158339146, east: 2.53494381704943, west: -0.60112285060716},
+	{north: -1.20305471830087, south: -1.52480158339146, east: -0.60112285060716, west: 2.53494381704943},
+}
+
+// validRangeBBox is the full valid latitude/longitude domain in degrees. It
+// guards the first-vertex containment check from out-of-range coordinates.
+var validRangeBBox = bbox{north: halfPiDeg, south: -halfPiDeg, east: piDeg, west: -piDeg}
+
+// validateContainmentMode reports an error if the mode is out of range or sets
+// flag bits outside the containment-mode field.
+func validateContainmentMode(mode ContainmentMode) error {
+	flags := uint32(mode)
+	if flags&^uint32(flagContainmentModeMask) != 0 || flags&flagContainmentModeMask >= uint32(ContainmentInvalid) {
+		return ErrOptionInvalid
+	}
+
+	return nil
+}
+
+// baseCellNumToCell returns the resolution-0 cell for a base cell number, or the
+// zero cell if the number is out of range.
+func baseCellNumToCell(baseCellNum int) Cell {
+	if baseCellNum < 0 || baseCellNum >= NumBaseCells {
+		return 0
+	}
+
+	return setH3Index(0, baseCellNum, centerDigit)
+}
+
+// toDegrees returns the bounding box with each coordinate converted from radians
+// to degrees.
+func (b bbox) toDegrees() bbox {
+	return bbox{
+		north: b.north * RadsToDegs,
+		south: b.south * RadsToDegs,
+		east:  b.east * RadsToDegs,
+		west:  b.west * RadsToDegs,
+	}
+}
+
+// cellToBBox returns the bounding box of a cell in degrees. When coverChildren
+// is true the box is guaranteed to contain the cell's children at any finer
+// resolution. The box is approximate and may carry a significant margin.
+func cellToBBox(cell Cell, coverChildren bool) bbox {
+	res := cell.Resolution()
+
+	var out bbox
+	if res == 0 {
+		out = res0BBoxesRads[cell.BaseCellNumber()].toDegrees()
+	} else {
+		// cell is always valid here (it comes from the hierarchy walk), so the
+		// center projection cannot fail.
+		center, _ := cell.LatLng()
+		edge := maxEdgeLengthRads[res] * RadsToDegs
+		lngRatio := 1 / math.Cos(center.Lat*DegsToRads)
+		out = bbox{
+			north: center.Lat + edge,
+			south: center.Lat - edge,
+			east:  center.Lng + edge*lngRatio,
+			west:  center.Lng - edge*lngRatio,
+		}
+	}
+
+	// Scale the box, which also normalizes it to the lat/lng domain.
+	scale := cellScaleFactor
+	if coverChildren {
+		scale = childScaleFactor
+	}
+
+	out = out.scaled(scale)
+
+	if cell == northPoleCells[res] {
+		out.north = halfPiDeg
+	}
+
+	if cell == southPoleCells[res] {
+		out.south = -halfPiDeg
+	}
+
+	// A box covering a pole spans the full longitude domain, making it a circle
+	// around the pole.
+	if out.north == halfPiDeg || out.south == -halfPiDeg {
+		out.east = piDeg
+		out.west = -piDeg
+	}
+
+	return out
+}
+
+// nextCell returns the next cell to visit in the depth-first traversal of the
+// global cell hierarchy: the next sibling, ascending to a parent's next sibling
+// when the current cell is the last sibling, or the next base cell at the top.
+func nextCell(cell Cell) Cell {
+	res := cell.Resolution()
+	for {
+		if res == 0 {
+			return baseCellNumToCell(cell.BaseCellNumber() + 1)
+		}
+
+		parent := cell.setResolution(res-1).setIndexDigit(res, digitMask)
+
+		digit := cell.indexDigit(res)
+		if digit < invalidDigit-1 {
+			step := 1
+			// Skip the missing center child of a pentagon.
+			if parent.IsPentagon() && digit == centerDigit {
+				step = 2
+			}
+
+			return cell.setIndexDigit(res, digit+step)
+		}
+
+		res--
+		cell = parent
+	}
+}
+
+// targetCellInPolygon reports whether a cell at the target resolution should be
+// included for the given containment mode.
+func targetCellInPolygon(cell Cell, polygon GeoPolygon, bboxes []bbox, mode ContainmentMode) bool {
+	if mode == ContainmentCenter || mode == ContainmentOverlapping || mode == ContainmentOverlappingBbox {
+		center, _ := cell.LatLng()
+		if pointInsidePolygon(polygon, bboxes, center) {
+			return true
+		}
+	}
+
+	if mode == ContainmentOverlapping || mode == ContainmentOverlappingBbox {
+		// If the polygon is wholly contained by the cell, its first vertex maps
+		// to this cell. Guard against out-of-range coordinates first.
+		firstVertex := polygon.GeoLoop[0]
+		if validRangeBBox.contains(firstVertex) {
+			polygonCell, _ := LatLngToCell(firstVertex, cell.Resolution())
+			if polygonCell == cell {
+				return true
+			}
+		}
+	}
+
+	if mode == ContainmentFull || mode == ContainmentOverlapping || mode == ContainmentOverlappingBbox {
+		if targetCellBoundaryInPolygon(cell, polygon, bboxes, mode) {
+			return true
+		}
+	}
+
+	if mode == ContainmentOverlappingBbox {
+		return cellBBoxOverlapsPolygon(cell, polygon, bboxes)
+	}
+
+	return false
+}
+
+// targetCellBoundaryInPolygon checks containment or crossing of a target-res
+// cell's exact boundary against the polygon, per the containment mode.
+func targetCellBoundaryInPolygon(cell Cell, polygon GeoPolygon, bboxes []bbox, mode ContainmentMode) bool {
+	boundary, _ := cell.Boundary()
+	box := cellToBBox(cell, false)
+
+	if (mode == ContainmentFull || mode == ContainmentOverlappingBbox) &&
+		cellBoundaryInsidePolygon(polygon, bboxes, boundary, box) {
+		return true
+	}
+
+	// Center inclusion was already checked, so for overlap only the boundary
+	// crossing remains.
+	if (mode == ContainmentOverlapping || mode == ContainmentOverlappingBbox) &&
+		cellBoundaryCrossesPolygon(polygon, bboxes, boundary, box) {
+		return true
+	}
+
+	return false
+}
+
+// cellBBoxOverlapsPolygon reports whether a child-covering cell bounding box
+// overlaps the polygon, used by the overlapping-bbox mode.
+func cellBBoxOverlapsPolygon(cell Cell, polygon GeoPolygon, bboxes []bbox) bool {
+	box := cellToBBox(cell, true)
+	if !bboxes[0].overlaps(box) {
+		return false
+	}
+
+	boxBoundary := box.toCellBoundary()
+
+	return box.containsBBox(bboxes[0]) ||
+		pointInsidePolygon(polygon, bboxes, boxBoundary[0]) ||
+		cellBoundaryCrossesPolygon(polygon, bboxes, boxBoundary, box)
+}
+
+// coarseCellInPolygon reports whether a coarser-than-target cell is wholly
+// contained by the polygon (so all of its children are included). It returns the
+// containment result and whether the traversal should recurse into the children.
+func coarseCellInPolygon(cell Cell, polygon GeoPolygon, bboxes []bbox) (contained, recurse bool) {
+	box := cellToBBox(cell, true)
+	if !bboxes[0].overlaps(box) {
+		return false, false
+	}
+
+	if bboxes[0].containsBBox(box) {
+		boxBoundary := box.toCellBoundary()
+		if cellBoundaryInsidePolygon(polygon, bboxes, boxBoundary, box) {
+			return true, false
+		}
+	}
+
+	return false, true
+}
+
+// polygonCompactCells yields the compact set of cells covering the polygon: each
+// cell is either at the target resolution or a coarser cell whose every child is
+// contained. Inclusion at the target resolution follows the containment mode.
+func polygonCompactCells(polygon GeoPolygon, bboxes []bbox, res int, mode ContainmentMode) iter.Seq[Cell] {
+	return func(yield func(Cell) bool) {
+		cell := baseCellNumToCell(0)
+		for cell != 0 {
+			cellRes := cell.Resolution()
+
+			if cellRes == res {
+				if targetCellInPolygon(cell, polygon, bboxes, mode) && !yield(cell) {
+					return
+				}
+
+				cell = nextCell(cell)
+
+				continue
+			}
+
+			contained, recurse := coarseCellInPolygon(cell, polygon, bboxes)
+			if contained {
+				if !yield(cell) {
+					return
+				}
+
+				cell = nextCell(cell)
+
+				continue
+			}
+
+			if recurse {
+				// cell is coarser than the target, so a center child always exists.
+				cell, _ = cell.CenterChild(cellRes + 1)
+
+				continue
+			}
+
+			cell = nextCell(cell)
+		}
+	}
+}
+
+// PolygonToCellsExperimental fills the polygon with cells of the given
+// resolution, including a cell according to the containment mode. The optional
+// maxNumCellsReturn caps the number of cells; exceeding it returns
+// ErrMemoryBounds. Output ordering is not significant.
+func PolygonToCellsExperimental(polygon GeoPolygon, res int, mode ContainmentMode, maxNumCellsReturn ...int64) ([]Cell, error) {
+	maxNumCells := int64(math.MaxInt64)
+	if len(maxNumCellsReturn) > 0 {
+		maxNumCells = maxNumCellsReturn[0]
+	}
+
+	if len(polygon.GeoLoop) == 0 {
+		return nil, nil
+	}
+
+	if res < 0 || res > MaxResolution {
+		return nil, ErrResolutionDomain
+	}
+
+	if err := validateContainmentMode(mode); err != nil {
+		return nil, err
+	}
+
+	bboxes := bboxesFromGeoPolygon(polygon)
+
+	var (
+		out   []Cell
+		count int64
+	)
+
+	for compact := range polygonCompactCells(polygon, bboxes, res, mode) {
+		for child := range compact.childCells(res) {
+			if count >= maxNumCells {
+				return nil, ErrMemoryBounds
+			}
+
+			out = append(out, child)
+			count++
+		}
+	}
+
+	return out, nil
+}
diff --git a/x/h3go/polyfill_test.go b/x/h3go/polyfill_test.go
new file mode 100644
index 0000000..accacb8
--- /dev/null
+++ b/x/h3go/polyfill_test.go
@@ -0,0 +1,565 @@
+/*
+ * Copyright 2026 Uber Technologies, Inc.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package h3go
+
+import (
+	"errors"
+	"testing"
+)
+
+// TestPolygonToCellsExperimentalWithHole runs a holed polygon through the
+// boundary-based modes, exercising the hole branches of the boundary checks.
+func TestPolygonToCellsExperimentalWithHole(t *testing.T) {
+	t.Parallel()
+
+	hole := GeoLoop{
+		{Lat: 37.78, Lng: -122.45},
+		{Lat: 37.78, Lng: -122.42},
+		{Lat: 37.76, Lng: -122.42},
+		{Lat: 37.76, Lng: -122.45},
+	}
+	holed := GeoPolygon{GeoLoop: sfSquareLoop, Holes: []GeoLoop{hole}}
+	solid := GeoPolygon{GeoLoop: sfSquareLoop}
+
+	for _, mode := range []ContainmentMode{ContainmentFull, ContainmentOverlapping} {
+		withHole, err := PolygonToCellsExperimental(holed, 9, mode)
+		if err != nil {
+			t.Fatalf("mode %d holed: %v", mode, err)
+		}
+
+		without, err := PolygonToCellsExperimental(solid, 9, mode)
+		if err != nil {
+			t.Fatalf("mode %d solid: %v", mode, err)
+		}
+
+		if len(withHole) >= len(without) {
+			t.Fatalf("mode %d: hole did not reduce count: with=%d without=%d", mode, len(withHole), len(without))
+		}
+	}
+}
+
+// TestBBoxScaled covers the latitude clamps and the longitude wrap branches of
+// the scaling helper.
+func TestBBoxScaled(t *testing.T) {
+	t.Parallel()
+
+	tests := map[string]struct {
+		box   bbox
+		scale float64
+	}{
+		"latitude_clamps":  {bbox{north: 89, south: -89, east: 10, west: -10}, 1.4},
+		"east_west_wrap":   {bbox{north: 10, south: -10, east: 179, west: -179}, 1.1},
+		"west_over_pi":     {bbox{north: 10, south: -10, east: 195, west: 188}, 2},
+		"east_under_negpi": {bbox{north: 10, south: -10, east: -185, west: -190}, 2},
+	}
+
+	for name, tt := range tests {
+		t.Run(name, func(t *testing.T) {
+			t.Parallel()
+
+			got := tt.box.scaled(tt.scale)
+			if got.north > halfPiDeg || got.south < -halfPiDeg {
+				t.Fatalf("latitude out of domain: %+v", got)
+			}
+
+			if got.east > piDeg || got.east < -piDeg || got.west > piDeg || got.west < -piDeg {
+				t.Fatalf("longitude out of domain: %+v", got)
+			}
+		})
+	}
+}
+
+// TestBBoxNormalizationEastTrend covers the eastward-default normalization of a
+// standard box paired with a far-east transmeridian box.
+func TestBBoxNormalizationEastTrend(t *testing.T) {
+	t.Parallel()
+
+	standardEast := bbox{north: 5, south: -5, east: 160, west: 150}
+	trans := bbox{north: 10, south: -10, east: -170, west: 170}
+
+	firstNorm, secondNorm := normalizationFor(standardEast, trans)
+	if firstNorm != normalizeNone {
+		t.Fatalf("first normalization: got %d, want normalizeNone", firstNorm)
+	}
+
+	if secondNorm != normalizeEast {
+		t.Fatalf("second normalization: got %d, want normalizeEast", secondNorm)
+	}
+}
+
+// transmeridianLoop is a square straddling the antimeridian, used to exercise
+// the longitude-normalization branches.
+var transmeridianLoop = GeoLoop{
+	{Lat: 10, Lng: 178},
+	{Lat: 10, Lng: -178},
+	{Lat: -10, Lng: -178},
+	{Lat: -10, Lng: 178},
+}
+
+// TestPolygonToCellsExperimentalModes runs every containment mode and checks the
+// output is non-empty and, for center mode, that every cell center is inside.
+func TestPolygonToCellsExperimentalModes(t *testing.T) {
+	t.Parallel()
+
+	polygon := GeoPolygon{GeoLoop: sfSquareLoop}
+	bboxes := bboxesFromGeoPolygon(polygon)
+
+	modes := []ContainmentMode{
+		ContainmentCenter,
+		ContainmentFull,
+		ContainmentOverlapping,
+		ContainmentOverlappingBbox,
+	}
+
+	for _, mode := range modes {
+		cells, err := PolygonToCellsExperimental(polygon, 8, mode)
+		if err != nil {
+			t.Fatalf("mode %d: %v", mode, err)
+		}
+
+		if len(cells) == 0 {
+			t.Fatalf("mode %d: got no cells", mode)
+		}
+
+		if mode != ContainmentCenter {
+			continue
+		}
+
+		for _, cell := range cells {
+			center, err := cell.LatLng()
+			if err != nil {
+				t.Fatalf("LatLng: %v", err)
+			}
+
+			if !pointInsidePolygon(polygon, bboxes, center) {
+				t.Fatalf("center mode: cell %015x center not inside polygon", uint64(cell))
+			}
+		}
+	}
+}
+
+// TestPolygonToCellsExperimentalTransmeridian fills a polygon crossing the
+// antimeridian, exercising the longitude-normalization paths in the algorithm.
+func TestPolygonToCellsExperimentalTransmeridian(t *testing.T) {
+	t.Parallel()
+
+	polygon := GeoPolygon{GeoLoop: transmeridianLoop}
+
+	cells, err := PolygonToCellsExperimental(polygon, 4, ContainmentOverlappingBbox)
+	if err != nil {
+		t.Fatalf("PolygonToCellsExperimental: %v", err)
+	}
+
+	if len(cells) == 0 {
+		t.Fatal("got no cells for transmeridian polygon")
+	}
+}
+
+// TestPolygonToCellsExperimentalEmpty covers the empty-loop short circuit.
+func TestPolygonToCellsExperimentalEmpty(t *testing.T) {
+	t.Parallel()
+
+	cells, err := PolygonToCellsExperimental(GeoPolygon{}, 7, ContainmentCenter)
+	if err != nil || cells != nil {
+		t.Fatalf("empty: got %v (%v), want nil, nil", cells, err)
+	}
+}
+
+// TestPolygonToCellsExperimentalErrors covers the resolution and mode validation
+// paths.
+func TestPolygonToCellsExperimentalErrors(t *testing.T) {
+	t.Parallel()
+
+	polygon := GeoPolygon{GeoLoop: sfSquareLoop}
+
+	for _, res := range []int{-1, MaxResolution + 1} {
+		if _, err := PolygonToCellsExperimental(polygon, res, ContainmentCenter); !errors.Is(err, ErrResolutionDomain) {
+			t.Fatalf("res %d: got %v, want ErrResolutionDomain", res, err)
+		}
+	}
+
+	for _, mode := range []ContainmentMode{ContainmentInvalid, ContainmentMode(16)} {
+		if _, err := PolygonToCellsExperimental(polygon, 6, mode); !errors.Is(err, ErrOptionInvalid) {
+			t.Fatalf("mode %d: got %v, want ErrOptionInvalid", mode, err)
+		}
+	}
+}
+
+// TestPolygonToCellsExperimentalBounds covers the cell-cap path, where exceeding
+// the maximum returns ErrMemoryBounds.
+func TestPolygonToCellsExperimentalBounds(t *testing.T) {
+	t.Parallel()
+
+	polygon := GeoPolygon{GeoLoop: sfSquareLoop}
+	if _, err := PolygonToCellsExperimental(polygon, 8, ContainmentCenter, 1); !errors.Is(err, ErrMemoryBounds) {
+		t.Fatalf("bounds: got %v, want ErrMemoryBounds", err)
+	}
+}
+
+// TestPolygonCompactCellsEarlyStop stops iteration on the first coarse cell and
+// on the first target cell, covering both early-return paths of the iterator.
+func TestPolygonCompactCellsEarlyStop(t *testing.T) {
+	t.Parallel()
+
+	polygon := GeoPolygon{GeoLoop: sfSquareLoop}
+	bboxes := bboxesFromGeoPolygon(polygon)
+
+	res := 9
+
+	sawCoarse := false
+
+	for cell := range polygonCompactCells(polygon, bboxes, res, ContainmentFull) {
+		if cell.Resolution() < res {
+			sawCoarse = true
+
+			break
+		}
+	}
+
+	if !sawCoarse {
+		t.Fatal("expected at least one coarse compact cell")
+	}
+
+	sawTarget := false
+
+	for cell := range polygonCompactCells(polygon, bboxes, res, ContainmentFull) {
+		if cell.Resolution() == res {
+			sawTarget = true
+
+			break
+		}
+	}
+
+	if !sawTarget {
+		t.Fatal("expected at least one target-resolution compact cell")
+	}
+}
+
+// TestCellToBBoxPoles covers the pole-cell branches, where the bounding box is
+// expanded to a full circle around the pole.
+func TestCellToBBoxPoles(t *testing.T) {
+	t.Parallel()
+
+	for res := 0; res <= 3; res++ {
+		north := cellToBBox(northPoleCells[res], false)
+		if north.north != halfPiDeg || north.east != piDeg || north.west != -piDeg {
+			t.Fatalf("north pole res %d: got %+v", res, north)
+		}
+
+		south := cellToBBox(southPoleCells[res], true)
+		if south.south != -halfPiDeg || south.east != piDeg || south.west != -piDeg {
+			t.Fatalf("south pole res %d: got %+v", res, south)
+		}
+	}
+}
+
+// TestNextCellPentagonSkip covers the missing-pentagon-child skip in nextCell.
+func TestNextCellPentagonSkip(t *testing.T) {
+	t.Parallel()
+
+	// Base cell 4 is a pentagon; its center child has digit 0 at resolution 1.
+	centerChild := setH3Index(1, 4, centerDigit)
+
+	next := nextCell(centerChild)
+	if got := next.indexDigit(1); got != 2 {
+		t.Fatalf("nextCell pentagon center child: digit %d, want 2 (skipped the deleted 1)", got)
+	}
+}
+
+// TestBaseCellNumToCellRange covers the in-range and out-of-range cases.
+func TestBaseCellNumToCellRange(t *testing.T) {
+	t.Parallel()
+
+	if got := baseCellNumToCell(0); got.BaseCellNumber() != 0 || got.Resolution() != 0 {
+		t.Fatalf("baseCellNumToCell(0): got %015x", uint64(got))
+	}
+
+	for _, num := range []int{-1, NumBaseCells} {
+		if got := baseCellNumToCell(num); got != 0 {
+			t.Fatalf("baseCellNumToCell(%d): got %015x, want 0", num, uint64(got))
+		}
+	}
+}
+
+// TestValidateContainmentMode covers the valid, out-of-range, and extra-bit
+// cases.
+func TestValidateContainmentMode(t *testing.T) {
+	t.Parallel()
+
+	for _, mode := range []ContainmentMode{ContainmentCenter, ContainmentOverlappingBbox} {
+		if err := validateContainmentMode(mode); err != nil {
+			t.Fatalf("validateContainmentMode(%d): %v", mode, err)
+		}
+	}
+
+	for _, mode := range []ContainmentMode{ContainmentInvalid, ContainmentMode(16)} {
+		if err := validateContainmentMode(mode); !errors.Is(err, ErrOptionInvalid) {
+			t.Fatalf("validateContainmentMode(%d): got %v, want ErrOptionInvalid", mode, err)
+		}
+	}
+}
+
+// TestBBoxNormalization covers the longitude-normalization helper across the
+// standard and transmeridian box-pair combinations.
+func TestBBoxNormalization(t *testing.T) {
+	t.Parallel()
+
+	standardA := bbox{north: 10, south: -10, east: 10, west: -10}
+	standardB := bbox{north: 5, south: -5, east: 20, west: 5}
+	transEast := bbox{north: 10, south: -10, east: -170, west: 170}
+	transOther := bbox{north: 8, south: -8, east: -160, west: 165}
+
+	tests := map[string]struct {
+		first, second  bbox
+		wantFirst      longitudeNormalization
+		wantSecondNone bool
+	}{
+		"both_standard": {standardA, standardB, normalizeNone, true},
+		"first_trans":   {transEast, standardB, normalizeEast, true},
+		"both_trans":    {transEast, transOther, normalizeEast, false},
+	}
+
+	for name, tt := range tests {
+		t.Run(name, func(t *testing.T) {
+			t.Parallel()
+
+			firstNorm, secondNorm := normalizationFor(tt.first, tt.second)
+			if firstNorm != tt.wantFirst {
+				t.Fatalf("first normalization: got %d, want %d", firstNorm, tt.wantFirst)
+			}
+
+			if (secondNorm == normalizeNone) != tt.wantSecondNone {
+				t.Fatalf("second normalization: got %d", secondNorm)
+			}
+		})
+	}
+}
+
+// TestBBoxNormalizationWestTrend covers the westward-trending normalization of a
+// transmeridian box paired with a far-west standard box.
+func TestBBoxNormalizationWestTrend(t *testing.T) {
+	t.Parallel()
+
+	trans := bbox{north: 10, south: -10, east: -170, west: 170}
+	standardWest := bbox{north: 5, south: -5, east: -150, west: -160}
+
+	firstNorm, secondNorm := normalizationFor(trans, standardWest)
+	if firstNorm != normalizeWest {
+		t.Fatalf("first normalization: got %d, want normalizeWest", firstNorm)
+	}
+
+	// Standard second box needs no normalization.
+	if secondNorm != normalizeNone {
+		t.Fatalf("second normalization: got %d, want normalizeNone", secondNorm)
+	}
+
+	// Pairing a standard first box with a transmeridian second exercises the
+	// second-box branches.
+	firstNorm, secondNorm = normalizationFor(standardWest, trans)
+	if firstNorm != normalizeNone {
+		t.Fatalf("first normalization (swapped): got %d, want normalizeNone", firstNorm)
+	}
+
+	if secondNorm == normalizeNone {
+		t.Fatal("second normalization (swapped): got normalizeNone, want a shift")
+	}
+}
+
+// TestApplyNormalization covers each normalization case and its guard.
+func TestApplyNormalization(t *testing.T) {
+	t.Parallel()
+
+	if got := applyNormalization(-170, normalizeEast); got != 190 {
+		t.Fatalf("east of -170: got %v, want 190", got)
+	}
+
+	if got := applyNormalization(10, normalizeEast); got != 10 {
+		t.Fatalf("east of 10: got %v, want 10", got)
+	}
+
+	if got := applyNormalization(170, normalizeWest); got != -190 {
+		t.Fatalf("west of 170: got %v, want -190", got)
+	}
+
+	if got := applyNormalization(-10, normalizeWest); got != -10 {
+		t.Fatalf("west of -10: got %v, want -10", got)
+	}
+
+	if got := applyNormalization(42, normalizeNone); got != 42 {
+		t.Fatalf("none: got %v, want 42", got)
+	}
+}
+
+// TestBBoxOverlapAndContains covers the overlap and containment predicates,
+// including transmeridian normalization and the early-reject branches.
+func TestBBoxOverlapAndContains(t *testing.T) {
+	t.Parallel()
+
+	outer := bbox{north: 20, south: -20, east: 20, west: -20}
+	inner := bbox{north: 10, south: -10, east: 10, west: -10}
+	disjointLat := bbox{north: 40, south: 30, east: 10, west: -10}
+	disjointLng := bbox{north: 10, south: -10, east: 40, west: 30}
+
+	if !outer.overlaps(inner) {
+		t.Fatal("outer should overlap inner")
+	}
+
+	if outer.overlaps(disjointLat) {
+		t.Fatal("latitude-disjoint boxes should not overlap")
+	}
+
+	if outer.overlaps(disjointLng) {
+		t.Fatal("longitude-disjoint boxes should not overlap")
+	}
+
+	if !outer.containsBBox(inner) {
+		t.Fatal("outer should contain inner")
+	}
+
+	if inner.containsBBox(outer) {
+		t.Fatal("inner should not contain outer")
+	}
+
+	if outer.containsBBox(disjointLat) {
+		t.Fatal("outer should not contain a latitude-disjoint box")
+	}
+}
+
+// TestLineCrossesLine covers the crossing, parallel, and out-of-range cases.
+func TestLineCrossesLine(t *testing.T) {
+	t.Parallel()
+
+	a1 := LatLng{Lat: 0, Lng: 0}
+	a2 := LatLng{Lat: 0, Lng: 10}
+	b1 := LatLng{Lat: -5, Lng: 5}
+	b2 := LatLng{Lat: 5, Lng: 5}
+
+	if !lineCrossesLine(a1, a2, b1, b2) {
+		t.Fatal("crossing segments should report true")
+	}
+
+	// Parallel segments never intersect (zero denominator).
+	if lineCrossesLine(a1, a2, LatLng{Lat: 1, Lng: 0}, LatLng{Lat: 1, Lng: 10}) {
+		t.Fatal("parallel segments should report false")
+	}
+
+	// Segment b is to the side, so the intersection parameter is out of range.
+	if lineCrossesLine(a1, a2, LatLng{Lat: -5, Lng: 50}, LatLng{Lat: 5, Lng: 50}) {
+		t.Fatal("non-overlapping segments should report false")
+	}
+}
+
+// TestPolygonToCellsExperimentalSFCounts ports the
+// testPolygonToCellsExperimental.c regression: exact res-9 cell counts for the
+// San Francisco polygon across all four containment modes. The counts encode the
+// ordering Full < Center < Overlapping < OverlappingBbox.
+func TestPolygonToCellsExperimentalSFCounts(t *testing.T) {
+	t.Parallel()
+
+	sf := GeoPolygon{GeoLoop: radLoop([][2]float64{
+		{0.659966917655, -2.1364398519396}, {0.6595011102219, -2.1359434279405},
+		{0.6583348114025, -2.1354884206045}, {0.6581220034068, -2.1382437718946},
+		{0.6594479998527, -2.1384597563896}, {0.6599990002976, -2.1376771158464},
+	})}
+
+	tests := map[string]struct {
+		giveMode ContainmentMode
+		want     int
+	}{
+		"center":           {giveMode: ContainmentCenter, want: 1253},
+		"full":             {giveMode: ContainmentFull, want: 1175},
+		"overlapping":      {giveMode: ContainmentOverlapping, want: 1334},
+		"overlapping_bbox": {giveMode: ContainmentOverlappingBbox, want: 1416},
+	}
+
+	for name, tt := range tests {
+		t.Run(name, func(t *testing.T) {
+			t.Parallel()
+
+			cells, err := PolygonToCellsExperimental(sf, 9, tt.giveMode)
+			if err != nil {
+				t.Fatalf("PolygonToCellsExperimental: %v", err)
+			}
+
+			if len(cells) != tt.want {
+				t.Fatalf("got %d cells, want %d", len(cells), tt.want)
+			}
+		})
+	}
+}
+
+// TestCellToBBoxContainsGeometry ports the testCellToBBoxExhaustive.c
+// correctness properties: a cell's bounding box contains all of its own boundary
+// vertices, and a parent's child-covering bounding box contains every boundary
+// vertex of its descendants.
+func TestCellToBBoxContainsGeometry(t *testing.T) {
+	t.Parallel()
+
+	res0, err := Res0Cells()
+	if err != nil {
+		t.Fatalf("Res0Cells: %v", err)
+	}
+
+	containsAll := func(t *testing.T, box bbox, cell Cell) {
+		t.Helper()
+
+		boundary, err := cell.Boundary()
+		if err != nil {
+			t.Fatalf("Boundary(%015x): %v", uint64(cell), err)
+		}
+
+		for _, vertex := range boundary {
+			if !box.contains(vertex) {
+				t.Fatalf("bbox does not contain vertex %v of cell %015x", vertex, uint64(cell))
+			}
+		}
+	}
+
+	t.Run("cell_bbox_bounds_self", func(t *testing.T) {
+		t.Parallel()
+
+		for _, parent := range res0 {
+			for res := 0; res <= 2; res++ {
+				cells, err := parent.Children(res)
+				if err != nil {
+					t.Fatalf("Children(%d): %v", res, err)
+				}
+
+				for _, cell := range cells {
+					containsAll(t, cellToBBox(cell, false), cell)
+				}
+			}
+		}
+	})
+
+	t.Run("parent_bbox_bounds_children", func(t *testing.T) {
+		t.Parallel()
+
+		for _, parent := range res0 {
+			box := cellToBBox(parent, true)
+
+			children, err := parent.Children(2)
+			if err != nil {
+				t.Fatalf("Children(2): %v", err)
+			}
+
+			for _, child := range children {
+				containsAll(t, box, child)
+			}
+		}
+	})
+}
diff --git a/x/h3go/polygon.go b/x/h3go/polygon.go
new file mode 100644
index 0000000..3d7afac
--- /dev/null
+++ b/x/h3go/polygon.go
@@ -0,0 +1,221 @@
+/*
+ * Copyright 2026 Uber Technologies, Inc.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package h3go
+
+// dblEpsilon is the difference between 1 and the least value greater than 1 that
+// is representable as a float64. It nudges the ray-casting test off exact vertex
+// latitudes and longitudes, matching the H3 C library's use of DBL_EPSILON.
+const dblEpsilon = 2.2204460492503131e-16
+
+// normalizeLng shifts a negative longitude (in radians) east by a full turn when
+// the loop is transmeridian, so a crossing loop's coordinates compare in one
+// frame.
+func normalizeLng(lng float64, isTransmeridian bool) float64 {
+	if isTransmeridian && lng < 0 {
+		return lng + twoPiRad
+	}
+
+	return lng
+}
+
+// pointInsideGeoLoop reports whether the loop contains the coordinate, using a
+// ray-casting test. It fails fast when the point is outside the loop's bounding
+// box.
+//
+// The bounding-box test runs in degrees, the units of the public GeoLoop, but
+// the ray casting converts to radians first: the dblEpsilon vertex nudges below
+// only break ties correctly at radian coordinate scale, where dblEpsilon exceeds
+// half a ULP. At degree scale (a coordinate ~57x larger) the nudge rounds away,
+// which would drop cells whose center latitude exactly matches a polygon vertex
+// (see uber/h3#595).
+func pointInsideGeoLoop(loop GeoLoop, box bbox, coord LatLng) bool {
+	if !box.contains(coord) {
+		return false
+	}
+
+	isTransmeridian := box.isTransmeridian()
+	contains := false
+
+	lat := coord.Lat * DegsToRads
+	lng := normalizeLng(coord.Lng*DegsToRads, isTransmeridian)
+
+	for i := range loop {
+		a := LatLng{Lat: loop[i].Lat * DegsToRads, Lng: loop[i].Lng * DegsToRads}
+		next := loop[(i+1)%len(loop)]
+		b := LatLng{Lat: next.Lat * DegsToRads, Lng: next.Lng * DegsToRads}
+
+		// Ray casting requires the second point to be the higher one.
+		if a.Lat > b.Lat {
+			a, b = b, a
+		}
+
+		// Nudge north off an exact vertex latitude to avoid counting the ray
+		// crossing the same vertex twice on successive segments.
+		if lat == a.Lat || lat == b.Lat {
+			lat += dblEpsilon
+		}
+
+		// Skip segments the horizontal ray cannot reach.
+		if lat < a.Lat || lat > b.Lat {
+			continue
+		}
+
+		aLng := normalizeLng(a.Lng, isTransmeridian)
+		bLng := normalizeLng(b.Lng, isTransmeridian)
+
+		// Bias westerly on an exact longitude match to break ties consistently.
+		if aLng == lng || bLng == lng {
+			lng -= dblEpsilon
+		}
+
+		// Longitude of the segment at the point's latitude.
+		ratio := (lat - a.Lat) / (b.Lat - a.Lat)
+		testLng := normalizeLng(aLng+(bLng-aLng)*ratio, isTransmeridian)
+
+		if testLng > lng {
+			contains = !contains
+		}
+	}
+
+	return contains
+}
+
+// pointInsidePolygon reports whether the polygon contains the coordinate: inside
+// the outer loop and outside every hole. bboxes holds the outer loop's bounding
+// box followed by one per hole.
+func pointInsidePolygon(polygon GeoPolygon, bboxes []bbox, coord LatLng) bool {
+	contains := pointInsideGeoLoop(polygon.GeoLoop, bboxes[0], coord)
+	if !contains {
+		return false
+	}
+
+	for i := range polygon.Holes {
+		if pointInsideGeoLoop(polygon.Holes[i], bboxes[i+1], coord) {
+			return false
+		}
+	}
+
+	return true
+}
+
+// lineCrossesLine reports whether segment a1→a2 intersects segment b1→b2. This
+// is a purely Cartesian test that ignores antimeridian wrapping and poles.
+func lineCrossesLine(a1, a2, b1, b2 LatLng) bool {
+	denom := (b2.Lng-b1.Lng)*(a2.Lat-a1.Lat) - (b2.Lat-b1.Lat)*(a2.Lng-a1.Lng)
+	if denom == 0 {
+		return false
+	}
+
+	test := ((b2.Lat-b1.Lat)*(a1.Lng-b1.Lng) - (b2.Lng-b1.Lng)*(a1.Lat-b1.Lat)) / denom
+	if test < 0 || test > 1 {
+		return false
+	}
+
+	test = ((a2.Lat-a1.Lat)*(a1.Lng-b1.Lng) - (a2.Lng-a1.Lng)*(a1.Lat-b1.Lat)) / denom
+
+	return test >= 0 && test <= 1
+}
+
+// cellBoundaryCrossesGeoLoop reports whether any segment of the cell boundary
+// intersects any segment of the loop. Crossing means line-segment intersection;
+// it does not include containment.
+func cellBoundaryCrossesGeoLoop(loop GeoLoop, loopBBox bbox, boundary CellBoundary, boundaryBBox bbox) bool {
+	if !loopBBox.overlaps(boundaryBBox) {
+		return false
+	}
+
+	loopNorm, boundaryNorm := normalizationFor(loopBBox, boundaryBBox)
+
+	normalBoundary := make([]LatLng, len(boundary))
+	for i := range boundary {
+		normalBoundary[i] = LatLng{Lat: boundary[i].Lat, Lng: applyNormalization(boundary[i].Lng, boundaryNorm)}
+	}
+
+	normalBoundaryBBox := bbox{
+		north: boundaryBBox.north,
+		south: boundaryBBox.south,
+		east:  applyNormalization(boundaryBBox.east, boundaryNorm),
+		west:  applyNormalization(boundaryBBox.west, boundaryNorm),
+	}
+
+	for i := range loop {
+		loop1 := LatLng{Lat: loop[i].Lat, Lng: applyNormalization(loop[i].Lng, loopNorm)}
+		next := loop[(i+1)%len(loop)]
+		loop2 := LatLng{Lat: next.Lat, Lng: applyNormalization(next.Lng, loopNorm)}
+
+		// Skip segments that cannot reach the boundary's bounding box.
+		if (loop1.Lat >= normalBoundaryBBox.north && loop2.Lat >= normalBoundaryBBox.north) ||
+			(loop1.Lat <= normalBoundaryBBox.south && loop2.Lat <= normalBoundaryBBox.south) ||
+			(loop1.Lng <= normalBoundaryBBox.west && loop2.Lng <= normalBoundaryBBox.west) ||
+			(loop1.Lng >= normalBoundaryBBox.east && loop2.Lng >= normalBoundaryBBox.east) {
+			continue
+		}
+
+		for j := range normalBoundary {
+			other := normalBoundary[(j+1)%len(normalBoundary)]
+			if lineCrossesLine(loop1, loop2, normalBoundary[j], other) {
+				return true
+			}
+		}
+	}
+
+	return false
+}
+
+// cellBoundaryInsidePolygon reports whether the cell boundary is completely
+// contained by the polygon: its first vertex is inside, it crosses neither the
+// outer loop nor any hole, and it contains no hole.
+func cellBoundaryInsidePolygon(polygon GeoPolygon, bboxes []bbox, boundary CellBoundary, boundaryBBox bbox) bool {
+	// Fails fast when the first vertex is outside the bounding box.
+	if !pointInsidePolygon(polygon, bboxes, boundary[0]) {
+		return false
+	}
+
+	if cellBoundaryCrossesGeoLoop(polygon.GeoLoop, bboxes[0], boundary, boundaryBBox) {
+		return false
+	}
+
+	boundaryLoop := GeoLoop(boundary)
+
+	for i := range polygon.Holes {
+		hole := polygon.Holes[i]
+		if len(hole) > 0 &&
+			(pointInsideGeoLoop(boundaryLoop, boundaryBBox, hole[0]) ||
+				cellBoundaryCrossesGeoLoop(hole, bboxes[i+1], boundary, boundaryBBox)) {
+			return false
+		}
+	}
+
+	return true
+}
+
+// cellBoundaryCrossesPolygon reports whether any part of the cell boundary
+// crosses the polygon's outer loop or any hole. Crossing means line-segment
+// intersection; it does not include containment.
+func cellBoundaryCrossesPolygon(polygon GeoPolygon, bboxes []bbox, boundary CellBoundary, boundaryBBox bbox) bool {
+	if cellBoundaryCrossesGeoLoop(polygon.GeoLoop, bboxes[0], boundary, boundaryBBox) {
+		return true
+	}
+
+	for i := range polygon.Holes {
+		if cellBoundaryCrossesGeoLoop(polygon.Holes[i], bboxes[i+1], boundary, boundaryBBox) {
+			return true
+		}
+	}
+
+	return false
+}
diff --git a/x/h3go/region.go b/x/h3go/region.go
new file mode 100644
index 0000000..57fdb9f
--- /dev/null
+++ b/x/h3go/region.go
@@ -0,0 +1,171 @@
+/*
+ * Copyright 2026 Uber Technologies, Inc.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package h3go
+
+// GeoLoop is an ordered list of geographic coordinates in degrees describing a
+// closed loop; the final point is implicitly connected back to the first.
+type GeoLoop []LatLng
+
+// GeoPolygon is a GeoLoop outer boundary with zero or more GeoLoop holes.
+type GeoPolygon struct {
+	GeoLoop GeoLoop
+	Holes   []GeoLoop
+}
+
+// polygonToCellsBuffer is added to the estimated cell count to cover small
+// polygons near icosahedron edges at odd resolutions, where line tracing needs a
+// little more room than the estimator provides.
+const polygonToCellsBuffer = 12
+
+// maxPolygonToCellsSize returns an upper bound on the number of cells that
+// PolygonToCells may produce for the polygon at the given resolution. It is the
+// larger of the bounding-box cell estimate and the total vertex count, plus a
+// small buffer.
+func maxPolygonToCellsSize(polygon GeoPolygon, res int) (int, error) {
+	numHexagons, err := bboxHexEstimate(bboxFromGeoLoop(polygon.GeoLoop), res)
+	if err != nil {
+		return 0, err
+	}
+
+	// The estimate usually exceeds the vertex count, but guard the rare case it
+	// does not.
+	totalVerts := len(polygon.GeoLoop)
+	for i := range polygon.Holes {
+		totalVerts += len(polygon.Holes[i])
+	}
+
+	numHexagons = max(numHexagons, totalVerts)
+
+	return numHexagons + polygonToCellsBuffer, nil
+}
+
+// getEdgeHexagons traces a loop with cells of the given resolution, adding every
+// cell whose center the loop passes through to the search set. These seed the
+// flood fill in PolygonToCells.
+func getEdgeHexagons(loop GeoLoop, res int, search map[Cell]bool) error {
+	for i := range loop {
+		origin := loop[i]
+		destination := loop[(i+1)%len(loop)]
+
+		numHexes, err := lineHexEstimate(origin, destination, res)
+		if err != nil {
+			return err
+		}
+
+		invNumHexes := 1.0 / float64(numHexes)
+		for j := range numHexes {
+			interpolate := LatLng{
+				Lat: origin.Lat*float64(numHexes-j)*invNumHexes + destination.Lat*float64(j)*invNumHexes,
+				Lng: origin.Lng*float64(numHexes-j)*invNumHexes + destination.Lng*float64(j)*invNumHexes,
+			}
+
+			// res and finiteness are already validated by lineHexEstimate above,
+			// so this conversion cannot fail here.
+			cell, _ := LatLngToCell(interpolate, res)
+			search[cell] = true
+		}
+	}
+
+	return nil
+}
+
+// PolygonToCells returns the cells of the given resolution whose centers fall
+// within the polygon. The polygon is considered in Cartesian (lat/lng) space:
+// the outer loop minus any holes. Output ordering is not significant.
+//
+// The algorithm traces the loops with cells, then flood-fills outward from those
+// seeds, keeping every neighbor whose center is contained, until no new cells are
+// found. This means two adjacent polygons with no overlap produce disjoint cell
+// sets.
+func PolygonToCells(polygon GeoPolygon, res int) ([]Cell, error) {
+	if len(polygon.GeoLoop) == 0 {
+		return nil, nil
+	}
+
+	bboxes := bboxesFromGeoPolygon(polygon)
+
+	// 1. Trace the outer loop and any holes to seed the search set. Tracing the
+	// first loop surfaces an invalid resolution.
+	search := make(map[Cell]bool)
+	for _, loop := range append([]GeoLoop{polygon.GeoLoop}, polygon.Holes...) {
+		if err := getEdgeHexagons(loop, res, search); err != nil {
+			return nil, err
+		}
+	}
+
+	// A degenerate bounding box (zero width or height) cannot be filled.
+	sizeHint, err := maxPolygonToCellsSize(polygon, res)
+	if err != nil {
+		return nil, err
+	}
+
+	// 2. Flood fill: from each search cell, test it and its neighbors for
+	// containment, and use the newly contained cells as the next search set.
+	found := make(map[Cell]bool, sizeHint)
+
+	searchList := make([]Cell, 0, len(search))
+	for cell := range search {
+		searchList = append(searchList, cell)
+	}
+
+	for len(searchList) > 0 {
+		searchList = polygonFloodStep(polygon, bboxes, searchList, found)
+	}
+
+	out := make([]Cell, 0, len(found))
+	for cell := range found {
+		out = append(out, cell)
+	}
+
+	return out, nil
+}
+
+// polygonFloodStep expands one generation of the polygon flood fill: for each
+// cell in searchList, it tests the cell and its neighbors and records those whose
+// center is inside the polygon, returning the newly found cells. Every cell here
+// comes from a prior LatLngToCell or grid-disk step, so it is always valid and
+// the projection calls cannot fail.
+func polygonFloodStep(polygon GeoPolygon, bboxes []bbox, searchList []Cell, found map[Cell]bool) []Cell {
+	var nextSearch []Cell
+
+	for _, searchHex := range searchList {
+		ring, _ := searchHex.GridDisk(1)
+
+		for _, hex := range ring {
+			if found[hex] {
+				continue
+			}
+
+			center, _ := hex.LatLng()
+			if !pointInsidePolygon(polygon, bboxes, center) {
+				continue
+			}
+
+			found[hex] = true
+
+			nextSearch = append(nextSearch, hex)
+		}
+	}
+
+	return nextSearch
+}
+
+// Cells returns the cells of the given resolution whose centers fall within the
+// polygon.
+func (p GeoPolygon) Cells(res int) ([]Cell, error) {
+	return PolygonToCells(p, res)
+}
diff --git a/x/h3go/region_test.go b/x/h3go/region_test.go
new file mode 100644
index 0000000..1e6df5d
--- /dev/null
+++ b/x/h3go/region_test.go
@@ -0,0 +1,565 @@
+/*
+ * Copyright 2026 Uber Technologies, Inc.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package h3go
+
+import (
+	"errors"
+	"math"
+	"testing"
+)
+
+// sfSquareLoop is a small square over San Francisco used across region tests.
+var sfSquareLoop = GeoLoop{
+	{Lat: 37.813, Lng: -122.513},
+	{Lat: 37.813, Lng: -122.345},
+	{Lat: 37.700, Lng: -122.345},
+	{Lat: 37.700, Lng: -122.513},
+}
+
+// radLoop converts a loop given in radians (as the C test fixtures are) into a
+// degrees GeoLoop.
+func radLoop(radVerts [][2]float64) GeoLoop {
+	loop := make(GeoLoop, len(radVerts))
+	for i, vert := range radVerts {
+		loop[i] = LatLng{Lat: vert[0] * RadsToDegs, Lng: vert[1] * RadsToDegs}
+	}
+
+	return loop
+}
+
+// TestPolygonToCellsReported ports the testPolygonToCellsReported.c regression
+// cases: the entire world split into two polygons, and several real-world
+// polygons with exact expected cell counts.
+func TestPolygonToCellsReported(t *testing.T) {
+	t.Parallel()
+
+	t.Run("entire_world", func(t *testing.T) {
+		t.Parallel()
+
+		world1 := GeoPolygon{GeoLoop: GeoLoop{
+			{Lat: -90, Lng: -180}, {Lat: 90, Lng: -180},
+			{Lat: 90, Lng: 0}, {Lat: -90, Lng: 0},
+		}}
+		world2 := GeoPolygon{GeoLoop: GeoLoop{
+			{Lat: -90, Lng: 0}, {Lat: 90, Lng: 0},
+			{Lat: 90, Lng: 180}, {Lat: -90, Lng: 180},
+		}}
+
+		for res := range 3 {
+			cells1, err := PolygonToCells(world1, res)
+			if err != nil {
+				t.Fatalf("PolygonToCells(world1, %d): %v", res, err)
+			}
+
+			cells2, err := PolygonToCells(world2, res)
+			if err != nil {
+				t.Fatalf("PolygonToCells(world2, %d): %v", res, err)
+			}
+
+			if got, want := len(cells1)+len(cells2), NumCells(res); got != want {
+				t.Fatalf("res %d: got %d cells, want %d (entire world)", res, got, want)
+			}
+
+			seen := make(map[Cell]bool, len(cells1))
+			for _, cell := range cells1 {
+				seen[cell] = true
+			}
+
+			for _, cell := range cells2 {
+				if seen[cell] {
+					t.Fatalf("res %d: cell %015x found in both halves", res, uint64(cell))
+				}
+			}
+		}
+	})
+
+	t.Run("exact_counts", func(t *testing.T) {
+		t.Parallel()
+
+		// https://github.com/uber/h3/issues/595: a vertex due east of the center
+		// at exactly the same latitude.
+		center595, err := Cell(0x85283473fffffff).LatLng()
+		if err != nil {
+			t.Fatalf("center LatLng: %v", err)
+		}
+
+		tests := map[string]struct {
+			giveLoop GeoLoop
+			giveRes  int
+			want     int
+		}{
+			"h3js_67": {
+				giveLoop: GeoLoop{
+					{Lat: -33.13755119234615, Lng: -56.25},
+					{Lat: -34.30714385628804, Lng: -56.25},
+					{Lat: -34.30714385628804, Lng: -57.65625},
+					{Lat: -33.13755119234615, Lng: -57.65625},
+				},
+				giveRes: 7,
+				want:    4499,
+			},
+			"h3js_67_2nd": {
+				giveLoop: GeoLoop{
+					{Lat: -34.30714385628804, Lng: -57.65625},
+					{Lat: -35.4606699514953, Lng: -57.65625},
+					{Lat: -35.4606699514953, Lng: -59.0625},
+					{Lat: -34.30714385628804, Lng: -59.0625},
+				},
+				giveRes: 7,
+				want:    4609,
+			},
+			"h3_136": {
+				giveLoop: radLoop([][2]float64{
+					{0.10068990369902957, 0.8920772174196191},
+					{0.10032914690616246, 0.8915914753447348},
+					{0.10033349237998787, 0.8915860128746426},
+					{0.10069496685903621, 0.8920742194546231},
+				}),
+				giveRes: 13,
+				want:    4353,
+			},
+			"issue_595": {
+				giveLoop: radLoop([][2]float64{
+					{center595.Lat * DegsToRads, -2.121207808248113},
+					{0.6565301558937859, -2.1281107217935986},
+					{0.6515463604919347, -2.1345342663428695},
+					{0.6466583305904194, -2.1276313527973842},
+				}),
+				giveRes: 5,
+				want:    8,
+			},
+		}
+
+		for name, tt := range tests {
+			t.Run(name, func(t *testing.T) {
+				t.Parallel()
+
+				cells, err := PolygonToCells(GeoPolygon{GeoLoop: tt.giveLoop}, tt.giveRes)
+				if err != nil {
+					t.Fatalf("PolygonToCells: %v", err)
+				}
+
+				if len(cells) != tt.want {
+					t.Fatalf("got %d cells, want %d", len(cells), tt.want)
+				}
+			})
+		}
+	})
+}
+
+// TestPolygonToCellsBasic checks a simple polygon fills with contained cells
+// whose centers are all inside, and that the method form agrees.
+func TestPolygonToCellsBasic(t *testing.T) {
+	t.Parallel()
+
+	polygon := GeoPolygon{GeoLoop: sfSquareLoop}
+
+	cells, err := PolygonToCells(polygon, 7)
+	if err != nil {
+		t.Fatalf("PolygonToCells: %v", err)
+	}
+
+	if len(cells) == 0 {
+		t.Fatal("PolygonToCells: got no cells, want some")
+	}
+
+	bboxes := bboxesFromGeoPolygon(polygon)
+
+	for _, cell := range cells {
+		center, err := cell.LatLng()
+		if err != nil {
+			t.Fatalf("LatLng: %v", err)
+		}
+
+		if !pointInsidePolygon(polygon, bboxes, center) {
+			t.Fatalf("cell %015x center not inside polygon", uint64(cell))
+		}
+	}
+
+	viaMethod, err := polygon.Cells(7)
+	if err != nil || len(viaMethod) != len(cells) {
+		t.Fatalf("Cells method: got %d (%v), want %d", len(viaMethod), err, len(cells))
+	}
+}
+
+// TestPolygonToCellsWithHole checks that cells whose centers are inside a hole
+// are excluded.
+func TestPolygonToCellsWithHole(t *testing.T) {
+	t.Parallel()
+
+	hole := GeoLoop{
+		{Lat: 37.78, Lng: -122.45},
+		{Lat: 37.78, Lng: -122.42},
+		{Lat: 37.76, Lng: -122.42},
+		{Lat: 37.76, Lng: -122.45},
+	}
+	polygon := GeoPolygon{GeoLoop: sfSquareLoop, Holes: []GeoLoop{hole}}
+
+	withHole, err := PolygonToCells(polygon, 8)
+	if err != nil {
+		t.Fatalf("PolygonToCells(with hole): %v", err)
+	}
+
+	withoutHole, err := PolygonToCells(GeoPolygon{GeoLoop: sfSquareLoop}, 8)
+	if err != nil {
+		t.Fatalf("PolygonToCells(no hole): %v", err)
+	}
+
+	if len(withHole) >= len(withoutHole) {
+		t.Fatalf("hole did not reduce cell count: with=%d without=%d", len(withHole), len(withoutHole))
+	}
+}
+
+// TestPolygonToCellsEmpty covers the empty-loop short circuit.
+func TestPolygonToCellsEmpty(t *testing.T) {
+	t.Parallel()
+
+	cells, err := PolygonToCells(GeoPolygon{}, 7)
+	if err != nil || cells != nil {
+		t.Fatalf("PolygonToCells(empty): got %v (%v), want nil, nil", cells, err)
+	}
+}
+
+// TestPolygonToCellsResolutionError covers the invalid-resolution path, surfaced
+// while tracing the loop edges.
+func TestPolygonToCellsResolutionError(t *testing.T) {
+	t.Parallel()
+
+	for _, res := range []int{-1, MaxResolution + 1} {
+		if _, err := PolygonToCells(GeoPolygon{GeoLoop: sfSquareLoop}, res); err == nil {
+			t.Fatalf("PolygonToCells(res %d): got nil error, want failure", res)
+		}
+	}
+}
+
+// TestPolygonToCellsDegenerate covers the degenerate bounding box path, where the
+// polygon has zero height (all vertices on one parallel).
+func TestPolygonToCellsDegenerate(t *testing.T) {
+	t.Parallel()
+
+	flat := GeoPolygon{GeoLoop: GeoLoop{
+		{Lat: 1, Lng: -1},
+		{Lat: 1, Lng: 0},
+		{Lat: 1, Lng: 1},
+	}}
+
+	if _, err := PolygonToCells(flat, 7); !errors.Is(err, ErrFailed) {
+		t.Fatalf("PolygonToCells(degenerate): got %v, want ErrFailed", err)
+	}
+}
+
+// TestGetEdgeHexagonsError covers the resolution error path of the loop tracer.
+func TestGetEdgeHexagonsError(t *testing.T) {
+	t.Parallel()
+
+	search := make(map[Cell]bool)
+	if err := getEdgeHexagons(sfSquareLoop, -1, search); err == nil {
+		t.Fatal("getEdgeHexagons(res -1): got nil error, want failure")
+	}
+}
+
+// TestPolygonFloodStepValidCells confirms the flood step expands a valid search
+// cell into contained neighbors.
+func TestPolygonFloodStepValidCells(t *testing.T) {
+	t.Parallel()
+
+	polygon := GeoPolygon{GeoLoop: sfSquareLoop}
+	bboxes := bboxesFromGeoPolygon(polygon)
+
+	center := LatLng{Lat: 37.76, Lng: -122.43}
+
+	seed, err := LatLngToCell(center, 8)
+	if err != nil {
+		t.Fatalf("LatLngToCell: %v", err)
+	}
+
+	found := map[Cell]bool{}
+
+	next := polygonFloodStep(polygon, bboxes, []Cell{seed}, found)
+	if len(found) == 0 || len(next) == 0 {
+		t.Fatalf("polygonFloodStep found %d cells, next %d; want some", len(found), len(next))
+	}
+}
+
+// TestMaxPolygonToCellsSizeVertexFloor covers the branch where the vertex count
+// exceeds the bounding-box estimate, so the vertex count is used.
+func TestMaxPolygonToCellsSizeVertexFloor(t *testing.T) {
+	t.Parallel()
+
+	// A tiny polygon at a coarse resolution estimates very few cells, so the
+	// vertex count dominates.
+	loop := GeoLoop{
+		{Lat: 0.0, Lng: 0.0},
+		{Lat: 0.0, Lng: 0.001},
+		{Lat: 0.001, Lng: 0.001},
+		{Lat: 0.001, Lng: 0.0},
+	}
+
+	size, err := maxPolygonToCellsSize(GeoPolygon{GeoLoop: loop}, 0)
+	if err != nil {
+		t.Fatalf("maxPolygonToCellsSize: %v", err)
+	}
+
+	if size < len(loop)+polygonToCellsBuffer {
+		t.Fatalf("size %d should be at least vertex count plus buffer", size)
+	}
+}
+
+// TestBBoxFromGeoLoop covers the empty, normal, and transmeridian cases.
+func TestBBoxFromGeoLoop(t *testing.T) {
+	t.Parallel()
+
+	if got := bboxFromGeoLoop(GeoLoop{}); got != (bbox{}) {
+		t.Fatalf("bboxFromGeoLoop(empty): got %+v, want zero", got)
+	}
+
+	normal := bboxFromGeoLoop(sfSquareLoop)
+	if normal.north <= normal.south || normal.east <= normal.west {
+		t.Fatalf("bboxFromGeoLoop(normal): unexpected %+v", normal)
+	}
+
+	if normal.isTransmeridian() {
+		t.Fatal("sf square should not be transmeridian")
+	}
+
+	trans := bboxFromGeoLoop(GeoLoop{
+		{Lat: 10, Lng: 178},
+		{Lat: 10, Lng: -178},
+		{Lat: -10, Lng: -178},
+		{Lat: -10, Lng: 178},
+	})
+	if !trans.isTransmeridian() {
+		t.Fatalf("expected transmeridian bbox, got %+v", trans)
+	}
+}
+
+// TestBBoxContains covers the standard and transmeridian containment branches.
+func TestBBoxContains(t *testing.T) {
+	t.Parallel()
+
+	standard := bbox{north: 10, south: -10, east: 10, west: -10}
+	if !standard.contains(LatLng{Lat: 0, Lng: 0}) {
+		t.Fatal("standard bbox should contain origin")
+	}
+
+	if standard.contains(LatLng{Lat: 20, Lng: 0}) {
+		t.Fatal("standard bbox should not contain out-of-range latitude")
+	}
+
+	if standard.contains(LatLng{Lat: 0, Lng: 20}) {
+		t.Fatal("standard bbox should not contain out-of-range longitude")
+	}
+
+	trans := bbox{north: 10, south: -10, east: -170, west: 170}
+	if !trans.contains(LatLng{Lat: 0, Lng: 179}) || !trans.contains(LatLng{Lat: 0, Lng: -179}) {
+		t.Fatal("transmeridian bbox should contain points on both sides")
+	}
+
+	if trans.contains(LatLng{Lat: 0, Lng: 0}) {
+		t.Fatal("transmeridian bbox should not contain the prime meridian")
+	}
+}
+
+// TestHexRadiusKm covers both the hexagon and pentagon boundary branches.
+func TestHexRadiusKm(t *testing.T) {
+	t.Parallel()
+
+	hexagon := CellFromString("8928308280fffff")
+	if !hexagon.IsValid() || hexagon.IsPentagon() {
+		t.Fatalf("fixture %015x should be a valid hexagon", uint64(hexagon))
+	}
+
+	if got := hexagon.hexRadiusKm(); got <= 0 {
+		t.Fatalf("hexRadiusKm(hexagon): got %v, want positive", got)
+	}
+
+	pentagons, err := Pentagons(9)
+	if err != nil {
+		t.Fatalf("Pentagons: %v", err)
+	}
+
+	if got := pentagons[0].hexRadiusKm(); got <= 0 {
+		t.Fatalf("hexRadiusKm(pentagon): got %v, want positive", got)
+	}
+}
+
+// TestBBoxEstimatesResolutionError covers the resolution error path of both
+// estimators.
+func TestBBoxEstimatesResolutionError(t *testing.T) {
+	t.Parallel()
+
+	box := bboxFromGeoLoop(sfSquareLoop)
+	if _, err := bboxHexEstimate(box, -1); err == nil {
+		t.Fatal("bboxHexEstimate(res -1): got nil error, want failure")
+	}
+
+	if _, err := lineHexEstimate(sfSquareLoop[0], sfSquareLoop[1], -1); err == nil {
+		t.Fatal("lineHexEstimate(res -1): got nil error, want failure")
+	}
+}
+
+// TestBBoxHexEstimateNonFinite covers the non-finite estimate guard, reached when
+// a bounding-box corner is NaN so the diagonal and area become NaN.
+func TestBBoxHexEstimateNonFinite(t *testing.T) {
+	t.Parallel()
+
+	box := bbox{north: math.NaN(), south: 0, east: 1, west: 0}
+	if _, err := bboxHexEstimate(box, 5); !errors.Is(err, ErrFailed) {
+		t.Fatalf("bboxHexEstimate(NaN): got %v, want ErrFailed", err)
+	}
+}
+
+// TestLineHexEstimateNonFinite covers the non-finite distance guard, reached when
+// an endpoint is NaN so the great-circle distance is NaN.
+func TestLineHexEstimateNonFinite(t *testing.T) {
+	t.Parallel()
+
+	origin := LatLng{Lat: math.NaN(), Lng: 0}
+	destination := LatLng{Lat: 1, Lng: 1}
+
+	if _, err := lineHexEstimate(origin, destination, 5); !errors.Is(err, ErrFailed) {
+		t.Fatalf("lineHexEstimate(NaN): got %v, want ErrFailed", err)
+	}
+}
+
+// TestPointInsideGeoLoopNudges covers the latitude and longitude nudge branches
+// of the ray-casting test, reached when the point lies exactly on a vertex
+// latitude or a segment-endpoint longitude.
+func TestPointInsideGeoLoopNudges(t *testing.T) {
+	t.Parallel()
+
+	loop := GeoLoop{
+		{Lat: 0, Lng: 0},
+		{Lat: 0, Lng: 2},
+		{Lat: 2, Lng: 2},
+		{Lat: 2, Lng: 0},
+	}
+	box := bboxFromGeoLoop(loop)
+
+	// Latitude exactly equal to a vertex latitude triggers the lat nudge; the
+	// longitude equal to a vertex longitude triggers the lng nudge.
+	onVertex := LatLng{Lat: 0, Lng: 0}
+	_ = pointInsideGeoLoop(loop, box, onVertex)
+
+	// A point clearly inside still reports as contained after the nudges.
+	inside := LatLng{Lat: 1, Lng: 1}
+	if !pointInsideGeoLoop(loop, box, inside) {
+		t.Fatal("interior point should be contained")
+	}
+}
+
+// TestPointInsidePolygonHole checks the point-in-polygon test excludes points in
+// holes and the normalizeLng helper.
+func TestPointInsidePolygonHole(t *testing.T) {
+	t.Parallel()
+
+	hole := GeoLoop{
+		{Lat: 37.78, Lng: -122.45},
+		{Lat: 37.78, Lng: -122.42},
+		{Lat: 37.76, Lng: -122.42},
+		{Lat: 37.76, Lng: -122.45},
+	}
+	polygon := GeoPolygon{GeoLoop: sfSquareLoop, Holes: []GeoLoop{hole}}
+	bboxes := bboxesFromGeoPolygon(polygon)
+
+	inHole := LatLng{Lat: 37.77, Lng: -122.435}
+	if pointInsidePolygon(polygon, bboxes, inHole) {
+		t.Fatal("point in hole should not be contained")
+	}
+
+	inPolygon := LatLng{Lat: 37.80, Lng: -122.40}
+	if !pointInsidePolygon(polygon, bboxes, inPolygon) {
+		t.Fatal("point in polygon (outside hole) should be contained")
+	}
+
+	outside := LatLng{Lat: 0, Lng: 0}
+	if pointInsidePolygon(polygon, bboxes, outside) {
+		t.Fatal("point outside should not be contained")
+	}
+}
+
+// TestNormalizeLng covers the transmeridian normalization branch.
+func TestNormalizeLng(t *testing.T) {
+	t.Parallel()
+
+	if got := normalizeLng(-1, true); got != -1+twoPiRad {
+		t.Fatalf("normalizeLng(-1, true): got %v, want %v", got, -1+twoPiRad)
+	}
+
+	if got := normalizeLng(-1, false); got != -1 {
+		t.Fatalf("normalizeLng(-1, false): got %v, want -1", got)
+	}
+
+	if got := normalizeLng(1, true); got != 1 {
+		t.Fatalf("normalizeLng(1, true): got %v, want 1", got)
+	}
+}
+
+// TestPolygonToCellsTransmeridian ports the testPolygonToCells.c transmeridian
+// regressions: a small prime-meridian box, the antimeridian-crossing box, and a
+// >4-vertex complex transmeridian polygon, each with an exact expected count.
+// The complex case guards the historical bug of using min/max longitude as the
+// transmeridian bounds.
+func TestPolygonToCellsTransmeridian(t *testing.T) {
+	t.Parallel()
+
+	tests := map[string]struct {
+		giveLoop GeoLoop
+		giveRes  int
+		want     int
+	}{
+		"prime_meridian": {
+			giveLoop: radLoop([][2]float64{
+				{0.01, 0.01}, {0.01, -0.01}, {-0.01, -0.01}, {-0.01, 0.01},
+			}),
+			giveRes: 7,
+			want:    4228,
+		},
+		"transmeridian": {
+			giveLoop: radLoop([][2]float64{
+				{0.01, -math.Pi + 0.01}, {0.01, math.Pi - 0.01},
+				{-0.01, math.Pi - 0.01}, {-0.01, -math.Pi + 0.01},
+			}),
+			giveRes: 7,
+			want:    4238,
+		},
+		"complex": {
+			giveLoop: radLoop([][2]float64{
+				{0.1, -math.Pi + 0.00001}, {0.1, math.Pi - 0.00001},
+				{0.05, math.Pi - 0.2}, {-0.1, math.Pi - 0.00001},
+				{-0.1, -math.Pi + 0.00001}, {-0.05, -math.Pi + 0.2},
+			}),
+			giveRes: 4,
+			want:    1204,
+		},
+	}
+
+	for name, tt := range tests {
+		t.Run(name, func(t *testing.T) {
+			t.Parallel()
+
+			cells, err := PolygonToCells(GeoPolygon{GeoLoop: tt.giveLoop}, tt.giveRes)
+			if err != nil {
+				t.Fatalf("PolygonToCells: %v", err)
+			}
+
+			if len(cells) != tt.want {
+				t.Fatalf("got %d cells, want %d", len(cells), tt.want)
+			}
+		})
+	}
+}
diff --git a/x/h3go/sets.go b/x/h3go/sets.go
index 527b001..062fcc1 100644
--- a/x/h3go/sets.go
+++ b/x/h3go/sets.go
@@ -34,9 +34,10 @@ func setH3Index(res, baseCell, initDigit int) Cell {
 }
 
 // Res0Cells returns all the cells at resolution 0 (the base cells). The error
-// return is always nil; it exists for parity with the cgo-backed h3 package.
+// return is always nil; it exists to match the signature of the other cell
+// enumerators.
 func Res0Cells() ([]Cell, error) {
-	out := make([]Cell, numBaseCells)
+	out := make([]Cell, NumBaseCells)
 	for bc := range out {
 		out[bc] = setH3Index(0, bc, centerDigit)
 	}
@@ -46,12 +47,12 @@ func Res0Cells() ([]Cell, error) {
 
 // Pentagons returns all the pentagons at the given resolution.
 func Pentagons(res int) ([]Cell, error) {
-	if res < 0 || res > maxResolution {
+	if res < 0 || res > MaxResolution {
 		return nil, ErrResolutionDomain
 	}
-	out := make([]Cell, 0, numPentagons)
+	out := make([]Cell, 0, NumPentagons)
 
-	for bc := range numBaseCells {
+	for bc := range NumBaseCells {
 		if h3core.IsBaseCellPentagon[bc] {
 			out = append(out, setH3Index(res, bc, centerDigit))
 		}
diff --git a/x/h3go/sets_test.go b/x/h3go/sets_test.go
index 6dda094..e387dc3 100644
--- a/x/h3go/sets_test.go
+++ b/x/h3go/sets_test.go
@@ -41,7 +41,7 @@ func TestSetH3IndexValue(t *testing.T) {
 		}
 	}
 
-	for r := 6; r <= maxResolution; r++ {
+	for r := 6; r <= MaxResolution; r++ {
 		if h.indexDigit(r) != invalidDigit {
 			t.Fatalf("digit %d = %d, want %d", r, h.indexDigit(r), invalidDigit)
 		}
@@ -62,8 +62,8 @@ func TestRes0Cells(t *testing.T) {
 		t.Fatalf("Res0Cells: %v", err)
 	}
 
-	if len(cells) != numBaseCells {
-		t.Fatalf("Res0Cells: got %d cells, want %d", len(cells), numBaseCells)
+	if len(cells) != NumBaseCells {
+		t.Fatalf("Res0Cells: got %d cells, want %d", len(cells), NumBaseCells)
 	}
 
 	for _, c := range cells {
@@ -71,6 +71,15 @@ func TestRes0Cells(t *testing.T) {
 			t.Fatalf("Res0Cells: %015x is not a valid res-0 cell", uint64(c))
 		}
 	}
+
+	// testBaseCells.c getRes0Cells regression: the first and last base cells.
+	if cells[0] != 0x8001fffffffffff {
+		t.Fatalf("first base cell: got %015x, want 8001fffffffffff", uint64(cells[0]))
+	}
+
+	if cells[121] != 0x80f3fffffffffff {
+		t.Fatalf("last base cell: got %015x, want 80f3fffffffffff", uint64(cells[121]))
+	}
 }
 
 // TestPentagons checks that Pentagons returns the 12 pentagons at each
@@ -78,14 +87,14 @@ func TestRes0Cells(t *testing.T) {
 func TestPentagons(t *testing.T) {
 	t.Parallel()
 
-	for res := 0; res <= maxResolution; res++ {
+	for res := 0; res <= MaxResolution; res++ {
 		pents, err := Pentagons(res)
 		if err != nil {
 			t.Fatalf("Pentagons(%d): %v", res, err)
 		}
 
-		if len(pents) != numPentagons {
-			t.Fatalf("Pentagons(%d): got %d, want %d", res, len(pents), numPentagons)
+		if len(pents) != NumPentagons {
+			t.Fatalf("Pentagons(%d): got %d, want %d", res, len(pents), NumPentagons)
 		}
 
 		for _, p := range pents {
diff --git a/x/h3go/vertex.go b/x/h3go/vertex.go
index 357f2dd..d448dd5 100644
--- a/x/h3go/vertex.go
+++ b/x/h3go/vertex.go
@@ -35,7 +35,7 @@ const directionIndexOffset = 2
 // pentagonDirectionFaces maps each pentagon base cell to the icosahedron faces
 // found in each axial direction, in order starting at the J axis. Generated by
 // the upstream generatePentagonDirectionFaces script.
-var pentagonDirectionFaces = [numPentagons]struct {
+var pentagonDirectionFaces = [NumPentagons]struct {
 	baseCell int
 	faces    [numPentVerts]int
 }{
@@ -366,7 +366,7 @@ func (v Vertex) Resolution() int {
 }
 
 // IndexDigit returns the indexing digit of the vertex at res, for res in
-// [1, maxResolution].
+// [1, MaxResolution].
 func (v Vertex) IndexDigit(res int) (int, error) {
 	return Cell(v).IndexDigit(res)
 }
diff --git a/x/h3go/vertex_test.go b/x/h3go/vertex_test.go
index 9ba27e0..0c9d3df 100644
--- a/x/h3go/vertex_test.go
+++ b/x/h3go/vertex_test.go
@@ -149,7 +149,7 @@ func TestDirectionForVertexNumInvalid(t *testing.T) {
 func TestVertexRotationsError(t *testing.T) {
 	t.Parallel()
 
-	corrupt := CellFromString("8928308280fffff").setBaseCell(numBaseCells)
+	corrupt := CellFromString("8928308280fffff").setBaseCell(NumBaseCells)
 
 	if _, err := corrupt.vertexRotations(); !errors.Is(err, ErrCellInvalid) {
 		t.Fatalf("vertexRotations(corrupt): got %v, want ErrCellInvalid", err)
@@ -171,7 +171,7 @@ func TestVertexOwnerFails(t *testing.T) {
 
 	// A res-0 cell with an out-of-range base cell enters the owner search (res 0)
 	// but cannot determine its vertex direction, so the lookup fails.
-	corrupt := Cell(h3Init).setMode(cellMode).setBaseCell(numBaseCells)
+	corrupt := Cell(h3Init).setMode(cellMode).setBaseCell(NumBaseCells)
 
 	if _, err := corrupt.Vertex(0); !errors.Is(err, ErrFailed) {
 		t.Fatalf("Vertex(corrupt res0): got %v, want ErrFailed", err)
@@ -320,7 +320,7 @@ func TestBaseCellToCCWrot60Invalid(t *testing.T) {
 		t.Fatalf("baseCellToCCWrot60(0, -1): got %d, want %d", got, invalidRotations)
 	}
 
-	if got := baseCellToCCWrot60(numBaseCells, 0); got != invalidRotations {
+	if got := baseCellToCCWrot60(NumBaseCells, 0); got != invalidRotations {
 		t.Fatalf("baseCellToCCWrot60(out of range, 0): got %d, want %d", got, invalidRotations)
 	}
 }
@@ -377,7 +377,7 @@ func TestVertexNeighborErrors(t *testing.T) {
 func TestVertexLatLngOwnerError(t *testing.T) {
 	t.Parallel()
 
-	corrupt := Cell(h3Init).setMode(vertexMode).setBaseCell(numBaseCells)
+	corrupt := Cell(h3Init).setMode(vertexMode).setBaseCell(NumBaseCells)
 
 	if _, err := Vertex(corrupt).LatLng(); !errors.Is(err, ErrCellInvalid) {
 		t.Fatalf("LatLng(corrupt owner): got %v, want ErrCellInvalid", err)
@@ -396,3 +396,38 @@ func TestVertexLatLngOwnerError(t *testing.T) {
 		t.Fatal("IsValidVertex(valid) should be true")
 	}
 }
+
+// TestCellToVertexInvalidLiterals ports the testVertex.c cellToVertex_invalid2
+// and invalid3 regressions: specific malformed indexes must fail with
+// ErrCellInvalid.
+func TestCellToVertexInvalidLiterals(t *testing.T) {
+	t.Parallel()
+
+	tests := map[string]struct {
+		giveCell      Cell
+		giveVertexNum int
+	}{
+		"invalid2": {giveCell: Cell(0x685b2396e900fff9), giveVertexNum: 2},
+		"invalid3": {giveCell: Cell(0x20ff20202020ff35), giveVertexNum: 0},
+	}
+
+	for name, tt := range tests {
+		t.Run(name, func(t *testing.T) {
+			t.Parallel()
+
+			if _, err := tt.giveCell.Vertex(tt.giveVertexNum); !errors.Is(err, ErrCellInvalid) {
+				t.Fatalf("Vertex: got %v, want ErrCellInvalid", err)
+			}
+		})
+	}
+}
+
+// TestIsValidVertexKnownLiteral ports the testVertex.c isValidVertex_hex
+// regression: a specific known-valid vertex index validates.
+func TestIsValidVertexKnownLiteral(t *testing.T) {
+	t.Parallel()
+
+	if !IsValidVertex(Vertex(0x2222597fffffffff)) {
+		t.Fatal("IsValidVertex(0x2222597fffffffff): got false, want true")
+	}
+}