pgLatLon

changeset 42:1b9cd45e9e48
Bugfix for type casts to ecluster; New "fair_distance" function
author: jbe
date: Tue Oct 25 18:44:43 2016 +0200 (2016-10-25)
parents: a90c0d750405
children: 3293906ea6a0
files: GNUmakefile README.html README.mkd latlon--0.9.sql latlon-v0008.c latlon.control
     1.1 --- a/GNUmakefile	Fri Oct 21 12:57:46 2016 +0200
     1.2 +++ b/GNUmakefile	Tue Oct 25 18:44:43 2016 +0200
     1.3 @@ -1,6 +1,6 @@
     1.4  EXTENSION = latlon
     1.5 -DATA = latlon--0.7--0.8.sql latlon--0.8.sql
     1.6 -MODULES = latlon-v0007
     1.7 +DATA = latlon--0.7--0.8.sql latlon--0.8.sql latlon--0.9.sql
     1.8 +MODULES = latlon-v0007 latlon-v0008
     1.9  
    1.10  PG_CONFIG = pg_config
    1.11  PGXS := $(shell $(PG_CONFIG) --pgxs)

     2.1 --- a/README.html	Fri Oct 21 12:57:46 2016 +0200
     2.2 +++ b/README.html	Tue Oct 25 18:44:43 2016 +0200
     2.3 @@ -1,5 +1,5 @@
     2.4 -<html><head><title>pgLatLon v0.8 documentation</title></head><body>
     2.5 -<h1>pgLatLon v0.8 documentation</h1>
     2.6 +<html><head><title>pgLatLon v0.9 documentation</title></head><body>
     2.7 +<h1>pgLatLon v0.9 documentation</h1>
     2.8  
     2.9  <p>pgLatLon is a spatial database extension for the PostgreSQL object-relational
    2.10  database management system providing geographic data types and spatial indexing
    2.11 @@ -40,8 +40,8 @@
    2.12  <p>It is also possible to compile and install the extension without GNU Make as
    2.13  follows:</p>
    2.14  
    2.15 -<pre><code>cc -Wall -O2 -fPIC -shared -I `pg_config --includedir-server` -o latlon-v0007.so latlon-v0007.c
    2.16 -cp latlon-v0007.so `pg_config --pkglibdir`
    2.17 +<pre><code>cc -Wall -O2 -fPIC -shared -I `pg_config --includedir-server` -o latlon-v0008.so latlon-v0008.c
    2.18 +cp latlon-v0008.so `pg_config --pkglibdir`
    2.19  cp latlon.control `pg_config --sharedir`/extension/
    2.20  cp latlon--*.sql `pg_config --sharedir`/extension/
    2.21  </code></pre>
    2.22 @@ -486,6 +486,40 @@
    2.23  
    2.24  <p>Same as <code>epoint(float8, float8)</code> but with arguments reversed.</p>
    2.25  
    2.26 +<h4><code>fair_distance(ecluster, epoint,</code> samples <code>int4 = 10000)</code></h4>
    2.27 +
    2.28 +<p>When working with user-generated content, users may be tempted to create
    2.29 +intentionally oversized objects in order to optimize search results in an
    2.30 +unfair manner. The <code>fair_distance</code> function aims to handle this by returning an
    2.31 +adjusted distance (i.e. distance increased by a penalty) if a geographic object
    2.32 +(the <code>ecluster</code>) consists of more than one point.</p>
    2.33 +
    2.34 +<p>The first argument to this function is an <code>ecluster</code>, the second argument is a
    2.35 +search point (<code>epoint</code>), and the third argument is an interger related to the
    2.36 +precision (higher precision will require more computation time).</p>
    2.37 +
    2.38 +<p>The penalty by which the returned distance is increased fulfills (at least) the
    2.39 +following properties:</p>
    2.40 +
    2.41 +<ul>
    2.42 +<li>For search points far away from the <code>ecluster</code> (i.e. distance approaching
    2.43 +infinity), the penalty approaches zero (i.e. <code>fair_distance</code> behaves the same
    2.44 +as <code>distance</code>).</li>
    2.45 +<li>If the <code>ecluster</code> consists of a set of points, the penalty for a search point
    2.46 +close to one of that points (closer than half of the minimum distance between
    2.47 +each pair of points in the <code>ecluster</code>) is chosen in such a way that the
    2.48 +adjusted distance is equal to the distance from the search point to the
    2.49 +closest point in the <code>ecluster</code> multiplied by the square root of the count of
    2.50 +points in the <code>ecluster</code>.</li>
    2.51 +</ul>
    2.52 +
    2.53 +<p>The function interally uses a Monte Carlo simulation to compute the result. The
    2.54 +third parameter (which defaults to 10000) can be used to adjust the number of
    2.55 +samples taken. It is ensured that the penalty is always positive, i.e. results
    2.56 +returned by the <code>fair_distance</code> function are always equal to or greater than
    2.57 +the results returned by the <code>distance</code> function regardless of stochastic
    2.58 +effects.</p>
    2.59 +
    2.60  <h4><code>GeoJSON_to_epoint(jsonb, text)</code></h4>
    2.61  
    2.62  <p>Maps a GeoJSON object of type "Point" or "Feature" (which contains a

     3.1 --- a/README.mkd	Fri Oct 21 12:57:46 2016 +0200
     3.2 +++ b/README.mkd	Tue Oct 25 18:44:43 2016 +0200
     3.3 @@ -1,4 +1,4 @@
     3.4 -pgLatLon v0.8 documentation
     3.5 +pgLatLon v0.9 documentation
     3.6  ===========================
     3.7  
     3.8  pgLatLon is a spatial database extension for the PostgreSQL object-relational
     3.9 @@ -39,8 +39,8 @@
    3.10  It is also possible to compile and install the extension without GNU Make as
    3.11  follows:
    3.12  
    3.13 -    cc -Wall -O2 -fPIC -shared -I `pg_config --includedir-server` -o latlon-v0007.so latlon-v0007.c
    3.14 -    cp latlon-v0007.so `pg_config --pkglibdir`
    3.15 +    cc -Wall -O2 -fPIC -shared -I `pg_config --includedir-server` -o latlon-v0008.so latlon-v0008.c
    3.16 +    cp latlon-v0008.so `pg_config --pkglibdir`
    3.17      cp latlon.control `pg_config --sharedir`/extension/
    3.18      cp latlon--*.sql `pg_config --sharedir`/extension/
    3.19  
    3.20 @@ -469,6 +469,38 @@
    3.21  
    3.22  Same as `epoint(float8, float8)` but with arguments reversed.
    3.23  
    3.24 +#### `fair_distance(ecluster, epoint,` samples `int4 = 10000)`
    3.25 +
    3.26 +When working with user-generated content, users may be tempted to create
    3.27 +intentionally oversized objects in order to optimize search results in an
    3.28 +unfair manner. The `fair_distance` function aims to handle this by returning an
    3.29 +adjusted distance (i.e. distance increased by a penalty) if a geographic object
    3.30 +(the `ecluster`) consists of more than one point.
    3.31 +
    3.32 +The first argument to this function is an `ecluster`, the second argument is a
    3.33 +search point (`epoint`), and the third argument is an interger related to the
    3.34 +precision (higher precision will require more computation time).
    3.35 +
    3.36 +The penalty by which the returned distance is increased fulfills (at least) the
    3.37 +following properties:
    3.38 +
    3.39 +* For search points far away from the `ecluster` (i.e. distance approaching
    3.40 +  infinity), the penalty approaches zero (i.e. `fair_distance` behaves the same
    3.41 +  as `distance`).
    3.42 +* If the `ecluster` consists of a set of points, the penalty for a search point
    3.43 +  close to one of that points (closer than half of the minimum distance between
    3.44 +  each pair of points in the `ecluster`) is chosen in such a way that the
    3.45 +  adjusted distance is equal to the distance from the search point to the
    3.46 +  closest point in the `ecluster` multiplied by the square root of the count of
    3.47 +  points in the `ecluster`.
    3.48 +
    3.49 +The function interally uses a Monte Carlo simulation to compute the result. The
    3.50 +third parameter (which defaults to 10000) can be used to adjust the number of
    3.51 +samples taken. It is ensured that the penalty is always positive, i.e. results
    3.52 +returned by the `fair_distance` function are always equal to or greater than
    3.53 +the results returned by the `distance` function regardless of stochastic
    3.54 +effects.
    3.55 +
    3.56  #### `GeoJSON_to_epoint(jsonb, text)`
    3.57  
    3.58  Maps a GeoJSON object of type "Point" or "Feature" (which contains a

     4.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     4.2 +++ b/latlon--0.9.sql	Tue Oct 25 18:44:43 2016 +0200
     4.3 @@ -0,0 +1,1692 @@
     4.4 +
     4.5 +----------------------------------------
     4.6 +-- forward declarations (shell types) --
     4.7 +----------------------------------------
     4.8 +
     4.9 +CREATE TYPE epoint;
    4.10 +CREATE TYPE ebox;
    4.11 +CREATE TYPE ecircle;
    4.12 +CREATE TYPE ecluster;
    4.13 +
    4.14 +
    4.15 +------------------------------------------------------------
    4.16 +-- dummy input/output functions for dummy index key types --
    4.17 +------------------------------------------------------------
    4.18 +
    4.19 +CREATE FUNCTION ekey_point_in_dummy(cstring)
    4.20 +  RETURNS ekey_point
    4.21 +  LANGUAGE C IMMUTABLE STRICT
    4.22 +  AS '$libdir/latlon-v0008', 'pgl_notimpl';
    4.23 +
    4.24 +CREATE FUNCTION ekey_point_out_dummy(ekey_point)
    4.25 +  RETURNS cstring
    4.26 +  LANGUAGE C IMMUTABLE STRICT
    4.27 +  AS '$libdir/latlon-v0008', 'pgl_notimpl';
    4.28 +
    4.29 +CREATE FUNCTION ekey_area_in_dummy(cstring)
    4.30 +  RETURNS ekey_area
    4.31 +  LANGUAGE C IMMUTABLE STRICT
    4.32 +  AS '$libdir/latlon-v0008', 'pgl_notimpl';
    4.33 +
    4.34 +CREATE FUNCTION ekey_area_out_dummy(ekey_area)
    4.35 +  RETURNS cstring
    4.36 +  LANGUAGE C IMMUTABLE STRICT
    4.37 +  AS '$libdir/latlon-v0008', 'pgl_notimpl';
    4.38 +
    4.39 +
    4.40 +--------------------------
    4.41 +-- text input functions --
    4.42 +--------------------------
    4.43 +
    4.44 +CREATE FUNCTION epoint_in(cstring)
    4.45 +  RETURNS epoint
    4.46 +  LANGUAGE C IMMUTABLE STRICT
    4.47 +  AS '$libdir/latlon-v0008', 'pgl_epoint_in';
    4.48 +
    4.49 +CREATE FUNCTION ebox_in(cstring)
    4.50 +  RETURNS ebox
    4.51 +  LANGUAGE C IMMUTABLE STRICT
    4.52 +  AS '$libdir/latlon-v0008', 'pgl_ebox_in';
    4.53 +
    4.54 +CREATE FUNCTION ecircle_in(cstring)
    4.55 +  RETURNS ecircle
    4.56 +  LANGUAGE C IMMUTABLE STRICT
    4.57 +  AS '$libdir/latlon-v0008', 'pgl_ecircle_in';
    4.58 +
    4.59 +CREATE FUNCTION ecluster_in(cstring)
    4.60 +  RETURNS ecluster
    4.61 +  LANGUAGE C IMMUTABLE STRICT
    4.62 +  AS '$libdir/latlon-v0008', 'pgl_ecluster_in';
    4.63 +
    4.64 +
    4.65 +---------------------------
    4.66 +-- text output functions --
    4.67 +---------------------------
    4.68 +
    4.69 +CREATE FUNCTION epoint_out(epoint)
    4.70 +  RETURNS cstring
    4.71 +  LANGUAGE C IMMUTABLE STRICT
    4.72 +  AS '$libdir/latlon-v0008', 'pgl_epoint_out';
    4.73 +
    4.74 +CREATE FUNCTION ebox_out(ebox)
    4.75 +  RETURNS cstring
    4.76 +  LANGUAGE C IMMUTABLE STRICT
    4.77 +  AS '$libdir/latlon-v0008', 'pgl_ebox_out';
    4.78 +
    4.79 +CREATE FUNCTION ecircle_out(ecircle)
    4.80 +  RETURNS cstring
    4.81 +  LANGUAGE C IMMUTABLE STRICT
    4.82 +  AS '$libdir/latlon-v0008', 'pgl_ecircle_out';
    4.83 +
    4.84 +CREATE FUNCTION ecluster_out(ecluster)
    4.85 +  RETURNS cstring
    4.86 +  LANGUAGE C IMMUTABLE STRICT
    4.87 +  AS '$libdir/latlon-v0008', 'pgl_ecluster_out';
    4.88 +
    4.89 +
    4.90 +--------------------------
    4.91 +-- binary I/O functions --
    4.92 +--------------------------
    4.93 +
    4.94 +CREATE FUNCTION epoint_recv(internal)
    4.95 +  RETURNS epoint
    4.96 +  LANGUAGE C IMMUTABLE STRICT
    4.97 +  AS '$libdir/latlon-v0008', 'pgl_epoint_recv';
    4.98 +
    4.99 +CREATE FUNCTION ebox_recv(internal)
   4.100 +  RETURNS ebox
   4.101 +  LANGUAGE C IMMUTABLE STRICT
   4.102 +  AS '$libdir/latlon-v0008', 'pgl_ebox_recv';
   4.103 +
   4.104 +CREATE FUNCTION ecircle_recv(internal)
   4.105 +  RETURNS ecircle
   4.106 +  LANGUAGE C IMMUTABLE STRICT
   4.107 +  AS '$libdir/latlon-v0008', 'pgl_ecircle_recv';
   4.108 +
   4.109 +CREATE FUNCTION epoint_send(epoint)
   4.110 +  RETURNS bytea
   4.111 +  LANGUAGE C IMMUTABLE STRICT
   4.112 +  AS '$libdir/latlon-v0008', 'pgl_epoint_send';
   4.113 +
   4.114 +CREATE FUNCTION ebox_send(ebox)
   4.115 +  RETURNS bytea
   4.116 +  LANGUAGE C IMMUTABLE STRICT
   4.117 +  AS '$libdir/latlon-v0008', 'pgl_ebox_send';
   4.118 +
   4.119 +CREATE FUNCTION ecircle_send(ecircle)
   4.120 +  RETURNS bytea
   4.121 +  LANGUAGE C IMMUTABLE STRICT
   4.122 +  AS '$libdir/latlon-v0008', 'pgl_ecircle_send';
   4.123 +
   4.124 +
   4.125 +-----------------------------------------------
   4.126 +-- type definitions of dummy index key types --
   4.127 +-----------------------------------------------
   4.128 +
   4.129 +CREATE TYPE ekey_point (
   4.130 +  internallength = 8,
   4.131 +  input = ekey_point_in_dummy,
   4.132 +  output = ekey_point_out_dummy,
   4.133 +  alignment = char );
   4.134 +
   4.135 +CREATE TYPE ekey_area (
   4.136 +  internallength = 9,
   4.137 +  input = ekey_area_in_dummy,
   4.138 +  output = ekey_area_out_dummy,
   4.139 +  alignment = char );
   4.140 +
   4.141 +
   4.142 +------------------------------------------
   4.143 +-- definitions of geographic data types --
   4.144 +------------------------------------------
   4.145 +
   4.146 +CREATE TYPE epoint (
   4.147 +  internallength = 16,
   4.148 +  input = epoint_in,
   4.149 +  output = epoint_out,
   4.150 +  receive = epoint_recv,
   4.151 +  send = epoint_send,
   4.152 +  alignment = double );
   4.153 +
   4.154 +CREATE TYPE ebox (
   4.155 +  internallength = 32,
   4.156 +  input = ebox_in,
   4.157 +  output = ebox_out,
   4.158 +  receive = ebox_recv,
   4.159 +  send = ebox_send,
   4.160 +  alignment = double );
   4.161 +
   4.162 +CREATE TYPE ecircle (
   4.163 +  internallength = 24,
   4.164 +  input = ecircle_in,
   4.165 +  output = ecircle_out,
   4.166 +  receive = ecircle_recv,
   4.167 +  send = ecircle_send,
   4.168 +  alignment = double );
   4.169 +
   4.170 +CREATE TYPE ecluster (
   4.171 +  internallength = VARIABLE,
   4.172 +  input = ecluster_in,
   4.173 +  output = ecluster_out,
   4.174 +  alignment = double,
   4.175 +  storage = external );
   4.176 +
   4.177 +
   4.178 +--------------------
   4.179 +-- B-tree support --
   4.180 +--------------------
   4.181 +
   4.182 +-- begin of B-tree support for epoint
   4.183 +
   4.184 +CREATE FUNCTION epoint_btree_lt(epoint, epoint)
   4.185 +  RETURNS boolean
   4.186 +  LANGUAGE C IMMUTABLE STRICT
   4.187 +  AS '$libdir/latlon-v0008', 'pgl_btree_epoint_lt';
   4.188 +
   4.189 +CREATE FUNCTION epoint_btree_le(epoint, epoint)
   4.190 +  RETURNS boolean
   4.191 +  LANGUAGE C IMMUTABLE STRICT
   4.192 +  AS '$libdir/latlon-v0008', 'pgl_btree_epoint_le';
   4.193 +
   4.194 +CREATE FUNCTION epoint_btree_eq(epoint, epoint)
   4.195 +  RETURNS boolean
   4.196 +  LANGUAGE C IMMUTABLE STRICT
   4.197 +  AS '$libdir/latlon-v0008', 'pgl_btree_epoint_eq';
   4.198 +
   4.199 +CREATE FUNCTION epoint_btree_ne(epoint, epoint)
   4.200 +  RETURNS boolean
   4.201 +  LANGUAGE C IMMUTABLE STRICT
   4.202 +  AS '$libdir/latlon-v0008', 'pgl_btree_epoint_ne';
   4.203 +
   4.204 +CREATE FUNCTION epoint_btree_ge(epoint, epoint)
   4.205 +  RETURNS boolean
   4.206 +  LANGUAGE C IMMUTABLE STRICT
   4.207 +  AS '$libdir/latlon-v0008', 'pgl_btree_epoint_ge';
   4.208 +
   4.209 +CREATE FUNCTION epoint_btree_gt(epoint, epoint)
   4.210 +  RETURNS boolean
   4.211 +  LANGUAGE C IMMUTABLE STRICT
   4.212 +  AS '$libdir/latlon-v0008', 'pgl_btree_epoint_gt';
   4.213 +
   4.214 +CREATE OPERATOR <<< (
   4.215 +  leftarg = epoint,
   4.216 +  rightarg = epoint,
   4.217 +  procedure = epoint_btree_lt,
   4.218 +  commutator = >>>,
   4.219 +  negator = >>>=,
   4.220 +  restrict = scalarltsel,
   4.221 +  join = scalarltjoinsel
   4.222 +);
   4.223 +
   4.224 +CREATE OPERATOR <<<= (
   4.225 +  leftarg = epoint,
   4.226 +  rightarg = epoint,
   4.227 +  procedure = epoint_btree_le,
   4.228 +  commutator = >>>=,
   4.229 +  negator = >>>,
   4.230 +  restrict = scalarltsel,
   4.231 +  join = scalarltjoinsel
   4.232 +);
   4.233 +
   4.234 +CREATE OPERATOR = (
   4.235 +  leftarg = epoint,
   4.236 +  rightarg = epoint,
   4.237 +  procedure = epoint_btree_eq,
   4.238 +  commutator = =,
   4.239 +  negator = <>,
   4.240 +  restrict = eqsel,
   4.241 +  join = eqjoinsel,
   4.242 +  merges
   4.243 +);
   4.244 +
   4.245 +CREATE OPERATOR <> (
   4.246 +  leftarg = epoint,
   4.247 +  rightarg = epoint,
   4.248 +  procedure = epoint_btree_eq,
   4.249 +  commutator = <>,
   4.250 +  negator = =,
   4.251 +  restrict = neqsel,
   4.252 +  join = neqjoinsel
   4.253 +);
   4.254 +
   4.255 +CREATE OPERATOR >>>= (
   4.256 +  leftarg = epoint,
   4.257 +  rightarg = epoint,
   4.258 +  procedure = epoint_btree_ge,
   4.259 +  commutator = <<<=,
   4.260 +  negator = <<<,
   4.261 +  restrict = scalargtsel,
   4.262 +  join = scalargtjoinsel
   4.263 +);
   4.264 +
   4.265 +CREATE OPERATOR >>> (
   4.266 +  leftarg = epoint,
   4.267 +  rightarg = epoint,
   4.268 +  procedure = epoint_btree_gt,
   4.269 +  commutator = <<<,
   4.270 +  negator = <<<=,
   4.271 +  restrict = scalargtsel,
   4.272 +  join = scalargtjoinsel
   4.273 +);
   4.274 +
   4.275 +CREATE FUNCTION epoint_btree_cmp(epoint, epoint)
   4.276 +  RETURNS int4
   4.277 +  LANGUAGE C IMMUTABLE STRICT
   4.278 +  AS '$libdir/latlon-v0008', 'pgl_btree_epoint_cmp';
   4.279 +
   4.280 +CREATE OPERATOR CLASS epoint_btree_ops
   4.281 +  DEFAULT FOR TYPE epoint USING btree AS
   4.282 +  OPERATOR 1 <<< ,
   4.283 +  OPERATOR 2 <<<= ,
   4.284 +  OPERATOR 3 = ,
   4.285 +  OPERATOR 4 >>>= ,
   4.286 +  OPERATOR 5 >>> ,
   4.287 +  FUNCTION 1 epoint_btree_cmp(epoint, epoint);
   4.288 +
   4.289 +-- end of B-tree support for epoint
   4.290 +
   4.291 +-- begin of B-tree support for ebox
   4.292 +
   4.293 +CREATE FUNCTION ebox_btree_lt(ebox, ebox)
   4.294 +  RETURNS boolean
   4.295 +  LANGUAGE C IMMUTABLE STRICT
   4.296 +  AS '$libdir/latlon-v0008', 'pgl_btree_ebox_lt';
   4.297 +
   4.298 +CREATE FUNCTION ebox_btree_le(ebox, ebox)
   4.299 +  RETURNS boolean
   4.300 +  LANGUAGE C IMMUTABLE STRICT
   4.301 +  AS '$libdir/latlon-v0008', 'pgl_btree_ebox_le';
   4.302 +
   4.303 +CREATE FUNCTION ebox_btree_eq(ebox, ebox)
   4.304 +  RETURNS boolean
   4.305 +  LANGUAGE C IMMUTABLE STRICT
   4.306 +  AS '$libdir/latlon-v0008', 'pgl_btree_ebox_eq';
   4.307 +
   4.308 +CREATE FUNCTION ebox_btree_ne(ebox, ebox)
   4.309 +  RETURNS boolean
   4.310 +  LANGUAGE C IMMUTABLE STRICT
   4.311 +  AS '$libdir/latlon-v0008', 'pgl_btree_ebox_ne';
   4.312 +
   4.313 +CREATE FUNCTION ebox_btree_ge(ebox, ebox)
   4.314 +  RETURNS boolean
   4.315 +  LANGUAGE C IMMUTABLE STRICT
   4.316 +  AS '$libdir/latlon-v0008', 'pgl_btree_ebox_ge';
   4.317 +
   4.318 +CREATE FUNCTION ebox_btree_gt(ebox, ebox)
   4.319 +  RETURNS boolean
   4.320 +  LANGUAGE C IMMUTABLE STRICT
   4.321 +  AS '$libdir/latlon-v0008', 'pgl_btree_ebox_gt';
   4.322 +
   4.323 +CREATE OPERATOR <<< (
   4.324 +  leftarg = ebox,
   4.325 +  rightarg = ebox,
   4.326 +  procedure = ebox_btree_lt,
   4.327 +  commutator = >>>,
   4.328 +  negator = >>>=,
   4.329 +  restrict = scalarltsel,
   4.330 +  join = scalarltjoinsel
   4.331 +);
   4.332 +
   4.333 +CREATE OPERATOR <<<= (
   4.334 +  leftarg = ebox,
   4.335 +  rightarg = ebox,
   4.336 +  procedure = ebox_btree_le,
   4.337 +  commutator = >>>=,
   4.338 +  negator = >>>,
   4.339 +  restrict = scalarltsel,
   4.340 +  join = scalarltjoinsel
   4.341 +);
   4.342 +
   4.343 +CREATE OPERATOR = (
   4.344 +  leftarg = ebox,
   4.345 +  rightarg = ebox,
   4.346 +  procedure = ebox_btree_eq,
   4.347 +  commutator = =,
   4.348 +  negator = <>,
   4.349 +  restrict = eqsel,
   4.350 +  join = eqjoinsel,
   4.351 +  merges
   4.352 +);
   4.353 +
   4.354 +CREATE OPERATOR <> (
   4.355 +  leftarg = ebox,
   4.356 +  rightarg = ebox,
   4.357 +  procedure = ebox_btree_eq,
   4.358 +  commutator = <>,
   4.359 +  negator = =,
   4.360 +  restrict = neqsel,
   4.361 +  join = neqjoinsel
   4.362 +);
   4.363 +
   4.364 +CREATE OPERATOR >>>= (
   4.365 +  leftarg = ebox,
   4.366 +  rightarg = ebox,
   4.367 +  procedure = ebox_btree_ge,
   4.368 +  commutator = <<<=,
   4.369 +  negator = <<<,
   4.370 +  restrict = scalargtsel,
   4.371 +  join = scalargtjoinsel
   4.372 +);
   4.373 +
   4.374 +CREATE OPERATOR >>> (
   4.375 +  leftarg = ebox,
   4.376 +  rightarg = ebox,
   4.377 +  procedure = ebox_btree_gt,
   4.378 +  commutator = <<<,
   4.379 +  negator = <<<=,
   4.380 +  restrict = scalargtsel,
   4.381 +  join = scalargtjoinsel
   4.382 +);
   4.383 +
   4.384 +CREATE FUNCTION ebox_btree_cmp(ebox, ebox)
   4.385 +  RETURNS int4
   4.386 +  LANGUAGE C IMMUTABLE STRICT
   4.387 +  AS '$libdir/latlon-v0008', 'pgl_btree_ebox_cmp';
   4.388 +
   4.389 +CREATE OPERATOR CLASS ebox_btree_ops
   4.390 +  DEFAULT FOR TYPE ebox USING btree AS
   4.391 +  OPERATOR 1 <<< ,
   4.392 +  OPERATOR 2 <<<= ,
   4.393 +  OPERATOR 3 = ,
   4.394 +  OPERATOR 4 >>>= ,
   4.395 +  OPERATOR 5 >>> ,
   4.396 +  FUNCTION 1 ebox_btree_cmp(ebox, ebox);
   4.397 +
   4.398 +-- end of B-tree support for ebox
   4.399 +
   4.400 +-- begin of B-tree support for ecircle
   4.401 +
   4.402 +CREATE FUNCTION ecircle_btree_lt(ecircle, ecircle)
   4.403 +  RETURNS boolean
   4.404 +  LANGUAGE C IMMUTABLE STRICT
   4.405 +  AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_lt';
   4.406 +
   4.407 +CREATE FUNCTION ecircle_btree_le(ecircle, ecircle)
   4.408 +  RETURNS boolean
   4.409 +  LANGUAGE C IMMUTABLE STRICT
   4.410 +  AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_le';
   4.411 +
   4.412 +CREATE FUNCTION ecircle_btree_eq(ecircle, ecircle)
   4.413 +  RETURNS boolean
   4.414 +  LANGUAGE C IMMUTABLE STRICT
   4.415 +  AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_eq';
   4.416 +
   4.417 +CREATE FUNCTION ecircle_btree_ne(ecircle, ecircle)
   4.418 +  RETURNS boolean
   4.419 +  LANGUAGE C IMMUTABLE STRICT
   4.420 +  AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_ne';
   4.421 +
   4.422 +CREATE FUNCTION ecircle_btree_ge(ecircle, ecircle)
   4.423 +  RETURNS boolean
   4.424 +  LANGUAGE C IMMUTABLE STRICT
   4.425 +  AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_ge';
   4.426 +
   4.427 +CREATE FUNCTION ecircle_btree_gt(ecircle, ecircle)
   4.428 +  RETURNS boolean
   4.429 +  LANGUAGE C IMMUTABLE STRICT
   4.430 +  AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_gt';
   4.431 +
   4.432 +CREATE OPERATOR <<< (
   4.433 +  leftarg = ecircle,
   4.434 +  rightarg = ecircle,
   4.435 +  procedure = ecircle_btree_lt,
   4.436 +  commutator = >>>,
   4.437 +  negator = >>>=,
   4.438 +  restrict = scalarltsel,
   4.439 +  join = scalarltjoinsel
   4.440 +);
   4.441 +
   4.442 +CREATE OPERATOR <<<= (
   4.443 +  leftarg = ecircle,
   4.444 +  rightarg = ecircle,
   4.445 +  procedure = ecircle_btree_le,
   4.446 +  commutator = >>>=,
   4.447 +  negator = >>>,
   4.448 +  restrict = scalarltsel,
   4.449 +  join = scalarltjoinsel
   4.450 +);
   4.451 +
   4.452 +CREATE OPERATOR = (
   4.453 +  leftarg = ecircle,
   4.454 +  rightarg = ecircle,
   4.455 +  procedure = ecircle_btree_eq,
   4.456 +  commutator = =,
   4.457 +  negator = <>,
   4.458 +  restrict = eqsel,
   4.459 +  join = eqjoinsel,
   4.460 +  merges
   4.461 +);
   4.462 +
   4.463 +CREATE OPERATOR <> (
   4.464 +  leftarg = ecircle,
   4.465 +  rightarg = ecircle,
   4.466 +  procedure = ecircle_btree_eq,
   4.467 +  commutator = <>,
   4.468 +  negator = =,
   4.469 +  restrict = neqsel,
   4.470 +  join = neqjoinsel
   4.471 +);
   4.472 +
   4.473 +CREATE OPERATOR >>>= (
   4.474 +  leftarg = ecircle,
   4.475 +  rightarg = ecircle,
   4.476 +  procedure = ecircle_btree_ge,
   4.477 +  commutator = <<<=,
   4.478 +  negator = <<<,
   4.479 +  restrict = scalargtsel,
   4.480 +  join = scalargtjoinsel
   4.481 +);
   4.482 +
   4.483 +CREATE OPERATOR >>> (
   4.484 +  leftarg = ecircle,
   4.485 +  rightarg = ecircle,
   4.486 +  procedure = ecircle_btree_gt,
   4.487 +  commutator = <<<,
   4.488 +  negator = <<<=,
   4.489 +  restrict = scalargtsel,
   4.490 +  join = scalargtjoinsel
   4.491 +);
   4.492 +
   4.493 +CREATE FUNCTION ecircle_btree_cmp(ecircle, ecircle)
   4.494 +  RETURNS int4
   4.495 +  LANGUAGE C IMMUTABLE STRICT
   4.496 +  AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_cmp';
   4.497 +
   4.498 +CREATE OPERATOR CLASS ecircle_btree_ops
   4.499 +  DEFAULT FOR TYPE ecircle USING btree AS
   4.500 +  OPERATOR 1 <<< ,
   4.501 +  OPERATOR 2 <<<= ,
   4.502 +  OPERATOR 3 = ,
   4.503 +  OPERATOR 4 >>>= ,
   4.504 +  OPERATOR 5 >>> ,
   4.505 +  FUNCTION 1 ecircle_btree_cmp(ecircle, ecircle);
   4.506 +
   4.507 +-- end of B-tree support for ecircle
   4.508 +
   4.509 +
   4.510 +----------------
   4.511 +-- type casts --
   4.512 +----------------
   4.513 +
   4.514 +CREATE FUNCTION cast_epoint_to_ebox(epoint)
   4.515 +  RETURNS ebox
   4.516 +  LANGUAGE C IMMUTABLE STRICT
   4.517 +  AS '$libdir/latlon-v0008', 'pgl_epoint_to_ebox';
   4.518 +
   4.519 +CREATE CAST (epoint AS ebox) WITH FUNCTION cast_epoint_to_ebox(epoint);
   4.520 +
   4.521 +CREATE FUNCTION cast_epoint_to_ecircle(epoint)
   4.522 +  RETURNS ecircle
   4.523 +  LANGUAGE C IMMUTABLE STRICT
   4.524 +  AS '$libdir/latlon-v0008', 'pgl_epoint_to_ecircle';
   4.525 +
   4.526 +CREATE CAST (epoint AS ecircle) WITH FUNCTION cast_epoint_to_ecircle(epoint);
   4.527 +
   4.528 +CREATE FUNCTION cast_epoint_to_ecluster(epoint)
   4.529 +  RETURNS ecluster
   4.530 +  LANGUAGE C IMMUTABLE STRICT
   4.531 +  AS '$libdir/latlon-v0008', 'pgl_epoint_to_ecluster';
   4.532 +
   4.533 +CREATE CAST (epoint AS ecluster) WITH FUNCTION cast_epoint_to_ecluster(epoint);
   4.534 +
   4.535 +CREATE FUNCTION cast_ebox_to_ecluster(ebox)
   4.536 +  RETURNS ecluster
   4.537 +  LANGUAGE C IMMUTABLE STRICT
   4.538 +  AS '$libdir/latlon-v0008', 'pgl_ebox_to_ecluster';
   4.539 +
   4.540 +CREATE CAST (ebox AS ecluster) WITH FUNCTION cast_ebox_to_ecluster(ebox);
   4.541 +
   4.542 +
   4.543 +---------------------------
   4.544 +-- constructor functions --
   4.545 +---------------------------
   4.546 +
   4.547 +CREATE FUNCTION epoint(float8, float8)
   4.548 +  RETURNS epoint
   4.549 +  LANGUAGE C IMMUTABLE STRICT
   4.550 +  AS '$libdir/latlon-v0008', 'pgl_create_epoint';
   4.551 +
   4.552 +CREATE FUNCTION epoint_latlon(float8, float8)
   4.553 +  RETURNS epoint
   4.554 +  LANGUAGE SQL IMMUTABLE STRICT AS $$
   4.555 +    SELECT epoint($1, $2)
   4.556 +  $$;
   4.557 +
   4.558 +CREATE FUNCTION epoint_lonlat(float8, float8)
   4.559 +  RETURNS epoint
   4.560 +  LANGUAGE SQL IMMUTABLE STRICT AS $$
   4.561 +    SELECT epoint($2, $1)
   4.562 +  $$;
   4.563 +
   4.564 +CREATE FUNCTION empty_ebox()
   4.565 +  RETURNS ebox
   4.566 +  LANGUAGE C IMMUTABLE STRICT
   4.567 +  AS '$libdir/latlon-v0008', 'pgl_create_empty_ebox';
   4.568 +
   4.569 +CREATE FUNCTION ebox(float8, float8, float8, float8)
   4.570 +  RETURNS ebox
   4.571 +  LANGUAGE C IMMUTABLE STRICT
   4.572 +  AS '$libdir/latlon-v0008', 'pgl_create_ebox';
   4.573 +
   4.574 +CREATE FUNCTION ebox(epoint, epoint)
   4.575 +  RETURNS ebox
   4.576 +  LANGUAGE C IMMUTABLE STRICT
   4.577 +  AS '$libdir/latlon-v0008', 'pgl_create_ebox_from_epoints';
   4.578 +
   4.579 +CREATE FUNCTION ecircle(float8, float8, float8)
   4.580 +  RETURNS ecircle
   4.581 +  LANGUAGE C IMMUTABLE STRICT
   4.582 +  AS '$libdir/latlon-v0008', 'pgl_create_ecircle';
   4.583 +
   4.584 +CREATE FUNCTION ecircle(epoint, float8)
   4.585 +  RETURNS ecircle
   4.586 +  LANGUAGE C IMMUTABLE STRICT
   4.587 +  AS '$libdir/latlon-v0008', 'pgl_create_ecircle_from_epoint';
   4.588 +
   4.589 +CREATE FUNCTION ecluster_concat(ecluster[])
   4.590 +  RETURNS ecluster
   4.591 +  LANGUAGE sql IMMUTABLE STRICT AS $$
   4.592 +    SELECT array_to_string($1, ' ')::ecluster
   4.593 +  $$;
   4.594 +
   4.595 +CREATE FUNCTION ecluster_concat(ecluster, ecluster)
   4.596 +  RETURNS ecluster
   4.597 +  LANGUAGE sql IMMUTABLE STRICT AS $$
   4.598 +    SELECT ($1::text || ' ' || $2::text)::ecluster
   4.599 +  $$;
   4.600 +
   4.601 +CREATE FUNCTION ecluster_create_multipoint(epoint[])
   4.602 +  RETURNS ecluster
   4.603 +  LANGUAGE sql IMMUTABLE STRICT AS $$
   4.604 +    SELECT
   4.605 +      array_to_string(array_agg('point (' || unnest || ')'), ' ')::ecluster
   4.606 +    FROM unnest($1)
   4.607 +  $$;
   4.608 +
   4.609 +CREATE FUNCTION ecluster_create_path(epoint[])
   4.610 +  RETURNS ecluster
   4.611 +  LANGUAGE sql IMMUTABLE STRICT AS $$
   4.612 +    SELECT CASE WHEN "str" = '' THEN 'empty'::ecluster ELSE
   4.613 +      ('path (' || array_to_string($1, ' ') || ')')::ecluster
   4.614 +    END
   4.615 +    FROM array_to_string($1, ' ') AS "str"
   4.616 +  $$;
   4.617 +
   4.618 +CREATE FUNCTION ecluster_create_outline(epoint[])
   4.619 +  RETURNS ecluster
   4.620 +  LANGUAGE sql IMMUTABLE STRICT AS $$
   4.621 +    SELECT CASE WHEN "str" = '' THEN 'empty'::ecluster ELSE
   4.622 +      ('outline (' || array_to_string($1, ' ') || ')')::ecluster
   4.623 +    END
   4.624 +    FROM array_to_string($1, ' ') AS "str"
   4.625 +  $$;
   4.626 +
   4.627 +CREATE FUNCTION ecluster_create_polygon(epoint[])
   4.628 +  RETURNS ecluster
   4.629 +  LANGUAGE sql IMMUTABLE STRICT AS $$
   4.630 +    SELECT CASE WHEN "str" = '' THEN 'empty'::ecluster ELSE
   4.631 +      ('polygon (' || array_to_string($1, ' ') || ')')::ecluster
   4.632 +    END
   4.633 +    FROM array_to_string($1, ' ') AS "str"
   4.634 +  $$;
   4.635 +
   4.636 +
   4.637 +----------------------
   4.638 +-- getter functions --
   4.639 +----------------------
   4.640 +
   4.641 +CREATE FUNCTION latitude(epoint)
   4.642 +  RETURNS float8
   4.643 +  LANGUAGE C IMMUTABLE STRICT
   4.644 +  AS '$libdir/latlon-v0008', 'pgl_epoint_lat';
   4.645 +
   4.646 +CREATE FUNCTION longitude(epoint)
   4.647 +  RETURNS float8
   4.648 +  LANGUAGE C IMMUTABLE STRICT
   4.649 +  AS '$libdir/latlon-v0008', 'pgl_epoint_lon';
   4.650 +
   4.651 +CREATE FUNCTION min_latitude(ebox)
   4.652 +  RETURNS float8
   4.653 +  LANGUAGE C IMMUTABLE STRICT
   4.654 +  AS '$libdir/latlon-v0008', 'pgl_ebox_lat_min';
   4.655 +
   4.656 +CREATE FUNCTION max_latitude(ebox)
   4.657 +  RETURNS float8
   4.658 +  LANGUAGE C IMMUTABLE STRICT
   4.659 +  AS '$libdir/latlon-v0008', 'pgl_ebox_lat_max';
   4.660 +
   4.661 +CREATE FUNCTION min_longitude(ebox)
   4.662 +  RETURNS float8
   4.663 +  LANGUAGE C IMMUTABLE STRICT
   4.664 +  AS '$libdir/latlon-v0008', 'pgl_ebox_lon_min';
   4.665 +
   4.666 +CREATE FUNCTION max_longitude(ebox)
   4.667 +  RETURNS float8
   4.668 +  LANGUAGE C IMMUTABLE STRICT
   4.669 +  AS '$libdir/latlon-v0008', 'pgl_ebox_lon_max';
   4.670 +
   4.671 +CREATE FUNCTION center(ecircle)
   4.672 +  RETURNS epoint
   4.673 +  LANGUAGE C IMMUTABLE STRICT
   4.674 +  AS '$libdir/latlon-v0008', 'pgl_ecircle_center';
   4.675 +
   4.676 +CREATE FUNCTION radius(ecircle)
   4.677 +  RETURNS float8
   4.678 +  LANGUAGE C IMMUTABLE STRICT
   4.679 +  AS '$libdir/latlon-v0008', 'pgl_ecircle_radius';
   4.680 +
   4.681 +CREATE FUNCTION ecluster_extract_points(ecluster)
   4.682 +  RETURNS SETOF epoint
   4.683 +  LANGUAGE sql IMMUTABLE STRICT AS $$
   4.684 +    SELECT "match"[2]::epoint
   4.685 +    FROM regexp_matches($1::text, e'(^| )point \\(([^)]+)\\)', 'g') AS "match"
   4.686 +  $$;
   4.687 +
   4.688 +CREATE FUNCTION ecluster_extract_paths(ecluster)
   4.689 +  RETURNS SETOF epoint[]
   4.690 +  LANGUAGE sql IMMUTABLE STRICT AS $$
   4.691 +    SELECT (
   4.692 +      SELECT array_agg("m2"[1]::epoint)
   4.693 +      FROM regexp_matches("m1"[2], e'[^ ]+ [^ ]+', 'g') AS "m2"
   4.694 +    )
   4.695 +    FROM regexp_matches($1::text, e'(^| )path \\(([^)]+)\\)', 'g') AS "m1"
   4.696 +  $$;
   4.697 +
   4.698 +CREATE FUNCTION ecluster_extract_outlines(ecluster)
   4.699 +  RETURNS SETOF epoint[]
   4.700 +  LANGUAGE sql IMMUTABLE STRICT AS $$
   4.701 +    SELECT (
   4.702 +      SELECT array_agg("m2"[1]::epoint)
   4.703 +      FROM regexp_matches("m1"[2], e'[^ ]+ [^ ]+', 'g') AS "m2"
   4.704 +    )
   4.705 +    FROM regexp_matches($1::text, e'(^| )outline \\(([^)]+)\\)', 'g') AS "m1"
   4.706 +  $$;
   4.707 +
   4.708 +CREATE FUNCTION ecluster_extract_polygons(ecluster)
   4.709 +  RETURNS SETOF epoint[]
   4.710 +  LANGUAGE sql IMMUTABLE STRICT AS $$
   4.711 +    SELECT (
   4.712 +      SELECT array_agg("m2"[1]::epoint)
   4.713 +      FROM regexp_matches("m1"[2], e'[^ ]+ [^ ]+', 'g') AS "m2"
   4.714 +    )
   4.715 +    FROM regexp_matches($1::text, e'(^| )polygon \\(([^)]+)\\)', 'g') AS "m1"
   4.716 +  $$;
   4.717 +
   4.718 +
   4.719 +---------------
   4.720 +-- operators --
   4.721 +---------------
   4.722 +
   4.723 +CREATE FUNCTION epoint_ebox_overlap_proc(epoint, ebox)
   4.724 +  RETURNS boolean
   4.725 +  LANGUAGE C IMMUTABLE STRICT
   4.726 +  AS '$libdir/latlon-v0008', 'pgl_epoint_ebox_overlap';
   4.727 +
   4.728 +CREATE FUNCTION epoint_ecircle_overlap_proc(epoint, ecircle)
   4.729 +  RETURNS boolean
   4.730 +  LANGUAGE C IMMUTABLE STRICT
   4.731 +  AS '$libdir/latlon-v0008', 'pgl_epoint_ecircle_overlap';
   4.732 +
   4.733 +CREATE FUNCTION epoint_ecluster_overlap_proc(epoint, ecluster)
   4.734 +  RETURNS boolean
   4.735 +  LANGUAGE C IMMUTABLE STRICT
   4.736 +  AS '$libdir/latlon-v0008', 'pgl_epoint_ecluster_overlap';
   4.737 +
   4.738 +CREATE FUNCTION epoint_ecluster_may_overlap_proc(epoint, ecluster)
   4.739 +  RETURNS boolean
   4.740 +  LANGUAGE C IMMUTABLE STRICT
   4.741 +  AS '$libdir/latlon-v0008', 'pgl_epoint_ecluster_may_overlap';
   4.742 +
   4.743 +CREATE FUNCTION ebox_overlap_proc(ebox, ebox)
   4.744 +  RETURNS boolean
   4.745 +  LANGUAGE C IMMUTABLE STRICT
   4.746 +  AS '$libdir/latlon-v0008', 'pgl_ebox_overlap';
   4.747 +
   4.748 +CREATE FUNCTION ebox_ecircle_may_overlap_proc(ebox, ecircle)
   4.749 +  RETURNS boolean
   4.750 +  LANGUAGE C IMMUTABLE STRICT
   4.751 +  AS '$libdir/latlon-v0008', 'pgl_ebox_ecircle_may_overlap';
   4.752 +
   4.753 +CREATE FUNCTION ebox_ecluster_may_overlap_proc(ebox, ecluster)
   4.754 +  RETURNS boolean
   4.755 +  LANGUAGE C IMMUTABLE STRICT
   4.756 +  AS '$libdir/latlon-v0008', 'pgl_ebox_ecluster_may_overlap';
   4.757 +
   4.758 +CREATE FUNCTION ecircle_overlap_proc(ecircle, ecircle)
   4.759 +  RETURNS boolean
   4.760 +  LANGUAGE C IMMUTABLE STRICT
   4.761 +  AS '$libdir/latlon-v0008', 'pgl_ecircle_overlap';
   4.762 +
   4.763 +CREATE FUNCTION ecircle_ecluster_overlap_proc(ecircle, ecluster)
   4.764 +  RETURNS boolean
   4.765 +  LANGUAGE C IMMUTABLE STRICT
   4.766 +  AS '$libdir/latlon-v0008', 'pgl_ecircle_ecluster_overlap';
   4.767 +
   4.768 +CREATE FUNCTION ecircle_ecluster_may_overlap_proc(ecircle, ecluster)
   4.769 +  RETURNS boolean
   4.770 +  LANGUAGE C IMMUTABLE STRICT
   4.771 +  AS '$libdir/latlon-v0008', 'pgl_ecircle_ecluster_may_overlap';
   4.772 +
   4.773 +CREATE FUNCTION ecluster_overlap_proc(ecluster, ecluster)
   4.774 +  RETURNS boolean
   4.775 +  LANGUAGE C IMMUTABLE STRICT
   4.776 +  AS '$libdir/latlon-v0008', 'pgl_ecluster_overlap';
   4.777 +
   4.778 +CREATE FUNCTION ecluster_may_overlap_proc(ecluster, ecluster)
   4.779 +  RETURNS boolean
   4.780 +  LANGUAGE C IMMUTABLE STRICT
   4.781 +  AS '$libdir/latlon-v0008', 'pgl_ecluster_may_overlap';
   4.782 +
   4.783 +CREATE FUNCTION ecluster_contains_proc(ecluster, ecluster)
   4.784 +  RETURNS boolean
   4.785 +  LANGUAGE C IMMUTABLE STRICT
   4.786 +  AS '$libdir/latlon-v0008', 'pgl_ecluster_contains';
   4.787 +
   4.788 +CREATE FUNCTION epoint_distance_proc(epoint, epoint)
   4.789 +  RETURNS float8
   4.790 +  LANGUAGE C IMMUTABLE STRICT
   4.791 +  AS '$libdir/latlon-v0008', 'pgl_epoint_distance';
   4.792 +
   4.793 +CREATE FUNCTION epoint_ecircle_distance_proc(epoint, ecircle)
   4.794 +  RETURNS float8
   4.795 +  LANGUAGE C IMMUTABLE STRICT
   4.796 +  AS '$libdir/latlon-v0008', 'pgl_epoint_ecircle_distance';
   4.797 +
   4.798 +CREATE FUNCTION epoint_ecluster_distance_proc(epoint, ecluster)
   4.799 +  RETURNS float8
   4.800 +  LANGUAGE C IMMUTABLE STRICT
   4.801 +  AS '$libdir/latlon-v0008', 'pgl_epoint_ecluster_distance';
   4.802 +
   4.803 +CREATE FUNCTION ecircle_distance_proc(ecircle, ecircle)
   4.804 +  RETURNS float8
   4.805 +  LANGUAGE C IMMUTABLE STRICT
   4.806 +  AS '$libdir/latlon-v0008', 'pgl_ecircle_distance';
   4.807 +
   4.808 +CREATE FUNCTION ecircle_ecluster_distance_proc(ecircle, ecluster)
   4.809 +  RETURNS float8
   4.810 +  LANGUAGE C IMMUTABLE STRICT
   4.811 +  AS '$libdir/latlon-v0008', 'pgl_ecircle_ecluster_distance';
   4.812 +
   4.813 +CREATE FUNCTION ecluster_distance_proc(ecluster, ecluster)
   4.814 +  RETURNS float8
   4.815 +  LANGUAGE C IMMUTABLE STRICT
   4.816 +  AS '$libdir/latlon-v0008', 'pgl_ecluster_distance';
   4.817 +
   4.818 +CREATE OPERATOR && (
   4.819 +  leftarg = epoint,
   4.820 +  rightarg = ebox,
   4.821 +  procedure = epoint_ebox_overlap_proc,
   4.822 +  commutator = &&,
   4.823 +  restrict = areasel,
   4.824 +  join = areajoinsel
   4.825 +);
   4.826 +
   4.827 +CREATE FUNCTION epoint_ebox_overlap_commutator(ebox, epoint)
   4.828 +  RETURNS boolean
   4.829 +  LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1';
   4.830 +
   4.831 +CREATE OPERATOR && (
   4.832 +  leftarg = ebox,
   4.833 +  rightarg = epoint,
   4.834 +  procedure = epoint_ebox_overlap_commutator,
   4.835 +  commutator = &&,
   4.836 +  restrict = areasel,
   4.837 +  join = areajoinsel
   4.838 +);
   4.839 +
   4.840 +CREATE OPERATOR && (
   4.841 +  leftarg = epoint,
   4.842 +  rightarg = ecircle,
   4.843 +  procedure = epoint_ecircle_overlap_proc,
   4.844 +  commutator = &&,
   4.845 +  restrict = areasel,
   4.846 +  join = areajoinsel
   4.847 +);
   4.848 +
   4.849 +CREATE FUNCTION epoint_ecircle_overlap_commutator(ecircle, epoint)
   4.850 +  RETURNS boolean
   4.851 +  LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1';
   4.852 +
   4.853 +CREATE OPERATOR && (
   4.854 +  leftarg = ecircle,
   4.855 +  rightarg = epoint,
   4.856 +  procedure = epoint_ecircle_overlap_commutator,
   4.857 +  commutator = &&,
   4.858 +  restrict = areasel,
   4.859 +  join = areajoinsel
   4.860 +);
   4.861 +
   4.862 +CREATE OPERATOR && (
   4.863 +  leftarg = epoint,
   4.864 +  rightarg = ecluster,
   4.865 +  procedure = epoint_ecluster_overlap_proc,
   4.866 +  commutator = &&,
   4.867 +  restrict = areasel,
   4.868 +  join = areajoinsel
   4.869 +);
   4.870 +
   4.871 +CREATE FUNCTION epoint_ecluster_overlap_commutator(ecluster, epoint)
   4.872 +  RETURNS boolean
   4.873 +  LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1';
   4.874 +
   4.875 +CREATE OPERATOR && (
   4.876 +  leftarg = ecluster,
   4.877 +  rightarg = epoint,
   4.878 +  procedure = epoint_ecluster_overlap_commutator,
   4.879 +  commutator = &&,
   4.880 +  restrict = areasel,
   4.881 +  join = areajoinsel
   4.882 +);
   4.883 +
   4.884 +CREATE OPERATOR && (
   4.885 +  leftarg = ebox,
   4.886 +  rightarg = ebox,
   4.887 +  procedure = ebox_overlap_proc,
   4.888 +  commutator = &&,
   4.889 +  restrict = areasel,
   4.890 +  join = areajoinsel
   4.891 +);
   4.892 +
   4.893 +CREATE OPERATOR && (
   4.894 +  leftarg = ecircle,
   4.895 +  rightarg = ecircle,
   4.896 +  procedure = ecircle_overlap_proc,
   4.897 +  commutator = &&,
   4.898 +  restrict = areasel,
   4.899 +  join = areajoinsel
   4.900 +);
   4.901 +
   4.902 +CREATE OPERATOR && (
   4.903 +  leftarg = ecircle,
   4.904 +  rightarg = ecluster,
   4.905 +  procedure = ecircle_ecluster_overlap_proc,
   4.906 +  commutator = &&,
   4.907 +  restrict = areasel,
   4.908 +  join = areajoinsel
   4.909 +);
   4.910 +
   4.911 +CREATE FUNCTION ecircle_ecluster_overlap_commutator(ecluster, ecircle)
   4.912 +  RETURNS boolean
   4.913 +  LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1';
   4.914 +
   4.915 +CREATE OPERATOR && (
   4.916 +  leftarg = ecluster,
   4.917 +  rightarg = ecircle,
   4.918 +  procedure = ecircle_ecluster_overlap_commutator,
   4.919 +  commutator = &&,
   4.920 +  restrict = areasel,
   4.921 +  join = areajoinsel
   4.922 +);
   4.923 +
   4.924 +CREATE OPERATOR && (
   4.925 +  leftarg = ecluster,
   4.926 +  rightarg = ecluster,
   4.927 +  procedure = ecluster_overlap_proc,
   4.928 +  commutator = &&,
   4.929 +  restrict = areasel,
   4.930 +  join = areajoinsel
   4.931 +);
   4.932 +
   4.933 +CREATE FUNCTION ebox_ecircle_overlap_castwrap(ebox, ecircle)
   4.934 +  RETURNS boolean
   4.935 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster && $2';
   4.936 +
   4.937 +CREATE OPERATOR && (
   4.938 +  leftarg = ebox,
   4.939 +  rightarg = ecircle,
   4.940 +  procedure = ebox_ecircle_overlap_castwrap,
   4.941 +  commutator = &&,
   4.942 +  restrict = areasel,
   4.943 +  join = areajoinsel
   4.944 +);
   4.945 +
   4.946 +CREATE FUNCTION ebox_ecircle_overlap_castwrap(ecircle, ebox)
   4.947 +  RETURNS boolean
   4.948 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1 && $2::ecluster';
   4.949 +
   4.950 +CREATE OPERATOR && (
   4.951 +  leftarg = ecircle,
   4.952 +  rightarg = ebox,
   4.953 +  procedure = ebox_ecircle_overlap_castwrap,
   4.954 +  commutator = &&,
   4.955 +  restrict = areasel,
   4.956 +  join = areajoinsel
   4.957 +);
   4.958 +
   4.959 +CREATE FUNCTION ebox_ecluster_overlap_castwrap(ebox, ecluster)
   4.960 +  RETURNS boolean
   4.961 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster && $2';
   4.962 +
   4.963 +CREATE OPERATOR && (
   4.964 +  leftarg = ebox,
   4.965 +  rightarg = ecluster,
   4.966 +  procedure = ebox_ecluster_overlap_castwrap,
   4.967 +  commutator = &&,
   4.968 +  restrict = areasel,
   4.969 +  join = areajoinsel
   4.970 +);
   4.971 +
   4.972 +CREATE FUNCTION ebox_ecluster_overlap_castwrap(ecluster, ebox)
   4.973 +  RETURNS boolean
   4.974 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1 && $2::ecluster';
   4.975 +
   4.976 +CREATE OPERATOR && (
   4.977 +  leftarg = ecluster,
   4.978 +  rightarg = ebox,
   4.979 +  procedure = ebox_ecluster_overlap_castwrap,
   4.980 +  commutator = &&,
   4.981 +  restrict = areasel,
   4.982 +  join = areajoinsel
   4.983 +);
   4.984 +
   4.985 +CREATE OPERATOR &&+ (
   4.986 +  leftarg = epoint,
   4.987 +  rightarg = ecluster,
   4.988 +  procedure = epoint_ecluster_may_overlap_proc,
   4.989 +  commutator = &&+,
   4.990 +  restrict = areasel,
   4.991 +  join = areajoinsel
   4.992 +);
   4.993 +
   4.994 +CREATE FUNCTION epoint_ecluster_may_overlap_commutator(ecluster, epoint)
   4.995 +  RETURNS boolean
   4.996 +  LANGUAGE sql IMMUTABLE AS 'SELECT $2 &&+ $1';
   4.997 +
   4.998 +CREATE OPERATOR &&+ (
   4.999 +  leftarg = ecluster,
  4.1000 +  rightarg = epoint,
  4.1001 +  procedure = epoint_ecluster_may_overlap_commutator,
  4.1002 +  commutator = &&+,
  4.1003 +  restrict = areasel,
  4.1004 +  join = areajoinsel
  4.1005 +);
  4.1006 +
  4.1007 +CREATE OPERATOR &&+ (
  4.1008 +  leftarg = ebox,
  4.1009 +  rightarg = ecircle,
  4.1010 +  procedure = ebox_ecircle_may_overlap_proc,
  4.1011 +  commutator = &&+,
  4.1012 +  restrict = areasel,
  4.1013 +  join = areajoinsel
  4.1014 +);
  4.1015 +
  4.1016 +CREATE FUNCTION ebox_ecircle_may_overlap_commutator(ecircle, ebox)
  4.1017 +  RETURNS boolean
  4.1018 +  LANGUAGE sql IMMUTABLE AS 'SELECT $2 &&+ $1';
  4.1019 +
  4.1020 +CREATE OPERATOR &&+ (
  4.1021 +  leftarg = ecircle,
  4.1022 +  rightarg = ebox,
  4.1023 +  procedure = ebox_ecircle_may_overlap_commutator,
  4.1024 +  commutator = &&+,
  4.1025 +  restrict = areasel,
  4.1026 +  join = areajoinsel
  4.1027 +);
  4.1028 +
  4.1029 +CREATE OPERATOR &&+ (
  4.1030 +  leftarg = ebox,
  4.1031 +  rightarg = ecluster,
  4.1032 +  procedure = ebox_ecluster_may_overlap_proc,
  4.1033 +  commutator = &&+,
  4.1034 +  restrict = areasel,
  4.1035 +  join = areajoinsel
  4.1036 +);
  4.1037 +
  4.1038 +CREATE FUNCTION ebox_ecluster_may_overlap_commutator(ecluster, ebox)
  4.1039 +  RETURNS boolean
  4.1040 +  LANGUAGE sql IMMUTABLE AS 'SELECT $2 &&+ $1';
  4.1041 +
  4.1042 +CREATE OPERATOR &&+ (
  4.1043 +  leftarg = ecluster,
  4.1044 +  rightarg = ebox,
  4.1045 +  procedure = ebox_ecluster_may_overlap_commutator,
  4.1046 +  commutator = &&+,
  4.1047 +  restrict = areasel,
  4.1048 +  join = areajoinsel
  4.1049 +);
  4.1050 +
  4.1051 +CREATE OPERATOR &&+ (
  4.1052 +  leftarg = ecircle,
  4.1053 +  rightarg = ecluster,
  4.1054 +  procedure = ecircle_ecluster_may_overlap_proc,
  4.1055 +  commutator = &&+,
  4.1056 +  restrict = areasel,
  4.1057 +  join = areajoinsel
  4.1058 +);
  4.1059 +
  4.1060 +CREATE FUNCTION ecircle_ecluster_may_overlap_commutator(ecluster, ecircle)
  4.1061 +  RETURNS boolean
  4.1062 +  LANGUAGE sql IMMUTABLE AS 'SELECT $2 &&+ $1';
  4.1063 +
  4.1064 +CREATE OPERATOR &&+ (
  4.1065 +  leftarg = ecluster,
  4.1066 +  rightarg = ecircle,
  4.1067 +  procedure = ecircle_ecluster_may_overlap_commutator,
  4.1068 +  commutator = &&+,
  4.1069 +  restrict = areasel,
  4.1070 +  join = areajoinsel
  4.1071 +);
  4.1072 +
  4.1073 +CREATE OPERATOR &&+ (
  4.1074 +  leftarg = ecluster,
  4.1075 +  rightarg = ecluster,
  4.1076 +  procedure = ecluster_may_overlap_proc,
  4.1077 +  commutator = &&+,
  4.1078 +  restrict = areasel,
  4.1079 +  join = areajoinsel
  4.1080 +);
  4.1081 +
  4.1082 +CREATE OPERATOR @> (
  4.1083 +  leftarg = ebox,
  4.1084 +  rightarg = epoint,
  4.1085 +  procedure = epoint_ebox_overlap_commutator,
  4.1086 +  commutator = <@,
  4.1087 +  restrict = areasel,
  4.1088 +  join = areajoinsel
  4.1089 +);
  4.1090 +
  4.1091 +CREATE OPERATOR <@ (
  4.1092 +  leftarg = epoint,
  4.1093 +  rightarg = ebox,
  4.1094 +  procedure = epoint_ebox_overlap_proc,
  4.1095 +  commutator = @>,
  4.1096 +  restrict = areasel,
  4.1097 +  join = areajoinsel
  4.1098 +);
  4.1099 +
  4.1100 +CREATE OPERATOR @> (
  4.1101 +  leftarg = ecluster,
  4.1102 +  rightarg = epoint,
  4.1103 +  procedure = epoint_ecluster_overlap_commutator,
  4.1104 +  commutator = <@,
  4.1105 +  restrict = areasel,
  4.1106 +  join = areajoinsel
  4.1107 +);
  4.1108 +
  4.1109 +CREATE OPERATOR <@ (
  4.1110 +  leftarg = epoint,
  4.1111 +  rightarg = ecluster,
  4.1112 +  procedure = epoint_ecluster_overlap_proc,
  4.1113 +  commutator = <@,
  4.1114 +  restrict = areasel,
  4.1115 +  join = areajoinsel
  4.1116 +);
  4.1117 +
  4.1118 +CREATE OPERATOR @> (
  4.1119 +  leftarg = ecluster,
  4.1120 +  rightarg = ecluster,
  4.1121 +  procedure = ecluster_contains_proc,
  4.1122 +  commutator = <@,
  4.1123 +  restrict = areasel,
  4.1124 +  join = areajoinsel
  4.1125 +);
  4.1126 +
  4.1127 +CREATE FUNCTION ecluster_contains_commutator(ecluster, ecluster)
  4.1128 +  RETURNS boolean
  4.1129 +  LANGUAGE sql IMMUTABLE AS 'SELECT $2 @> $1';
  4.1130 +
  4.1131 +CREATE OPERATOR <@ (
  4.1132 +  leftarg = ecluster,
  4.1133 +  rightarg = ecluster,
  4.1134 +  procedure = ecluster_contains_commutator,
  4.1135 +  commutator = @>,
  4.1136 +  restrict = areasel,
  4.1137 +  join = areajoinsel
  4.1138 +);
  4.1139 +
  4.1140 +CREATE FUNCTION ebox_contains_castwrap(ebox, ebox)
  4.1141 +  RETURNS boolean
  4.1142 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster @> $2::ecluster';
  4.1143 +
  4.1144 +CREATE OPERATOR @> (
  4.1145 +  leftarg = ebox,
  4.1146 +  rightarg = ebox,
  4.1147 +  procedure = ebox_contains_castwrap,
  4.1148 +  commutator = <@,
  4.1149 +  restrict = areasel,
  4.1150 +  join = areajoinsel
  4.1151 +);
  4.1152 +
  4.1153 +CREATE FUNCTION ebox_contains_swapped_castwrap(ebox, ebox)
  4.1154 +  RETURNS boolean
  4.1155 +  LANGUAGE sql IMMUTABLE AS 'SELECT $2::ecluster @> $1::ecluster';
  4.1156 +
  4.1157 +CREATE OPERATOR <@ (
  4.1158 +  leftarg = ebox,
  4.1159 +  rightarg = ebox,
  4.1160 +  procedure = ebox_contains_swapped_castwrap,
  4.1161 +  commutator = @>,
  4.1162 +  restrict = areasel,
  4.1163 +  join = areajoinsel
  4.1164 +);
  4.1165 +
  4.1166 +CREATE FUNCTION ebox_ecluster_contains_castwrap(ebox, ecluster)
  4.1167 +  RETURNS boolean
  4.1168 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster @> $2';
  4.1169 +
  4.1170 +CREATE OPERATOR @> (
  4.1171 +  leftarg = ebox,
  4.1172 +  rightarg = ecluster,
  4.1173 +  procedure = ebox_ecluster_contains_castwrap,
  4.1174 +  commutator = <@,
  4.1175 +  restrict = areasel,
  4.1176 +  join = areajoinsel
  4.1177 +);
  4.1178 +
  4.1179 +CREATE FUNCTION ebox_ecluster_contains_castwrap(ecluster, ebox)
  4.1180 +  RETURNS boolean
  4.1181 +  LANGUAGE sql IMMUTABLE AS 'SELECT $2::ecluster @> $1';
  4.1182 +
  4.1183 +CREATE OPERATOR <@ (
  4.1184 +  leftarg = ecluster,
  4.1185 +  rightarg = ebox,
  4.1186 +  procedure = ebox_ecluster_contains_castwrap,
  4.1187 +  commutator = @>,
  4.1188 +  restrict = areasel,
  4.1189 +  join = areajoinsel
  4.1190 +);
  4.1191 +
  4.1192 +CREATE FUNCTION ecluster_ebox_contains_castwrap(ecluster, ebox)
  4.1193 +  RETURNS boolean
  4.1194 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1 @> $2::ecluster';
  4.1195 +
  4.1196 +CREATE OPERATOR @> (
  4.1197 +  leftarg = ecluster,
  4.1198 +  rightarg = ebox,
  4.1199 +  procedure = ecluster_ebox_contains_castwrap,
  4.1200 +  commutator = <@,
  4.1201 +  restrict = areasel,
  4.1202 +  join = areajoinsel
  4.1203 +);
  4.1204 +
  4.1205 +CREATE FUNCTION ecluster_ebox_contains_castwrap(ebox, ecluster)
  4.1206 +  RETURNS boolean
  4.1207 +  LANGUAGE sql IMMUTABLE AS 'SELECT $2 @> $1::ecluster';
  4.1208 +
  4.1209 +CREATE OPERATOR <@ (
  4.1210 +  leftarg = ebox,
  4.1211 +  rightarg = ecluster,
  4.1212 +  procedure = ecluster_ebox_contains_castwrap,
  4.1213 +  commutator = @>,
  4.1214 +  restrict = areasel,
  4.1215 +  join = areajoinsel
  4.1216 +);
  4.1217 +
  4.1218 +CREATE OPERATOR <-> (
  4.1219 +  leftarg = epoint,
  4.1220 +  rightarg = epoint,
  4.1221 +  procedure = epoint_distance_proc,
  4.1222 +  commutator = <->
  4.1223 +);
  4.1224 +
  4.1225 +CREATE OPERATOR <-> (
  4.1226 +  leftarg = epoint,
  4.1227 +  rightarg = ecircle,
  4.1228 +  procedure = epoint_ecircle_distance_proc,
  4.1229 +  commutator = <->
  4.1230 +);
  4.1231 +
  4.1232 +CREATE FUNCTION epoint_ecircle_distance_commutator(ecircle, epoint)
  4.1233 +  RETURNS float8
  4.1234 +  LANGUAGE sql IMMUTABLE AS 'SELECT $2 <-> $1';
  4.1235 +
  4.1236 +CREATE OPERATOR <-> (
  4.1237 +  leftarg = ecircle,
  4.1238 +  rightarg = epoint,
  4.1239 +  procedure = epoint_ecircle_distance_commutator,
  4.1240 +  commutator = <->
  4.1241 +);
  4.1242 +
  4.1243 +CREATE OPERATOR <-> (
  4.1244 +  leftarg = epoint,
  4.1245 +  rightarg = ecluster,
  4.1246 +  procedure = epoint_ecluster_distance_proc,
  4.1247 +  commutator = <->
  4.1248 +);
  4.1249 +
  4.1250 +CREATE FUNCTION epoint_ecluster_distance_commutator(ecluster, epoint)
  4.1251 +  RETURNS float8
  4.1252 +  LANGUAGE sql IMMUTABLE AS 'SELECT $2 <-> $1';
  4.1253 +
  4.1254 +CREATE OPERATOR <-> (
  4.1255 +  leftarg = ecluster,
  4.1256 +  rightarg = epoint,
  4.1257 +  procedure = epoint_ecluster_distance_commutator,
  4.1258 +  commutator = <->
  4.1259 +);
  4.1260 +
  4.1261 +CREATE OPERATOR <-> (
  4.1262 +  leftarg = ecircle,
  4.1263 +  rightarg = ecircle,
  4.1264 +  procedure = ecircle_distance_proc,
  4.1265 +  commutator = <->
  4.1266 +);
  4.1267 +
  4.1268 +CREATE OPERATOR <-> (
  4.1269 +  leftarg = ecircle,
  4.1270 +  rightarg = ecluster,
  4.1271 +  procedure = ecircle_ecluster_distance_proc,
  4.1272 +  commutator = <->
  4.1273 +);
  4.1274 +
  4.1275 +CREATE FUNCTION ecircle_ecluster_distance_commutator(ecluster, ecircle)
  4.1276 +  RETURNS float8
  4.1277 +  LANGUAGE sql IMMUTABLE AS 'SELECT $2 <-> $1';
  4.1278 +
  4.1279 +CREATE OPERATOR <-> (
  4.1280 +  leftarg = ecluster,
  4.1281 +  rightarg = ecircle,
  4.1282 +  procedure = ecircle_ecluster_distance_commutator,
  4.1283 +  commutator = <->
  4.1284 +);
  4.1285 +
  4.1286 +CREATE OPERATOR <-> (
  4.1287 +  leftarg = ecluster,
  4.1288 +  rightarg = ecluster,
  4.1289 +  procedure = ecluster_distance_proc,
  4.1290 +  commutator = <->
  4.1291 +);
  4.1292 +
  4.1293 +CREATE FUNCTION epoint_ebox_distance_castwrap(epoint, ebox)
  4.1294 +  RETURNS float8
  4.1295 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2::ecluster';
  4.1296 +
  4.1297 +CREATE OPERATOR <-> (
  4.1298 +  leftarg = epoint,
  4.1299 +  rightarg = ebox,
  4.1300 +  procedure = epoint_ebox_distance_castwrap,
  4.1301 +  commutator = <->
  4.1302 +);
  4.1303 +
  4.1304 +CREATE FUNCTION epoint_ebox_distance_castwrap(ebox, epoint)
  4.1305 +  RETURNS float8
  4.1306 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster <-> $2';
  4.1307 +
  4.1308 +CREATE OPERATOR <-> (
  4.1309 +  leftarg = ebox,
  4.1310 +  rightarg = epoint,
  4.1311 +  procedure = epoint_ebox_distance_castwrap,
  4.1312 +  commutator = <->
  4.1313 +);
  4.1314 +
  4.1315 +CREATE FUNCTION ebox_distance_castwrap(ebox, ebox)
  4.1316 +  RETURNS float8
  4.1317 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster <-> $2::ecluster';
  4.1318 +
  4.1319 +CREATE OPERATOR <-> (
  4.1320 +  leftarg = ebox,
  4.1321 +  rightarg = ebox,
  4.1322 +  procedure = ebox_distance_castwrap,
  4.1323 +  commutator = <->
  4.1324 +);
  4.1325 +
  4.1326 +CREATE FUNCTION ebox_ecircle_distance_castwrap(ebox, ecircle)
  4.1327 +  RETURNS float8
  4.1328 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster <-> $2';
  4.1329 +
  4.1330 +CREATE OPERATOR <-> (
  4.1331 +  leftarg = ebox,
  4.1332 +  rightarg = ecircle,
  4.1333 +  procedure = ebox_ecircle_distance_castwrap,
  4.1334 +  commutator = <->
  4.1335 +);
  4.1336 +
  4.1337 +CREATE FUNCTION ebox_ecircle_distance_castwrap(ecircle, ebox)
  4.1338 +  RETURNS float8
  4.1339 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2::ecluster';
  4.1340 +
  4.1341 +CREATE OPERATOR <-> (
  4.1342 +  leftarg = ecircle,
  4.1343 +  rightarg = ebox,
  4.1344 +  procedure = ebox_ecircle_distance_castwrap,
  4.1345 +  commutator = <->
  4.1346 +);
  4.1347 +
  4.1348 +CREATE FUNCTION ebox_ecluster_distance_castwrap(ebox, ecluster)
  4.1349 +  RETURNS float8
  4.1350 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster <-> $2';
  4.1351 +
  4.1352 +CREATE OPERATOR <-> (
  4.1353 +  leftarg = ebox,
  4.1354 +  rightarg = ecluster,
  4.1355 +  procedure = ebox_ecluster_distance_castwrap,
  4.1356 +  commutator = <->
  4.1357 +);
  4.1358 +
  4.1359 +CREATE FUNCTION ebox_ecluster_distance_castwrap(ecluster, ebox)
  4.1360 +  RETURNS float8
  4.1361 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2::ecluster';
  4.1362 +
  4.1363 +CREATE OPERATOR <-> (
  4.1364 +  leftarg = ecluster,
  4.1365 +  rightarg = ebox,
  4.1366 +  procedure = ebox_ecluster_distance_castwrap,
  4.1367 +  commutator = <->
  4.1368 +);
  4.1369 +
  4.1370 +
  4.1371 +---------------------
  4.1372 +-- other functions --
  4.1373 +---------------------
  4.1374 +
  4.1375 +CREATE FUNCTION monte_carlo_area(ecluster, int4 = 10000)
  4.1376 +  RETURNS float8
  4.1377 +  LANGUAGE C STRICT
  4.1378 +  AS '$libdir/latlon-v0008', 'pgl_ecluster_monte_carlo_area';
  4.1379 +
  4.1380 +CREATE FUNCTION fair_distance(ecluster, epoint, int4 = 10000)
  4.1381 +  RETURNS float8
  4.1382 +  LANGUAGE C STRICT
  4.1383 +  AS '$libdir/latlon-v0008', 'pgl_ecluster_epoint_fair_distance';
  4.1384 +
  4.1385 +
  4.1386 +----------------
  4.1387 +-- GiST index --
  4.1388 +----------------
  4.1389 +
  4.1390 +CREATE FUNCTION pgl_gist_consistent(internal, internal, smallint, oid, internal)
  4.1391 +  RETURNS boolean
  4.1392 +  LANGUAGE C STRICT
  4.1393 +  AS '$libdir/latlon-v0008', 'pgl_gist_consistent';
  4.1394 +
  4.1395 +CREATE FUNCTION pgl_gist_union(internal, internal)
  4.1396 +  RETURNS internal
  4.1397 +  LANGUAGE C STRICT
  4.1398 +  AS '$libdir/latlon-v0008', 'pgl_gist_union';
  4.1399 +
  4.1400 +CREATE FUNCTION pgl_gist_compress_epoint(internal)
  4.1401 +  RETURNS internal
  4.1402 +  LANGUAGE C STRICT
  4.1403 +  AS '$libdir/latlon-v0008', 'pgl_gist_compress_epoint';
  4.1404 +
  4.1405 +CREATE FUNCTION pgl_gist_compress_ecircle(internal)
  4.1406 +  RETURNS internal
  4.1407 +  LANGUAGE C STRICT
  4.1408 +  AS '$libdir/latlon-v0008', 'pgl_gist_compress_ecircle';
  4.1409 +
  4.1410 +CREATE FUNCTION pgl_gist_compress_ecluster(internal)
  4.1411 +  RETURNS internal
  4.1412 +  LANGUAGE C STRICT
  4.1413 +  AS '$libdir/latlon-v0008', 'pgl_gist_compress_ecluster';
  4.1414 +
  4.1415 +CREATE FUNCTION pgl_gist_decompress(internal)
  4.1416 +  RETURNS internal
  4.1417 +  LANGUAGE C STRICT
  4.1418 +  AS '$libdir/latlon-v0008', 'pgl_gist_decompress';
  4.1419 +
  4.1420 +CREATE FUNCTION pgl_gist_penalty(internal, internal, internal)
  4.1421 +  RETURNS internal
  4.1422 +  LANGUAGE C STRICT
  4.1423 +  AS '$libdir/latlon-v0008', 'pgl_gist_penalty';
  4.1424 +
  4.1425 +CREATE FUNCTION pgl_gist_picksplit(internal, internal)
  4.1426 +  RETURNS internal
  4.1427 +  LANGUAGE C STRICT
  4.1428 +  AS '$libdir/latlon-v0008', 'pgl_gist_picksplit';
  4.1429 +
  4.1430 +CREATE FUNCTION pgl_gist_same(internal, internal, internal)
  4.1431 +  RETURNS internal
  4.1432 +  LANGUAGE C STRICT
  4.1433 +  AS '$libdir/latlon-v0008', 'pgl_gist_same';
  4.1434 +
  4.1435 +CREATE FUNCTION pgl_gist_distance(internal, internal, smallint, oid)
  4.1436 +  RETURNS internal
  4.1437 +  LANGUAGE C STRICT
  4.1438 +  AS '$libdir/latlon-v0008', 'pgl_gist_distance';
  4.1439 +
  4.1440 +CREATE OPERATOR CLASS epoint_ops
  4.1441 +  DEFAULT FOR TYPE epoint USING gist AS
  4.1442 +  OPERATOR  11 = ,
  4.1443 +  OPERATOR  22 &&  (epoint, ebox),
  4.1444 +  OPERATOR 222 <@  (epoint, ebox),
  4.1445 +  OPERATOR  23 &&  (epoint, ecircle),
  4.1446 +  OPERATOR  24 &&  (epoint, ecluster),
  4.1447 +  OPERATOR 124 &&+ (epoint, ecluster),
  4.1448 +  OPERATOR 224 <@  (epoint, ecluster),
  4.1449 +  OPERATOR  31 <-> (epoint, epoint) FOR ORDER BY float_ops,
  4.1450 +  OPERATOR  32 <-> (epoint, ebox) FOR ORDER BY float_ops,
  4.1451 +  OPERATOR  33 <-> (epoint, ecircle) FOR ORDER BY float_ops,
  4.1452 +  OPERATOR  34 <-> (epoint, ecluster) FOR ORDER BY float_ops,
  4.1453 +  FUNCTION 1 pgl_gist_consistent(internal, internal, smallint, oid, internal),
  4.1454 +  FUNCTION 2 pgl_gist_union(internal, internal),
  4.1455 +  FUNCTION 3 pgl_gist_compress_epoint(internal),
  4.1456 +  FUNCTION 4 pgl_gist_decompress(internal),
  4.1457 +  FUNCTION 5 pgl_gist_penalty(internal, internal, internal),
  4.1458 +  FUNCTION 6 pgl_gist_picksplit(internal, internal),
  4.1459 +  FUNCTION 7 pgl_gist_same(internal, internal, internal),
  4.1460 +  FUNCTION 8 pgl_gist_distance(internal, internal, smallint, oid),
  4.1461 +  STORAGE ekey_point;
  4.1462 +
  4.1463 +CREATE OPERATOR CLASS ecircle_ops
  4.1464 +  DEFAULT FOR TYPE ecircle USING gist AS
  4.1465 +  OPERATOR  13 = ,
  4.1466 +  OPERATOR  21 &&  (ecircle, epoint),
  4.1467 +  OPERATOR  22 &&  (ecircle, ebox),
  4.1468 +  OPERATOR 122 &&+ (ecircle, ebox),
  4.1469 +  OPERATOR  23 &&  (ecircle, ecircle),
  4.1470 +  OPERATOR  24 &&  (ecircle, ecluster),
  4.1471 +  OPERATOR 124 &&+ (ecircle, ecluster),
  4.1472 +  OPERATOR  31 <-> (ecircle, epoint) FOR ORDER BY float_ops,
  4.1473 +  OPERATOR  32 <-> (ecircle, ebox) FOR ORDER BY float_ops,
  4.1474 +  OPERATOR  33 <-> (ecircle, ecircle) FOR ORDER BY float_ops,
  4.1475 +  OPERATOR  34 <-> (ecircle, ecluster) FOR ORDER BY float_ops,
  4.1476 +  FUNCTION 1 pgl_gist_consistent(internal, internal, smallint, oid, internal),
  4.1477 +  FUNCTION 2 pgl_gist_union(internal, internal),
  4.1478 +  FUNCTION 3 pgl_gist_compress_ecircle(internal),
  4.1479 +  FUNCTION 4 pgl_gist_decompress(internal),
  4.1480 +  FUNCTION 5 pgl_gist_penalty(internal, internal, internal),
  4.1481 +  FUNCTION 6 pgl_gist_picksplit(internal, internal),
  4.1482 +  FUNCTION 7 pgl_gist_same(internal, internal, internal),
  4.1483 +  FUNCTION 8 pgl_gist_distance(internal, internal, smallint, oid),
  4.1484 +  STORAGE ekey_area;
  4.1485 +
  4.1486 +CREATE OPERATOR CLASS ecluster_ops
  4.1487 +  DEFAULT FOR TYPE ecluster USING gist AS
  4.1488 +  OPERATOR  21 &&  (ecluster, epoint),
  4.1489 +  OPERATOR 121 &&+ (ecluster, epoint),
  4.1490 +  OPERATOR 221 @>  (ecluster, epoint),
  4.1491 +  OPERATOR  22 &&  (ecluster, ebox),
  4.1492 +  OPERATOR 122 &&+ (ecluster, ebox),
  4.1493 +  OPERATOR 222 @>  (ecluster, ebox),
  4.1494 +  OPERATOR 322 <@  (ecluster, ebox),
  4.1495 +  OPERATOR  23 &&  (ecluster, ecircle),
  4.1496 +  OPERATOR 123 &&+ (ecluster, ecircle),
  4.1497 +  OPERATOR  24 &&  (ecluster, ecluster),
  4.1498 +  OPERATOR 124 &&+ (ecluster, ecluster),
  4.1499 +  OPERATOR 224 @>  (ecluster, ecluster),
  4.1500 +  OPERATOR 324 <@  (ecluster, ecluster),
  4.1501 +  OPERATOR  31 <-> (ecluster, epoint) FOR ORDER BY float_ops,
  4.1502 +  OPERATOR  32 <-> (ecluster, ebox) FOR ORDER BY float_ops,
  4.1503 +  OPERATOR  33 <-> (ecluster, ecircle) FOR ORDER BY float_ops,
  4.1504 +  OPERATOR  34 <-> (ecluster, ecluster) FOR ORDER BY float_ops,
  4.1505 +  FUNCTION 1 pgl_gist_consistent(internal, internal, smallint, oid, internal),
  4.1506 +  FUNCTION 2 pgl_gist_union(internal, internal),
  4.1507 +  FUNCTION 3 pgl_gist_compress_ecluster(internal),
  4.1508 +  FUNCTION 4 pgl_gist_decompress(internal),
  4.1509 +  FUNCTION 5 pgl_gist_penalty(internal, internal, internal),
  4.1510 +  FUNCTION 6 pgl_gist_picksplit(internal, internal),
  4.1511 +  FUNCTION 7 pgl_gist_same(internal, internal, internal),
  4.1512 +  FUNCTION 8 pgl_gist_distance(internal, internal, smallint, oid),
  4.1513 +  STORAGE ekey_area;
  4.1514 +
  4.1515 +
  4.1516 +---------------------
  4.1517 +-- alias functions --
  4.1518 +---------------------
  4.1519 +
  4.1520 +CREATE FUNCTION distance(epoint, epoint)
  4.1521 +  RETURNS float8
  4.1522 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2';
  4.1523 +
  4.1524 +CREATE FUNCTION distance(ecluster, epoint)
  4.1525 +  RETURNS float8
  4.1526 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2';
  4.1527 +
  4.1528 +CREATE FUNCTION distance_within(epoint, epoint, float8)
  4.1529 +  RETURNS boolean
  4.1530 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1 && ecircle($2, $3)';
  4.1531 +
  4.1532 +CREATE FUNCTION distance_within(ecluster, epoint, float8)
  4.1533 +  RETURNS boolean
  4.1534 +  LANGUAGE sql IMMUTABLE AS 'SELECT $1 && ecircle($2, $3)';
  4.1535 +
  4.1536 +
  4.1537 +--------------------------------
  4.1538 +-- other data storage formats --
  4.1539 +--------------------------------
  4.1540 +
  4.1541 +CREATE FUNCTION coords_to_epoint(float8, float8, text = 'epoint')
  4.1542 +  RETURNS epoint
  4.1543 +  LANGUAGE plpgsql IMMUTABLE STRICT AS $$
  4.1544 +    DECLARE
  4.1545 +      "result" epoint;
  4.1546 +    BEGIN
  4.1547 +      IF $3 = 'epoint_lonlat' THEN
  4.1548 +        -- avoid dynamic command execution for better performance
  4.1549 +        RETURN epoint($2, $1);
  4.1550 +      END IF;
  4.1551 +      IF $3 = 'epoint' OR $3 = 'epoint_latlon' THEN
  4.1552 +        -- avoid dynamic command execution for better performance
  4.1553 +        RETURN epoint($1, $2);
  4.1554 +      END IF;
  4.1555 +      EXECUTE 'SELECT ' || $3 || '($1, $2)' INTO STRICT "result" USING $1, $2;
  4.1556 +      RETURN "result";
  4.1557 +    END;
  4.1558 +  $$;
  4.1559 +
  4.1560 +CREATE FUNCTION GeoJSON_LinearRing_vertices(jsonb, text = 'epoint_lonlat')
  4.1561 +  RETURNS SETOF jsonb
  4.1562 +  LANGUAGE sql IMMUTABLE STRICT AS $$
  4.1563 +    SELECT "result" FROM
  4.1564 +      ( SELECT jsonb_array_length($1) - 1 ) AS "lastindex_row" ("lastindex")
  4.1565 +      CROSS JOIN LATERAL jsonb_array_elements(
  4.1566 +        CASE WHEN
  4.1567 +          coords_to_epoint(
  4.1568 +            ($1->0->>0)::float8,
  4.1569 +            ($1->0->>1)::float8,
  4.1570 +            $2
  4.1571 +          ) = coords_to_epoint(
  4.1572 +            ($1->"lastindex"->>0)::float8,
  4.1573 +            ($1->"lastindex"->>1)::float8,
  4.1574 +            $2
  4.1575 +          )
  4.1576 +        THEN
  4.1577 +          $1 - "lastindex"
  4.1578 +        ELSE
  4.1579 +          $1
  4.1580 +        END
  4.1581 +      ) AS "result_row" ("result")
  4.1582 +  $$;
  4.1583 +
  4.1584 +CREATE FUNCTION GeoJSON_to_epoint(jsonb, text = 'epoint_lonlat')
  4.1585 +  RETURNS epoint
  4.1586 +  LANGUAGE sql IMMUTABLE STRICT AS $$
  4.1587 +    SELECT CASE
  4.1588 +    WHEN $1->>'type' = 'Point' THEN
  4.1589 +      coords_to_epoint(
  4.1590 +        ($1->'coordinates'->>0)::float8,
  4.1591 +        ($1->'coordinates'->>1)::float8,
  4.1592 +        $2
  4.1593 +      )
  4.1594 +    WHEN $1->>'type' = 'Feature' THEN
  4.1595 +      GeoJSON_to_epoint($1->'geometry', $2)
  4.1596 +    ELSE
  4.1597 +      NULL
  4.1598 +    END
  4.1599 +  $$;
  4.1600 +
  4.1601 +CREATE FUNCTION GeoJSON_to_ecluster(jsonb, text = 'epoint_lonlat')
  4.1602 +  RETURNS ecluster
  4.1603 +  LANGUAGE sql IMMUTABLE STRICT AS $$
  4.1604 +    SELECT CASE $1->>'type'
  4.1605 +    WHEN 'Point' THEN
  4.1606 +      coords_to_epoint(
  4.1607 +        ($1->'coordinates'->>0)::float8,
  4.1608 +        ($1->'coordinates'->>1)::float8,
  4.1609 +        $2
  4.1610 +      )::ecluster
  4.1611 +    WHEN 'MultiPoint' THEN
  4.1612 +      ( SELECT ecluster_create_multipoint(array_agg(
  4.1613 +          coords_to_epoint(
  4.1614 +            ("coord"->>0)::float8,
  4.1615 +            ("coord"->>1)::float8,
  4.1616 +            $2
  4.1617 +          )
  4.1618 +        ))
  4.1619 +        FROM jsonb_array_elements($1->'coordinates') AS "coord"
  4.1620 +      )
  4.1621 +    WHEN 'LineString' THEN
  4.1622 +      ( SELECT ecluster_create_path(array_agg(
  4.1623 +          coords_to_epoint(
  4.1624 +            ("coord"->>0)::float8,
  4.1625 +            ("coord"->>1)::float8,
  4.1626 +            $2
  4.1627 +          )
  4.1628 +        ))
  4.1629 +        FROM jsonb_array_elements($1->'coordinates') AS "coord"
  4.1630 +      )
  4.1631 +    WHEN 'MultiLineString' THEN
  4.1632 +      ( SELECT ecluster_concat(array_agg(
  4.1633 +          ( SELECT ecluster_create_path(array_agg(
  4.1634 +              coords_to_epoint(
  4.1635 +                ("coord"->>0)::float8,
  4.1636 +                ("coord"->>1)::float8,
  4.1637 +                $2
  4.1638 +              )
  4.1639 +            ))
  4.1640 +            FROM jsonb_array_elements("coord_array") AS "coord"
  4.1641 +          )
  4.1642 +        ))
  4.1643 +        FROM jsonb_array_elements($1->'coordinates') AS "coord_array"
  4.1644 +      )
  4.1645 +    WHEN 'Polygon' THEN
  4.1646 +      ( SELECT ecluster_concat(array_agg(
  4.1647 +          ( SELECT ecluster_create_polygon(array_agg(
  4.1648 +              coords_to_epoint(
  4.1649 +                ("coord"->>0)::float8,
  4.1650 +                ("coord"->>1)::float8,
  4.1651 +                $2
  4.1652 +              )
  4.1653 +            ))
  4.1654 +            FROM GeoJSON_LinearRing_vertices("coord_array", $2) AS "coord"
  4.1655 +          )
  4.1656 +        ))
  4.1657 +        FROM jsonb_array_elements($1->'coordinates') AS "coord_array"
  4.1658 +      )
  4.1659 +    WHEN 'MultiPolygon' THEN
  4.1660 +      ( SELECT ecluster_concat(array_agg(
  4.1661 +          ( SELECT ecluster_concat(array_agg(
  4.1662 +              ( SELECT ecluster_create_polygon(array_agg(
  4.1663 +                  coords_to_epoint(
  4.1664 +                    ("coord"->>0)::float8,
  4.1665 +                    ("coord"->>1)::float8,
  4.1666 +                    $2
  4.1667 +                  )
  4.1668 +                ))
  4.1669 +                FROM GeoJSON_LinearRing_vertices("coord_array", $2) AS "coord"
  4.1670 +              )
  4.1671 +            ))
  4.1672 +            FROM jsonb_array_elements("coord_array_array") AS "coord_array"
  4.1673 +          )
  4.1674 +        ))
  4.1675 +        FROM jsonb_array_elements($1->'coordinates') AS "coord_array_array"
  4.1676 +      )
  4.1677 +    WHEN 'GeometryCollection' THEN
  4.1678 +      ( SELECT ecluster_concat(array_agg(
  4.1679 +          GeoJSON_to_ecluster("geometry", $2)
  4.1680 +        ))
  4.1681 +        FROM jsonb_array_elements($1->'geometries') AS "geometry"
  4.1682 +      )
  4.1683 +    WHEN 'Feature' THEN
  4.1684 +      GeoJSON_to_ecluster($1->'geometry', $2)
  4.1685 +    WHEN 'FeatureCollection' THEN
  4.1686 +      ( SELECT ecluster_concat(array_agg(
  4.1687 +          GeoJSON_to_ecluster("feature", $2)
  4.1688 +        ))
  4.1689 +        FROM jsonb_array_elements($1->'features') AS "feature"
  4.1690 +      )
  4.1691 +    ELSE
  4.1692 +      NULL
  4.1693 +    END
  4.1694 +  $$;
  4.1695 +

     5.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     5.2 +++ b/latlon-v0008.c	Tue Oct 25 18:44:43 2016 +0200
     5.3 @@ -0,0 +1,3285 @@
     5.4 +
     5.5 +/*-------------*
     5.6 + *  C prelude  *
     5.7 + *-------------*/
     5.8 +
     5.9 +#include "postgres.h"
    5.10 +#include "fmgr.h"
    5.11 +#include "libpq/pqformat.h"
    5.12 +#include "access/gist.h"
    5.13 +#include "access/stratnum.h"
    5.14 +#include "utils/array.h"
    5.15 +#include <math.h>
    5.16 +
    5.17 +#ifdef PG_MODULE_MAGIC
    5.18 +PG_MODULE_MAGIC;
    5.19 +#endif
    5.20 +
    5.21 +#if INT_MAX < 2147483647
    5.22 +#error Expected int type to be at least 32 bit wide
    5.23 +#endif
    5.24 +
    5.25 +
    5.26 +/*---------------------------------*
    5.27 + *  distance calculation on earth  *
    5.28 + *  (using WGS-84 spheroid)        *
    5.29 + *---------------------------------*/
    5.30 +
    5.31 +/*  WGS-84 spheroid with following parameters:
    5.32 +    semi-major axis  a = 6378137
    5.33 +    semi-minor axis  b = a * (1 - 1/298.257223563)
    5.34 +    estimated diameter = 2 * (2*a+b)/3
    5.35 +*/
    5.36 +#define PGL_SPHEROID_A 6378137.0            /* semi major axis */
    5.37 +#define PGL_SPHEROID_F (1.0/298.257223563)  /* flattening */
    5.38 +#define PGL_SPHEROID_B (PGL_SPHEROID_A * (1.0-PGL_SPHEROID_F))
    5.39 +#define PGL_EPS2       ( ( PGL_SPHEROID_A * PGL_SPHEROID_A - \
    5.40 +                           PGL_SPHEROID_B * PGL_SPHEROID_B ) / \
    5.41 +                         ( PGL_SPHEROID_A * PGL_SPHEROID_A ) )
    5.42 +#define PGL_SUBEPS2    (1.0-PGL_EPS2)
    5.43 +#define PGL_RADIUS     ((2.0*PGL_SPHEROID_A + PGL_SPHEROID_B) / 3.0)
    5.44 +#define PGL_DIAMETER   (2.0 * PGL_RADIUS)
    5.45 +#define PGL_SCALE      (PGL_SPHEROID_A / PGL_DIAMETER)  /* semi-major ref. */
    5.46 +#define PGL_MAXDIST    (PGL_RADIUS * M_PI)              /* maximum distance */
    5.47 +#define PGL_FADELIMIT  (PGL_MAXDIST / 3.0)              /* 1/6 circumference */
    5.48 +
    5.49 +/* calculate distance between two points on earth (given in degrees) */
    5.50 +static inline double pgl_distance(
    5.51 +  double lat1, double lon1, double lat2, double lon2
    5.52 +) {
    5.53 +  float8 lat1cos, lat1sin, lat2cos, lat2sin, lon2cos, lon2sin;
    5.54 +  float8 nphi1, nphi2, x1, z1, x2, y2, z2, g, s, t;
    5.55 +  /* normalize delta longitude (lon2 > 0 && lon1 = 0) */
    5.56 +  /* lon1 = 0 (not used anymore) */
    5.57 +  lon2 = fabs(lon2-lon1);
    5.58 +  /* convert to radians (first divide, then multiply) */
    5.59 +  lat1 = (lat1 / 180.0) * M_PI;
    5.60 +  lat2 = (lat2 / 180.0) * M_PI;
    5.61 +  lon2 = (lon2 / 180.0) * M_PI;
    5.62 +  /* make lat2 >= lat1 to ensure reversal-symmetry despite floating point
    5.63 +     operations (lon2 >= lon1 is already ensured in a previous step) */
    5.64 +  if (lat2 < lat1) { float8 swap = lat1; lat1 = lat2; lat2 = swap; }
    5.65 +  /* calculate 3d coordinates on scaled ellipsoid which has an average diameter
    5.66 +     of 1.0 */
    5.67 +  lat1cos = cos(lat1); lat1sin = sin(lat1);
    5.68 +  lat2cos = cos(lat2); lat2sin = sin(lat2);
    5.69 +  lon2cos = cos(lon2); lon2sin = sin(lon2);
    5.70 +  nphi1 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat1sin * lat1sin);
    5.71 +  nphi2 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat2sin * lat2sin);
    5.72 +  x1 = nphi1 * lat1cos;
    5.73 +  z1 = nphi1 * PGL_SUBEPS2 * lat1sin;
    5.74 +  x2 = nphi2 * lat2cos * lon2cos;
    5.75 +  y2 = nphi2 * lat2cos * lon2sin;
    5.76 +  z2 = nphi2 * PGL_SUBEPS2 * lat2sin;
    5.77 +  /* calculate tunnel distance through scaled (diameter 1.0) ellipsoid */
    5.78 +  g = sqrt((x2-x1)*(x2-x1) + y2*y2 + (z2-z1)*(z2-z1));
    5.79 +  /* convert tunnel distance through scaled ellipsoid to approximated surface
    5.80 +     distance on original ellipsoid */
    5.81 +  if (g > 1.0) g = 1.0;
    5.82 +  s = PGL_DIAMETER * asin(g);
    5.83 +  /* return result only if small enough to be precise (less than 1/3 of
    5.84 +     maximum possible distance) */
    5.85 +  if (s <= PGL_FADELIMIT) return s;
    5.86 +  /* calculate tunnel distance to antipodal point through scaled ellipsoid */
    5.87 +  g = sqrt((x2+x1)*(x2+x1) + y2*y2 + (z2+z1)*(z2+z1));
    5.88 +  /* convert tunnel distance to antipodal point through scaled ellipsoid to
    5.89 +     approximated surface distance to antipodal point on original ellipsoid */
    5.90 +  if (g > 1.0) g = 1.0;
    5.91 +  t = PGL_DIAMETER * asin(g);
    5.92 +  /* surface distance between original points can now be approximated by
    5.93 +     substracting antipodal distance from maximum possible distance;
    5.94 +     return result only if small enough (less than 1/3 of maximum possible
    5.95 +     distance) */
    5.96 +  if (t <= PGL_FADELIMIT) return PGL_MAXDIST-t;
    5.97 +  /* otherwise crossfade direct and antipodal result to ensure monotonicity */
    5.98 +  return (
    5.99 +    (s * (t-PGL_FADELIMIT) + (PGL_MAXDIST-t) * (s-PGL_FADELIMIT)) /
   5.100 +    (s + t - 2*PGL_FADELIMIT)
   5.101 +  );
   5.102 +}
   5.103 +
   5.104 +/* finite distance that can not be reached on earth */
   5.105 +#define PGL_ULTRA_DISTANCE (3 * PGL_MAXDIST)
   5.106 +
   5.107 +
   5.108 +/*--------------------------------*
   5.109 + *  simple geographic data types  *
   5.110 + *--------------------------------*/
   5.111 +
   5.112 +/* point on earth given by latitude and longitude in degrees */
   5.113 +/* (type "epoint" in SQL) */
   5.114 +typedef struct {
   5.115 +  double lat;  /* between  -90 and  90 (both inclusive) */
   5.116 +  double lon;  /* between -180 and 180 (both inclusive) */
   5.117 +} pgl_point;
   5.118 +
   5.119 +/* box delimited by two parallels and two meridians (all in degrees) */
   5.120 +/* (type "ebox" in SQL) */
   5.121 +typedef struct {
   5.122 +  double lat_min;  /* between  -90 and  90 (both inclusive) */
   5.123 +  double lat_max;  /* between  -90 and  90 (both inclusive) */
   5.124 +  double lon_min;  /* between -180 and 180 (both inclusive) */
   5.125 +  double lon_max;  /* between -180 and 180 (both inclusive) */
   5.126 +  /* if lat_min > lat_max, then box is empty */
   5.127 +  /* if lon_min > lon_max, then 180th meridian is crossed */
   5.128 +} pgl_box;
   5.129 +
   5.130 +/* circle on earth surface (for radial searches with fixed radius) */
   5.131 +/* (type "ecircle" in SQL) */
   5.132 +typedef struct {
   5.133 +  pgl_point center;
   5.134 +  double radius; /* positive (including +0 but excluding -0), or -INFINITY */
   5.135 +  /* A negative radius (i.e. -INFINITY) denotes nothing (i.e. no point),
   5.136 +     zero radius (0) denotes a single point,
   5.137 +     a finite radius (0 < radius < INFINITY) denotes a filled circle, and
   5.138 +     a radius of INFINITY is valid and means complete coverage of earth. */
   5.139 +} pgl_circle;
   5.140 +
   5.141 +
   5.142 +/*----------------------------------*
   5.143 + *  geographic "cluster" data type  *
   5.144 + *----------------------------------*/
   5.145 +
   5.146 +/* A cluster is a collection of points, paths, outlines, and polygons. If two
   5.147 +   polygons in a cluster overlap, the area covered by both polygons does not
   5.148 +   belong to the cluster. This way, a cluster can be used to describe complex
   5.149 +   shapes like polygons with holes. Outlines are non-filled polygons. Paths are
   5.150 +   open by default (i.e. the last point in the list is not connected with the
   5.151 +   first point in the list). Note that each outline or polygon in a cluster
   5.152 +   must cover a longitude range of less than 180 degrees to avoid ambiguities.
   5.153 +   Areas which are larger may be split into multiple polygons. */
   5.154 +
   5.155 +/* maximum number of points in a cluster */
   5.156 +/* (limited to avoid integer overflows, e.g. when allocating memory) */
   5.157 +#define PGL_CLUSTER_MAXPOINTS 16777216
   5.158 +
   5.159 +/* types of cluster entries */
   5.160 +#define PGL_ENTRY_POINT   1  /* a point */
   5.161 +#define PGL_ENTRY_PATH    2  /* a path from first point to last point */
   5.162 +#define PGL_ENTRY_OUTLINE 3  /* a non-filled polygon with given vertices */
   5.163 +#define PGL_ENTRY_POLYGON 4  /* a filled polygon with given vertices */
   5.164 +
   5.165 +/* Entries of a cluster are described by two different structs: pgl_newentry
   5.166 +   and pgl_entry. The first is used only during construction of a cluster, the
   5.167 +   second is used in all other cases (e.g. when reading clusters from the
   5.168 +   database, performing operations, etc). */
   5.169 +
   5.170 +/* entry for new geographic cluster during construction of that cluster */
   5.171 +typedef struct {
   5.172 +  int32_t entrytype;
   5.173 +  int32_t npoints;
   5.174 +  pgl_point *points;  /* pointer to an array of points (pgl_point) */
   5.175 +} pgl_newentry;
   5.176 +
   5.177 +/* entry of geographic cluster */
   5.178 +typedef struct {
   5.179 +  int32_t entrytype;  /* type of entry: point, path, outline, polygon */
   5.180 +  int32_t npoints;    /* number of stored points (set to 1 for point entry) */
   5.181 +  int32_t offset;     /* offset of pgl_point array from cluster base address */
   5.182 +  /* use macro PGL_ENTRY_POINTS to obtain a pointer to the array of points */
   5.183 +} pgl_entry;
   5.184 +
   5.185 +/* geographic cluster which is a collection of points, (open) paths, polygons,
   5.186 +   and outlines (non-filled polygons) */
   5.187 +typedef struct {
   5.188 +  char header[VARHDRSZ];  /* PostgreSQL header for variable size data types */
   5.189 +  int32_t nentries;       /* number of stored points */
   5.190 +  pgl_circle bounding;    /* bounding circle */
   5.191 +  /* Note: bounding circle ensures alignment of pgl_cluster for points */
   5.192 +  pgl_entry entries[FLEXIBLE_ARRAY_MEMBER];  /* var-length data */
   5.193 +} pgl_cluster;
   5.194 +
   5.195 +/* macro to determine memory alignment of points */
   5.196 +/* (needed to store pgl_point array after entries in pgl_cluster) */
   5.197 +typedef struct { char dummy; pgl_point aligned; } pgl_point_alignment;
   5.198 +#define PGL_POINT_ALIGNMENT offsetof(pgl_point_alignment, aligned)
   5.199 +
   5.200 +/* macro to extract a pointer to the array of points of a cluster entry */
   5.201 +#define PGL_ENTRY_POINTS(cluster, idx) \
   5.202 +  ((pgl_point *)(((intptr_t)cluster)+(cluster)->entries[idx].offset))
   5.203 +
   5.204 +/* convert pgl_newentry array to pgl_cluster */
   5.205 +/* NOTE: requires pgl_finalize_cluster to be called to finalize result */
   5.206 +static pgl_cluster *pgl_new_cluster(int nentries, pgl_newentry *entries) {
   5.207 +  int i;              /* index of current entry */
   5.208 +  int npoints = 0;    /* number of points in whole cluster */
   5.209 +  int entry_npoints;  /* number of points in current entry */
   5.210 +  int points_offset = PGL_POINT_ALIGNMENT * (
   5.211 +    ( offsetof(pgl_cluster, entries) +
   5.212 +      nentries * sizeof(pgl_entry) +
   5.213 +      PGL_POINT_ALIGNMENT - 1
   5.214 +    ) / PGL_POINT_ALIGNMENT
   5.215 +  );  /* offset of pgl_point array from base address (considering alignment) */
   5.216 +  pgl_cluster *cluster;  /* new cluster to be returned */
   5.217 +  /* determine total number of points */
   5.218 +  for (i=0; i<nentries; i++) npoints += entries[i].npoints;
   5.219 +  /* allocate memory for cluster (including entries and points) */
   5.220 +  cluster = palloc(points_offset + npoints * sizeof(pgl_point));
   5.221 +  /* re-count total number of points to determine offset for each entry */
   5.222 +  npoints = 0;
   5.223 +  /* copy entries and points */
   5.224 +  for (i=0; i<nentries; i++) {
   5.225 +    /* determine number of points in entry */
   5.226 +    entry_npoints = entries[i].npoints;
   5.227 +    /* copy entry */
   5.228 +    cluster->entries[i].entrytype = entries[i].entrytype;
   5.229 +    cluster->entries[i].npoints = entry_npoints;
   5.230 +    /* calculate offset (in bytes) of pgl_point array */
   5.231 +    cluster->entries[i].offset = points_offset + npoints * sizeof(pgl_point);
   5.232 +    /* copy points */
   5.233 +    memcpy(
   5.234 +      PGL_ENTRY_POINTS(cluster, i),
   5.235 +      entries[i].points,
   5.236 +      entry_npoints * sizeof(pgl_point)
   5.237 +    );
   5.238 +    /* update total number of points processed */
   5.239 +    npoints += entry_npoints;
   5.240 +  }
   5.241 +  /* set number of entries in cluster */
   5.242 +  cluster->nentries = nentries;
   5.243 +  /* set PostgreSQL header for variable sized data */
   5.244 +  SET_VARSIZE(cluster, points_offset + npoints * sizeof(pgl_point));
   5.245 +  /* return newly created cluster */
   5.246 +  return cluster;
   5.247 +}
   5.248 +
   5.249 +
   5.250 +/*----------------------------------------*
   5.251 + *  C functions on geographic data types  *
   5.252 + *----------------------------------------*/
   5.253 +
   5.254 +/* round latitude or longitude to 12 digits after decimal point */
   5.255 +static inline double pgl_round(double val) {
   5.256 +  return round(val * 1e12) / 1e12;
   5.257 +}
   5.258 +
   5.259 +/* compare two points */
   5.260 +/* (equality when same point on earth is described, otherwise an arbitrary
   5.261 +   linear order) */
   5.262 +static int pgl_point_cmp(pgl_point *point1, pgl_point *point2) {
   5.263 +  double lon1, lon2;  /* modified longitudes for special cases */
   5.264 +  /* use latitude as first ordering criterion */
   5.265 +  if (point1->lat < point2->lat) return -1;
   5.266 +  if (point1->lat > point2->lat) return 1;
   5.267 +  /* determine modified longitudes (considering special case of poles and
   5.268 +     180th meridian which can be described as W180 or E180) */
   5.269 +  if (point1->lat == -90 || point1->lat == 90) lon1 = 0;
   5.270 +  else if (point1->lon == 180) lon1 = -180;
   5.271 +  else lon1 = point1->lon;
   5.272 +  if (point2->lat == -90 || point2->lat == 90) lon2 = 0;
   5.273 +  else if (point2->lon == 180) lon2 = -180;
   5.274 +  else lon2 = point2->lon;
   5.275 +  /* use (modified) longitude as secondary ordering criterion */
   5.276 +  if (lon1 < lon2) return -1;
   5.277 +  if (lon1 > lon2) return 1;
   5.278 +  /* no difference found, points are equal */
   5.279 +  return 0;
   5.280 +}
   5.281 +
   5.282 +/* compare two boxes */
   5.283 +/* (equality when same box on earth is described, otherwise an arbitrary linear
   5.284 +   order) */
   5.285 +static int pgl_box_cmp(pgl_box *box1, pgl_box *box2) {
   5.286 +  /* two empty boxes are equal, and an empty box is always considered "less
   5.287 +     than" a non-empty box */
   5.288 +  if (box1->lat_min> box1->lat_max && box2->lat_min<=box2->lat_max) return -1;
   5.289 +  if (box1->lat_min> box1->lat_max && box2->lat_min> box2->lat_max) return 0;
   5.290 +  if (box1->lat_min<=box1->lat_max && box2->lat_min> box2->lat_max) return 1;
   5.291 +  /* use southern border as first ordering criterion */
   5.292 +  if (box1->lat_min < box2->lat_min) return -1;
   5.293 +  if (box1->lat_min > box2->lat_min) return 1;
   5.294 +  /* use northern border as second ordering criterion */
   5.295 +  if (box1->lat_max < box2->lat_max) return -1;
   5.296 +  if (box1->lat_max > box2->lat_max) return 1;
   5.297 +  /* use western border as third ordering criterion */
   5.298 +  if (box1->lon_min < box2->lon_min) return -1;
   5.299 +  if (box1->lon_min > box2->lon_min) return 1;
   5.300 +  /* use eastern border as fourth ordering criterion */
   5.301 +  if (box1->lon_max < box2->lon_max) return -1;
   5.302 +  if (box1->lon_max > box2->lon_max) return 1;
   5.303 +  /* no difference found, boxes are equal */
   5.304 +  return 0;
   5.305 +}
   5.306 +
   5.307 +/* compare two circles */
   5.308 +/* (equality when same circle on earth is described, otherwise an arbitrary
   5.309 +   linear order) */
   5.310 +static int pgl_circle_cmp(pgl_circle *circle1, pgl_circle *circle2) {
   5.311 +  /* two circles with same infinite radius (positive or negative infinity) are
   5.312 +     considered equal independently of center point */
   5.313 +  if (
   5.314 +    !isfinite(circle1->radius) && !isfinite(circle2->radius) &&
   5.315 +    circle1->radius == circle2->radius
   5.316 +  ) return 0;
   5.317 +  /* use radius as first ordering criterion */
   5.318 +  if (circle1->radius < circle2->radius) return -1;
   5.319 +  if (circle1->radius > circle2->radius) return 1;
   5.320 +  /* use center point as secondary ordering criterion */
   5.321 +  return pgl_point_cmp(&(circle1->center), &(circle2->center));
   5.322 +}
   5.323 +
   5.324 +/* set box to empty box*/
   5.325 +static void pgl_box_set_empty(pgl_box *box) {
   5.326 +  box->lat_min = INFINITY;
   5.327 +  box->lat_max = -INFINITY;
   5.328 +  box->lon_min = 0;
   5.329 +  box->lon_max = 0;
   5.330 +}
   5.331 +
   5.332 +/* check if point is inside a box */
   5.333 +static bool pgl_point_in_box(pgl_point *point, pgl_box *box) {
   5.334 +  return (
   5.335 +    point->lat >= box->lat_min && point->lat <= box->lat_max && (
   5.336 +      (box->lon_min > box->lon_max) ? (
   5.337 +        /* box crosses 180th meridian */
   5.338 +        point->lon >= box->lon_min || point->lon <= box->lon_max
   5.339 +      ) : (
   5.340 +        /* box does not cross the 180th meridian */
   5.341 +        point->lon >= box->lon_min && point->lon <= box->lon_max
   5.342 +      )
   5.343 +    )
   5.344 +  );
   5.345 +}
   5.346 +
   5.347 +/* check if two boxes overlap */
   5.348 +static bool pgl_boxes_overlap(pgl_box *box1, pgl_box *box2) {
   5.349 +  return (
   5.350 +    box2->lat_max >= box2->lat_min &&  /* ensure box2 is not empty */
   5.351 +    ( box2->lat_min >= box1->lat_min || box2->lat_max >= box1->lat_min ) &&
   5.352 +    ( box2->lat_min <= box1->lat_max || box2->lat_max <= box1->lat_max ) && (
   5.353 +      (
   5.354 +        /* check if one and only one box crosses the 180th meridian */
   5.355 +        ((box1->lon_min > box1->lon_max) ? 1 : 0) ^
   5.356 +        ((box2->lon_min > box2->lon_max) ? 1 : 0)
   5.357 +      ) ? (
   5.358 +        /* exactly one box crosses the 180th meridian */
   5.359 +        box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min ||
   5.360 +        box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max
   5.361 +      ) : (
   5.362 +        /* no box or both boxes cross the 180th meridian */
   5.363 +        (
   5.364 +          (box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min) &&
   5.365 +          (box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max)
   5.366 +        ) ||
   5.367 +        /* handle W180 == E180 */
   5.368 +        ( box1->lon_min == -180 && box2->lon_max == 180 ) ||
   5.369 +        ( box2->lon_min == -180 && box1->lon_max == 180 )
   5.370 +      )
   5.371 +    )
   5.372 +  );
   5.373 +}
   5.374 +
   5.375 +/* check unambiguousness of east/west orientation of cluster entries and set
   5.376 +   bounding circle of cluster */
   5.377 +static bool pgl_finalize_cluster(pgl_cluster *cluster) {
   5.378 +  int i, j;                 /* i: index of entry, j: index of point in entry */
   5.379 +  int npoints;              /* number of points in entry */
   5.380 +  int total_npoints = 0;    /* total number of points in cluster */
   5.381 +  pgl_point *points;        /* points in entry */
   5.382 +  int lon_dir;              /* first point of entry west (-1) or east (+1) */
   5.383 +  double lon_break = 0;     /* antipodal longitude of first point in entry */
   5.384 +  double lon_min, lon_max;  /* covered longitude range of entry */
   5.385 +  double value;             /* temporary variable */
   5.386 +  /* reset bounding circle center to empty circle at 0/0 coordinates */
   5.387 +  cluster->bounding.center.lat = 0;
   5.388 +  cluster->bounding.center.lon = 0;
   5.389 +  cluster->bounding.radius = -INFINITY;
   5.390 +  /* if cluster is not empty */
   5.391 +  if (cluster->nentries != 0) {
   5.392 +    /* iterate over all cluster entries and ensure they each cover a longitude
   5.393 +       range less than 180 degrees */
   5.394 +    for (i=0; i<cluster->nentries; i++) {
   5.395 +      /* get properties of entry */
   5.396 +      npoints = cluster->entries[i].npoints;
   5.397 +      points = PGL_ENTRY_POINTS(cluster, i);
   5.398 +      /* get longitude of first point of entry */
   5.399 +      value = points[0].lon;
   5.400 +      /* initialize lon_min and lon_max with longitude of first point */
   5.401 +      lon_min = value;
   5.402 +      lon_max = value;
   5.403 +      /* determine east/west orientation of first point and calculate antipodal
   5.404 +         longitude (Note: rounding required here) */
   5.405 +      if      (value < 0) { lon_dir = -1; lon_break = pgl_round(value + 180); }
   5.406 +      else if (value > 0) { lon_dir =  1; lon_break = pgl_round(value - 180); }
   5.407 +      else lon_dir = 0;
   5.408 +      /* iterate over all other points in entry */
   5.409 +      for (j=1; j<npoints; j++) {
   5.410 +        /* consider longitude wrap-around */
   5.411 +        value = points[j].lon;
   5.412 +        if      (lon_dir<0 && value>lon_break) value = pgl_round(value - 360);
   5.413 +        else if (lon_dir>0 && value<lon_break) value = pgl_round(value + 360);
   5.414 +        /* update lon_min and lon_max */
   5.415 +        if      (value < lon_min) lon_min = value;
   5.416 +        else if (value > lon_max) lon_max = value;
   5.417 +        /* return false if 180 degrees or more are covered */
   5.418 +        if (lon_max - lon_min >= 180) return false;
   5.419 +      }
   5.420 +    }
   5.421 +    /* iterate over all points of all entries and calculate arbitrary center
   5.422 +       point for bounding circle (best if center point minimizes the radius,
   5.423 +       but some error is allowed here) */
   5.424 +    for (i=0; i<cluster->nentries; i++) {
   5.425 +      /* get properties of entry */
   5.426 +      npoints = cluster->entries[i].npoints;
   5.427 +      points = PGL_ENTRY_POINTS(cluster, i);
   5.428 +      /* check if first entry */
   5.429 +      if (i==0) {
   5.430 +        /* get longitude of first point of first entry in whole cluster */
   5.431 +        value = points[0].lon;
   5.432 +        /* initialize lon_min and lon_max with longitude of first point of
   5.433 +           first entry in whole cluster (used to determine if whole cluster
   5.434 +           covers a longitude range of 180 degrees or more) */
   5.435 +        lon_min = value;
   5.436 +        lon_max = value;
   5.437 +        /* determine east/west orientation of first point and calculate
   5.438 +           antipodal longitude (Note: rounding not necessary here) */
   5.439 +        if      (value < 0) { lon_dir = -1; lon_break = value + 180; }
   5.440 +        else if (value > 0) { lon_dir =  1; lon_break = value - 180; }
   5.441 +        else lon_dir = 0;
   5.442 +      }
   5.443 +      /* iterate over all points in entry */
   5.444 +      for (j=0; j<npoints; j++) {
   5.445 +        /* longitude wrap-around (Note: rounding not necessary here) */
   5.446 +        value = points[j].lon;
   5.447 +        if      (lon_dir < 0 && value > lon_break) value -= 360;
   5.448 +        else if (lon_dir > 0 && value < lon_break) value += 360;
   5.449 +        if      (value < lon_min) lon_min = value;
   5.450 +        else if (value > lon_max) lon_max = value;
   5.451 +        /* set bounding circle to cover whole earth if more than 180 degrees
   5.452 +           are covered */
   5.453 +        if (lon_max - lon_min >= 180) {
   5.454 +          cluster->bounding.center.lat = 0;
   5.455 +          cluster->bounding.center.lon = 0;
   5.456 +          cluster->bounding.radius = INFINITY;
   5.457 +          return true;
   5.458 +        }
   5.459 +        /* add point to bounding circle center (for average calculation) */
   5.460 +        cluster->bounding.center.lat += points[j].lat;
   5.461 +        cluster->bounding.center.lon += value;
   5.462 +      }
   5.463 +      /* count total number of points */
   5.464 +      total_npoints += npoints;
   5.465 +    }
   5.466 +    /* determine average latitude and longitude of cluster */
   5.467 +    cluster->bounding.center.lat /= total_npoints;
   5.468 +    cluster->bounding.center.lon /= total_npoints;
   5.469 +    /* normalize longitude of center of cluster bounding circle */
   5.470 +    if (cluster->bounding.center.lon < -180) {
   5.471 +      cluster->bounding.center.lon += 360;
   5.472 +    }
   5.473 +    else if (cluster->bounding.center.lon > 180) {
   5.474 +      cluster->bounding.center.lon -= 360;
   5.475 +    }
   5.476 +    /* round bounding circle center (useful if it is used by other functions) */
   5.477 +    cluster->bounding.center.lat = pgl_round(cluster->bounding.center.lat);
   5.478 +    cluster->bounding.center.lon = pgl_round(cluster->bounding.center.lon);
   5.479 +    /* calculate radius of bounding circle */
   5.480 +    for (i=0; i<cluster->nentries; i++) {
   5.481 +      npoints = cluster->entries[i].npoints;
   5.482 +      points = PGL_ENTRY_POINTS(cluster, i);
   5.483 +      for (j=0; j<npoints; j++) {
   5.484 +        value = pgl_distance(
   5.485 +          cluster->bounding.center.lat, cluster->bounding.center.lon,
   5.486 +          points[j].lat, points[j].lon
   5.487 +        );
   5.488 +        if (value > cluster->bounding.radius) cluster->bounding.radius = value;
   5.489 +      }
   5.490 +    }
   5.491 +  }
   5.492 +  /* return true (east/west orientation is unambiguous) */
   5.493 +  return true;
   5.494 +}
   5.495 +
   5.496 +/* check if point is inside cluster */
   5.497 +/* (if point is on perimeter, then true is returned if and only if
   5.498 +   strict == false) */
   5.499 +static bool pgl_point_in_cluster(
   5.500 +  pgl_point *point,
   5.501 +  pgl_cluster *cluster,
   5.502 +  bool strict
   5.503 +) {
   5.504 +  int i, j, k;  /* i: entry, j: point in entry, k: next point in entry */
   5.505 +  int entrytype;         /* type of entry */
   5.506 +  int npoints;           /* number of points in entry */
   5.507 +  pgl_point *points;     /* array of points in entry */
   5.508 +  int lon_dir = 0;       /* first vertex west (-1) or east (+1) */
   5.509 +  double lon_break = 0;  /* antipodal longitude of first vertex */
   5.510 +  double lat0 = point->lat;  /* latitude of point */
   5.511 +  double lon0;           /* (adjusted) longitude of point */
   5.512 +  double lat1, lon1;     /* latitude and (adjusted) longitude of vertex */
   5.513 +  double lat2, lon2;     /* latitude and (adjusted) longitude of next vertex */
   5.514 +  double lon;            /* longitude of intersection */
   5.515 +  int counter = 0;       /* counter for intersections east of point */
   5.516 +  /* iterate over all entries */
   5.517 +  for (i=0; i<cluster->nentries; i++) {
   5.518 +    /* get type of entry */
   5.519 +    entrytype = cluster->entries[i].entrytype;
   5.520 +    /* skip all entries but polygons if perimeters are excluded */
   5.521 +    if (strict && entrytype != PGL_ENTRY_POLYGON) continue;
   5.522 +    /* get points of entry */
   5.523 +    npoints = cluster->entries[i].npoints;
   5.524 +    points = PGL_ENTRY_POINTS(cluster, i);
   5.525 +    /* determine east/west orientation of first point of entry and calculate
   5.526 +       antipodal longitude */
   5.527 +    lon_break = points[0].lon;
   5.528 +    if      (lon_break < 0) { lon_dir = -1; lon_break += 180; }
   5.529 +    else if (lon_break > 0) { lon_dir =  1; lon_break -= 180; }
   5.530 +    else lon_dir = 0;
   5.531 +    /* get longitude of point */
   5.532 +    lon0 = point->lon;
   5.533 +    /* consider longitude wrap-around for point */
   5.534 +    if      (lon_dir < 0 && lon0 > lon_break) lon0 = pgl_round(lon0 - 360);
   5.535 +    else if (lon_dir > 0 && lon0 < lon_break) lon0 = pgl_round(lon0 + 360);
   5.536 +    /* iterate over all edges and vertices */
   5.537 +    for (j=0; j<npoints; j++) {
   5.538 +      /* return if point is on vertex of polygon */
   5.539 +      if (pgl_point_cmp(point, &(points[j])) == 0) return !strict;
   5.540 +      /* calculate index of next vertex */
   5.541 +      k = (j+1) % npoints;
   5.542 +      /* skip last edge unless entry is (closed) outline or polygon */
   5.543 +      if (
   5.544 +        k == 0 &&
   5.545 +        entrytype != PGL_ENTRY_OUTLINE &&
   5.546 +        entrytype != PGL_ENTRY_POLYGON
   5.547 +      ) continue;
   5.548 +      /* use previously calculated values for lat1 and lon1 if possible */
   5.549 +      if (j) {
   5.550 +        lat1 = lat2;
   5.551 +        lon1 = lon2;
   5.552 +      } else {
   5.553 +        /* otherwise get latitude and longitude values of first vertex */
   5.554 +        lat1 = points[0].lat;
   5.555 +        lon1 = points[0].lon;
   5.556 +        /* and consider longitude wrap-around for first vertex */
   5.557 +        if      (lon_dir < 0 && lon1 > lon_break) lon1 = pgl_round(lon1 - 360);
   5.558 +        else if (lon_dir > 0 && lon1 < lon_break) lon1 = pgl_round(lon1 + 360);
   5.559 +      }
   5.560 +      /* get latitude and longitude of next vertex */
   5.561 +      lat2 = points[k].lat;
   5.562 +      lon2 = points[k].lon;
   5.563 +      /* consider longitude wrap-around for next vertex */
   5.564 +      if      (lon_dir < 0 && lon2 > lon_break) lon2 = pgl_round(lon2 - 360);
   5.565 +      else if (lon_dir > 0 && lon2 < lon_break) lon2 = pgl_round(lon2 + 360);
   5.566 +      /* return if point is on horizontal (west to east) edge of polygon */
   5.567 +      if (
   5.568 +        lat0 == lat1 && lat0 == lat2 &&
   5.569 +        ( (lon0 >= lon1 && lon0 <= lon2) || (lon0 >= lon2 && lon0 <= lon1) )
   5.570 +      ) return !strict;
   5.571 +      /* check if edge crosses east/west line of point */
   5.572 +      if ((lat1 < lat0 && lat2 >= lat0) || (lat2 < lat0 && lat1 >= lat0)) {
   5.573 +        /* calculate longitude of intersection */
   5.574 +        lon = (lon1 * (lat2-lat0) + lon2 * (lat0-lat1)) / (lat2-lat1);
   5.575 +        /* return if intersection goes (approximately) through point */
   5.576 +        if (pgl_round(lon) == lon0) return !strict;
   5.577 +        /* count intersection if east of point and entry is polygon*/
   5.578 +        if (entrytype == PGL_ENTRY_POLYGON && lon > lon0) counter++;
   5.579 +      }
   5.580 +    }
   5.581 +  }
   5.582 +  /* return true if number of intersections is odd */
   5.583 +  return counter & 1;
   5.584 +}
   5.585 +
   5.586 +/* check if all points of the second cluster are strictly inside the first
   5.587 +   cluster */
   5.588 +static inline bool pgl_all_cluster_points_strictly_in_cluster(
   5.589 +  pgl_cluster *outer, pgl_cluster *inner
   5.590 +) {
   5.591 +  int i, j;           /* i: entry, j: point in entry */
   5.592 +  int npoints;        /* number of points in entry */
   5.593 +  pgl_point *points;  /* array of points in entry */
   5.594 +  /* iterate over all entries of "inner" cluster */
   5.595 +  for (i=0; i<inner->nentries; i++) {
   5.596 +    /* get properties of entry */
   5.597 +    npoints = inner->entries[i].npoints;
   5.598 +    points = PGL_ENTRY_POINTS(inner, i);
   5.599 +    /* iterate over all points in entry of "inner" cluster */
   5.600 +    for (j=0; j<npoints; j++) {
   5.601 +      /* return false if one point of inner cluster is not in outer cluster */
   5.602 +      if (!pgl_point_in_cluster(points+j, outer, true)) return false;
   5.603 +    }
   5.604 +  }
   5.605 +  /* otherwise return true */
   5.606 +  return true;
   5.607 +}
   5.608 +
   5.609 +/* check if any point the second cluster is inside the first cluster */
   5.610 +static inline bool pgl_any_cluster_points_in_cluster(
   5.611 +  pgl_cluster *outer, pgl_cluster *inner
   5.612 +) {
   5.613 +  int i, j;           /* i: entry, j: point in entry */
   5.614 +  int npoints;        /* number of points in entry */
   5.615 +  pgl_point *points;  /* array of points in entry */
   5.616 +  /* iterate over all entries of "inner" cluster */
   5.617 +  for (i=0; i<inner->nentries; i++) {
   5.618 +    /* get properties of entry */
   5.619 +    npoints = inner->entries[i].npoints;
   5.620 +    points = PGL_ENTRY_POINTS(inner, i);
   5.621 +    /* iterate over all points in entry of "inner" cluster */
   5.622 +    for (j=0; j<npoints; j++) {
   5.623 +      /* return true if one point of inner cluster is in outer cluster */
   5.624 +      if (pgl_point_in_cluster(points+j, outer, false)) return true;
   5.625 +    }
   5.626 +  }
   5.627 +  /* otherwise return false */
   5.628 +  return false;
   5.629 +}
   5.630 +
   5.631 +/* check if line segment strictly crosses line (not just touching) */
   5.632 +static inline bool pgl_lseg_crosses_line(
   5.633 +  double seg_x1,  double seg_y1,  double seg_x2,  double seg_y2,
   5.634 +  double line_x1, double line_y1, double line_x2, double line_y2
   5.635 +) {
   5.636 +  return (
   5.637 +    (
   5.638 +      (seg_x1-line_x1) * (line_y2-line_y1) -
   5.639 +      (seg_y1-line_y1) * (line_x2-line_x1)
   5.640 +    ) * (
   5.641 +      (seg_x2-line_x1) * (line_y2-line_y1) -
   5.642 +      (seg_y2-line_y1) * (line_x2-line_x1)
   5.643 +    )
   5.644 +  ) < 0;
   5.645 +}
   5.646 +
   5.647 +/* check if paths and outlines of two clusters strictly overlap (not just
   5.648 +   touching) */
   5.649 +static bool pgl_outlines_overlap(
   5.650 +  pgl_cluster *cluster1, pgl_cluster *cluster2
   5.651 +) {
   5.652 +  int i1, j1, k1;  /* i: entry, j: point in entry, k: next point in entry */
   5.653 +  int i2, j2, k2;
   5.654 +  int entrytype1, entrytype2;     /* type of entry */
   5.655 +  int npoints1, npoints2;         /* number of points in entry */
   5.656 +  pgl_point *points1;             /* array of points in entry of cluster1 */
   5.657 +  pgl_point *points2;             /* array of points in entry of cluster2 */
   5.658 +  int lon_dir1, lon_dir2;         /* first vertex west (-1) or east (+1) */
   5.659 +  double lon_break1, lon_break2;  /* antipodal longitude of first vertex */
   5.660 +  double lat11, lon11;  /* latitude and (adjusted) longitude of vertex */
   5.661 +  double lat12, lon12;  /* latitude and (adjusted) longitude of next vertex */
   5.662 +  double lat21, lon21;  /* latitude and (adjusted) longitudes for cluster2 */
   5.663 +  double lat22, lon22;
   5.664 +  double wrapvalue;     /* temporary helper value to adjust wrap-around */  
   5.665 +  /* iterate over all entries of cluster1 */
   5.666 +  for (i1=0; i1<cluster1->nentries; i1++) {
   5.667 +    /* get properties of entry in cluster1 and skip points */
   5.668 +    npoints1 = cluster1->entries[i1].npoints;
   5.669 +    if (npoints1 < 2) continue;
   5.670 +    entrytype1 = cluster1->entries[i1].entrytype;
   5.671 +    points1 = PGL_ENTRY_POINTS(cluster1, i1);
   5.672 +    /* determine east/west orientation of first point and calculate antipodal
   5.673 +       longitude */
   5.674 +    lon_break1 = points1[0].lon;
   5.675 +    if (lon_break1 < 0) {
   5.676 +      lon_dir1   = -1;
   5.677 +      lon_break1 = pgl_round(lon_break1 + 180);
   5.678 +    } else if (lon_break1 > 0) {
   5.679 +      lon_dir1   = 1;
   5.680 +      lon_break1 = pgl_round(lon_break1 - 180);
   5.681 +    } else lon_dir1 = 0;
   5.682 +    /* iterate over all edges and vertices in cluster1 */
   5.683 +    for (j1=0; j1<npoints1; j1++) {
   5.684 +      /* calculate index of next vertex */
   5.685 +      k1 = (j1+1) % npoints1;
   5.686 +      /* skip last edge unless entry is (closed) outline or polygon */
   5.687 +      if (
   5.688 +        k1 == 0 &&
   5.689 +        entrytype1 != PGL_ENTRY_OUTLINE &&
   5.690 +        entrytype1 != PGL_ENTRY_POLYGON
   5.691 +      ) continue;
   5.692 +      /* use previously calculated values for lat1 and lon1 if possible */
   5.693 +      if (j1) {
   5.694 +        lat11 = lat12;
   5.695 +        lon11 = lon12;
   5.696 +      } else {
   5.697 +        /* otherwise get latitude and longitude values of first vertex */
   5.698 +        lat11 = points1[0].lat;
   5.699 +        lon11 = points1[0].lon;
   5.700 +        /* and consider longitude wrap-around for first vertex */
   5.701 +        if      (lon_dir1<0 && lon11>lon_break1) lon11 = pgl_round(lon11-360);
   5.702 +        else if (lon_dir1>0 && lon11<lon_break1) lon11 = pgl_round(lon11+360);
   5.703 +      }
   5.704 +      /* get latitude and longitude of next vertex */
   5.705 +      lat12 = points1[k1].lat;
   5.706 +      lon12 = points1[k1].lon;
   5.707 +      /* consider longitude wrap-around for next vertex */
   5.708 +      if      (lon_dir1<0 && lon12>lon_break1) lon12 = pgl_round(lon12-360);
   5.709 +      else if (lon_dir1>0 && lon12<lon_break1) lon12 = pgl_round(lon12+360);
   5.710 +      /* skip degenerated edges */
   5.711 +      if (lat11 == lat12 && lon11 == lon12) continue;
   5.712 +      /* iterate over all entries of cluster2 */
   5.713 +      for (i2=0; i2<cluster2->nentries; i2++) {
   5.714 +        /* get points and number of points of entry in cluster2 */
   5.715 +        npoints2 = cluster2->entries[i2].npoints;
   5.716 +        if (npoints2 < 2) continue;
   5.717 +        entrytype2 = cluster2->entries[i2].entrytype;
   5.718 +        points2 = PGL_ENTRY_POINTS(cluster2, i2);
   5.719 +        /* determine east/west orientation of first point and calculate antipodal
   5.720 +           longitude */
   5.721 +        lon_break2 = points2[0].lon;
   5.722 +        if (lon_break2 < 0) {
   5.723 +          lon_dir2   = -1;
   5.724 +          lon_break2 = pgl_round(lon_break2 + 180);
   5.725 +        } else if (lon_break2 > 0) {
   5.726 +          lon_dir2   = 1;
   5.727 +          lon_break2 = pgl_round(lon_break2 - 180);
   5.728 +        } else lon_dir2 = 0;
   5.729 +        /* iterate over all edges and vertices in cluster2 */
   5.730 +        for (j2=0; j2<npoints2; j2++) {
   5.731 +          /* calculate index of next vertex */
   5.732 +          k2 = (j2+1) % npoints2;
   5.733 +          /* skip last edge unless entry is (closed) outline or polygon */
   5.734 +          if (
   5.735 +            k2 == 0 &&
   5.736 +            entrytype2 != PGL_ENTRY_OUTLINE &&
   5.737 +            entrytype2 != PGL_ENTRY_POLYGON
   5.738 +          ) continue;
   5.739 +          /* use previously calculated values for lat1 and lon1 if possible */
   5.740 +          if (j2) {
   5.741 +            lat21 = lat22;
   5.742 +            lon21 = lon22;
   5.743 +          } else {
   5.744 +            /* otherwise get latitude and longitude values of first vertex */
   5.745 +            lat21 = points2[0].lat;
   5.746 +            lon21 = points2[0].lon;
   5.747 +            /* and consider longitude wrap-around for first vertex */
   5.748 +            if      (lon_dir2<0 && lon21>lon_break2) lon21 = pgl_round(lon21-360);
   5.749 +            else if (lon_dir2>0 && lon21<lon_break2) lon21 = pgl_round(lon21+360);
   5.750 +          }
   5.751 +          /* get latitude and longitude of next vertex */
   5.752 +          lat22 = points2[k2].lat;
   5.753 +          lon22 = points2[k2].lon;
   5.754 +          /* consider longitude wrap-around for next vertex */
   5.755 +          if      (lon_dir2<0 && lon22>lon_break2) lon22 = pgl_round(lon22-360);
   5.756 +          else if (lon_dir2>0 && lon22<lon_break2) lon22 = pgl_round(lon22+360);
   5.757 +          /* skip degenerated edges */
   5.758 +          if (lat21 == lat22 && lon21 == lon22) continue;
   5.759 +          /* perform another wrap-around where necessary */
   5.760 +          /* TODO: improve performance of whole wrap-around mechanism */
   5.761 +          wrapvalue = (lon21 + lon22) - (lon11 + lon12);
   5.762 +          if (wrapvalue > 360) {
   5.763 +            lon21 = pgl_round(lon21 - 360);
   5.764 +            lon22 = pgl_round(lon22 - 360);
   5.765 +          } else if (wrapvalue < -360) {
   5.766 +            lon21 = pgl_round(lon21 + 360);
   5.767 +            lon22 = pgl_round(lon22 + 360);
   5.768 +          }
   5.769 +          /* return true if segments overlap */
   5.770 +          if (
   5.771 +            pgl_lseg_crosses_line(
   5.772 +              lat11, lon11, lat12, lon12,
   5.773 +              lat21, lon21, lat22, lon22
   5.774 +            ) && pgl_lseg_crosses_line(
   5.775 +              lat21, lon21, lat22, lon22,
   5.776 +              lat11, lon11, lat12, lon12
   5.777 +            )
   5.778 +          ) {
   5.779 +            return true;
   5.780 +          }
   5.781 +        }
   5.782 +      }
   5.783 +    }
   5.784 +  }
   5.785 +  /* otherwise return false */
   5.786 +  return false;
   5.787 +}
   5.788 +
   5.789 +/* check if second cluster is completely contained in first cluster */
   5.790 +static bool pgl_cluster_in_cluster(pgl_cluster *outer, pgl_cluster *inner) {
   5.791 +  if (!pgl_all_cluster_points_strictly_in_cluster(outer, inner)) return false;
   5.792 +  if (pgl_any_cluster_points_in_cluster(inner, outer)) return false;
   5.793 +  if (pgl_outlines_overlap(outer, inner)) return false;
   5.794 +  return true;
   5.795 +}
   5.796 +
   5.797 +/* check if two clusters overlap */
   5.798 +static bool pgl_clusters_overlap(
   5.799 +  pgl_cluster *cluster1, pgl_cluster *cluster2
   5.800 +) {
   5.801 +  if (pgl_any_cluster_points_in_cluster(cluster1, cluster2)) return true;
   5.802 +  if (pgl_any_cluster_points_in_cluster(cluster2, cluster1)) return true;
   5.803 +  if (pgl_outlines_overlap(cluster1, cluster2)) return true;
   5.804 +  return false;
   5.805 +}
   5.806 +
   5.807 +
   5.808 +/* calculate (approximate) distance between point and cluster */
   5.809 +static double pgl_point_cluster_distance(pgl_point *point, pgl_cluster *cluster) {
   5.810 +  double comp;           /* square of compression of meridians */
   5.811 +  int i, j, k;  /* i: entry, j: point in entry, k: next point in entry */
   5.812 +  int entrytype;         /* type of entry */
   5.813 +  int npoints;           /* number of points in entry */
   5.814 +  pgl_point *points;     /* array of points in entry */
   5.815 +  int lon_dir = 0;       /* first vertex west (-1) or east (+1) */
   5.816 +  double lon_break = 0;  /* antipodal longitude of first vertex */
   5.817 +  double lon_min = 0;    /* minimum (adjusted) longitude of entry vertices */
   5.818 +  double lon_max = 0;    /* maximum (adjusted) longitude of entry vertices */
   5.819 +  double lat0 = point->lat;  /* latitude of point */
   5.820 +  double lon0;           /* (adjusted) longitude of point */
   5.821 +  double lat1, lon1;     /* latitude and (adjusted) longitude of vertex */
   5.822 +  double lat2, lon2;     /* latitude and (adjusted) longitude of next vertex */
   5.823 +  double s;              /* scalar for vector calculations */
   5.824 +  double dist;           /* distance calculated in one step */
   5.825 +  double min_dist = INFINITY;   /* minimum distance */
   5.826 +  /* distance is zero if point is contained in cluster */
   5.827 +  if (pgl_point_in_cluster(point, cluster, false)) return 0;
   5.828 +  /* calculate approximate square compression of meridians */
   5.829 +  comp = cos((lat0 / 180.0) * M_PI);
   5.830 +  comp *= comp;
   5.831 +  /* calculate exact square compression of meridians */
   5.832 +  comp *= (
   5.833 +    (1.0 - PGL_EPS2 * (1.0-comp)) *
   5.834 +    (1.0 - PGL_EPS2 * (1.0-comp)) /
   5.835 +    (PGL_SUBEPS2 * PGL_SUBEPS2)
   5.836 +  );
   5.837 +  /* iterate over all entries */
   5.838 +  for (i=0; i<cluster->nentries; i++) {
   5.839 +    /* get properties of entry */
   5.840 +    entrytype = cluster->entries[i].entrytype;
   5.841 +    npoints = cluster->entries[i].npoints;
   5.842 +    points = PGL_ENTRY_POINTS(cluster, i);
   5.843 +    /* determine east/west orientation of first point of entry and calculate
   5.844 +       antipodal longitude */
   5.845 +    lon_break = points[0].lon;
   5.846 +    if      (lon_break < 0) { lon_dir = -1; lon_break += 180; }
   5.847 +    else if (lon_break > 0) { lon_dir =  1; lon_break -= 180; }
   5.848 +    else lon_dir = 0;
   5.849 +    /* determine covered longitude range */
   5.850 +    for (j=0; j<npoints; j++) {
   5.851 +      /* get longitude of vertex */
   5.852 +      lon1 = points[j].lon;
   5.853 +      /* adjust longitude to fix potential wrap-around */
   5.854 +      if      (lon_dir < 0 && lon1 > lon_break) lon1 -= 360;
   5.855 +      else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
   5.856 +      /* update minimum and maximum longitude of polygon */
   5.857 +      if (j == 0 || lon1 < lon_min) lon_min = lon1;
   5.858 +      if (j == 0 || lon1 > lon_max) lon_max = lon1;
   5.859 +    }
   5.860 +    /* adjust longitude wrap-around according to full longitude range */
   5.861 +    lon_break = (lon_max + lon_min) / 2;
   5.862 +    if      (lon_break < 0) { lon_dir = -1; lon_break += 180; }
   5.863 +    else if (lon_break > 0) { lon_dir =  1; lon_break -= 180; }
   5.864 +    /* get longitude of point */
   5.865 +    lon0 = point->lon;
   5.866 +    /* consider longitude wrap-around for point */
   5.867 +    if      (lon_dir < 0 && lon0 > lon_break) lon0 -= 360;
   5.868 +    else if (lon_dir > 0 && lon0 < lon_break) lon0 += 360;
   5.869 +    /* iterate over all edges and vertices */
   5.870 +    for (j=0; j<npoints; j++) {
   5.871 +      /* use previously calculated values for lat1 and lon1 if possible */
   5.872 +      if (j) {
   5.873 +        lat1 = lat2;
   5.874 +        lon1 = lon2;
   5.875 +      } else {
   5.876 +        /* otherwise get latitude and longitude values of first vertex */
   5.877 +        lat1 = points[0].lat;
   5.878 +        lon1 = points[0].lon;
   5.879 +        /* and consider longitude wrap-around for first vertex */
   5.880 +        if      (lon_dir < 0 && lon1 > lon_break) lon1 -= 360;
   5.881 +        else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
   5.882 +      }
   5.883 +      /* calculate distance to vertex */
   5.884 +      dist = pgl_distance(lat0, lon0, lat1, lon1);
   5.885 +      /* store calculated distance if smallest */
   5.886 +      if (dist < min_dist) min_dist = dist;
   5.887 +      /* calculate index of next vertex */
   5.888 +      k = (j+1) % npoints;
   5.889 +      /* skip last edge unless entry is (closed) outline or polygon */
   5.890 +      if (
   5.891 +        k == 0 &&
   5.892 +        entrytype != PGL_ENTRY_OUTLINE &&
   5.893 +        entrytype != PGL_ENTRY_POLYGON
   5.894 +      ) continue;
   5.895 +      /* get latitude and longitude of next vertex */
   5.896 +      lat2 = points[k].lat;
   5.897 +      lon2 = points[k].lon;
   5.898 +      /* consider longitude wrap-around for next vertex */
   5.899 +      if      (lon_dir < 0 && lon2 > lon_break) lon2 -= 360;
   5.900 +      else if (lon_dir > 0 && lon2 < lon_break) lon2 += 360;
   5.901 +      /* go to next vertex and edge if edge is degenerated */
   5.902 +      if (lat1 == lat2 && lon1 == lon2) continue;
   5.903 +      /* otherwise test if point can be projected onto edge of polygon */
   5.904 +      s = (
   5.905 +        ((lat0-lat1) * (lat2-lat1) + comp * (lon0-lon1) * (lon2-lon1)) /
   5.906 +        ((lat2-lat1) * (lat2-lat1) + comp * (lon2-lon1) * (lon2-lon1))
   5.907 +      );
   5.908 +      /* go to next vertex and edge if point cannot be projected */
   5.909 +      if (!(s > 0 && s < 1)) continue;
   5.910 +      /* calculate distance from original point to projected point */
   5.911 +      dist = pgl_distance(
   5.912 +        lat0, lon0,
   5.913 +        lat1 + s * (lat2-lat1),
   5.914 +        lon1 + s * (lon2-lon1)
   5.915 +      );
   5.916 +      /* store calculated distance if smallest */
   5.917 +      if (dist < min_dist) min_dist = dist;
   5.918 +    }
   5.919 +  }
   5.920 +  /* return minimum distance */
   5.921 +  return min_dist;
   5.922 +}
   5.923 +
   5.924 +/* calculate (approximate) distance between two clusters */
   5.925 +static double pgl_cluster_distance(pgl_cluster *cluster1, pgl_cluster *cluster2) {
   5.926 +  int i, j;                    /* i: entry, j: point in entry */
   5.927 +  int npoints;                 /* number of points in entry */
   5.928 +  pgl_point *points;           /* array of points in entry */
   5.929 +  double dist;                 /* distance calculated in one step */
   5.930 +  double min_dist = INFINITY;  /* minimum distance */
   5.931 +  /* consider distance from each point in one cluster to the whole other */
   5.932 +  for (i=0; i<cluster1->nentries; i++) {
   5.933 +    npoints = cluster1->entries[i].npoints;
   5.934 +    points = PGL_ENTRY_POINTS(cluster1, i);
   5.935 +    for (j=0; j<npoints; j++) {
   5.936 +      dist = pgl_point_cluster_distance(points+j, cluster2);
   5.937 +      if (dist == 0) return dist;
   5.938 +      if (dist < min_dist) min_dist = dist;
   5.939 +    }
   5.940 +  }
   5.941 +  /* consider distance from each point in other cluster to the first cluster */
   5.942 +  for (i=0; i<cluster2->nentries; i++) {
   5.943 +    npoints = cluster2->entries[i].npoints;
   5.944 +    points = PGL_ENTRY_POINTS(cluster2, i);
   5.945 +    for (j=0; j<npoints; j++) {
   5.946 +      dist = pgl_point_cluster_distance(points+j, cluster1);
   5.947 +      if (dist == 0) return dist;
   5.948 +      if (dist < min_dist) min_dist = dist;
   5.949 +    }
   5.950 +  }
   5.951 +  return min_dist;
   5.952 +}
   5.953 +
   5.954 +/* estimator function for distance between box and point */
   5.955 +/* always returns a smaller value than actually correct or zero */
   5.956 +static double pgl_estimate_point_box_distance(pgl_point *point, pgl_box *box) {
   5.957 +  double dlon;      /* longitude range of box (delta longitude) */
   5.958 +  double distance;  /* return value */
   5.959 +  /* return infinity if box is empty */
   5.960 +  if (box->lat_min > box->lat_max) return INFINITY;
   5.961 +  /* return zero if point is inside box */
   5.962 +  if (pgl_point_in_box(point, box)) return 0;
   5.963 +  /* calculate delta longitude */
   5.964 +  dlon = box->lon_max - box->lon_min;
   5.965 +  if (dlon < 0) dlon += 360;  /* 180th meridian crossed */
   5.966 +  /* if delta longitude is greater than 150 degrees, perform safe fall-back */
   5.967 +  if (dlon > 150) return 0;
   5.968 +  /* calculate lower limit for distance (formula below requires dlon <= 150) */
   5.969 +  /* TODO: provide better estimation function to improve performance */
   5.970 +  distance = (
   5.971 +    (1.0-1e-14) * pgl_distance(
   5.972 +      point->lat,
   5.973 +      point->lon,
   5.974 +      (box->lat_min + box->lat_max) / 2,
   5.975 +      box->lon_min + dlon/2
   5.976 +    ) - pgl_distance(
   5.977 +      box->lat_min, box->lon_min,
   5.978 +      box->lat_max, box->lon_max
   5.979 +    )
   5.980 +  );
   5.981 +  /* truncate negative results to zero */
   5.982 +  if (distance <= 0) distance = 0;
   5.983 +  /* return result */
   5.984 +  return distance;
   5.985 +}
   5.986 +
   5.987 +
   5.988 +/*-------------------------------*
   5.989 + *  Monte Carlo based functions  *
   5.990 + *-------------------------------*/
   5.991 +
   5.992 +/* half of (spherical) earth's surface area */
   5.993 +#define PGL_HALF_SURFACE (PGL_RADIUS * PGL_DIAMETER * M_PI)
   5.994 +
   5.995 +/* golden angle in radians */
   5.996 +#define PGL_GOLDEN_ANGLE (M_PI * (sqrt(5) - 1.0))
   5.997 +
   5.998 +/* create a list of sample points covering a bounding circle 
   5.999 +   and return covered area */
  5.1000 +static double pgl_sample_points(
  5.1001 +  pgl_point *center,  /* center of bounding circle */
  5.1002 +  double radius,      /* radius of bounding circle */
  5.1003 +  int samples,        /* number of sample points */
  5.1004 +  pgl_point *result   /* pointer to result array */
  5.1005 +) {
  5.1006 +  double double_share = 2.0;  /* double of covered share of earth's surface */
  5.1007 +  double double_share_div_samples;  /* double_share divided by sample count */
  5.1008 +  int i;
  5.1009 +  double t;  /* parameter of spiral laid on (spherical) earth's surface */
  5.1010 +  double x, y, z;  /* normalized coordinates of point on non-rotated spiral */
  5.1011 +  double sin_phi;  /* sine of sph. coordinate of point of non-rotated spiral */
  5.1012 +  double lambda;   /* other sph. coordinate of point of non-rotated spiral */
  5.1013 +  double rot = (0.5 - center->lat / 180.0) * M_PI;  /* needed rot. (in rad) */
  5.1014 +  double cos_rot = cos(rot);  /* cosine of rotation by latitude */
  5.1015 +  double sin_rot = sin(rot);  /* sine of rotation by latitude */
  5.1016 +  double x_rot, z_rot;  /* normalized coordinates of point on rotated spiral */
  5.1017 +  double center_lon = center->lon;  /* second rotation in degree */
  5.1018 +  /* add safety margin to bounding circle because of spherical approximation */
  5.1019 +  radius *= PGL_SPHEROID_A / PGL_RADIUS;
  5.1020 +  /* if whole earth is covered, use initialized value, otherwise calculate
  5.1021 +     share of covered area (multiplied by 2) */
  5.1022 +  if (radius < PGL_MAXDIST) double_share = 1.0 - cos(radius / PGL_RADIUS);
  5.1023 +  /* divide double_share by sample count for later calculations */
  5.1024 +  double_share_div_samples = double_share / samples;
  5.1025 +  /* generate sample points */
  5.1026 +  for (i=0; i<samples; i++) {
  5.1027 +    /* use an offset of 1/2 to avoid too dense clustering at spiral center */
  5.1028 +    t = 0.5 + i;
  5.1029 +    /* calculate normalized coordinates of point on non-rotated spiral */
  5.1030 +    z = 1.0 - double_share_div_samples * t;
  5.1031 +    sin_phi = sqrt(1.0 - z*z);
  5.1032 +    lambda = t * PGL_GOLDEN_ANGLE;
  5.1033 +    x = sin_phi * cos(lambda);
  5.1034 +    y = sin_phi * sin(lambda);
  5.1035 +    /* rotate spiral by latitude value of bounding circle */
  5.1036 +    x_rot = cos_rot * x + sin_rot * z;
  5.1037 +    z_rot = cos_rot * z - sin_rot * x;
  5.1038 +    /* set resulting sample point in result array */
  5.1039 +    /* (while performing second rotation by bounding circle longitude) */
  5.1040 +    result[i].lat = 180.0 * (atan(z_rot / fabs(x_rot)) / M_PI);
  5.1041 +    result[i].lon = center_lon + 180.0 * (atan2(y, x_rot) / M_PI);
  5.1042 +  }
  5.1043 +  /* return covered area */
  5.1044 +  return PGL_HALF_SURFACE * double_share;
  5.1045 +}
  5.1046 +
  5.1047 +/* calculate covered area using Monte Carlo method */
  5.1048 +/* TODO: inefficient, should be replaced by different method */
  5.1049 +static double pgl_monte_carlo_area(pgl_cluster *cluster, int samples) {
  5.1050 +  pgl_point *points;
  5.1051 +  double area;
  5.1052 +  int i;
  5.1053 +  int matches = 0;
  5.1054 +  if (samples > PGL_CLUSTER_MAXPOINTS) {
  5.1055 +    ereport(ERROR, (
  5.1056 +      errcode(ERRCODE_DATA_EXCEPTION),
  5.1057 +      errmsg(
  5.1058 +        "too many sample points for Monte Carlo method (maximum %i)",
  5.1059 +        PGL_CLUSTER_MAXPOINTS
  5.1060 +      )
  5.1061 +    ));
  5.1062 +  }
  5.1063 +  if (!(cluster->bounding.radius > 0)) return 0;
  5.1064 +  for (i=0; ; i++) {
  5.1065 +    if (i == cluster->nentries) return 0;
  5.1066 +    if (cluster->entries[i].entrytype == PGL_ENTRY_POLYGON) break;
  5.1067 +  }
  5.1068 +  points = palloc(samples * sizeof(pgl_point));
  5.1069 +  area = pgl_sample_points(
  5.1070 +    &cluster->bounding.center,
  5.1071 +    cluster->bounding.radius,
  5.1072 +    samples,
  5.1073 +    points
  5.1074 +  );
  5.1075 +  for (i=0; i<samples; i++) {
  5.1076 +    if (pgl_point_in_cluster(points+i, cluster, true)) matches++;
  5.1077 +  }
  5.1078 +  pfree(points);
  5.1079 +  return area * ((double)matches / samples);
  5.1080 +}
  5.1081 +
  5.1082 +/* fair distance between point and cluster (see README file for explanation) */
  5.1083 +static double pgl_fair_distance(
  5.1084 +  pgl_point *point, pgl_cluster *cluster, int samples
  5.1085 +) {
  5.1086 +  double distance;       /* shortest distance from point to cluster */
  5.1087 +  pgl_point *points;     /* sample points for Monte Carlo method */
  5.1088 +  double area;           /* area covered by sample points */
  5.1089 +  int i;
  5.1090 +  int matches = 0;       /* number of points within enlarged area (multiplied
  5.1091 +                            by two to be able to store half matches) */
  5.1092 +  double result;         /* result determined by Monte Carlo method */
  5.1093 +  /* limit sample count to avoid integer overflows on memory allocation */
  5.1094 +  if (samples > PGL_CLUSTER_MAXPOINTS) {
  5.1095 +    ereport(ERROR, (
  5.1096 +      errcode(ERRCODE_DATA_EXCEPTION),
  5.1097 +      errmsg(
  5.1098 +        "too many sample points for Monte Carlo method (maximum %i)",
  5.1099 +        PGL_CLUSTER_MAXPOINTS
  5.1100 +      )
  5.1101 +    ));
  5.1102 +  }
  5.1103 +  /* calculate shortest distance from point to cluster */
  5.1104 +  distance = pgl_point_cluster_distance(point, cluster);
  5.1105 +  /* if cluster consists of a single point or has no bounding circle with
  5.1106 +      positive radius, simply return distance */
  5.1107 +  if (
  5.1108 +    (cluster->nentries==1 && cluster->entries[0].entrytype==PGL_ENTRY_POINT) ||
  5.1109 +    !(cluster->bounding.radius > 0)
  5.1110 +  ) return distance;
  5.1111 +  /* if cluster consists of two points which are twice as far apart, return
  5.1112 +     distance between point and cluster multiplied by square root of two */
  5.1113 +  if (
  5.1114 +    cluster->nentries == 2 &&
  5.1115 +    cluster->entries[0].entrytype == PGL_ENTRY_POINT &&
  5.1116 +    cluster->entries[1].entrytype == PGL_ENTRY_POINT &&
  5.1117 +    pgl_distance(
  5.1118 +      PGL_ENTRY_POINTS(cluster, 0)[0].lat,
  5.1119 +      PGL_ENTRY_POINTS(cluster, 0)[0].lon,
  5.1120 +      PGL_ENTRY_POINTS(cluster, 1)[0].lat,
  5.1121 +      PGL_ENTRY_POINTS(cluster, 1)[0].lon
  5.1122 +    ) >= 2.0 * distance
  5.1123 +  ) {
  5.1124 +    return distance * M_SQRT2;
  5.1125 +  }
  5.1126 +  /* otherwise create sample points for Monte Carlo method and determine area
  5.1127 +     covered by sample points */
  5.1128 +  points = palloc(samples * sizeof(pgl_point));
  5.1129 +  area = pgl_sample_points(
  5.1130 +    &cluster->bounding.center,
  5.1131 +    cluster->bounding.radius + distance,  /* pad bounding circle by distance */
  5.1132 +    samples,
  5.1133 +    points
  5.1134 +  );
  5.1135 +  /* perform Monte Carlo method */
  5.1136 +  if (distance > 0) {
  5.1137 +    /* point is outside cluster */
  5.1138 +    for (i=0; i<samples; i++) {
  5.1139 +      /* count sample poitns within cluster as half match */
  5.1140 +      if (pgl_point_in_cluster(points+i, cluster, true)) matches += 1;
  5.1141 +      /* count sample points outside of cluster but within cluster grown by
  5.1142 +         distance as full match */
  5.1143 +      else if (
  5.1144 +        pgl_point_cluster_distance(points+i, cluster) < distance
  5.1145 +      ) matches += 2;
  5.1146 +    }
  5.1147 +  } else {
  5.1148 +    /* if point is within cluster,
  5.1149 +       just count sample points within cluster as half match */
  5.1150 +    for (i=0; i<samples; i++) {
  5.1151 +      if (pgl_point_in_cluster(points+i, cluster, true)) matches += 1;
  5.1152 +    }
  5.1153 +  }
  5.1154 +  /* release memory for sample points needed by Monte Carlo method */
  5.1155 +  pfree(points);
  5.1156 +  /* convert area determined by Monte Carlo method into distance
  5.1157 +     (i.e. radius of circle with the same area) */
  5.1158 +  result = sqrt(area * ((double)matches / 2.0 / samples) / M_PI);
  5.1159 +  /* return result only if it is greater than the distance between point and
  5.1160 +     cluster to avoid unexpected results because of errors due to limited
  5.1161 +     precision */
  5.1162 +  if (result > distance) return result;
  5.1163 +  /* otherwise return distance between point and cluster */
  5.1164 +  else return distance;
  5.1165 +}
  5.1166 +
  5.1167 +
  5.1168 +/*-------------------------------------------------*
  5.1169 + *  geographic index based on space-filling curve  *
  5.1170 + *-------------------------------------------------*/
  5.1171 +
  5.1172 +/* number of bytes used for geographic (center) position in keys */
  5.1173 +#define PGL_KEY_LATLON_BYTELEN 7
  5.1174 +
  5.1175 +/* maximum reference value for logarithmic size of geographic objects */
  5.1176 +#define PGL_AREAKEY_REFOBJSIZE (PGL_DIAMETER/3.0)  /* can be tweaked */
  5.1177 +
  5.1178 +/* pointer to index key (either pgl_pointkey or pgl_areakey) */
  5.1179 +typedef unsigned char *pgl_keyptr;
  5.1180 +
  5.1181 +/* index key for points (objects with zero area) on the spheroid */
  5.1182 +/* bit  0..55: interspersed bits of latitude and longitude,
  5.1183 +   bit 56..57: always zero,
  5.1184 +   bit 58..63: node depth in hypothetic (full) tree from 0 to 56 (incl.) */
  5.1185 +typedef unsigned char pgl_pointkey[PGL_KEY_LATLON_BYTELEN+1];
  5.1186 +
  5.1187 +/* index key for geographic objects on spheroid with area greater than zero */
  5.1188 +/* bit  0..55: interspersed bits of latitude and longitude of center point,
  5.1189 +   bit     56: always set to 1,
  5.1190 +   bit 57..63: node depth in hypothetic (full) tree from 0 to (2*56)+1 (incl.),
  5.1191 +   bit 64..71: logarithmic object size from 0 to 56+1 = 57 (incl.), but set to
  5.1192 +               PGL_KEY_OBJSIZE_EMPTY (with interspersed bits = 0 and node depth
  5.1193 +               = 113) for empty objects, and set to PGL_KEY_OBJSIZE_UNIVERSAL
  5.1194 +               (with interspersed bits = 0 and node depth = 0) for keys which
  5.1195 +               cover both empty and non-empty objects */
  5.1196 +
  5.1197 +typedef unsigned char pgl_areakey[PGL_KEY_LATLON_BYTELEN+2];
  5.1198 +
  5.1199 +/* helper macros for reading/writing index keys */
  5.1200 +#define PGL_KEY_NODEDEPTH_OFFSET  PGL_KEY_LATLON_BYTELEN
  5.1201 +#define PGL_KEY_OBJSIZE_OFFSET    (PGL_KEY_NODEDEPTH_OFFSET+1)
  5.1202 +#define PGL_POINTKEY_MAXDEPTH     (PGL_KEY_LATLON_BYTELEN*8)
  5.1203 +#define PGL_AREAKEY_MAXDEPTH      (2*PGL_POINTKEY_MAXDEPTH+1)
  5.1204 +#define PGL_AREAKEY_MAXOBJSIZE    (PGL_POINTKEY_MAXDEPTH+1)
  5.1205 +#define PGL_AREAKEY_TYPEMASK      0x80
  5.1206 +#define PGL_KEY_LATLONBIT(key, n) ((key)[(n)/8] & (0x80 >> ((n)%8)))
  5.1207 +#define PGL_KEY_LATLONBIT_DIFF(key1, key2, n) \
  5.1208 +                                  ( PGL_KEY_LATLONBIT(key1, n) ^ \
  5.1209 +                                    PGL_KEY_LATLONBIT(key2, n) )
  5.1210 +#define PGL_KEY_IS_AREAKEY(key)   ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \
  5.1211 +                                    PGL_AREAKEY_TYPEMASK)
  5.1212 +#define PGL_KEY_NODEDEPTH(key)    ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \
  5.1213 +                                    (PGL_AREAKEY_TYPEMASK-1))
  5.1214 +#define PGL_KEY_OBJSIZE(key)      ((key)[PGL_KEY_OBJSIZE_OFFSET])
  5.1215 +#define PGL_KEY_OBJSIZE_EMPTY     126
  5.1216 +#define PGL_KEY_OBJSIZE_UNIVERSAL 127
  5.1217 +#define PGL_KEY_IS_EMPTY(key)     ( PGL_KEY_IS_AREAKEY(key) && \
  5.1218 +                                    (key)[PGL_KEY_OBJSIZE_OFFSET] == \
  5.1219 +                                    PGL_KEY_OBJSIZE_EMPTY )
  5.1220 +#define PGL_KEY_IS_UNIVERSAL(key) ( PGL_KEY_IS_AREAKEY(key) && \
  5.1221 +                                    (key)[PGL_KEY_OBJSIZE_OFFSET] == \
  5.1222 +                                    PGL_KEY_OBJSIZE_UNIVERSAL )
  5.1223 +
  5.1224 +/* set area key to match empty objects only */
  5.1225 +static void pgl_key_set_empty(pgl_keyptr key) {
  5.1226 +  memset(key, 0, sizeof(pgl_areakey));
  5.1227 +  /* Note: setting node depth to maximum is required for picksplit function */
  5.1228 +  key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH;
  5.1229 +  key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_EMPTY;
  5.1230 +}
  5.1231 +
  5.1232 +/* set area key to match any object (including empty objects) */
  5.1233 +static void pgl_key_set_universal(pgl_keyptr key) {
  5.1234 +  memset(key, 0, sizeof(pgl_areakey));
  5.1235 +  key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK;
  5.1236 +  key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_UNIVERSAL;
  5.1237 +}
  5.1238 +
  5.1239 +/* convert a point on earth into a max-depth key to be used in index */
  5.1240 +static void pgl_point_to_key(pgl_point *point, pgl_keyptr key) {
  5.1241 +  double lat = point->lat;
  5.1242 +  double lon = point->lon;
  5.1243 +  int i;
  5.1244 +  /* clear latitude and longitude bits */
  5.1245 +  memset(key, 0, PGL_KEY_LATLON_BYTELEN);
  5.1246 +  /* set node depth to maximum and type bit to zero */
  5.1247 +  key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_POINTKEY_MAXDEPTH;
  5.1248 +  /* iterate over all latitude/longitude bit pairs */
  5.1249 +  for (i=0; i<PGL_POINTKEY_MAXDEPTH/2; i++) {
  5.1250 +    /* determine latitude bit */
  5.1251 +    if (lat >= 0) {
  5.1252 +      key[i/4] |= 0x80 >> (2*(i%4));
  5.1253 +      lat *= 2; lat -= 90;
  5.1254 +    } else {
  5.1255 +      lat *= 2; lat += 90;
  5.1256 +    }
  5.1257 +    /* determine longitude bit */
  5.1258 +    if (lon >= 0) {
  5.1259 +      key[i/4] |= 0x80 >> (2*(i%4)+1);
  5.1260 +      lon *= 2; lon -= 180;
  5.1261 +    } else {
  5.1262 +      lon *= 2; lon += 180;
  5.1263 +    }
  5.1264 +  }
  5.1265 +}
  5.1266 +
  5.1267 +/* convert a circle on earth into a max-depth key to be used in an index */
  5.1268 +static void pgl_circle_to_key(pgl_circle *circle, pgl_keyptr key) {
  5.1269 +  /* handle special case of empty circle */
  5.1270 +  if (circle->radius < 0) {
  5.1271 +    pgl_key_set_empty(key);
  5.1272 +    return;
  5.1273 +  }
  5.1274 +  /* perform same action as for point keys */
  5.1275 +  pgl_point_to_key(&(circle->center), key);
  5.1276 +  /* but overwrite type and node depth to fit area index key */
  5.1277 +  key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH;
  5.1278 +  /* check if radius is greater than (or equal to) reference size */
  5.1279 +  /* (treat equal values as greater values for numerical safety) */
  5.1280 +  if (circle->radius >= PGL_AREAKEY_REFOBJSIZE) {
  5.1281 +    /* if yes, set logarithmic size to zero */
  5.1282 +    key[PGL_KEY_OBJSIZE_OFFSET] = 0;
  5.1283 +  } else {
  5.1284 +    /* otherwise, determine logarithmic size iteratively */
  5.1285 +    /* (one step is equivalent to a factor of sqrt(2)) */
  5.1286 +    double reference = PGL_AREAKEY_REFOBJSIZE / M_SQRT2;
  5.1287 +    int objsize = 1;
  5.1288 +    while (objsize < PGL_AREAKEY_MAXOBJSIZE) {
  5.1289 +      /* stop when radius is greater than (or equal to) adjusted reference */
  5.1290 +      /* (treat equal values as greater values for numerical safety) */
  5.1291 +      if (circle->radius >= reference) break;
  5.1292 +      reference /= M_SQRT2;
  5.1293 +      objsize++;
  5.1294 +    }
  5.1295 +    /* set logarithmic size to determined value */
  5.1296 +    key[PGL_KEY_OBJSIZE_OFFSET] = objsize;
  5.1297 +  }
  5.1298 +}
  5.1299 +
  5.1300 +/* check if one key is subkey of another key or vice versa */
  5.1301 +static bool pgl_keys_overlap(pgl_keyptr key1, pgl_keyptr key2) {
  5.1302 +  int i;  /* key bit offset (includes both lat/lon and log. obj. size bits) */
  5.1303 +  /* determine smallest depth */
  5.1304 +  int depth1 = PGL_KEY_NODEDEPTH(key1);
  5.1305 +  int depth2 = PGL_KEY_NODEDEPTH(key2);
  5.1306 +  int depth = (depth1 < depth2) ? depth1 : depth2;
  5.1307 +  /* check if keys are area keys (assuming that both keys have same type) */
  5.1308 +  if (PGL_KEY_IS_AREAKEY(key1)) {
  5.1309 +    int j = 0;  /* bit offset for logarithmic object size bits */
  5.1310 +    int k = 0;  /* bit offset for latitude and longitude */
  5.1311 +    /* fetch logarithmic object size information */
  5.1312 +    int objsize1 = PGL_KEY_OBJSIZE(key1);
  5.1313 +    int objsize2 = PGL_KEY_OBJSIZE(key2);
  5.1314 +    /* handle special cases for empty objects (universal and empty keys) */
  5.1315 +    if (
  5.1316 +      objsize1 == PGL_KEY_OBJSIZE_UNIVERSAL ||
  5.1317 +      objsize2 == PGL_KEY_OBJSIZE_UNIVERSAL
  5.1318 +    ) return true;
  5.1319 +    if (
  5.1320 +      objsize1 == PGL_KEY_OBJSIZE_EMPTY ||
  5.1321 +      objsize2 == PGL_KEY_OBJSIZE_EMPTY
  5.1322 +    ) return objsize1 == objsize2;
  5.1323 +    /* iterate through key bits */
  5.1324 +    for (i=0; i<depth; i++) {
  5.1325 +      /* every second bit is a bit describing the object size */
  5.1326 +      if (i%2 == 0) {
  5.1327 +        /* check if object size bit is different in both keys (objsize1 and
  5.1328 +           objsize2 describe the minimum index when object size bit is set) */
  5.1329 +        if (
  5.1330 +          (objsize1 <= j && objsize2 > j) ||
  5.1331 +          (objsize2 <= j && objsize1 > j)
  5.1332 +        ) {
  5.1333 +          /* bit differs, therefore keys are in separate branches */
  5.1334 +          return false;
  5.1335 +        }
  5.1336 +        /* increase bit counter for object size bits */
  5.1337 +        j++;
  5.1338 +      }
  5.1339 +      /* all other bits describe latitude and longitude */
  5.1340 +      else {
  5.1341 +        /* check if bit differs in both keys */
  5.1342 +        if (PGL_KEY_LATLONBIT_DIFF(key1, key2, k)) {
  5.1343 +          /* bit differs, therefore keys are in separate branches */
  5.1344 +          return false;
  5.1345 +        }
  5.1346 +        /* increase bit counter for latitude/longitude bits */
  5.1347 +        k++;
  5.1348 +      }
  5.1349 +    }
  5.1350 +  }
  5.1351 +  /* if not, keys are point keys */
  5.1352 +  else {
  5.1353 +    /* iterate through key bits */
  5.1354 +    for (i=0; i<depth; i++) {
  5.1355 +      /* check if bit differs in both keys */
  5.1356 +      if (PGL_KEY_LATLONBIT_DIFF(key1, key2, i)) {
  5.1357 +        /* bit differs, therefore keys are in separate branches */
  5.1358 +        return false;
  5.1359 +      }
  5.1360 +    }
  5.1361 +  }
  5.1362 +  /* return true because keys are in the same branch */
  5.1363 +  return true;
  5.1364 +}
  5.1365 +
  5.1366 +/* combine two keys into new key which covers both original keys */
  5.1367 +/* (result stored in first argument) */
  5.1368 +static void pgl_unite_keys(pgl_keyptr dst, pgl_keyptr src) {
  5.1369 +  int i;  /* key bit offset (includes both lat/lon and log. obj. size bits) */
  5.1370 +  /* determine smallest depth */
  5.1371 +  int depth1 = PGL_KEY_NODEDEPTH(dst);
  5.1372 +  int depth2 = PGL_KEY_NODEDEPTH(src);
  5.1373 +  int depth = (depth1 < depth2) ? depth1 : depth2;
  5.1374 +  /* check if keys are area keys (assuming that both keys have same type) */
  5.1375 +  if (PGL_KEY_IS_AREAKEY(dst)) {
  5.1376 +    pgl_areakey dstbuf = { 0, };  /* destination buffer (cleared) */
  5.1377 +    int j = 0;  /* bit offset for logarithmic object size bits */
  5.1378 +    int k = 0;  /* bit offset for latitude and longitude */
  5.1379 +    /* fetch logarithmic object size information */
  5.1380 +    int objsize1 = PGL_KEY_OBJSIZE(dst);
  5.1381 +    int objsize2 = PGL_KEY_OBJSIZE(src);
  5.1382 +    /* handle special cases for empty objects (universal and empty keys) */
  5.1383 +    if (
  5.1384 +      objsize1 > PGL_AREAKEY_MAXOBJSIZE ||
  5.1385 +      objsize2 > PGL_AREAKEY_MAXOBJSIZE
  5.1386 +    ) {
  5.1387 +      if (
  5.1388 +        objsize1 == PGL_KEY_OBJSIZE_EMPTY &&
  5.1389 +        objsize2 == PGL_KEY_OBJSIZE_EMPTY
  5.1390 +      ) pgl_key_set_empty(dst);
  5.1391 +      else pgl_key_set_universal(dst);
  5.1392 +      return;
  5.1393 +    }
  5.1394 +    /* iterate through key bits */
  5.1395 +    for (i=0; i<depth; i++) {
  5.1396 +      /* every second bit is a bit describing the object size */
  5.1397 +      if (i%2 == 0) {
  5.1398 +        /* increase bit counter for object size bits first */
  5.1399 +        /* (handy when setting objsize variable) */
  5.1400 +        j++;
  5.1401 +        /* check if object size bit is set in neither key */
  5.1402 +        if (objsize1 >= j && objsize2 >= j) {
  5.1403 +          /* set objsize in destination buffer to indicate that size bit is
  5.1404 +             unset in destination buffer at the current bit position */
  5.1405 +          dstbuf[PGL_KEY_OBJSIZE_OFFSET] = j;
  5.1406 +        }
  5.1407 +        /* break if object size bit is set in one key only */
  5.1408 +        else if (objsize1 >= j || objsize2 >= j) break;
  5.1409 +      }
  5.1410 +      /* all other bits describe latitude and longitude */
  5.1411 +      else {
  5.1412 +        /* break if bit differs in both keys */
  5.1413 +        if (PGL_KEY_LATLONBIT(dst, k)) {
  5.1414 +          if (!PGL_KEY_LATLONBIT(src, k)) break;
  5.1415 +          /* but set bit in destination buffer if bit is set in both keys */
  5.1416 +          dstbuf[k/8] |= 0x80 >> (k%8);
  5.1417 +        } else if (PGL_KEY_LATLONBIT(src, k)) break;
  5.1418 +        /* increase bit counter for latitude/longitude bits */
  5.1419 +        k++;
  5.1420 +      }
  5.1421 +    }
  5.1422 +    /* set common node depth and type bit (type bit = 1) */
  5.1423 +    dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | i;
  5.1424 +    /* copy contents of destination buffer to first key */
  5.1425 +    memcpy(dst, dstbuf, sizeof(pgl_areakey));
  5.1426 +  }
  5.1427 +  /* if not, keys are point keys */
  5.1428 +  else {
  5.1429 +    pgl_pointkey dstbuf = { 0, };  /* destination buffer (cleared) */
  5.1430 +    /* iterate through key bits */
  5.1431 +    for (i=0; i<depth; i++) {
  5.1432 +      /* break if bit differs in both keys */
  5.1433 +      if (PGL_KEY_LATLONBIT(dst, i)) {
  5.1434 +        if (!PGL_KEY_LATLONBIT(src, i)) break;
  5.1435 +        /* but set bit in destination buffer if bit is set in both keys */
  5.1436 +        dstbuf[i/8] |= 0x80 >> (i%8);
  5.1437 +      } else if (PGL_KEY_LATLONBIT(src, i)) break;
  5.1438 +    }
  5.1439 +    /* set common node depth (type bit = 0) */
  5.1440 +    dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = i;
  5.1441 +    /* copy contents of destination buffer to first key */
  5.1442 +    memcpy(dst, dstbuf, sizeof(pgl_pointkey));
  5.1443 +  }
  5.1444 +}
  5.1445 +
  5.1446 +/* determine center(!) boundaries and radius estimation of index key */
  5.1447 +static double pgl_key_to_box(pgl_keyptr key, pgl_box *box) {
  5.1448 +  int i;
  5.1449 +  /* determine node depth */
  5.1450 +  int depth = PGL_KEY_NODEDEPTH(key);
  5.1451 +  /* center point of possible result */
  5.1452 +  double lat = 0;
  5.1453 +  double lon = 0;
  5.1454 +  /* maximum distance of real center point from key center */
  5.1455 +  double dlat = 90;
  5.1456 +  double dlon = 180;
  5.1457 +  /* maximum radius of contained objects */
  5.1458 +  double radius = 0;  /* always return zero for point index keys */
  5.1459 +  /* check if key is area key */
  5.1460 +  if (PGL_KEY_IS_AREAKEY(key)) {
  5.1461 +    /* get logarithmic object size */
  5.1462 +    int objsize = PGL_KEY_OBJSIZE(key);
  5.1463 +    /* handle special cases for empty objects (universal and empty keys) */
  5.1464 +    if (objsize == PGL_KEY_OBJSIZE_EMPTY) {
  5.1465 +      pgl_box_set_empty(box);
  5.1466 +      return 0;
  5.1467 +    } else if (objsize == PGL_KEY_OBJSIZE_UNIVERSAL) {
  5.1468 +      box->lat_min = -90;
  5.1469 +      box->lat_max =  90;
  5.1470 +      box->lon_min = -180;
  5.1471 +      box->lon_max =  180;
  5.1472 +      return 0;  /* any value >= 0 would do */
  5.1473 +    }
  5.1474 +    /* calculate maximum possible radius of objects covered by the given key */
  5.1475 +    if (objsize == 0) radius = INFINITY;
  5.1476 +    else {
  5.1477 +      radius = PGL_AREAKEY_REFOBJSIZE;
  5.1478 +      while (--objsize) radius /= M_SQRT2;
  5.1479 +    }
  5.1480 +    /* iterate over latitude and longitude bits in key */
  5.1481 +    /* (every second bit is a latitude or longitude bit) */
  5.1482 +    for (i=0; i<depth/2; i++) {
  5.1483 +      /* check if latitude bit */
  5.1484 +      if (i%2 == 0) {
  5.1485 +        /* cut latitude dimension in half */
  5.1486 +        dlat /= 2;
  5.1487 +        /* increase center latitude if bit is 1, otherwise decrease */
  5.1488 +        if (PGL_KEY_LATLONBIT(key, i)) lat += dlat;
  5.1489 +        else lat -= dlat;
  5.1490 +      }
  5.1491 +      /* otherwise longitude bit */
  5.1492 +      else {
  5.1493 +        /* cut longitude dimension in half */
  5.1494 +        dlon /= 2;
  5.1495 +        /* increase center longitude if bit is 1, otherwise decrease */
  5.1496 +        if (PGL_KEY_LATLONBIT(key, i)) lon += dlon;
  5.1497 +        else lon -= dlon;
  5.1498 +      }
  5.1499 +    }
  5.1500 +  }
  5.1501 +  /* if not, keys are point keys */
  5.1502 +  else {
  5.1503 +    /* iterate over all bits in key */
  5.1504 +    for (i=0; i<depth; i++) {
  5.1505 +      /* check if latitude bit */
  5.1506 +      if (i%2 == 0) {
  5.1507 +        /* cut latitude dimension in half */
  5.1508 +        dlat /= 2;
  5.1509 +        /* increase center latitude if bit is 1, otherwise decrease */
  5.1510 +        if (PGL_KEY_LATLONBIT(key, i)) lat += dlat;
  5.1511 +        else lat -= dlat;
  5.1512 +      }
  5.1513 +      /* otherwise longitude bit */
  5.1514 +      else {
  5.1515 +        /* cut longitude dimension in half */
  5.1516 +        dlon /= 2;
  5.1517 +        /* increase center longitude if bit is 1, otherwise decrease */
  5.1518 +        if (PGL_KEY_LATLONBIT(key, i)) lon += dlon;
  5.1519 +        else lon -= dlon;
  5.1520 +      }
  5.1521 +    }
  5.1522 +  }
  5.1523 +  /* calculate boundaries from center point and remaining dlat and dlon */
  5.1524 +  /* (return values through pointer to box) */
  5.1525 +  box->lat_min = lat - dlat;
  5.1526 +  box->lat_max = lat + dlat;
  5.1527 +  box->lon_min = lon - dlon;
  5.1528 +  box->lon_max = lon + dlon;
  5.1529 +  /* return radius (as a function return value) */
  5.1530 +  return radius;
  5.1531 +}
  5.1532 +
  5.1533 +/* estimator function for distance between point and index key */
  5.1534 +/* always returns a smaller value than actually correct or zero */
  5.1535 +static double pgl_estimate_key_distance(pgl_keyptr key, pgl_point *point) {
  5.1536 +  pgl_box box;  /* center(!) bounding box of area index key */
  5.1537 +  /* calculate center(!) bounding box and maximum radius of objects covered
  5.1538 +     by area index key (radius is zero for point index keys) */
  5.1539 +  double distance = pgl_key_to_box(key, &box);
  5.1540 +  /* calculate estimated distance between bounding box of center point of
  5.1541 +     indexed object and point passed as second argument, then substract maximum
  5.1542 +     radius of objects covered by index key */
  5.1543 +  distance = pgl_estimate_point_box_distance(point, &box) - distance;
  5.1544 +  /* truncate negative results to zero */
  5.1545 +  if (distance <= 0) distance = 0;
  5.1546 +  /* return result */
  5.1547 +  return distance;
  5.1548 +}
  5.1549 +
  5.1550 +
  5.1551 +/*---------------------------------*
  5.1552 + *  helper functions for text I/O  *
  5.1553 + *---------------------------------*/
  5.1554 +
  5.1555 +#define PGL_NUMBUFLEN 64  /* buffer size for number to string conversion */
  5.1556 +
  5.1557 +/* convert floating point number to string (round-trip safe) */
  5.1558 +static void pgl_print_float(char *buf, double flt) {
  5.1559 +  /* check if number is integral */
  5.1560 +  if (trunc(flt) == flt) {
  5.1561 +    /* for integral floats use maximum precision */
  5.1562 +    snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt);
  5.1563 +  } else {
  5.1564 +    /* otherwise check if 15, 16, or 17 digits needed (round-trip safety) */
  5.1565 +    snprintf(buf, PGL_NUMBUFLEN, "%.15g", flt);
  5.1566 +    if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.16g", flt);
  5.1567 +    if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt);
  5.1568 +  }
  5.1569 +}
  5.1570 +
  5.1571 +/* convert latitude floating point number (in degrees) to string */
  5.1572 +static void pgl_print_lat(char *buf, double lat) {
  5.1573 +  if (signbit(lat)) {
  5.1574 +    /* treat negative latitudes (including -0) as south */
  5.1575 +    snprintf(buf, PGL_NUMBUFLEN, "S%015.12f", -lat);
  5.1576 +  } else {
  5.1577 +    /* treat positive latitudes (including +0) as north */
  5.1578 +    snprintf(buf, PGL_NUMBUFLEN, "N%015.12f", lat);
  5.1579 +  }
  5.1580 +}
  5.1581 +
  5.1582 +/* convert longitude floating point number (in degrees) to string */
  5.1583 +static void pgl_print_lon(char *buf, double lon) {
  5.1584 +  if (signbit(lon)) {
  5.1585 +    /* treat negative longitudes (including -0) as west */
  5.1586 +    snprintf(buf, PGL_NUMBUFLEN, "W%016.12f", -lon);
  5.1587 +  } else {
  5.1588 +    /* treat positive longitudes (including +0) as east */
  5.1589 +    snprintf(buf, PGL_NUMBUFLEN, "E%016.12f", lon);
  5.1590 +  }
  5.1591 +}
  5.1592 +
  5.1593 +/* bit masks used as return value of pgl_scan() function */
  5.1594 +#define PGL_SCAN_NONE 0      /* no value has been parsed */
  5.1595 +#define PGL_SCAN_LAT (1<<0)  /* latitude has been parsed */
  5.1596 +#define PGL_SCAN_LON (1<<1)  /* longitude has been parsed */
  5.1597 +#define PGL_SCAN_LATLON (PGL_SCAN_LAT | PGL_SCAN_LON)  /* bitwise OR of both */
  5.1598 +
  5.1599 +/* parse a coordinate (can be latitude or longitude) */
  5.1600 +static int pgl_scan(char **str, double *lat, double *lon) {
  5.1601 +  double val;
  5.1602 +  int len;
  5.1603 +  if (
  5.1604 +    sscanf(*str, " N %lf %n", &val, &len) ||
  5.1605 +    sscanf(*str, " n %lf %n", &val, &len)
  5.1606 +  ) {
  5.1607 +    *str += len; *lat = val; return PGL_SCAN_LAT;
  5.1608 +  }
  5.1609 +  if (
  5.1610 +    sscanf(*str, " S %lf %n", &val, &len) ||
  5.1611 +    sscanf(*str, " s %lf %n", &val, &len)
  5.1612 +  ) {
  5.1613 +    *str += len; *lat = -val; return PGL_SCAN_LAT;
  5.1614 +  }
  5.1615 +  if (
  5.1616 +    sscanf(*str, " E %lf %n", &val, &len) ||
  5.1617 +    sscanf(*str, " e %lf %n", &val, &len)
  5.1618 +  ) {
  5.1619 +    *str += len; *lon = val; return PGL_SCAN_LON;
  5.1620 +  }
  5.1621 +  if (
  5.1622 +    sscanf(*str, " W %lf %n", &val, &len) ||
  5.1623 +    sscanf(*str, " w %lf %n", &val, &len)
  5.1624 +  ) {
  5.1625 +    *str += len; *lon = -val; return PGL_SCAN_LON;
  5.1626 +  }
  5.1627 +  return PGL_SCAN_NONE;
  5.1628 +}
  5.1629 +
  5.1630 +
  5.1631 +/*-----------------*
  5.1632 + *  SQL functions  *
  5.1633 + *-----------------*/
  5.1634 +
  5.1635 +/* Note: These function names use "epoint", "ebox", etc. notation here instead
  5.1636 +   of "point", "box", etc. in order to distinguish them from any previously
  5.1637 +   defined functions. */
  5.1638 +
  5.1639 +/* function needed for dummy types and/or not implemented features */
  5.1640 +PG_FUNCTION_INFO_V1(pgl_notimpl);
  5.1641 +Datum pgl_notimpl(PG_FUNCTION_ARGS) {
  5.1642 +  ereport(ERROR, (errmsg("not implemented by pgLatLon")));
  5.1643 +}
  5.1644 +
  5.1645 +/* set point to latitude and longitude (including checks) */
  5.1646 +static void pgl_epoint_set_latlon(pgl_point *point, double lat, double lon) {
  5.1647 +  /* reject infinite or NaN values */
  5.1648 +  if (!isfinite(lat) || !isfinite(lon)) {
  5.1649 +    ereport(ERROR, (
  5.1650 +      errcode(ERRCODE_DATA_EXCEPTION),
  5.1651 +      errmsg("epoint requires finite coordinates")
  5.1652 +    ));
  5.1653 +  }
  5.1654 +  /* check latitude bounds */
  5.1655 +  if (lat < -90) {
  5.1656 +    ereport(WARNING, (errmsg("latitude exceeds south pole")));
  5.1657 +    lat = -90;
  5.1658 +  } else if (lat > 90) {
  5.1659 +    ereport(WARNING, (errmsg("latitude exceeds north pole")));
  5.1660 +    lat = 90;
  5.1661 +  }
  5.1662 +  /* check longitude bounds */
  5.1663 +  if (lon < -180) {
  5.1664 +    ereport(NOTICE, (errmsg("longitude west of 180th meridian normalized")));
  5.1665 +    lon += 360 - trunc(lon / 360) * 360;
  5.1666 +  } else if (lon > 180) {
  5.1667 +    ereport(NOTICE, (errmsg("longitude east of 180th meridian normalized")));
  5.1668 +    lon -= 360 + trunc(lon / 360) * 360;
  5.1669 +  }
  5.1670 +  /* store rounded latitude/longitude values for round-trip safety */
  5.1671 +  point->lat = pgl_round(lat);
  5.1672 +  point->lon = pgl_round(lon);
  5.1673 +}
  5.1674 +
  5.1675 +/* create point ("epoint" in SQL) from latitude and longitude */
  5.1676 +PG_FUNCTION_INFO_V1(pgl_create_epoint);
  5.1677 +Datum pgl_create_epoint(PG_FUNCTION_ARGS) {
  5.1678 +  pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point));
  5.1679 +  pgl_epoint_set_latlon(point, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1));
  5.1680 +  PG_RETURN_POINTER(point);
  5.1681 +}
  5.1682 +
  5.1683 +/* parse point ("epoint" in SQL) */
  5.1684 +/* format: '[NS]<float> [EW]<float>' */
  5.1685 +PG_FUNCTION_INFO_V1(pgl_epoint_in);
  5.1686 +Datum pgl_epoint_in(PG_FUNCTION_ARGS) {
  5.1687 +  char *str = PG_GETARG_CSTRING(0);  /* input string */
  5.1688 +  char *strptr = str;  /* current position within string */
  5.1689 +  int done = 0;        /* bit mask storing if latitude or longitude was read */
  5.1690 +  double lat, lon;     /* parsed values as double precision floats */
  5.1691 +  pgl_point *point;    /* return value (to be palloc'ed) */
  5.1692 +  /* parse two floats (each latitude or longitude) separated by white-space */
  5.1693 +  done |= pgl_scan(&strptr, &lat, &lon);
  5.1694 +  if (strptr != str && isspace(strptr[-1])) {
  5.1695 +    done |= pgl_scan(&strptr, &lat, &lon);
  5.1696 +  }
  5.1697 +  /* require end of string, and latitude and longitude parsed successfully */
  5.1698 +  if (strptr[0] || done != PGL_SCAN_LATLON) {
  5.1699 +    ereport(ERROR, (
  5.1700 +      errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  5.1701 +      errmsg("invalid input syntax for type epoint: \"%s\"", str)
  5.1702 +    ));
  5.1703 +  }
  5.1704 +  /* allocate memory for result */
  5.1705 +  point = (pgl_point *)palloc(sizeof(pgl_point));
  5.1706 +  /* set latitude and longitude (and perform checks) */
  5.1707 +  pgl_epoint_set_latlon(point, lat, lon);
  5.1708 +  /* return result */
  5.1709 +  PG_RETURN_POINTER(point);
  5.1710 +}
  5.1711 +
  5.1712 +/* create box ("ebox" in SQL) that is empty */
  5.1713 +PG_FUNCTION_INFO_V1(pgl_create_empty_ebox);
  5.1714 +Datum pgl_create_empty_ebox(PG_FUNCTION_ARGS) {
  5.1715 +  pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
  5.1716 +  pgl_box_set_empty(box);
  5.1717 +  PG_RETURN_POINTER(box);
  5.1718 +}
  5.1719 +
  5.1720 +/* set box to given boundaries (including checks) */
  5.1721 +static void pgl_ebox_set_boundaries(
  5.1722 +  pgl_box *box,
  5.1723 +  double lat_min, double lat_max, double lon_min, double lon_max
  5.1724 +) {
  5.1725 +  /* if minimum latitude is greater than maximum latitude, return empty box */
  5.1726 +  if (lat_min > lat_max) {
  5.1727 +    pgl_box_set_empty(box);
  5.1728 +    return;
  5.1729 +  }
  5.1730 +  /* otherwise reject infinite or NaN values */
  5.1731 +  if (
  5.1732 +    !isfinite(lat_min) || !isfinite(lat_max) ||
  5.1733 +    !isfinite(lon_min) || !isfinite(lon_max)
  5.1734 +  ) {
  5.1735 +    ereport(ERROR, (
  5.1736 +      errcode(ERRCODE_DATA_EXCEPTION),
  5.1737 +      errmsg("ebox requires finite coordinates")
  5.1738 +    ));
  5.1739 +  }
  5.1740 +  /* check latitude bounds */
  5.1741 +  if (lat_max < -90) {
  5.1742 +    ereport(WARNING, (errmsg("northern latitude exceeds south pole")));
  5.1743 +    lat_max = -90;
  5.1744 +  } else if (lat_max > 90) {
  5.1745 +    ereport(WARNING, (errmsg("northern latitude exceeds north pole")));
  5.1746 +    lat_max = 90;
  5.1747 +  }
  5.1748 +  if (lat_min < -90) {
  5.1749 +    ereport(WARNING, (errmsg("southern latitude exceeds south pole")));
  5.1750 +    lat_min = -90;
  5.1751 +  } else if (lat_min > 90) {
  5.1752 +    ereport(WARNING, (errmsg("southern latitude exceeds north pole")));
  5.1753 +    lat_min = 90;
  5.1754 +  }
  5.1755 +  /* check if all longitudes are included */
  5.1756 +  if (lon_max - lon_min >= 360) {
  5.1757 +    if (lon_max - lon_min > 360) ereport(WARNING, (
  5.1758 +      errmsg("longitude coverage greater than 360 degrees")
  5.1759 +    ));
  5.1760 +    lon_min = -180;
  5.1761 +    lon_max = 180;
  5.1762 +  } else {
  5.1763 +    /* normalize longitude bounds */
  5.1764 +    if      (lon_min < -180) lon_min += 360 - trunc(lon_min / 360) * 360;
  5.1765 +    else if (lon_min >  180) lon_min -= 360 + trunc(lon_min / 360) * 360;
  5.1766 +    if      (lon_max < -180) lon_max += 360 - trunc(lon_max / 360) * 360;
  5.1767 +    else if (lon_max >  180) lon_max -= 360 + trunc(lon_max / 360) * 360;
  5.1768 +  }
  5.1769 +  /* store rounded latitude/longitude values for round-trip safety */
  5.1770 +  box->lat_min = pgl_round(lat_min);
  5.1771 +  box->lat_max = pgl_round(lat_max);
  5.1772 +  box->lon_min = pgl_round(lon_min);
  5.1773 +  box->lon_max = pgl_round(lon_max);
  5.1774 +  /* ensure that rounding does not change orientation */
  5.1775 +  if (lon_min > lon_max && box->lon_min == box->lon_max) {
  5.1776 +    box->lon_min = -180;
  5.1777 +    box->lon_max = 180;
  5.1778 +  }
  5.1779 +}
  5.1780 +
  5.1781 +/* create box ("ebox" in SQL) from min/max latitude and min/max longitude */
  5.1782 +PG_FUNCTION_INFO_V1(pgl_create_ebox);
  5.1783 +Datum pgl_create_ebox(PG_FUNCTION_ARGS) {
  5.1784 +  pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
  5.1785 +  pgl_ebox_set_boundaries(
  5.1786 +    box,
  5.1787 +    PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1),
  5.1788 +    PG_GETARG_FLOAT8(2), PG_GETARG_FLOAT8(3)
  5.1789 +  );
  5.1790 +  PG_RETURN_POINTER(box);
  5.1791 +}
  5.1792 +
  5.1793 +/* create box ("ebox" in SQL) from two points ("epoint"s) */
  5.1794 +/* (can not be used to cover a longitude range of more than 120 degrees) */
  5.1795 +PG_FUNCTION_INFO_V1(pgl_create_ebox_from_epoints);
  5.1796 +Datum pgl_create_ebox_from_epoints(PG_FUNCTION_ARGS) {
  5.1797 +  pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0);
  5.1798 +  pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1);
  5.1799 +  pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
  5.1800 +  double lat_min, lat_max, lon_min, lon_max;
  5.1801 +  double dlon;  /* longitude range (delta longitude) */
  5.1802 +  /* order latitude and longitude boundaries */
  5.1803 +  if (point2->lat < point1->lat) {
  5.1804 +    lat_min = point2->lat;
  5.1805 +    lat_max = point1->lat;
  5.1806 +  } else {
  5.1807 +    lat_min = point1->lat;
  5.1808 +    lat_max = point2->lat;
  5.1809 +  }
  5.1810 +  if (point2->lon < point1->lon) {
  5.1811 +    lon_min = point2->lon;
  5.1812 +    lon_max = point1->lon;
  5.1813 +  } else {
  5.1814 +    lon_min = point1->lon;
  5.1815 +    lon_max = point2->lon;
  5.1816 +  }
  5.1817 +  /* calculate longitude range (round to avoid floating point errors) */
  5.1818 +  dlon = pgl_round(lon_max - lon_min);
  5.1819 +  /* determine east-west direction */
  5.1820 +  if (dlon >= 240) {
  5.1821 +    /* assume that 180th meridian is crossed and swap min/max longitude */
  5.1822 +    double swap = lon_min; lon_min = lon_max; lon_max = swap;
  5.1823 +  } else if (dlon > 120) {
  5.1824 +    /* unclear orientation since delta longitude > 120 */
  5.1825 +    ereport(ERROR, (
  5.1826 +      errcode(ERRCODE_DATA_EXCEPTION),
  5.1827 +      errmsg("can not determine east/west orientation for ebox")
  5.1828 +    ));
  5.1829 +  }
  5.1830 +  /* use boundaries to setup box (and perform checks) */
  5.1831 +  pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max);
  5.1832 +  /* return result */
  5.1833 +  PG_RETURN_POINTER(box);
  5.1834 +}
  5.1835 +
  5.1836 +/* parse box ("ebox" in SQL) */
  5.1837 +/* format: '[NS]<float> [EW]<float> [NS]<float> [EW]<float>'
  5.1838 +       or: '[NS]<float> [NS]<float> [EW]<float> [EW]<float>' */
  5.1839 +PG_FUNCTION_INFO_V1(pgl_ebox_in);
  5.1840 +Datum pgl_ebox_in(PG_FUNCTION_ARGS) {
  5.1841 +  char *str = PG_GETARG_CSTRING(0);  /* input string */
  5.1842 +  char *str_lower;     /* lower case version of input string */
  5.1843 +  char *strptr;        /* current position within string */
  5.1844 +  int valid;           /* number of valid chars */
  5.1845 +  int done;            /* specifies if latitude or longitude was read */
  5.1846 +  double val;          /* temporary variable */
  5.1847 +  int lat_count = 0;   /* count of latitude values parsed */
  5.1848 +  int lon_count = 0;   /* count of longitufde values parsed */
  5.1849 +  double lat_min, lat_max, lon_min, lon_max;  /* see pgl_box struct */
  5.1850 +  pgl_box *box;        /* return value (to be palloc'ed) */
  5.1851 +  /* lowercase input */
  5.1852 +  str_lower = psprintf("%s", str);
  5.1853 +  for (strptr=str_lower; *strptr; strptr++) {
  5.1854 +    if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A';
  5.1855 +  }
  5.1856 +  /* reset reading position to start of (lowercase) string */
  5.1857 +  strptr = str_lower;
  5.1858 +  /* check if empty box */
  5.1859 +  valid = 0;
  5.1860 +  sscanf(strptr, " empty %n", &valid);
  5.1861 +  if (valid && strptr[valid] == 0) {
  5.1862 +    /* allocate and return empty box */
  5.1863 +    box = (pgl_box *)palloc(sizeof(pgl_box));
  5.1864 +    pgl_box_set_empty(box);
  5.1865 +    PG_RETURN_POINTER(box);
  5.1866 +  }
  5.1867 +  /* demand four blocks separated by whitespace */
  5.1868 +  valid = 0;
  5.1869 +  sscanf(strptr, " %*s %*s %*s %*s %n", &valid);
  5.1870 +  /* if four blocks separated by whitespace exist, parse those blocks */
  5.1871 +  if (strptr[valid] == 0) while (strptr[0]) {
  5.1872 +    /* parse either latitude or longitude (whichever found in input string) */
  5.1873 +    done = pgl_scan(&strptr, &val, &val);
  5.1874 +    /* store latitude or longitude in lat_min, lat_max, lon_min, or lon_max */
  5.1875 +    if (done == PGL_SCAN_LAT) {
  5.1876 +      if (!lat_count) lat_min = val; else lat_max = val;
  5.1877 +      lat_count++;
  5.1878 +    } else if (done == PGL_SCAN_LON) {
  5.1879 +      if (!lon_count) lon_min = val; else lon_max = val;
  5.1880 +      lon_count++;
  5.1881 +    } else {
  5.1882 +      break;
  5.1883 +    }
  5.1884 +  }
  5.1885 +  /* require end of string, and two latitude and two longitude values */
  5.1886 +  if (strptr[0] || lat_count != 2 || lon_count != 2) {
  5.1887 +    ereport(ERROR, (
  5.1888 +      errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  5.1889 +      errmsg("invalid input syntax for type ebox: \"%s\"", str)
  5.1890 +    ));
  5.1891 +  }
  5.1892 +  /* free lower case string */
  5.1893 +  pfree(str_lower);
  5.1894 +  /* order boundaries (maximum greater than minimum) */
  5.1895 +  if (lat_min > lat_max) { val = lat_min; lat_min = lat_max; lat_max = val; }
  5.1896 +  if (lon_min > lon_max) { val = lon_min; lon_min = lon_max; lon_max = val; }
  5.1897 +  /* allocate memory for result */
  5.1898 +  box = (pgl_box *)palloc(sizeof(pgl_box));
  5.1899 +  /* set boundaries (and perform checks) */
  5.1900 +  pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max);
  5.1901 +  /* return result */
  5.1902 +  PG_RETURN_POINTER(box);
  5.1903 +}
  5.1904 +
  5.1905 +/* set circle to given latitude, longitude, and radius (including checks) */
  5.1906 +static void pgl_ecircle_set_latlon_radius(
  5.1907 +  pgl_circle *circle, double lat, double lon, double radius
  5.1908 +) {
  5.1909 +  /* set center point (including checks) */
  5.1910 +  pgl_epoint_set_latlon(&(circle->center), lat, lon);
  5.1911 +  /* handle non-positive radius */
  5.1912 +  if (isnan(radius)) {
  5.1913 +    ereport(ERROR, (
  5.1914 +      errcode(ERRCODE_DATA_EXCEPTION),
  5.1915 +      errmsg("invalid radius for ecircle")
  5.1916 +    ));
  5.1917 +  }
  5.1918 +  if (radius == 0) radius = 0;  /* avoids -0 */
  5.1919 +  else if (radius < 0) {
  5.1920 +    if (isfinite(radius)) {
  5.1921 +      ereport(NOTICE, (errmsg("negative radius converted to minus infinity")));
  5.1922 +    }
  5.1923 +    radius = -INFINITY;
  5.1924 +  }
  5.1925 +  /* store radius (round-trip safety is ensured by pgl_print_float) */
  5.1926 +  circle->radius = radius;
  5.1927 +}
  5.1928 +
  5.1929 +/* create circle ("ecircle" in SQL) from latitude, longitude, and radius */
  5.1930 +PG_FUNCTION_INFO_V1(pgl_create_ecircle);
  5.1931 +Datum pgl_create_ecircle(PG_FUNCTION_ARGS) {
  5.1932 +  pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
  5.1933 +  pgl_ecircle_set_latlon_radius(
  5.1934 +    circle, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1), PG_GETARG_FLOAT8(2)
  5.1935 +  );
  5.1936 +  PG_RETURN_POINTER(circle);
  5.1937 +}
  5.1938 +
  5.1939 +/* create circle ("ecircle" in SQL) from point ("epoint"), and radius */
  5.1940 +PG_FUNCTION_INFO_V1(pgl_create_ecircle_from_epoint);
  5.1941 +Datum pgl_create_ecircle_from_epoint(PG_FUNCTION_ARGS) {
  5.1942 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  5.1943 +  double radius = PG_GETARG_FLOAT8(1);
  5.1944 +  pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
  5.1945 +  /* set latitude, longitude, radius (and perform checks) */
  5.1946 +  pgl_ecircle_set_latlon_radius(circle, point->lat, point->lon, radius);
  5.1947 +  /* return result */
  5.1948 +  PG_RETURN_POINTER(circle);
  5.1949 +}
  5.1950 +
  5.1951 +/* parse circle ("ecircle" in SQL) */
  5.1952 +/* format: '[NS]<float> [EW]<float> <float>' */
  5.1953 +PG_FUNCTION_INFO_V1(pgl_ecircle_in);
  5.1954 +Datum pgl_ecircle_in(PG_FUNCTION_ARGS) {
  5.1955 +  char *str = PG_GETARG_CSTRING(0);  /* input string */
  5.1956 +  char *strptr = str;       /* current position within string */
  5.1957 +  double lat, lon, radius;  /* parsed values as double precision flaots */
  5.1958 +  int valid = 0;            /* number of valid chars */
  5.1959 +  int done = 0;             /* stores if latitude and/or longitude was read */
  5.1960 +  pgl_circle *circle;       /* return value (to be palloc'ed) */
  5.1961 +  /* demand three blocks separated by whitespace */
  5.1962 +  sscanf(strptr, " %*s %*s %*s %n", &valid);
  5.1963 +  /* if three blocks separated by whitespace exist, parse those blocks */
  5.1964 +  if (strptr[valid] == 0) {
  5.1965 +    /* parse latitude and longitude */
  5.1966 +    done |= pgl_scan(&strptr, &lat, &lon);
  5.1967 +    done |= pgl_scan(&strptr, &lat, &lon);
  5.1968 +    /* parse radius (while incrementing strptr by number of bytes parsed) */
  5.1969 +    valid = 0;
  5.1970 +    if (sscanf(strptr, " %lf %n", &radius, &valid) == 1) strptr += valid;
  5.1971 +  }
  5.1972 +  /* require end of string and both latitude and longitude being parsed */
  5.1973 +  if (strptr[0] || done != PGL_SCAN_LATLON) {
  5.1974 +    ereport(ERROR, (
  5.1975 +      errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  5.1976 +      errmsg("invalid input syntax for type ecircle: \"%s\"", str)
  5.1977 +    ));
  5.1978 +  }
  5.1979 +  /* allocate memory for result */
  5.1980 +  circle = (pgl_circle *)palloc(sizeof(pgl_circle));
  5.1981 +  /* set latitude, longitude, radius (and perform checks) */
  5.1982 +  pgl_ecircle_set_latlon_radius(circle, lat, lon, radius);
  5.1983 +  /* return result */
  5.1984 +  PG_RETURN_POINTER(circle);
  5.1985 +}
  5.1986 +
  5.1987 +/* parse cluster ("ecluster" in SQL) */
  5.1988 +PG_FUNCTION_INFO_V1(pgl_ecluster_in);
  5.1989 +Datum pgl_ecluster_in(PG_FUNCTION_ARGS) {
  5.1990 +  int i;
  5.1991 +  char *str = PG_GETARG_CSTRING(0);  /* input string */
  5.1992 +  char *str_lower;         /* lower case version of input string */
  5.1993 +  char *strptr;            /* pointer to current reading position of input */
  5.1994 +  int npoints_total = 0;   /* total number of points in cluster */
  5.1995 +  int nentries = 0;        /* total number of entries */
  5.1996 +  pgl_newentry *entries;   /* array of pgl_newentry to create pgl_cluster */
  5.1997 +  int entries_buflen = 4;  /* maximum number of elements in entries array */
  5.1998 +  int valid;               /* number of valid chars processed */
  5.1999 +  double lat, lon;         /* latitude and longitude of parsed point */
  5.2000 +  int entrytype;           /* current entry type */
  5.2001 +  int npoints;             /* number of points in current entry */
  5.2002 +  pgl_point *points;       /* array of pgl_point for pgl_newentry */
  5.2003 +  int points_buflen;       /* maximum number of elements in points array */
  5.2004 +  int done;                /* return value of pgl_scan function */
  5.2005 +  pgl_cluster *cluster;    /* created cluster */
  5.2006 +  /* lowercase input */
  5.2007 +  str_lower = psprintf("%s", str);
  5.2008 +  for (strptr=str_lower; *strptr; strptr++) {
  5.2009 +    if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A';
  5.2010 +  }
  5.2011 +  /* reset reading position to start of (lowercase) string */
  5.2012 +  strptr = str_lower;
  5.2013 +  /* allocate initial buffer for entries */
  5.2014 +  entries = palloc(entries_buflen * sizeof(pgl_newentry));
  5.2015 +  /* parse until end of string */
  5.2016 +  while (strptr[0]) {
  5.2017 +    /* require previous white-space or closing parenthesis before next token */
  5.2018 +    if (strptr != str_lower && !isspace(strptr[-1]) && strptr[-1] != ')') {
  5.2019 +      goto pgl_ecluster_in_error;
  5.2020 +    }
  5.2021 +    /* ignore token "empty" */
  5.2022 +    valid = 0; sscanf(strptr, " empty %n", &valid);
  5.2023 +    if (valid) { strptr += valid; continue; }
  5.2024 +    /* test for "point" token */
  5.2025 +    valid = 0; sscanf(strptr, " point ( %n", &valid);
  5.2026 +    if (valid) {
  5.2027 +      strptr += valid;
  5.2028 +      entrytype = PGL_ENTRY_POINT;
  5.2029 +      goto pgl_ecluster_in_type_ok;
  5.2030 +    }
  5.2031 +    /* test for "path" token */
  5.2032 +    valid = 0; sscanf(strptr, " path ( %n", &valid);
  5.2033 +    if (valid) {
  5.2034 +      strptr += valid;
  5.2035 +      entrytype = PGL_ENTRY_PATH;
  5.2036 +      goto pgl_ecluster_in_type_ok;
  5.2037 +    }
  5.2038 +    /* test for "outline" token */
  5.2039 +    valid = 0; sscanf(strptr, " outline ( %n", &valid);
  5.2040 +    if (valid) {
  5.2041 +      strptr += valid;
  5.2042 +      entrytype = PGL_ENTRY_OUTLINE;
  5.2043 +      goto pgl_ecluster_in_type_ok;
  5.2044 +    }
  5.2045 +    /* test for "polygon" token */
  5.2046 +    valid = 0; sscanf(strptr, " polygon ( %n", &valid);
  5.2047 +    if (valid) {
  5.2048 +      strptr += valid;
  5.2049 +      entrytype = PGL_ENTRY_POLYGON;
  5.2050 +      goto pgl_ecluster_in_type_ok;
  5.2051 +    }
  5.2052 +    /* error if no valid token found */
  5.2053 +    goto pgl_ecluster_in_error;
  5.2054 +    pgl_ecluster_in_type_ok:
  5.2055 +    /* check if pgl_newentry array needs to grow */
  5.2056 +    if (nentries == entries_buflen) {
  5.2057 +      pgl_newentry *newbuf;
  5.2058 +      entries_buflen *= 2;
  5.2059 +      newbuf = palloc(entries_buflen * sizeof(pgl_newentry));
  5.2060 +      memcpy(newbuf, entries, nentries * sizeof(pgl_newentry));
  5.2061 +      pfree(entries);
  5.2062 +      entries = newbuf;
  5.2063 +    }
  5.2064 +    /* reset number of points for current entry */
  5.2065 +    npoints = 0;
  5.2066 +    /* allocate array for points */
  5.2067 +    points_buflen = 4;
  5.2068 +    points = palloc(points_buflen * sizeof(pgl_point));
  5.2069 +    /* parse until closing parenthesis */
  5.2070 +    while (strptr[0] != ')') {
  5.2071 +      /* error on unexpected end of string */
  5.2072 +      if (strptr[0] == 0) goto pgl_ecluster_in_error;
  5.2073 +      /* mark neither latitude nor longitude as read */
  5.2074 +      done = PGL_SCAN_NONE;
  5.2075 +      /* require white-space before second, third, etc. point */
  5.2076 +      if (npoints != 0 && !isspace(strptr[-1])) goto pgl_ecluster_in_error;
  5.2077 +      /* scan latitude (or longitude) */
  5.2078 +      done |= pgl_scan(&strptr, &lat, &lon);
  5.2079 +      /* require white-space before second coordinate */
  5.2080 +      if (strptr != str && !isspace(strptr[-1])) goto pgl_ecluster_in_error;
  5.2081 +      /* scan longitude (or latitude) */
  5.2082 +      done |= pgl_scan(&strptr, &lat, &lon);
  5.2083 +      /* error unless both latitude and longitude were parsed */
  5.2084 +      if (done != PGL_SCAN_LATLON) goto pgl_ecluster_in_error;
  5.2085 +      /* throw error if number of points is too high */
  5.2086 +      if (npoints_total == PGL_CLUSTER_MAXPOINTS) {
  5.2087 +        ereport(ERROR, (
  5.2088 +          errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  5.2089 +          errmsg(
  5.2090 +            "too many points for ecluster entry (maximum %i)",
  5.2091 +            PGL_CLUSTER_MAXPOINTS
  5.2092 +          )
  5.2093 +        ));
  5.2094 +      }
  5.2095 +      /* check if pgl_point array needs to grow */
  5.2096 +      if (npoints == points_buflen) {
  5.2097 +        pgl_point *newbuf;
  5.2098 +        points_buflen *= 2;
  5.2099 +        newbuf = palloc(points_buflen * sizeof(pgl_point));
  5.2100 +        memcpy(newbuf, points, npoints * sizeof(pgl_point));
  5.2101 +        pfree(points);
  5.2102 +        points = newbuf;
  5.2103 +      }
  5.2104 +      /* append point to pgl_point array (includes checks) */
  5.2105 +      pgl_epoint_set_latlon(&(points[npoints++]), lat, lon);
  5.2106 +      /* increase total number of points */
  5.2107 +      npoints_total++;
  5.2108 +    }
  5.2109 +    /* error if entry has no points */
  5.2110 +    if (!npoints) goto pgl_ecluster_in_error;
  5.2111 +    /* entries with one point are automatically of type "point" */
  5.2112 +    if (npoints == 1) entrytype = PGL_ENTRY_POINT;
  5.2113 +    /* if entries have more than one point */
  5.2114 +    else {
  5.2115 +      /* throw error if entry type is "point" */
  5.2116 +      if (entrytype == PGL_ENTRY_POINT) {
  5.2117 +        ereport(ERROR, (
  5.2118 +          errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  5.2119 +          errmsg("invalid input syntax for type ecluster (point entry with more than one point)")
  5.2120 +        ));
  5.2121 +      }
  5.2122 +      /* coerce outlines and polygons with more than 2 points to be a path */
  5.2123 +      if (npoints == 2) entrytype = PGL_ENTRY_PATH;
  5.2124 +    }
  5.2125 +    /* append entry to pgl_newentry array */
  5.2126 +    entries[nentries].entrytype = entrytype;
  5.2127 +    entries[nentries].npoints = npoints;
  5.2128 +    entries[nentries].points = points;
  5.2129 +    nentries++;
  5.2130 +    /* consume closing parenthesis */
  5.2131 +    strptr++;
  5.2132 +    /* consume white-space */
  5.2133 +    while (isspace(strptr[0])) strptr++;
  5.2134 +  }
  5.2135 +  /* free lower case string */
  5.2136 +  pfree(str_lower);
  5.2137 +  /* create cluster from pgl_newentry array */
  5.2138 +  cluster = pgl_new_cluster(nentries, entries);
  5.2139 +  /* free pgl_newentry array */
  5.2140 +  for (i=0; i<nentries; i++) pfree(entries[i].points);
  5.2141 +  pfree(entries);
  5.2142 +  /* set bounding circle of cluster and check east/west orientation */
  5.2143 +  if (!pgl_finalize_cluster(cluster)) {
  5.2144 +    ereport(ERROR, (
  5.2145 +      errcode(ERRCODE_DATA_EXCEPTION),
  5.2146 +      errmsg("can not determine east/west orientation for ecluster"),
  5.2147 +      errhint("Ensure that each entry has a longitude span of less than 180 degrees.")
  5.2148 +    ));
  5.2149 +  }
  5.2150 +  /* return cluster */
  5.2151 +  PG_RETURN_POINTER(cluster);
  5.2152 +  /* code to throw error */
  5.2153 +  pgl_ecluster_in_error:
  5.2154 +  ereport(ERROR, (
  5.2155 +    errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  5.2156 +    errmsg("invalid input syntax for type ecluster: \"%s\"", str)
  5.2157 +  ));
  5.2158 +}
  5.2159 +
  5.2160 +/* convert point ("epoint") to string representation */
  5.2161 +PG_FUNCTION_INFO_V1(pgl_epoint_out);
  5.2162 +Datum pgl_epoint_out(PG_FUNCTION_ARGS) {
  5.2163 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  5.2164 +  char latstr[PGL_NUMBUFLEN];
  5.2165 +  char lonstr[PGL_NUMBUFLEN];
  5.2166 +  pgl_print_lat(latstr, point->lat);
  5.2167 +  pgl_print_lon(lonstr, point->lon);
  5.2168 +  PG_RETURN_CSTRING(psprintf("%s %s", latstr, lonstr));
  5.2169 +}
  5.2170 +
  5.2171 +/* convert box ("ebox") to string representation */
  5.2172 +PG_FUNCTION_INFO_V1(pgl_ebox_out);
  5.2173 +Datum pgl_ebox_out(PG_FUNCTION_ARGS) {
  5.2174 +  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
  5.2175 +  double lon_min = box->lon_min;
  5.2176 +  double lon_max = box->lon_max;
  5.2177 +  char lat_min_str[PGL_NUMBUFLEN];
  5.2178 +  char lat_max_str[PGL_NUMBUFLEN];
  5.2179 +  char lon_min_str[PGL_NUMBUFLEN];
  5.2180 +  char lon_max_str[PGL_NUMBUFLEN];
  5.2181 +  /* return string "empty" if box is set to be empty */
  5.2182 +  if (box->lat_min > box->lat_max) PG_RETURN_CSTRING("empty");
  5.2183 +  /* use boundaries exceeding W180 or E180 if 180th meridian is enclosed */
  5.2184 +  /* (required since pgl_box_in orders the longitude boundaries) */
  5.2185 +  if (lon_min > lon_max) {
  5.2186 +    if (lon_min + lon_max >= 0) lon_min -= 360;
  5.2187 +    else lon_max += 360;
  5.2188 +  }
  5.2189 +  /* format and return result */
  5.2190 +  pgl_print_lat(lat_min_str, box->lat_min);
  5.2191 +  pgl_print_lat(lat_max_str, box->lat_max);
  5.2192 +  pgl_print_lon(lon_min_str, lon_min);
  5.2193 +  pgl_print_lon(lon_max_str, lon_max);
  5.2194 +  PG_RETURN_CSTRING(psprintf(
  5.2195 +    "%s %s %s %s",
  5.2196 +    lat_min_str, lon_min_str, lat_max_str, lon_max_str
  5.2197 +  ));
  5.2198 +}
  5.2199 +
  5.2200 +/* convert circle ("ecircle") to string representation */
  5.2201 +PG_FUNCTION_INFO_V1(pgl_ecircle_out);
  5.2202 +Datum pgl_ecircle_out(PG_FUNCTION_ARGS) {
  5.2203 +  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
  5.2204 +  char latstr[PGL_NUMBUFLEN];
  5.2205 +  char lonstr[PGL_NUMBUFLEN];
  5.2206 +  char radstr[PGL_NUMBUFLEN];
  5.2207 +  pgl_print_lat(latstr, circle->center.lat);
  5.2208 +  pgl_print_lon(lonstr, circle->center.lon);
  5.2209 +  pgl_print_float(radstr, circle->radius);
  5.2210 +  PG_RETURN_CSTRING(psprintf("%s %s %s", latstr, lonstr, radstr));
  5.2211 +}
  5.2212 +
  5.2213 +/* convert cluster ("ecluster") to string representation */
  5.2214 +PG_FUNCTION_INFO_V1(pgl_ecluster_out);
  5.2215 +Datum pgl_ecluster_out(PG_FUNCTION_ARGS) {
  5.2216 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  5.2217 +  char latstr[PGL_NUMBUFLEN];  /* string buffer for latitude */
  5.2218 +  char lonstr[PGL_NUMBUFLEN];  /* string buffer for longitude */
  5.2219 +  char ***strings;     /* array of array of strings */
  5.2220 +  char *string;        /* string of current token */
  5.2221 +  char *res, *resptr;  /* result and pointer to current write position */
  5.2222 +  size_t reslen = 1;   /* length of result (init with 1 for terminator) */
  5.2223 +  int npoints;         /* number of points of current entry */
  5.2224 +  int i, j;            /* i: entry, j: point in entry */
  5.2225 +  /* handle empty clusters */
  5.2226 +  if (cluster->nentries == 0) {
  5.2227 +    /* free detoasted cluster (if copy) */
  5.2228 +    PG_FREE_IF_COPY(cluster, 0);
  5.2229 +    /* return static result */
  5.2230 +    PG_RETURN_CSTRING("empty");
  5.2231 +  }
  5.2232 +  /* allocate array of array of strings */
  5.2233 +  strings = palloc(cluster->nentries * sizeof(char **));
  5.2234 +  /* iterate over all entries in cluster */
  5.2235 +  for (i=0; i<cluster->nentries; i++) {
  5.2236 +    /* get number of points in entry */
  5.2237 +    npoints = cluster->entries[i].npoints;
  5.2238 +    /* allocate array of strings (one string for each point plus two extra) */
  5.2239 +    strings[i] = palloc((2 + npoints) * sizeof(char *));
  5.2240 +    /* determine opening string */
  5.2241 +    switch (cluster->entries[i].entrytype) {
  5.2242 +      case PGL_ENTRY_POINT:   string = (i==0)?"point ("  :" point (";   break;
  5.2243 +      case PGL_ENTRY_PATH:    string = (i==0)?"path ("   :" path (";    break;
  5.2244 +      case PGL_ENTRY_OUTLINE: string = (i==0)?"outline (":" outline ("; break;
  5.2245 +      case PGL_ENTRY_POLYGON: string = (i==0)?"polygon (":" polygon ("; break;
  5.2246 +      default:                string = (i==0)?"unknown"  :" unknown";
  5.2247 +    }
  5.2248 +    /* use opening string as first string in array */
  5.2249 +    strings[i][0] = string;
  5.2250 +    /* update result length (for allocating result string later) */
  5.2251 +    reslen += strlen(string);
  5.2252 +    /* iterate over all points */
  5.2253 +    for (j=0; j<npoints; j++) {
  5.2254 +      /* create string representation of point */
  5.2255 +      pgl_print_lat(latstr, PGL_ENTRY_POINTS(cluster, i)[j].lat);
  5.2256 +      pgl_print_lon(lonstr, PGL_ENTRY_POINTS(cluster, i)[j].lon);
  5.2257 +      string = psprintf((j == 0) ? "%s %s" : " %s %s", latstr, lonstr);
  5.2258 +      /* copy string pointer to string array */
  5.2259 +      strings[i][j+1] = string;
  5.2260 +      /* update result length (for allocating result string later) */
  5.2261 +      reslen += strlen(string);
  5.2262 +    }
  5.2263 +    /* use closing parenthesis as last string in array */
  5.2264 +    strings[i][npoints+1] = ")";
  5.2265 +    /* update result length (for allocating result string later) */
  5.2266 +    reslen++;
  5.2267 +  }
  5.2268 +  /* allocate result string */
  5.2269 +  res = palloc(reslen);
  5.2270 +  /* set write pointer to begin of result string */
  5.2271 +  resptr = res;
  5.2272 +  /* copy strings into result string */
  5.2273 +  for (i=0; i<cluster->nentries; i++) {
  5.2274 +    npoints = cluster->entries[i].npoints;
  5.2275 +    for (j=0; j<npoints+2; j++) {
  5.2276 +      string = strings[i][j];
  5.2277 +      strcpy(resptr, string);
  5.2278 +      resptr += strlen(string);
  5.2279 +      /* free strings allocated by psprintf */
  5.2280 +      if (j != 0 && j != npoints+1) pfree(string);
  5.2281 +    }
  5.2282 +    /* free array of strings */
  5.2283 +    pfree(strings[i]);
  5.2284 +  }
  5.2285 +  /* free array of array of strings */
  5.2286 +  pfree(strings);
  5.2287 +  /* free detoasted cluster (if copy) */
  5.2288 +  PG_FREE_IF_COPY(cluster, 0);
  5.2289 +  /* return result */
  5.2290 +  PG_RETURN_CSTRING(res);
  5.2291 +}
  5.2292 +
  5.2293 +/* binary input function for point ("epoint") */
  5.2294 +PG_FUNCTION_INFO_V1(pgl_epoint_recv);
  5.2295 +Datum pgl_epoint_recv(PG_FUNCTION_ARGS) {
  5.2296 +  StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
  5.2297 +  pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point));
  5.2298 +  point->lat = pq_getmsgfloat8(buf);
  5.2299 +  point->lon = pq_getmsgfloat8(buf);
  5.2300 +  PG_RETURN_POINTER(point);
  5.2301 +}
  5.2302 +
  5.2303 +/* binary input function for box ("ebox") */
  5.2304 +PG_FUNCTION_INFO_V1(pgl_ebox_recv);
  5.2305 +Datum pgl_ebox_recv(PG_FUNCTION_ARGS) {
  5.2306 +  StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
  5.2307 +  pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
  5.2308 +  box->lat_min = pq_getmsgfloat8(buf);
  5.2309 +  box->lat_max = pq_getmsgfloat8(buf);
  5.2310 +  box->lon_min = pq_getmsgfloat8(buf);
  5.2311 +  box->lon_max = pq_getmsgfloat8(buf);
  5.2312 +  PG_RETURN_POINTER(box);
  5.2313 +}
  5.2314 +
  5.2315 +/* binary input function for circle ("ecircle") */
  5.2316 +PG_FUNCTION_INFO_V1(pgl_ecircle_recv);
  5.2317 +Datum pgl_ecircle_recv(PG_FUNCTION_ARGS) {
  5.2318 +  StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
  5.2319 +  pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
  5.2320 +  circle->center.lat = pq_getmsgfloat8(buf);
  5.2321 +  circle->center.lon = pq_getmsgfloat8(buf);
  5.2322 +  circle->radius = pq_getmsgfloat8(buf);
  5.2323 +  PG_RETURN_POINTER(circle);
  5.2324 +}
  5.2325 +
  5.2326 +/* TODO: binary receive function for cluster */
  5.2327 +
  5.2328 +/* binary output function for point ("epoint") */
  5.2329 +PG_FUNCTION_INFO_V1(pgl_epoint_send);
  5.2330 +Datum pgl_epoint_send(PG_FUNCTION_ARGS) {
  5.2331 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  5.2332 +  StringInfoData buf;
  5.2333 +  pq_begintypsend(&buf);
  5.2334 +  pq_sendfloat8(&buf, point->lat);
  5.2335 +  pq_sendfloat8(&buf, point->lon);
  5.2336 +  PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
  5.2337 +}
  5.2338 +
  5.2339 +/* binary output function for box ("ebox") */
  5.2340 +PG_FUNCTION_INFO_V1(pgl_ebox_send);
  5.2341 +Datum pgl_ebox_send(PG_FUNCTION_ARGS) {
  5.2342 +  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
  5.2343 +  StringInfoData buf;
  5.2344 +  pq_begintypsend(&buf);
  5.2345 +  pq_sendfloat8(&buf, box->lat_min);
  5.2346 +  pq_sendfloat8(&buf, box->lat_max);
  5.2347 +  pq_sendfloat8(&buf, box->lon_min);
  5.2348 +  pq_sendfloat8(&buf, box->lon_max);
  5.2349 +  PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
  5.2350 +}
  5.2351 +
  5.2352 +/* binary output function for circle ("ecircle") */
  5.2353 +PG_FUNCTION_INFO_V1(pgl_ecircle_send);
  5.2354 +Datum pgl_ecircle_send(PG_FUNCTION_ARGS) {
  5.2355 +  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
  5.2356 +  StringInfoData buf;
  5.2357 +  pq_begintypsend(&buf);
  5.2358 +  pq_sendfloat8(&buf, circle->center.lat);
  5.2359 +  pq_sendfloat8(&buf, circle->center.lon);
  5.2360 +  pq_sendfloat8(&buf, circle->radius);
  5.2361 +  PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
  5.2362 +}
  5.2363 +
  5.2364 +/* TODO: binary send functions for cluster */
  5.2365 +
  5.2366 +/* cast point ("epoint") to box ("ebox") */
  5.2367 +PG_FUNCTION_INFO_V1(pgl_epoint_to_ebox);
  5.2368 +Datum pgl_epoint_to_ebox(PG_FUNCTION_ARGS) {
  5.2369 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  5.2370 +  pgl_box *box = palloc(sizeof(pgl_box));
  5.2371 +  box->lat_min = point->lat;
  5.2372 +  box->lat_max = point->lat;
  5.2373 +  box->lon_min = point->lon;
  5.2374 +  box->lon_max = point->lon;
  5.2375 +  PG_RETURN_POINTER(box);
  5.2376 +}
  5.2377 +
  5.2378 +/* cast point ("epoint") to circle ("ecircle") */
  5.2379 +PG_FUNCTION_INFO_V1(pgl_epoint_to_ecircle);
  5.2380 +Datum pgl_epoint_to_ecircle(PG_FUNCTION_ARGS) {
  5.2381 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  5.2382 +  pgl_circle *circle = palloc(sizeof(pgl_box));
  5.2383 +  circle->center = *point;
  5.2384 +  circle->radius = 0;
  5.2385 +  PG_RETURN_POINTER(circle);
  5.2386 +}
  5.2387 +
  5.2388 +/* cast point ("epoint") to cluster ("ecluster") */
  5.2389 +PG_FUNCTION_INFO_V1(pgl_epoint_to_ecluster);
  5.2390 +Datum pgl_epoint_to_ecluster(PG_FUNCTION_ARGS) {
  5.2391 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  5.2392 +  pgl_newentry entry;
  5.2393 +  pgl_cluster *cluster;
  5.2394 +  entry.entrytype = PGL_ENTRY_POINT;
  5.2395 +  entry.npoints = 1;
  5.2396 +  entry.points = point;
  5.2397 +  cluster = pgl_new_cluster(1, &entry);
  5.2398 +  pgl_finalize_cluster(cluster);  /* NOTE: should not fail */
  5.2399 +  PG_RETURN_POINTER(cluster);
  5.2400 +}
  5.2401 +
  5.2402 +/* cast box ("ebox") to cluster ("ecluster") */
  5.2403 +#define pgl_ebox_to_ecluster_macro(i, a, b) \
  5.2404 +  entries[i].entrytype = PGL_ENTRY_POLYGON; \
  5.2405 +  entries[i].npoints = 4; \
  5.2406 +  entries[i].points = points[i]; \
  5.2407 +  points[i][0].lat = box->lat_min; \
  5.2408 +  points[i][0].lon = (a); \
  5.2409 +  points[i][1].lat = box->lat_min; \
  5.2410 +  points[i][1].lon = (b); \
  5.2411 +  points[i][2].lat = box->lat_max; \
  5.2412 +  points[i][2].lon = (b); \
  5.2413 +  points[i][3].lat = box->lat_max; \
  5.2414 +  points[i][3].lon = (a);
  5.2415 +PG_FUNCTION_INFO_V1(pgl_ebox_to_ecluster);
  5.2416 +Datum pgl_ebox_to_ecluster(PG_FUNCTION_ARGS) {
  5.2417 +  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
  5.2418 +  double lon, dlon;
  5.2419 +  int nentries;
  5.2420 +  pgl_newentry entries[3];
  5.2421 +  pgl_point points[3][4];
  5.2422 +  pgl_cluster *cluster;
  5.2423 +  if (box->lat_min > box->lat_max) {
  5.2424 +    nentries = 0;
  5.2425 +  } else if (box->lon_min > box->lon_max) {
  5.2426 +    if (box->lon_min < 0) {
  5.2427 +      lon = pgl_round((box->lon_min + 180) / 2.0);
  5.2428 +      nentries = 3;
  5.2429 +      pgl_ebox_to_ecluster_macro(0, box->lon_min, lon);
  5.2430 +      pgl_ebox_to_ecluster_macro(1, lon, 180);
  5.2431 +      pgl_ebox_to_ecluster_macro(2, -180, box->lon_max);
  5.2432 +    } else if (box->lon_max > 0) {
  5.2433 +      lon = pgl_round((box->lon_max - 180) / 2.0);
  5.2434 +      nentries = 3;
  5.2435 +      pgl_ebox_to_ecluster_macro(0, box->lon_min, 180);
  5.2436 +      pgl_ebox_to_ecluster_macro(1, -180, lon);
  5.2437 +      pgl_ebox_to_ecluster_macro(2, lon, box->lon_max);
  5.2438 +    } else {
  5.2439 +      nentries = 2;
  5.2440 +      pgl_ebox_to_ecluster_macro(0, box->lon_min, 180);
  5.2441 +      pgl_ebox_to_ecluster_macro(1, -180, box->lon_max);
  5.2442 +    }
  5.2443 +  } else {
  5.2444 +    dlon = pgl_round(box->lon_max - box->lon_min);
  5.2445 +    if (dlon < 180) {
  5.2446 +      nentries = 1;
  5.2447 +      pgl_ebox_to_ecluster_macro(0, box->lon_min, box->lon_max);
  5.2448 +    } else {
  5.2449 +      lon = pgl_round((box->lon_min + box->lon_max) / 2.0);
  5.2450 +      if (
  5.2451 +        pgl_round(lon - box->lon_min) < 180 &&
  5.2452 +        pgl_round(box->lon_max - lon) < 180
  5.2453 +      ) {
  5.2454 +        nentries = 2;
  5.2455 +        pgl_ebox_to_ecluster_macro(0, box->lon_min, lon);
  5.2456 +        pgl_ebox_to_ecluster_macro(1, lon, box->lon_max);
  5.2457 +      } else {
  5.2458 +        nentries = 3;
  5.2459 +        pgl_ebox_to_ecluster_macro(0, box->lon_min, -60);
  5.2460 +        pgl_ebox_to_ecluster_macro(1, -60, 60);
  5.2461 +        pgl_ebox_to_ecluster_macro(2, 60, box->lon_max);
  5.2462 +      }
  5.2463 +    }
  5.2464 +  }
  5.2465 +  cluster = pgl_new_cluster(nentries, entries);
  5.2466 +  pgl_finalize_cluster(cluster);  /* NOTE: should not fail */
  5.2467 +  PG_RETURN_POINTER(cluster);
  5.2468 +}
  5.2469 +
  5.2470 +/* extract latitude from point ("epoint") */
  5.2471 +PG_FUNCTION_INFO_V1(pgl_epoint_lat);
  5.2472 +Datum pgl_epoint_lat(PG_FUNCTION_ARGS) {
  5.2473 +  PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lat);
  5.2474 +}
  5.2475 +
  5.2476 +/* extract longitude from point ("epoint") */
  5.2477 +PG_FUNCTION_INFO_V1(pgl_epoint_lon);
  5.2478 +Datum pgl_epoint_lon(PG_FUNCTION_ARGS) {
  5.2479 +  PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lon);
  5.2480 +}
  5.2481 +
  5.2482 +/* extract minimum latitude from box ("ebox") */
  5.2483 +PG_FUNCTION_INFO_V1(pgl_ebox_lat_min);
  5.2484 +Datum pgl_ebox_lat_min(PG_FUNCTION_ARGS) {
  5.2485 +  PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_min);
  5.2486 +}
  5.2487 +
  5.2488 +/* extract maximum latitude from box ("ebox") */
  5.2489 +PG_FUNCTION_INFO_V1(pgl_ebox_lat_max);
  5.2490 +Datum pgl_ebox_lat_max(PG_FUNCTION_ARGS) {
  5.2491 +  PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_max);
  5.2492 +}
  5.2493 +
  5.2494 +/* extract minimum longitude from box ("ebox") */
  5.2495 +PG_FUNCTION_INFO_V1(pgl_ebox_lon_min);
  5.2496 +Datum pgl_ebox_lon_min(PG_FUNCTION_ARGS) {
  5.2497 +  PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_min);
  5.2498 +}
  5.2499 +
  5.2500 +/* extract maximum longitude from box ("ebox") */
  5.2501 +PG_FUNCTION_INFO_V1(pgl_ebox_lon_max);
  5.2502 +Datum pgl_ebox_lon_max(PG_FUNCTION_ARGS) {
  5.2503 +  PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_max);
  5.2504 +}
  5.2505 +
  5.2506 +/* extract center point from circle ("ecircle") */
  5.2507 +PG_FUNCTION_INFO_V1(pgl_ecircle_center);
  5.2508 +Datum pgl_ecircle_center(PG_FUNCTION_ARGS) {
  5.2509 +  PG_RETURN_POINTER(&(((pgl_circle *)PG_GETARG_POINTER(0))->center));
  5.2510 +}
  5.2511 +
  5.2512 +/* extract radius from circle ("ecircle") */
  5.2513 +PG_FUNCTION_INFO_V1(pgl_ecircle_radius);
  5.2514 +Datum pgl_ecircle_radius(PG_FUNCTION_ARGS) {
  5.2515 +  PG_RETURN_FLOAT8(((pgl_circle *)PG_GETARG_POINTER(0))->radius);
  5.2516 +}
  5.2517 +
  5.2518 +/* check if point is inside box (overlap operator "&&") in SQL */
  5.2519 +PG_FUNCTION_INFO_V1(pgl_epoint_ebox_overlap);
  5.2520 +Datum pgl_epoint_ebox_overlap(PG_FUNCTION_ARGS) {
  5.2521 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  5.2522 +  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(1);
  5.2523 +  PG_RETURN_BOOL(pgl_point_in_box(point, box));
  5.2524 +}
  5.2525 +
  5.2526 +/* check if point is inside circle (overlap operator "&&") in SQL */
  5.2527 +PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_overlap);
  5.2528 +Datum pgl_epoint_ecircle_overlap(PG_FUNCTION_ARGS) {
  5.2529 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  5.2530 +  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
  5.2531 +  PG_RETURN_BOOL(
  5.2532 +    pgl_distance(
  5.2533 +      point->lat, point->lon,
  5.2534 +      circle->center.lat, circle->center.lon
  5.2535 +    ) <= circle->radius
  5.2536 +  );
  5.2537 +}
  5.2538 +
  5.2539 +/* check if point is inside cluster (overlap operator "&&") in SQL */
  5.2540 +PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_overlap);
  5.2541 +Datum pgl_epoint_ecluster_overlap(PG_FUNCTION_ARGS) {
  5.2542 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  5.2543 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  5.2544 +  bool retval;
  5.2545 +  /* points outside bounding circle are always assumed to be non-overlapping
  5.2546 +     (necessary for consistent table and index scans) */
  5.2547 +  if (
  5.2548 +    pgl_distance(
  5.2549 +      point->lat, point->lon,
  5.2550 +      cluster->bounding.center.lat, cluster->bounding.center.lon
  5.2551 +    ) > cluster->bounding.radius
  5.2552 +  ) retval = false;
  5.2553 +  else retval = pgl_point_in_cluster(point, cluster, false);
  5.2554 +  PG_FREE_IF_COPY(cluster, 1);
  5.2555 +  PG_RETURN_BOOL(retval);
  5.2556 +}
  5.2557 +
  5.2558 +/* check if point may be inside cluster (lossy overl. operator "&&+") in SQL */
  5.2559 +PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_may_overlap);
  5.2560 +Datum pgl_epoint_ecluster_may_overlap(PG_FUNCTION_ARGS) {
  5.2561 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  5.2562 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  5.2563 +  bool retval = pgl_distance(
  5.2564 +    point->lat, point->lon,
  5.2565 +    cluster->bounding.center.lat, cluster->bounding.center.lon
  5.2566 +  ) <= cluster->bounding.radius;
  5.2567 +  PG_FREE_IF_COPY(cluster, 1);
  5.2568 +  PG_RETURN_BOOL(retval);
  5.2569 +}
  5.2570 +
  5.2571 +/* check if two boxes overlap (overlap operator "&&") in SQL */
  5.2572 +PG_FUNCTION_INFO_V1(pgl_ebox_overlap);
  5.2573 +Datum pgl_ebox_overlap(PG_FUNCTION_ARGS) {
  5.2574 +  pgl_box *box1 = (pgl_box *)PG_GETARG_POINTER(0);
  5.2575 +  pgl_box *box2 = (pgl_box *)PG_GETARG_POINTER(1);
  5.2576 +  PG_RETURN_BOOL(pgl_boxes_overlap(box1, box2));
  5.2577 +}
  5.2578 +
  5.2579 +/* check if box and circle may overlap (lossy overl. operator "&&+") in SQL */
  5.2580 +PG_FUNCTION_INFO_V1(pgl_ebox_ecircle_may_overlap);
  5.2581 +Datum pgl_ebox_ecircle_may_overlap(PG_FUNCTION_ARGS) {
  5.2582 +  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
  5.2583 +  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
  5.2584 +  PG_RETURN_BOOL(
  5.2585 +    pgl_estimate_point_box_distance(&circle->center, box) <= circle->radius
  5.2586 +  );
  5.2587 +}
  5.2588 +
  5.2589 +/* check if box and cluster may overlap (lossy overl. operator "&&+") in SQL */
  5.2590 +PG_FUNCTION_INFO_V1(pgl_ebox_ecluster_may_overlap);
  5.2591 +Datum pgl_ebox_ecluster_may_overlap(PG_FUNCTION_ARGS) {
  5.2592 +  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
  5.2593 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  5.2594 +  bool retval = pgl_estimate_point_box_distance(
  5.2595 +    &cluster->bounding.center,
  5.2596 +    box
  5.2597 +  ) <= cluster->bounding.radius;
  5.2598 +  PG_FREE_IF_COPY(cluster, 1);
  5.2599 +  PG_RETURN_BOOL(retval);
  5.2600 +}
  5.2601 +
  5.2602 +/* check if two circles overlap (overlap operator "&&") in SQL */
  5.2603 +PG_FUNCTION_INFO_V1(pgl_ecircle_overlap);
  5.2604 +Datum pgl_ecircle_overlap(PG_FUNCTION_ARGS) {
  5.2605 +  pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0);
  5.2606 +  pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1);
  5.2607 +  PG_RETURN_BOOL(
  5.2608 +    pgl_distance(
  5.2609 +      circle1->center.lat, circle1->center.lon,
  5.2610 +      circle2->center.lat, circle2->center.lon
  5.2611 +    ) <= circle1->radius + circle2->radius
  5.2612 +  );
  5.2613 +}
  5.2614 +
  5.2615 +/* check if circle and cluster overlap (overlap operator "&&") in SQL */
  5.2616 +PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_overlap);
  5.2617 +Datum pgl_ecircle_ecluster_overlap(PG_FUNCTION_ARGS) {
  5.2618 +  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
  5.2619 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  5.2620 +  bool retval = (
  5.2621 +    pgl_point_cluster_distance(&(circle->center), cluster) <= circle->radius
  5.2622 +  );
  5.2623 +  PG_FREE_IF_COPY(cluster, 1);
  5.2624 +  PG_RETURN_BOOL(retval);
  5.2625 +}
  5.2626 +
  5.2627 +/* check if circle and cluster may overlap (l. ov. operator "&&+") in SQL */
  5.2628 +PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_may_overlap);
  5.2629 +Datum pgl_ecircle_ecluster_may_overlap(PG_FUNCTION_ARGS) {
  5.2630 +  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
  5.2631 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  5.2632 +  bool retval = pgl_distance(
  5.2633 +    circle->center.lat, circle->center.lon,
  5.2634 +    cluster->bounding.center.lat, cluster->bounding.center.lon
  5.2635 +  ) <= circle->radius + cluster->bounding.radius;
  5.2636 +  PG_FREE_IF_COPY(cluster, 1);
  5.2637 +  PG_RETURN_BOOL(retval);
  5.2638 +}
  5.2639 +
  5.2640 +/* check if two clusters overlap (overlap operator "&&") in SQL */
  5.2641 +PG_FUNCTION_INFO_V1(pgl_ecluster_overlap);
  5.2642 +Datum pgl_ecluster_overlap(PG_FUNCTION_ARGS) {
  5.2643 +  pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  5.2644 +  pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  5.2645 +  bool retval;
  5.2646 +  /* clusters with non-touching bounding circles are always assumed to be
  5.2647 +     non-overlapping (improves performance and is necessary for consistent
  5.2648 +     table and index scans) */
  5.2649 +  if (
  5.2650 +    pgl_distance(
  5.2651 +      cluster1->bounding.center.lat, cluster1->bounding.center.lon,
  5.2652 +      cluster2->bounding.center.lat, cluster2->bounding.center.lon
  5.2653 +    ) > cluster1->bounding.radius + cluster2->bounding.radius
  5.2654 +  ) retval = false;
  5.2655 +  else retval = pgl_clusters_overlap(cluster1, cluster2);
  5.2656 +  PG_FREE_IF_COPY(cluster1, 0);
  5.2657 +  PG_FREE_IF_COPY(cluster2, 1);
  5.2658 +  PG_RETURN_BOOL(retval);
  5.2659 +}
  5.2660 +
  5.2661 +/* check if two clusters may overlap (lossy overlap operator "&&+") in SQL */
  5.2662 +PG_FUNCTION_INFO_V1(pgl_ecluster_may_overlap);
  5.2663 +Datum pgl_ecluster_may_overlap(PG_FUNCTION_ARGS) {
  5.2664 +  pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  5.2665 +  pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  5.2666 +  bool retval = pgl_distance(
  5.2667 +    cluster1->bounding.center.lat, cluster1->bounding.center.lon,
  5.2668 +    cluster2->bounding.center.lat, cluster2->bounding.center.lon
  5.2669 +  ) <= cluster1->bounding.radius + cluster2->bounding.radius;
  5.2670 +  PG_FREE_IF_COPY(cluster1, 0);
  5.2671 +  PG_FREE_IF_COPY(cluster2, 1);
  5.2672 +  PG_RETURN_BOOL(retval);
  5.2673 +}
  5.2674 +
  5.2675 +/* check if second cluster is in first cluster (cont. operator "@>) in SQL */
  5.2676 +PG_FUNCTION_INFO_V1(pgl_ecluster_contains);
  5.2677 +Datum pgl_ecluster_contains(PG_FUNCTION_ARGS) {
  5.2678 +  pgl_cluster *outer = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  5.2679 +  pgl_cluster *inner = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  5.2680 +  bool retval;
  5.2681 +  /* clusters with non-touching bounding circles are always assumed to be
  5.2682 +     non-overlapping (improves performance and is necessary for consistent
  5.2683 +     table and index scans) */
  5.2684 +  if (
  5.2685 +    pgl_distance(
  5.2686 +      outer->bounding.center.lat, outer->bounding.center.lon,
  5.2687 +      inner->bounding.center.lat, inner->bounding.center.lon
  5.2688 +    ) > outer->bounding.radius + inner->bounding.radius
  5.2689 +  ) retval = false;
  5.2690 +  else retval = pgl_cluster_in_cluster(outer, inner);
  5.2691 +  PG_FREE_IF_COPY(outer, 0);
  5.2692 +  PG_FREE_IF_COPY(inner, 1);
  5.2693 +  PG_RETURN_BOOL(retval);
  5.2694 +}
  5.2695 +
  5.2696 +/* calculate distance between two points ("<->" operator) in SQL */
  5.2697 +PG_FUNCTION_INFO_V1(pgl_epoint_distance);
  5.2698 +Datum pgl_epoint_distance(PG_FUNCTION_ARGS) {
  5.2699 +  pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0);
  5.2700 +  pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1);
  5.2701 +  PG_RETURN_FLOAT8(pgl_distance(
  5.2702 +    point1->lat, point1->lon, point2->lat, point2->lon
  5.2703 +  ));
  5.2704 +}
  5.2705 +
  5.2706 +/* calculate point to circle distance ("<->" operator) in SQL */
  5.2707 +PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_distance);
  5.2708 +Datum pgl_epoint_ecircle_distance(PG_FUNCTION_ARGS) {
  5.2709 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  5.2710 +  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
  5.2711 +  double distance = pgl_distance(
  5.2712 +    point->lat, point->lon, circle->center.lat, circle->center.lon
  5.2713 +  ) - circle->radius;
  5.2714 +  PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
  5.2715 +}
  5.2716 +
  5.2717 +/* calculate point to cluster distance ("<->" operator) in SQL */
  5.2718 +PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_distance);
  5.2719 +Datum pgl_epoint_ecluster_distance(PG_FUNCTION_ARGS) {
  5.2720 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  5.2721 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  5.2722 +  double distance = pgl_point_cluster_distance(point, cluster);
  5.2723 +  PG_FREE_IF_COPY(cluster, 1);
  5.2724 +  PG_RETURN_FLOAT8(distance);
  5.2725 +}
  5.2726 +
  5.2727 +/* calculate distance between two circles ("<->" operator) in SQL */
  5.2728 +PG_FUNCTION_INFO_V1(pgl_ecircle_distance);
  5.2729 +Datum pgl_ecircle_distance(PG_FUNCTION_ARGS) {
  5.2730 +  pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0);
  5.2731 +  pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1);
  5.2732 +  double distance = pgl_distance(
  5.2733 +    circle1->center.lat, circle1->center.lon,
  5.2734 +    circle2->center.lat, circle2->center.lon
  5.2735 +  ) - (circle1->radius + circle2->radius);
  5.2736 +  PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
  5.2737 +}
  5.2738 +
  5.2739 +/* calculate circle to cluster distance ("<->" operator) in SQL */
  5.2740 +PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_distance);
  5.2741 +Datum pgl_ecircle_ecluster_distance(PG_FUNCTION_ARGS) {
  5.2742 +  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
  5.2743 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  5.2744 +  double distance = (
  5.2745 +    pgl_point_cluster_distance(&(circle->center), cluster) - circle->radius
  5.2746 +  );
  5.2747 +  PG_FREE_IF_COPY(cluster, 1);
  5.2748 +  PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
  5.2749 +}
  5.2750 +
  5.2751 +/* calculate distance between two clusters ("<->" operator) in SQL */
  5.2752 +PG_FUNCTION_INFO_V1(pgl_ecluster_distance);
  5.2753 +Datum pgl_ecluster_distance(PG_FUNCTION_ARGS) {
  5.2754 +  pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  5.2755 +  pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  5.2756 +  double retval = pgl_cluster_distance(cluster1, cluster2);
  5.2757 +  PG_FREE_IF_COPY(cluster1, 0);
  5.2758 +  PG_FREE_IF_COPY(cluster2, 1);
  5.2759 +  PG_RETURN_FLOAT8(retval);
  5.2760 +}
  5.2761 +
  5.2762 +PG_FUNCTION_INFO_V1(pgl_ecluster_monte_carlo_area);
  5.2763 +Datum pgl_ecluster_monte_carlo_area(PG_FUNCTION_ARGS) {
  5.2764 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  5.2765 +  int samples = PG_GETARG_INT32(1);
  5.2766 +  double retval = pgl_monte_carlo_area(cluster, samples);
  5.2767 +  PG_FREE_IF_COPY(cluster, 0);
  5.2768 +  PG_RETURN_FLOAT8(retval);
  5.2769 +}
  5.2770 +
  5.2771 +PG_FUNCTION_INFO_V1(pgl_ecluster_epoint_fair_distance);
  5.2772 +Datum pgl_ecluster_epoint_fair_distance(PG_FUNCTION_ARGS) {
  5.2773 +  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  5.2774 +  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(1);
  5.2775 +  int samples = PG_GETARG_INT32(2);
  5.2776 +  double retval = pgl_fair_distance(point, cluster, samples);
  5.2777 +  PG_FREE_IF_COPY(cluster, 0);
  5.2778 +  PG_RETURN_FLOAT8(retval);
  5.2779 +}
  5.2780 +
  5.2781 +
  5.2782 +/*-----------------------------------------------------------*
  5.2783 + *  B-tree comparison operators and index support functions  *
  5.2784 + *-----------------------------------------------------------*/
  5.2785 +
  5.2786 +/* macro for a B-tree operator (without detoasting) */
  5.2787 +#define PGL_BTREE_OPER(func, type, cmpfunc, oper) \
  5.2788 +  PG_FUNCTION_INFO_V1(func); \
  5.2789 +  Datum func(PG_FUNCTION_ARGS) { \
  5.2790 +    type *a = (type *)PG_GETARG_POINTER(0); \
  5.2791 +    type *b = (type *)PG_GETARG_POINTER(1); \
  5.2792 +    PG_RETURN_BOOL(cmpfunc(a, b) oper 0); \
  5.2793 +  }
  5.2794 +
  5.2795 +/* macro for a B-tree comparison function (without detoasting) */
  5.2796 +#define PGL_BTREE_CMP(func, type, cmpfunc) \
  5.2797 +  PG_FUNCTION_INFO_V1(func); \
  5.2798 +  Datum func(PG_FUNCTION_ARGS) { \
  5.2799 +    type *a = (type *)PG_GETARG_POINTER(0); \
  5.2800 +    type *b = (type *)PG_GETARG_POINTER(1); \
  5.2801 +    PG_RETURN_INT32(cmpfunc(a, b)); \
  5.2802 +  }
  5.2803 +
  5.2804 +/* macro for a B-tree operator (with detoasting) */
  5.2805 +#define PGL_BTREE_OPER_DETOAST(func, type, cmpfunc, oper) \
  5.2806 +  PG_FUNCTION_INFO_V1(func); \
  5.2807 +  Datum func(PG_FUNCTION_ARGS) { \
  5.2808 +    bool res; \
  5.2809 +    type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \
  5.2810 +    type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \
  5.2811 +    res = cmpfunc(a, b) oper 0; \
  5.2812 +    PG_FREE_IF_COPY(a, 0); \
  5.2813 +    PG_FREE_IF_COPY(b, 1); \
  5.2814 +    PG_RETURN_BOOL(res); \
  5.2815 +  }
  5.2816 +
  5.2817 +/* macro for a B-tree comparison function (with detoasting) */
  5.2818 +#define PGL_BTREE_CMP_DETOAST(func, type, cmpfunc) \
  5.2819 +  PG_FUNCTION_INFO_V1(func); \
  5.2820 +  Datum func(PG_FUNCTION_ARGS) { \
  5.2821 +    int32_t res; \
  5.2822 +    type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \
  5.2823 +    type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \
  5.2824 +    res = cmpfunc(a, b); \
  5.2825 +    PG_FREE_IF_COPY(a, 0); \
  5.2826 +    PG_FREE_IF_COPY(b, 1); \
  5.2827 +    PG_RETURN_INT32(res); \
  5.2828 +  }
  5.2829 +
  5.2830 +/* B-tree operators and comparison function for point */
  5.2831 +PGL_BTREE_OPER(pgl_btree_epoint_lt, pgl_point, pgl_point_cmp, <)
  5.2832 +PGL_BTREE_OPER(pgl_btree_epoint_le, pgl_point, pgl_point_cmp, <=)
  5.2833 +PGL_BTREE_OPER(pgl_btree_epoint_eq, pgl_point, pgl_point_cmp, ==)
  5.2834 +PGL_BTREE_OPER(pgl_btree_epoint_ne, pgl_point, pgl_point_cmp, !=)
  5.2835 +PGL_BTREE_OPER(pgl_btree_epoint_ge, pgl_point, pgl_point_cmp, >=)
  5.2836 +PGL_BTREE_OPER(pgl_btree_epoint_gt, pgl_point, pgl_point_cmp, >)
  5.2837 +PGL_BTREE_CMP(pgl_btree_epoint_cmp, pgl_point, pgl_point_cmp)
  5.2838 +
  5.2839 +/* B-tree operators and comparison function for box */
  5.2840 +PGL_BTREE_OPER(pgl_btree_ebox_lt, pgl_box, pgl_box_cmp, <)
  5.2841 +PGL_BTREE_OPER(pgl_btree_ebox_le, pgl_box, pgl_box_cmp, <=)
  5.2842 +PGL_BTREE_OPER(pgl_btree_ebox_eq, pgl_box, pgl_box_cmp, ==)
  5.2843 +PGL_BTREE_OPER(pgl_btree_ebox_ne, pgl_box, pgl_box_cmp, !=)
  5.2844 +PGL_BTREE_OPER(pgl_btree_ebox_ge, pgl_box, pgl_box_cmp, >=)
  5.2845 +PGL_BTREE_OPER(pgl_btree_ebox_gt, pgl_box, pgl_box_cmp, >)
  5.2846 +PGL_BTREE_CMP(pgl_btree_ebox_cmp, pgl_box, pgl_box_cmp)
  5.2847 +
  5.2848 +/* B-tree operators and comparison function for circle */
  5.2849 +PGL_BTREE_OPER(pgl_btree_ecircle_lt, pgl_circle, pgl_circle_cmp, <)
  5.2850 +PGL_BTREE_OPER(pgl_btree_ecircle_le, pgl_circle, pgl_circle_cmp, <=)
  5.2851 +PGL_BTREE_OPER(pgl_btree_ecircle_eq, pgl_circle, pgl_circle_cmp, ==)
  5.2852 +PGL_BTREE_OPER(pgl_btree_ecircle_ne, pgl_circle, pgl_circle_cmp, !=)
  5.2853 +PGL_BTREE_OPER(pgl_btree_ecircle_ge, pgl_circle, pgl_circle_cmp, >=)
  5.2854 +PGL_BTREE_OPER(pgl_btree_ecircle_gt, pgl_circle, pgl_circle_cmp, >)
  5.2855 +PGL_BTREE_CMP(pgl_btree_ecircle_cmp, pgl_circle, pgl_circle_cmp)
  5.2856 +
  5.2857 +
  5.2858 +/*--------------------------------*
  5.2859 + *  GiST index support functions  *
  5.2860 + *--------------------------------*/
  5.2861 +
  5.2862 +/* GiST "consistent" support function */
  5.2863 +PG_FUNCTION_INFO_V1(pgl_gist_consistent);
  5.2864 +Datum pgl_gist_consistent(PG_FUNCTION_ARGS) {
  5.2865 +  GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  5.2866 +  pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key);
  5.2867 +  StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2);
  5.2868 +  bool *recheck = (bool *)PG_GETARG_POINTER(4);
  5.2869 +  /* demand recheck because index and query methods are lossy */
  5.2870 +  *recheck = true;
  5.2871 +  /* strategy number aliases for different operators using the same strategy */
  5.2872 +  strategy %= 100;
  5.2873 +  /* strategy number 11: equality of two points */
  5.2874 +  if (strategy == 11) {
  5.2875 +    /* query datum is another point */
  5.2876 +    pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
  5.2877 +    /* convert other point to key */
  5.2878 +    pgl_pointkey querykey;
  5.2879 +    pgl_point_to_key(query, querykey);
  5.2880 +    /* return true if both keys overlap */
  5.2881 +    PG_RETURN_BOOL(pgl_keys_overlap(key, querykey));
  5.2882 +  }
  5.2883 +  /* strategy number 13: equality of two circles */
  5.2884 +  if (strategy == 13) {
  5.2885 +    /* query datum is another circle */
  5.2886 +    pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
  5.2887 +    /* convert other circle to key */
  5.2888 +    pgl_areakey querykey;
  5.2889 +    pgl_circle_to_key(query, querykey);
  5.2890 +    /* return true if both keys overlap */
  5.2891 +    PG_RETURN_BOOL(pgl_keys_overlap(key, querykey));
  5.2892 +  }
  5.2893 +  /* for all remaining strategies, keys on empty objects produce no match */
  5.2894 +  /* (check necessary because query radius may be infinite) */
  5.2895 +  if (PGL_KEY_IS_EMPTY(key)) PG_RETURN_BOOL(false);
  5.2896 +  /* strategy number 21: overlapping with point */
  5.2897 +  if (strategy == 21) {
  5.2898 +    /* query datum is a point */
  5.2899 +    pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
  5.2900 +    /* return true if estimated distance (allowed to be smaller than real
  5.2901 +       distance) between index key and point is zero */
  5.2902 +    PG_RETURN_BOOL(pgl_estimate_key_distance(key, query) == 0);
  5.2903 +  }
  5.2904 +  /* strategy number 22: (point) overlapping with box */
  5.2905 +  if (strategy == 22) {
  5.2906 +    /* query datum is a box */
  5.2907 +    pgl_box *query = (pgl_box *)PG_GETARG_POINTER(1);
  5.2908 +    /* determine bounding box of indexed key */
  5.2909 +    pgl_box keybox;
  5.2910 +    pgl_key_to_box(key, &keybox);
  5.2911 +    /* return true if query box overlaps with bounding box of indexed key */
  5.2912 +    PG_RETURN_BOOL(pgl_boxes_overlap(query, &keybox));
  5.2913 +  }
  5.2914 +  /* strategy number 23: overlapping with circle */
  5.2915 +  if (strategy == 23) {
  5.2916 +    /* query datum is a circle */
  5.2917 +    pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
  5.2918 +    /* return true if estimated distance (allowed to be smaller than real
  5.2919 +       distance) between index key and circle center is smaller than radius */
  5.2920 +    PG_RETURN_BOOL(
  5.2921 +      pgl_estimate_key_distance(key, &(query->center)) <= query->radius
  5.2922 +    );
  5.2923 +  }
  5.2924 +  /* strategy number 24: overlapping with cluster */
  5.2925 +  if (strategy == 24) {
  5.2926 +    bool retval;  /* return value */
  5.2927 +    /* query datum is a cluster */
  5.2928 +    pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  5.2929 +    /* return true if estimated distance (allowed to be smaller than real
  5.2930 +       distance) between index key and circle center is smaller than radius */
  5.2931 +    retval = (
  5.2932 +      pgl_estimate_key_distance(key, &(query->bounding.center)) <=
  5.2933 +      query->bounding.radius
  5.2934 +    );
  5.2935 +    PG_FREE_IF_COPY(query, 1);  /* free detoasted cluster (if copy) */
  5.2936 +    PG_RETURN_BOOL(retval);
  5.2937 +  }
  5.2938 +  /* throw error for any unknown strategy number */
  5.2939 +  elog(ERROR, "unrecognized strategy number: %d", strategy);
  5.2940 +}
  5.2941 +
  5.2942 +/* GiST "union" support function */
  5.2943 +PG_FUNCTION_INFO_V1(pgl_gist_union);
  5.2944 +Datum pgl_gist_union(PG_FUNCTION_ARGS) {
  5.2945 +  GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0);
  5.2946 +  pgl_keyptr out;  /* return value (to be palloc'ed) */
  5.2947 +  int i;
  5.2948 +  /* determine key size */
  5.2949 +  size_t keysize = PGL_KEY_IS_AREAKEY(
  5.2950 +    (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key)
  5.2951 +  ) ? sizeof (pgl_areakey) : sizeof(pgl_pointkey);
  5.2952 +  /* begin with first key as result */
  5.2953 +  out = palloc(keysize);
  5.2954 +  memcpy(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key), keysize);
  5.2955 +  /* unite current result with second, third, etc. key */
  5.2956 +  for (i=1; i<entryvec->n; i++) {
  5.2957 +    pgl_unite_keys(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key));
  5.2958 +  }
  5.2959 +  /* return result */
  5.2960 +  PG_RETURN_POINTER(out);
  5.2961 +}
  5.2962 +
  5.2963 +/* GiST "compress" support function for indicis on points */
  5.2964 +PG_FUNCTION_INFO_V1(pgl_gist_compress_epoint);
  5.2965 +Datum pgl_gist_compress_epoint(PG_FUNCTION_ARGS) {
  5.2966 +  GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  5.2967 +  GISTENTRY *retval;  /* return value (to be palloc'ed unless set to entry) */
  5.2968 +  /* only transform new leaves */
  5.2969 +  if (entry->leafkey) {
  5.2970 +    /* get point to be transformed */
  5.2971 +    pgl_point *point = (pgl_point *)DatumGetPointer(entry->key);
  5.2972 +    /* allocate memory for key */
  5.2973 +    pgl_keyptr key = palloc(sizeof(pgl_pointkey));
  5.2974 +    /* transform point to key */
  5.2975 +    pgl_point_to_key(point, key);
  5.2976 +    /* create new GISTENTRY structure as return value */
  5.2977 +    retval = palloc(sizeof(GISTENTRY));
  5.2978 +    gistentryinit(
  5.2979 +      *retval, PointerGetDatum(key),
  5.2980 +      entry->rel, entry->page, entry->offset, FALSE
  5.2981 +    );
  5.2982 +  } else {
  5.2983 +    /* inner nodes have already been transformed */
  5.2984 +    retval = entry;
  5.2985 +  }
  5.2986 +  /* return pointer to old or new GISTENTRY structure */
  5.2987 +  PG_RETURN_POINTER(retval);
  5.2988 +}
  5.2989 +
  5.2990 +/* GiST "compress" support function for indicis on circles */
  5.2991 +PG_FUNCTION_INFO_V1(pgl_gist_compress_ecircle);
  5.2992 +Datum pgl_gist_compress_ecircle(PG_FUNCTION_ARGS) {
  5.2993 +  GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  5.2994 +  GISTENTRY *retval;  /* return value (to be palloc'ed unless set to entry) */
  5.2995 +  /* only transform new leaves */
  5.2996 +  if (entry->leafkey) {
  5.2997 +    /* get circle to be transformed */
  5.2998 +    pgl_circle *circle = (pgl_circle *)DatumGetPointer(entry->key);
  5.2999 +    /* allocate memory for key */
  5.3000 +    pgl_keyptr key = palloc(sizeof(pgl_areakey));
  5.3001 +    /* transform circle to key */
  5.3002 +    pgl_circle_to_key(circle, key);
  5.3003 +    /* create new GISTENTRY structure as return value */
  5.3004 +    retval = palloc(sizeof(GISTENTRY));
  5.3005 +    gistentryinit(
  5.3006 +      *retval, PointerGetDatum(key),
  5.3007 +      entry->rel, entry->page, entry->offset, FALSE
  5.3008 +    );
  5.3009 +  } else {
  5.3010 +    /* inner nodes have already been transformed */
  5.3011 +    retval = entry;
  5.3012 +  }
  5.3013 +  /* return pointer to old or new GISTENTRY structure */
  5.3014 +  PG_RETURN_POINTER(retval);
  5.3015 +}
  5.3016 +
  5.3017 +/* GiST "compress" support function for indices on clusters */
  5.3018 +PG_FUNCTION_INFO_V1(pgl_gist_compress_ecluster);
  5.3019 +Datum pgl_gist_compress_ecluster(PG_FUNCTION_ARGS) {
  5.3020 +  GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  5.3021 +  GISTENTRY *retval;  /* return value (to be palloc'ed unless set to entry) */
  5.3022 +  /* only transform new leaves */
  5.3023 +  if (entry->leafkey) {
  5.3024 +    /* get cluster to be transformed (detoasting necessary!) */
  5.3025 +    pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(entry->key);
  5.3026 +    /* allocate memory for key */
  5.3027 +    pgl_keyptr key = palloc(sizeof(pgl_areakey));
  5.3028 +    /* transform cluster to key */
  5.3029 +    pgl_circle_to_key(&(cluster->bounding), key);
  5.3030 +    /* create new GISTENTRY structure as return value */
  5.3031 +    retval = palloc(sizeof(GISTENTRY));
  5.3032 +    gistentryinit(
  5.3033 +      *retval, PointerGetDatum(key),
  5.3034 +      entry->rel, entry->page, entry->offset, FALSE
  5.3035 +    );
  5.3036 +    /* free detoasted datum */
  5.3037 +    if ((void *)cluster != (void *)DatumGetPointer(entry->key)) pfree(cluster);
  5.3038 +  } else {
  5.3039 +    /* inner nodes have already been transformed */
  5.3040 +    retval = entry;
  5.3041 +  }
  5.3042 +  /* return pointer to old or new GISTENTRY structure */
  5.3043 +  PG_RETURN_POINTER(retval);
  5.3044 +}
  5.3045 +
  5.3046 +/* GiST "decompress" support function for indices */
  5.3047 +PG_FUNCTION_INFO_V1(pgl_gist_decompress);
  5.3048 +Datum pgl_gist_decompress(PG_FUNCTION_ARGS) {
  5.3049 +  /* return passed pointer without transformation */
  5.3050 +  PG_RETURN_POINTER(PG_GETARG_POINTER(0));
  5.3051 +}
  5.3052 +
  5.3053 +/* GiST "penalty" support function */
  5.3054 +PG_FUNCTION_INFO_V1(pgl_gist_penalty);
  5.3055 +Datum pgl_gist_penalty(PG_FUNCTION_ARGS) {
  5.3056 +  GISTENTRY *origentry = (GISTENTRY *)PG_GETARG_POINTER(0);
  5.3057 +  GISTENTRY *newentry = (GISTENTRY *)PG_GETARG_POINTER(1);
  5.3058 +  float *penalty = (float *)PG_GETARG_POINTER(2);
  5.3059 +  /* get original key and key to insert */
  5.3060 +  pgl_keyptr orig = (pgl_keyptr)DatumGetPointer(origentry->key);
  5.3061 +  pgl_keyptr new = (pgl_keyptr)DatumGetPointer(newentry->key);
  5.3062 +  /* copy original key */
  5.3063 +  union { pgl_pointkey pointkey; pgl_areakey areakey; } union_key;
  5.3064 +  if (PGL_KEY_IS_AREAKEY(orig)) {
  5.3065 +    memcpy(union_key.areakey, orig, sizeof(union_key.areakey));
  5.3066 +  } else {
  5.3067 +    memcpy(union_key.pointkey, orig, sizeof(union_key.pointkey));
  5.3068 +  }
  5.3069 +  /* calculate union of both keys */
  5.3070 +  pgl_unite_keys((pgl_keyptr)&union_key, new);
  5.3071 +  /* penalty equal to reduction of key length (logarithm of added area) */
  5.3072 +  /* (return value by setting referenced value and returning pointer) */
  5.3073 +  *penalty = (
  5.3074 +    PGL_KEY_NODEDEPTH(orig) - PGL_KEY_NODEDEPTH((pgl_keyptr)&union_key)
  5.3075 +  );
  5.3076 +  PG_RETURN_POINTER(penalty);
  5.3077 +}
  5.3078 +
  5.3079 +/* GiST "picksplit" support function */
  5.3080 +PG_FUNCTION_INFO_V1(pgl_gist_picksplit);
  5.3081 +Datum pgl_gist_picksplit(PG_FUNCTION_ARGS) {
  5.3082 +  GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0);
  5.3083 +  GIST_SPLITVEC *v = (GIST_SPLITVEC *)PG_GETARG_POINTER(1);
  5.3084 +  OffsetNumber i;  /* between FirstOffsetNumber and entryvec->n (inclusive) */
  5.3085 +  union {
  5.3086 +    pgl_pointkey pointkey;
  5.3087 +    pgl_areakey areakey;
  5.3088 +  } union_all;  /* union of all keys (to be calculated from scratch)
  5.3089 +                   (later cut in half) */
  5.3090 +  int is_areakey = PGL_KEY_IS_AREAKEY(
  5.3091 +    (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key)
  5.3092 +  );
  5.3093 +  int keysize = is_areakey ? sizeof(pgl_areakey) : sizeof(pgl_pointkey);
  5.3094 +  pgl_keyptr unionL = palloc(keysize);  /* union of keys that go left */
  5.3095 +  pgl_keyptr unionR = palloc(keysize);  /* union of keys that go right */
  5.3096 +  pgl_keyptr key;  /* current key to be processed */
  5.3097 +  /* allocate memory for array of left and right keys, set counts to zero */
  5.3098 +  v->spl_left = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber));
  5.3099 +  v->spl_nleft = 0;
  5.3100 +  v->spl_right = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber));
  5.3101 +  v->spl_nright = 0;
  5.3102 +  /* calculate union of all keys from scratch */
  5.3103 +  memcpy(
  5.3104 +    (pgl_keyptr)&union_all,
  5.3105 +    (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key),
  5.3106 +    keysize
  5.3107 +  );
  5.3108 +  for (i=FirstOffsetNumber+1; i<entryvec->n; i=OffsetNumberNext(i)) {
  5.3109 +    pgl_unite_keys(
  5.3110 +      (pgl_keyptr)&union_all,
  5.3111 +      (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key)
  5.3112 +    );
  5.3113 +  }
  5.3114 +  /* check if trivial split is necessary due to exhausted key length */
  5.3115 +  /* (Note: keys for empty objects must have node depth set to maximum) */
  5.3116 +  if (PGL_KEY_NODEDEPTH((pgl_keyptr)&union_all) == (
  5.3117 +    is_areakey ? PGL_AREAKEY_MAXDEPTH : PGL_POINTKEY_MAXDEPTH
  5.3118 +  )) {
  5.3119 +    /* half of all keys go left */
  5.3120 +    for (
  5.3121 +      i=FirstOffsetNumber;
  5.3122 +      i<FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2;
  5.3123 +      i=OffsetNumberNext(i)
  5.3124 +    ) {
  5.3125 +      /* pointer to current key */
  5.3126 +      key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
  5.3127 +      /* update unionL */
  5.3128 +      /* check if key is first key that goes left */
  5.3129 +      if (!v->spl_nleft) {
  5.3130 +        /* first key that goes left is just copied to unionL */
  5.3131 +        memcpy(unionL, key, keysize);
  5.3132 +      } else {
  5.3133 +        /* unite current value and next key */
  5.3134 +        pgl_unite_keys(unionL, key);
  5.3135 +      }
  5.3136 +      /* append offset number to list of keys that go left */
  5.3137 +      v->spl_left[v->spl_nleft++] = i;
  5.3138 +    }
  5.3139 +    /* other half goes right */
  5.3140 +    for (
  5.3141 +      i=FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2;
  5.3142 +      i<entryvec->n;
  5.3143 +      i=OffsetNumberNext(i)
  5.3144 +    ) {
  5.3145 +      /* pointer to current key */
  5.3146 +      key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
  5.3147 +      /* update unionR */
  5.3148 +      /* check if key is first key that goes right */
  5.3149 +      if (!v->spl_nright) {
  5.3150 +        /* first key that goes right is just copied to unionR */
  5.3151 +        memcpy(unionR, key, keysize);
  5.3152 +      } else {
  5.3153 +        /* unite current value and next key */
  5.3154 +        pgl_unite_keys(unionR, key);
  5.3155 +      }
  5.3156 +      /* append offset number to list of keys that go right */
  5.3157 +      v->spl_right[v->spl_nright++] = i;
  5.3158 +    }
  5.3159 +  }
  5.3160 +  /* otherwise, a non-trivial split is possible */
  5.3161 +  else {
  5.3162 +    /* cut covered area in half */
  5.3163 +    /* (union_all then refers to area of keys that go left) */
  5.3164 +    /* check if union of all keys covers empty and non-empty objects */
  5.3165 +    if (PGL_KEY_IS_UNIVERSAL((pgl_keyptr)&union_all)) {
  5.3166 +      /* if yes, split into empty and non-empty objects */
  5.3167 +      pgl_key_set_empty((pgl_keyptr)&union_all);
  5.3168 +    } else {
  5.3169 +      /* otherwise split by next bit */
  5.3170 +      ((pgl_keyptr)&union_all)[PGL_KEY_NODEDEPTH_OFFSET]++;
  5.3171 +      /* NOTE: type bit conserved */
  5.3172 +    }
  5.3173 +    /* determine for each key if it goes left or right */
  5.3174 +    for (i=FirstOffsetNumber; i<entryvec->n; i=OffsetNumberNext(i)) {
  5.3175 +      /* pointer to current key */
  5.3176 +      key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
  5.3177 +      /* keys within one half of the area go left */
  5.3178 +      if (pgl_keys_overlap((pgl_keyptr)&union_all, key)) {
  5.3179 +        /* update unionL */
  5.3180 +        /* check if key is first key that goes left */
  5.3181 +        if (!v->spl_nleft) {
  5.3182 +          /* first key that goes left is just copied to unionL */
  5.3183 +          memcpy(unionL, key, keysize);
  5.3184 +        } else {
  5.3185 +          /* unite current value of unionL and processed key */
  5.3186 +          pgl_unite_keys(unionL, key);
  5.3187 +        }
  5.3188 +        /* append offset number to list of keys that go left */
  5.3189 +        v->spl_left[v->spl_nleft++] = i;
  5.3190 +      }
  5.3191 +      /* the other keys go right */
  5.3192 +      else {
  5.3193 +        /* update unionR */
  5.3194 +        /* check if key is first key that goes right */
  5.3195 +        if (!v->spl_nright) {
  5.3196 +          /* first key that goes right is just copied to unionR */
  5.3197 +          memcpy(unionR, key, keysize);
  5.3198 +        } else {
  5.3199 +          /* unite current value of unionR and processed key */
  5.3200 +          pgl_unite_keys(unionR, key);
  5.3201 +        }
  5.3202 +        /* append offset number to list of keys that go right */
  5.3203 +        v->spl_right[v->spl_nright++] = i;
  5.3204 +      }
  5.3205 +    }
  5.3206 +  }
  5.3207 +  /* store unions in return value */
  5.3208 +  v->spl_ldatum = PointerGetDatum(unionL);
  5.3209 +  v->spl_rdatum = PointerGetDatum(unionR);
  5.3210 +  /* return all results */
  5.3211 +  PG_RETURN_POINTER(v);
  5.3212 +}
  5.3213 +
  5.3214 +/* GiST "same"/"equal" support function */
  5.3215 +PG_FUNCTION_INFO_V1(pgl_gist_same);
  5.3216 +Datum pgl_gist_same(PG_FUNCTION_ARGS) {
  5.3217 +  pgl_keyptr key1 = (pgl_keyptr)PG_GETARG_POINTER(0);
  5.3218 +  pgl_keyptr key2 = (pgl_keyptr)PG_GETARG_POINTER(1);
  5.3219 +  bool *boolptr = (bool *)PG_GETARG_POINTER(2);
  5.3220 +  /* two keys are equal if they are binary equal */
  5.3221 +  /* (return result by setting referenced boolean and returning pointer) */
  5.3222 +  *boolptr = !memcmp(
  5.3223 +    key1,
  5.3224 +    key2,
  5.3225 +    PGL_KEY_IS_AREAKEY(key1) ? sizeof(pgl_areakey) : sizeof(pgl_pointkey)
  5.3226 +  );
  5.3227 +  PG_RETURN_POINTER(boolptr);
  5.3228 +}
  5.3229 +
  5.3230 +/* GiST "distance" support function */
  5.3231 +PG_FUNCTION_INFO_V1(pgl_gist_distance);
  5.3232 +Datum pgl_gist_distance(PG_FUNCTION_ARGS) {
  5.3233 +  GISTENTRY *entry = (GISTENTRY *)PG_GETARG_POINTER(0);
  5.3234 +  pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key);
  5.3235 +  StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2);
  5.3236 +  bool *recheck = (bool *)PG_GETARG_POINTER(4);
  5.3237 +  double distance;  /* return value */
  5.3238 +  /* demand recheck because distance is just an estimation */
  5.3239 +  /* (real distance may be bigger) */
  5.3240 +  *recheck = true;
  5.3241 +  /* strategy number aliases for different operators using the same strategy */
  5.3242 +  strategy %= 100;
  5.3243 +  /* strategy number 31: distance to point */
  5.3244 +  if (strategy == 31) {
  5.3245 +    /* query datum is a point */
  5.3246 +    pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
  5.3247 +    /* use pgl_estimate_pointkey_distance() function to compute result */
  5.3248 +    distance = pgl_estimate_key_distance(key, query);
  5.3249 +    /* avoid infinity (reserved!) */
  5.3250 +    if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
  5.3251 +    /* return result */
  5.3252 +    PG_RETURN_FLOAT8(distance);
  5.3253 +  }
  5.3254 +  /* strategy number 33: distance to circle */
  5.3255 +  if (strategy == 33) {
  5.3256 +    /* query datum is a circle */
  5.3257 +    pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
  5.3258 +    /* estimate distance to circle center and substract circle radius */
  5.3259 +    distance = (
  5.3260 +      pgl_estimate_key_distance(key, &(query->center)) - query->radius
  5.3261 +    );
  5.3262 +    /* convert non-positive values to zero and avoid infinity (reserved!) */
  5.3263 +    if (distance <= 0) distance = 0;
  5.3264 +    else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
  5.3265 +    /* return result */
  5.3266 +    PG_RETURN_FLOAT8(distance);
  5.3267 +  }
  5.3268 +  /* strategy number 34: distance to cluster */
  5.3269 +  if (strategy == 34) {
  5.3270 +    /* query datum is a cluster */
  5.3271 +    pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  5.3272 +    /* estimate distance to bounding center and substract bounding radius */
  5.3273 +    distance = (
  5.3274 +      pgl_estimate_key_distance(key, &(query->bounding.center)) -
  5.3275 +      query->bounding.radius
  5.3276 +    );
  5.3277 +    /* convert non-positive values to zero and avoid infinity (reserved!) */
  5.3278 +    if (distance <= 0) distance = 0;
  5.3279 +    else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
  5.3280 +    /* free detoasted cluster (if copy) */
  5.3281 +    PG_FREE_IF_COPY(query, 1);
  5.3282 +    /* return result */
  5.3283 +    PG_RETURN_FLOAT8(distance);
  5.3284 +  }
  5.3285 +  /* throw error for any unknown strategy number */
  5.3286 +  elog(ERROR, "unrecognized strategy number: %d", strategy);
  5.3287 +}
  5.3288 +

     6.1 --- a/latlon.control	Fri Oct 21 12:57:46 2016 +0200
     6.2 +++ b/latlon.control	Tue Oct 25 18:44:43 2016 +0200
     6.3 @@ -1,3 +1,3 @@
     6.4  comment = 'Geographic data types and spatial indexing for the WGS-84 spheroid'
     6.5 -default_version = '0.8'
     6.6 +default_version = '0.9'
     6.7  relocatable = true
author	jbe
date	Tue Oct 25 18:44:43 2016 +0200 (2016-10-25)
parents	a90c0d750405
children	3293906ea6a0
files	GNUmakefile README.html README.mkd latlon--0.9.sql latlon-v0008.c latlon.control