feat: Cutover to new non-selfintersecting curve

benchmarks/curve.bench.ts

@@ -3,7 +3,6 @@
 // Copyright (c) A5 contributors
 
 // Benchmarks for the space-filling curve: s -> cell decode, cell -> s encode,
-// and fractional-point location (IJToS).
 //
 // CI runs these same files against both the PR and its merge-base, so they
 // must run on either side of the L-system migration: the adapters below pick
@@ -105,19 +104,3 @@ describe('tripleToS', () => {
     );
   }
 });
-
-describe('IJToS', () => {
-  const values = sampleS(15, N);
-  const ijs: IJ[] = new Array(N);
-  for (let i = 0; i < N; i++) {
-    ijs[i] = centroidIJ(tripleOf(values[i], 15, 'uv'));
-  }
-  let i = 0;
-  bench(
-    'IJToS res 15',
-    () => {
-      lattice.IJToS(ijs[i++ & (N - 1)], 15, 'uv');
-    },
-    BENCH_OPTS
-  );
-});

modules/core/origin.ts

@@ -114,6 +114,38 @@ export function segmentToQuintant(segment: number, origin: Origin): {quintant: n
   return {quintant, orientation};
 }
 
+// Lookup tables for the two mappings above, built once at startup — there are
+// only 60 (origin, quintant) pairs. Indexed by origin.id * 5 + quintant
+// (resp. + segment, the global quintant number as encoded in serialized cell
+// ids). Call sites read these directly: a pair of flat array loads that V8
+// cannot pessimize, unlike the object-returning functions (a cold second call
+// site to quintantToSegment measurably degraded the optimized hot path).
+export const QUINTANT_TO_SEGMENT = new Uint8Array(60);
+export const QUINTANT_TO_ORIENTATION: Orientation[] = new Array(60);
+export const SEGMENT_TO_QUINTANT = new Uint8Array(60);
+export const SEGMENT_TO_ORIENTATION: Orientation[] = new Array(60);
+for (const origin of origins) {
+  for (let i = 0; i < 5; i++) {
+    const {segment, orientation} = quintantToSegment(i, origin);
+    QUINTANT_TO_SEGMENT[origin.id * 5 + i] = segment;
+    QUINTANT_TO_ORIENTATION[origin.id * 5 + i] = orientation;
+    const s2q = segmentToQuintant(i, origin);
+    SEGMENT_TO_QUINTANT[origin.id * 5 + i] = s2q.quintant;
+    SEGMENT_TO_ORIENTATION[origin.id * 5 + i] = s2q.orientation;
+  }
+}
+
+/**
+ * The `count` origins nearest to a point, by haversine distance, nearest first.
+ * Used by the boundary resolution in sphericalToCell: a point on (or within
+ * float noise of) a face seam or dodecahedron vertex may belong to a cell of
+ * the 2nd- or 3rd-nearest face.
+ */
+export function findNearestOrigins(point: Spherical, count: number): Origin[] {
+  const sorted = [...origins].sort((a, b) => haversine(point, a.axis) - haversine(point, b.axis));
+  return sorted.slice(0, count);
+}
+
 /**
  * Find the nearest origin to a point on the sphere
  * Uses haversine formula to calculate great-circle distance
@@ -136,25 +168,6 @@ export function isNearestOrigin(point: Spherical, origin: Origin): boolean {
   return haversine(point, origin.axis) > 0.49999999;
 }
 
-/**
- * Same as `findNearestOrigin` but takes a Cartesian unit vector. The
- * argmin of `1 − a·b` matches the argmin of haversine, so this returns
- * the same origin without any spherical-trig conversions.
- */
-export function findNearestOriginCartesian(c: Cartesian): Origin {
-  let minDistance = Infinity;
-  let nearest = origins[0];
-  for (const origin of origins) {
-    const ax = origin.axisCartesian;
-    const distance = 1 - (c[0] * ax[0] + c[1] * ax[1] + c[2] * ax[2]);
-    if (distance < minDistance) {
-      minDistance = distance;
-      nearest = origin;
-    }
-  }
-  return nearest;
-}
-
 /**
  * Modified haversine formula to calculate great-circle distance.
  * Retruns the "angle" between the two points. We need to minimize this to find the nearest origin

modules/core/tiling.ts

@@ -36,6 +36,51 @@ const FLAVOR_CENTERS = [0, 1, 2, 3].map(flavor => {
   return p.getCenter();
 });
 
+// The base PENTAGON under each flavor's orientation ops, flattened to
+// [x0,y0,...,x4,y4] for the allocation-free containment test below.
+const FLAVOR_PENTAGONS = [0, 1, 2, 3].map(flavor => {
+  const p = PENTAGON.clone();
+  if (flavor & 1) p.rotate180();
+  if (flavor & 2) p.reflectY();
+  const verts = p.getVertices();
+  const flat = new Float64Array(10);
+  for (let i = 0; i < 5; i++) {
+    flat[i * 2] = verts[i][0];
+    flat[i * 2 + 1] = verts[i][1];
+  }
+  return flat;
+});
+
+/**
+ * Signed containment margin of a point in the pentagon of (triple, flavor)
+ * (> 0 ⇔ strictly inside; the most-violated-edge cross product otherwise),
+ * tested in the SCALED quintant-0 frame (face coords rotated into quintant 0 and scaled
+ * by 2^resolution — the frame `_faceToEstimate` works in). In this frame the
+ * cell's pentagon is the flavor-oriented base pentagon translated by
+ * BASIS·(x+y, -x+(flavor&1)), so the test needs no curve decode, no
+ * re-projection, and — pentagons being unit-size here — stays well-conditioned
+ * at every resolution.
+ */
+export function cellMarginScaled(px: number, py: number, x: number, y: number, flavor: number): number {
+  const rx = x + y;
+  const ry = -x + (flavor & 1);
+  const tx = BASIS[0] * rx + BASIS[2] * ry;
+  const ty = BASIS[1] * rx + BASIS[3] * ry;
+  const pent = FLAVOR_PENTAGONS[flavor];
+  let margin = Infinity;
+  for (let i = 0; i < 5; i++) {
+    const j = i === 4 ? 0 : i + 1;
+    const v1x = pent[i * 2] + tx;
+    const v1y = pent[i * 2 + 1] + ty;
+    const v2x = pent[j * 2] + tx;
+    const v2y = pent[j * 2 + 1] + ty;
+    // (v1 - v2) × (p - v1): < 0 ⇒ strictly outside this edge
+    const cross = (v1x - v2x) * (py - v1y) - (v1y - v2y) * (px - v1x);
+    if (cross < margin) margin = cross;
+  }
+  return margin;
+}
+
 /**
  * Get pentagon vertices for a cell.
  *

modules/lattice/compat.ts

@@ -228,7 +228,7 @@ export function compatTripleToS(t: Triple, resolution: number, orientation: Orie
     raw = {x: raw.z, y: raw.y, z: raw.x};
   }
   const ab = tripleToAB(raw);
-  const sGeo = axiomTargetToS(ORIGINAL, ab.a, ab.b, resolution, AXIOM_W, true)[0];
+  const sGeo = axiomTargetToS(ORIGINAL, ab.a, ab.b, resolution, AXIOM_W)[0];
   const digits = digitsOf(sGeo, resolution);
   inverseShift(digits, rec.invertJ, rec.flipIJ);
   const v = packDigits(digits);

modules/lattice/curve.ts

@@ -7,17 +7,43 @@
 // (sToCell / sToTriple) and triple.ts (tripleToS).
 
 import type {IJ} from '../core/coordinate-systems';
-import type {Orientation} from './types';
-import {sumPointToS} from './lsystem';
+import type {Orientation, Triple} from './types';
+import {tripleToSLattice} from './lsystem';
+import {POW2} from './lsystem/tables';
 
 /**
- * Fractional IJ point -> curve position `s` of the containing cell, by direct
- * L-system descent. The IJ plane maps onto the L-system's corner-sum frame by
- * the exact affine map target = (12*(i+j), -12*j) — derived by matching cell
- * centroids across the two frames (parity-0 cell (x,y,z): centroid (x+y+1/3,
- * -x+1/3) in IJ, corner sum (12y+8, 12x-4); parity-1: (x+y-1/3, -x+2/3) and
- * (12y+4, 12x-8); both fit sum = (12(i+j), -12j) exactly), and validated
- * against the old-engine discretization over all resolutions and orientations.
+ * Locate the lattice triangle containing a fractional IJ point, as a triple.
+ *
+ * The triples tile the IJ plane as triangles: the unit square (m, n) =
+ * (floor(i), floor(j)) splits along the diagonal u+v = 1 into a lower triangle
+ * (the parity-0 cell (-n, m+n, -m), centroid (m+1/3, n+1/3)) and an upper
+ * triangle (the parity-1 cell (-n, m+n+1, -m), centroid (m+2/3, n+2/3)) — the
+ * centroid correspondences were derived from the exact IJ <-> corner-sum
+ * affine map target = (12*(i+j), -12*j) and validated against the old-engine
+ * discretization over all resolutions and orientations. Point location is two floors + one
+ * diagonal comparison. Points exactly on a triangle edge have no unique cell;
+ * the >= tie-break below is the fixed convention.
+ *
+ * The result is clamped into quintant bounds (m >= 0, n >= 0, m+n+parity <=
+ * maxRow, equivalent to tripleInBounds): a point slightly outside the quintant
+ * (as the estimate path can produce near quintant edges) must still map to a
+ * valid cell for the exact encode.
  */
-export const IJToS = (ij: IJ, resolution: number, orientation: Orientation = 'uv'): bigint =>
-  sumPointToS(12 * (ij[0] + ij[1]), -12 * ij[1], resolution, orientation);
+export function roundToTriple(ij: IJ, resolution: number): Triple {
+  const maxRow = POW2[resolution] - 1;
+  let m = Math.floor(ij[0]);
+  let n = Math.floor(ij[1]);
+  let parity = ij[0] - m + (ij[1] - n) >= 1 ? 1 : 0;
+  if (m < 0) m = 0;
+  if (n < 0) n = 0;
+  if (m + n + parity > maxRow) {
+    parity = 0;
+    if (m + n > maxRow) {
+      const over = m + n - maxRow;
+      const dm = Math.min(m, over);
+      m -= dm;
+      n -= over - dm;
+    }
+  }
+  return {x: -n, y: m + n + parity, z: -m};
+}

modules/lattice/index.ts

@@ -2,28 +2,21 @@
 // SPDX-License-Identifier: Apache-2.0
 // Copyright (c) A5 contributors
 
-// The canonical A5 curve is currently the ORIGINAL construction (compat.ts):
-// the two-motif quaternary L-system with the shiftDigits recode on top, so
-// cell IDs remain bit-identical to previous releases. The non-self-intersecting
-// L-system curve (lsystem.ts / curve.ts) powers the machinery underneath and is
-// fully implemented and pinned by fixtures (tests/lattice/lsystem.test.ts);
-// making it canonical is a planned follow-up — a breaking change of all cell
-// IDs that swaps the exports below to lsystem.ts/curve.ts and regenerates the
-// fixtures.
+// The canonical A5 curve is the non-self-intersecting L-system curve
+// (lsystem/ + curve.ts): point location via roundToTriple, s <-> cell mappings via
+// sToCell / sToTriple / tripleToS. This is a breaking change from previous
+// releases — cell IDs differ from the original construction. The original curve
+// remains available bit-for-bit via the compat* exports below for migration.
 
 export type {Orientation} from './types';
 
-export {
-  compatSToCell as sToCell,
-  compatSToTriple as sToTriple,
-  compatTripleToS as tripleToS,
-  compatIJToS as IJToS
-} from './compat';
+export {roundToTriple} from './curve';
+export {sToCell, sToTriple} from './lsystem';
 export type {Cell} from './lsystem';
 
 export type {Triple} from './triple';
-export {tripleParity, tripleInBounds} from './triple';
+export {tripleParity, tripleInBounds, tripleFlavor, tripleToS} from './triple';
 
-// Also exported under their own names, so the old-curve behavior stays pinned
-// explicitly (tests/lattice/compat.test.ts) across the future canonical swap.
+// The ORIGINAL (pre-L-system) curve, bit-for-bit, for the migration path —
+// same cells, same pentagon flavors, old visiting order (tests/lattice/compat.test.ts).
 export {compatSToCell, compatSToTriple, compatTripleToS, compatIJToS} from './compat';