diff --git a/CMakeLists.txt b/CMakeLists.txt index 9426b868b..92f29201d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -184,6 +184,7 @@ set(LIB_SOURCE_FILES src/h3lib/lib/localij.c src/h3lib/lib/latLng.c src/h3lib/lib/directedEdge.c + src/h3lib/lib/iterGosper.c src/h3lib/lib/mathExtensions.c src/h3lib/lib/iterators.c src/h3lib/lib/faceijk.c @@ -281,6 +282,7 @@ set(OTHER_SOURCE_FILES src/apps/testapps/testH3IteratorsInternal.c src/apps/testapps/testMathExtensionsInternal.c src/apps/testapps/testDescribeH3Error.c + src/apps/testapps/testGosperIter.c src/apps/testapps/testGeoLoopArea.c src/apps/miscapps/cellToBoundaryHier.c src/apps/miscapps/cellToLatLngHier.c @@ -320,6 +322,7 @@ set(OTHER_SOURCE_FILES src/apps/benchmarks/benchmarkGridDiskCells.c src/apps/benchmarks/benchmarkGridPathCells.c src/apps/benchmarks/benchmarkDirectedEdge.c + src/apps/benchmarks/benchmarkGosperIter.c src/apps/benchmarks/benchmarkVertex.c src/apps/benchmarks/benchmarkIsValidCell.c src/apps/benchmarks/benchmarkH3Api.c @@ -672,6 +675,8 @@ if(BUILD_BENCHMARKS) src/apps/benchmarks/benchmarkGridPathCells.c) add_h3_benchmark(benchmarkDirectedEdge src/apps/benchmarks/benchmarkDirectedEdge.c) + add_h3_benchmark(benchmarkGosperIter + src/apps/benchmarks/benchmarkGosperIter.c) add_h3_benchmark(benchmarkVertex src/apps/benchmarks/benchmarkVertex.c) add_h3_benchmark(benchmarkIsValidCell src/apps/benchmarks/benchmarkIsValidCell.c) diff --git a/CMakeTests.cmake b/CMakeTests.cmake index 7500b4729..f48c8f931 100644 --- a/CMakeTests.cmake +++ b/CMakeTests.cmake @@ -260,6 +260,7 @@ add_h3_test(testH3IteratorsInternal src/apps/testapps/testH3IteratorsInternal.c) add_h3_test(testMathExtensionsInternal src/apps/testapps/testMathExtensionsInternal.c) add_h3_test(testDescribeH3Error src/apps/testapps/testDescribeH3Error.c) +add_h3_test(testGosperIter src/apps/testapps/testGosperIter.c) add_h3_test(testGeoLoopArea src/apps/testapps/testGeoLoopArea.c) add_h3_test_with_arg(testH3NeighborRotations diff --git a/src/apps/benchmarks/benchmarkGosperIter.c b/src/apps/benchmarks/benchmarkGosperIter.c new file mode 100644 index 000000000..98232b708 --- /dev/null +++ b/src/apps/benchmarks/benchmarkGosperIter.c @@ -0,0 +1,53 @@ +/* + * 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. + */ +#include "benchmark.h" +#include "iterators.h" + +H3Index hex2 = 0x820887fffffffff; // res 2 hexagon +H3Index pent2 = 0x820807fffffffff; // res 2 pentagon + +void exhaustGosper(H3Index cell, int childRes) { + IterEdgesGosper iter = iterInitGosper(cell, childRes); + while (iter.e) iterStepGosper(&iter); +} + +BEGIN_BENCHMARKS(); + +// res: 2 -> 2 +BENCHMARK(hex_plus0, 50000000, { exhaustGosper(hex2, 2); }); +BENCHMARK(pent_plus0, 50000000, { exhaustGosper(pent2, 2); }); + +// res: 2 -> 3 +BENCHMARK(hex_plus1, 500000, { exhaustGosper(hex2, 3); }); +BENCHMARK(pent_plus1, 500000, { exhaustGosper(pent2, 3); }); + +// res: 2 -> 4 +BENCHMARK(hex_plus2, 100000, { exhaustGosper(hex2, 4); }); +BENCHMARK(pent_plus2, 100000, { exhaustGosper(pent2, 4); }); + +// res: 2 -> 7 +BENCHMARK(hex_plus5, 1000, { exhaustGosper(hex2, 7); }); +BENCHMARK(pent_plus5, 1000, { exhaustGosper(pent2, 7); }); + +// res: 2 -> 10 +BENCHMARK(hex_plus8, 100, { exhaustGosper(hex2, 10); }); +BENCHMARK(pent_plus8, 100, { exhaustGosper(pent2, 10); }); + +// res: 2 -> 13 +BENCHMARK(hex_plus11, 10, { exhaustGosper(hex2, 13); }); +BENCHMARK(pent_plus11, 10, { exhaustGosper(pent2, 13); }); + +END_BENCHMARKS(); diff --git a/src/apps/testapps/testGosperIter.c b/src/apps/testapps/testGosperIter.c new file mode 100644 index 000000000..59a09d896 --- /dev/null +++ b/src/apps/testapps/testGosperIter.c @@ -0,0 +1,320 @@ +/* + * 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. + */ + +#include + +#include "iterators.h" +#include "mathExtensions.h" +#include "test.h" +#include "utility.h" + +// Check iterator mechanics: +// 1. produces correct number of edges +// 2. each step produces an edge that's different from the previous edge +// 3. stays H3_NULL after finishing +void check_iterator_mechanics(H3Index h, int childRes) { + bool isPent = H3_EXPORT(isPentagon)(h); + int parentRes = H3_EXPORT(getResolution)(h); + int numSides = isPent ? 5 : 6; + int64_t expectedEdges = numSides * _ipow(3, childRes - parentRes); + + IterEdgesGosper iter = iterInitGosper(h, childRes); + t_assert(iter.remaining == expectedEdges, "correct number of edges"); + + int64_t i; + for (i = 0; iter.e; i++) { + H3Index prev = iter.e; + iterStepGosper(&iter); + t_assert(iter.e != prev, "edge should advance"); + } + t_assert(i == expectedEdges, "correct number of edges"); + + // iterator should stay exhausted + for (int j = 0; j < 100; j++) { + t_assert(iter.remaining == 0, "num edges to zero"); + t_assert(iter.e == H3_NULL, "iterator exhausted"); + iterStepGosper(&iter); + } +} + +// Per-edge correctness: +// 1. each edge is a valid directed edge +// 2. each edge is at childRes +// 3. confirm the edge is on the Gosper island boundary: +// origin cell is a child of parent cell, destination cell is not +void check_edge_correctness(H3Index h, int childRes) { + int parentRes = H3_EXPORT(getResolution)(h); + IterEdgesGosper iter = iterInitGosper(h, childRes); + + while (iter.e) { + t_assert(H3_EXPORT(isValidDirectedEdge)(iter.e), "edge is valid"); + t_assert(H3_EXPORT(getResolution)(iter.e) == childRes, + "correct resolution"); + + H3Index s; // source/origin cell + H3Index d; // destination cell + H3Index ps; // parent of source + H3Index pd; // parent of destination + + t_assertSuccess(H3_EXPORT(getDirectedEdgeOrigin)(iter.e, &s)); + t_assertSuccess(H3_EXPORT(getDirectedEdgeDestination)(iter.e, &d)); + t_assertSuccess(H3_EXPORT(cellToParent)(s, parentRes, &ps)); + t_assertSuccess(H3_EXPORT(cellToParent)(d, parentRes, &pd)); + + t_assert(ps == h, "origin *is* child of parent cell"); + t_assert(pd != h, "destination *is not* child of parent cell"); + + iterStepGosper(&iter); + } +} + +// Check that edge_a's last vertex matches edge_b's first vertex, +// i.e. the edges connect end-to-start. Tolerance is relative to +// the shorter edge length to handle varying resolutions. +// TODO: Replace floating-point distance check with exact vertex comparison +// once an edgeToVertexes function exists. +bool do_edges_connect(H3Index edge_a, H3Index edge_b) { + double len_a, len_b; + H3_EXPORT(edgeLengthRads)(edge_a, &len_a); + H3_EXPORT(edgeLengthRads)(edge_b, &len_b); + double tol = MIN(len_a, len_b) / 1000.0; + + CellBoundary bd_a, bd_b; + t_assertSuccess(H3_EXPORT(directedEdgeToBoundary)(edge_a, &bd_a)); + t_assertSuccess(H3_EXPORT(directedEdgeToBoundary)(edge_b, &bd_b)); + + LatLng end_a = bd_a.verts[bd_a.numVerts - 1]; + LatLng start_b = bd_b.verts[0]; + double dist = H3_EXPORT(greatCircleDistanceRads)(&end_a, &start_b); + + return dist < tol; +} + +// Loop structure: +// 1. consecutive edges share an endpoint (floating-point test) +// 2. loop closes (returns to where it started) +void check_loop_valid(H3Index h, int childRes) { + IterEdgesGosper iter = iterInitGosper(h, childRes); + H3Index first_edge = iter.e; + H3Index prev_edge = iter.e; + iterStepGosper(&iter); + + while (iter.e) { + t_assert(do_edges_connect(prev_edge, iter.e), "edges should connect"); + + prev_edge = iter.e; + iterStepGosper(&iter); + } + + t_assert(do_edges_connect(prev_edge, first_edge), "loop should close"); +} + +void check_gosper_island(H3Index h, int childRes) { + check_iterator_mechanics(h, childRes); + check_edge_correctness(h, childRes); + check_loop_valid(h, childRes); +} + +// Verify iterator produces some cyclic rotation of the expected edges +void check_expected_edges(H3Index h, int childRes, H3Index *expected, + int numExpected) { + IterEdgesGosper iter = iterInitGosper(h, childRes); + t_assert(iter.remaining == numExpected, "correct number of edges"); + + // Find rotation offset: where does the iterator's first edge appear? + int offset = -1; + for (int i = 0; i < numExpected; i++) { + if (expected[i] == iter.e) { + offset = i; + break; + } + } + t_assert(offset >= 0, "first edge found in expected array"); + + for (int i = 0; iter.e; i++) { + t_assert(iter.e == expected[(offset + i) % numExpected], + "correct edge"); + iterStepGosper(&iter); + } +} + +SUITE(gosper_iter) { + // Exhaustive test of all cells, resolutions 0,1,2,3,4. + // Expanding to Gosper islands up to res 4. + TEST(coarse_resolutions_exhaustive) { + for (int childRes = 0; childRes <= 4; childRes++) { + for (int parentRes = 0; parentRes <= childRes; parentRes++) { + IterCellsResolution cell_iter = iterInitRes(parentRes); + while (cell_iter.h) { + check_gosper_island(cell_iter.h, childRes); + iterStepRes(&cell_iter); + } + } + } + } + + // Since an exhaustive test of all resolutions would take too long, + // we test a few specific cells at finer resolutions, including + // expanding up to res 15. + TEST(finer_resolutions) { + H3Index cells[] = { + 0x8508000ffffffff, // res 5 hexagon + 0x85080003fffffff, // res 5 pentagon + + 0x8808000009fffff, // res 8 hexagon + 0x8808000001fffff, // res 8 pentagon + + 0x8a0800000017fff, // res 10 hexagon + 0x8a0800000007fff, // res 10 pentagon + + 0x8e754e64992d6c7, // res 14 hexagon + 0x8e0800000000007, // res 14 pentagon + + 0x8f754e64992d6d8, // res 15 hexagon + 0x8f0800000000000, // res 15 pentagon + }; + for (int i = 0; i < (int)ARRAY_SIZE(cells); i++) { + int parentRes = H3_EXPORT(getResolution)(cells[i]); + for (int gap = 0; gap <= 5; gap++) { + int childRes = parentRes + gap; + if (childRes > 15) continue; + check_gosper_island(cells[i], childRes); + } + } + } + + // `check_expected_edges` tests below check for a specific sequence + // edges for a given cell and resolution input. Any cyclic rotation + // should pass the test. + + TEST(just_hex_2) { + H3Index h = 0x820887fffffffff; + int res = 2; + H3Index expected[] = { + 0x1320887fffffffff, 0x1120887fffffffff, 0x1520887fffffffff, + 0x1420887fffffffff, 0x1620887fffffffff, 0x1220887fffffffff, + }; + check_expected_edges(h, res, expected, ARRAY_SIZE(expected)); + } + + TEST(just_pent_2) { + H3Index h = 0x820807fffffffff; + int res = 2; + H3Index expected[] = { + 0x1320807fffffffff, 0x1520807fffffffff, 0x1420807fffffffff, + 0x1620807fffffffff, 0x1220807fffffffff, + }; + check_expected_edges(h, res, expected, ARRAY_SIZE(expected)); + } + + TEST(pent_edges_2_to_4) { + H3Index h = 0x820807fffffffff; + int res = 4; + H3Index expected[] = { + 0x11408053ffffffff, 0x15408053ffffffff, 0x1140805bffffffff, + 0x1540805bffffffff, 0x1440805bffffffff, 0x15408059ffffffff, + 0x14408059ffffffff, 0x16408059ffffffff, 0x1440805dffffffff, + 0x1540804bffffffff, 0x1440804bffffffff, 0x15408049ffffffff, + 0x14408049ffffffff, 0x16408049ffffffff, 0x1440804dffffffff, + 0x1640804dffffffff, 0x1240804dffffffff, 0x16408045ffffffff, + 0x14408069ffffffff, 0x16408069ffffffff, 0x1440806dffffffff, + 0x1640806dffffffff, 0x1240806dffffffff, 0x16408065ffffffff, + 0x12408065ffffffff, 0x13408065ffffffff, 0x12408067ffffffff, + 0x1640802dffffffff, 0x1240802dffffffff, 0x16408025ffffffff, + 0x12408025ffffffff, 0x13408025ffffffff, 0x12408027ffffffff, + 0x13408027ffffffff, 0x11408027ffffffff, 0x13408023ffffffff, + 0x12408035ffffffff, 0x13408035ffffffff, 0x12408037ffffffff, + 0x13408037ffffffff, 0x11408037ffffffff, 0x13408033ffffffff, + 0x11408033ffffffff, 0x15408033ffffffff, 0x1140803bffffffff, + }; + check_expected_edges(h, res, expected, ARRAY_SIZE(expected)); + } + + TEST(pent_edges_1_to_3) { + H3Index h = 0x81083ffffffffff; + int res = 3; + H3Index expected[] = { + 0x113082bfffffffff, 0x1330829fffffffff, 0x1130829fffffffff, + 0x1530829fffffffff, 0x113082dfffffffff, 0x153082dfffffffff, + 0x143082dfffffffff, 0x153082cfffffffff, 0x143082cfffffffff, + 0x1530821fffffffff, 0x1130825fffffffff, 0x1530825fffffffff, + 0x1430825fffffffff, 0x1530824fffffffff, 0x1430824fffffffff, + 0x1630824fffffffff, 0x1430826fffffffff, 0x1630826fffffffff, + 0x1430835fffffffff, 0x1530834fffffffff, 0x1430834fffffffff, + 0x1630834fffffffff, 0x1430836fffffffff, 0x1630836fffffffff, + 0x1230836fffffffff, 0x1630832fffffffff, 0x1230832fffffffff, + 0x1630814fffffffff, 0x1430816fffffffff, 0x1630816fffffffff, + 0x1230816fffffffff, 0x1630812fffffffff, 0x1230812fffffffff, + 0x1330812fffffffff, 0x1230813fffffffff, 0x1330813fffffffff, + 0x123081efffffffff, 0x163081afffffffff, 0x123081afffffffff, + 0x133081afffffffff, 0x123081bfffffffff, 0x133081bfffffffff, + 0x113081bfffffffff, 0x1330819fffffffff, 0x1130819fffffffff, + }; + check_expected_edges(h, res, expected, ARRAY_SIZE(expected)); + } + + TEST(hex_edges_1_to_3) { + H3Index h = 0x81c67ffffffffff; + int res = 3; + H3Index expected[] = { + 0x133c64afffffffff, 0x123c64bfffffffff, 0x133c64bfffffffff, + 0x113c64bfffffffff, 0x133c649fffffffff, 0x113c649fffffffff, + 0x153c649fffffffff, 0x113c64dfffffffff, 0x153c64dfffffffff, + 0x113c66bfffffffff, 0x133c669fffffffff, 0x113c669fffffffff, + 0x153c669fffffffff, 0x113c66dfffffffff, 0x153c66dfffffffff, + 0x143c66dfffffffff, 0x153c66cfffffffff, 0x143c66cfffffffff, + 0x153c661fffffffff, 0x113c665fffffffff, 0x153c665fffffffff, + 0x143c665fffffffff, 0x153c664fffffffff, 0x143c664fffffffff, + 0x163c664fffffffff, 0x143c666fffffffff, 0x163c666fffffffff, + 0x143c675fffffffff, 0x153c674fffffffff, 0x143c674fffffffff, + 0x163c674fffffffff, 0x143c676fffffffff, 0x163c676fffffffff, + 0x123c676fffffffff, 0x163c672fffffffff, 0x123c672fffffffff, + 0x163c654fffffffff, 0x143c656fffffffff, 0x163c656fffffffff, + 0x123c656fffffffff, 0x163c652fffffffff, 0x123c652fffffffff, + 0x133c652fffffffff, 0x123c653fffffffff, 0x133c653fffffffff, + 0x123c65efffffffff, 0x163c65afffffffff, 0x123c65afffffffff, + 0x133c65afffffffff, 0x123c65bfffffffff, 0x133c65bfffffffff, + 0x113c65bfffffffff, 0x133c659fffffffff, 0x113c659fffffffff, + }; + check_expected_edges(h, res, expected, ARRAY_SIZE(expected)); + } + + TEST(hex_edges_2_to_4) { + H3Index h = 0x82c64ffffffffff; + int res = 4; + H3Index expected[] = { + 0x134c6497ffffffff, 0x114c6497ffffffff, 0x134c6493ffffffff, + 0x114c6493ffffffff, 0x154c6493ffffffff, 0x114c649bffffffff, + 0x154c649bffffffff, 0x144c649bffffffff, 0x154c6499ffffffff, + 0x114c64d3ffffffff, 0x154c64d3ffffffff, 0x114c64dbffffffff, + 0x154c64dbffffffff, 0x144c64dbffffffff, 0x154c64d9ffffffff, + 0x144c64d9ffffffff, 0x164c64d9ffffffff, 0x144c64ddffffffff, + 0x154c64cbffffffff, 0x144c64cbffffffff, 0x154c64c9ffffffff, + 0x144c64c9ffffffff, 0x164c64c9ffffffff, 0x144c64cdffffffff, + 0x164c64cdffffffff, 0x124c64cdffffffff, 0x164c64c5ffffffff, + 0x144c64e9ffffffff, 0x164c64e9ffffffff, 0x144c64edffffffff, + 0x164c64edffffffff, 0x124c64edffffffff, 0x164c64e5ffffffff, + 0x124c64e5ffffffff, 0x134c64e5ffffffff, 0x124c64e7ffffffff, + 0x164c64adffffffff, 0x124c64adffffffff, 0x164c64a5ffffffff, + 0x124c64a5ffffffff, 0x134c64a5ffffffff, 0x124c64a7ffffffff, + 0x134c64a7ffffffff, 0x114c64a7ffffffff, 0x134c64a3ffffffff, + 0x124c64b5ffffffff, 0x134c64b5ffffffff, 0x124c64b7ffffffff, + 0x134c64b7ffffffff, 0x114c64b7ffffffff, 0x134c64b3ffffffff, + 0x114c64b3ffffffff, 0x154c64b3ffffffff, 0x114c64bbffffffff, + }; + check_expected_edges(h, res, expected, ARRAY_SIZE(expected)); + } +} diff --git a/src/h3lib/include/iterators.h b/src/h3lib/include/iterators.h index 9e1de9bfd..899e08c9e 100644 --- a/src/h3lib/include/iterators.h +++ b/src/h3lib/include/iterators.h @@ -22,8 +22,10 @@ #ifndef ITERATORS_H #define ITERATORS_H +#include #include +#include "h3Index.h" #include "h3api.h" /** @@ -84,4 +86,45 @@ typedef struct { DECLSPEC IterCellsResolution iterInitRes(int res); DECLSPEC void iterStepRes(IterCellsResolution *iter); +/** + * IterEdgesGosper: iterator for the directed edges on the boundary of a cell's + * child set (the Gosper island outline) at a given resolution. + * + * Each yielded edge is a directed edge whose origin cell is a child of the + * parent cell and whose destination cell is not. The edges form a closed loop + * in counter-clockwise order around the Gosper island. The starting edge is + * not guaranteed; any cyclic rotation of the sequence is a valid output. + * + * Constructor: + * + * Initialize with `iterInitGosper(cell, childRes)`. + * Requires `childRes >= getResolution(cell)` and a valid cell. + * + * Iteration: + * + * Step iterator with `iterStepGosper`. + * The current edge is accessed via the `.e` member. + * The `.remaining` member gives the number of edges left to yield + * (including the current edge). + * When the iterator is exhausted, `.e` will be `H3_NULL` even after + * further calls to `iterStepGosper`. + */ +typedef struct { + H3Index e; // Current directed edge (H3_NULL when exhausted) + int64_t remaining; // Number of edges left to yield (including current) + + // Internal state — cyclic position trackers: + int8_t _edgePos; // Position in edge_dir cycle (period 6) + int8_t _walkPos[MAX_H3_RES + 1]; // Position in walk_digit cycle + // (period 18) per resolution level + + // Additional internal state + int8_t _parentRes; + int8_t _childRes; + bool _isPentagon; +} IterEdgesGosper; + +DECLSPEC IterEdgesGosper iterInitGosper(H3Index h, int childRes); +DECLSPEC void iterStepGosper(IterEdgesGosper *iter); + #endif diff --git a/src/h3lib/include/mathExtensions.h b/src/h3lib/include/mathExtensions.h index 1ac6cfdba..c5715b2c3 100644 --- a/src/h3lib/include/mathExtensions.h +++ b/src/h3lib/include/mathExtensions.h @@ -28,6 +28,11 @@ */ #define MAX(a, b) (((a) > (b)) ? (a) : (b)) +/** + * MIN returns the minimum of two values. + */ +#define MIN(a, b) (((a) < (b)) ? (a) : (b)) + /** Evaluates to true if a + b would overflow for int32 */ static inline bool ADD_INT32S_OVERFLOWS(int32_t a, int32_t b) { if (a > 0) { diff --git a/src/h3lib/lib/iterGosper.c b/src/h3lib/lib/iterGosper.c new file mode 100644 index 000000000..b78b91e2c --- /dev/null +++ b/src/h3lib/lib/iterGosper.c @@ -0,0 +1,231 @@ +/* + * 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. + */ +/** @file iterGosper.c + * @brief Iterator for the directed edges forming the "Gosper island" outline + * of a cell's child set at a given resolution. + * + * For a cell with resolution `r`, the boundary of the cell's children at + * resolution `R` forms a closed loop of directed edges (the Gosper island). + * For a hexagon this is 6 * 3^(R - r) edges; for a pentagon, 5 * 3^(R - r). + * We refer to `R - r` as the resolution delta. + * + * ## Geometric intuition + * + * At delta 0, we just return the cell's own edges in counter-clockwise + * order, cycling through the `edge_dir` array. + * + * At delta 1 (for a hexagon), each side of the parent subdivides into 3 edges + * at the finer resolution, giving 6 * 3 = 18 edges total. The origin cells + * share the parent's digit path, extended by one digit into the child + * resolution. The `walk_digit` array maps these WALK_CYCLE_LEN walk positions + * to H3 digit values: the 6 non-center child digits in counter-clockwise order, + * each repeated 3x for the subdivision. + * + * ## Pentagon handling + * + * Pentagons have 5 sides instead of 6: the K-axis (H3 digit 1) has no + * child cell (the deleted subsequence). To handle pentagonal input, we build + * the digit sequence as if it were a hexagon, then skip over the deleted + * subsequences (edges that land on the missing K-axis side). + */ + +#include +#include + +#include "coordijk.h" +#include "h3Assert.h" +#include "h3Index.h" +#include "iterators.h" +#include "mathExtensions.h" + +// H3 digit at each of the WALK_CYCLE_LEN boundary walk positions. +// The 6 non-center digits {1,5,4,6,2,3} in counter-clockwise order around +// the hexagon, each repeated 3x for the fractal subdivision at each side. +// Note: walk_digit[i] == edge_dir[(i/3 + 1) % 6]; the same cyclic sequence +// of digits, rotated by 1 and each element repeated 3x. +static const Direction walk_digit[] = { + K_AXES_DIGIT, K_AXES_DIGIT, K_AXES_DIGIT, // 1,1,1 + IK_AXES_DIGIT, IK_AXES_DIGIT, IK_AXES_DIGIT, // 5,5,5 + I_AXES_DIGIT, I_AXES_DIGIT, I_AXES_DIGIT, // 4,4,4 + IJ_AXES_DIGIT, IJ_AXES_DIGIT, IJ_AXES_DIGIT, // 6,6,6 + J_AXES_DIGIT, J_AXES_DIGIT, J_AXES_DIGIT, // 2,2,2 + JK_AXES_DIGIT, JK_AXES_DIGIT, JK_AXES_DIGIT}; // 3,3,3 + +#define WALK_CYCLE_LEN ((int)(sizeof(walk_digit) / sizeof(walk_digit[0]))) + +// H3 edge direction at each edge index, in counter-clockwise order around the +// hexagon. Stored in the directed edge's reserved bits. +// Numeric values: {3, 1, 5, 4, 6, 2} +static const Direction edge_dir[] = {JK_AXES_DIGIT, K_AXES_DIGIT, + IK_AXES_DIGIT, I_AXES_DIGIT, + IJ_AXES_DIGIT, J_AXES_DIGIT}; + +/** + * Advance the walk along origin cells on the Gosper island boundary, + * updating the child digits. + * + * Each resolution level cycles through WALK_CYCLE_LEN walk positions (6 sides * + * 3 segments per side). Across resolutions, each step at a coarser level maps + * to 3 steps at the next finer level. When the coarser level moves to the next + * cell, the finer level shifts back by 6 positions. + * + * The point where each level moves to the next cell depends on + * resolution parity: Class II (even r) at walkPos % 3 == 0, + * Class III (odd r) at walkPos % 3 == 1. + * + * @return true if the origin cell changed. + */ +static bool advanceOriginCell(int8_t *walkPos, H3Index *h, int8_t r, + int8_t parentRes) { + int prevDigit = H3_GET_INDEX_DIGIT(*h, r); + + // When moving to the next cell, recurse to advance the parent resolution. + if ((r > parentRes + 1) && (walkPos[r] % 3 == r % 2)) { + bool parentChanged = advanceOriginCell(walkPos, h, r - 1, parentRes); + if (parentChanged) { + // New parent cell: shift back by 6 positions. + walkPos[r] -= 6; + } + } + + // Advance to next position: Increase the cyclic index by one, but + // we use + (WALK_CYCLE_LEN + 1) == +1 (mod WALK_CYCLE_LEN), so that it + // stays positive, even after potentially subtracting by 6 above. + walkPos[r] = (walkPos[r] + WALK_CYCLE_LEN + 1) % WALK_CYCLE_LEN; + + // Update the child digit + Direction newDigit = walk_digit[walkPos[r]]; + H3_SET_INDEX_DIGIT(*h, r, newDigit); + + // A change in the finest digit is sufficient to detect a cell change. + return newDigit != prevDigit; +} + +/** + * Advance to the next boundary edge (origin cell + edge direction). + */ +static void advanceEdge(IterEdgesGosper *iter) { + bool cellChanged = advanceOriginCell(iter->_walkPos, &(iter->e), + iter->_childRes, iter->_parentRes); + + // If we move to a new origin cell, shift the _edgePos back by 2 + if (cellChanged) { + iter->_edgePos -= 2; + } + + // Advance to next edge. Increase the cyclic index by one, but + // we use +7 == +1 (mod 6) so that it stays positive, + // even after potentially subtracting by 2 above. + iter->_edgePos = (iter->_edgePos + 7) % 6; + H3_SET_RESERVED_BITS(iter->e, edge_dir[iter->_edgePos]); +} + +// Skip edges corresponding to deleted subsequences of a pentagon. +static void skipPentagonEdges(IterEdgesGosper *iter) { + while (iter->_isPentagon && ALWAYS(iter->e != H3_NULL) && + H3_GET_INDEX_DIGIT(iter->e, iter->_parentRes + 1) == 1) { + advanceEdge(iter); + } +} + +void iterStepGosper(IterEdgesGosper *iter) { + if (iter->e == H3_NULL) return; + + iter->remaining--; + if (iter->remaining <= 0) { + iter->e = H3_NULL; + return; + } + + if (iter->_parentRes == iter->_childRes) { + // Easy case of 0 resolution delta: + // Cycle through edge indices, skipping the + // deleted subsequence for pentagons. + iter->_edgePos++; + if (iter->_isPentagon && edge_dir[iter->_edgePos] == K_AXES_DIGIT) { + iter->_edgePos++; + } + + H3_SET_RESERVED_BITS(iter->e, edge_dir[iter->_edgePos]); + } else { + // Resolution delta > 0 + advanceEdge(iter); + + // Skip deleted subsequences for pentagons. + skipPentagonEdges(iter); + } +} + +// Note: No input validation here. This is an internal function where +// validation is expected to be done upstream. If this is ever added to the +// public API, we'll need to validate h, childRes >= parentRes, +// and childRes <= MAX_H3_RES. +IterEdgesGosper iterInitGosper(H3Index h, int childRes) { + int parentRes = H3_GET_RESOLUTION(h); + bool isPent = H3_EXPORT(isPentagon)(h); + + IterEdgesGosper iter = { + .e = h, + .remaining = 0, + + ._edgePos = 0, + ._walkPos = {0}, + + ._parentRes = parentRes, + ._childRes = childRes, + ._isPentagon = isPent, + }; + + // Set number of edges for the iterator to yield + int numSides = isPent ? 5 : 6; + iter.remaining = numSides * _ipow(3, childRes - parentRes); + + // Same-resolution fast path: no digit/resolution setup needed + if (parentRes == childRes) { + H3_SET_MODE(iter.e, H3_DIRECTEDEDGE_MODE); + // Skip the deleted subsequence for pentagons (if edge_dir happens + // to start there) + if (NEVER(isPent && edge_dir[iter._edgePos] == K_AXES_DIGIT)) { + iter._edgePos++; + } + H3_SET_RESERVED_BITS(iter.e, edge_dir[iter._edgePos]); + return iter; + } + + // Set resolution to child level and initialize boundary walk. + // The first child level starts at walk position 0. + // Subsequent levels start at either 14 or 16, depending on + // if those resolutions are Class II or III. + H3_SET_RESOLUTION(iter.e, childRes); + for (int r = parentRes + 1; r <= childRes; r++) { + if (r == parentRes + 1) { + iter._walkPos[r] = 0; + } else { + iter._walkPos[r] = (r % 2) ? 14 : 16; + } + H3_SET_INDEX_DIGIT(iter.e, r, walk_digit[iter._walkPos[r]]); + } + + // Set edge mode and initial direction + H3_SET_MODE(iter.e, H3_DIRECTEDEDGE_MODE); + H3_SET_RESERVED_BITS(iter.e, edge_dir[iter._edgePos]); + + // For pentagons, the walk may start on a deleted subsequence. + // Advance past those positions to reach the first valid edge. + skipPentagonEdges(&iter); + + return iter; +}