From 28a3446aea2c1a3b4ed2f1123b7c1e1ff385712d 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/4] geo-float-v04

Use the built-in float datatype to implement geometric types

This will provide:

* Check for underflow and overflow
* Check for division by zero
* Handle NaNs consistently

The patch also replaces all occurrences of "double" as "float8".  They
are the same, but were randomly spread on the same file.
---
 src/backend/access/gist/gistproc.c | 156 +++++------
 src/backend/utils/adt/geo_ops.c    | 533 ++++++++++++++++++++-----------------
 src/backend/utils/adt/geo_spgist.c |  36 +--
 src/include/utils/geo_decls.h      |  41 ++-
 4 files changed, 401 insertions(+), 365 deletions(-)

diff --git a/src/backend/access/gist/gistproc.c b/src/backend/access/gist/gistproc.c
index 13524b4da0..6fa97602c5 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,52 +364,54 @@ 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)
+	if (float4_gt(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.
 			 */
-			if (overlap < context->overlap ||
-				(overlap == context->overlap && ratio > context->ratio))
+			if (float4_lt(overlap, context->overlap) ||
+				(float4_eq(overlap, context->overlap) &&
+				 float4_gt(ratio, context->ratio)))
 				selectthis = true;
 		}
 		else
 		{
 			/*
 			 * Across dimensions, choose the new split if it has a smaller
 			 * *non-negative* overlap, or same *non-negative* overlap but
 			 * bigger range. This condition differs from the one described in
 			 * the article. On the datasets where leaf MBRs don't overlap
 			 * themselves, non-overlapping splits (i.e. splits which have zero
@@ -416,55 +420,49 @@ g_box_consider_split(ConsiderSplitContext *context, int dimNum,
 			 * non-overlapping splits (i.e. having lowest negative overlap)
 			 * appears to be in the same dimension as in the previous split.
 			 * Therefore MBRs appear to be very prolonged along another
 			 * dimension, which leads to bad search performance. Using range
 			 * as the second split criteria makes MBRs more quadratic. Using
 			 * *non-negative* overlap instead of overlap as the first split
 			 * criteria gives to range criteria a chance to matter, because
 			 * non-overlapping splits are equivalent in this criteria.
 			 */
 			if (non_negative(overlap) < non_negative(context->overlap) ||
-				(range > context->range &&
+				(float4_gt(range, context->range) &&
 				 non_negative(overlap) <= non_negative(context->overlap)))
 				selectthis = true;
 		}
 
 		if (selectthis)
 		{
 			/* save information about selected split */
 			context->first = false;
 			context->ratio = ratio;
 			context->range = range;
 			context->overlap = overlap;
 			context->rightLower = rightLower;
 			context->leftUpper = leftUpper;
 			context->dim = dimNum;
 		}
 	}
 }
 
 /*
  * Compare common entries by their deltas.
- * (We assume the deltas can't be NaN.)
  */
 static int
 common_entry_cmp(const void *i1, const void *i2)
 {
-	double		delta1 = ((const CommonEntry *) i1)->delta,
+	float8		delta1 = ((const CommonEntry *) i1)->delta,
 				delta2 = ((const CommonEntry *) i2)->delta;
 
-	if (delta1 < delta2)
-		return -1;
-	else if (delta1 > delta2)
-		return 1;
-	else
-		return 0;
+	return float8_cmp_internal(delta1, delta2);
 }
 
 /*
  * --------------------------------------------------------------------------
  * Double sorting split algorithm. This is used for both boxes and points.
  *
  * The algorithm finds split of boxes by considering splits along each axis.
  * Each entry is first projected as an interval on the X-axis, and different
  * ways to split the intervals into two groups are considered, trying to
  * minimize the overlap of the groups. Then the same is repeated for the
@@ -524,21 +522,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 +719,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 +774,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 +1098,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 +1139,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,80 +1207,84 @@ 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)
+	else if (float8_le(point->x, box->high.x) &&
+			 float8_ge(point->x, box->low.x) &&
+			 float8_le(point->y, box->high.y) &&
+			 float8_ge(point->y, box->low.y))
 	{
 		/* point inside the box */
 		result = 0.0;
 	}
-	else if (point->x <= box->high.x && point->x >= box->low.x)
+	else if (float8_le(point->x, box->high.x) &&
+			 float8_ge(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;
-		else if (point->y < box->low.y)
-			result = box->low.y - point->y;
+		Assert(float8_le(box->low.y, box->high.y));
+		if (float8_gt(point->y, box->high.y))
+			result = float8_mi(point->y, box->high.y);
+		else if (float8_lt(point->y, box->low.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)
+	else if (float8_le(point->y, box->high.y) &&
+			 float8_ge(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;
-		else if (point->x < box->low.x)
-			result = box->low.x - point->x;
+		Assert(float8_lt(box->low.x, box->high.x));
+		if (float8_gt(point->x, box->high.x))
+			result = float8_mi(point->x, box->high.x);
+		else if (float8_lt(point->x, box->low.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)
+		if (float8_gt(result, subresult))
 			result = subresult;
 
 		p.x = box->low.x;
 		p.y = box->high.y;
 		subresult = point_point_distance(point, &p);
-		if (result > subresult)
+		if (float8_gt(result, subresult))
 			result = subresult;
 
 		p.x = box->high.x;
 		p.y = box->low.y;
 		subresult = point_point_distance(point, &p);
-		if (result > subresult)
+		if (float8_gt(result, subresult))
 			result = subresult;
 	}
 
 	return result;
 }
 
 static bool
 gist_point_consistent_internal(StrategyNumber strategy,
 							   bool isLeaf, BOX *key, Point *query)
 {
@@ -1363,24 +1365,24 @@ gist_point_consistent(PG_FUNCTION_ARGS)
 				 * Instead we write a non-fuzzy overlap test.  The same code
 				 * will also serve for leaf-page tests, since leaf keys have
 				 * high == low.
 				 */
 				BOX		   *query,
 						   *key;
 
 				query = PG_GETARG_BOX_P(1);
 				key = DatumGetBoxP(entry->key);
 
-				result = (key->high.x >= query->low.x &&
-						  key->low.x <= query->high.x &&
-						  key->high.y >= query->low.y &&
-						  key->low.y <= query->high.y);
+				result = (float8_ge(key->high.x, query->low.x) &&
+						  float8_le(key->low.x, query->high.x) &&
+						  float8_ge(key->high.y, query->low.y) &&
+						  float8_le(key->low.y, query->high.y));
 				*recheck = false;
 			}
 			break;
 		case PolygonStrategyNumberGroup:
 			{
 				POLYGON    *query = PG_GETARG_POLYGON_P(1);
 
 				result = DatumGetBool(DirectFunctionCall5(
 														  gist_poly_consistent,
 														  PointerGetDatum(entry),
@@ -1389,22 +1391,22 @@ gist_point_consistent(PG_FUNCTION_ARGS)
 														  0, PointerGetDatum(recheck)));
 
 				if (GIST_LEAF(entry) && result)
 				{
 					/*
 					 * We are on leaf page and quick check shows overlapping
 					 * of polygon's bounding box and point
 					 */
 					BOX		   *box = DatumGetBoxP(entry->key);
 
-					Assert(box->high.x == box->low.x
-						   && box->high.y == box->low.y);
+					Assert(float8_eq(box->high.x, box->low.x) &&
+						   float8_eq(box->high.y, box->low.y));
 					result = DatumGetBool(DirectFunctionCall2(
 															  poly_contain_pt,
 															  PolygonPGetDatum(query),
 															  PointPGetDatum(&box->high)));
 					*recheck = false;
 				}
 			}
 			break;
 		case CircleStrategyNumberGroup:
 			{
@@ -1418,22 +1420,22 @@ gist_point_consistent(PG_FUNCTION_ARGS)
 														  0, PointerGetDatum(recheck)));
 
 				if (GIST_LEAF(entry) && result)
 				{
 					/*
 					 * We are on leaf page and quick check shows overlapping
 					 * of polygon's bounding box and point
 					 */
 					BOX		   *box = DatumGetBoxP(entry->key);
 
-					Assert(box->high.x == box->low.x
-						   && box->high.y == box->low.y);
+					Assert(float8_eq(box->high.x, box->low.x) &&
+						   float8_eq(box->high.y, box->low.y));
 					result = DatumGetBool(DirectFunctionCall2(
 															  circle_contain_pt,
 															  CirclePGetDatum(query),
 															  PointPGetDatum(&box->high)));
 					*recheck = false;
 				}
 			}
 			break;
 		default:
 			elog(ERROR, "unrecognized strategy number: %d", strategy);
@@ -1442,21 +1444,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 +1473,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 +1507,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 567554735f..dbf1de0ca2 100644
--- a/src/backend/utils/adt/geo_ops.c
+++ b/src/backend/utils/adt/geo_ops.c
@@ -29,76 +29,77 @@
  * Internal routines
  */
 
 enum path_delim
 {
 	PATH_NONE, PATH_OPEN, PATH_CLOSED
 };
 
 /* Routines for float8 */
 static inline float8 float8_slope(float8 x1, float8 x2, float8 y1, float8 y2);
+static inline float8 float8_hypot(float8 x, float8 y);
 
 /* Routines for two-dimensional points */
 static inline void point_construct(Point *result, float8 x, float8 y);
 static inline void point_add_internal(Point *result, Point *pt1, Point *pt2);
 static inline void point_sub_internal(Point *result, Point *pt1, Point *pt2);
 static inline void point_mul_internal(Point *result, Point *pt1, Point *pt2);
 static inline void point_div_internal(Point *result, Point *pt1, Point *pt2);
 static inline bool point_eq_internal(Point *pt1, Point *pt2);
 static float8 point_dt(Point *pt1, Point *pt2);
 static int	point_inside(Point *p, int npts, Point *plist);
 
 /* Routines for two-dimensional 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);
 
 /* Routines for two-dimensional circles */
-static double circle_ar(CIRCLE *circle);
+static float8 circle_ar(CIRCLE *circle);
 
 /* Routines for two-dimensional lines */
 static inline void line_construct_pm(LINE *result, Point *pt, float8 m);
 static inline void line_construct_pts(LINE *result, Point *pt1, Point *pt2);
 static bool line_interpt_internal(Point *result, LINE *l1, LINE *l2);
 static inline float8 line_calculate_point(LINE *line, Point *pt);
 
 /* Routines for two-dimensional line segments */
 static inline void statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2);
 static bool lseg_interpt_line_internal(Point *result, LSEG *lseg, LINE *line);
 static bool lseg_interpt_internal(Point *result, LSEG *l1, LSEG *l2);
 static bool lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start);
 static float8 lseg_dt(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 on_ps_internal(Point *pt, LSEG *lseg);
 
 /* Routines for two-dimensional polygons */
 static void make_bound_box(POLYGON *poly);
 
 /* Unsorted */
 static bool plist_same(int npts, Point *p1, Point *p2);
-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);
-static double dist_pl_internal(Point *pt, LINE *line);
-static double dist_ps_internal(Point *pt, LSEG *lseg);
-static double dist_ppoly_internal(Point *pt, POLYGON *poly);
+static float8 dist_pl_internal(Point *pt, LINE *line);
+static float8 dist_ps_internal(Point *pt, LSEG *lseg);
+static float8 dist_ppoly_internal(Point *pt, POLYGON *poly);
 
 
 /*
  * Delimiters for input and output strings.
  * LDELIM, RDELIM, and DELIM are left, right, and separator delimiters, respectively.
  * LDELIM_EP, RDELIM_EP are left and right delimiters for paths with endpoints.
  */
 
 #define LDELIM			'('
 #define RDELIM			')'
@@ -127,38 +128,38 @@ static double dist_ppoly_internal(Point *pt, POLYGON *poly);
  *
  * 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);
@@ -358,21 +359,42 @@ pair_count(char *s, char delim)
  * to calculate inverse slope.  To achieve that pass the values in
  * (y1, y2, x2, x1) order.
  */
 static inline float8
 float8_slope(float8 x1, float8 x2, float8 y1, float8 y2)
 {
 	if (FPeq(x1, x2))
 		return DBL_MAX;
 	if (FPeq(y1, y2))
 		return 0.0;
-	return (y1 - y2) / (x1 - x2);
+	return float8_div(float8_mi(y1, y2), float8_mi(x1, x2));
+}
+
+
+/*
+ * Wrapper around libc hypot()
+ */
+static inline float8
+float8_hypot(float8 x, float8 y)
+{
+	float8		result;
+
+	errno = 0;
+	result = hypot(x, y);
+	if (errno == ERANGE)
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("value out of range: overflow")));
+	check_float8_val(result, isinf(x) || isinf(y), x == 0.0 && y == 0.0);
+	Assert(result >= 0.0);
+
+	return result;
 }
 
 
 /***********************************************************************
  **
  **		Routines for two-dimensional boxes.
  **
  ***********************************************************************/
 
 /*----------------------------------------------------------
@@ -384,33 +406,33 @@ float8_slope(float8 x1, float8 x2, float8 y1, float8 y2)
  *		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.
@@ -424,38 +446,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);
 }
 
 /*
@@ -474,31 +496,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;
 	}
 }
@@ -810,54 +832,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.
@@ -867,24 +889,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
@@ -1015,30 +1037,30 @@ line_send(PG_FUNCTION_ARGS)
 
 /*
  * Fill already-allocated LINE struct from the point and the slope
  */
 static inline void
 line_construct_pm(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));
 	}
 }
 
 
 /*
  * Fill already-allocated LINE struct from two points on the line
  */
 static inline void
 line_construct_pts(LINE *result, Point *pt1, Point *pt2)
 {
@@ -1066,21 +1088,22 @@ line_construct_pp(PG_FUNCTION_ARGS)
 
 /*
  * Calculate the line equation for a point
  *
  * This returns the result of the line equation Ax + By + C.  The result
  * needs to be 0 for the point to be on the line.
  */
 static inline float8
 line_calculate_point(LINE *line, Point *pt)
 {
-	return line->A * pt->x + line->B * pt->y + line->C;
+	return float8_pl(float8_pl(float8_mul(line->A, pt->x),
+							   float8_mul(line->B, pt->y)), line->C);
 }
 
 
 /*----------------------------------------------------------
  *	Relative position routines.
  *---------------------------------------------------------*/
 
 Datum
 line_intersect(PG_FUNCTION_ARGS)
 {
@@ -1103,21 +1126,22 @@ 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));
 }
 
@@ -1127,34 +1151,34 @@ line_horizontal(PG_FUNCTION_ARGS)
 	LINE	   *line = PG_GETARG_LINE_P(0);
 
 	PG_RETURN_BOOL(FPzero(line->A));
 }
 
 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;
+		ratio = float8_div(l1->A, l2->A);
 	else if (!FPzero(l2->B))
-		k = l1->B / l2->B;
+		ratio = float8_div(l1->B, l2->B);
 	else if (!FPzero(l2->C))
-		k = l1->C / 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)));
 }
 
 
 /*----------------------------------------------------------
  *	Line arithmetic routines.
  *---------------------------------------------------------*/
 
 /* line_distance()
  * Distance between two lines.
  */
@@ -1162,21 +1186,21 @@ 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_internal(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 = dist_pl_internal(&tmp, l2);
 	PG_RETURN_FLOAT8(result);
 }
 
 /* line_interpt()
  * Point where two lines l1, l2 intersect (if any)
  */
 Datum
 line_interpt(PG_FUNCTION_ARGS)
@@ -1199,41 +1223,41 @@ line_interpt(PG_FUNCTION_ARGS)
  * 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.
  */
 static bool
 line_interpt_internal(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 (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;
 
 		if (FPzero(l2->B))
 			x = l2->C;
 		else
-			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);
 
 #ifdef GEODEBUG
 	printf("line_interpt- lines are A=%.*g, B=%.*g, C=%.*g, A=%.*g, B=%.*g, C=%.*g\n",
 		   DBL_DIG, l1->A, DBL_DIG, l1->B, DBL_DIG, l1->C, DBL_DIG, l2->A, DBL_DIG, l2->B, DBL_DIG, l2->C);
 	printf("line_interpt- lines intersect at (%.*g,%.*g)\n", DBL_DIG, x, DBL_DIG, y);
 #endif
 
@@ -1260,36 +1284,35 @@ line_interpt_internal(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;
@@ -1506,31 +1529,31 @@ path_npoints(PG_FUNCTION_ARGS)
 
 	PG_RETURN_INT32(path->npts);
 }
 
 
 Datum
 path_close(PG_FUNCTION_ARGS)
 {
 	PATH	   *path = PG_GETARG_PATH_P_COPY(0);
 
-	path->closed = TRUE;
+	path->closed = true;
 
 	PG_RETURN_PATH_P(path);
 }
 
 Datum
 path_open(PG_FUNCTION_ARGS)
 {
 	PATH	   *path = PG_GETARG_PATH_P_COPY(0);
 
-	path->closed = FALSE;
+	path->closed = false;
 
 	PG_RETURN_PATH_P(path);
 }
 
 
 /* path_inter -
  *		Does p1 intersect p2 at any point?
  *		Use bounding boxes for a quick (O(n)) check, then do a
  *		O(n^2) iterative edge check.
  */
@@ -1546,33 +1569,33 @@ path_inter(PG_FUNCTION_ARGS)
 	LSEG		seg1,
 				seg2;
 
 	if (p1->npts <= 0 || p2->npts <= 0)
 		PG_RETURN_BOOL(false);
 
 	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)
@@ -1648,21 +1671,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_dt(&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();
 
@@ -1687,21 +1710,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.
  **
  ***********************************************************************/
@@ -1875,21 +1898,21 @@ point_distance(PG_FUNCTION_ARGS)
 	Point	   *pt2 = PG_GETARG_POINT_P(1);
 
 	PG_RETURN_FLOAT8(point_dt(pt1, pt2));
 }
 
 static float8
 point_dt(Point *pt1, Point *pt2)
 {
 	float8		result;
 
-	result = hypot(pt1->x - pt2->x, pt1->y - pt2->y);
+	result = float8_hypot(float8_mi(pt1->x, pt2->x), float8_mi(pt1->y, pt2->y));
 
 #ifdef GEODEBUG
 	printf("point_dt- segment (%f,%f),(%f,%f) length is %f\n",
 		   pt1->x, pt1->y, pt2->x, pt2->y, result);
 #endif
 
 	return result;
 }
 
 Datum
@@ -2049,35 +2072,35 @@ lseg_parallel(PG_FUNCTION_ARGS)
  *	'((0,0),(0,1))'::lseg ?-| '((0,0),(1,0))'::lseg
  * So, modified it to check explicitly for slope of vertical line
  *	returned by float8_slope() and the results seem better.
  * - thomas 1998-01-31
  */
 Datum
 lseg_perp(PG_FUNCTION_ARGS)
 {
 	LSEG	   *l1 = PG_GETARG_LSEG_P(0);
 	LSEG	   *l2 = PG_GETARG_LSEG_P(1);
-	double		m1,
+	float8		m1,
 				m2;
 
 	m1 = float8_slope(l1->p[0].x, l1->p[1].x, l1->p[0].y, l1->p[1].y);
 	m2 = float8_slope(l2->p[0].x, l2->p[1].x, l2->p[0].y, l2->p[1].y);
 
 #ifdef GEODEBUG
 	printf("lseg_perp- slopes are %g and %g\n", m1, m2);
 #endif
 	if (FPzero(m1))
 		PG_RETURN_BOOL(FPeq(m2, DBL_MAX));
 	else if (FPzero(m2))
 		PG_RETURN_BOOL(FPeq(m1, DBL_MAX));
 
-	PG_RETURN_BOOL(FPeq(m1 / m2, -1.0));
+	PG_RETURN_BOOL(FPeq(float8_div(m1, m2), -1.0));
 }
 
 Datum
 lseg_vertical(PG_FUNCTION_ARGS)
 {
 	LSEG	   *lseg = PG_GETARG_LSEG_P(0);
 
 	PG_RETURN_BOOL(FPeq(lseg->p[0].x, lseg->p[1].x));
 }
 
@@ -2167,52 +2190,47 @@ lseg_distance(PG_FUNCTION_ARGS)
 	LSEG	   *l2 = PG_GETARG_LSEG_P(1);
 
 	PG_RETURN_FLOAT8(lseg_dt(l1, l2));
 }
 
 /* lseg_dt()
  * Distance between two line segments.
  * Must check both sets of endpoints to ensure minimum distance is found.
  * - thomas 1998-02-01
  */
-static double
+static float8
 lseg_dt(LSEG *l1, LSEG *l2)
 {
-	double		result,
-				d;
+	float8		result;
 
 	if (lseg_interpt_internal(NULL, l1, l2))
 		return 0.0;
 
-	d = dist_ps_internal(&l1->p[0], l2);
-	result = d;
-	d = dist_ps_internal(&l1->p[1], l2);
-	result = Min(result, d);
-	d = dist_ps_internal(&l2->p[0], l1);
-	result = Min(result, d);
-	d = dist_ps_internal(&l2->p[1], l1);
-	result = Min(result, d);
+	result = dist_ps_internal(&l1->p[0], l2);
+	result = float8_min(result, dist_ps_internal(&l1->p[1], l2));
+	result = float8_min(result, dist_ps_internal(&l2->p[0], l1));
+	result = float8_min(result, dist_ps_internal(&l2->p[1], l1));
 
 	return result;
 }
 
 
 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);
 }
 
 static bool
 lseg_interpt_internal(Point *result, LSEG *l1, LSEG *l2)
 {
 	Point		interpt;
 	LINE		tmp1,
 				tmp2;
@@ -2296,45 +2314,44 @@ lseg_interpt(PG_FUNCTION_ARGS)
  */
 Datum
 dist_pl(PG_FUNCTION_ARGS)
 {
 	Point	   *pt = PG_GETARG_POINT_P(0);
 	LINE	   *line = PG_GETARG_LINE_P(1);
 
 	PG_RETURN_FLOAT8(dist_pl_internal(pt, line));
 }
 
-static double
+static float8
 dist_pl_internal(Point *pt, LINE *line)
 {
-	return fabs(line_calculate_point(line, pt) /
-				hypot(line->A, line->B));
+	return float8_div(fabs(line_calculate_point(line, pt)),
+					  float8_hypot(line->A, line->B));
 }
 
 /*
  * Distance from a point to a lseg
  */
 Datum
 dist_ps(PG_FUNCTION_ARGS)
 {
 	Point	   *pt = PG_GETARG_POINT_P(0);
 	LSEG	   *lseg = PG_GETARG_LSEG_P(1);
 
 	PG_RETURN_FLOAT8(dist_ps_internal(pt, lseg));
 }
 
-static double
+static float8
 dist_ps_internal(Point *pt, LSEG *lseg)
 {
-	double		m;				/* slope of perp. */
-	double		result,
-				tmpdist;
+	float8		m;				/* slope of perp. */
+	float8		result;
 	Point		interpt;
 	LINE		ln;
 
 	/*
 	 * Construct a line perpendicular to the input segment and through the
 	 * input point
 	 */
 	m = float8_slope(lseg->p[0].y, lseg->p[1].y, lseg->p[1].x, lseg->p[0].x);
 	line_construct_pm(&ln, pt, m);
 
@@ -2354,24 +2371,22 @@ dist_ps_internal(Point *pt, LSEG *lseg)
 		/* yes, so use distance to the intersection point */
 		result = point_dt(pt, &interpt);
 #ifdef GEODEBUG
 		printf("dist_ps- distance is %f to intersection point is (%f,%f)\n",
 			   result, tmp.x, tmp.y);
 #endif
 	}
 	else
 	{
 		/* no, so use distance to the nearer endpoint */
-		result = point_dt(pt, &lseg->p[0]);
-		tmpdist = point_dt(pt, &lseg->p[1]);
-		if (tmpdist < result)
-			result = tmpdist;
+		result = float8_min(point_dt(pt, &lseg->p[0]),
+							point_dt(pt, &lseg->p[1]));
 	}
 
 	return result;
 }
 
 /*
  * Distance from a point to a path
  */
 Datum
 dist_ppath(PG_FUNCTION_ARGS)
@@ -2409,21 +2424,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 = dist_ps_internal(pt, &lseg);
-				if (!have_min || tmp < result)
+				if (!have_min || float8_lt(tmp, result))
 				{
 					result = tmp;
 					have_min = true;
 				}
 			}
 			break;
 	}
 	PG_RETURN_FLOAT8(result);
 }
 
@@ -2447,33 +2462,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_internal(NULL, lseg, line))
 		result = 0.0;
 	else
-	{
-		result = dist_pl_internal(&lseg->p[0], line);
-		d2 = dist_pl_internal(&lseg->p[1], line);
 		/* XXX shouldn't we take the min not max? */
-		if (d2 > result)
-			result = d2;
-	}
+		result = float8_max(dist_pl_internal(&lseg->p[0], line),
+							dist_pl_internal(&lseg->p[1], line));
 
 	PG_RETURN_FLOAT8(result);
 }
 
 /*
  * Distance from a lseg to a box
  */
 Datum
 dist_sb(PG_FUNCTION_ARGS)
 {
@@ -2515,25 +2525,22 @@ 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_min(float8_mi(dist_ppoly_internal(&circle->center, poly),
+								  circle->radius), 0.0);
 
 	PG_RETURN_FLOAT8(result);
 }
 
 /*
  * Distance from a point to a polygon
  */
 Datum
 dist_ppoly(PG_FUNCTION_ARGS)
 {
@@ -2551,21 +2558,21 @@ dist_polyp(PG_FUNCTION_ARGS)
 {
 	POLYGON    *poly = PG_GETARG_POLYGON_P(0);
 	Point	   *point = PG_GETARG_POINT_P(1);
 	float8		result;
 
 	result = dist_ppoly_internal(point, poly);
 
 	PG_RETURN_FLOAT8(result);
 }
 
-static double
+static float8
 dist_ppoly_internal(Point *pt, POLYGON *poly)
 {
 	float8		result;
 	float8		d;
 	int			i;
 	LSEG		seg;
 
 	if (point_inside(pt, poly->npts, poly->p) != 0)
 	{
 #ifdef GEODEBUG
@@ -2588,21 +2595,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 = dist_ps_internal(pt, &seg);
 #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.
@@ -2656,41 +2663,41 @@ lseg_interpt_line_internal(Point *result, LSEG *lseg, LINE *line)
 /* close_pl -
  *		The intersection point of a perpendicular of the line
  *		through the point.
  */
 Datum
 close_pl(PG_FUNCTION_ARGS)
 {
 	Point	   *pt = PG_GETARG_POINT_P(0);
 	LINE	   *line = PG_GETARG_LINE_P(1);
 	Point	   *result;
-	double		invm;
+	float8		invm;
 	LINE		tmp;
 
 	result = (Point *) palloc(sizeof(Point));
 
 	if (FPzero(line->B))		/* vertical? */
 	{
 		result->x = line->C;
 		result->y = pt->y;
 		PG_RETURN_POINT_P(result);
 	}
 	if (FPzero(line->A))		/* horizontal? */
 	{
 		result->x = pt->x;
 		result->y = line->C;
 		PG_RETURN_POINT_P(result);
 	}
 	/* drop a perpendicular and find the intersection point */
 
 	/* invert and flip the sign on the slope to get a perpendicular */
-	invm = line->B / line->A;
+	invm = float8_div(line->B, line->A);
 	line_construct_pm(&tmp, pt, invm);
 	line_interpt_internal(result, &tmp, line);
 	PG_RETURN_POINT_P(result);
 }
 
 
 /* close_ps()
  * Closest point on line segment to specified point.
  * Take the closest endpoint if the point is left, right,
  *	above, or below the segment, otherwise find the intersection
@@ -2699,21 +2706,21 @@ close_pl(PG_FUNCTION_ARGS)
  * Some tricky code here, relying on boolean expressions
  *	evaluating to only zero or one to use as an array index.
  *		bug fixes by gthaker@atl.lmco.com; May 1, 1998
  */
 Datum
 close_ps(PG_FUNCTION_ARGS)
 {
 	Point	   *pt = PG_GETARG_POINT_P(0);
 	LSEG	   *lseg = PG_GETARG_LSEG_P(1);
 	Point	   *result;
-	double		invm;
+	float8		invm;
 	int			xh,
 				yh;
 	LINE		tmp;
 
 	result = (Point *) palloc(sizeof(Point));
 
 #ifdef GEODEBUG
 	printf("close_sp:pt->x %f pt->y %f\nlseg(0).x %f lseg(0).y %f  lseg(1).x %f lseg(1).y %f\n",
 		   pt->x, pt->y, lseg->p[0].x, lseg->p[0].y,
 		   lseg->p[1].x, lseg->p[1].y);
@@ -2758,32 +2765,32 @@ close_ps(PG_FUNCTION_ARGS)
 	}
 
 	/*
 	 * vert. and horiz. cases are down, now check if the closest point is one
 	 * of the end points or someplace on the lseg.
 	 */
 
 	invm = float8_slope(lseg->p[0].y, lseg->p[1].y, lseg->p[1].x, lseg->p[0].x);
 	line_construct_pm(&tmp, &lseg->p[!yh], invm);	/* lower edge of the
 													 * "band" */
-	if (pt->y < (tmp.A * pt->x + tmp.C))
+	if (pt->y < float8_pl(float8_mul(tmp.A, pt->x), tmp.C))
 	{							/* we are below the lower edge */
 		*result = lseg->p[!yh]; /* below the lseg, take lower end pt */
 #ifdef GEODEBUG
 		printf("close_ps below: tmp A %f  B %f   C %f\n",
 			   tmp->A, tmp->B, tmp->C);
 #endif
 		PG_RETURN_POINT_P(result);
 	}
 	line_construct_pm(&tmp, &lseg->p[yh], invm);	/* upper edge of the
 													 * "band" */
-	if (pt->y > (tmp.A * pt->x + tmp.C))
+	if (pt->y > float8_pl(float8_mul(tmp.A, pt->x), tmp.C))
 	{							/* we are below the lower edge */
 		*result = lseg->p[yh];	/* above the lseg, take higher end pt */
 #ifdef GEODEBUG
 		printf("close_ps above: tmp A %f  B %f   C %f\n",
 			   tmp->A, tmp->B, tmp->C);
 #endif
 		PG_RETURN_POINT_P(result);
 	}
 
 	/*
@@ -2819,32 +2826,32 @@ close_lseg(PG_FUNCTION_ARGS)
 {
 	LSEG	   *l1 = PG_GETARG_LSEG_P(0);
 	LSEG	   *l2 = PG_GETARG_LSEG_P(1);
 	Point	   *result;
 	double		dist,
 				dist0,
 				dist1;
 
 	dist0 = dist_ps_internal(&l1->p[0], l2);
 	dist1 = dist_ps_internal(&l1->p[1], l2);
-	dist = Min(dist0, dist1);
+	dist = float8_min(dist0, dist1);
 
-	if (dist_ps_internal(&l2->p[0], l1) < dist)
+	if (float8_lt(dist_ps_internal(&l2->p[0], l1), dist))
 	{
 		result = DatumGetPointP(DirectFunctionCall2(close_ps,
 													PointPGetDatum(&l2->p[0]),
 													LsegPGetDatum(l1)));
 		result = DatumGetPointP(DirectFunctionCall2(close_ps,
 													PointPGetDatum(result),
 													LsegPGetDatum(l2)));
 	}
-	else if (dist_ps_internal(&l2->p[1], l1) < dist)
+	else if (float8_lt(dist_ps_internal(&l2->p[1], l1), dist))
 	{
 		result = DatumGetPointP(DirectFunctionCall2(close_ps,
 													PointPGetDatum(&l2->p[1]),
 													LsegPGetDatum(l1)));
 		result = DatumGetPointP(DirectFunctionCall2(close_ps,
 													PointPGetDatum(result),
 													LsegPGetDatum(l2)));
 	}
 	else
 	{
@@ -2859,52 +2866,55 @@ close_lseg(PG_FUNCTION_ARGS)
  * Closest point on or in box to specified point.
  */
 Datum
 close_pb(PG_FUNCTION_ARGS)
 {
 	Point	   *pt = PG_GETARG_POINT_P(0);
 	BOX		   *box = PG_GETARG_BOX_P(1);
 	LSEG		lseg,
 				seg;
 	Point		point;
-	double		dist,
+	float8		dist,
 				d;
 
 	if (DatumGetBool(DirectFunctionCall2(on_pb,
 										 PointPGetDatum(pt),
 										 BoxPGetDatum(box))))
 		PG_RETURN_POINT_P(pt);
 
 	/* pairwise check lseg distances */
 	point.x = box->low.x;
 	point.y = box->high.y;
 	statlseg_construct(&lseg, &box->low, &point);
 	dist = dist_ps_internal(pt, &lseg);
 
 	statlseg_construct(&seg, &box->high, &point);
-	if ((d = dist_ps_internal(pt, &seg)) < dist)
+	d = dist_ps_internal(pt, &seg);
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		memcpy(&lseg, &seg, sizeof(lseg));
 	}
 
 	point.x = box->high.x;
 	point.y = box->low.y;
 	statlseg_construct(&seg, &box->low, &point);
-	if ((d = dist_ps_internal(pt, &seg)) < dist)
+	d = dist_ps_internal(pt, &seg);
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		memcpy(&lseg, &seg, sizeof(lseg));
 	}
 
 	statlseg_construct(&seg, &box->high, &point);
-	if ((d = dist_ps_internal(pt, &seg)) < dist)
+	d = dist_ps_internal(pt, &seg);
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		memcpy(&lseg, &seg, sizeof(lseg));
 	}
 
 	PG_RETURN_DATUM(DirectFunctionCall2(close_ps,
 										PointPGetDatum(pt),
 										LsegPGetDatum(&lseg)));
 }
 
@@ -2927,21 +2937,21 @@ close_sl(PG_FUNCTION_ARGS)
 	float8		d1,
 				d2;
 
 	result = (Point *) palloc(sizeof(Point));
 
 	if (lseg_interpt_line_internal(result, lseg, line))
 		PG_RETURN_POINT_P(result);
 
 	d1 = dist_pl_internal(&lseg->p[0], line);
 	d2 = dist_pl_internal(&lseg->p[1], line);
-	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")));
@@ -2961,77 +2971,80 @@ close_ls(PG_FUNCTION_ARGS)
 	float8		d1,
 				d2;
 
 	result = (Point *) palloc(sizeof(Point));
 
 	if (lseg_interpt_line_internal(result, lseg, line))
 		PG_RETURN_POINT_P(result);
 
 	d1 = dist_pl_internal(&lseg->p[0], line);
 	d2 = dist_pl_internal(&lseg->p[1], line);
-	if (d1 < d2)
+	if (float8_lt(d1, d2))
 		*result = lseg->p[0];
 	else
 		*result = lseg->p[1];
 
 	PG_RETURN_POINT_P(result);
 }
 
 /* close_sb()
  * Closest point on or in box to line segment.
  */
 Datum
 close_sb(PG_FUNCTION_ARGS)
 {
 	LSEG	   *lseg = PG_GETARG_LSEG_P(0);
 	BOX		   *box = PG_GETARG_BOX_P(1);
 	Point		point;
 	LSEG		bseg,
 				seg;
-	double		dist,
+	float8		dist,
 				d;
 
 	/* segment intersects box? then just return closest point to center */
 	if (DatumGetBool(DirectFunctionCall2(inter_sb,
 										 LsegPGetDatum(lseg),
 										 BoxPGetDatum(box))))
 	{
 		box_cn(&point, box);
 		PG_RETURN_DATUM(DirectFunctionCall2(close_ps,
 											PointPGetDatum(&point),
 											LsegPGetDatum(lseg)));
 	}
 
 	/* pairwise check lseg distances */
 	point.x = box->low.x;
 	point.y = box->high.y;
 	statlseg_construct(&bseg, &box->low, &point);
 	dist = lseg_dt(lseg, &bseg);
 
 	statlseg_construct(&seg, &box->high, &point);
-	if ((d = lseg_dt(lseg, &seg)) < dist)
+	d = lseg_dt(lseg, &seg);
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		memcpy(&bseg, &seg, sizeof(bseg));
 	}
 
 	point.x = box->high.x;
 	point.y = box->low.y;
 	statlseg_construct(&seg, &box->low, &point);
-	if ((d = lseg_dt(lseg, &seg)) < dist)
+	d = lseg_dt(lseg, &seg);
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		memcpy(&bseg, &seg, sizeof(bseg));
 	}
 
 	statlseg_construct(&seg, &box->high, &point);
-	if ((d = lseg_dt(lseg, &seg)) < dist)
+	d = lseg_dt(lseg, &seg);
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		memcpy(&bseg, &seg, sizeof(bseg));
 	}
 
 	/* OK, we now have the closest line segment on the box boundary */
 	PG_RETURN_DATUM(DirectFunctionCall2(close_lseg,
 										LsegPGetDatum(lseg),
 										LsegPGetDatum(&bseg)));
 }
@@ -3079,50 +3092,55 @@ on_ps(PG_FUNCTION_ARGS)
 {
 	Point	   *pt = PG_GETARG_POINT_P(0);
 	LSEG	   *lseg = PG_GETARG_LSEG_P(1);
 
 	PG_RETURN_BOOL(on_ps_internal(pt, lseg));
 }
 
 static bool
 on_ps_internal(Point *pt, LSEG *lseg)
 {
-	return FPeq(point_dt(pt, &lseg->p[0]) + point_dt(pt, &lseg->p[1]),
+	return FPeq(float8_pl(point_dt(pt, &lseg->p[0]),
+						  point_dt(pt, &lseg->p[1])),
 				point_dt(&lseg->p[0], &lseg->p[1]));
 }
 
 
 /*
  * Check whether the point is in the box or on its border
  */
 Datum
 on_pb(PG_FUNCTION_ARGS)
 {
 	Point	   *pt = PG_GETARG_POINT_P(0);
 	BOX		   *box = PG_GETARG_BOX_P(1);
 
-	PG_RETURN_BOOL(pt->x <= box->high.x && pt->x >= box->low.x &&
-				   pt->y <= box->high.y && pt->y >= box->low.y);
+	PG_RETURN_BOOL(float8_le(pt->x, box->high.x) &&
+				   float8_ge(pt->x, box->low.x) &&
+				   float8_le(pt->y, box->high.y) &&
+				   float8_ge(pt->y, box->low.y));
 }
 
 
 /*
  * Commutator of on_pb()
  */
 Datum
 box_contain_pt(PG_FUNCTION_ARGS)
 {
 	BOX		   *box = PG_GETARG_BOX_P(0);
 	Point	   *pt = PG_GETARG_POINT_P(1);
 
-	PG_RETURN_BOOL(pt->x <= box->high.x && pt->x >= box->low.x &&
-				   pt->y <= box->high.y && pt->y >= box->low.y);
+	PG_RETURN_BOOL(float8_le(pt->x, box->high.x) &&
+				   float8_ge(pt->x, box->low.x) &&
+				   float8_le(pt->y, box->high.y) &&
+				   float8_ge(pt->y, box->low.y));
 }
 
 
 /* on_ppath -
  *		Whether a point lies within (on) a polyline.
  *		If open, we have to (groan) check each segment.
  * (uses same algorithm as for point intersecting segment - tgl 1997-07-09)
  *		If closed, we use the old O(n) ray method for point-in-polygon.
  *				The ray is horizontal, from pt out to the right.
  *				Each segment that crosses the ray counts as an
@@ -3130,33 +3148,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);
 }
 
@@ -3226,24 +3243,24 @@ inter_sl(PG_FUNCTION_ARGS)
  */
 Datum
 inter_sb(PG_FUNCTION_ARGS)
 {
 	LSEG	   *lseg = PG_GETARG_LSEG_P(0);
 	BOX		   *box = PG_GETARG_BOX_P(1);
 	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))
 		PG_RETURN_BOOL(false);
 
 	/* an endpoint of segment is inside box? then clearly intersects */
 	if (DatumGetBool(DirectFunctionCall2(on_pb,
 										 PointPGetDatum(&lseg->p[0]),
 										 BoxPGetDatum(box))) ||
 		DatumGetBool(DirectFunctionCall2(on_pb,
@@ -3332,27 +3349,27 @@ make_bound_box(POLYGON *poly)
 {
 	int			i;
 	BOX		   *boundbox = &poly->boundbox;
 
 	if (poly->npts > 0)
 	{
 		boundbox->low.x = boundbox->high.x = poly->p[0].x;
 		boundbox->low.y = boundbox->high.y = poly->p[0].y;
 		for (i = 1; i < poly->npts; i++)
 		{
-			if (poly->p[i].x < boundbox->low.x)
+			if (float8_lt(poly->p[i].x, boundbox->low.x))
 				boundbox->low.x = poly->p[i].x;
-			if (poly->p[i].x > boundbox->high.x)
+			if (float8_gt(poly->p[i].x, boundbox->high.x))
 				boundbox->high.x = poly->p[i].x;
-			if (poly->p[i].y < boundbox->low.y)
+			if (float8_lt(poly->p[i].y, boundbox->low.y))
 				boundbox->low.y = poly->p[i].y;
-			if (poly->p[i].y > boundbox->high.y)
+			if (float8_gt(poly->p[i].y, boundbox->high.y))
 				boundbox->high.y = poly->p[i].y;
 		}
 	}
 	else
 		ereport(ERROR,
 				(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
 				 errmsg("cannot create bounding box for empty polygon")));
 }
 
 /*------------------------------------------------------------------
@@ -3478,21 +3495,21 @@ poly_send(PG_FUNCTION_ARGS)
  * the right most point of A left of the left most point
  * of B?
  *-------------------------------------------------------*/
 Datum
 poly_left(PG_FUNCTION_ARGS)
 {
 	POLYGON    *polya = PG_GETARG_POLYGON_P(0);
 	POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
 	bool		result;
 
-	result = polya->boundbox.high.x < polyb->boundbox.low.x;
+	result = float8_lt(polya->boundbox.high.x, polyb->boundbox.low.x);
 
 	/*
 	 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
 	 */
 	PG_FREE_IF_COPY(polya, 0);
 	PG_FREE_IF_COPY(polyb, 1);
 
 	PG_RETURN_BOOL(result);
 }
 
@@ -3501,21 +3518,21 @@ poly_left(PG_FUNCTION_ARGS)
  * the right most point of A at or left of the right most point
  * of B?
  *-------------------------------------------------------*/
 Datum
 poly_overleft(PG_FUNCTION_ARGS)
 {
 	POLYGON    *polya = PG_GETARG_POLYGON_P(0);
 	POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
 	bool		result;
 
-	result = polya->boundbox.high.x <= polyb->boundbox.high.x;
+	result = float8_le(polya->boundbox.high.x, polyb->boundbox.high.x);
 
 	/*
 	 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
 	 */
 	PG_FREE_IF_COPY(polya, 0);
 	PG_FREE_IF_COPY(polyb, 1);
 
 	PG_RETURN_BOOL(result);
 }
 
@@ -3524,21 +3541,21 @@ poly_overleft(PG_FUNCTION_ARGS)
  * the left most point of A right of the right most point
  * of B?
  *-------------------------------------------------------*/
 Datum
 poly_right(PG_FUNCTION_ARGS)
 {
 	POLYGON    *polya = PG_GETARG_POLYGON_P(0);
 	POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
 	bool		result;
 
-	result = polya->boundbox.low.x > polyb->boundbox.high.x;
+	result = float8_gt(polya->boundbox.low.x, polyb->boundbox.high.x);
 
 	/*
 	 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
 	 */
 	PG_FREE_IF_COPY(polya, 0);
 	PG_FREE_IF_COPY(polyb, 1);
 
 	PG_RETURN_BOOL(result);
 }
 
@@ -3547,21 +3564,21 @@ poly_right(PG_FUNCTION_ARGS)
  * the left most point of A at or right of the left most point
  * of B?
  *-------------------------------------------------------*/
 Datum
 poly_overright(PG_FUNCTION_ARGS)
 {
 	POLYGON    *polya = PG_GETARG_POLYGON_P(0);
 	POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
 	bool		result;
 
-	result = polya->boundbox.low.x >= polyb->boundbox.low.x;
+	result = float8_ge(polya->boundbox.low.x, polyb->boundbox.low.x);
 
 	/*
 	 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
 	 */
 	PG_FREE_IF_COPY(polya, 0);
 	PG_FREE_IF_COPY(polyb, 1);
 
 	PG_RETURN_BOOL(result);
 }
 
@@ -3570,21 +3587,21 @@ poly_overright(PG_FUNCTION_ARGS)
  * the upper most point of A below the lower most point
  * of B?
  *-------------------------------------------------------*/
 Datum
 poly_below(PG_FUNCTION_ARGS)
 {
 	POLYGON    *polya = PG_GETARG_POLYGON_P(0);
 	POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
 	bool		result;
 
-	result = polya->boundbox.high.y < polyb->boundbox.low.y;
+	result = float8_lt(polya->boundbox.high.y, polyb->boundbox.low.y);
 
 	/*
 	 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
 	 */
 	PG_FREE_IF_COPY(polya, 0);
 	PG_FREE_IF_COPY(polyb, 1);
 
 	PG_RETURN_BOOL(result);
 }
 
@@ -3593,21 +3610,21 @@ poly_below(PG_FUNCTION_ARGS)
  * the upper most point of A at or below the upper most point
  * of B?
  *-------------------------------------------------------*/
 Datum
 poly_overbelow(PG_FUNCTION_ARGS)
 {
 	POLYGON    *polya = PG_GETARG_POLYGON_P(0);
 	POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
 	bool		result;
 
-	result = polya->boundbox.high.y <= polyb->boundbox.high.y;
+	result = float8_le(polya->boundbox.high.y, polyb->boundbox.high.y);
 
 	/*
 	 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
 	 */
 	PG_FREE_IF_COPY(polya, 0);
 	PG_FREE_IF_COPY(polyb, 1);
 
 	PG_RETURN_BOOL(result);
 }
 
@@ -3616,21 +3633,21 @@ poly_overbelow(PG_FUNCTION_ARGS)
  * the lower most point of A above the upper most point
  * of B?
  *-------------------------------------------------------*/
 Datum
 poly_above(PG_FUNCTION_ARGS)
 {
 	POLYGON    *polya = PG_GETARG_POLYGON_P(0);
 	POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
 	bool		result;
 
-	result = polya->boundbox.low.y > polyb->boundbox.high.y;
+	result = float8_gt(polya->boundbox.low.y, polyb->boundbox.high.y);
 
 	/*
 	 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
 	 */
 	PG_FREE_IF_COPY(polya, 0);
 	PG_FREE_IF_COPY(polyb, 1);
 
 	PG_RETURN_BOOL(result);
 }
 
@@ -3639,21 +3656,21 @@ poly_above(PG_FUNCTION_ARGS)
  * the lower most point of A at or above the lower most point
  * of B?
  *-------------------------------------------------------*/
 Datum
 poly_overabove(PG_FUNCTION_ARGS)
 {
 	POLYGON    *polya = PG_GETARG_POLYGON_P(0);
 	POLYGON    *polyb = PG_GETARG_POLYGON_P(1);
 	bool		result;
 
-	result = polya->boundbox.low.y >= polyb->boundbox.low.y;
+	result = float8_ge(polya->boundbox.low.y, polyb->boundbox.low.y);
 
 	/*
 	 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
 	 */
 	PG_FREE_IF_COPY(polya, 0);
 	PG_FREE_IF_COPY(polyb, 1);
 
 	PG_RETURN_BOOL(result);
 }
 
@@ -3850,22 +3867,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.
  *-----------------------------------------------------------------*/
@@ -3992,22 +4009,22 @@ point_add(PG_FUNCTION_ARGS)
 
 	point_add_internal(result, p1, p2);
 
 	PG_RETURN_POINT_P(result);
 }
 
 static inline void
 point_add_internal(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_sub(PG_FUNCTION_ARGS)
 {
 	Point	   *p1 = PG_GETARG_POINT_P(0);
 	Point	   *p2 = PG_GETARG_POINT_P(1);
 	Point	   *result;
 
@@ -4015,22 +4032,22 @@ point_sub(PG_FUNCTION_ARGS)
 
 	point_sub_internal(result, p1, p2);
 
 	PG_RETURN_POINT_P(result);
 }
 
 static inline void
 point_sub_internal(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_mul(PG_FUNCTION_ARGS)
 {
 	Point	   *p1 = PG_GETARG_POINT_P(0);
 	Point	   *p2 = PG_GETARG_POINT_P(1);
 	Point	   *result;
 
@@ -4038,22 +4055,24 @@ point_mul(PG_FUNCTION_ARGS)
 
 	point_mul_internal(result, p1, p2);
 
 	PG_RETURN_POINT_P(result);
 }
 
 static inline void
 point_mul_internal(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_div(PG_FUNCTION_ARGS)
 {
 	Point	   *p1 = PG_GETARG_POINT_P(0);
 	Point	   *p2 = PG_GETARG_POINT_P(1);
 	Point	   *result;
 
@@ -4062,24 +4081,26 @@ point_div(PG_FUNCTION_ARGS)
 	point_div_internal(result, p1, p2);
 
 	PG_RETURN_POINT_P(result);
 }
 
 static inline void
 point_div_internal(Point *result, Point *pt1, Point *pt2)
 {
 	float8		div;
 
-	div = (pt2->x * pt2->x) + (pt2->y * pt2->y);
+	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));
 }
 
 
 /***********************************************************************
  **
  **		Routines for 2D boxes.
  **
  ***********************************************************************/
 
 Datum
@@ -4188,24 +4209,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.
  **
  ***********************************************************************/
@@ -4514,21 +4535,21 @@ 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)
+	if (float8_lt(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--;
@@ -4581,21 +4602,21 @@ 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)
+	if (float8_lt(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
  */
@@ -4632,144 +4653,146 @@ circle_same(PG_FUNCTION_ARGS)
 
 /*		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) + circle1->radius), circle2->radius));
+	PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
+						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) + circle2->radius), circle1->radius));
+	PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
+						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);
@@ -4868,36 +4891,38 @@ 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_internal(&result->center, &circle->center, point);
-	result->radius *= hypot(point->x, point->y);
+	result->radius = float8_mul(circle->radius,
+								float8_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_internal(&result->center, &circle->center, point);
-	result->radius /= hypot(point->x, point->y);
+	result->radius = float8_div(circle->radius,
+								float8_hypot(point->x, point->y));
 
 	PG_RETURN_CIRCLE_P(result);
 }
 
 
 /*		circle_area		-		returns the area of the circle.
  */
 Datum
 circle_area(PG_FUNCTION_ARGS)
 {
@@ -4907,21 +4932,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);
 
@@ -4932,81 +4957,84 @@ 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 (float8_lt(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);
+	PG_RETURN_BOOL(float8_le(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);
+	PG_RETURN_BOOL(float8_le(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 (float8_lt(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 (float8_lt(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;
@@ -5014,24 +5042,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)
 {
@@ -5046,65 +5074,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")));
@@ -5115,27 +5143,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);
 }
 
 /*		poly_circle		- convert polygon to circle
  *
  * XXX This algorithm should use weighted means of line segments
@@ -5154,29 +5185,30 @@ poly_circle(PG_FUNCTION_ARGS)
 				 errmsg("cannot convert empty polygon to circle")));
 
 	circle = (CIRCLE *) palloc(sizeof(CIRCLE));
 
 	circle->center.x = 0;
 	circle->center.y = 0;
 	circle->radius = 0;
 
 	for (i = 0; i < poly->npts; i++)
 	{
-		circle->center.x += poly->p[i].x;
-		circle->center.y += poly->p[i].y;
+		circle->center.x = float8_pl(circle->center.x, poly->p[i].x);
+		circle->center.y = float8_pl(circle->center.y, poly->p[i].y);
 	}
-	circle->center.x /= poly->npts;
-	circle->center.y /= poly->npts;
+	circle->center.x = float8_div(circle->center.x, poly->npts);
+	circle->center.y = float8_div(circle->center.y, poly->npts);
 
 	for (i = 0; i < poly->npts; i++)
-		circle->radius += point_dt(&poly->p[i], &circle->center);
-	circle->radius /= poly->npts;
+		circle->radius = float8_pl(circle->radius, point_dt(&poly->p[i],
+															&circle->center));
+	circle->radius = float8_div(circle->radius, poly->npts);
 
 	PG_RETURN_CIRCLE_P(circle);
 }
 
 
 /***********************************************************************
  **
  **		Private routines for multiple types.
  **
  ***********************************************************************/
@@ -5190,45 +5222,45 @@ 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;
 
 	if (npts <= 0)
 		return 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;
 	}
 
@@ -5246,69 +5278,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,
diff --git a/src/backend/utils/adt/geo_spgist.c b/src/backend/utils/adt/geo_spgist.c
index c800bb1338..f62be1061e 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
 {
@@ -121,30 +121,30 @@ typedef struct
  * The quadrant is 8 bit unsigned integer with 4 least bits in use.
  * This function accepts BOXes as input.  They are not casted to
  * RangeBoxes, yet.  All 4 bits are set by comparing a corner of the box.
  * This makes 16 quadrants in total.
  */
 static uint8
 getQuadrant(BOX *centroid, BOX *inBox)
 {
 	uint8		quadrant = 0;
 
-	if (inBox->low.x > centroid->low.x)
+	if (float8_gt(inBox->low.x, centroid->low.x))
 		quadrant |= 0x8;
 
-	if (inBox->high.x > centroid->high.x)
+	if (float8_gt(inBox->high.x, centroid->high.x))
 		quadrant |= 0x4;
 
-	if (inBox->low.y > centroid->low.y)
+	if (float8_gt(inBox->low.y, centroid->low.y))
 		quadrant |= 0x2;
 
-	if (inBox->high.y > centroid->high.y)
+	if (float8_gt(inBox->high.y, centroid->high.y))
 		quadrant |= 0x1;
 
 	return quadrant;
 }
 
 /*
  * Get RangeBox using BOX
  *
  * We are turning the BOX to our structures to emphasize their function
  * of representing points in 4D space.  It also is more convenient to
@@ -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 22c68f5b62..eee0ceff11 100644
--- a/src/include/utils/geo_decls.h
+++ b/src/include/utils/geo_decls.h
@@ -1,23 +1,20 @@
 /*-------------------------------------------------------------------------
  *
  * geo_decls.h - Declarations for various 2D constructs.
  *
  *
  * Portions Copyright (c) 1996-2017, 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
 
@@ -28,43 +25,43 @@
 
 /*--------------------------------------------------------------------
  * Useful floating point utilities and constants.
  *-------------------------------------------------------------------*/
 
 
 #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)
-#define FPge(A,B)				((B) - (A) <= EPSILON)
+#define FPeq(A, B)				(float8_le(fabs(float8_mi(A, B)), EPSILON))
+#define FPne(A, B)				(float8_gt(fabs(float8_mi(A, B)), EPSILON))
+#define FPlt(A, B)				(float8_gt(float8_mi(B, A), EPSILON))
+#define FPle(A, B)				(float8_le(float8_mi(A, B), EPSILON))
+#define FPgt(A, B)				(float8_gt(float8_mi(A, B), EPSILON))
+#define FPge(A, B)				(float8_le(float8_mi(B, A), EPSILON))
 #else
-#define FPzero(A)				((A) == 0)
-#define FPeq(A,B)				((A) == (B))
-#define FPne(A,B)				((A) != (B))
-#define FPlt(A,B)				((A) < (B))
-#define FPle(A,B)				((A) <= (B))
-#define FPgt(A,B)				((A) > (B))
-#define FPge(A,B)				((A) >= (B))
+#define FPzero(A)				((A) == 0.0)
+#define FPeq(A, B)				(float8_eq(A, B))
+#define FPne(A, B)				(float8_ne(A, B))
+#define FPlt(A, B)				(float8_lt(A, B))
+#define FPle(A, B)				(float8_le(A, B))
+#define FPgt(A, B)				(float8_gt(A, B))
+#define FPge(A, B)				(float8_ge(A, B))
 #endif
 
 
 /*---------------------------------------------------------------------
  * Point - (x,y)
  *-------------------------------------------------------------------*/
 typedef struct
 {
-	double		x,
+	float8		x,
 				y;
 } Point;
 
 
 /*---------------------------------------------------------------------
  * LSEG - A straight line, specified by endpoints.
  *-------------------------------------------------------------------*/
 typedef struct
 {
 	Point		p[2];
@@ -72,66 +69,66 @@ typedef struct
 
 
 /*---------------------------------------------------------------------
  * PATH - Specified by vertex points.
  *-------------------------------------------------------------------*/
 typedef struct
 {
 	int32		vl_len_;		/* varlena header (do not touch directly!) */
 	int32		npts;
 	int32		closed;			/* is this a closed polygon? */
-	int32		dummy;			/* padding to make it double align */
+	int32		dummy;			/* padding to make it float8 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
 {
 	Point		high,
 				low;			/* corner POINTs */
 } BOX;
 
 /*---------------------------------------------------------------------
- * POLYGON - Specified by an array of doubles defining the points,
+ * POLYGON - Specified by an array of float8s defining the points,
  *		keeping the number of points and the bounding box for
  *		speed purposes.
  *-------------------------------------------------------------------*/
 typedef struct
 {
 	int32		vl_len_;		/* varlena header (do not touch directly!) */
 	int32		npts;
 	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))
@@ -166,11 +163,11 @@ typedef struct
 #define PolygonPGetDatum(X)			PointerGetDatum(X)
 #define PG_GETARG_POLYGON_P(n)		DatumGetPolygonP(PG_GETARG_DATUM(n))
 #define PG_GETARG_POLYGON_P_COPY(n) DatumGetPolygonPCopy(PG_GETARG_DATUM(n))
 #define PG_RETURN_POLYGON_P(x)		return PolygonPGetDatum(x)
 
 #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)
 
-#endif   /* GEO_DECLS_H */
+#endif							/* GEO_DECLS_H */
-- 
2.13.5 (Apple Git-94)

