1.1 --- a/GNUmakefile	Mon Oct 31 13:06:31 2016 +0100
     1.2 +++ b/GNUmakefile	Mon Oct 31 13:10:25 2016 +0100
     1.3 @@ -1,6 +1,6 @@
     1.4  EXTENSION = latlon
     1.5 -DATA = latlon--0.9.sql latlon--0.9--0.10.sql latlon--0.10.sql
     1.6 -MODULES = latlon-v0008 latlon-v0009
     1.7 +DATA = latlon--0.9--0.10.sql latlon--0.10.sql
     1.8 +MODULES = latlon-v0009
     1.9  
    1.10  PG_CONFIG = pg_config
    1.11  PGXS := $(shell $(PG_CONFIG) --pgxs)

     2.1 --- a/latlon--0.9.sql	Mon Oct 31 13:06:31 2016 +0100
     2.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.3 @@ -1,1692 +0,0 @@
     2.4 -
     2.5 -----------------------------------------
     2.6 --- forward declarations (shell types) --
     2.7 -----------------------------------------
     2.8 -
     2.9 -CREATE TYPE epoint;
    2.10 -CREATE TYPE ebox;
    2.11 -CREATE TYPE ecircle;
    2.12 -CREATE TYPE ecluster;
    2.13 -
    2.14 -
    2.15 -------------------------------------------------------------
    2.16 --- dummy input/output functions for dummy index key types --
    2.17 -------------------------------------------------------------
    2.18 -
    2.19 -CREATE FUNCTION ekey_point_in_dummy(cstring)
    2.20 -  RETURNS ekey_point
    2.21 -  LANGUAGE C IMMUTABLE STRICT
    2.22 -  AS '$libdir/latlon-v0008', 'pgl_notimpl';
    2.23 -
    2.24 -CREATE FUNCTION ekey_point_out_dummy(ekey_point)
    2.25 -  RETURNS cstring
    2.26 -  LANGUAGE C IMMUTABLE STRICT
    2.27 -  AS '$libdir/latlon-v0008', 'pgl_notimpl';
    2.28 -
    2.29 -CREATE FUNCTION ekey_area_in_dummy(cstring)
    2.30 -  RETURNS ekey_area
    2.31 -  LANGUAGE C IMMUTABLE STRICT
    2.32 -  AS '$libdir/latlon-v0008', 'pgl_notimpl';
    2.33 -
    2.34 -CREATE FUNCTION ekey_area_out_dummy(ekey_area)
    2.35 -  RETURNS cstring
    2.36 -  LANGUAGE C IMMUTABLE STRICT
    2.37 -  AS '$libdir/latlon-v0008', 'pgl_notimpl';
    2.38 -
    2.39 -
    2.40 ---------------------------
    2.41 --- text input functions --
    2.42 ---------------------------
    2.43 -
    2.44 -CREATE FUNCTION epoint_in(cstring)
    2.45 -  RETURNS epoint
    2.46 -  LANGUAGE C IMMUTABLE STRICT
    2.47 -  AS '$libdir/latlon-v0008', 'pgl_epoint_in';
    2.48 -
    2.49 -CREATE FUNCTION ebox_in(cstring)
    2.50 -  RETURNS ebox
    2.51 -  LANGUAGE C IMMUTABLE STRICT
    2.52 -  AS '$libdir/latlon-v0008', 'pgl_ebox_in';
    2.53 -
    2.54 -CREATE FUNCTION ecircle_in(cstring)
    2.55 -  RETURNS ecircle
    2.56 -  LANGUAGE C IMMUTABLE STRICT
    2.57 -  AS '$libdir/latlon-v0008', 'pgl_ecircle_in';
    2.58 -
    2.59 -CREATE FUNCTION ecluster_in(cstring)
    2.60 -  RETURNS ecluster
    2.61 -  LANGUAGE C IMMUTABLE STRICT
    2.62 -  AS '$libdir/latlon-v0008', 'pgl_ecluster_in';
    2.63 -
    2.64 -
    2.65 ----------------------------
    2.66 --- text output functions --
    2.67 ----------------------------
    2.68 -
    2.69 -CREATE FUNCTION epoint_out(epoint)
    2.70 -  RETURNS cstring
    2.71 -  LANGUAGE C IMMUTABLE STRICT
    2.72 -  AS '$libdir/latlon-v0008', 'pgl_epoint_out';
    2.73 -
    2.74 -CREATE FUNCTION ebox_out(ebox)
    2.75 -  RETURNS cstring
    2.76 -  LANGUAGE C IMMUTABLE STRICT
    2.77 -  AS '$libdir/latlon-v0008', 'pgl_ebox_out';
    2.78 -
    2.79 -CREATE FUNCTION ecircle_out(ecircle)
    2.80 -  RETURNS cstring
    2.81 -  LANGUAGE C IMMUTABLE STRICT
    2.82 -  AS '$libdir/latlon-v0008', 'pgl_ecircle_out';
    2.83 -
    2.84 -CREATE FUNCTION ecluster_out(ecluster)
    2.85 -  RETURNS cstring
    2.86 -  LANGUAGE C IMMUTABLE STRICT
    2.87 -  AS '$libdir/latlon-v0008', 'pgl_ecluster_out';
    2.88 -
    2.89 -
    2.90 ---------------------------
    2.91 --- binary I/O functions --
    2.92 ---------------------------
    2.93 -
    2.94 -CREATE FUNCTION epoint_recv(internal)
    2.95 -  RETURNS epoint
    2.96 -  LANGUAGE C IMMUTABLE STRICT
    2.97 -  AS '$libdir/latlon-v0008', 'pgl_epoint_recv';
    2.98 -
    2.99 -CREATE FUNCTION ebox_recv(internal)
   2.100 -  RETURNS ebox
   2.101 -  LANGUAGE C IMMUTABLE STRICT
   2.102 -  AS '$libdir/latlon-v0008', 'pgl_ebox_recv';
   2.103 -
   2.104 -CREATE FUNCTION ecircle_recv(internal)
   2.105 -  RETURNS ecircle
   2.106 -  LANGUAGE C IMMUTABLE STRICT
   2.107 -  AS '$libdir/latlon-v0008', 'pgl_ecircle_recv';
   2.108 -
   2.109 -CREATE FUNCTION epoint_send(epoint)
   2.110 -  RETURNS bytea
   2.111 -  LANGUAGE C IMMUTABLE STRICT
   2.112 -  AS '$libdir/latlon-v0008', 'pgl_epoint_send';
   2.113 -
   2.114 -CREATE FUNCTION ebox_send(ebox)
   2.115 -  RETURNS bytea
   2.116 -  LANGUAGE C IMMUTABLE STRICT
   2.117 -  AS '$libdir/latlon-v0008', 'pgl_ebox_send';
   2.118 -
   2.119 -CREATE FUNCTION ecircle_send(ecircle)
   2.120 -  RETURNS bytea
   2.121 -  LANGUAGE C IMMUTABLE STRICT
   2.122 -  AS '$libdir/latlon-v0008', 'pgl_ecircle_send';
   2.123 -
   2.124 -
   2.125 ------------------------------------------------
   2.126 --- type definitions of dummy index key types --
   2.127 ------------------------------------------------
   2.128 -
   2.129 -CREATE TYPE ekey_point (
   2.130 -  internallength = 8,
   2.131 -  input = ekey_point_in_dummy,
   2.132 -  output = ekey_point_out_dummy,
   2.133 -  alignment = char );
   2.134 -
   2.135 -CREATE TYPE ekey_area (
   2.136 -  internallength = 9,
   2.137 -  input = ekey_area_in_dummy,
   2.138 -  output = ekey_area_out_dummy,
   2.139 -  alignment = char );
   2.140 -
   2.141 -
   2.142 -------------------------------------------
   2.143 --- definitions of geographic data types --
   2.144 -------------------------------------------
   2.145 -
   2.146 -CREATE TYPE epoint (
   2.147 -  internallength = 16,
   2.148 -  input = epoint_in,
   2.149 -  output = epoint_out,
   2.150 -  receive = epoint_recv,
   2.151 -  send = epoint_send,
   2.152 -  alignment = double );
   2.153 -
   2.154 -CREATE TYPE ebox (
   2.155 -  internallength = 32,
   2.156 -  input = ebox_in,
   2.157 -  output = ebox_out,
   2.158 -  receive = ebox_recv,
   2.159 -  send = ebox_send,
   2.160 -  alignment = double );
   2.161 -
   2.162 -CREATE TYPE ecircle (
   2.163 -  internallength = 24,
   2.164 -  input = ecircle_in,
   2.165 -  output = ecircle_out,
   2.166 -  receive = ecircle_recv,
   2.167 -  send = ecircle_send,
   2.168 -  alignment = double );
   2.169 -
   2.170 -CREATE TYPE ecluster (
   2.171 -  internallength = VARIABLE,
   2.172 -  input = ecluster_in,
   2.173 -  output = ecluster_out,
   2.174 -  alignment = double,
   2.175 -  storage = external );
   2.176 -
   2.177 -
   2.178 ---------------------
   2.179 --- B-tree support --
   2.180 ---------------------
   2.181 -
   2.182 --- begin of B-tree support for epoint
   2.183 -
   2.184 -CREATE FUNCTION epoint_btree_lt(epoint, epoint)
   2.185 -  RETURNS boolean
   2.186 -  LANGUAGE C IMMUTABLE STRICT
   2.187 -  AS '$libdir/latlon-v0008', 'pgl_btree_epoint_lt';
   2.188 -
   2.189 -CREATE FUNCTION epoint_btree_le(epoint, epoint)
   2.190 -  RETURNS boolean
   2.191 -  LANGUAGE C IMMUTABLE STRICT
   2.192 -  AS '$libdir/latlon-v0008', 'pgl_btree_epoint_le';
   2.193 -
   2.194 -CREATE FUNCTION epoint_btree_eq(epoint, epoint)
   2.195 -  RETURNS boolean
   2.196 -  LANGUAGE C IMMUTABLE STRICT
   2.197 -  AS '$libdir/latlon-v0008', 'pgl_btree_epoint_eq';
   2.198 -
   2.199 -CREATE FUNCTION epoint_btree_ne(epoint, epoint)
   2.200 -  RETURNS boolean
   2.201 -  LANGUAGE C IMMUTABLE STRICT
   2.202 -  AS '$libdir/latlon-v0008', 'pgl_btree_epoint_ne';
   2.203 -
   2.204 -CREATE FUNCTION epoint_btree_ge(epoint, epoint)
   2.205 -  RETURNS boolean
   2.206 -  LANGUAGE C IMMUTABLE STRICT
   2.207 -  AS '$libdir/latlon-v0008', 'pgl_btree_epoint_ge';
   2.208 -
   2.209 -CREATE FUNCTION epoint_btree_gt(epoint, epoint)
   2.210 -  RETURNS boolean
   2.211 -  LANGUAGE C IMMUTABLE STRICT
   2.212 -  AS '$libdir/latlon-v0008', 'pgl_btree_epoint_gt';
   2.213 -
   2.214 -CREATE OPERATOR <<< (
   2.215 -  leftarg = epoint,
   2.216 -  rightarg = epoint,
   2.217 -  procedure = epoint_btree_lt,
   2.218 -  commutator = >>>,
   2.219 -  negator = >>>=,
   2.220 -  restrict = scalarltsel,
   2.221 -  join = scalarltjoinsel
   2.222 -);
   2.223 -
   2.224 -CREATE OPERATOR <<<= (
   2.225 -  leftarg = epoint,
   2.226 -  rightarg = epoint,
   2.227 -  procedure = epoint_btree_le,
   2.228 -  commutator = >>>=,
   2.229 -  negator = >>>,
   2.230 -  restrict = scalarltsel,
   2.231 -  join = scalarltjoinsel
   2.232 -);
   2.233 -
   2.234 -CREATE OPERATOR = (
   2.235 -  leftarg = epoint,
   2.236 -  rightarg = epoint,
   2.237 -  procedure = epoint_btree_eq,
   2.238 -  commutator = =,
   2.239 -  negator = <>,
   2.240 -  restrict = eqsel,
   2.241 -  join = eqjoinsel,
   2.242 -  merges
   2.243 -);
   2.244 -
   2.245 -CREATE OPERATOR <> (
   2.246 -  leftarg = epoint,
   2.247 -  rightarg = epoint,
   2.248 -  procedure = epoint_btree_eq,
   2.249 -  commutator = <>,
   2.250 -  negator = =,
   2.251 -  restrict = neqsel,
   2.252 -  join = neqjoinsel
   2.253 -);
   2.254 -
   2.255 -CREATE OPERATOR >>>= (
   2.256 -  leftarg = epoint,
   2.257 -  rightarg = epoint,
   2.258 -  procedure = epoint_btree_ge,
   2.259 -  commutator = <<<=,
   2.260 -  negator = <<<,
   2.261 -  restrict = scalargtsel,
   2.262 -  join = scalargtjoinsel
   2.263 -);
   2.264 -
   2.265 -CREATE OPERATOR >>> (
   2.266 -  leftarg = epoint,
   2.267 -  rightarg = epoint,
   2.268 -  procedure = epoint_btree_gt,
   2.269 -  commutator = <<<,
   2.270 -  negator = <<<=,
   2.271 -  restrict = scalargtsel,
   2.272 -  join = scalargtjoinsel
   2.273 -);
   2.274 -
   2.275 -CREATE FUNCTION epoint_btree_cmp(epoint, epoint)
   2.276 -  RETURNS int4
   2.277 -  LANGUAGE C IMMUTABLE STRICT
   2.278 -  AS '$libdir/latlon-v0008', 'pgl_btree_epoint_cmp';
   2.279 -
   2.280 -CREATE OPERATOR CLASS epoint_btree_ops
   2.281 -  DEFAULT FOR TYPE epoint USING btree AS
   2.282 -  OPERATOR 1 <<< ,
   2.283 -  OPERATOR 2 <<<= ,
   2.284 -  OPERATOR 3 = ,
   2.285 -  OPERATOR 4 >>>= ,
   2.286 -  OPERATOR 5 >>> ,
   2.287 -  FUNCTION 1 epoint_btree_cmp(epoint, epoint);
   2.288 -
   2.289 --- end of B-tree support for epoint
   2.290 -
   2.291 --- begin of B-tree support for ebox
   2.292 -
   2.293 -CREATE FUNCTION ebox_btree_lt(ebox, ebox)
   2.294 -  RETURNS boolean
   2.295 -  LANGUAGE C IMMUTABLE STRICT
   2.296 -  AS '$libdir/latlon-v0008', 'pgl_btree_ebox_lt';
   2.297 -
   2.298 -CREATE FUNCTION ebox_btree_le(ebox, ebox)
   2.299 -  RETURNS boolean
   2.300 -  LANGUAGE C IMMUTABLE STRICT
   2.301 -  AS '$libdir/latlon-v0008', 'pgl_btree_ebox_le';
   2.302 -
   2.303 -CREATE FUNCTION ebox_btree_eq(ebox, ebox)
   2.304 -  RETURNS boolean
   2.305 -  LANGUAGE C IMMUTABLE STRICT
   2.306 -  AS '$libdir/latlon-v0008', 'pgl_btree_ebox_eq';
   2.307 -
   2.308 -CREATE FUNCTION ebox_btree_ne(ebox, ebox)
   2.309 -  RETURNS boolean
   2.310 -  LANGUAGE C IMMUTABLE STRICT
   2.311 -  AS '$libdir/latlon-v0008', 'pgl_btree_ebox_ne';
   2.312 -
   2.313 -CREATE FUNCTION ebox_btree_ge(ebox, ebox)
   2.314 -  RETURNS boolean
   2.315 -  LANGUAGE C IMMUTABLE STRICT
   2.316 -  AS '$libdir/latlon-v0008', 'pgl_btree_ebox_ge';
   2.317 -
   2.318 -CREATE FUNCTION ebox_btree_gt(ebox, ebox)
   2.319 -  RETURNS boolean
   2.320 -  LANGUAGE C IMMUTABLE STRICT
   2.321 -  AS '$libdir/latlon-v0008', 'pgl_btree_ebox_gt';
   2.322 -
   2.323 -CREATE OPERATOR <<< (
   2.324 -  leftarg = ebox,
   2.325 -  rightarg = ebox,
   2.326 -  procedure = ebox_btree_lt,
   2.327 -  commutator = >>>,
   2.328 -  negator = >>>=,
   2.329 -  restrict = scalarltsel,
   2.330 -  join = scalarltjoinsel
   2.331 -);
   2.332 -
   2.333 -CREATE OPERATOR <<<= (
   2.334 -  leftarg = ebox,
   2.335 -  rightarg = ebox,
   2.336 -  procedure = ebox_btree_le,
   2.337 -  commutator = >>>=,
   2.338 -  negator = >>>,
   2.339 -  restrict = scalarltsel,
   2.340 -  join = scalarltjoinsel
   2.341 -);
   2.342 -
   2.343 -CREATE OPERATOR = (
   2.344 -  leftarg = ebox,
   2.345 -  rightarg = ebox,
   2.346 -  procedure = ebox_btree_eq,
   2.347 -  commutator = =,
   2.348 -  negator = <>,
   2.349 -  restrict = eqsel,
   2.350 -  join = eqjoinsel,
   2.351 -  merges
   2.352 -);
   2.353 -
   2.354 -CREATE OPERATOR <> (
   2.355 -  leftarg = ebox,
   2.356 -  rightarg = ebox,
   2.357 -  procedure = ebox_btree_eq,
   2.358 -  commutator = <>,
   2.359 -  negator = =,
   2.360 -  restrict = neqsel,
   2.361 -  join = neqjoinsel
   2.362 -);
   2.363 -
   2.364 -CREATE OPERATOR >>>= (
   2.365 -  leftarg = ebox,
   2.366 -  rightarg = ebox,
   2.367 -  procedure = ebox_btree_ge,
   2.368 -  commutator = <<<=,
   2.369 -  negator = <<<,
   2.370 -  restrict = scalargtsel,
   2.371 -  join = scalargtjoinsel
   2.372 -);
   2.373 -
   2.374 -CREATE OPERATOR >>> (
   2.375 -  leftarg = ebox,
   2.376 -  rightarg = ebox,
   2.377 -  procedure = ebox_btree_gt,
   2.378 -  commutator = <<<,
   2.379 -  negator = <<<=,
   2.380 -  restrict = scalargtsel,
   2.381 -  join = scalargtjoinsel
   2.382 -);
   2.383 -
   2.384 -CREATE FUNCTION ebox_btree_cmp(ebox, ebox)
   2.385 -  RETURNS int4
   2.386 -  LANGUAGE C IMMUTABLE STRICT
   2.387 -  AS '$libdir/latlon-v0008', 'pgl_btree_ebox_cmp';
   2.388 -
   2.389 -CREATE OPERATOR CLASS ebox_btree_ops
   2.390 -  DEFAULT FOR TYPE ebox USING btree AS
   2.391 -  OPERATOR 1 <<< ,
   2.392 -  OPERATOR 2 <<<= ,
   2.393 -  OPERATOR 3 = ,
   2.394 -  OPERATOR 4 >>>= ,
   2.395 -  OPERATOR 5 >>> ,
   2.396 -  FUNCTION 1 ebox_btree_cmp(ebox, ebox);
   2.397 -
   2.398 --- end of B-tree support for ebox
   2.399 -
   2.400 --- begin of B-tree support for ecircle
   2.401 -
   2.402 -CREATE FUNCTION ecircle_btree_lt(ecircle, ecircle)
   2.403 -  RETURNS boolean
   2.404 -  LANGUAGE C IMMUTABLE STRICT
   2.405 -  AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_lt';
   2.406 -
   2.407 -CREATE FUNCTION ecircle_btree_le(ecircle, ecircle)
   2.408 -  RETURNS boolean
   2.409 -  LANGUAGE C IMMUTABLE STRICT
   2.410 -  AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_le';
   2.411 -
   2.412 -CREATE FUNCTION ecircle_btree_eq(ecircle, ecircle)
   2.413 -  RETURNS boolean
   2.414 -  LANGUAGE C IMMUTABLE STRICT
   2.415 -  AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_eq';
   2.416 -
   2.417 -CREATE FUNCTION ecircle_btree_ne(ecircle, ecircle)
   2.418 -  RETURNS boolean
   2.419 -  LANGUAGE C IMMUTABLE STRICT
   2.420 -  AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_ne';
   2.421 -
   2.422 -CREATE FUNCTION ecircle_btree_ge(ecircle, ecircle)
   2.423 -  RETURNS boolean
   2.424 -  LANGUAGE C IMMUTABLE STRICT
   2.425 -  AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_ge';
   2.426 -
   2.427 -CREATE FUNCTION ecircle_btree_gt(ecircle, ecircle)
   2.428 -  RETURNS boolean
   2.429 -  LANGUAGE C IMMUTABLE STRICT
   2.430 -  AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_gt';
   2.431 -
   2.432 -CREATE OPERATOR <<< (
   2.433 -  leftarg = ecircle,
   2.434 -  rightarg = ecircle,
   2.435 -  procedure = ecircle_btree_lt,
   2.436 -  commutator = >>>,
   2.437 -  negator = >>>=,
   2.438 -  restrict = scalarltsel,
   2.439 -  join = scalarltjoinsel
   2.440 -);
   2.441 -
   2.442 -CREATE OPERATOR <<<= (
   2.443 -  leftarg = ecircle,
   2.444 -  rightarg = ecircle,
   2.445 -  procedure = ecircle_btree_le,
   2.446 -  commutator = >>>=,
   2.447 -  negator = >>>,
   2.448 -  restrict = scalarltsel,
   2.449 -  join = scalarltjoinsel
   2.450 -);
   2.451 -
   2.452 -CREATE OPERATOR = (
   2.453 -  leftarg = ecircle,
   2.454 -  rightarg = ecircle,
   2.455 -  procedure = ecircle_btree_eq,
   2.456 -  commutator = =,
   2.457 -  negator = <>,
   2.458 -  restrict = eqsel,
   2.459 -  join = eqjoinsel,
   2.460 -  merges
   2.461 -);
   2.462 -
   2.463 -CREATE OPERATOR <> (
   2.464 -  leftarg = ecircle,
   2.465 -  rightarg = ecircle,
   2.466 -  procedure = ecircle_btree_eq,
   2.467 -  commutator = <>,
   2.468 -  negator = =,
   2.469 -  restrict = neqsel,
   2.470 -  join = neqjoinsel
   2.471 -);
   2.472 -
   2.473 -CREATE OPERATOR >>>= (
   2.474 -  leftarg = ecircle,
   2.475 -  rightarg = ecircle,
   2.476 -  procedure = ecircle_btree_ge,
   2.477 -  commutator = <<<=,
   2.478 -  negator = <<<,
   2.479 -  restrict = scalargtsel,
   2.480 -  join = scalargtjoinsel
   2.481 -);
   2.482 -
   2.483 -CREATE OPERATOR >>> (
   2.484 -  leftarg = ecircle,
   2.485 -  rightarg = ecircle,
   2.486 -  procedure = ecircle_btree_gt,
   2.487 -  commutator = <<<,
   2.488 -  negator = <<<=,
   2.489 -  restrict = scalargtsel,
   2.490 -  join = scalargtjoinsel
   2.491 -);
   2.492 -
   2.493 -CREATE FUNCTION ecircle_btree_cmp(ecircle, ecircle)
   2.494 -  RETURNS int4
   2.495 -  LANGUAGE C IMMUTABLE STRICT
   2.496 -  AS '$libdir/latlon-v0008', 'pgl_btree_ecircle_cmp';
   2.497 -
   2.498 -CREATE OPERATOR CLASS ecircle_btree_ops
   2.499 -  DEFAULT FOR TYPE ecircle USING btree AS
   2.500 -  OPERATOR 1 <<< ,
   2.501 -  OPERATOR 2 <<<= ,
   2.502 -  OPERATOR 3 = ,
   2.503 -  OPERATOR 4 >>>= ,
   2.504 -  OPERATOR 5 >>> ,
   2.505 -  FUNCTION 1 ecircle_btree_cmp(ecircle, ecircle);
   2.506 -
   2.507 --- end of B-tree support for ecircle
   2.508 -
   2.509 -
   2.510 -----------------
   2.511 --- type casts --
   2.512 -----------------
   2.513 -
   2.514 -CREATE FUNCTION cast_epoint_to_ebox(epoint)
   2.515 -  RETURNS ebox
   2.516 -  LANGUAGE C IMMUTABLE STRICT
   2.517 -  AS '$libdir/latlon-v0008', 'pgl_epoint_to_ebox';
   2.518 -
   2.519 -CREATE CAST (epoint AS ebox) WITH FUNCTION cast_epoint_to_ebox(epoint);
   2.520 -
   2.521 -CREATE FUNCTION cast_epoint_to_ecircle(epoint)
   2.522 -  RETURNS ecircle
   2.523 -  LANGUAGE C IMMUTABLE STRICT
   2.524 -  AS '$libdir/latlon-v0008', 'pgl_epoint_to_ecircle';
   2.525 -
   2.526 -CREATE CAST (epoint AS ecircle) WITH FUNCTION cast_epoint_to_ecircle(epoint);
   2.527 -
   2.528 -CREATE FUNCTION cast_epoint_to_ecluster(epoint)
   2.529 -  RETURNS ecluster
   2.530 -  LANGUAGE C IMMUTABLE STRICT
   2.531 -  AS '$libdir/latlon-v0008', 'pgl_epoint_to_ecluster';
   2.532 -
   2.533 -CREATE CAST (epoint AS ecluster) WITH FUNCTION cast_epoint_to_ecluster(epoint);
   2.534 -
   2.535 -CREATE FUNCTION cast_ebox_to_ecluster(ebox)
   2.536 -  RETURNS ecluster
   2.537 -  LANGUAGE C IMMUTABLE STRICT
   2.538 -  AS '$libdir/latlon-v0008', 'pgl_ebox_to_ecluster';
   2.539 -
   2.540 -CREATE CAST (ebox AS ecluster) WITH FUNCTION cast_ebox_to_ecluster(ebox);
   2.541 -
   2.542 -
   2.543 ----------------------------
   2.544 --- constructor functions --
   2.545 ----------------------------
   2.546 -
   2.547 -CREATE FUNCTION epoint(float8, float8)
   2.548 -  RETURNS epoint
   2.549 -  LANGUAGE C IMMUTABLE STRICT
   2.550 -  AS '$libdir/latlon-v0008', 'pgl_create_epoint';
   2.551 -
   2.552 -CREATE FUNCTION epoint_latlon(float8, float8)
   2.553 -  RETURNS epoint
   2.554 -  LANGUAGE SQL IMMUTABLE STRICT AS $$
   2.555 -    SELECT epoint($1, $2)
   2.556 -  $$;
   2.557 -
   2.558 -CREATE FUNCTION epoint_lonlat(float8, float8)
   2.559 -  RETURNS epoint
   2.560 -  LANGUAGE SQL IMMUTABLE STRICT AS $$
   2.561 -    SELECT epoint($2, $1)
   2.562 -  $$;
   2.563 -
   2.564 -CREATE FUNCTION empty_ebox()
   2.565 -  RETURNS ebox
   2.566 -  LANGUAGE C IMMUTABLE STRICT
   2.567 -  AS '$libdir/latlon-v0008', 'pgl_create_empty_ebox';
   2.568 -
   2.569 -CREATE FUNCTION ebox(float8, float8, float8, float8)
   2.570 -  RETURNS ebox
   2.571 -  LANGUAGE C IMMUTABLE STRICT
   2.572 -  AS '$libdir/latlon-v0008', 'pgl_create_ebox';
   2.573 -
   2.574 -CREATE FUNCTION ebox(epoint, epoint)
   2.575 -  RETURNS ebox
   2.576 -  LANGUAGE C IMMUTABLE STRICT
   2.577 -  AS '$libdir/latlon-v0008', 'pgl_create_ebox_from_epoints';
   2.578 -
   2.579 -CREATE FUNCTION ecircle(float8, float8, float8)
   2.580 -  RETURNS ecircle
   2.581 -  LANGUAGE C IMMUTABLE STRICT
   2.582 -  AS '$libdir/latlon-v0008', 'pgl_create_ecircle';
   2.583 -
   2.584 -CREATE FUNCTION ecircle(epoint, float8)
   2.585 -  RETURNS ecircle
   2.586 -  LANGUAGE C IMMUTABLE STRICT
   2.587 -  AS '$libdir/latlon-v0008', 'pgl_create_ecircle_from_epoint';
   2.588 -
   2.589 -CREATE FUNCTION ecluster_concat(ecluster[])
   2.590 -  RETURNS ecluster
   2.591 -  LANGUAGE sql IMMUTABLE STRICT AS $$
   2.592 -    SELECT array_to_string($1, ' ')::ecluster
   2.593 -  $$;
   2.594 -
   2.595 -CREATE FUNCTION ecluster_concat(ecluster, ecluster)
   2.596 -  RETURNS ecluster
   2.597 -  LANGUAGE sql IMMUTABLE STRICT AS $$
   2.598 -    SELECT ($1::text || ' ' || $2::text)::ecluster
   2.599 -  $$;
   2.600 -
   2.601 -CREATE FUNCTION ecluster_create_multipoint(epoint[])
   2.602 -  RETURNS ecluster
   2.603 -  LANGUAGE sql IMMUTABLE STRICT AS $$
   2.604 -    SELECT
   2.605 -      array_to_string(array_agg('point (' || unnest || ')'), ' ')::ecluster
   2.606 -    FROM unnest($1)
   2.607 -  $$;
   2.608 -
   2.609 -CREATE FUNCTION ecluster_create_path(epoint[])
   2.610 -  RETURNS ecluster
   2.611 -  LANGUAGE sql IMMUTABLE STRICT AS $$
   2.612 -    SELECT CASE WHEN "str" = '' THEN 'empty'::ecluster ELSE
   2.613 -      ('path (' || array_to_string($1, ' ') || ')')::ecluster
   2.614 -    END
   2.615 -    FROM array_to_string($1, ' ') AS "str"
   2.616 -  $$;
   2.617 -
   2.618 -CREATE FUNCTION ecluster_create_outline(epoint[])
   2.619 -  RETURNS ecluster
   2.620 -  LANGUAGE sql IMMUTABLE STRICT AS $$
   2.621 -    SELECT CASE WHEN "str" = '' THEN 'empty'::ecluster ELSE
   2.622 -      ('outline (' || array_to_string($1, ' ') || ')')::ecluster
   2.623 -    END
   2.624 -    FROM array_to_string($1, ' ') AS "str"
   2.625 -  $$;
   2.626 -
   2.627 -CREATE FUNCTION ecluster_create_polygon(epoint[])
   2.628 -  RETURNS ecluster
   2.629 -  LANGUAGE sql IMMUTABLE STRICT AS $$
   2.630 -    SELECT CASE WHEN "str" = '' THEN 'empty'::ecluster ELSE
   2.631 -      ('polygon (' || array_to_string($1, ' ') || ')')::ecluster
   2.632 -    END
   2.633 -    FROM array_to_string($1, ' ') AS "str"
   2.634 -  $$;
   2.635 -
   2.636 -
   2.637 -----------------------
   2.638 --- getter functions --
   2.639 -----------------------
   2.640 -
   2.641 -CREATE FUNCTION latitude(epoint)
   2.642 -  RETURNS float8
   2.643 -  LANGUAGE C IMMUTABLE STRICT
   2.644 -  AS '$libdir/latlon-v0008', 'pgl_epoint_lat';
   2.645 -
   2.646 -CREATE FUNCTION longitude(epoint)
   2.647 -  RETURNS float8
   2.648 -  LANGUAGE C IMMUTABLE STRICT
   2.649 -  AS '$libdir/latlon-v0008', 'pgl_epoint_lon';
   2.650 -
   2.651 -CREATE FUNCTION min_latitude(ebox)
   2.652 -  RETURNS float8
   2.653 -  LANGUAGE C IMMUTABLE STRICT
   2.654 -  AS '$libdir/latlon-v0008', 'pgl_ebox_lat_min';
   2.655 -
   2.656 -CREATE FUNCTION max_latitude(ebox)
   2.657 -  RETURNS float8
   2.658 -  LANGUAGE C IMMUTABLE STRICT
   2.659 -  AS '$libdir/latlon-v0008', 'pgl_ebox_lat_max';
   2.660 -
   2.661 -CREATE FUNCTION min_longitude(ebox)
   2.662 -  RETURNS float8
   2.663 -  LANGUAGE C IMMUTABLE STRICT
   2.664 -  AS '$libdir/latlon-v0008', 'pgl_ebox_lon_min';
   2.665 -
   2.666 -CREATE FUNCTION max_longitude(ebox)
   2.667 -  RETURNS float8
   2.668 -  LANGUAGE C IMMUTABLE STRICT
   2.669 -  AS '$libdir/latlon-v0008', 'pgl_ebox_lon_max';
   2.670 -
   2.671 -CREATE FUNCTION center(ecircle)
   2.672 -  RETURNS epoint
   2.673 -  LANGUAGE C IMMUTABLE STRICT
   2.674 -  AS '$libdir/latlon-v0008', 'pgl_ecircle_center';
   2.675 -
   2.676 -CREATE FUNCTION radius(ecircle)
   2.677 -  RETURNS float8
   2.678 -  LANGUAGE C IMMUTABLE STRICT
   2.679 -  AS '$libdir/latlon-v0008', 'pgl_ecircle_radius';
   2.680 -
   2.681 -CREATE FUNCTION ecluster_extract_points(ecluster)
   2.682 -  RETURNS SETOF epoint
   2.683 -  LANGUAGE sql IMMUTABLE STRICT AS $$
   2.684 -    SELECT "match"[2]::epoint
   2.685 -    FROM regexp_matches($1::text, e'(^| )point \\(([^)]+)\\)', 'g') AS "match"
   2.686 -  $$;
   2.687 -
   2.688 -CREATE FUNCTION ecluster_extract_paths(ecluster)
   2.689 -  RETURNS SETOF epoint[]
   2.690 -  LANGUAGE sql IMMUTABLE STRICT AS $$
   2.691 -    SELECT (
   2.692 -      SELECT array_agg("m2"[1]::epoint)
   2.693 -      FROM regexp_matches("m1"[2], e'[^ ]+ [^ ]+', 'g') AS "m2"
   2.694 -    )
   2.695 -    FROM regexp_matches($1::text, e'(^| )path \\(([^)]+)\\)', 'g') AS "m1"
   2.696 -  $$;
   2.697 -
   2.698 -CREATE FUNCTION ecluster_extract_outlines(ecluster)
   2.699 -  RETURNS SETOF epoint[]
   2.700 -  LANGUAGE sql IMMUTABLE STRICT AS $$
   2.701 -    SELECT (
   2.702 -      SELECT array_agg("m2"[1]::epoint)
   2.703 -      FROM regexp_matches("m1"[2], e'[^ ]+ [^ ]+', 'g') AS "m2"
   2.704 -    )
   2.705 -    FROM regexp_matches($1::text, e'(^| )outline \\(([^)]+)\\)', 'g') AS "m1"
   2.706 -  $$;
   2.707 -
   2.708 -CREATE FUNCTION ecluster_extract_polygons(ecluster)
   2.709 -  RETURNS SETOF epoint[]
   2.710 -  LANGUAGE sql IMMUTABLE STRICT AS $$
   2.711 -    SELECT (
   2.712 -      SELECT array_agg("m2"[1]::epoint)
   2.713 -      FROM regexp_matches("m1"[2], e'[^ ]+ [^ ]+', 'g') AS "m2"
   2.714 -    )
   2.715 -    FROM regexp_matches($1::text, e'(^| )polygon \\(([^)]+)\\)', 'g') AS "m1"
   2.716 -  $$;
   2.717 -
   2.718 -
   2.719 ----------------
   2.720 --- operators --
   2.721 ----------------
   2.722 -
   2.723 -CREATE FUNCTION epoint_ebox_overlap_proc(epoint, ebox)
   2.724 -  RETURNS boolean
   2.725 -  LANGUAGE C IMMUTABLE STRICT
   2.726 -  AS '$libdir/latlon-v0008', 'pgl_epoint_ebox_overlap';
   2.727 -
   2.728 -CREATE FUNCTION epoint_ecircle_overlap_proc(epoint, ecircle)
   2.729 -  RETURNS boolean
   2.730 -  LANGUAGE C IMMUTABLE STRICT
   2.731 -  AS '$libdir/latlon-v0008', 'pgl_epoint_ecircle_overlap';
   2.732 -
   2.733 -CREATE FUNCTION epoint_ecluster_overlap_proc(epoint, ecluster)
   2.734 -  RETURNS boolean
   2.735 -  LANGUAGE C IMMUTABLE STRICT
   2.736 -  AS '$libdir/latlon-v0008', 'pgl_epoint_ecluster_overlap';
   2.737 -
   2.738 -CREATE FUNCTION epoint_ecluster_may_overlap_proc(epoint, ecluster)
   2.739 -  RETURNS boolean
   2.740 -  LANGUAGE C IMMUTABLE STRICT
   2.741 -  AS '$libdir/latlon-v0008', 'pgl_epoint_ecluster_may_overlap';
   2.742 -
   2.743 -CREATE FUNCTION ebox_overlap_proc(ebox, ebox)
   2.744 -  RETURNS boolean
   2.745 -  LANGUAGE C IMMUTABLE STRICT
   2.746 -  AS '$libdir/latlon-v0008', 'pgl_ebox_overlap';
   2.747 -
   2.748 -CREATE FUNCTION ebox_ecircle_may_overlap_proc(ebox, ecircle)
   2.749 -  RETURNS boolean
   2.750 -  LANGUAGE C IMMUTABLE STRICT
   2.751 -  AS '$libdir/latlon-v0008', 'pgl_ebox_ecircle_may_overlap';
   2.752 -
   2.753 -CREATE FUNCTION ebox_ecluster_may_overlap_proc(ebox, ecluster)
   2.754 -  RETURNS boolean
   2.755 -  LANGUAGE C IMMUTABLE STRICT
   2.756 -  AS '$libdir/latlon-v0008', 'pgl_ebox_ecluster_may_overlap';
   2.757 -
   2.758 -CREATE FUNCTION ecircle_overlap_proc(ecircle, ecircle)
   2.759 -  RETURNS boolean
   2.760 -  LANGUAGE C IMMUTABLE STRICT
   2.761 -  AS '$libdir/latlon-v0008', 'pgl_ecircle_overlap';
   2.762 -
   2.763 -CREATE FUNCTION ecircle_ecluster_overlap_proc(ecircle, ecluster)
   2.764 -  RETURNS boolean
   2.765 -  LANGUAGE C IMMUTABLE STRICT
   2.766 -  AS '$libdir/latlon-v0008', 'pgl_ecircle_ecluster_overlap';
   2.767 -
   2.768 -CREATE FUNCTION ecircle_ecluster_may_overlap_proc(ecircle, ecluster)
   2.769 -  RETURNS boolean
   2.770 -  LANGUAGE C IMMUTABLE STRICT
   2.771 -  AS '$libdir/latlon-v0008', 'pgl_ecircle_ecluster_may_overlap';
   2.772 -
   2.773 -CREATE FUNCTION ecluster_overlap_proc(ecluster, ecluster)
   2.774 -  RETURNS boolean
   2.775 -  LANGUAGE C IMMUTABLE STRICT
   2.776 -  AS '$libdir/latlon-v0008', 'pgl_ecluster_overlap';
   2.777 -
   2.778 -CREATE FUNCTION ecluster_may_overlap_proc(ecluster, ecluster)
   2.779 -  RETURNS boolean
   2.780 -  LANGUAGE C IMMUTABLE STRICT
   2.781 -  AS '$libdir/latlon-v0008', 'pgl_ecluster_may_overlap';
   2.782 -
   2.783 -CREATE FUNCTION ecluster_contains_proc(ecluster, ecluster)
   2.784 -  RETURNS boolean
   2.785 -  LANGUAGE C IMMUTABLE STRICT
   2.786 -  AS '$libdir/latlon-v0008', 'pgl_ecluster_contains';
   2.787 -
   2.788 -CREATE FUNCTION epoint_distance_proc(epoint, epoint)
   2.789 -  RETURNS float8
   2.790 -  LANGUAGE C IMMUTABLE STRICT
   2.791 -  AS '$libdir/latlon-v0008', 'pgl_epoint_distance';
   2.792 -
   2.793 -CREATE FUNCTION epoint_ecircle_distance_proc(epoint, ecircle)
   2.794 -  RETURNS float8
   2.795 -  LANGUAGE C IMMUTABLE STRICT
   2.796 -  AS '$libdir/latlon-v0008', 'pgl_epoint_ecircle_distance';
   2.797 -
   2.798 -CREATE FUNCTION epoint_ecluster_distance_proc(epoint, ecluster)
   2.799 -  RETURNS float8
   2.800 -  LANGUAGE C IMMUTABLE STRICT
   2.801 -  AS '$libdir/latlon-v0008', 'pgl_epoint_ecluster_distance';
   2.802 -
   2.803 -CREATE FUNCTION ecircle_distance_proc(ecircle, ecircle)
   2.804 -  RETURNS float8
   2.805 -  LANGUAGE C IMMUTABLE STRICT
   2.806 -  AS '$libdir/latlon-v0008', 'pgl_ecircle_distance';
   2.807 -
   2.808 -CREATE FUNCTION ecircle_ecluster_distance_proc(ecircle, ecluster)
   2.809 -  RETURNS float8
   2.810 -  LANGUAGE C IMMUTABLE STRICT
   2.811 -  AS '$libdir/latlon-v0008', 'pgl_ecircle_ecluster_distance';
   2.812 -
   2.813 -CREATE FUNCTION ecluster_distance_proc(ecluster, ecluster)
   2.814 -  RETURNS float8
   2.815 -  LANGUAGE C IMMUTABLE STRICT
   2.816 -  AS '$libdir/latlon-v0008', 'pgl_ecluster_distance';
   2.817 -
   2.818 -CREATE OPERATOR && (
   2.819 -  leftarg = epoint,
   2.820 -  rightarg = ebox,
   2.821 -  procedure = epoint_ebox_overlap_proc,
   2.822 -  commutator = &&,
   2.823 -  restrict = areasel,
   2.824 -  join = areajoinsel
   2.825 -);
   2.826 -
   2.827 -CREATE FUNCTION epoint_ebox_overlap_commutator(ebox, epoint)
   2.828 -  RETURNS boolean
   2.829 -  LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1';
   2.830 -
   2.831 -CREATE OPERATOR && (
   2.832 -  leftarg = ebox,
   2.833 -  rightarg = epoint,
   2.834 -  procedure = epoint_ebox_overlap_commutator,
   2.835 -  commutator = &&,
   2.836 -  restrict = areasel,
   2.837 -  join = areajoinsel
   2.838 -);
   2.839 -
   2.840 -CREATE OPERATOR && (
   2.841 -  leftarg = epoint,
   2.842 -  rightarg = ecircle,
   2.843 -  procedure = epoint_ecircle_overlap_proc,
   2.844 -  commutator = &&,
   2.845 -  restrict = areasel,
   2.846 -  join = areajoinsel
   2.847 -);
   2.848 -
   2.849 -CREATE FUNCTION epoint_ecircle_overlap_commutator(ecircle, epoint)
   2.850 -  RETURNS boolean
   2.851 -  LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1';
   2.852 -
   2.853 -CREATE OPERATOR && (
   2.854 -  leftarg = ecircle,
   2.855 -  rightarg = epoint,
   2.856 -  procedure = epoint_ecircle_overlap_commutator,
   2.857 -  commutator = &&,
   2.858 -  restrict = areasel,
   2.859 -  join = areajoinsel
   2.860 -);
   2.861 -
   2.862 -CREATE OPERATOR && (
   2.863 -  leftarg = epoint,
   2.864 -  rightarg = ecluster,
   2.865 -  procedure = epoint_ecluster_overlap_proc,
   2.866 -  commutator = &&,
   2.867 -  restrict = areasel,
   2.868 -  join = areajoinsel
   2.869 -);
   2.870 -
   2.871 -CREATE FUNCTION epoint_ecluster_overlap_commutator(ecluster, epoint)
   2.872 -  RETURNS boolean
   2.873 -  LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1';
   2.874 -
   2.875 -CREATE OPERATOR && (
   2.876 -  leftarg = ecluster,
   2.877 -  rightarg = epoint,
   2.878 -  procedure = epoint_ecluster_overlap_commutator,
   2.879 -  commutator = &&,
   2.880 -  restrict = areasel,
   2.881 -  join = areajoinsel
   2.882 -);
   2.883 -
   2.884 -CREATE OPERATOR && (
   2.885 -  leftarg = ebox,
   2.886 -  rightarg = ebox,
   2.887 -  procedure = ebox_overlap_proc,
   2.888 -  commutator = &&,
   2.889 -  restrict = areasel,
   2.890 -  join = areajoinsel
   2.891 -);
   2.892 -
   2.893 -CREATE OPERATOR && (
   2.894 -  leftarg = ecircle,
   2.895 -  rightarg = ecircle,
   2.896 -  procedure = ecircle_overlap_proc,
   2.897 -  commutator = &&,
   2.898 -  restrict = areasel,
   2.899 -  join = areajoinsel
   2.900 -);
   2.901 -
   2.902 -CREATE OPERATOR && (
   2.903 -  leftarg = ecircle,
   2.904 -  rightarg = ecluster,
   2.905 -  procedure = ecircle_ecluster_overlap_proc,
   2.906 -  commutator = &&,
   2.907 -  restrict = areasel,
   2.908 -  join = areajoinsel
   2.909 -);
   2.910 -
   2.911 -CREATE FUNCTION ecircle_ecluster_overlap_commutator(ecluster, ecircle)
   2.912 -  RETURNS boolean
   2.913 -  LANGUAGE sql IMMUTABLE AS 'SELECT $2 && $1';
   2.914 -
   2.915 -CREATE OPERATOR && (
   2.916 -  leftarg = ecluster,
   2.917 -  rightarg = ecircle,
   2.918 -  procedure = ecircle_ecluster_overlap_commutator,
   2.919 -  commutator = &&,
   2.920 -  restrict = areasel,
   2.921 -  join = areajoinsel
   2.922 -);
   2.923 -
   2.924 -CREATE OPERATOR && (
   2.925 -  leftarg = ecluster,
   2.926 -  rightarg = ecluster,
   2.927 -  procedure = ecluster_overlap_proc,
   2.928 -  commutator = &&,
   2.929 -  restrict = areasel,
   2.930 -  join = areajoinsel
   2.931 -);
   2.932 -
   2.933 -CREATE FUNCTION ebox_ecircle_overlap_castwrap(ebox, ecircle)
   2.934 -  RETURNS boolean
   2.935 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster && $2';
   2.936 -
   2.937 -CREATE OPERATOR && (
   2.938 -  leftarg = ebox,
   2.939 -  rightarg = ecircle,
   2.940 -  procedure = ebox_ecircle_overlap_castwrap,
   2.941 -  commutator = &&,
   2.942 -  restrict = areasel,
   2.943 -  join = areajoinsel
   2.944 -);
   2.945 -
   2.946 -CREATE FUNCTION ebox_ecircle_overlap_castwrap(ecircle, ebox)
   2.947 -  RETURNS boolean
   2.948 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1 && $2::ecluster';
   2.949 -
   2.950 -CREATE OPERATOR && (
   2.951 -  leftarg = ecircle,
   2.952 -  rightarg = ebox,
   2.953 -  procedure = ebox_ecircle_overlap_castwrap,
   2.954 -  commutator = &&,
   2.955 -  restrict = areasel,
   2.956 -  join = areajoinsel
   2.957 -);
   2.958 -
   2.959 -CREATE FUNCTION ebox_ecluster_overlap_castwrap(ebox, ecluster)
   2.960 -  RETURNS boolean
   2.961 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster && $2';
   2.962 -
   2.963 -CREATE OPERATOR && (
   2.964 -  leftarg = ebox,
   2.965 -  rightarg = ecluster,
   2.966 -  procedure = ebox_ecluster_overlap_castwrap,
   2.967 -  commutator = &&,
   2.968 -  restrict = areasel,
   2.969 -  join = areajoinsel
   2.970 -);
   2.971 -
   2.972 -CREATE FUNCTION ebox_ecluster_overlap_castwrap(ecluster, ebox)
   2.973 -  RETURNS boolean
   2.974 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1 && $2::ecluster';
   2.975 -
   2.976 -CREATE OPERATOR && (
   2.977 -  leftarg = ecluster,
   2.978 -  rightarg = ebox,
   2.979 -  procedure = ebox_ecluster_overlap_castwrap,
   2.980 -  commutator = &&,
   2.981 -  restrict = areasel,
   2.982 -  join = areajoinsel
   2.983 -);
   2.984 -
   2.985 -CREATE OPERATOR &&+ (
   2.986 -  leftarg = epoint,
   2.987 -  rightarg = ecluster,
   2.988 -  procedure = epoint_ecluster_may_overlap_proc,
   2.989 -  commutator = &&+,
   2.990 -  restrict = areasel,
   2.991 -  join = areajoinsel
   2.992 -);
   2.993 -
   2.994 -CREATE FUNCTION epoint_ecluster_may_overlap_commutator(ecluster, epoint)
   2.995 -  RETURNS boolean
   2.996 -  LANGUAGE sql IMMUTABLE AS 'SELECT $2 &&+ $1';
   2.997 -
   2.998 -CREATE OPERATOR &&+ (
   2.999 -  leftarg = ecluster,
  2.1000 -  rightarg = epoint,
  2.1001 -  procedure = epoint_ecluster_may_overlap_commutator,
  2.1002 -  commutator = &&+,
  2.1003 -  restrict = areasel,
  2.1004 -  join = areajoinsel
  2.1005 -);
  2.1006 -
  2.1007 -CREATE OPERATOR &&+ (
  2.1008 -  leftarg = ebox,
  2.1009 -  rightarg = ecircle,
  2.1010 -  procedure = ebox_ecircle_may_overlap_proc,
  2.1011 -  commutator = &&+,
  2.1012 -  restrict = areasel,
  2.1013 -  join = areajoinsel
  2.1014 -);
  2.1015 -
  2.1016 -CREATE FUNCTION ebox_ecircle_may_overlap_commutator(ecircle, ebox)
  2.1017 -  RETURNS boolean
  2.1018 -  LANGUAGE sql IMMUTABLE AS 'SELECT $2 &&+ $1';
  2.1019 -
  2.1020 -CREATE OPERATOR &&+ (
  2.1021 -  leftarg = ecircle,
  2.1022 -  rightarg = ebox,
  2.1023 -  procedure = ebox_ecircle_may_overlap_commutator,
  2.1024 -  commutator = &&+,
  2.1025 -  restrict = areasel,
  2.1026 -  join = areajoinsel
  2.1027 -);
  2.1028 -
  2.1029 -CREATE OPERATOR &&+ (
  2.1030 -  leftarg = ebox,
  2.1031 -  rightarg = ecluster,
  2.1032 -  procedure = ebox_ecluster_may_overlap_proc,
  2.1033 -  commutator = &&+,
  2.1034 -  restrict = areasel,
  2.1035 -  join = areajoinsel
  2.1036 -);
  2.1037 -
  2.1038 -CREATE FUNCTION ebox_ecluster_may_overlap_commutator(ecluster, ebox)
  2.1039 -  RETURNS boolean
  2.1040 -  LANGUAGE sql IMMUTABLE AS 'SELECT $2 &&+ $1';
  2.1041 -
  2.1042 -CREATE OPERATOR &&+ (
  2.1043 -  leftarg = ecluster,
  2.1044 -  rightarg = ebox,
  2.1045 -  procedure = ebox_ecluster_may_overlap_commutator,
  2.1046 -  commutator = &&+,
  2.1047 -  restrict = areasel,
  2.1048 -  join = areajoinsel
  2.1049 -);
  2.1050 -
  2.1051 -CREATE OPERATOR &&+ (
  2.1052 -  leftarg = ecircle,
  2.1053 -  rightarg = ecluster,
  2.1054 -  procedure = ecircle_ecluster_may_overlap_proc,
  2.1055 -  commutator = &&+,
  2.1056 -  restrict = areasel,
  2.1057 -  join = areajoinsel
  2.1058 -);
  2.1059 -
  2.1060 -CREATE FUNCTION ecircle_ecluster_may_overlap_commutator(ecluster, ecircle)
  2.1061 -  RETURNS boolean
  2.1062 -  LANGUAGE sql IMMUTABLE AS 'SELECT $2 &&+ $1';
  2.1063 -
  2.1064 -CREATE OPERATOR &&+ (
  2.1065 -  leftarg = ecluster,
  2.1066 -  rightarg = ecircle,
  2.1067 -  procedure = ecircle_ecluster_may_overlap_commutator,
  2.1068 -  commutator = &&+,
  2.1069 -  restrict = areasel,
  2.1070 -  join = areajoinsel
  2.1071 -);
  2.1072 -
  2.1073 -CREATE OPERATOR &&+ (
  2.1074 -  leftarg = ecluster,
  2.1075 -  rightarg = ecluster,
  2.1076 -  procedure = ecluster_may_overlap_proc,
  2.1077 -  commutator = &&+,
  2.1078 -  restrict = areasel,
  2.1079 -  join = areajoinsel
  2.1080 -);
  2.1081 -
  2.1082 -CREATE OPERATOR @> (
  2.1083 -  leftarg = ebox,
  2.1084 -  rightarg = epoint,
  2.1085 -  procedure = epoint_ebox_overlap_commutator,
  2.1086 -  commutator = <@,
  2.1087 -  restrict = areasel,
  2.1088 -  join = areajoinsel
  2.1089 -);
  2.1090 -
  2.1091 -CREATE OPERATOR <@ (
  2.1092 -  leftarg = epoint,
  2.1093 -  rightarg = ebox,
  2.1094 -  procedure = epoint_ebox_overlap_proc,
  2.1095 -  commutator = @>,
  2.1096 -  restrict = areasel,
  2.1097 -  join = areajoinsel
  2.1098 -);
  2.1099 -
  2.1100 -CREATE OPERATOR @> (
  2.1101 -  leftarg = ecluster,
  2.1102 -  rightarg = epoint,
  2.1103 -  procedure = epoint_ecluster_overlap_commutator,
  2.1104 -  commutator = <@,
  2.1105 -  restrict = areasel,
  2.1106 -  join = areajoinsel
  2.1107 -);
  2.1108 -
  2.1109 -CREATE OPERATOR <@ (
  2.1110 -  leftarg = epoint,
  2.1111 -  rightarg = ecluster,
  2.1112 -  procedure = epoint_ecluster_overlap_proc,
  2.1113 -  commutator = <@,
  2.1114 -  restrict = areasel,
  2.1115 -  join = areajoinsel
  2.1116 -);
  2.1117 -
  2.1118 -CREATE OPERATOR @> (
  2.1119 -  leftarg = ecluster,
  2.1120 -  rightarg = ecluster,
  2.1121 -  procedure = ecluster_contains_proc,
  2.1122 -  commutator = <@,
  2.1123 -  restrict = areasel,
  2.1124 -  join = areajoinsel
  2.1125 -);
  2.1126 -
  2.1127 -CREATE FUNCTION ecluster_contains_commutator(ecluster, ecluster)
  2.1128 -  RETURNS boolean
  2.1129 -  LANGUAGE sql IMMUTABLE AS 'SELECT $2 @> $1';
  2.1130 -
  2.1131 -CREATE OPERATOR <@ (
  2.1132 -  leftarg = ecluster,
  2.1133 -  rightarg = ecluster,
  2.1134 -  procedure = ecluster_contains_commutator,
  2.1135 -  commutator = @>,
  2.1136 -  restrict = areasel,
  2.1137 -  join = areajoinsel
  2.1138 -);
  2.1139 -
  2.1140 -CREATE FUNCTION ebox_contains_castwrap(ebox, ebox)
  2.1141 -  RETURNS boolean
  2.1142 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster @> $2::ecluster';
  2.1143 -
  2.1144 -CREATE OPERATOR @> (
  2.1145 -  leftarg = ebox,
  2.1146 -  rightarg = ebox,
  2.1147 -  procedure = ebox_contains_castwrap,
  2.1148 -  commutator = <@,
  2.1149 -  restrict = areasel,
  2.1150 -  join = areajoinsel
  2.1151 -);
  2.1152 -
  2.1153 -CREATE FUNCTION ebox_contains_swapped_castwrap(ebox, ebox)
  2.1154 -  RETURNS boolean
  2.1155 -  LANGUAGE sql IMMUTABLE AS 'SELECT $2::ecluster @> $1::ecluster';
  2.1156 -
  2.1157 -CREATE OPERATOR <@ (
  2.1158 -  leftarg = ebox,
  2.1159 -  rightarg = ebox,
  2.1160 -  procedure = ebox_contains_swapped_castwrap,
  2.1161 -  commutator = @>,
  2.1162 -  restrict = areasel,
  2.1163 -  join = areajoinsel
  2.1164 -);
  2.1165 -
  2.1166 -CREATE FUNCTION ebox_ecluster_contains_castwrap(ebox, ecluster)
  2.1167 -  RETURNS boolean
  2.1168 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster @> $2';
  2.1169 -
  2.1170 -CREATE OPERATOR @> (
  2.1171 -  leftarg = ebox,
  2.1172 -  rightarg = ecluster,
  2.1173 -  procedure = ebox_ecluster_contains_castwrap,
  2.1174 -  commutator = <@,
  2.1175 -  restrict = areasel,
  2.1176 -  join = areajoinsel
  2.1177 -);
  2.1178 -
  2.1179 -CREATE FUNCTION ebox_ecluster_contains_castwrap(ecluster, ebox)
  2.1180 -  RETURNS boolean
  2.1181 -  LANGUAGE sql IMMUTABLE AS 'SELECT $2::ecluster @> $1';
  2.1182 -
  2.1183 -CREATE OPERATOR <@ (
  2.1184 -  leftarg = ecluster,
  2.1185 -  rightarg = ebox,
  2.1186 -  procedure = ebox_ecluster_contains_castwrap,
  2.1187 -  commutator = @>,
  2.1188 -  restrict = areasel,
  2.1189 -  join = areajoinsel
  2.1190 -);
  2.1191 -
  2.1192 -CREATE FUNCTION ecluster_ebox_contains_castwrap(ecluster, ebox)
  2.1193 -  RETURNS boolean
  2.1194 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1 @> $2::ecluster';
  2.1195 -
  2.1196 -CREATE OPERATOR @> (
  2.1197 -  leftarg = ecluster,
  2.1198 -  rightarg = ebox,
  2.1199 -  procedure = ecluster_ebox_contains_castwrap,
  2.1200 -  commutator = <@,
  2.1201 -  restrict = areasel,
  2.1202 -  join = areajoinsel
  2.1203 -);
  2.1204 -
  2.1205 -CREATE FUNCTION ecluster_ebox_contains_castwrap(ebox, ecluster)
  2.1206 -  RETURNS boolean
  2.1207 -  LANGUAGE sql IMMUTABLE AS 'SELECT $2 @> $1::ecluster';
  2.1208 -
  2.1209 -CREATE OPERATOR <@ (
  2.1210 -  leftarg = ebox,
  2.1211 -  rightarg = ecluster,
  2.1212 -  procedure = ecluster_ebox_contains_castwrap,
  2.1213 -  commutator = @>,
  2.1214 -  restrict = areasel,
  2.1215 -  join = areajoinsel
  2.1216 -);
  2.1217 -
  2.1218 -CREATE OPERATOR <-> (
  2.1219 -  leftarg = epoint,
  2.1220 -  rightarg = epoint,
  2.1221 -  procedure = epoint_distance_proc,
  2.1222 -  commutator = <->
  2.1223 -);
  2.1224 -
  2.1225 -CREATE OPERATOR <-> (
  2.1226 -  leftarg = epoint,
  2.1227 -  rightarg = ecircle,
  2.1228 -  procedure = epoint_ecircle_distance_proc,
  2.1229 -  commutator = <->
  2.1230 -);
  2.1231 -
  2.1232 -CREATE FUNCTION epoint_ecircle_distance_commutator(ecircle, epoint)
  2.1233 -  RETURNS float8
  2.1234 -  LANGUAGE sql IMMUTABLE AS 'SELECT $2 <-> $1';
  2.1235 -
  2.1236 -CREATE OPERATOR <-> (
  2.1237 -  leftarg = ecircle,
  2.1238 -  rightarg = epoint,
  2.1239 -  procedure = epoint_ecircle_distance_commutator,
  2.1240 -  commutator = <->
  2.1241 -);
  2.1242 -
  2.1243 -CREATE OPERATOR <-> (
  2.1244 -  leftarg = epoint,
  2.1245 -  rightarg = ecluster,
  2.1246 -  procedure = epoint_ecluster_distance_proc,
  2.1247 -  commutator = <->
  2.1248 -);
  2.1249 -
  2.1250 -CREATE FUNCTION epoint_ecluster_distance_commutator(ecluster, epoint)
  2.1251 -  RETURNS float8
  2.1252 -  LANGUAGE sql IMMUTABLE AS 'SELECT $2 <-> $1';
  2.1253 -
  2.1254 -CREATE OPERATOR <-> (
  2.1255 -  leftarg = ecluster,
  2.1256 -  rightarg = epoint,
  2.1257 -  procedure = epoint_ecluster_distance_commutator,
  2.1258 -  commutator = <->
  2.1259 -);
  2.1260 -
  2.1261 -CREATE OPERATOR <-> (
  2.1262 -  leftarg = ecircle,
  2.1263 -  rightarg = ecircle,
  2.1264 -  procedure = ecircle_distance_proc,
  2.1265 -  commutator = <->
  2.1266 -);
  2.1267 -
  2.1268 -CREATE OPERATOR <-> (
  2.1269 -  leftarg = ecircle,
  2.1270 -  rightarg = ecluster,
  2.1271 -  procedure = ecircle_ecluster_distance_proc,
  2.1272 -  commutator = <->
  2.1273 -);
  2.1274 -
  2.1275 -CREATE FUNCTION ecircle_ecluster_distance_commutator(ecluster, ecircle)
  2.1276 -  RETURNS float8
  2.1277 -  LANGUAGE sql IMMUTABLE AS 'SELECT $2 <-> $1';
  2.1278 -
  2.1279 -CREATE OPERATOR <-> (
  2.1280 -  leftarg = ecluster,
  2.1281 -  rightarg = ecircle,
  2.1282 -  procedure = ecircle_ecluster_distance_commutator,
  2.1283 -  commutator = <->
  2.1284 -);
  2.1285 -
  2.1286 -CREATE OPERATOR <-> (
  2.1287 -  leftarg = ecluster,
  2.1288 -  rightarg = ecluster,
  2.1289 -  procedure = ecluster_distance_proc,
  2.1290 -  commutator = <->
  2.1291 -);
  2.1292 -
  2.1293 -CREATE FUNCTION epoint_ebox_distance_castwrap(epoint, ebox)
  2.1294 -  RETURNS float8
  2.1295 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2::ecluster';
  2.1296 -
  2.1297 -CREATE OPERATOR <-> (
  2.1298 -  leftarg = epoint,
  2.1299 -  rightarg = ebox,
  2.1300 -  procedure = epoint_ebox_distance_castwrap,
  2.1301 -  commutator = <->
  2.1302 -);
  2.1303 -
  2.1304 -CREATE FUNCTION epoint_ebox_distance_castwrap(ebox, epoint)
  2.1305 -  RETURNS float8
  2.1306 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster <-> $2';
  2.1307 -
  2.1308 -CREATE OPERATOR <-> (
  2.1309 -  leftarg = ebox,
  2.1310 -  rightarg = epoint,
  2.1311 -  procedure = epoint_ebox_distance_castwrap,
  2.1312 -  commutator = <->
  2.1313 -);
  2.1314 -
  2.1315 -CREATE FUNCTION ebox_distance_castwrap(ebox, ebox)
  2.1316 -  RETURNS float8
  2.1317 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster <-> $2::ecluster';
  2.1318 -
  2.1319 -CREATE OPERATOR <-> (
  2.1320 -  leftarg = ebox,
  2.1321 -  rightarg = ebox,
  2.1322 -  procedure = ebox_distance_castwrap,
  2.1323 -  commutator = <->
  2.1324 -);
  2.1325 -
  2.1326 -CREATE FUNCTION ebox_ecircle_distance_castwrap(ebox, ecircle)
  2.1327 -  RETURNS float8
  2.1328 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster <-> $2';
  2.1329 -
  2.1330 -CREATE OPERATOR <-> (
  2.1331 -  leftarg = ebox,
  2.1332 -  rightarg = ecircle,
  2.1333 -  procedure = ebox_ecircle_distance_castwrap,
  2.1334 -  commutator = <->
  2.1335 -);
  2.1336 -
  2.1337 -CREATE FUNCTION ebox_ecircle_distance_castwrap(ecircle, ebox)
  2.1338 -  RETURNS float8
  2.1339 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2::ecluster';
  2.1340 -
  2.1341 -CREATE OPERATOR <-> (
  2.1342 -  leftarg = ecircle,
  2.1343 -  rightarg = ebox,
  2.1344 -  procedure = ebox_ecircle_distance_castwrap,
  2.1345 -  commutator = <->
  2.1346 -);
  2.1347 -
  2.1348 -CREATE FUNCTION ebox_ecluster_distance_castwrap(ebox, ecluster)
  2.1349 -  RETURNS float8
  2.1350 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1::ecluster <-> $2';
  2.1351 -
  2.1352 -CREATE OPERATOR <-> (
  2.1353 -  leftarg = ebox,
  2.1354 -  rightarg = ecluster,
  2.1355 -  procedure = ebox_ecluster_distance_castwrap,
  2.1356 -  commutator = <->
  2.1357 -);
  2.1358 -
  2.1359 -CREATE FUNCTION ebox_ecluster_distance_castwrap(ecluster, ebox)
  2.1360 -  RETURNS float8
  2.1361 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2::ecluster';
  2.1362 -
  2.1363 -CREATE OPERATOR <-> (
  2.1364 -  leftarg = ecluster,
  2.1365 -  rightarg = ebox,
  2.1366 -  procedure = ebox_ecluster_distance_castwrap,
  2.1367 -  commutator = <->
  2.1368 -);
  2.1369 -
  2.1370 -
  2.1371 ----------------------
  2.1372 --- other functions --
  2.1373 ----------------------
  2.1374 -
  2.1375 -CREATE FUNCTION monte_carlo_area(ecluster, int4 = 10000)
  2.1376 -  RETURNS float8
  2.1377 -  LANGUAGE C STRICT
  2.1378 -  AS '$libdir/latlon-v0008', 'pgl_ecluster_monte_carlo_area';
  2.1379 -
  2.1380 -CREATE FUNCTION fair_distance(ecluster, epoint, int4 = 10000)
  2.1381 -  RETURNS float8
  2.1382 -  LANGUAGE C STRICT
  2.1383 -  AS '$libdir/latlon-v0008', 'pgl_ecluster_epoint_fair_distance';
  2.1384 -
  2.1385 -
  2.1386 -----------------
  2.1387 --- GiST index --
  2.1388 -----------------
  2.1389 -
  2.1390 -CREATE FUNCTION pgl_gist_consistent(internal, internal, smallint, oid, internal)
  2.1391 -  RETURNS boolean
  2.1392 -  LANGUAGE C STRICT
  2.1393 -  AS '$libdir/latlon-v0008', 'pgl_gist_consistent';
  2.1394 -
  2.1395 -CREATE FUNCTION pgl_gist_union(internal, internal)
  2.1396 -  RETURNS internal
  2.1397 -  LANGUAGE C STRICT
  2.1398 -  AS '$libdir/latlon-v0008', 'pgl_gist_union';
  2.1399 -
  2.1400 -CREATE FUNCTION pgl_gist_compress_epoint(internal)
  2.1401 -  RETURNS internal
  2.1402 -  LANGUAGE C STRICT
  2.1403 -  AS '$libdir/latlon-v0008', 'pgl_gist_compress_epoint';
  2.1404 -
  2.1405 -CREATE FUNCTION pgl_gist_compress_ecircle(internal)
  2.1406 -  RETURNS internal
  2.1407 -  LANGUAGE C STRICT
  2.1408 -  AS '$libdir/latlon-v0008', 'pgl_gist_compress_ecircle';
  2.1409 -
  2.1410 -CREATE FUNCTION pgl_gist_compress_ecluster(internal)
  2.1411 -  RETURNS internal
  2.1412 -  LANGUAGE C STRICT
  2.1413 -  AS '$libdir/latlon-v0008', 'pgl_gist_compress_ecluster';
  2.1414 -
  2.1415 -CREATE FUNCTION pgl_gist_decompress(internal)
  2.1416 -  RETURNS internal
  2.1417 -  LANGUAGE C STRICT
  2.1418 -  AS '$libdir/latlon-v0008', 'pgl_gist_decompress';
  2.1419 -
  2.1420 -CREATE FUNCTION pgl_gist_penalty(internal, internal, internal)
  2.1421 -  RETURNS internal
  2.1422 -  LANGUAGE C STRICT
  2.1423 -  AS '$libdir/latlon-v0008', 'pgl_gist_penalty';
  2.1424 -
  2.1425 -CREATE FUNCTION pgl_gist_picksplit(internal, internal)
  2.1426 -  RETURNS internal
  2.1427 -  LANGUAGE C STRICT
  2.1428 -  AS '$libdir/latlon-v0008', 'pgl_gist_picksplit';
  2.1429 -
  2.1430 -CREATE FUNCTION pgl_gist_same(internal, internal, internal)
  2.1431 -  RETURNS internal
  2.1432 -  LANGUAGE C STRICT
  2.1433 -  AS '$libdir/latlon-v0008', 'pgl_gist_same';
  2.1434 -
  2.1435 -CREATE FUNCTION pgl_gist_distance(internal, internal, smallint, oid)
  2.1436 -  RETURNS internal
  2.1437 -  LANGUAGE C STRICT
  2.1438 -  AS '$libdir/latlon-v0008', 'pgl_gist_distance';
  2.1439 -
  2.1440 -CREATE OPERATOR CLASS epoint_ops
  2.1441 -  DEFAULT FOR TYPE epoint USING gist AS
  2.1442 -  OPERATOR  11 = ,
  2.1443 -  OPERATOR  22 &&  (epoint, ebox),
  2.1444 -  OPERATOR 222 <@  (epoint, ebox),
  2.1445 -  OPERATOR  23 &&  (epoint, ecircle),
  2.1446 -  OPERATOR  24 &&  (epoint, ecluster),
  2.1447 -  OPERATOR 124 &&+ (epoint, ecluster),
  2.1448 -  OPERATOR 224 <@  (epoint, ecluster),
  2.1449 -  OPERATOR  31 <-> (epoint, epoint) FOR ORDER BY float_ops,
  2.1450 -  OPERATOR  32 <-> (epoint, ebox) FOR ORDER BY float_ops,
  2.1451 -  OPERATOR  33 <-> (epoint, ecircle) FOR ORDER BY float_ops,
  2.1452 -  OPERATOR  34 <-> (epoint, ecluster) FOR ORDER BY float_ops,
  2.1453 -  FUNCTION 1 pgl_gist_consistent(internal, internal, smallint, oid, internal),
  2.1454 -  FUNCTION 2 pgl_gist_union(internal, internal),
  2.1455 -  FUNCTION 3 pgl_gist_compress_epoint(internal),
  2.1456 -  FUNCTION 4 pgl_gist_decompress(internal),
  2.1457 -  FUNCTION 5 pgl_gist_penalty(internal, internal, internal),
  2.1458 -  FUNCTION 6 pgl_gist_picksplit(internal, internal),
  2.1459 -  FUNCTION 7 pgl_gist_same(internal, internal, internal),
  2.1460 -  FUNCTION 8 pgl_gist_distance(internal, internal, smallint, oid),
  2.1461 -  STORAGE ekey_point;
  2.1462 -
  2.1463 -CREATE OPERATOR CLASS ecircle_ops
  2.1464 -  DEFAULT FOR TYPE ecircle USING gist AS
  2.1465 -  OPERATOR  13 = ,
  2.1466 -  OPERATOR  21 &&  (ecircle, epoint),
  2.1467 -  OPERATOR  22 &&  (ecircle, ebox),
  2.1468 -  OPERATOR 122 &&+ (ecircle, ebox),
  2.1469 -  OPERATOR  23 &&  (ecircle, ecircle),
  2.1470 -  OPERATOR  24 &&  (ecircle, ecluster),
  2.1471 -  OPERATOR 124 &&+ (ecircle, ecluster),
  2.1472 -  OPERATOR  31 <-> (ecircle, epoint) FOR ORDER BY float_ops,
  2.1473 -  OPERATOR  32 <-> (ecircle, ebox) FOR ORDER BY float_ops,
  2.1474 -  OPERATOR  33 <-> (ecircle, ecircle) FOR ORDER BY float_ops,
  2.1475 -  OPERATOR  34 <-> (ecircle, ecluster) FOR ORDER BY float_ops,
  2.1476 -  FUNCTION 1 pgl_gist_consistent(internal, internal, smallint, oid, internal),
  2.1477 -  FUNCTION 2 pgl_gist_union(internal, internal),
  2.1478 -  FUNCTION 3 pgl_gist_compress_ecircle(internal),
  2.1479 -  FUNCTION 4 pgl_gist_decompress(internal),
  2.1480 -  FUNCTION 5 pgl_gist_penalty(internal, internal, internal),
  2.1481 -  FUNCTION 6 pgl_gist_picksplit(internal, internal),
  2.1482 -  FUNCTION 7 pgl_gist_same(internal, internal, internal),
  2.1483 -  FUNCTION 8 pgl_gist_distance(internal, internal, smallint, oid),
  2.1484 -  STORAGE ekey_area;
  2.1485 -
  2.1486 -CREATE OPERATOR CLASS ecluster_ops
  2.1487 -  DEFAULT FOR TYPE ecluster USING gist AS
  2.1488 -  OPERATOR  21 &&  (ecluster, epoint),
  2.1489 -  OPERATOR 121 &&+ (ecluster, epoint),
  2.1490 -  OPERATOR 221 @>  (ecluster, epoint),
  2.1491 -  OPERATOR  22 &&  (ecluster, ebox),
  2.1492 -  OPERATOR 122 &&+ (ecluster, ebox),
  2.1493 -  OPERATOR 222 @>  (ecluster, ebox),
  2.1494 -  OPERATOR 322 <@  (ecluster, ebox),
  2.1495 -  OPERATOR  23 &&  (ecluster, ecircle),
  2.1496 -  OPERATOR 123 &&+ (ecluster, ecircle),
  2.1497 -  OPERATOR  24 &&  (ecluster, ecluster),
  2.1498 -  OPERATOR 124 &&+ (ecluster, ecluster),
  2.1499 -  OPERATOR 224 @>  (ecluster, ecluster),
  2.1500 -  OPERATOR 324 <@  (ecluster, ecluster),
  2.1501 -  OPERATOR  31 <-> (ecluster, epoint) FOR ORDER BY float_ops,
  2.1502 -  OPERATOR  32 <-> (ecluster, ebox) FOR ORDER BY float_ops,
  2.1503 -  OPERATOR  33 <-> (ecluster, ecircle) FOR ORDER BY float_ops,
  2.1504 -  OPERATOR  34 <-> (ecluster, ecluster) FOR ORDER BY float_ops,
  2.1505 -  FUNCTION 1 pgl_gist_consistent(internal, internal, smallint, oid, internal),
  2.1506 -  FUNCTION 2 pgl_gist_union(internal, internal),
  2.1507 -  FUNCTION 3 pgl_gist_compress_ecluster(internal),
  2.1508 -  FUNCTION 4 pgl_gist_decompress(internal),
  2.1509 -  FUNCTION 5 pgl_gist_penalty(internal, internal, internal),
  2.1510 -  FUNCTION 6 pgl_gist_picksplit(internal, internal),
  2.1511 -  FUNCTION 7 pgl_gist_same(internal, internal, internal),
  2.1512 -  FUNCTION 8 pgl_gist_distance(internal, internal, smallint, oid),
  2.1513 -  STORAGE ekey_area;
  2.1514 -
  2.1515 -
  2.1516 ----------------------
  2.1517 --- alias functions --
  2.1518 ----------------------
  2.1519 -
  2.1520 -CREATE FUNCTION distance(epoint, epoint)
  2.1521 -  RETURNS float8
  2.1522 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2';
  2.1523 -
  2.1524 -CREATE FUNCTION distance(ecluster, epoint)
  2.1525 -  RETURNS float8
  2.1526 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1 <-> $2';
  2.1527 -
  2.1528 -CREATE FUNCTION distance_within(epoint, epoint, float8)
  2.1529 -  RETURNS boolean
  2.1530 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1 && ecircle($2, $3)';
  2.1531 -
  2.1532 -CREATE FUNCTION distance_within(ecluster, epoint, float8)
  2.1533 -  RETURNS boolean
  2.1534 -  LANGUAGE sql IMMUTABLE AS 'SELECT $1 && ecircle($2, $3)';
  2.1535 -
  2.1536 -
  2.1537 ---------------------------------
  2.1538 --- other data storage formats --
  2.1539 ---------------------------------
  2.1540 -
  2.1541 -CREATE FUNCTION coords_to_epoint(float8, float8, text = 'epoint')
  2.1542 -  RETURNS epoint
  2.1543 -  LANGUAGE plpgsql IMMUTABLE STRICT AS $$
  2.1544 -    DECLARE
  2.1545 -      "result" epoint;
  2.1546 -    BEGIN
  2.1547 -      IF $3 = 'epoint_lonlat' THEN
  2.1548 -        -- avoid dynamic command execution for better performance
  2.1549 -        RETURN epoint($2, $1);
  2.1550 -      END IF;
  2.1551 -      IF $3 = 'epoint' OR $3 = 'epoint_latlon' THEN
  2.1552 -        -- avoid dynamic command execution for better performance
  2.1553 -        RETURN epoint($1, $2);
  2.1554 -      END IF;
  2.1555 -      EXECUTE 'SELECT ' || $3 || '($1, $2)' INTO STRICT "result" USING $1, $2;
  2.1556 -      RETURN "result";
  2.1557 -    END;
  2.1558 -  $$;
  2.1559 -
  2.1560 -CREATE FUNCTION GeoJSON_LinearRing_vertices(jsonb, text = 'epoint_lonlat')
  2.1561 -  RETURNS SETOF jsonb
  2.1562 -  LANGUAGE sql IMMUTABLE STRICT AS $$
  2.1563 -    SELECT "result" FROM
  2.1564 -      ( SELECT jsonb_array_length($1) - 1 ) AS "lastindex_row" ("lastindex")
  2.1565 -      CROSS JOIN LATERAL jsonb_array_elements(
  2.1566 -        CASE WHEN
  2.1567 -          coords_to_epoint(
  2.1568 -            ($1->0->>0)::float8,
  2.1569 -            ($1->0->>1)::float8,
  2.1570 -            $2
  2.1571 -          ) = coords_to_epoint(
  2.1572 -            ($1->"lastindex"->>0)::float8,
  2.1573 -            ($1->"lastindex"->>1)::float8,
  2.1574 -            $2
  2.1575 -          )
  2.1576 -        THEN
  2.1577 -          $1 - "lastindex"
  2.1578 -        ELSE
  2.1579 -          $1
  2.1580 -        END
  2.1581 -      ) AS "result_row" ("result")
  2.1582 -  $$;
  2.1583 -
  2.1584 -CREATE FUNCTION GeoJSON_to_epoint(jsonb, text = 'epoint_lonlat')
  2.1585 -  RETURNS epoint
  2.1586 -  LANGUAGE sql IMMUTABLE STRICT AS $$
  2.1587 -    SELECT CASE
  2.1588 -    WHEN $1->>'type' = 'Point' THEN
  2.1589 -      coords_to_epoint(
  2.1590 -        ($1->'coordinates'->>0)::float8,
  2.1591 -        ($1->'coordinates'->>1)::float8,
  2.1592 -        $2
  2.1593 -      )
  2.1594 -    WHEN $1->>'type' = 'Feature' THEN
  2.1595 -      GeoJSON_to_epoint($1->'geometry', $2)
  2.1596 -    ELSE
  2.1597 -      NULL
  2.1598 -    END
  2.1599 -  $$;
  2.1600 -
  2.1601 -CREATE FUNCTION GeoJSON_to_ecluster(jsonb, text = 'epoint_lonlat')
  2.1602 -  RETURNS ecluster
  2.1603 -  LANGUAGE sql IMMUTABLE STRICT AS $$
  2.1604 -    SELECT CASE $1->>'type'
  2.1605 -    WHEN 'Point' THEN
  2.1606 -      coords_to_epoint(
  2.1607 -        ($1->'coordinates'->>0)::float8,
  2.1608 -        ($1->'coordinates'->>1)::float8,
  2.1609 -        $2
  2.1610 -      )::ecluster
  2.1611 -    WHEN 'MultiPoint' THEN
  2.1612 -      ( SELECT ecluster_create_multipoint(array_agg(
  2.1613 -          coords_to_epoint(
  2.1614 -            ("coord"->>0)::float8,
  2.1615 -            ("coord"->>1)::float8,
  2.1616 -            $2
  2.1617 -          )
  2.1618 -        ))
  2.1619 -        FROM jsonb_array_elements($1->'coordinates') AS "coord"
  2.1620 -      )
  2.1621 -    WHEN 'LineString' THEN
  2.1622 -      ( SELECT ecluster_create_path(array_agg(
  2.1623 -          coords_to_epoint(
  2.1624 -            ("coord"->>0)::float8,
  2.1625 -            ("coord"->>1)::float8,
  2.1626 -            $2
  2.1627 -          )
  2.1628 -        ))
  2.1629 -        FROM jsonb_array_elements($1->'coordinates') AS "coord"
  2.1630 -      )
  2.1631 -    WHEN 'MultiLineString' THEN
  2.1632 -      ( SELECT ecluster_concat(array_agg(
  2.1633 -          ( SELECT ecluster_create_path(array_agg(
  2.1634 -              coords_to_epoint(
  2.1635 -                ("coord"->>0)::float8,
  2.1636 -                ("coord"->>1)::float8,
  2.1637 -                $2
  2.1638 -              )
  2.1639 -            ))
  2.1640 -            FROM jsonb_array_elements("coord_array") AS "coord"
  2.1641 -          )
  2.1642 -        ))
  2.1643 -        FROM jsonb_array_elements($1->'coordinates') AS "coord_array"
  2.1644 -      )
  2.1645 -    WHEN 'Polygon' THEN
  2.1646 -      ( SELECT ecluster_concat(array_agg(
  2.1647 -          ( SELECT ecluster_create_polygon(array_agg(
  2.1648 -              coords_to_epoint(
  2.1649 -                ("coord"->>0)::float8,
  2.1650 -                ("coord"->>1)::float8,
  2.1651 -                $2
  2.1652 -              )
  2.1653 -            ))
  2.1654 -            FROM GeoJSON_LinearRing_vertices("coord_array", $2) AS "coord"
  2.1655 -          )
  2.1656 -        ))
  2.1657 -        FROM jsonb_array_elements($1->'coordinates') AS "coord_array"
  2.1658 -      )
  2.1659 -    WHEN 'MultiPolygon' THEN
  2.1660 -      ( SELECT ecluster_concat(array_agg(
  2.1661 -          ( SELECT ecluster_concat(array_agg(
  2.1662 -              ( SELECT ecluster_create_polygon(array_agg(
  2.1663 -                  coords_to_epoint(
  2.1664 -                    ("coord"->>0)::float8,
  2.1665 -                    ("coord"->>1)::float8,
  2.1666 -                    $2
  2.1667 -                  )
  2.1668 -                ))
  2.1669 -                FROM GeoJSON_LinearRing_vertices("coord_array", $2) AS "coord"
  2.1670 -              )
  2.1671 -            ))
  2.1672 -            FROM jsonb_array_elements("coord_array_array") AS "coord_array"
  2.1673 -          )
  2.1674 -        ))
  2.1675 -        FROM jsonb_array_elements($1->'coordinates') AS "coord_array_array"
  2.1676 -      )
  2.1677 -    WHEN 'GeometryCollection' THEN
  2.1678 -      ( SELECT ecluster_concat(array_agg(
  2.1679 -          GeoJSON_to_ecluster("geometry", $2)
  2.1680 -        ))
  2.1681 -        FROM jsonb_array_elements($1->'geometries') AS "geometry"
  2.1682 -      )
  2.1683 -    WHEN 'Feature' THEN
  2.1684 -      GeoJSON_to_ecluster($1->'geometry', $2)
  2.1685 -    WHEN 'FeatureCollection' THEN
  2.1686 -      ( SELECT ecluster_concat(array_agg(
  2.1687 -          GeoJSON_to_ecluster("feature", $2)
  2.1688 -        ))
  2.1689 -        FROM jsonb_array_elements($1->'features') AS "feature"
  2.1690 -      )
  2.1691 -    ELSE
  2.1692 -      NULL
  2.1693 -    END
  2.1694 -  $$;
  2.1695 -

     3.1 --- a/latlon-v0008.c	Mon Oct 31 13:06:31 2016 +0100
     3.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
     3.3 @@ -1,3285 +0,0 @@
     3.4 -
     3.5 -/*-------------*
     3.6 - *  C prelude  *
     3.7 - *-------------*/
     3.8 -
     3.9 -#include "postgres.h"
    3.10 -#include "fmgr.h"
    3.11 -#include "libpq/pqformat.h"
    3.12 -#include "access/gist.h"
    3.13 -#include "access/stratnum.h"
    3.14 -#include "utils/array.h"
    3.15 -#include <math.h>
    3.16 -
    3.17 -#ifdef PG_MODULE_MAGIC
    3.18 -PG_MODULE_MAGIC;
    3.19 -#endif
    3.20 -
    3.21 -#if INT_MAX < 2147483647
    3.22 -#error Expected int type to be at least 32 bit wide
    3.23 -#endif
    3.24 -
    3.25 -
    3.26 -/*---------------------------------*
    3.27 - *  distance calculation on earth  *
    3.28 - *  (using WGS-84 spheroid)        *
    3.29 - *---------------------------------*/
    3.30 -
    3.31 -/*  WGS-84 spheroid with following parameters:
    3.32 -    semi-major axis  a = 6378137
    3.33 -    semi-minor axis  b = a * (1 - 1/298.257223563)
    3.34 -    estimated diameter = 2 * (2*a+b)/3
    3.35 -*/
    3.36 -#define PGL_SPHEROID_A 6378137.0            /* semi major axis */
    3.37 -#define PGL_SPHEROID_F (1.0/298.257223563)  /* flattening */
    3.38 -#define PGL_SPHEROID_B (PGL_SPHEROID_A * (1.0-PGL_SPHEROID_F))
    3.39 -#define PGL_EPS2       ( ( PGL_SPHEROID_A * PGL_SPHEROID_A - \
    3.40 -                           PGL_SPHEROID_B * PGL_SPHEROID_B ) / \
    3.41 -                         ( PGL_SPHEROID_A * PGL_SPHEROID_A ) )
    3.42 -#define PGL_SUBEPS2    (1.0-PGL_EPS2)
    3.43 -#define PGL_RADIUS     ((2.0*PGL_SPHEROID_A + PGL_SPHEROID_B) / 3.0)
    3.44 -#define PGL_DIAMETER   (2.0 * PGL_RADIUS)
    3.45 -#define PGL_SCALE      (PGL_SPHEROID_A / PGL_DIAMETER)  /* semi-major ref. */
    3.46 -#define PGL_MAXDIST    (PGL_RADIUS * M_PI)              /* maximum distance */
    3.47 -#define PGL_FADELIMIT  (PGL_MAXDIST / 3.0)              /* 1/6 circumference */
    3.48 -
    3.49 -/* calculate distance between two points on earth (given in degrees) */
    3.50 -static inline double pgl_distance(
    3.51 -  double lat1, double lon1, double lat2, double lon2
    3.52 -) {
    3.53 -  float8 lat1cos, lat1sin, lat2cos, lat2sin, lon2cos, lon2sin;
    3.54 -  float8 nphi1, nphi2, x1, z1, x2, y2, z2, g, s, t;
    3.55 -  /* normalize delta longitude (lon2 > 0 && lon1 = 0) */
    3.56 -  /* lon1 = 0 (not used anymore) */
    3.57 -  lon2 = fabs(lon2-lon1);
    3.58 -  /* convert to radians (first divide, then multiply) */
    3.59 -  lat1 = (lat1 / 180.0) * M_PI;
    3.60 -  lat2 = (lat2 / 180.0) * M_PI;
    3.61 -  lon2 = (lon2 / 180.0) * M_PI;
    3.62 -  /* make lat2 >= lat1 to ensure reversal-symmetry despite floating point
    3.63 -     operations (lon2 >= lon1 is already ensured in a previous step) */
    3.64 -  if (lat2 < lat1) { float8 swap = lat1; lat1 = lat2; lat2 = swap; }
    3.65 -  /* calculate 3d coordinates on scaled ellipsoid which has an average diameter
    3.66 -     of 1.0 */
    3.67 -  lat1cos = cos(lat1); lat1sin = sin(lat1);
    3.68 -  lat2cos = cos(lat2); lat2sin = sin(lat2);
    3.69 -  lon2cos = cos(lon2); lon2sin = sin(lon2);
    3.70 -  nphi1 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat1sin * lat1sin);
    3.71 -  nphi2 = PGL_SCALE / sqrt(1 - PGL_EPS2 * lat2sin * lat2sin);
    3.72 -  x1 = nphi1 * lat1cos;
    3.73 -  z1 = nphi1 * PGL_SUBEPS2 * lat1sin;
    3.74 -  x2 = nphi2 * lat2cos * lon2cos;
    3.75 -  y2 = nphi2 * lat2cos * lon2sin;
    3.76 -  z2 = nphi2 * PGL_SUBEPS2 * lat2sin;
    3.77 -  /* calculate tunnel distance through scaled (diameter 1.0) ellipsoid */
    3.78 -  g = sqrt((x2-x1)*(x2-x1) + y2*y2 + (z2-z1)*(z2-z1));
    3.79 -  /* convert tunnel distance through scaled ellipsoid to approximated surface
    3.80 -     distance on original ellipsoid */
    3.81 -  if (g > 1.0) g = 1.0;
    3.82 -  s = PGL_DIAMETER * asin(g);
    3.83 -  /* return result only if small enough to be precise (less than 1/3 of
    3.84 -     maximum possible distance) */
    3.85 -  if (s <= PGL_FADELIMIT) return s;
    3.86 -  /* calculate tunnel distance to antipodal point through scaled ellipsoid */
    3.87 -  g = sqrt((x2+x1)*(x2+x1) + y2*y2 + (z2+z1)*(z2+z1));
    3.88 -  /* convert tunnel distance to antipodal point through scaled ellipsoid to
    3.89 -     approximated surface distance to antipodal point on original ellipsoid */
    3.90 -  if (g > 1.0) g = 1.0;
    3.91 -  t = PGL_DIAMETER * asin(g);
    3.92 -  /* surface distance between original points can now be approximated by
    3.93 -     substracting antipodal distance from maximum possible distance;
    3.94 -     return result only if small enough (less than 1/3 of maximum possible
    3.95 -     distance) */
    3.96 -  if (t <= PGL_FADELIMIT) return PGL_MAXDIST-t;
    3.97 -  /* otherwise crossfade direct and antipodal result to ensure monotonicity */
    3.98 -  return (
    3.99 -    (s * (t-PGL_FADELIMIT) + (PGL_MAXDIST-t) * (s-PGL_FADELIMIT)) /
   3.100 -    (s + t - 2*PGL_FADELIMIT)
   3.101 -  );
   3.102 -}
   3.103 -
   3.104 -/* finite distance that can not be reached on earth */
   3.105 -#define PGL_ULTRA_DISTANCE (3 * PGL_MAXDIST)
   3.106 -
   3.107 -
   3.108 -/*--------------------------------*
   3.109 - *  simple geographic data types  *
   3.110 - *--------------------------------*/
   3.111 -
   3.112 -/* point on earth given by latitude and longitude in degrees */
   3.113 -/* (type "epoint" in SQL) */
   3.114 -typedef struct {
   3.115 -  double lat;  /* between  -90 and  90 (both inclusive) */
   3.116 -  double lon;  /* between -180 and 180 (both inclusive) */
   3.117 -} pgl_point;
   3.118 -
   3.119 -/* box delimited by two parallels and two meridians (all in degrees) */
   3.120 -/* (type "ebox" in SQL) */
   3.121 -typedef struct {
   3.122 -  double lat_min;  /* between  -90 and  90 (both inclusive) */
   3.123 -  double lat_max;  /* between  -90 and  90 (both inclusive) */
   3.124 -  double lon_min;  /* between -180 and 180 (both inclusive) */
   3.125 -  double lon_max;  /* between -180 and 180 (both inclusive) */
   3.126 -  /* if lat_min > lat_max, then box is empty */
   3.127 -  /* if lon_min > lon_max, then 180th meridian is crossed */
   3.128 -} pgl_box;
   3.129 -
   3.130 -/* circle on earth surface (for radial searches with fixed radius) */
   3.131 -/* (type "ecircle" in SQL) */
   3.132 -typedef struct {
   3.133 -  pgl_point center;
   3.134 -  double radius; /* positive (including +0 but excluding -0), or -INFINITY */
   3.135 -  /* A negative radius (i.e. -INFINITY) denotes nothing (i.e. no point),
   3.136 -     zero radius (0) denotes a single point,
   3.137 -     a finite radius (0 < radius < INFINITY) denotes a filled circle, and
   3.138 -     a radius of INFINITY is valid and means complete coverage of earth. */
   3.139 -} pgl_circle;
   3.140 -
   3.141 -
   3.142 -/*----------------------------------*
   3.143 - *  geographic "cluster" data type  *
   3.144 - *----------------------------------*/
   3.145 -
   3.146 -/* A cluster is a collection of points, paths, outlines, and polygons. If two
   3.147 -   polygons in a cluster overlap, the area covered by both polygons does not
   3.148 -   belong to the cluster. This way, a cluster can be used to describe complex
   3.149 -   shapes like polygons with holes. Outlines are non-filled polygons. Paths are
   3.150 -   open by default (i.e. the last point in the list is not connected with the
   3.151 -   first point in the list). Note that each outline or polygon in a cluster
   3.152 -   must cover a longitude range of less than 180 degrees to avoid ambiguities.
   3.153 -   Areas which are larger may be split into multiple polygons. */
   3.154 -
   3.155 -/* maximum number of points in a cluster */
   3.156 -/* (limited to avoid integer overflows, e.g. when allocating memory) */
   3.157 -#define PGL_CLUSTER_MAXPOINTS 16777216
   3.158 -
   3.159 -/* types of cluster entries */
   3.160 -#define PGL_ENTRY_POINT   1  /* a point */
   3.161 -#define PGL_ENTRY_PATH    2  /* a path from first point to last point */
   3.162 -#define PGL_ENTRY_OUTLINE 3  /* a non-filled polygon with given vertices */
   3.163 -#define PGL_ENTRY_POLYGON 4  /* a filled polygon with given vertices */
   3.164 -
   3.165 -/* Entries of a cluster are described by two different structs: pgl_newentry
   3.166 -   and pgl_entry. The first is used only during construction of a cluster, the
   3.167 -   second is used in all other cases (e.g. when reading clusters from the
   3.168 -   database, performing operations, etc). */
   3.169 -
   3.170 -/* entry for new geographic cluster during construction of that cluster */
   3.171 -typedef struct {
   3.172 -  int32_t entrytype;
   3.173 -  int32_t npoints;
   3.174 -  pgl_point *points;  /* pointer to an array of points (pgl_point) */
   3.175 -} pgl_newentry;
   3.176 -
   3.177 -/* entry of geographic cluster */
   3.178 -typedef struct {
   3.179 -  int32_t entrytype;  /* type of entry: point, path, outline, polygon */
   3.180 -  int32_t npoints;    /* number of stored points (set to 1 for point entry) */
   3.181 -  int32_t offset;     /* offset of pgl_point array from cluster base address */
   3.182 -  /* use macro PGL_ENTRY_POINTS to obtain a pointer to the array of points */
   3.183 -} pgl_entry;
   3.184 -
   3.185 -/* geographic cluster which is a collection of points, (open) paths, polygons,
   3.186 -   and outlines (non-filled polygons) */
   3.187 -typedef struct {
   3.188 -  char header[VARHDRSZ];  /* PostgreSQL header for variable size data types */
   3.189 -  int32_t nentries;       /* number of stored points */
   3.190 -  pgl_circle bounding;    /* bounding circle */
   3.191 -  /* Note: bounding circle ensures alignment of pgl_cluster for points */
   3.192 -  pgl_entry entries[FLEXIBLE_ARRAY_MEMBER];  /* var-length data */
   3.193 -} pgl_cluster;
   3.194 -
   3.195 -/* macro to determine memory alignment of points */
   3.196 -/* (needed to store pgl_point array after entries in pgl_cluster) */
   3.197 -typedef struct { char dummy; pgl_point aligned; } pgl_point_alignment;
   3.198 -#define PGL_POINT_ALIGNMENT offsetof(pgl_point_alignment, aligned)
   3.199 -
   3.200 -/* macro to extract a pointer to the array of points of a cluster entry */
   3.201 -#define PGL_ENTRY_POINTS(cluster, idx) \
   3.202 -  ((pgl_point *)(((intptr_t)cluster)+(cluster)->entries[idx].offset))
   3.203 -
   3.204 -/* convert pgl_newentry array to pgl_cluster */
   3.205 -/* NOTE: requires pgl_finalize_cluster to be called to finalize result */
   3.206 -static pgl_cluster *pgl_new_cluster(int nentries, pgl_newentry *entries) {
   3.207 -  int i;              /* index of current entry */
   3.208 -  int npoints = 0;    /* number of points in whole cluster */
   3.209 -  int entry_npoints;  /* number of points in current entry */
   3.210 -  int points_offset = PGL_POINT_ALIGNMENT * (
   3.211 -    ( offsetof(pgl_cluster, entries) +
   3.212 -      nentries * sizeof(pgl_entry) +
   3.213 -      PGL_POINT_ALIGNMENT - 1
   3.214 -    ) / PGL_POINT_ALIGNMENT
   3.215 -  );  /* offset of pgl_point array from base address (considering alignment) */
   3.216 -  pgl_cluster *cluster;  /* new cluster to be returned */
   3.217 -  /* determine total number of points */
   3.218 -  for (i=0; i<nentries; i++) npoints += entries[i].npoints;
   3.219 -  /* allocate memory for cluster (including entries and points) */
   3.220 -  cluster = palloc(points_offset + npoints * sizeof(pgl_point));
   3.221 -  /* re-count total number of points to determine offset for each entry */
   3.222 -  npoints = 0;
   3.223 -  /* copy entries and points */
   3.224 -  for (i=0; i<nentries; i++) {
   3.225 -    /* determine number of points in entry */
   3.226 -    entry_npoints = entries[i].npoints;
   3.227 -    /* copy entry */
   3.228 -    cluster->entries[i].entrytype = entries[i].entrytype;
   3.229 -    cluster->entries[i].npoints = entry_npoints;
   3.230 -    /* calculate offset (in bytes) of pgl_point array */
   3.231 -    cluster->entries[i].offset = points_offset + npoints * sizeof(pgl_point);
   3.232 -    /* copy points */
   3.233 -    memcpy(
   3.234 -      PGL_ENTRY_POINTS(cluster, i),
   3.235 -      entries[i].points,
   3.236 -      entry_npoints * sizeof(pgl_point)
   3.237 -    );
   3.238 -    /* update total number of points processed */
   3.239 -    npoints += entry_npoints;
   3.240 -  }
   3.241 -  /* set number of entries in cluster */
   3.242 -  cluster->nentries = nentries;
   3.243 -  /* set PostgreSQL header for variable sized data */
   3.244 -  SET_VARSIZE(cluster, points_offset + npoints * sizeof(pgl_point));
   3.245 -  /* return newly created cluster */
   3.246 -  return cluster;
   3.247 -}
   3.248 -
   3.249 -
   3.250 -/*----------------------------------------*
   3.251 - *  C functions on geographic data types  *
   3.252 - *----------------------------------------*/
   3.253 -
   3.254 -/* round latitude or longitude to 12 digits after decimal point */
   3.255 -static inline double pgl_round(double val) {
   3.256 -  return round(val * 1e12) / 1e12;
   3.257 -}
   3.258 -
   3.259 -/* compare two points */
   3.260 -/* (equality when same point on earth is described, otherwise an arbitrary
   3.261 -   linear order) */
   3.262 -static int pgl_point_cmp(pgl_point *point1, pgl_point *point2) {
   3.263 -  double lon1, lon2;  /* modified longitudes for special cases */
   3.264 -  /* use latitude as first ordering criterion */
   3.265 -  if (point1->lat < point2->lat) return -1;
   3.266 -  if (point1->lat > point2->lat) return 1;
   3.267 -  /* determine modified longitudes (considering special case of poles and
   3.268 -     180th meridian which can be described as W180 or E180) */
   3.269 -  if (point1->lat == -90 || point1->lat == 90) lon1 = 0;
   3.270 -  else if (point1->lon == 180) lon1 = -180;
   3.271 -  else lon1 = point1->lon;
   3.272 -  if (point2->lat == -90 || point2->lat == 90) lon2 = 0;
   3.273 -  else if (point2->lon == 180) lon2 = -180;
   3.274 -  else lon2 = point2->lon;
   3.275 -  /* use (modified) longitude as secondary ordering criterion */
   3.276 -  if (lon1 < lon2) return -1;
   3.277 -  if (lon1 > lon2) return 1;
   3.278 -  /* no difference found, points are equal */
   3.279 -  return 0;
   3.280 -}
   3.281 -
   3.282 -/* compare two boxes */
   3.283 -/* (equality when same box on earth is described, otherwise an arbitrary linear
   3.284 -   order) */
   3.285 -static int pgl_box_cmp(pgl_box *box1, pgl_box *box2) {
   3.286 -  /* two empty boxes are equal, and an empty box is always considered "less
   3.287 -     than" a non-empty box */
   3.288 -  if (box1->lat_min> box1->lat_max && box2->lat_min<=box2->lat_max) return -1;
   3.289 -  if (box1->lat_min> box1->lat_max && box2->lat_min> box2->lat_max) return 0;
   3.290 -  if (box1->lat_min<=box1->lat_max && box2->lat_min> box2->lat_max) return 1;
   3.291 -  /* use southern border as first ordering criterion */
   3.292 -  if (box1->lat_min < box2->lat_min) return -1;
   3.293 -  if (box1->lat_min > box2->lat_min) return 1;
   3.294 -  /* use northern border as second ordering criterion */
   3.295 -  if (box1->lat_max < box2->lat_max) return -1;
   3.296 -  if (box1->lat_max > box2->lat_max) return 1;
   3.297 -  /* use western border as third ordering criterion */
   3.298 -  if (box1->lon_min < box2->lon_min) return -1;
   3.299 -  if (box1->lon_min > box2->lon_min) return 1;
   3.300 -  /* use eastern border as fourth ordering criterion */
   3.301 -  if (box1->lon_max < box2->lon_max) return -1;
   3.302 -  if (box1->lon_max > box2->lon_max) return 1;
   3.303 -  /* no difference found, boxes are equal */
   3.304 -  return 0;
   3.305 -}
   3.306 -
   3.307 -/* compare two circles */
   3.308 -/* (equality when same circle on earth is described, otherwise an arbitrary
   3.309 -   linear order) */
   3.310 -static int pgl_circle_cmp(pgl_circle *circle1, pgl_circle *circle2) {
   3.311 -  /* two circles with same infinite radius (positive or negative infinity) are
   3.312 -     considered equal independently of center point */
   3.313 -  if (
   3.314 -    !isfinite(circle1->radius) && !isfinite(circle2->radius) &&
   3.315 -    circle1->radius == circle2->radius
   3.316 -  ) return 0;
   3.317 -  /* use radius as first ordering criterion */
   3.318 -  if (circle1->radius < circle2->radius) return -1;
   3.319 -  if (circle1->radius > circle2->radius) return 1;
   3.320 -  /* use center point as secondary ordering criterion */
   3.321 -  return pgl_point_cmp(&(circle1->center), &(circle2->center));
   3.322 -}
   3.323 -
   3.324 -/* set box to empty box*/
   3.325 -static void pgl_box_set_empty(pgl_box *box) {
   3.326 -  box->lat_min = INFINITY;
   3.327 -  box->lat_max = -INFINITY;
   3.328 -  box->lon_min = 0;
   3.329 -  box->lon_max = 0;
   3.330 -}
   3.331 -
   3.332 -/* check if point is inside a box */
   3.333 -static bool pgl_point_in_box(pgl_point *point, pgl_box *box) {
   3.334 -  return (
   3.335 -    point->lat >= box->lat_min && point->lat <= box->lat_max && (
   3.336 -      (box->lon_min > box->lon_max) ? (
   3.337 -        /* box crosses 180th meridian */
   3.338 -        point->lon >= box->lon_min || point->lon <= box->lon_max
   3.339 -      ) : (
   3.340 -        /* box does not cross the 180th meridian */
   3.341 -        point->lon >= box->lon_min && point->lon <= box->lon_max
   3.342 -      )
   3.343 -    )
   3.344 -  );
   3.345 -}
   3.346 -
   3.347 -/* check if two boxes overlap */
   3.348 -static bool pgl_boxes_overlap(pgl_box *box1, pgl_box *box2) {
   3.349 -  return (
   3.350 -    box2->lat_max >= box2->lat_min &&  /* ensure box2 is not empty */
   3.351 -    ( box2->lat_min >= box1->lat_min || box2->lat_max >= box1->lat_min ) &&
   3.352 -    ( box2->lat_min <= box1->lat_max || box2->lat_max <= box1->lat_max ) && (
   3.353 -      (
   3.354 -        /* check if one and only one box crosses the 180th meridian */
   3.355 -        ((box1->lon_min > box1->lon_max) ? 1 : 0) ^
   3.356 -        ((box2->lon_min > box2->lon_max) ? 1 : 0)
   3.357 -      ) ? (
   3.358 -        /* exactly one box crosses the 180th meridian */
   3.359 -        box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min ||
   3.360 -        box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max
   3.361 -      ) : (
   3.362 -        /* no box or both boxes cross the 180th meridian */
   3.363 -        (
   3.364 -          (box2->lon_min >= box1->lon_min || box2->lon_max >= box1->lon_min) &&
   3.365 -          (box2->lon_min <= box1->lon_max || box2->lon_max <= box1->lon_max)
   3.366 -        ) ||
   3.367 -        /* handle W180 == E180 */
   3.368 -        ( box1->lon_min == -180 && box2->lon_max == 180 ) ||
   3.369 -        ( box2->lon_min == -180 && box1->lon_max == 180 )
   3.370 -      )
   3.371 -    )
   3.372 -  );
   3.373 -}
   3.374 -
   3.375 -/* check unambiguousness of east/west orientation of cluster entries and set
   3.376 -   bounding circle of cluster */
   3.377 -static bool pgl_finalize_cluster(pgl_cluster *cluster) {
   3.378 -  int i, j;                 /* i: index of entry, j: index of point in entry */
   3.379 -  int npoints;              /* number of points in entry */
   3.380 -  int total_npoints = 0;    /* total number of points in cluster */
   3.381 -  pgl_point *points;        /* points in entry */
   3.382 -  int lon_dir;              /* first point of entry west (-1) or east (+1) */
   3.383 -  double lon_break = 0;     /* antipodal longitude of first point in entry */
   3.384 -  double lon_min, lon_max;  /* covered longitude range of entry */
   3.385 -  double value;             /* temporary variable */
   3.386 -  /* reset bounding circle center to empty circle at 0/0 coordinates */
   3.387 -  cluster->bounding.center.lat = 0;
   3.388 -  cluster->bounding.center.lon = 0;
   3.389 -  cluster->bounding.radius = -INFINITY;
   3.390 -  /* if cluster is not empty */
   3.391 -  if (cluster->nentries != 0) {
   3.392 -    /* iterate over all cluster entries and ensure they each cover a longitude
   3.393 -       range less than 180 degrees */
   3.394 -    for (i=0; i<cluster->nentries; i++) {
   3.395 -      /* get properties of entry */
   3.396 -      npoints = cluster->entries[i].npoints;
   3.397 -      points = PGL_ENTRY_POINTS(cluster, i);
   3.398 -      /* get longitude of first point of entry */
   3.399 -      value = points[0].lon;
   3.400 -      /* initialize lon_min and lon_max with longitude of first point */
   3.401 -      lon_min = value;
   3.402 -      lon_max = value;
   3.403 -      /* determine east/west orientation of first point and calculate antipodal
   3.404 -         longitude (Note: rounding required here) */
   3.405 -      if      (value < 0) { lon_dir = -1; lon_break = pgl_round(value + 180); }
   3.406 -      else if (value > 0) { lon_dir =  1; lon_break = pgl_round(value - 180); }
   3.407 -      else lon_dir = 0;
   3.408 -      /* iterate over all other points in entry */
   3.409 -      for (j=1; j<npoints; j++) {
   3.410 -        /* consider longitude wrap-around */
   3.411 -        value = points[j].lon;
   3.412 -        if      (lon_dir<0 && value>lon_break) value = pgl_round(value - 360);
   3.413 -        else if (lon_dir>0 && value<lon_break) value = pgl_round(value + 360);
   3.414 -        /* update lon_min and lon_max */
   3.415 -        if      (value < lon_min) lon_min = value;
   3.416 -        else if (value > lon_max) lon_max = value;
   3.417 -        /* return false if 180 degrees or more are covered */
   3.418 -        if (lon_max - lon_min >= 180) return false;
   3.419 -      }
   3.420 -    }
   3.421 -    /* iterate over all points of all entries and calculate arbitrary center
   3.422 -       point for bounding circle (best if center point minimizes the radius,
   3.423 -       but some error is allowed here) */
   3.424 -    for (i=0; i<cluster->nentries; i++) {
   3.425 -      /* get properties of entry */
   3.426 -      npoints = cluster->entries[i].npoints;
   3.427 -      points = PGL_ENTRY_POINTS(cluster, i);
   3.428 -      /* check if first entry */
   3.429 -      if (i==0) {
   3.430 -        /* get longitude of first point of first entry in whole cluster */
   3.431 -        value = points[0].lon;
   3.432 -        /* initialize lon_min and lon_max with longitude of first point of
   3.433 -           first entry in whole cluster (used to determine if whole cluster
   3.434 -           covers a longitude range of 180 degrees or more) */
   3.435 -        lon_min = value;
   3.436 -        lon_max = value;
   3.437 -        /* determine east/west orientation of first point and calculate
   3.438 -           antipodal longitude (Note: rounding not necessary here) */
   3.439 -        if      (value < 0) { lon_dir = -1; lon_break = value + 180; }
   3.440 -        else if (value > 0) { lon_dir =  1; lon_break = value - 180; }
   3.441 -        else lon_dir = 0;
   3.442 -      }
   3.443 -      /* iterate over all points in entry */
   3.444 -      for (j=0; j<npoints; j++) {
   3.445 -        /* longitude wrap-around (Note: rounding not necessary here) */
   3.446 -        value = points[j].lon;
   3.447 -        if      (lon_dir < 0 && value > lon_break) value -= 360;
   3.448 -        else if (lon_dir > 0 && value < lon_break) value += 360;
   3.449 -        if      (value < lon_min) lon_min = value;
   3.450 -        else if (value > lon_max) lon_max = value;
   3.451 -        /* set bounding circle to cover whole earth if more than 180 degrees
   3.452 -           are covered */
   3.453 -        if (lon_max - lon_min >= 180) {
   3.454 -          cluster->bounding.center.lat = 0;
   3.455 -          cluster->bounding.center.lon = 0;
   3.456 -          cluster->bounding.radius = INFINITY;
   3.457 -          return true;
   3.458 -        }
   3.459 -        /* add point to bounding circle center (for average calculation) */
   3.460 -        cluster->bounding.center.lat += points[j].lat;
   3.461 -        cluster->bounding.center.lon += value;
   3.462 -      }
   3.463 -      /* count total number of points */
   3.464 -      total_npoints += npoints;
   3.465 -    }
   3.466 -    /* determine average latitude and longitude of cluster */
   3.467 -    cluster->bounding.center.lat /= total_npoints;
   3.468 -    cluster->bounding.center.lon /= total_npoints;
   3.469 -    /* normalize longitude of center of cluster bounding circle */
   3.470 -    if (cluster->bounding.center.lon < -180) {
   3.471 -      cluster->bounding.center.lon += 360;
   3.472 -    }
   3.473 -    else if (cluster->bounding.center.lon > 180) {
   3.474 -      cluster->bounding.center.lon -= 360;
   3.475 -    }
   3.476 -    /* round bounding circle center (useful if it is used by other functions) */
   3.477 -    cluster->bounding.center.lat = pgl_round(cluster->bounding.center.lat);
   3.478 -    cluster->bounding.center.lon = pgl_round(cluster->bounding.center.lon);
   3.479 -    /* calculate radius of bounding circle */
   3.480 -    for (i=0; i<cluster->nentries; i++) {
   3.481 -      npoints = cluster->entries[i].npoints;
   3.482 -      points = PGL_ENTRY_POINTS(cluster, i);
   3.483 -      for (j=0; j<npoints; j++) {
   3.484 -        value = pgl_distance(
   3.485 -          cluster->bounding.center.lat, cluster->bounding.center.lon,
   3.486 -          points[j].lat, points[j].lon
   3.487 -        );
   3.488 -        if (value > cluster->bounding.radius) cluster->bounding.radius = value;
   3.489 -      }
   3.490 -    }
   3.491 -  }
   3.492 -  /* return true (east/west orientation is unambiguous) */
   3.493 -  return true;
   3.494 -}
   3.495 -
   3.496 -/* check if point is inside cluster */
   3.497 -/* (if point is on perimeter, then true is returned if and only if
   3.498 -   strict == false) */
   3.499 -static bool pgl_point_in_cluster(
   3.500 -  pgl_point *point,
   3.501 -  pgl_cluster *cluster,
   3.502 -  bool strict
   3.503 -) {
   3.504 -  int i, j, k;  /* i: entry, j: point in entry, k: next point in entry */
   3.505 -  int entrytype;         /* type of entry */
   3.506 -  int npoints;           /* number of points in entry */
   3.507 -  pgl_point *points;     /* array of points in entry */
   3.508 -  int lon_dir = 0;       /* first vertex west (-1) or east (+1) */
   3.509 -  double lon_break = 0;  /* antipodal longitude of first vertex */
   3.510 -  double lat0 = point->lat;  /* latitude of point */
   3.511 -  double lon0;           /* (adjusted) longitude of point */
   3.512 -  double lat1, lon1;     /* latitude and (adjusted) longitude of vertex */
   3.513 -  double lat2, lon2;     /* latitude and (adjusted) longitude of next vertex */
   3.514 -  double lon;            /* longitude of intersection */
   3.515 -  int counter = 0;       /* counter for intersections east of point */
   3.516 -  /* iterate over all entries */
   3.517 -  for (i=0; i<cluster->nentries; i++) {
   3.518 -    /* get type of entry */
   3.519 -    entrytype = cluster->entries[i].entrytype;
   3.520 -    /* skip all entries but polygons if perimeters are excluded */
   3.521 -    if (strict && entrytype != PGL_ENTRY_POLYGON) continue;
   3.522 -    /* get points of entry */
   3.523 -    npoints = cluster->entries[i].npoints;
   3.524 -    points = PGL_ENTRY_POINTS(cluster, i);
   3.525 -    /* determine east/west orientation of first point of entry and calculate
   3.526 -       antipodal longitude */
   3.527 -    lon_break = points[0].lon;
   3.528 -    if      (lon_break < 0) { lon_dir = -1; lon_break += 180; }
   3.529 -    else if (lon_break > 0) { lon_dir =  1; lon_break -= 180; }
   3.530 -    else lon_dir = 0;
   3.531 -    /* get longitude of point */
   3.532 -    lon0 = point->lon;
   3.533 -    /* consider longitude wrap-around for point */
   3.534 -    if      (lon_dir < 0 && lon0 > lon_break) lon0 = pgl_round(lon0 - 360);
   3.535 -    else if (lon_dir > 0 && lon0 < lon_break) lon0 = pgl_round(lon0 + 360);
   3.536 -    /* iterate over all edges and vertices */
   3.537 -    for (j=0; j<npoints; j++) {
   3.538 -      /* return if point is on vertex of polygon */
   3.539 -      if (pgl_point_cmp(point, &(points[j])) == 0) return !strict;
   3.540 -      /* calculate index of next vertex */
   3.541 -      k = (j+1) % npoints;
   3.542 -      /* skip last edge unless entry is (closed) outline or polygon */
   3.543 -      if (
   3.544 -        k == 0 &&
   3.545 -        entrytype != PGL_ENTRY_OUTLINE &&
   3.546 -        entrytype != PGL_ENTRY_POLYGON
   3.547 -      ) continue;
   3.548 -      /* use previously calculated values for lat1 and lon1 if possible */
   3.549 -      if (j) {
   3.550 -        lat1 = lat2;
   3.551 -        lon1 = lon2;
   3.552 -      } else {
   3.553 -        /* otherwise get latitude and longitude values of first vertex */
   3.554 -        lat1 = points[0].lat;
   3.555 -        lon1 = points[0].lon;
   3.556 -        /* and consider longitude wrap-around for first vertex */
   3.557 -        if      (lon_dir < 0 && lon1 > lon_break) lon1 = pgl_round(lon1 - 360);
   3.558 -        else if (lon_dir > 0 && lon1 < lon_break) lon1 = pgl_round(lon1 + 360);
   3.559 -      }
   3.560 -      /* get latitude and longitude of next vertex */
   3.561 -      lat2 = points[k].lat;
   3.562 -      lon2 = points[k].lon;
   3.563 -      /* consider longitude wrap-around for next vertex */
   3.564 -      if      (lon_dir < 0 && lon2 > lon_break) lon2 = pgl_round(lon2 - 360);
   3.565 -      else if (lon_dir > 0 && lon2 < lon_break) lon2 = pgl_round(lon2 + 360);
   3.566 -      /* return if point is on horizontal (west to east) edge of polygon */
   3.567 -      if (
   3.568 -        lat0 == lat1 && lat0 == lat2 &&
   3.569 -        ( (lon0 >= lon1 && lon0 <= lon2) || (lon0 >= lon2 && lon0 <= lon1) )
   3.570 -      ) return !strict;
   3.571 -      /* check if edge crosses east/west line of point */
   3.572 -      if ((lat1 < lat0 && lat2 >= lat0) || (lat2 < lat0 && lat1 >= lat0)) {
   3.573 -        /* calculate longitude of intersection */
   3.574 -        lon = (lon1 * (lat2-lat0) + lon2 * (lat0-lat1)) / (lat2-lat1);
   3.575 -        /* return if intersection goes (approximately) through point */
   3.576 -        if (pgl_round(lon) == lon0) return !strict;
   3.577 -        /* count intersection if east of point and entry is polygon*/
   3.578 -        if (entrytype == PGL_ENTRY_POLYGON && lon > lon0) counter++;
   3.579 -      }
   3.580 -    }
   3.581 -  }
   3.582 -  /* return true if number of intersections is odd */
   3.583 -  return counter & 1;
   3.584 -}
   3.585 -
   3.586 -/* check if all points of the second cluster are strictly inside the first
   3.587 -   cluster */
   3.588 -static inline bool pgl_all_cluster_points_strictly_in_cluster(
   3.589 -  pgl_cluster *outer, pgl_cluster *inner
   3.590 -) {
   3.591 -  int i, j;           /* i: entry, j: point in entry */
   3.592 -  int npoints;        /* number of points in entry */
   3.593 -  pgl_point *points;  /* array of points in entry */
   3.594 -  /* iterate over all entries of "inner" cluster */
   3.595 -  for (i=0; i<inner->nentries; i++) {
   3.596 -    /* get properties of entry */
   3.597 -    npoints = inner->entries[i].npoints;
   3.598 -    points = PGL_ENTRY_POINTS(inner, i);
   3.599 -    /* iterate over all points in entry of "inner" cluster */
   3.600 -    for (j=0; j<npoints; j++) {
   3.601 -      /* return false if one point of inner cluster is not in outer cluster */
   3.602 -      if (!pgl_point_in_cluster(points+j, outer, true)) return false;
   3.603 -    }
   3.604 -  }
   3.605 -  /* otherwise return true */
   3.606 -  return true;
   3.607 -}
   3.608 -
   3.609 -/* check if any point the second cluster is inside the first cluster */
   3.610 -static inline bool pgl_any_cluster_points_in_cluster(
   3.611 -  pgl_cluster *outer, pgl_cluster *inner
   3.612 -) {
   3.613 -  int i, j;           /* i: entry, j: point in entry */
   3.614 -  int npoints;        /* number of points in entry */
   3.615 -  pgl_point *points;  /* array of points in entry */
   3.616 -  /* iterate over all entries of "inner" cluster */
   3.617 -  for (i=0; i<inner->nentries; i++) {
   3.618 -    /* get properties of entry */
   3.619 -    npoints = inner->entries[i].npoints;
   3.620 -    points = PGL_ENTRY_POINTS(inner, i);
   3.621 -    /* iterate over all points in entry of "inner" cluster */
   3.622 -    for (j=0; j<npoints; j++) {
   3.623 -      /* return true if one point of inner cluster is in outer cluster */
   3.624 -      if (pgl_point_in_cluster(points+j, outer, false)) return true;
   3.625 -    }
   3.626 -  }
   3.627 -  /* otherwise return false */
   3.628 -  return false;
   3.629 -}
   3.630 -
   3.631 -/* check if line segment strictly crosses line (not just touching) */
   3.632 -static inline bool pgl_lseg_crosses_line(
   3.633 -  double seg_x1,  double seg_y1,  double seg_x2,  double seg_y2,
   3.634 -  double line_x1, double line_y1, double line_x2, double line_y2
   3.635 -) {
   3.636 -  return (
   3.637 -    (
   3.638 -      (seg_x1-line_x1) * (line_y2-line_y1) -
   3.639 -      (seg_y1-line_y1) * (line_x2-line_x1)
   3.640 -    ) * (
   3.641 -      (seg_x2-line_x1) * (line_y2-line_y1) -
   3.642 -      (seg_y2-line_y1) * (line_x2-line_x1)
   3.643 -    )
   3.644 -  ) < 0;
   3.645 -}
   3.646 -
   3.647 -/* check if paths and outlines of two clusters strictly overlap (not just
   3.648 -   touching) */
   3.649 -static bool pgl_outlines_overlap(
   3.650 -  pgl_cluster *cluster1, pgl_cluster *cluster2
   3.651 -) {
   3.652 -  int i1, j1, k1;  /* i: entry, j: point in entry, k: next point in entry */
   3.653 -  int i2, j2, k2;
   3.654 -  int entrytype1, entrytype2;     /* type of entry */
   3.655 -  int npoints1, npoints2;         /* number of points in entry */
   3.656 -  pgl_point *points1;             /* array of points in entry of cluster1 */
   3.657 -  pgl_point *points2;             /* array of points in entry of cluster2 */
   3.658 -  int lon_dir1, lon_dir2;         /* first vertex west (-1) or east (+1) */
   3.659 -  double lon_break1, lon_break2;  /* antipodal longitude of first vertex */
   3.660 -  double lat11, lon11;  /* latitude and (adjusted) longitude of vertex */
   3.661 -  double lat12, lon12;  /* latitude and (adjusted) longitude of next vertex */
   3.662 -  double lat21, lon21;  /* latitude and (adjusted) longitudes for cluster2 */
   3.663 -  double lat22, lon22;
   3.664 -  double wrapvalue;     /* temporary helper value to adjust wrap-around */  
   3.665 -  /* iterate over all entries of cluster1 */
   3.666 -  for (i1=0; i1<cluster1->nentries; i1++) {
   3.667 -    /* get properties of entry in cluster1 and skip points */
   3.668 -    npoints1 = cluster1->entries[i1].npoints;
   3.669 -    if (npoints1 < 2) continue;
   3.670 -    entrytype1 = cluster1->entries[i1].entrytype;
   3.671 -    points1 = PGL_ENTRY_POINTS(cluster1, i1);
   3.672 -    /* determine east/west orientation of first point and calculate antipodal
   3.673 -       longitude */
   3.674 -    lon_break1 = points1[0].lon;
   3.675 -    if (lon_break1 < 0) {
   3.676 -      lon_dir1   = -1;
   3.677 -      lon_break1 = pgl_round(lon_break1 + 180);
   3.678 -    } else if (lon_break1 > 0) {
   3.679 -      lon_dir1   = 1;
   3.680 -      lon_break1 = pgl_round(lon_break1 - 180);
   3.681 -    } else lon_dir1 = 0;
   3.682 -    /* iterate over all edges and vertices in cluster1 */
   3.683 -    for (j1=0; j1<npoints1; j1++) {
   3.684 -      /* calculate index of next vertex */
   3.685 -      k1 = (j1+1) % npoints1;
   3.686 -      /* skip last edge unless entry is (closed) outline or polygon */
   3.687 -      if (
   3.688 -        k1 == 0 &&
   3.689 -        entrytype1 != PGL_ENTRY_OUTLINE &&
   3.690 -        entrytype1 != PGL_ENTRY_POLYGON
   3.691 -      ) continue;
   3.692 -      /* use previously calculated values for lat1 and lon1 if possible */
   3.693 -      if (j1) {
   3.694 -        lat11 = lat12;
   3.695 -        lon11 = lon12;
   3.696 -      } else {
   3.697 -        /* otherwise get latitude and longitude values of first vertex */
   3.698 -        lat11 = points1[0].lat;
   3.699 -        lon11 = points1[0].lon;
   3.700 -        /* and consider longitude wrap-around for first vertex */
   3.701 -        if      (lon_dir1<0 && lon11>lon_break1) lon11 = pgl_round(lon11-360);
   3.702 -        else if (lon_dir1>0 && lon11<lon_break1) lon11 = pgl_round(lon11+360);
   3.703 -      }
   3.704 -      /* get latitude and longitude of next vertex */
   3.705 -      lat12 = points1[k1].lat;
   3.706 -      lon12 = points1[k1].lon;
   3.707 -      /* consider longitude wrap-around for next vertex */
   3.708 -      if      (lon_dir1<0 && lon12>lon_break1) lon12 = pgl_round(lon12-360);
   3.709 -      else if (lon_dir1>0 && lon12<lon_break1) lon12 = pgl_round(lon12+360);
   3.710 -      /* skip degenerated edges */
   3.711 -      if (lat11 == lat12 && lon11 == lon12) continue;
   3.712 -      /* iterate over all entries of cluster2 */
   3.713 -      for (i2=0; i2<cluster2->nentries; i2++) {
   3.714 -        /* get points and number of points of entry in cluster2 */
   3.715 -        npoints2 = cluster2->entries[i2].npoints;
   3.716 -        if (npoints2 < 2) continue;
   3.717 -        entrytype2 = cluster2->entries[i2].entrytype;
   3.718 -        points2 = PGL_ENTRY_POINTS(cluster2, i2);
   3.719 -        /* determine east/west orientation of first point and calculate antipodal
   3.720 -           longitude */
   3.721 -        lon_break2 = points2[0].lon;
   3.722 -        if (lon_break2 < 0) {
   3.723 -          lon_dir2   = -1;
   3.724 -          lon_break2 = pgl_round(lon_break2 + 180);
   3.725 -        } else if (lon_break2 > 0) {
   3.726 -          lon_dir2   = 1;
   3.727 -          lon_break2 = pgl_round(lon_break2 - 180);
   3.728 -        } else lon_dir2 = 0;
   3.729 -        /* iterate over all edges and vertices in cluster2 */
   3.730 -        for (j2=0; j2<npoints2; j2++) {
   3.731 -          /* calculate index of next vertex */
   3.732 -          k2 = (j2+1) % npoints2;
   3.733 -          /* skip last edge unless entry is (closed) outline or polygon */
   3.734 -          if (
   3.735 -            k2 == 0 &&
   3.736 -            entrytype2 != PGL_ENTRY_OUTLINE &&
   3.737 -            entrytype2 != PGL_ENTRY_POLYGON
   3.738 -          ) continue;
   3.739 -          /* use previously calculated values for lat1 and lon1 if possible */
   3.740 -          if (j2) {
   3.741 -            lat21 = lat22;
   3.742 -            lon21 = lon22;
   3.743 -          } else {
   3.744 -            /* otherwise get latitude and longitude values of first vertex */
   3.745 -            lat21 = points2[0].lat;
   3.746 -            lon21 = points2[0].lon;
   3.747 -            /* and consider longitude wrap-around for first vertex */
   3.748 -            if      (lon_dir2<0 && lon21>lon_break2) lon21 = pgl_round(lon21-360);
   3.749 -            else if (lon_dir2>0 && lon21<lon_break2) lon21 = pgl_round(lon21+360);
   3.750 -          }
   3.751 -          /* get latitude and longitude of next vertex */
   3.752 -          lat22 = points2[k2].lat;
   3.753 -          lon22 = points2[k2].lon;
   3.754 -          /* consider longitude wrap-around for next vertex */
   3.755 -          if      (lon_dir2<0 && lon22>lon_break2) lon22 = pgl_round(lon22-360);
   3.756 -          else if (lon_dir2>0 && lon22<lon_break2) lon22 = pgl_round(lon22+360);
   3.757 -          /* skip degenerated edges */
   3.758 -          if (lat21 == lat22 && lon21 == lon22) continue;
   3.759 -          /* perform another wrap-around where necessary */
   3.760 -          /* TODO: improve performance of whole wrap-around mechanism */
   3.761 -          wrapvalue = (lon21 + lon22) - (lon11 + lon12);
   3.762 -          if (wrapvalue > 360) {
   3.763 -            lon21 = pgl_round(lon21 - 360);
   3.764 -            lon22 = pgl_round(lon22 - 360);
   3.765 -          } else if (wrapvalue < -360) {
   3.766 -            lon21 = pgl_round(lon21 + 360);
   3.767 -            lon22 = pgl_round(lon22 + 360);
   3.768 -          }
   3.769 -          /* return true if segments overlap */
   3.770 -          if (
   3.771 -            pgl_lseg_crosses_line(
   3.772 -              lat11, lon11, lat12, lon12,
   3.773 -              lat21, lon21, lat22, lon22
   3.774 -            ) && pgl_lseg_crosses_line(
   3.775 -              lat21, lon21, lat22, lon22,
   3.776 -              lat11, lon11, lat12, lon12
   3.777 -            )
   3.778 -          ) {
   3.779 -            return true;
   3.780 -          }
   3.781 -        }
   3.782 -      }
   3.783 -    }
   3.784 -  }
   3.785 -  /* otherwise return false */
   3.786 -  return false;
   3.787 -}
   3.788 -
   3.789 -/* check if second cluster is completely contained in first cluster */
   3.790 -static bool pgl_cluster_in_cluster(pgl_cluster *outer, pgl_cluster *inner) {
   3.791 -  if (!pgl_all_cluster_points_strictly_in_cluster(outer, inner)) return false;
   3.792 -  if (pgl_any_cluster_points_in_cluster(inner, outer)) return false;
   3.793 -  if (pgl_outlines_overlap(outer, inner)) return false;
   3.794 -  return true;
   3.795 -}
   3.796 -
   3.797 -/* check if two clusters overlap */
   3.798 -static bool pgl_clusters_overlap(
   3.799 -  pgl_cluster *cluster1, pgl_cluster *cluster2
   3.800 -) {
   3.801 -  if (pgl_any_cluster_points_in_cluster(cluster1, cluster2)) return true;
   3.802 -  if (pgl_any_cluster_points_in_cluster(cluster2, cluster1)) return true;
   3.803 -  if (pgl_outlines_overlap(cluster1, cluster2)) return true;
   3.804 -  return false;
   3.805 -}
   3.806 -
   3.807 -
   3.808 -/* calculate (approximate) distance between point and cluster */
   3.809 -static double pgl_point_cluster_distance(pgl_point *point, pgl_cluster *cluster) {
   3.810 -  double comp;           /* square of compression of meridians */
   3.811 -  int i, j, k;  /* i: entry, j: point in entry, k: next point in entry */
   3.812 -  int entrytype;         /* type of entry */
   3.813 -  int npoints;           /* number of points in entry */
   3.814 -  pgl_point *points;     /* array of points in entry */
   3.815 -  int lon_dir = 0;       /* first vertex west (-1) or east (+1) */
   3.816 -  double lon_break = 0;  /* antipodal longitude of first vertex */
   3.817 -  double lon_min = 0;    /* minimum (adjusted) longitude of entry vertices */
   3.818 -  double lon_max = 0;    /* maximum (adjusted) longitude of entry vertices */
   3.819 -  double lat0 = point->lat;  /* latitude of point */
   3.820 -  double lon0;           /* (adjusted) longitude of point */
   3.821 -  double lat1, lon1;     /* latitude and (adjusted) longitude of vertex */
   3.822 -  double lat2, lon2;     /* latitude and (adjusted) longitude of next vertex */
   3.823 -  double s;              /* scalar for vector calculations */
   3.824 -  double dist;           /* distance calculated in one step */
   3.825 -  double min_dist = INFINITY;   /* minimum distance */
   3.826 -  /* distance is zero if point is contained in cluster */
   3.827 -  if (pgl_point_in_cluster(point, cluster, false)) return 0;
   3.828 -  /* calculate approximate square compression of meridians */
   3.829 -  comp = cos((lat0 / 180.0) * M_PI);
   3.830 -  comp *= comp;
   3.831 -  /* calculate exact square compression of meridians */
   3.832 -  comp *= (
   3.833 -    (1.0 - PGL_EPS2 * (1.0-comp)) *
   3.834 -    (1.0 - PGL_EPS2 * (1.0-comp)) /
   3.835 -    (PGL_SUBEPS2 * PGL_SUBEPS2)
   3.836 -  );
   3.837 -  /* iterate over all entries */
   3.838 -  for (i=0; i<cluster->nentries; i++) {
   3.839 -    /* get properties of entry */
   3.840 -    entrytype = cluster->entries[i].entrytype;
   3.841 -    npoints = cluster->entries[i].npoints;
   3.842 -    points = PGL_ENTRY_POINTS(cluster, i);
   3.843 -    /* determine east/west orientation of first point of entry and calculate
   3.844 -       antipodal longitude */
   3.845 -    lon_break = points[0].lon;
   3.846 -    if      (lon_break < 0) { lon_dir = -1; lon_break += 180; }
   3.847 -    else if (lon_break > 0) { lon_dir =  1; lon_break -= 180; }
   3.848 -    else lon_dir = 0;
   3.849 -    /* determine covered longitude range */
   3.850 -    for (j=0; j<npoints; j++) {
   3.851 -      /* get longitude of vertex */
   3.852 -      lon1 = points[j].lon;
   3.853 -      /* adjust longitude to fix potential wrap-around */
   3.854 -      if      (lon_dir < 0 && lon1 > lon_break) lon1 -= 360;
   3.855 -      else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
   3.856 -      /* update minimum and maximum longitude of polygon */
   3.857 -      if (j == 0 || lon1 < lon_min) lon_min = lon1;
   3.858 -      if (j == 0 || lon1 > lon_max) lon_max = lon1;
   3.859 -    }
   3.860 -    /* adjust longitude wrap-around according to full longitude range */
   3.861 -    lon_break = (lon_max + lon_min) / 2;
   3.862 -    if      (lon_break < 0) { lon_dir = -1; lon_break += 180; }
   3.863 -    else if (lon_break > 0) { lon_dir =  1; lon_break -= 180; }
   3.864 -    /* get longitude of point */
   3.865 -    lon0 = point->lon;
   3.866 -    /* consider longitude wrap-around for point */
   3.867 -    if      (lon_dir < 0 && lon0 > lon_break) lon0 -= 360;
   3.868 -    else if (lon_dir > 0 && lon0 < lon_break) lon0 += 360;
   3.869 -    /* iterate over all edges and vertices */
   3.870 -    for (j=0; j<npoints; j++) {
   3.871 -      /* use previously calculated values for lat1 and lon1 if possible */
   3.872 -      if (j) {
   3.873 -        lat1 = lat2;
   3.874 -        lon1 = lon2;
   3.875 -      } else {
   3.876 -        /* otherwise get latitude and longitude values of first vertex */
   3.877 -        lat1 = points[0].lat;
   3.878 -        lon1 = points[0].lon;
   3.879 -        /* and consider longitude wrap-around for first vertex */
   3.880 -        if      (lon_dir < 0 && lon1 > lon_break) lon1 -= 360;
   3.881 -        else if (lon_dir > 0 && lon1 < lon_break) lon1 += 360;
   3.882 -      }
   3.883 -      /* calculate distance to vertex */
   3.884 -      dist = pgl_distance(lat0, lon0, lat1, lon1);
   3.885 -      /* store calculated distance if smallest */
   3.886 -      if (dist < min_dist) min_dist = dist;
   3.887 -      /* calculate index of next vertex */
   3.888 -      k = (j+1) % npoints;
   3.889 -      /* skip last edge unless entry is (closed) outline or polygon */
   3.890 -      if (
   3.891 -        k == 0 &&
   3.892 -        entrytype != PGL_ENTRY_OUTLINE &&
   3.893 -        entrytype != PGL_ENTRY_POLYGON
   3.894 -      ) continue;
   3.895 -      /* get latitude and longitude of next vertex */
   3.896 -      lat2 = points[k].lat;
   3.897 -      lon2 = points[k].lon;
   3.898 -      /* consider longitude wrap-around for next vertex */
   3.899 -      if      (lon_dir < 0 && lon2 > lon_break) lon2 -= 360;
   3.900 -      else if (lon_dir > 0 && lon2 < lon_break) lon2 += 360;
   3.901 -      /* go to next vertex and edge if edge is degenerated */
   3.902 -      if (lat1 == lat2 && lon1 == lon2) continue;
   3.903 -      /* otherwise test if point can be projected onto edge of polygon */
   3.904 -      s = (
   3.905 -        ((lat0-lat1) * (lat2-lat1) + comp * (lon0-lon1) * (lon2-lon1)) /
   3.906 -        ((lat2-lat1) * (lat2-lat1) + comp * (lon2-lon1) * (lon2-lon1))
   3.907 -      );
   3.908 -      /* go to next vertex and edge if point cannot be projected */
   3.909 -      if (!(s > 0 && s < 1)) continue;
   3.910 -      /* calculate distance from original point to projected point */
   3.911 -      dist = pgl_distance(
   3.912 -        lat0, lon0,
   3.913 -        lat1 + s * (lat2-lat1),
   3.914 -        lon1 + s * (lon2-lon1)
   3.915 -      );
   3.916 -      /* store calculated distance if smallest */
   3.917 -      if (dist < min_dist) min_dist = dist;
   3.918 -    }
   3.919 -  }
   3.920 -  /* return minimum distance */
   3.921 -  return min_dist;
   3.922 -}
   3.923 -
   3.924 -/* calculate (approximate) distance between two clusters */
   3.925 -static double pgl_cluster_distance(pgl_cluster *cluster1, pgl_cluster *cluster2) {
   3.926 -  int i, j;                    /* i: entry, j: point in entry */
   3.927 -  int npoints;                 /* number of points in entry */
   3.928 -  pgl_point *points;           /* array of points in entry */
   3.929 -  double dist;                 /* distance calculated in one step */
   3.930 -  double min_dist = INFINITY;  /* minimum distance */
   3.931 -  /* consider distance from each point in one cluster to the whole other */
   3.932 -  for (i=0; i<cluster1->nentries; i++) {
   3.933 -    npoints = cluster1->entries[i].npoints;
   3.934 -    points = PGL_ENTRY_POINTS(cluster1, i);
   3.935 -    for (j=0; j<npoints; j++) {
   3.936 -      dist = pgl_point_cluster_distance(points+j, cluster2);
   3.937 -      if (dist == 0) return dist;
   3.938 -      if (dist < min_dist) min_dist = dist;
   3.939 -    }
   3.940 -  }
   3.941 -  /* consider distance from each point in other cluster to the first cluster */
   3.942 -  for (i=0; i<cluster2->nentries; i++) {
   3.943 -    npoints = cluster2->entries[i].npoints;
   3.944 -    points = PGL_ENTRY_POINTS(cluster2, i);
   3.945 -    for (j=0; j<npoints; j++) {
   3.946 -      dist = pgl_point_cluster_distance(points+j, cluster1);
   3.947 -      if (dist == 0) return dist;
   3.948 -      if (dist < min_dist) min_dist = dist;
   3.949 -    }
   3.950 -  }
   3.951 -  return min_dist;
   3.952 -}
   3.953 -
   3.954 -/* estimator function for distance between box and point */
   3.955 -/* always returns a smaller value than actually correct or zero */
   3.956 -static double pgl_estimate_point_box_distance(pgl_point *point, pgl_box *box) {
   3.957 -  double dlon;      /* longitude range of box (delta longitude) */
   3.958 -  double distance;  /* return value */
   3.959 -  /* return infinity if box is empty */
   3.960 -  if (box->lat_min > box->lat_max) return INFINITY;
   3.961 -  /* return zero if point is inside box */
   3.962 -  if (pgl_point_in_box(point, box)) return 0;
   3.963 -  /* calculate delta longitude */
   3.964 -  dlon = box->lon_max - box->lon_min;
   3.965 -  if (dlon < 0) dlon += 360;  /* 180th meridian crossed */
   3.966 -  /* if delta longitude is greater than 150 degrees, perform safe fall-back */
   3.967 -  if (dlon > 150) return 0;
   3.968 -  /* calculate lower limit for distance (formula below requires dlon <= 150) */
   3.969 -  /* TODO: provide better estimation function to improve performance */
   3.970 -  distance = (
   3.971 -    (1.0-1e-14) * pgl_distance(
   3.972 -      point->lat,
   3.973 -      point->lon,
   3.974 -      (box->lat_min + box->lat_max) / 2,
   3.975 -      box->lon_min + dlon/2
   3.976 -    ) - pgl_distance(
   3.977 -      box->lat_min, box->lon_min,
   3.978 -      box->lat_max, box->lon_max
   3.979 -    )
   3.980 -  );
   3.981 -  /* truncate negative results to zero */
   3.982 -  if (distance <= 0) distance = 0;
   3.983 -  /* return result */
   3.984 -  return distance;
   3.985 -}
   3.986 -
   3.987 -
   3.988 -/*------------------------------------------------------------*
   3.989 - *  Functions using numerical integration (Monte Carlo like)  *
   3.990 - *------------------------------------------------------------*/
   3.991 -
   3.992 -/* half of (spherical) earth's surface area */
   3.993 -#define PGL_HALF_SURFACE (PGL_RADIUS * PGL_DIAMETER * M_PI)
   3.994 -
   3.995 -/* golden angle in radians */
   3.996 -#define PGL_GOLDEN_ANGLE (M_PI * (sqrt(5) - 1.0))
   3.997 -
   3.998 -/* create a list of sample points covering a bounding circle 
   3.999 -   and return covered area */
  3.1000 -static double pgl_sample_points(
  3.1001 -  pgl_point *center,  /* center of bounding circle */
  3.1002 -  double radius,      /* radius of bounding circle */
  3.1003 -  int samples,        /* number of sample points */
  3.1004 -  pgl_point *result   /* pointer to result array */
  3.1005 -) {
  3.1006 -  double double_share = 2.0;  /* double of covered share of earth's surface */
  3.1007 -  double double_share_div_samples;  /* double_share divided by sample count */
  3.1008 -  int i;
  3.1009 -  double t;  /* parameter of spiral laid on (spherical) earth's surface */
  3.1010 -  double x, y, z;  /* normalized coordinates of point on non-rotated spiral */
  3.1011 -  double sin_phi;  /* sine of sph. coordinate of point of non-rotated spiral */
  3.1012 -  double lambda;   /* other sph. coordinate of point of non-rotated spiral */
  3.1013 -  double rot = (0.5 - center->lat / 180.0) * M_PI;  /* needed rot. (in rad) */
  3.1014 -  double cos_rot = cos(rot);  /* cosine of rotation by latitude */
  3.1015 -  double sin_rot = sin(rot);  /* sine of rotation by latitude */
  3.1016 -  double x_rot, z_rot;  /* normalized coordinates of point on rotated spiral */
  3.1017 -  double center_lon = center->lon;  /* second rotation in degree */
  3.1018 -  /* add safety margin to bounding circle because of spherical approximation */
  3.1019 -  radius *= PGL_SPHEROID_A / PGL_RADIUS;
  3.1020 -  /* if whole earth is covered, use initialized value, otherwise calculate
  3.1021 -     share of covered area (multiplied by 2) */
  3.1022 -  if (radius < PGL_MAXDIST) double_share = 1.0 - cos(radius / PGL_RADIUS);
  3.1023 -  /* divide double_share by sample count for later calculations */
  3.1024 -  double_share_div_samples = double_share / samples;
  3.1025 -  /* generate sample points */
  3.1026 -  for (i=0; i<samples; i++) {
  3.1027 -    /* use an offset of 1/2 to avoid too dense clustering at spiral center */
  3.1028 -    t = 0.5 + i;
  3.1029 -    /* calculate normalized coordinates of point on non-rotated spiral */
  3.1030 -    z = 1.0 - double_share_div_samples * t;
  3.1031 -    sin_phi = sqrt(1.0 - z*z);
  3.1032 -    lambda = t * PGL_GOLDEN_ANGLE;
  3.1033 -    x = sin_phi * cos(lambda);
  3.1034 -    y = sin_phi * sin(lambda);
  3.1035 -    /* rotate spiral by latitude value of bounding circle */
  3.1036 -    x_rot = cos_rot * x + sin_rot * z;
  3.1037 -    z_rot = cos_rot * z - sin_rot * x;
  3.1038 -    /* set resulting sample point in result array */
  3.1039 -    /* (while performing second rotation by bounding circle longitude) */
  3.1040 -    result[i].lat = 180.0 * (atan(z_rot / fabs(x_rot)) / M_PI);
  3.1041 -    result[i].lon = center_lon + 180.0 * (atan2(y, x_rot) / M_PI);
  3.1042 -  }
  3.1043 -  /* return covered area */
  3.1044 -  return PGL_HALF_SURFACE * double_share;
  3.1045 -}
  3.1046 -
  3.1047 -/* calculate covered area using Monte Carlo like method */
  3.1048 -/* TODO: inefficient, should be replaced by different method */
  3.1049 -static double pgl_monte_carlo_area(pgl_cluster *cluster, int samples) {
  3.1050 -  pgl_point *points;
  3.1051 -  double area;
  3.1052 -  int i;
  3.1053 -  int matches = 0;
  3.1054 -  if (samples > PGL_CLUSTER_MAXPOINTS) {
  3.1055 -    ereport(ERROR, (
  3.1056 -      errcode(ERRCODE_DATA_EXCEPTION),
  3.1057 -      errmsg(
  3.1058 -        "too many sample points for Monte Carlo method (maximum %i)",
  3.1059 -        PGL_CLUSTER_MAXPOINTS
  3.1060 -      )
  3.1061 -    ));
  3.1062 -  }
  3.1063 -  if (!(cluster->bounding.radius > 0)) return 0;
  3.1064 -  for (i=0; ; i++) {
  3.1065 -    if (i == cluster->nentries) return 0;
  3.1066 -    if (cluster->entries[i].entrytype == PGL_ENTRY_POLYGON) break;
  3.1067 -  }
  3.1068 -  points = palloc(samples * sizeof(pgl_point));
  3.1069 -  area = pgl_sample_points(
  3.1070 -    &cluster->bounding.center,
  3.1071 -    cluster->bounding.radius,
  3.1072 -    samples,
  3.1073 -    points
  3.1074 -  );
  3.1075 -  for (i=0; i<samples; i++) {
  3.1076 -    if (pgl_point_in_cluster(points+i, cluster, true)) matches++;
  3.1077 -  }
  3.1078 -  pfree(points);
  3.1079 -  return area * ((double)matches / samples);
  3.1080 -}
  3.1081 -
  3.1082 -/* fair distance between point and cluster (see README file for explanation) */
  3.1083 -static double pgl_fair_distance(
  3.1084 -  pgl_point *point, pgl_cluster *cluster, int samples
  3.1085 -) {
  3.1086 -  double distance;       /* shortest distance from point to cluster */
  3.1087 -  pgl_point *points;     /* sample points for Monte Carlo method */
  3.1088 -  double area;           /* area covered by sample points */
  3.1089 -  int i;
  3.1090 -  int matches = 0;       /* number of points within enlarged area (multiplied
  3.1091 -                            by two to be able to store half matches) */
  3.1092 -  double result;         /* result determined by Monte Carlo method */
  3.1093 -  /* limit sample count to avoid integer overflows on memory allocation */
  3.1094 -  if (samples > PGL_CLUSTER_MAXPOINTS) {
  3.1095 -    ereport(ERROR, (
  3.1096 -      errcode(ERRCODE_DATA_EXCEPTION),
  3.1097 -      errmsg(
  3.1098 -        "too many sample points for Monte Carlo method (maximum %i)",
  3.1099 -        PGL_CLUSTER_MAXPOINTS
  3.1100 -      )
  3.1101 -    ));
  3.1102 -  }
  3.1103 -  /* calculate shortest distance from point to cluster */
  3.1104 -  distance = pgl_point_cluster_distance(point, cluster);
  3.1105 -  /* if cluster consists of a single point or has no bounding circle with
  3.1106 -      positive radius, simply return distance */
  3.1107 -  if (
  3.1108 -    (cluster->nentries==1 && cluster->entries[0].entrytype==PGL_ENTRY_POINT) ||
  3.1109 -    !(cluster->bounding.radius > 0)
  3.1110 -  ) return distance;
  3.1111 -  /* if cluster consists of two points which are twice as far apart, return
  3.1112 -     distance between point and cluster multiplied by square root of two */
  3.1113 -  if (
  3.1114 -    cluster->nentries == 2 &&
  3.1115 -    cluster->entries[0].entrytype == PGL_ENTRY_POINT &&
  3.1116 -    cluster->entries[1].entrytype == PGL_ENTRY_POINT &&
  3.1117 -    pgl_distance(
  3.1118 -      PGL_ENTRY_POINTS(cluster, 0)[0].lat,
  3.1119 -      PGL_ENTRY_POINTS(cluster, 0)[0].lon,
  3.1120 -      PGL_ENTRY_POINTS(cluster, 1)[0].lat,
  3.1121 -      PGL_ENTRY_POINTS(cluster, 1)[0].lon
  3.1122 -    ) >= 2.0 * distance
  3.1123 -  ) {
  3.1124 -    return distance * M_SQRT2;
  3.1125 -  }
  3.1126 -  /* otherwise create sample points for Monte Carlo method and determine area
  3.1127 -     covered by sample points */
  3.1128 -  points = palloc(samples * sizeof(pgl_point));
  3.1129 -  area = pgl_sample_points(
  3.1130 -    &cluster->bounding.center,
  3.1131 -    cluster->bounding.radius + distance,  /* pad bounding circle by distance */
  3.1132 -    samples,
  3.1133 -    points
  3.1134 -  );
  3.1135 -  /* perform Monte Carlo method */
  3.1136 -  if (distance > 0) {
  3.1137 -    /* point is outside cluster */
  3.1138 -    for (i=0; i<samples; i++) {
  3.1139 -      /* count sample poitns within cluster as half match */
  3.1140 -      if (pgl_point_in_cluster(points+i, cluster, true)) matches += 1;
  3.1141 -      /* count sample points outside of cluster but within cluster grown by
  3.1142 -         distance as full match */
  3.1143 -      else if (
  3.1144 -        pgl_point_cluster_distance(points+i, cluster) < distance
  3.1145 -      ) matches += 2;
  3.1146 -    }
  3.1147 -  } else {
  3.1148 -    /* if point is within cluster,
  3.1149 -       just count sample points within cluster as half match */
  3.1150 -    for (i=0; i<samples; i++) {
  3.1151 -      if (pgl_point_in_cluster(points+i, cluster, true)) matches += 1;
  3.1152 -    }
  3.1153 -  }
  3.1154 -  /* release memory for sample points needed by Monte Carlo method */
  3.1155 -  pfree(points);
  3.1156 -  /* convert area determined by Monte Carlo method into distance
  3.1157 -     (i.e. radius of circle with the same area) */
  3.1158 -  result = sqrt(area * ((double)matches / 2.0 / samples) / M_PI);
  3.1159 -  /* return result only if it is greater than the distance between point and
  3.1160 -     cluster to avoid unexpected results because of errors due to limited
  3.1161 -     precision */
  3.1162 -  if (result > distance) return result;
  3.1163 -  /* otherwise return distance between point and cluster */
  3.1164 -  else return distance;
  3.1165 -}
  3.1166 -
  3.1167 -
  3.1168 -/*-------------------------------------------------*
  3.1169 - *  geographic index based on space-filling curve  *
  3.1170 - *-------------------------------------------------*/
  3.1171 -
  3.1172 -/* number of bytes used for geographic (center) position in keys */
  3.1173 -#define PGL_KEY_LATLON_BYTELEN 7
  3.1174 -
  3.1175 -/* maximum reference value for logarithmic size of geographic objects */
  3.1176 -#define PGL_AREAKEY_REFOBJSIZE (PGL_DIAMETER/3.0)  /* can be tweaked */
  3.1177 -
  3.1178 -/* pointer to index key (either pgl_pointkey or pgl_areakey) */
  3.1179 -typedef unsigned char *pgl_keyptr;
  3.1180 -
  3.1181 -/* index key for points (objects with zero area) on the spheroid */
  3.1182 -/* bit  0..55: interspersed bits of latitude and longitude,
  3.1183 -   bit 56..57: always zero,
  3.1184 -   bit 58..63: node depth in hypothetic (full) tree from 0 to 56 (incl.) */
  3.1185 -typedef unsigned char pgl_pointkey[PGL_KEY_LATLON_BYTELEN+1];
  3.1186 -
  3.1187 -/* index key for geographic objects on spheroid with area greater than zero */
  3.1188 -/* bit  0..55: interspersed bits of latitude and longitude of center point,
  3.1189 -   bit     56: always set to 1,
  3.1190 -   bit 57..63: node depth in hypothetic (full) tree from 0 to (2*56)+1 (incl.),
  3.1191 -   bit 64..71: logarithmic object size from 0 to 56+1 = 57 (incl.), but set to
  3.1192 -               PGL_KEY_OBJSIZE_EMPTY (with interspersed bits = 0 and node depth
  3.1193 -               = 113) for empty objects, and set to PGL_KEY_OBJSIZE_UNIVERSAL
  3.1194 -               (with interspersed bits = 0 and node depth = 0) for keys which
  3.1195 -               cover both empty and non-empty objects */
  3.1196 -
  3.1197 -typedef unsigned char pgl_areakey[PGL_KEY_LATLON_BYTELEN+2];
  3.1198 -
  3.1199 -/* helper macros for reading/writing index keys */
  3.1200 -#define PGL_KEY_NODEDEPTH_OFFSET  PGL_KEY_LATLON_BYTELEN
  3.1201 -#define PGL_KEY_OBJSIZE_OFFSET    (PGL_KEY_NODEDEPTH_OFFSET+1)
  3.1202 -#define PGL_POINTKEY_MAXDEPTH     (PGL_KEY_LATLON_BYTELEN*8)
  3.1203 -#define PGL_AREAKEY_MAXDEPTH      (2*PGL_POINTKEY_MAXDEPTH+1)
  3.1204 -#define PGL_AREAKEY_MAXOBJSIZE    (PGL_POINTKEY_MAXDEPTH+1)
  3.1205 -#define PGL_AREAKEY_TYPEMASK      0x80
  3.1206 -#define PGL_KEY_LATLONBIT(key, n) ((key)[(n)/8] & (0x80 >> ((n)%8)))
  3.1207 -#define PGL_KEY_LATLONBIT_DIFF(key1, key2, n) \
  3.1208 -                                  ( PGL_KEY_LATLONBIT(key1, n) ^ \
  3.1209 -                                    PGL_KEY_LATLONBIT(key2, n) )
  3.1210 -#define PGL_KEY_IS_AREAKEY(key)   ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \
  3.1211 -                                    PGL_AREAKEY_TYPEMASK)
  3.1212 -#define PGL_KEY_NODEDEPTH(key)    ((key)[PGL_KEY_NODEDEPTH_OFFSET] & \
  3.1213 -                                    (PGL_AREAKEY_TYPEMASK-1))
  3.1214 -#define PGL_KEY_OBJSIZE(key)      ((key)[PGL_KEY_OBJSIZE_OFFSET])
  3.1215 -#define PGL_KEY_OBJSIZE_EMPTY     126
  3.1216 -#define PGL_KEY_OBJSIZE_UNIVERSAL 127
  3.1217 -#define PGL_KEY_IS_EMPTY(key)     ( PGL_KEY_IS_AREAKEY(key) && \
  3.1218 -                                    (key)[PGL_KEY_OBJSIZE_OFFSET] == \
  3.1219 -                                    PGL_KEY_OBJSIZE_EMPTY )
  3.1220 -#define PGL_KEY_IS_UNIVERSAL(key) ( PGL_KEY_IS_AREAKEY(key) && \
  3.1221 -                                    (key)[PGL_KEY_OBJSIZE_OFFSET] == \
  3.1222 -                                    PGL_KEY_OBJSIZE_UNIVERSAL )
  3.1223 -
  3.1224 -/* set area key to match empty objects only */
  3.1225 -static void pgl_key_set_empty(pgl_keyptr key) {
  3.1226 -  memset(key, 0, sizeof(pgl_areakey));
  3.1227 -  /* Note: setting node depth to maximum is required for picksplit function */
  3.1228 -  key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH;
  3.1229 -  key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_EMPTY;
  3.1230 -}
  3.1231 -
  3.1232 -/* set area key to match any object (including empty objects) */
  3.1233 -static void pgl_key_set_universal(pgl_keyptr key) {
  3.1234 -  memset(key, 0, sizeof(pgl_areakey));
  3.1235 -  key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK;
  3.1236 -  key[PGL_KEY_OBJSIZE_OFFSET] = PGL_KEY_OBJSIZE_UNIVERSAL;
  3.1237 -}
  3.1238 -
  3.1239 -/* convert a point on earth into a max-depth key to be used in index */
  3.1240 -static void pgl_point_to_key(pgl_point *point, pgl_keyptr key) {
  3.1241 -  double lat = point->lat;
  3.1242 -  double lon = point->lon;
  3.1243 -  int i;
  3.1244 -  /* clear latitude and longitude bits */
  3.1245 -  memset(key, 0, PGL_KEY_LATLON_BYTELEN);
  3.1246 -  /* set node depth to maximum and type bit to zero */
  3.1247 -  key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_POINTKEY_MAXDEPTH;
  3.1248 -  /* iterate over all latitude/longitude bit pairs */
  3.1249 -  for (i=0; i<PGL_POINTKEY_MAXDEPTH/2; i++) {
  3.1250 -    /* determine latitude bit */
  3.1251 -    if (lat >= 0) {
  3.1252 -      key[i/4] |= 0x80 >> (2*(i%4));
  3.1253 -      lat *= 2; lat -= 90;
  3.1254 -    } else {
  3.1255 -      lat *= 2; lat += 90;
  3.1256 -    }
  3.1257 -    /* determine longitude bit */
  3.1258 -    if (lon >= 0) {
  3.1259 -      key[i/4] |= 0x80 >> (2*(i%4)+1);
  3.1260 -      lon *= 2; lon -= 180;
  3.1261 -    } else {
  3.1262 -      lon *= 2; lon += 180;
  3.1263 -    }
  3.1264 -  }
  3.1265 -}
  3.1266 -
  3.1267 -/* convert a circle on earth into a max-depth key to be used in an index */
  3.1268 -static void pgl_circle_to_key(pgl_circle *circle, pgl_keyptr key) {
  3.1269 -  /* handle special case of empty circle */
  3.1270 -  if (circle->radius < 0) {
  3.1271 -    pgl_key_set_empty(key);
  3.1272 -    return;
  3.1273 -  }
  3.1274 -  /* perform same action as for point keys */
  3.1275 -  pgl_point_to_key(&(circle->center), key);
  3.1276 -  /* but overwrite type and node depth to fit area index key */
  3.1277 -  key[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | PGL_AREAKEY_MAXDEPTH;
  3.1278 -  /* check if radius is greater than (or equal to) reference size */
  3.1279 -  /* (treat equal values as greater values for numerical safety) */
  3.1280 -  if (circle->radius >= PGL_AREAKEY_REFOBJSIZE) {
  3.1281 -    /* if yes, set logarithmic size to zero */
  3.1282 -    key[PGL_KEY_OBJSIZE_OFFSET] = 0;
  3.1283 -  } else {
  3.1284 -    /* otherwise, determine logarithmic size iteratively */
  3.1285 -    /* (one step is equivalent to a factor of sqrt(2)) */
  3.1286 -    double reference = PGL_AREAKEY_REFOBJSIZE / M_SQRT2;
  3.1287 -    int objsize = 1;
  3.1288 -    while (objsize < PGL_AREAKEY_MAXOBJSIZE) {
  3.1289 -      /* stop when radius is greater than (or equal to) adjusted reference */
  3.1290 -      /* (treat equal values as greater values for numerical safety) */
  3.1291 -      if (circle->radius >= reference) break;
  3.1292 -      reference /= M_SQRT2;
  3.1293 -      objsize++;
  3.1294 -    }
  3.1295 -    /* set logarithmic size to determined value */
  3.1296 -    key[PGL_KEY_OBJSIZE_OFFSET] = objsize;
  3.1297 -  }
  3.1298 -}
  3.1299 -
  3.1300 -/* check if one key is subkey of another key or vice versa */
  3.1301 -static bool pgl_keys_overlap(pgl_keyptr key1, pgl_keyptr key2) {
  3.1302 -  int i;  /* key bit offset (includes both lat/lon and log. obj. size bits) */
  3.1303 -  /* determine smallest depth */
  3.1304 -  int depth1 = PGL_KEY_NODEDEPTH(key1);
  3.1305 -  int depth2 = PGL_KEY_NODEDEPTH(key2);
  3.1306 -  int depth = (depth1 < depth2) ? depth1 : depth2;
  3.1307 -  /* check if keys are area keys (assuming that both keys have same type) */
  3.1308 -  if (PGL_KEY_IS_AREAKEY(key1)) {
  3.1309 -    int j = 0;  /* bit offset for logarithmic object size bits */
  3.1310 -    int k = 0;  /* bit offset for latitude and longitude */
  3.1311 -    /* fetch logarithmic object size information */
  3.1312 -    int objsize1 = PGL_KEY_OBJSIZE(key1);
  3.1313 -    int objsize2 = PGL_KEY_OBJSIZE(key2);
  3.1314 -    /* handle special cases for empty objects (universal and empty keys) */
  3.1315 -    if (
  3.1316 -      objsize1 == PGL_KEY_OBJSIZE_UNIVERSAL ||
  3.1317 -      objsize2 == PGL_KEY_OBJSIZE_UNIVERSAL
  3.1318 -    ) return true;
  3.1319 -    if (
  3.1320 -      objsize1 == PGL_KEY_OBJSIZE_EMPTY ||
  3.1321 -      objsize2 == PGL_KEY_OBJSIZE_EMPTY
  3.1322 -    ) return objsize1 == objsize2;
  3.1323 -    /* iterate through key bits */
  3.1324 -    for (i=0; i<depth; i++) {
  3.1325 -      /* every second bit is a bit describing the object size */
  3.1326 -      if (i%2 == 0) {
  3.1327 -        /* check if object size bit is different in both keys (objsize1 and
  3.1328 -           objsize2 describe the minimum index when object size bit is set) */
  3.1329 -        if (
  3.1330 -          (objsize1 <= j && objsize2 > j) ||
  3.1331 -          (objsize2 <= j && objsize1 > j)
  3.1332 -        ) {
  3.1333 -          /* bit differs, therefore keys are in separate branches */
  3.1334 -          return false;
  3.1335 -        }
  3.1336 -        /* increase bit counter for object size bits */
  3.1337 -        j++;
  3.1338 -      }
  3.1339 -      /* all other bits describe latitude and longitude */
  3.1340 -      else {
  3.1341 -        /* check if bit differs in both keys */
  3.1342 -        if (PGL_KEY_LATLONBIT_DIFF(key1, key2, k)) {
  3.1343 -          /* bit differs, therefore keys are in separate branches */
  3.1344 -          return false;
  3.1345 -        }
  3.1346 -        /* increase bit counter for latitude/longitude bits */
  3.1347 -        k++;
  3.1348 -      }
  3.1349 -    }
  3.1350 -  }
  3.1351 -  /* if not, keys are point keys */
  3.1352 -  else {
  3.1353 -    /* iterate through key bits */
  3.1354 -    for (i=0; i<depth; i++) {
  3.1355 -      /* check if bit differs in both keys */
  3.1356 -      if (PGL_KEY_LATLONBIT_DIFF(key1, key2, i)) {
  3.1357 -        /* bit differs, therefore keys are in separate branches */
  3.1358 -        return false;
  3.1359 -      }
  3.1360 -    }
  3.1361 -  }
  3.1362 -  /* return true because keys are in the same branch */
  3.1363 -  return true;
  3.1364 -}
  3.1365 -
  3.1366 -/* combine two keys into new key which covers both original keys */
  3.1367 -/* (result stored in first argument) */
  3.1368 -static void pgl_unite_keys(pgl_keyptr dst, pgl_keyptr src) {
  3.1369 -  int i;  /* key bit offset (includes both lat/lon and log. obj. size bits) */
  3.1370 -  /* determine smallest depth */
  3.1371 -  int depth1 = PGL_KEY_NODEDEPTH(dst);
  3.1372 -  int depth2 = PGL_KEY_NODEDEPTH(src);
  3.1373 -  int depth = (depth1 < depth2) ? depth1 : depth2;
  3.1374 -  /* check if keys are area keys (assuming that both keys have same type) */
  3.1375 -  if (PGL_KEY_IS_AREAKEY(dst)) {
  3.1376 -    pgl_areakey dstbuf = { 0, };  /* destination buffer (cleared) */
  3.1377 -    int j = 0;  /* bit offset for logarithmic object size bits */
  3.1378 -    int k = 0;  /* bit offset for latitude and longitude */
  3.1379 -    /* fetch logarithmic object size information */
  3.1380 -    int objsize1 = PGL_KEY_OBJSIZE(dst);
  3.1381 -    int objsize2 = PGL_KEY_OBJSIZE(src);
  3.1382 -    /* handle special cases for empty objects (universal and empty keys) */
  3.1383 -    if (
  3.1384 -      objsize1 > PGL_AREAKEY_MAXOBJSIZE ||
  3.1385 -      objsize2 > PGL_AREAKEY_MAXOBJSIZE
  3.1386 -    ) {
  3.1387 -      if (
  3.1388 -        objsize1 == PGL_KEY_OBJSIZE_EMPTY &&
  3.1389 -        objsize2 == PGL_KEY_OBJSIZE_EMPTY
  3.1390 -      ) pgl_key_set_empty(dst);
  3.1391 -      else pgl_key_set_universal(dst);
  3.1392 -      return;
  3.1393 -    }
  3.1394 -    /* iterate through key bits */
  3.1395 -    for (i=0; i<depth; i++) {
  3.1396 -      /* every second bit is a bit describing the object size */
  3.1397 -      if (i%2 == 0) {
  3.1398 -        /* increase bit counter for object size bits first */
  3.1399 -        /* (handy when setting objsize variable) */
  3.1400 -        j++;
  3.1401 -        /* check if object size bit is set in neither key */
  3.1402 -        if (objsize1 >= j && objsize2 >= j) {
  3.1403 -          /* set objsize in destination buffer to indicate that size bit is
  3.1404 -             unset in destination buffer at the current bit position */
  3.1405 -          dstbuf[PGL_KEY_OBJSIZE_OFFSET] = j;
  3.1406 -        }
  3.1407 -        /* break if object size bit is set in one key only */
  3.1408 -        else if (objsize1 >= j || objsize2 >= j) break;
  3.1409 -      }
  3.1410 -      /* all other bits describe latitude and longitude */
  3.1411 -      else {
  3.1412 -        /* break if bit differs in both keys */
  3.1413 -        if (PGL_KEY_LATLONBIT(dst, k)) {
  3.1414 -          if (!PGL_KEY_LATLONBIT(src, k)) break;
  3.1415 -          /* but set bit in destination buffer if bit is set in both keys */
  3.1416 -          dstbuf[k/8] |= 0x80 >> (k%8);
  3.1417 -        } else if (PGL_KEY_LATLONBIT(src, k)) break;
  3.1418 -        /* increase bit counter for latitude/longitude bits */
  3.1419 -        k++;
  3.1420 -      }
  3.1421 -    }
  3.1422 -    /* set common node depth and type bit (type bit = 1) */
  3.1423 -    dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = PGL_AREAKEY_TYPEMASK | i;
  3.1424 -    /* copy contents of destination buffer to first key */
  3.1425 -    memcpy(dst, dstbuf, sizeof(pgl_areakey));
  3.1426 -  }
  3.1427 -  /* if not, keys are point keys */
  3.1428 -  else {
  3.1429 -    pgl_pointkey dstbuf = { 0, };  /* destination buffer (cleared) */
  3.1430 -    /* iterate through key bits */
  3.1431 -    for (i=0; i<depth; i++) {
  3.1432 -      /* break if bit differs in both keys */
  3.1433 -      if (PGL_KEY_LATLONBIT(dst, i)) {
  3.1434 -        if (!PGL_KEY_LATLONBIT(src, i)) break;
  3.1435 -        /* but set bit in destination buffer if bit is set in both keys */
  3.1436 -        dstbuf[i/8] |= 0x80 >> (i%8);
  3.1437 -      } else if (PGL_KEY_LATLONBIT(src, i)) break;
  3.1438 -    }
  3.1439 -    /* set common node depth (type bit = 0) */
  3.1440 -    dstbuf[PGL_KEY_NODEDEPTH_OFFSET] = i;
  3.1441 -    /* copy contents of destination buffer to first key */
  3.1442 -    memcpy(dst, dstbuf, sizeof(pgl_pointkey));
  3.1443 -  }
  3.1444 -}
  3.1445 -
  3.1446 -/* determine center(!) boundaries and radius estimation of index key */
  3.1447 -static double pgl_key_to_box(pgl_keyptr key, pgl_box *box) {
  3.1448 -  int i;
  3.1449 -  /* determine node depth */
  3.1450 -  int depth = PGL_KEY_NODEDEPTH(key);
  3.1451 -  /* center point of possible result */
  3.1452 -  double lat = 0;
  3.1453 -  double lon = 0;
  3.1454 -  /* maximum distance of real center point from key center */
  3.1455 -  double dlat = 90;
  3.1456 -  double dlon = 180;
  3.1457 -  /* maximum radius of contained objects */
  3.1458 -  double radius = 0;  /* always return zero for point index keys */
  3.1459 -  /* check if key is area key */
  3.1460 -  if (PGL_KEY_IS_AREAKEY(key)) {
  3.1461 -    /* get logarithmic object size */
  3.1462 -    int objsize = PGL_KEY_OBJSIZE(key);
  3.1463 -    /* handle special cases for empty objects (universal and empty keys) */
  3.1464 -    if (objsize == PGL_KEY_OBJSIZE_EMPTY) {
  3.1465 -      pgl_box_set_empty(box);
  3.1466 -      return 0;
  3.1467 -    } else if (objsize == PGL_KEY_OBJSIZE_UNIVERSAL) {
  3.1468 -      box->lat_min = -90;
  3.1469 -      box->lat_max =  90;
  3.1470 -      box->lon_min = -180;
  3.1471 -      box->lon_max =  180;
  3.1472 -      return 0;  /* any value >= 0 would do */
  3.1473 -    }
  3.1474 -    /* calculate maximum possible radius of objects covered by the given key */
  3.1475 -    if (objsize == 0) radius = INFINITY;
  3.1476 -    else {
  3.1477 -      radius = PGL_AREAKEY_REFOBJSIZE;
  3.1478 -      while (--objsize) radius /= M_SQRT2;
  3.1479 -    }
  3.1480 -    /* iterate over latitude and longitude bits in key */
  3.1481 -    /* (every second bit is a latitude or longitude bit) */
  3.1482 -    for (i=0; i<depth/2; i++) {
  3.1483 -      /* check if latitude bit */
  3.1484 -      if (i%2 == 0) {
  3.1485 -        /* cut latitude dimension in half */
  3.1486 -        dlat /= 2;
  3.1487 -        /* increase center latitude if bit is 1, otherwise decrease */
  3.1488 -        if (PGL_KEY_LATLONBIT(key, i)) lat += dlat;
  3.1489 -        else lat -= dlat;
  3.1490 -      }
  3.1491 -      /* otherwise longitude bit */
  3.1492 -      else {
  3.1493 -        /* cut longitude dimension in half */
  3.1494 -        dlon /= 2;
  3.1495 -        /* increase center longitude if bit is 1, otherwise decrease */
  3.1496 -        if (PGL_KEY_LATLONBIT(key, i)) lon += dlon;
  3.1497 -        else lon -= dlon;
  3.1498 -      }
  3.1499 -    }
  3.1500 -  }
  3.1501 -  /* if not, keys are point keys */
  3.1502 -  else {
  3.1503 -    /* iterate over all bits in key */
  3.1504 -    for (i=0; i<depth; i++) {
  3.1505 -      /* check if latitude bit */
  3.1506 -      if (i%2 == 0) {
  3.1507 -        /* cut latitude dimension in half */
  3.1508 -        dlat /= 2;
  3.1509 -        /* increase center latitude if bit is 1, otherwise decrease */
  3.1510 -        if (PGL_KEY_LATLONBIT(key, i)) lat += dlat;
  3.1511 -        else lat -= dlat;
  3.1512 -      }
  3.1513 -      /* otherwise longitude bit */
  3.1514 -      else {
  3.1515 -        /* cut longitude dimension in half */
  3.1516 -        dlon /= 2;
  3.1517 -        /* increase center longitude if bit is 1, otherwise decrease */
  3.1518 -        if (PGL_KEY_LATLONBIT(key, i)) lon += dlon;
  3.1519 -        else lon -= dlon;
  3.1520 -      }
  3.1521 -    }
  3.1522 -  }
  3.1523 -  /* calculate boundaries from center point and remaining dlat and dlon */
  3.1524 -  /* (return values through pointer to box) */
  3.1525 -  box->lat_min = lat - dlat;
  3.1526 -  box->lat_max = lat + dlat;
  3.1527 -  box->lon_min = lon - dlon;
  3.1528 -  box->lon_max = lon + dlon;
  3.1529 -  /* return radius (as a function return value) */
  3.1530 -  return radius;
  3.1531 -}
  3.1532 -
  3.1533 -/* estimator function for distance between point and index key */
  3.1534 -/* always returns a smaller value than actually correct or zero */
  3.1535 -static double pgl_estimate_key_distance(pgl_keyptr key, pgl_point *point) {
  3.1536 -  pgl_box box;  /* center(!) bounding box of area index key */
  3.1537 -  /* calculate center(!) bounding box and maximum radius of objects covered
  3.1538 -     by area index key (radius is zero for point index keys) */
  3.1539 -  double distance = pgl_key_to_box(key, &box);
  3.1540 -  /* calculate estimated distance between bounding box of center point of
  3.1541 -     indexed object and point passed as second argument, then substract maximum
  3.1542 -     radius of objects covered by index key */
  3.1543 -  distance = pgl_estimate_point_box_distance(point, &box) - distance;
  3.1544 -  /* truncate negative results to zero */
  3.1545 -  if (distance <= 0) distance = 0;
  3.1546 -  /* return result */
  3.1547 -  return distance;
  3.1548 -}
  3.1549 -
  3.1550 -
  3.1551 -/*---------------------------------*
  3.1552 - *  helper functions for text I/O  *
  3.1553 - *---------------------------------*/
  3.1554 -
  3.1555 -#define PGL_NUMBUFLEN 64  /* buffer size for number to string conversion */
  3.1556 -
  3.1557 -/* convert floating point number to string (round-trip safe) */
  3.1558 -static void pgl_print_float(char *buf, double flt) {
  3.1559 -  /* check if number is integral */
  3.1560 -  if (trunc(flt) == flt) {
  3.1561 -    /* for integral floats use maximum precision */
  3.1562 -    snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt);
  3.1563 -  } else {
  3.1564 -    /* otherwise check if 15, 16, or 17 digits needed (round-trip safety) */
  3.1565 -    snprintf(buf, PGL_NUMBUFLEN, "%.15g", flt);
  3.1566 -    if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.16g", flt);
  3.1567 -    if (strtod(buf, NULL) != flt) snprintf(buf, PGL_NUMBUFLEN, "%.17g", flt);
  3.1568 -  }
  3.1569 -}
  3.1570 -
  3.1571 -/* convert latitude floating point number (in degrees) to string */
  3.1572 -static void pgl_print_lat(char *buf, double lat) {
  3.1573 -  if (signbit(lat)) {
  3.1574 -    /* treat negative latitudes (including -0) as south */
  3.1575 -    snprintf(buf, PGL_NUMBUFLEN, "S%015.12f", -lat);
  3.1576 -  } else {
  3.1577 -    /* treat positive latitudes (including +0) as north */
  3.1578 -    snprintf(buf, PGL_NUMBUFLEN, "N%015.12f", lat);
  3.1579 -  }
  3.1580 -}
  3.1581 -
  3.1582 -/* convert longitude floating point number (in degrees) to string */
  3.1583 -static void pgl_print_lon(char *buf, double lon) {
  3.1584 -  if (signbit(lon)) {
  3.1585 -    /* treat negative longitudes (including -0) as west */
  3.1586 -    snprintf(buf, PGL_NUMBUFLEN, "W%016.12f", -lon);
  3.1587 -  } else {
  3.1588 -    /* treat positive longitudes (including +0) as east */
  3.1589 -    snprintf(buf, PGL_NUMBUFLEN, "E%016.12f", lon);
  3.1590 -  }
  3.1591 -}
  3.1592 -
  3.1593 -/* bit masks used as return value of pgl_scan() function */
  3.1594 -#define PGL_SCAN_NONE 0      /* no value has been parsed */
  3.1595 -#define PGL_SCAN_LAT (1<<0)  /* latitude has been parsed */
  3.1596 -#define PGL_SCAN_LON (1<<1)  /* longitude has been parsed */
  3.1597 -#define PGL_SCAN_LATLON (PGL_SCAN_LAT | PGL_SCAN_LON)  /* bitwise OR of both */
  3.1598 -
  3.1599 -/* parse a coordinate (can be latitude or longitude) */
  3.1600 -static int pgl_scan(char **str, double *lat, double *lon) {
  3.1601 -  double val;
  3.1602 -  int len;
  3.1603 -  if (
  3.1604 -    sscanf(*str, " N %lf %n", &val, &len) ||
  3.1605 -    sscanf(*str, " n %lf %n", &val, &len)
  3.1606 -  ) {
  3.1607 -    *str += len; *lat = val; return PGL_SCAN_LAT;
  3.1608 -  }
  3.1609 -  if (
  3.1610 -    sscanf(*str, " S %lf %n", &val, &len) ||
  3.1611 -    sscanf(*str, " s %lf %n", &val, &len)
  3.1612 -  ) {
  3.1613 -    *str += len; *lat = -val; return PGL_SCAN_LAT;
  3.1614 -  }
  3.1615 -  if (
  3.1616 -    sscanf(*str, " E %lf %n", &val, &len) ||
  3.1617 -    sscanf(*str, " e %lf %n", &val, &len)
  3.1618 -  ) {
  3.1619 -    *str += len; *lon = val; return PGL_SCAN_LON;
  3.1620 -  }
  3.1621 -  if (
  3.1622 -    sscanf(*str, " W %lf %n", &val, &len) ||
  3.1623 -    sscanf(*str, " w %lf %n", &val, &len)
  3.1624 -  ) {
  3.1625 -    *str += len; *lon = -val; return PGL_SCAN_LON;
  3.1626 -  }
  3.1627 -  return PGL_SCAN_NONE;
  3.1628 -}
  3.1629 -
  3.1630 -
  3.1631 -/*-----------------*
  3.1632 - *  SQL functions  *
  3.1633 - *-----------------*/
  3.1634 -
  3.1635 -/* Note: These function names use "epoint", "ebox", etc. notation here instead
  3.1636 -   of "point", "box", etc. in order to distinguish them from any previously
  3.1637 -   defined functions. */
  3.1638 -
  3.1639 -/* function needed for dummy types and/or not implemented features */
  3.1640 -PG_FUNCTION_INFO_V1(pgl_notimpl);
  3.1641 -Datum pgl_notimpl(PG_FUNCTION_ARGS) {
  3.1642 -  ereport(ERROR, (errmsg("not implemented by pgLatLon")));
  3.1643 -}
  3.1644 -
  3.1645 -/* set point to latitude and longitude (including checks) */
  3.1646 -static void pgl_epoint_set_latlon(pgl_point *point, double lat, double lon) {
  3.1647 -  /* reject infinite or NaN values */
  3.1648 -  if (!isfinite(lat) || !isfinite(lon)) {
  3.1649 -    ereport(ERROR, (
  3.1650 -      errcode(ERRCODE_DATA_EXCEPTION),
  3.1651 -      errmsg("epoint requires finite coordinates")
  3.1652 -    ));
  3.1653 -  }
  3.1654 -  /* check latitude bounds */
  3.1655 -  if (lat < -90) {
  3.1656 -    ereport(WARNING, (errmsg("latitude exceeds south pole")));
  3.1657 -    lat = -90;
  3.1658 -  } else if (lat > 90) {
  3.1659 -    ereport(WARNING, (errmsg("latitude exceeds north pole")));
  3.1660 -    lat = 90;
  3.1661 -  }
  3.1662 -  /* check longitude bounds */
  3.1663 -  if (lon < -180) {
  3.1664 -    ereport(NOTICE, (errmsg("longitude west of 180th meridian normalized")));
  3.1665 -    lon += 360 - trunc(lon / 360) * 360;
  3.1666 -  } else if (lon > 180) {
  3.1667 -    ereport(NOTICE, (errmsg("longitude east of 180th meridian normalized")));
  3.1668 -    lon -= 360 + trunc(lon / 360) * 360;
  3.1669 -  }
  3.1670 -  /* store rounded latitude/longitude values for round-trip safety */
  3.1671 -  point->lat = pgl_round(lat);
  3.1672 -  point->lon = pgl_round(lon);
  3.1673 -}
  3.1674 -
  3.1675 -/* create point ("epoint" in SQL) from latitude and longitude */
  3.1676 -PG_FUNCTION_INFO_V1(pgl_create_epoint);
  3.1677 -Datum pgl_create_epoint(PG_FUNCTION_ARGS) {
  3.1678 -  pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point));
  3.1679 -  pgl_epoint_set_latlon(point, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1));
  3.1680 -  PG_RETURN_POINTER(point);
  3.1681 -}
  3.1682 -
  3.1683 -/* parse point ("epoint" in SQL) */
  3.1684 -/* format: '[NS]<float> [EW]<float>' */
  3.1685 -PG_FUNCTION_INFO_V1(pgl_epoint_in);
  3.1686 -Datum pgl_epoint_in(PG_FUNCTION_ARGS) {
  3.1687 -  char *str = PG_GETARG_CSTRING(0);  /* input string */
  3.1688 -  char *strptr = str;  /* current position within string */
  3.1689 -  int done = 0;        /* bit mask storing if latitude or longitude was read */
  3.1690 -  double lat, lon;     /* parsed values as double precision floats */
  3.1691 -  pgl_point *point;    /* return value (to be palloc'ed) */
  3.1692 -  /* parse two floats (each latitude or longitude) separated by white-space */
  3.1693 -  done |= pgl_scan(&strptr, &lat, &lon);
  3.1694 -  if (strptr != str && isspace(strptr[-1])) {
  3.1695 -    done |= pgl_scan(&strptr, &lat, &lon);
  3.1696 -  }
  3.1697 -  /* require end of string, and latitude and longitude parsed successfully */
  3.1698 -  if (strptr[0] || done != PGL_SCAN_LATLON) {
  3.1699 -    ereport(ERROR, (
  3.1700 -      errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  3.1701 -      errmsg("invalid input syntax for type epoint: \"%s\"", str)
  3.1702 -    ));
  3.1703 -  }
  3.1704 -  /* allocate memory for result */
  3.1705 -  point = (pgl_point *)palloc(sizeof(pgl_point));
  3.1706 -  /* set latitude and longitude (and perform checks) */
  3.1707 -  pgl_epoint_set_latlon(point, lat, lon);
  3.1708 -  /* return result */
  3.1709 -  PG_RETURN_POINTER(point);
  3.1710 -}
  3.1711 -
  3.1712 -/* create box ("ebox" in SQL) that is empty */
  3.1713 -PG_FUNCTION_INFO_V1(pgl_create_empty_ebox);
  3.1714 -Datum pgl_create_empty_ebox(PG_FUNCTION_ARGS) {
  3.1715 -  pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
  3.1716 -  pgl_box_set_empty(box);
  3.1717 -  PG_RETURN_POINTER(box);
  3.1718 -}
  3.1719 -
  3.1720 -/* set box to given boundaries (including checks) */
  3.1721 -static void pgl_ebox_set_boundaries(
  3.1722 -  pgl_box *box,
  3.1723 -  double lat_min, double lat_max, double lon_min, double lon_max
  3.1724 -) {
  3.1725 -  /* if minimum latitude is greater than maximum latitude, return empty box */
  3.1726 -  if (lat_min > lat_max) {
  3.1727 -    pgl_box_set_empty(box);
  3.1728 -    return;
  3.1729 -  }
  3.1730 -  /* otherwise reject infinite or NaN values */
  3.1731 -  if (
  3.1732 -    !isfinite(lat_min) || !isfinite(lat_max) ||
  3.1733 -    !isfinite(lon_min) || !isfinite(lon_max)
  3.1734 -  ) {
  3.1735 -    ereport(ERROR, (
  3.1736 -      errcode(ERRCODE_DATA_EXCEPTION),
  3.1737 -      errmsg("ebox requires finite coordinates")
  3.1738 -    ));
  3.1739 -  }
  3.1740 -  /* check latitude bounds */
  3.1741 -  if (lat_max < -90) {
  3.1742 -    ereport(WARNING, (errmsg("northern latitude exceeds south pole")));
  3.1743 -    lat_max = -90;
  3.1744 -  } else if (lat_max > 90) {
  3.1745 -    ereport(WARNING, (errmsg("northern latitude exceeds north pole")));
  3.1746 -    lat_max = 90;
  3.1747 -  }
  3.1748 -  if (lat_min < -90) {
  3.1749 -    ereport(WARNING, (errmsg("southern latitude exceeds south pole")));
  3.1750 -    lat_min = -90;
  3.1751 -  } else if (lat_min > 90) {
  3.1752 -    ereport(WARNING, (errmsg("southern latitude exceeds north pole")));
  3.1753 -    lat_min = 90;
  3.1754 -  }
  3.1755 -  /* check if all longitudes are included */
  3.1756 -  if (lon_max - lon_min >= 360) {
  3.1757 -    if (lon_max - lon_min > 360) ereport(WARNING, (
  3.1758 -      errmsg("longitude coverage greater than 360 degrees")
  3.1759 -    ));
  3.1760 -    lon_min = -180;
  3.1761 -    lon_max = 180;
  3.1762 -  } else {
  3.1763 -    /* normalize longitude bounds */
  3.1764 -    if      (lon_min < -180) lon_min += 360 - trunc(lon_min / 360) * 360;
  3.1765 -    else if (lon_min >  180) lon_min -= 360 + trunc(lon_min / 360) * 360;
  3.1766 -    if      (lon_max < -180) lon_max += 360 - trunc(lon_max / 360) * 360;
  3.1767 -    else if (lon_max >  180) lon_max -= 360 + trunc(lon_max / 360) * 360;
  3.1768 -  }
  3.1769 -  /* store rounded latitude/longitude values for round-trip safety */
  3.1770 -  box->lat_min = pgl_round(lat_min);
  3.1771 -  box->lat_max = pgl_round(lat_max);
  3.1772 -  box->lon_min = pgl_round(lon_min);
  3.1773 -  box->lon_max = pgl_round(lon_max);
  3.1774 -  /* ensure that rounding does not change orientation */
  3.1775 -  if (lon_min > lon_max && box->lon_min == box->lon_max) {
  3.1776 -    box->lon_min = -180;
  3.1777 -    box->lon_max = 180;
  3.1778 -  }
  3.1779 -}
  3.1780 -
  3.1781 -/* create box ("ebox" in SQL) from min/max latitude and min/max longitude */
  3.1782 -PG_FUNCTION_INFO_V1(pgl_create_ebox);
  3.1783 -Datum pgl_create_ebox(PG_FUNCTION_ARGS) {
  3.1784 -  pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
  3.1785 -  pgl_ebox_set_boundaries(
  3.1786 -    box,
  3.1787 -    PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1),
  3.1788 -    PG_GETARG_FLOAT8(2), PG_GETARG_FLOAT8(3)
  3.1789 -  );
  3.1790 -  PG_RETURN_POINTER(box);
  3.1791 -}
  3.1792 -
  3.1793 -/* create box ("ebox" in SQL) from two points ("epoint"s) */
  3.1794 -/* (can not be used to cover a longitude range of more than 120 degrees) */
  3.1795 -PG_FUNCTION_INFO_V1(pgl_create_ebox_from_epoints);
  3.1796 -Datum pgl_create_ebox_from_epoints(PG_FUNCTION_ARGS) {
  3.1797 -  pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0);
  3.1798 -  pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1);
  3.1799 -  pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
  3.1800 -  double lat_min, lat_max, lon_min, lon_max;
  3.1801 -  double dlon;  /* longitude range (delta longitude) */
  3.1802 -  /* order latitude and longitude boundaries */
  3.1803 -  if (point2->lat < point1->lat) {
  3.1804 -    lat_min = point2->lat;
  3.1805 -    lat_max = point1->lat;
  3.1806 -  } else {
  3.1807 -    lat_min = point1->lat;
  3.1808 -    lat_max = point2->lat;
  3.1809 -  }
  3.1810 -  if (point2->lon < point1->lon) {
  3.1811 -    lon_min = point2->lon;
  3.1812 -    lon_max = point1->lon;
  3.1813 -  } else {
  3.1814 -    lon_min = point1->lon;
  3.1815 -    lon_max = point2->lon;
  3.1816 -  }
  3.1817 -  /* calculate longitude range (round to avoid floating point errors) */
  3.1818 -  dlon = pgl_round(lon_max - lon_min);
  3.1819 -  /* determine east-west direction */
  3.1820 -  if (dlon >= 240) {
  3.1821 -    /* assume that 180th meridian is crossed and swap min/max longitude */
  3.1822 -    double swap = lon_min; lon_min = lon_max; lon_max = swap;
  3.1823 -  } else if (dlon > 120) {
  3.1824 -    /* unclear orientation since delta longitude > 120 */
  3.1825 -    ereport(ERROR, (
  3.1826 -      errcode(ERRCODE_DATA_EXCEPTION),
  3.1827 -      errmsg("can not determine east/west orientation for ebox")
  3.1828 -    ));
  3.1829 -  }
  3.1830 -  /* use boundaries to setup box (and perform checks) */
  3.1831 -  pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max);
  3.1832 -  /* return result */
  3.1833 -  PG_RETURN_POINTER(box);
  3.1834 -}
  3.1835 -
  3.1836 -/* parse box ("ebox" in SQL) */
  3.1837 -/* format: '[NS]<float> [EW]<float> [NS]<float> [EW]<float>'
  3.1838 -       or: '[NS]<float> [NS]<float> [EW]<float> [EW]<float>' */
  3.1839 -PG_FUNCTION_INFO_V1(pgl_ebox_in);
  3.1840 -Datum pgl_ebox_in(PG_FUNCTION_ARGS) {
  3.1841 -  char *str = PG_GETARG_CSTRING(0);  /* input string */
  3.1842 -  char *str_lower;     /* lower case version of input string */
  3.1843 -  char *strptr;        /* current position within string */
  3.1844 -  int valid;           /* number of valid chars */
  3.1845 -  int done;            /* specifies if latitude or longitude was read */
  3.1846 -  double val;          /* temporary variable */
  3.1847 -  int lat_count = 0;   /* count of latitude values parsed */
  3.1848 -  int lon_count = 0;   /* count of longitufde values parsed */
  3.1849 -  double lat_min, lat_max, lon_min, lon_max;  /* see pgl_box struct */
  3.1850 -  pgl_box *box;        /* return value (to be palloc'ed) */
  3.1851 -  /* lowercase input */
  3.1852 -  str_lower = psprintf("%s", str);
  3.1853 -  for (strptr=str_lower; *strptr; strptr++) {
  3.1854 -    if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A';
  3.1855 -  }
  3.1856 -  /* reset reading position to start of (lowercase) string */
  3.1857 -  strptr = str_lower;
  3.1858 -  /* check if empty box */
  3.1859 -  valid = 0;
  3.1860 -  sscanf(strptr, " empty %n", &valid);
  3.1861 -  if (valid && strptr[valid] == 0) {
  3.1862 -    /* allocate and return empty box */
  3.1863 -    box = (pgl_box *)palloc(sizeof(pgl_box));
  3.1864 -    pgl_box_set_empty(box);
  3.1865 -    PG_RETURN_POINTER(box);
  3.1866 -  }
  3.1867 -  /* demand four blocks separated by whitespace */
  3.1868 -  valid = 0;
  3.1869 -  sscanf(strptr, " %*s %*s %*s %*s %n", &valid);
  3.1870 -  /* if four blocks separated by whitespace exist, parse those blocks */
  3.1871 -  if (strptr[valid] == 0) while (strptr[0]) {
  3.1872 -    /* parse either latitude or longitude (whichever found in input string) */
  3.1873 -    done = pgl_scan(&strptr, &val, &val);
  3.1874 -    /* store latitude or longitude in lat_min, lat_max, lon_min, or lon_max */
  3.1875 -    if (done == PGL_SCAN_LAT) {
  3.1876 -      if (!lat_count) lat_min = val; else lat_max = val;
  3.1877 -      lat_count++;
  3.1878 -    } else if (done == PGL_SCAN_LON) {
  3.1879 -      if (!lon_count) lon_min = val; else lon_max = val;
  3.1880 -      lon_count++;
  3.1881 -    } else {
  3.1882 -      break;
  3.1883 -    }
  3.1884 -  }
  3.1885 -  /* require end of string, and two latitude and two longitude values */
  3.1886 -  if (strptr[0] || lat_count != 2 || lon_count != 2) {
  3.1887 -    ereport(ERROR, (
  3.1888 -      errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  3.1889 -      errmsg("invalid input syntax for type ebox: \"%s\"", str)
  3.1890 -    ));
  3.1891 -  }
  3.1892 -  /* free lower case string */
  3.1893 -  pfree(str_lower);
  3.1894 -  /* order boundaries (maximum greater than minimum) */
  3.1895 -  if (lat_min > lat_max) { val = lat_min; lat_min = lat_max; lat_max = val; }
  3.1896 -  if (lon_min > lon_max) { val = lon_min; lon_min = lon_max; lon_max = val; }
  3.1897 -  /* allocate memory for result */
  3.1898 -  box = (pgl_box *)palloc(sizeof(pgl_box));
  3.1899 -  /* set boundaries (and perform checks) */
  3.1900 -  pgl_ebox_set_boundaries(box, lat_min, lat_max, lon_min, lon_max);
  3.1901 -  /* return result */
  3.1902 -  PG_RETURN_POINTER(box);
  3.1903 -}
  3.1904 -
  3.1905 -/* set circle to given latitude, longitude, and radius (including checks) */
  3.1906 -static void pgl_ecircle_set_latlon_radius(
  3.1907 -  pgl_circle *circle, double lat, double lon, double radius
  3.1908 -) {
  3.1909 -  /* set center point (including checks) */
  3.1910 -  pgl_epoint_set_latlon(&(circle->center), lat, lon);
  3.1911 -  /* handle non-positive radius */
  3.1912 -  if (isnan(radius)) {
  3.1913 -    ereport(ERROR, (
  3.1914 -      errcode(ERRCODE_DATA_EXCEPTION),
  3.1915 -      errmsg("invalid radius for ecircle")
  3.1916 -    ));
  3.1917 -  }
  3.1918 -  if (radius == 0) radius = 0;  /* avoids -0 */
  3.1919 -  else if (radius < 0) {
  3.1920 -    if (isfinite(radius)) {
  3.1921 -      ereport(NOTICE, (errmsg("negative radius converted to minus infinity")));
  3.1922 -    }
  3.1923 -    radius = -INFINITY;
  3.1924 -  }
  3.1925 -  /* store radius (round-trip safety is ensured by pgl_print_float) */
  3.1926 -  circle->radius = radius;
  3.1927 -}
  3.1928 -
  3.1929 -/* create circle ("ecircle" in SQL) from latitude, longitude, and radius */
  3.1930 -PG_FUNCTION_INFO_V1(pgl_create_ecircle);
  3.1931 -Datum pgl_create_ecircle(PG_FUNCTION_ARGS) {
  3.1932 -  pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
  3.1933 -  pgl_ecircle_set_latlon_radius(
  3.1934 -    circle, PG_GETARG_FLOAT8(0), PG_GETARG_FLOAT8(1), PG_GETARG_FLOAT8(2)
  3.1935 -  );
  3.1936 -  PG_RETURN_POINTER(circle);
  3.1937 -}
  3.1938 -
  3.1939 -/* create circle ("ecircle" in SQL) from point ("epoint"), and radius */
  3.1940 -PG_FUNCTION_INFO_V1(pgl_create_ecircle_from_epoint);
  3.1941 -Datum pgl_create_ecircle_from_epoint(PG_FUNCTION_ARGS) {
  3.1942 -  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  3.1943 -  double radius = PG_GETARG_FLOAT8(1);
  3.1944 -  pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
  3.1945 -  /* set latitude, longitude, radius (and perform checks) */
  3.1946 -  pgl_ecircle_set_latlon_radius(circle, point->lat, point->lon, radius);
  3.1947 -  /* return result */
  3.1948 -  PG_RETURN_POINTER(circle);
  3.1949 -}
  3.1950 -
  3.1951 -/* parse circle ("ecircle" in SQL) */
  3.1952 -/* format: '[NS]<float> [EW]<float> <float>' */
  3.1953 -PG_FUNCTION_INFO_V1(pgl_ecircle_in);
  3.1954 -Datum pgl_ecircle_in(PG_FUNCTION_ARGS) {
  3.1955 -  char *str = PG_GETARG_CSTRING(0);  /* input string */
  3.1956 -  char *strptr = str;       /* current position within string */
  3.1957 -  double lat, lon, radius;  /* parsed values as double precision flaots */
  3.1958 -  int valid = 0;            /* number of valid chars */
  3.1959 -  int done = 0;             /* stores if latitude and/or longitude was read */
  3.1960 -  pgl_circle *circle;       /* return value (to be palloc'ed) */
  3.1961 -  /* demand three blocks separated by whitespace */
  3.1962 -  sscanf(strptr, " %*s %*s %*s %n", &valid);
  3.1963 -  /* if three blocks separated by whitespace exist, parse those blocks */
  3.1964 -  if (strptr[valid] == 0) {
  3.1965 -    /* parse latitude and longitude */
  3.1966 -    done |= pgl_scan(&strptr, &lat, &lon);
  3.1967 -    done |= pgl_scan(&strptr, &lat, &lon);
  3.1968 -    /* parse radius (while incrementing strptr by number of bytes parsed) */
  3.1969 -    valid = 0;
  3.1970 -    if (sscanf(strptr, " %lf %n", &radius, &valid) == 1) strptr += valid;
  3.1971 -  }
  3.1972 -  /* require end of string and both latitude and longitude being parsed */
  3.1973 -  if (strptr[0] || done != PGL_SCAN_LATLON) {
  3.1974 -    ereport(ERROR, (
  3.1975 -      errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  3.1976 -      errmsg("invalid input syntax for type ecircle: \"%s\"", str)
  3.1977 -    ));
  3.1978 -  }
  3.1979 -  /* allocate memory for result */
  3.1980 -  circle = (pgl_circle *)palloc(sizeof(pgl_circle));
  3.1981 -  /* set latitude, longitude, radius (and perform checks) */
  3.1982 -  pgl_ecircle_set_latlon_radius(circle, lat, lon, radius);
  3.1983 -  /* return result */
  3.1984 -  PG_RETURN_POINTER(circle);
  3.1985 -}
  3.1986 -
  3.1987 -/* parse cluster ("ecluster" in SQL) */
  3.1988 -PG_FUNCTION_INFO_V1(pgl_ecluster_in);
  3.1989 -Datum pgl_ecluster_in(PG_FUNCTION_ARGS) {
  3.1990 -  int i;
  3.1991 -  char *str = PG_GETARG_CSTRING(0);  /* input string */
  3.1992 -  char *str_lower;         /* lower case version of input string */
  3.1993 -  char *strptr;            /* pointer to current reading position of input */
  3.1994 -  int npoints_total = 0;   /* total number of points in cluster */
  3.1995 -  int nentries = 0;        /* total number of entries */
  3.1996 -  pgl_newentry *entries;   /* array of pgl_newentry to create pgl_cluster */
  3.1997 -  int entries_buflen = 4;  /* maximum number of elements in entries array */
  3.1998 -  int valid;               /* number of valid chars processed */
  3.1999 -  double lat, lon;         /* latitude and longitude of parsed point */
  3.2000 -  int entrytype;           /* current entry type */
  3.2001 -  int npoints;             /* number of points in current entry */
  3.2002 -  pgl_point *points;       /* array of pgl_point for pgl_newentry */
  3.2003 -  int points_buflen;       /* maximum number of elements in points array */
  3.2004 -  int done;                /* return value of pgl_scan function */
  3.2005 -  pgl_cluster *cluster;    /* created cluster */
  3.2006 -  /* lowercase input */
  3.2007 -  str_lower = psprintf("%s", str);
  3.2008 -  for (strptr=str_lower; *strptr; strptr++) {
  3.2009 -    if (*strptr >= 'A' && *strptr <= 'Z') *strptr += 'a' - 'A';
  3.2010 -  }
  3.2011 -  /* reset reading position to start of (lowercase) string */
  3.2012 -  strptr = str_lower;
  3.2013 -  /* allocate initial buffer for entries */
  3.2014 -  entries = palloc(entries_buflen * sizeof(pgl_newentry));
  3.2015 -  /* parse until end of string */
  3.2016 -  while (strptr[0]) {
  3.2017 -    /* require previous white-space or closing parenthesis before next token */
  3.2018 -    if (strptr != str_lower && !isspace(strptr[-1]) && strptr[-1] != ')') {
  3.2019 -      goto pgl_ecluster_in_error;
  3.2020 -    }
  3.2021 -    /* ignore token "empty" */
  3.2022 -    valid = 0; sscanf(strptr, " empty %n", &valid);
  3.2023 -    if (valid) { strptr += valid; continue; }
  3.2024 -    /* test for "point" token */
  3.2025 -    valid = 0; sscanf(strptr, " point ( %n", &valid);
  3.2026 -    if (valid) {
  3.2027 -      strptr += valid;
  3.2028 -      entrytype = PGL_ENTRY_POINT;
  3.2029 -      goto pgl_ecluster_in_type_ok;
  3.2030 -    }
  3.2031 -    /* test for "path" token */
  3.2032 -    valid = 0; sscanf(strptr, " path ( %n", &valid);
  3.2033 -    if (valid) {
  3.2034 -      strptr += valid;
  3.2035 -      entrytype = PGL_ENTRY_PATH;
  3.2036 -      goto pgl_ecluster_in_type_ok;
  3.2037 -    }
  3.2038 -    /* test for "outline" token */
  3.2039 -    valid = 0; sscanf(strptr, " outline ( %n", &valid);
  3.2040 -    if (valid) {
  3.2041 -      strptr += valid;
  3.2042 -      entrytype = PGL_ENTRY_OUTLINE;
  3.2043 -      goto pgl_ecluster_in_type_ok;
  3.2044 -    }
  3.2045 -    /* test for "polygon" token */
  3.2046 -    valid = 0; sscanf(strptr, " polygon ( %n", &valid);
  3.2047 -    if (valid) {
  3.2048 -      strptr += valid;
  3.2049 -      entrytype = PGL_ENTRY_POLYGON;
  3.2050 -      goto pgl_ecluster_in_type_ok;
  3.2051 -    }
  3.2052 -    /* error if no valid token found */
  3.2053 -    goto pgl_ecluster_in_error;
  3.2054 -    pgl_ecluster_in_type_ok:
  3.2055 -    /* check if pgl_newentry array needs to grow */
  3.2056 -    if (nentries == entries_buflen) {
  3.2057 -      pgl_newentry *newbuf;
  3.2058 -      entries_buflen *= 2;
  3.2059 -      newbuf = palloc(entries_buflen * sizeof(pgl_newentry));
  3.2060 -      memcpy(newbuf, entries, nentries * sizeof(pgl_newentry));
  3.2061 -      pfree(entries);
  3.2062 -      entries = newbuf;
  3.2063 -    }
  3.2064 -    /* reset number of points for current entry */
  3.2065 -    npoints = 0;
  3.2066 -    /* allocate array for points */
  3.2067 -    points_buflen = 4;
  3.2068 -    points = palloc(points_buflen * sizeof(pgl_point));
  3.2069 -    /* parse until closing parenthesis */
  3.2070 -    while (strptr[0] != ')') {
  3.2071 -      /* error on unexpected end of string */
  3.2072 -      if (strptr[0] == 0) goto pgl_ecluster_in_error;
  3.2073 -      /* mark neither latitude nor longitude as read */
  3.2074 -      done = PGL_SCAN_NONE;
  3.2075 -      /* require white-space before second, third, etc. point */
  3.2076 -      if (npoints != 0 && !isspace(strptr[-1])) goto pgl_ecluster_in_error;
  3.2077 -      /* scan latitude (or longitude) */
  3.2078 -      done |= pgl_scan(&strptr, &lat, &lon);
  3.2079 -      /* require white-space before second coordinate */
  3.2080 -      if (strptr != str && !isspace(strptr[-1])) goto pgl_ecluster_in_error;
  3.2081 -      /* scan longitude (or latitude) */
  3.2082 -      done |= pgl_scan(&strptr, &lat, &lon);
  3.2083 -      /* error unless both latitude and longitude were parsed */
  3.2084 -      if (done != PGL_SCAN_LATLON) goto pgl_ecluster_in_error;
  3.2085 -      /* throw error if number of points is too high */
  3.2086 -      if (npoints_total == PGL_CLUSTER_MAXPOINTS) {
  3.2087 -        ereport(ERROR, (
  3.2088 -          errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  3.2089 -          errmsg(
  3.2090 -            "too many points for ecluster entry (maximum %i)",
  3.2091 -            PGL_CLUSTER_MAXPOINTS
  3.2092 -          )
  3.2093 -        ));
  3.2094 -      }
  3.2095 -      /* check if pgl_point array needs to grow */
  3.2096 -      if (npoints == points_buflen) {
  3.2097 -        pgl_point *newbuf;
  3.2098 -        points_buflen *= 2;
  3.2099 -        newbuf = palloc(points_buflen * sizeof(pgl_point));
  3.2100 -        memcpy(newbuf, points, npoints * sizeof(pgl_point));
  3.2101 -        pfree(points);
  3.2102 -        points = newbuf;
  3.2103 -      }
  3.2104 -      /* append point to pgl_point array (includes checks) */
  3.2105 -      pgl_epoint_set_latlon(&(points[npoints++]), lat, lon);
  3.2106 -      /* increase total number of points */
  3.2107 -      npoints_total++;
  3.2108 -    }
  3.2109 -    /* error if entry has no points */
  3.2110 -    if (!npoints) goto pgl_ecluster_in_error;
  3.2111 -    /* entries with one point are automatically of type "point" */
  3.2112 -    if (npoints == 1) entrytype = PGL_ENTRY_POINT;
  3.2113 -    /* if entries have more than one point */
  3.2114 -    else {
  3.2115 -      /* throw error if entry type is "point" */
  3.2116 -      if (entrytype == PGL_ENTRY_POINT) {
  3.2117 -        ereport(ERROR, (
  3.2118 -          errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  3.2119 -          errmsg("invalid input syntax for type ecluster (point entry with more than one point)")
  3.2120 -        ));
  3.2121 -      }
  3.2122 -      /* coerce outlines and polygons with more than 2 points to be a path */
  3.2123 -      if (npoints == 2) entrytype = PGL_ENTRY_PATH;
  3.2124 -    }
  3.2125 -    /* append entry to pgl_newentry array */
  3.2126 -    entries[nentries].entrytype = entrytype;
  3.2127 -    entries[nentries].npoints = npoints;
  3.2128 -    entries[nentries].points = points;
  3.2129 -    nentries++;
  3.2130 -    /* consume closing parenthesis */
  3.2131 -    strptr++;
  3.2132 -    /* consume white-space */
  3.2133 -    while (isspace(strptr[0])) strptr++;
  3.2134 -  }
  3.2135 -  /* free lower case string */
  3.2136 -  pfree(str_lower);
  3.2137 -  /* create cluster from pgl_newentry array */
  3.2138 -  cluster = pgl_new_cluster(nentries, entries);
  3.2139 -  /* free pgl_newentry array */
  3.2140 -  for (i=0; i<nentries; i++) pfree(entries[i].points);
  3.2141 -  pfree(entries);
  3.2142 -  /* set bounding circle of cluster and check east/west orientation */
  3.2143 -  if (!pgl_finalize_cluster(cluster)) {
  3.2144 -    ereport(ERROR, (
  3.2145 -      errcode(ERRCODE_DATA_EXCEPTION),
  3.2146 -      errmsg("can not determine east/west orientation for ecluster"),
  3.2147 -      errhint("Ensure that each entry has a longitude span of less than 180 degrees.")
  3.2148 -    ));
  3.2149 -  }
  3.2150 -  /* return cluster */
  3.2151 -  PG_RETURN_POINTER(cluster);
  3.2152 -  /* code to throw error */
  3.2153 -  pgl_ecluster_in_error:
  3.2154 -  ereport(ERROR, (
  3.2155 -    errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
  3.2156 -    errmsg("invalid input syntax for type ecluster: \"%s\"", str)
  3.2157 -  ));
  3.2158 -}
  3.2159 -
  3.2160 -/* convert point ("epoint") to string representation */
  3.2161 -PG_FUNCTION_INFO_V1(pgl_epoint_out);
  3.2162 -Datum pgl_epoint_out(PG_FUNCTION_ARGS) {
  3.2163 -  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  3.2164 -  char latstr[PGL_NUMBUFLEN];
  3.2165 -  char lonstr[PGL_NUMBUFLEN];
  3.2166 -  pgl_print_lat(latstr, point->lat);
  3.2167 -  pgl_print_lon(lonstr, point->lon);
  3.2168 -  PG_RETURN_CSTRING(psprintf("%s %s", latstr, lonstr));
  3.2169 -}
  3.2170 -
  3.2171 -/* convert box ("ebox") to string representation */
  3.2172 -PG_FUNCTION_INFO_V1(pgl_ebox_out);
  3.2173 -Datum pgl_ebox_out(PG_FUNCTION_ARGS) {
  3.2174 -  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
  3.2175 -  double lon_min = box->lon_min;
  3.2176 -  double lon_max = box->lon_max;
  3.2177 -  char lat_min_str[PGL_NUMBUFLEN];
  3.2178 -  char lat_max_str[PGL_NUMBUFLEN];
  3.2179 -  char lon_min_str[PGL_NUMBUFLEN];
  3.2180 -  char lon_max_str[PGL_NUMBUFLEN];
  3.2181 -  /* return string "empty" if box is set to be empty */
  3.2182 -  if (box->lat_min > box->lat_max) PG_RETURN_CSTRING("empty");
  3.2183 -  /* use boundaries exceeding W180 or E180 if 180th meridian is enclosed */
  3.2184 -  /* (required since pgl_box_in orders the longitude boundaries) */
  3.2185 -  if (lon_min > lon_max) {
  3.2186 -    if (lon_min + lon_max >= 0) lon_min -= 360;
  3.2187 -    else lon_max += 360;
  3.2188 -  }
  3.2189 -  /* format and return result */
  3.2190 -  pgl_print_lat(lat_min_str, box->lat_min);
  3.2191 -  pgl_print_lat(lat_max_str, box->lat_max);
  3.2192 -  pgl_print_lon(lon_min_str, lon_min);
  3.2193 -  pgl_print_lon(lon_max_str, lon_max);
  3.2194 -  PG_RETURN_CSTRING(psprintf(
  3.2195 -    "%s %s %s %s",
  3.2196 -    lat_min_str, lon_min_str, lat_max_str, lon_max_str
  3.2197 -  ));
  3.2198 -}
  3.2199 -
  3.2200 -/* convert circle ("ecircle") to string representation */
  3.2201 -PG_FUNCTION_INFO_V1(pgl_ecircle_out);
  3.2202 -Datum pgl_ecircle_out(PG_FUNCTION_ARGS) {
  3.2203 -  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
  3.2204 -  char latstr[PGL_NUMBUFLEN];
  3.2205 -  char lonstr[PGL_NUMBUFLEN];
  3.2206 -  char radstr[PGL_NUMBUFLEN];
  3.2207 -  pgl_print_lat(latstr, circle->center.lat);
  3.2208 -  pgl_print_lon(lonstr, circle->center.lon);
  3.2209 -  pgl_print_float(radstr, circle->radius);
  3.2210 -  PG_RETURN_CSTRING(psprintf("%s %s %s", latstr, lonstr, radstr));
  3.2211 -}
  3.2212 -
  3.2213 -/* convert cluster ("ecluster") to string representation */
  3.2214 -PG_FUNCTION_INFO_V1(pgl_ecluster_out);
  3.2215 -Datum pgl_ecluster_out(PG_FUNCTION_ARGS) {
  3.2216 -  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  3.2217 -  char latstr[PGL_NUMBUFLEN];  /* string buffer for latitude */
  3.2218 -  char lonstr[PGL_NUMBUFLEN];  /* string buffer for longitude */
  3.2219 -  char ***strings;     /* array of array of strings */
  3.2220 -  char *string;        /* string of current token */
  3.2221 -  char *res, *resptr;  /* result and pointer to current write position */
  3.2222 -  size_t reslen = 1;   /* length of result (init with 1 for terminator) */
  3.2223 -  int npoints;         /* number of points of current entry */
  3.2224 -  int i, j;            /* i: entry, j: point in entry */
  3.2225 -  /* handle empty clusters */
  3.2226 -  if (cluster->nentries == 0) {
  3.2227 -    /* free detoasted cluster (if copy) */
  3.2228 -    PG_FREE_IF_COPY(cluster, 0);
  3.2229 -    /* return static result */
  3.2230 -    PG_RETURN_CSTRING("empty");
  3.2231 -  }
  3.2232 -  /* allocate array of array of strings */
  3.2233 -  strings = palloc(cluster->nentries * sizeof(char **));
  3.2234 -  /* iterate over all entries in cluster */
  3.2235 -  for (i=0; i<cluster->nentries; i++) {
  3.2236 -    /* get number of points in entry */
  3.2237 -    npoints = cluster->entries[i].npoints;
  3.2238 -    /* allocate array of strings (one string for each point plus two extra) */
  3.2239 -    strings[i] = palloc((2 + npoints) * sizeof(char *));
  3.2240 -    /* determine opening string */
  3.2241 -    switch (cluster->entries[i].entrytype) {
  3.2242 -      case PGL_ENTRY_POINT:   string = (i==0)?"point ("  :" point (";   break;
  3.2243 -      case PGL_ENTRY_PATH:    string = (i==0)?"path ("   :" path (";    break;
  3.2244 -      case PGL_ENTRY_OUTLINE: string = (i==0)?"outline (":" outline ("; break;
  3.2245 -      case PGL_ENTRY_POLYGON: string = (i==0)?"polygon (":" polygon ("; break;
  3.2246 -      default:                string = (i==0)?"unknown"  :" unknown";
  3.2247 -    }
  3.2248 -    /* use opening string as first string in array */
  3.2249 -    strings[i][0] = string;
  3.2250 -    /* update result length (for allocating result string later) */
  3.2251 -    reslen += strlen(string);
  3.2252 -    /* iterate over all points */
  3.2253 -    for (j=0; j<npoints; j++) {
  3.2254 -      /* create string representation of point */
  3.2255 -      pgl_print_lat(latstr, PGL_ENTRY_POINTS(cluster, i)[j].lat);
  3.2256 -      pgl_print_lon(lonstr, PGL_ENTRY_POINTS(cluster, i)[j].lon);
  3.2257 -      string = psprintf((j == 0) ? "%s %s" : " %s %s", latstr, lonstr);
  3.2258 -      /* copy string pointer to string array */
  3.2259 -      strings[i][j+1] = string;
  3.2260 -      /* update result length (for allocating result string later) */
  3.2261 -      reslen += strlen(string);
  3.2262 -    }
  3.2263 -    /* use closing parenthesis as last string in array */
  3.2264 -    strings[i][npoints+1] = ")";
  3.2265 -    /* update result length (for allocating result string later) */
  3.2266 -    reslen++;
  3.2267 -  }
  3.2268 -  /* allocate result string */
  3.2269 -  res = palloc(reslen);
  3.2270 -  /* set write pointer to begin of result string */
  3.2271 -  resptr = res;
  3.2272 -  /* copy strings into result string */
  3.2273 -  for (i=0; i<cluster->nentries; i++) {
  3.2274 -    npoints = cluster->entries[i].npoints;
  3.2275 -    for (j=0; j<npoints+2; j++) {
  3.2276 -      string = strings[i][j];
  3.2277 -      strcpy(resptr, string);
  3.2278 -      resptr += strlen(string);
  3.2279 -      /* free strings allocated by psprintf */
  3.2280 -      if (j != 0 && j != npoints+1) pfree(string);
  3.2281 -    }
  3.2282 -    /* free array of strings */
  3.2283 -    pfree(strings[i]);
  3.2284 -  }
  3.2285 -  /* free array of array of strings */
  3.2286 -  pfree(strings);
  3.2287 -  /* free detoasted cluster (if copy) */
  3.2288 -  PG_FREE_IF_COPY(cluster, 0);
  3.2289 -  /* return result */
  3.2290 -  PG_RETURN_CSTRING(res);
  3.2291 -}
  3.2292 -
  3.2293 -/* binary input function for point ("epoint") */
  3.2294 -PG_FUNCTION_INFO_V1(pgl_epoint_recv);
  3.2295 -Datum pgl_epoint_recv(PG_FUNCTION_ARGS) {
  3.2296 -  StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
  3.2297 -  pgl_point *point = (pgl_point *)palloc(sizeof(pgl_point));
  3.2298 -  point->lat = pq_getmsgfloat8(buf);
  3.2299 -  point->lon = pq_getmsgfloat8(buf);
  3.2300 -  PG_RETURN_POINTER(point);
  3.2301 -}
  3.2302 -
  3.2303 -/* binary input function for box ("ebox") */
  3.2304 -PG_FUNCTION_INFO_V1(pgl_ebox_recv);
  3.2305 -Datum pgl_ebox_recv(PG_FUNCTION_ARGS) {
  3.2306 -  StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
  3.2307 -  pgl_box *box = (pgl_box *)palloc(sizeof(pgl_box));
  3.2308 -  box->lat_min = pq_getmsgfloat8(buf);
  3.2309 -  box->lat_max = pq_getmsgfloat8(buf);
  3.2310 -  box->lon_min = pq_getmsgfloat8(buf);
  3.2311 -  box->lon_max = pq_getmsgfloat8(buf);
  3.2312 -  PG_RETURN_POINTER(box);
  3.2313 -}
  3.2314 -
  3.2315 -/* binary input function for circle ("ecircle") */
  3.2316 -PG_FUNCTION_INFO_V1(pgl_ecircle_recv);
  3.2317 -Datum pgl_ecircle_recv(PG_FUNCTION_ARGS) {
  3.2318 -  StringInfo buf = (StringInfo)PG_GETARG_POINTER(0);
  3.2319 -  pgl_circle *circle = (pgl_circle *)palloc(sizeof(pgl_circle));
  3.2320 -  circle->center.lat = pq_getmsgfloat8(buf);
  3.2321 -  circle->center.lon = pq_getmsgfloat8(buf);
  3.2322 -  circle->radius = pq_getmsgfloat8(buf);
  3.2323 -  PG_RETURN_POINTER(circle);
  3.2324 -}
  3.2325 -
  3.2326 -/* TODO: binary receive function for cluster */
  3.2327 -
  3.2328 -/* binary output function for point ("epoint") */
  3.2329 -PG_FUNCTION_INFO_V1(pgl_epoint_send);
  3.2330 -Datum pgl_epoint_send(PG_FUNCTION_ARGS) {
  3.2331 -  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  3.2332 -  StringInfoData buf;
  3.2333 -  pq_begintypsend(&buf);
  3.2334 -  pq_sendfloat8(&buf, point->lat);
  3.2335 -  pq_sendfloat8(&buf, point->lon);
  3.2336 -  PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
  3.2337 -}
  3.2338 -
  3.2339 -/* binary output function for box ("ebox") */
  3.2340 -PG_FUNCTION_INFO_V1(pgl_ebox_send);
  3.2341 -Datum pgl_ebox_send(PG_FUNCTION_ARGS) {
  3.2342 -  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
  3.2343 -  StringInfoData buf;
  3.2344 -  pq_begintypsend(&buf);
  3.2345 -  pq_sendfloat8(&buf, box->lat_min);
  3.2346 -  pq_sendfloat8(&buf, box->lat_max);
  3.2347 -  pq_sendfloat8(&buf, box->lon_min);
  3.2348 -  pq_sendfloat8(&buf, box->lon_max);
  3.2349 -  PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
  3.2350 -}
  3.2351 -
  3.2352 -/* binary output function for circle ("ecircle") */
  3.2353 -PG_FUNCTION_INFO_V1(pgl_ecircle_send);
  3.2354 -Datum pgl_ecircle_send(PG_FUNCTION_ARGS) {
  3.2355 -  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
  3.2356 -  StringInfoData buf;
  3.2357 -  pq_begintypsend(&buf);
  3.2358 -  pq_sendfloat8(&buf, circle->center.lat);
  3.2359 -  pq_sendfloat8(&buf, circle->center.lon);
  3.2360 -  pq_sendfloat8(&buf, circle->radius);
  3.2361 -  PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
  3.2362 -}
  3.2363 -
  3.2364 -/* TODO: binary send functions for cluster */
  3.2365 -
  3.2366 -/* cast point ("epoint") to box ("ebox") */
  3.2367 -PG_FUNCTION_INFO_V1(pgl_epoint_to_ebox);
  3.2368 -Datum pgl_epoint_to_ebox(PG_FUNCTION_ARGS) {
  3.2369 -  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  3.2370 -  pgl_box *box = palloc(sizeof(pgl_box));
  3.2371 -  box->lat_min = point->lat;
  3.2372 -  box->lat_max = point->lat;
  3.2373 -  box->lon_min = point->lon;
  3.2374 -  box->lon_max = point->lon;
  3.2375 -  PG_RETURN_POINTER(box);
  3.2376 -}
  3.2377 -
  3.2378 -/* cast point ("epoint") to circle ("ecircle") */
  3.2379 -PG_FUNCTION_INFO_V1(pgl_epoint_to_ecircle);
  3.2380 -Datum pgl_epoint_to_ecircle(PG_FUNCTION_ARGS) {
  3.2381 -  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  3.2382 -  pgl_circle *circle = palloc(sizeof(pgl_box));
  3.2383 -  circle->center = *point;
  3.2384 -  circle->radius = 0;
  3.2385 -  PG_RETURN_POINTER(circle);
  3.2386 -}
  3.2387 -
  3.2388 -/* cast point ("epoint") to cluster ("ecluster") */
  3.2389 -PG_FUNCTION_INFO_V1(pgl_epoint_to_ecluster);
  3.2390 -Datum pgl_epoint_to_ecluster(PG_FUNCTION_ARGS) {
  3.2391 -  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  3.2392 -  pgl_newentry entry;
  3.2393 -  pgl_cluster *cluster;
  3.2394 -  entry.entrytype = PGL_ENTRY_POINT;
  3.2395 -  entry.npoints = 1;
  3.2396 -  entry.points = point;
  3.2397 -  cluster = pgl_new_cluster(1, &entry);
  3.2398 -  pgl_finalize_cluster(cluster);  /* NOTE: should not fail */
  3.2399 -  PG_RETURN_POINTER(cluster);
  3.2400 -}
  3.2401 -
  3.2402 -/* cast box ("ebox") to cluster ("ecluster") */
  3.2403 -#define pgl_ebox_to_ecluster_macro(i, a, b) \
  3.2404 -  entries[i].entrytype = PGL_ENTRY_POLYGON; \
  3.2405 -  entries[i].npoints = 4; \
  3.2406 -  entries[i].points = points[i]; \
  3.2407 -  points[i][0].lat = box->lat_min; \
  3.2408 -  points[i][0].lon = (a); \
  3.2409 -  points[i][1].lat = box->lat_min; \
  3.2410 -  points[i][1].lon = (b); \
  3.2411 -  points[i][2].lat = box->lat_max; \
  3.2412 -  points[i][2].lon = (b); \
  3.2413 -  points[i][3].lat = box->lat_max; \
  3.2414 -  points[i][3].lon = (a);
  3.2415 -PG_FUNCTION_INFO_V1(pgl_ebox_to_ecluster);
  3.2416 -Datum pgl_ebox_to_ecluster(PG_FUNCTION_ARGS) {
  3.2417 -  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
  3.2418 -  double lon, dlon;
  3.2419 -  int nentries;
  3.2420 -  pgl_newentry entries[3];
  3.2421 -  pgl_point points[3][4];
  3.2422 -  pgl_cluster *cluster;
  3.2423 -  if (box->lat_min > box->lat_max) {
  3.2424 -    nentries = 0;
  3.2425 -  } else if (box->lon_min > box->lon_max) {
  3.2426 -    if (box->lon_min < 0) {
  3.2427 -      lon = pgl_round((box->lon_min + 180) / 2.0);
  3.2428 -      nentries = 3;
  3.2429 -      pgl_ebox_to_ecluster_macro(0, box->lon_min, lon);
  3.2430 -      pgl_ebox_to_ecluster_macro(1, lon, 180);
  3.2431 -      pgl_ebox_to_ecluster_macro(2, -180, box->lon_max);
  3.2432 -    } else if (box->lon_max > 0) {
  3.2433 -      lon = pgl_round((box->lon_max - 180) / 2.0);
  3.2434 -      nentries = 3;
  3.2435 -      pgl_ebox_to_ecluster_macro(0, box->lon_min, 180);
  3.2436 -      pgl_ebox_to_ecluster_macro(1, -180, lon);
  3.2437 -      pgl_ebox_to_ecluster_macro(2, lon, box->lon_max);
  3.2438 -    } else {
  3.2439 -      nentries = 2;
  3.2440 -      pgl_ebox_to_ecluster_macro(0, box->lon_min, 180);
  3.2441 -      pgl_ebox_to_ecluster_macro(1, -180, box->lon_max);
  3.2442 -    }
  3.2443 -  } else {
  3.2444 -    dlon = pgl_round(box->lon_max - box->lon_min);
  3.2445 -    if (dlon < 180) {
  3.2446 -      nentries = 1;
  3.2447 -      pgl_ebox_to_ecluster_macro(0, box->lon_min, box->lon_max);
  3.2448 -    } else {
  3.2449 -      lon = pgl_round((box->lon_min + box->lon_max) / 2.0);
  3.2450 -      if (
  3.2451 -        pgl_round(lon - box->lon_min) < 180 &&
  3.2452 -        pgl_round(box->lon_max - lon) < 180
  3.2453 -      ) {
  3.2454 -        nentries = 2;
  3.2455 -        pgl_ebox_to_ecluster_macro(0, box->lon_min, lon);
  3.2456 -        pgl_ebox_to_ecluster_macro(1, lon, box->lon_max);
  3.2457 -      } else {
  3.2458 -        nentries = 3;
  3.2459 -        pgl_ebox_to_ecluster_macro(0, box->lon_min, -60);
  3.2460 -        pgl_ebox_to_ecluster_macro(1, -60, 60);
  3.2461 -        pgl_ebox_to_ecluster_macro(2, 60, box->lon_max);
  3.2462 -      }
  3.2463 -    }
  3.2464 -  }
  3.2465 -  cluster = pgl_new_cluster(nentries, entries);
  3.2466 -  pgl_finalize_cluster(cluster);  /* NOTE: should not fail */
  3.2467 -  PG_RETURN_POINTER(cluster);
  3.2468 -}
  3.2469 -
  3.2470 -/* extract latitude from point ("epoint") */
  3.2471 -PG_FUNCTION_INFO_V1(pgl_epoint_lat);
  3.2472 -Datum pgl_epoint_lat(PG_FUNCTION_ARGS) {
  3.2473 -  PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lat);
  3.2474 -}
  3.2475 -
  3.2476 -/* extract longitude from point ("epoint") */
  3.2477 -PG_FUNCTION_INFO_V1(pgl_epoint_lon);
  3.2478 -Datum pgl_epoint_lon(PG_FUNCTION_ARGS) {
  3.2479 -  PG_RETURN_FLOAT8(((pgl_point *)PG_GETARG_POINTER(0))->lon);
  3.2480 -}
  3.2481 -
  3.2482 -/* extract minimum latitude from box ("ebox") */
  3.2483 -PG_FUNCTION_INFO_V1(pgl_ebox_lat_min);
  3.2484 -Datum pgl_ebox_lat_min(PG_FUNCTION_ARGS) {
  3.2485 -  PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_min);
  3.2486 -}
  3.2487 -
  3.2488 -/* extract maximum latitude from box ("ebox") */
  3.2489 -PG_FUNCTION_INFO_V1(pgl_ebox_lat_max);
  3.2490 -Datum pgl_ebox_lat_max(PG_FUNCTION_ARGS) {
  3.2491 -  PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lat_max);
  3.2492 -}
  3.2493 -
  3.2494 -/* extract minimum longitude from box ("ebox") */
  3.2495 -PG_FUNCTION_INFO_V1(pgl_ebox_lon_min);
  3.2496 -Datum pgl_ebox_lon_min(PG_FUNCTION_ARGS) {
  3.2497 -  PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_min);
  3.2498 -}
  3.2499 -
  3.2500 -/* extract maximum longitude from box ("ebox") */
  3.2501 -PG_FUNCTION_INFO_V1(pgl_ebox_lon_max);
  3.2502 -Datum pgl_ebox_lon_max(PG_FUNCTION_ARGS) {
  3.2503 -  PG_RETURN_FLOAT8(((pgl_box *)PG_GETARG_POINTER(0))->lon_max);
  3.2504 -}
  3.2505 -
  3.2506 -/* extract center point from circle ("ecircle") */
  3.2507 -PG_FUNCTION_INFO_V1(pgl_ecircle_center);
  3.2508 -Datum pgl_ecircle_center(PG_FUNCTION_ARGS) {
  3.2509 -  PG_RETURN_POINTER(&(((pgl_circle *)PG_GETARG_POINTER(0))->center));
  3.2510 -}
  3.2511 -
  3.2512 -/* extract radius from circle ("ecircle") */
  3.2513 -PG_FUNCTION_INFO_V1(pgl_ecircle_radius);
  3.2514 -Datum pgl_ecircle_radius(PG_FUNCTION_ARGS) {
  3.2515 -  PG_RETURN_FLOAT8(((pgl_circle *)PG_GETARG_POINTER(0))->radius);
  3.2516 -}
  3.2517 -
  3.2518 -/* check if point is inside box (overlap operator "&&") in SQL */
  3.2519 -PG_FUNCTION_INFO_V1(pgl_epoint_ebox_overlap);
  3.2520 -Datum pgl_epoint_ebox_overlap(PG_FUNCTION_ARGS) {
  3.2521 -  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  3.2522 -  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(1);
  3.2523 -  PG_RETURN_BOOL(pgl_point_in_box(point, box));
  3.2524 -}
  3.2525 -
  3.2526 -/* check if point is inside circle (overlap operator "&&") in SQL */
  3.2527 -PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_overlap);
  3.2528 -Datum pgl_epoint_ecircle_overlap(PG_FUNCTION_ARGS) {
  3.2529 -  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  3.2530 -  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
  3.2531 -  PG_RETURN_BOOL(
  3.2532 -    pgl_distance(
  3.2533 -      point->lat, point->lon,
  3.2534 -      circle->center.lat, circle->center.lon
  3.2535 -    ) <= circle->radius
  3.2536 -  );
  3.2537 -}
  3.2538 -
  3.2539 -/* check if point is inside cluster (overlap operator "&&") in SQL */
  3.2540 -PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_overlap);
  3.2541 -Datum pgl_epoint_ecluster_overlap(PG_FUNCTION_ARGS) {
  3.2542 -  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  3.2543 -  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  3.2544 -  bool retval;
  3.2545 -  /* points outside bounding circle are always assumed to be non-overlapping
  3.2546 -     (necessary for consistent table and index scans) */
  3.2547 -  if (
  3.2548 -    pgl_distance(
  3.2549 -      point->lat, point->lon,
  3.2550 -      cluster->bounding.center.lat, cluster->bounding.center.lon
  3.2551 -    ) > cluster->bounding.radius
  3.2552 -  ) retval = false;
  3.2553 -  else retval = pgl_point_in_cluster(point, cluster, false);
  3.2554 -  PG_FREE_IF_COPY(cluster, 1);
  3.2555 -  PG_RETURN_BOOL(retval);
  3.2556 -}
  3.2557 -
  3.2558 -/* check if point may be inside cluster (lossy overl. operator "&&+") in SQL */
  3.2559 -PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_may_overlap);
  3.2560 -Datum pgl_epoint_ecluster_may_overlap(PG_FUNCTION_ARGS) {
  3.2561 -  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  3.2562 -  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  3.2563 -  bool retval = pgl_distance(
  3.2564 -    point->lat, point->lon,
  3.2565 -    cluster->bounding.center.lat, cluster->bounding.center.lon
  3.2566 -  ) <= cluster->bounding.radius;
  3.2567 -  PG_FREE_IF_COPY(cluster, 1);
  3.2568 -  PG_RETURN_BOOL(retval);
  3.2569 -}
  3.2570 -
  3.2571 -/* check if two boxes overlap (overlap operator "&&") in SQL */
  3.2572 -PG_FUNCTION_INFO_V1(pgl_ebox_overlap);
  3.2573 -Datum pgl_ebox_overlap(PG_FUNCTION_ARGS) {
  3.2574 -  pgl_box *box1 = (pgl_box *)PG_GETARG_POINTER(0);
  3.2575 -  pgl_box *box2 = (pgl_box *)PG_GETARG_POINTER(1);
  3.2576 -  PG_RETURN_BOOL(pgl_boxes_overlap(box1, box2));
  3.2577 -}
  3.2578 -
  3.2579 -/* check if box and circle may overlap (lossy overl. operator "&&+") in SQL */
  3.2580 -PG_FUNCTION_INFO_V1(pgl_ebox_ecircle_may_overlap);
  3.2581 -Datum pgl_ebox_ecircle_may_overlap(PG_FUNCTION_ARGS) {
  3.2582 -  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
  3.2583 -  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
  3.2584 -  PG_RETURN_BOOL(
  3.2585 -    pgl_estimate_point_box_distance(&circle->center, box) <= circle->radius
  3.2586 -  );
  3.2587 -}
  3.2588 -
  3.2589 -/* check if box and cluster may overlap (lossy overl. operator "&&+") in SQL */
  3.2590 -PG_FUNCTION_INFO_V1(pgl_ebox_ecluster_may_overlap);
  3.2591 -Datum pgl_ebox_ecluster_may_overlap(PG_FUNCTION_ARGS) {
  3.2592 -  pgl_box *box = (pgl_box *)PG_GETARG_POINTER(0);
  3.2593 -  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  3.2594 -  bool retval = pgl_estimate_point_box_distance(
  3.2595 -    &cluster->bounding.center,
  3.2596 -    box
  3.2597 -  ) <= cluster->bounding.radius;
  3.2598 -  PG_FREE_IF_COPY(cluster, 1);
  3.2599 -  PG_RETURN_BOOL(retval);
  3.2600 -}
  3.2601 -
  3.2602 -/* check if two circles overlap (overlap operator "&&") in SQL */
  3.2603 -PG_FUNCTION_INFO_V1(pgl_ecircle_overlap);
  3.2604 -Datum pgl_ecircle_overlap(PG_FUNCTION_ARGS) {
  3.2605 -  pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0);
  3.2606 -  pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1);
  3.2607 -  PG_RETURN_BOOL(
  3.2608 -    pgl_distance(
  3.2609 -      circle1->center.lat, circle1->center.lon,
  3.2610 -      circle2->center.lat, circle2->center.lon
  3.2611 -    ) <= circle1->radius + circle2->radius
  3.2612 -  );
  3.2613 -}
  3.2614 -
  3.2615 -/* check if circle and cluster overlap (overlap operator "&&") in SQL */
  3.2616 -PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_overlap);
  3.2617 -Datum pgl_ecircle_ecluster_overlap(PG_FUNCTION_ARGS) {
  3.2618 -  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
  3.2619 -  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  3.2620 -  bool retval = (
  3.2621 -    pgl_point_cluster_distance(&(circle->center), cluster) <= circle->radius
  3.2622 -  );
  3.2623 -  PG_FREE_IF_COPY(cluster, 1);
  3.2624 -  PG_RETURN_BOOL(retval);
  3.2625 -}
  3.2626 -
  3.2627 -/* check if circle and cluster may overlap (l. ov. operator "&&+") in SQL */
  3.2628 -PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_may_overlap);
  3.2629 -Datum pgl_ecircle_ecluster_may_overlap(PG_FUNCTION_ARGS) {
  3.2630 -  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
  3.2631 -  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  3.2632 -  bool retval = pgl_distance(
  3.2633 -    circle->center.lat, circle->center.lon,
  3.2634 -    cluster->bounding.center.lat, cluster->bounding.center.lon
  3.2635 -  ) <= circle->radius + cluster->bounding.radius;
  3.2636 -  PG_FREE_IF_COPY(cluster, 1);
  3.2637 -  PG_RETURN_BOOL(retval);
  3.2638 -}
  3.2639 -
  3.2640 -/* check if two clusters overlap (overlap operator "&&") in SQL */
  3.2641 -PG_FUNCTION_INFO_V1(pgl_ecluster_overlap);
  3.2642 -Datum pgl_ecluster_overlap(PG_FUNCTION_ARGS) {
  3.2643 -  pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  3.2644 -  pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  3.2645 -  bool retval;
  3.2646 -  /* clusters with non-touching bounding circles are always assumed to be
  3.2647 -     non-overlapping (improves performance and is necessary for consistent
  3.2648 -     table and index scans) */
  3.2649 -  if (
  3.2650 -    pgl_distance(
  3.2651 -      cluster1->bounding.center.lat, cluster1->bounding.center.lon,
  3.2652 -      cluster2->bounding.center.lat, cluster2->bounding.center.lon
  3.2653 -    ) > cluster1->bounding.radius + cluster2->bounding.radius
  3.2654 -  ) retval = false;
  3.2655 -  else retval = pgl_clusters_overlap(cluster1, cluster2);
  3.2656 -  PG_FREE_IF_COPY(cluster1, 0);
  3.2657 -  PG_FREE_IF_COPY(cluster2, 1);
  3.2658 -  PG_RETURN_BOOL(retval);
  3.2659 -}
  3.2660 -
  3.2661 -/* check if two clusters may overlap (lossy overlap operator "&&+") in SQL */
  3.2662 -PG_FUNCTION_INFO_V1(pgl_ecluster_may_overlap);
  3.2663 -Datum pgl_ecluster_may_overlap(PG_FUNCTION_ARGS) {
  3.2664 -  pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  3.2665 -  pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  3.2666 -  bool retval = pgl_distance(
  3.2667 -    cluster1->bounding.center.lat, cluster1->bounding.center.lon,
  3.2668 -    cluster2->bounding.center.lat, cluster2->bounding.center.lon
  3.2669 -  ) <= cluster1->bounding.radius + cluster2->bounding.radius;
  3.2670 -  PG_FREE_IF_COPY(cluster1, 0);
  3.2671 -  PG_FREE_IF_COPY(cluster2, 1);
  3.2672 -  PG_RETURN_BOOL(retval);
  3.2673 -}
  3.2674 -
  3.2675 -/* check if second cluster is in first cluster (cont. operator "@>) in SQL */
  3.2676 -PG_FUNCTION_INFO_V1(pgl_ecluster_contains);
  3.2677 -Datum pgl_ecluster_contains(PG_FUNCTION_ARGS) {
  3.2678 -  pgl_cluster *outer = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  3.2679 -  pgl_cluster *inner = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  3.2680 -  bool retval;
  3.2681 -  /* clusters with non-touching bounding circles are always assumed to be
  3.2682 -     non-overlapping (improves performance and is necessary for consistent
  3.2683 -     table and index scans) */
  3.2684 -  if (
  3.2685 -    pgl_distance(
  3.2686 -      outer->bounding.center.lat, outer->bounding.center.lon,
  3.2687 -      inner->bounding.center.lat, inner->bounding.center.lon
  3.2688 -    ) > outer->bounding.radius + inner->bounding.radius
  3.2689 -  ) retval = false;
  3.2690 -  else retval = pgl_cluster_in_cluster(outer, inner);
  3.2691 -  PG_FREE_IF_COPY(outer, 0);
  3.2692 -  PG_FREE_IF_COPY(inner, 1);
  3.2693 -  PG_RETURN_BOOL(retval);
  3.2694 -}
  3.2695 -
  3.2696 -/* calculate distance between two points ("<->" operator) in SQL */
  3.2697 -PG_FUNCTION_INFO_V1(pgl_epoint_distance);
  3.2698 -Datum pgl_epoint_distance(PG_FUNCTION_ARGS) {
  3.2699 -  pgl_point *point1 = (pgl_point *)PG_GETARG_POINTER(0);
  3.2700 -  pgl_point *point2 = (pgl_point *)PG_GETARG_POINTER(1);
  3.2701 -  PG_RETURN_FLOAT8(pgl_distance(
  3.2702 -    point1->lat, point1->lon, point2->lat, point2->lon
  3.2703 -  ));
  3.2704 -}
  3.2705 -
  3.2706 -/* calculate point to circle distance ("<->" operator) in SQL */
  3.2707 -PG_FUNCTION_INFO_V1(pgl_epoint_ecircle_distance);
  3.2708 -Datum pgl_epoint_ecircle_distance(PG_FUNCTION_ARGS) {
  3.2709 -  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  3.2710 -  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(1);
  3.2711 -  double distance = pgl_distance(
  3.2712 -    point->lat, point->lon, circle->center.lat, circle->center.lon
  3.2713 -  ) - circle->radius;
  3.2714 -  PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
  3.2715 -}
  3.2716 -
  3.2717 -/* calculate point to cluster distance ("<->" operator) in SQL */
  3.2718 -PG_FUNCTION_INFO_V1(pgl_epoint_ecluster_distance);
  3.2719 -Datum pgl_epoint_ecluster_distance(PG_FUNCTION_ARGS) {
  3.2720 -  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(0);
  3.2721 -  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  3.2722 -  double distance = pgl_point_cluster_distance(point, cluster);
  3.2723 -  PG_FREE_IF_COPY(cluster, 1);
  3.2724 -  PG_RETURN_FLOAT8(distance);
  3.2725 -}
  3.2726 -
  3.2727 -/* calculate distance between two circles ("<->" operator) in SQL */
  3.2728 -PG_FUNCTION_INFO_V1(pgl_ecircle_distance);
  3.2729 -Datum pgl_ecircle_distance(PG_FUNCTION_ARGS) {
  3.2730 -  pgl_circle *circle1 = (pgl_circle *)PG_GETARG_POINTER(0);
  3.2731 -  pgl_circle *circle2 = (pgl_circle *)PG_GETARG_POINTER(1);
  3.2732 -  double distance = pgl_distance(
  3.2733 -    circle1->center.lat, circle1->center.lon,
  3.2734 -    circle2->center.lat, circle2->center.lon
  3.2735 -  ) - (circle1->radius + circle2->radius);
  3.2736 -  PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
  3.2737 -}
  3.2738 -
  3.2739 -/* calculate circle to cluster distance ("<->" operator) in SQL */
  3.2740 -PG_FUNCTION_INFO_V1(pgl_ecircle_ecluster_distance);
  3.2741 -Datum pgl_ecircle_ecluster_distance(PG_FUNCTION_ARGS) {
  3.2742 -  pgl_circle *circle = (pgl_circle *)PG_GETARG_POINTER(0);
  3.2743 -  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  3.2744 -  double distance = (
  3.2745 -    pgl_point_cluster_distance(&(circle->center), cluster) - circle->radius
  3.2746 -  );
  3.2747 -  PG_FREE_IF_COPY(cluster, 1);
  3.2748 -  PG_RETURN_FLOAT8((distance <= 0) ? 0 : distance);
  3.2749 -}
  3.2750 -
  3.2751 -/* calculate distance between two clusters ("<->" operator) in SQL */
  3.2752 -PG_FUNCTION_INFO_V1(pgl_ecluster_distance);
  3.2753 -Datum pgl_ecluster_distance(PG_FUNCTION_ARGS) {
  3.2754 -  pgl_cluster *cluster1 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  3.2755 -  pgl_cluster *cluster2 = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  3.2756 -  double retval = pgl_cluster_distance(cluster1, cluster2);
  3.2757 -  PG_FREE_IF_COPY(cluster1, 0);
  3.2758 -  PG_FREE_IF_COPY(cluster2, 1);
  3.2759 -  PG_RETURN_FLOAT8(retval);
  3.2760 -}
  3.2761 -
  3.2762 -PG_FUNCTION_INFO_V1(pgl_ecluster_monte_carlo_area);
  3.2763 -Datum pgl_ecluster_monte_carlo_area(PG_FUNCTION_ARGS) {
  3.2764 -  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  3.2765 -  int samples = PG_GETARG_INT32(1);
  3.2766 -  double retval = pgl_monte_carlo_area(cluster, samples);
  3.2767 -  PG_FREE_IF_COPY(cluster, 0);
  3.2768 -  PG_RETURN_FLOAT8(retval);
  3.2769 -}
  3.2770 -
  3.2771 -PG_FUNCTION_INFO_V1(pgl_ecluster_epoint_fair_distance);
  3.2772 -Datum pgl_ecluster_epoint_fair_distance(PG_FUNCTION_ARGS) {
  3.2773 -  pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
  3.2774 -  pgl_point *point = (pgl_point *)PG_GETARG_POINTER(1);
  3.2775 -  int samples = PG_GETARG_INT32(2);
  3.2776 -  double retval = pgl_fair_distance(point, cluster, samples);
  3.2777 -  PG_FREE_IF_COPY(cluster, 0);
  3.2778 -  PG_RETURN_FLOAT8(retval);
  3.2779 -}
  3.2780 -
  3.2781 -
  3.2782 -/*-----------------------------------------------------------*
  3.2783 - *  B-tree comparison operators and index support functions  *
  3.2784 - *-----------------------------------------------------------*/
  3.2785 -
  3.2786 -/* macro for a B-tree operator (without detoasting) */
  3.2787 -#define PGL_BTREE_OPER(func, type, cmpfunc, oper) \
  3.2788 -  PG_FUNCTION_INFO_V1(func); \
  3.2789 -  Datum func(PG_FUNCTION_ARGS) { \
  3.2790 -    type *a = (type *)PG_GETARG_POINTER(0); \
  3.2791 -    type *b = (type *)PG_GETARG_POINTER(1); \
  3.2792 -    PG_RETURN_BOOL(cmpfunc(a, b) oper 0); \
  3.2793 -  }
  3.2794 -
  3.2795 -/* macro for a B-tree comparison function (without detoasting) */
  3.2796 -#define PGL_BTREE_CMP(func, type, cmpfunc) \
  3.2797 -  PG_FUNCTION_INFO_V1(func); \
  3.2798 -  Datum func(PG_FUNCTION_ARGS) { \
  3.2799 -    type *a = (type *)PG_GETARG_POINTER(0); \
  3.2800 -    type *b = (type *)PG_GETARG_POINTER(1); \
  3.2801 -    PG_RETURN_INT32(cmpfunc(a, b)); \
  3.2802 -  }
  3.2803 -
  3.2804 -/* macro for a B-tree operator (with detoasting) */
  3.2805 -#define PGL_BTREE_OPER_DETOAST(func, type, cmpfunc, oper) \
  3.2806 -  PG_FUNCTION_INFO_V1(func); \
  3.2807 -  Datum func(PG_FUNCTION_ARGS) { \
  3.2808 -    bool res; \
  3.2809 -    type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \
  3.2810 -    type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \
  3.2811 -    res = cmpfunc(a, b) oper 0; \
  3.2812 -    PG_FREE_IF_COPY(a, 0); \
  3.2813 -    PG_FREE_IF_COPY(b, 1); \
  3.2814 -    PG_RETURN_BOOL(res); \
  3.2815 -  }
  3.2816 -
  3.2817 -/* macro for a B-tree comparison function (with detoasting) */
  3.2818 -#define PGL_BTREE_CMP_DETOAST(func, type, cmpfunc) \
  3.2819 -  PG_FUNCTION_INFO_V1(func); \
  3.2820 -  Datum func(PG_FUNCTION_ARGS) { \
  3.2821 -    int32_t res; \
  3.2822 -    type *a = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); \
  3.2823 -    type *b = (type *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); \
  3.2824 -    res = cmpfunc(a, b); \
  3.2825 -    PG_FREE_IF_COPY(a, 0); \
  3.2826 -    PG_FREE_IF_COPY(b, 1); \
  3.2827 -    PG_RETURN_INT32(res); \
  3.2828 -  }
  3.2829 -
  3.2830 -/* B-tree operators and comparison function for point */
  3.2831 -PGL_BTREE_OPER(pgl_btree_epoint_lt, pgl_point, pgl_point_cmp, <)
  3.2832 -PGL_BTREE_OPER(pgl_btree_epoint_le, pgl_point, pgl_point_cmp, <=)
  3.2833 -PGL_BTREE_OPER(pgl_btree_epoint_eq, pgl_point, pgl_point_cmp, ==)
  3.2834 -PGL_BTREE_OPER(pgl_btree_epoint_ne, pgl_point, pgl_point_cmp, !=)
  3.2835 -PGL_BTREE_OPER(pgl_btree_epoint_ge, pgl_point, pgl_point_cmp, >=)
  3.2836 -PGL_BTREE_OPER(pgl_btree_epoint_gt, pgl_point, pgl_point_cmp, >)
  3.2837 -PGL_BTREE_CMP(pgl_btree_epoint_cmp, pgl_point, pgl_point_cmp)
  3.2838 -
  3.2839 -/* B-tree operators and comparison function for box */
  3.2840 -PGL_BTREE_OPER(pgl_btree_ebox_lt, pgl_box, pgl_box_cmp, <)
  3.2841 -PGL_BTREE_OPER(pgl_btree_ebox_le, pgl_box, pgl_box_cmp, <=)
  3.2842 -PGL_BTREE_OPER(pgl_btree_ebox_eq, pgl_box, pgl_box_cmp, ==)
  3.2843 -PGL_BTREE_OPER(pgl_btree_ebox_ne, pgl_box, pgl_box_cmp, !=)
  3.2844 -PGL_BTREE_OPER(pgl_btree_ebox_ge, pgl_box, pgl_box_cmp, >=)
  3.2845 -PGL_BTREE_OPER(pgl_btree_ebox_gt, pgl_box, pgl_box_cmp, >)
  3.2846 -PGL_BTREE_CMP(pgl_btree_ebox_cmp, pgl_box, pgl_box_cmp)
  3.2847 -
  3.2848 -/* B-tree operators and comparison function for circle */
  3.2849 -PGL_BTREE_OPER(pgl_btree_ecircle_lt, pgl_circle, pgl_circle_cmp, <)
  3.2850 -PGL_BTREE_OPER(pgl_btree_ecircle_le, pgl_circle, pgl_circle_cmp, <=)
  3.2851 -PGL_BTREE_OPER(pgl_btree_ecircle_eq, pgl_circle, pgl_circle_cmp, ==)
  3.2852 -PGL_BTREE_OPER(pgl_btree_ecircle_ne, pgl_circle, pgl_circle_cmp, !=)
  3.2853 -PGL_BTREE_OPER(pgl_btree_ecircle_ge, pgl_circle, pgl_circle_cmp, >=)
  3.2854 -PGL_BTREE_OPER(pgl_btree_ecircle_gt, pgl_circle, pgl_circle_cmp, >)
  3.2855 -PGL_BTREE_CMP(pgl_btree_ecircle_cmp, pgl_circle, pgl_circle_cmp)
  3.2856 -
  3.2857 -
  3.2858 -/*--------------------------------*
  3.2859 - *  GiST index support functions  *
  3.2860 - *--------------------------------*/
  3.2861 -
  3.2862 -/* GiST "consistent" support function */
  3.2863 -PG_FUNCTION_INFO_V1(pgl_gist_consistent);
  3.2864 -Datum pgl_gist_consistent(PG_FUNCTION_ARGS) {
  3.2865 -  GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  3.2866 -  pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key);
  3.2867 -  StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2);
  3.2868 -  bool *recheck = (bool *)PG_GETARG_POINTER(4);
  3.2869 -  /* demand recheck because index and query methods are lossy */
  3.2870 -  *recheck = true;
  3.2871 -  /* strategy number aliases for different operators using the same strategy */
  3.2872 -  strategy %= 100;
  3.2873 -  /* strategy number 11: equality of two points */
  3.2874 -  if (strategy == 11) {
  3.2875 -    /* query datum is another point */
  3.2876 -    pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
  3.2877 -    /* convert other point to key */
  3.2878 -    pgl_pointkey querykey;
  3.2879 -    pgl_point_to_key(query, querykey);
  3.2880 -    /* return true if both keys overlap */
  3.2881 -    PG_RETURN_BOOL(pgl_keys_overlap(key, querykey));
  3.2882 -  }
  3.2883 -  /* strategy number 13: equality of two circles */
  3.2884 -  if (strategy == 13) {
  3.2885 -    /* query datum is another circle */
  3.2886 -    pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
  3.2887 -    /* convert other circle to key */
  3.2888 -    pgl_areakey querykey;
  3.2889 -    pgl_circle_to_key(query, querykey);
  3.2890 -    /* return true if both keys overlap */
  3.2891 -    PG_RETURN_BOOL(pgl_keys_overlap(key, querykey));
  3.2892 -  }
  3.2893 -  /* for all remaining strategies, keys on empty objects produce no match */
  3.2894 -  /* (check necessary because query radius may be infinite) */
  3.2895 -  if (PGL_KEY_IS_EMPTY(key)) PG_RETURN_BOOL(false);
  3.2896 -  /* strategy number 21: overlapping with point */
  3.2897 -  if (strategy == 21) {
  3.2898 -    /* query datum is a point */
  3.2899 -    pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
  3.2900 -    /* return true if estimated distance (allowed to be smaller than real
  3.2901 -       distance) between index key and point is zero */
  3.2902 -    PG_RETURN_BOOL(pgl_estimate_key_distance(key, query) == 0);
  3.2903 -  }
  3.2904 -  /* strategy number 22: (point) overlapping with box */
  3.2905 -  if (strategy == 22) {
  3.2906 -    /* query datum is a box */
  3.2907 -    pgl_box *query = (pgl_box *)PG_GETARG_POINTER(1);
  3.2908 -    /* determine bounding box of indexed key */
  3.2909 -    pgl_box keybox;
  3.2910 -    pgl_key_to_box(key, &keybox);
  3.2911 -    /* return true if query box overlaps with bounding box of indexed key */
  3.2912 -    PG_RETURN_BOOL(pgl_boxes_overlap(query, &keybox));
  3.2913 -  }
  3.2914 -  /* strategy number 23: overlapping with circle */
  3.2915 -  if (strategy == 23) {
  3.2916 -    /* query datum is a circle */
  3.2917 -    pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
  3.2918 -    /* return true if estimated distance (allowed to be smaller than real
  3.2919 -       distance) between index key and circle center is smaller than radius */
  3.2920 -    PG_RETURN_BOOL(
  3.2921 -      pgl_estimate_key_distance(key, &(query->center)) <= query->radius
  3.2922 -    );
  3.2923 -  }
  3.2924 -  /* strategy number 24: overlapping with cluster */
  3.2925 -  if (strategy == 24) {
  3.2926 -    bool retval;  /* return value */
  3.2927 -    /* query datum is a cluster */
  3.2928 -    pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  3.2929 -    /* return true if estimated distance (allowed to be smaller than real
  3.2930 -       distance) between index key and circle center is smaller than radius */
  3.2931 -    retval = (
  3.2932 -      pgl_estimate_key_distance(key, &(query->bounding.center)) <=
  3.2933 -      query->bounding.radius
  3.2934 -    );
  3.2935 -    PG_FREE_IF_COPY(query, 1);  /* free detoasted cluster (if copy) */
  3.2936 -    PG_RETURN_BOOL(retval);
  3.2937 -  }
  3.2938 -  /* throw error for any unknown strategy number */
  3.2939 -  elog(ERROR, "unrecognized strategy number: %d", strategy);
  3.2940 -}
  3.2941 -
  3.2942 -/* GiST "union" support function */
  3.2943 -PG_FUNCTION_INFO_V1(pgl_gist_union);
  3.2944 -Datum pgl_gist_union(PG_FUNCTION_ARGS) {
  3.2945 -  GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0);
  3.2946 -  pgl_keyptr out;  /* return value (to be palloc'ed) */
  3.2947 -  int i;
  3.2948 -  /* determine key size */
  3.2949 -  size_t keysize = PGL_KEY_IS_AREAKEY(
  3.2950 -    (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key)
  3.2951 -  ) ? sizeof (pgl_areakey) : sizeof(pgl_pointkey);
  3.2952 -  /* begin with first key as result */
  3.2953 -  out = palloc(keysize);
  3.2954 -  memcpy(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[0].key), keysize);
  3.2955 -  /* unite current result with second, third, etc. key */
  3.2956 -  for (i=1; i<entryvec->n; i++) {
  3.2957 -    pgl_unite_keys(out, (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key));
  3.2958 -  }
  3.2959 -  /* return result */
  3.2960 -  PG_RETURN_POINTER(out);
  3.2961 -}
  3.2962 -
  3.2963 -/* GiST "compress" support function for indicis on points */
  3.2964 -PG_FUNCTION_INFO_V1(pgl_gist_compress_epoint);
  3.2965 -Datum pgl_gist_compress_epoint(PG_FUNCTION_ARGS) {
  3.2966 -  GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  3.2967 -  GISTENTRY *retval;  /* return value (to be palloc'ed unless set to entry) */
  3.2968 -  /* only transform new leaves */
  3.2969 -  if (entry->leafkey) {
  3.2970 -    /* get point to be transformed */
  3.2971 -    pgl_point *point = (pgl_point *)DatumGetPointer(entry->key);
  3.2972 -    /* allocate memory for key */
  3.2973 -    pgl_keyptr key = palloc(sizeof(pgl_pointkey));
  3.2974 -    /* transform point to key */
  3.2975 -    pgl_point_to_key(point, key);
  3.2976 -    /* create new GISTENTRY structure as return value */
  3.2977 -    retval = palloc(sizeof(GISTENTRY));
  3.2978 -    gistentryinit(
  3.2979 -      *retval, PointerGetDatum(key),
  3.2980 -      entry->rel, entry->page, entry->offset, FALSE
  3.2981 -    );
  3.2982 -  } else {
  3.2983 -    /* inner nodes have already been transformed */
  3.2984 -    retval = entry;
  3.2985 -  }
  3.2986 -  /* return pointer to old or new GISTENTRY structure */
  3.2987 -  PG_RETURN_POINTER(retval);
  3.2988 -}
  3.2989 -
  3.2990 -/* GiST "compress" support function for indicis on circles */
  3.2991 -PG_FUNCTION_INFO_V1(pgl_gist_compress_ecircle);
  3.2992 -Datum pgl_gist_compress_ecircle(PG_FUNCTION_ARGS) {
  3.2993 -  GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  3.2994 -  GISTENTRY *retval;  /* return value (to be palloc'ed unless set to entry) */
  3.2995 -  /* only transform new leaves */
  3.2996 -  if (entry->leafkey) {
  3.2997 -    /* get circle to be transformed */
  3.2998 -    pgl_circle *circle = (pgl_circle *)DatumGetPointer(entry->key);
  3.2999 -    /* allocate memory for key */
  3.3000 -    pgl_keyptr key = palloc(sizeof(pgl_areakey));
  3.3001 -    /* transform circle to key */
  3.3002 -    pgl_circle_to_key(circle, key);
  3.3003 -    /* create new GISTENTRY structure as return value */
  3.3004 -    retval = palloc(sizeof(GISTENTRY));
  3.3005 -    gistentryinit(
  3.3006 -      *retval, PointerGetDatum(key),
  3.3007 -      entry->rel, entry->page, entry->offset, FALSE
  3.3008 -    );
  3.3009 -  } else {
  3.3010 -    /* inner nodes have already been transformed */
  3.3011 -    retval = entry;
  3.3012 -  }
  3.3013 -  /* return pointer to old or new GISTENTRY structure */
  3.3014 -  PG_RETURN_POINTER(retval);
  3.3015 -}
  3.3016 -
  3.3017 -/* GiST "compress" support function for indices on clusters */
  3.3018 -PG_FUNCTION_INFO_V1(pgl_gist_compress_ecluster);
  3.3019 -Datum pgl_gist_compress_ecluster(PG_FUNCTION_ARGS) {
  3.3020 -  GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  3.3021 -  GISTENTRY *retval;  /* return value (to be palloc'ed unless set to entry) */
  3.3022 -  /* only transform new leaves */
  3.3023 -  if (entry->leafkey) {
  3.3024 -    /* get cluster to be transformed (detoasting necessary!) */
  3.3025 -    pgl_cluster *cluster = (pgl_cluster *)PG_DETOAST_DATUM(entry->key);
  3.3026 -    /* allocate memory for key */
  3.3027 -    pgl_keyptr key = palloc(sizeof(pgl_areakey));
  3.3028 -    /* transform cluster to key */
  3.3029 -    pgl_circle_to_key(&(cluster->bounding), key);
  3.3030 -    /* create new GISTENTRY structure as return value */
  3.3031 -    retval = palloc(sizeof(GISTENTRY));
  3.3032 -    gistentryinit(
  3.3033 -      *retval, PointerGetDatum(key),
  3.3034 -      entry->rel, entry->page, entry->offset, FALSE
  3.3035 -    );
  3.3036 -    /* free detoasted datum */
  3.3037 -    if ((void *)cluster != (void *)DatumGetPointer(entry->key)) pfree(cluster);
  3.3038 -  } else {
  3.3039 -    /* inner nodes have already been transformed */
  3.3040 -    retval = entry;
  3.3041 -  }
  3.3042 -  /* return pointer to old or new GISTENTRY structure */
  3.3043 -  PG_RETURN_POINTER(retval);
  3.3044 -}
  3.3045 -
  3.3046 -/* GiST "decompress" support function for indices */
  3.3047 -PG_FUNCTION_INFO_V1(pgl_gist_decompress);
  3.3048 -Datum pgl_gist_decompress(PG_FUNCTION_ARGS) {
  3.3049 -  /* return passed pointer without transformation */
  3.3050 -  PG_RETURN_POINTER(PG_GETARG_POINTER(0));
  3.3051 -}
  3.3052 -
  3.3053 -/* GiST "penalty" support function */
  3.3054 -PG_FUNCTION_INFO_V1(pgl_gist_penalty);
  3.3055 -Datum pgl_gist_penalty(PG_FUNCTION_ARGS) {
  3.3056 -  GISTENTRY *origentry = (GISTENTRY *)PG_GETARG_POINTER(0);
  3.3057 -  GISTENTRY *newentry = (GISTENTRY *)PG_GETARG_POINTER(1);
  3.3058 -  float *penalty = (float *)PG_GETARG_POINTER(2);
  3.3059 -  /* get original key and key to insert */
  3.3060 -  pgl_keyptr orig = (pgl_keyptr)DatumGetPointer(origentry->key);
  3.3061 -  pgl_keyptr new = (pgl_keyptr)DatumGetPointer(newentry->key);
  3.3062 -  /* copy original key */
  3.3063 -  union { pgl_pointkey pointkey; pgl_areakey areakey; } union_key;
  3.3064 -  if (PGL_KEY_IS_AREAKEY(orig)) {
  3.3065 -    memcpy(union_key.areakey, orig, sizeof(union_key.areakey));
  3.3066 -  } else {
  3.3067 -    memcpy(union_key.pointkey, orig, sizeof(union_key.pointkey));
  3.3068 -  }
  3.3069 -  /* calculate union of both keys */
  3.3070 -  pgl_unite_keys((pgl_keyptr)&union_key, new);
  3.3071 -  /* penalty equal to reduction of key length (logarithm of added area) */
  3.3072 -  /* (return value by setting referenced value and returning pointer) */
  3.3073 -  *penalty = (
  3.3074 -    PGL_KEY_NODEDEPTH(orig) - PGL_KEY_NODEDEPTH((pgl_keyptr)&union_key)
  3.3075 -  );
  3.3076 -  PG_RETURN_POINTER(penalty);
  3.3077 -}
  3.3078 -
  3.3079 -/* GiST "picksplit" support function */
  3.3080 -PG_FUNCTION_INFO_V1(pgl_gist_picksplit);
  3.3081 -Datum pgl_gist_picksplit(PG_FUNCTION_ARGS) {
  3.3082 -  GistEntryVector *entryvec = (GistEntryVector *)PG_GETARG_POINTER(0);
  3.3083 -  GIST_SPLITVEC *v = (GIST_SPLITVEC *)PG_GETARG_POINTER(1);
  3.3084 -  OffsetNumber i;  /* between FirstOffsetNumber and entryvec->n (inclusive) */
  3.3085 -  union {
  3.3086 -    pgl_pointkey pointkey;
  3.3087 -    pgl_areakey areakey;
  3.3088 -  } union_all;  /* union of all keys (to be calculated from scratch)
  3.3089 -                   (later cut in half) */
  3.3090 -  int is_areakey = PGL_KEY_IS_AREAKEY(
  3.3091 -    (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key)
  3.3092 -  );
  3.3093 -  int keysize = is_areakey ? sizeof(pgl_areakey) : sizeof(pgl_pointkey);
  3.3094 -  pgl_keyptr unionL = palloc(keysize);  /* union of keys that go left */
  3.3095 -  pgl_keyptr unionR = palloc(keysize);  /* union of keys that go right */
  3.3096 -  pgl_keyptr key;  /* current key to be processed */
  3.3097 -  /* allocate memory for array of left and right keys, set counts to zero */
  3.3098 -  v->spl_left = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber));
  3.3099 -  v->spl_nleft = 0;
  3.3100 -  v->spl_right = (OffsetNumber *)palloc(entryvec->n * sizeof(OffsetNumber));
  3.3101 -  v->spl_nright = 0;
  3.3102 -  /* calculate union of all keys from scratch */
  3.3103 -  memcpy(
  3.3104 -    (pgl_keyptr)&union_all,
  3.3105 -    (pgl_keyptr)DatumGetPointer(entryvec->vector[FirstOffsetNumber].key),
  3.3106 -    keysize
  3.3107 -  );
  3.3108 -  for (i=FirstOffsetNumber+1; i<entryvec->n; i=OffsetNumberNext(i)) {
  3.3109 -    pgl_unite_keys(
  3.3110 -      (pgl_keyptr)&union_all,
  3.3111 -      (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key)
  3.3112 -    );
  3.3113 -  }
  3.3114 -  /* check if trivial split is necessary due to exhausted key length */
  3.3115 -  /* (Note: keys for empty objects must have node depth set to maximum) */
  3.3116 -  if (PGL_KEY_NODEDEPTH((pgl_keyptr)&union_all) == (
  3.3117 -    is_areakey ? PGL_AREAKEY_MAXDEPTH : PGL_POINTKEY_MAXDEPTH
  3.3118 -  )) {
  3.3119 -    /* half of all keys go left */
  3.3120 -    for (
  3.3121 -      i=FirstOffsetNumber;
  3.3122 -      i<FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2;
  3.3123 -      i=OffsetNumberNext(i)
  3.3124 -    ) {
  3.3125 -      /* pointer to current key */
  3.3126 -      key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
  3.3127 -      /* update unionL */
  3.3128 -      /* check if key is first key that goes left */
  3.3129 -      if (!v->spl_nleft) {
  3.3130 -        /* first key that goes left is just copied to unionL */
  3.3131 -        memcpy(unionL, key, keysize);
  3.3132 -      } else {
  3.3133 -        /* unite current value and next key */
  3.3134 -        pgl_unite_keys(unionL, key);
  3.3135 -      }
  3.3136 -      /* append offset number to list of keys that go left */
  3.3137 -      v->spl_left[v->spl_nleft++] = i;
  3.3138 -    }
  3.3139 -    /* other half goes right */
  3.3140 -    for (
  3.3141 -      i=FirstOffsetNumber+(entryvec->n - FirstOffsetNumber)/2;
  3.3142 -      i<entryvec->n;
  3.3143 -      i=OffsetNumberNext(i)
  3.3144 -    ) {
  3.3145 -      /* pointer to current key */
  3.3146 -      key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
  3.3147 -      /* update unionR */
  3.3148 -      /* check if key is first key that goes right */
  3.3149 -      if (!v->spl_nright) {
  3.3150 -        /* first key that goes right is just copied to unionR */
  3.3151 -        memcpy(unionR, key, keysize);
  3.3152 -      } else {
  3.3153 -        /* unite current value and next key */
  3.3154 -        pgl_unite_keys(unionR, key);
  3.3155 -      }
  3.3156 -      /* append offset number to list of keys that go right */
  3.3157 -      v->spl_right[v->spl_nright++] = i;
  3.3158 -    }
  3.3159 -  }
  3.3160 -  /* otherwise, a non-trivial split is possible */
  3.3161 -  else {
  3.3162 -    /* cut covered area in half */
  3.3163 -    /* (union_all then refers to area of keys that go left) */
  3.3164 -    /* check if union of all keys covers empty and non-empty objects */
  3.3165 -    if (PGL_KEY_IS_UNIVERSAL((pgl_keyptr)&union_all)) {
  3.3166 -      /* if yes, split into empty and non-empty objects */
  3.3167 -      pgl_key_set_empty((pgl_keyptr)&union_all);
  3.3168 -    } else {
  3.3169 -      /* otherwise split by next bit */
  3.3170 -      ((pgl_keyptr)&union_all)[PGL_KEY_NODEDEPTH_OFFSET]++;
  3.3171 -      /* NOTE: type bit conserved */
  3.3172 -    }
  3.3173 -    /* determine for each key if it goes left or right */
  3.3174 -    for (i=FirstOffsetNumber; i<entryvec->n; i=OffsetNumberNext(i)) {
  3.3175 -      /* pointer to current key */
  3.3176 -      key = (pgl_keyptr)DatumGetPointer(entryvec->vector[i].key);
  3.3177 -      /* keys within one half of the area go left */
  3.3178 -      if (pgl_keys_overlap((pgl_keyptr)&union_all, key)) {
  3.3179 -        /* update unionL */
  3.3180 -        /* check if key is first key that goes left */
  3.3181 -        if (!v->spl_nleft) {
  3.3182 -          /* first key that goes left is just copied to unionL */
  3.3183 -          memcpy(unionL, key, keysize);
  3.3184 -        } else {
  3.3185 -          /* unite current value of unionL and processed key */
  3.3186 -          pgl_unite_keys(unionL, key);
  3.3187 -        }
  3.3188 -        /* append offset number to list of keys that go left */
  3.3189 -        v->spl_left[v->spl_nleft++] = i;
  3.3190 -      }
  3.3191 -      /* the other keys go right */
  3.3192 -      else {
  3.3193 -        /* update unionR */
  3.3194 -        /* check if key is first key that goes right */
  3.3195 -        if (!v->spl_nright) {
  3.3196 -          /* first key that goes right is just copied to unionR */
  3.3197 -          memcpy(unionR, key, keysize);
  3.3198 -        } else {
  3.3199 -          /* unite current value of unionR and processed key */
  3.3200 -          pgl_unite_keys(unionR, key);
  3.3201 -        }
  3.3202 -        /* append offset number to list of keys that go right */
  3.3203 -        v->spl_right[v->spl_nright++] = i;
  3.3204 -      }
  3.3205 -    }
  3.3206 -  }
  3.3207 -  /* store unions in return value */
  3.3208 -  v->spl_ldatum = PointerGetDatum(unionL);
  3.3209 -  v->spl_rdatum = PointerGetDatum(unionR);
  3.3210 -  /* return all results */
  3.3211 -  PG_RETURN_POINTER(v);
  3.3212 -}
  3.3213 -
  3.3214 -/* GiST "same"/"equal" support function */
  3.3215 -PG_FUNCTION_INFO_V1(pgl_gist_same);
  3.3216 -Datum pgl_gist_same(PG_FUNCTION_ARGS) {
  3.3217 -  pgl_keyptr key1 = (pgl_keyptr)PG_GETARG_POINTER(0);
  3.3218 -  pgl_keyptr key2 = (pgl_keyptr)PG_GETARG_POINTER(1);
  3.3219 -  bool *boolptr = (bool *)PG_GETARG_POINTER(2);
  3.3220 -  /* two keys are equal if they are binary equal */
  3.3221 -  /* (return result by setting referenced boolean and returning pointer) */
  3.3222 -  *boolptr = !memcmp(
  3.3223 -    key1,
  3.3224 -    key2,
  3.3225 -    PGL_KEY_IS_AREAKEY(key1) ? sizeof(pgl_areakey) : sizeof(pgl_pointkey)
  3.3226 -  );
  3.3227 -  PG_RETURN_POINTER(boolptr);
  3.3228 -}
  3.3229 -
  3.3230 -/* GiST "distance" support function */
  3.3231 -PG_FUNCTION_INFO_V1(pgl_gist_distance);
  3.3232 -Datum pgl_gist_distance(PG_FUNCTION_ARGS) {
  3.3233 -  GISTENTRY *entry = (GISTENTRY *)PG_GETARG_POINTER(0);
  3.3234 -  pgl_keyptr key = (pgl_keyptr)DatumGetPointer(entry->key);
  3.3235 -  StrategyNumber strategy = (StrategyNumber)PG_GETARG_UINT16(2);
  3.3236 -  bool *recheck = (bool *)PG_GETARG_POINTER(4);
  3.3237 -  double distance;  /* return value */
  3.3238 -  /* demand recheck because distance is just an estimation */
  3.3239 -  /* (real distance may be bigger) */
  3.3240 -  *recheck = true;
  3.3241 -  /* strategy number aliases for different operators using the same strategy */
  3.3242 -  strategy %= 100;
  3.3243 -  /* strategy number 31: distance to point */
  3.3244 -  if (strategy == 31) {
  3.3245 -    /* query datum is a point */
  3.3246 -    pgl_point *query = (pgl_point *)PG_GETARG_POINTER(1);
  3.3247 -    /* use pgl_estimate_pointkey_distance() function to compute result */
  3.3248 -    distance = pgl_estimate_key_distance(key, query);
  3.3249 -    /* avoid infinity (reserved!) */
  3.3250 -    if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
  3.3251 -    /* return result */
  3.3252 -    PG_RETURN_FLOAT8(distance);
  3.3253 -  }
  3.3254 -  /* strategy number 33: distance to circle */
  3.3255 -  if (strategy == 33) {
  3.3256 -    /* query datum is a circle */
  3.3257 -    pgl_circle *query = (pgl_circle *)PG_GETARG_POINTER(1);
  3.3258 -    /* estimate distance to circle center and substract circle radius */
  3.3259 -    distance = (
  3.3260 -      pgl_estimate_key_distance(key, &(query->center)) - query->radius
  3.3261 -    );
  3.3262 -    /* convert non-positive values to zero and avoid infinity (reserved!) */
  3.3263 -    if (distance <= 0) distance = 0;
  3.3264 -    else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
  3.3265 -    /* return result */
  3.3266 -    PG_RETURN_FLOAT8(distance);
  3.3267 -  }
  3.3268 -  /* strategy number 34: distance to cluster */
  3.3269 -  if (strategy == 34) {
  3.3270 -    /* query datum is a cluster */
  3.3271 -    pgl_cluster *query = (pgl_cluster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
  3.3272 -    /* estimate distance to bounding center and substract bounding radius */
  3.3273 -    distance = (
  3.3274 -      pgl_estimate_key_distance(key, &(query->bounding.center)) -
  3.3275 -      query->bounding.radius
  3.3276 -    );
  3.3277 -    /* convert non-positive values to zero and avoid infinity (reserved!) */
  3.3278 -    if (distance <= 0) distance = 0;
  3.3279 -    else if (!isfinite(distance)) distance = PGL_ULTRA_DISTANCE;
  3.3280 -    /* free detoasted cluster (if copy) */
  3.3281 -    PG_FREE_IF_COPY(query, 1);
  3.3282 -    /* return result */
  3.3283 -    PG_RETURN_FLOAT8(distance);
  3.3284 -  }
  3.3285 -  /* throw error for any unknown strategy number */
  3.3286 -  elog(ERROR, "unrecognized strategy number: %d", strategy);
  3.3287 -}
  3.3288 -