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 "utils/array.h"
16 #include "utils/float.h"
17
18 #include "cubedata.h"
19
20 PG_MODULE_MAGIC;
21
22 /*
23 * Taken from the intarray contrib header
24 */
25 #define ARRPTR(x) ( (double *) ARR_DATA_PTR(x) )
26 #define ARRNELEMS(x) ArrayGetNItems( ARR_NDIM(x), ARR_DIMS(x))
27
28 /*
29 ** Input/Output routines
30 */
31 PG_FUNCTION_INFO_V1(cube_in);
32 PG_FUNCTION_INFO_V1(cube_a_f8_f8);
33 PG_FUNCTION_INFO_V1(cube_a_f8);
34 PG_FUNCTION_INFO_V1(cube_out);
35 PG_FUNCTION_INFO_V1(cube_f8);
36 PG_FUNCTION_INFO_V1(cube_f8_f8);
37 PG_FUNCTION_INFO_V1(cube_c_f8);
38 PG_FUNCTION_INFO_V1(cube_c_f8_f8);
39 PG_FUNCTION_INFO_V1(cube_dim);
40 PG_FUNCTION_INFO_V1(cube_ll_coord);
41 PG_FUNCTION_INFO_V1(cube_ur_coord);
42 PG_FUNCTION_INFO_V1(cube_coord);
43 PG_FUNCTION_INFO_V1(cube_coord_llur);
44 PG_FUNCTION_INFO_V1(cube_subset);
45
46 /*
47 ** GiST support methods
48 */
49
50 PG_FUNCTION_INFO_V1(g_cube_consistent);
51 PG_FUNCTION_INFO_V1(g_cube_compress);
52 PG_FUNCTION_INFO_V1(g_cube_decompress);
53 PG_FUNCTION_INFO_V1(g_cube_penalty);
54 PG_FUNCTION_INFO_V1(g_cube_picksplit);
55 PG_FUNCTION_INFO_V1(g_cube_union);
56 PG_FUNCTION_INFO_V1(g_cube_same);
57 PG_FUNCTION_INFO_V1(g_cube_distance);
58
59 /*
60 ** B-tree support functions
61 */
62 PG_FUNCTION_INFO_V1(cube_eq);
63 PG_FUNCTION_INFO_V1(cube_ne);
64 PG_FUNCTION_INFO_V1(cube_lt);
65 PG_FUNCTION_INFO_V1(cube_gt);
66 PG_FUNCTION_INFO_V1(cube_le);
67 PG_FUNCTION_INFO_V1(cube_ge);
68 PG_FUNCTION_INFO_V1(cube_cmp);
69
70 /*
71 ** R-tree support functions
72 */
73
74 PG_FUNCTION_INFO_V1(cube_contains);
75 PG_FUNCTION_INFO_V1(cube_contained);
76 PG_FUNCTION_INFO_V1(cube_overlap);
77 PG_FUNCTION_INFO_V1(cube_union);
78 PG_FUNCTION_INFO_V1(cube_inter);
79 PG_FUNCTION_INFO_V1(cube_size);
80
81 /*
82 ** miscellaneous
83 */
84 PG_FUNCTION_INFO_V1(distance_taxicab);
85 PG_FUNCTION_INFO_V1(cube_distance);
86 PG_FUNCTION_INFO_V1(distance_chebyshev);
87 PG_FUNCTION_INFO_V1(cube_is_point);
88 PG_FUNCTION_INFO_V1(cube_enlarge);
89
90 /*
91 ** For internal use only
92 */
93 int32 cube_cmp_v0(NDBOX *a, NDBOX *b);
94 bool cube_contains_v0(NDBOX *a, NDBOX *b);
95 bool cube_overlap_v0(NDBOX *a, NDBOX *b);
96 NDBOX *cube_union_v0(NDBOX *a, NDBOX *b);
97 void rt_cube_size(NDBOX *a, double *sz);
98 NDBOX *g_cube_binary_union(NDBOX *r1, NDBOX *r2, int *sizep);
99 bool g_cube_leaf_consistent(NDBOX *key, NDBOX *query, StrategyNumber strategy);
100 bool g_cube_internal_consistent(NDBOX *key, NDBOX *query, StrategyNumber strategy);
101
102 /*
103 ** Auxiliary funxtions
104 */
105 static double distance_1D(double a1, double a2, double b1, double b2);
106 static bool cube_is_point_internal(NDBOX *cube);
107
108
109 /*****************************************************************************
110 * Input/Output functions
111 *****************************************************************************/
112
113 /* NdBox = [(lowerleft),(upperright)] */
114 /* [(xLL(1)...xLL(N)),(xUR(1)...xUR(n))] */
115 Datum
cube_in(PG_FUNCTION_ARGS)116 cube_in(PG_FUNCTION_ARGS)
117 {
118 char *str = PG_GETARG_CSTRING(0);
119 NDBOX *result;
120
121 cube_scanner_init(str);
122
123 if (cube_yyparse(&result) != 0)
124 cube_yyerror(&result, "cube parser failed");
125
126 cube_scanner_finish();
127
128 PG_RETURN_NDBOX_P(result);
129 }
130
131
132 /*
133 ** Allows the construction of a cube from 2 float[]'s
134 */
135 Datum
cube_a_f8_f8(PG_FUNCTION_ARGS)136 cube_a_f8_f8(PG_FUNCTION_ARGS)
137 {
138 ArrayType *ur = PG_GETARG_ARRAYTYPE_P(0);
139 ArrayType *ll = PG_GETARG_ARRAYTYPE_P(1);
140 NDBOX *result;
141 int i;
142 int dim;
143 int size;
144 bool point;
145 double *dur,
146 *dll;
147
148 if (array_contains_nulls(ur) || array_contains_nulls(ll))
149 ereport(ERROR,
150 (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
151 errmsg("cannot work with arrays containing NULLs")));
152
153 dim = ARRNELEMS(ur);
154 if (dim > CUBE_MAX_DIM)
155 ereport(ERROR,
156 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
157 errmsg("can't extend cube"),
158 errdetail("A cube cannot have more than %d dimensions.",
159 CUBE_MAX_DIM)));
160
161 if (ARRNELEMS(ll) != dim)
162 ereport(ERROR,
163 (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
164 errmsg("UR and LL arrays must be of same length")));
165
166 dur = ARRPTR(ur);
167 dll = ARRPTR(ll);
168
169 /* Check if it's a point */
170 point = true;
171 for (i = 0; i < dim; i++)
172 {
173 if (dur[i] != dll[i])
174 {
175 point = false;
176 break;
177 }
178 }
179
180 size = point ? POINT_SIZE(dim) : CUBE_SIZE(dim);
181 result = (NDBOX *) palloc0(size);
182 SET_VARSIZE(result, size);
183 SET_DIM(result, dim);
184
185 for (i = 0; i < dim; i++)
186 result->x[i] = dur[i];
187
188 if (!point)
189 {
190 for (i = 0; i < dim; i++)
191 result->x[i + dim] = dll[i];
192 }
193 else
194 SET_POINT_BIT(result);
195
196 PG_RETURN_NDBOX_P(result);
197 }
198
199 /*
200 ** Allows the construction of a zero-volume cube from a float[]
201 */
202 Datum
cube_a_f8(PG_FUNCTION_ARGS)203 cube_a_f8(PG_FUNCTION_ARGS)
204 {
205 ArrayType *ur = PG_GETARG_ARRAYTYPE_P(0);
206 NDBOX *result;
207 int i;
208 int dim;
209 int size;
210 double *dur;
211
212 if (array_contains_nulls(ur))
213 ereport(ERROR,
214 (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
215 errmsg("cannot work with arrays containing NULLs")));
216
217 dim = ARRNELEMS(ur);
218 if (dim > CUBE_MAX_DIM)
219 ereport(ERROR,
220 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
221 errmsg("array is too long"),
222 errdetail("A cube cannot have more than %d dimensions.",
223 CUBE_MAX_DIM)));
224
225 dur = ARRPTR(ur);
226
227 size = POINT_SIZE(dim);
228 result = (NDBOX *) palloc0(size);
229 SET_VARSIZE(result, size);
230 SET_DIM(result, dim);
231 SET_POINT_BIT(result);
232
233 for (i = 0; i < dim; i++)
234 result->x[i] = dur[i];
235
236 PG_RETURN_NDBOX_P(result);
237 }
238
239 Datum
cube_subset(PG_FUNCTION_ARGS)240 cube_subset(PG_FUNCTION_ARGS)
241 {
242 NDBOX *c = PG_GETARG_NDBOX_P(0);
243 ArrayType *idx = PG_GETARG_ARRAYTYPE_P(1);
244 NDBOX *result;
245 int size,
246 dim,
247 i;
248 int *dx;
249
250 if (array_contains_nulls(idx))
251 ereport(ERROR,
252 (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
253 errmsg("cannot work with arrays containing NULLs")));
254
255 dx = (int32 *) ARR_DATA_PTR(idx);
256
257 dim = ARRNELEMS(idx);
258 if (dim > CUBE_MAX_DIM)
259 ereport(ERROR,
260 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
261 errmsg("array is too long"),
262 errdetail("A cube cannot have more than %d dimensions.",
263 CUBE_MAX_DIM)));
264
265 size = IS_POINT(c) ? POINT_SIZE(dim) : CUBE_SIZE(dim);
266 result = (NDBOX *) palloc0(size);
267 SET_VARSIZE(result, size);
268 SET_DIM(result, dim);
269
270 if (IS_POINT(c))
271 SET_POINT_BIT(result);
272
273 for (i = 0; i < dim; i++)
274 {
275 if ((dx[i] <= 0) || (dx[i] > DIM(c)))
276 ereport(ERROR,
277 (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
278 errmsg("Index out of bounds")));
279 result->x[i] = c->x[dx[i] - 1];
280 if (!IS_POINT(c))
281 result->x[i + dim] = c->x[dx[i] + DIM(c) - 1];
282 }
283
284 PG_FREE_IF_COPY(c, 0);
285 PG_RETURN_NDBOX_P(result);
286 }
287
288 Datum
cube_out(PG_FUNCTION_ARGS)289 cube_out(PG_FUNCTION_ARGS)
290 {
291 NDBOX *cube = PG_GETARG_NDBOX_P(0);
292 StringInfoData buf;
293 int dim = DIM(cube);
294 int i;
295
296 initStringInfo(&buf);
297
298 appendStringInfoChar(&buf, '(');
299 for (i = 0; i < dim; i++)
300 {
301 if (i > 0)
302 appendStringInfoString(&buf, ", ");
303 appendStringInfoString(&buf, float8out_internal(LL_COORD(cube, i)));
304 }
305 appendStringInfoChar(&buf, ')');
306
307 if (!cube_is_point_internal(cube))
308 {
309 appendStringInfoString(&buf, ",(");
310 for (i = 0; i < dim; i++)
311 {
312 if (i > 0)
313 appendStringInfoString(&buf, ", ");
314 appendStringInfoString(&buf, float8out_internal(UR_COORD(cube, i)));
315 }
316 appendStringInfoChar(&buf, ')');
317 }
318
319 PG_FREE_IF_COPY(cube, 0);
320 PG_RETURN_CSTRING(buf.data);
321 }
322
323
324 /*****************************************************************************
325 * GiST functions
326 *****************************************************************************/
327
328 /*
329 ** The GiST Consistent method for boxes
330 ** Should return false if for all data items x below entry,
331 ** the predicate x op query == false, where op is the oper
332 ** corresponding to strategy in the pg_amop table.
333 */
334 Datum
g_cube_consistent(PG_FUNCTION_ARGS)335 g_cube_consistent(PG_FUNCTION_ARGS)
336 {
337 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
338 NDBOX *query = PG_GETARG_NDBOX_P(1);
339 StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
340
341 /* Oid subtype = PG_GETARG_OID(3); */
342 bool *recheck = (bool *) PG_GETARG_POINTER(4);
343 bool res;
344
345 /* All cases served by this function are exact */
346 *recheck = false;
347
348 /*
349 * if entry is not leaf, use g_cube_internal_consistent, else use
350 * g_cube_leaf_consistent
351 */
352 if (GIST_LEAF(entry))
353 res = g_cube_leaf_consistent(DatumGetNDBOXP(entry->key),
354 query, strategy);
355 else
356 res = g_cube_internal_consistent(DatumGetNDBOXP(entry->key),
357 query, strategy);
358
359 PG_FREE_IF_COPY(query, 1);
360 PG_RETURN_BOOL(res);
361 }
362
363
364 /*
365 ** The GiST Union method for boxes
366 ** returns the minimal bounding box that encloses all the entries in entryvec
367 */
368 Datum
g_cube_union(PG_FUNCTION_ARGS)369 g_cube_union(PG_FUNCTION_ARGS)
370 {
371 GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
372 int *sizep = (int *) PG_GETARG_POINTER(1);
373 NDBOX *out = (NDBOX *) NULL;
374 NDBOX *tmp;
375 int i;
376
377 tmp = DatumGetNDBOXP(entryvec->vector[0].key);
378
379 /*
380 * sizep = sizeof(NDBOX); -- NDBOX has variable size
381 */
382 *sizep = VARSIZE(tmp);
383
384 for (i = 1; i < entryvec->n; i++)
385 {
386 out = g_cube_binary_union(tmp,
387 DatumGetNDBOXP(entryvec->vector[i].key),
388 sizep);
389 tmp = out;
390 }
391
392 PG_RETURN_POINTER(out);
393 }
394
395 /*
396 ** GiST Compress and Decompress methods for boxes
397 ** do not do anything.
398 */
399
400 Datum
g_cube_compress(PG_FUNCTION_ARGS)401 g_cube_compress(PG_FUNCTION_ARGS)
402 {
403 PG_RETURN_DATUM(PG_GETARG_DATUM(0));
404 }
405
406 Datum
g_cube_decompress(PG_FUNCTION_ARGS)407 g_cube_decompress(PG_FUNCTION_ARGS)
408 {
409 GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
410 NDBOX *key = DatumGetNDBOXP(entry->key);
411
412 if (key != DatumGetNDBOXP(entry->key))
413 {
414 GISTENTRY *retval = (GISTENTRY *) palloc(sizeof(GISTENTRY));
415
416 gistentryinit(*retval, PointerGetDatum(key),
417 entry->rel, entry->page,
418 entry->offset, false);
419 PG_RETURN_POINTER(retval);
420 }
421 PG_RETURN_POINTER(entry);
422 }
423
424
425 /*
426 ** The GiST Penalty method for boxes
427 ** As in the R-tree paper, we use change in area as our penalty metric
428 */
429 Datum
g_cube_penalty(PG_FUNCTION_ARGS)430 g_cube_penalty(PG_FUNCTION_ARGS)
431 {
432 GISTENTRY *origentry = (GISTENTRY *) PG_GETARG_POINTER(0);
433 GISTENTRY *newentry = (GISTENTRY *) PG_GETARG_POINTER(1);
434 float *result = (float *) PG_GETARG_POINTER(2);
435 NDBOX *ud;
436 double tmp1,
437 tmp2;
438
439 ud = cube_union_v0(DatumGetNDBOXP(origentry->key),
440 DatumGetNDBOXP(newentry->key));
441 rt_cube_size(ud, &tmp1);
442 rt_cube_size(DatumGetNDBOXP(origentry->key), &tmp2);
443 *result = (float) (tmp1 - tmp2);
444
445 PG_RETURN_FLOAT8(*result);
446 }
447
448
449
450 /*
451 ** The GiST PickSplit method for boxes
452 ** We use Guttman's poly time split algorithm
453 */
454 Datum
g_cube_picksplit(PG_FUNCTION_ARGS)455 g_cube_picksplit(PG_FUNCTION_ARGS)
456 {
457 GistEntryVector *entryvec = (GistEntryVector *) PG_GETARG_POINTER(0);
458 GIST_SPLITVEC *v = (GIST_SPLITVEC *) PG_GETARG_POINTER(1);
459 OffsetNumber i,
460 j;
461 NDBOX *datum_alpha,
462 *datum_beta;
463 NDBOX *datum_l,
464 *datum_r;
465 NDBOX *union_d,
466 *union_dl,
467 *union_dr;
468 NDBOX *inter_d;
469 bool firsttime;
470 double size_alpha,
471 size_beta,
472 size_union,
473 size_inter;
474 double size_waste,
475 waste;
476 double size_l,
477 size_r;
478 int nbytes;
479 OffsetNumber seed_1 = 1,
480 seed_2 = 2;
481 OffsetNumber *left,
482 *right;
483 OffsetNumber maxoff;
484
485 maxoff = entryvec->n - 2;
486 nbytes = (maxoff + 2) * sizeof(OffsetNumber);
487 v->spl_left = (OffsetNumber *) palloc(nbytes);
488 v->spl_right = (OffsetNumber *) palloc(nbytes);
489
490 firsttime = true;
491 waste = 0.0;
492
493 for (i = FirstOffsetNumber; i < maxoff; i = OffsetNumberNext(i))
494 {
495 datum_alpha = DatumGetNDBOXP(entryvec->vector[i].key);
496 for (j = OffsetNumberNext(i); j <= maxoff; j = OffsetNumberNext(j))
497 {
498 datum_beta = DatumGetNDBOXP(entryvec->vector[j].key);
499
500 /* compute the wasted space by unioning these guys */
501 /* size_waste = size_union - size_inter; */
502 union_d = cube_union_v0(datum_alpha, datum_beta);
503 rt_cube_size(union_d, &size_union);
504 inter_d = DatumGetNDBOXP(DirectFunctionCall2(cube_inter,
505 entryvec->vector[i].key,
506 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 = DatumGetNDBOXP(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 = DatumGetNDBOXP(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 = DatumGetNDBOXP(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_P(0);
608 NDBOX *b2 = PG_GETARG_NDBOX_P(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_P(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 = cube_overlap_v0(key, query);
633 break;
634 case RTSameStrategyNumber:
635 retval = (cube_cmp_v0(key, query) == 0);
636 break;
637 case RTContainsStrategyNumber:
638 case RTOldContainsStrategyNumber:
639 retval = cube_contains_v0(key, query);
640 break;
641 case RTContainedByStrategyNumber:
642 case RTOldContainedByStrategyNumber:
643 retval = 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_P(0);
759 NDBOX *b = PG_GETARG_NDBOX_P(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_P(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_P(0);
774 NDBOX *b = PG_GETARG_NDBOX_P(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_P(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_P(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_P(0),
972 *b = PG_GETARG_NDBOX_P(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_P(0),
987 *b = PG_GETARG_NDBOX_P(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_P(0),
1002 *b = PG_GETARG_NDBOX_P(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_P(0),
1017 *b = PG_GETARG_NDBOX_P(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_P(0),
1032 *b = PG_GETARG_NDBOX_P(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_P(0),
1047 *b = PG_GETARG_NDBOX_P(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_P(0),
1062 *b = PG_GETARG_NDBOX_P(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_P(0),
1117 *b = PG_GETARG_NDBOX_P(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_P(0),
1133 *b = PG_GETARG_NDBOX_P(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_P(0),
1188 *b = PG_GETARG_NDBOX_P(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_P(0),
1208 *b = PG_GETARG_NDBOX_P(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_P(0),
1257 *b = PG_GETARG_NDBOX_P(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_P(0),
1301 *b = PG_GETARG_NDBOX_P(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 = DatumGetNDBOXP(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 bool inverse = false;
1366
1367 /* 0 is the only unsupported coordinate value */
1368 if (coord == 0)
1369 ereport(ERROR,
1370 (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
1371 errmsg("zero cube index is not defined")));
1372
1373 /* Return inversed value for negative coordinate */
1374 if (coord < 0)
1375 {
1376 coord = -coord;
1377 inverse = true;
1378 }
1379
1380 if (coord <= 2 * DIM(cube))
1381 {
1382 /* dimension index */
1383 int index = (coord - 1) / 2;
1384
1385 /* whether this is upper bound (lower bound otherwise) */
1386 bool upper = ((coord - 1) % 2 == 1);
1387
1388 if (IS_POINT(cube))
1389 {
1390 retval = cube->x[index];
1391 }
1392 else
1393 {
1394 if (isLeaf)
1395 {
1396 /* For leaf just return required upper/lower bound */
1397 if (upper)
1398 retval = Max(cube->x[index], cube->x[index + DIM(cube)]);
1399 else
1400 retval = Min(cube->x[index], cube->x[index + DIM(cube)]);
1401 }
1402 else
1403 {
1404 /*
1405 * For non-leaf we should always return lower bound,
1406 * because even upper bound of a child in the subtree can
1407 * be as small as our lower bound. For inversed case we
1408 * return upper bound because it becomes lower bound for
1409 * inversed value.
1410 */
1411 if (!inverse)
1412 retval = Min(cube->x[index], cube->x[index + DIM(cube)]);
1413 else
1414 retval = Max(cube->x[index], cube->x[index + DIM(cube)]);
1415 }
1416 }
1417 }
1418 else
1419 {
1420 retval = 0.0;
1421 }
1422
1423 /* Inverse return value if needed */
1424 if (inverse)
1425 retval = -retval;
1426 }
1427 else
1428 {
1429 NDBOX *query = PG_GETARG_NDBOX_P(1);
1430
1431 switch (strategy)
1432 {
1433 case CubeKNNDistanceTaxicab:
1434 retval = DatumGetFloat8(DirectFunctionCall2(distance_taxicab,
1435 PointerGetDatum(cube), PointerGetDatum(query)));
1436 break;
1437 case CubeKNNDistanceEuclid:
1438 retval = DatumGetFloat8(DirectFunctionCall2(cube_distance,
1439 PointerGetDatum(cube), PointerGetDatum(query)));
1440 break;
1441 case CubeKNNDistanceChebyshev:
1442 retval = DatumGetFloat8(DirectFunctionCall2(distance_chebyshev,
1443 PointerGetDatum(cube), PointerGetDatum(query)));
1444 break;
1445 default:
1446 elog(ERROR, "unrecognized cube strategy number: %d", strategy);
1447 retval = 0; /* keep compiler quiet */
1448 break;
1449 }
1450 }
1451 PG_RETURN_FLOAT8(retval);
1452 }
1453
1454 static double
distance_1D(double a1,double a2,double b1,double b2)1455 distance_1D(double a1, double a2, double b1, double b2)
1456 {
1457 /* interval (a) is entirely on the left of (b) */
1458 if ((a1 <= b1) && (a2 <= b1) && (a1 <= b2) && (a2 <= b2))
1459 return (Min(b1, b2) - Max(a1, a2));
1460
1461 /* interval (a) is entirely on the right of (b) */
1462 if ((a1 > b1) && (a2 > b1) && (a1 > b2) && (a2 > b2))
1463 return (Min(a1, a2) - Max(b1, b2));
1464
1465 /* the rest are all sorts of intersections */
1466 return 0.0;
1467 }
1468
1469 /* Test if a box is also a point */
1470 Datum
cube_is_point(PG_FUNCTION_ARGS)1471 cube_is_point(PG_FUNCTION_ARGS)
1472 {
1473 NDBOX *cube = PG_GETARG_NDBOX_P(0);
1474 bool result;
1475
1476 result = cube_is_point_internal(cube);
1477 PG_FREE_IF_COPY(cube, 0);
1478 PG_RETURN_BOOL(result);
1479 }
1480
1481 static bool
cube_is_point_internal(NDBOX * cube)1482 cube_is_point_internal(NDBOX *cube)
1483 {
1484 int i;
1485
1486 if (IS_POINT(cube))
1487 return true;
1488
1489 /*
1490 * Even if the point-flag is not set, all the lower-left coordinates might
1491 * match the upper-right coordinates, so that the value is in fact a
1492 * point. Such values don't arise with current code - the point flag is
1493 * always set if appropriate - but they might be present on-disk in
1494 * clusters upgraded from pre-9.4 versions.
1495 */
1496 for (i = 0; i < DIM(cube); i++)
1497 {
1498 if (LL_COORD(cube, i) != UR_COORD(cube, i))
1499 return false;
1500 }
1501 return true;
1502 }
1503
1504 /* Return dimensions in use in the data structure */
1505 Datum
cube_dim(PG_FUNCTION_ARGS)1506 cube_dim(PG_FUNCTION_ARGS)
1507 {
1508 NDBOX *c = PG_GETARG_NDBOX_P(0);
1509 int dim = DIM(c);
1510
1511 PG_FREE_IF_COPY(c, 0);
1512 PG_RETURN_INT32(dim);
1513 }
1514
1515 /* Return a specific normalized LL coordinate */
1516 Datum
cube_ll_coord(PG_FUNCTION_ARGS)1517 cube_ll_coord(PG_FUNCTION_ARGS)
1518 {
1519 NDBOX *c = PG_GETARG_NDBOX_P(0);
1520 int n = PG_GETARG_INT32(1);
1521 double result;
1522
1523 if (DIM(c) >= n && n > 0)
1524 result = Min(LL_COORD(c, n - 1), UR_COORD(c, n - 1));
1525 else
1526 result = 0;
1527
1528 PG_FREE_IF_COPY(c, 0);
1529 PG_RETURN_FLOAT8(result);
1530 }
1531
1532 /* Return a specific normalized UR coordinate */
1533 Datum
cube_ur_coord(PG_FUNCTION_ARGS)1534 cube_ur_coord(PG_FUNCTION_ARGS)
1535 {
1536 NDBOX *c = PG_GETARG_NDBOX_P(0);
1537 int n = PG_GETARG_INT32(1);
1538 double result;
1539
1540 if (DIM(c) >= n && n > 0)
1541 result = Max(LL_COORD(c, n - 1), UR_COORD(c, n - 1));
1542 else
1543 result = 0;
1544
1545 PG_FREE_IF_COPY(c, 0);
1546 PG_RETURN_FLOAT8(result);
1547 }
1548
1549 /*
1550 * Function returns cube coordinate.
1551 * Numbers from 1 to DIM denotes first corner coordinates.
1552 * Numbers from DIM+1 to 2*DIM denotes second corner coordinates.
1553 */
1554 Datum
cube_coord(PG_FUNCTION_ARGS)1555 cube_coord(PG_FUNCTION_ARGS)
1556 {
1557 NDBOX *cube = PG_GETARG_NDBOX_P(0);
1558 int coord = PG_GETARG_INT32(1);
1559
1560 if (coord <= 0 || coord > 2 * DIM(cube))
1561 ereport(ERROR,
1562 (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
1563 errmsg("cube index %d is out of bounds", coord)));
1564
1565 if (IS_POINT(cube))
1566 PG_RETURN_FLOAT8(cube->x[(coord - 1) % DIM(cube)]);
1567 else
1568 PG_RETURN_FLOAT8(cube->x[coord - 1]);
1569 }
1570
1571
1572 /*----
1573 * This function works like cube_coord(), but rearranges coordinates in the
1574 * way suitable to support coordinate ordering using KNN-GiST. For historical
1575 * reasons this extension allows us to create cubes in form ((2,1),(1,2)) and
1576 * instead of normalizing such cube to ((1,1),(2,2)) it stores cube in original
1577 * way. But in order to get cubes ordered by one of dimensions from the index
1578 * without explicit sort step we need this representation-independent coordinate
1579 * getter. Moreover, indexed dataset may contain cubes of different dimensions
1580 * number. Accordingly, this coordinate getter should be able to return
1581 * lower/upper bound for particular dimension independently on number of cube
1582 * dimensions. Also, KNN-GiST supports only ascending sorting. In order to
1583 * support descending sorting, this function returns inverse of value when
1584 * negative coordinate is given.
1585 *
1586 * Long story short, this function uses following meaning of coordinates:
1587 * # (2 * N - 1) -- lower bound of Nth dimension,
1588 * # (2 * N) -- upper bound of Nth dimension,
1589 * # - (2 * N - 1) -- negative of lower bound of Nth dimension,
1590 * # - (2 * N) -- negative of upper bound of Nth dimension.
1591 *
1592 * When given coordinate exceeds number of cube dimensions, then 0 returned
1593 * (reproducing logic of GiST indexing of variable-length cubes).
1594 */
1595 Datum
cube_coord_llur(PG_FUNCTION_ARGS)1596 cube_coord_llur(PG_FUNCTION_ARGS)
1597 {
1598 NDBOX *cube = PG_GETARG_NDBOX_P(0);
1599 int coord = PG_GETARG_INT32(1);
1600 bool inverse = false;
1601 float8 result;
1602
1603 /* 0 is the only unsupported coordinate value */
1604 if (coord == 0)
1605 ereport(ERROR,
1606 (errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
1607 errmsg("zero cube index is not defined")));
1608
1609 /* Return inversed value for negative coordinate */
1610 if (coord < 0)
1611 {
1612 coord = -coord;
1613 inverse = true;
1614 }
1615
1616 if (coord <= 2 * DIM(cube))
1617 {
1618 /* dimension index */
1619 int index = (coord - 1) / 2;
1620
1621 /* whether this is upper bound (lower bound otherwise) */
1622 bool upper = ((coord - 1) % 2 == 1);
1623
1624 if (IS_POINT(cube))
1625 {
1626 result = cube->x[index];
1627 }
1628 else
1629 {
1630 if (upper)
1631 result = Max(cube->x[index], cube->x[index + DIM(cube)]);
1632 else
1633 result = Min(cube->x[index], cube->x[index + DIM(cube)]);
1634 }
1635 }
1636 else
1637 {
1638 /*
1639 * Return zero if coordinate is out of bound. That reproduces logic
1640 * of how cubes with low dimension number are expanded during GiST
1641 * indexing.
1642 */
1643 result = 0.0;
1644 }
1645
1646 /* Inverse value if needed */
1647 if (inverse)
1648 result = -result;
1649
1650 PG_RETURN_FLOAT8(result);
1651 }
1652
1653 /* Increase or decrease box size by a radius in at least n dimensions. */
1654 Datum
cube_enlarge(PG_FUNCTION_ARGS)1655 cube_enlarge(PG_FUNCTION_ARGS)
1656 {
1657 NDBOX *a = PG_GETARG_NDBOX_P(0);
1658 double r = PG_GETARG_FLOAT8(1);
1659 int32 n = PG_GETARG_INT32(2);
1660 NDBOX *result;
1661 int dim = 0;
1662 int size;
1663 int i,
1664 j;
1665
1666 if (n > CUBE_MAX_DIM)
1667 n = CUBE_MAX_DIM;
1668 if (r > 0 && n > 0)
1669 dim = n;
1670 if (DIM(a) > dim)
1671 dim = DIM(a);
1672
1673 size = CUBE_SIZE(dim);
1674 result = (NDBOX *) palloc0(size);
1675 SET_VARSIZE(result, size);
1676 SET_DIM(result, dim);
1677
1678 for (i = 0, j = dim; i < DIM(a); i++, j++)
1679 {
1680 if (LL_COORD(a, i) >= UR_COORD(a, i))
1681 {
1682 result->x[i] = UR_COORD(a, i) - r;
1683 result->x[j] = LL_COORD(a, i) + r;
1684 }
1685 else
1686 {
1687 result->x[i] = LL_COORD(a, i) - r;
1688 result->x[j] = UR_COORD(a, i) + r;
1689 }
1690 if (result->x[i] > result->x[j])
1691 {
1692 result->x[i] = (result->x[i] + result->x[j]) / 2;
1693 result->x[j] = result->x[i];
1694 }
1695 }
1696 /* dim > a->dim only if r > 0 */
1697 for (; i < dim; i++, j++)
1698 {
1699 result->x[i] = -r;
1700 result->x[j] = r;
1701 }
1702
1703 /*
1704 * Check if the result was in fact a point, and set the flag in the datum
1705 * accordingly. (we don't bother to repalloc it smaller)
1706 */
1707 if (cube_is_point_internal(result))
1708 {
1709 size = POINT_SIZE(dim);
1710 SET_VARSIZE(result, size);
1711 SET_POINT_BIT(result);
1712 }
1713
1714 PG_FREE_IF_COPY(a, 0);
1715 PG_RETURN_NDBOX_P(result);
1716 }
1717
1718 /* Create a one dimensional box with identical upper and lower coordinates */
1719 Datum
cube_f8(PG_FUNCTION_ARGS)1720 cube_f8(PG_FUNCTION_ARGS)
1721 {
1722 double x = PG_GETARG_FLOAT8(0);
1723 NDBOX *result;
1724 int size;
1725
1726 size = POINT_SIZE(1);
1727 result = (NDBOX *) palloc0(size);
1728 SET_VARSIZE(result, size);
1729 SET_DIM(result, 1);
1730 SET_POINT_BIT(result);
1731 result->x[0] = x;
1732
1733 PG_RETURN_NDBOX_P(result);
1734 }
1735
1736 /* Create a one dimensional box */
1737 Datum
cube_f8_f8(PG_FUNCTION_ARGS)1738 cube_f8_f8(PG_FUNCTION_ARGS)
1739 {
1740 double x0 = PG_GETARG_FLOAT8(0);
1741 double x1 = PG_GETARG_FLOAT8(1);
1742 NDBOX *result;
1743 int size;
1744
1745 if (x0 == x1)
1746 {
1747 size = POINT_SIZE(1);
1748 result = (NDBOX *) palloc0(size);
1749 SET_VARSIZE(result, size);
1750 SET_DIM(result, 1);
1751 SET_POINT_BIT(result);
1752 result->x[0] = x0;
1753 }
1754 else
1755 {
1756 size = CUBE_SIZE(1);
1757 result = (NDBOX *) palloc0(size);
1758 SET_VARSIZE(result, size);
1759 SET_DIM(result, 1);
1760 result->x[0] = x0;
1761 result->x[1] = x1;
1762 }
1763
1764 PG_RETURN_NDBOX_P(result);
1765 }
1766
1767 /* Add a dimension to an existing cube with the same values for the new
1768 coordinate */
1769 Datum
cube_c_f8(PG_FUNCTION_ARGS)1770 cube_c_f8(PG_FUNCTION_ARGS)
1771 {
1772 NDBOX *cube = PG_GETARG_NDBOX_P(0);
1773 double x = PG_GETARG_FLOAT8(1);
1774 NDBOX *result;
1775 int size;
1776 int i;
1777
1778 if (DIM(cube) + 1 > CUBE_MAX_DIM)
1779 ereport(ERROR,
1780 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
1781 errmsg("can't extend cube"),
1782 errdetail("A cube cannot have more than %d dimensions.",
1783 CUBE_MAX_DIM)));
1784
1785 if (IS_POINT(cube))
1786 {
1787 size = POINT_SIZE((DIM(cube) + 1));
1788 result = (NDBOX *) palloc0(size);
1789 SET_VARSIZE(result, size);
1790 SET_DIM(result, DIM(cube) + 1);
1791 SET_POINT_BIT(result);
1792 for (i = 0; i < DIM(cube); i++)
1793 result->x[i] = cube->x[i];
1794 result->x[DIM(result) - 1] = x;
1795 }
1796 else
1797 {
1798 size = CUBE_SIZE((DIM(cube) + 1));
1799 result = (NDBOX *) palloc0(size);
1800 SET_VARSIZE(result, size);
1801 SET_DIM(result, DIM(cube) + 1);
1802 for (i = 0; i < DIM(cube); i++)
1803 {
1804 result->x[i] = cube->x[i];
1805 result->x[DIM(result) + i] = cube->x[DIM(cube) + i];
1806 }
1807 result->x[DIM(result) - 1] = x;
1808 result->x[2 * DIM(result) - 1] = x;
1809 }
1810
1811 PG_FREE_IF_COPY(cube, 0);
1812 PG_RETURN_NDBOX_P(result);
1813 }
1814
1815 /* Add a dimension to an existing cube */
1816 Datum
cube_c_f8_f8(PG_FUNCTION_ARGS)1817 cube_c_f8_f8(PG_FUNCTION_ARGS)
1818 {
1819 NDBOX *cube = PG_GETARG_NDBOX_P(0);
1820 double x1 = PG_GETARG_FLOAT8(1);
1821 double x2 = PG_GETARG_FLOAT8(2);
1822 NDBOX *result;
1823 int size;
1824 int i;
1825
1826 if (DIM(cube) + 1 > CUBE_MAX_DIM)
1827 ereport(ERROR,
1828 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
1829 errmsg("can't extend cube"),
1830 errdetail("A cube cannot have more than %d dimensions.",
1831 CUBE_MAX_DIM)));
1832
1833 if (IS_POINT(cube) && (x1 == x2))
1834 {
1835 size = POINT_SIZE((DIM(cube) + 1));
1836 result = (NDBOX *) palloc0(size);
1837 SET_VARSIZE(result, size);
1838 SET_DIM(result, DIM(cube) + 1);
1839 SET_POINT_BIT(result);
1840 for (i = 0; i < DIM(cube); i++)
1841 result->x[i] = cube->x[i];
1842 result->x[DIM(result) - 1] = x1;
1843 }
1844 else
1845 {
1846 size = CUBE_SIZE((DIM(cube) + 1));
1847 result = (NDBOX *) palloc0(size);
1848 SET_VARSIZE(result, size);
1849 SET_DIM(result, DIM(cube) + 1);
1850 for (i = 0; i < DIM(cube); i++)
1851 {
1852 result->x[i] = LL_COORD(cube, i);
1853 result->x[DIM(result) + i] = UR_COORD(cube, i);
1854 }
1855 result->x[DIM(result) - 1] = x1;
1856 result->x[2 * DIM(result) - 1] = x2;
1857 }
1858
1859 PG_FREE_IF_COPY(cube, 0);
1860 PG_RETURN_NDBOX_P(result);
1861 }
1862