diff --git a/contrib/cube/cube--1.0.sql b/contrib/cube/cube--1.0.sql
new file mode 100644
index 0307811..053c3dc
*** a/contrib/cube/cube--1.0.sql
--- b/contrib/cube/cube--1.0.sql
*************** CREATE FUNCTION cube(cube, float8, float
*** 173,178 ****
--- 173,182 ----
  AS 'MODULE_PATHNAME', 'cube_c_f8_f8'
  LANGUAGE C IMMUTABLE STRICT;
  
+ CREATE FUNCTION cube_sort_by(cube, int4) RETURNS float8
+ AS 'MODULE_PATHNAME'
+ LANGUAGE C IMMUTABLE STRICT;
+ 
  -- Test if cube is also a point
  
  CREATE FUNCTION cube_is_point(cube)
*************** CREATE OPERATOR <@ (
*** 246,251 ****
--- 250,260 ----
  	RESTRICT = contsel, JOIN = contjoinsel
  );
  
+ CREATE OPERATOR <*> (
+ 	LEFTARG = cube, RIGHTARG = int4, PROCEDURE = cube_sort_by,
+ 	COMMUTATOR = '<*>'
+ );
+ 
  -- these are obsolete/deprecated:
  CREATE OPERATOR @ (
  	LEFTARG = cube, RIGHTARG = cube, PROCEDURE = cube_contains,
*************** RETURNS internal
*** 296,301 ****
--- 305,315 ----
  AS 'MODULE_PATHNAME'
  LANGUAGE C IMMUTABLE STRICT;
  
+ CREATE FUNCTION g_cube_distance(internal, cube, smallint, oid)
+ RETURNS float8
+ AS 'MODULE_PATHNAME'
+ LANGUAGE C STRICT;
+ 
  
  -- Create the operator classes for indexing
  
*************** CREATE OPERATOR CLASS gist_cube_ops
*** 316,325 ****
  	OPERATOR	8	<@ ,
  	OPERATOR	13	@ ,
  	OPERATOR	14	~ ,
  	FUNCTION	1	g_cube_consistent (internal, cube, int, oid, internal),
  	FUNCTION	2	g_cube_union (internal, internal),
  	FUNCTION	3	g_cube_compress (internal),
  	FUNCTION	4	g_cube_decompress (internal),
  	FUNCTION	5	g_cube_penalty (internal, internal, internal),
  	FUNCTION	6	g_cube_picksplit (internal, internal),
! 	FUNCTION	7	g_cube_same (cube, cube, internal);
--- 330,341 ----
  	OPERATOR	8	<@ ,
  	OPERATOR	13	@ ,
  	OPERATOR	14	~ ,
+ 	OPERATOR	15	<*> (cube, int) FOR ORDER BY float_ops,
  	FUNCTION	1	g_cube_consistent (internal, cube, int, oid, internal),
  	FUNCTION	2	g_cube_union (internal, internal),
  	FUNCTION	3	g_cube_compress (internal),
  	FUNCTION	4	g_cube_decompress (internal),
  	FUNCTION	5	g_cube_penalty (internal, internal, internal),
  	FUNCTION	6	g_cube_picksplit (internal, internal),
! 	FUNCTION	7	g_cube_same (cube, cube, internal),
! 	FUNCTION	8	g_cube_distance (internal, cube, smallint, oid);
diff --git a/contrib/cube/cube--unpackaged--1.0.sql b/contrib/cube/cube--unpackaged--1.0.sql
new file mode 100644
index 6859682..e14986b
*** a/contrib/cube/cube--unpackaged--1.0.sql
--- b/contrib/cube/cube--unpackaged--1.0.sql
*************** ALTER EXTENSION cube ADD operator <@(cub
*** 43,48 ****
--- 43,49 ----
  ALTER EXTENSION cube ADD operator @>(cube,cube);
  ALTER EXTENSION cube ADD operator ~(cube,cube);
  ALTER EXTENSION cube ADD operator @(cube,cube);
+ ALTER EXTENSION cube ADD operator <*>(cube,integer);
  ALTER EXTENSION cube ADD function g_cube_consistent(internal,cube,integer,oid,internal);
  ALTER EXTENSION cube ADD function g_cube_compress(internal);
  ALTER EXTENSION cube ADD function g_cube_decompress(internal);
*************** ALTER EXTENSION cube ADD function g_cube
*** 50,55 ****
--- 51,57 ----
  ALTER EXTENSION cube ADD function g_cube_picksplit(internal,internal);
  ALTER EXTENSION cube ADD function g_cube_union(internal,internal);
  ALTER EXTENSION cube ADD function g_cube_same(cube,cube,internal);
+ ALTER EXTENSION cube ADD function g_cube_distance(internal, cube, smallint, oid);
  ALTER EXTENSION cube ADD operator family cube_ops using btree;
  ALTER EXTENSION cube ADD operator class cube_ops using btree;
  ALTER EXTENSION cube ADD operator family gist_cube_ops using gist;
diff --git a/contrib/cube/cube.c b/contrib/cube/cube.c
new file mode 100644
index b591133..a5a1213
*** a/contrib/cube/cube.c
--- b/contrib/cube/cube.c
*************** PG_FUNCTION_INFO_V1(cube_dim);
*** 47,52 ****
--- 47,53 ----
  PG_FUNCTION_INFO_V1(cube_ll_coord);
  PG_FUNCTION_INFO_V1(cube_ur_coord);
  PG_FUNCTION_INFO_V1(cube_subset);
+ PG_FUNCTION_INFO_V1(cube_sort_by);
  
  Datum		cube_in(PG_FUNCTION_ARGS);
  Datum		cube(PG_FUNCTION_ARGS);
*************** Datum		cube_dim(PG_FUNCTION_ARGS);
*** 61,71 ****
  Datum		cube_ll_coord(PG_FUNCTION_ARGS);
  Datum		cube_ur_coord(PG_FUNCTION_ARGS);
  Datum		cube_subset(PG_FUNCTION_ARGS);
  
  /*
  ** GiST support methods
  */
- 
  PG_FUNCTION_INFO_V1(g_cube_consistent);
  PG_FUNCTION_INFO_V1(g_cube_compress);
  PG_FUNCTION_INFO_V1(g_cube_decompress);
--- 62,72 ----
  Datum		cube_ll_coord(PG_FUNCTION_ARGS);
  Datum		cube_ur_coord(PG_FUNCTION_ARGS);
  Datum		cube_subset(PG_FUNCTION_ARGS);
+ Datum		cube_sort_by(PG_FUNCTION_ARGS);
  
  /*
  ** GiST support methods
  */
  PG_FUNCTION_INFO_V1(g_cube_consistent);
  PG_FUNCTION_INFO_V1(g_cube_compress);
  PG_FUNCTION_INFO_V1(g_cube_decompress);
*************** PG_FUNCTION_INFO_V1(g_cube_penalty);
*** 73,78 ****
--- 74,80 ----
  PG_FUNCTION_INFO_V1(g_cube_picksplit);
  PG_FUNCTION_INFO_V1(g_cube_union);
  PG_FUNCTION_INFO_V1(g_cube_same);
+ PG_FUNCTION_INFO_V1(g_cube_distance);
  
  Datum		g_cube_consistent(PG_FUNCTION_ARGS);
  Datum		g_cube_compress(PG_FUNCTION_ARGS);
*************** Datum		g_cube_penalty(PG_FUNCTION_ARGS);
*** 81,86 ****
--- 83,89 ----
  Datum		g_cube_picksplit(PG_FUNCTION_ARGS);
  Datum		g_cube_union(PG_FUNCTION_ARGS);
  Datum		g_cube_same(PG_FUNCTION_ARGS);
+ Datum		g_cube_distance(PG_FUNCTION_ARGS);
  
  /*
  ** B-tree support functions
*************** Datum		cube_cmp(PG_FUNCTION_ARGS);
*** 104,110 ****
  /*
  ** R-tree support functions
  */
- 
  PG_FUNCTION_INFO_V1(cube_contains);
  PG_FUNCTION_INFO_V1(cube_contained);
  PG_FUNCTION_INFO_V1(cube_overlap);
--- 107,112 ----
*************** Datum		cube_enlarge(PG_FUNCTION_ARGS);
*** 136,146 ****
  int32		cube_cmp_v0(NDBOX *a, NDBOX *b);
  bool		cube_contains_v0(NDBOX *a, NDBOX *b);
  bool		cube_overlap_v0(NDBOX *a, NDBOX *b);
! NDBOX	   *cube_union_v0(NDBOX *a, NDBOX *b);
  void		rt_cube_size(NDBOX *a, double *sz);
! NDBOX	   *g_cube_binary_union(NDBOX *r1, NDBOX *r2, int *sizep);
  bool		g_cube_leaf_consistent(NDBOX *key, NDBOX *query, StrategyNumber strategy);
  bool		g_cube_internal_consistent(NDBOX *key, NDBOX *query, StrategyNumber strategy);
  
  /*
  ** Auxiliary funxtions
--- 138,149 ----
  int32		cube_cmp_v0(NDBOX *a, NDBOX *b);
  bool		cube_contains_v0(NDBOX *a, NDBOX *b);
  bool		cube_overlap_v0(NDBOX *a, NDBOX *b);
! NDBOX*		cube_union_v0(NDBOX *a, NDBOX *b);
  void		rt_cube_size(NDBOX *a, double *sz);
! NDBOX*		g_cube_binary_union(NDBOX *r1, NDBOX *r2, int *sizep);
  bool		g_cube_leaf_consistent(NDBOX *key, NDBOX *query, StrategyNumber strategy);
  bool		g_cube_internal_consistent(NDBOX *key, NDBOX *query, StrategyNumber strategy);
+ float8		cube_sort_by_v0(NDBOX *cube, int d);
  
  /*
  ** Auxiliary funxtions
*************** cube_a_f8_f8(PG_FUNCTION_ARGS)
*** 183,188 ****
--- 186,192 ----
  	int			i;
  	int			dim;
  	int			size;
+ 	bool		point = true;
  	double	   *dur,
  			   *dll;
  
*************** cube_a_f8_f8(PG_FUNCTION_ARGS)
*** 200,219 ****
  	dur = ARRPTR(ur);
  	dll = ARRPTR(ll);
  
! 	size = offsetof(NDBOX, x[0]) +sizeof(double) * 2 * dim;
  	result = (NDBOX *) palloc0(size);
  	SET_VARSIZE(result, size);
! 	result->dim = dim;
  
  	for (i = 0; i < dim; i++)
  	{
  		result->x[i] = dur[i];
  		result->x[i + dim] = dll[i];
  	}
  
  	PG_RETURN_NDBOX(result);
  }
  
  /*
  ** Allows the construction of a zero-volume cube from a float[]
  */
--- 204,234 ----
  	dur = ARRPTR(ur);
  	dll = ARRPTR(ll);
  
! 	size = CUBE_SIZE(dim);
  	result = (NDBOX *) palloc0(size);
  	SET_VARSIZE(result, size);
! 	SET_DIM(result, dim);
  
  	for (i = 0; i < dim; i++)
  	{
  		result->x[i] = dur[i];
  		result->x[i + dim] = dll[i];
+ 		if (dur[i] != dll[i])
+ 			point = false;
+ 	}
+ 
+ 	if (point)
+ 	{
+ 		size = POINT_SIZE(dim);
+ 		result = repalloc(result, size);
+ 		SET_VARSIZE(result, size);
+ 		SET_POINT_BIT(result);
  	}
  
  	PG_RETURN_NDBOX(result);
  }
  
+ 
  /*
  ** Allows the construction of a zero-volume cube from a float[]
  */
*************** cube_a_f8(PG_FUNCTION_ARGS)
*** 233,255 ****
  				 errmsg("cannot work with arrays containing NULLs")));
  
  	dim = ARRNELEMS(ur);
- 
  	dur = ARRPTR(ur);
  
! 	size = offsetof(NDBOX, x[0]) +sizeof(double) * 2 * dim;
  	result = (NDBOX *) palloc0(size);
  	SET_VARSIZE(result, size);
! 	result->dim = dim;
  
  	for (i = 0; i < dim; i++)
- 	{
  		result->x[i] = dur[i];
- 		result->x[i + dim] = dur[i];
- 	}
  
  	PG_RETURN_NDBOX(result);
  }
  
  Datum
  cube_subset(PG_FUNCTION_ARGS)
  {
--- 248,268 ----
  				 errmsg("cannot work with arrays containing NULLs")));
  
  	dim = ARRNELEMS(ur);
  	dur = ARRPTR(ur);
  
! 	size = POINT_SIZE(dim);
  	result = (NDBOX *) palloc0(size);
  	SET_VARSIZE(result, size);
! 	SET_DIM(result, dim);
! 	SET_POINT_BIT(result);
  
  	for (i = 0; i < dim; i++)
  		result->x[i] = dur[i];
  
  	PG_RETURN_NDBOX(result);
  }
  
+ 
  Datum
  cube_subset(PG_FUNCTION_ARGS)
  {
*************** cube_subset(PG_FUNCTION_ARGS)
*** 269,282 ****
  	dx = (int32 *) ARR_DATA_PTR(idx);
  
  	dim = ARRNELEMS(idx);
! 	size = offsetof(NDBOX, x[0]) +sizeof(double) * 2 * dim;
  	result = (NDBOX *) palloc0(size);
  	SET_VARSIZE(result, size);
! 	result->dim = dim;
  
  	for (i = 0; i < dim; i++)
  	{
! 		if ((dx[i] <= 0) || (dx[i] > c->dim))
  		{
  			pfree(result);
  			ereport(ERROR,
--- 282,298 ----
  	dx = (int32 *) ARR_DATA_PTR(idx);
  
  	dim = ARRNELEMS(idx);
! 	size = IS_POINT(c) ? POINT_SIZE(dim) : CUBE_SIZE(dim);
  	result = (NDBOX *) palloc0(size);
  	SET_VARSIZE(result, size);
! 	SET_DIM(result, dim);
! 
! 	if (IS_POINT(c))
! 		SET_POINT_BIT(result);
  
  	for (i = 0; i < dim; i++)
  	{
! 		if ((dx[i] <= 0) || (dx[i] > DIM(c)))
  		{
  			pfree(result);
  			ereport(ERROR,
*************** cube_subset(PG_FUNCTION_ARGS)
*** 284,303 ****
  					 errmsg("Index out of bounds")));
  		}
  		result->x[i] = c->x[dx[i] - 1];
! 		result->x[i + dim] = c->x[dx[i] + c->dim - 1];
  	}
  
  	PG_FREE_IF_COPY(c, 0);
  	PG_RETURN_NDBOX(result);
  }
  
  Datum
  cube_out(PG_FUNCTION_ARGS)
  {
  	NDBOX	   *cube = PG_GETARG_NDBOX(0);
  	StringInfoData buf;
! 	int			dim = cube->dim;
! 	bool		equal = true;
  	int			i;
  	int			ndig;
  
--- 300,329 ----
  					 errmsg("Index out of bounds")));
  		}
  		result->x[i] = c->x[dx[i] - 1];
! 		if (!IS_POINT(c))
! 			result->x[i + dim] = c->x[dx[i] + DIM(c) - 1];
  	}
  
  	PG_FREE_IF_COPY(c, 0);
  	PG_RETURN_NDBOX(result);
  }
  
+ 
+ Datum
+ cube_sort_by(PG_FUNCTION_ARGS)
+ {
+ 	NDBOX *key = PG_GETARG_NDBOX(0);
+ 	int d = PG_GETARG_INT32(1);
+ 	PG_RETURN_FLOAT8(cube_sort_by_v0(key, d));
+ }
+ 
+ 
  Datum
  cube_out(PG_FUNCTION_ARGS)
  {
  	NDBOX	   *cube = PG_GETARG_NDBOX(0);
  	StringInfoData buf;
! 	int			dim = DIM(cube);
  	int			i;
  	int			ndig;
  
*************** cube_out(PG_FUNCTION_ARGS)
*** 319,338 ****
  	{
  		if (i > 0)
  			appendStringInfo(&buf, ", ");
! 		appendStringInfo(&buf, "%.*g", ndig, cube->x[i]);
! 		if (cube->x[i] != cube->x[i + dim])
! 			equal = false;
  	}
  	appendStringInfoChar(&buf, ')');
  
! 	if (!equal)
  	{
  		appendStringInfo(&buf, ",(");
  		for (i = 0; i < dim; i++)
  		{
  			if (i > 0)
  				appendStringInfo(&buf, ", ");
! 			appendStringInfo(&buf, "%.*g", ndig, cube->x[i + dim]);
  		}
  		appendStringInfoChar(&buf, ')');
  	}
--- 345,362 ----
  	{
  		if (i > 0)
  			appendStringInfo(&buf, ", ");
! 		appendStringInfo(&buf, "%.*g", ndig, LL_COORD(cube,i));
  	}
  	appendStringInfoChar(&buf, ')');
  
! 	if (!IS_POINT(cube))
  	{
  		appendStringInfo(&buf, ",(");
  		for (i = 0; i < dim; i++)
  		{
  			if (i > 0)
  				appendStringInfo(&buf, ", ");
! 			appendStringInfo(&buf, "%.*g", ndig, UR_COORD(cube, i));
  		}
  		appendStringInfoChar(&buf, ')');
  	}
*************** g_cube_union(PG_FUNCTION_ARGS)
*** 395,403 ****
  	NDBOX	   *tmp;
  	int			i;
  
- 	/*
- 	 * fprintf(stderr, "union\n");
- 	 */
  	tmp = DatumGetNDBOX(entryvec->vector[0].key);
  
  	/*
--- 419,424 ----
*************** g_cube_union(PG_FUNCTION_ARGS)
*** 416,447 ****
  	PG_RETURN_POINTER(out);
  }
  
  /*
  ** GiST Compress and Decompress methods for boxes
  ** do not do anything.
  */
- 
  Datum
  g_cube_compress(PG_FUNCTION_ARGS)
  {
! 	PG_RETURN_DATUM(PG_GETARG_DATUM(0));
  }
  
  Datum
  g_cube_decompress(PG_FUNCTION_ARGS)
  {
! 	GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
! 	NDBOX	   *key = DatumGetNDBOX(PG_DETOAST_DATUM(entry->key));
! 
! 	if (key != DatumGetNDBOX(entry->key))
! 	{
! 		GISTENTRY  *retval = (GISTENTRY *) palloc(sizeof(GISTENTRY));
! 
! 		gistentryinit(*retval, PointerGetDatum(key),
! 					  entry->rel, entry->page,
! 					  entry->offset, FALSE);
! 		PG_RETURN_POINTER(retval);
! 	}
  	PG_RETURN_POINTER(entry);
  }
  
--- 437,458 ----
  	PG_RETURN_POINTER(out);
  }
  
+ 
  /*
  ** GiST Compress and Decompress methods for boxes
  ** do not do anything.
  */
  Datum
  g_cube_compress(PG_FUNCTION_ARGS)
  {
! 	GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
! 	PG_RETURN_POINTER(entry);
  }
  
  Datum
  g_cube_decompress(PG_FUNCTION_ARGS)
  {
! 	GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
  	PG_RETURN_POINTER(entry);
  }
  
*************** g_cube_penalty(PG_FUNCTION_ARGS)
*** 466,479 ****
  	rt_cube_size(DatumGetNDBOX(origentry->key), &tmp2);
  	*result = (float) (tmp1 - tmp2);
  
- 	/*
- 	 * fprintf(stderr, "penalty\n"); fprintf(stderr, "\t%g\n", *result);
- 	 */
  	PG_RETURN_FLOAT8(*result);
  }
  
  
- 
  /*
  ** The GiST PickSplit method for boxes
  ** We use Guttman's poly time split algorithm
--- 477,486 ----
*************** g_cube_picksplit(PG_FUNCTION_ARGS)
*** 508,517 ****
  	OffsetNumber *left,
  			   *right;
  	OffsetNumber maxoff;
! 
! 	/*
! 	 * fprintf(stderr, "picksplit\n");
! 	 */
  	maxoff = entryvec->n - 2;
  	nbytes = (maxoff + 2) * sizeof(OffsetNumber);
  	v->spl_left = (OffsetNumber *) palloc(nbytes);
--- 515,521 ----
  	OffsetNumber *left,
  			   *right;
  	OffsetNumber maxoff;
! 	 
  	maxoff = entryvec->n - 2;
  	nbytes = (maxoff + 2) * sizeof(OffsetNumber);
  	v->spl_left = (OffsetNumber *) palloc(nbytes);
*************** g_cube_picksplit(PG_FUNCTION_ARGS)
*** 627,632 ****
--- 631,637 ----
  	PG_RETURN_POINTER(v);
  }
  
+ 
  /*
  ** Equality method
  */
*************** g_cube_same(PG_FUNCTION_ARGS)
*** 642,653 ****
  	else
  		*result = FALSE;
  
- 	/*
- 	 * fprintf(stderr, "same: %s\n", (*result ? "TRUE" : "FALSE" ));
- 	 */
  	PG_RETURN_NDBOX(result);
  }
  
  /*
  ** SUPPORT ROUTINES
  */
--- 647,670 ----
  	else
  		*result = FALSE;
  
  	PG_RETURN_NDBOX(result);
  }
  
+ 
+ /*
+ ** kNN sorting support
+ */
+ Datum
+ g_cube_distance(PG_FUNCTION_ARGS)
+ {
+ 	GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
+ 	int			d = PG_GETARG_INT32(1);
+ 	NDBOX		*key = DatumGetNDBOX(entry->key);
+ 
+ 	PG_RETURN_FLOAT8(cube_sort_by_v0(key, d));
+ }
+ 
+ 
  /*
  ** SUPPORT ROUTINES
  */
*************** g_cube_leaf_consistent(NDBOX *key,
*** 658,666 ****
  {
  	bool		retval;
  
- 	/*
- 	 * fprintf(stderr, "leaf_consistent, %d\n", strategy);
- 	 */
  	switch (strategy)
  	{
  		case RTOverlapStrategyNumber:
--- 675,680 ----
*************** g_cube_leaf_consistent(NDBOX *key,
*** 683,688 ****
--- 697,703 ----
  	return (retval);
  }
  
+ 
  bool
  g_cube_internal_consistent(NDBOX *key,
  						   NDBOX *query,
*************** g_cube_internal_consistent(NDBOX *key,
*** 690,698 ****
  {
  	bool		retval;
  
- 	/*
- 	 * fprintf(stderr, "internal_consistent, %d\n", strategy);
- 	 */
  	switch (strategy)
  	{
  		case RTOverlapStrategyNumber:
--- 705,710 ----
*************** g_cube_internal_consistent(NDBOX *key,
*** 713,718 ****
--- 725,747 ----
  	return (retval);
  }
  
+ 
+ float8
+ cube_sort_by_v0(NDBOX *cube, int d)
+ {
+ 	bool is_min;
+ 
+ 	is_min = ((d % 2) == 0);
+ 	d = d/2;
+ 	if (DIM(cube) <= d)
+ 		return 0.0;
+ 	if (is_min)
+ 		return LL_COORD(cube,d);
+ 	else
+ 		return 1.0 / UR_COORD(cube,d);
+ }
+ 
+ 
  NDBOX *
  g_cube_binary_union(NDBOX *r1, NDBOX *r2, int *sizep)
  {
*************** NDBOX *
*** 730,786 ****
  cube_union_v0(NDBOX *a, NDBOX *b)
  {
  	int			i;
  	NDBOX	   *result;
  
! 	if (a->dim >= b->dim)
! 	{
! 		result = palloc0(VARSIZE(a));
! 		SET_VARSIZE(result, VARSIZE(a));
! 		result->dim = a->dim;
! 	}
! 	else
! 	{
! 		result = palloc0(VARSIZE(b));
! 		SET_VARSIZE(result, VARSIZE(b));
! 		result->dim = b->dim;
! 	}
  
  	/* swap the box pointers if needed */
! 	if (a->dim < b->dim)
  	{
  		NDBOX	   *tmp = b;
- 
  		b = a;
  		a = tmp;
  	}
  
! 	/*
! 	 * use the potentially smaller of the two boxes (b) to fill in the result,
! 	 * padding absent dimensions with zeroes
! 	 */
! 	for (i = 0; i < b->dim; i++)
  	{
! 		result->x[i] = Min(b->x[i], b->x[i + b->dim]);
! 		result->x[i + a->dim] = Max(b->x[i], b->x[i + b->dim]);
  	}
! 	for (i = b->dim; i < a->dim; i++)
  	{
! 		result->x[i] = 0;
! 		result->x[i + a->dim] = 0;
  	}
  
! 	/* compute the union */
! 	for (i = 0; i < a->dim; i++)
  	{
! 		result->x[i] =
! 			Min(Min(a->x[i], a->x[i + a->dim]), result->x[i]);
! 		result->x[i + a->dim] = Max(Max(a->x[i],
! 								   a->x[i + a->dim]), result->x[i + a->dim]);
  	}
  
  	return (result);
  }
  
  Datum
  cube_union(PG_FUNCTION_ARGS)
  {
--- 759,820 ----
  cube_union_v0(NDBOX *a, NDBOX *b)
  {
  	int			i;
+ 	bool		point_result = true;
  	NDBOX	   *result;
  
! 	/* let's try to guess result for same pointers */
! 	if (a == b)
! 		return a;
  
  	/* swap the box pointers if needed */
! 	if (DIM(a) < DIM(b))
  	{
  		NDBOX	   *tmp = b;
  		b = a;
  		a = tmp;
  	}
  
! 	result = palloc0(CUBE_SIZE(DIM(a)));
! 	SET_VARSIZE(result, CUBE_SIZE(DIM(a)));
! 	SET_DIM(result, DIM(a));
! 
! 	 /* compute the union */
! 	for (i = 0; i < DIM(b); i++)
  	{
! 		result->x[i] = Min(
! 			Min(LL_COORD(a,i), UR_COORD(a,i)),
! 			Min(LL_COORD(b,i), UR_COORD(b,i))
! 		);
! 		result->x[i + DIM(a)] = Max(
! 			Max(LL_COORD(a,i), UR_COORD(a,i)), 
! 			Max(LL_COORD(b,i), UR_COORD(b,i))
! 		);
! 		if (result->x[i] != result->x[i + DIM(a)])
! 			point_result = false;
  	}
! 	for (i = DIM(b); i < DIM(a); i++)
  	{
! 		result->x[i] = Min(0,
! 			Min(LL_COORD(a,i), UR_COORD(a,i))
! 		);
! 		result->x[i + DIM(a)] = Max(0,
! 			Max(LL_COORD(a,i), UR_COORD(a,i))
! 		);
! 		if (result->x[i] != result->x[i + DIM(a)])
! 			point_result = false;
  	}
  
! 	if (point_result)
  	{
! 		result = repalloc(result, POINT_SIZE(DIM(a)));
! 		SET_VARSIZE(result, POINT_SIZE(DIM(a)));
! 		SET_POINT_BIT(result);
  	}
  
  	return (result);
  }
  
+ 
  Datum
  cube_union(PG_FUNCTION_ARGS)
  {
*************** cube_union(PG_FUNCTION_ARGS)
*** 795,800 ****
--- 829,835 ----
  	PG_RETURN_NDBOX(res);
  }
  
+ 
  /* cube_inter */
  Datum
  cube_inter(PG_FUNCTION_ARGS)
*************** cube_inter(PG_FUNCTION_ARGS)
*** 802,825 ****
  	NDBOX	   *a = PG_GETARG_NDBOX(0);
  	NDBOX	   *b = PG_GETARG_NDBOX(1);
  	NDBOX	   *result;
! 	bool		swapped = false;
  	int			i;
  
- 	if (a->dim >= b->dim)
- 	{
- 		result = palloc0(VARSIZE(a));
- 		SET_VARSIZE(result, VARSIZE(a));
- 		result->dim = a->dim;
- 	}
- 	else
- 	{
- 		result = palloc0(VARSIZE(b));
- 		SET_VARSIZE(result, VARSIZE(b));
- 		result->dim = b->dim;
- 	}
- 
  	/* swap the box pointers if needed */
! 	if (a->dim < b->dim)
  	{
  		NDBOX	   *tmp = b;
  
--- 837,848 ----
  	NDBOX	   *a = PG_GETARG_NDBOX(0);
  	NDBOX	   *b = PG_GETARG_NDBOX(1);
  	NDBOX	   *result;
! 	bool		swapped = false,
! 				point_result = true;
  	int			i;
  
  	/* swap the box pointers if needed */
! 	if (DIM(a) < DIM(b))
  	{
  		NDBOX	   *tmp = b;
  
*************** cube_inter(PG_FUNCTION_ARGS)
*** 828,855 ****
  		swapped = true;
  	}
  
! 	/*
! 	 * use the potentially	smaller of the two boxes (b) to fill in the
! 	 * result, padding absent dimensions with zeroes
! 	 */
! 	for (i = 0; i < b->dim; i++)
  	{
! 		result->x[i] = Min(b->x[i], b->x[i + b->dim]);
! 		result->x[i + a->dim] = Max(b->x[i], b->x[i + b->dim]);
  	}
! 	for (i = b->dim; i < a->dim; i++)
  	{
! 		result->x[i] = 0;
! 		result->x[i + a->dim] = 0;
  	}
  
! 	/* compute the intersection */
! 	for (i = 0; i < a->dim; i++)
  	{
! 		result->x[i] =
! 			Max(Min(a->x[i], a->x[i + a->dim]), result->x[i]);
! 		result->x[i + a->dim] = Min(Max(a->x[i],
! 								   a->x[i + a->dim]), result->x[i + a->dim]);
  	}
  
  	if (swapped)
--- 851,890 ----
  		swapped = true;
  	}
  
! 	result = (NDBOX *) palloc0(CUBE_SIZE(DIM(a)));
! 	SET_VARSIZE(result, CUBE_SIZE(DIM(a)));
! 	SET_DIM(result, DIM(a));
! 
! 	for (i = 0; i < DIM(b); i++)
  	{
! 		result->x[i] = Max(
! 			Min(LL_COORD(a,i), UR_COORD(a,i)),
! 			Min(LL_COORD(b,i), UR_COORD(b,i))
! 		);
! 		result->x[i + DIM(a)] = Min(
! 			Max(LL_COORD(a,i), UR_COORD(a,i)),
! 			Max(LL_COORD(b,i), UR_COORD(b,i))
! 		);
! 		if (result->x[i] != result->x[i + DIM(a)])
! 			point_result = false;
  	}
! 	for (i = DIM(b); i < DIM(a); i++)
  	{
! 		result->x[i] = Max(0,
! 			Min(LL_COORD(a,i), UR_COORD(a,i))
! 		);
! 		result->x[i + DIM(a)] = Min(0,
! 			Max(LL_COORD(a,i), UR_COORD(a,i))
! 		);
! 		if (result->x[i] != result->x[i + DIM(a)])
! 			point_result = false;
  	}
  
! 	if (point_result)
  	{
! 		result = repalloc(result, POINT_SIZE(DIM(a)));
! 		SET_VARSIZE(result, POINT_SIZE(DIM(a)));
! 		SET_POINT_BIT(result);
  	}
  
  	if (swapped)
*************** cube_inter(PG_FUNCTION_ARGS)
*** 869,908 ****
  	PG_RETURN_NDBOX(result);
  }
  
  /* cube_size */
  Datum
  cube_size(PG_FUNCTION_ARGS)
  {
  	NDBOX	   *a = PG_GETARG_NDBOX(0);
  	double		result;
! 	int			i,
! 				j;
  
  	result = 1.0;
! 	for (i = 0, j = a->dim; i < a->dim; i++, j++)
! 		result = result * Abs((a->x[j] - a->x[i]));
  
  	PG_FREE_IF_COPY(a, 0);
  	PG_RETURN_FLOAT8(result);
  }
  
  void
  rt_cube_size(NDBOX *a, double *size)
  {
! 	int			i,
! 				j;
  
  	if (a == (NDBOX *) NULL)
  		*size = 0.0;
  	else
  	{
  		*size = 1.0;
! 		for (i = 0, j = a->dim; i < a->dim; i++, j++)
! 			*size = (*size) * Abs((a->x[j] - a->x[i]));
  	}
  	return;
  }
  
  /* make up a metric in which one box will be 'lower' than the other
     -- this can be useful for sorting and to determine uniqueness */
  int32
--- 904,944 ----
  	PG_RETURN_NDBOX(result);
  }
  
+ 
  /* cube_size */
  Datum
  cube_size(PG_FUNCTION_ARGS)
  {
  	NDBOX	   *a = PG_GETARG_NDBOX(0);
  	double		result;
! 	int			i;
  
  	result = 1.0;
! 	for (i = 0; i < DIM(a); i++)
! 		result = result * Abs((LL_COORD(a,i) - UR_COORD(a,i)));
  
  	PG_FREE_IF_COPY(a, 0);
  	PG_RETURN_FLOAT8(result);
  }
  
+ 
  void
  rt_cube_size(NDBOX *a, double *size)
  {
! 	int i;
  
  	if (a == (NDBOX *) NULL)
  		*size = 0.0;
  	else
  	{
  		*size = 1.0;
! 		for (i = 0; i < DIM(a); i++)
! 			*size = (*size) * Abs(UR_COORD(a,i) - LL_COORD(a,i));
  	}
  	return;
  }
  
+ 
  /* make up a metric in which one box will be 'lower' than the other
     -- this can be useful for sorting and to determine uniqueness */
  int32
*************** cube_cmp_v0(NDBOX *a, NDBOX *b)
*** 911,953 ****
  	int			i;
  	int			dim;
  
! 	dim = Min(a->dim, b->dim);
  
  	/* compare the common dimensions */
  	for (i = 0; i < dim; i++)
  	{
! 		if (Min(a->x[i], a->x[a->dim + i]) >
! 			Min(b->x[i], b->x[b->dim + i]))
  			return 1;
! 		if (Min(a->x[i], a->x[a->dim + i]) <
! 			Min(b->x[i], b->x[b->dim + i]))
  			return -1;
  	}
  	for (i = 0; i < dim; i++)
  	{
! 		if (Max(a->x[i], a->x[a->dim + i]) >
! 			Max(b->x[i], b->x[b->dim + i]))
  			return 1;
! 		if (Max(a->x[i], a->x[a->dim + i]) <
! 			Max(b->x[i], b->x[b->dim + i]))
  			return -1;
  	}
  
  	/* compare extra dimensions to zero */
! 	if (a->dim > b->dim)
  	{
! 		for (i = dim; i < a->dim; i++)
  		{
! 			if (Min(a->x[i], a->x[a->dim + i]) > 0)
  				return 1;
! 			if (Min(a->x[i], a->x[a->dim + i]) < 0)
  				return -1;
  		}
! 		for (i = dim; i < a->dim; i++)
  		{
! 			if (Max(a->x[i], a->x[a->dim + i]) > 0)
  				return 1;
! 			if (Max(a->x[i], a->x[a->dim + i]) < 0)
  				return -1;
  		}
  
--- 947,989 ----
  	int			i;
  	int			dim;
  
! 	dim = Min(DIM(a), DIM(b));
  
  	/* compare the common dimensions */
  	for (i = 0; i < dim; i++)
  	{
! 		if (Min(LL_COORD(a, i), UR_COORD(a, i)) >
! 			Min(LL_COORD(b, i), UR_COORD(b, i)))
  			return 1;
! 		if (Min(LL_COORD(a, i), UR_COORD(a, i)) <
! 			Min(LL_COORD(b, i), UR_COORD(b, i)))
  			return -1;
  	}
  	for (i = 0; i < dim; i++)
  	{
! 		if (Max(LL_COORD(a, i), UR_COORD(a, i)) >
! 			Max(LL_COORD(b, i), UR_COORD(b, i)))
  			return 1;
! 		if (Max(LL_COORD(a, i), UR_COORD(a, i)) <
! 			Max(LL_COORD(b, i), UR_COORD(b, i)))
  			return -1;
  	}
  
  	/* compare extra dimensions to zero */
! 	if (DIM(a) > DIM(b))
  	{
! 		for (i = dim; i < DIM(a); i++)
  		{
! 			if (Min(LL_COORD(a, i), UR_COORD(a, i)) > 0)
  				return 1;
! 			if (Min(LL_COORD(a, i), UR_COORD(a, i)) < 0)
  				return -1;
  		}
! 		for (i = dim; i < DIM(a); i++)
  		{
! 			if (Max(LL_COORD(a, i), UR_COORD(a, i)) > 0)
  				return 1;
! 			if (Max(LL_COORD(a, i), UR_COORD(a, i)) < 0)
  				return -1;
  		}
  
*************** cube_cmp_v0(NDBOX *a, NDBOX *b)
*** 957,976 ****
  		 */
  		return 1;
  	}
! 	if (a->dim < b->dim)
  	{
! 		for (i = dim; i < b->dim; i++)
  		{
! 			if (Min(b->x[i], b->x[b->dim + i]) > 0)
  				return -1;
! 			if (Min(b->x[i], b->x[b->dim + i]) < 0)
  				return 1;
  		}
! 		for (i = dim; i < b->dim; i++)
  		{
! 			if (Max(b->x[i], b->x[b->dim + i]) > 0)
  				return -1;
! 			if (Max(b->x[i], b->x[b->dim + i]) < 0)
  				return 1;
  		}
  
--- 993,1012 ----
  		 */
  		return 1;
  	}
! 	if (DIM(a) < DIM(b))
  	{
! 		for (i = dim; i < DIM(b); i++)
  		{
! 			if (Min(LL_COORD(b, i), UR_COORD(b, i)) > 0)
  				return -1;
! 			if (Min(LL_COORD(b, i), UR_COORD(b, i)) < 0)
  				return 1;
  		}
! 		for (i = dim; i < DIM(b); i++)
  		{
! 			if (Max(LL_COORD(b, i), UR_COORD(b, i)) > 0)
  				return -1;
! 			if (Max(LL_COORD(b, i), UR_COORD(b, i)) < 0)
  				return 1;
  		}
  
*************** cube_cmp_v0(NDBOX *a, NDBOX *b)
*** 985,990 ****
--- 1021,1027 ----
  	return 0;
  }
  
+ 
  Datum
  cube_cmp(PG_FUNCTION_ARGS)
  {
*************** cube_contains_v0(NDBOX *a, NDBOX *b)
*** 1100,1135 ****
  	if ((a == NULL) || (b == NULL))
  		return (FALSE);
  
! 	if (a->dim < b->dim)
  	{
  		/*
  		 * the further comparisons will make sense if the excess dimensions of
  		 * (b) were zeroes Since both UL and UR coordinates must be zero, we
  		 * can check them all without worrying about which is which.
  		 */
! 		for (i = a->dim; i < b->dim; i++)
  		{
! 			if (b->x[i] != 0)
  				return (FALSE);
! 			if (b->x[i + b->dim] != 0)
  				return (FALSE);
  		}
  	}
  
  	/* Can't care less about the excess dimensions of (a), if any */
! 	for (i = 0; i < Min(a->dim, b->dim); i++)
  	{
! 		if (Min(a->x[i], a->x[a->dim + i]) >
! 			Min(b->x[i], b->x[b->dim + i]))
  			return (FALSE);
! 		if (Max(a->x[i], a->x[a->dim + i]) <
! 			Max(b->x[i], b->x[b->dim + i]))
  			return (FALSE);
  	}
  
  	return (TRUE);
  }
  
  Datum
  cube_contains(PG_FUNCTION_ARGS)
  {
--- 1137,1173 ----
  	if ((a == NULL) || (b == NULL))
  		return (FALSE);
  
! 	if (DIM(a) < DIM(b))
  	{
  		/*
  		 * the further comparisons will make sense if the excess dimensions of
  		 * (b) were zeroes Since both UL and UR coordinates must be zero, we
  		 * can check them all without worrying about which is which.
  		 */
! 		for (i = DIM(a); i < DIM(b); i++)
  		{
! 			if (LL_COORD(b,i) != 0)
  				return (FALSE);
! 			if (UR_COORD(b, i) != 0)
  				return (FALSE);
  		}
  	}
  
  	/* Can't care less about the excess dimensions of (a), if any */
! 	for (i = 0; i < Min(DIM(a), DIM(b)); i++)
  	{
! 		if (Min(LL_COORD(a,i), UR_COORD(a, i)) >
! 			Min(LL_COORD(b,i), UR_COORD(b, i)))
  			return (FALSE);
! 		if (Max(LL_COORD(a,i), UR_COORD(a, i)) <
! 			Max(LL_COORD(b,i), UR_COORD(b, i)))
  			return (FALSE);
  	}
  
  	return (TRUE);
  }
  
+ 
  Datum
  cube_contains(PG_FUNCTION_ARGS)
  {
*************** cube_contains(PG_FUNCTION_ARGS)
*** 1144,1149 ****
--- 1182,1188 ----
  	PG_RETURN_BOOL(res);
  }
  
+ 
  /* Contained */
  /* Box(A) Contained by Box(B) IFF Box(B) Contains Box(A) */
  Datum
*************** cube_contained(PG_FUNCTION_ARGS)
*** 1160,1165 ****
--- 1199,1205 ----
  	PG_RETURN_BOOL(res);
  }
  
+ 
  /* Overlap */
  /* Box(A) Overlap Box(B) IFF (pt(a)LL < pt(B)UR) && (pt(b)LL < pt(a)UR) */
  bool
*************** cube_overlap_v0(NDBOX *a, NDBOX *b)
*** 1175,1181 ****
  		return (FALSE);
  
  	/* swap the box pointers if needed */
! 	if (a->dim < b->dim)
  	{
  		NDBOX	   *tmp = b;
  
--- 1215,1221 ----
  		return (FALSE);
  
  	/* swap the box pointers if needed */
! 	if (DIM(a) < DIM(b))
  	{
  		NDBOX	   *tmp = b;
  
*************** cube_overlap_v0(NDBOX *a, NDBOX *b)
*** 1184,1205 ****
  	}
  
  	/* compare within the dimensions of (b) */
! 	for (i = 0; i < b->dim; i++)
  	{
! 		if (Min(a->x[i], a->x[a->dim + i]) >
! 			Max(b->x[i], b->x[b->dim + i]))
  			return (FALSE);
! 		if (Max(a->x[i], a->x[a->dim + i]) <
! 			Min(b->x[i], b->x[b->dim + i]))
  			return (FALSE);
  	}
  
  	/* compare to zero those dimensions in (a) absent in (b) */
! 	for (i = b->dim; i < a->dim; i++)
  	{
! 		if (Min(a->x[i], a->x[a->dim + i]) > 0)
  			return (FALSE);
! 		if (Max(a->x[i], a->x[a->dim + i]) < 0)
  			return (FALSE);
  	}
  
--- 1224,1243 ----
  	}
  
  	/* compare within the dimensions of (b) */
! 	for (i = 0; i < DIM(b); i++)
  	{
! 		if (Min(LL_COORD(a,i), UR_COORD(a,i)) > Max(LL_COORD(b,i), UR_COORD(b,i)))
  			return (FALSE);
! 		if (Max(LL_COORD(a,i), UR_COORD(a,i)) < Min(LL_COORD(b,i), UR_COORD(b,i)))
  			return (FALSE);
  	}
  
  	/* compare to zero those dimensions in (a) absent in (b) */
! 	for (i = DIM(b); i < DIM(a); i++)
  	{
! 		if (Min(LL_COORD(a,i), UR_COORD(a,i)) > 0)
  			return (FALSE);
! 		if (Max(LL_COORD(a,i), UR_COORD(a,i)) < 0)
  			return (FALSE);
  	}
  
*************** cube_distance(PG_FUNCTION_ARGS)
*** 1238,1244 ****
  	int			i;
  
  	/* swap the box pointers if needed */
! 	if (a->dim < b->dim)
  	{
  		NDBOX	   *tmp = b;
  
--- 1276,1282 ----
  	int			i;
  
  	/* swap the box pointers if needed */
! 	if (DIM(a) < DIM(b))
  	{
  		NDBOX	   *tmp = b;
  
*************** cube_distance(PG_FUNCTION_ARGS)
*** 1249,1264 ****
  
  	distance = 0.0;
  	/* compute within the dimensions of (b) */
! 	for (i = 0; i < b->dim; i++)
  	{
! 		d = distance_1D(a->x[i], a->x[i + a->dim], b->x[i], b->x[i + b->dim]);
  		distance += d * d;
  	}
  
  	/* compute distance to zero for those dimensions in (a) absent in (b) */
! 	for (i = b->dim; i < a->dim; i++)
  	{
! 		d = distance_1D(a->x[i], a->x[i + a->dim], 0.0, 0.0);
  		distance += d * d;
  	}
  
--- 1287,1302 ----
  
  	distance = 0.0;
  	/* compute within the dimensions of (b) */
! 	for (i = 0; i < DIM(b); i++)
  	{
! 		d = distance_1D(LL_COORD(a,i), UR_COORD(a,i), LL_COORD(b,i), UR_COORD(b,i));
  		distance += d * d;
  	}
  
  	/* compute distance to zero for those dimensions in (a) absent in (b) */
! 	for (i = DIM(b); i < DIM(a); i++)
  	{
! 		d = distance_1D(LL_COORD(a,i), UR_COORD(a,i), 0.0, 0.0);
  		distance += d * d;
  	}
  
*************** cube_distance(PG_FUNCTION_ARGS)
*** 1276,1281 ****
--- 1314,1320 ----
  	PG_RETURN_FLOAT8(sqrt(distance));
  }
  
+ 
  static double
  distance_1D(double a1, double a2, double b1, double b2)
  {
*************** distance_1D(double a1, double a2, double
*** 1291,1325 ****
  	return (0.0);
  }
  
  /* Test if a box is also a point */
  Datum
  cube_is_point(PG_FUNCTION_ARGS)
  {
! 	NDBOX	   *a = PG_GETARG_NDBOX(0);
! 	int			i,
! 				j;
! 
! 	for (i = 0, j = a->dim; i < a->dim; i++, j++)
! 	{
! 		if (a->x[i] != a->x[j])
! 			PG_RETURN_BOOL(FALSE);
! 	}
! 
! 	PG_FREE_IF_COPY(a, 0);
! 	PG_RETURN_BOOL(TRUE);
  }
  
  /* Return dimensions in use in the data structure */
  Datum
  cube_dim(PG_FUNCTION_ARGS)
  {
  	NDBOX	   *c = PG_GETARG_NDBOX(0);
! 	int			dim = c->dim;
! 
  	PG_FREE_IF_COPY(c, 0);
  	PG_RETURN_INT32(dim);
  }
  
  /* Return a specific normalized LL coordinate */
  Datum
  cube_ll_coord(PG_FUNCTION_ARGS)
--- 1330,1359 ----
  	return (0.0);
  }
  
+ 
  /* Test if a box is also a point */
  Datum
  cube_is_point(PG_FUNCTION_ARGS)
  {
! 	NDBOX		*cube = PG_GETARG_NDBOX(0);
! 	if (IS_POINT(cube))
! 		PG_RETURN_BOOL(TRUE);
! 	else
! 		PG_RETURN_BOOL(FALSE);
  }
  
+ 
  /* Return dimensions in use in the data structure */
  Datum
  cube_dim(PG_FUNCTION_ARGS)
  {
  	NDBOX	   *c = PG_GETARG_NDBOX(0);
! 	int			dim = DIM(c);
  	PG_FREE_IF_COPY(c, 0);
  	PG_RETURN_INT32(dim);
  }
  
+ 
  /* Return a specific normalized LL coordinate */
  Datum
  cube_ll_coord(PG_FUNCTION_ARGS)
*************** cube_ll_coord(PG_FUNCTION_ARGS)
*** 1328,1335 ****
  	int			n = PG_GETARG_INT16(1);
  	double		result;
  
! 	if (c->dim >= n && n > 0)
! 		result = Min(c->x[n - 1], c->x[c->dim + n - 1]);
  	else
  		result = 0;
  
--- 1362,1369 ----
  	int			n = PG_GETARG_INT16(1);
  	double		result;
  
! 	if (DIM(c) >= n && n > 0)
! 		result = Min(LL_COORD(c, n-1), UR_COORD(c, n-1));
  	else
  		result = 0;
  
*************** cube_ll_coord(PG_FUNCTION_ARGS)
*** 1337,1342 ****
--- 1371,1377 ----
  	PG_RETURN_FLOAT8(result);
  }
  
+ 
  /* Return a specific normalized UR coordinate */
  Datum
  cube_ur_coord(PG_FUNCTION_ARGS)
*************** cube_ur_coord(PG_FUNCTION_ARGS)
*** 1345,1352 ****
  	int			n = PG_GETARG_INT16(1);
  	double		result;
  
! 	if (c->dim >= n && n > 0)
! 		result = Max(c->x[n - 1], c->x[c->dim + n - 1]);
  	else
  		result = 0;
  
--- 1380,1387 ----
  	int			n = PG_GETARG_INT16(1);
  	double		result;
  
! 	if (DIM(c) >= n && n > 0)
! 		result = Max(LL_COORD(c, n-1), UR_COORD(c, n-1));
  	else
  		result = 0;
  
*************** cube_ur_coord(PG_FUNCTION_ARGS)
*** 1354,1359 ****
--- 1389,1395 ----
  	PG_RETURN_FLOAT8(result);
  }
  
+ 
  /* Increase or decrease box size by a radius in at least n dimensions. */
  Datum
  cube_enlarge(PG_FUNCTION_ARGS)
*************** cube_enlarge(PG_FUNCTION_ARGS)
*** 1366,1400 ****
  	int			size;
  	int			i,
  				j,
! 				k;
  
  	if (n > CUBE_MAX_DIM)
  		n = CUBE_MAX_DIM;
  	if (r > 0 && n > 0)
  		dim = n;
! 	if (a->dim > dim)
! 		dim = a->dim;
! 	size = offsetof(NDBOX, x[0]) +sizeof(double) * dim * 2;
  	result = (NDBOX *) palloc0(size);
  	SET_VARSIZE(result, size);
! 	result->dim = dim;
! 	for (i = 0, j = dim, k = a->dim; i < a->dim; i++, j++, k++)
  	{
! 		if (a->x[i] >= a->x[k])
  		{
! 			result->x[i] = a->x[k] - r;
! 			result->x[j] = a->x[i] + r;
  		}
  		else
  		{
! 			result->x[i] = a->x[i] - r;
! 			result->x[j] = a->x[k] + r;
  		}
  		if (result->x[i] > result->x[j])
  		{
  			result->x[i] = (result->x[i] + result->x[j]) / 2;
  			result->x[j] = result->x[i];
  		}
  	}
  	/* dim > a->dim only if r > 0 */
  	for (; i < dim; i++, j++)
--- 1402,1441 ----
  	int			size;
  	int			i,
  				j,
! 				shrunk_coordinates = 0;
  
  	if (n > CUBE_MAX_DIM)
  		n = CUBE_MAX_DIM;
  	if (r > 0 && n > 0)
  		dim = n;
! 	if (DIM(a) > dim)
! 		dim = DIM(a);
! 
! 	size = CUBE_SIZE(dim);
  	result = (NDBOX *) palloc0(size);
  	SET_VARSIZE(result, size);
! 	SET_DIM(result, dim);
! 	
! 	for (i = 0, j = dim; i < DIM(a); i++, j++)
  	{
! 		if (LL_COORD(a,i) >= UR_COORD(a,i))
  		{
! 			result->x[i] = UR_COORD(a,i) - r;
! 			result->x[j] = LL_COORD(a,i) + r;
  		}
  		else
  		{
! 			result->x[i] = LL_COORD(a,i) - r;
! 			result->x[j] = UR_COORD(a,i) + r;
  		}
  		if (result->x[i] > result->x[j])
  		{
  			result->x[i] = (result->x[i] + result->x[j]) / 2;
  			result->x[j] = result->x[i];
+ 			shrunk_coordinates++;
  		}
+ 		else if (result->x[i] == result->x[j])
+ 			shrunk_coordinates++;
  	}
  	/* dim > a->dim only if r > 0 */
  	for (; i < dim; i++, j++)
*************** cube_enlarge(PG_FUNCTION_ARGS)
*** 1403,1412 ****
--- 1444,1464 ----
  		result->x[j] = r;
  	}
  
+ 	/* Point can arise in two cases:
+ 	   1) When argument is point and r == 0
+ 	   2) When all coordinates was set to their averages */
+ 	if ( (IS_POINT(a) && r == 0) || (shrunk_coordinates == dim) ){
+ 		size = POINT_SIZE(dim);
+ 		result = (NDBOX *) repalloc(result, size);
+ 		SET_VARSIZE(result, size);
+ 		SET_POINT_BIT(result);
+ 	}
+ 
  	PG_FREE_IF_COPY(a, 0);
  	PG_RETURN_NDBOX(result);
  }
  
+ 
  /* Create a one dimensional box with identical upper and lower coordinates */
  Datum
  cube_f8(PG_FUNCTION_ARGS)
*************** cube_f8(PG_FUNCTION_ARGS)
*** 1415,1498 ****
  	NDBOX	   *result;
  	int			size;
  
! 	size = offsetof(NDBOX, x[0]) +sizeof(double) * 2;
  	result = (NDBOX *) palloc0(size);
  	SET_VARSIZE(result, size);
! 	result->dim = 1;
! 	result->x[0] = result->x[1] = x;
  
  	PG_RETURN_NDBOX(result);
  }
  
  /* Create a one dimensional box */
  Datum
  cube_f8_f8(PG_FUNCTION_ARGS)
  {
  	double		x0 = PG_GETARG_FLOAT8(0);
  	double		x1 = PG_GETARG_FLOAT8(1);
! 	NDBOX	   *result;
  	int			size;
  
! 	size = offsetof(NDBOX, x[0]) +sizeof(double) * 2;
! 	result = (NDBOX *) palloc0(size);
! 	SET_VARSIZE(result, size);
! 	result->dim = 1;
! 	result->x[0] = x0;
! 	result->x[1] = x1;
  
  	PG_RETURN_NDBOX(result);
  }
  
  /* Add a dimension to an existing cube with the same values for the new
     coordinate */
  Datum
  cube_c_f8(PG_FUNCTION_ARGS)
  {
! 	NDBOX	   *c = PG_GETARG_NDBOX(0);
  	double		x = PG_GETARG_FLOAT8(1);
! 	NDBOX	   *result;
  	int			size;
  	int			i;
  
! 	size = offsetof(NDBOX, x[0]) +sizeof(double) * (c->dim + 1) *2;
! 	result = (NDBOX *) palloc0(size);
! 	SET_VARSIZE(result, size);
! 	result->dim = c->dim + 1;
! 	for (i = 0; i < c->dim; i++)
  	{
! 		result->x[i] = c->x[i];
! 		result->x[result->dim + i] = c->x[c->dim + i];
  	}
- 	result->x[result->dim - 1] = x;
- 	result->x[2 * result->dim - 1] = x;
  
! 	PG_FREE_IF_COPY(c, 0);
  	PG_RETURN_NDBOX(result);
  }
  
  /* Add a dimension to an existing cube */
  Datum
  cube_c_f8_f8(PG_FUNCTION_ARGS)
  {
! 	NDBOX	   *c = PG_GETARG_NDBOX(0);
  	double		x1 = PG_GETARG_FLOAT8(1);
  	double		x2 = PG_GETARG_FLOAT8(2);
  	NDBOX	   *result;
  	int			size;
  	int			i;
  
! 	size = offsetof(NDBOX, x[0]) +sizeof(double) * (c->dim + 1) *2;
! 	result = (NDBOX *) palloc0(size);
! 	SET_VARSIZE(result, size);
! 	result->dim = c->dim + 1;
! 	for (i = 0; i < c->dim; i++)
  	{
! 		result->x[i] = c->x[i];
! 		result->x[result->dim + i] = c->x[c->dim + i];
  	}
- 	result->x[result->dim - 1] = x1;
- 	result->x[2 * result->dim - 1] = x2;
  
! 	PG_FREE_IF_COPY(c, 0);
  	PG_RETURN_NDBOX(result);
  }
--- 1467,1595 ----
  	NDBOX	   *result;
  	int			size;
  
! 	size = POINT_SIZE(1);
  	result = (NDBOX *) palloc0(size);
  	SET_VARSIZE(result, size);
! 	SET_DIM(result, 1);
! 	SET_POINT_BIT(result);
! 
! 	result->x[0] = x;
  
  	PG_RETURN_NDBOX(result);
  }
  
+ 
  /* Create a one dimensional box */
  Datum
  cube_f8_f8(PG_FUNCTION_ARGS)
  {
  	double		x0 = PG_GETARG_FLOAT8(0);
  	double		x1 = PG_GETARG_FLOAT8(1);
! 	NDBOX		*result;
  	int			size;
  
! 	if (x0 == x1)
! 	{
! 		size = POINT_SIZE(1);
! 		result = (NDBOX *) palloc0(size);
! 		SET_VARSIZE(result, size);
! 		SET_DIM(result, 1);
! 		SET_POINT_BIT(result);
! 		result->x[0] = x0;
! 	}
! 	else
! 	{
! 		size = CUBE_SIZE(1);
! 		result = (NDBOX *) palloc0(size);
! 		SET_VARSIZE(result, size);
! 		SET_DIM(result, 1);
! 		result->x[0] = x0;
! 		result->x[1] = x1;
! 	}
  
  	PG_RETURN_NDBOX(result);
  }
  
+ 
  /* Add a dimension to an existing cube with the same values for the new
     coordinate */
  Datum
  cube_c_f8(PG_FUNCTION_ARGS)
  {
! 	NDBOX		*cube = PG_GETARG_NDBOX(0);
  	double		x = PG_GETARG_FLOAT8(1);
! 	NDBOX		*result;
  	int			size;
  	int			i;
  
! 	if (IS_POINT(cube))
  	{
! 		size = POINT_SIZE((DIM(cube) + 1));
! 		result = (NDBOX *) palloc0(size);
! 		SET_VARSIZE(result, size);
! 		SET_DIM(result, DIM(cube) + 1);
! 		SET_POINT_BIT(result);
! 		for (i = 0; i < DIM(cube); i++)
! 			result->x[i] = cube->x[i];
! 		result->x[DIM(result) - 1] = x;
! 	}
! 	else
! 	{
! 		size = CUBE_SIZE((DIM(cube) + 1));
! 		result = (NDBOX *) palloc0(size);
! 		SET_VARSIZE(result, size);
! 		SET_DIM(result, DIM(cube) + 1);
! 		for (i = 0; i < DIM(cube); i++)
! 		{
! 			result->x[i] = cube->x[i];
! 			result->x[DIM(result) + i] = cube->x[DIM(cube) + i];
! 		}
! 		result->x[DIM(result) - 1] = x;
! 		result->x[2*DIM(result) - 1] = x;
  	}
  
! 	PG_FREE_IF_COPY(cube, 0);
  	PG_RETURN_NDBOX(result);
  }
  
+ 
  /* Add a dimension to an existing cube */
  Datum
  cube_c_f8_f8(PG_FUNCTION_ARGS)
  {
! 	NDBOX	   *cube = PG_GETARG_NDBOX(0);
  	double		x1 = PG_GETARG_FLOAT8(1);
  	double		x2 = PG_GETARG_FLOAT8(2);
  	NDBOX	   *result;
  	int			size;
  	int			i;
  
! 	if (IS_POINT(cube) && (x1 == x2)){
! 		size = POINT_SIZE((DIM(cube) + 1));
! 		result = (NDBOX *) palloc0(size);
! 		SET_VARSIZE(result, size);
! 		SET_DIM(result, DIM(cube) + 1);
! 		SET_POINT_BIT(result);
! 		for (i = 0; i < DIM(cube); i++)
! 			result->x[i] = cube->x[i];
! 		result->x[DIM(result) - 1] = x1;
! 	}
! 	else
  	{
! 		size = CUBE_SIZE((DIM(cube) + 1));
! 		result = (NDBOX *) palloc0(size);
! 		SET_VARSIZE(result, size);
! 		SET_DIM(result, DIM(cube) + 1);
! 		for (i = 0; i < DIM(cube); i++)
! 		{
! 			result->x[i] = LL_COORD(cube, i);
! 			result->x[DIM(result) + i] = UR_COORD(cube, i);
! 		}
! 		result->x[DIM(result) - 1] = x1;
! 		result->x[2 * DIM(result) - 1] = x2;
  	}
  
! 	PG_FREE_IF_COPY(cube, 0);
  	PG_RETURN_NDBOX(result);
  }
+ 
diff --git a/contrib/cube/cubedata.h b/contrib/cube/cubedata.h
new file mode 100644
index fd0c26a..9ef297d
*** a/contrib/cube/cubedata.h
--- b/contrib/cube/cubedata.h
***************
*** 4,14 ****
  
  typedef struct NDBOX
  {
! 	int32		vl_len_;		/* varlena header (do not touch directly!) */
! 	unsigned int dim;
! 	double		x[1];
  } NDBOX;
  
! #define DatumGetNDBOX(x)	((NDBOX*)DatumGetPointer(x))
! #define PG_GETARG_NDBOX(x)	DatumGetNDBOX( PG_DETOAST_DATUM(PG_GETARG_DATUM(x)) )
! #define PG_RETURN_NDBOX(x)	PG_RETURN_POINTER(x)
--- 4,39 ----
  
  typedef struct NDBOX
  {
! 	/* varlena header (do not touch directly!) */
! 	int32 vl_len_;
! 
! 	/* 
! 	 * Header contains info about NDBOX. For binary
! 	 * compatibility with old versions it is defined
! 	 * as uint32.
! 	 * 
! 	 * Following information is stored:
! 	 *
! 	 *  bits 0-7  : number of cube dimensions;
! 	 *  bits 8-30 : not used;
! 	 *  bit  31   : point flag. If set, then NDBOX stores
! 	 *             n dimensions instead of 2*n;
! 	 */
! 	unsigned int header;
! 	double x[1];
  } NDBOX;
  
! #define DatumGetNDBOX(x)   ((NDBOX*)DatumGetPointer(x))
! #define PG_GETARG_NDBOX(x) DatumGetNDBOX( PG_DETOAST_DATUM(PG_GETARG_DATUM(x)) )
! #define PG_RETURN_NDBOX(x) PG_RETURN_POINTER(x)
! 
! #define IS_POINT(cube)		( cube->header >> 31 )
! #define SET_POINT_BIT(cube)	( cube->header = cube->header + 0x80000000 )
! #define DIM(cube)			( cube->header & 0x7fffffff )
! #define SET_DIM(cube, _dim)	( cube->header = _dim )
! 
! #define LL_COORD(cube, i) ( cube->x[i] )
! #define UR_COORD(cube, i) ( IS_POINT(cube) ? cube->x[i] : cube->x[i + DIM(cube)] )
! 
! #define POINT_SIZE(_dim) (offsetof(NDBOX, x[0]) + sizeof(double)*_dim)
! #define CUBE_SIZE(_dim) (offsetof(NDBOX, x[0]) + sizeof(double)*_dim*2)
diff --git a/contrib/cube/cubeparse.y b/contrib/cube/cubeparse.y
new file mode 100644
index 21fe537..abc0790
*** a/contrib/cube/cubeparse.y
--- b/contrib/cube/cubeparse.y
*************** write_box(unsigned int dim, char *str1, 
*** 175,185 ****
  	NDBOX	   *bp;
  	char	   *s;
  	int			i;
! 	int			size = offsetof(NDBOX, x[0]) + sizeof(double) * dim * 2;
  
  	bp = palloc0(size);
  	SET_VARSIZE(bp, size);
! 	bp->dim = dim;
  
  	s = str1;
  	bp->x[i=0] = strtod(s, NULL);
--- 175,186 ----
  	NDBOX	   *bp;
  	char	   *s;
  	int			i;
! 	int			size = CUBE_SIZE(dim);
! 	bool		point = true;
  
  	bp = palloc0(size);
  	SET_VARSIZE(bp, size);
! 	SET_DIM(bp, dim);
  
  	s = str1;
  	bp->x[i=0] = strtod(s, NULL);
*************** write_box(unsigned int dim, char *str1, 
*** 195,200 ****
--- 196,213 ----
  	{
  		s++; i++;
  		bp->x[i] = strtod(s, NULL);
+ 		if (bp->x[i] != bp->x[i-dim])
+ 			point = false;
+ 	}
+ 
+ 	if (bp->x[0] != bp->x[dim])
+ 		point = false;
+ 	if (point)
+ 	{
+ 		size = POINT_SIZE(dim);
+ 		bp = repalloc(bp, size);
+ 		SET_VARSIZE(bp, size);
+ 		SET_POINT_BIT(bp);
  	}
  
  	return(bp);
*************** write_box(unsigned int dim, char *str1, 
*** 203,233 ****
  static NDBOX *
  write_point_as_box(char *str, int dim)
  {
!   NDBOX		   *bp;
!   int			i,
  				size;
!   double		x;
!   char		   *s = str;
! 
!   size = offsetof(NDBOX, x[0]) + sizeof(double) * dim * 2;
  
!   bp = palloc0(size);
!   SET_VARSIZE(bp, size);
!   bp->dim = dim;
  
!   i = 0;
!   x = strtod(s, NULL);
!   bp->x[0] = x;
!   bp->x[dim] = x;
!   while ((s = strchr(s, ',')) != NULL)
!   {
! 	  s++; i++;
! 	  x = strtod(s, NULL);
! 	  bp->x[i] = x;
! 	  bp->x[i+dim] = x;
!   }
  
!   return(bp);
  }
  
  #include "cubescan.c"
--- 216,244 ----
  static NDBOX *
  write_point_as_box(char *str, int dim)
  {
! 	NDBOX		*bp;
! 	int			i,
  				size;
! 	double		x;
! 	char		*s = str;
  
! 	size = POINT_SIZE(dim);
! 	bp = palloc0(size);
! 	SET_VARSIZE(bp, size);
! 	SET_DIM(bp, dim);
! 	SET_POINT_BIT(bp);
  
! 	i = 0;
! 	x = strtod(s, NULL);
! 	bp->x[0] = x;
! 	while ((s = strchr(s, ',')) != NULL)
! 	{
! 		s++; i++;
! 		x = strtod(s, NULL);
! 		bp->x[i] = x;
! 	}
  
! 	return(bp);
  }
  
  #include "cubescan.c"
diff --git a/contrib/cube/expected/cube_1.out b/contrib/cube/expected/cube_1.out
new file mode 100644
index fefebf5..c07d61d
*** a/contrib/cube/expected/cube_1.out
--- b/contrib/cube/expected/cube_1.out
*************** SELECT cube_subset(cube('(1,3,5),(6,7,8)
*** 473,480 ****
--- 473,557 ----
   (5, 3, 1, 1),(8, 7, 6, 6)
  (1 row)
  
+ SELECT cube_subset(cube('(1,3,5),(1,3,5)'), ARRAY[3,2,1,1]);
+  cube_subset  
+ --------------
+  (5, 3, 1, 1)
+ (1 row)
+ 
  SELECT cube_subset(cube('(1,3,5),(6,7,8)'), ARRAY[4,0]);
  ERROR:  Index out of bounds
+ SELECT cube_subset(cube('(6,7,8),(6,7,8)'), ARRAY[4,0]);
+ ERROR:  Index out of bounds
+ --
+ -- Test point processing
+ --
+ SELECT cube('(1,2),(1,2)'); -- cube_in
+   cube  
+ --------
+  (1, 2)
+ (1 row)
+ 
+ SELECT cube('{0,1,2}'::float[], '{0,1,2}'::float[]); -- cube_a_f8_f8
+    cube    
+ -----------
+  (0, 1, 2)
+ (1 row)
+ 
+ SELECT cube('{5,6,7,8}'::float[]); -- cube_a_f8
+      cube     
+ --------------
+  (5, 6, 7, 8)
+ (1 row)
+ 
+ SELECT cube(1.37); -- cube_f8
+   cube  
+ --------
+  (1.37)
+ (1 row)
+ 
+ SELECT cube(1.37, 1.37); -- cube_f8_f8
+   cube  
+ --------
+  (1.37)
+ (1 row)
+ 
+ SELECT cube(cube(1,1), 42); -- cube_c_f8
+   cube   
+ ---------
+  (1, 42)
+ (1 row)
+ 
+ SELECT cube(cube(1,2), 42); -- cube_c_f8
+       cube       
+ -----------------
+  (1, 42),(2, 42)
+ (1 row)
+ 
+ SELECT cube(cube(1,1), 42, 42); -- cube_c_f8_f8
+   cube   
+ ---------
+  (1, 42)
+ (1 row)
+ 
+ SELECT cube(cube(1,1), 42, 24); -- cube_c_f8_f8
+       cube       
+ -----------------
+  (1, 42),(1, 24)
+ (1 row)
+ 
+ SELECT cube(cube(1,2), 42, 42); -- cube_c_f8_f8
+       cube       
+ -----------------
+  (1, 42),(2, 42)
+ (1 row)
+ 
+ SELECT cube(cube(1,2), 42, 24); -- cube_c_f8_f8
+       cube       
+ -----------------
+  (1, 42),(2, 24)
+ (1 row)
+ 
  --
  -- Testing limit of CUBE_MAX_DIM dimensions check in cube_in.
  --
*************** SELECT cube_distance('(0)'::cube,'(.3,.4
*** 878,883 ****
--- 955,978 ----
             0.5
  (1 row)
  
+ SELECT cube_distance('(2,3,4)'::cube,'(2,3,4)'::cube);
+  cube_distance 
+ ---------------
+              0
+ (1 row)
+ 
+ SELECT cube_distance('(42,42,42,42)'::cube,'(137,137,137,137)'::cube);
+  cube_distance 
+ ---------------
+            190
+ (1 row)
+ 
+ SELECT cube_distance('(42,42,42)'::cube,'(137,137)'::cube);
+   cube_distance   
+ ------------------
+  140.762210837994
+ (1 row)
+ 
  -- Test of cube function (text to cube)
  --
  SELECT cube('(1,1.2)'::text);
*************** SELECT cube_dim('(0,0,0)'::cube);
*** 912,917 ****
--- 1007,1024 ----
          3
  (1 row)
  
+ SELECT cube_dim('(42,42,42),(42,42,42)'::cube);
+  cube_dim 
+ ----------
+         3
+ (1 row)
+ 
+ SELECT cube_dim('(4,8,15,16,23),(4,8,15,16,23)'::cube);
+  cube_dim 
+ ----------
+         5
+ (1 row)
+ 
  -- Test of cube_ll_coord function (retrieves LL coodinate values)
  --
  SELECT cube_ll_coord('(-1,1),(2,-2)'::cube, 1);
*************** SELECT cube_ll_coord('(-1,1),(2,-2)'::cu
*** 932,937 ****
--- 1039,1080 ----
               0
  (1 row)
  
+ SELECT cube_ll_coord('(1,2),(1,2)'::cube, 1);
+  cube_ll_coord 
+ ---------------
+              1
+ (1 row)
+ 
+ SELECT cube_ll_coord('(1,2),(1,2)'::cube, 2);
+  cube_ll_coord 
+ ---------------
+              2
+ (1 row)
+ 
+ SELECT cube_ll_coord('(1,2),(1,2)'::cube, 3);
+  cube_ll_coord 
+ ---------------
+              0
+ (1 row)
+ 
+ SELECT cube_ll_coord('(42,137)'::cube, 1);
+  cube_ll_coord 
+ ---------------
+             42
+ (1 row)
+ 
+ SELECT cube_ll_coord('(42,137)'::cube, 2);
+  cube_ll_coord 
+ ---------------
+            137
+ (1 row)
+ 
+ SELECT cube_ll_coord('(42,137)'::cube, 3);
+  cube_ll_coord 
+ ---------------
+              0
+ (1 row)
+ 
  -- Test of cube_ur_coord function (retrieves UR coodinate values)
  --
  SELECT cube_ur_coord('(-1,1),(2,-2)'::cube, 1);
*************** SELECT cube_ur_coord('(-1,1),(2,-2)'::cu
*** 952,957 ****
--- 1095,1136 ----
               0
  (1 row)
  
+ SELECT cube_ur_coord('(1,2),(1,2)'::cube, 1);
+  cube_ur_coord 
+ ---------------
+              1
+ (1 row)
+ 
+ SELECT cube_ur_coord('(1,2),(1,2)'::cube, 2);
+  cube_ur_coord 
+ ---------------
+              2
+ (1 row)
+ 
+ SELECT cube_ur_coord('(1,2),(1,2)'::cube, 3);
+  cube_ur_coord 
+ ---------------
+              0
+ (1 row)
+ 
+ SELECT cube_ur_coord('(42,137)'::cube, 1);
+  cube_ur_coord 
+ ---------------
+             42
+ (1 row)
+ 
+ SELECT cube_ur_coord('(42,137)'::cube, 2);
+  cube_ur_coord 
+ ---------------
+            137
+ (1 row)
+ 
+ SELECT cube_ur_coord('(42,137)'::cube, 3);
+  cube_ur_coord 
+ ---------------
+              0
+ (1 row)
+ 
  -- Test of cube_is_point
  --
  SELECT cube_is_point('(0)'::cube);
*************** SELECT cube_enlarge('(2,-2),(-3,7)'::cub
*** 1100,1105 ****
--- 1279,1386 ----
   (-0.5, 1),(-0.5, 4)
  (1 row)
  
+ SELECT cube_enlarge('(42,-23,-23),(42,23,23)'::cube, -23, 5);
+  cube_enlarge 
+ --------------
+  (42, 0, 0)
+ (1 row)
+ 
+ SELECT cube_enlarge('(42,-23,-23),(42,23,23)'::cube, -24, 5);
+  cube_enlarge 
+ --------------
+  (42, 0, 0)
+ (1 row)
+ 
+ -- Test of cube_union (MBR for two cubes)
+ --
+ SELECT cube_union('(1,2),(3,4)'::cube, '(5,6,7),(8,9,10)'::cube);
+       cube_union      
+ ----------------------
+  (1, 2, 0),(8, 9, 10)
+ (1 row)
+ 
+ SELECT cube_union('(1,2)'::cube, '(4,2,0,0)'::cube);
+         cube_union         
+ ---------------------------
+  (1, 2, 0, 0),(4, 2, 0, 0)
+ (1 row)
+ 
+ SELECT cube_union('(1,2),(1,2)'::cube, '(4,2),(4,2)'::cube);
+   cube_union   
+ ---------------
+  (1, 2),(4, 2)
+ (1 row)
+ 
+ SELECT cube_union('(1,2),(1,2)'::cube, '(1,2),(1,2)'::cube);
+  cube_union 
+ ------------
+  (1, 2)
+ (1 row)
+ 
+ SELECT cube_union('(1,2),(1,2)'::cube, '(1,2,0),(1,2,0)'::cube);
+  cube_union 
+ ------------
+  (1, 2, 0)
+ (1 row)
+ 
+ -- Test of cube_inter
+ --
+ SELECT cube_inter('(1,2),(10,11)'::cube, '(3,4), (16,15)'::cube); -- intersects
+    cube_inter    
+ -----------------
+  (3, 4),(10, 11)
+ (1 row)
+ 
+ SELECT cube_inter('(1,2),(10,11)'::cube, '(3,4), (6,5)'::cube); -- includes
+   cube_inter   
+ ---------------
+  (3, 4),(6, 5)
+ (1 row)
+ 
+ SELECT cube_inter('(1,2),(10,11)'::cube, '(13,14), (16,15)'::cube); -- no intersection
+     cube_inter     
+ -------------------
+  (13, 14),(10, 11)
+ (1 row)
+ 
+ SELECT cube_inter('(1,2),(10,11)'::cube, '(3,14), (16,15)'::cube); -- no intersection, but one dimension intersects
+     cube_inter    
+ ------------------
+  (3, 14),(10, 11)
+ (1 row)
+ 
+ SELECT cube_inter('(1,2),(10,11)'::cube, '(10,11), (16,15)'::cube); -- point intersection
+  cube_inter 
+ ------------
+  (10, 11)
+ (1 row)
+ 
+ SELECT cube_inter('(1,2,3)'::cube, '(1,2,3)'::cube); -- point args
+  cube_inter 
+ ------------
+  (1, 2, 3)
+ (1 row)
+ 
+ SELECT cube_inter('(1,2,3)'::cube, '(5,6,3)'::cube); -- point args
+      cube_inter      
+ ---------------------
+  (5, 6, 3),(1, 2, 3)
+ (1 row)
+ 
+ -- Test of cube_size
+ --
+ SELECT cube_size('(4,8),(15,16)'::cube);
+  cube_size 
+ -----------
+         88
+ (1 row)
+ 
+ SELECT cube_size('(42,137)'::cube);
+  cube_size 
+ -----------
+          0
+ (1 row)
+ 
  -- Load some example data and build the index
  --
  CREATE TABLE test_cube (c cube);
diff --git a/contrib/cube/sql/cube.sql b/contrib/cube/sql/cube.sql
new file mode 100644
index 02e068e..d58974c
*** a/contrib/cube/sql/cube.sql
--- b/contrib/cube/sql/cube.sql
*************** SELECT cube('{0,1,2}'::float[], '{3}'::f
*** 112,118 ****
--- 112,135 ----
  SELECT cube(NULL::float[], '{3}'::float[]);
  SELECT cube('{0,1,2}'::float[]);
  SELECT cube_subset(cube('(1,3,5),(6,7,8)'), ARRAY[3,2,1,1]);
+ SELECT cube_subset(cube('(1,3,5),(1,3,5)'), ARRAY[3,2,1,1]);
  SELECT cube_subset(cube('(1,3,5),(6,7,8)'), ARRAY[4,0]);
+ SELECT cube_subset(cube('(6,7,8),(6,7,8)'), ARRAY[4,0]);
+ 
+ --
+ -- Test point processing
+ --
+ SELECT cube('(1,2),(1,2)'); -- cube_in
+ SELECT cube('{0,1,2}'::float[], '{0,1,2}'::float[]); -- cube_a_f8_f8
+ SELECT cube('{5,6,7,8}'::float[]); -- cube_a_f8
+ SELECT cube(1.37); -- cube_f8
+ SELECT cube(1.37, 1.37); -- cube_f8_f8
+ SELECT cube(cube(1,1), 42); -- cube_c_f8
+ SELECT cube(cube(1,2), 42); -- cube_c_f8
+ SELECT cube(cube(1,1), 42, 42); -- cube_c_f8_f8
+ SELECT cube(cube(1,1), 42, 24); -- cube_c_f8_f8
+ SELECT cube(cube(1,2), 42, 42); -- cube_c_f8_f8
+ SELECT cube(cube(1,2), 42, 24); -- cube_c_f8_f8
  
  --
  -- Testing limit of CUBE_MAX_DIM dimensions check in cube_in.
*************** SELECT '(-1,-1),(1,1)'::cube            
*** 212,217 ****
--- 229,237 ----
  --
  SELECT cube_distance('(0)'::cube,'(2,2,2,2)'::cube);
  SELECT cube_distance('(0)'::cube,'(.3,.4)'::cube);
+ SELECT cube_distance('(2,3,4)'::cube,'(2,3,4)'::cube);
+ SELECT cube_distance('(42,42,42,42)'::cube,'(137,137,137,137)'::cube);
+ SELECT cube_distance('(42,42,42)'::cube,'(137,137)'::cube);
  
  -- Test of cube function (text to cube)
  --
*************** SELECT cube(NULL);
*** 223,240 ****
--- 243,274 ----
  SELECT cube_dim('(0)'::cube);
  SELECT cube_dim('(0,0)'::cube);
  SELECT cube_dim('(0,0,0)'::cube);
+ SELECT cube_dim('(42,42,42),(42,42,42)'::cube);
+ SELECT cube_dim('(4,8,15,16,23),(4,8,15,16,23)'::cube);
  
  -- Test of cube_ll_coord function (retrieves LL coodinate values)
  --
  SELECT cube_ll_coord('(-1,1),(2,-2)'::cube, 1);
  SELECT cube_ll_coord('(-1,1),(2,-2)'::cube, 2);
  SELECT cube_ll_coord('(-1,1),(2,-2)'::cube, 3);
+ SELECT cube_ll_coord('(1,2),(1,2)'::cube, 1);
+ SELECT cube_ll_coord('(1,2),(1,2)'::cube, 2);
+ SELECT cube_ll_coord('(1,2),(1,2)'::cube, 3);
+ SELECT cube_ll_coord('(42,137)'::cube, 1);
+ SELECT cube_ll_coord('(42,137)'::cube, 2);
+ SELECT cube_ll_coord('(42,137)'::cube, 3);
  
  -- Test of cube_ur_coord function (retrieves UR coodinate values)
  --
  SELECT cube_ur_coord('(-1,1),(2,-2)'::cube, 1);
  SELECT cube_ur_coord('(-1,1),(2,-2)'::cube, 2);
  SELECT cube_ur_coord('(-1,1),(2,-2)'::cube, 3);
+ SELECT cube_ur_coord('(1,2),(1,2)'::cube, 1);
+ SELECT cube_ur_coord('(1,2),(1,2)'::cube, 2);
+ SELECT cube_ur_coord('(1,2),(1,2)'::cube, 3);
+ SELECT cube_ur_coord('(42,137)'::cube, 1);
+ SELECT cube_ur_coord('(42,137)'::cube, 2);
+ SELECT cube_ur_coord('(42,137)'::cube, 3);
  
  -- Test of cube_is_point
  --
*************** SELECT cube_enlarge('(2,-2),(-3,7)'::cub
*** 265,270 ****
--- 299,329 ----
  SELECT cube_enlarge('(2,-2),(-3,7)'::cube, 3, 2);
  SELECT cube_enlarge('(2,-2),(-3,7)'::cube, -1, 2);
  SELECT cube_enlarge('(2,-2),(-3,7)'::cube, -3, 2);
+ SELECT cube_enlarge('(42,-23,-23),(42,23,23)'::cube, -23, 5);
+ SELECT cube_enlarge('(42,-23,-23),(42,23,23)'::cube, -24, 5);
+ 
+ -- Test of cube_union (MBR for two cubes)
+ --
+ SELECT cube_union('(1,2),(3,4)'::cube, '(5,6,7),(8,9,10)'::cube);
+ SELECT cube_union('(1,2)'::cube, '(4,2,0,0)'::cube);
+ SELECT cube_union('(1,2),(1,2)'::cube, '(4,2),(4,2)'::cube);
+ SELECT cube_union('(1,2),(1,2)'::cube, '(1,2),(1,2)'::cube);
+ SELECT cube_union('(1,2),(1,2)'::cube, '(1,2,0),(1,2,0)'::cube);
+ 
+ -- Test of cube_inter
+ --
+ SELECT cube_inter('(1,2),(10,11)'::cube, '(3,4), (16,15)'::cube); -- intersects
+ SELECT cube_inter('(1,2),(10,11)'::cube, '(3,4), (6,5)'::cube); -- includes
+ SELECT cube_inter('(1,2),(10,11)'::cube, '(13,14), (16,15)'::cube); -- no intersection
+ SELECT cube_inter('(1,2),(10,11)'::cube, '(3,14), (16,15)'::cube); -- no intersection, but one dimension intersects
+ SELECT cube_inter('(1,2),(10,11)'::cube, '(10,11), (16,15)'::cube); -- point intersection
+ SELECT cube_inter('(1,2,3)'::cube, '(1,2,3)'::cube); -- point args
+ SELECT cube_inter('(1,2,3)'::cube, '(5,6,3)'::cube); -- point args
+ 
+ -- Test of cube_size
+ --
+ SELECT cube_size('(4,8),(15,16)'::cube);
+ SELECT cube_size('(42,137)'::cube);
  
  -- Load some example data and build the index
  --
