1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/latlon-v0008.c	Tue Oct 25 18:44:43 2016 +0200
     1.3 @@ -0,0 +1,3285 @@
     1.4 +
     1.5 +/*-------------*
     1.6 + *  C prelude  *
     1.7 + *-------------*/
     1.8 +
     1.9 +#include "postgres.h"
    1.10 +#include "fmgr.h"
    1.11 +#include "libpq/pqformat.h"
    1.12 +#include "access/gist.h"
    1.13 +#include "access/stratnum.h"
    1.14 +#include "utils/array.h"
    1.15 +#include <math.h>
    1.16 +
    1.17 +#ifdef PG_MODULE_MAGIC
    1.18 +PG_MODULE_MAGIC;
    1.19 +#endif
    1.20 +
    1.21 +#if INT_MAX < 2147483647
    1.22 +#error Expected int type to be at least 32 bit wide
    1.23 +#endif
    1.24 +
    1.25 +
    1.26 +/*---------------------------------*
    1.27 + *  distance calculation on earth  *
    1.28 + *  (using WGS-84 spheroid)        *
    1.29 + *---------------------------------*/
    1.30 +
    1.31 +/*  WGS-84 spheroid with following parameters:
    1.32 +    semi-major axis  a = 6378137
    1.33 +    semi-minor axis  b = a * (1 - 1/298.257223563)
    1.34 +    estimated diameter = 2 * (2*a+b)/3
    1.35 +*/
    1.36 +#define PGL_SPHEROID_A 6378137.0            /* semi major axis */
    1.37 +#define PGL_SPHEROID_F (1.0/298.257223563)  /* flattening */
    1.38 +#define PGL_SPHEROID_B (PGL_SPHEROID_A * (1.0-PGL_SPHEROID_F))
    1.39 +#define PGL_EPS2       ( ( PGL_SPHEROID_A * PGL_SPHEROID_A - \
    1.40 +                           PGL_SPHEROID_B * PGL_SPHEROID_B ) / \
    1.41 +                         ( PGL_SPHEROID_A * PGL_SPHEROID_A ) )
    1.42 +#define PGL_SUBEPS2    (1.0-PGL_EPS2)
    1.43 +#define PGL_RADIUS     ((2.0*PGL_SPHEROID_A + PGL_SPHEROID_B) / 3.0)
    1.44 +#define PGL_DIAMETER   (2.0 * PGL_RADIUS)
    1.45 +#define PGL_SCALE      (PGL_SPHEROID_A / PGL_DIAMETER)  /* semi-major ref. */
    1.46 +#define PGL_MAXDIST    (PGL_RADIUS * M_PI)              /* maximum distance */
    1.47 +#define PGL_FADELIMIT  (PGL_MAXDIST / 3.0)              /* 1/6 circumference */
    1.48 +
    1.49 +/* calculate distance between two points on earth (given in degrees) */
    1.50 +static inline double pgl_distance(
    1.51 +  double lat1, double lon1, double lat2, double lon2
    1.52 +) {
    1.53 +  float8 lat1cos, lat1sin, lat2cos, lat2sin, lon2cos, lon2sin;
    1.54 +  float8 nphi1, nphi2, x1, z1, x2, y2, z2, g, s, t;
    1.55 +  /* normalize delta longitude (lon2 > 0 && lon1 = 0) */
    1.56 +  /* lon1 = 0 (not used anymore) */
    1.57 +  lon2 = fabs(lon2-lon1);
    1.58 +  /* convert to radians (first divide, then multiply) */
    1.59 +  lat1 = (lat1 / 180.0) * M_PI;
    1.60 +  lat2 = (lat2 / 180.0) * M_PI;
    1.61 +  lon2 = (lon2 / 180.0) * M_PI;
    1.62 +  /* make lat2 >= lat1 to ensure reversal-symmetry despite floating point
    1.63 +     operations (lon2 >= lon1 is already ensured in a previous step) */
    1.64 +  if (lat2 < lat1) { float8 swap = lat1; lat1 = lat2; lat2 = swap; }
    1.65 +  /* calculate 3d coordinates on scaled ellipsoid which has an average diameter
    1.66 +     of 1.0 */
    1.67 +  lat1cos = cos(lat1); lat1sin = sin(lat1);
    1.68 +  lat2cos = cos(lat2); lat2sin = sin(lat2);
    1.69 +  lon2cos = cos(lon2); lon2sin = sin(lon2);
    1.70 +  nphi1 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat1sin * lat1sin);
    1.71 +  nphi2 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat2sin * lat2sin);
    1.72 +  x1 = nphi1 * lat1cos;
    1.73 +  z1 = nphi1 * PGL_SUBEPS2 * lat1sin;
    1.74 +  x2 = nphi2 * lat2cos * lon2cos;
    1.75 +  y2 = nphi2 * lat2cos * lon2sin;
    1.76 +  z2 = nphi2 * PGL_SUBEPS2 * lat2sin;
    1.77 +  /* calculate tunnel distance through scaled (diameter 1.0) ellipsoid */
    1.78 +  g = sqrt((x2-x1)*(x2-x1) + y2*y2 + (z2-z1)*(z2-z1));
    1.79 +  /* convert tunnel distance through scaled ellipsoid to approximated surface
    1.80 +     distance on original ellipsoid */
    1.81 +  if (g > 1.0) g = 1.0;
    1.82 +  s = PGL_DIAMETER * asin(g);
    1.83 +  /* return result only if small enough to be precise (less than 1/3 of
    1.84 +     maximum possible distance) */
    1.85 +  if (s <= PGL_FADELIMIT) return s;
    1.86 +  /* calculate tunnel distance to antipodal point through scaled ellipsoid */
    1.87 +  g = sqrt((x2+x1)*(x2+x1) + y2*y2 + (z2+z1)*(z2+z1));
    1.88 +  /* convert tunnel distance to antipodal point through scaled ellipsoid to
    1.89 +     approximated surface distance to antipodal point on original ellipsoid */
    1.90 +  if (g > 1.0) g = 1.0;
    1.91 +  t = PGL_DIAMETER * asin(g);
    1.92 +  /* surface distance between original points can now be approximated by
    1.93 +     substracting antipodal distance from maximum possible distance;
    1.94 +     return result only if small enough (less than 1/3 of maximum possible
    1.95 +     distance) */
    1.96 +  if (t <= PGL_FADELIMIT) return PGL_MAXDIST-t;
    1.97 +  /* otherwise crossfade direct and antipodal result to ensure monotonicity */
    1.98 +  return (
    1.99 +    (s * (t-PGL_FADELIMIT) + (PGL_MAXDIST-t) * (s-PGL_FADELIMIT)) /
   1.100 +    (s + t - 2*PGL_FADELIMIT)
   1.101 +  );
   1.102 +}
   1.103 +
   1.104 +/* finite distance that can not be reached on earth */
   1.105 +#define PGL_ULTRA_DISTANCE (3 * PGL_MAXDIST)
   1.106 +
   1.107 +
   1.108 +/*--------------------------------*
   1.109 + *  simple geographic data types  *
   1.110 + *--------------------------------*/
   1.111 +
   1.112 +/* point on earth given by latitude and longitude in degrees */
   1.113 +/* (type "epoint" in SQL) */
   1.114 +typedef struct {
   1.115 +  double lat;  /* between  -90 and  90 (both inclusive) */
   1.116 +  double lon;  /* between -180 and 180 (both inclusive) */
   1.117 +} pgl_point;
   1.118 +
   1.119 +/* box delimited by two parallels and two meridians (all in degrees) */
   1.120 +/* (type "ebox" in SQL) */
   1.121 +typedef struct {
   1.122 +  double lat_min;  /* between  -90 and  90 (both inclusive) */
   1.123 +  double lat_max;  /* between  -90 and  90 (both inclusive) */
   1.124 +  double lon_min;  /* between -180 and 180 (both inclusive) */
   1.125 +  double lon_max;  /* between -180 and 180 (both inclusive) */
   1.126 +  /* if lat_min > lat_max, then box is empty */
   1.127 +  /* if lon_min > lon_max, then 180th meridian is crossed */
   1.128 +} pgl_box;
   1.129 +
   1.130 +/* circle on earth surface (for radial searches with fixed radius) */
   1.131 +/* (type "ecircle" in SQL) */
   1.132 +typedef struct {
   1.133 +  pgl_point center;
   1.134 +  double radius; /* positive (including +0 but excluding -0), or -INFINITY */
   1.135 +  /* A negative radius (i.e. -INFINITY) denotes nothing (i.e. no point),
   1.136 +     zero radius (0) denotes a single point,
   1.137 +     a finite radius (0 < radius < INFINITY) denotes a filled circle, and
   1.138 +     a radius of INFINITY is valid and means complete coverage of earth. */
   1.139 +} pgl_circle;
   1.140 +
   1.141 +
   1.142 +/*----------------------------------*
   1.143 + *  geographic "cluster" data type  *
   1.144 + *----------------------------------*/
   1.145 +
   1.146 +/* A cluster is a collection of points, paths, outlines, and polygons. If two
   1.147 +   polygons in a cluster overlap, the area covered by both polygons does not
   1.148 +   belong to the cluster. This way, a cluster can be used to describe complex
   1.149 +   shapes like polygons with holes. Outlines are non-filled polygons. Paths are
   1.150 +   open by default (i.e. the last point in the list is not connected with the
   1.151 +   first point in the list). Note that each outline or polygon in a cluster
   1.152 +   must cover a longitude range of less than 180 degrees to avoid ambiguities.
   1.153 +   Areas which are larger may be split into multiple polygons. */
   1.154 +
   1.155 +/* maximum number of points in a cluster */
   1.156 +/* (limited to avoid integer overflows, e.g. when allocating memory) */
   1.157 +#define PGL_CLUSTER_MAXPOINTS 16777216
   1.158 +
   1.159 +/* types of cluster entries */
   1.160 +#define PGL_ENTRY_POINT   1  /* a point */
   1.161 +#define PGL_ENTRY_PATH    2  /* a path from first point to last point */
   1.162 +#define PGL_ENTRY_OUTLINE 3  /* a non-filled polygon with given vertices */
   1.163 +#define PGL_ENTRY_POLYGON 4  /* a filled polygon with given vertices */
   1.164 +
   1.165 +/* Entries of a cluster are described by two different structs: pgl_newentry
   1.166 +   and pgl_entry. The first is used only during construction of a cluster, the
   1.167 +   second is used in all other cases (e.g. when reading clusters from the
   1.168 +   database, performing operations, etc). */
   1.169 +
   1.170 +/* entry for new geographic cluster during construction of that cluster */
   1.171 +typedef struct {
   1.172 +  int32_t entrytype;
   1.173 +  int32_t npoints;
   1.174 +  pgl_point *points;  /* pointer to an array of points (pgl_point) */
   1.175 +} pgl_newentry;
   1.176 +
   1.177 +/* entry of geographic cluster */
   1.178 +typedef struct {
   1.179 +  int32_t entrytype;  /* type of entry: point, path, outline, polygon */
   1.180 +  int32_t npoints;    /* number of stored points (set to 1 for point entry) */
   1.181 +  int32_t offset;     /* offset of pgl_point array from cluster base address */
   1.182 +  /* use macro PGL_ENTRY_POINTS to obtain a pointer to the array of points */
   1.183 +} pgl_entry;
   1.184 +
   1.185 +/* geographic cluster which is a collection of points, (open) paths, polygons,
   1.186 +   and outlines (non-filled polygons) */
   1.187 +typedef struct {
   1.188 +  char header[VARHDRSZ];  /* PostgreSQL header for variable size data types */
   1.189 +  int32_t nentries;       /* number of stored points */
   1.190 +  pgl_circle bounding;    /* bounding circle */
   1.191 +  /* Note: bounding circle ensures alignment of pgl_cluster for points */
   1.192 +  pgl_entry entries[FLEXIBLE_ARRAY_MEMBER];  /* var-length data */
   1.193 +} pgl_cluster;
   1.194 +
   1.195 +/* macro to determine memory alignment of points */
   1.196 +/* (needed to store pgl_point array after entries in pgl_cluster) */
   1.197 +typedef struct { char dummy; pgl_point aligned; } pgl_point_alignment;
   1.198 +#define PGL_POINT_ALIGNMENT offsetof(pgl_point_alignment, aligned)
   1.199 +
   1.200 +/* macro to extract a pointer to the array of points of a cluster entry */
   1.201 +#define PGL_ENTRY_POINTS(cluster, idx) \
   1.202 +  ((pgl_point *)(((intptr_t)cluster)+(cluster)->entries[idx].offset))
   1.203 +
   1.204 +/* convert pgl_newentry array to pgl_cluster */
   1.205 +/* NOTE: requires pgl_finalize_cluster to be called to finalize result */
   1.206 +static pgl_cluster *pgl_new_cluster(int nentries, pgl_newentry *entries) {
   1.207 +  int i;              /* index of current entry */
   1.208 +  int npoints = 0;    /* number of points in whole cluster */
   1.209 +  int entry_npoints;  /* number of points in current entry */
   1.210 +  int points_offset = PGL_POINT_ALIGNMENT * (
   1.211 +    ( offsetof(pgl_cluster, entries) +
   1.212 +      nentries * sizeof(pgl_entry) +
   1.213 +      PGL_POINT_ALIGNMENT - 1
   1.214 +    ) / PGL_POINT_ALIGNMENT
   1.215 +  );  /* offset of pgl_point array from base address (considering alignment) */
   1.216 +  pgl_cluster *cluster;  /* new cluster to be returned */
   1.217 +  /* determine total number of points */
   1.218 +  for (i=0; i<nentries; i++) npoints += entries[i].npoints;
   1.219 +  /* allocate memory for cluster (including entries and points) */
   1.220 +  cluster = palloc(points_offset + npoints * sizeof(pgl_point));
   1.221 +  /* re-count total number of points to determine offset for each entry */
   1.222 +  npoints = 0;
   1.223 +  /* copy entries and points */
   1.224 +  for (i=0; i<nentries; i++) {
   1.225 +    /* determine number of points in entry */
   1.226 +    entry_npoints = entries[i].npoints;
   1.227 +    /* copy entry */
   1.228 +    cluster->entries[i].entrytype = entries[i].entrytype;
   1.229 +    cluster->entries[i].npoints = entry_npoints;
   1.230 +    /* calculate offset (in bytes) of pgl_point array */
   1.231 +    cluster->entries[i].offset = points_offset + npoints * sizeof(pgl_point);
   1.232 +    /* copy points */
   1.233 +    memcpy(
   1.234 +      PGL_ENTRY_POINTS(cluster, i),
   1.235 +      entries[i].points,
   1.236 +      entry_npoints * sizeof(pgl_point)
   1.237 +    );
   1.238 +    /* update total number of points processed */
   1.239 +    npoints += entry_npoints;
   1.240 +  }
   1.241 +  /* set number of entries in cluster */
   1.242 +  cluster->nentries = nentries;
   1.243 +  /* set PostgreSQL header for variable sized data */
   1.244 +  SET_VARSIZE(cluster, points_offset + npoints * sizeof(pgl_point));
   1.245 +  /* return newly created cluster */
   1.246 +  return cluster;
   1.247 +}
   1.248 +
   1.249 +
   1.250 +/*----------------------------------------*
   1.251 + *  C functions on geographic data types  *
   1.252 + *----------------------------------------*/
   1.253 +
   1.254 +/* round latitude or longitude to 12 digits after decimal point */
   1.255 +static inline double pgl_round(double val) {
   1.256 +  return round(val * 1e12) / 1e12;
   1.257 +}
   1.258 +
   1.259 +/* compare two points */
   1.260 +/* (equality when same point on earth is described, otherwise an arbitrary
   1.261 +   linear order) */
   1.262 +static int pgl_point_cmp(pgl_point *point1, pgl_point *point2) {
   1.263 +  double lon1, lon2;  /* modified longitudes for special cases */
   1.264 +  /* use latitude as first ordering criterion */
   1.265 +  if (point1->lat < point2->lat) return -1;
   1.266 +  if (point1->lat > point2->lat) return 1;
   1.267 +  /* determine modified longitudes (considering special case of poles and
   1.268 +     180th meridian which can be described as W180 or E180) */
   1.269 +  if (point1->lat == -90 || point1->lat == 90) lon1 = 0;
   1.270 +  else if (point1->lon == 180) lon1 = -180;
   1.271 +  else lon1 = point1->lon;
   1.272 +  if (point2->lat == -90 || point2->lat == 90) lon2 = 0;
   1.273 +  else if (point2->lon == 180) lon2 = -180;
   1.274 +  else lon2 = point2->lon;
   1.275 +  /* use (modified) longitude as secondary ordering criterion */
   1.276 +  if (lon1 < lon2) return -1;
   1.277 +  if (lon1 > lon2) return 1;
   1.278 +  /* no difference found, points are equal */
   1.279 +  return 0;
   1.280 +}
   1.281 +
   1.282 +/* compare two boxes */
   1.283 +/* (equality when same box on earth is described, otherwise an arbitrary linear
   1.284 +   order) */
   1.285 +static int pgl_box_cmp(pgl_box *box1, pgl_box *box2) {
   1.286 +  /* two empty boxes are equal, and an empty box is always considered "less
   1.287 +     than" a non-empty box */
   1.288 +  if (box1->lat_min> box1->lat_max && box2->lat_min<=box2->lat_max) return -1;
   1.289 +  if (box1->lat_min> box1->lat_max && box2->lat_min> box2->lat_max) return 0;
   1.290 +  if (box1->lat_min<=box1->lat_max && box2->lat_min> box2->lat_max) return 1;
   1.291 +  /* use southern border as first ordering criterion */
   1.292 +  if (box1->lat_min < box2->lat_min) return -1;
   1.293 +  if (box1->lat_min > box2->lat_min) return 1;
   1.294 +  /* use northern border as second ordering criterion */
   1.295 +  if (box1->lat_max < box2->lat_max) return -1;
   1.296 +  if (box1->lat_max > box2->lat_max) return 1;
   1.297 +  /* use western border as third ordering criterion */
   1.298 +  if (box1->lon_min < box2->lon_min) return -1;
   1.299 +  if (box1->lon_min > box2->lon_min) return 1;
   1.300 +  /* use eastern border as fourth ordering criterion */
   1.301 +  if (box1->lon_max < box2->lon_max) return -1;
   1.302 +  if (box1->lon_max > box2->lon_max) return 1;
   1.303 +  /* no difference found, boxes are equal */
   1.304 +  return 0;
   1.305 +}
   1.306 +
   1.307 +/* compare two circles */
   1.308 +/* (equality when same circle on earth is described, otherwise an arbitrary
   1.309 +   linear order) */
   1.310 +static int pgl_circle_cmp(pgl_circle *circle1, pgl_circle *circle2) {
   1.311 +  /* two circles with same infinite radius (positive or negative infinity) are
   1.312 +     considered equal independently of center point */
   1.313 +  if (
   1.314 +    !isfinite(circle1->radius) && !isfinite(circle2->radius) &&
   1.315 +    circle1->radius == circle2->radius
   1.316 +  ) return 0;
   1.317 +  /* use radius as first ordering criterion */
   1.318 +  if (circle1->radius < circle2->radius) return -1;
   1.319 +  if (circle1->radius > circle2->radius) return 1;
   1.320 +  /* use center point as secondary ordering criterion */
   1.321 +  return pgl_point_cmp(&(circle1->center), &(circle2->center));
   1.322 +}
   1.323 +
   1.324 +/* set box to empty box*/
   1.325 +static void pgl_box_set_empty(pgl_box *box) {
   1.326 +  box->lat_min = INFINITY;
   1.327 +  box->lat_max = -INFINITY;
   1.328 +  box->lon_min = 0;
   1.329 +  box->lon_max = 0;
   1.330 +}
   1.331 +
   1.332 +/* check if point is inside a box */
   1.333 +static bool pgl_point_in_box(pgl_point *point, pgl_box *box) {
   1.334 +  return (
   1.335 +    point->lat >= box->lat_min && point->lat <= box->lat_max && (
   1.336 +      (box->lon_min > box->lon_max) ? (
   1.337 +        /* box crosses 180th meridian */
   1.338 +        point->lon >= box->lon_min || point->lon <= box->lon_max
   1.339 +      ) : (
   1.340 +        /* box does not cross the 180th meridian */
   1.341 +        point->lon >= box->lon_min && point->lon <= box->lon_max
   1.342 +      )
   1.343 +    )
   1.344 +  );
   1.345 +}
   1.346 +
   1.347 +/* check if two boxes overlap */
   1.348 +static bool pgl_boxes_overlap(pgl_box *box1, pgl_box *box2) {
   1.349 +  return (
   1.350 +    box2->lat_max >= box2->lat_min &&  /* ensure box2 is not empty */
   1.351 +    ( box2->lat_min >= box1->lat_min || box2->lat_max >= box1->lat_min ) &&
   1.352 +    ( box2->lat_min <= box1->lat_max || box2->lat_max <= box1->lat_max ) && (
   1.353 +      (
   1.354 +        /* check if one and only one box crosses the 180th meridian */
   1.355 +        ((box1->lon_min > box1->lon_max) ? 1 : 0) ^
   1.356 +        ((box2->lon_min > box2->lon_max) ? 1 : 0)
   1.357 +      ) ? (
   1.358 +        /* exactly one box crosses the 180th meridian */
   1.359 +        box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min ||
   1.360 +        box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max
   1.361 +      ) : (
   1.362 +        /* no box or both boxes cross the 180th meridian */
   1.363 +        (
   1.364 +          (box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min) &&
   1.365 +          (box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max)
   1.366 +        ) ||
   1.367 +        /* handle W180 == E180 */
   1.368 +        ( box1->lon_min == -180 && box2->lon_max == 180 ) ||
   1.369 +        ( box2->lon_min == -180 && box1->lon_max == 180 )
   1.370 +      )
   1.371 +    )
   1.372 +  );
   1.373 +}
   1.374 +
   1.375 +/* check unambiguousness of east/west orientation of cluster entries and set
   1.376 +   bounding circle of cluster */
   1.377 +static bool pgl_finalize_cluster(pgl_cluster *cluster) {
   1.378 +  int i, j;                 /* i: index of entry, j: index of point in entry */
   1.379 +  int npoints;              /* number of points in entry */
   1.380 +  int total_npoints = 0;    /* total number of points in cluster */
   1.381 +  pgl_point *points;        /* points in entry */
   1.382 +  int lon_dir;              /* first point of entry west (-1) or east (+1) */
   1.383 +  double lon_break = 0;     /* antipodal longitude of first point in entry */
   1.384 +  double lon_min, lon_max;  /* covered longitude range of entry */
   1.385 +  double value;             /* temporary variable */
   1.386 +  /* reset bounding circle center to empty circle at 0/0 coordinates */
   1.387 +  cluster->bounding.center.lat = 0;
   1.388 +  cluster->bounding.center.lon = 0;
   1.389 +  cluster->bounding.radius = -INFINITY;
   1.390 +  /* if cluster is not empty */
   1.391 +  if (cluster->nentries != 0) {
   1.392 +    /* iterate over all cluster entries and ensure they each cover a longitude
   1.393 +       range less than 180 degrees */
   1.394 +    for (i=0; i<cluster->nentries; i++) {
   1.395 +      /* get properties of entry */
   1.396 +      npoints = cluster->entries[i].npoints;
   1.397 +      points = PGL_ENTRY_POINTS(cluster, i);
   1.398 +      /* get longitude of first point of entry */
   1.399 +      value = points[0].lon;
   1.400 +      /* initialize lon_min and lon_max with longitude of first point */
   1.401 +      lon_min = value;
   1.402 +      lon_max = value;
   1.403 +      /* determine east/west orientation of first point and calculate antipodal
   1.404 +         longitude (Note: rounding required here) */
   1.405 +      if      (value < 0) { lon_dir = -1; lon_break = pgl_round(value + 180); }
   1.406 +      else if (value > 0) { lon_dir =  1; lon_break = pgl_round(value - 180); }
   1.407 +      else lon_dir = 0;
   1.408 +      /* iterate over all other points in entry */
   1.409 +      for (j=1; j<npoints; j++) {
   1.410 +        /* consider longitude wrap-around */
   1.411 +        value = points[j].lon;
   1.412 +        if      (lon_dir<0 && value>lon_break) value = pgl_round(value - 360);
   1.413 +        else if (lon_dir>0 && value<lon_break) value = pgl_round(value + 360);
   1.414 +        /* update lon_min and lon_max */
   1.415 +        if      (value < lon_min) lon_min = value;
   1.416 +        else if (value > lon_max) lon_max = value;
   1.417 +        /* return false if 180 degrees or more are covered */
   1.418 +        if (lon_max - lon_min >= 180) return false;
   1.419 +      }
   1.420 +    }
   1.421 +    /* iterate over all points of all entries and calculate arbitrary center
   1.422 +       point for bounding circle (best if center point minimizes the radius,
   1.423 +       but some error is allowed here) */
   1.424 +    for (i=0; i<cluster->nentries; i++) {
   1.425 +      /* get properties of entry */
   1.426 +      npoints = cluster->entries[i].npoints;
   1.427 +      points = PGL_ENTRY_POINTS(cluster, i);
   1.428 +      /* check if first entry */
   1.429 +      if (i==0) {
   1.430 +        /* get longitude of first point of first entry in whole cluster */
   1.431 +        value = points[0].lon;
   1.432 +        /* initialize lon_min and lon_max with longitude of first point of
   1.433 +           first entry in whole cluster (used to determine if whole cluster
   1.434 +           covers a longitude range of 180 degrees or more) */
   1.435 +        lon_min = value;
   1.436 +        lon_max = value;
   1.437 +        /* determine east/west orientation of first point and calculate
   1.438 +           antipodal longitude (Note: rounding not necessary here) */
   1.439 +        if      (value < 0) { lon_dir = -1; lon_break = value + 180; }
   1.440 +        else if (value > 0) { lon_dir =  1; lon_break = value - 180; }
   1.441 +        else lon_dir = 0;
   1.442 +      }
   1.443 +      /* iterate over all points in entry */
   1.444 +      for (j=0; j<npoints; j++) {
   1.445 +        /* longitude wrap-around (Note: rounding not necessary here) */
   1.446 +        value = points[j].lon;
   1.447 +        if      (lon_dir < 0 && value > lon_break) value -= 360;
   1.448 +        else if (lon_dir > 0 && value < lon_break) value += 360;
   1.449 +        if      (value < lon_min) lon_min = value;
   1.450 +        else if (value > lon_max) lon_max = value;
   1.451 +        /* set bounding circle to cover whole earth if more than 180 degrees
   1.452 +           are covered */
   1.453 +        if (lon_max - lon_min >= 180) {
   1.454 +          cluster->bounding.center.lat = 0;
   1.455 +          cluster->bounding.center.lon = 0;
   1.456 +          cluster->bounding.radius = INFINITY;
   1.457 +          return true;
   1.458 +        }
   1.459 +        /* add point to bounding circle center (for average calculation) */
   1.460 +        cluster->bounding.center.lat += points[j].lat;
   1.461 +        cluster->bounding.center.lon += value;
   1.462 +      }
   1.463 +      /* count total number of points */
   1.464 +      total_npoints += npoints;
   1.465 +    }
   1.466 +    /* determine average latitude and longitude of cluster */
   1.467 +    cluster->bounding.center.lat /= total_npoints;
   1.468 +    cluster->bounding.center.lon /= total_npoints;
   1.469 +    /* normalize longitude of center of cluster bounding circle */
   1.470 +    if (cluster->bounding.center.lon < -180) {
   1.471 +      cluster->bounding.center.lon += 360;
   1.472 +    }
   1.473 +    else if (cluster->bounding.center.lon > 180) {
   1.474 +      cluster->bounding.center.lon -= 360;
   1.475 +    }
   1.476 +    /* round bounding circle center (useful if it is used by other functions) */
   1.477 +    cluster->bounding.center.lat = pgl_round(cluster->bounding.center.lat);
   1.478 +    cluster->bounding.center.lon = pgl_round(cluster->bounding.center.lon);
   1.479 +    /* calculate radius of bounding circle */
   1.480 +    for (i=0; i<cluster->nentries; i++) {
   1.481 +      npoints = cluster->entries[i].npoints;
   1.482 +      points = PGL_ENTRY_POINTS(cluster, i);
   1.483 +      for (j=0; j<npoints; j++) {
   1.484 +        value = pgl_distance(
   1.485 +          cluster->bounding.center.lat, cluster->bounding.center.lon,
   1.486 +          points[j].lat, points[j].lon
   1.487 +        );
   1.488 +        if (value > cluster->bounding.radius) cluster->bounding.radius = value;
   1.489 +      }
   1.490 +    }
   1.491 +  }
   1.492 +  /* return true (east/west orientation is unambiguous) */
   1.493 +  return true;
   1.494 +}
   1.495 +
   1.496 +/* check if point is inside cluster */
   1.497 +/* (if point is on perimeter, then true is returned if and only if
   1.498 +   strict == false) */
   1.499 +static bool pgl_point_in_cluster(
   1.500 +  pgl_point *point,
   1.501 +  pgl_cluster *cluster,
   1.502 +  bool strict
   1.503 +) {
   1.504 +  int i, j, k;  /* i: entry, j: point in entry, k: next point in entry */
   1.505 +  int entrytype;         /* type of entry */
   1.506 +  int npoints;           /* number of points in entry */
   1.507 +  pgl_point *points;     /* array of points in entry */
   1.508 +  int lon_dir = 0;       /* first vertex west (-1) or east (+1) */
   1.509 +  double lon_break = 0;  /* antipodal longitude of first vertex */
   1.510 +  double lat0 = point->lat;  /* latitude of point */
   1.511 +  double lon0;           /* (adjusted) longitude of point */
   1.512 +  double lat1, lon1;     /* latitude and (adjusted) longitude of vertex */
   1.513 +  double lat2, lon2;     /* latitude and (adjusted) longitude of next vertex */
   1.514 +  double lon;            /* longitude of intersection */
   1.515 +  int counter = 0;       /* counter for intersections east of point */
   1.516 +  /* iterate over all entries */
   1.517 +  for (i=0; i<cluster->nentries; i++) {
   1.518 +    /* get type of entry */
   1.519 +    entrytype = cluster->entries[i].entrytype;
   1.520 +    /* skip all entries but polygons if perimeters are excluded */
   1.521 +    if (strict && entrytype != PGL_ENTRY_POLYGON) continue;
   1.522 +    /* get points of entry */
   1.523 +    npoints = cluster->entries[i].npoints;
   1.524 +    points = PGL_ENTRY_POINTS(cluster, i);
   1.525 +    /* determine east/west orientation of first point of entry and calculate
   1.526 +       antipodal longitude */
   1.527 +    lon_break = points[0].lon;
   1.528 +    if      (lon_break < 0) { lon_dir = -1; lon_break += 180; }
   1.529 +    else if (lon_break > 0) { lon_dir =  1; lon_break -= 180; }
   1.530 +    else lon_dir = 0;
   1.531 +    /* get longitude of point */
   1.532 +    lon0 = point->lon;
   1.533 +    /* consider longitude wrap-around for point */
   1.534 +    if      (lon_dir < 0 && lon0 > lon_break) lon0 = pgl_round(lon0 - 360);
   1.535 +    else if (lon_dir > 0 && lon0 < lon_break) lon0 = pgl_round(lon0 + 360);
   1.536 +    /* iterate over all edges and vertices */
   1.537 +    for (j=0; j<npoints; j++) {
   1.538 +      /* return if point is on vertex of polygon */
   1.539 +      if (pgl_point_cmp(point, &(points[j])) == 0) return !strict;
   1.540 +      /* calculate index of next vertex */
   1.541 +      k = (j+1) % npoints;
   1.542 +      /* skip last edge unless entry is (closed) outline or polygon */
   1.543 +      if (
   1.544 +        k == 0 &&
   1.545 +        entrytype != PGL_ENTRY_OUTLINE &&
   1.546 +        entrytype != PGL_ENTRY_POLYGON
   1.547 +      ) continue;
   1.548 +      /* use previously calculated values for lat1 and lon1 if possible */
   1.549 +      if (j) {
   1.550 +        lat1 = lat2;
   1.551 +        lon1 = lon2;
   1.552 +      } else {
   1.553 +        /* otherwise get latitude and longitude values of first vertex */
   1.554 +        lat1 = points[0].lat;
   1.555 +        lon1 = points[0].lon;
   1.556 +        /* and consider longitude wrap-around for first vertex */
   1.557 +        if      (lon_dir < 0 && lon1 > lon_break) lon1 = pgl_round(lon1 - 360);
   1.558 +        else if (lon_dir > 0 && lon1 < lon_break) lon1 = pgl_round(lon1 + 360);
   1.559 +      }
   1.560 +      /* get latitude and longitude of next vertex */
   1.561 +      lat2 = points[k].lat;
   1.562 +      lon2 = points[k].lon;
   1.563 +      /* consider longitude wrap-around for next vertex */
   1.564 +      if      (lon_dir < 0 && lon2 > lon_break) lon2 = pgl_round(lon2 - 360);
   1.565 +      else if (lon_dir > 0 && lon2 < lon_break) lon2 = pgl_round(lon2 + 360);
   1.566 +      /* return if point is on horizontal (west to east) edge of polygon */
   1.567 +      if (
   1.568 +        lat0 == lat1 && lat0 == lat2 &&
   1.569 +        ( (lon0 >= lon1 && lon0 <= lon2) || (lon0 >= lon2 && lon0 <= lon1) )
   1.570 +      ) return !strict;
   1.571 +      /* check if edge crosses east/west line of point */
   1.572 +      if ((lat1 < lat0 && lat2 >= lat0) || (lat2 < lat0 && lat1 >= lat0)) {
   1.573 +        /* calculate longitude of intersection */
   1.574 +        lon = (lon1 * (lat2-lat0) + lon2 * (lat0-lat1)) / (lat2-lat1);
   1.575 +        /* return if intersection goes (approximately) through point */
   1.576 +        if (pgl_round(lon) == lon0) return !strict;
   1.577 +        /* count intersection if east of point and entry is polygon*/
   1.578 +        if (entrytype == PGL_ENTRY_POLYGON && lon > lon0) counter++;
   1.579 +      }
   1.580 +    }
   1.581 +  }
   1.582 +  /* return true if number of intersections is odd */
   1.583 +  return counter & 1;
   1.584 +}
   1.585 +
   1.586 +/* check if all points of the second cluster are strictly inside the first
   1.587 +   cluster */
   1.588 +static inline bool pgl_all_cluster_points_strictly_in_cluster(
   1.589 +  pgl_cluster *outer, pgl_cluster *inner
   1.590 +) {
   1.591 +  int i, j;           /* i: entry, j: point in entry */
   1.592 +  int npoints;        /* number of points in entry */
   1.593 +  pgl_point *points;  /* array of points in entry */
   1.594 +  /* iterate over all entries of "inner" cluster */
   1.595 +  for (i=0; i<inner->nentries; i++) {
   1.596 +    /* get properties of entry */
   1.597 +    npoints = inner->entries[i].npoints;
   1.598 +    points = PGL_ENTRY_POINTS(inner, i);
   1.599 +    /* iterate over all points in entry of "inner" cluster */
   1.600 +    for (j=0; j<npoints; j++) {
   1.601 +      /* return false if one point of inner cluster is not in outer cluster */
   1.602 +      if (!pgl_point_in_cluster(points+j, outer, true)) return false;
   1.603 +    }
   1.604 +  }
   1.605 +  /* otherwise return true */
   1.606 +  return true;
   1.607 +}
   1.608 +
   1.609 +/* check if any point the second cluster is inside the first cluster */
   1.610 +static inline bool pgl_any_cluster_points_in_cluster(
   1.611 +  pgl_cluster *outer, pgl_cluster *inner
   1.612 +) {
   1.613 +  int i, j;           /* i: entry, j: point in entry */
   1.614 +  int npoints;        /* number of points in entry */
   1.615 +  pgl_point *points;  /* array of points in entry */
   1.616 +  /* iterate over all entries of "inner" cluster */
   1.617 +  for (i=0; i<inner->nentries; i++) {
   1.618 +    /* get properties of entry */
   1.619 +    npoints = inner->entries[i].npoints;
   1.620 +    points = PGL_ENTRY_POINTS(inner, i);
   1.621 +    /* iterate over all points in entry of "inner" cluster */
   1.622 +    for (j=0; j<npoints; j++) {
   1.623 +      /* return true if one point of inner cluster is in outer cluster */
   1.624 +      if (pgl_point_in_cluster(points+j, outer, false)) return true;
   1.625 +    }
   1.626 +  }
   1.627 +  /* otherwise return false */
   1.628 +  return false;
   1.629 +}
   1.630 +
   1.631 +/* check if line segment strictly crosses line (not just touching) */
   1.632 +static inline bool pgl_lseg_crosses_line(
   1.633 +  double seg_x1,  double seg_y1,  double seg_x2,  double seg_y2,
   1.634 +  double line_x1, double line_y1, double line_x2, double line_y2
   1.635 +) {
   1.636 +  return (
   1.637 +    (
   1.638 +      (seg_x1-line_x1) * (line_y2-line_y1) -
   1.639 +      (seg_y1-line_y1) * (line_x2-line_x1)
   1.640 +    ) * (
   1.641 +      (seg_x2-line_x1) * (line_y2-line_y1) -
   1.642 +      (seg_y2-line_y1) * (line_x2-line_x1)
   1.643 +    )
   1.644 +  ) < 0;
   1.645 +}
   1.646 +
   1.647 +/* check if paths and outlines of two clusters strictly overlap (not just
   1.648 +   touching) */
   1.649 +static bool pgl_outlines_overlap(
   1.650 +  pgl_cluster *cluster1, pgl_cluster *cluster2
   1.651 +) {
   1.652 +  int i1, j1, k1;  /* i: entry, j: point in entry, k: next point in entry */
   1.653 +  int i2, j2, k2;
   1.654 +  int entrytype1, entrytype2;     /* type of entry */
   1.655 +  int npoints1, npoints2;         /* number of points in entry */
   1.656 +  pgl_point *points1;             /* array of points in entry of cluster1 */
   1.657 +  pgl_point *points2;             /* array of points in entry of cluster2 */
   1.658 +  int lon_dir1, lon_dir2;         /* first vertex west (-1) or east (+1) */
   1.659 +  double lon_break1, lon_break2;  /* antipodal longitude of first vertex */
   1.660 +  double lat11, lon11;  /* latitude and (adjusted) longitude of vertex */
   1.661 +  double lat12, lon12;  /* latitude and (adjusted) longitude of next vertex */
   1.662 +  double lat21, lon21;  /* latitude and (adjusted) longitudes for cluster2 */
   1.663 +  double lat22, lon22;
   1.664 +  double wrapvalue;     /* temporary helper value to adjust wrap-around */  
   1.665 +  /* iterate over all entries of cluster1 */
   1.666 +  for (i1=0; i1<cluster1->nentries; i1++) {
   1.667 +    /* get properties of entry in cluster1 and skip points */
   1.668 +    npoints1 = cluster1->entries[i1].npoints;
   1.669 +    if (npoints1 < 2) continue;
   1.670 +    entrytype1 = cluster1->entries[i1].entrytype;
   1.671 +    points1 = PGL_ENTRY_POINTS(cluster1, i1);
   1.672 +    /* determine east/west orientation of first point and calculate antipodal
   1.673 +       longitude */
   1.674 +    lon_break1 = points1[0].lon;
   1.675 +    if (lon_break1 < 0) {
   1.676 +      lon_dir1   = -1;
   1.677 +      lon_break1 = pgl_round(lon_break1 + 180);
   1.678 +    } else if (lon_break1 > 0) {
   1.679 +      lon_dir1   = 1;
   1.680 +      lon_break1 = pgl_round(lon_break1 - 180);
   1.681 +    } else lon_dir1 = 0;
   1.682 +    /* iterate over all edges and vertices in cluster1 */
   1.683 +    for (j1=0; j1<npoints1; j1++) {
   1.684 +      /* calculate index of next vertex */
   1.685 +      k1 = (j1+1) % npoints1;
   1.686 +      /* skip last edge unless entry is (closed) outline or polygon */
   1.687 +      if (
   1.688 +        k1 == 0 &&
   1.689 +        entrytype1 != PGL_ENTRY_OUTLINE &&
   1.690 +        entrytype1 != PGL_ENTRY_POLYGON
   1.691 +      ) continue;
   1.692 +      /* use previously calculated values for lat1 and lon1 if possible */
   1.693 +      if (j1) {
   1.694 +        lat11 = lat12;
   1.695 +        lon11 = lon12;
   1.696 +      } else {
   1.697 +        /* otherwise get latitude and longitude values of first vertex */
   1.698 +        lat11 = points1[0].lat;
   1.699 +        lon11 = points1[0].lon;
   1.700 +        /* and consider longitude wrap-around for first vertex */
   1.701 +        if      (lon_dir1<0 && lon11>lon_break1) lon11 = pgl_round(lon11-360);
   1.702 +        else if (lon_dir1>0 && lon11<lon_break1) lon11 = pgl_round(lon11+360);
   1.703 +      }
   1.704 +      /* get latitude and longitude of next vertex */
   1.705 +      lat12 = points1[k1].lat;
   1.706 +      lon12 = points1[k1].lon;
   1.707 +      /* consider longitude wrap-around for next vertex */
   1.708 +      if      (lon_dir1<0 && lon12>lon_break1) lon12 = pgl_round(lon12-360);
   1.709 +      else if (lon_dir1>0 && lon12<lon_break1) lon12 = pgl_round(lon12+360);
   1.710 +      /* skip degenerated edges */
   1.711 +      if (lat11 == lat12 && lon11 == lon12) continue;
   1.712 +      /* iterate over all entries of cluster2 */
   1.713 +      for (i2=0; i2<cluster2->nentries; i2++) {
   1.714 +        /* get points and number of points of entry in cluster2 */
   1.715 +        npoints2 = cluster2->entries[i2].npoints;
   1.716 +        if (npoints2 < 2) continue;
   1.717 +        entrytype2 = cluster2->entries[i2].entrytype;
   1.718 +        points2 = PGL_ENTRY_POINTS(cluster2, i2);
   1.719 +        /* determine east/west orientation of first point and calculate antipodal
   1.720 +           longitude */
   1.721 +        lon_break2 = points2[0].lon;
   1.722 +        if (lon_break2 < 0) {
   1.723 +          lon_dir2   = -1;
   1.724 +          lon_break2 = pgl_round(lon_break2 + 180);
   1.725 +        } else if (lon_break2 > 0) {
   1.726 +          lon_dir2   = 1;
   1.727 +          lon_break2 = pgl_round(lon_break2 - 180);
   1.728 +        } else lon_dir2 = 0;
   1.729 +        /* iterate over all edges and vertices in cluster2 */
   1.730 +        for (j2=0; j2<npoints2; j2++) {
   1.731 +          /* calculate index of next vertex */
   1.732 +          k2 = (j2+1) % npoints2;
   1.733 +          /* skip last edge unless entry is (closed) outline or polygon */
   1.734 +          if (
   1.735 +            k2 == 0 &&
   1.736 +            entrytype2 != PGL_ENTRY_OUTLINE &&
   1.737 +            entrytype2 != PGL_ENTRY_POLYGON
   1.738 +          ) continue;
   1.739 +          /* use previously calculated values for lat1 and lon1 if possible */
   1.740 +          if (j2) {
   1.741 +            lat21 = lat22;
   1.742 +            lon21 = lon22;
   1.743 +          } else {
   1.744 +            /* otherwise get latitude and longitude values of first vertex */
   1.745 +            lat21 = points2[0].lat;
   1.746 +            lon21 = points2[0].lon;
   1.747 +            /* and consider longitude wrap-around for first vertex */
   1.748 +            if      (lon_dir2<0 && lon21>lon_break2) lon21 = pgl_round(lon21-360);
   1.749 +            else if (lon_dir2>0 && lon21<lon_break2) lon21 = pgl_round(lon21+360);
   1.750 +          }
   1.751 +          /* get latitude and longitude of next vertex */
   1.752 +          lat22 = points2[k2].lat;
   1.753 +          lon22 = points2[k2].lon;
   1.754 +          /* consider longitude wrap-around for next vertex */
   1.755 +          if      (lon_dir2<0 && lon22>lon_break2) lon22 = pgl_round(lon22-360);
   1.756 +          else if (lon_dir2>0 && lon22<lon_break2) lon22 = pgl_round(lon22+360);
   1.757 +          /* skip degenerated edges */
   1.758 +          if (lat21 == lat22 && lon21 == lon22) continue;
   1.759 +          /* perform another wrap-around where necessary */
   1.760 +          /* TODO: improve performance of whole wrap-around mechanism */
   1.761 +          wrapvalue = (lon21 + lon22) - (lon11 + lon12);
   1.762 +          if (wrapvalue > 360) {
   1.763 +            lon21 = pgl_round(lon21 - 360);
   1.764 +            lon22 = pgl_round(lon22 - 360);
   1.765 +          } else if (wrapvalue < -360) {
   1.766 +            lon21 = pgl_round(lon21 + 360);
   1.767 +            lon22 = pgl_round(lon22 + 360);
   1.768 +          }
   1.769 +          /* return true if segments overlap */
   1.770 +          if (
   1.771 +            pgl_lseg_crosses_line(
   1.772 +              lat11, lon11, lat12, lon12,
   1.773 +              lat21, lon21, lat22, lon22
   1.774 +            ) && pgl_lseg_crosses_line(
   1.775 +              lat21, lon21, lat22, lon22,
   1.776 +              lat11, lon11, lat12, lon12
   1.777 +            )
   1.778 +          ) {
   1.779 +            return true;
   1.780 +          }
   1.781 +        }
   1.782 +      }
   1.783 +    }
   1.784 +  }
   1.785 +  /* otherwise return false */
   1.786 +  return false;
   1.787 +}
   1.788 +
   1.789 +/* check if second cluster is completely contained in first cluster */
   1.790 +static bool pgl_cluster_in_cluster(pgl_cluster *outer, pgl_cluster *inner) {
   1.791 +  if (!pgl_all_cluster_points_strictly_in_cluster(outer, inner)) return false;
   1.792 +  if (pgl_any_cluster_points_in_cluster(inner, outer)) return false;
   1.793 +  if (pgl_outlines_overlap(outer, inner)) return false;
   1.794 +  return true;
   1.795 +}
   1.796 +
   1.797 +/* check if two clusters overlap */
   1.798 +static bool pgl_clusters_overlap(
   1.799 +  pgl_cluster *cluster1, pgl_cluster *cluster2
   1.800 +) {
   1.801 +  if (pgl_any_cluster_points_in_cluster(cluster1, cluster2)) return true;
   1.802 +  if (pgl_any_cluster_points_in_cluster(cluster2, cluster1)) return true;
   1.803 +  if (pgl_outlines_overlap(cluster1, cluster2)) return true;
   1.804 +  return false;
   1.805 +}
   1.806 +
   1.807 +
   1.808 +/* calculate (approximate) distance between point and cluster */
   1.809 +static double pgl_point_cluster_distance(pgl_point *point, pgl_cluster *cluster) {
   1.810 +  double comp;           /* square of compression of meridians */
   1.811 +  int i, j, k;  /* i: entry, j: point in entry, k: next point in entry */
   1.812 +  int entrytype;         /* type of entry */
   1.813 +  int npoints;           /* number of points in entry */
   1.814 +  pgl_point *points;     /* array of points in entry */
   1.815 +  int lon_dir = 0;       /* first vertex west (-1) or east (+1) */
   1.816 +  double lon_break = 0;  /* antipodal longitude of first vertex */
   1.817 +  double lon_min = 0;    /* minimum (adjusted) longitude of entry vertices */
   1.818 +  double lon_max = 0;    /* maximum (adjusted) longitude of entry vertices */
   1.819 +  double lat0 = point->lat;  /* latitude of point */
   1.820 +  double lon0;           /* (adjusted) longitude of point */
   1.821 +  double lat1, lon1;     /* latitude and (adjusted) longitude of vertex */
   1.822 +  double lat2, lon2;     /* latitude and (adjusted) longitude of next vertex */
   1.823 +  double s;              /* scalar for vector calculations */
   1.824 +  double dist;           /* distance calculated in one step */
   1.825 +  double min_dist = INFINITY;   /* minimum distance */
   1.826 +  /* distance is zero if point is contained in cluster */
   1.827 +  if (pgl_point_in_cluster(point, cluster, false)) return 0;
   1.828 +  /* calculate approximate square compression of meridians */
   1.829 +  comp = cos((lat0 / 180.0) * M_PI);
   1.830 +  comp *= comp;
   1.831 +  /* calculate exact square compression of meridians */
   1.832 +  comp *= (
   1.833 +    (1.0 - PGL_EPS2 * (1.0-comp)) *
   1.834 +    (1.0 - PGL_EPS2 * (1.0-comp)) /
   1.835 +    (PGL_SUBEPS2 * PGL_SUBEPS2)
   1.836 +  );
   1.837 +  /* iterate over all entries */
   1.838 +  for (i=0; i<cluster->nentries; i++) {
   1.839 +    /* get properties of entry */
   1.840 +    entrytype = cluster->entries[i].entrytype;
   1.841 +    npoints = cluster->entries[i].npoints;
   1.842 +    points = PGL_ENTRY_POINTS(cluster, i);
   1.843 +    /* determine east/west orientation of first point of entry and calculate
   1.844 +       antipodal longitude */
   1.845 +    lon_break = points[0].lon;
   1.846 +    if      (lon_break < 0) { lon_dir = -1; lon_break += 180; }
   1.847 +    else if (lon_break > 0) { lon_dir =  1; lon_break -= 180; }
   1.848 +    else lon_dir = 0;
   1.849 +    /* determine covered longitude range */
   1.850 +    for (j=0; j<npoints; j++) {
   1.851 +      /* get longitude of vertex */
   1.852 +      lon1 = points[j].lon;
   1.853 +      /* adjust longitude to fix potential wrap-around */
   1.854 +      if      (lon_dir < 0 && lon1 > lon_break) lon1 -= 360;
   1.855 +      else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
   1.856 +      /* update minimum and maximum longitude of polygon */
   1.857 +      if (j == 0 || lon1 < lon_min) lon_min = lon1;
   1.858 +      if (j == 0 || lon1 > lon_max) lon_max = lon1;
   1.859 +    }
   1.860 +    /* adjust longitude wrap-around according to full longitude range */
   1.861 +    lon_break = (lon_max + lon_min) / 2;
   1.862 +    if      (lon_break < 0) { lon_dir = -1; lon_break += 180; }
   1.863 +    else if (lon_break > 0) { lon_dir =  1; lon_break -= 180; }
   1.864 +    /* get longitude of point */
   1.865 +    lon0 = point->lon;
   1.866 +    /* consider longitude wrap-around for point */
   1.867 +    if      (lon_dir < 0 && lon0 > lon_break) lon0 -= 360;
   1.868 +    else if (lon_dir > 0 && lon0 < lon_break) lon0 += 360;
   1.869 +    /* iterate over all edges and vertices */
   1.870 +    for (j=0; j<npoints; j++) {
   1.871 +      /* use previously calculated values for lat1 and lon1 if possible */
   1.872 +      if (j) {
   1.873 +        lat1 = lat2;
   1.874 +        lon1 = lon2;
   1.875 +      } else {
   1.876 +        /* otherwise get latitude and longitude values of first vertex */
   1.877 +        lat1 = points[0].lat;
   1.878 +        lon1 = points[0].lon;
   1.879 +        /* and consider longitude wrap-around for first vertex */
   1.880 +        if      (lon_dir < 0 && lon1 > lon_break) lon1 -= 360;
   1.881 +        else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
   1.882 +      }
   1.883 +      /* calculate distance to vertex */
   1.884 +      dist = pgl_distance(lat0, lon0, lat1, lon1);
   1.885 +      /* store calculated distance if smallest */
   1.886 +      if (dist < min_dist) min_dist = dist;
   1.887 +      /* calculate index of next vertex */
   1.888 +      k = (j+1) % npoints;
   1.889 +      /* skip last edge unless entry is (closed) outline or polygon */
   1.890 +      if (
   1.891 +        k == 0 &&
   1.892 +        entrytype != PGL_ENTRY_OUTLINE &&
   1.893 +        entrytype != PGL_ENTRY_POLYGON
   1.894 +      ) continue;
   1.895 +      /* get latitude and longitude of next vertex */
   1.896 +      lat2 = points[k].lat;
   1.897 +      lon2 = points[k].lon;
   1.898 +      /* consider longitude wrap-around for next vertex */
   1.899 +      if      (lon_dir < 0 && lon2 > lon_break) lon2 -= 360;
   1.900 +      else if (lon_dir > 0 && lon2 < lon_break) lon2 += 360;
   1.901 +      /* go to next vertex and edge if edge is degenerated */
   1.902 +      if (lat1 == lat2 && lon1 == lon2) continue;
   1.903 +      /* otherwise test if point can be projected onto edge of polygon */
   1.904 +      s = (
   1.905 +        ((lat0-lat1) * (lat2-lat1) + comp * (lon0-lon1) * (lon2-lon1)) /
   1.906 +        ((lat2-lat1) * (lat2-lat1) + comp * (lon2-lon1) * (lon2-lon1))
   1.907 +      );
   1.908 +      /* go to next vertex and edge if point cannot be projected */
   1.909 +      if (!(s > 0 && s < 1)) continue;
   1.910 +      /* calculate distance from original point to projected point */
   1.911 +      dist = pgl_distance(
   1.912 +        lat0, lon0,
   1.913 +        lat1 + s * (lat2-lat1),
   1.914 +        lon1 + s * (lon2-lon1)
   1.915 +      );
   1.916 +      /* store calculated distance if smallest */
   1.917 +      if (dist < min_dist) min_dist = dist;
   1.918 +    }
   1.919 +  }
   1.920 +  /* return minimum distance */
   1.921 +  return min_dist;
   1.922 +}
   1.923 +
   1.924 +/* calculate (approximate) distance between two clusters */
   1.925 +static double pgl_cluster_distance(pgl_cluster *cluster1, pgl_cluster *cluster2) {
   1.926 +  int i, j;                    /* i: entry, j: point in entry */
   1.927 +  int npoints;                 /* number of points in entry */
   1.928 +  pgl_point *points;           /* array of points in entry */
   1.929 +  double dist;                 /* distance calculated in one step */
   1.930 +  double min_dist = INFINITY;  /* minimum distance */
   1.931 +  /* consider distance from each point in one cluster to the whole other */
   1.932 +  for (i=0; i<cluster1->nentries; i++) {
   1.933 +    npoints = cluster1->entries[i].npoints;
   1.934 +    points = PGL_ENTRY_POINTS(cluster1, i);
   1.935 +    for (j=0; j<npoints; j++) {
   1.936 +      dist = pgl_point_cluster_distance(points+j, cluster2);
   1.937 +      if (dist == 0) return dist;
   1.938 +      if (dist < min_dist) min_dist = dist;
   1.939 +    }
   1.940 +  }
   1.941 +  /* consider distance from each point in other cluster to the first cluster */
   1.942 +  for (i=0; i<cluster2->nentries; i++) {
   1.943 +    npoints = cluster2->entries[i].npoints;
   1.944 +    points = PGL_ENTRY_POINTS(cluster2, i);
   1.945 +    for (j=0; j<npoints; j++) {
   1.946 +      dist = pgl_point_cluster_distance(points+j, cluster1);
   1.947 +      if (dist == 0) return dist;
   1.948 +      if (dist < min_dist) min_dist = dist;
   1.949 +    }
   1.950 +  }
   1.951 +  return min_dist;
   1.952 +}
   1.953 +
   1.954 +/* estimator function for distance between box and point */
   1.955 +/* always returns a smaller value than actually correct or zero */
   1.956 +static double pgl_estimate_point_box_distance(pgl_point *point, pgl_box *box) {
   1.957 +  double dlon;      /* longitude range of box (delta longitude) */
   1.958 +  double distance;  /* return value */
   1.959 +  /* return infinity if box is empty */
   1.960 +  if (box->lat_min > box->lat_max) return INFINITY;
   1.961 +  /* return zero if point is inside box */
   1.962 +  if (pgl_point_in_box(point, box)) return 0;
   1.963 +  /* calculate delta longitude */
   1.964 +  dlon = box->lon_max - box->lon_min;
   1.965 +  if (dlon < 0) dlon += 360;  /* 180th meridian crossed */
   1.966 +  /* if delta longitude is greater than 150 degrees, perform safe fall-back */
   1.967 +  if (dlon > 150) return 0;
   1.968 +  /* calculate lower limit for distance (formula below requires dlon <= 150) */
   1.969 +  /* TODO: provide better estimation function to improve performance */
   1.970 +  distance = (
   1.971 +    (1.0-1e-14) * pgl_distance(
   1.972 +      point->lat,
   1.973 +      point->lon,
   1.974 +      (box->lat_min + box->lat_max) / 2,
   1.975 +      box->lon_min + dlon/2
   1.976 +    ) - pgl_distance(
   1.977 +      box->lat_min, box->lon_min,
   1.978 +      box->lat_max, box->lon_max
   1.979 +    )
   1.980 +  );
   1.981 +  /* truncate negative results to zero */
   1.982 +  if (distance <= 0) distance = 0;
   1.983 +  /* return result */
   1.984 +  return distance;
   1.985 +}
   1.986 +
   1.987 +
   1.988 +/*-------------------------------*
   1.989 + *  Monte Carlo based functions  *
   1.990 + *-------------------------------*/
   1.991 +
   1.992 +/* half of (spherical) earth's surface area */
   1.993 +#define PGL_HALF_SURFACE (PGL_RADIUS * PGL_DIAMETER * M_PI)
   1.994 +
   1.995 +/* golden angle in radians */
   1.996 +#define PGL_GOLDEN_ANGLE (M_PI * (sqrt(5) - 1.0))
   1.997 +
   1.998 +/* create a list of sample points covering a bounding circle 
   1.999 +   and return covered area */
  1.1000 +static double pgl_sample_points(
  1.1001 +  pgl_point *center,  /* center of bounding circle */
  1.1002 +  double radius,      /* radius of bounding circle */
  1.1003 +  int samples,        /* number of sample points */
  1.1004 +  pgl_point *result   /* pointer to result array */
  1.1005 +) {
  1.1006 +  double double_share = 2.0;  /* double of covered share of earth's surface */
  1.1007 +  double double_share_div_samples;  /* double_share divided by sample count */
  1.1008 +  int i;
  1.1009 +  double t;  /* parameter of spiral laid on (spherical) earth's surface */
  1.1010 +  double x, y, z;  /* normalized coordinates of point on non-rotated spiral */
  1.1011 +  double sin_phi;  /* sine of sph. coordinate of point of non-rotated spiral */
  1.1012 +  double lambda;   /* other sph. coordinate of point of non-rotated spiral */
  1.1013 +  double rot = (0.5 - center->lat / 180.0) * M_PI;  /* needed rot. (in rad) */
  1.1014 +  double cos_rot = cos(rot);  /* cosine of rotation by latitude */
  1.1015 +  double sin_rot = sin(rot);  /* sine of rotation by latitude */
  1.1016 +  double x_rot, z_rot;  /* normalized coordinates of point on rotated spiral */
  1.1017 +  double center_lon = center->lon;  /* second rotation in degree */
  1.1018 +  /* add safety margin to bounding circle because of spherical approximation */
  1.1019 +  radius *= PGL_SPHEROID_A / PGL_RADIUS;
  1.1020 +  /* if whole earth is covered, use initialized value, otherwise calculate
  1.1021 +     share of covered area (multiplied by 2) */
  1.1022 +  if (radius < PGL_MAXDIST) double_share = 1.0 - cos(radius / PGL_RADIUS);
  1.1023 +  /* divide double_share by sample count for later calculations */
  1.1024 +  double_share_div_samples = double_share / samples;
  1.1025 +  /* generate sample points */
  1.1026 +  for (i=0; i<samples; i++) {
  1.1027 +    /* use an offset of 1/2 to avoid too dense clustering at spiral center */
  1.1028 +    t = 0.5 + i;
  1.1029 +    /* calculate normalized coordinates of point on non-rotated spiral */
  1.1030 +    z = 1.0 - double_share_div_samples * t;
  1.1031 +    sin_phi = sqrt(1.0 - z*z);
  1.1032 +    lambda = t * PGL_GOLDEN_ANGLE;
  1.1033 +    x = sin_phi * cos(lambda);
  1.1034 +    y = sin_phi * sin(lambda);
  1.1035 +    /* rotate spiral by latitude value of bounding circle */
  1.1036 +    x_rot = cos_rot * x + sin_rot * z;
  1.1037 +    z_rot = cos_rot * z - sin_rot * x;
  1.1038 +    /* set resulting sample point in result array */
  1.1039 +    /* (while performing second rotation by bounding circle longitude) */
  1.1040 +    result[i].lat = 180.0 * (atan(z_rot / fabs(x_rot)) / M_PI);
  1.1041 +    result[i].lon = center_lon + 180.0 * (atan2(y, x_rot) / M_PI);
  1.1042 +  }
  1.1043 +  /* return covered area */
  1.1044 +  return PGL_HALF_SURFACE * double_share;
  1.1045 +}
  1.1046 +
  1.1047 +/* calculate covered area using Monte Carlo method */
  1.1048 +/* TODO: inefficient, should be replaced by different method */
  1.1049 +static double pgl_monte_carlo_area(pgl_cluster *cluster, int samples) {
  1.1050 +  pgl_point *points;
  1.1051 +  double area;
  1.1052 +  int i;
  1.1053 +  int matches = 0;
  1.1054 +  if (samples > PGL_CLUSTER_MAXPOINTS) {
  1.1055 +    ereport(ERROR, (
  1.1056 +      errcode(ERRCODE_DATA_EXCEPTION),
  1.1057 +      errmsg(
  1.1058 +        "too many sample points for Monte Carlo method (maximum %i)",
  1.1059 +        PGL_CLUSTER_MAXPOINTS
  1.1060 +      )
  1.1061 +    ));
  1.1062 +  }
  1.1063 +  if (!(cluster->bounding.radius > 0)) return 0;
  1.1064 +  for (i=0; ; i++) {
  1.1065 +    if (i == cluster->nentries) return 0;
  1.1066 +    if (cluster->entries[i].entrytype == PGL_ENTRY_POLYGON) break;
  1.1067 +  }
  1.1068 +  points = palloc(samples * sizeof(pgl_point));
  1.1069 +  area = pgl_sample_points(
  1.1070 +    &cluster->bounding.center,
  1.1071 +    cluster->bounding.radius,
  1.1072 +    samples,
  1.1073 +    points
  1.1074 +  );
  1.1075 +  for (i=0; i<samples; i++) {
  1.1076 +    if (pgl_point_in_cluster(points+i, cluster, true)) matches++;
  1.1077 +  }
  1.1078 +  pfree(points);
  1.1079 +  return area * ((double)matches / samples);
  1.1080 +}
  1.1081 +
  1.1082 +/* fair distance between point and cluster (see README file for explanation) */
  1.1083 +static double pgl_fair_distance(
  1.1084 +  pgl_point *point, pgl_cluster *cluster, int samples
  1.1085 +) {
  1.1086 +  double distance;       /* shortest distance from point to cluster */
  1.1087 +  pgl_point *points;     /* sample points for Monte Carlo method */
  1.1088 +  double area;           /* area covered by sample points */
  1.1089 +  int i;
  1.1090 +  int matches = 0;       /* number of points within enlarged area (multiplied
  1.1091 +                            by two to be able to store half matches) */
  1.1092 +  double result;         /* result determined by Monte Carlo method */
  1.1093 +  /* limit sample count to avoid integer overflows on memory allocation */
  1.1094 +  if (samples > PGL_CLUSTER_MAXPOINTS) {
  1.1095 +    ereport(ERROR, (
  1.1096 +      errcode(ERRCODE_DATA_EXCEPTION),
  1.1097 +      errmsg(
  1.1098 +        "too many sample points for Monte Carlo method (maximum %i)",
  1.1099 +        PGL_CLUSTER_MAXPOINTS
  1.1100 +      )
  1.1101 +    ));
  1.1102 +  }
  1.1103 +  /* calculate shortest distance from point to cluster */
  1.1104 +  distance = pgl_point_cluster_distance(point, cluster);
  1.1105 +  /* if cluster consists of a single point or has no bounding circle with
  1.1106 +      positive radius, simply return distance */
  1.1107 +  if (
  1.1108 +    (cluster->nentries==1 && cluster->entries[0].entrytype==PGL_ENTRY_POINT) ||
  1.1109 +    !(cluster->bounding.radius > 0)
  1.1110 +  ) return distance;
  1.1111 +  /* if cluster consists of two points which are twice as far apart, return
  1.1112 +     distance between point and cluster multiplied by square root of two */
  1.1113 +  if (
  1.1114 +    cluster->nentries == 2 &&
  1.1115 +    cluster->entries[0].entrytype == PGL_ENTRY_POINT &&
  1.1116 +    cluster->entries[1].entrytype == PGL_ENTRY_POINT &&
  1.1117 +    pgl_distance(
  1.1118 +      PGL_ENTRY_POINTS(cluster, 0)[0].lat,
  1.1119 +      PGL_ENTRY_POINTS(cluster, 0)[0].lon,
  1.1120 +      PGL_ENTRY_POINTS(cluster, 1)[0].lat,
  1.1121 +      PGL_ENTRY_POINTS(cluster, 1)[0].lon
  1.1122 +    ) >= 2.0 * distance
  1.1123 +  ) {
  1.1124 +    return distance * M_SQRT2;
  1.1125 +  }
  1.1126 +  /* otherwise create sample points for Monte Carlo method and determine area
  1.1127 +     covered by sample points */
  1.1128 +  points = palloc(samples * sizeof(pgl_point));
  1.1129 +  area = pgl_sample_points(
  1.1130 +    &cluster->bounding.center,
  1.1131 +    cluster->bounding.radius + distance,  /* pad bounding circle by distance */
  1.1132 +    samples,
  1.1133 +    points
  1.1134 +  );
  1.1135 +  /* perform Monte Carlo method */
  1.1136 +  if (distance > 0) {
  1.1137 +    /* point is outside cluster */
  1.1138 +    for (i=0; i<samples; i++) {
  1.1139 +      /* count sample poitns within cluster as half match */
  1.1140 +      if (pgl_point_in_cluster(points+i, cluster, true)) matches += 1;
  1.1141 +      /* count sample points outside of cluster but within cluster grown by
  1.1142 +         distance as full match */
  1.1143 +      else if (
  1.1144 +        pgl_point_cluster_distance(points+i, cluster) < distance
  1.1145 +      ) matches += 2;
  1.1146 +    }
  1.1147 +  } else {
  1.1148 +    /* if point is within cluster,
  1.1149 +       just count sample points within cluster as half match */
  1.1150 +    for (i=0; i<samples; i++) {
  1.1151 +      if (pgl_point_in_cluster(points+i, cluster, true)) matches += 1;
  1.1152 +    }
  1.1153 +  }
  1.1154 +  /* release memory for sample points needed by Monte Carlo method */
  1.1155 +  pfree(points);
  1.1156 +  /* convert area determined by Monte Carlo method into distance
  1.1157 +     (i.e. radius of circle with the same area) */
  1.1158 +  result = sqrt(area * ((double)matches / 2.0 / samples) / M_PI);
  1.1159 +  /* return result only if it is greater than the distance between point and
  1.1160 +     cluster to avoid unexpected results because of errors due to limited
  1.1161 +     precision */
  1.1162 +  if (result > distance) return result;
  1.1163 +  /* otherwise return distance between point and cluster */
  1.1164 +  else return distance;
  1.1165 +}
  1.1166 +
  1.1167 +
  1.1168 +/*-------------------------------------------------*
  1.1169 + *  geographic index based on space-filling curve  *
  1.1170 + *-------------------------------------------------*/
  1.1171 +
  1.1172 +/* number of bytes used for geographic (center) position in keys */
  1.1173 +#define PGL_KEY_LATLON_BYTELEN 7
  1.1174 +
  1.1175 +/* maximum reference value for logarithmic size of geographic objects */
  1.1176 +#define PGL_AREAKEY_REFOBJSIZE (PGL_DIAMETER/3.0)  /* can be tweaked */
  1.1177 +
  1.1178 +/* pointer to index key (either pgl_pointkey or pgl_areakey) */
  1.1179 +typedef unsigned char *pgl_keyptr;
  1.1180 +
  1.1181 +/* index key for points (objects with zero area) on the spheroid */
  1.1182 +/* bit  0..55: interspersed bits of latitude and longitude,
  1.1183 +   bit 56..57: always zero,
  1.1184 +   bit 58..63: node depth in hypothetic (full) tree from 0 to 56 (incl.) */
  1.1185 +typedef unsigned char pgl_pointkey[PGL_KEY_LATLON_BYTELEN+1];
  1.1186 +
  1.1187 +/* index key for geographic objects on spheroid with area greater than zero */
  1.1188 +/* bit  0..55: interspersed bits of latitude and longitude of center point,
  1.1189 +   bit     56: always set to 1,
  1.1190 +   bit 57..63: node depth in hypothetic (full) tree from 0 to (2*56)+1 (incl.),
  1.1191 +   bit 64..71: logarithmic object size from 0 to 56+1 = 57 (incl.), but set to
  1.1192 +               PGL_KEY_OBJSIZE_EMPTY (with interspersed bits = 0 and node depth
  1.1193 +               = 113) for empty objects, and set to PGL_KEY_OBJSIZE_UNIVERSAL
  1.1194 +               (with interspersed bits = 0 and node depth = 0) for keys which
  1.1195 +               cover both empty and non-empty objects */
  1.1196 +
  1.1197 +typedef unsigned char pgl_areakey[PGL_KEY_LATLON_BYTELEN+2];
  1.1198 +
  1.1199 +/* helper macros for reading/writing index keys */
  1.1200 +#define PGL_KEY_NODEDEPTH_OFFSET  PGL_KEY_LATLON_BYTELEN
  1.1201 +#define PGL_KEY_OBJSIZE_OFFSET    (PGL_KEY_NODEDEPTH_OFFSET+1)
  1.1202 +#define PGL_POINTKEY_MAXDEPTH     (PGL_KEY_LATLON_BYTELEN*8)
  1.1203 +#define PGL_AREAKEY_MAXDEPTH      (2*PGL_POINTKEY_MAXDEPTH+1)
  1.1204 +#define PGL_AREAKEY_MAXOBJSIZE    (PGL_POINTKEY_MAXDEPTH+1)
  1.1205 +#define PGL_AREAKEY_TYPEMASK      0x80
  1.1206 +#define PGL_KEY_LATLONBIT(key, n) ((key)[(n)/8] & (0x80 >> ((n)%8)))
  1.1207 +#define PGL_KEY_LATLONBIT_DIFF(key1, key2, n) \
  1.1208 +                                  ( PGL_KEY_LATLONBIT(key1, n) ^ \
  1.1209 +                                    PGL_KEY_LATLONBIT(key2, n) )
  1.1210 +#define PGL_KEY_IS_AREAKEY(key)   ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \
  1.1211 +                                    PGL_AREAKEY_TYPEMASK)
  1.1212 +#define PGL_KEY_NODEDEPTH(key)    ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \
  1.1213 +                                    (PGL_AREAKEY_TYPEMASK-1))
  1.1214 +#define PGL_KEY_OBJSIZE(key)      ((key)[PGL_KEY_OBJSIZE_OFFSET])
  1.1215 +#define PGL_KEY_OBJSIZE_EMPTY     126
  1.1216 +#define PGL_KEY_OBJSIZE_UNIVERSAL 127
  1.1217 +#define PGL_KEY_IS_EMPTY(key)     ( PGL_KEY_IS_AREAKEY(key) && \
  1.1218 +                                    (key)[PGL_KEY_OBJSIZE_OFFSET] == \
  1.1219 +                                    PGL_KEY_OBJSIZE_EMPTY )
  1.1220 +#define PGL_KEY_IS_UNIVERSAL(key) ( PGL_KEY_IS_AREAKEY(key) && \
  1.1221 +                                    (key)[PGL_KEY_OBJSIZE_OFFSET] == \
  1.1222 +                                    PGL_KEY_OBJSIZE_UNIVERSAL )
  1.1223 +
  1.1224 +/* set area key to match empty objects only */
  1.1225 +static void pgl_key_set_empty(pgl_keyptr key) {
  1.1226 +  memset(key, 0, sizeof(pgl_areakey));
  1.1227 +  /* Note: setting node depth to maximum is required for picksplit function */
  1.1228 +  key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH;
  1.1229 +  key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_EMPTY;
  1.1230 +}
  1.1231 +
  1.1232 +/* set area key to match any object (including empty objects) */
  1.1233 +static void pgl_key_set_universal(pgl_keyptr key) {
  1.1234 +  memset(key, 0, sizeof(pgl_areakey));
  1.1235 +  key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK;
  1.1236 +  key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_UNIVERSAL;
  1.1237 +}
  1.1238 +
  1.1239 +/* convert a point on earth into a max-depth key to be used in index */
  1.1240 +static void pgl_point_to_key(pgl_point *point, pgl_keyptr key) {
  1.1241 +  double lat = point->lat;
  1.1242 +  double lon = point->lon;
  1.1243 +  int i;
  1.1244 +  /* clear latitude and longitude bits */
  1.1245 +  memset(key, 0, PGL_KEY_LATLON_BYTELEN);
  1.1246 +  /* set node depth to maximum and type bit to zero */
  1.1247 +  key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_POINTKEY_MAXDEPTH;
  1.1248 +  /* iterate over all latitude/longitude bit pairs */
  1.1249 +  for (i=0; i<PGL_POINTKEY_MAXDEPTH/2; i++) {
  1.1250 +    /* determine latitude bit */
  1.1251 +    if (lat >= 0) {
  1.1252 +      key[i/4] |= 0x80 >> (2*(i%4));
  1.1253 +      lat *= 2; lat -= 90;
  1.1254 +    } else {
  1.1255 +      lat *= 2; lat += 90;
  1.1256 +    }
  1.1257 +    /* determine longitude bit */
  1.1258 +    if (lon >= 0) {
  1.1259 +      key[i/4] |= 0x80 >> (2*(i%4)+1);
  1.1260 +      lon *= 2; lon -= 180;
  1.1261 +    } else {
  1.1262 +      lon *= 2; lon += 180;
  1.1263 +    }
  1.1264 +  }
  1.1265 +}
  1.1266 +
  1.1267 +/* convert a circle on earth into a max-depth key to be used in an index */
  1.1268 +static void pgl_circle_to_key(pgl_circle *circle, pgl_keyptr key) {
  1.1269 +  /* handle special case of empty circle */
  1.1270 +  if (circle->radius < 0) {
  1.1271 +    pgl_key_set_empty(key);
  1.1272 +    return;
  1.1273 +  }
  1.1274 +  /* perform same action as for point keys */
  1.1275 +  pgl_point_to_key(&(circle->center), key);
  1.1276 +  /* but overwrite type and node depth to fit area index key */
  1.1277 +  key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH;
  1.1278 +  /* check if radius is greater than (or equal to) reference size */
  1.1279 +  /* (treat equal values as greater values for numerical safety) */
  1.1280 +  if (circle->radius >= PGL_AREAKEY_REFOBJSIZE) {
  1.1281 +    /* if yes, set logarithmic size to zero */
  1.1282 +    key[PGL_KEY_OBJSIZE_OFFSET] = 0;
  1.1283 +  } else {
  1.1284 +    /* otherwise, determine logarithmic size iteratively */
  1.1285 +    /* (one step is equivalent to a factor of sqrt(2)) */
  1.1286 +    double reference = PGL_AREAKEY_REFOBJSIZE / M_SQRT2;
  1.1287 +    int objsize = 1;
  1.1288 +    while (objsize < PGL_AREAKEY_MAXOBJSIZE) {
  1.1289 +      /* stop when radius is greater than (or equal to) adjusted reference */
  1.1290 +      /* (treat equal values as greater values for numerical safety) */
  1.1291 +      if (circle->radius >= reference) break;
  1.1292 +      reference /= M_SQRT2;
  1.1293 +      objsize++;
  1.1294 +    }
  1.1295 +    /* set logarithmic size to determined value */
  1.1296 +    key[PGL_KEY_OBJSIZE_OFFSET] = objsize;
  1.1297 +  }
  1.1298 +}
  1.1299 +
  1.1300 +/* check if one key is subkey of another key or vice versa */
  1.1301 +static bool pgl_keys_overlap(pgl_keyptr key1, pgl_keyptr key2) {
  1.1302 +  int i;  /* key bit offset (includes both lat/lon and log. obj. size bits) */
  1.1303 +  /* determine smallest depth */
  1.1304 +  int depth1 = PGL_KEY_NODEDEPTH(key1);
  1.1305 +  int depth2 = PGL_KEY_NODEDEPTH(key2);
  1.1306 +  int depth = (depth1 < depth2) ? depth1 : depth2;
  1.1307 +  /* check if keys are area keys (assuming that both keys have same type) */
  1.1308 +  if (PGL_KEY_IS_AREAKEY(key1)) {
  1.1309 +    int j = 0;  /* bit offset for logarithmic object size bits */
  1.1310 +    int k = 0;  /* bit offset for latitude and longitude */
  1.1311 +    /* fetch logarithmic object size information */
  1.1312 +    int objsize1 = PGL_KEY_OBJSIZE(key1);
  1.1313 +    int objsize2 = PGL_KEY_OBJSIZE(key2);
  1.1314 +    /* handle special cases for empty objects (universal and empty keys) */
  1.1315 +    if (
  1.1316 +      objsize1 == PGL_KEY_OBJSIZE_UNIVERSAL ||
  1.1317 +      objsize2 == PGL_KEY_OBJSIZE_UNIVERSAL
  1.1318 +    ) return true;
  1.1319 +    if (
  1.1320 +      objsize1 == PGL_KEY_OBJSIZE_EMPTY ||
  1.1321 +      objsize2 == PGL_KEY_OBJSIZE_EMPTY
  1.1322 +    ) return objsize1 == objsize2;
  1.1323 +    /* iterate through key bits */
  1.1324 +    for (i=0; i<depth; i++) {
  1.1325 +      /* every second bit is a bit describing the object size */
  1.1326 +      if (i%2 == 0) {
  1.1327 +        /* check if object size bit is different in both keys (objsize1 and
  1.1328 +           objsize2 describe the minimum index when object size bit is set) */
  1.1329 +        if (
  1.1330 +          (objsize1 <= j && objsize2 > j) ||
  1.1331 +          (objsize2 <= j && objsize1 > j)
  1.1332 +        ) {
  1.1333 +          /* bit differs, therefore keys are in separate branches */
  1.1334 +          return false;
  1.1335 +        }
  1.1336 +        /* increase bit counter for object size bits */
  1.1337 +        j++;
  1.1338 +      }
  1.1339 +      /* all other bits describe latitude and longitude */
  1.1340 +      else {
  1.1341 +        /* check if bit differs in both keys */
  1.1342 +        if (PGL_KEY_LATLONBIT_DIFF(key1, key2, k)) {
  1.1343 +          /* bit differs, therefore keys are in separate branches */
  1.1344 +          return false;
  1.1345 +        }
  1.1346 +        /* increase bit counter for latitude/longitude bits */
  1.1347 +        k++;
  1.1348 +      }
  1.1349 +    }
  1.1350 +  }
  1.1351 +  /* if not, keys are point keys */
  1.1352 +  else {
  1.1353 +    /* iterate through key bits */
  1.1354 +    for (i=0; i<depth; i++) {
  1.1355 +      /* check if bit differs in both keys */
  1.1356 +      if (PGL_KEY_LATLONBIT_DIFF(key1, key2, i)) {
  1.1357 +        /* bit differs, therefore keys are in separate branches */
  1.1358 +        return false;
  1.1359 +      }
  1.1360 +    }
  1.1361 +  }
  1.1362 +  /* return true because keys are in the same branch */
  1.1363 +  return true;
  1.1364 +}
  1.1365 +
  1.1366 +/* combine two keys into new key which covers both original keys */
  1.1367 +/* (result stored in first argument) */
  1.1368 +static void pgl_unite_keys(pgl_keyptr dst, pgl_keyptr src) {
  1.1369 +  int i;  /* key bit offset (includes both lat/lon and log. obj. size bits) */
  1.1370 +  /* determine smallest depth */
  1.1371 +  int depth1 = PGL_KEY_NODEDEPTH(dst);
  1.1372 +  int depth2 = PGL_KEY_NODEDEPTH(src);
  1.1373 +  int depth = (depth1 < depth2) ? depth1 : depth2;
  1.1374 +  /* check if keys are area keys (assuming that both keys have same type) */
  1.1375 +  if (PGL_KEY_IS_AREAKEY(dst)) {
  1.1376 +    pgl_areakey dstbuf = { 0, };  /* destination buffer (cleared) */
  1.1377 +    int j = 0;  /* bit offset for logarithmic object size bits */
  1.1378 +    int k = 0;  /* bit offset for latitude and longitude */
  1.1379 +    /* fetch logarithmic object size information */
  1.1380 +    int objsize1 = PGL_KEY_OBJSIZE(dst);
  1.1381 +    int objsize2 = PGL_KEY_OBJSIZE(src);
  1.1382 +    /* handle special cases for empty objects (universal and empty keys) */
  1.1383 +    if (
  1.1384 +      objsize1 > PGL_AREAKEY_MAXOBJSIZE ||
  1.1385 +      objsize2 > PGL_AREAKEY_MAXOBJSIZE
  1.1386 +    ) {
  1.1387 +      if (
  1.1388 +        objsize1 == PGL_KEY_OBJSIZE_EMPTY &&
  1.1389 +        objsize2 == PGL_KEY_OBJSIZE_EMPTY
  1.1390 +      ) pgl_key_set_empty(dst);
  1.1391 +      else pgl_key_set_universal(dst);
  1.1392 +      return;
  1.1393 +    }
  1.1394 +    /* iterate through key bits */
  1.1395 +    for (i=0; i<depth; i++) {
  1.1396 +      /* every second bit is a bit describing the object size */
  1.1397 +      if (i%2 == 0) {
  1.1398 +        /* increase bit counter for object size bits first */
  1.1399 +        /* (handy when setting objsize variable) */
  1.1400 +        j++;
  1.1401 +        /* check if object size bit is set in neither key */
  1.1402 +        if (objsize1 >= j && objsize2 >= j) {
  1.1403 +          /* set objsize in destination buffer to indicate that size bit is
  1.1404 +             unset in destination buffer at the current bit position */
  1.1405 +          dstbuf[PGL_KEY_OBJSIZE_OFFSET] = j;
  1.1406 +        }
  1.1407 +        /* break if object size bit is set in one key only */
  1.1408 +        else if (objsize1 >= j || objsize2 >= j) break;
  1.1409 +      }
  1.1410 +      /* all other bits describe latitude and longitude */
  1.1411 +      else {
  1.1412 +        /* break if bit differs in both keys */
  1.1413 +        if (PGL_KEY_LATLONBIT(dst, k)) {
  1.1414 +          if (!PGL_KEY_LATLONBIT(src, k)) break;
  1.1415 +          /* but set bit in destination buffer if bit is set in both keys */
  1.1416 +          dstbuf[k/8] |= 0x80 >> (k%8);
  1.1417 +        } else if (PGL_KEY_LATLONBIT(src, k)) break;
  1.1418 +        /* increase bit counter for latitude/longitude bits */
  1.1419 +        k++;
  1.1420 +      }
  1.1421 +    }
  1.1422 +    /* set common node depth and type bit (type bit = 1) */
  1.1423 +    dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | i;
  1.1424 +    /* copy contents of destination buffer to first key */
  1.1425 +    memcpy(dst, dstbuf, sizeof(pgl_areakey));
  1.1426 +  }
  1.1427 +  /* if not, keys are point keys */
  1.1428 +  else {
  1.1429 +    pgl_pointkey dstbuf = { 0, };  /* destination buffer (cleared) */
  1.1430 +    /* iterate through key bits */
  1.1431 +    for (i=0; i<depth; i++) {
  1.1432 +      /* break if bit differs in both keys */
  1.1433 +      if (PGL_KEY_LATLONBIT(dst, i)) {
  1.1434 +        if (!PGL_KEY_LATLONBIT(src, i)) break;
  1.1435 +        /* but set bit in destination buffer if bit is set in both keys */
  1.1436 +        dstbuf[i/8] |= 0x80 >> (i%8);
  1.1437 +      } else if (PGL_KEY_LATLONBIT(src, i)) break;
  1.1438 +    }
  1.1439 +    /* set common node depth (type bit = 0) */
  1.1440 +    dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = i;
  1.1441 +    /* copy contents of destination buffer to first key */
  1.1442 +    memcpy(dst, dstbuf, sizeof(pgl_pointkey));
  1.1443 +  }
  1.1444 +}
  1.1445 +
  1.1446 +/* determine center(!) boundaries and radius estimation of index key */
  1.1447 +static double pgl_key_to_box(pgl_keyptr key, pgl_box *box) {
  1.1448 +  int i;
  1.1449 +  /* determine node depth */
  1.1450 +  int depth = PGL_KEY_NODEDEPTH(key);
  1.1451 +  /* center point of possible result */
  1.1452 +  double lat = 0;
  1.1453 +  double lon = 0;
  1.1454 +  /* maximum distance of real center point from key center */
  1.1455 +  double dlat = 90;
  1.1456 +  double dlon = 180;
  1.1457 +  /* maximum radius of contained objects */
  1.1458 +  double radius = 0;  /* always return zero for point index keys */
  1.1459 +  /* check if key is area key */
  1.1460 +  if (PGL_KEY_IS_AREAKEY(key)) {
  1.1461 +    /* get logarithmic object size */
  1.1462 +    int objsize = PGL_KEY_OBJSIZE(key);
  1.1463 +    /* handle special cases for empty objects (universal and empty keys) */
  1.1464 +    if (objsize == PGL_KEY_OBJSIZE_EMPTY) {
  1.1465 +      pgl_box_set_empty(box);
  1.1466 +      return 0;
  1.1467 +    } else if (objsize == PGL_KEY_OBJSIZE_UNIVERSAL) {
  1.1468 +      box->lat_min = -90;
  1.1469 +      box->lat_max =  90;
  1.1470 +      box->lon_min = -180;
  1.1471 +      box->lon_max =  180;
  1.1472 +      return 0;  /* any value >= 0 would do */
  1.1473 +    }
  1.1474 +    /* calculate maximum possible radius of objects covered by the given key */
  1.1475 +    if (objsize == 0) radius = INFINITY;
  1.1476 +    else {
  1.1477 +      radius = PGL_AREAKEY_REFOBJSIZE;
  1.1478 +      while (--objsize) radius /= M_SQRT2;
  1.1479 +    }
  1.1480 +    /* iterate over latitude and longitude bits in key */
  1.1481 +    /* (every second bit is a latitude or longitude bit) */
  1.1482 +    for (i=0; i<depth/2; i++) {
  1.1483 +      /* check if latitude bit */
  1.1484 +      if (i%2 == 0) {
  1.1485 +        /* cut latitude dimension in half */
  1.1486 +        dlat /= 2;
  1.1487 +        /* increase center latitude if bit is 1, otherwise decrease */
  1.1488 +        if (PGL_KEY_LATLONBIT(key, i)) lat += dlat;
  1.1489 +        else lat -= dlat;
  1.1490 +      }
  1.1491 +      /* otherwise longitude bit */
  1.1492 +      else {
  1.1493 +        /* cut longitude dimension in half */
  1.1494 +        dlon /= 2;
  1.1495 +        /* increase center longitude if bit is 1, otherwise decrease */
  1.1496 +        if (PGL_KEY_LATLONBIT(key, i)) lon += dlon;
  1.1497 +        else lon -= dlon;
  1.1498 +      }
  1.1499 +    }
  1.1500 +  }
  1.1501 +  /* if not, keys are point keys */
  1.1502 +  else {
  1.1503 +    /* iterate over all bits in key */
  1.1504 +    for (i=0; i<depth; i++) {
  1.1505 +      /* check if latitude bit */
  1.1506 +      if (i%2 == 0) {
  1.1507 +        /* cut latitude dimension in half */
  1.1508 +        dlat /= 2;
  1.1509 +        /* increase center latitude if bit is 1, otherwise decrease */
  1.1510 +        if (PGL_KEY_LATLONBIT(key, i)) lat += dlat;
  1.1511 +        else lat -= dlat;
  1.1512 +      }
  1.1513 +      /* otherwise longitude bit */
  1.1514 +      else {
  1.1515 +        /* cut longitude dimension in half */
  1.1516 +        dlon /= 2;
  1.1517 +        /* increase center longitude if bit is 1, otherwise decrease */
  1.1518 +        if (PGL_KEY_LATLONBIT(key, i)) lon += dlon;
  1.1519 +        else lon -= dlon;
  1.1520 +      }
  1.1521 +    }
  1.1522 +  }
  1.1523 +  /* calculate boundaries from center point and remaining dlat and dlon */
  1.1524 +  /* (return values through pointer to box) */
  1.1525 +  box->lat_min = lat - dlat;
  1.1526 +  box->lat_max = lat + dlat;
  1.1527 +  box->lon_min = lon - dlon;
  1.1528 +  box->lon_max = lon + dlon;
  1.1529 +  /* return radius (as a function return value) */
  1.1530 +  return radius;
  1.1531 +}
  1.1532 +
  1.1533 +/* estimator function for distance between point and index key */
  1.1534 +/* always returns a smaller value than actually correct or zero */
  1.1535 +static double pgl_estimate_key_distance(pgl_keyptr key, pgl_point *point) {
  1.1536 +  pgl_box box;  /* center(!) bounding box of area index key */
  1.1537 +  /* calculate center(!) bounding box and maximum radius of objects covered
  1.1538 +     by area index key (radius is zero for point index keys) */
  1.1539 +  double distance = pgl_key_to_box(key, &box);
  1.1540 +  /* calculate estimated distance between bounding box of center point of
  1.1541 +     indexed object and point passed as second argument, then substract maximum
  1.1542 +     radius of objects covered by index key */
  1.1543 +  distance = pgl_estimate_point_box_distance(point, &box) - distance;
  1.1544 +  /* truncate negative results to zero */
  1.1545 +  if (distance <= 0) distance = 0;
  1.1546 +  /* return result */
  1.1547 +  return distance;
  1.1548 +}
  1.1549 +
  1.1550 +
  1.1551 +/*---------------------------------*
  1.1552 + *  helper functions for text I/O  *
  1.1553 + *---------------------------------*/
  1.1554 +
  1.1555 +#define PGL_NUMBUFLEN 64  /* buffer size for number to string conversion */
  1.1556 +
  1.1557 +/* convert floating point number to string (round-trip safe) */
  1.1558 +static void pgl_print_float(char *buf, double flt) {
  1.1559 +  /* check if number is integral */
  1.1560 +  if (trunc(flt) == flt) {
  1.1561 +    /* for integral floats use maximum precision */
  1.1562 +    snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt);
  1.1563 +  } else {
  1.1564 +    /* otherwise check if 15, 16, or 17 digits needed (round-trip safety) */
  1.1565 +    snprintf(buf, PGL_NUMBUFLEN, "%.15g", flt);
  1.1566 +    if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.16g", flt);
  1.1567 +    if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt);
  1.1568 +  }
  1.1569 +}
  1.1570 +
  1.1571 +/* convert latitude floating point number (in degrees) to string */
  1.1572 +static void pgl_print_lat(char *buf, double lat) {
  1.1573 +  if (signbit(lat)) {
  1.1574 +    /* treat negative latitudes (including -0) as south */
  1.1575 +    snprintf(buf, PGL_NUMBUFLEN, "S%015.12f", -lat);
  1.1576 +  } else {
  1.1577 +    /* treat positive latitudes (including +0) as north */
  1.1578 +    snprintf(buf, PGL_NUMBUFLEN, "N%015.12f", lat);
  1.1579 +  }
  1.1580 +}
  1.1581 +
  1.1582 +/* convert longitude floating point number (in degrees) to string */
  1.1583 +static void pgl_print_lon(char *buf, double lon) {
  1.1584 +  if (signbit(lon)) {
  1.1585 +    /* treat negative longitudes (including -0) as west */
  1.1586 +    snprintf(buf, PGL_NUMBUFLEN, "W%016.12f", -lon);
  1.1587 +  } else {
  1.1588 +    /* treat positive longitudes (including +0) as east */
  1.1589 +    snprintf(buf, PGL_NUMBUFLEN, "E%016.12f", lon);
  1.1590 +  }
  1.1591 +}
  1.1592 +
  1.1593 +/* bit masks used as return value of pgl_scan() function */
  1.1594 +#define PGL_SCAN_NONE 0      /* no value has been parsed */
  1.1595 +#define PGL_SCAN_LAT (1<<0)  /* latitude has been parsed */
  1.1596 +#define PGL_SCAN_LON (1<<1)  /* longitude has been parsed */
  1.1597 +#define PGL_SCAN_LATLON (PGL_SCAN_LAT | PGL_SCAN_LON)  /* bitwise OR of both */
  1.1598 +
  1.1599 +/* parse a coordinate (can be latitude or longitude) */
  1.1600 +static int pgl_scan(char **str, double *lat, double *lon) {
  1.1601 +  double val;
  1.1602 +  int len;
  1.1603 +  if (
  1.1604 +    sscanf(*str, " N %lf %n", &val, &len) ||
  1.1605 +    sscanf(*str, " n %lf %n", &val, &len)
  1.1606 +  ) {
  1.1607 +    *str += len; *lat = val; return PGL_SCAN_LAT;
  1.1608 +  }
  1.1609 +  if (
  1.1610 +    sscanf(*str, " S %lf %n", &val, &len) ||
  1.1611 +    sscanf(*str, " s %lf %n", &val, &len)
  1.1612 +  ) {
  1.1613 +    *str += len; *lat = -val; return PGL_SCAN_LAT;
  1.1614 +  }
  1.1615 +  if (
  1.1616 +    sscanf(*str, " E %lf %n", &val, &len) ||
  1.1617 +    sscanf(*str, " e %lf %n", &val, &len)
  1.1618 +  ) {
  1.1619 +    *str += len; *lon = val; return PGL_SCAN_LON;
  1.1620 +  }
  1.1621 +  if (
  1.1622 +    sscanf(*str, " W %lf %n", &val, &len) ||
  1.1623 +    sscanf(*str, " w %lf %n", &val, &len)
  1.1624 +  ) {
  1.1625 +    *str += len; *lon = -val; return PGL_SCAN_LON;
  1.1626 +  }
  1.1627 +  return PGL_SCAN_NONE;
  1.1628 +}
  1.1629 +
  1.1630 +
  1.1631 +/*-----------------*
  1.1632 + *  SQL functions  *
  1.1633 + *-----------------*/
  1.1634 +
  1.1635 +/* Note: These function names use "epoint", "ebox", etc. notation here instead
  1.1636 +   of "point", "box", etc. in order to distinguish them from any previously
  1.1637 +   defined functions. */
  1.1638 +
  1.1639 +/* function needed for dummy types and/or not implemented features */
  1.1640 +PG_FUNCTION_INFO_V1(pgl_notimpl);
  1.1641 +Datum pgl_notimpl(PG_FUNCTION_ARGS) {
  1.1642 +  ereport(ERROR, (errmsg("not implemented by pgLatLon")));
  1.1643 +}
  1.1644 +
  1.1645 +/* set point to latitude and longitude (including checks) */
  1.1646 +static void pgl_epoint_set_latlon(pgl_point *point, double lat, double lon) {
  1.1647 +  /* reject infinite or NaN values */
  1.1648 +  if (!isfinite(lat) || !isfinite(lon)) {
  1.1649 +    ereport(ERROR, (
  1.1650 +      errcode(ERRCODE_DATA_EXCEPTION),
  1.1651 +      errmsg("epoint requires finite coordinates")
  1.1652 +    ));
  1.1653 +  }
  1.1654 +  /* check latitude bounds */
  1.1655 +  if (lat < -90) {
  1.1656 +    ereport(WARNING, (errmsg("latitude exceeds south pole")));
  1.1657 +    lat = -90;
  1.1658 +  } else if (lat > 90) {
  1.1659 +    ereport(WARNING, (errmsg("latitude exceeds north pole")));
  1.1660 +    lat = 90;
  1.1661 +  }
  1.1662 +  /* check longitude bounds */
  1.1663 +  if (lon < -180) {
  1.1664 +    ereport(NOTICE, (errmsg("longitude west of 180th meridian normalized")));
  1.1665 +    lon += 360 - trunc(lon / 360) * 360;
  1.1666 +  } else if (lon > 180) {
  1.1667 +    ereport(NOTICE, (errmsg("longitude east of 180th meridian normalized")));
  1.1668 +    lon -= 360 + trunc(lon / 360) * 360;
  1.1669 +  }
  1.1670 +  /* store rounded latitude/longitude values for round-trip safety */
  1.1671 +  point->lat = pgl_round(lat);
  1.1672 +  point->lon = pgl_round(lon);
  1.1673 +}
  1.1674 +
  1.1675 +/* create point ("epoint" in SQL) from latitude and longitude */
  1.1676 +PG_FUNCTION_INFO_V1(pgl_create_epoint);
  1.1677 +Datum pgl_create_epoint(PG_FUNCTION_ARGS) {
  1.1678 +  pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point));
  1.1679 +  pgl_epoint_set_latlon(point, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1));
  1.1680 +  PG_RETURN_POINTER(point);
  1.1681 +}
  1.1682 +
  1.1683 +/* parse point ("epoint" in SQL) */
  1.1684 +/* format: '[NS]<float> [EW]<float>' */
  1.1685 +PG_FUNCTION_INFO_V1(pgl_epoint_in);
  1.1686 +Datum pgl_epoint_in(PG_FUNCTION_ARGS) {
  1.1687 +  char *str = PG_GETARG_CSTRING(0);  /* input string */
  1.1688 +  char *strptr = str;  /* current position within string */
  1.1689 +  int done = 0;        /* bit mask storing if latitude or longitude was read */
  1.1690 +  double lat, lon;     /* parsed values as double precision floats */
  1.1691 +  pgl_point *point;    /* return value (to be palloc'ed) */
  1.1692 +  /* parse two floats (each latitude or longitude) separated by white-space */
  1.1693 +  done |= pgl_scan(&strptr, &lat, &lon);
  1.1694 +  if (strptr != str && isspace(strptr[-1])) {
  1.1695 +    done |= pgl_scan(&strptr, &lat, &lon);
  1.1696 +  }
  1.1697 +  /* require end of string, and latitude and longitude parsed successfully */
  1.1698 +  if (strptr[0] || done != PGL_SCAN_LATLON) {
  1.1699 +    ereport(ERROR, (
  1.1700 +      errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  1.1701 +      errmsg("invalid input syntax for type epoint: \"%s\"", str)
  1.1702 +    ));
  1.1703 +  }
  1.1704 +  /* allocate memory for result */
  1.1705 +  point = (pgl_point *)palloc(sizeof(pgl_point));
  1.1706 +  /* set latitude and longitude (and perform checks) */
  1.1707 +  pgl_epoint_set_latlon(point, lat, lon);
  1.1708 +  /* return result */
  1.1709 +  PG_RETURN_POINTER(point);
  1.1710 +}
  1.1711 +
  1.1712 +/* create box ("ebox" in SQL) that is empty */
  1.1713 +PG_FUNCTION_INFO_V1(pgl_create_empty_ebox);
  1.1714 +Datum pgl_create_empty_ebox(PG_FUNCTION_ARGS) {
  1.1715 +  pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
  1.1716 +  pgl_box_set_empty(box);
  1.1717 +  PG_RETURN_POINTER(box);
  1.1718 +}
  1.1719 +
  1.1720 +/* set box to given boundaries (including checks) */
  1.1721 +static void pgl_ebox_set_boundaries(
  1.1722 +  pgl_box *box,
  1.1723 +  double lat_min, double lat_max, double lon_min, double lon_max
  1.1724 +) {
  1.1725 +  /* if minimum latitude is greater than maximum latitude, return empty box */
  1.1726 +  if (lat_min > lat_max) {
  1.1727 +    pgl_box_set_empty(box);
  1.1728 +    return;
  1.1729 +  }
  1.1730 +  /* otherwise reject infinite or NaN values */
  1.1731 +  if (
  1.1732 +    !isfinite(lat_min) || !isfinite(lat_max) ||
  1.1733 +    !isfinite(lon_min) || !isfinite(lon_max)
  1.1734 +  ) {
  1.1735 +    ereport(ERROR, (
  1.1736 +      errcode(ERRCODE_DATA_EXCEPTION),
  1.1737 +      errmsg("ebox requires finite coordinates")
  1.1738 +    ));
  1.1739 +  }
  1.1740 +  /* check latitude bounds */
  1.1741 +  if (lat_max < -90) {
  1.1742 +    ereport(WARNING, (errmsg("northern latitude exceeds south pole")));
  1.1743 +    lat_max = -90;
  1.1744 +  } else if (lat_max > 90) {
  1.1745 +    ereport(WARNING, (errmsg("northern latitude exceeds north pole")));
  1.1746 +    lat_max = 90;
  1.1747 +  }
  1.1748 +  if (lat_min < -90) {
  1.1749 +    ereport(WARNING, (errmsg("southern latitude exceeds south pole")));
  1.1750 +    lat_min = -90;
  1.1751 +  } else if (lat_min > 90) {
  1.1752 +    ereport(WARNING, (errmsg("southern latitude exceeds north pole")));
  1.1753 +    lat_min = 90;
  1.1754 +  }
  1.1755 +  /* check if all longitudes are included */
  1.1756 +  if (lon_max - lon_min >= 360) {
  1.1757 +    if (lon_max - lon_min > 360) ereport(WARNING, (
  1.1758 +      errmsg("longitude coverage greater than 360 degrees")
  1.1759 +    ));
  1.1760 +    lon_min = -180;
  1.1761 +    lon_max = 180;
  1.1762 +  } else {
  1.1763 +    /* normalize longitude bounds */
  1.1764 +    if      (lon_min < -180) lon_min += 360 - trunc(lon_min / 360) * 360;
  1.1765 +    else if (lon_min >  180) lon_min -= 360 + trunc(lon_min / 360) * 360;
  1.1766 +    if      (lon_max < -180) lon_max += 360 - trunc(lon_max / 360) * 360;
  1.1767 +    else if (lon_max >  180) lon_max -= 360 + trunc(lon_max / 360) * 360;
  1.1768 +  }
  1.1769 +  /* store rounded latitude/longitude values for round-trip safety */
  1.1770 +  box->lat_min = pgl_round(lat_min);
  1.1771 +  box->lat_max = pgl_round(lat_max);
  1.1772 +  box->lon_min = pgl_round(lon_min);
  1.1773 +  box->lon_max = pgl_round(lon_max);
  1.1774 +  /* ensure that rounding does not change orientation */
  1.1775 +  if (lon_min > lon_max && box->lon_min == box->lon_max) {
  1.1776 +    box->lon_min = -180;
  1.1777 +    box->lon_max = 180;
  1.1778 +  }
  1.1779 +}
  1.1780 +
  1.1781 +/* create box ("ebox" in SQL) from min/max latitude and min/max longitude */
  1.1782 +PG_FUNCTION_INFO_V1(pgl_create_ebox);
  1.1783 +Datum pgl_create_ebox(PG_FUNCTION_ARGS) {
  1.1784 +  pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
  1.1785 +  pgl_ebox_set_boundaries(
  1.1786 +    box,
  1.1787 +    PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1),
  1.1788 +    PG_GETARG_FLOAT8(2), PG_GETARG_FLOAT8(3)
  1.1789 +  );
  1.1790 +  PG_RETURN_POINTER(box);
  1.1791 +}
  1.1792 +
  1.1793 +/* create box ("ebox" in SQL) from two points ("epoint"s) */
  1.1794 +/* (can not be used to cover a longitude range of more than 120 degrees) */
  1.1795 +PG_FUNCTION_INFO_V1(pgl_create_ebox_from_epoints);
  1.1796 +Datum pgl_create_ebox_from_epoints(PG_FUNCTION_ARGS) {
  1.1797 +  pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0);
  1.1798 +  pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1);
  1.1799 +  pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
  1.1800 +  double lat_min, lat_max, lon_min, lon_max;
  1.1801 +  double dlon;  /* longitude range (delta longitude) */
  1.1802 +  /* order latitude and longitude boundaries */
  1.1803 +  if (point2->lat < point1->lat) {
  1.1804 +    lat_min = point2->lat;
  1.1805 +    lat_max = point1->lat;
  1.1806 +  } else {
  1.1807 +    lat_min = point1->lat;
  1.1808 +    lat_max = point2->lat;
  1.1809 +  }
  1.1810 +  if (point2->lon < point1->lon) {
  1.1811 +    lon_min = point2->lon;
  1.1812 +    lon_max = point1->lon;
  1.1813 +  } else {
  1.1814 +    lon_min = point1->lon;
  1.1815 +    lon_max = point2->lon;
  1.1816 +  }
  1.1817 +  /* calculate longitude range (round to avoid floating point errors) */
  1.1818 +  dlon = pgl_round(lon_max - lon_min);
  1.1819 +  /* determine east-west direction */
  1.1820 +  if (dlon >= 240) {
  1.1821 +    /* assume that 180th meridian is crossed and swap min/max longitude */
  1.1822 +    double swap = lon_min; lon_min = lon_max; lon_max = swap;
  1.1823 +  } else if (dlon > 120) {
  1.1824 +    /* unclear orientation since delta longitude > 120 */
  1.1825 +    ereport(ERROR, (
  1.1826 +      errcode(ERRCODE_DATA_EXCEPTION),
  1.1827 +      errmsg("can not determine east/west orientation for ebox")
  1.1828 +    ));
  1.1829 +  }
  1.1830 +  /* use boundaries to setup box (and perform checks) */
  1.1831 +  pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max);
  1.1832 +  /* return result */
  1.1833 +  PG_RETURN_POINTER(box);
  1.1834 +}
  1.1835 +
  1.1836 +/* parse box ("ebox" in SQL) */
  1.1837 +/* format: '[NS]<float> [EW]<float> [NS]<float> [EW]<float>'
  1.1838 +       or: '[NS]<float> [NS]<float> [EW]<float> [EW]<float>' */
  1.1839 +PG_FUNCTION_INFO_V1(pgl_ebox_in);
  1.1840 +Datum pgl_ebox_in(PG_FUNCTION_ARGS) {
  1.1841 +  char *str = PG_GETARG_CSTRING(0);  /* input string */
  1.1842 +  char *str_lower;     /* lower case version of input string */
  1.1843 +  char *strptr;        /* current position within string */
  1.1844 +  int valid;           /* number of valid chars */
  1.1845 +  int done;            /* specifies if latitude or longitude was read */
  1.1846 +  double val;          /* temporary variable */
  1.1847 +  int lat_count = 0;   /* count of latitude values parsed */
  1.1848 +  int lon_count = 0;   /* count of longitufde values parsed */
  1.1849 +  double lat_min, lat_max, lon_min, lon_max;  /* see pgl_box struct */
  1.1850 +  pgl_box *box;        /* return value (to be palloc'ed) */
  1.1851 +  /* lowercase input */
  1.1852 +  str_lower = psprintf("%s", str);
  1.1853 +  for (strptr=str_lower; *strptr; strptr++) {
  1.1854 +    if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A';
  1.1855 +  }
  1.1856 +  /* reset reading position to start of (lowercase) string */
  1.1857 +  strptr = str_lower;
  1.1858 +  /* check if empty box */
  1.1859 +  valid = 0;
  1.1860 +  sscanf(strptr, " empty %n", &valid);
  1.1861 +  if (valid && strptr[valid] == 0) {
  1.1862 +    /* allocate and return empty box */
  1.1863 +    box = (pgl_box *)palloc(sizeof(pgl_box));
  1.1864 +    pgl_box_set_empty(box);
  1.1865 +    PG_RETURN_POINTER(box);
  1.1866 +  }
  1.1867 +  /* demand four blocks separated by whitespace */
  1.1868 +  valid = 0;
  1.1869 +  sscanf(strptr, " %*s %*s %*s %*s %n", &valid);
  1.1870 +  /* if four blocks separated by whitespace exist, parse those blocks */
  1.1871 +  if (strptr[valid] == 0) while (strptr[0]) {
  1.1872 +    /* parse either latitude or longitude (whichever found in input string) */
  1.1873 +    done = pgl_scan(&strptr, &val, &val);
  1.1874 +    /* store latitude or longitude in lat_min, lat_max, lon_min, or lon_max */
  1.1875 +    if (done == PGL_SCAN_LAT) {
  1.1876 +      if (!lat_count) lat_min = val; else lat_max = val;
  1.1877 +      lat_count++;
  1.1878 +    } else if (done == PGL_SCAN_LON) {
  1.1879 +      if (!lon_count) lon_min = val; else lon_max = val;
  1.1880 +      lon_count++;
  1.1881 +    } else {
  1.1882 +      break;
  1.1883 +    }
  1.1884 +  }
  1.1885 +  /* require end of string, and two latitude and two longitude values */
  1.1886 +  if (strptr[0] || lat_count != 2 || lon_count != 2) {
  1.1887 +    ereport(ERROR, (
  1.1888 +      errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  1.1889 +      errmsg("invalid input syntax for type ebox: \"%s\"", str)
  1.1890 +    ));
  1.1891 +  }
  1.1892 +  /* free lower case string */
  1.1893 +  pfree(str_lower);
  1.1894 +  /* order boundaries (maximum greater than minimum) */
  1.1895 +  if (lat_min > lat_max) { val = lat_min; lat_min = lat_max; lat_max = val; }
  1.1896 +  if (lon_min > lon_max) { val = lon_min; lon_min = lon_max; lon_max = val; }
  1.1897 +  /* allocate memory for result */
  1.1898 +  box = (pgl_box *)palloc(sizeof(pgl_box));
  1.1899 +  /* set boundaries (and perform checks) */
  1.1900 +  pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max);
  1.1901 +  /* return result */
  1.1902 +  PG_RETURN_POINTER(box);
  1.1903 +}
  1.1904 +
  1.1905 +/* set circle to given latitude, longitude, and radius (including checks) */
  1.1906 +static void pgl_ecircle_set_latlon_radius(
  1.1907 +  pgl_circle *circle, double lat, double lon, double radius
  1.1908 +) {
  1.1909 +  /* set center point (including checks) */
  1.1910 +  pgl_epoint_set_latlon(&(circle->center), lat, lon);
  1.1911 +  /* handle non-positive radius */
  1.1912 +  if (isnan(radius)) {
  1.1913 +    ereport(ERROR, (
  1.1914 +      errcode(ERRCODE_DATA_EXCEPTION),
  1.1915 +      errmsg("invalid radius for ecircle")
  1.1916 +    ));
  1.1917 +  }
  1.1918 +  if (radius == 0) radius = 0;  /* avoids -0 */
  1.1919 +  else if (radius < 0) {
  1.1920 +    if (isfinite(radius)) {
  1.1921 +      ereport(NOTICE, (errmsg("negative radius converted to minus infinity")));
  1.1922 +    }
  1.1923 +    radius = -INFINITY;
  1.1924 +  }
  1.1925 +  /* store radius (round-trip safety is ensured by pgl_print_float) */
  1.1926 +  circle->radius = radius;
  1.1927 +}
  1.1928 +
  1.1929 +/* create circle ("ecircle" in SQL) from latitude, longitude, and radius */
  1.1930 +PG_FUNCTION_INFO_V1(pgl_create_ecircle);
  1.1931 +Datum pgl_create_ecircle(PG_FUNCTION_ARGS) {
  1.1932 +  pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
  1.1933 +  pgl_ecircle_set_latlon_radius(
  1.1934 +    circle, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1), PG_GETARG_FLOAT8(2)
  1.1935 +  );
  1.1936 +  PG_RETURN_POINTER(circle);
  1.1937 +}
  1.1938 +
  1.1939 +/* create circle ("ecircle" in SQL) from point ("epoint"), and radius */
  1.1940 +PG_FUNCTION_INFO_V1(pgl_create_ecircle_from_epoint);
  1.1941 +Datum pgl_create_ecircle_from_epoint(PG_FUNCTION_ARGS) {
  1.1942 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  1.1943 +  double radius = PG_GETARG_FLOAT8(1);
  1.1944 +  pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
  1.1945 +  /* set latitude, longitude, radius (and perform checks) */
  1.1946 +  pgl_ecircle_set_latlon_radius(circle, point->lat, point->lon, radius);
  1.1947 +  /* return result */
  1.1948 +  PG_RETURN_POINTER(circle);
  1.1949 +}
  1.1950 +
  1.1951 +/* parse circle ("ecircle" in SQL) */
  1.1952 +/* format: '[NS]<float> [EW]<float> <float>' */
  1.1953 +PG_FUNCTION_INFO_V1(pgl_ecircle_in);
  1.1954 +Datum pgl_ecircle_in(PG_FUNCTION_ARGS) {
  1.1955 +  char *str = PG_GETARG_CSTRING(0);  /* input string */
  1.1956 +  char *strptr = str;       /* current position within string */
  1.1957 +  double lat, lon, radius;  /* parsed values as double precision flaots */
  1.1958 +  int valid = 0;            /* number of valid chars */
  1.1959 +  int done = 0;             /* stores if latitude and/or longitude was read */
  1.1960 +  pgl_circle *circle;       /* return value (to be palloc'ed) */
  1.1961 +  /* demand three blocks separated by whitespace */
  1.1962 +  sscanf(strptr, " %*s %*s %*s %n", &valid);
  1.1963 +  /* if three blocks separated by whitespace exist, parse those blocks */
  1.1964 +  if (strptr[valid] == 0) {
  1.1965 +    /* parse latitude and longitude */
  1.1966 +    done |= pgl_scan(&strptr, &lat, &lon);
  1.1967 +    done |= pgl_scan(&strptr, &lat, &lon);
  1.1968 +    /* parse radius (while incrementing strptr by number of bytes parsed) */
  1.1969 +    valid = 0;
  1.1970 +    if (sscanf(strptr, " %lf %n", &radius, &valid) == 1) strptr += valid;
  1.1971 +  }
  1.1972 +  /* require end of string and both latitude and longitude being parsed */
  1.1973 +  if (strptr[0] || done != PGL_SCAN_LATLON) {
  1.1974 +    ereport(ERROR, (
  1.1975 +      errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  1.1976 +      errmsg("invalid input syntax for type ecircle: \"%s\"", str)
  1.1977 +    ));
  1.1978 +  }
  1.1979 +  /* allocate memory for result */
  1.1980 +  circle = (pgl_circle *)palloc(sizeof(pgl_circle));
  1.1981 +  /* set latitude, longitude, radius (and perform checks) */
  1.1982 +  pgl_ecircle_set_latlon_radius(circle, lat, lon, radius);
  1.1983 +  /* return result */
  1.1984 +  PG_RETURN_POINTER(circle);
  1.1985 +}
  1.1986 +
  1.1987 +/* parse cluster ("ecluster" in SQL) */
  1.1988 +PG_FUNCTION_INFO_V1(pgl_ecluster_in);
  1.1989 +Datum pgl_ecluster_in(PG_FUNCTION_ARGS) {
  1.1990 +  int i;
  1.1991 +  char *str = PG_GETARG_CSTRING(0);  /* input string */
  1.1992 +  char *str_lower;         /* lower case version of input string */
  1.1993 +  char *strptr;            /* pointer to current reading position of input */
  1.1994 +  int npoints_total = 0;   /* total number of points in cluster */
  1.1995 +  int nentries = 0;        /* total number of entries */
  1.1996 +  pgl_newentry *entries;   /* array of pgl_newentry to create pgl_cluster */
  1.1997 +  int entries_buflen = 4;  /* maximum number of elements in entries array */
  1.1998 +  int valid;               /* number of valid chars processed */
  1.1999 +  double lat, lon;         /* latitude and longitude of parsed point */
  1.2000 +  int entrytype;           /* current entry type */
  1.2001 +  int npoints;             /* number of points in current entry */
  1.2002 +  pgl_point *points;       /* array of pgl_point for pgl_newentry */
  1.2003 +  int points_buflen;       /* maximum number of elements in points array */
  1.2004 +  int done;                /* return value of pgl_scan function */
  1.2005 +  pgl_cluster *cluster;    /* created cluster */
  1.2006 +  /* lowercase input */
  1.2007 +  str_lower = psprintf("%s", str);
  1.2008 +  for (strptr=str_lower; *strptr; strptr++) {
  1.2009 +    if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A';
  1.2010 +  }
  1.2011 +  /* reset reading position to start of (lowercase) string */
  1.2012 +  strptr = str_lower;
  1.2013 +  /* allocate initial buffer for entries */
  1.2014 +  entries = palloc(entries_buflen * sizeof(pgl_newentry));
  1.2015 +  /* parse until end of string */
  1.2016 +  while (strptr[0]) {
  1.2017 +    /* require previous white-space or closing parenthesis before next token */
  1.2018 +    if (strptr != str_lower && !isspace(strptr[-1]) && strptr[-1] != ')') {
  1.2019 +      goto pgl_ecluster_in_error;
  1.2020 +    }
  1.2021 +    /* ignore token "empty" */
  1.2022 +    valid = 0; sscanf(strptr, " empty %n", &valid);
  1.2023 +    if (valid) { strptr += valid; continue; }
  1.2024 +    /* test for "point" token */
  1.2025 +    valid = 0; sscanf(strptr, " point ( %n", &valid);
  1.2026 +    if (valid) {
  1.2027 +      strptr += valid;
  1.2028 +      entrytype = PGL_ENTRY_POINT;
  1.2029 +      goto pgl_ecluster_in_type_ok;
  1.2030 +    }
  1.2031 +    /* test for "path" token */
  1.2032 +    valid = 0; sscanf(strptr, " path ( %n", &valid);
  1.2033 +    if (valid) {
  1.2034 +      strptr += valid;
  1.2035 +      entrytype = PGL_ENTRY_PATH;
  1.2036 +      goto pgl_ecluster_in_type_ok;
  1.2037 +    }
  1.2038 +    /* test for "outline" token */
  1.2039 +    valid = 0; sscanf(strptr, " outline ( %n", &valid);
  1.2040 +    if (valid) {
  1.2041 +      strptr += valid;
  1.2042 +      entrytype = PGL_ENTRY_OUTLINE;
  1.2043 +      goto pgl_ecluster_in_type_ok;
  1.2044 +    }
  1.2045 +    /* test for "polygon" token */
  1.2046 +    valid = 0; sscanf(strptr, " polygon ( %n", &valid);
  1.2047 +    if (valid) {
  1.2048 +      strptr += valid;
  1.2049 +      entrytype = PGL_ENTRY_POLYGON;
  1.2050 +      goto pgl_ecluster_in_type_ok;
  1.2051 +    }
  1.2052 +    /* error if no valid token found */
  1.2053 +    goto pgl_ecluster_in_error;
  1.2054 +    pgl_ecluster_in_type_ok:
  1.2055 +    /* check if pgl_newentry array needs to grow */
  1.2056 +    if (nentries == entries_buflen) {
  1.2057 +      pgl_newentry *newbuf;
  1.2058 +      entries_buflen *= 2;
  1.2059 +      newbuf = palloc(entries_buflen * sizeof(pgl_newentry));
  1.2060 +      memcpy(newbuf, entries, nentries * sizeof(pgl_newentry));
  1.2061 +      pfree(entries);
  1.2062 +      entries = newbuf;
  1.2063 +    }
  1.2064 +    /* reset number of points for current entry */
  1.2065 +    npoints = 0;
  1.2066 +    /* allocate array for points */
  1.2067 +    points_buflen = 4;
  1.2068 +    points = palloc(points_buflen * sizeof(pgl_point));
  1.2069 +    /* parse until closing parenthesis */
  1.2070 +    while (strptr[0] != ')') {
  1.2071 +      /* error on unexpected end of string */
  1.2072 +      if (strptr[0] == 0) goto pgl_ecluster_in_error;
  1.2073 +      /* mark neither latitude nor longitude as read */
  1.2074 +      done = PGL_SCAN_NONE;
  1.2075 +      /* require white-space before second, third, etc. point */
  1.2076 +      if (npoints != 0 && !isspace(strptr[-1])) goto pgl_ecluster_in_error;
  1.2077 +      /* scan latitude (or longitude) */
  1.2078 +      done |= pgl_scan(&strptr, &lat, &lon);
  1.2079 +      /* require white-space before second coordinate */
  1.2080 +      if (strptr != str && !isspace(strptr[-1])) goto pgl_ecluster_in_error;
  1.2081 +      /* scan longitude (or latitude) */
  1.2082 +      done |= pgl_scan(&strptr, &lat, &lon);
  1.2083 +      /* error unless both latitude and longitude were parsed */
  1.2084 +      if (done != PGL_SCAN_LATLON) goto pgl_ecluster_in_error;
  1.2085 +      /* throw error if number of points is too high */
  1.2086 +      if (npoints_total == PGL_CLUSTER_MAXPOINTS) {
  1.2087 +        ereport(ERROR, (
  1.2088 +          errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  1.2089 +          errmsg(
  1.2090 +            "too many points for ecluster entry (maximum %i)",
  1.2091 +            PGL_CLUSTER_MAXPOINTS
  1.2092 +          )
  1.2093 +        ));
  1.2094 +      }
  1.2095 +      /* check if pgl_point array needs to grow */
  1.2096 +      if (npoints == points_buflen) {
  1.2097 +        pgl_point *newbuf;
  1.2098 +        points_buflen *= 2;
  1.2099 +        newbuf = palloc(points_buflen * sizeof(pgl_point));
  1.2100 +        memcpy(newbuf, points, npoints * sizeof(pgl_point));
  1.2101 +        pfree(points);
  1.2102 +        points = newbuf;
  1.2103 +      }
  1.2104 +      /* append point to pgl_point array (includes checks) */
  1.2105 +      pgl_epoint_set_latlon(&(points[npoints++]), lat, lon);
  1.2106 +      /* increase total number of points */
  1.2107 +      npoints_total++;
  1.2108 +    }
  1.2109 +    /* error if entry has no points */
  1.2110 +    if (!npoints) goto pgl_ecluster_in_error;
  1.2111 +    /* entries with one point are automatically of type "point" */
  1.2112 +    if (npoints == 1) entrytype = PGL_ENTRY_POINT;
  1.2113 +    /* if entries have more than one point */
  1.2114 +    else {
  1.2115 +      /* throw error if entry type is "point" */
  1.2116 +      if (entrytype == PGL_ENTRY_POINT) {
  1.2117 +        ereport(ERROR, (
  1.2118 +          errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  1.2119 +          errmsg("invalid input syntax for type ecluster (point entry with more than one point)")
  1.2120 +        ));
  1.2121 +      }
  1.2122 +      /* coerce outlines and polygons with more than 2 points to be a path */
  1.2123 +      if (npoints == 2) entrytype = PGL_ENTRY_PATH;
  1.2124 +    }
  1.2125 +    /* append entry to pgl_newentry array */
  1.2126 +    entries[nentries].entrytype = entrytype;
  1.2127 +    entries[nentries].npoints = npoints;
  1.2128 +    entries[nentries].points = points;
  1.2129 +    nentries++;
  1.2130 +    /* consume closing parenthesis */
  1.2131 +    strptr++;
  1.2132 +    /* consume white-space */
  1.2133 +    while (isspace(strptr[0])) strptr++;
  1.2134 +  }
  1.2135 +  /* free lower case string */
  1.2136 +  pfree(str_lower);
  1.2137 +  /* create cluster from pgl_newentry array */
  1.2138 +  cluster = pgl_new_cluster(nentries, entries);
  1.2139 +  /* free pgl_newentry array */
  1.2140 +  for (i=0; i<nentries; i++) pfree(entries[i].points);
  1.2141 +  pfree(entries);
  1.2142 +  /* set bounding circle of cluster and check east/west orientation */
  1.2143 +  if (!pgl_finalize_cluster(cluster)) {
  1.2144 +    ereport(ERROR, (
  1.2145 +      errcode(ERRCODE_DATA_EXCEPTION),
  1.2146 +      errmsg("can not determine east/west orientation for ecluster"),
  1.2147 +      errhint("Ensure that each entry has a longitude span of less than 180 degrees.")
  1.2148 +    ));
  1.2149 +  }
  1.2150 +  /* return cluster */
  1.2151 +  PG_RETURN_POINTER(cluster);
  1.2152 +  /* code to throw error */
  1.2153 +  pgl_ecluster_in_error:
  1.2154 +  ereport(ERROR, (
  1.2155 +    errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  1.2156 +    errmsg("invalid input syntax for type ecluster: \"%s\"", str)
  1.2157 +  ));
  1.2158 +}
  1.2159 +
  1.2160 +/* convert point ("epoint") to string representation */
  1.2161 +PG_FUNCTION_INFO_V1(pgl_epoint_out);
  1.2162 +Datum pgl_epoint_out(PG_FUNCTION_ARGS) {
  1.2163 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  1.2164 +  char latstr[PGL_NUMBUFLEN];
  1.2165 +  char lonstr[PGL_NUMBUFLEN];
  1.2166 +  pgl_print_lat(latstr, point->lat);
  1.2167 +  pgl_print_lon(lonstr, point->lon);
  1.2168 +  PG_RETURN_CSTRING(psprintf("%s %s", latstr, lonstr));
  1.2169 +}
  1.2170 +
  1.2171 +/* convert box ("ebox") to string representation */
  1.2172 +PG_FUNCTION_INFO_V1(pgl_ebox_out);
  1.2173 +Datum pgl_ebox_out(PG_FUNCTION_ARGS) {
  1.2174 +  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
  1.2175 +  double lon_min = box->lon_min;
  1.2176 +  double lon_max = box->lon_max;
  1.2177 +  char lat_min_str[PGL_NUMBUFLEN];
  1.2178 +  char lat_max_str[PGL_NUMBUFLEN];
  1.2179 +  char lon_min_str[PGL_NUMBUFLEN];
  1.2180 +  char lon_max_str[PGL_NUMBUFLEN];
  1.2181 +  /* return string "empty" if box is set to be empty */
  1.2182 +  if (box->lat_min > box->lat_max) PG_RETURN_CSTRING("empty");
  1.2183 +  /* use boundaries exceeding W180 or E180 if 180th meridian is enclosed */
  1.2184 +  /* (required since pgl_box_in orders the longitude boundaries) */
  1.2185 +  if (lon_min > lon_max) {
  1.2186 +    if (lon_min + lon_max >= 0) lon_min -= 360;
  1.2187 +    else lon_max += 360;
  1.2188 +  }
  1.2189 +  /* format and return result */
  1.2190 +  pgl_print_lat(lat_min_str, box->lat_min);
  1.2191 +  pgl_print_lat(lat_max_str, box->lat_max);
  1.2192 +  pgl_print_lon(lon_min_str, lon_min);
  1.2193 +  pgl_print_lon(lon_max_str, lon_max);
  1.2194 +  PG_RETURN_CSTRING(psprintf(
  1.2195 +    "%s %s %s %s",
  1.2196 +    lat_min_str, lon_min_str, lat_max_str, lon_max_str
  1.2197 +  ));
  1.2198 +}
  1.2199 +
  1.2200 +/* convert circle ("ecircle") to string representation */
  1.2201 +PG_FUNCTION_INFO_V1(pgl_ecircle_out);
  1.2202 +Datum pgl_ecircle_out(PG_FUNCTION_ARGS) {
  1.2203 +  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
  1.2204 +  char latstr[PGL_NUMBUFLEN];
  1.2205 +  char lonstr[PGL_NUMBUFLEN];
  1.2206 +  char radstr[PGL_NUMBUFLEN];
  1.2207 +  pgl_print_lat(latstr, circle->center.lat);
  1.2208 +  pgl_print_lon(lonstr, circle->center.lon);
  1.2209 +  pgl_print_float(radstr, circle->radius);
  1.2210 +  PG_RETURN_CSTRING(psprintf("%s %s %s", latstr, lonstr, radstr));
  1.2211 +}
  1.2212 +
  1.2213 +/* convert cluster ("ecluster") to string representation */
  1.2214 +PG_FUNCTION_INFO_V1(pgl_ecluster_out);
  1.2215 +Datum pgl_ecluster_out(PG_FUNCTION_ARGS) {
  1.2216 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  1.2217 +  char latstr[PGL_NUMBUFLEN];  /* string buffer for latitude */
  1.2218 +  char lonstr[PGL_NUMBUFLEN];  /* string buffer for longitude */
  1.2219 +  char ***strings;     /* array of array of strings */
  1.2220 +  char *string;        /* string of current token */
  1.2221 +  char *res, *resptr;  /* result and pointer to current write position */
  1.2222 +  size_t reslen = 1;   /* length of result (init with 1 for terminator) */
  1.2223 +  int npoints;         /* number of points of current entry */
  1.2224 +  int i, j;            /* i: entry, j: point in entry */
  1.2225 +  /* handle empty clusters */
  1.2226 +  if (cluster->nentries == 0) {
  1.2227 +    /* free detoasted cluster (if copy) */
  1.2228 +    PG_FREE_IF_COPY(cluster, 0);
  1.2229 +    /* return static result */
  1.2230 +    PG_RETURN_CSTRING("empty");
  1.2231 +  }
  1.2232 +  /* allocate array of array of strings */
  1.2233 +  strings = palloc(cluster->nentries * sizeof(char **));
  1.2234 +  /* iterate over all entries in cluster */
  1.2235 +  for (i=0; i<cluster->nentries; i++) {
  1.2236 +    /* get number of points in entry */
  1.2237 +    npoints = cluster->entries[i].npoints;
  1.2238 +    /* allocate array of strings (one string for each point plus two extra) */
  1.2239 +    strings[i] = palloc((2 + npoints) * sizeof(char *));
  1.2240 +    /* determine opening string */
  1.2241 +    switch (cluster->entries[i].entrytype) {
  1.2242 +      case PGL_ENTRY_POINT:   string = (i==0)?"point ("  :" point (";   break;
  1.2243 +      case PGL_ENTRY_PATH:    string = (i==0)?"path ("   :" path (";    break;
  1.2244 +      case PGL_ENTRY_OUTLINE: string = (i==0)?"outline (":" outline ("; break;
  1.2245 +      case PGL_ENTRY_POLYGON: string = (i==0)?"polygon (":" polygon ("; break;
  1.2246 +      default:                string = (i==0)?"unknown"  :" unknown";
  1.2247 +    }
  1.2248 +    /* use opening string as first string in array */
  1.2249 +    strings[i][0] = string;
  1.2250 +    /* update result length (for allocating result string later) */
  1.2251 +    reslen += strlen(string);
  1.2252 +    /* iterate over all points */
  1.2253 +    for (j=0; j<npoints; j++) {
  1.2254 +      /* create string representation of point */
  1.2255 +      pgl_print_lat(latstr, PGL_ENTRY_POINTS(cluster, i)[j].lat);
  1.2256 +      pgl_print_lon(lonstr, PGL_ENTRY_POINTS(cluster, i)[j].lon);
  1.2257 +      string = psprintf((j == 0) ? "%s %s" : " %s %s", latstr, lonstr);
  1.2258 +      /* copy string pointer to string array */
  1.2259 +      strings[i][j+1] = string;
  1.2260 +      /* update result length (for allocating result string later) */
  1.2261 +      reslen += strlen(string);
  1.2262 +    }
  1.2263 +    /* use closing parenthesis as last string in array */
  1.2264 +    strings[i][npoints+1] = ")";
  1.2265 +    /* update result length (for allocating result string later) */
  1.2266 +    reslen++;
  1.2267 +  }
  1.2268 +  /* allocate result string */
  1.2269 +  res = palloc(reslen);
  1.2270 +  /* set write pointer to begin of result string */
  1.2271 +  resptr = res;
  1.2272 +  /* copy strings into result string */
  1.2273 +  for (i=0; i<cluster->nentries; i++) {
  1.2274 +    npoints = cluster->entries[i].npoints;
  1.2275 +    for (j=0; j<npoints+2; j++) {
  1.2276 +      string = strings[i][j];
  1.2277 +      strcpy(resptr, string);
  1.2278 +      resptr += strlen(string);
  1.2279 +      /* free strings allocated by psprintf */
  1.2280 +      if (j != 0 && j != npoints+1) pfree(string);
  1.2281 +    }
  1.2282 +    /* free array of strings */
  1.2283 +    pfree(strings[i]);
  1.2284 +  }
  1.2285 +  /* free array of array of strings */
  1.2286 +  pfree(strings);
  1.2287 +  /* free detoasted cluster (if copy) */
  1.2288 +  PG_FREE_IF_COPY(cluster, 0);
  1.2289 +  /* return result */
  1.2290 +  PG_RETURN_CSTRING(res);
  1.2291 +}
  1.2292 +
  1.2293 +/* binary input function for point ("epoint") */
  1.2294 +PG_FUNCTION_INFO_V1(pgl_epoint_recv);
  1.2295 +Datum pgl_epoint_recv(PG_FUNCTION_ARGS) {
  1.2296 +  StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
  1.2297 +  pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point));
  1.2298 +  point->lat = pq_getmsgfloat8(buf);
  1.2299 +  point->lon = pq_getmsgfloat8(buf);
  1.2300 +  PG_RETURN_POINTER(point);
  1.2301 +}
  1.2302 +
  1.2303 +/* binary input function for box ("ebox") */
  1.2304 +PG_FUNCTION_INFO_V1(pgl_ebox_recv);
  1.2305 +Datum pgl_ebox_recv(PG_FUNCTION_ARGS) {
  1.2306 +  StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
  1.2307 +  pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
  1.2308 +  box->lat_min = pq_getmsgfloat8(buf);
  1.2309 +  box->lat_max = pq_getmsgfloat8(buf);
  1.2310 +  box->lon_min = pq_getmsgfloat8(buf);
  1.2311 +  box->lon_max = pq_getmsgfloat8(buf);
  1.2312 +  PG_RETURN_POINTER(box);
  1.2313 +}
  1.2314 +
  1.2315 +/* binary input function for circle ("ecircle") */
  1.2316 +PG_FUNCTION_INFO_V1(pgl_ecircle_recv);
  1.2317 +Datum pgl_ecircle_recv(PG_FUNCTION_ARGS) {
  1.2318 +  StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
  1.2319 +  pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
  1.2320 +  circle->center.lat = pq_getmsgfloat8(buf);
  1.2321 +  circle->center.lon = pq_getmsgfloat8(buf);
  1.2322 +  circle->radius = pq_getmsgfloat8(buf);
  1.2323 +  PG_RETURN_POINTER(circle);
  1.2324 +}
  1.2325 +
  1.2326 +/* TODO: binary receive function for cluster */
  1.2327 +
  1.2328 +/* binary output function for point ("epoint") */
  1.2329 +PG_FUNCTION_INFO_V1(pgl_epoint_send);
  1.2330 +Datum pgl_epoint_send(PG_FUNCTION_ARGS) {
  1.2331 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  1.2332 +  StringInfoData buf;
  1.2333 +  pq_begintypsend(&buf);
  1.2334 +  pq_sendfloat8(&buf, point->lat);
  1.2335 +  pq_sendfloat8(&buf, point->lon);
  1.2336 +  PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
  1.2337 +}
  1.2338 +
  1.2339 +/* binary output function for box ("ebox") */
  1.2340 +PG_FUNCTION_INFO_V1(pgl_ebox_send);
  1.2341 +Datum pgl_ebox_send(PG_FUNCTION_ARGS) {
  1.2342 +  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
  1.2343 +  StringInfoData buf;
  1.2344 +  pq_begintypsend(&buf);
  1.2345 +  pq_sendfloat8(&buf, box->lat_min);
  1.2346 +  pq_sendfloat8(&buf, box->lat_max);
  1.2347 +  pq_sendfloat8(&buf, box->lon_min);
  1.2348 +  pq_sendfloat8(&buf, box->lon_max);
  1.2349 +  PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
  1.2350 +}
  1.2351 +
  1.2352 +/* binary output function for circle ("ecircle") */
  1.2353 +PG_FUNCTION_INFO_V1(pgl_ecircle_send);
  1.2354 +Datum pgl_ecircle_send(PG_FUNCTION_ARGS) {
  1.2355 +  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
  1.2356 +  StringInfoData buf;
  1.2357 +  pq_begintypsend(&buf);
  1.2358 +  pq_sendfloat8(&buf, circle->center.lat);
  1.2359 +  pq_sendfloat8(&buf, circle->center.lon);
  1.2360 +  pq_sendfloat8(&buf, circle->radius);
  1.2361 +  PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
  1.2362 +}
  1.2363 +
  1.2364 +/* TODO: binary send functions for cluster */
  1.2365 +
  1.2366 +/* cast point ("epoint") to box ("ebox") */
  1.2367 +PG_FUNCTION_INFO_V1(pgl_epoint_to_ebox);
  1.2368 +Datum pgl_epoint_to_ebox(PG_FUNCTION_ARGS) {
  1.2369 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  1.2370 +  pgl_box *box = palloc(sizeof(pgl_box));
  1.2371 +  box->lat_min = point->lat;
  1.2372 +  box->lat_max = point->lat;
  1.2373 +  box->lon_min = point->lon;
  1.2374 +  box->lon_max = point->lon;
  1.2375 +  PG_RETURN_POINTER(box);
  1.2376 +}
  1.2377 +
  1.2378 +/* cast point ("epoint") to circle ("ecircle") */
  1.2379 +PG_FUNCTION_INFO_V1(pgl_epoint_to_ecircle);
  1.2380 +Datum pgl_epoint_to_ecircle(PG_FUNCTION_ARGS) {
  1.2381 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  1.2382 +  pgl_circle *circle = palloc(sizeof(pgl_box));
  1.2383 +  circle->center = *point;
  1.2384 +  circle->radius = 0;
  1.2385 +  PG_RETURN_POINTER(circle);
  1.2386 +}
  1.2387 +
  1.2388 +/* cast point ("epoint") to cluster ("ecluster") */
  1.2389 +PG_FUNCTION_INFO_V1(pgl_epoint_to_ecluster);
  1.2390 +Datum pgl_epoint_to_ecluster(PG_FUNCTION_ARGS) {
  1.2391 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  1.2392 +  pgl_newentry entry;
  1.2393 +  pgl_cluster *cluster;
  1.2394 +  entry.entrytype = PGL_ENTRY_POINT;
  1.2395 +  entry.npoints = 1;
  1.2396 +  entry.points = point;
  1.2397 +  cluster = pgl_new_cluster(1, &entry);
  1.2398 +  pgl_finalize_cluster(cluster);  /* NOTE: should not fail */
  1.2399 +  PG_RETURN_POINTER(cluster);
  1.2400 +}
  1.2401 +
  1.2402 +/* cast box ("ebox") to cluster ("ecluster") */
  1.2403 +#define pgl_ebox_to_ecluster_macro(i, a, b) \
  1.2404 +  entries[i].entrytype = PGL_ENTRY_POLYGON; \
  1.2405 +  entries[i].npoints = 4; \
  1.2406 +  entries[i].points = points[i]; \
  1.2407 +  points[i][0].lat = box->lat_min; \
  1.2408 +  points[i][0].lon = (a); \
  1.2409 +  points[i][1].lat = box->lat_min; \
  1.2410 +  points[i][1].lon = (b); \
  1.2411 +  points[i][2].lat = box->lat_max; \
  1.2412 +  points[i][2].lon = (b); \
  1.2413 +  points[i][3].lat = box->lat_max; \
  1.2414 +  points[i][3].lon = (a);
  1.2415 +PG_FUNCTION_INFO_V1(pgl_ebox_to_ecluster);
  1.2416 +Datum pgl_ebox_to_ecluster(PG_FUNCTION_ARGS) {
  1.2417 +  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
  1.2418 +  double lon, dlon;
  1.2419 +  int nentries;
  1.2420 +  pgl_newentry entries[3];
  1.2421 +  pgl_point points[3][4];
  1.2422 +  pgl_cluster *cluster;
  1.2423 +  if (box->lat_min > box->lat_max) {
  1.2424 +    nentries = 0;
  1.2425 +  } else if (box->lon_min > box->lon_max) {
  1.2426 +    if (box->lon_min < 0) {
  1.2427 +      lon = pgl_round((box->lon_min + 180) / 2.0);
  1.2428 +      nentries = 3;
  1.2429 +      pgl_ebox_to_ecluster_macro(0, box->lon_min, lon);
  1.2430 +      pgl_ebox_to_ecluster_macro(1, lon, 180);
  1.2431 +      pgl_ebox_to_ecluster_macro(2, -180, box->lon_max);
  1.2432 +    } else if (box->lon_max > 0) {
  1.2433 +      lon = pgl_round((box->lon_max - 180) / 2.0);
  1.2434 +      nentries = 3;
  1.2435 +      pgl_ebox_to_ecluster_macro(0, box->lon_min, 180);
  1.2436 +      pgl_ebox_to_ecluster_macro(1, -180, lon);
  1.2437 +      pgl_ebox_to_ecluster_macro(2, lon, box->lon_max);
  1.2438 +    } else {
  1.2439 +      nentries = 2;
  1.2440 +      pgl_ebox_to_ecluster_macro(0, box->lon_min, 180);
  1.2441 +      pgl_ebox_to_ecluster_macro(1, -180, box->lon_max);
  1.2442 +    }
  1.2443 +  } else {
  1.2444 +    dlon = pgl_round(box->lon_max - box->lon_min);
  1.2445 +    if (dlon < 180) {
  1.2446 +      nentries = 1;
  1.2447 +      pgl_ebox_to_ecluster_macro(0, box->lon_min, box->lon_max);
  1.2448 +    } else {
  1.2449 +      lon = pgl_round((box->lon_min + box->lon_max) / 2.0);
  1.2450 +      if (
  1.2451 +        pgl_round(lon - box->lon_min) < 180 &&
  1.2452 +        pgl_round(box->lon_max - lon) < 180
  1.2453 +      ) {
  1.2454 +        nentries = 2;
  1.2455 +        pgl_ebox_to_ecluster_macro(0, box->lon_min, lon);
  1.2456 +        pgl_ebox_to_ecluster_macro(1, lon, box->lon_max);
  1.2457 +      } else {
  1.2458 +        nentries = 3;
  1.2459 +        pgl_ebox_to_ecluster_macro(0, box->lon_min, -60);
  1.2460 +        pgl_ebox_to_ecluster_macro(1, -60, 60);
  1.2461 +        pgl_ebox_to_ecluster_macro(2, 60, box->lon_max);
  1.2462 +      }
  1.2463 +    }
  1.2464 +  }
  1.2465 +  cluster = pgl_new_cluster(nentries, entries);
  1.2466 +  pgl_finalize_cluster(cluster);  /* NOTE: should not fail */
  1.2467 +  PG_RETURN_POINTER(cluster);
  1.2468 +}
  1.2469 +
  1.2470 +/* extract latitude from point ("epoint") */
  1.2471 +PG_FUNCTION_INFO_V1(pgl_epoint_lat);
  1.2472 +Datum pgl_epoint_lat(PG_FUNCTION_ARGS) {
  1.2473 +  PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lat);
  1.2474 +}
  1.2475 +
  1.2476 +/* extract longitude from point ("epoint") */
  1.2477 +PG_FUNCTION_INFO_V1(pgl_epoint_lon);
  1.2478 +Datum pgl_epoint_lon(PG_FUNCTION_ARGS) {
  1.2479 +  PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lon);
  1.2480 +}
  1.2481 +
  1.2482 +/* extract minimum latitude from box ("ebox") */
  1.2483 +PG_FUNCTION_INFO_V1(pgl_ebox_lat_min);
  1.2484 +Datum pgl_ebox_lat_min(PG_FUNCTION_ARGS) {
  1.2485 +  PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_min);
  1.2486 +}
  1.2487 +
  1.2488 +/* extract maximum latitude from box ("ebox") */
  1.2489 +PG_FUNCTION_INFO_V1(pgl_ebox_lat_max);
  1.2490 +Datum pgl_ebox_lat_max(PG_FUNCTION_ARGS) {
  1.2491 +  PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_max);
  1.2492 +}
  1.2493 +
  1.2494 +/* extract minimum longitude from box ("ebox") */
  1.2495 +PG_FUNCTION_INFO_V1(pgl_ebox_lon_min);
  1.2496 +Datum pgl_ebox_lon_min(PG_FUNCTION_ARGS) {
  1.2497 +  PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_min);
  1.2498 +}
  1.2499 +
  1.2500 +/* extract maximum longitude from box ("ebox") */
  1.2501 +PG_FUNCTION_INFO_V1(pgl_ebox_lon_max);
  1.2502 +Datum pgl_ebox_lon_max(PG_FUNCTION_ARGS) {
  1.2503 +  PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_max);
  1.2504 +}
  1.2505 +
  1.2506 +/* extract center point from circle ("ecircle") */
  1.2507 +PG_FUNCTION_INFO_V1(pgl_ecircle_center);
  1.2508 +Datum pgl_ecircle_center(PG_FUNCTION_ARGS) {
  1.2509 +  PG_RETURN_POINTER(&(((pgl_circle *)PG_GETARG_POINTER(0))->center));
  1.2510 +}
  1.2511 +
  1.2512 +/* extract radius from circle ("ecircle") */
  1.2513 +PG_FUNCTION_INFO_V1(pgl_ecircle_radius);
  1.2514 +Datum pgl_ecircle_radius(PG_FUNCTION_ARGS) {
  1.2515 +  PG_RETURN_FLOAT8(((pgl_circle *)PG_GETARG_POINTER(0))->radius);
  1.2516 +}
  1.2517 +
  1.2518 +/* check if point is inside box (overlap operator "&&") in SQL */
  1.2519 +PG_FUNCTION_INFO_V1(pgl_epoint_ebox_overlap);
  1.2520 +Datum pgl_epoint_ebox_overlap(PG_FUNCTION_ARGS) {
  1.2521 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  1.2522 +  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(1);
  1.2523 +  PG_RETURN_BOOL(pgl_point_in_box(point, box));
  1.2524 +}
  1.2525 +
  1.2526 +/* check if point is inside circle (overlap operator "&&") in SQL */
  1.2527 +PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_overlap);
  1.2528 +Datum pgl_epoint_ecircle_overlap(PG_FUNCTION_ARGS) {
  1.2529 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  1.2530 +  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
  1.2531 +  PG_RETURN_BOOL(
  1.2532 +    pgl_distance(
  1.2533 +      point->lat, point->lon,
  1.2534 +      circle->center.lat, circle->center.lon
  1.2535 +    ) <= circle->radius
  1.2536 +  );
  1.2537 +}
  1.2538 +
  1.2539 +/* check if point is inside cluster (overlap operator "&&") in SQL */
  1.2540 +PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_overlap);
  1.2541 +Datum pgl_epoint_ecluster_overlap(PG_FUNCTION_ARGS) {
  1.2542 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  1.2543 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  1.2544 +  bool retval;
  1.2545 +  /* points outside bounding circle are always assumed to be non-overlapping
  1.2546 +     (necessary for consistent table and index scans) */
  1.2547 +  if (
  1.2548 +    pgl_distance(
  1.2549 +      point->lat, point->lon,
  1.2550 +      cluster->bounding.center.lat, cluster->bounding.center.lon
  1.2551 +    ) > cluster->bounding.radius
  1.2552 +  ) retval = false;
  1.2553 +  else retval = pgl_point_in_cluster(point, cluster, false);
  1.2554 +  PG_FREE_IF_COPY(cluster, 1);
  1.2555 +  PG_RETURN_BOOL(retval);
  1.2556 +}
  1.2557 +
  1.2558 +/* check if point may be inside cluster (lossy overl. operator "&&+") in SQL */
  1.2559 +PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_may_overlap);
  1.2560 +Datum pgl_epoint_ecluster_may_overlap(PG_FUNCTION_ARGS) {
  1.2561 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  1.2562 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  1.2563 +  bool retval = pgl_distance(
  1.2564 +    point->lat, point->lon,
  1.2565 +    cluster->bounding.center.lat, cluster->bounding.center.lon
  1.2566 +  ) <= cluster->bounding.radius;
  1.2567 +  PG_FREE_IF_COPY(cluster, 1);
  1.2568 +  PG_RETURN_BOOL(retval);
  1.2569 +}
  1.2570 +
  1.2571 +/* check if two boxes overlap (overlap operator "&&") in SQL */
  1.2572 +PG_FUNCTION_INFO_V1(pgl_ebox_overlap);
  1.2573 +Datum pgl_ebox_overlap(PG_FUNCTION_ARGS) {
  1.2574 +  pgl_box *box1 = (pgl_box *)PG_GETARG_POINTER(0);
  1.2575 +  pgl_box *box2 = (pgl_box *)PG_GETARG_POINTER(1);
  1.2576 +  PG_RETURN_BOOL(pgl_boxes_overlap(box1, box2));
  1.2577 +}
  1.2578 +
  1.2579 +/* check if box and circle may overlap (lossy overl. operator "&&+") in SQL */
  1.2580 +PG_FUNCTION_INFO_V1(pgl_ebox_ecircle_may_overlap);
  1.2581 +Datum pgl_ebox_ecircle_may_overlap(PG_FUNCTION_ARGS) {
  1.2582 +  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
  1.2583 +  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
  1.2584 +  PG_RETURN_BOOL(
  1.2585 +    pgl_estimate_point_box_distance(&circle->center, box) <= circle->radius
  1.2586 +  );
  1.2587 +}
  1.2588 +
  1.2589 +/* check if box and cluster may overlap (lossy overl. operator "&&+") in SQL */
  1.2590 +PG_FUNCTION_INFO_V1(pgl_ebox_ecluster_may_overlap);
  1.2591 +Datum pgl_ebox_ecluster_may_overlap(PG_FUNCTION_ARGS) {
  1.2592 +  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
  1.2593 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  1.2594 +  bool retval = pgl_estimate_point_box_distance(
  1.2595 +    &cluster->bounding.center,
  1.2596 +    box
  1.2597 +  ) <= cluster->bounding.radius;
  1.2598 +  PG_FREE_IF_COPY(cluster, 1);
  1.2599 +  PG_RETURN_BOOL(retval);
  1.2600 +}
  1.2601 +
  1.2602 +/* check if two circles overlap (overlap operator "&&") in SQL */
  1.2603 +PG_FUNCTION_INFO_V1(pgl_ecircle_overlap);
  1.2604 +Datum pgl_ecircle_overlap(PG_FUNCTION_ARGS) {
  1.2605 +  pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0);
  1.2606 +  pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1);
  1.2607 +  PG_RETURN_BOOL(
  1.2608 +    pgl_distance(
  1.2609 +      circle1->center.lat, circle1->center.lon,
  1.2610 +      circle2->center.lat, circle2->center.lon
  1.2611 +    ) <= circle1->radius + circle2->radius
  1.2612 +  );
  1.2613 +}
  1.2614 +
  1.2615 +/* check if circle and cluster overlap (overlap operator "&&") in SQL */
  1.2616 +PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_overlap);
  1.2617 +Datum pgl_ecircle_ecluster_overlap(PG_FUNCTION_ARGS) {
  1.2618 +  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
  1.2619 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  1.2620 +  bool retval = (
  1.2621 +    pgl_point_cluster_distance(&(circle->center), cluster) <= circle->radius
  1.2622 +  );
  1.2623 +  PG_FREE_IF_COPY(cluster, 1);
  1.2624 +  PG_RETURN_BOOL(retval);
  1.2625 +}
  1.2626 +
  1.2627 +/* check if circle and cluster may overlap (l. ov. operator "&&+") in SQL */
  1.2628 +PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_may_overlap);
  1.2629 +Datum pgl_ecircle_ecluster_may_overlap(PG_FUNCTION_ARGS) {
  1.2630 +  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
  1.2631 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  1.2632 +  bool retval = pgl_distance(
  1.2633 +    circle->center.lat, circle->center.lon,
  1.2634 +    cluster->bounding.center.lat, cluster->bounding.center.lon
  1.2635 +  ) <= circle->radius + cluster->bounding.radius;
  1.2636 +  PG_FREE_IF_COPY(cluster, 1);
  1.2637 +  PG_RETURN_BOOL(retval);
  1.2638 +}
  1.2639 +
  1.2640 +/* check if two clusters overlap (overlap operator "&&") in SQL */
  1.2641 +PG_FUNCTION_INFO_V1(pgl_ecluster_overlap);
  1.2642 +Datum pgl_ecluster_overlap(PG_FUNCTION_ARGS) {
  1.2643 +  pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  1.2644 +  pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  1.2645 +  bool retval;
  1.2646 +  /* clusters with non-touching bounding circles are always assumed to be
  1.2647 +     non-overlapping (improves performance and is necessary for consistent
  1.2648 +     table and index scans) */
  1.2649 +  if (
  1.2650 +    pgl_distance(
  1.2651 +      cluster1->bounding.center.lat, cluster1->bounding.center.lon,
  1.2652 +      cluster2->bounding.center.lat, cluster2->bounding.center.lon
  1.2653 +    ) > cluster1->bounding.radius + cluster2->bounding.radius
  1.2654 +  ) retval = false;
  1.2655 +  else retval = pgl_clusters_overlap(cluster1, cluster2);
  1.2656 +  PG_FREE_IF_COPY(cluster1, 0);
  1.2657 +  PG_FREE_IF_COPY(cluster2, 1);
  1.2658 +  PG_RETURN_BOOL(retval);
  1.2659 +}
  1.2660 +
  1.2661 +/* check if two clusters may overlap (lossy overlap operator "&&+") in SQL */
  1.2662 +PG_FUNCTION_INFO_V1(pgl_ecluster_may_overlap);
  1.2663 +Datum pgl_ecluster_may_overlap(PG_FUNCTION_ARGS) {
  1.2664 +  pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  1.2665 +  pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  1.2666 +  bool retval = pgl_distance(
  1.2667 +    cluster1->bounding.center.lat, cluster1->bounding.center.lon,
  1.2668 +    cluster2->bounding.center.lat, cluster2->bounding.center.lon
  1.2669 +  ) <= cluster1->bounding.radius + cluster2->bounding.radius;
  1.2670 +  PG_FREE_IF_COPY(cluster1, 0);
  1.2671 +  PG_FREE_IF_COPY(cluster2, 1);
  1.2672 +  PG_RETURN_BOOL(retval);
  1.2673 +}
  1.2674 +
  1.2675 +/* check if second cluster is in first cluster (cont. operator "@>) in SQL */
  1.2676 +PG_FUNCTION_INFO_V1(pgl_ecluster_contains);
  1.2677 +Datum pgl_ecluster_contains(PG_FUNCTION_ARGS) {
  1.2678 +  pgl_cluster *outer = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  1.2679 +  pgl_cluster *inner = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  1.2680 +  bool retval;
  1.2681 +  /* clusters with non-touching bounding circles are always assumed to be
  1.2682 +     non-overlapping (improves performance and is necessary for consistent
  1.2683 +     table and index scans) */
  1.2684 +  if (
  1.2685 +    pgl_distance(
  1.2686 +      outer->bounding.center.lat, outer->bounding.center.lon,
  1.2687 +      inner->bounding.center.lat, inner->bounding.center.lon
  1.2688 +    ) > outer->bounding.radius + inner->bounding.radius
  1.2689 +  ) retval = false;
  1.2690 +  else retval = pgl_cluster_in_cluster(outer, inner);
  1.2691 +  PG_FREE_IF_COPY(outer, 0);
  1.2692 +  PG_FREE_IF_COPY(inner, 1);
  1.2693 +  PG_RETURN_BOOL(retval);
  1.2694 +}
  1.2695 +
  1.2696 +/* calculate distance between two points ("<->" operator) in SQL */
  1.2697 +PG_FUNCTION_INFO_V1(pgl_epoint_distance);
  1.2698 +Datum pgl_epoint_distance(PG_FUNCTION_ARGS) {
  1.2699 +  pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0);
  1.2700 +  pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1);
  1.2701 +  PG_RETURN_FLOAT8(pgl_distance(
  1.2702 +    point1->lat, point1->lon, point2->lat, point2->lon
  1.2703 +  ));
  1.2704 +}
  1.2705 +
  1.2706 +/* calculate point to circle distance ("<->" operator) in SQL */
  1.2707 +PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_distance);
  1.2708 +Datum pgl_epoint_ecircle_distance(PG_FUNCTION_ARGS) {
  1.2709 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  1.2710 +  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
  1.2711 +  double distance = pgl_distance(
  1.2712 +    point->lat, point->lon, circle->center.lat, circle->center.lon
  1.2713 +  ) - circle->radius;
  1.2714 +  PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
  1.2715 +}
  1.2716 +
  1.2717 +/* calculate point to cluster distance ("<->" operator) in SQL */
  1.2718 +PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_distance);
  1.2719 +Datum pgl_epoint_ecluster_distance(PG_FUNCTION_ARGS) {
  1.2720 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  1.2721 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  1.2722 +  double distance = pgl_point_cluster_distance(point, cluster);
  1.2723 +  PG_FREE_IF_COPY(cluster, 1);
  1.2724 +  PG_RETURN_FLOAT8(distance);
  1.2725 +}
  1.2726 +
  1.2727 +/* calculate distance between two circles ("<->" operator) in SQL */
  1.2728 +PG_FUNCTION_INFO_V1(pgl_ecircle_distance);
  1.2729 +Datum pgl_ecircle_distance(PG_FUNCTION_ARGS) {
  1.2730 +  pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0);
  1.2731 +  pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1);
  1.2732 +  double distance = pgl_distance(
  1.2733 +    circle1->center.lat, circle1->center.lon,
  1.2734 +    circle2->center.lat, circle2->center.lon
  1.2735 +  ) - (circle1->radius + circle2->radius);
  1.2736 +  PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
  1.2737 +}
  1.2738 +
  1.2739 +/* calculate circle to cluster distance ("<->" operator) in SQL */
  1.2740 +PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_distance);
  1.2741 +Datum pgl_ecircle_ecluster_distance(PG_FUNCTION_ARGS) {
  1.2742 +  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
  1.2743 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  1.2744 +  double distance = (
  1.2745 +    pgl_point_cluster_distance(&(circle->center), cluster) - circle->radius
  1.2746 +  );
  1.2747 +  PG_FREE_IF_COPY(cluster, 1);
  1.2748 +  PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
  1.2749 +}
  1.2750 +
  1.2751 +/* calculate distance between two clusters ("<->" operator) in SQL */
  1.2752 +PG_FUNCTION_INFO_V1(pgl_ecluster_distance);
  1.2753 +Datum pgl_ecluster_distance(PG_FUNCTION_ARGS) {
  1.2754 +  pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  1.2755 +  pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  1.2756 +  double retval = pgl_cluster_distance(cluster1, cluster2);
  1.2757 +  PG_FREE_IF_COPY(cluster1, 0);
  1.2758 +  PG_FREE_IF_COPY(cluster2, 1);
  1.2759 +  PG_RETURN_FLOAT8(retval);
  1.2760 +}
  1.2761 +
  1.2762 +PG_FUNCTION_INFO_V1(pgl_ecluster_monte_carlo_area);
  1.2763 +Datum pgl_ecluster_monte_carlo_area(PG_FUNCTION_ARGS) {
  1.2764 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  1.2765 +  int samples = PG_GETARG_INT32(1);
  1.2766 +  double retval = pgl_monte_carlo_area(cluster, samples);
  1.2767 +  PG_FREE_IF_COPY(cluster, 0);
  1.2768 +  PG_RETURN_FLOAT8(retval);
  1.2769 +}
  1.2770 +
  1.2771 +PG_FUNCTION_INFO_V1(pgl_ecluster_epoint_fair_distance);
  1.2772 +Datum pgl_ecluster_epoint_fair_distance(PG_FUNCTION_ARGS) {
  1.2773 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  1.2774 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(1);
  1.2775 +  int samples = PG_GETARG_INT32(2);
  1.2776 +  double retval = pgl_fair_distance(point, cluster, samples);
  1.2777 +  PG_FREE_IF_COPY(cluster, 0);
  1.2778 +  PG_RETURN_FLOAT8(retval);
  1.2779 +}
  1.2780 +
  1.2781 +
  1.2782 +/*-----------------------------------------------------------*
  1.2783 + *  B-tree comparison operators and index support functions  *
  1.2784 + *-----------------------------------------------------------*/
  1.2785 +
  1.2786 +/* macro for a B-tree operator (without detoasting) */
  1.2787 +#define PGL_BTREE_OPER(func, type, cmpfunc, oper) \
  1.2788 +  PG_FUNCTION_INFO_V1(func); \
  1.2789 +  Datum func(PG_FUNCTION_ARGS) { \
  1.2790 +    type *a = (type *)PG_GETARG_POINTER(0); \
  1.2791 +    type *b = (type *)PG_GETARG_POINTER(1); \
  1.2792 +    PG_RETURN_BOOL(cmpfunc(a, b) oper 0); \
  1.2793 +  }
  1.2794 +
  1.2795 +/* macro for a B-tree comparison function (without detoasting) */
  1.2796 +#define PGL_BTREE_CMP(func, type, cmpfunc) \
  1.2797 +  PG_FUNCTION_INFO_V1(func); \
  1.2798 +  Datum func(PG_FUNCTION_ARGS) { \
  1.2799 +    type *a = (type *)PG_GETARG_POINTER(0); \
  1.2800 +    type *b = (type *)PG_GETARG_POINTER(1); \
  1.2801 +    PG_RETURN_INT32(cmpfunc(a, b)); \
  1.2802 +  }
  1.2803 +
  1.2804 +/* macro for a B-tree operator (with detoasting) */
  1.2805 +#define PGL_BTREE_OPER_DETOAST(func, type, cmpfunc, oper) \
  1.2806 +  PG_FUNCTION_INFO_V1(func); \
  1.2807 +  Datum func(PG_FUNCTION_ARGS) { \
  1.2808 +    bool res; \
  1.2809 +    type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \
  1.2810 +    type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \
  1.2811 +    res = cmpfunc(a, b) oper 0; \
  1.2812 +    PG_FREE_IF_COPY(a, 0); \
  1.2813 +    PG_FREE_IF_COPY(b, 1); \
  1.2814 +    PG_RETURN_BOOL(res); \
  1.2815 +  }
  1.2816 +
  1.2817 +/* macro for a B-tree comparison function (with detoasting) */
  1.2818 +#define PGL_BTREE_CMP_DETOAST(func, type, cmpfunc) \
  1.2819 +  PG_FUNCTION_INFO_V1(func); \
  1.2820 +  Datum func(PG_FUNCTION_ARGS) { \
  1.2821 +    int32_t res; \
  1.2822 +    type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \
  1.2823 +    type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \
  1.2824 +    res = cmpfunc(a, b); \
  1.2825 +    PG_FREE_IF_COPY(a, 0); \
  1.2826 +    PG_FREE_IF_COPY(b, 1); \
  1.2827 +    PG_RETURN_INT32(res); \
  1.2828 +  }
  1.2829 +
  1.2830 +/* B-tree operators and comparison function for point */
  1.2831 +PGL_BTREE_OPER(pgl_btree_epoint_lt, pgl_point, pgl_point_cmp, <)
  1.2832 +PGL_BTREE_OPER(pgl_btree_epoint_le, pgl_point, pgl_point_cmp, <=)
  1.2833 +PGL_BTREE_OPER(pgl_btree_epoint_eq, pgl_point, pgl_point_cmp, ==)
  1.2834 +PGL_BTREE_OPER(pgl_btree_epoint_ne, pgl_point, pgl_point_cmp, !=)
  1.2835 +PGL_BTREE_OPER(pgl_btree_epoint_ge, pgl_point, pgl_point_cmp, >=)
  1.2836 +PGL_BTREE_OPER(pgl_btree_epoint_gt, pgl_point, pgl_point_cmp, >)
  1.2837 +PGL_BTREE_CMP(pgl_btree_epoint_cmp, pgl_point, pgl_point_cmp)
  1.2838 +
  1.2839 +/* B-tree operators and comparison function for box */
  1.2840 +PGL_BTREE_OPER(pgl_btree_ebox_lt, pgl_box, pgl_box_cmp, <)
  1.2841 +PGL_BTREE_OPER(pgl_btree_ebox_le, pgl_box, pgl_box_cmp, <=)
  1.2842 +PGL_BTREE_OPER(pgl_btree_ebox_eq, pgl_box, pgl_box_cmp, ==)
  1.2843 +PGL_BTREE_OPER(pgl_btree_ebox_ne, pgl_box, pgl_box_cmp, !=)
  1.2844 +PGL_BTREE_OPER(pgl_btree_ebox_ge, pgl_box, pgl_box_cmp, >=)
  1.2845 +PGL_BTREE_OPER(pgl_btree_ebox_gt, pgl_box, pgl_box_cmp, >)
  1.2846 +PGL_BTREE_CMP(pgl_btree_ebox_cmp, pgl_box, pgl_box_cmp)
  1.2847 +
  1.2848 +/* B-tree operators and comparison function for circle */
  1.2849 +PGL_BTREE_OPER(pgl_btree_ecircle_lt, pgl_circle, pgl_circle_cmp, <)
  1.2850 +PGL_BTREE_OPER(pgl_btree_ecircle_le, pgl_circle, pgl_circle_cmp, <=)
  1.2851 +PGL_BTREE_OPER(pgl_btree_ecircle_eq, pgl_circle, pgl_circle_cmp, ==)
  1.2852 +PGL_BTREE_OPER(pgl_btree_ecircle_ne, pgl_circle, pgl_circle_cmp, !=)
  1.2853 +PGL_BTREE_OPER(pgl_btree_ecircle_ge, pgl_circle, pgl_circle_cmp, >=)
  1.2854 +PGL_BTREE_OPER(pgl_btree_ecircle_gt, pgl_circle, pgl_circle_cmp, >)
  1.2855 +PGL_BTREE_CMP(pgl_btree_ecircle_cmp, pgl_circle, pgl_circle_cmp)
  1.2856 +
  1.2857 +
  1.2858 +/*--------------------------------*
  1.2859 + *  GiST index support functions  *
  1.2860 + *--------------------------------*/
  1.2861 +
  1.2862 +/* GiST "consistent" support function */
  1.2863 +PG_FUNCTION_INFO_V1(pgl_gist_consistent);
  1.2864 +Datum pgl_gist_consistent(PG_FUNCTION_ARGS) {
  1.2865 +  GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  1.2866 +  pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key);
  1.2867 +  StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2);
  1.2868 +  bool *recheck = (bool *)PG_GETARG_POINTER(4);
  1.2869 +  /* demand recheck because index and query methods are lossy */
  1.2870 +  *recheck = true;
  1.2871 +  /* strategy number aliases for different operators using the same strategy */
  1.2872 +  strategy %= 100;
  1.2873 +  /* strategy number 11: equality of two points */
  1.2874 +  if (strategy == 11) {
  1.2875 +    /* query datum is another point */
  1.2876 +    pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
  1.2877 +    /* convert other point to key */
  1.2878 +    pgl_pointkey querykey;
  1.2879 +    pgl_point_to_key(query, querykey);
  1.2880 +    /* return true if both keys overlap */
  1.2881 +    PG_RETURN_BOOL(pgl_keys_overlap(key, querykey));
  1.2882 +  }
  1.2883 +  /* strategy number 13: equality of two circles */
  1.2884 +  if (strategy == 13) {
  1.2885 +    /* query datum is another circle */
  1.2886 +    pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
  1.2887 +    /* convert other circle to key */
  1.2888 +    pgl_areakey querykey;
  1.2889 +    pgl_circle_to_key(query, querykey);
  1.2890 +    /* return true if both keys overlap */
  1.2891 +    PG_RETURN_BOOL(pgl_keys_overlap(key, querykey));
  1.2892 +  }
  1.2893 +  /* for all remaining strategies, keys on empty objects produce no match */
  1.2894 +  /* (check necessary because query radius may be infinite) */
  1.2895 +  if (PGL_KEY_IS_EMPTY(key)) PG_RETURN_BOOL(false);
  1.2896 +  /* strategy number 21: overlapping with point */
  1.2897 +  if (strategy == 21) {
  1.2898 +    /* query datum is a point */
  1.2899 +    pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
  1.2900 +    /* return true if estimated distance (allowed to be smaller than real
  1.2901 +       distance) between index key and point is zero */
  1.2902 +    PG_RETURN_BOOL(pgl_estimate_key_distance(key, query) == 0);
  1.2903 +  }
  1.2904 +  /* strategy number 22: (point) overlapping with box */
  1.2905 +  if (strategy == 22) {
  1.2906 +    /* query datum is a box */
  1.2907 +    pgl_box *query = (pgl_box *)PG_GETARG_POINTER(1);
  1.2908 +    /* determine bounding box of indexed key */
  1.2909 +    pgl_box keybox;
  1.2910 +    pgl_key_to_box(key, &keybox);
  1.2911 +    /* return true if query box overlaps with bounding box of indexed key */
  1.2912 +    PG_RETURN_BOOL(pgl_boxes_overlap(query, &keybox));
  1.2913 +  }
  1.2914 +  /* strategy number 23: overlapping with circle */
  1.2915 +  if (strategy == 23) {
  1.2916 +    /* query datum is a circle */
  1.2917 +    pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
  1.2918 +    /* return true if estimated distance (allowed to be smaller than real
  1.2919 +       distance) between index key and circle center is smaller than radius */
  1.2920 +    PG_RETURN_BOOL(
  1.2921 +      pgl_estimate_key_distance(key, &(query->center)) <= query->radius
  1.2922 +    );
  1.2923 +  }
  1.2924 +  /* strategy number 24: overlapping with cluster */
  1.2925 +  if (strategy == 24) {
  1.2926 +    bool retval;  /* return value */
  1.2927 +    /* query datum is a cluster */
  1.2928 +    pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  1.2929 +    /* return true if estimated distance (allowed to be smaller than real
  1.2930 +       distance) between index key and circle center is smaller than radius */
  1.2931 +    retval = (
  1.2932 +      pgl_estimate_key_distance(key, &(query->bounding.center)) <=
  1.2933 +      query->bounding.radius
  1.2934 +    );
  1.2935 +    PG_FREE_IF_COPY(query, 1);  /* free detoasted cluster (if copy) */
  1.2936 +    PG_RETURN_BOOL(retval);
  1.2937 +  }
  1.2938 +  /* throw error for any unknown strategy number */
  1.2939 +  elog(ERROR, "unrecognized strategy number: %d", strategy);
  1.2940 +}
  1.2941 +
  1.2942 +/* GiST "union" support function */
  1.2943 +PG_FUNCTION_INFO_V1(pgl_gist_union);
  1.2944 +Datum pgl_gist_union(PG_FUNCTION_ARGS) {
  1.2945 +  GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0);
  1.2946 +  pgl_keyptr out;  /* return value (to be palloc'ed) */
  1.2947 +  int i;
  1.2948 +  /* determine key size */
  1.2949 +  size_t keysize = PGL_KEY_IS_AREAKEY(
  1.2950 +    (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key)
  1.2951 +  ) ? sizeof (pgl_areakey) : sizeof(pgl_pointkey);
  1.2952 +  /* begin with first key as result */
  1.2953 +  out = palloc(keysize);
  1.2954 +  memcpy(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key), keysize);
  1.2955 +  /* unite current result with second, third, etc. key */
  1.2956 +  for (i=1; i<entryvec->n; i++) {
  1.2957 +    pgl_unite_keys(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key));
  1.2958 +  }
  1.2959 +  /* return result */
  1.2960 +  PG_RETURN_POINTER(out);
  1.2961 +}
  1.2962 +
  1.2963 +/* GiST "compress" support function for indicis on points */
  1.2964 +PG_FUNCTION_INFO_V1(pgl_gist_compress_epoint);
  1.2965 +Datum pgl_gist_compress_epoint(PG_FUNCTION_ARGS) {
  1.2966 +  GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  1.2967 +  GISTENTRY *retval;  /* return value (to be palloc'ed unless set to entry) */
  1.2968 +  /* only transform new leaves */
  1.2969 +  if (entry->leafkey) {
  1.2970 +    /* get point to be transformed */
  1.2971 +    pgl_point *point = (pgl_point *)DatumGetPointer(entry->key);
  1.2972 +    /* allocate memory for key */
  1.2973 +    pgl_keyptr key = palloc(sizeof(pgl_pointkey));
  1.2974 +    /* transform point to key */
  1.2975 +    pgl_point_to_key(point, key);
  1.2976 +    /* create new GISTENTRY structure as return value */
  1.2977 +    retval = palloc(sizeof(GISTENTRY));
  1.2978 +    gistentryinit(
  1.2979 +      *retval, PointerGetDatum(key),
  1.2980 +      entry->rel, entry->page, entry->offset, FALSE
  1.2981 +    );
  1.2982 +  } else {
  1.2983 +    /* inner nodes have already been transformed */
  1.2984 +    retval = entry;
  1.2985 +  }
  1.2986 +  /* return pointer to old or new GISTENTRY structure */
  1.2987 +  PG_RETURN_POINTER(retval);
  1.2988 +}
  1.2989 +
  1.2990 +/* GiST "compress" support function for indicis on circles */
  1.2991 +PG_FUNCTION_INFO_V1(pgl_gist_compress_ecircle);
  1.2992 +Datum pgl_gist_compress_ecircle(PG_FUNCTION_ARGS) {
  1.2993 +  GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  1.2994 +  GISTENTRY *retval;  /* return value (to be palloc'ed unless set to entry) */
  1.2995 +  /* only transform new leaves */
  1.2996 +  if (entry->leafkey) {
  1.2997 +    /* get circle to be transformed */
  1.2998 +    pgl_circle *circle = (pgl_circle *)DatumGetPointer(entry->key);
  1.2999 +    /* allocate memory for key */
  1.3000 +    pgl_keyptr key = palloc(sizeof(pgl_areakey));
  1.3001 +    /* transform circle to key */
  1.3002 +    pgl_circle_to_key(circle, key);
  1.3003 +    /* create new GISTENTRY structure as return value */
  1.3004 +    retval = palloc(sizeof(GISTENTRY));
  1.3005 +    gistentryinit(
  1.3006 +      *retval, PointerGetDatum(key),
  1.3007 +      entry->rel, entry->page, entry->offset, FALSE
  1.3008 +    );
  1.3009 +  } else {
  1.3010 +    /* inner nodes have already been transformed */
  1.3011 +    retval = entry;
  1.3012 +  }
  1.3013 +  /* return pointer to old or new GISTENTRY structure */
  1.3014 +  PG_RETURN_POINTER(retval);
  1.3015 +}
  1.3016 +
  1.3017 +/* GiST "compress" support function for indices on clusters */
  1.3018 +PG_FUNCTION_INFO_V1(pgl_gist_compress_ecluster);
  1.3019 +Datum pgl_gist_compress_ecluster(PG_FUNCTION_ARGS) {
  1.3020 +  GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  1.3021 +  GISTENTRY *retval;  /* return value (to be palloc'ed unless set to entry) */
  1.3022 +  /* only transform new leaves */
  1.3023 +  if (entry->leafkey) {
  1.3024 +    /* get cluster to be transformed (detoasting necessary!) */
  1.3025 +    pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(entry->key);
  1.3026 +    /* allocate memory for key */
  1.3027 +    pgl_keyptr key = palloc(sizeof(pgl_areakey));
  1.3028 +    /* transform cluster to key */
  1.3029 +    pgl_circle_to_key(&(cluster->bounding), key);
  1.3030 +    /* create new GISTENTRY structure as return value */
  1.3031 +    retval = palloc(sizeof(GISTENTRY));
  1.3032 +    gistentryinit(
  1.3033 +      *retval, PointerGetDatum(key),
  1.3034 +      entry->rel, entry->page, entry->offset, FALSE
  1.3035 +    );
  1.3036 +    /* free detoasted datum */
  1.3037 +    if ((void *)cluster != (void *)DatumGetPointer(entry->key)) pfree(cluster);
  1.3038 +  } else {
  1.3039 +    /* inner nodes have already been transformed */
  1.3040 +    retval = entry;
  1.3041 +  }
  1.3042 +  /* return pointer to old or new GISTENTRY structure */
  1.3043 +  PG_RETURN_POINTER(retval);
  1.3044 +}
  1.3045 +
  1.3046 +/* GiST "decompress" support function for indices */
  1.3047 +PG_FUNCTION_INFO_V1(pgl_gist_decompress);
  1.3048 +Datum pgl_gist_decompress(PG_FUNCTION_ARGS) {
  1.3049 +  /* return passed pointer without transformation */
  1.3050 +  PG_RETURN_POINTER(PG_GETARG_POINTER(0));
  1.3051 +}
  1.3052 +
  1.3053 +/* GiST "penalty" support function */
  1.3054 +PG_FUNCTION_INFO_V1(pgl_gist_penalty);
  1.3055 +Datum pgl_gist_penalty(PG_FUNCTION_ARGS) {
  1.3056 +  GISTENTRY *origentry = (GISTENTRY *)PG_GETARG_POINTER(0);
  1.3057 +  GISTENTRY *newentry = (GISTENTRY *)PG_GETARG_POINTER(1);
  1.3058 +  float *penalty = (float *)PG_GETARG_POINTER(2);
  1.3059 +  /* get original key and key to insert */
  1.3060 +  pgl_keyptr orig = (pgl_keyptr)DatumGetPointer(origentry->key);
  1.3061 +  pgl_keyptr new = (pgl_keyptr)DatumGetPointer(newentry->key);
  1.3062 +  /* copy original key */
  1.3063 +  union { pgl_pointkey pointkey; pgl_areakey areakey; } union_key;
  1.3064 +  if (PGL_KEY_IS_AREAKEY(orig)) {
  1.3065 +    memcpy(union_key.areakey, orig, sizeof(union_key.areakey));
  1.3066 +  } else {
  1.3067 +    memcpy(union_key.pointkey, orig, sizeof(union_key.pointkey));
  1.3068 +  }
  1.3069 +  /* calculate union of both keys */
  1.3070 +  pgl_unite_keys((pgl_keyptr)&union_key, new);
  1.3071 +  /* penalty equal to reduction of key length (logarithm of added area) */
  1.3072 +  /* (return value by setting referenced value and returning pointer) */
  1.3073 +  *penalty = (
  1.3074 +    PGL_KEY_NODEDEPTH(orig) - PGL_KEY_NODEDEPTH((pgl_keyptr)&union_key)
  1.3075 +  );
  1.3076 +  PG_RETURN_POINTER(penalty);
  1.3077 +}
  1.3078 +
  1.3079 +/* GiST "picksplit" support function */
  1.3080 +PG_FUNCTION_INFO_V1(pgl_gist_picksplit);
  1.3081 +Datum pgl_gist_picksplit(PG_FUNCTION_ARGS) {
  1.3082 +  GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0);
  1.3083 +  GIST_SPLITVEC *v = (GIST_SPLITVEC *)PG_GETARG_POINTER(1);
  1.3084 +  OffsetNumber i;  /* between FirstOffsetNumber and entryvec->n (inclusive) */
  1.3085 +  union {
  1.3086 +    pgl_pointkey pointkey;
  1.3087 +    pgl_areakey areakey;
  1.3088 +  } union_all;  /* union of all keys (to be calculated from scratch)
  1.3089 +                   (later cut in half) */
  1.3090 +  int is_areakey = PGL_KEY_IS_AREAKEY(
  1.3091 +    (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key)
  1.3092 +  );
  1.3093 +  int keysize = is_areakey ? sizeof(pgl_areakey) : sizeof(pgl_pointkey);
  1.3094 +  pgl_keyptr unionL = palloc(keysize);  /* union of keys that go left */
  1.3095 +  pgl_keyptr unionR = palloc(keysize);  /* union of keys that go right */
  1.3096 +  pgl_keyptr key;  /* current key to be processed */
  1.3097 +  /* allocate memory for array of left and right keys, set counts to zero */
  1.3098 +  v->spl_left = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber));
  1.3099 +  v->spl_nleft = 0;
  1.3100 +  v->spl_right = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber));
  1.3101 +  v->spl_nright = 0;
  1.3102 +  /* calculate union of all keys from scratch */
  1.3103 +  memcpy(
  1.3104 +    (pgl_keyptr)&union_all,
  1.3105 +    (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key),
  1.3106 +    keysize
  1.3107 +  );
  1.3108 +  for (i=FirstOffsetNumber+1; i<entryvec->n; i=OffsetNumberNext(i)) {
  1.3109 +    pgl_unite_keys(
  1.3110 +      (pgl_keyptr)&union_all,
  1.3111 +      (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key)
  1.3112 +    );
  1.3113 +  }
  1.3114 +  /* check if trivial split is necessary due to exhausted key length */
  1.3115 +  /* (Note: keys for empty objects must have node depth set to maximum) */
  1.3116 +  if (PGL_KEY_NODEDEPTH((pgl_keyptr)&union_all) == (
  1.3117 +    is_areakey ? PGL_AREAKEY_MAXDEPTH : PGL_POINTKEY_MAXDEPTH
  1.3118 +  )) {
  1.3119 +    /* half of all keys go left */
  1.3120 +    for (
  1.3121 +      i=FirstOffsetNumber;
  1.3122 +      i<FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2;
  1.3123 +      i=OffsetNumberNext(i)
  1.3124 +    ) {
  1.3125 +      /* pointer to current key */
  1.3126 +      key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
  1.3127 +      /* update unionL */
  1.3128 +      /* check if key is first key that goes left */
  1.3129 +      if (!v->spl_nleft) {
  1.3130 +        /* first key that goes left is just copied to unionL */
  1.3131 +        memcpy(unionL, key, keysize);
  1.3132 +      } else {
  1.3133 +        /* unite current value and next key */
  1.3134 +        pgl_unite_keys(unionL, key);
  1.3135 +      }
  1.3136 +      /* append offset number to list of keys that go left */
  1.3137 +      v->spl_left[v->spl_nleft++] = i;
  1.3138 +    }
  1.3139 +    /* other half goes right */
  1.3140 +    for (
  1.3141 +      i=FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2;
  1.3142 +      i<entryvec->n;
  1.3143 +      i=OffsetNumberNext(i)
  1.3144 +    ) {
  1.3145 +      /* pointer to current key */
  1.3146 +      key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
  1.3147 +      /* update unionR */
  1.3148 +      /* check if key is first key that goes right */
  1.3149 +      if (!v->spl_nright) {
  1.3150 +        /* first key that goes right is just copied to unionR */
  1.3151 +        memcpy(unionR, key, keysize);
  1.3152 +      } else {
  1.3153 +        /* unite current value and next key */
  1.3154 +        pgl_unite_keys(unionR, key);
  1.3155 +      }
  1.3156 +      /* append offset number to list of keys that go right */
  1.3157 +      v->spl_right[v->spl_nright++] = i;
  1.3158 +    }
  1.3159 +  }
  1.3160 +  /* otherwise, a non-trivial split is possible */
  1.3161 +  else {
  1.3162 +    /* cut covered area in half */
  1.3163 +    /* (union_all then refers to area of keys that go left) */
  1.3164 +    /* check if union of all keys covers empty and non-empty objects */
  1.3165 +    if (PGL_KEY_IS_UNIVERSAL((pgl_keyptr)&union_all)) {
  1.3166 +      /* if yes, split into empty and non-empty objects */
  1.3167 +      pgl_key_set_empty((pgl_keyptr)&union_all);
  1.3168 +    } else {
  1.3169 +      /* otherwise split by next bit */
  1.3170 +      ((pgl_keyptr)&union_all)[PGL_KEY_NODEDEPTH_OFFSET]++;
  1.3171 +      /* NOTE: type bit conserved */
  1.3172 +    }
  1.3173 +    /* determine for each key if it goes left or right */
  1.3174 +    for (i=FirstOffsetNumber; i<entryvec->n; i=OffsetNumberNext(i)) {
  1.3175 +      /* pointer to current key */
  1.3176 +      key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
  1.3177 +      /* keys within one half of the area go left */
  1.3178 +      if (pgl_keys_overlap((pgl_keyptr)&union_all, key)) {
  1.3179 +        /* update unionL */
  1.3180 +        /* check if key is first key that goes left */
  1.3181 +        if (!v->spl_nleft) {
  1.3182 +          /* first key that goes left is just copied to unionL */
  1.3183 +          memcpy(unionL, key, keysize);
  1.3184 +        } else {
  1.3185 +          /* unite current value of unionL and processed key */
  1.3186 +          pgl_unite_keys(unionL, key);
  1.3187 +        }
  1.3188 +        /* append offset number to list of keys that go left */
  1.3189 +        v->spl_left[v->spl_nleft++] = i;
  1.3190 +      }
  1.3191 +      /* the other keys go right */
  1.3192 +      else {
  1.3193 +        /* update unionR */
  1.3194 +        /* check if key is first key that goes right */
  1.3195 +        if (!v->spl_nright) {
  1.3196 +          /* first key that goes right is just copied to unionR */
  1.3197 +          memcpy(unionR, key, keysize);
  1.3198 +        } else {
  1.3199 +          /* unite current value of unionR and processed key */
  1.3200 +          pgl_unite_keys(unionR, key);
  1.3201 +        }
  1.3202 +        /* append offset number to list of keys that go right */
  1.3203 +        v->spl_right[v->spl_nright++] = i;
  1.3204 +      }
  1.3205 +    }
  1.3206 +  }
  1.3207 +  /* store unions in return value */
  1.3208 +  v->spl_ldatum = PointerGetDatum(unionL);
  1.3209 +  v->spl_rdatum = PointerGetDatum(unionR);
  1.3210 +  /* return all results */
  1.3211 +  PG_RETURN_POINTER(v);
  1.3212 +}
  1.3213 +
  1.3214 +/* GiST "same"/"equal" support function */
  1.3215 +PG_FUNCTION_INFO_V1(pgl_gist_same);
  1.3216 +Datum pgl_gist_same(PG_FUNCTION_ARGS) {
  1.3217 +  pgl_keyptr key1 = (pgl_keyptr)PG_GETARG_POINTER(0);
  1.3218 +  pgl_keyptr key2 = (pgl_keyptr)PG_GETARG_POINTER(1);
  1.3219 +  bool *boolptr = (bool *)PG_GETARG_POINTER(2);
  1.3220 +  /* two keys are equal if they are binary equal */
  1.3221 +  /* (return result by setting referenced boolean and returning pointer) */
  1.3222 +  *boolptr = !memcmp(
  1.3223 +    key1,
  1.3224 +    key2,
  1.3225 +    PGL_KEY_IS_AREAKEY(key1) ? sizeof(pgl_areakey) : sizeof(pgl_pointkey)
  1.3226 +  );
  1.3227 +  PG_RETURN_POINTER(boolptr);
  1.3228 +}
  1.3229 +
  1.3230 +/* GiST "distance" support function */
  1.3231 +PG_FUNCTION_INFO_V1(pgl_gist_distance);
  1.3232 +Datum pgl_gist_distance(PG_FUNCTION_ARGS) {
  1.3233 +  GISTENTRY *entry = (GISTENTRY *)PG_GETARG_POINTER(0);
  1.3234 +  pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key);
  1.3235 +  StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2);
  1.3236 +  bool *recheck = (bool *)PG_GETARG_POINTER(4);
  1.3237 +  double distance;  /* return value */
  1.3238 +  /* demand recheck because distance is just an estimation */
  1.3239 +  /* (real distance may be bigger) */
  1.3240 +  *recheck = true;
  1.3241 +  /* strategy number aliases for different operators using the same strategy */
  1.3242 +  strategy %= 100;
  1.3243 +  /* strategy number 31: distance to point */
  1.3244 +  if (strategy == 31) {
  1.3245 +    /* query datum is a point */
  1.3246 +    pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
  1.3247 +    /* use pgl_estimate_pointkey_distance() function to compute result */
  1.3248 +    distance = pgl_estimate_key_distance(key, query);
  1.3249 +    /* avoid infinity (reserved!) */
  1.3250 +    if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
  1.3251 +    /* return result */
  1.3252 +    PG_RETURN_FLOAT8(distance);
  1.3253 +  }
  1.3254 +  /* strategy number 33: distance to circle */
  1.3255 +  if (strategy == 33) {
  1.3256 +    /* query datum is a circle */
  1.3257 +    pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
  1.3258 +    /* estimate distance to circle center and substract circle radius */
  1.3259 +    distance = (
  1.3260 +      pgl_estimate_key_distance(key, &(query->center)) - query->radius
  1.3261 +    );
  1.3262 +    /* convert non-positive values to zero and avoid infinity (reserved!) */
  1.3263 +    if (distance <= 0) distance = 0;
  1.3264 +    else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
  1.3265 +    /* return result */
  1.3266 +    PG_RETURN_FLOAT8(distance);
  1.3267 +  }
  1.3268 +  /* strategy number 34: distance to cluster */
  1.3269 +  if (strategy == 34) {
  1.3270 +    /* query datum is a cluster */
  1.3271 +    pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  1.3272 +    /* estimate distance to bounding center and substract bounding radius */
  1.3273 +    distance = (
  1.3274 +      pgl_estimate_key_distance(key, &(query->bounding.center)) -
  1.3275 +      query->bounding.radius
  1.3276 +    );
  1.3277 +    /* convert non-positive values to zero and avoid infinity (reserved!) */
  1.3278 +    if (distance <= 0) distance = 0;
  1.3279 +    else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
  1.3280 +    /* free detoasted cluster (if copy) */
  1.3281 +    PG_FREE_IF_COPY(query, 1);
  1.3282 +    /* return result */
  1.3283 +    PG_RETURN_FLOAT8(distance);
  1.3284 +  }
  1.3285 +  /* throw error for any unknown strategy number */
  1.3286 +  elog(ERROR, "unrecognized strategy number: %d", strategy);
  1.3287 +}
  1.3288 +