/src/postgres/src/backend/access/gist/gistproc.c

Source
/*-------------------------------------------------------------------------
 *
 * gistproc.c
 *    Support procedures for GiSTs over 2-D objects (boxes, polygons, circles,
 *    points).
 *
 * This gives R-tree behavior, with Guttman's poly-time split algorithm.
 *
 *
 * Portions Copyright (c) 1996-2026, PostgreSQL Global Development Group
 * Portions Copyright (c) 1994, Regents of the University of California
 *
 * IDENTIFICATION
 *  src/backend/access/gist/gistproc.c
 *
 *-------------------------------------------------------------------------
 */
#include "postgres.h"

#include <math.h>

#include "access/gist.h"
#include "access/stratnum.h"
#include "utils/float.h"
#include "utils/fmgrprotos.h"
#include "utils/geo_decls.h"
#include "utils/sortsupport.h"


static bool gist_box_leaf_consistent(BOX *key, BOX *query,
                   StrategyNumber strategy);
static bool rtree_internal_consistent(BOX *key, BOX *query,
                    StrategyNumber strategy);

static uint64 point_zorder_internal(float4 x, float4 y);
static uint64 part_bits32_by2(uint32 x);
static uint32 ieee_float32_to_uint32(float f);
static int  gist_bbox_zorder_cmp(Datum a, Datum b, SortSupport ssup);
static Datum gist_bbox_zorder_abbrev_convert(Datum original, SortSupport ssup);
static bool gist_bbox_zorder_abbrev_abort(int memtupcount, SortSupport ssup);


/* Minimum accepted ratio of split */
#define LIMIT_RATIO 0.3


/**************************************************
 * Box ops
 **************************************************/

/*
 * Calculates union of two boxes, a and b. The result is stored in *n.
 */
static void
rt_box_union(BOX *n, const BOX *a, const BOX *b)
{
  n->high.x = float8_max(a->high.x, b->high.x);
  n->high.y = float8_max(a->high.y, b->high.y);
  n->low.x = float8_min(a->low.x, b->low.x);
  n->low.y = float8_min(a->low.y, b->low.y);
}

/*
 * Size of a BOX for penalty-calculation purposes.
 * The result can be +Infinity, but not NaN.
 */
static float8
size_box(const BOX *box)
{
  /*
   * Check for zero-width cases.  Note that we define the size of a zero-
   * by-infinity box as zero.  It's important to special-case this somehow,
   * as naively multiplying infinity by zero will produce NaN.
   *
   * The less-than cases should not happen, but if they do, say "zero".
   */
  if (float8_le(box->high.x, box->low.x) ||
    float8_le(box->high.y, box->low.y))
    return 0.0;

  /*
   * We treat NaN as larger than +Infinity, so any distance involving a NaN
   * and a non-NaN is infinite.  Note the previous check eliminated the
   * possibility that the low fields are NaNs.
   */
  if (isnan(box->high.x) || isnan(box->high.y))
    return get_float8_infinity();
  return float8_mul(float8_mi(box->high.x, box->low.x),
            float8_mi(box->high.y, box->low.y));
}

/*
 * Return amount by which the union of the two boxes is larger than
 * the original BOX's area.  The result can be +Infinity, but not NaN.
 */
static float8
box_penalty(const BOX *original, const BOX *new)
{
  BOX     unionbox;

  rt_box_union(&unionbox, original, new);
  return float8_mi(size_box(&unionbox), size_box(original));
}

/*
 * The GiST Consistent method for boxes
 *
 * Should return false if for all data items x below entry,
 * the predicate x op query must be false, where op is the oper
 * corresponding to strategy in the pg_amop table.
 */
Datum
gist_box_consistent(PG_FUNCTION_ARGS)
{
  GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  BOX      *query = PG_GETARG_BOX_P(1);
  StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
#ifdef NOT_USED
  Oid     subtype = PG_GETARG_OID(3);
#endif
  bool     *recheck = (bool *) PG_GETARG_POINTER(4);

  /* All cases served by this function are exact */
  *recheck = false;

  if (DatumGetBoxP(entry->key) == NULL || query == NULL)
    PG_RETURN_BOOL(false);

  /*
   * if entry is not leaf, use rtree_internal_consistent, else use
   * gist_box_leaf_consistent
   */
  if (GIST_LEAF(entry))
    PG_RETURN_BOOL(gist_box_leaf_consistent(DatumGetBoxP(entry->key),
                        query,
                        strategy));
  else
    PG_RETURN_BOOL(rtree_internal_consistent(DatumGetBoxP(entry->key),
                         query,
                         strategy));
}

/*
 * Increase BOX b to include addon.
 */
static void
adjustBox(BOX *b, const BOX *addon)
{
  if (float8_lt(b->high.x, addon->high.x))
    b->high.x = addon->high.x;
  if (float8_gt(b->low.x, addon->low.x))
    b->low.x = addon->low.x;
  if (float8_lt(b->high.y, addon->high.y))
    b->high.y = addon->high.y;
  if (float8_gt(b->low.y, addon->low.y))
    b->low.y = addon->low.y;
}

/*
 * The GiST Union method for boxes
 *
 * returns the minimal bounding box that encloses all the entries in entryvec
 */
Datum
gist_box_union(PG_FUNCTION_ARGS)
{
  GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
  int      *sizep = (int *) PG_GETARG_POINTER(1);
  int     numranges,
        i;
  BOX      *cur,
         *pageunion;

  numranges = entryvec->n;
  pageunion = palloc_object(BOX);
  cur = DatumGetBoxP(entryvec->vector[0].key);
  memcpy(pageunion, cur, sizeof(BOX));

  for (i = 1; i < numranges; i++)
  {
    cur = DatumGetBoxP(entryvec->vector[i].key);
    adjustBox(pageunion, cur);
  }
  *sizep = sizeof(BOX);

  PG_RETURN_POINTER(pageunion);
}

/*
 * We store boxes as boxes in GiST indexes, so we do not need
 * compress, decompress, or fetch functions.
 */

/*
 * The GiST Penalty method for boxes (also used for points)
 *
 * As in the R-tree paper, we use change in area as our penalty metric
 */
Datum
gist_box_penalty(PG_FUNCTION_ARGS)
{
  GISTENTRY  *origentry = (GISTENTRY *) PG_GETARG_POINTER(0);
  GISTENTRY  *newentry = (GISTENTRY *) PG_GETARG_POINTER(1);
  float    *result = (float *) PG_GETARG_POINTER(2);
  BOX      *origbox = DatumGetBoxP(origentry->key);
  BOX      *newbox = DatumGetBoxP(newentry->key);

  *result = (float) box_penalty(origbox, newbox);
  PG_RETURN_POINTER(result);
}

/*
 * Trivial split: half of entries will be placed on one page
 * and another half - to another
 */
static void
fallbackSplit(GistEntryVector *entryvec, GIST_SPLITVEC *v)
{
  OffsetNumber i,
        maxoff;
  BOX      *unionL = NULL,
         *unionR = NULL;
  int     nbytes;

  maxoff = entryvec->n - 1;

  nbytes = (maxoff + 2) * sizeof(OffsetNumber);
  v->spl_left = (OffsetNumber *) palloc(nbytes);
  v->spl_right = (OffsetNumber *) palloc(nbytes);
  v->spl_nleft = v->spl_nright = 0;

  for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
  {
    BOX      *cur = DatumGetBoxP(entryvec->vector[i].key);

    if (i <= (maxoff - FirstOffsetNumber + 1) / 2)
    {
      v->spl_left[v->spl_nleft] = i;
      if (unionL == NULL)
      {
        unionL = palloc_object(BOX);
        *unionL = *cur;
      }
      else
        adjustBox(unionL, cur);

      v->spl_nleft++;
    }
    else
    {
      v->spl_right[v->spl_nright] = i;
      if (unionR == NULL)
      {
        unionR = palloc_object(BOX);
        *unionR = *cur;
      }
      else
        adjustBox(unionR, cur);

      v->spl_nright++;
    }
  }

  v->spl_ldatum = BoxPGetDatum(unionL);
  v->spl_rdatum = BoxPGetDatum(unionR);
}

/*
 * Represents information about an entry that can be placed to either group
 * without affecting overlap over selected axis ("common entry").
 */
typedef struct
{
  /* Index of entry in the initial array */
  int     index;
  /* Delta between penalties of entry insertion into different groups */
  float8    delta;
} CommonEntry;

/*
 * Context for g_box_consider_split. Contains information about currently
 * selected split and some general information.
 */
typedef struct
{
  int     entriesCount; /* total number of entries being split */
  BOX     boundingBox;  /* minimum bounding box across all entries */

  /* Information about currently selected split follows */

  bool    first;      /* true if no split was selected yet */

  float8    leftUpper;    /* upper bound of left interval */
  float8    rightLower;   /* lower bound of right interval */

  float4    ratio;
  float4    overlap;
  int     dim;      /* axis of this split */
  float8    range;      /* width of general MBR projection to the
                 * selected axis */
} ConsiderSplitContext;

/*
 * Interval represents projection of box to axis.
 */
typedef struct
{
  float8    lower,
        upper;
} SplitInterval;

/*
 * Interval comparison function by lower bound of the interval;
 */
static int
interval_cmp_lower(const void *i1, const void *i2)
{
  float8    lower1 = ((const SplitInterval *) i1)->lower,
        lower2 = ((const SplitInterval *) i2)->lower;

  return float8_cmp_internal(lower1, lower2);
}

/*
 * Interval comparison function by upper bound of the interval;
 */
static int
interval_cmp_upper(const void *i1, const void *i2)
{
  float8    upper1 = ((const SplitInterval *) i1)->upper,
        upper2 = ((const SplitInterval *) i2)->upper;

  return float8_cmp_internal(upper1, upper2);
}

/*
 * Replace negative (or NaN) value with zero.
 */
static inline float
non_negative(float val)
{
  if (val >= 0.0f)
    return val;
  else
    return 0.0f;
}

/*
 * Consider replacement of currently selected split with the better one.
 */
static inline void
g_box_consider_split(ConsiderSplitContext *context, int dimNum,
           float8 rightLower, int minLeftCount,
           float8 leftUpper, int maxLeftCount)
{
  int     leftCount,
        rightCount;
  float4    ratio,
        overlap;
  float8    range;

  /*
   * Calculate entries distribution ratio assuming most uniform distribution
   * of common entries.
   */
  if (minLeftCount >= (context->entriesCount + 1) / 2)
  {
    leftCount = minLeftCount;
  }
  else
  {
    if (maxLeftCount <= context->entriesCount / 2)
      leftCount = maxLeftCount;
    else
      leftCount = context->entriesCount / 2;
  }
  rightCount = context->entriesCount - leftCount;

  /*
   * Ratio of split - quotient between size of lesser group and total
   * entries count.
   */
  ratio = float4_div(Min(leftCount, rightCount), context->entriesCount);

  if (ratio > LIMIT_RATIO)
  {
    bool    selectthis = false;

    /*
     * The ratio is acceptable, so compare current split with previously
     * selected one. Between splits of one dimension we search for minimal
     * overlap (allowing negative values) and minimal ration (between same
     * overlaps. We switch dimension if find less overlap (non-negative)
     * or less range with same overlap.
     */
    if (dimNum == 0)
      range = float8_mi(context->boundingBox.high.x,
                context->boundingBox.low.x);
    else
      range = float8_mi(context->boundingBox.high.y,
                context->boundingBox.low.y);

    overlap = float8_div(float8_mi(leftUpper, rightLower), range);

    /* If there is no previous selection, select this */
    if (context->first)
      selectthis = true;
    else if (context->dim == dimNum)
    {
      /*
       * Within the same dimension, choose the new split if it has a
       * smaller overlap, or same overlap but better ratio.
       */
      if (overlap < context->overlap ||
        (overlap == context->overlap && ratio > context->ratio))
        selectthis = true;
    }
    else
    {
      /*
       * Across dimensions, choose the new split if it has a smaller
       * *non-negative* overlap, or same *non-negative* overlap but
       * bigger range. This condition differs from the one described in
       * the article. On the datasets where leaf MBRs don't overlap
       * themselves, non-overlapping splits (i.e. splits which have zero
       * *non-negative* overlap) are frequently possible. In this case
       * splits tends to be along one dimension, because most distant
       * non-overlapping splits (i.e. having lowest negative overlap)
       * appears to be in the same dimension as in the previous split.
       * Therefore MBRs appear to be very prolonged along another
       * dimension, which leads to bad search performance. Using range
       * as the second split criteria makes MBRs more quadratic. Using
       * *non-negative* overlap instead of overlap as the first split
       * criteria gives to range criteria a chance to matter, because
       * non-overlapping splits are equivalent in this criteria.
       */
      if (non_negative(overlap) < non_negative(context->overlap) ||
        (range > context->range &&
         non_negative(overlap) <= non_negative(context->overlap)))
        selectthis = true;
    }

    if (selectthis)
    {
      /* save information about selected split */
      context->first = false;
      context->ratio = ratio;
      context->range = range;
      context->overlap = overlap;
      context->rightLower = rightLower;
      context->leftUpper = leftUpper;
      context->dim = dimNum;
    }
  }
}

/*
 * Compare common entries by their deltas.
 */
static int
common_entry_cmp(const void *i1, const void *i2)
{
  float8    delta1 = ((const CommonEntry *) i1)->delta,
        delta2 = ((const CommonEntry *) i2)->delta;

  return float8_cmp_internal(delta1, delta2);
}

/*
 * --------------------------------------------------------------------------
 * Double sorting split algorithm. This is used for both boxes and points.
 *
 * The algorithm finds split of boxes by considering splits along each axis.
 * Each entry is first projected as an interval on the X-axis, and different
 * ways to split the intervals into two groups are considered, trying to
 * minimize the overlap of the groups. Then the same is repeated for the
 * Y-axis, and the overall best split is chosen. The quality of a split is
 * determined by overlap along that axis and some other criteria (see
 * g_box_consider_split).
 *
 * After that, all the entries are divided into three groups:
 *
 * 1) Entries which should be placed to the left group
 * 2) Entries which should be placed to the right group
 * 3) "Common entries" which can be placed to any of groups without affecting
 *    of overlap along selected axis.
 *
 * The common entries are distributed by minimizing penalty.
 *
 * For details see:
 * "A new double sorting-based node splitting algorithm for R-tree", A. Korotkov
 * http://syrcose.ispras.ru/2011/files/SYRCoSE2011_Proceedings.pdf#page=36
 * --------------------------------------------------------------------------
 */
Datum
gist_box_picksplit(PG_FUNCTION_ARGS)
{
  GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
  GIST_SPLITVEC *v = (GIST_SPLITVEC *) PG_GETARG_POINTER(1);
  OffsetNumber i,
        maxoff;
  ConsiderSplitContext context;
  BOX      *box,
         *leftBox,
         *rightBox;
  int     dim,
        commonEntriesCount;
  SplitInterval *intervalsLower,
         *intervalsUpper;
  CommonEntry *commonEntries;
  int     nentries;

  memset(&context, 0, sizeof(ConsiderSplitContext));

  maxoff = entryvec->n - 1;
  nentries = context.entriesCount = maxoff - FirstOffsetNumber + 1;

  /* Allocate arrays for intervals along axes */
  intervalsLower = palloc_array(SplitInterval, nentries);
  intervalsUpper = palloc_array(SplitInterval, nentries);

  /*
   * Calculate the overall minimum bounding box over all the entries.
   */
  for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
  {
    box = DatumGetBoxP(entryvec->vector[i].key);
    if (i == FirstOffsetNumber)
      context.boundingBox = *box;
    else
      adjustBox(&context.boundingBox, box);
  }

  /*
   * Iterate over axes for optimal split searching.
   */
  context.first = true;   /* nothing selected yet */
  for (dim = 0; dim < 2; dim++)
  {
    float8    leftUpper,
          rightLower;
    int     i1,
          i2;

    /* Project each entry as an interval on the selected axis. */
    for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
    {
      box = DatumGetBoxP(entryvec->vector[i].key);
      if (dim == 0)
      {
        intervalsLower[i - FirstOffsetNumber].lower = box->low.x;
        intervalsLower[i - FirstOffsetNumber].upper = box->high.x;
      }
      else
      {
        intervalsLower[i - FirstOffsetNumber].lower = box->low.y;
        intervalsLower[i - FirstOffsetNumber].upper = box->high.y;
      }
    }

    /*
     * Make two arrays of intervals: one sorted by lower bound and another
     * sorted by upper bound.
     */
    memcpy(intervalsUpper, intervalsLower,
         sizeof(SplitInterval) * nentries);
    qsort(intervalsLower, nentries, sizeof(SplitInterval),
        interval_cmp_lower);
    qsort(intervalsUpper, nentries, sizeof(SplitInterval),
        interval_cmp_upper);

    /*----
     * The goal is to form a left and right interval, so that every entry
     * interval is contained by either left or right interval (or both).
     *
     * For example, with the intervals (0,1), (1,3), (2,3), (2,4):
     *
     * 0 1 2 3 4
     * +-+
     *   +---+
     *     +-+
     *     +---+
     *
     * The left and right intervals are of the form (0,a) and (b,4).
     * We first consider splits where b is the lower bound of an entry.
     * We iterate through all entries, and for each b, calculate the
     * smallest possible a. Then we consider splits where a is the
     * upper bound of an entry, and for each a, calculate the greatest
     * possible b.
     *
     * In the above example, the first loop would consider splits:
     * b=0: (0,1)-(0,4)
     * b=1: (0,1)-(1,4)
     * b=2: (0,3)-(2,4)
     *
     * And the second loop:
     * a=1: (0,1)-(1,4)
     * a=3: (0,3)-(2,4)
     * a=4: (0,4)-(2,4)
     */

    /*
     * Iterate over lower bound of right group, finding smallest possible
     * upper bound of left group.
     */
    i1 = 0;
    i2 = 0;
    rightLower = intervalsLower[i1].lower;
    leftUpper = intervalsUpper[i2].lower;
    while (true)
    {
      /*
       * Find next lower bound of right group.
       */
      while (i1 < nentries &&
           float8_eq(rightLower, intervalsLower[i1].lower))
      {
        if (float8_lt(leftUpper, intervalsLower[i1].upper))
          leftUpper = intervalsLower[i1].upper;
        i1++;
      }
      if (i1 >= nentries)
        break;
      rightLower = intervalsLower[i1].lower;

      /*
       * Find count of intervals which anyway should be placed to the
       * left group.
       */
      while (i2 < nentries &&
           float8_le(intervalsUpper[i2].upper, leftUpper))
        i2++;

      /*
       * Consider found split.
       */
      g_box_consider_split(&context, dim, rightLower, i1, leftUpper, i2);
    }

    /*
     * Iterate over upper bound of left group finding greatest possible
     * lower bound of right group.
     */
    i1 = nentries - 1;
    i2 = nentries - 1;
    rightLower = intervalsLower[i1].upper;
    leftUpper = intervalsUpper[i2].upper;
    while (true)
    {
      /*
       * Find next upper bound of left group.
       */
      while (i2 >= 0 && float8_eq(leftUpper, intervalsUpper[i2].upper))
      {
        if (float8_gt(rightLower, intervalsUpper[i2].lower))
          rightLower = intervalsUpper[i2].lower;
        i2--;
      }
      if (i2 < 0)
        break;
      leftUpper = intervalsUpper[i2].upper;

      /*
       * Find count of intervals which anyway should be placed to the
       * right group.
       */
      while (i1 >= 0 && float8_ge(intervalsLower[i1].lower, rightLower))
        i1--;

      /*
       * Consider found split.
       */
      g_box_consider_split(&context, dim,
                 rightLower, i1 + 1, leftUpper, i2 + 1);
    }
  }

  /*
   * If we failed to find any acceptable splits, use trivial split.
   */
  if (context.first)
  {
    fallbackSplit(entryvec, v);
    PG_RETURN_POINTER(v);
  }

  /*
   * Ok, we have now selected the split across one axis.
   *
   * While considering the splits, we already determined that there will be
   * enough entries in both groups to reach the desired ratio, but we did
   * not memorize which entries go to which group. So determine that now.
   */

  /* Allocate vectors for results */
  v->spl_left = palloc_array(OffsetNumber, nentries);
  v->spl_right = palloc_array(OffsetNumber, nentries);
  v->spl_nleft = 0;
  v->spl_nright = 0;

  /* Allocate bounding boxes of left and right groups */
  leftBox = palloc0_object(BOX);
  rightBox = palloc0_object(BOX);

  /*
   * Allocate an array for "common entries" - entries which can be placed to
   * either group without affecting overlap along selected axis.
   */
  commonEntriesCount = 0;
  commonEntries = palloc_array(CommonEntry, nentries);

  /* Helper macros to place an entry in the left or right group */
#define PLACE_LEFT(box, off)          \
  do {                   \
    if (v->spl_nleft > 0)         \
      adjustBox(leftBox, box);     \
    else                  \
      *leftBox = *(box);         \
    v->spl_left[v->spl_nleft++] = off;    \
  } while(0)

#define PLACE_RIGHT(box, off)         \
  do {                   \
    if (v->spl_nright > 0)         \
      adjustBox(rightBox, box);     \
    else                  \
      *rightBox = *(box);         \
    v->spl_right[v->spl_nright++] = off;  \
  } while(0)

  /*
   * Distribute entries which can be distributed unambiguously, and collect
   * common entries.
   */
  for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
  {
    float8    lower,
          upper;

    /*
     * Get upper and lower bounds along selected axis.
     */
    box = DatumGetBoxP(entryvec->vector[i].key);
    if (context.dim == 0)
    {
      lower = box->low.x;
      upper = box->high.x;
    }
    else
    {
      lower = box->low.y;
      upper = box->high.y;
    }

    if (float8_le(upper, context.leftUpper))
    {
      /* Fits to the left group */
      if (float8_ge(lower, context.rightLower))
      {
        /* Fits also to the right group, so "common entry" */
        commonEntries[commonEntriesCount++].index = i;
      }
      else
      {
        /* Doesn't fit to the right group, so join to the left group */
        PLACE_LEFT(box, i);
      }
    }
    else
    {
      /*
       * Each entry should fit on either left or right group. Since this
       * entry didn't fit on the left group, it better fit in the right
       * group.
       */
      Assert(float8_ge(lower, context.rightLower));

      /* Doesn't fit to the left group, so join to the right group */
      PLACE_RIGHT(box, i);
    }
  }

  /*
   * Distribute "common entries", if any.
   */
  if (commonEntriesCount > 0)
  {
    /*
     * Calculate minimum number of entries that must be placed in both
     * groups, to reach LIMIT_RATIO.
     */
    int     m = ceil(LIMIT_RATIO * nentries);

    /*
     * Calculate delta between penalties of join "common entries" to
     * different groups.
     */
    for (i = 0; i < commonEntriesCount; i++)
    {
      box = DatumGetBoxP(entryvec->vector[commonEntries[i].index].key);
      commonEntries[i].delta = fabs(float8_mi(box_penalty(leftBox, box),
                          box_penalty(rightBox, box)));
    }

    /*
     * Sort "common entries" by calculated deltas in order to distribute
     * the most ambiguous entries first.
     */
    qsort(commonEntries, commonEntriesCount, sizeof(CommonEntry), common_entry_cmp);

    /*
     * Distribute "common entries" between groups.
     */
    for (i = 0; i < commonEntriesCount; i++)
    {
      box = DatumGetBoxP(entryvec->vector[commonEntries[i].index].key);

      /*
       * Check if we have to place this entry in either group to achieve
       * LIMIT_RATIO.
       */
      if (v->spl_nleft + (commonEntriesCount - i) <= m)
        PLACE_LEFT(box, commonEntries[i].index);
      else if (v->spl_nright + (commonEntriesCount - i) <= m)
        PLACE_RIGHT(box, commonEntries[i].index);
      else
      {
        /* Otherwise select the group by minimal penalty */
        if (box_penalty(leftBox, box) < box_penalty(rightBox, box))
          PLACE_LEFT(box, commonEntries[i].index);
        else
          PLACE_RIGHT(box, commonEntries[i].index);
      }
    }
  }

  v->spl_ldatum = PointerGetDatum(leftBox);
  v->spl_rdatum = PointerGetDatum(rightBox);
  PG_RETURN_POINTER(v);
}

/*
 * Equality method
 *
 * This is used for boxes, points, circles, and polygons, all of which store
 * boxes as GiST index entries.
 *
 * Returns true only when boxes are exactly the same.  We can't use fuzzy
 * comparisons here without breaking index consistency; therefore, this isn't
 * equivalent to box_same().
 */
Datum
gist_box_same(PG_FUNCTION_ARGS)
{
  BOX      *b1 = PG_GETARG_BOX_P(0);
  BOX      *b2 = PG_GETARG_BOX_P(1);
  bool     *result = (bool *) PG_GETARG_POINTER(2);

  if (b1 && b2)
    *result = (float8_eq(b1->low.x, b2->low.x) &&
           float8_eq(b1->low.y, b2->low.y) &&
           float8_eq(b1->high.x, b2->high.x) &&
           float8_eq(b1->high.y, b2->high.y));
  else
    *result = (b1 == NULL && b2 == NULL);
  PG_RETURN_POINTER(result);
}

/*
 * Leaf-level consistency for boxes: just apply the query operator
 */
static bool
gist_box_leaf_consistent(BOX *key, BOX *query, StrategyNumber strategy)
{
  bool    retval;

  switch (strategy)
  {
    case RTLeftStrategyNumber:
      retval = DatumGetBool(DirectFunctionCall2(box_left,
                            PointerGetDatum(key),
                            PointerGetDatum(query)));
      break;
    case RTOverLeftStrategyNumber:
      retval = DatumGetBool(DirectFunctionCall2(box_overleft,
                            PointerGetDatum(key),
                            PointerGetDatum(query)));
      break;
    case RTOverlapStrategyNumber:
      retval = DatumGetBool(DirectFunctionCall2(box_overlap,
                            PointerGetDatum(key),
                            PointerGetDatum(query)));
      break;
    case RTOverRightStrategyNumber:
      retval = DatumGetBool(DirectFunctionCall2(box_overright,
                            PointerGetDatum(key),
                            PointerGetDatum(query)));
      break;
    case RTRightStrategyNumber:
      retval = DatumGetBool(DirectFunctionCall2(box_right,
                            PointerGetDatum(key),
                            PointerGetDatum(query)));
      break;
    case RTSameStrategyNumber:
      retval = DatumGetBool(DirectFunctionCall2(box_same,
                            PointerGetDatum(key),
                            PointerGetDatum(query)));
      break;
    case RTContainsStrategyNumber:
      retval = DatumGetBool(DirectFunctionCall2(box_contain,
                            PointerGetDatum(key),
                            PointerGetDatum(query)));
      break;
    case RTContainedByStrategyNumber:
      retval = DatumGetBool(DirectFunctionCall2(box_contained,
                            PointerGetDatum(key),
                            PointerGetDatum(query)));
      break;
    case RTOverBelowStrategyNumber:
      retval = DatumGetBool(DirectFunctionCall2(box_overbelow,
                            PointerGetDatum(key),
                            PointerGetDatum(query)));
      break;
    case RTBelowStrategyNumber:
      retval = DatumGetBool(DirectFunctionCall2(box_below,
                            PointerGetDatum(key),
                            PointerGetDatum(query)));
      break;
    case RTAboveStrategyNumber:
      retval = DatumGetBool(DirectFunctionCall2(box_above,
                            PointerGetDatum(key),
                            PointerGetDatum(query)));
      break;
    case RTOverAboveStrategyNumber:
      retval = DatumGetBool(DirectFunctionCall2(box_overabove,
                            PointerGetDatum(key),
                            PointerGetDatum(query)));
      break;
    default:
      elog(ERROR, "unrecognized strategy number: %d", strategy);
      retval = false;    /* keep compiler quiet */
      break;
  }
  return retval;
}

/*****************************************
 * Common rtree functions (for boxes, polygons, and circles)
 *****************************************/

/*
 * Internal-page consistency for all these types
 *
 * We can use the same function since all types use bounding boxes as the
 * internal-page representation.
 */
static bool
rtree_internal_consistent(BOX *key, BOX *query, StrategyNumber strategy)
{
  bool    retval;

  switch (strategy)
  {
    case RTLeftStrategyNumber:
      retval = !DatumGetBool(DirectFunctionCall2(box_overright,
                             PointerGetDatum(key),
                             PointerGetDatum(query)));
      break;
    case RTOverLeftStrategyNumber:
      retval = !DatumGetBool(DirectFunctionCall2(box_right,
                             PointerGetDatum(key),
                             PointerGetDatum(query)));
      break;
    case RTOverlapStrategyNumber:
      retval = DatumGetBool(DirectFunctionCall2(box_overlap,
                            PointerGetDatum(key),
                            PointerGetDatum(query)));
      break;
    case RTOverRightStrategyNumber:
      retval = !DatumGetBool(DirectFunctionCall2(box_left,
                             PointerGetDatum(key),
                             PointerGetDatum(query)));
      break;
    case RTRightStrategyNumber:
      retval = !DatumGetBool(DirectFunctionCall2(box_overleft,
                             PointerGetDatum(key),
                             PointerGetDatum(query)));
      break;
    case RTSameStrategyNumber:
    case RTContainsStrategyNumber:
      retval = DatumGetBool(DirectFunctionCall2(box_contain,
                            PointerGetDatum(key),
                            PointerGetDatum(query)));
      break;
    case RTContainedByStrategyNumber:
      retval = DatumGetBool(DirectFunctionCall2(box_overlap,
                            PointerGetDatum(key),
                            PointerGetDatum(query)));
      break;
    case RTOverBelowStrategyNumber:
      retval = !DatumGetBool(DirectFunctionCall2(box_above,
                             PointerGetDatum(key),
                             PointerGetDatum(query)));
      break;
    case RTBelowStrategyNumber:
      retval = !DatumGetBool(DirectFunctionCall2(box_overabove,
                             PointerGetDatum(key),
                             PointerGetDatum(query)));
      break;
    case RTAboveStrategyNumber:
      retval = !DatumGetBool(DirectFunctionCall2(box_overbelow,
                             PointerGetDatum(key),
                             PointerGetDatum(query)));
      break;
    case RTOverAboveStrategyNumber:
      retval = !DatumGetBool(DirectFunctionCall2(box_below,
                             PointerGetDatum(key),
                             PointerGetDatum(query)));
      break;
    default:
      elog(ERROR, "unrecognized strategy number: %d", strategy);
      retval = false;    /* keep compiler quiet */
      break;
  }
  return retval;
}

/**************************************************
 * Polygon ops
 **************************************************/

/*
 * GiST compress for polygons: represent a polygon by its bounding box
 */
Datum
gist_poly_compress(PG_FUNCTION_ARGS)
{
  GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  GISTENTRY  *retval;

  if (entry->leafkey)
  {
    POLYGON    *in = DatumGetPolygonP(entry->key);
    BOX      *r;

    r = palloc_object(BOX);
    memcpy(r, &(in->boundbox), sizeof(BOX));

    retval = palloc_object(GISTENTRY);
    gistentryinit(*retval, PointerGetDatum(r),
            entry->rel, entry->page,
            entry->offset, false);
  }
  else
    retval = entry;
  PG_RETURN_POINTER(retval);
}

/*
 * The GiST Consistent method for polygons
 */
Datum
gist_poly_consistent(PG_FUNCTION_ARGS)
{
  GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  POLYGON    *query = PG_GETARG_POLYGON_P(1);
  StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
#ifdef NOT_USED
  Oid     subtype = PG_GETARG_OID(3);
#endif
  bool     *recheck = (bool *) PG_GETARG_POINTER(4);
  bool    result;

  /* All cases served by this function are inexact */
  *recheck = true;

  if (DatumGetBoxP(entry->key) == NULL || query == NULL)
    PG_RETURN_BOOL(false);

  /*
   * Since the operators require recheck anyway, we can just use
   * rtree_internal_consistent even at leaf nodes.  (This works in part
   * because the index entries are bounding boxes not polygons.)
   */
  result = rtree_internal_consistent(DatumGetBoxP(entry->key),
                     &(query->boundbox), strategy);

  /* Avoid memory leak if supplied poly is toasted */
  PG_FREE_IF_COPY(query, 1);

  PG_RETURN_BOOL(result);
}

/**************************************************
 * Circle ops
 **************************************************/

/*
 * GiST compress for circles: represent a circle by its bounding box
 */
Datum
gist_circle_compress(PG_FUNCTION_ARGS)
{
  GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  GISTENTRY  *retval;

  if (entry->leafkey)
  {
    CIRCLE     *in = DatumGetCircleP(entry->key);
    BOX      *r;

    r = palloc_object(BOX);
    r->high.x = float8_pl(in->center.x, in->radius);
    r->low.x = float8_mi(in->center.x, in->radius);
    r->high.y = float8_pl(in->center.y, in->radius);
    r->low.y = float8_mi(in->center.y, in->radius);

    retval = palloc_object(GISTENTRY);
    gistentryinit(*retval, PointerGetDatum(r),
            entry->rel, entry->page,
            entry->offset, false);
  }
  else
    retval = entry;
  PG_RETURN_POINTER(retval);
}

/*
 * The GiST Consistent method for circles
 */
Datum
gist_circle_consistent(PG_FUNCTION_ARGS)
{
  GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  CIRCLE     *query = PG_GETARG_CIRCLE_P(1);
  StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
#ifdef NOT_USED
  Oid     subtype = PG_GETARG_OID(3);
#endif
  bool     *recheck = (bool *) PG_GETARG_POINTER(4);
  BOX     bbox;
  bool    result;

  /* All cases served by this function are inexact */
  *recheck = true;

  if (DatumGetBoxP(entry->key) == NULL || query == NULL)
    PG_RETURN_BOOL(false);

  /*
   * Since the operators require recheck anyway, we can just use
   * rtree_internal_consistent even at leaf nodes.  (This works in part
   * because the index entries are bounding boxes not circles.)
   */
  bbox.high.x = float8_pl(query->center.x, query->radius);
  bbox.low.x = float8_mi(query->center.x, query->radius);
  bbox.high.y = float8_pl(query->center.y, query->radius);
  bbox.low.y = float8_mi(query->center.y, query->radius);

  result = rtree_internal_consistent(DatumGetBoxP(entry->key),
                     &bbox, strategy);

  PG_RETURN_BOOL(result);
}

/**************************************************
 * Point ops
 **************************************************/

Datum
gist_point_compress(PG_FUNCTION_ARGS)
{
  GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);

  if (entry->leafkey)     /* Point, actually */
  {
    BOX      *box = palloc_object(BOX);
    Point    *point = DatumGetPointP(entry->key);
    GISTENTRY  *retval = palloc_object(GISTENTRY);

    box->high = box->low = *point;

    gistentryinit(*retval, BoxPGetDatum(box),
            entry->rel, entry->page, entry->offset, false);

    PG_RETURN_POINTER(retval);
  }

  PG_RETURN_POINTER(entry);
}

/*
 * GiST Fetch method for point
 *
 * Get point coordinates from its bounding box coordinates and form new
 * gistentry.
 */
Datum
gist_point_fetch(PG_FUNCTION_ARGS)
{
  GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  BOX      *in = DatumGetBoxP(entry->key);
  Point    *r;
  GISTENTRY  *retval;

  retval = palloc_object(GISTENTRY);

  r = palloc_object(Point);
  r->x = in->high.x;
  r->y = in->high.y;
  gistentryinit(*retval, PointerGetDatum(r),
          entry->rel, entry->page,
          entry->offset, false);

  PG_RETURN_POINTER(retval);
}


#define point_point_distance(p1,p2) \
  DatumGetFloat8(DirectFunctionCall2(point_distance, \
                     PointPGetDatum(p1), PointPGetDatum(p2)))

static float8
computeDistance(bool isLeaf, BOX *box, Point *point)
{
  float8    result = 0.0;

  if (isLeaf)
  {
    /* simple point to point distance */
    result = point_point_distance(point, &box->low);
  }
  else if (point->x <= box->high.x && point->x >= box->low.x &&
       point->y <= box->high.y && point->y >= box->low.y)
  {
    /* point inside the box */
    result = 0.0;
  }
  else if (point->x <= box->high.x && point->x >= box->low.x)
  {
    /* point is over or below box */
    Assert(box->low.y <= box->high.y);
    if (point->y > box->high.y)
      result = float8_mi(point->y, box->high.y);
    else if (point->y < box->low.y)
      result = float8_mi(box->low.y, point->y);
    else
      elog(ERROR, "inconsistent point values");
  }
  else if (point->y <= box->high.y && point->y >= box->low.y)
  {
    /* point is to left or right of box */
    Assert(box->low.x <= box->high.x);
    if (point->x > box->high.x)
      result = float8_mi(point->x, box->high.x);
    else if (point->x < box->low.x)
      result = float8_mi(box->low.x, point->x);
    else
      elog(ERROR, "inconsistent point values");
  }
  else
  {
    /* closest point will be a vertex */
    Point   p;
    float8    subresult;

    result = point_point_distance(point, &box->low);

    subresult = point_point_distance(point, &box->high);
    if (result > subresult)
      result = subresult;

    p.x = box->low.x;
    p.y = box->high.y;
    subresult = point_point_distance(point, &p);
    if (result > subresult)
      result = subresult;

    p.x = box->high.x;
    p.y = box->low.y;
    subresult = point_point_distance(point, &p);
    if (result > subresult)
      result = subresult;
  }

  return result;
}

static bool
gist_point_consistent_internal(StrategyNumber strategy,
                 bool isLeaf, BOX *key, Point *query)
{
  bool    result = false;

  switch (strategy)
  {
    case RTLeftStrategyNumber:
      result = FPlt(key->low.x, query->x);
      break;
    case RTRightStrategyNumber:
      result = FPgt(key->high.x, query->x);
      break;
    case RTAboveStrategyNumber:
      result = FPgt(key->high.y, query->y);
      break;
    case RTBelowStrategyNumber:
      result = FPlt(key->low.y, query->y);
      break;
    case RTSameStrategyNumber:
      if (isLeaf)
      {
        /* key.high must equal key.low, so we can disregard it */
        result = (FPeq(key->low.x, query->x) &&
              FPeq(key->low.y, query->y));
      }
      else
      {
        result = (FPle(query->x, key->high.x) &&
              FPge(query->x, key->low.x) &&
              FPle(query->y, key->high.y) &&
              FPge(query->y, key->low.y));
      }
      break;
    default:
      elog(ERROR, "unrecognized strategy number: %d", strategy);
      result = false;    /* keep compiler quiet */
      break;
  }

  return result;
}

#define GeoStrategyNumberOffset   20
#define PointStrategyNumberGroup  0
#define BoxStrategyNumberGroup    1
#define PolygonStrategyNumberGroup  2
#define CircleStrategyNumberGroup 3

Datum
gist_point_consistent(PG_FUNCTION_ARGS)
{
  GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
  bool     *recheck = (bool *) PG_GETARG_POINTER(4);
  bool    result;
  StrategyNumber strategyGroup;

  /*
   * We have to remap these strategy numbers to get this klugy
   * classification logic to work.
   */
  if (strategy == RTOldBelowStrategyNumber)
    strategy = RTBelowStrategyNumber;
  else if (strategy == RTOldAboveStrategyNumber)
    strategy = RTAboveStrategyNumber;

  strategyGroup = strategy / GeoStrategyNumberOffset;
  switch (strategyGroup)
  {
    case PointStrategyNumberGroup:
      result = gist_point_consistent_internal(strategy % GeoStrategyNumberOffset,
                          GIST_LEAF(entry),
                          DatumGetBoxP(entry->key),
                          PG_GETARG_POINT_P(1));
      *recheck = false;
      break;
    case BoxStrategyNumberGroup:
      {
        /*
         * The only operator in this group is point <@ box (on_pb), so
         * we needn't examine strategy again.
         *
         * For historical reasons, on_pb uses exact rather than fuzzy
         * comparisons.  We could use box_overlap when at an internal
         * page, but that would lead to possibly visiting child pages
         * uselessly, because box_overlap uses fuzzy comparisons.
         * Instead we write a non-fuzzy overlap test.  The same code
         * will also serve for leaf-page tests, since leaf keys have
         * high == low.
         */
        BOX      *query,
               *key;

        query = PG_GETARG_BOX_P(1);
        key = DatumGetBoxP(entry->key);

        result = (key->high.x >= query->low.x &&
              key->low.x <= query->high.x &&
              key->high.y >= query->low.y &&
              key->low.y <= query->high.y);
        *recheck = false;
      }
      break;
    case PolygonStrategyNumberGroup:
      {
        POLYGON    *query = PG_GETARG_POLYGON_P(1);

        result = DatumGetBool(DirectFunctionCall5(gist_poly_consistent,
                              PointerGetDatum(entry),
                              PolygonPGetDatum(query),
                              Int16GetDatum(RTOverlapStrategyNumber),
                              0, PointerGetDatum(recheck)));

        if (GIST_LEAF(entry) && result)
        {
          /*
           * We are on leaf page and quick check shows overlapping
           * of polygon's bounding box and point
           */
          BOX      *box = DatumGetBoxP(entry->key);

          Assert(box->high.x == box->low.x
               && box->high.y == box->low.y);
          result = DatumGetBool(DirectFunctionCall2(poly_contain_pt,
                                PolygonPGetDatum(query),
                                PointPGetDatum(&box->high)));
          *recheck = false;
        }
      }
      break;
    case CircleStrategyNumberGroup:
      {
        CIRCLE     *query = PG_GETARG_CIRCLE_P(1);

        result = DatumGetBool(DirectFunctionCall5(gist_circle_consistent,
                              PointerGetDatum(entry),
                              CirclePGetDatum(query),
                              Int16GetDatum(RTOverlapStrategyNumber),
                              0, PointerGetDatum(recheck)));

        if (GIST_LEAF(entry) && result)
        {
          /*
           * We are on leaf page and quick check shows overlapping
           * of polygon's bounding box and point
           */
          BOX      *box = DatumGetBoxP(entry->key);

          Assert(box->high.x == box->low.x
               && box->high.y == box->low.y);
          result = DatumGetBool(DirectFunctionCall2(circle_contain_pt,
                                CirclePGetDatum(query),
                                PointPGetDatum(&box->high)));
          *recheck = false;
        }
      }
      break;
    default:
      elog(ERROR, "unrecognized strategy number: %d", strategy);
      result = false;    /* keep compiler quiet */
      break;
  }

  PG_RETURN_BOOL(result);
}

Datum
gist_point_distance(PG_FUNCTION_ARGS)
{
  GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
  float8    distance;
  StrategyNumber strategyGroup = strategy / GeoStrategyNumberOffset;

  switch (strategyGroup)
  {
    case PointStrategyNumberGroup:
      distance = computeDistance(GIST_LEAF(entry),
                     DatumGetBoxP(entry->key),
                     PG_GETARG_POINT_P(1));
      break;
    default:
      elog(ERROR, "unrecognized strategy number: %d", strategy);
      distance = 0.0;   /* keep compiler quiet */
      break;
  }

  PG_RETURN_FLOAT8(distance);
}

static float8
gist_bbox_distance(GISTENTRY *entry, Datum query, StrategyNumber strategy)
{
  float8    distance;
  StrategyNumber strategyGroup = strategy / GeoStrategyNumberOffset;

  switch (strategyGroup)
  {
    case PointStrategyNumberGroup:
      distance = computeDistance(false,
                     DatumGetBoxP(entry->key),
                     DatumGetPointP(query));
      break;
    default:
      elog(ERROR, "unrecognized strategy number: %d", strategy);
      distance = 0.0;   /* keep compiler quiet */
  }

  return distance;
}

Datum
gist_box_distance(PG_FUNCTION_ARGS)
{
  GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  Datum   query = PG_GETARG_DATUM(1);
  StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
#ifdef NOT_USED
  Oid     subtype = PG_GETARG_OID(3);
  bool     *recheck = (bool *) PG_GETARG_POINTER(4);
#endif
  float8    distance;

  distance = gist_bbox_distance(entry, query, strategy);

  PG_RETURN_FLOAT8(distance);
}

/*
 * The inexact GiST distance methods for geometric types that store bounding
 * boxes.
 *
 * Compute lossy distance from point to index entries.  The result is inexact
 * because index entries are bounding boxes, not the exact shapes of the
 * indexed geometric types.  We use distance from point to MBR of index entry.
 * This is a lower bound estimate of distance from point to indexed geometric
 * type.
 */
Datum
gist_circle_distance(PG_FUNCTION_ARGS)
{
  GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  Datum   query = PG_GETARG_DATUM(1);
  StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
#ifdef NOT_USED
  Oid     subtype = PG_GETARG_OID(3);
#endif
  bool     *recheck = (bool *) PG_GETARG_POINTER(4);
  float8    distance;

  distance = gist_bbox_distance(entry, query, strategy);
  *recheck = true;

  PG_RETURN_FLOAT8(distance);
}

Datum
gist_poly_distance(PG_FUNCTION_ARGS)
{
  GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  Datum   query = PG_GETARG_DATUM(1);
  StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
#ifdef NOT_USED
  Oid     subtype = PG_GETARG_OID(3);
#endif
  bool     *recheck = (bool *) PG_GETARG_POINTER(4);
  float8    distance;

  distance = gist_bbox_distance(entry, query, strategy);
  *recheck = true;

  PG_RETURN_FLOAT8(distance);
}

/*
 * Z-order routines for fast index build
 */

/*
 * Compute Z-value of a point
 *
 * Z-order (also known as Morton Code) maps a two-dimensional point to a
 * single integer, in a way that preserves locality. Points that are close in
 * the two-dimensional 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.
 *
 * 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.
 */
static uint64
point_zorder_internal(float4 x, float4 y)
{
  uint32    ix = ieee_float32_to_uint32(x);
  uint32    iy = ieee_float32_to_uint32(y);

  /* Interleave the bits */
  return part_bits32_by2(ix) | (part_bits32_by2(iy) << 1);
}

/* Interleave 32 bits with zeroes */
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);

  return n;
}

/*
 * Convert a 32-bit IEEE float to uint32 in a way that preserves the ordering
 */
static uint32
ieee_float32_to_uint32(float f)
{
  /*----
   *
   * 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 all
   * zeros. 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 0-7FFFFFFF
   * x = 0 to value 8000000 (both positive and negative zero)
   * x > 0 to range 8000001-FFFFFFFF
   *
   * 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 any non-NaN values that would be
   * mapped to FFFFFFFF. In fact, there is a range of unused values on both
   * ends of the uint32 space.
   */
  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 0-7FFFFFFF. This flips the sign
       * bit to 0 in the same instruction.
       */
      Assert(f <= 0);   /* can be -0 */
      u.i ^= 0xFFFFFFFF;
    }
    else
    {
      /* Map the positive value (or 0) to range 80000000-FFFFFFFF */
      u.i |= 0x80000000;
    }

    return u.i;
  }
}

/*
 * Compare the Z-order of points
 */
static int
gist_bbox_zorder_cmp(Datum a, Datum b, SortSupport ssup)
{
  Point    *p1 = &(DatumGetBoxP(a)->low);
  Point    *p2 = &(DatumGetBoxP(b)->low);
  uint64    z1;
  uint64    z2;

  /*
   * Do a quick check for equality first. It's not clear if this is worth it
   * in general, but certainly is when used as tie-breaker with abbreviated
   * keys,
   */
  if (p1->x == p2->x && p1->y == p2->y)
    return 0;

  z1 = point_zorder_internal(p1->x, p1->y);
  z2 = point_zorder_internal(p2->x, p2->y);
  if (z1 > z2)
    return 1;
  else if (z1 < z2)
    return -1;
  else
    return 0;
}

/*
 * Abbreviated version of Z-order comparison
 *
 * The abbreviated format is a Z-order value computed from the two 32-bit
 * floats.  Now that sizeof(Datum) is always 8, the 64-bit Z-order value
 * always fits fully in the abbreviated Datum.
 */
static Datum
gist_bbox_zorder_abbrev_convert(Datum original, SortSupport ssup)
{
  Point    *p = &(DatumGetBoxP(original)->low);
  uint64    z;

  z = point_zorder_internal(p->x, p->y);

  return UInt64GetDatum(z);
}

/*
 * We never consider aborting the abbreviation.
 *
 * On 64-bit systems, the abbreviation is not lossy so it is always
 * worthwhile. (Perhaps it's not on 32-bit systems, but we don't bother
 * with logic to decide.)
 */
static bool
gist_bbox_zorder_abbrev_abort(int memtupcount, SortSupport ssup)
{
  return false;
}

/*
 * Sort support routine for fast GiST index build by sorting.
 */
Datum
gist_point_sortsupport(PG_FUNCTION_ARGS)
{
  SortSupport ssup = (SortSupport) PG_GETARG_POINTER(0);

  if (ssup->abbreviate)
  {
    ssup->comparator = ssup_datum_uint64_cmp;
    ssup->abbrev_converter = gist_bbox_zorder_abbrev_convert;
    ssup->abbrev_abort = gist_bbox_zorder_abbrev_abort;
    ssup->abbrev_full_comparator = gist_bbox_zorder_cmp;
  }
  else
  {
    ssup->comparator = gist_bbox_zorder_cmp;
  }
  PG_RETURN_VOID();
}

Coverage Report

Created: 2026-10-16 06:48