1 /******************************************************************************
2 contrib/cube/cube.c
3
4 This file contains routines that can be bound to a Postgres backend and
5 called by the backend in the process of processing queries. The calling
6 format for these routines is dictated by Postgres architecture.
7 ******************************************************************************/
8
9 #include "postgres.h"
10
11 #include <float.h>
12 #include <math.h>
13
14 #include "access/gist.h"
15 #include "access/stratnum.h"
16 #include "utils/array.h"
17 #include "utils/builtins.h"
18
19 #include "cubedata.h"
20
21 PG_MODULE_MAGIC;
22
23 /*
24 * Taken from the intarray contrib header
25 */
26 #define ARRPTR(x) ( (double *) ARR_DATA_PTR(x) )
27 #define ARRNELEMS(x) ArrayGetNItems( ARR_NDIM(x), ARR_DIMS(x))
28
29 /*
30 ** Input/Output routines
31 */
32 PG_FUNCTION_INFO_V1(cube_in);
33 PG_FUNCTION_INFO_V1(cube_a_f8_f8);
34 PG_FUNCTION_INFO_V1(cube_a_f8);
35 PG_FUNCTION_INFO_V1(cube_out);
36 PG_FUNCTION_INFO_V1(cube_f8);
37 PG_FUNCTION_INFO_V1(cube_f8_f8);
38 PG_FUNCTION_INFO_V1(cube_c_f8);
39 PG_FUNCTION_INFO_V1(cube_c_f8_f8);
40 PG_FUNCTION_INFO_V1(cube_dim);
41 PG_FUNCTION_INFO_V1(cube_ll_coord);
42 PG_FUNCTION_INFO_V1(cube_ur_coord);
43 PG_FUNCTION_INFO_V1(cube_coord);
44 PG_FUNCTION_INFO_V1(cube_coord_llur);
45 PG_FUNCTION_INFO_V1(cube_subset);
46
47 /*
48 ** GiST support methods
49 */
50
51 PG_FUNCTION_INFO_V1(g_cube_consistent);
52 PG_FUNCTION_INFO_V1(g_cube_compress);
53 PG_FUNCTION_INFO_V1(g_cube_decompress);
54 PG_FUNCTION_INFO_V1(g_cube_penalty);
55 PG_FUNCTION_INFO_V1(g_cube_picksplit);
56 PG_FUNCTION_INFO_V1(g_cube_union);
57 PG_FUNCTION_INFO_V1(g_cube_same);
58 PG_FUNCTION_INFO_V1(g_cube_distance);
59
60 /*
61 ** B-tree support functions
62 */
63 PG_FUNCTION_INFO_V1(cube_eq);
64 PG_FUNCTION_INFO_V1(cube_ne);
65 PG_FUNCTION_INFO_V1(cube_lt);
66 PG_FUNCTION_INFO_V1(cube_gt);
67 PG_FUNCTION_INFO_V1(cube_le);
68 PG_FUNCTION_INFO_V1(cube_ge);
69 PG_FUNCTION_INFO_V1(cube_cmp);
70
71 /*
72 ** R-tree support functions
73 */
74
75 PG_FUNCTION_INFO_V1(cube_contains);
76 PG_FUNCTION_INFO_V1(cube_contained);
77 PG_FUNCTION_INFO_V1(cube_overlap);
78 PG_FUNCTION_INFO_V1(cube_union);
79 PG_FUNCTION_INFO_V1(cube_inter);
80 PG_FUNCTION_INFO_V1(cube_size);
81
82 /*
83 ** miscellaneous
84 */
85 PG_FUNCTION_INFO_V1(distance_taxicab);
86 PG_FUNCTION_INFO_V1(cube_distance);
87 PG_FUNCTION_INFO_V1(distance_chebyshev);
88 PG_FUNCTION_INFO_V1(cube_is_point);
89 PG_FUNCTION_INFO_V1(cube_enlarge);
90
91 /*
92 ** For internal use only
93 */
94 int32 cube_cmp_v0(NDBOX *a, NDBOX *b);
95 bool cube_contains_v0(NDBOX *a, NDBOX *b);
96 bool cube_overlap_v0(NDBOX *a, NDBOX *b);
97 NDBOX *cube_union_v0(NDBOX *a, NDBOX *b);
98 void rt_cube_size(NDBOX *a, double *sz);
99 NDBOX *g_cube_binary_union(NDBOX *r1, NDBOX *r2, int *sizep);
100 bool g_cube_leaf_consistent(NDBOX *key, NDBOX *query, StrategyNumber strategy);
101 bool g_cube_internal_consistent(NDBOX *key, NDBOX *query, StrategyNumber strategy);
102
103 /*
104 ** Auxiliary funxtions
105 */
106 static double distance_1D(double a1, double a2, double b1, double b2);
107 static bool cube_is_point_internal(NDBOX *cube);
108
109
110 /*****************************************************************************
111 * Input/Output functions
112 *****************************************************************************/
113
114 /* NdBox = [(lowerleft),(upperright)] */
115 /* [(xLL(1)...xLL(N)),(xUR(1)...xUR(n))] */
116 Datum
cube_in(PG_FUNCTION_ARGS)117 cube_in(PG_FUNCTION_ARGS)
118 {
119 char *str = PG_GETARG_CSTRING(0);
120 NDBOX *result;
121
122 cube_scanner_init(str);
123
124 if (cube_yyparse(&result) != 0)
125 cube_yyerror(&result, "cube parser failed");
126
127 cube_scanner_finish();
128
129 PG_RETURN_NDBOX(result);
130 }
131
132
133 /*
134 ** Allows the construction of a cube from 2 float[]'s
135 */
136 Datum
cube_a_f8_f8(PG_FUNCTION_ARGS)137 cube_a_f8_f8(PG_FUNCTION_ARGS)
138 {
139 ArrayType *ur = PG_GETARG_ARRAYTYPE_P(0);
140 ArrayType *ll = PG_GETARG_ARRAYTYPE_P(1);
141 NDBOX *result;
142 int i;
143 int dim;
144 int size;
145 bool point;
146 double *dur,
147 *dll;
148
149 if (array_contains_nulls(ur) || array_contains_nulls(ll))
150 ereport(ERROR,
151 (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
152 errmsg("cannot work with arrays containing NULLs")));
153
154 dim = ARRNELEMS(ur);
155 if (dim > CUBE_MAX_DIM)
156 ereport(ERROR,
157 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
158 errmsg("can't extend cube"),
159 errdetail("A cube cannot have more than %d dimensions.",
160 CUBE_MAX_DIM)));
161
162 if (ARRNELEMS(ll) != dim)
163 ereport(ERROR,
164 (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
165 errmsg("UR and LL arrays must be of same length")));
166
167 dur = ARRPTR(ur);
168 dll = ARRPTR(ll);
169
170 /* Check if it's a point */
171 point = true;
172 for (i = 0; i < dim; i++)
173 {
174 if (dur[i] != dll[i])
175 {
176 point = false;
177 break;
178 }
179 }
180
181 size = point ? POINT_SIZE(dim) : CUBE_SIZE(dim);
182 result = (NDBOX *) palloc0(size);
183 SET_VARSIZE(result, size);
184 SET_DIM(result, dim);
185
186 for (i = 0; i < dim; i++)
187 result->x[i] = dur[i];
188
189 if (!point)
190 {
191 for (i = 0; i < dim; i++)
192 result->x[i + dim] = dll[i];
193 }
194 else
195 SET_POINT_BIT(result);
196
197 PG_RETURN_NDBOX(result);
198 }
199
200 /*
201 ** Allows the construction of a zero-volume cube from a float[]
202 */
203 Datum
cube_a_f8(PG_FUNCTION_ARGS)204 cube_a_f8(PG_FUNCTION_ARGS)
205 {
206 ArrayType *ur = PG_GETARG_ARRAYTYPE_P(0);
207 NDBOX *result;
208 int i;
209 int dim;
210 int size;
211 double *dur;
212
213 if (array_contains_nulls(ur))
214 ereport(ERROR,
215 (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
216 errmsg("cannot work with arrays containing NULLs")));
217
218 dim = ARRNELEMS(ur);
219 if (dim > CUBE_MAX_DIM)
220 ereport(ERROR,
221 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
222 errmsg("array is too long"),
223 errdetail("A cube cannot have more than %d dimensions.",
224 CUBE_MAX_DIM)));
225
226 dur = ARRPTR(ur);
227
228 size = POINT_SIZE(dim);
229 result = (NDBOX *) palloc0(size);
230 SET_VARSIZE(result, size);
231 SET_DIM(result, dim);
232 SET_POINT_BIT(result);
233
234 for (i = 0; i < dim; i++)
235 result->x[i] = dur[i];
236
237 PG_RETURN_NDBOX(result);
238 }
239
240 Datum
cube_subset(PG_FUNCTION_ARGS)241 cube_subset(PG_FUNCTION_ARGS)
242 {
243 NDBOX *c = PG_GETARG_NDBOX(0);
244 ArrayType *idx = PG_GETARG_ARRAYTYPE_P(1);
245 NDBOX *result;
246 int size,
247 dim,
248 i;
249 int *dx;
250
251 if (array_contains_nulls(idx))
252 ereport(ERROR,
253 (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
254 errmsg("cannot work with arrays containing NULLs")));
255
256 dx = (int32 *) ARR_DATA_PTR(idx);
257
258 dim = ARRNELEMS(idx);
259 if (dim > CUBE_MAX_DIM)
260 ereport(ERROR,
261 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
262 errmsg("array is too long"),
263 errdetail("A cube cannot have more than %d dimensions.",
264 CUBE_MAX_DIM)));
265
266 size = IS_POINT(c) ? POINT_SIZE(dim) : CUBE_SIZE(dim);
267 result = (NDBOX *) palloc0(size);
268 SET_VARSIZE(result, size);
269 SET_DIM(result, dim);
270
271 if (IS_POINT(c))
272 SET_POINT_BIT(result);
273
274 for (i = 0; i < dim; i++)
275 {
276 if ((dx[i] <= 0) || (dx[i] > DIM(c)))
277 ereport(ERROR,
278 (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
279 errmsg("Index out of bounds")));
280 result->x[i] = c->x[dx[i] - 1];
281 if (!IS_POINT(c))
282 result->x[i + dim] = c->x[dx[i] + DIM(c) - 1];
283 }
284
285 PG_FREE_IF_COPY(c, 0);
286 PG_RETURN_NDBOX(result);
287 }
288
289 Datum
cube_out(PG_FUNCTION_ARGS)290 cube_out(PG_FUNCTION_ARGS)
291 {
292 NDBOX *cube = PG_GETARG_NDBOX(0);
293 StringInfoData buf;
294 int dim = DIM(cube);
295 int i;
296
297 initStringInfo(&buf);
298
299 appendStringInfoChar(&buf, '(');
300 for (i = 0; i < dim; i++)
301 {
302 if (i > 0)
303 appendStringInfoString(&buf, ", ");
304 appendStringInfoString(&buf, float8out_internal(LL_COORD(cube, i)));
305 }
306 appendStringInfoChar(&buf, ')');
307
308 if (!cube_is_point_internal(cube))
309 {
310 appendStringInfoString(&buf, ",(");
311 for (i = 0; i < dim; i++)
312 {
313 if (i > 0)
314 appendStringInfoString(&buf, ", ");
315 appendStringInfoString(&buf, float8out_internal(UR_COORD(cube, i)));
316 }
317 appendStringInfoChar(&buf, ')');
318 }
319
320 PG_FREE_IF_COPY(cube, 0);
321 PG_RETURN_CSTRING(buf.data);
322 }
323
324
325 /*****************************************************************************
326 * GiST functions
327 *****************************************************************************/
328
329 /*
330 ** The GiST Consistent method for boxes
331 ** Should return false if for all data items x below entry,
332 ** the predicate x op query == FALSE, where op is the oper
333 ** corresponding to strategy in the pg_amop table.
334 */
335 Datum
g_cube_consistent(PG_FUNCTION_ARGS)336 g_cube_consistent(PG_FUNCTION_ARGS)
337 {
338 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
339 NDBOX *query = PG_GETARG_NDBOX(1);
340 StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
341
342 /* Oid subtype = PG_GETARG_OID(3); */
343 bool *recheck = (bool *) PG_GETARG_POINTER(4);
344 bool res;
345
346 /* All cases served by this function are exact */
347 *recheck = false;
348
349 /*
350 * if entry is not leaf, use g_cube_internal_consistent, else use
351 * g_cube_leaf_consistent
352 */
353 if (GIST_LEAF(entry))
354 res = g_cube_leaf_consistent(DatumGetNDBOX(entry->key),
355 query, strategy);
356 else
357 res = g_cube_internal_consistent(DatumGetNDBOX(entry->key),
358 query, strategy);
359
360 PG_FREE_IF_COPY(query, 1);
361 PG_RETURN_BOOL(res);
362 }
363
364
365 /*
366 ** The GiST Union method for boxes
367 ** returns the minimal bounding box that encloses all the entries in entryvec
368 */
369 Datum
g_cube_union(PG_FUNCTION_ARGS)370 g_cube_union(PG_FUNCTION_ARGS)
371 {
372 GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
373 int *sizep = (int *) PG_GETARG_POINTER(1);
374 NDBOX *out = (NDBOX *) NULL;
375 NDBOX *tmp;
376 int i;
377
378 tmp = DatumGetNDBOX(entryvec->vector[0].key);
379
380 /*
381 * sizep = sizeof(NDBOX); -- NDBOX has variable size
382 */
383 *sizep = VARSIZE(tmp);
384
385 for (i = 1; i < entryvec->n; i++)
386 {
387 out = g_cube_binary_union(tmp,
388 DatumGetNDBOX(entryvec->vector[i].key),
389 sizep);
390 tmp = out;
391 }
392
393 PG_RETURN_POINTER(out);
394 }
395
396 /*
397 ** GiST Compress and Decompress methods for boxes
398 ** do not do anything.
399 */
400
401 Datum
g_cube_compress(PG_FUNCTION_ARGS)402 g_cube_compress(PG_FUNCTION_ARGS)
403 {
404 PG_RETURN_DATUM(PG_GETARG_DATUM(0));
405 }
406
407 Datum
g_cube_decompress(PG_FUNCTION_ARGS)408 g_cube_decompress(PG_FUNCTION_ARGS)
409 {
410 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
411 NDBOX *key = DatumGetNDBOX(PG_DETOAST_DATUM(entry->key));
412
413 if (key != DatumGetNDBOX(entry->key))
414 {
415 GISTENTRY *retval = (GISTENTRY *) palloc(sizeof(GISTENTRY));
416
417 gistentryinit(*retval, PointerGetDatum(key),
418 entry->rel, entry->page,
419 entry->offset, FALSE);
420 PG_RETURN_POINTER(retval);
421 }
422 PG_RETURN_POINTER(entry);
423 }
424
425
426 /*
427 ** The GiST Penalty method for boxes
428 ** As in the R-tree paper, we use change in area as our penalty metric
429 */
430 Datum
g_cube_penalty(PG_FUNCTION_ARGS)431 g_cube_penalty(PG_FUNCTION_ARGS)
432 {
433 GISTENTRY *origentry = (GISTENTRY *) PG_GETARG_POINTER(0);
434 GISTENTRY *newentry = (GISTENTRY *) PG_GETARG_POINTER(1);
435 float *result = (float *) PG_GETARG_POINTER(2);
436 NDBOX *ud;
437 double tmp1,
438 tmp2;
439
440 ud = cube_union_v0(DatumGetNDBOX(origentry->key),
441 DatumGetNDBOX(newentry->key));
442 rt_cube_size(ud, &tmp1);
443 rt_cube_size(DatumGetNDBOX(origentry->key), &tmp2);
444 *result = (float) (tmp1 - tmp2);
445
446 PG_RETURN_FLOAT8(*result);
447 }
448
449
450
451 /*
452 ** The GiST PickSplit method for boxes
453 ** We use Guttman's poly time split algorithm
454 */
455 Datum
g_cube_picksplit(PG_FUNCTION_ARGS)456 g_cube_picksplit(PG_FUNCTION_ARGS)
457 {
458 GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
459 GIST_SPLITVEC *v = (GIST_SPLITVEC *) PG_GETARG_POINTER(1);
460 OffsetNumber i,
461 j;
462 NDBOX *datum_alpha,
463 *datum_beta;
464 NDBOX *datum_l,
465 *datum_r;
466 NDBOX *union_d,
467 *union_dl,
468 *union_dr;
469 NDBOX *inter_d;
470 bool firsttime;
471 double size_alpha,
472 size_beta,
473 size_union,
474 size_inter;
475 double size_waste,
476 waste;
477 double size_l,
478 size_r;
479 int nbytes;
480 OffsetNumber seed_1 = 1,
481 seed_2 = 2;
482 OffsetNumber *left,
483 *right;
484 OffsetNumber maxoff;
485
486 maxoff = entryvec->n - 2;
487 nbytes = (maxoff + 2) * sizeof(OffsetNumber);
488 v->spl_left = (OffsetNumber *) palloc(nbytes);
489 v->spl_right = (OffsetNumber *) palloc(nbytes);
490
491 firsttime = true;
492 waste = 0.0;
493
494 for (i = FirstOffsetNumber; i < maxoff; i = OffsetNumberNext(i))
495 {
496 datum_alpha = DatumGetNDBOX(entryvec->vector[i].key);
497 for (j = OffsetNumberNext(i); j <= maxoff; j = OffsetNumberNext(j))
498 {
499 datum_beta = DatumGetNDBOX(entryvec->vector[j].key);
500
501 /* compute the wasted space by unioning these guys */
502 /* size_waste = size_union - size_inter; */
503 union_d = cube_union_v0(datum_alpha, datum_beta);
504 rt_cube_size(union_d, &size_union);
505 inter_d = DatumGetNDBOX(DirectFunctionCall2(cube_inter,
506 entryvec->vector[i].key, entryvec->vector[j].key));
507 rt_cube_size(inter_d, &size_inter);
508 size_waste = size_union - size_inter;
509
510 /*
511 * are these a more promising split than what we've already seen?
512 */
513
514 if (size_waste > waste || firsttime)
515 {
516 waste = size_waste;
517 seed_1 = i;
518 seed_2 = j;
519 firsttime = false;
520 }
521 }
522 }
523
524 left = v->spl_left;
525 v->spl_nleft = 0;
526 right = v->spl_right;
527 v->spl_nright = 0;
528
529 datum_alpha = DatumGetNDBOX(entryvec->vector[seed_1].key);
530 datum_l = cube_union_v0(datum_alpha, datum_alpha);
531 rt_cube_size(datum_l, &size_l);
532 datum_beta = DatumGetNDBOX(entryvec->vector[seed_2].key);
533 datum_r = cube_union_v0(datum_beta, datum_beta);
534 rt_cube_size(datum_r, &size_r);
535
536 /*
537 * Now split up the regions between the two seeds. An important property
538 * of this split algorithm is that the split vector v has the indices of
539 * items to be split in order in its left and right vectors. We exploit
540 * this property by doing a merge in the code that actually splits the
541 * page.
542 *
543 * For efficiency, we also place the new index tuple in this loop. This is
544 * handled at the very end, when we have placed all the existing tuples
545 * and i == maxoff + 1.
546 */
547
548 maxoff = OffsetNumberNext(maxoff);
549 for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
550 {
551 /*
552 * If we've already decided where to place this item, just put it on
553 * the right list. Otherwise, we need to figure out which page needs
554 * the least enlargement in order to store the item.
555 */
556
557 if (i == seed_1)
558 {
559 *left++ = i;
560 v->spl_nleft++;
561 continue;
562 }
563 else if (i == seed_2)
564 {
565 *right++ = i;
566 v->spl_nright++;
567 continue;
568 }
569
570 /* okay, which page needs least enlargement? */
571 datum_alpha = DatumGetNDBOX(entryvec->vector[i].key);
572 union_dl = cube_union_v0(datum_l, datum_alpha);
573 union_dr = cube_union_v0(datum_r, datum_alpha);
574 rt_cube_size(union_dl, &size_alpha);
575 rt_cube_size(union_dr, &size_beta);
576
577 /* pick which page to add it to */
578 if (size_alpha - size_l < size_beta - size_r)
579 {
580 datum_l = union_dl;
581 size_l = size_alpha;
582 *left++ = i;
583 v->spl_nleft++;
584 }
585 else
586 {
587 datum_r = union_dr;
588 size_r = size_beta;
589 *right++ = i;
590 v->spl_nright++;
591 }
592 }
593 *left = *right = FirstOffsetNumber; /* sentinel value, see dosplit() */
594
595 v->spl_ldatum = PointerGetDatum(datum_l);
596 v->spl_rdatum = PointerGetDatum(datum_r);
597
598 PG_RETURN_POINTER(v);
599 }
600
601 /*
602 ** Equality method
603 */
604 Datum
g_cube_same(PG_FUNCTION_ARGS)605 g_cube_same(PG_FUNCTION_ARGS)
606 {
607 NDBOX *b1 = PG_GETARG_NDBOX(0);
608 NDBOX *b2 = PG_GETARG_NDBOX(1);
609 bool *result = (bool *) PG_GETARG_POINTER(2);
610
611 if (cube_cmp_v0(b1, b2) == 0)
612 *result = TRUE;
613 else
614 *result = FALSE;
615
616 PG_RETURN_NDBOX(result);
617 }
618
619 /*
620 ** SUPPORT ROUTINES
621 */
622 bool
g_cube_leaf_consistent(NDBOX * key,NDBOX * query,StrategyNumber strategy)623 g_cube_leaf_consistent(NDBOX *key,
624 NDBOX *query,
625 StrategyNumber strategy)
626 {
627 bool retval;
628
629 switch (strategy)
630 {
631 case RTOverlapStrategyNumber:
632 retval = (bool) cube_overlap_v0(key, query);
633 break;
634 case RTSameStrategyNumber:
635 retval = (bool) (cube_cmp_v0(key, query) == 0);
636 break;
637 case RTContainsStrategyNumber:
638 case RTOldContainsStrategyNumber:
639 retval = (bool) cube_contains_v0(key, query);
640 break;
641 case RTContainedByStrategyNumber:
642 case RTOldContainedByStrategyNumber:
643 retval = (bool) cube_contains_v0(query, key);
644 break;
645 default:
646 retval = FALSE;
647 }
648 return (retval);
649 }
650
651 bool
g_cube_internal_consistent(NDBOX * key,NDBOX * query,StrategyNumber strategy)652 g_cube_internal_consistent(NDBOX *key,
653 NDBOX *query,
654 StrategyNumber strategy)
655 {
656 bool retval;
657
658 switch (strategy)
659 {
660 case RTOverlapStrategyNumber:
661 retval = (bool) cube_overlap_v0(key, query);
662 break;
663 case RTSameStrategyNumber:
664 case RTContainsStrategyNumber:
665 case RTOldContainsStrategyNumber:
666 retval = (bool) cube_contains_v0(key, query);
667 break;
668 case RTContainedByStrategyNumber:
669 case RTOldContainedByStrategyNumber:
670 retval = (bool) cube_overlap_v0(key, query);
671 break;
672 default:
673 retval = FALSE;
674 }
675 return (retval);
676 }
677
678 NDBOX *
g_cube_binary_union(NDBOX * r1,NDBOX * r2,int * sizep)679 g_cube_binary_union(NDBOX *r1, NDBOX *r2, int *sizep)
680 {
681 NDBOX *retval;
682
683 retval = cube_union_v0(r1, r2);
684 *sizep = VARSIZE(retval);
685
686 return (retval);
687 }
688
689
690 /* cube_union_v0 */
691 NDBOX *
cube_union_v0(NDBOX * a,NDBOX * b)692 cube_union_v0(NDBOX *a, NDBOX *b)
693 {
694 int i;
695 NDBOX *result;
696 int dim;
697 int size;
698
699 /* trivial case */
700 if (a == b)
701 return a;
702
703 /* swap the arguments if needed, so that 'a' is always larger than 'b' */
704 if (DIM(a) < DIM(b))
705 {
706 NDBOX *tmp = b;
707
708 b = a;
709 a = tmp;
710 }
711 dim = DIM(a);
712
713 size = CUBE_SIZE(dim);
714 result = palloc0(size);
715 SET_VARSIZE(result, size);
716 SET_DIM(result, dim);
717
718 /* First compute the union of the dimensions present in both args */
719 for (i = 0; i < DIM(b); i++)
720 {
721 result->x[i] = Min(
722 Min(LL_COORD(a, i), UR_COORD(a, i)),
723 Min(LL_COORD(b, i), UR_COORD(b, i))
724 );
725 result->x[i + DIM(a)] = Max(
726 Max(LL_COORD(a, i), UR_COORD(a, i)),
727 Max(LL_COORD(b, i), UR_COORD(b, i))
728 );
729 }
730 /* continue on the higher dimensions only present in 'a' */
731 for (; i < DIM(a); i++)
732 {
733 result->x[i] = Min(0,
734 Min(LL_COORD(a, i), UR_COORD(a, i))
735 );
736 result->x[i + dim] = Max(0,
737 Max(LL_COORD(a, i), UR_COORD(a, i))
738 );
739 }
740
741 /*
742 * Check if the result was in fact a point, and set the flag in the datum
743 * accordingly. (we don't bother to repalloc it smaller)
744 */
745 if (cube_is_point_internal(result))
746 {
747 size = POINT_SIZE(dim);
748 SET_VARSIZE(result, size);
749 SET_POINT_BIT(result);
750 }
751
752 return (result);
753 }
754
755 Datum
cube_union(PG_FUNCTION_ARGS)756 cube_union(PG_FUNCTION_ARGS)
757 {
758 NDBOX *a = PG_GETARG_NDBOX(0);
759 NDBOX *b = PG_GETARG_NDBOX(1);
760 NDBOX *res;
761
762 res = cube_union_v0(a, b);
763
764 PG_FREE_IF_COPY(a, 0);
765 PG_FREE_IF_COPY(b, 1);
766 PG_RETURN_NDBOX(res);
767 }
768
769 /* cube_inter */
770 Datum
cube_inter(PG_FUNCTION_ARGS)771 cube_inter(PG_FUNCTION_ARGS)
772 {
773 NDBOX *a = PG_GETARG_NDBOX(0);
774 NDBOX *b = PG_GETARG_NDBOX(1);
775 NDBOX *result;
776 bool swapped = false;
777 int i;
778 int dim;
779 int size;
780
781 /* swap the arguments if needed, so that 'a' is always larger than 'b' */
782 if (DIM(a) < DIM(b))
783 {
784 NDBOX *tmp = b;
785
786 b = a;
787 a = tmp;
788 swapped = true;
789 }
790 dim = DIM(a);
791
792 size = CUBE_SIZE(dim);
793 result = (NDBOX *) palloc0(size);
794 SET_VARSIZE(result, size);
795 SET_DIM(result, dim);
796
797 /* First compute intersection of the dimensions present in both args */
798 for (i = 0; i < DIM(b); i++)
799 {
800 result->x[i] = Max(
801 Min(LL_COORD(a, i), UR_COORD(a, i)),
802 Min(LL_COORD(b, i), UR_COORD(b, i))
803 );
804 result->x[i + DIM(a)] = Min(
805 Max(LL_COORD(a, i), UR_COORD(a, i)),
806 Max(LL_COORD(b, i), UR_COORD(b, i))
807 );
808 }
809 /* continue on the higher dimensions only present in 'a' */
810 for (; i < DIM(a); i++)
811 {
812 result->x[i] = Max(0,
813 Min(LL_COORD(a, i), UR_COORD(a, i))
814 );
815 result->x[i + DIM(a)] = Min(0,
816 Max(LL_COORD(a, i), UR_COORD(a, i))
817 );
818 }
819
820 /*
821 * Check if the result was in fact a point, and set the flag in the datum
822 * accordingly. (we don't bother to repalloc it smaller)
823 */
824 if (cube_is_point_internal(result))
825 {
826 size = POINT_SIZE(dim);
827 result = repalloc(result, size);
828 SET_VARSIZE(result, size);
829 SET_POINT_BIT(result);
830 }
831
832 if (swapped)
833 {
834 PG_FREE_IF_COPY(b, 0);
835 PG_FREE_IF_COPY(a, 1);
836 }
837 else
838 {
839 PG_FREE_IF_COPY(a, 0);
840 PG_FREE_IF_COPY(b, 1);
841 }
842
843 /*
844 * Is it OK to return a non-null intersection for non-overlapping boxes?
845 */
846 PG_RETURN_NDBOX(result);
847 }
848
849 /* cube_size */
850 Datum
cube_size(PG_FUNCTION_ARGS)851 cube_size(PG_FUNCTION_ARGS)
852 {
853 NDBOX *a = PG_GETARG_NDBOX(0);
854 double result;
855
856 rt_cube_size(a, &result);
857 PG_FREE_IF_COPY(a, 0);
858 PG_RETURN_FLOAT8(result);
859 }
860
861 void
rt_cube_size(NDBOX * a,double * size)862 rt_cube_size(NDBOX *a, double *size)
863 {
864 double result;
865 int i;
866
867 if (a == (NDBOX *) NULL)
868 {
869 /* special case for GiST */
870 result = 0.0;
871 }
872 else if (IS_POINT(a) || DIM(a) == 0)
873 {
874 /* necessarily has zero size */
875 result = 0.0;
876 }
877 else
878 {
879 result = 1.0;
880 for (i = 0; i < DIM(a); i++)
881 result *= Abs(UR_COORD(a, i) - LL_COORD(a, i));
882 }
883 *size = result;
884 }
885
886 /* make up a metric in which one box will be 'lower' than the other
887 -- this can be useful for sorting and to determine uniqueness */
888 int32
cube_cmp_v0(NDBOX * a,NDBOX * b)889 cube_cmp_v0(NDBOX *a, NDBOX *b)
890 {
891 int i;
892 int dim;
893
894 dim = Min(DIM(a), DIM(b));
895
896 /* compare the common dimensions */
897 for (i = 0; i < dim; i++)
898 {
899 if (Min(LL_COORD(a, i), UR_COORD(a, i)) >
900 Min(LL_COORD(b, i), UR_COORD(b, i)))
901 return 1;
902 if (Min(LL_COORD(a, i), UR_COORD(a, i)) <
903 Min(LL_COORD(b, i), UR_COORD(b, i)))
904 return -1;
905 }
906 for (i = 0; i < dim; i++)
907 {
908 if (Max(LL_COORD(a, i), UR_COORD(a, i)) >
909 Max(LL_COORD(b, i), UR_COORD(b, i)))
910 return 1;
911 if (Max(LL_COORD(a, i), UR_COORD(a, i)) <
912 Max(LL_COORD(b, i), UR_COORD(b, i)))
913 return -1;
914 }
915
916 /* compare extra dimensions to zero */
917 if (DIM(a) > DIM(b))
918 {
919 for (i = dim; i < DIM(a); i++)
920 {
921 if (Min(LL_COORD(a, i), UR_COORD(a, i)) > 0)
922 return 1;
923 if (Min(LL_COORD(a, i), UR_COORD(a, i)) < 0)
924 return -1;
925 }
926 for (i = dim; i < DIM(a); i++)
927 {
928 if (Max(LL_COORD(a, i), UR_COORD(a, i)) > 0)
929 return 1;
930 if (Max(LL_COORD(a, i), UR_COORD(a, i)) < 0)
931 return -1;
932 }
933
934 /*
935 * if all common dimensions are equal, the cube with more dimensions
936 * wins
937 */
938 return 1;
939 }
940 if (DIM(a) < DIM(b))
941 {
942 for (i = dim; i < DIM(b); i++)
943 {
944 if (Min(LL_COORD(b, i), UR_COORD(b, i)) > 0)
945 return -1;
946 if (Min(LL_COORD(b, i), UR_COORD(b, i)) < 0)
947 return 1;
948 }
949 for (i = dim; i < DIM(b); i++)
950 {
951 if (Max(LL_COORD(b, i), UR_COORD(b, i)) > 0)
952 return -1;
953 if (Max(LL_COORD(b, i), UR_COORD(b, i)) < 0)
954 return 1;
955 }
956
957 /*
958 * if all common dimensions are equal, the cube with more dimensions
959 * wins
960 */
961 return -1;
962 }
963
964 /* They're really equal */
965 return 0;
966 }
967
968 Datum
cube_cmp(PG_FUNCTION_ARGS)969 cube_cmp(PG_FUNCTION_ARGS)
970 {
971 NDBOX *a = PG_GETARG_NDBOX(0),
972 *b = PG_GETARG_NDBOX(1);
973 int32 res;
974
975 res = cube_cmp_v0(a, b);
976
977 PG_FREE_IF_COPY(a, 0);
978 PG_FREE_IF_COPY(b, 1);
979 PG_RETURN_INT32(res);
980 }
981
982
983 Datum
cube_eq(PG_FUNCTION_ARGS)984 cube_eq(PG_FUNCTION_ARGS)
985 {
986 NDBOX *a = PG_GETARG_NDBOX(0),
987 *b = PG_GETARG_NDBOX(1);
988 int32 res;
989
990 res = cube_cmp_v0(a, b);
991
992 PG_FREE_IF_COPY(a, 0);
993 PG_FREE_IF_COPY(b, 1);
994 PG_RETURN_BOOL(res == 0);
995 }
996
997
998 Datum
cube_ne(PG_FUNCTION_ARGS)999 cube_ne(PG_FUNCTION_ARGS)
1000 {
1001 NDBOX *a = PG_GETARG_NDBOX(0),
1002 *b = PG_GETARG_NDBOX(1);
1003 int32 res;
1004
1005 res = cube_cmp_v0(a, b);
1006
1007 PG_FREE_IF_COPY(a, 0);
1008 PG_FREE_IF_COPY(b, 1);
1009 PG_RETURN_BOOL(res != 0);
1010 }
1011
1012
1013 Datum
cube_lt(PG_FUNCTION_ARGS)1014 cube_lt(PG_FUNCTION_ARGS)
1015 {
1016 NDBOX *a = PG_GETARG_NDBOX(0),
1017 *b = PG_GETARG_NDBOX(1);
1018 int32 res;
1019
1020 res = cube_cmp_v0(a, b);
1021
1022 PG_FREE_IF_COPY(a, 0);
1023 PG_FREE_IF_COPY(b, 1);
1024 PG_RETURN_BOOL(res < 0);
1025 }
1026
1027
1028 Datum
cube_gt(PG_FUNCTION_ARGS)1029 cube_gt(PG_FUNCTION_ARGS)
1030 {
1031 NDBOX *a = PG_GETARG_NDBOX(0),
1032 *b = PG_GETARG_NDBOX(1);
1033 int32 res;
1034
1035 res = cube_cmp_v0(a, b);
1036
1037 PG_FREE_IF_COPY(a, 0);
1038 PG_FREE_IF_COPY(b, 1);
1039 PG_RETURN_BOOL(res > 0);
1040 }
1041
1042
1043 Datum
cube_le(PG_FUNCTION_ARGS)1044 cube_le(PG_FUNCTION_ARGS)
1045 {
1046 NDBOX *a = PG_GETARG_NDBOX(0),
1047 *b = PG_GETARG_NDBOX(1);
1048 int32 res;
1049
1050 res = cube_cmp_v0(a, b);
1051
1052 PG_FREE_IF_COPY(a, 0);
1053 PG_FREE_IF_COPY(b, 1);
1054 PG_RETURN_BOOL(res <= 0);
1055 }
1056
1057
1058 Datum
cube_ge(PG_FUNCTION_ARGS)1059 cube_ge(PG_FUNCTION_ARGS)
1060 {
1061 NDBOX *a = PG_GETARG_NDBOX(0),
1062 *b = PG_GETARG_NDBOX(1);
1063 int32 res;
1064
1065 res = cube_cmp_v0(a, b);
1066
1067 PG_FREE_IF_COPY(a, 0);
1068 PG_FREE_IF_COPY(b, 1);
1069 PG_RETURN_BOOL(res >= 0);
1070 }
1071
1072
1073 /* Contains */
1074 /* Box(A) CONTAINS Box(B) IFF pt(A) < pt(B) */
1075 bool
cube_contains_v0(NDBOX * a,NDBOX * b)1076 cube_contains_v0(NDBOX *a, NDBOX *b)
1077 {
1078 int i;
1079
1080 if ((a == NULL) || (b == NULL))
1081 return (FALSE);
1082
1083 if (DIM(a) < DIM(b))
1084 {
1085 /*
1086 * the further comparisons will make sense if the excess dimensions of
1087 * (b) were zeroes Since both UL and UR coordinates must be zero, we
1088 * can check them all without worrying about which is which.
1089 */
1090 for (i = DIM(a); i < DIM(b); i++)
1091 {
1092 if (LL_COORD(b, i) != 0)
1093 return (FALSE);
1094 if (UR_COORD(b, i) != 0)
1095 return (FALSE);
1096 }
1097 }
1098
1099 /* Can't care less about the excess dimensions of (a), if any */
1100 for (i = 0; i < Min(DIM(a), DIM(b)); i++)
1101 {
1102 if (Min(LL_COORD(a, i), UR_COORD(a, i)) >
1103 Min(LL_COORD(b, i), UR_COORD(b, i)))
1104 return (FALSE);
1105 if (Max(LL_COORD(a, i), UR_COORD(a, i)) <
1106 Max(LL_COORD(b, i), UR_COORD(b, i)))
1107 return (FALSE);
1108 }
1109
1110 return (TRUE);
1111 }
1112
1113 Datum
cube_contains(PG_FUNCTION_ARGS)1114 cube_contains(PG_FUNCTION_ARGS)
1115 {
1116 NDBOX *a = PG_GETARG_NDBOX(0),
1117 *b = PG_GETARG_NDBOX(1);
1118 bool res;
1119
1120 res = cube_contains_v0(a, b);
1121
1122 PG_FREE_IF_COPY(a, 0);
1123 PG_FREE_IF_COPY(b, 1);
1124 PG_RETURN_BOOL(res);
1125 }
1126
1127 /* Contained */
1128 /* Box(A) Contained by Box(B) IFF Box(B) Contains Box(A) */
1129 Datum
cube_contained(PG_FUNCTION_ARGS)1130 cube_contained(PG_FUNCTION_ARGS)
1131 {
1132 NDBOX *a = PG_GETARG_NDBOX(0),
1133 *b = PG_GETARG_NDBOX(1);
1134 bool res;
1135
1136 res = cube_contains_v0(b, a);
1137
1138 PG_FREE_IF_COPY(a, 0);
1139 PG_FREE_IF_COPY(b, 1);
1140 PG_RETURN_BOOL(res);
1141 }
1142
1143 /* Overlap */
1144 /* Box(A) Overlap Box(B) IFF (pt(a)LL < pt(B)UR) && (pt(b)LL < pt(a)UR) */
1145 bool
cube_overlap_v0(NDBOX * a,NDBOX * b)1146 cube_overlap_v0(NDBOX *a, NDBOX *b)
1147 {
1148 int i;
1149
1150 if ((a == NULL) || (b == NULL))
1151 return (FALSE);
1152
1153 /* swap the box pointers if needed */
1154 if (DIM(a) < DIM(b))
1155 {
1156 NDBOX *tmp = b;
1157
1158 b = a;
1159 a = tmp;
1160 }
1161
1162 /* compare within the dimensions of (b) */
1163 for (i = 0; i < DIM(b); i++)
1164 {
1165 if (Min(LL_COORD(a, i), UR_COORD(a, i)) > Max(LL_COORD(b, i), UR_COORD(b, i)))
1166 return (FALSE);
1167 if (Max(LL_COORD(a, i), UR_COORD(a, i)) < Min(LL_COORD(b, i), UR_COORD(b, i)))
1168 return (FALSE);
1169 }
1170
1171 /* compare to zero those dimensions in (a) absent in (b) */
1172 for (i = DIM(b); i < DIM(a); i++)
1173 {
1174 if (Min(LL_COORD(a, i), UR_COORD(a, i)) > 0)
1175 return (FALSE);
1176 if (Max(LL_COORD(a, i), UR_COORD(a, i)) < 0)
1177 return (FALSE);
1178 }
1179
1180 return (TRUE);
1181 }
1182
1183
1184 Datum
cube_overlap(PG_FUNCTION_ARGS)1185 cube_overlap(PG_FUNCTION_ARGS)
1186 {
1187 NDBOX *a = PG_GETARG_NDBOX(0),
1188 *b = PG_GETARG_NDBOX(1);
1189 bool res;
1190
1191 res = cube_overlap_v0(a, b);
1192
1193 PG_FREE_IF_COPY(a, 0);
1194 PG_FREE_IF_COPY(b, 1);
1195 PG_RETURN_BOOL(res);
1196 }
1197
1198
1199 /* Distance */
1200 /* The distance is computed as a per axis sum of the squared distances
1201 between 1D projections of the boxes onto Cartesian axes. Assuming zero
1202 distance between overlapping projections, this metric coincides with the
1203 "common sense" geometric distance */
1204 Datum
cube_distance(PG_FUNCTION_ARGS)1205 cube_distance(PG_FUNCTION_ARGS)
1206 {
1207 NDBOX *a = PG_GETARG_NDBOX(0),
1208 *b = PG_GETARG_NDBOX(1);
1209 bool swapped = false;
1210 double d,
1211 distance;
1212 int i;
1213
1214 /* swap the box pointers if needed */
1215 if (DIM(a) < DIM(b))
1216 {
1217 NDBOX *tmp = b;
1218
1219 b = a;
1220 a = tmp;
1221 swapped = true;
1222 }
1223
1224 distance = 0.0;
1225 /* compute within the dimensions of (b) */
1226 for (i = 0; i < DIM(b); i++)
1227 {
1228 d = distance_1D(LL_COORD(a, i), UR_COORD(a, i), LL_COORD(b, i), UR_COORD(b, i));
1229 distance += d * d;
1230 }
1231
1232 /* compute distance to zero for those dimensions in (a) absent in (b) */
1233 for (i = DIM(b); i < DIM(a); i++)
1234 {
1235 d = distance_1D(LL_COORD(a, i), UR_COORD(a, i), 0.0, 0.0);
1236 distance += d * d;
1237 }
1238
1239 if (swapped)
1240 {
1241 PG_FREE_IF_COPY(b, 0);
1242 PG_FREE_IF_COPY(a, 1);
1243 }
1244 else
1245 {
1246 PG_FREE_IF_COPY(a, 0);
1247 PG_FREE_IF_COPY(b, 1);
1248 }
1249
1250 PG_RETURN_FLOAT8(sqrt(distance));
1251 }
1252
1253 Datum
distance_taxicab(PG_FUNCTION_ARGS)1254 distance_taxicab(PG_FUNCTION_ARGS)
1255 {
1256 NDBOX *a = PG_GETARG_NDBOX(0),
1257 *b = PG_GETARG_NDBOX(1);
1258 bool swapped = false;
1259 double distance;
1260 int i;
1261
1262 /* swap the box pointers if needed */
1263 if (DIM(a) < DIM(b))
1264 {
1265 NDBOX *tmp = b;
1266
1267 b = a;
1268 a = tmp;
1269 swapped = true;
1270 }
1271
1272 distance = 0.0;
1273 /* compute within the dimensions of (b) */
1274 for (i = 0; i < DIM(b); i++)
1275 distance += fabs(distance_1D(LL_COORD(a, i), UR_COORD(a, i),
1276 LL_COORD(b, i), UR_COORD(b, i)));
1277
1278 /* compute distance to zero for those dimensions in (a) absent in (b) */
1279 for (i = DIM(b); i < DIM(a); i++)
1280 distance += fabs(distance_1D(LL_COORD(a, i), UR_COORD(a, i),
1281 0.0, 0.0));
1282
1283 if (swapped)
1284 {
1285 PG_FREE_IF_COPY(b, 0);
1286 PG_FREE_IF_COPY(a, 1);
1287 }
1288 else
1289 {
1290 PG_FREE_IF_COPY(a, 0);
1291 PG_FREE_IF_COPY(b, 1);
1292 }
1293
1294 PG_RETURN_FLOAT8(distance);
1295 }
1296
1297 Datum
distance_chebyshev(PG_FUNCTION_ARGS)1298 distance_chebyshev(PG_FUNCTION_ARGS)
1299 {
1300 NDBOX *a = PG_GETARG_NDBOX(0),
1301 *b = PG_GETARG_NDBOX(1);
1302 bool swapped = false;
1303 double d,
1304 distance;
1305 int i;
1306
1307 /* swap the box pointers if needed */
1308 if (DIM(a) < DIM(b))
1309 {
1310 NDBOX *tmp = b;
1311
1312 b = a;
1313 a = tmp;
1314 swapped = true;
1315 }
1316
1317 distance = 0.0;
1318 /* compute within the dimensions of (b) */
1319 for (i = 0; i < DIM(b); i++)
1320 {
1321 d = fabs(distance_1D(LL_COORD(a, i), UR_COORD(a, i),
1322 LL_COORD(b, i), UR_COORD(b, i)));
1323 if (d > distance)
1324 distance = d;
1325 }
1326
1327 /* compute distance to zero for those dimensions in (a) absent in (b) */
1328 for (i = DIM(b); i < DIM(a); i++)
1329 {
1330 d = fabs(distance_1D(LL_COORD(a, i), UR_COORD(a, i), 0.0, 0.0));
1331 if (d > distance)
1332 distance = d;
1333 }
1334
1335 if (swapped)
1336 {
1337 PG_FREE_IF_COPY(b, 0);
1338 PG_FREE_IF_COPY(a, 1);
1339 }
1340 else
1341 {
1342 PG_FREE_IF_COPY(a, 0);
1343 PG_FREE_IF_COPY(b, 1);
1344 }
1345
1346 PG_RETURN_FLOAT8(distance);
1347 }
1348
1349 Datum
g_cube_distance(PG_FUNCTION_ARGS)1350 g_cube_distance(PG_FUNCTION_ARGS)
1351 {
1352 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
1353 StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
1354 NDBOX *cube = DatumGetNDBOX(entry->key);
1355 double retval;
1356
1357 if (strategy == CubeKNNDistanceCoord)
1358 {
1359 /*
1360 * Handle ordering by ~> operator. See comments of cube_coord_llur()
1361 * for details
1362 */
1363 int coord = PG_GETARG_INT32(1);
1364 bool isLeaf = GistPageIsLeaf(entry->page);
1365
1366 /* 0 is the only unsupported coordinate value */
1367 if (coord <= 0)
1368 ereport(ERROR,
1369 (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
1370 errmsg("cube index %d is out of bounds", coord)));
1371
1372 if (coord <= 2 * DIM(cube))
1373 {
1374 /* dimension index */
1375 int index = (coord - 1) / 2;
1376 /* whether this is upper bound (lower bound otherwise) */
1377 bool upper = ((coord - 1) % 2 == 1);
1378
1379 if (IS_POINT(cube))
1380 {
1381 retval = cube->x[index];
1382 }
1383 else
1384 {
1385 if (isLeaf)
1386 {
1387 /* For leaf just return required upper/lower bound */
1388 if (upper)
1389 retval = Max(cube->x[index], cube->x[index + DIM(cube)]);
1390 else
1391 retval = Min(cube->x[index], cube->x[index + DIM(cube)]);
1392 }
1393 else
1394 {
1395 /*
1396 * For non-leaf we should always return lower bound,
1397 * because even upper bound of a child in the subtree can
1398 * be as small as our lower bound.
1399 */
1400 retval = Min(cube->x[index], cube->x[index + DIM(cube)]);
1401 }
1402 }
1403 }
1404 else
1405 {
1406 retval = 0.0;
1407 }
1408 }
1409 else
1410 {
1411 NDBOX *query = PG_GETARG_NDBOX(1);
1412
1413 switch (strategy)
1414 {
1415 case CubeKNNDistanceTaxicab:
1416 retval = DatumGetFloat8(DirectFunctionCall2(distance_taxicab,
1417 PointerGetDatum(cube), PointerGetDatum(query)));
1418 break;
1419 case CubeKNNDistanceEuclid:
1420 retval = DatumGetFloat8(DirectFunctionCall2(cube_distance,
1421 PointerGetDatum(cube), PointerGetDatum(query)));
1422 break;
1423 case CubeKNNDistanceChebyshev:
1424 retval = DatumGetFloat8(DirectFunctionCall2(distance_chebyshev,
1425 PointerGetDatum(cube), PointerGetDatum(query)));
1426 break;
1427 default:
1428 elog(ERROR, "unrecognized cube strategy number: %d", strategy);
1429 retval = 0; /* keep compiler quiet */
1430 break;
1431 }
1432 }
1433 PG_RETURN_FLOAT8(retval);
1434 }
1435
1436 static double
distance_1D(double a1,double a2,double b1,double b2)1437 distance_1D(double a1, double a2, double b1, double b2)
1438 {
1439 /* interval (a) is entirely on the left of (b) */
1440 if ((a1 <= b1) && (a2 <= b1) && (a1 <= b2) && (a2 <= b2))
1441 return (Min(b1, b2) - Max(a1, a2));
1442
1443 /* interval (a) is entirely on the right of (b) */
1444 if ((a1 > b1) && (a2 > b1) && (a1 > b2) && (a2 > b2))
1445 return (Min(a1, a2) - Max(b1, b2));
1446
1447 /* the rest are all sorts of intersections */
1448 return (0.0);
1449 }
1450
1451 /* Test if a box is also a point */
1452 Datum
cube_is_point(PG_FUNCTION_ARGS)1453 cube_is_point(PG_FUNCTION_ARGS)
1454 {
1455 NDBOX *cube = PG_GETARG_NDBOX(0);
1456 bool result;
1457
1458 result = cube_is_point_internal(cube);
1459 PG_FREE_IF_COPY(cube, 0);
1460 PG_RETURN_BOOL(result);
1461 }
1462
1463 static bool
cube_is_point_internal(NDBOX * cube)1464 cube_is_point_internal(NDBOX *cube)
1465 {
1466 int i;
1467
1468 if (IS_POINT(cube))
1469 return true;
1470
1471 /*
1472 * Even if the point-flag is not set, all the lower-left coordinates might
1473 * match the upper-right coordinates, so that the value is in fact a
1474 * point. Such values don't arise with current code - the point flag is
1475 * always set if appropriate - but they might be present on-disk in
1476 * clusters upgraded from pre-9.4 versions.
1477 */
1478 for (i = 0; i < DIM(cube); i++)
1479 {
1480 if (LL_COORD(cube, i) != UR_COORD(cube, i))
1481 return false;
1482 }
1483 return true;
1484 }
1485
1486 /* Return dimensions in use in the data structure */
1487 Datum
cube_dim(PG_FUNCTION_ARGS)1488 cube_dim(PG_FUNCTION_ARGS)
1489 {
1490 NDBOX *c = PG_GETARG_NDBOX(0);
1491 int dim = DIM(c);
1492
1493 PG_FREE_IF_COPY(c, 0);
1494 PG_RETURN_INT32(dim);
1495 }
1496
1497 /* Return a specific normalized LL coordinate */
1498 Datum
cube_ll_coord(PG_FUNCTION_ARGS)1499 cube_ll_coord(PG_FUNCTION_ARGS)
1500 {
1501 NDBOX *c = PG_GETARG_NDBOX(0);
1502 int n = PG_GETARG_INT32(1);
1503 double result;
1504
1505 if (DIM(c) >= n && n > 0)
1506 result = Min(LL_COORD(c, n - 1), UR_COORD(c, n - 1));
1507 else
1508 result = 0;
1509
1510 PG_FREE_IF_COPY(c, 0);
1511 PG_RETURN_FLOAT8(result);
1512 }
1513
1514 /* Return a specific normalized UR coordinate */
1515 Datum
cube_ur_coord(PG_FUNCTION_ARGS)1516 cube_ur_coord(PG_FUNCTION_ARGS)
1517 {
1518 NDBOX *c = PG_GETARG_NDBOX(0);
1519 int n = PG_GETARG_INT32(1);
1520 double result;
1521
1522 if (DIM(c) >= n && n > 0)
1523 result = Max(LL_COORD(c, n - 1), UR_COORD(c, n - 1));
1524 else
1525 result = 0;
1526
1527 PG_FREE_IF_COPY(c, 0);
1528 PG_RETURN_FLOAT8(result);
1529 }
1530
1531 /*
1532 * Function returns cube coordinate.
1533 * Numbers from 1 to DIM denotes first corner coordinates.
1534 * Numbers from DIM+1 to 2*DIM denotes second corner coordinates.
1535 */
1536 Datum
cube_coord(PG_FUNCTION_ARGS)1537 cube_coord(PG_FUNCTION_ARGS)
1538 {
1539 NDBOX *cube = PG_GETARG_NDBOX(0);
1540 int coord = PG_GETARG_INT32(1);
1541
1542 if (coord <= 0 || coord > 2 * DIM(cube))
1543 ereport(ERROR,
1544 (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
1545 errmsg("cube index %d is out of bounds", coord)));
1546
1547 if (IS_POINT(cube))
1548 PG_RETURN_FLOAT8(cube->x[(coord - 1) % DIM(cube)]);
1549 else
1550 PG_RETURN_FLOAT8(cube->x[coord - 1]);
1551 }
1552
1553
1554 /*----
1555 * This function works like cube_coord(), but rearranges coordinates in the
1556 * way suitable to support coordinate ordering using KNN-GiST. For historical
1557 * reasons this extension allows us to create cubes in form ((2,1),(1,2)) and
1558 * instead of normalizing such cube to ((1,1),(2,2)) it stores cube in original
1559 * way. But in order to get cubes ordered by one of dimensions from the index
1560 * without explicit sort step we need this representation-independent coordinate
1561 * getter. Moreover, indexed dataset may contain cubes of different dimensions
1562 * number. Accordingly, this coordinate getter should be able to return
1563 * lower/upper bound for particular dimension independently on number of cube
1564 * dimensions.
1565 *
1566 * Long story short, this function uses following meaning of coordinates:
1567 * # (2 * N - 1) -- lower bound of Nth dimension,
1568 * # (2 * N) -- upper bound of Nth dimension.
1569 *
1570 * When given coordinate exceeds number of cube dimensions, then 0 returned
1571 * (reproducing logic of GiST indexing of variable-length cubes).
1572 */
1573 Datum
cube_coord_llur(PG_FUNCTION_ARGS)1574 cube_coord_llur(PG_FUNCTION_ARGS)
1575 {
1576 NDBOX *cube = PG_GETARG_NDBOX(0);
1577 int coord = PG_GETARG_INT32(1);
1578 float8 result;
1579
1580 /* 0 is the only unsupported coordinate value */
1581 if (coord <= 0)
1582 ereport(ERROR,
1583 (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
1584 errmsg("cube index %d is out of bounds", coord)));
1585
1586 if (coord <= 2 * DIM(cube))
1587 {
1588 /* dimension index */
1589 int index = (coord - 1) / 2;
1590 /* whether this is upper bound (lower bound otherwise) */
1591 bool upper = ((coord - 1) % 2 == 1);
1592
1593 if (IS_POINT(cube))
1594 {
1595 result = cube->x[index];
1596 }
1597 else
1598 {
1599 if (upper)
1600 result = Max(cube->x[index], cube->x[index + DIM(cube)]);
1601 else
1602 result = Min(cube->x[index], cube->x[index + DIM(cube)]);
1603 }
1604 }
1605 else
1606 {
1607 /*
1608 * Return zero if coordinate is out of bound. That reproduces logic of
1609 * how cubes with low dimension number are expanded during GiST
1610 * indexing.
1611 */
1612 result = 0.0;
1613 }
1614
1615 PG_RETURN_FLOAT8(result);
1616 }
1617
1618 /* Increase or decrease box size by a radius in at least n dimensions. */
1619 Datum
cube_enlarge(PG_FUNCTION_ARGS)1620 cube_enlarge(PG_FUNCTION_ARGS)
1621 {
1622 NDBOX *a = PG_GETARG_NDBOX(0);
1623 double r = PG_GETARG_FLOAT8(1);
1624 int32 n = PG_GETARG_INT32(2);
1625 NDBOX *result;
1626 int dim = 0;
1627 int size;
1628 int i,
1629 j;
1630
1631 if (n > CUBE_MAX_DIM)
1632 n = CUBE_MAX_DIM;
1633 if (r > 0 && n > 0)
1634 dim = n;
1635 if (DIM(a) > dim)
1636 dim = DIM(a);
1637
1638 size = CUBE_SIZE(dim);
1639 result = (NDBOX *) palloc0(size);
1640 SET_VARSIZE(result, size);
1641 SET_DIM(result, dim);
1642
1643 for (i = 0, j = dim; i < DIM(a); i++, j++)
1644 {
1645 if (LL_COORD(a, i) >= UR_COORD(a, i))
1646 {
1647 result->x[i] = UR_COORD(a, i) - r;
1648 result->x[j] = LL_COORD(a, i) + r;
1649 }
1650 else
1651 {
1652 result->x[i] = LL_COORD(a, i) - r;
1653 result->x[j] = UR_COORD(a, i) + r;
1654 }
1655 if (result->x[i] > result->x[j])
1656 {
1657 result->x[i] = (result->x[i] + result->x[j]) / 2;
1658 result->x[j] = result->x[i];
1659 }
1660 }
1661 /* dim > a->dim only if r > 0 */
1662 for (; i < dim; i++, j++)
1663 {
1664 result->x[i] = -r;
1665 result->x[j] = r;
1666 }
1667
1668 /*
1669 * Check if the result was in fact a point, and set the flag in the datum
1670 * accordingly. (we don't bother to repalloc it smaller)
1671 */
1672 if (cube_is_point_internal(result))
1673 {
1674 size = POINT_SIZE(dim);
1675 SET_VARSIZE(result, size);
1676 SET_POINT_BIT(result);
1677 }
1678
1679 PG_FREE_IF_COPY(a, 0);
1680 PG_RETURN_NDBOX(result);
1681 }
1682
1683 /* Create a one dimensional box with identical upper and lower coordinates */
1684 Datum
cube_f8(PG_FUNCTION_ARGS)1685 cube_f8(PG_FUNCTION_ARGS)
1686 {
1687 double x = PG_GETARG_FLOAT8(0);
1688 NDBOX *result;
1689 int size;
1690
1691 size = POINT_SIZE(1);
1692 result = (NDBOX *) palloc0(size);
1693 SET_VARSIZE(result, size);
1694 SET_DIM(result, 1);
1695 SET_POINT_BIT(result);
1696 result->x[0] = x;
1697
1698 PG_RETURN_NDBOX(result);
1699 }
1700
1701 /* Create a one dimensional box */
1702 Datum
cube_f8_f8(PG_FUNCTION_ARGS)1703 cube_f8_f8(PG_FUNCTION_ARGS)
1704 {
1705 double x0 = PG_GETARG_FLOAT8(0);
1706 double x1 = PG_GETARG_FLOAT8(1);
1707 NDBOX *result;
1708 int size;
1709
1710 if (x0 == x1)
1711 {
1712 size = POINT_SIZE(1);
1713 result = (NDBOX *) palloc0(size);
1714 SET_VARSIZE(result, size);
1715 SET_DIM(result, 1);
1716 SET_POINT_BIT(result);
1717 result->x[0] = x0;
1718 }
1719 else
1720 {
1721 size = CUBE_SIZE(1);
1722 result = (NDBOX *) palloc0(size);
1723 SET_VARSIZE(result, size);
1724 SET_DIM(result, 1);
1725 result->x[0] = x0;
1726 result->x[1] = x1;
1727 }
1728
1729 PG_RETURN_NDBOX(result);
1730 }
1731
1732 /* Add a dimension to an existing cube with the same values for the new
1733 coordinate */
1734 Datum
cube_c_f8(PG_FUNCTION_ARGS)1735 cube_c_f8(PG_FUNCTION_ARGS)
1736 {
1737 NDBOX *cube = PG_GETARG_NDBOX(0);
1738 double x = PG_GETARG_FLOAT8(1);
1739 NDBOX *result;
1740 int size;
1741 int i;
1742
1743 if (DIM(cube) + 1 > CUBE_MAX_DIM)
1744 ereport(ERROR,
1745 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
1746 errmsg("can't extend cube"),
1747 errdetail("A cube cannot have more than %d dimensions.",
1748 CUBE_MAX_DIM)));
1749
1750 if (IS_POINT(cube))
1751 {
1752 size = POINT_SIZE((DIM(cube) + 1));
1753 result = (NDBOX *) palloc0(size);
1754 SET_VARSIZE(result, size);
1755 SET_DIM(result, DIM(cube) + 1);
1756 SET_POINT_BIT(result);
1757 for (i = 0; i < DIM(cube); i++)
1758 result->x[i] = cube->x[i];
1759 result->x[DIM(result) - 1] = x;
1760 }
1761 else
1762 {
1763 size = CUBE_SIZE((DIM(cube) + 1));
1764 result = (NDBOX *) palloc0(size);
1765 SET_VARSIZE(result, size);
1766 SET_DIM(result, DIM(cube) + 1);
1767 for (i = 0; i < DIM(cube); i++)
1768 {
1769 result->x[i] = cube->x[i];
1770 result->x[DIM(result) + i] = cube->x[DIM(cube) + i];
1771 }
1772 result->x[DIM(result) - 1] = x;
1773 result->x[2 * DIM(result) - 1] = x;
1774 }
1775
1776 PG_FREE_IF_COPY(cube, 0);
1777 PG_RETURN_NDBOX(result);
1778 }
1779
1780 /* Add a dimension to an existing cube */
1781 Datum
cube_c_f8_f8(PG_FUNCTION_ARGS)1782 cube_c_f8_f8(PG_FUNCTION_ARGS)
1783 {
1784 NDBOX *cube = PG_GETARG_NDBOX(0);
1785 double x1 = PG_GETARG_FLOAT8(1);
1786 double x2 = PG_GETARG_FLOAT8(2);
1787 NDBOX *result;
1788 int size;
1789 int i;
1790
1791 if (DIM(cube) + 1 > CUBE_MAX_DIM)
1792 ereport(ERROR,
1793 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
1794 errmsg("can't extend cube"),
1795 errdetail("A cube cannot have more than %d dimensions.",
1796 CUBE_MAX_DIM)));
1797
1798 if (IS_POINT(cube) && (x1 == x2))
1799 {
1800 size = POINT_SIZE((DIM(cube) + 1));
1801 result = (NDBOX *) palloc0(size);
1802 SET_VARSIZE(result, size);
1803 SET_DIM(result, DIM(cube) + 1);
1804 SET_POINT_BIT(result);
1805 for (i = 0; i < DIM(cube); i++)
1806 result->x[i] = cube->x[i];
1807 result->x[DIM(result) - 1] = x1;
1808 }
1809 else
1810 {
1811 size = CUBE_SIZE((DIM(cube) + 1));
1812 result = (NDBOX *) palloc0(size);
1813 SET_VARSIZE(result, size);
1814 SET_DIM(result, DIM(cube) + 1);
1815 for (i = 0; i < DIM(cube); i++)
1816 {
1817 result->x[i] = LL_COORD(cube, i);
1818 result->x[DIM(result) + i] = UR_COORD(cube, i);
1819 }
1820 result->x[DIM(result) - 1] = x1;
1821 result->x[2 * DIM(result) - 1] = x2;
1822 }
1823
1824 PG_FREE_IF_COPY(cube, 0);
1825 PG_RETURN_NDBOX(result);
1826 }
1827