diff --git a/contrib/cube/cube--1.0.sql b/contrib/cube/cube--1.0.sql
index 0307811..9c821de 100644
--- a/contrib/cube/cube--1.0.sql
+++ b/contrib/cube/cube--1.0.sql
@@ -135,7 +135,17 @@ LANGUAGE C IMMUTABLE STRICT;
 
 -- proximity routines
 
-CREATE FUNCTION cube_distance(cube, cube)
+CREATE FUNCTION distance_taxicab(cube, cube)
+RETURNS float8
+AS 'MODULE_PATHNAME'
+LANGUAGE C IMMUTABLE STRICT;
+
+CREATE FUNCTION distance_euclid(cube, cube)
+RETURNS float8
+AS 'MODULE_PATHNAME'
+LANGUAGE C IMMUTABLE STRICT;
+
+CREATE FUNCTION distance_chebyshev(cube, cube)
 RETURNS float8
 AS 'MODULE_PATHNAME'
 LANGUAGE C IMMUTABLE STRICT;
@@ -157,6 +167,11 @@ RETURNS float8
 AS 'MODULE_PATHNAME'
 LANGUAGE C IMMUTABLE STRICT;
 
+CREATE FUNCTION cube_coord(cube, int4)
+RETURNS float8
+AS 'MODULE_PATHNAME'
+LANGUAGE C IMMUTABLE STRICT;
+
 CREATE FUNCTION cube(float8) RETURNS cube
 AS 'MODULE_PATHNAME', 'cube_f8'
 LANGUAGE C IMMUTABLE STRICT;
@@ -246,6 +261,25 @@ CREATE OPERATOR <@ (
 	RESTRICT = contsel, JOIN = contjoinsel
 );
 
+CREATE OPERATOR -> (
+	LEFTARG = cube, RIGHTARG = int, PROCEDURE = cube_coord
+);
+
+CREATE OPERATOR <#> (
+	LEFTARG = cube, RIGHTARG = cube, PROCEDURE = distance_taxicab,
+	COMMUTATOR = '<#>'
+);
+
+CREATE OPERATOR <-> (
+	LEFTARG = cube, RIGHTARG = cube, PROCEDURE = distance_euclid,
+	COMMUTATOR = '<->'
+);
+
+CREATE OPERATOR <=> (
+	LEFTARG = cube, RIGHTARG = cube, PROCEDURE = distance_chebyshev,
+	COMMUTATOR = '<=>'
+);
+
 -- these are obsolete/deprecated:
 CREATE OPERATOR @ (
 	LEFTARG = cube, RIGHTARG = cube, PROCEDURE = cube_contains,
@@ -296,6 +330,10 @@ RETURNS internal
 AS 'MODULE_PATHNAME'
 LANGUAGE C IMMUTABLE STRICT;
 
+CREATE FUNCTION g_cube_distance (internal, cube, smallint, oid)
+RETURNS internal
+AS 'MODULE_PATHNAME'
+LANGUAGE C IMMUTABLE STRICT;
 
 -- Create the operator classes for indexing
 
@@ -316,10 +354,17 @@ CREATE OPERATOR CLASS gist_cube_ops
 	OPERATOR	8	<@ ,
 	OPERATOR	13	@ ,
 	OPERATOR	14	~ ,
+	OPERATOR	15	-> (cube, int) FOR ORDER BY float_ops,
+	OPERATOR	16	<#> (cube, cube) FOR ORDER BY float_ops,
+	OPERATOR	17	<-> (cube, cube) FOR ORDER BY float_ops,
+	OPERATOR	18	<=> (cube, cube) 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	7	g_cube_same (cube, cube, internal),
+	FUNCTION	8	g_cube_distance (internal, cube, smallint, oid);
+
diff --git a/contrib/cube/cube.c b/contrib/cube/cube.c
index dab0e6e..0702fac 100644
--- a/contrib/cube/cube.c
+++ b/contrib/cube/cube.c
@@ -45,6 +45,7 @@ PG_FUNCTION_INFO_V1(cube_c_f8_f8);
 PG_FUNCTION_INFO_V1(cube_dim);
 PG_FUNCTION_INFO_V1(cube_ll_coord);
 PG_FUNCTION_INFO_V1(cube_ur_coord);
+PG_FUNCTION_INFO_V1(cube_coord);
 PG_FUNCTION_INFO_V1(cube_subset);
 
 Datum		cube_in(PG_FUNCTION_ARGS);
@@ -58,6 +59,7 @@ Datum		cube_c_f8_f8(PG_FUNCTION_ARGS);
 Datum		cube_dim(PG_FUNCTION_ARGS);
 Datum		cube_ll_coord(PG_FUNCTION_ARGS);
 Datum		cube_ur_coord(PG_FUNCTION_ARGS);
+Datum		cube_coord(PG_FUNCTION_ARGS);
 Datum		cube_subset(PG_FUNCTION_ARGS);
 
 /*
@@ -71,6 +73,7 @@ PG_FUNCTION_INFO_V1(g_cube_penalty);
 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);
@@ -79,6 +82,7 @@ Datum		g_cube_penalty(PG_FUNCTION_ARGS);
 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
@@ -120,11 +124,15 @@ Datum		cube_size(PG_FUNCTION_ARGS);
 /*
 ** miscellaneous
 */
-PG_FUNCTION_INFO_V1(cube_distance);
+PG_FUNCTION_INFO_V1(distance_taxicab);
+PG_FUNCTION_INFO_V1(distance_euclid);
+PG_FUNCTION_INFO_V1(distance_chebyshev);
 PG_FUNCTION_INFO_V1(cube_is_point);
 PG_FUNCTION_INFO_V1(cube_enlarge);
 
-Datum		cube_distance(PG_FUNCTION_ARGS);
+Datum		distance_taxicab(PG_FUNCTION_ARGS);
+Datum		distance_euclid(PG_FUNCTION_ARGS);
+Datum		distance_chebyshev(PG_FUNCTION_ARGS);
 Datum		cube_is_point(PG_FUNCTION_ARGS);
 Datum		cube_enlarge(PG_FUNCTION_ARGS);
 
@@ -1219,6 +1227,61 @@ cube_overlap(PG_FUNCTION_ARGS)
 	PG_RETURN_BOOL(res);
 }
 
+static double
+distance_1D(double a1, double a2, double b1, double b2)
+{
+	/* interval (a) is entirely on the left of (b) */
+	if ((a1 <= b1) && (a2 <= b1) && (a1 <= b2) && (a2 <= b2))
+		return (Min(b1, b2) - Max(a1, a2));
+
+	/* interval (a) is entirely on the right of (b) */
+	if ((a1 > b1) && (a2 > b1) && (a1 > b2) && (a2 > b2))
+		return (Min(a1, a2) - Max(b1, b2));
+
+	/* the rest are all sorts of intersections */
+	return (0.0);
+}
+
+Datum
+distance_taxicab(PG_FUNCTION_ARGS)
+{
+	NDBOX	   *a = PG_GETARG_NDBOX(0),
+			   *b = PG_GETARG_NDBOX(1);
+	bool		swapped = false;
+	double		distance;
+	int			i;
+
+	/* swap the box pointers if needed */
+	if (a->dim < b->dim)
+	{
+		NDBOX	   *tmp = b;
+		b = a;
+		a = tmp;
+		swapped = true;
+	}
+
+	distance = 0.0;
+	/* compute within the dimensions of (b) */
+	for (i = 0; i < b->dim; i++)
+		distance += abs(distance_1D(a->x[i], a->x[i + a->dim], b->x[i], b->x[i + b->dim]));
+
+	/* compute distance to zero for those dimensions in (a) absent in (b) */
+	for (i = b->dim; i < a->dim; i++)
+		distance += abs(distance_1D(a->x[i], a->x[i + a->dim], 0.0, 0.0));
+
+	if (swapped)
+	{
+		PG_FREE_IF_COPY(b, 0);
+		PG_FREE_IF_COPY(a, 1);
+	}
+	else
+	{
+		PG_FREE_IF_COPY(a, 0);
+		PG_FREE_IF_COPY(b, 1);
+	}
+
+	PG_RETURN_FLOAT8(distance);
+}
 
 /* Distance */
 /* The distance is computed as a per axis sum of the squared distances
@@ -1226,7 +1289,7 @@ cube_overlap(PG_FUNCTION_ARGS)
    distance between overlapping projections, this metric coincides with the
    "common sense" geometric distance */
 Datum
-cube_distance(PG_FUNCTION_ARGS)
+distance_euclid(PG_FUNCTION_ARGS)
 {
 	NDBOX	   *a = PG_GETARG_NDBOX(0),
 			   *b = PG_GETARG_NDBOX(1);
@@ -1239,7 +1302,6 @@ cube_distance(PG_FUNCTION_ARGS)
 	if (a->dim < b->dim)
 	{
 		NDBOX	   *tmp = b;
-
 		b = a;
 		a = tmp;
 		swapped = true;
@@ -1250,16 +1312,18 @@ cube_distance(PG_FUNCTION_ARGS)
 	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;
+		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;
+		distance += d*d;
 	}
 
+	distance = sqrt(distance);
+
 	if (swapped)
 	{
 		PG_FREE_IF_COPY(b, 0);
@@ -1271,22 +1335,96 @@ cube_distance(PG_FUNCTION_ARGS)
 		PG_FREE_IF_COPY(b, 1);
 	}
 
-	PG_RETURN_FLOAT8(sqrt(distance));
+	PG_RETURN_FLOAT8(distance);
 }
 
-static double
-distance_1D(double a1, double a2, double b1, double b2)
+/* Distance */
+/* The distance is computed as a per axis sum of the squared distances
+   between 1D projections of the boxes onto Cartesian axes. Assuming zero
+   distance between overlapping projections, this metric coincides with the
+   "common sense" geometric distance */
+Datum
+distance_chebyshev(PG_FUNCTION_ARGS)
 {
-	/* interval (a) is entirely on the left of (b) */
-	if ((a1 <= b1) && (a2 <= b1) && (a1 <= b2) && (a2 <= b2))
-		return (Min(b1, b2) - Max(a1, a2));
+	NDBOX	   *a = PG_GETARG_NDBOX(0),
+			   *b = PG_GETARG_NDBOX(1);
+	bool		swapped = false;
+	double		d, distance;
+	int			i;
 
-	/* interval (a) is entirely on the right of (b) */
-	if ((a1 > b1) && (a2 > b1) && (a1 > b2) && (a2 > b2))
-		return (Min(a1, a2) - Max(b1, b2));
+	/* swap the box pointers if needed */
+	if (a->dim < b->dim)
+	{
+		NDBOX	   *tmp = b;
+		b = a;
+		a = tmp;
+		swapped = true;
+	}
 
-	/* the rest are all sorts of intersections */
-	return (0.0);
+	distance = 0.0;
+	/* compute within the dimensions of (b) */
+	for (i = 0; i < b->dim; i++)
+	{
+		d = abs(distance_1D(a->x[i], a->x[i + a->dim], b->x[i], b->x[i + b->dim]));
+		if (d > distance)
+			distance = d;
+	}
+
+	/* compute distance to zero for those dimensions in (a) absent in (b) */
+	for (i = b->dim; i < a->dim; i++)
+	{
+		d = abs(distance_1D(a->x[i], a->x[i + a->dim], 0.0, 0.0));
+		if (d > distance)
+			distance = d;
+	}
+
+	if (swapped)
+	{
+		PG_FREE_IF_COPY(b, 0);
+		PG_FREE_IF_COPY(a, 1);
+	}
+	else
+	{
+		PG_FREE_IF_COPY(a, 0);
+		PG_FREE_IF_COPY(b, 1);
+	}
+
+	PG_RETURN_FLOAT8(distance);
+}
+
+Datum
+g_cube_distance(PG_FUNCTION_ARGS)
+{
+	GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
+	StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
+	NDBOX      *cube = DatumGetNDBOX(entry->key);
+	double      retval;
+
+	if (strategy == 15)
+	{
+		int coord = PG_GETARG_INT32(1);
+		if(coord > 0)
+			retval =  cube->x[coord-1];
+		else
+			retval = -cube->x[-coord-1+cube->dim];
+	}
+	else
+	{
+		NDBOX *query = PG_GETARG_NDBOX(1);
+		switch(strategy)
+		{
+		case 16:
+			retval = DatumGetFloat8(DirectFunctionCall2(distance_taxicab, PointerGetDatum(cube), PointerGetDatum(query)));
+			break;
+		case 17:
+			retval = DatumGetFloat8(DirectFunctionCall2(distance_euclid, PointerGetDatum(cube), PointerGetDatum(query)));
+			break;
+		case 18:
+			retval = DatumGetFloat8(DirectFunctionCall2(distance_chebyshev, PointerGetDatum(cube), PointerGetDatum(query)));
+			break;
+		}
+	}
+	PG_RETURN_FLOAT8(retval);
 }
 
 /* Test if a box is also a point */
@@ -1352,6 +1490,22 @@ cube_ur_coord(PG_FUNCTION_ARGS)
 	PG_RETURN_FLOAT8(result);
 }
 
+Datum
+cube_coord(PG_FUNCTION_ARGS)
+{
+	NDBOX	   *cube = PG_GETARG_NDBOX(0);
+	int			coord = PG_GETARG_INT16(1);
+
+	if ((-2*cube->dim <= coord) && (coord < 0))
+		PG_RETURN_FLOAT8(-cube->x[-coord - 1]);
+	else if ((0 < coord) && (coord <= 2*cube->dim))
+		PG_RETURN_FLOAT8(cube->x[coord - 1]);
+	else
+		ereport(ERROR,
+					(errcode(ERRCODE_ARRAY_ELEMENT_ERROR),
+					 errmsg("Index out of bounds")));
+}
+
 /* Increase or decrease box size by a radius in at least n dimensions. */
 Datum
 cube_enlarge(PG_FUNCTION_ARGS)
