BUG #19597: getQuadrant: impossible case is reachable

The following bug has been logged on the website:

Bug reference: 19597
Logged by: Michael Malis
Email address: malis(at)pgrust(dot)com
PostgreSQL version: 18.3
Operating system: Debian 18.3-1.pgdg13+1
Description:

getQuadrant has four arms covering the quadrants and an elog(ERROR,
"getQuadrant: impossible case") fallthrough. The arms are exhaustive only if
the above/below/horizontal predicates are exhaustive. They are not: the
three comparisons are written in three algebraically different forms, and
near a power-of-two boundary — where the gap between adjacent doubles
doubles — adding EPSILON rounds up on one side and vanishes on the other. A
point can then be neither above, below, nor level with the centroid, and the
"impossible" branch executes.

Reproducer (runnable against stock PostgreSQL 18.3)
---------------------------------------------------
CREATE TABLE gq(p point);
-- 4000 identical points: their mean (the quad centroid) is exactly this
value,
-- which is nextafter(2^34, -inf)
INSERT INTO gq SELECT point('17179869183.999998','17179869183.999998')
FROM generate_series(1,4000);
CREATE INDEX gqx ON gq USING spgist(p);
-- probe at 2^34, one representable step away
INSERT INTO gq SELECT point('17179869184.0','17179869184.0')
FROM generate_series(1,400);
ERROR: getQuadrant: impossible case

The identical-point fill is what forces the computed centroid onto the
chosen value — spg_quad_picksplit uses the mean of the split set, so the
centroid cannot simply be inserted.

Expected vs. actual
-------------------
- Expected: the insert succeeds, or fails with a meaningful user-facing
error. A branch labelled "impossible case" should not be reachable from
well-formed finite input.
- Actual: ERROR: getQuadrant: impossible case. The backend survives (this is
a catchable error, not a crash), but the index operation fails and the
message is an internal invariant leaking to the user.

Mechanism, with file:line into the 18.3 source
----------------------------------------------
src/backend/access/spgist/spgquadtreeproc.c, getQuadrant (function at 55):
57: if ((SPTEST(point_above, tst, centroid) || SPTEST(point_horiz, tst,
centroid)) && ...
63: if (SPTEST(point_below, tst, centroid) && ...
68: if ((SPTEST(point_below, tst, centroid) || SPTEST(point_horiz, tst,
centroid)) && ...
73: if (SPTEST(point_above, tst, centroid) && SPTEST(point_left, tst,
centroid))
77: elog(ERROR, "getQuadrant: impossible case");

point_above/point_below/point_horiz reduce to FPgt/FPlt/FPeq on the y
coordinate. src/include/utils/geo_decls.h:
41: #define EPSILON 1.0E-06
47: FPeq(A,B) { return A == B || fabs(A - B) <= EPSILON; }
59: FPlt(A,B) { return A + EPSILON < B; }
71: FPgt(A,B) { return A > B + EPSILON; }

FPeq tests the rounded difference; FPlt and FPgt test rounded sums. Over
exact reals these partition the line; in floating point they do not. Just
below 2^34 the spacing between doubles is 2^-19 ≈ 1.907e-6 (> EPSILON, so
adding EPSILON rounds up a step); just above it is 2^-18 ≈ 3.815e-6 (>
2·EPSILON, so adding EPSILON rounds back down). For a = 2^34, b =
nextafter(2^34, -inf) all three are false — verified on 18.3's own
arithmetic:

expression value
----------- --------------------------------------------
a - b 1.9073486328125e-06 (> EPSILON ⇒ FPeq false)
a > b + eps false (b + eps rounds up to exactly a)
a + eps < b false (a + eps rounds back to a)

With the y trichotomy failing, no arm matches regardless of the x result,
and line 77 executes.