From a7e0237d0bfd4909c0c4e0237013efe0cd071dd4 Mon Sep 17 00:00:00 2001
From: Heikki Linnakangas
Date: Mon, 7 Sep 2020 12:33:36 +0300
Subject: [PATCH v15 4/4] Map negative values better.

src/backend/access/gist/gistproc.c  169 +++++++++++++++++++++
1 file changed, 121 insertions(+), 48 deletions()
diff git a/src/backend/access/gist/gistproc.c b/src/backend/access/gist/gistproc.c
index 4ed9f46c9bb..7ca7eda84b0 100644
 a/src/backend/access/gist/gistproc.c
+++ b/src/backend/access/gist/gistproc.c
@@ 31,10 +31,11 @@ static bool gist_box_leaf_consistent(BOX *key, BOX *query,
StrategyNumber strategy);
static bool rtree_internal_consistent(BOX *key, BOX *query,
StrategyNumber strategy);
static int64 part_bits32_by2(uint32 x);
static int64 interleave_bits32(uint32 x, uint32 y);
static inline uint64 point_zorder_internal(Point *p);
static int gist_bbox_fastcmp(Datum x, Datum y, SortSupport ssup);
+static uint64 part_bits32_by2(uint32 x);
+static uint64 interleave_bits32(uint32 x, uint32 y);
+static uint32 ieee_float32_to_uint32(float f);
+static uint64 point_zorder_internal(Point *p);
+static int gist_bbox_fastcmp(Datum x, Datum y, SortSupport ssup);
/* Minimum accepted ratio of split */
@@ 1549,75 +1550,147 @@ gist_poly_distance(PG_FUNCTION_ARGS)
/* Zorder routines */
/* Interleave 32 bits with zeroes */
static int64
+static uint64
part_bits32_by2(uint32 x)
{
uint64 n = x;
n = (n  (n << 16)) & UINT64CONST(0x0000FFFF0000FFFF);
 n = (n  (n << 8)) & UINT64CONST(0x00FF00FF00FF00FF);
 n = (n  (n << 4)) & UINT64CONST(0x0F0F0F0F0F0F0F0F);
 n = (n  (n << 2)) & UINT64CONST(0x3333333333333333);
 n = (n  (n << 1)) & UINT64CONST(0x5555555555555555);
+ n = (n  (n << 8)) & UINT64CONST(0x00FF00FF00FF00FF);
+ n = (n  (n << 4)) & UINT64CONST(0x0F0F0F0F0F0F0F0F);
+ n = (n  (n << 2)) & UINT64CONST(0x3333333333333333);
+ n = (n  (n << 1)) & UINT64CONST(0x5555555555555555);
return n;
}
/*
 * Compute Zorder for integers. Also called Morton code.
+ * Compute Zorder for a point
+ *
+ * Map a twodimensional point to a single integer, in a way that preserves
+ * locality. Points that are close in the twodimensional space are mapped to
+ * integer that are not far from each other. We do that by interleaving the
+ * bits in the X and Y components, this is called a Zorder or Morton Code.
+ *
+ * A Morton Code is normally defined only for integers, but the X and Y values
+ * of a point are floating point. We expect floats to be in IEEE format, and
+ * the sort order of IEEE floats is mostly correlated to the binary sort order
+ * of the bits reinterpreted as an int. It isn't in some special cases, but
+ * for this use case we don't really care about that, we're just trying to
+ * encourage locality.
*/
static int64
+
+static uint64
+point_zorder_internal(Point *p)
+{
+ return interleave_bits32(ieee_float32_to_uint32(p>x),
+ ieee_float32_to_uint32(p>y));
+}
+
+static uint64
interleave_bits32(uint32 x, uint32 y)
{
return part_bits32_by2(x)  (part_bits32_by2(y) << 1);
}
/* Compute Zorder for Point */
static inline uint64
point_zorder_internal(Point *p)
+/*
+ * Convert a 32bit IEEE float to uint32 in a way that preserves the ordering.
+ */
+static uint32
+ieee_float32_to_uint32(float f)
{
 /*
 * In this function we need to compute Morton codes for floating point
 * components p>x and p>y. But Morton codes are defined only for
 * integers.
 * We expect floats to be in IEEE format, and the sort order of IEEE
 * floats is mostly correlated to the binary sort order of the bits
 * reinterpreted as an int. It isn't in some special cases, but for this
 * use case we don't really care about that, we're just trying to
 * encourage locality.
 * There is a big jump in integer value (whether signed or
 * unsigned) as you cross from positive to negative floats, and then the
 * sort order is reversed. This can have negative effect on searches when
 * query window touches many quadrants simultaneously. In worst case this
 * searches can be x4 more costly.
 * We generate a Morton code that interleaves the bits of N integers
 * to produce a single integer that preserves locality: things that were
 * close in the N dimensional space are close in the resulting integer.
+ /*
+ *
+ * IEEE 754 floating point format
+ * 
+ *
+ * IEEE 754 floating point numbers have this format:
+ *
+ * exponent (8 bits)
+ * 
+ * s eeeeeeee mmmmmmmmmmmmmmmmmmmmmmm
+ *  
+ * sign mantissa (23 bits)
+ *
+ * Infinity has all bits in the exponent set and the mantissa is
+ * allzeros. Negative infinity is the same but with the sign bit set.
+ *
+ * NaNs are represented with all bits in the exponent set, and the least
+ * significant bit in the mantissa also set. The rest of the mantissa bits
+ * can be used to distinguish different kinds of NaNs.
+ *
+ * The IEEE format has the nice property that when you take the bit
+ * representation and interpret it as an integer, the order is preserved,
+ * except for the sign. That holds for the +Infinity values too.
+ *
+ * Mapping to uint32
+ * 
+ *
+ * In order to have a smooth transition from negative to positive numbers,
+ * we map floats to unsigned integers like this:
+ *
+ * x < 0 to range 07FFFFFFF
+ * x = 0 to value 8000000 (both positive and negative zero)
+ * x > 0 to range 8000001FFFFFFFF
+ *
+ * We don't care to distinguish different kind of NaNs, so they are all
+ * mapped to the same arbitrary value, FFFFFFFF. Because of the IEEE bit
+ * representation of NaNs, there aren't actually any nonNaN values that
+ * would be mapped to FFFFFFFF, there is actually a range of unused values
+ * on both ends of the uint32 space.
*/
 union {
 float f;
 uint32 i;
 } a,b;
 a.f = p>x;
 b.f = p>y;
 if (isnan(a.f))
 a.i = PG_INT32_MAX;
 if (isnan(b.f))
 b.i = PG_INT32_MAX;
 return interleave_bits32(a.i, b.i);
+ if (isnan(f))
+ return 0xFFFFFFFF;
+ else
+ {
+ union
+ {
+ float f;
+ uint32 i;
+ } u;
+
+ u.f = f;
+
+ /* Check the sign bit */
+ if ((u.i & 0x80000000) != 0)
+ {
+ /*
+ * Map the negative value to range 07FFFFFFF. This flips the sign
+ * bit to 0 in the same instruction.
+ */
+ Assert(f < 0);
+
+ u.i ^= 0xFFFFFFFF;
+ }
+ else
+ {
+ /* Map the positive value (or 0) to range 80000000FFFFFFFF */
+ u.i = 0x80000000;
+ }
+
+ return u.i;
+ }
}
static int
gist_bbox_fastcmp(Datum x, Datum y, SortSupport ssup)
{
 Point *p1 = &(DatumGetBoxP(x)>low);
 Point *p2 = &(DatumGetBoxP(y)>low);
 uint64 z1 = point_zorder_internal(p1);
 uint64 z2 = point_zorder_internal(p2);

 return z1 == z2 ? 0 : z1 > z2 ? 1 : 1;
+ Point *p1 = &(DatumGetBoxP(x)>low);
+ Point *p2 = &(DatumGetBoxP(y)>low);
+ uint64 z1 = point_zorder_internal(p1);
+ uint64 z2 = point_zorder_internal(p2);
+
+ if (z1 > z2)
+ return 1;
+ else if (z1 < z2)
+ return 1;
+ else
+ return 0;
}
+/*
+ * Sort support routine for fast GiST index build by sorting.
+ */
Datum
gist_point_sortsupport(PG_FUNCTION_ARGS)
{

2.20.1