diff --git a/Makefile b/Makefile
index c9c225f..14350ec 100644
--- a/Makefile
+++ b/Makefile
@@ -34,12 +34,13 @@ DATA_built  = $(RELEASE_SQL) \
 			  pg_sphere--1.3.0--1.3.1.sql \
 			  pg_sphere--1.3.1--1.4.0.sql \
 			  pg_sphere--1.4.0--1.4.1.sql \
-			  pg_sphere--1.4.1--1.4.2.sql
+			  pg_sphere--1.4.1--1.4.2.sql \
+			  pg_sphere--1.4.2--1.5.0.sql
 
 DOCS        = README.pg_sphere COPYRIGHT.pg_sphere
 TESTS       = version tables points euler circle line ellipse poly path box \
               index contains_ops contains_ops_compat bounding_box_gist gnomo \
-              epochprop contains overlaps spoint_brin sbox_brin selectivity
+              epochprop contains overlaps spoint_brin sbox_brin selectivity knn
 REGRESS     = init $(TESTS)
 
 PG_CFLAGS	+= -DPGSPHERE_VERSION=$(PGSPHERE_VERSION)
@@ -209,6 +210,9 @@ pg_sphere--1.3.1--1.4.0.sql: pgs_circle_sel.sql.in pgs_hash.sql.in
 pg_sphere--1.4.0--1.4.1.sql pg_sphere--1.4.1--1.4.2.sql:
 	@echo "-- Nothing to upgrade in the schema" > $@
 
+pg_sphere--1.4.2--1.5.0.sql:
+	cat upgrade_scripts/$@.in $^ > $@
+
 # end of local stuff
 
 src/sscan.o : src/sparse.c
diff --git a/Makefile.common.mk b/Makefile.common.mk
index 9a73710..7764d84 100644
--- a/Makefile.common.mk
+++ b/Makefile.common.mk
@@ -5,4 +5,4 @@
 #----------------------------------------------------------------------------
 
 EXTENSION        	:= pg_sphere
-PGSPHERE_VERSION 	:= 1.4.2
+PGSPHERE_VERSION 	:= 1.5.0
diff --git a/doc/indices.sgm b/doc/indices.sgm
index 998b546..6b874b0 100644
--- a/doc/indices.sgm
+++ b/doc/indices.sgm
@@ -68,6 +68,11 @@
             </para>
           </listitem>
         </itemizedlist>
+        <para>
+          A GiST index can be also used for quickly finding the points closest to the given one
+          when ordering by an expression with the <literal>&lt;-&gt;</literal> operator,
+          as shown in an example below.
+        </para>
         <para>
           BRIN indexing supports just spherical points (<type>spoint</type>)
           and spherical coordinates range (<type>sbox</type>) at the moment.
@@ -82,6 +87,13 @@
 <![CDATA[CREATE INDEX test_pos_idx ON test USING GIST (pos);]]>
 <![CDATA[VACUUM ANALYZE test;]]>
           </programlisting>
+          <para>
+           To find the points closest to a given spherical position, use the <literal>&lt;-&gt;</literal> operator:
+          </para>
+          <programlisting>
+<![CDATA[SELECT * FROM test ORDER BY pos <-> spoint (0.2, 0.3) LIMIT 10 ]]>
+          </programlisting>
+
           <para>
             BRIN index can be created through the following syntax:
           </para>
diff --git a/expected/knn.out b/expected/knn.out
new file mode 100644
index 0000000..49c8903
--- /dev/null
+++ b/expected/knn.out
@@ -0,0 +1,122 @@
+CREATE TABLE points (id int, p spoint, pos int);
+INSERT INTO points (id, p) SELECT x, spoint(random()*6.28, (2*random()-1)*1.57) FROM generate_series(1,314159) x;
+CREATE INDEX i ON points USING gist (p);
+SET enable_indexscan = true;
+EXPLAIN (costs off) SELECT p <-> spoint (0.2, 0.3) FROM points ORDER BY 1 LIMIT 100;
+                   QUERY PLAN                    
+-------------------------------------------------
+ Limit
+   ->  Index Scan using i on points
+         Order By: (p <-> '(0.2 , 0.3)'::spoint)
+(3 rows)
+
+UPDATE points SET pos = n FROM
+  (SELECT id, row_number() OVER (ORDER BY p <-> spoint (0.2, 0.3)) n  FROM points ORDER BY p <-> spoint (0.2, 0.3) LIMIT 100) sel 
+  WHERE points.id = sel.id;
+SET enable_indexscan = false;
+SELECT pos, row_number() OVER (ORDER BY p <-> spoint (0.2, 0.3)) n  FROM points ORDER BY p <-> spoint (0.2, 0.3) LIMIT 100;
+ pos |  n  
+-----+-----
+   1 |   1
+   2 |   2
+   3 |   3
+   4 |   4
+   5 |   5
+   6 |   6
+   7 |   7
+   8 |   8
+   9 |   9
+  10 |  10
+  11 |  11
+  12 |  12
+  13 |  13
+  14 |  14
+  15 |  15
+  16 |  16
+  17 |  17
+  18 |  18
+  19 |  19
+  20 |  20
+  21 |  21
+  22 |  22
+  23 |  23
+  24 |  24
+  25 |  25
+  26 |  26
+  27 |  27
+  28 |  28
+  29 |  29
+  30 |  30
+  31 |  31
+  32 |  32
+  33 |  33
+  34 |  34
+  35 |  35
+  36 |  36
+  37 |  37
+  38 |  38
+  39 |  39
+  40 |  40
+  41 |  41
+  42 |  42
+  43 |  43
+  44 |  44
+  45 |  45
+  46 |  46
+  47 |  47
+  48 |  48
+  49 |  49
+  50 |  50
+  51 |  51
+  52 |  52
+  53 |  53
+  54 |  54
+  55 |  55
+  56 |  56
+  57 |  57
+  58 |  58
+  59 |  59
+  60 |  60
+  61 |  61
+  62 |  62
+  63 |  63
+  64 |  64
+  65 |  65
+  66 |  66
+  67 |  67
+  68 |  68
+  69 |  69
+  70 |  70
+  71 |  71
+  72 |  72
+  73 |  73
+  74 |  74
+  75 |  75
+  76 |  76
+  77 |  77
+  78 |  78
+  79 |  79
+  80 |  80
+  81 |  81
+  82 |  82
+  83 |  83
+  84 |  84
+  85 |  85
+  86 |  86
+  87 |  87
+  88 |  88
+  89 |  89
+  90 |  90
+  91 |  91
+  92 |  92
+  93 |  93
+  94 |  94
+  95 |  95
+  96 |  96
+  97 |  97
+  98 |  98
+  99 |  99
+ 100 | 100
+(100 rows)
+
+DROP TABLE points;
diff --git a/expected/version.out b/expected/version.out
index 5e361bf..1d82d75 100644
--- a/expected/version.out
+++ b/expected/version.out
@@ -2,6 +2,6 @@
 SELECT pg_sphere_version();
  pg_sphere_version 
 -------------------
- 1.4.2
+ 1.5.0
 (1 row)
 
diff --git a/pg_sphere.control b/pg_sphere.control
index 69ca772..4ecf4cd 100644
--- a/pg_sphere.control
+++ b/pg_sphere.control
@@ -1,5 +1,5 @@
 # pg_sphere extension
 comment = 'spherical objects with useful functions, operators and index support'
-default_version = '1.4.2'
+default_version = '1.5.0'
 module_pathname = '$libdir/pg_sphere'
 relocatable = true
diff --git a/pgs_gist.sql.in b/pgs_gist.sql.in
index e8ac119..fcd4f11 100644
--- a/pgs_gist.sql.in
+++ b/pgs_gist.sql.in
@@ -98,12 +98,15 @@ CREATE FUNCTION g_spoint_compress(internal)
    AS 'MODULE_PATHNAME', 'g_spoint_compress'
    LANGUAGE 'c';
 
-
 CREATE FUNCTION g_spoint_consistent(internal, internal, int4, oid, internal)
    RETURNS internal
    AS 'MODULE_PATHNAME', 'g_spoint_consistent'
    LANGUAGE 'c';
 
+CREATE FUNCTION g_spoint_distance(internal, spoint, smallint, oid, internal)
+   RETURNS internal
+   AS 'MODULE_PATHNAME', 'g_spoint_distance'
+   LANGUAGE 'c';
 
 CREATE OPERATOR CLASS spoint
    DEFAULT FOR TYPE spoint USING gist AS
@@ -114,6 +117,7 @@ CREATE OPERATOR CLASS spoint
    OPERATOR  14 @ (spoint, spoly),
    OPERATOR  15 @ (spoint, sellipse),
    OPERATOR  16 @ (spoint, sbox),
+   OPERATOR  17 <-> (spoint, spoint) FOR ORDER BY float_ops,
    OPERATOR  37 <@ (spoint, scircle),
    OPERATOR  38 <@ (spoint, sline),
    OPERATOR  39 <@ (spoint, spath),
@@ -127,6 +131,7 @@ CREATE OPERATOR CLASS spoint
    FUNCTION  5 g_spherekey_penalty (internal, internal, internal),
    FUNCTION  6 g_spherekey_picksplit (internal, internal),
    FUNCTION  7 g_spherekey_same (spherekey, spherekey, internal),
+   FUNCTION  8 g_spoint_distance (internal, spoint, smallint, oid, internal),
    STORAGE   spherekey;
 
 
diff --git a/sql/knn.sql b/sql/knn.sql
new file mode 100644
index 0000000..0207c98
--- /dev/null
+++ b/sql/knn.sql
@@ -0,0 +1,12 @@
+CREATE TABLE points (id int, p spoint, pos int);
+INSERT INTO points (id, p) SELECT x, spoint(random()*6.28, (2*random()-1)*1.57) FROM generate_series(1,314159) x;
+CREATE INDEX i ON points USING gist (p);
+SET enable_indexscan = true;
+EXPLAIN (costs off) SELECT p <-> spoint (0.2, 0.3) FROM points ORDER BY 1 LIMIT 100;
+UPDATE points SET pos = n FROM
+  (SELECT id, row_number() OVER (ORDER BY p <-> spoint (0.2, 0.3)) n  FROM points ORDER BY p <-> spoint (0.2, 0.3) LIMIT 100) sel 
+  WHERE points.id = sel.id;
+SET enable_indexscan = false;
+SELECT pos, row_number() OVER (ORDER BY p <-> spoint (0.2, 0.3)) n  FROM points ORDER BY p <-> spoint (0.2, 0.3) LIMIT 100;
+DROP TABLE points;
+
diff --git a/src/gist.c b/src/gist.c
index b3223a5..0fef2ac 100644
--- a/src/gist.c
+++ b/src/gist.c
@@ -36,6 +36,7 @@ PG_FUNCTION_INFO_V1(g_spherekey_penalty);
 PG_FUNCTION_INFO_V1(g_spherekey_picksplit);
 PG_FUNCTION_INFO_V1(g_spoint3_penalty);
 PG_FUNCTION_INFO_V1(g_spoint3_picksplit);
+PG_FUNCTION_INFO_V1(g_spoint_distance);
 PG_FUNCTION_INFO_V1(g_spoint3_distance);
 PG_FUNCTION_INFO_V1(g_spoint3_fetch);
 
@@ -510,7 +511,7 @@ g_spoint_consistent(PG_FUNCTION_ARGS)
 				SCK_INTERLEAVE(SELLIPSE, sphereellipse_gen_key, 0);
 				break;
 			case 42:
-				SCK_INTERLEAVE(SBOX, spherebox_gen_key , 0);
+				SCK_INTERLEAVE(SBOX, spherebox_gen_key, 0);
 				break;
 		}
 
@@ -681,6 +682,13 @@ g_spoint3_consistent(PG_FUNCTION_ARGS)
 	PG_RETURN_BOOL(false);
 }
 
+static double
+distance_vector_point_3d(Vector3D *v, double x, double y, double z)
+{
+	/* as v has length = 1 by design */
+	return acos((v->x * x + v->y * y + v->z * z) / sqrt(x * x + y * y + z * z));
+}
+
 Datum
 g_spoint3_distance(PG_FUNCTION_ARGS)
 {
@@ -1672,6 +1680,224 @@ fallbackSplit(Box3D *boxes, OffsetNumber maxoff, GIST_SPLITVEC *v)
 	v->spl_ldatum_exists = v->spl_rdatum_exists = false;
 }
 
+
+Datum
+g_spoint_distance(PG_FUNCTION_ARGS)
+{
+	GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
+	StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
+	Box3D	   *box = (Box3D *) DatumGetPointer(entry->key);
+	double		retval;
+	SPoint	   *point = (SPoint *) PG_GETARG_POINTER(1);
+	Vector3D	v_point,
+				v_low,
+				v_high;
+
+	switch (strategy)
+	{
+		case 17:
+			/* Prepare data for calculation */
+			spoint_vector3d(&v_point, point);
+			v_low.x = (double) box->low.coord[0] / MAXCVALUE;
+			v_low.y = (double) box->low.coord[1] / MAXCVALUE;
+			v_low.z = (double) box->low.coord[2] / MAXCVALUE;
+			v_high.x = (double) box->high.coord[0] / MAXCVALUE;
+			v_high.y = (double) box->high.coord[1] / MAXCVALUE;
+			v_high.z = (double) box->high.coord[2] / MAXCVALUE;
+
+			/*
+			 * a box splits space into 27 subspaces (6+12+8+1) with different
+			 * distance calculation
+			 */
+			if (v_point.x < v_low.x)
+			{
+				if (v_point.y < v_low.y)
+				{
+					if (v_point.z < v_low.z)
+					{
+						/* point2point distance */
+						retval = distance_vector_point_3d(&v_point, v_low.x, v_low.y, v_low.z);
+					}
+					else if (v_point.z < v_high.z)
+					{
+						/* point2line distance */
+						retval = distance_vector_point_3d(&v_point, v_low.x, v_low.y, v_point.z);
+					}
+					else
+					{
+						/* point2point distance */
+						retval = distance_vector_point_3d(&v_point, v_low.x, v_low.y, v_high.z);
+					}
+				}
+				else if (v_point.y < v_high.y)
+				{
+					if (v_point.z < v_low.z)
+					{
+						/* point2line distance */
+						retval = distance_vector_point_3d(&v_point, v_low.x, v_point.y, v_low.z);
+					}
+					else if (v_point.z < v_high.z)
+					{
+						/* point2plane distance */
+						retval = distance_vector_point_3d(&v_point, v_low.x, v_point.y, v_point.z);
+					}
+					else
+					{
+						/* point2line distance */
+						retval = distance_vector_point_3d(&v_point, v_low.x, v_point.y, v_high.z);
+					}
+				}
+				else
+				{
+					if (v_point.z < v_low.z)
+					{
+						/* point2point distance */
+						retval = distance_vector_point_3d(&v_point, v_low.x, v_high.y, v_low.z);
+					}
+					else if (v_point.z < v_high.z)
+					{
+						/* point2line distance */
+						retval = distance_vector_point_3d(&v_point, v_low.x, v_high.y, v_point.z);
+					}
+					else
+					{
+						/* point2point distance */
+						retval = distance_vector_point_3d(&v_point, v_low.x, v_high.y, v_high.z);
+					}
+				}
+			}
+			else if (v_point.x < v_high.x)
+			{
+				if (v_point.y < v_low.y)
+				{
+					if (v_point.z < v_low.z)
+					{
+						/* p2line distance */
+						retval = distance_vector_point_3d(&v_point, v_point.x, v_low.y, v_low.z);
+					}
+					else if (v_point.z < v_high.z)
+					{
+						/* point2plane distance */
+						retval = distance_vector_point_3d(&v_point, v_point.x, v_low.y, v_point.z);
+					}
+					else
+					{
+						/* point2line distance */
+						retval = distance_vector_point_3d(&v_point, v_point.x, v_low.y, v_high.z);
+					}
+				}
+				else if (v_point.y < v_high.y)
+				{
+					if (v_point.z < v_low.z)
+					{
+						/* point2plane distance */
+						retval = distance_vector_point_3d(&v_point, v_point.x, v_point.y, v_low.z);
+					}
+					else if (v_point.z < v_high.z)
+					{
+						/* inside cube */
+						retval = 0;
+					}
+					else
+					{
+						/* point2plane distance */
+						retval = distance_vector_point_3d(&v_point, v_point.x, v_point.y, v_high.z);
+					}
+				}
+				else
+				{
+					if (v_point.z < v_low.z)
+					{
+						/* point2line distance */
+						retval = distance_vector_point_3d(&v_point, v_point.x, v_high.y, v_low.z);
+					}
+					else if (v_point.z < v_high.z)
+					{
+						/* point2plane distance */
+						retval = distance_vector_point_3d(&v_point, v_point.x, v_high.y, v_point.z);
+					}
+					else
+					{
+						/* point2line distance */
+						retval = distance_vector_point_3d(&v_point, v_point.x, v_high.y, v_high.z);
+					}
+				}
+			}
+			else
+			{
+				if (v_point.y < v_low.y)
+				{
+					if (v_point.z < v_low.z)
+					{
+						/* p2p distance */
+						retval = distance_vector_point_3d(&v_point, v_high.x, v_low.y, v_low.z);
+					}
+					else if (v_point.z < v_high.z)
+					{
+						/* point2line distance */
+						retval = distance_vector_point_3d(&v_point, v_high.x, v_low.y, v_point.z);
+					}
+					else
+					{
+						/* point2point distance */
+						retval = distance_vector_point_3d(&v_point, v_high.x, v_low.y, v_high.z);
+					}
+				}
+				else if (v_point.y < v_high.y)
+				{
+					if (v_point.z < v_low.z)
+					{
+						/* point2line distance */
+						retval = distance_vector_point_3d(&v_point, v_high.x, v_point.y, v_low.z);
+					}
+					else if (v_point.z < v_high.z)
+					{
+						/* point2plane distance */
+						retval = distance_vector_point_3d(&v_point, v_high.x, v_point.y, v_point.z);
+					}
+					else
+					{
+						/* point2line distance */
+						retval = distance_vector_point_3d(&v_point, v_high.x, v_point.y, v_high.z);
+					}
+				}
+				else
+				{
+					if (v_point.z < v_low.z)
+					{
+						/* point2point distance */
+						retval = distance_vector_point_3d(&v_point, v_high.x, v_high.y, v_low.z);
+					}
+					else if (v_point.z < v_high.z)
+					{
+						/* point2line distance */
+						retval = distance_vector_point_3d(&v_point, v_high.x, v_high.y, v_point.z);
+					}
+					else
+					{
+						/* point2point distance */
+						retval = distance_vector_point_3d(&v_point, v_high.x, v_high.y, v_high.z);
+					}
+				}
+			}
+
+			elog(DEBUG1, "distance (%lg,%lg,%lg %lg,%lg,%lg) <-> (%lg,%lg) = %lg",
+				 v_low.x, v_low.y, v_low.z,
+				 v_high.x, v_high.y, v_high.z,
+				 point->lng, point->lat,
+				 retval
+				);
+			break;
+
+		default:
+			elog(ERROR, "unrecognized cube strategy number: %d", strategy);
+			retval = 0;			/* keep compiler quiet */
+			break;
+	}
+
+	PG_RETURN_FLOAT8(retval);
+}
+
 /*
  * Represents information about an entry that can be placed to either group
  * without affecting overlap over selected axis ("common entry").
@@ -2203,7 +2429,7 @@ do_picksplit(Box3D *boxes, OffsetNumber maxoff, GIST_SPLITVEC *v)
 		{
 			box = &boxes[i];
 			commonEntries[i].delta = fabs((unionSizeBox3D(leftBox, box) - leftBoxSize) -
-							 (unionSizeBox3D(rightBox, box) - rightBoxSize));
+										  (unionSizeBox3D(rightBox, box) - rightBoxSize));
 		}
 
 		/*
diff --git a/upgrade_scripts/pg_sphere--1.4.2--1.5.0.sql.in b/upgrade_scripts/pg_sphere--1.4.2--1.5.0.sql.in
new file mode 100644
index 0000000..35c964e
--- /dev/null
+++ b/upgrade_scripts/pg_sphere--1.4.2--1.5.0.sql.in
@@ -0,0 +1,16 @@
+-- Upgrade: 1.4.2 -> 1.5.0
+
+CREATE OR REPLACE FUNCTION g_spoint_distance(internal, spoint, smallint, oid, internal)
+   RETURNS internal
+   AS 'MODULE_PATHNAME', 'g_spoint_distance'
+   LANGUAGE 'c';
+
+DO $$
+BEGIN
+   ALTER OPERATOR FAMILY spoint USING gist ADD
+      FUNCTION  8 (spoint, spoint) g_spoint_distance (internal, spoint, smallint, oid, internal);
+EXCEPTION
+   WHEN duplicate_object THEN NULL;
+   WHEN OTHERS THEN RAISE;
+END;
+$$;