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