Precision and determinism

A geo. operation gives the same digits on the interpreter, a native build and a wasm32 build, and a query gives the same records the code would keep. Floats appear only where PostGIS stores and searches:

No floats in Jadpo's own code

Karney's method needs additions, multiplications, divisions, square roots, sin, cos and atan2. All of it runs on whole numbers far wider than the result, and is rounded once at the end. A float's sin differs between platforms' maths libraries, so a float version could give different digits on wasm32.

Each result is rounded half to even from the exact binary value. It isn't correctly rounded, since the series are finite, but it is the same digits on every target. The checker and the runtime share the code, so a polygon literal is read by the same code the runtime runs.

Coordinates are doubles' shortest decimals

A Latitude or a Longitude, as a value, a field or a point's part, is rounded to the shortest decimal of the nearest IEEE double as it is built, before its range is checked, as a Vector's items are. That is the same on every target. Every coordinate with up to 15 significant digits, which is all a person writes and every GPS gives, is unchanged.

Rounding can't take a value in range out of it, since ±90 and ±180 are doubles. A double whose shortest decimal has more than 38 places (below 10⁻²³ of a degree) is rounded half to even to 38. A literal is checked as written.

Checked against reference answers

geo.distance gives the reference distance of Karney's published test set of 10,000 geodesics, rounded half to even to the micrometre. Before rounding, every distance is within 75 femtometres of the reference. geo.bearing and geo.area are checked against the same set, and every predicate, centroid and box against PostGIS's answers on 169 shapes: small, continental and pole-holding polygons, polygons across the antimeridian and the equator, holes, concave shapes, lines and multi-part shapes. Every yes or no is PostGIS's, except where ST_Covers parts from geo.contains (contains, within and intersects). The gate runs both checks natively and under WebAssembly, which give the same answers. The design notes give the method and its tolerances.

Predicates are exact

Whether a point is inside a shape, or two shapes meet, is decided exactly, by whole numbers, with no tolerance. So a point exactly on an edge, a corner, a meridian or the equator is found on it, and inside includes the edge. Where three points lie exactly on one great circle, the tie is broken by a fixed rule (Simulation of Simplicity, as Google's S2 library does). No float decides anything.

What it costs

These figures are measured natively on an Apple M1 Pro (cargo run --release -p jadpo-value --example geodesics -- bench).

• Measures: a distance along a meridian takes 9 µs; 2 km across London 31 µs; London to Paris 42 µs; London to Sydney 58 µs; nearly antipodal 64 µs. A bearing takes 47 µs, an 11-point line's length 315 µs, and 71 µs for the 1° square's area (meridians and the equator are quick). A 100-point patch's area takes 4.5 ms, about 45 µs an edge.
• Under WebAssembly: checking 500 of the test set's geodesics takes 0.30 s natively and 0.83 s under Node's WebAssembly, about 2.8 times, so a distance there takes about 110 µs, and 10 ms of CPU measures about 90 distances.
• Predicates and shapes: a point in a 100-point patch takes 0.44 ms, in a 1,000-point one 4.3 ms. A 100-point patch in another takes 30 ms, since each piece of every edge is classified, about n² in the points. Two 100-point patches intersecting take 0.9 ms, within(distance:) 33 µs, a 100-point patch's centroid 0.37 ms and its box 1.2 ms. Turning each corner into a direction (two sines and cosines) is most of a small polygon's cost.

A request measures a handful of things. A batch measuring millions, or the area of a 10,000-point outline (about half a second), would notice. No time limit is enforced, as none is elsewhere; the figures are the baseline that tuning measures against.

Where code and the database could disagree

They could disagree on a point on, or within a hair of, a polygon's edge, and on a distance within a hair of a radius: PostGIS decides in doubles with its own tolerances, and Jadpo exactly. ST_Covers also parts from geo.contains on four kinds of polygon. So a query's condition uses PostGIS only to find the rows, through the index, a millimetre wider than the condition, and the runtime decides each row with the same code a body runs (Queries). The index is used, and the answer is the code's.

The cost is the code's test on each row found: about 40 µs a distance natively, and 0.44 ms for a point in a 100-point patch. So a radius holding a thousand bakeries spends some 40 ms sorting them by distance.