From 11eec2d91d58875f749c5800c2e85738b2edc932 Mon Sep 17 00:00:00 2001
From: Emre Hasegeli <emre@hasegeli.com>
Date: Sat, 27 May 2017 16:17:42 -0400
Subject: [PATCH 3/6] geo-float-v09

Use the built-in float datatype to implement geometric types

This patch makes the geometric operators and functions use the exported
function of the float datatype.  The main reason of doing so is to check
for underflow and overflow, and to handle NaNs consciously.

The float datatypes consider NaNs to be equal and greater than any
non-NaN.  This change considers NaNs equal only for the equality
operators.  The placement operators, contains, overlaps, left/right of
etc. continues to return false when NaNs are involved.  We don't need
to worry about them being considered greater than any-NaN because there
aren't any basic comparison operators like less/greater than for the
geometric datatypes.

All changes can be summarised as:

* Check for underflow and overflow
* Check for division by zero
* Let the objects with NaN values to be the same
* Return NULL when the distance is NaN for all closest point operators
* Don't round the slope to DBL_MAX
* Favour not-NaN over NaN where it makes sense

The patch also replaces all occurrences of "double" as "float8".  They
are the same, but were randomly spread on the same file.

Discussion: https://www.postgresql.org/message-id/CAE2gYzxF7-5djV6-cEvqQu-fNsnt%3DEqbOURx7ZDg%2BVv6ZMTWbg%40mail.gmail.com

Parallel discussions:

* https://www.postgresql.org/message-id/28685.1468246504%40sss.pgh.pa.us
---
 src/backend/access/gist/gistproc.c |  91 +++----
 src/backend/utils/adt/geo_ops.c    | 497 ++++++++++++++++++++-----------------
 src/backend/utils/adt/geo_spgist.c |  28 +--
 src/include/utils/geo_decls.h      |  17 +-
 4 files changed, 340 insertions(+), 293 deletions(-)

diff --git a/src/backend/access/gist/gistproc.c b/src/backend/access/gist/gistproc.c
index 4bbfdadf9f..4b6233f212 100644
--- a/src/backend/access/gist/gistproc.c
+++ b/src/backend/access/gist/gistproc.c
@@ -16,20 +16,21 @@
  *-------------------------------------------------------------------------
  */
 #include "postgres.h"
 
 #include <float.h>
 #include <math.h>
 
 #include "access/gist.h"
 #include "access/stratnum.h"
 #include "utils/builtins.h"
+#include "utils/float.h"
 #include "utils/geo_decls.h"
 
 
 static bool gist_box_leaf_consistent(BOX *key, BOX *query,
 						 StrategyNumber strategy);
 static bool rtree_internal_consistent(BOX *key, BOX *query,
 						  StrategyNumber strategy);
 
 /* Minimum accepted ratio of split */
 #define LIMIT_RATIO 0.3
@@ -48,55 +49,56 @@ rt_box_union(BOX *n, const BOX *a, const BOX *b)
 	n->high.x = float8_max(a->high.x, b->high.x);
 	n->high.y = float8_max(a->high.y, b->high.y);
 	n->low.x = float8_min(a->low.x, b->low.x);
 	n->low.y = float8_min(a->low.y, b->low.y);
 }
 
 /*
  * Size of a BOX for penalty-calculation purposes.
  * The result can be +Infinity, but not NaN.
  */
-static double
+static float8
 size_box(const BOX *box)
 {
 	/*
 	 * Check for zero-width cases.  Note that we define the size of a zero-
 	 * by-infinity box as zero.  It's important to special-case this somehow,
 	 * as naively multiplying infinity by zero will produce NaN.
 	 *
 	 * The less-than cases should not happen, but if they do, say "zero".
 	 */
 	if (float8_le(box->high.x, box->low.x) ||
 		float8_le(box->high.y, box->low.y))
 		return 0.0;
 
 	/*
 	 * We treat NaN as larger than +Infinity, so any distance involving a NaN
 	 * and a non-NaN is infinite.  Note the previous check eliminated the
 	 * possibility that the low fields are NaNs.
 	 */
 	if (isnan(box->high.x) || isnan(box->high.y))
 		return get_float8_infinity();
-	return (box->high.x - box->low.x) * (box->high.y - box->low.y);
+	return float8_mul(float8_mi(box->high.x, box->low.x),
+					  float8_mi(box->high.y, box->low.y));
 }
 
 /*
  * Return amount by which the union of the two boxes is larger than
  * the original BOX's area.  The result can be +Infinity, but not NaN.
  */
-static double
+static float8
 box_penalty(const BOX *original, const BOX *new)
 {
 	BOX			unionbox;
 
 	rt_box_union(&unionbox, original, new);
-	return size_box(&unionbox) - size_box(original);
+	return float8_mi(size_box(&unionbox), size_box(original));
 }
 
 /*
  * The GiST Consistent method for boxes
  *
  * Should return false if for all data items x below entry,
  * the predicate x op query must be false, where op is the oper
  * corresponding to strategy in the pg_amop table.
  */
 Datum
@@ -256,74 +258,74 @@ fallbackSplit(GistEntryVector *entryvec, GIST_SPLITVEC *v)
 
 /*
  * Represents information about an entry that can be placed to either group
  * without affecting overlap over selected axis ("common entry").
  */
 typedef struct
 {
 	/* Index of entry in the initial array */
 	int			index;
 	/* Delta between penalties of entry insertion into different groups */
-	double		delta;
+	float8		delta;
 } CommonEntry;
 
 /*
  * Context for g_box_consider_split. Contains information about currently
  * selected split and some general information.
  */
 typedef struct
 {
 	int			entriesCount;	/* total number of entries being split */
 	BOX			boundingBox;	/* minimum bounding box across all entries */
 
 	/* Information about currently selected split follows */
 
 	bool		first;			/* true if no split was selected yet */
 
-	double		leftUpper;		/* upper bound of left interval */
-	double		rightLower;		/* lower bound of right interval */
+	float8		leftUpper;		/* upper bound of left interval */
+	float8		rightLower;		/* lower bound of right interval */
 
 	float4		ratio;
 	float4		overlap;
 	int			dim;			/* axis of this split */
-	double		range;			/* width of general MBR projection to the
+	float8		range;			/* width of general MBR projection to the
 								 * selected axis */
 } ConsiderSplitContext;
 
 /*
  * Interval represents projection of box to axis.
  */
 typedef struct
 {
-	double		lower,
+	float8		lower,
 				upper;
 } SplitInterval;
 
 /*
  * Interval comparison function by lower bound of the interval;
  */
 static int
 interval_cmp_lower(const void *i1, const void *i2)
 {
-	double		lower1 = ((const SplitInterval *) i1)->lower,
+	float8		lower1 = ((const SplitInterval *) i1)->lower,
 				lower2 = ((const SplitInterval *) i2)->lower;
 
 	return float8_cmp_internal(lower1, lower2);
 }
 
 /*
  * Interval comparison function by upper bound of the interval;
  */
 static int
 interval_cmp_upper(const void *i1, const void *i2)
 {
-	double		upper1 = ((const SplitInterval *) i1)->upper,
+	float8		upper1 = ((const SplitInterval *) i1)->upper,
 				upper2 = ((const SplitInterval *) i2)->upper;
 
 	return float8_cmp_internal(upper1, upper2);
 }
 
 /*
  * Replace negative (or NaN) value with zero.
  */
 static inline float
 non_negative(float val)
@@ -332,28 +334,28 @@ non_negative(float val)
 		return val;
 	else
 		return 0.0f;
 }
 
 /*
  * Consider replacement of currently selected split with the better one.
  */
 static inline void
 g_box_consider_split(ConsiderSplitContext *context, int dimNum,
-					 double rightLower, int minLeftCount,
-					 double leftUpper, int maxLeftCount)
+					 float8 rightLower, int minLeftCount,
+					 float8 leftUpper, int maxLeftCount)
 {
 	int			leftCount,
 				rightCount;
 	float4		ratio,
 				overlap;
-	double		range;
+	float8		range;
 
 	/*
 	 * Calculate entries distribution ratio assuming most uniform distribution
 	 * of common entries.
 	 */
 	if (minLeftCount >= (context->entriesCount + 1) / 2)
 	{
 		leftCount = minLeftCount;
 	}
 	else
@@ -362,40 +364,41 @@ g_box_consider_split(ConsiderSplitContext *context, int dimNum,
 			leftCount = maxLeftCount;
 		else
 			leftCount = context->entriesCount / 2;
 	}
 	rightCount = context->entriesCount - leftCount;
 
 	/*
 	 * Ratio of split - quotient between size of lesser group and total
 	 * entries count.
 	 */
-	ratio = ((float4) Min(leftCount, rightCount)) /
-		((float4) context->entriesCount);
+	ratio = float4_div(Min(leftCount, rightCount), context->entriesCount);
 
 	if (ratio > LIMIT_RATIO)
 	{
 		bool		selectthis = false;
 
 		/*
 		 * The ratio is acceptable, so compare current split with previously
 		 * selected one. Between splits of one dimension we search for minimal
 		 * overlap (allowing negative values) and minimal ration (between same
 		 * overlaps. We switch dimension if find less overlap (non-negative)
 		 * or less range with same overlap.
 		 */
 		if (dimNum == 0)
-			range = context->boundingBox.high.x - context->boundingBox.low.x;
+			range = float8_mi(context->boundingBox.high.x,
+							  context->boundingBox.low.x);
 		else
-			range = context->boundingBox.high.y - context->boundingBox.low.y;
+			range = float8_mi(context->boundingBox.high.y,
+							  context->boundingBox.low.y);
 
-		overlap = (leftUpper - rightLower) / range;
+		overlap = float8_div(float8_mi(leftUpper, rightLower), range);
 
 		/* If there is no previous selection, select this */
 		if (context->first)
 			selectthis = true;
 		else if (context->dim == dimNum)
 		{
 			/*
 			 * Within the same dimension, choose the new split if it has a
 			 * smaller overlap, or same overlap but better ratio.
 			 */
@@ -524,21 +527,21 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
 		else
 			adjustBox(&context.boundingBox, box);
 	}
 
 	/*
 	 * Iterate over axes for optimal split searching.
 	 */
 	context.first = true;		/* nothing selected yet */
 	for (dim = 0; dim < 2; dim++)
 	{
-		double		leftUpper,
+		float8		leftUpper,
 					rightLower;
 		int			i1,
 					i2;
 
 		/* Project each entry as an interval on the selected axis. */
 		for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
 		{
 			box = DatumGetBoxP(entryvec->vector[i].key);
 			if (dim == 0)
 			{
@@ -721,21 +724,21 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
 			*rightBox = *(box);					\
 		v->spl_right[v->spl_nright++] = off;	\
 	} while(0)
 
 	/*
 	 * Distribute entries which can be distributed unambiguously, and collect
 	 * common entries.
 	 */
 	for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
 	{
-		double		lower,
+		float8		lower,
 					upper;
 
 		/*
 		 * Get upper and lower bounds along selected axis.
 		 */
 		box = DatumGetBoxP(entryvec->vector[i].key);
 		if (context.dim == 0)
 		{
 			lower = box->low.x;
 			upper = box->high.x;
@@ -776,31 +779,31 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
 
 	/*
 	 * Distribute "common entries", if any.
 	 */
 	if (commonEntriesCount > 0)
 	{
 		/*
 		 * Calculate minimum number of entries that must be placed in both
 		 * groups, to reach LIMIT_RATIO.
 		 */
-		int			m = ceil(LIMIT_RATIO * (double) nentries);
+		int			m = ceil(LIMIT_RATIO * (float8) nentries);
 
 		/*
 		 * Calculate delta between penalties of join "common entries" to
 		 * different groups.
 		 */
 		for (i = 0; i < commonEntriesCount; i++)
 		{
 			box = DatumGetBoxP(entryvec->vector[commonEntries[i].index].key);
-			commonEntries[i].delta = Abs(box_penalty(leftBox, box) -
-										 box_penalty(rightBox, box));
+			commonEntries[i].delta = Abs(float8_mi(box_penalty(leftBox, box),
+												   box_penalty(rightBox, box)));
 		}
 
 		/*
 		 * Sort "common entries" by calculated deltas in order to distribute
 		 * the most ambiguous entries first.
 		 */
 		qsort(commonEntries, commonEntriesCount, sizeof(CommonEntry), common_entry_cmp);
 
 		/*
 		 * Distribute "common entries" between groups.
@@ -1100,24 +1103,24 @@ gist_circle_compress(PG_FUNCTION_ARGS)
 {
 	GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
 	GISTENTRY  *retval;
 
 	if (entry->leafkey)
 	{
 		CIRCLE	   *in = DatumGetCircleP(entry->key);
 		BOX		   *r;
 
 		r = (BOX *) palloc(sizeof(BOX));
-		r->high.x = in->center.x + in->radius;
-		r->low.x = in->center.x - in->radius;
-		r->high.y = in->center.y + in->radius;
-		r->low.y = in->center.y - in->radius;
+		r->high.x = float8_pl(in->center.x, in->radius);
+		r->low.x = float8_mi(in->center.x, in->radius);
+		r->high.y = float8_pl(in->center.y, in->radius);
+		r->low.y = float8_mi(in->center.y, in->radius);
 
 		retval = (GISTENTRY *) palloc(sizeof(GISTENTRY));
 		gistentryinit(*retval, PointerGetDatum(r),
 					  entry->rel, entry->page,
 					  entry->offset, false);
 	}
 	else
 		retval = entry;
 	PG_RETURN_POINTER(retval);
 }
@@ -1141,24 +1144,24 @@ gist_circle_consistent(PG_FUNCTION_ARGS)
 	*recheck = true;
 
 	if (DatumGetBoxP(entry->key) == NULL || query == NULL)
 		PG_RETURN_BOOL(false);
 
 	/*
 	 * Since the operators require recheck anyway, we can just use
 	 * rtree_internal_consistent even at leaf nodes.  (This works in part
 	 * because the index entries are bounding boxes not circles.)
 	 */
-	bbox.high.x = query->center.x + query->radius;
-	bbox.low.x = query->center.x - query->radius;
-	bbox.high.y = query->center.y + query->radius;
-	bbox.low.y = query->center.y - query->radius;
+	bbox.high.x = float8_pl(query->center.x, query->radius);
+	bbox.low.x = float8_mi(query->center.x, query->radius);
+	bbox.high.y = float8_pl(query->center.y, query->radius);
+	bbox.low.y = float8_mi(query->center.y, query->radius);
 
 	result = rtree_internal_consistent(DatumGetBoxP(entry->key),
 									   &bbox, strategy);
 
 	PG_RETURN_BOOL(result);
 }
 
 /**************************************************
  * Point ops
  **************************************************/
@@ -1209,63 +1212,63 @@ gist_point_fetch(PG_FUNCTION_ARGS)
 				  entry->offset, false);
 
 	PG_RETURN_POINTER(retval);
 }
 
 
 #define point_point_distance(p1,p2) \
 	DatumGetFloat8(DirectFunctionCall2(point_distance, \
 									   PointPGetDatum(p1), PointPGetDatum(p2)))
 
-static double
+static float8
 computeDistance(bool isLeaf, BOX *box, Point *point)
 {
-	double		result = 0.0;
+	float8		result = 0.0;
 
 	if (isLeaf)
 	{
 		/* simple point to point distance */
 		result = point_point_distance(point, &box->low);
 	}
 	else if (point->x <= box->high.x && point->x >= box->low.x &&
 			 point->y <= box->high.y && point->y >= box->low.y)
 	{
 		/* point inside the box */
 		result = 0.0;
 	}
 	else if (point->x <= box->high.x && point->x >= box->low.x)
 	{
 		/* point is over or below box */
 		Assert(box->low.y <= box->high.y);
 		if (point->y > box->high.y)
-			result = point->y - box->high.y;
+			result = float8_mi(point->y, box->high.y);
 		else if (point->y < box->low.y)
-			result = box->low.y - point->y;
+			result = float8_mi(box->low.y, point->y);
 		else
 			elog(ERROR, "inconsistent point values");
 	}
 	else if (point->y <= box->high.y && point->y >= box->low.y)
 	{
 		/* point is to left or right of box */
 		Assert(box->low.x <= box->high.x);
 		if (point->x > box->high.x)
-			result = point->x - box->high.x;
+			result = float8_mi(point->x, box->high.x);
 		else if (point->x < box->low.x)
-			result = box->low.x - point->x;
+			result = float8_mi(box->low.x, point->x);
 		else
 			elog(ERROR, "inconsistent point values");
 	}
 	else
 	{
 		/* closest point will be a vertex */
 		Point		p;
-		double		subresult;
+		float8		subresult;
 
 		result = point_point_distance(point, &box->low);
 
 		subresult = point_point_distance(point, &box->high);
 		if (result > subresult)
 			result = subresult;
 
 		p.x = box->low.x;
 		p.y = box->high.y;
 		subresult = point_point_distance(point, &p);
@@ -1442,21 +1445,21 @@ gist_point_consistent(PG_FUNCTION_ARGS)
 	}
 
 	PG_RETURN_BOOL(result);
 }
 
 Datum
 gist_point_distance(PG_FUNCTION_ARGS)
 {
 	GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
 	StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
-	double		distance;
+	float8		distance;
 	StrategyNumber strategyGroup = strategy / GeoStrategyNumberOffset;
 
 	switch (strategyGroup)
 	{
 		case PointStrategyNumberGroup:
 			distance = computeDistance(GIST_LEAF(entry),
 									   DatumGetBoxP(entry->key),
 									   PG_GETARG_POINT_P(1));
 			break;
 		default:
@@ -1471,25 +1474,25 @@ gist_point_distance(PG_FUNCTION_ARGS)
 /*
  * The inexact GiST distance method for geometric types that store bounding
  * boxes.
  *
  * Compute lossy distance from point to index entries.  The result is inexact
  * because index entries are bounding boxes, not the exact shapes of the
  * indexed geometric types.  We use distance from point to MBR of index entry.
  * This is a lower bound estimate of distance from point to indexed geometric
  * type.
  */
-static double
+static float8
 gist_bbox_distance(GISTENTRY *entry, Datum query,
 				   StrategyNumber strategy, bool *recheck)
 {
-	double		distance;
+	float8		distance;
 	StrategyNumber strategyGroup = strategy / GeoStrategyNumberOffset;
 
 	/* Bounding box distance is always inexact. */
 	*recheck = true;
 
 	switch (strategyGroup)
 	{
 		case PointStrategyNumberGroup:
 			distance = computeDistance(false,
 									   DatumGetBoxP(entry->key),
@@ -1505,32 +1508,32 @@ gist_bbox_distance(GISTENTRY *entry, Datum query,
 
 Datum
 gist_circle_distance(PG_FUNCTION_ARGS)
 {
 	GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
 	Datum		query = PG_GETARG_DATUM(1);
 	StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
 
 	/* Oid subtype = PG_GETARG_OID(3); */
 	bool	   *recheck = (bool *) PG_GETARG_POINTER(4);
-	double		distance;
+	float8		distance;
 
 	distance = gist_bbox_distance(entry, query, strategy, recheck);
 
 	PG_RETURN_FLOAT8(distance);
 }
 
 Datum
 gist_poly_distance(PG_FUNCTION_ARGS)
 {
 	GISTENTRY  *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
 	Datum		query = PG_GETARG_DATUM(1);
 	StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
 
 	/* Oid subtype = PG_GETARG_OID(3); */
 	bool	   *recheck = (bool *) PG_GETARG_POINTER(4);
-	double		distance;
+	float8		distance;
 
 	distance = gist_bbox_distance(entry, query, strategy, recheck);
 
 	PG_RETURN_FLOAT8(distance);
 }
diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index cdea309eeb..49667b17f7 100644
--- a/src/backend/utils/adt/geo_ops.c
+++ b/src/backend/utils/adt/geo_ops.c
@@ -51,56 +51,56 @@ static inline float8 line_invsl(LINE *line);
 static bool	line_interpt_line(Point *result, LINE *l1, LINE *l2);
 static bool line_contain_point(LINE *line, Point *point);
 static float8 line_closept_point(Point *result, LINE *line, Point *pt);
 
 /* Routines for line segments */
 static inline void statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2);
 static inline float8 lseg_sl(LSEG *lseg);
 static inline float8 lseg_invsl(LSEG *lseg);
 static bool	lseg_interpt_line(Point *result, LSEG *lseg, LINE *line);
 static bool	lseg_interpt_lseg(Point *result, LSEG *l1, LSEG *l2);
-static int	lseg_crossing(double x, double y, double px, double py);
+static int	lseg_crossing(float8 x, float8 y, float8 px, float8 py);
 static bool	lseg_contain_point(LSEG *lseg, Point *point);
 static float8 lseg_closept_point(Point *result, LSEG *lseg, Point *pt);
 static float8 lseg_closept_line(Point *result, LSEG *lseg, LINE *line);
 static float8 lseg_closept_lseg(Point *result, LSEG *l1, LSEG *l2);
 
 /* Routines for boxes */
 static inline void box_construct(BOX *result, Point *pt1, Point *pt2);
 static void box_cn(Point *center, BOX *box);
 static bool box_ov(BOX *box1, BOX *box2);
-static double box_ar(BOX *box);
-static double box_ht(BOX *box);
-static double box_wd(BOX *box);
+static float8 box_ar(BOX *box);
+static float8 box_ht(BOX *box);
+static float8 box_wd(BOX *box);
 static bool box_contain_point(BOX *box, Point *point);
 static bool box_contain_box(BOX *box1, BOX *box2);
 static bool box_contain_lseg(BOX *box, LSEG *lseg);
 static bool box_interpt_lseg(Point *result, BOX *box, LSEG *lseg);
 static float8 box_closept_point(Point *result, BOX *box, Point *point);
 static float8 box_closept_lseg(Point *result, BOX *box, LSEG *lseg);
 
 /* Routines for circles */
-static double circle_ar(CIRCLE *circle);
+static float8 circle_ar(CIRCLE *circle);
 
 /* Routines for polygons */
 static void make_bound_box(POLYGON *poly);
 static void poly_to_circle(CIRCLE *result, POLYGON *poly);
 static bool lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start);
 static bool poly_contain_poly(POLYGON *polya, POLYGON *polyb);
 static bool plist_same(int npts, Point *p1, Point *p2);
 static float8 dist_ppoly_internal(Point *pt, POLYGON *poly);
 
 /* Routines for encoding and decoding */
-static double single_decode(char *num, char **endptr_p,
+static float8 single_decode(char *num, char **endptr_p,
 			  const char *type_name, const char *orig_string);
 static void single_encode(float8 x, StringInfo str);
-static void pair_decode(char *str, double *x, double *y, char **endptr_p,
+static void pair_decode(char *str, float8 *x, float8 *y, char **endptr_p,
 			const char *type_name, const char *orig_string);
 static void pair_encode(float8 x, float8 y, StringInfo str);
 static int	pair_count(char *s, char delim);
 static void path_decode(char *str, bool opentype, int npts, Point *p,
 			bool *isopen, char **endptr_p,
 			const char *type_name, const char *orig_string);
 static char *path_encode(enum path_delim path_delim, int npts, Point *pt);
 
 
 /*
@@ -138,38 +138,38 @@ static char *path_encode(enum path_delim path_delim, int npts, Point *pt);
  *
  * For boxes, the points are opposite corners with the first point at the top right.
  * For closed paths and polygons, the points should be reordered to allow
  *	fast and correct equality comparisons.
  *
  * XXX perhaps points in complex shapes should be reordered internally
  *	to allow faster internal operations, but should keep track of input order
  *	and restore that order for text output - tgl 97/01/16
  */
 
-static double
+static float8
 single_decode(char *num, char **endptr_p,
 			  const char *type_name, const char *orig_string)
 {
 	return float8in_internal(num, endptr_p, type_name, orig_string);
 }								/* single_decode() */
 
 static void
 single_encode(float8 x, StringInfo str)
 {
 	char	   *xstr = float8out_internal(x);
 
 	appendStringInfoString(str, xstr);
 	pfree(xstr);
 }								/* single_encode() */
 
 static void
-pair_decode(char *str, double *x, double *y, char **endptr_p,
+pair_decode(char *str, float8 *x, float8 *y, char **endptr_p,
 			const char *type_name, const char *orig_string)
 {
 	bool		has_delim;
 
 	while (isspace((unsigned char) *str))
 		str++;
 	if ((has_delim = (*str == LDELIM)))
 		str++;
 
 	*x = float8in_internal(str, &str, type_name, orig_string);
@@ -368,33 +368,33 @@ pair_count(char *s, char delim)
  *		External format: (two corners of box)
  *				"(f8, f8), (f8, f8)"
  *				also supports the older style "(f8, f8, f8, f8)"
  */
 Datum
 box_in(PG_FUNCTION_ARGS)
 {
 	char	   *str = PG_GETARG_CSTRING(0);
 	BOX		   *box = (BOX *) palloc(sizeof(BOX));
 	bool		isopen;
-	double		x,
+	float8		x,
 				y;
 
 	path_decode(str, false, 2, &(box->high), &isopen, NULL, "box", str);
 
 	/* reorder corners if necessary... */
-	if (box->high.x < box->low.x)
+	if (float8_lt(box->high.x, box->low.x))
 	{
 		x = box->high.x;
 		box->high.x = box->low.x;
 		box->low.x = x;
 	}
-	if (box->high.y < box->low.y)
+	if (float8_lt(box->high.y, box->low.y))
 	{
 		y = box->high.y;
 		box->high.y = box->low.y;
 		box->low.y = y;
 	}
 
 	PG_RETURN_BOX_P(box);
 }
 
 /*		box_out -		convert a box to external form.
@@ -408,38 +408,38 @@ box_out(PG_FUNCTION_ARGS)
 }
 
 /*
  *		box_recv			- converts external binary format to box
  */
 Datum
 box_recv(PG_FUNCTION_ARGS)
 {
 	StringInfo	buf = (StringInfo) PG_GETARG_POINTER(0);
 	BOX		   *box;
-	double		x,
+	float8		x,
 				y;
 
 	box = (BOX *) palloc(sizeof(BOX));
 
 	box->high.x = pq_getmsgfloat8(buf);
 	box->high.y = pq_getmsgfloat8(buf);
 	box->low.x = pq_getmsgfloat8(buf);
 	box->low.y = pq_getmsgfloat8(buf);
 
 	/* reorder corners if necessary... */
-	if (box->high.x < box->low.x)
+	if (float8_lt(box->high.x, box->low.x))
 	{
 		x = box->high.x;
 		box->high.x = box->low.x;
 		box->low.x = x;
 	}
-	if (box->high.y < box->low.y)
+	if (float8_lt(box->high.y, box->low.y))
 	{
 		y = box->high.y;
 		box->high.y = box->low.y;
 		box->low.y = y;
 	}
 
 	PG_RETURN_BOX_P(box);
 }
 
 /*
@@ -458,31 +458,31 @@ box_send(PG_FUNCTION_ARGS)
 	pq_sendfloat8(&buf, box->low.y);
 	PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
 }
 
 
 /*		box_construct	-		fill in a new box.
  */
 static inline void
 box_construct(BOX *result, Point *pt1, Point *pt2)
 {
-	if (pt1->x > pt2->x)
+	if (float8_gt(pt1->x, pt2->x))
 	{
 		result->high.x = pt1->x;
 		result->low.x = pt2->x;
 	}
 	else
 	{
 		result->high.x = pt2->x;
 		result->low.x = pt1->x;
 	}
-	if (pt1->y > pt2->y)
+	if (float8_gt(pt1->y, pt2->y))
 	{
 		result->high.y = pt1->y;
 		result->low.y = pt2->y;
 	}
 	else
 	{
 		result->high.y = pt2->y;
 		result->low.y = pt1->y;
 	}
 }
@@ -800,54 +800,54 @@ box_center(PG_FUNCTION_ARGS)
 	Point	   *result = (Point *) palloc(sizeof(Point));
 
 	box_cn(result, box);
 
 	PG_RETURN_POINT_P(result);
 }
 
 
 /*		box_ar	-		returns the area of the box.
  */
-static double
+static float8
 box_ar(BOX *box)
 {
-	return box_wd(box) * box_ht(box);
+	return float8_mul(box_wd(box), box_ht(box));
 }
 
 
 /*		box_cn	-		stores the centerpoint of the box into *center.
  */
 static void
 box_cn(Point *center, BOX *box)
 {
-	center->x = (box->high.x + box->low.x) / 2.0;
-	center->y = (box->high.y + box->low.y) / 2.0;
+	center->x = float8_div(float8_pl(box->high.x, box->low.x), 2.0);
+	center->y = float8_div(float8_pl(box->high.y, box->low.y), 2.0);
 }
 
 
 /*		box_wd	-		returns the width (length) of the box
  *								  (horizontal magnitude).
  */
-static double
+static float8
 box_wd(BOX *box)
 {
-	return box->high.x - box->low.x;
+	return float8_mi(box->high.x, box->low.x);
 }
 
 
 /*		box_ht	-		returns the height of the box
  *								  (vertical magnitude).
  */
-static double
+static float8
 box_ht(BOX *box)
 {
-	return box->high.y - box->low.y;
+	return float8_mi(box->high.y, box->low.y);
 }
 
 
 /*----------------------------------------------------------
  *	Funky operations.
  *---------------------------------------------------------*/
 
 /*		box_intersect	-
  *				returns the overlapping portion of two boxes,
  *				  or NULL if they do not intersect.
@@ -857,24 +857,24 @@ box_intersect(PG_FUNCTION_ARGS)
 {
 	BOX		   *box1 = PG_GETARG_BOX_P(0);
 	BOX		   *box2 = PG_GETARG_BOX_P(1);
 	BOX		   *result;
 
 	if (!box_ov(box1, box2))
 		PG_RETURN_NULL();
 
 	result = (BOX *) palloc(sizeof(BOX));
 
-	result->high.x = Min(box1->high.x, box2->high.x);
-	result->low.x = Max(box1->low.x, box2->low.x);
-	result->high.y = Min(box1->high.y, box2->high.y);
-	result->low.y = Max(box1->low.y, box2->low.y);
+	result->high.x = float8_min(box1->high.x, box2->high.x);
+	result->low.x = float8_max(box1->low.x, box2->low.x);
+	result->high.y = float8_min(box1->high.y, box2->high.y);
+	result->low.y = float8_max(box1->low.y, box2->low.y);
 
 	PG_RETURN_BOX_P(result);
 }
 
 
 /*		box_diagonal	-
  *				returns a line segment which happens to be the
  *				  positive-slope diagonal of "box".
  */
 Datum
@@ -1006,30 +1006,30 @@ line_send(PG_FUNCTION_ARGS)
 
 /*
  * Fill already-allocated LINE struct from the point and the slope
  */
 static inline void
 line_construct(LINE *result, Point *pt, float8 m)
 {
 	if (m == DBL_MAX)
 	{
 		/* vertical - use "x = C" */
-		result->A = -1;
-		result->B = 0;
+		result->A = -1.0;
+		result->B = 0.0;
 		result->C = pt->x;
 	}
 	else
 	{
 		/* use "mx - y + yinter = 0" */
 		result->A = m;
 		result->B = -1.0;
-		result->C = pt->y - m * pt->x;
+		result->C = float8_mi(pt->y, float8_mul(m, pt->x));
 	}
 }
 
 /* line_construct_pp()
  * two points
  */
 Datum
 line_construct_pp(PG_FUNCTION_ARGS)
 {
 	Point	   *pt1 = PG_GETARG_POINT_P(0);
@@ -1068,94 +1068,105 @@ Datum
 line_perp(PG_FUNCTION_ARGS)
 {
 	LINE	   *l1 = PG_GETARG_LINE_P(0);
 	LINE	   *l2 = PG_GETARG_LINE_P(1);
 
 	if (FPzero(l1->A))
 		PG_RETURN_BOOL(FPzero(l2->B));
 	else if (FPzero(l1->B))
 		PG_RETURN_BOOL(FPzero(l2->A));
 
-	PG_RETURN_BOOL(FPeq(((l1->A * l2->B) / (l1->B * l2->A)), -1.0));
+	PG_RETURN_BOOL(FPeq(float8_div(float8_mul(l1->A, l2->B),
+								   float8_mul(l1->B, l2->A)), -1.0));
 }
 
 Datum
 line_vertical(PG_FUNCTION_ARGS)
 {
 	LINE	   *line = PG_GETARG_LINE_P(0);
 
 	PG_RETURN_BOOL(FPzero(line->B));
 }
 
 Datum
 line_horizontal(PG_FUNCTION_ARGS)
 {
 	LINE	   *line = PG_GETARG_LINE_P(0);
 
 	PG_RETURN_BOOL(FPzero(line->A));
 }
 
+
+/*
+ * Check whether the two lines are the same
+ *
+ * We consider NaNs values to be equal to each other to let those lines
+ * to be found.
+ */
 Datum
 line_eq(PG_FUNCTION_ARGS)
 {
 	LINE	   *l1 = PG_GETARG_LINE_P(0);
 	LINE	   *l2 = PG_GETARG_LINE_P(1);
-	double		k;
+	float8		ratio;
 
-	if (!FPzero(l2->A))
-		k = l1->A / l2->A;
-	else if (!FPzero(l2->B))
-		k = l1->B / l2->B;
-	else if (!FPzero(l2->C))
-		k = l1->C / l2->C;
+	if (!FPzero(l2->A) && !isnan(l2->A))
+		ratio = float8_div(l1->A, l2->A);
+	else if (!FPzero(l2->B) && !isnan(l2->B))
+		ratio = float8_div(l1->B, l2->B);
+	else if (!FPzero(l2->C) && !isnan(l2->C))
+		ratio = float8_div(l1->C, l2->C);
 	else
-		k = 1.0;
+		ratio = 1.0;
 
-	PG_RETURN_BOOL(FPeq(l1->A, k * l2->A) &&
-				   FPeq(l1->B, k * l2->B) &&
-				   FPeq(l1->C, k * l2->C));
+	PG_RETURN_BOOL((FPeq(l1->A, float8_mul(ratio, l2->A)) &&
+					FPeq(l1->B, float8_mul(ratio, l2->B)) &&
+					FPeq(l1->C, float8_mul(ratio, l2->C))) ||
+				   (float8_eq(l1->A, l2->A) &&
+					float8_eq(l1->B, l2->B) &&
+					float8_eq(l1->C, l2->C)));
 }
 
 
 /*----------------------------------------------------------
  *	Line arithmetic routines.
  *---------------------------------------------------------*/
 
 /*
  * Return inverse slope of the line
  */
 static inline float8
 line_invsl(LINE *line)
 {
 	if (FPzero(line->A))
 		return DBL_MAX;
 	if (FPzero(line->B))
 		return 0.0;
-	return line->B / line->A;
+	return float8_div(line->B, line->A);
 }
 
 
 /* line_distance()
  * Distance between two lines.
  */
 Datum
 line_distance(PG_FUNCTION_ARGS)
 {
 	LINE	   *l1 = PG_GETARG_LINE_P(0);
 	LINE	   *l2 = PG_GETARG_LINE_P(1);
 	float8		result;
 	Point		tmp;
 
 	if (line_interpt_line(NULL, l1, l2))
 		PG_RETURN_FLOAT8(0.0);
 	if (FPzero(l1->B))			/* vertical? */
-		PG_RETURN_FLOAT8(fabs(l1->C - l2->C));
+		PG_RETURN_FLOAT8(fabs(float8_mi(l1->C, l2->C)));
 	point_construct(&tmp, 0.0, l1->C);
 	result = line_closept_point(NULL, l2, &tmp);
 	PG_RETURN_FLOAT8(result);
 }
 
 /* line_interpt()
  * Point where two lines l1, l2 intersect (if any)
  */
 Datum
 line_interpt(PG_FUNCTION_ARGS)
@@ -1174,47 +1185,51 @@ line_interpt(PG_FUNCTION_ARGS)
 /*
  * Internal version of line_interpt
  *
  * This returns true if two lines intersect (they do, if they are not
  * parallel), false if they do not.  This also sets the intersection point
  * to *result, if it is not NULL.
  *
  * NOTE If the lines are identical then we will find they are parallel
  * and report "no intersection".  This is a little weird, but since
  * there's no *unique* intersection, maybe it's appropriate behavior.
+ *
+ * If the lines have NaN constants, we will return true, and the intersection
+ * point would have NaN coordinates.  We shouldn't return false in this case
+ * because that would mean the lines are parallel.
  */
 static bool
 line_interpt_line(Point *result, LINE *l1, LINE *l2)
 {
-	double		x,
+	float8		x,
 				y;
 
 	if (FPzero(l1->B))			/* l1 vertical? */
 	{
 		if (FPzero(l2->B))		/* l2 vertical? */
 			return false;
 
 		x = l1->C;
-		y = (l2->A * x + l2->C);
+		y = float8_pl(float8_mul(l2->A, x), l2->C);
 	}
 	else if (FPzero(l2->B))		/* l2 vertical? */
 	{
 		x = l2->C;
-		y = (l1->A * x + l1->C);
+		y = float8_pl(float8_mul(l1->A, x), l1->C);
 	}
 	else
 	{
-		if (FPeq(l2->A, l1->A * (l2->B / l1->B)))
+		if (FPeq(l2->A, float8_mul(l1->A, float8_div(l2->B, l1->B))))
 			return false;
 
-		x = (l1->C - l2->C) / (l2->A - l1->A);
-		y = (l1->A * x + l1->C);
+		x = float8_div(float8_mi(l1->C, l2->C), float8_mi(l2->A, l1->A));
+		y = float8_pl(float8_mul(l1->A, x), l1->C);
 	}
 	if (result != NULL)
 		point_construct(result, x, y);
 
 	return true;
 }
 
 
 /***********************************************************************
  **
@@ -1235,36 +1250,35 @@ line_interpt_line(Point *result, LINE *l1, LINE *l2)
  *				"(xcoord, ycoord),... "
  *				"[xcoord, ycoord,... ]"
  *		Also support older format:
  *				"(closed, npts, xcoord, ycoord,... )"
  *---------------------------------------------------------*/
 
 Datum
 path_area(PG_FUNCTION_ARGS)
 {
 	PATH	   *path = PG_GETARG_PATH_P(0);
-	double		area = 0.0;
+	float8		area = 0.0;
 	int			i,
 				j;
 
 	if (!path->closed)
 		PG_RETURN_NULL();
 
 	for (i = 0; i < path->npts; i++)
 	{
 		j = (i + 1) % path->npts;
-		area += path->p[i].x * path->p[j].y;
-		area -= path->p[i].y * path->p[j].x;
+		area = float8_pl(area, float8_mul(path->p[i].x, path->p[j].y));
+		area = float8_mi(area, float8_mul(path->p[i].y, path->p[j].x));
 	}
 
-	area *= 0.5;
-	PG_RETURN_FLOAT8(area < 0.0 ? -area : area);
+	PG_RETURN_FLOAT8(float8_div(fabs(area), 2.0));
 }
 
 
 Datum
 path_in(PG_FUNCTION_ARGS)
 {
 	char	   *str = PG_GETARG_CSTRING(0);
 	PATH	   *path;
 	bool		isopen;
 	char	   *s;
@@ -1520,33 +1534,33 @@ path_inter(PG_FUNCTION_ARGS)
 				j;
 	LSEG		seg1,
 				seg2;
 
 	Assert(p1->npts > 0 && p2->npts > 0);
 
 	b1.high.x = b1.low.x = p1->p[0].x;
 	b1.high.y = b1.low.y = p1->p[0].y;
 	for (i = 1; i < p1->npts; i++)
 	{
-		b1.high.x = Max(p1->p[i].x, b1.high.x);
-		b1.high.y = Max(p1->p[i].y, b1.high.y);
-		b1.low.x = Min(p1->p[i].x, b1.low.x);
-		b1.low.y = Min(p1->p[i].y, b1.low.y);
+		b1.high.x = float8_max(p1->p[i].x, b1.high.x);
+		b1.high.y = float8_max(p1->p[i].y, b1.high.y);
+		b1.low.x = float8_min(p1->p[i].x, b1.low.x);
+		b1.low.y = float8_min(p1->p[i].y, b1.low.y);
 	}
 	b2.high.x = b2.low.x = p2->p[0].x;
 	b2.high.y = b2.low.y = p2->p[0].y;
 	for (i = 1; i < p2->npts; i++)
 	{
-		b2.high.x = Max(p2->p[i].x, b2.high.x);
-		b2.high.y = Max(p2->p[i].y, b2.high.y);
-		b2.low.x = Min(p2->p[i].x, b2.low.x);
-		b2.low.y = Min(p2->p[i].y, b2.low.y);
+		b2.high.x = float8_max(p2->p[i].x, b2.high.x);
+		b2.high.y = float8_max(p2->p[i].y, b2.high.y);
+		b2.low.x = float8_min(p2->p[i].x, b2.low.x);
+		b2.low.y = float8_min(p2->p[i].y, b2.low.y);
 	}
 	if (!box_ov(&b1, &b2))
 		PG_RETURN_BOOL(false);
 
 	/* pairwise check lseg intersections */
 	for (i = 0; i < p1->npts; i++)
 	{
 		int			iprev;
 
 		if (i > 0)
@@ -1622,21 +1636,21 @@ path_distance(PG_FUNCTION_ARGS)
 			{
 				if (!p2->closed)
 					continue;
 				jprev = p2->npts - 1;	/* include the closure segment */
 			}
 
 			statlseg_construct(&seg1, &p1->p[iprev], &p1->p[i]);
 			statlseg_construct(&seg2, &p2->p[jprev], &p2->p[j]);
 
 			tmp = lseg_closept_lseg(NULL, &seg1, &seg2);
-			if (!have_min || tmp < min)
+			if (!have_min || float8_lt(tmp, min))
 			{
 				min = tmp;
 				have_min = true;
 			}
 		}
 	}
 
 	if (!have_min)
 		PG_RETURN_NULL();
 
@@ -1661,21 +1675,21 @@ path_length(PG_FUNCTION_ARGS)
 
 		if (i > 0)
 			iprev = i - 1;
 		else
 		{
 			if (!path->closed)
 				continue;
 			iprev = path->npts - 1; /* include the closure segment */
 		}
 
-		result += point_dt(&path->p[iprev], &path->p[i]);
+		result = float8_pl(result, point_dt(&path->p[iprev], &path->p[i]));
 	}
 
 	PG_RETURN_FLOAT8(result);
 }
 
 /***********************************************************************
  **
  **		Routines for 2D points.
  **
  ***********************************************************************/
@@ -1824,44 +1838,52 @@ point_eq(PG_FUNCTION_ARGS)
 
 Datum
 point_ne(PG_FUNCTION_ARGS)
 {
 	Point	   *pt1 = PG_GETARG_POINT_P(0);
 	Point	   *pt2 = PG_GETARG_POINT_P(1);
 
 	PG_RETURN_BOOL(!point_eq_point(pt1, pt2));
 }
 
+
+/*
+ * Check whether the two points are the same
+ *
+ * We consider NaNs coordinates to be equal to each other to let those points
+ * to be found.
+ */
 static inline bool
 point_eq_point(Point *pt1, Point *pt2)
 {
-	return FPeq(pt1->x, pt2->x) && FPeq(pt1->y, pt2->y);
+	return ((FPeq(pt1->x, pt2->x) && FPeq(pt1->y, pt2->y)) ||
+			(float8_eq(pt1->x, pt2->x) && float8_eq(pt1->y, pt2->y)));
 }
 
 
 /*----------------------------------------------------------
  *	"Arithmetic" operators on points.
  *---------------------------------------------------------*/
 
 Datum
 point_distance(PG_FUNCTION_ARGS)
 {
 	Point	   *pt1 = PG_GETARG_POINT_P(0);
 	Point	   *pt2 = PG_GETARG_POINT_P(1);
 
 	PG_RETURN_FLOAT8(point_dt(pt1, pt2));
 }
 
 static inline float8
 point_dt(Point *pt1, Point *pt2)
 {
-	return HYPOT(pt1->x - pt2->x, pt1->y - pt2->y);
+	return HYPOT(float8_mi(pt1->x, pt2->x), float8_mi(pt1->y, pt2->y));
 }
 
 Datum
 point_slope(PG_FUNCTION_ARGS)
 {
 	Point	   *pt1 = PG_GETARG_POINT_P(0);
 	Point	   *pt2 = PG_GETARG_POINT_P(1);
 
 	PG_RETURN_FLOAT8(point_sl(pt1, pt2));
 }
@@ -1872,37 +1894,37 @@ point_slope(PG_FUNCTION_ARGS)
  *
  * Note that this function returns DBL_MAX when the points are the same.
  */
 static inline float8
 point_sl(Point *pt1, Point *pt2)
 {
 	if (FPeq(pt1->x, pt2->x))
 		return DBL_MAX;
 	if (FPeq(pt1->y, pt2->y))
 		return 0.0;
-	return (pt1->y - pt2->y) / (pt1->x - pt2->x);
+	return float8_div(float8_mi(pt1->y, pt2->y), float8_mi(pt1->x, pt2->x));
 }
 
 
 /*
  * Return inverse slope of two points
  *
  * Note that this function returns 0.0 when the points are the same.
  */
 static inline float8
 point_invsl(Point *pt1, Point *pt2)
 {
 	if (FPeq(pt1->x, pt2->x))
 		return 0.0;
 	if (FPeq(pt1->y, pt2->y))
 		return DBL_MAX;
-	return (pt1->x - pt2->x) / (pt2->y - pt1->y);
+	return float8_div(float8_mi(pt1->x, pt2->x), float8_mi(pt2->y, pt1->y));
 }
 
 
 /***********************************************************************
  **
  **		Routines for 2D line segments.
  **
  ***********************************************************************/
 
 /*----------------------------------------------------------
@@ -2162,22 +2184,22 @@ lseg_distance(PG_FUNCTION_ARGS)
 
 
 Datum
 lseg_center(PG_FUNCTION_ARGS)
 {
 	LSEG	   *lseg = PG_GETARG_LSEG_P(0);
 	Point	   *result;
 
 	result = (Point *) palloc(sizeof(Point));
 
-	result->x = (lseg->p[0].x + lseg->p[1].x) / 2.0;
-	result->y = (lseg->p[0].y + lseg->p[1].y) / 2.0;
+	result->x = float8_div(float8_pl(lseg->p[0].x, lseg->p[1].x), 2.0);
+	result->y = float8_div(float8_pl(lseg->p[0].y, lseg->p[1].y), 2.0);
 
 	PG_RETURN_POINT_P(result);
 }
 
 
 /*
  *		Find the intersection point of two segments (if any).
  *
  * This returns true if two line segments intersect, false if they do not.
  * This also sets the intersection point to *result, if it is not NULL.
@@ -2286,21 +2308,21 @@ dist_ppath(PG_FUNCTION_ARGS)
 			iprev = i - 1;
 		else
 		{
 			if (!path->closed)
 				continue;
 			iprev = path->npts - 1; /* Include the closure segment */
 		}
 
 		statlseg_construct(&lseg, &path->p[iprev], &path->p[i]);
 		tmp = lseg_closept_point(NULL, &lseg, pt);
-		if (!have_min || tmp < result)
+		if (!have_min || float8_lt(tmp, result))
 		{
 			result = tmp;
 			have_min = true;
 		}
 	}
 
 	PG_RETURN_FLOAT8(result);
 }
 
 /*
@@ -2316,33 +2338,28 @@ dist_pb(PG_FUNCTION_ARGS)
 }
 
 /*
  * Distance from a lseg to a line
  */
 Datum
 dist_sl(PG_FUNCTION_ARGS)
 {
 	LSEG	   *lseg = PG_GETARG_LSEG_P(0);
 	LINE	   *line = PG_GETARG_LINE_P(1);
-	float8		result,
-				d2;
+	float8		result;
 
 	if (lseg_interpt_line(NULL, lseg, line))
 		result = 0.0;
 	else
-	{
-		result = line_closept_point(NULL, line, &lseg->p[0]);
-		d2 = line_closept_point(NULL, line, &lseg->p[1]);
 		/* XXX shouldn't we take the min not max? */
-		if (d2 > result)
-			result = d2;
-	}
+		result = float8_max(line_closept_point(NULL, line, &lseg->p[0]),
+							line_closept_point(NULL, line, &lseg->p[1]));
 
 	PG_RETURN_FLOAT8(result);
 }
 
 /*
  * Distance from a lseg to a box
  */
 Datum
 dist_sb(PG_FUNCTION_ARGS)
 {
@@ -2375,25 +2392,24 @@ dist_lb(PG_FUNCTION_ARGS)
  * Distance from a circle to a polygon
  */
 Datum
 dist_cpoly(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle = PG_GETARG_CIRCLE_P(0);
 	POLYGON    *poly = PG_GETARG_POLYGON_P(1);
 	float8		result;
 
 	/* calculate distance to center, and subtract radius */
-	result = dist_ppoly_internal(&circle->center, poly);
-
-	result -= circle->radius;
-	if (result < 0)
-		result = 0;
+	result = float8_mi(dist_ppoly_internal(&circle->center, poly),
+					   circle->radius);
+	if (result < 0.0)
+		result = 0.0;
 
 	PG_RETURN_FLOAT8(result);
 }
 
 /*
  * Distance from a point to a polygon
  */
 Datum
 dist_ppoly(PG_FUNCTION_ARGS)
 {
@@ -2442,21 +2458,21 @@ dist_ppoly_internal(Point *pt, POLYGON *poly)
 	for (i = 0; i < poly->npts - 1; i++)
 	{
 		seg.p[0].x = poly->p[i].x;
 		seg.p[0].y = poly->p[i].y;
 		seg.p[1].x = poly->p[i + 1].x;
 		seg.p[1].y = poly->p[i + 1].y;
 		d = lseg_closept_point(NULL, &seg, pt);
 #ifdef GEODEBUG
 		printf("dist_ppoly_internal- segment %d distance is %f\n", (i + 1), d);
 #endif
-		if (d < result)
+		if (float8_lt(d, result))
 			result = d;
 	}
 
 	return result;
 }
 
 
 /*---------------------------------------------------------------------
  *		interpt_
  *				Intersection point of objects.
@@ -2535,21 +2551,22 @@ line_closept_point(Point *result, LINE *line, Point *point)
 
 Datum
 close_pl(PG_FUNCTION_ARGS)
 {
 	Point	   *pt = PG_GETARG_POINT_P(0);
 	LINE	   *line = PG_GETARG_LINE_P(1);
 	Point	   *result;
 
 	result = (Point *) palloc(sizeof(Point));
 
-	line_closept_point(result, line, pt);
+	if (isnan(line_closept_point(result, line, pt)))
+		PG_RETURN_NULL();
 
 	PG_RETURN_POINT_P(result);
 }
 
 
 /*
  * Closest point on line segment to specified point.
  *
  * This sets the closest point to the *result if it is not NULL and returns
  * the distance to the closest point.
@@ -2595,134 +2612,136 @@ close_ps(PG_FUNCTION_ARGS)
  * This sets the closest point to the *result if it is not NULL and returns
  * the distance to the closest point.
  *
  * XXX This function is wrong.  If must never set the *result to a point on
  * the second segment.
  */
 static float8
 lseg_closept_lseg(Point *result, LSEG *l1, LSEG *l2)
 {
 	Point		point;
-	double		dist;
-	double		d;
+	float8		dist,
+				d;
 
 	d = lseg_closept_point(NULL, l1, &l2->p[0]);
 	dist = d;
 	if (result != NULL)
 		*result = l2->p[0];
 
 	d = lseg_closept_point(NULL, l1, &l2->p[1]);
-	if (d < dist)
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		if (result != NULL)
 			*result = l2->p[1];
 	}
 
-	if (lseg_closept_point(&point, l2, &l1->p[0]) < dist)
+	if (float8_lt(lseg_closept_point(&point, l2, &l1->p[0]), dist))
 		d = lseg_closept_point(result, l1, &point);
 
-	if (lseg_closept_point(&point, l2, &l1->p[1]) < dist)
+	if (float8_lt(lseg_closept_point(&point, l2, &l1->p[1]), dist))
 		d = lseg_closept_point(result, l1, &point);
 
-	if (d < dist)
+	if (float8_lt(d, dist))
 		dist = d;
 
 	return dist;
 }
 
 Datum
 close_lseg(PG_FUNCTION_ARGS)
 {
 	LSEG	   *l1 = PG_GETARG_LSEG_P(0);
 	LSEG	   *l2 = PG_GETARG_LSEG_P(1);
 	Point	   *result;
 
 	result = (Point *) palloc(sizeof(Point));
 
-	lseg_closept_lseg(result, l2, l1);
+	if (isnan(lseg_closept_lseg(result, l2, l1)))
+		PG_RETURN_NULL();
 
 	PG_RETURN_POINT_P(result);
 }
 
 
 /*
  * Closest point on or in box to specified point.
  *
  * This sets the closest point to the *result if it is not NULL and returns
  * the distance to the closest point.
  */
 static float8
 box_closept_point(Point *result, BOX *box, Point *pt)
 {
-	LSEG		lseg;
+	float8		dist,
+				d;
 	Point		point,
 				closept;
-	double		dist,
-				d;
+	LSEG		lseg;
 
 	if (box_contain_point(box, pt))
 	{
 		if (result != NULL)
 			*result = *pt;
 
 		return 0.0;
 	}
 
 	/* pairwise check lseg distances */
 	point.x = box->low.x;
 	point.y = box->high.y;
 	statlseg_construct(&lseg, &box->low, &point);
 	dist = lseg_closept_point(result, &lseg, pt);
 
 	statlseg_construct(&lseg, &box->high, &point);
 	d = lseg_closept_point(&closept, &lseg, pt);
-	if (d < dist)
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		if (result != NULL)
 			*result = closept;
 	}
 
 	point.x = box->high.x;
 	point.y = box->low.y;
 	statlseg_construct(&lseg, &box->low, &point);
 	d = lseg_closept_point(&closept, &lseg, pt);
-	if (d < dist)
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		if (result != NULL)
 			*result = closept;
 	}
 
 	statlseg_construct(&lseg, &box->high, &point);
 	d = lseg_closept_point(&closept, &lseg, pt);
-	if (d < dist)
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		if (result != NULL)
 			*result = closept;
 	}
 
 	return dist;
 }
 
 Datum
 close_pb(PG_FUNCTION_ARGS)
 {
 	Point	   *pt = PG_GETARG_POINT_P(0);
 	BOX		   *box = PG_GETARG_BOX_P(1);
 	Point	   *result;
 
 	result = (Point *) palloc(sizeof(Point));
 
-	box_closept_point(result, box, pt);
+	if (isnan(box_closept_point(result, box, pt)))
+		PG_RETURN_NULL();
 
 	PG_RETURN_POINT_P(result);
 }
 
 
 /* close_sl()
  * Closest point on line to line segment.
  *
  * XXX THIS CODE IS WRONG
  * The code is actually calculating the point on the line segment
@@ -2740,21 +2759,21 @@ close_sl(PG_FUNCTION_ARGS)
 	float8		d1,
 				d2;
 
 	result = (Point *) palloc(sizeof(Point));
 
 	if (lseg_interpt_line(result, lseg, line))
 		PG_RETURN_POINT_P(result);
 
 	d1 = line_closept_point(NULL, line, &lseg->p[0]);
 	d2 = line_closept_point(NULL, line, &lseg->p[1]);
-	if (d1 < d2)
+	if (float8_lt(d1, d2))
 		*result = lseg->p[0];
 	else
 		*result = lseg->p[1];
 
 	PG_RETURN_POINT_P(result);
 #endif
 
 	ereport(ERROR,
 			(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
 			 errmsg("function \"close_sl\" not implemented")));
@@ -2804,92 +2823,94 @@ lseg_closept_line(Point *result, LSEG *lseg, LINE *line)
 
 Datum
 close_ls(PG_FUNCTION_ARGS)
 {
 	LINE	   *line = PG_GETARG_LINE_P(0);
 	LSEG	   *lseg = PG_GETARG_LSEG_P(1);
 	Point	   *result;
 
 	result = (Point *) palloc(sizeof(Point));
 
-	lseg_closept_line(result, lseg, line);
+	if (isnan(lseg_closept_line(result, lseg, line)))
+		PG_RETURN_NULL();
 
 	PG_RETURN_POINT_P(result);
 }
 
 
 /*
  * Closest point on or in box to line segment.
  *
  * This sets the closest point to the *result if it is not NULL and returns
  * the distance to the closest point.
  */
 static float8
 box_closept_lseg(Point *result, BOX *box, LSEG *lseg)
 {
+	float8		dist,
+				d;
 	Point		point,
 				closept;
 	LSEG		bseg;
-	double		dist,
-				d;
 
 	if (box_interpt_lseg(result, box, lseg))
 		return 0.0;
 
 	/* pairwise check lseg distances */
 	point.x = box->low.x;
 	point.y = box->high.y;
 	statlseg_construct(&bseg, &box->low, &point);
 	dist = lseg_closept_lseg(result, &bseg, lseg);
 
 	statlseg_construct(&bseg, &box->high, &point);
 	d = lseg_closept_lseg(&closept, &bseg, lseg);
-	if (d < dist)
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		if (result != NULL)
 			*result = closept;
 	}
 
 	point.x = box->high.x;
 	point.y = box->low.y;
 	statlseg_construct(&bseg, &box->low, &point);
 	d = lseg_closept_lseg(&closept, &bseg, lseg);
-	if (d < dist)
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		if (result != NULL)
 			*result = closept;
 	}
 
 	statlseg_construct(&bseg, &box->high, &point);
 	d = lseg_closept_lseg(&closept, &bseg, lseg);
-	if (d < dist)
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		if (result != NULL)
 			*result = closept;
 	}
 
 	return dist;
 }
 
 Datum
 close_sb(PG_FUNCTION_ARGS)
 {
 	LSEG	   *lseg = PG_GETARG_LSEG_P(0);
 	BOX		   *box = PG_GETARG_BOX_P(1);
 	Point	   *result;
 
 	result = (Point *) palloc(sizeof(Point));
 
-	box_closept_lseg(result, box, lseg);
+	if (isnan(box_closept_lseg(result, box, lseg)))
+		PG_RETURN_NULL();
 
 	PG_RETURN_POINT_P(result);
 }
 
 
 Datum
 close_lb(PG_FUNCTION_ARGS)
 {
 #ifdef NOT_USED
 	LINE	   *line = PG_GETARG_LINE_P(0);
@@ -2908,21 +2929,23 @@ close_lb(PG_FUNCTION_ARGS)
  *		on_
  *				Whether one object lies completely within another.
  *-------------------------------------------------------------------*/
 
 /*
  *		Does the point satisfy the equation?
  */
 static bool
 line_contain_point(LINE *line, Point *point)
 {
-	return FPzero(line->A * point->x + line->B * point->y + line->C);
+	return FPzero(float8_pl(float8_pl(float8_mul(line->A, point->x),
+									  float8_mul(line->B, point->y)),
+							line->C));
 }
 
 Datum
 on_pl(PG_FUNCTION_ARGS)
 {
 	Point	   *pt = PG_GETARG_POINT_P(0);
 	LINE	   *line = PG_GETARG_LINE_P(1);
 
 	PG_RETURN_BOOL(line_contain_point(line, pt));
 }
@@ -2989,33 +3012,32 @@ box_contain_pt(PG_FUNCTION_ARGS)
  *				but not cross.
  *				(we can do p-in-p in lg(n), but it takes preprocessing)
  */
 Datum
 on_ppath(PG_FUNCTION_ARGS)
 {
 	Point	   *pt = PG_GETARG_POINT_P(0);
 	PATH	   *path = PG_GETARG_PATH_P(1);
 	int			i,
 				n;
-	double		a,
+	float8		a,
 				b;
 
 	/*-- OPEN --*/
 	if (!path->closed)
 	{
 		n = path->npts - 1;
 		a = point_dt(pt, &path->p[0]);
 		for (i = 0; i < n; i++)
 		{
 			b = point_dt(pt, &path->p[i + 1]);
-			if (FPeq(a + b,
-					 point_dt(&path->p[i], &path->p[i + 1])))
+			if (FPeq(float8_pl(a, b), point_dt(&path->p[i], &path->p[i + 1])))
 				PG_RETURN_BOOL(true);
 			a = b;
 		}
 		PG_RETURN_BOOL(false);
 	}
 
 	/*-- CLOSED --*/
 	PG_RETURN_BOOL(point_inside(pt, path->npts, path->p) != 0);
 }
 
@@ -3087,24 +3109,24 @@ inter_sl(PG_FUNCTION_ARGS)
  * Optimize for non-intersection by checking for box intersection first.
  * - thomas 1998-01-30
  */
 static bool
 box_interpt_lseg(Point *result, BOX *box, LSEG *lseg)
 {
 	BOX			lbox;
 	LSEG		bseg;
 	Point		point;
 
-	lbox.low.x = Min(lseg->p[0].x, lseg->p[1].x);
-	lbox.low.y = Min(lseg->p[0].y, lseg->p[1].y);
-	lbox.high.x = Max(lseg->p[0].x, lseg->p[1].x);
-	lbox.high.y = Max(lseg->p[0].y, lseg->p[1].y);
+	lbox.low.x = float8_min(lseg->p[0].x, lseg->p[1].x);
+	lbox.low.y = float8_min(lseg->p[0].y, lseg->p[1].y);
+	lbox.high.x = float8_max(lseg->p[0].x, lseg->p[1].x);
+	lbox.high.y = float8_max(lseg->p[0].y, lseg->p[1].y);
 
 	/* nothing close to overlap? then not going to intersect */
 	if (!box_ov(&lbox, box))
 		return false;
 
 	if (result != NULL)
 	{
 		box_cn(&point, box);
 		lseg_closept_point(result, lseg, &point);
 	}
@@ -3197,38 +3219,38 @@ inter_lb(PG_FUNCTION_ARGS)
  * make_bound_box - create the bounding box for the input polygon
  *------------------------------------------------------------------*/
 
 /*---------------------------------------------------------------------
  * Make the smallest bounding box for the given polygon.
  *---------------------------------------------------------------------*/
 static void
 make_bound_box(POLYGON *poly)
 {
 	int			i;
-	double		x1,
+	float8		x1,
 				y1,
 				x2,
 				y2;
 
 	Assert(poly->npts > 0);
 
 	x1 = x2 = poly->p[0].x;
 	y2 = y1 = poly->p[0].y;
 	for (i = 1; i < poly->npts; i++)
 	{
-		if (poly->p[i].x < x1)
+		if (float8_lt(poly->p[i].x, x1))
 			x1 = poly->p[i].x;
-		if (poly->p[i].x > x2)
+		if (float8_gt(poly->p[i].x, x2))
 			x2 = poly->p[i].x;
-		if (poly->p[i].y < y1)
+		if (float8_lt(poly->p[i].y, y1))
 			y1 = poly->p[i].y;
-		if (poly->p[i].y > y2)
+		if (float8_gt(poly->p[i].y, y2))
 			y2 = poly->p[i].y;
 	}
 
 	poly->boundbox.low.x = x1;
 	poly->boundbox.high.x = x2;
 	poly->boundbox.low.y = y1;
 	poly->boundbox.high.y = y2;
 }
 
 /*------------------------------------------------------------------
@@ -3727,22 +3749,22 @@ lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start)
 	}
 
 	if (res && !intersection)
 	{
 		Point		p;
 
 		/*
 		 * if X-intersection wasn't found  then check central point of tested
 		 * segment. In opposite case we already check all subsegments
 		 */
-		p.x = (t.p[0].x + t.p[1].x) / 2.0;
-		p.y = (t.p[0].y + t.p[1].y) / 2.0;
+		p.x = float8_div(float8_pl(t.p[0].x, t.p[1].x), 2.0);
+		p.y = float8_div(float8_pl(t.p[0].y, t.p[1].y), 2.0);
 
 		res = point_inside(&p, poly->npts, poly->p);
 	}
 
 	return res;
 }
 
 /*-----------------------------------------------------------------
  * Determine if polygon A contains polygon B.
  *-----------------------------------------------------------------*/
@@ -3868,22 +3890,22 @@ construct_point(PG_FUNCTION_ARGS)
 	point_construct(result, x, y);
 
 	PG_RETURN_POINT_P(result);
 }
 
 
 static inline void
 point_add_point(Point *result, Point *pt1, Point *pt2)
 {
 	point_construct(result,
-					pt1->x + pt2->x,
-					pt1->y + pt2->y);
+					float8_pl(pt1->x, pt2->x),
+					float8_pl(pt1->y, pt2->y));
 }
 
 Datum
 point_add(PG_FUNCTION_ARGS)
 {
 	Point	   *p1 = PG_GETARG_POINT_P(0);
 	Point	   *p2 = PG_GETARG_POINT_P(1);
 	Point	   *result;
 
 	result = (Point *) palloc(sizeof(Point));
@@ -3891,22 +3913,22 @@ point_add(PG_FUNCTION_ARGS)
 	point_add_point(result, p1, p2);
 
 	PG_RETURN_POINT_P(result);
 }
 
 
 static inline void
 point_sub_point(Point *result, Point *pt1, Point *pt2)
 {
 	point_construct(result,
-					pt1->x - pt2->x,
-					pt1->y - pt2->y);
+					float8_mi(pt1->x, pt2->x),
+					float8_mi(pt1->y, pt2->y));
 }
 
 Datum
 point_sub(PG_FUNCTION_ARGS)
 {
 	Point	   *p1 = PG_GETARG_POINT_P(0);
 	Point	   *p2 = PG_GETARG_POINT_P(1);
 	Point	   *result;
 
 	result = (Point *) palloc(sizeof(Point));
@@ -3914,54 +3936,53 @@ point_sub(PG_FUNCTION_ARGS)
 	point_sub_point(result, p1, p2);
 
 	PG_RETURN_POINT_P(result);
 }
 
 
 static inline void
 point_mul_point(Point *result, Point *pt1, Point *pt2)
 {
 	point_construct(result,
-					(pt1->x * pt2->x) - (pt1->y * pt2->y),
-					(pt1->x * pt2->y) + (pt1->y * pt2->x));
+					float8_mi(float8_mul(pt1->x, pt2->x),
+							  float8_mul(pt1->y, pt2->y)),
+					float8_pl(float8_mul(pt1->x, pt2->y),
+							  float8_mul(pt1->y, pt2->x)));
 }
 
 Datum
 point_mul(PG_FUNCTION_ARGS)
 {
 	Point	   *p1 = PG_GETARG_POINT_P(0);
 	Point	   *p2 = PG_GETARG_POINT_P(1);
 	Point	   *result;
 
 	result = (Point *) palloc(sizeof(Point));
 
 	point_mul_point(result, p1, p2);
 
 	PG_RETURN_POINT_P(result);
 }
 
 
 static inline void
 point_div_point(Point *result, Point *pt1, Point *pt2)
 {
-	double		div;
+	float8		div;
 
-	div = (pt2->x * pt2->x) + (pt2->y * pt2->y);
-
-	if (div == 0.0)
-		ereport(ERROR,
-				(errcode(ERRCODE_DIVISION_BY_ZERO),
-				 errmsg("division by zero")));
+	div = float8_pl(float8_mul(pt2->x, pt2->x), float8_mul(pt2->y, pt2->y));
 
 	point_construct(result,
-					((pt1->x * pt2->x) + (pt1->y * pt2->y)) / div,
-					((pt2->x * pt1->y) - (pt2->y * pt1->x)) / div);
+					float8_div(float8_pl(float8_mul(pt1->x, pt2->x),
+										 float8_mul(pt1->y, pt2->y)), div),
+					float8_div(float8_mi(float8_mul(pt1->y, pt2->x),
+										 float8_mul(pt1->x, pt2->y)), div));
 }
 
 Datum
 point_div(PG_FUNCTION_ARGS)
 {
 	Point	   *p1 = PG_GETARG_POINT_P(0);
 	Point	   *p2 = PG_GETARG_POINT_P(1);
 	Point	   *result;
 
 	result = (Point *) palloc(sizeof(Point));
@@ -4084,24 +4105,24 @@ point_box(PG_FUNCTION_ARGS)
  */
 Datum
 boxes_bound_box(PG_FUNCTION_ARGS)
 {
 	BOX		   *box1 = PG_GETARG_BOX_P(0),
 			   *box2 = PG_GETARG_BOX_P(1),
 			   *container;
 
 	container = (BOX *) palloc(sizeof(BOX));
 
-	container->high.x = Max(box1->high.x, box2->high.x);
-	container->low.x = Min(box1->low.x, box2->low.x);
-	container->high.y = Max(box1->high.y, box2->high.y);
-	container->low.y = Min(box1->low.y, box2->low.y);
+	container->high.x = float8_max(box1->high.x, box2->high.x);
+	container->low.x = float8_min(box1->low.x, box2->low.x);
+	container->high.y = float8_max(box1->high.y, box2->high.y);
+	container->low.y = float8_min(box1->low.y, box2->low.y);
 
 	PG_RETURN_BOX_P(container);
 }
 
 
 /***********************************************************************
  **
  **		Routines for 2D paths.
  **
  ***********************************************************************/
@@ -4407,21 +4428,22 @@ circle_in(PG_FUNCTION_ARGS)
 		if (*cp == LDELIM)
 			s = cp;
 	}
 
 	pair_decode(s, &circle->center.x, &circle->center.y, &s, "circle", str);
 
 	if (*s == DELIM)
 		s++;
 
 	circle->radius = single_decode(s, &s, "circle", str);
-	if (circle->radius < 0)
+	/* We have to accept NaN. */
+	if (circle->radius < 0.0)
 		ereport(ERROR,
 				(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
 				 errmsg("invalid input syntax for type %s: \"%s\"",
 						"circle", str)));
 
 	while (depth > 0)
 	{
 		if ((*s == RDELIM) || ((*s == RDELIM_C) && (depth == 1)))
 		{
 			depth--;
@@ -4474,21 +4496,22 @@ circle_recv(PG_FUNCTION_ARGS)
 {
 	StringInfo	buf = (StringInfo) PG_GETARG_POINTER(0);
 	CIRCLE	   *circle;
 
 	circle = (CIRCLE *) palloc(sizeof(CIRCLE));
 
 	circle->center.x = pq_getmsgfloat8(buf);
 	circle->center.y = pq_getmsgfloat8(buf);
 	circle->radius = pq_getmsgfloat8(buf);
 
-	if (circle->radius < 0)
+	/* We have to accept NaN. */
+	if (circle->radius < 0.0)
 		ereport(ERROR,
 				(errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
 				 errmsg("invalid radius in external \"circle\" value")));
 
 	PG_RETURN_CIRCLE_P(circle);
 }
 
 /*
  *		circle_send			- converts circle to binary format
  */
@@ -4505,166 +4528,170 @@ circle_send(PG_FUNCTION_ARGS)
 	PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
 }
 
 
 /*----------------------------------------------------------
  *	Relational operators for CIRCLEs.
  *		<, >, <=, >=, and == are based on circle area.
  *---------------------------------------------------------*/
 
 /*		circles identical?
+ *
+ * We consider NaNs values to be equal to each other to let those circles
+ * to be found.
  */
 Datum
 circle_same(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle1 = PG_GETARG_CIRCLE_P(0);
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 
-	PG_RETURN_BOOL(FPeq(circle1->radius, circle2->radius) &&
+	PG_RETURN_BOOL(((isnan(circle1->radius) && isnan(circle1->radius)) ||
+					FPeq(circle1->radius, circle2->radius)) &&
 				   point_eq_point(&circle1->center, &circle2->center));
 }
 
 /*		circle_overlap	-		does circle1 overlap circle2?
  */
 Datum
 circle_overlap(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle1 = PG_GETARG_CIRCLE_P(0);
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 
 	PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
-						circle1->radius + circle2->radius));
+						float8_pl(circle1->radius, circle2->radius)));
 }
 
 /*		circle_overleft -		is the right edge of circle1 at or left of
  *								the right edge of circle2?
  */
 Datum
 circle_overleft(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle1 = PG_GETARG_CIRCLE_P(0);
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 
-	PG_RETURN_BOOL(FPle((circle1->center.x + circle1->radius),
-						(circle2->center.x + circle2->radius)));
+	PG_RETURN_BOOL(FPle(float8_pl(circle1->center.x, circle1->radius),
+						float8_pl(circle2->center.x, circle2->radius)));
 }
 
 /*		circle_left		-		is circle1 strictly left of circle2?
  */
 Datum
 circle_left(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle1 = PG_GETARG_CIRCLE_P(0);
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 
-	PG_RETURN_BOOL(FPlt((circle1->center.x + circle1->radius),
-						(circle2->center.x - circle2->radius)));
+	PG_RETURN_BOOL(FPlt(float8_pl(circle1->center.x, circle1->radius),
+						float8_mi(circle2->center.x, circle2->radius)));
 }
 
 /*		circle_right	-		is circle1 strictly right of circle2?
  */
 Datum
 circle_right(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle1 = PG_GETARG_CIRCLE_P(0);
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 
-	PG_RETURN_BOOL(FPgt((circle1->center.x - circle1->radius),
-						(circle2->center.x + circle2->radius)));
+	PG_RETURN_BOOL(FPgt(float8_mi(circle1->center.x, circle1->radius),
+						float8_pl(circle2->center.x, circle2->radius)));
 }
 
 /*		circle_overright	-	is the left edge of circle1 at or right of
  *								the left edge of circle2?
  */
 Datum
 circle_overright(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle1 = PG_GETARG_CIRCLE_P(0);
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 
-	PG_RETURN_BOOL(FPge((circle1->center.x - circle1->radius),
-						(circle2->center.x - circle2->radius)));
+	PG_RETURN_BOOL(FPge(float8_mi(circle1->center.x, circle1->radius),
+						float8_mi(circle2->center.x, circle2->radius)));
 }
 
 /*		circle_contained		-		is circle1 contained by circle2?
  */
 Datum
 circle_contained(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle1 = PG_GETARG_CIRCLE_P(0);
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 
 	PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
-						circle2->radius - circle1->radius));
+						float8_mi(circle2->radius, circle1->radius)));
 }
 
 /*		circle_contain	-		does circle1 contain circle2?
  */
 Datum
 circle_contain(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle1 = PG_GETARG_CIRCLE_P(0);
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 
 	PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
-						circle1->radius - circle2->radius));
+						float8_mi(circle1->radius, circle2->radius)));
 }
 
 
 /*		circle_below		-		is circle1 strictly below circle2?
  */
 Datum
 circle_below(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle1 = PG_GETARG_CIRCLE_P(0);
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 
-	PG_RETURN_BOOL(FPlt((circle1->center.y + circle1->radius),
-						(circle2->center.y - circle2->radius)));
+	PG_RETURN_BOOL(FPlt(float8_pl(circle1->center.y, circle1->radius),
+						float8_mi(circle2->center.y, circle2->radius)));
 }
 
 /*		circle_above	-		is circle1 strictly above circle2?
  */
 Datum
 circle_above(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle1 = PG_GETARG_CIRCLE_P(0);
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 
-	PG_RETURN_BOOL(FPgt((circle1->center.y - circle1->radius),
-						(circle2->center.y + circle2->radius)));
+	PG_RETURN_BOOL(FPgt(float8_mi(circle1->center.y, circle1->radius),
+						float8_pl(circle2->center.y, circle2->radius)));
 }
 
 /*		circle_overbelow -		is the upper edge of circle1 at or below
  *								the upper edge of circle2?
  */
 Datum
 circle_overbelow(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle1 = PG_GETARG_CIRCLE_P(0);
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 
-	PG_RETURN_BOOL(FPle((circle1->center.y + circle1->radius),
-						(circle2->center.y + circle2->radius)));
+	PG_RETURN_BOOL(FPle(float8_pl(circle1->center.y, circle1->radius),
+						float8_pl(circle2->center.y, circle2->radius)));
 }
 
 /*		circle_overabove	-	is the lower edge of circle1 at or above
  *								the lower edge of circle2?
  */
 Datum
 circle_overabove(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle1 = PG_GETARG_CIRCLE_P(0);
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 
-	PG_RETURN_BOOL(FPge((circle1->center.y - circle1->radius),
-						(circle2->center.y - circle2->radius)));
+	PG_RETURN_BOOL(FPge(float8_mi(circle1->center.y, circle1->radius),
+						float8_mi(circle2->center.y, circle2->radius)));
 }
 
 
 /*		circle_relop	-		is area(circle1) relop area(circle2), within
  *								our accuracy constraint?
  */
 Datum
 circle_eq(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle1 = PG_GETARG_CIRCLE_P(0);
@@ -4763,36 +4790,36 @@ circle_sub_pt(PG_FUNCTION_ARGS)
 Datum
 circle_mul_pt(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle = PG_GETARG_CIRCLE_P(0);
 	Point	   *point = PG_GETARG_POINT_P(1);
 	CIRCLE	   *result;
 
 	result = (CIRCLE *) palloc(sizeof(CIRCLE));
 
 	point_mul_point(&result->center, &circle->center, point);
-	result->radius = circle->radius * HYPOT(point->x, point->y);
+	result->radius = float8_mul(circle->radius, HYPOT(point->x, point->y));
 
 	PG_RETURN_CIRCLE_P(result);
 }
 
 Datum
 circle_div_pt(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle = PG_GETARG_CIRCLE_P(0);
 	Point	   *point = PG_GETARG_POINT_P(1);
 	CIRCLE	   *result;
 
 	result = (CIRCLE *) palloc(sizeof(CIRCLE));
 
 	point_div_point(&result->center, &circle->center, point);
-	result->radius = circle->radius / HYPOT(point->x, point->y);
+	result->radius = float8_div(circle->radius, HYPOT(point->x, point->y));
 
 	PG_RETURN_CIRCLE_P(result);
 }
 
 
 /*		circle_area		-		returns the area of the circle.
  */
 Datum
 circle_area(PG_FUNCTION_ARGS)
 {
@@ -4802,21 +4829,21 @@ circle_area(PG_FUNCTION_ARGS)
 }
 
 
 /*		circle_diameter -		returns the diameter of the circle.
  */
 Datum
 circle_diameter(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle = PG_GETARG_CIRCLE_P(0);
 
-	PG_RETURN_FLOAT8(2 * circle->radius);
+	PG_RETURN_FLOAT8(float8_mul(circle->radius, 2.0));
 }
 
 
 /*		circle_radius	-		returns the radius of the circle.
  */
 Datum
 circle_radius(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle = PG_GETARG_CIRCLE_P(0);
 
@@ -4827,81 +4854,85 @@ circle_radius(PG_FUNCTION_ARGS)
 /*		circle_distance -		returns the distance between
  *								  two circles.
  */
 Datum
 circle_distance(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle1 = PG_GETARG_CIRCLE_P(0);
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 	float8		result;
 
-	result = point_dt(&circle1->center, &circle2->center) -
-		(circle1->radius + circle2->radius);
-	if (result < 0)
-		result = 0;
+	result = float8_mi(point_dt(&circle1->center, &circle2->center),
+					   float8_pl(circle1->radius, circle2->radius));
+	if (result < 0.0)
+		result = 0.0;
+
 	PG_RETURN_FLOAT8(result);
 }
 
 
 Datum
 circle_contain_pt(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle = PG_GETARG_CIRCLE_P(0);
 	Point	   *point = PG_GETARG_POINT_P(1);
-	double		d;
+	float8		d;
 
 	d = point_dt(&circle->center, point);
 	PG_RETURN_BOOL(d <= circle->radius);
 }
 
 
 Datum
 pt_contained_circle(PG_FUNCTION_ARGS)
 {
 	Point	   *point = PG_GETARG_POINT_P(0);
 	CIRCLE	   *circle = PG_GETARG_CIRCLE_P(1);
-	double		d;
+	float8		d;
 
 	d = point_dt(&circle->center, point);
 	PG_RETURN_BOOL(d <= circle->radius);
 }
 
 
 /*		dist_pc -		returns the distance between
  *						  a point and a circle.
  */
 Datum
 dist_pc(PG_FUNCTION_ARGS)
 {
 	Point	   *point = PG_GETARG_POINT_P(0);
 	CIRCLE	   *circle = PG_GETARG_CIRCLE_P(1);
 	float8		result;
 
-	result = point_dt(point, &circle->center) - circle->radius;
-	if (result < 0)
-		result = 0;
+	result = float8_mi(point_dt(point, &circle->center),
+					   circle->radius);
+	if (result < 0.0)
+		result = 0.0;
+
 	PG_RETURN_FLOAT8(result);
 }
 
 /*
  * Distance from a circle to a point
  */
 Datum
 dist_cpoint(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle = PG_GETARG_CIRCLE_P(0);
 	Point	   *point = PG_GETARG_POINT_P(1);
 	float8		result;
 
-	result = point_dt(point, &circle->center) - circle->radius;
-	if (result < 0)
-		result = 0;
+	result = float8_mi(point_dt(point, &circle->center), circle->radius);
+	if (result < 0.0)
+		result = 0.0;
+
 	PG_RETURN_FLOAT8(result);
 }
 
 /*		circle_center	-		returns the center point of the circle.
  */
 Datum
 circle_center(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle = PG_GETARG_CIRCLE_P(0);
 	Point	   *result;
@@ -4909,24 +4940,24 @@ circle_center(PG_FUNCTION_ARGS)
 	result = (Point *) palloc(sizeof(Point));
 	result->x = circle->center.x;
 	result->y = circle->center.y;
 
 	PG_RETURN_POINT_P(result);
 }
 
 
 /*		circle_ar		-		returns the area of the circle.
  */
-static double
+static float8
 circle_ar(CIRCLE *circle)
 {
-	return M_PI * (circle->radius * circle->radius);
+	return float8_mul(float8_mul(circle->radius, circle->radius), M_PI);
 }
 
 
 /*----------------------------------------------------------
  *	Conversion operators.
  *---------------------------------------------------------*/
 
 Datum
 cr_circle(PG_FUNCTION_ARGS)
 {
@@ -4941,65 +4972,65 @@ cr_circle(PG_FUNCTION_ARGS)
 	result->radius = radius;
 
 	PG_RETURN_CIRCLE_P(result);
 }
 
 Datum
 circle_box(PG_FUNCTION_ARGS)
 {
 	CIRCLE	   *circle = PG_GETARG_CIRCLE_P(0);
 	BOX		   *box;
-	double		delta;
+	float8		delta;
 
 	box = (BOX *) palloc(sizeof(BOX));
 
-	delta = circle->radius / sqrt(2.0);
+	delta = float8_div(circle->radius, sqrt(2.0));
 
-	box->high.x = circle->center.x + delta;
-	box->low.x = circle->center.x - delta;
-	box->high.y = circle->center.y + delta;
-	box->low.y = circle->center.y - delta;
+	box->high.x = float8_pl(circle->center.x, delta);
+	box->low.x = float8_mi(circle->center.x, delta);
+	box->high.y = float8_pl(circle->center.y, delta);
+	box->low.y = float8_mi(circle->center.y, delta);
 
 	PG_RETURN_BOX_P(box);
 }
 
 /* box_circle()
  * Convert a box to a circle.
  */
 Datum
 box_circle(PG_FUNCTION_ARGS)
 {
 	BOX		   *box = PG_GETARG_BOX_P(0);
 	CIRCLE	   *circle;
 
 	circle = (CIRCLE *) palloc(sizeof(CIRCLE));
 
-	circle->center.x = (box->high.x + box->low.x) / 2;
-	circle->center.y = (box->high.y + box->low.y) / 2;
+	circle->center.x = float8_div(float8_pl(box->high.x, box->low.x), 2.0);
+	circle->center.y = float8_div(float8_pl(box->high.y, box->low.y), 2.0);
 
 	circle->radius = point_dt(&circle->center, &box->high);
 
 	PG_RETURN_CIRCLE_P(circle);
 }
 
 
 Datum
 circle_poly(PG_FUNCTION_ARGS)
 {
 	int32		npts = PG_GETARG_INT32(0);
 	CIRCLE	   *circle = PG_GETARG_CIRCLE_P(1);
 	POLYGON    *poly;
 	int			base_size,
 				size;
 	int			i;
-	double		angle;
-	double		anglestep;
+	float8		angle;
+	float8		anglestep;
 
 	if (FPzero(circle->radius))
 		ereport(ERROR,
 				(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
 				 errmsg("cannot convert circle with radius zero to polygon")));
 
 	if (npts < 2)
 		ereport(ERROR,
 				(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
 				 errmsg("must request at least 2 points")));
@@ -5010,27 +5041,30 @@ circle_poly(PG_FUNCTION_ARGS)
 	/* Check for integer overflow */
 	if (base_size / npts != sizeof(poly->p[0]) || size <= base_size)
 		ereport(ERROR,
 				(errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
 				 errmsg("too many points requested")));
 
 	poly = (POLYGON *) palloc0(size);	/* zero any holes */
 	SET_VARSIZE(poly, size);
 	poly->npts = npts;
 
-	anglestep = (2.0 * M_PI) / npts;
+	anglestep = float8_div(2.0 * M_PI, npts);
 
 	for (i = 0; i < npts; i++)
 	{
-		angle = i * anglestep;
-		poly->p[i].x = circle->center.x - (circle->radius * cos(angle));
-		poly->p[i].y = circle->center.y + (circle->radius * sin(angle));
+		angle = float8_mul(anglestep, i);
+
+		poly->p[i].x = float8_mi(circle->center.x,
+								 float8_mul(circle->radius, cos(angle)));
+		poly->p[i].y = float8_pl(circle->center.y,
+								 float8_mul(circle->radius, sin(angle)));
 	}
 
 	make_bound_box(poly);
 
 	PG_RETURN_POLYGON_P(poly);
 }
 
 /*
  * Convert polygon to circle
  *
@@ -5045,26 +5079,27 @@ poly_to_circle(CIRCLE *result, POLYGON *poly)
 	int			i;
 
 	Assert(poly->npts > 0);
 
 	result->center.x = 0;
 	result->center.y = 0;
 	result->radius = 0;
 
 	for (i = 0; i < poly->npts; i++)
 		point_add_point(&result->center, &result->center, &poly->p[i]);
-	result->center.x /= poly->npts;
-	result->center.y /= poly->npts;
+	result->center.x = float8_div(result->center.x, poly->npts);
+	result->center.y = float8_div(result->center.y, poly->npts);
 
 	for (i = 0; i < poly->npts; i++)
-		result->radius += point_dt(&poly->p[i], &result->center);
-	result->radius /= poly->npts;
+		result->radius = float8_pl(result->radius,
+								   point_dt(&poly->p[i], &result->center));
+	result->radius = float8_div(result->radius, poly->npts);
 }
 
 Datum
 poly_circle(PG_FUNCTION_ARGS)
 {
 	POLYGON    *poly = PG_GETARG_POLYGON_P(0);
 	CIRCLE	   *result;
 
 	result = (CIRCLE *) palloc(sizeof(CIRCLE));
 
@@ -5089,44 +5124,44 @@ poly_circle(PG_FUNCTION_ARGS)
  *	http://hopf.math.northwestern.edu/index.html
  *	Description of algorithm:  http://www.linuxjournal.com/article/2197
  *							   http://www.linuxjournal.com/article/2029
  */
 
 #define POINT_ON_POLYGON INT_MAX
 
 static int
 point_inside(Point *p, int npts, Point *plist)
 {
-	double		x0,
+	float8		x0,
 				y0;
-	double		prev_x,
+	float8		prev_x,
 				prev_y;
 	int			i = 0;
-	double		x,
+	float8		x,
 				y;
 	int			cross,
 				total_cross = 0;
 
 	Assert(npts > 0);
 
 	/* compute first polygon point relative to single point */
-	x0 = plist[0].x - p->x;
-	y0 = plist[0].y - p->y;
+	x0 = float8_mi(plist[0].x, p->x);
+	y0 = float8_mi(plist[0].y, p->y);
 
 	prev_x = x0;
 	prev_y = y0;
 	/* loop over polygon points and aggregate total_cross */
 	for (i = 1; i < npts; i++)
 	{
 		/* compute next polygon point relative to single point */
-		x = plist[i].x - p->x;
-		y = plist[i].y - p->y;
+		x = float8_mi(plist[i].x, p->x);
+		y = float8_mi(plist[i].y, p->y);
 
 		/* compute previous to current point crossing */
 		if ((cross = lseg_crossing(x, y, prev_x, prev_y)) == POINT_ON_POLYGON)
 			return 2;
 		total_cross += cross;
 
 		prev_x = x;
 		prev_y = y;
 	}
 
@@ -5144,69 +5179,74 @@ point_inside(Point *p, int npts, Point *plist)
 /* lseg_crossing()
  * Returns +/-2 if line segment crosses the positive X-axis in a +/- direction.
  * Returns +/-1 if one point is on the positive X-axis.
  * Returns 0 if both points are on the positive X-axis, or there is no crossing.
  * Returns POINT_ON_POLYGON if the segment contains (0,0).
  * Wow, that is one confusing API, but it is used above, and when summed,
  * can tell is if a point is in a polygon.
  */
 
 static int
-lseg_crossing(double x, double y, double prev_x, double prev_y)
+lseg_crossing(float8 x, float8 y, float8 prev_x, float8 prev_y)
 {
-	double		z;
+	float8		z;
 	int			y_sign;
 
 	if (FPzero(y))
 	{							/* y == 0, on X axis */
 		if (FPzero(x))			/* (x,y) is (0,0)? */
 			return POINT_ON_POLYGON;
 		else if (FPgt(x, 0))
 		{						/* x > 0 */
 			if (FPzero(prev_y)) /* y and prev_y are zero */
 				/* prev_x > 0? */
-				return FPgt(prev_x, 0) ? 0 : POINT_ON_POLYGON;
-			return FPlt(prev_y, 0) ? 1 : -1;
+				return FPgt(prev_x, 0.0) ? 0 : POINT_ON_POLYGON;
+			return FPlt(prev_y, 0.0) ? 1 : -1;
 		}
 		else
 		{						/* x < 0, x not on positive X axis */
 			if (FPzero(prev_y))
 				/* prev_x < 0? */
-				return FPlt(prev_x, 0) ? 0 : POINT_ON_POLYGON;
+				return FPlt(prev_x, 0.0) ? 0 : POINT_ON_POLYGON;
 			return 0;
 		}
 	}
 	else
 	{							/* y != 0 */
 		/* compute y crossing direction from previous point */
-		y_sign = FPgt(y, 0) ? 1 : -1;
+		y_sign = FPgt(y, 0.0) ? 1 : -1;
 
 		if (FPzero(prev_y))
 			/* previous point was on X axis, so new point is either off or on */
-			return FPlt(prev_x, 0) ? 0 : y_sign;
-		else if (FPgt(y_sign * prev_y, 0))
+			return FPlt(prev_x, 0.0) ? 0 : y_sign;
+		else if ((y_sign < 0 && FPlt(prev_y, 0.0)) ||
+				 (y_sign > 0 && FPgt(prev_y, 0.0)))
 			/* both above or below X axis */
 			return 0;			/* same sign */
 		else
 		{						/* y and prev_y cross X-axis */
-			if (FPge(x, 0) && FPgt(prev_x, 0))
+			if (FPge(x, 0.0) && FPgt(prev_x, 0.0))
 				/* both non-negative so cross positive X-axis */
 				return 2 * y_sign;
-			if (FPlt(x, 0) && FPle(prev_x, 0))
+			if (FPlt(x, 0.0) && FPle(prev_x, 0.0))
 				/* both non-positive so do not cross positive X-axis */
 				return 0;
 
 			/* x and y cross axises, see URL above point_inside() */
-			z = (x - prev_x) * y - (y - prev_y) * x;
+			z = float8_mi(float8_mul(float8_mi(x, prev_x), y),
+						  float8_mul(float8_mi(y, prev_y), x));
 			if (FPzero(z))
 				return POINT_ON_POLYGON;
-			return FPgt((y_sign * z), 0) ? 0 : 2 * y_sign;
+			if ((y_sign < 0 && FPlt(z, 0.0)) ||
+				(y_sign > 0 && FPgt(z, 0.0)))
+				return 0;
+			return 2 * y_sign;
 		}
 	}
 }
 
 
 static bool
 plist_same(int npts, Point *p1, Point *p2)
 {
 	int			i,
 				ii,
@@ -5276,47 +5316,52 @@ plist_same(int npts, Point *p1, Point *p2)
  *					 = x * sqrt( 1 + y^2/x^2 )
  *					 = x * sqrt( 1 + y/x * y/x )
  *
  * It is expected that this routine will eventually be replaced with the
  * C99 hypot() function.
  *
  * This implementation conforms to IEEE Std 1003.1 and GLIBC, in that the
  * case of hypot(inf,nan) results in INF, and not NAN.
  *-----------------------------------------------------------------------
  */
-double
-pg_hypot(double x, double y)
+float8
+pg_hypot(float8 x, float8 y)
 {
-	double		yx;
+	float8		yx,
+				result;
 
 	/* Handle INF and NaN properly */
 	if (isinf(x) || isinf(y))
 		return get_float8_infinity();
 
 	if (isnan(x) || isnan(y))
 		return get_float8_nan();
 
 	/* Else, drop any minus signs */
 	x = fabs(x);
 	y = fabs(y);
 
 	/* Swap x and y if needed to make x the larger one */
 	if (x < y)
 	{
-		double		temp = x;
+		float8		temp = x;
 
 		x = y;
 		y = temp;
 	}
 
 	/*
 	 * If y is zero, the hypotenuse is x.  This test saves a few cycles in
 	 * such cases, but more importantly it also protects against
 	 * divide-by-zero errors, since now x >= y.
 	 */
 	if (y == 0.0)
 		return x;
 
 	/* Determine the hypotenuse */
 	yx = y / x;
-	return x * sqrt(1.0 + (yx * yx));
+	result = x * sqrt(1.0 + (yx * yx));
+
+	check_float8_val(result, false, false);
+
+	return result;
 }
diff --git a/src/backend/utils/adt/geo_spgist.c b/src/backend/utils/adt/geo_spgist.c
index 096e59f817..b4008d7423 100644
--- a/src/backend/utils/adt/geo_spgist.c
+++ b/src/backend/utils/adt/geo_spgist.c
@@ -76,38 +76,38 @@
 #include "access/spgist.h"
 #include "access/stratnum.h"
 #include "catalog/pg_type.h"
 #include "utils/fmgrprotos.h"
 #include "utils/geo_decls.h"
 
 /*
  * Comparator for qsort
  *
  * We don't need to use the floating point macros in here, because this
- * is going only going to be used in a place to effect the performance
+ * is only going to be used in a place to effect the performance
  * of the index, not the correctness.
  */
 static int
 compareDoubles(const void *a, const void *b)
 {
-	double		x = *(double *) a;
-	double		y = *(double *) b;
+	float8		x = *(float8 *) a;
+	float8		y = *(float8 *) b;
 
 	if (x == y)
 		return 0;
 	return (x > y) ? 1 : -1;
 }
 
 typedef struct
 {
-	double		low;
-	double		high;
+	float8		low;
+	float8		high;
 } Range;
 
 typedef struct
 {
 	Range		left;
 	Range		right;
 } RangeBox;
 
 typedef struct
 {
@@ -167,21 +167,21 @@ getRangeBox(BOX *box)
 /*
  * Initialize the traversal value
  *
  * In the beginning, we don't have any restrictions.  We have to
  * initialize the struct to cover the whole 4D space.
  */
 static RectBox *
 initRectBox(void)
 {
 	RectBox    *rect_box = (RectBox *) palloc(sizeof(RectBox));
-	double		infinity = get_float8_infinity();
+	float8		infinity = get_float8_infinity();
 
 	rect_box->range_box_x.left.low = -infinity;
 	rect_box->range_box_x.left.high = infinity;
 
 	rect_box->range_box_x.right.low = -infinity;
 	rect_box->range_box_x.right.high = infinity;
 
 	rect_box->range_box_y.left.low = -infinity;
 	rect_box->range_box_y.left.high = infinity;
 
@@ -410,40 +410,40 @@ spg_box_quad_choose(PG_FUNCTION_ARGS)
  * point as the median of the coordinates of the boxes.
  */
 Datum
 spg_box_quad_picksplit(PG_FUNCTION_ARGS)
 {
 	spgPickSplitIn *in = (spgPickSplitIn *) PG_GETARG_POINTER(0);
 	spgPickSplitOut *out = (spgPickSplitOut *) PG_GETARG_POINTER(1);
 	BOX		   *centroid;
 	int			median,
 				i;
-	double	   *lowXs = palloc(sizeof(double) * in->nTuples);
-	double	   *highXs = palloc(sizeof(double) * in->nTuples);
-	double	   *lowYs = palloc(sizeof(double) * in->nTuples);
-	double	   *highYs = palloc(sizeof(double) * in->nTuples);
+	float8	   *lowXs = palloc(sizeof(float8) * in->nTuples);
+	float8	   *highXs = palloc(sizeof(float8) * in->nTuples);
+	float8	   *lowYs = palloc(sizeof(float8) * in->nTuples);
+	float8	   *highYs = palloc(sizeof(float8) * in->nTuples);
 
 	/* Calculate median of all 4D coordinates */
 	for (i = 0; i < in->nTuples; i++)
 	{
 		BOX		   *box = DatumGetBoxP(in->datums[i]);
 
 		lowXs[i] = box->low.x;
 		highXs[i] = box->high.x;
 		lowYs[i] = box->low.y;
 		highYs[i] = box->high.y;
 	}
 
-	qsort(lowXs, in->nTuples, sizeof(double), compareDoubles);
-	qsort(highXs, in->nTuples, sizeof(double), compareDoubles);
-	qsort(lowYs, in->nTuples, sizeof(double), compareDoubles);
-	qsort(highYs, in->nTuples, sizeof(double), compareDoubles);
+	qsort(lowXs, in->nTuples, sizeof(float8), compareDoubles);
+	qsort(highXs, in->nTuples, sizeof(float8), compareDoubles);
+	qsort(lowYs, in->nTuples, sizeof(float8), compareDoubles);
+	qsort(highYs, in->nTuples, sizeof(float8), compareDoubles);
 
 	median = in->nTuples / 2;
 
 	centroid = palloc(sizeof(BOX));
 
 	centroid->low.x = lowXs[median];
 	centroid->high.x = highXs[median];
 	centroid->low.y = lowYs[median];
 	centroid->high.y = highYs[median];
 
diff --git a/src/include/utils/geo_decls.h b/src/include/utils/geo_decls.h
index aeb31515e4..d14e8c8f72 100644
--- a/src/include/utils/geo_decls.h
+++ b/src/include/utils/geo_decls.h
@@ -1,42 +1,41 @@
 /*-------------------------------------------------------------------------
  *
  * geo_decls.h - Declarations for various 2D constructs.
  *
  *
  * Portions Copyright (c) 1996-2018, PostgreSQL Global Development Group
  * Portions Copyright (c) 1994, Regents of the University of California
  *
  * src/include/utils/geo_decls.h
  *
- * NOTE
- *	  These routines do *not* use the float types from adt/.
- *
  *	  XXX These routines were not written by a numerical analyst.
  *
  *	  XXX I have made some attempt to flesh out the operators
  *		and data types. There are still some more to do. - tgl 97/04/19
  *
  *-------------------------------------------------------------------------
  */
 #ifndef GEO_DECLS_H
 #define GEO_DECLS_H
 
 #include <math.h>
 
 #include "fmgr.h"
 #include "utils/float.h"
 
 /*--------------------------------------------------------------------
  * Useful floating point utilities and constants.
- *-------------------------------------------------------------------*/
-
+ *-------------------------------------------------------------------
+ *
+ * XXX They are not NaN-aware.
+ */
 
 #define EPSILON					1.0E-06
 
 #ifdef EPSILON
 #define FPzero(A)				(fabs(A) <= EPSILON)
 #define FPeq(A,B)				(fabs((A) - (B)) <= EPSILON)
 #define FPne(A,B)				(fabs((A) - (B)) > EPSILON)
 #define FPlt(A,B)				((B) - (A) > EPSILON)
 #define FPle(A,B)				((A) - (B) <= EPSILON)
 #define FPgt(A,B)				((A) - (B) > EPSILON)
@@ -51,21 +50,21 @@
 #define FPge(A,B)				((A) >= (B))
 #endif
 
 #define HYPOT(A, B)				pg_hypot(A, B)
 
 /*---------------------------------------------------------------------
  * Point - (x,y)
  *-------------------------------------------------------------------*/
 typedef struct
 {
-	double		x,
+	float8		x,
 				y;
 } Point;
 
 
 /*---------------------------------------------------------------------
  * LSEG - A straight line, specified by endpoints.
  *-------------------------------------------------------------------*/
 typedef struct
 {
 	Point		p[2];
@@ -83,21 +82,21 @@ typedef struct
 	int32		dummy;			/* padding to make it double align */
 	Point		p[FLEXIBLE_ARRAY_MEMBER];
 } PATH;
 
 
 /*---------------------------------------------------------------------
  * LINE - Specified by its general equation (Ax+By+C=0).
  *-------------------------------------------------------------------*/
 typedef struct
 {
-	double		A,
+	float8		A,
 				B,
 				C;
 } LINE;
 
 
 /*---------------------------------------------------------------------
  * BOX	- Specified by two corner points, which are
  *		 sorted to save calculation time later.
  *-------------------------------------------------------------------*/
 typedef struct
@@ -118,21 +117,21 @@ typedef struct
 	BOX			boundbox;
 	Point		p[FLEXIBLE_ARRAY_MEMBER];
 } POLYGON;
 
 /*---------------------------------------------------------------------
  * CIRCLE - Specified by a center point and radius.
  *-------------------------------------------------------------------*/
 typedef struct
 {
 	Point		center;
-	double		radius;
+	float8		radius;
 } CIRCLE;
 
 /*
  * fmgr interface macros
  *
  * Path and Polygon are toastable varlena types, the others are just
  * fixed-size pass-by-reference types.
  */
 
 #define DatumGetPointP(X)	 ((Point *) DatumGetPointer(X))
@@ -172,13 +171,13 @@ typedef struct
 #define DatumGetCircleP(X)	  ((CIRCLE *) DatumGetPointer(X))
 #define CirclePGetDatum(X)	  PointerGetDatum(X)
 #define PG_GETARG_CIRCLE_P(n) DatumGetCircleP(PG_GETARG_DATUM(n))
 #define PG_RETURN_CIRCLE_P(x) return CirclePGetDatum(x)
 
 
 /*
  * in geo_ops.c
  */
 
-extern double pg_hypot(double x, double y);
+extern float8 pg_hypot(float8 x, float8 y);
 
 #endif							/* GEO_DECLS_H */
-- 
2.14.3 (Apple Git-98)

