# 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:

| Where | Arithmetic | Same on every target? |
|---|---|---|
| `geo.` operations in code: interpreter, native, wasm32 | Jadpo's own whole-number fixed point, 256 bits with 136 after the point (about 41 decimal places), its `sqrt` a whole square root and its `sin`, `cos` and `atan2` the series of the standard library's decimal maths ([Numbers](numbers.md)); results rounded half to even to the result type's scale | Yes: integer arithmetic only, no platform maths library |
| Reading GeoJSON | The double JSON gives, held as its shortest decimal | Yes |
| Coordinates built in code | Rounded to the nearest double's shortest decimal when built | Yes |
| Storage | PostGIS doubles; exact, given the row above | Yes |
| Query conditions and order | PostGIS's doubles find the rows, through the index and a millimetre wider than the condition; Jadpo's own code decides each and sorts them | Yes: the code's answer |

## 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`](geometry-operations.md#contains-within-and-intersects)).
The gate runs both checks natively and under WebAssembly, which give the
same answers. The [design notes](https://github.com/jadpo-dev/jadpo/blob/main/docs/language/notes/20-geospatial.md#how-the-geodesics-are-checked)
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](geospatial-queries.md)). 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.
