From 20c1d6645eff82fc54cfb905795044f59ab7916c Mon Sep 17 00:00:00 2001
From: Tomas Vondra <tomas@2ndquadrant.com>
Date: Thu, 26 Jul 2018 15:14:36 +0200
Subject: [PATCH 3/6] Use the built-in float datatype to implement geometric
 types

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

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

All changes can be summarised as:

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

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

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

Parallel discussions: https://www.postgresql.org/message-id/28685.1468246504%40sss.pgh.pa.us
---
 src/backend/access/gist/gistproc.c |  90 +++----
 src/backend/utils/adt/geo_ops.c    | 501 ++++++++++++++++++++-----------------
 src/backend/utils/adt/geo_spgist.c |  28 +--
 src/include/utils/geo_decls.h      |  19 +-
 4 files changed, 341 insertions(+), 297 deletions(-)

diff --git a/src/backend/access/gist/gistproc.c b/src/backend/access/gist/gistproc.c
index d6ce5ccf6c..895cf2d473 100644
--- a/src/backend/access/gist/gistproc.c
+++ b/src/backend/access/gist/gistproc.c
@@ -55,7 +55,7 @@ rt_box_union(BOX *n, const BOX *a, const BOX *b)
  * 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)
 {
 	/*
@@ -76,20 +76,21 @@ size_box(const BOX *box)
 	 */
 	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));
 }
 
 /*
@@ -263,7 +264,7 @@ 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;
 
 /*
@@ -279,13 +280,13 @@ typedef struct
 
 	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;
 
@@ -294,7 +295,7 @@ typedef struct
  */
 typedef struct
 {
-	double		lower,
+	float8		lower,
 				upper;
 } SplitInterval;
 
@@ -304,7 +305,7 @@ typedef struct
 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);
@@ -316,7 +317,7 @@ interval_cmp_lower(const void *i1, const void *i2)
 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);
@@ -339,14 +340,14 @@ non_negative(float val)
  */
 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
@@ -369,8 +370,7 @@ g_box_consider_split(ConsiderSplitContext *context, int dimNum,
 	 * 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)
 	{
@@ -384,11 +384,13 @@ g_box_consider_split(ConsiderSplitContext *context, int dimNum,
 		 * 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)
@@ -531,7 +533,7 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
 	context.first = true;		/* nothing selected yet */
 	for (dim = 0; dim < 2; dim++)
 	{
-		double		leftUpper,
+		float8		leftUpper,
 					rightLower;
 		int			i1,
 					i2;
@@ -728,7 +730,7 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
 	 */
 	for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i))
 	{
-		double		lower,
+		float8		lower,
 					upper;
 
 		/*
@@ -783,7 +785,7 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
 		 * 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
@@ -792,8 +794,8 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
 		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)));
 		}
 
 		/*
@@ -1107,10 +1109,10 @@ gist_circle_compress(PG_FUNCTION_ARGS)
 		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),
@@ -1148,10 +1150,10 @@ gist_circle_consistent(PG_FUNCTION_ARGS)
 	 * 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);
@@ -1216,10 +1218,10 @@ gist_point_fetch(PG_FUNCTION_ARGS)
 	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)
 	{
@@ -1237,9 +1239,9 @@ computeDistance(bool isLeaf, BOX *box, Point *point)
 		/* point is over or below box */
 		Assert(box->low.y <= box->high.y);
 		if (point->y > box->high.y)
-			result = point->y - box->high.y;
+			result = float8_mi(point->y, box->high.y);
 		else if (point->y < box->low.y)
-			result = box->low.y - point->y;
+			result = float8_mi(box->low.y, point->y);
 		else
 			elog(ERROR, "inconsistent point values");
 	}
@@ -1248,9 +1250,9 @@ computeDistance(bool isLeaf, BOX *box, Point *point)
 		/* point is to left or right of box */
 		Assert(box->low.x <= box->high.x);
 		if (point->x > box->high.x)
-			result = point->x - box->high.x;
+			result = float8_mi(point->x, box->high.x);
 		else if (point->x < box->low.x)
-			result = box->low.x - point->x;
+			result = float8_mi(box->low.x, point->x);
 		else
 			elog(ERROR, "inconsistent point values");
 	}
@@ -1258,7 +1260,7 @@ computeDistance(bool isLeaf, BOX *box, Point *point)
 	{
 		/* closest point will be a vertex */
 		Point		p;
-		double		subresult;
+		float8		subresult;
 
 		result = point_point_distance(point, &box->low);
 
@@ -1449,7 +1451,7 @@ 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)
@@ -1478,11 +1480,11 @@ gist_point_distance(PG_FUNCTION_ARGS)
  * 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. */
@@ -1512,7 +1514,7 @@ gist_circle_distance(PG_FUNCTION_ARGS)
 
 	/* 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);
 
@@ -1528,7 +1530,7 @@ gist_poly_distance(PG_FUNCTION_ARGS)
 
 	/* 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);
 
diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index 1fb2ff2603..314e442f2e 100644
--- a/src/backend/utils/adt/geo_ops.c
+++ b/src/backend/utils/adt/geo_ops.c
@@ -59,7 +59,7 @@ static inline float8 lseg_sl(LSEG *lseg);
 static inline float8 lseg_invsl(LSEG *lseg);
 static bool	lseg_interpt_line(Point *result, LSEG *lseg, LINE *line);
 static bool	lseg_interpt_lseg(Point *result, LSEG *l1, LSEG *l2);
-static int	lseg_crossing(double x, double y, double px, double py);
+static int	lseg_crossing(float8 x, float8 y, float8 px, float8 py);
 static bool	lseg_contain_point(LSEG *lseg, Point *point);
 static float8 lseg_closept_point(Point *result, LSEG *lseg, Point *pt);
 static float8 lseg_closept_line(Point *result, LSEG *lseg, LINE *line);
@@ -69,9 +69,9 @@ static float8 lseg_closept_lseg(Point *result, LSEG *l1, LSEG *l2);
 static inline void box_construct(BOX *result, Point *pt1, Point *pt2);
 static void box_cn(Point *center, BOX *box);
 static bool box_ov(BOX *box1, BOX *box2);
-static double box_ar(BOX *box);
-static double box_ht(BOX *box);
-static double box_wd(BOX *box);
+static float8 box_ar(BOX *box);
+static float8 box_ht(BOX *box);
+static float8 box_wd(BOX *box);
 static bool box_contain_point(BOX *box, Point *point);
 static bool box_contain_box(BOX *box1, BOX *box2);
 static bool box_contain_lseg(BOX *box, LSEG *lseg);
@@ -80,7 +80,7 @@ static float8 box_closept_point(Point *result, BOX *box, Point *point);
 static float8 box_closept_lseg(Point *result, BOX *box, LSEG *lseg);
 
 /* Routines for circles */
-static double circle_ar(CIRCLE *circle);
+static float8 circle_ar(CIRCLE *circle);
 
 /* Routines for polygons */
 static void make_bound_box(POLYGON *poly);
@@ -91,10 +91,10 @@ static bool plist_same(int npts, Point *p1, Point *p2);
 static float8 dist_ppoly_internal(Point *pt, POLYGON *poly);
 
 /* Routines for encoding and decoding */
-static double single_decode(char *num, char **endptr_p,
+static float8 single_decode(char *num, char **endptr_p,
 			  const char *type_name, const char *orig_string);
 static void single_encode(float8 x, StringInfo str);
-static void pair_decode(char *str, double *x, double *y, char **endptr_p,
+static void pair_decode(char *str, float8 *x, float8 *y, char **endptr_p,
 			const char *type_name, const char *orig_string);
 static void pair_encode(float8 x, float8 y, StringInfo str);
 static int	pair_count(char *s, char delim);
@@ -146,7 +146,7 @@ static char *path_encode(enum path_delim path_delim, int npts, Point *pt);
  *	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)
 {
@@ -163,7 +163,7 @@ single_encode(float8 x, StringInfo str)
 }								/* 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;
@@ -376,19 +376,19 @@ 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;
@@ -416,7 +416,7 @@ box_recv(PG_FUNCTION_ARGS)
 {
 	StringInfo	buf = (StringInfo) PG_GETARG_POINTER(0);
 	BOX		   *box;
-	double		x,
+	float8		x,
 				y;
 
 	box = (BOX *) palloc(sizeof(BOX));
@@ -427,13 +427,13 @@ box_recv(PG_FUNCTION_ARGS)
 	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;
@@ -466,7 +466,7 @@ box_send(PG_FUNCTION_ARGS)
 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;
@@ -476,7 +476,7 @@ box_construct(BOX *result, Point *pt1, Point *pt2)
 		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;
@@ -808,10 +808,10 @@ box_center(PG_FUNCTION_ARGS)
 
 /*		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));
 }
 
 
@@ -820,28 +820,28 @@ box_ar(BOX *box)
 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);
 }
 
 
@@ -865,10 +865,10 @@ box_intersect(PG_FUNCTION_ARGS)
 
 	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);
 }
@@ -1014,8 +1014,8 @@ line_construct(LINE *result, Point *pt, float8 m)
 	if (m == DBL_MAX)
 	{
 		/* vertical - use "x = C" */
-		result->A = -1;
-		result->B = 0;
+		result->A = -1.0;
+		result->B = 0.0;
 		result->C = pt->x;
 	}
 	else
@@ -1023,7 +1023,7 @@ line_construct(LINE *result, Point *pt, float8 m)
 		/* 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));
 	}
 }
 
@@ -1076,7 +1076,8 @@ line_perp(PG_FUNCTION_ARGS)
 	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
@@ -1095,25 +1096,35 @@ line_horizontal(PG_FUNCTION_ARGS)
 	PG_RETURN_BOOL(FPzero(line->A));
 }
 
+
+/*
+ * Check whether the two lines are the same
+ *
+ * We consider NaNs values to be equal to each other to let those lines
+ * to be found.
+ */
 Datum
 line_eq(PG_FUNCTION_ARGS)
 {
 	LINE	   *l1 = PG_GETARG_LINE_P(0);
 	LINE	   *l2 = PG_GETARG_LINE_P(1);
-	double		k;
-
-	if (!FPzero(l2->A))
-		k = l1->A / l2->A;
-	else if (!FPzero(l2->B))
-		k = l1->B / l2->B;
-	else if (!FPzero(l2->C))
-		k = l1->C / l2->C;
+	float8		ratio;
+
+	if (!FPzero(l2->A) && !isnan(l2->A))
+		ratio = float8_div(l1->A, l2->A);
+	else if (!FPzero(l2->B) && !isnan(l2->B))
+		ratio = float8_div(l1->B, l2->B);
+	else if (!FPzero(l2->C) && !isnan(l2->C))
+		ratio = float8_div(l1->C, l2->C);
 	else
-		k = 1.0;
+		ratio = 1.0;
 
-	PG_RETURN_BOOL(FPeq(l1->A, k * l2->A) &&
-				   FPeq(l1->B, k * l2->B) &&
-				   FPeq(l1->C, k * l2->C));
+	PG_RETURN_BOOL((FPeq(l1->A, float8_mul(ratio, l2->A)) &&
+					FPeq(l1->B, float8_mul(ratio, l2->B)) &&
+					FPeq(l1->C, float8_mul(ratio, l2->C))) ||
+				   (float8_eq(l1->A, l2->A) &&
+					float8_eq(l1->B, l2->B) &&
+					float8_eq(l1->C, l2->C)));
 }
 
 
@@ -1131,7 +1142,7 @@ line_invsl(LINE *line)
 		return DBL_MAX;
 	if (FPzero(line->B))
 		return 0.0;
-	return line->B / line->A;
+	return float8_div(line->B, line->A);
 }
 
 
@@ -1149,7 +1160,7 @@ line_distance(PG_FUNCTION_ARGS)
 	if (line_interpt_line(NULL, l1, l2))	/* intersecting? */
 		PG_RETURN_FLOAT8(0.0);
 	if (FPzero(l1->B))			/* vertical? */
-		PG_RETURN_FLOAT8(fabs(l1->C - l2->C));
+		PG_RETURN_FLOAT8(fabs(float8_mi(l1->C, l2->C)));
 	point_construct(&tmp, 0.0, l1->C);
 	result = line_closept_point(NULL, l2, &tmp);
 	PG_RETURN_FLOAT8(result);
@@ -1182,11 +1193,15 @@ line_interpt(PG_FUNCTION_ARGS)
  * NOTE: If the lines are identical then we will find they are parallel
  * and report "no intersection".  This is a little weird, but since
  * there's no *unique* intersection, maybe it's appropriate behavior.
+ *
+ * If the lines have NaN constants, we will return true, and the intersection
+ * point would have NaN coordinates.  We shouldn't return false in this case
+ * because that would mean the lines are parallel.
  */
 static bool
 line_interpt_line(Point *result, LINE *l1, LINE *l2)
 {
-	double		x,
+	float8		x,
 				y;
 
 	if (FPzero(l1->B))			/* l1 vertical? */
@@ -1195,20 +1210,20 @@ line_interpt_line(Point *result, LINE *l1, LINE *l2)
 			return false;
 
 		x = l1->C;
-		y = (l2->A * x + l2->C);
+		y = float8_pl(float8_mul(l2->A, x), l2->C);
 	}
 	else if (FPzero(l2->B))		/* l2 vertical? */
 	{
 		x = l2->C;
-		y = (l1->A * x + l1->C);
+		y = float8_pl(float8_mul(l1->A, x), l1->C);
 	}
 	else
 	{
-		if (FPeq(l2->A, l1->A * (l2->B / l1->B)))
+		if (FPeq(l2->A, float8_mul(l1->A, float8_div(l2->B, l1->B))))
 			return false;
 
-		x = (l1->C - l2->C) / (l2->A - l1->A);
-		y = (l1->A * x + l1->C);
+		x = float8_div(float8_mi(l1->C, l2->C), float8_mi(l2->A, l1->A));
+		y = float8_pl(float8_mul(l1->A, x), l1->C);
 	}
 
 	if (result != NULL)
@@ -1244,7 +1259,7 @@ Datum
 path_area(PG_FUNCTION_ARGS)
 {
 	PATH	   *path = PG_GETARG_PATH_P(0);
-	double		area = 0.0;
+	float8		area = 0.0;
 	int			i,
 				j;
 
@@ -1254,12 +1269,11 @@ path_area(PG_FUNCTION_ARGS)
 	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));
 }
 
 
@@ -1529,19 +1543,19 @@ path_inter(PG_FUNCTION_ARGS)
 	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);
@@ -1631,7 +1645,7 @@ path_distance(PG_FUNCTION_ARGS)
 			statlseg_construct(&seg2, &p2->p[jprev], &p2->p[j]);
 
 			tmp = lseg_closept_lseg(NULL, &seg1, &seg2);
-			if (!have_min || tmp < min)
+			if (!have_min || float8_lt(tmp, min))
 			{
 				min = tmp;
 				have_min = true;
@@ -1670,7 +1684,7 @@ path_length(PG_FUNCTION_ARGS)
 			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);
@@ -1830,10 +1844,18 @@ point_ne(PG_FUNCTION_ARGS)
 	PG_RETURN_BOOL(!point_eq_point(pt1, pt2));
 }
 
+
+/*
+ * Check whether the two points are the same
+ *
+ * We consider NaNs coordinates to be equal to each other to let those points
+ * to be found.
+ */
 static inline bool
 point_eq_point(Point *pt1, Point *pt2)
 {
-	return FPeq(pt1->x, pt2->x) && FPeq(pt1->y, pt2->y);
+	return ((FPeq(pt1->x, pt2->x) && FPeq(pt1->y, pt2->y)) ||
+			(float8_eq(pt1->x, pt2->x) && float8_eq(pt1->y, pt2->y)));
 }
 
 
@@ -1853,7 +1875,7 @@ point_distance(PG_FUNCTION_ARGS)
 static inline float8
 point_dt(Point *pt1, Point *pt2)
 {
-	return HYPOT(pt1->x - pt2->x, pt1->y - pt2->y);
+	return HYPOT(float8_mi(pt1->x, pt2->x), float8_mi(pt1->y, pt2->y));
 }
 
 Datum
@@ -1878,7 +1900,7 @@ point_sl(Point *pt1, Point *pt2)
 		return DBL_MAX;
 	if (FPeq(pt1->y, pt2->y))
 		return 0.0;
-	return (pt1->y - pt2->y) / (pt1->x - pt2->x);
+	return float8_div(float8_mi(pt1->y, pt2->y), float8_mi(pt1->x, pt2->x));
 }
 
 
@@ -1894,7 +1916,7 @@ point_invsl(Point *pt1, Point *pt2)
 		return 0.0;
 	if (FPeq(pt1->y, pt2->y))
 		return DBL_MAX;
-	return (pt1->x - pt2->x) / (pt2->y - pt1->y);
+	return float8_div(float8_mi(pt1->x, pt2->x), float8_mi(pt2->y, pt1->y));
 }
 
 
@@ -2168,8 +2190,8 @@ lseg_center(PG_FUNCTION_ARGS)
 
 	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);
 }
@@ -2292,7 +2314,7 @@ dist_ppath(PG_FUNCTION_ARGS)
 
 		statlseg_construct(&lseg, &path->p[iprev], &path->p[i]);
 		tmp = lseg_closept_point(NULL, &lseg, pt);
-		if (!have_min || tmp < result)
+		if (!have_min || float8_lt(tmp, result))
 		{
 			result = tmp;
 			have_min = true;
@@ -2322,19 +2344,14 @@ dist_sl(PG_FUNCTION_ARGS)
 {
 	LSEG	   *lseg = PG_GETARG_LSEG_P(0);
 	LINE	   *line = PG_GETARG_LINE_P(1);
-	float8		result,
-				d2;
+	float8		result;
 
 	if (lseg_interpt_line(NULL, lseg, line))
 		result = 0.0;
 	else
-	{
-		result = line_closept_point(NULL, line, &lseg->p[0]);
-		d2 = line_closept_point(NULL, line, &lseg->p[1]);
 		/* XXX shouldn't we take the min not max? */
-		if (d2 > result)
-			result = d2;
-	}
+		result = float8_max(line_closept_point(NULL, line, &lseg->p[0]),
+							line_closept_point(NULL, line, &lseg->p[1]));
 
 	PG_RETURN_FLOAT8(result);
 }
@@ -2381,11 +2398,10 @@ dist_cpoly(PG_FUNCTION_ARGS)
 	float8		result;
 
 	/* calculate distance to center, and subtract radius */
-	result = dist_ppoly_internal(&circle->center, poly);
-
-	result -= circle->radius;
-	if (result < 0)
-		result = 0;
+	result = float8_mi(dist_ppoly_internal(&circle->center, poly),
+					   circle->radius);
+	if (result < 0.0)
+		result = 0.0;
 
 	PG_RETURN_FLOAT8(result);
 }
@@ -2411,7 +2427,7 @@ dist_polyp(PG_FUNCTION_ARGS)
 	PG_RETURN_FLOAT8(dist_ppoly_internal(point, poly));
 }
 
-static double
+static float8
 dist_ppoly_internal(Point *pt, POLYGON *poly)
 {
 	float8		result;
@@ -2448,7 +2464,7 @@ dist_ppoly_internal(Point *pt, POLYGON *poly)
 #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;
 	}
 
@@ -2550,7 +2566,8 @@ close_pl(PG_FUNCTION_ARGS)
 
 	result = (Point *) palloc(sizeof(Point));
 
-	line_closept_point(result, line, pt);
+	if (isnan(line_closept_point(result, line, pt)))
+		PG_RETURN_NULL();
 
 	PG_RETURN_POINT_P(result);
 }
@@ -2607,8 +2624,8 @@ static float8
 lseg_closept_lseg(Point *result, LSEG *l1, LSEG *l2)
 {
 	Point		point;
-	double		dist;
-	double		d;
+	float8		dist,
+				d;
 
 	d = lseg_closept_point(NULL, l1, &l2->p[0]);
 	dist = d;
@@ -2616,20 +2633,20 @@ lseg_closept_lseg(Point *result, LSEG *l1, LSEG *l2)
 		*result = l2->p[0];
 
 	d = lseg_closept_point(NULL, l1, &l2->p[1]);
-	if (d < dist)
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		if (result != NULL)
 			*result = l2->p[1];
 	}
 
-	if (lseg_closept_point(&point, l2, &l1->p[0]) < dist)
+	if (float8_lt(lseg_closept_point(&point, l2, &l1->p[0]), dist))
 		d = lseg_closept_point(result, l1, &point);
 
-	if (lseg_closept_point(&point, l2, &l1->p[1]) < dist)
+	if (float8_lt(lseg_closept_point(&point, l2, &l1->p[1]), dist))
 		d = lseg_closept_point(result, l1, &point);
 
-	if (d < dist)
+	if (float8_lt(d, dist))
 		dist = d;
 
 	return dist;
@@ -2644,7 +2661,8 @@ close_lseg(PG_FUNCTION_ARGS)
 
 	result = (Point *) palloc(sizeof(Point));
 
-	lseg_closept_lseg(result, l2, l1);
+	if (isnan(lseg_closept_lseg(result, l2, l1)))
+		PG_RETURN_NULL();
 
 	PG_RETURN_POINT_P(result);
 }
@@ -2659,11 +2677,11 @@ close_lseg(PG_FUNCTION_ARGS)
 static float8
 box_closept_point(Point *result, BOX *box, Point *pt)
 {
-	LSEG		lseg;
+	float8		dist,
+				d;
 	Point		point,
 				closept;
-	double		dist,
-				d;
+	LSEG		lseg;
 
 	if (box_contain_point(box, pt))
 	{
@@ -2681,7 +2699,7 @@ box_closept_point(Point *result, BOX *box, Point *pt)
 
 	statlseg_construct(&lseg, &box->high, &point);
 	d = lseg_closept_point(&closept, &lseg, pt);
-	if (d < dist)
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		if (result != NULL)
@@ -2692,7 +2710,7 @@ box_closept_point(Point *result, BOX *box, Point *pt)
 	point.y = box->low.y;
 	statlseg_construct(&lseg, &box->low, &point);
 	d = lseg_closept_point(&closept, &lseg, pt);
-	if (d < dist)
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		if (result != NULL)
@@ -2701,7 +2719,7 @@ box_closept_point(Point *result, BOX *box, Point *pt)
 
 	statlseg_construct(&lseg, &box->high, &point);
 	d = lseg_closept_point(&closept, &lseg, pt);
-	if (d < dist)
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		if (result != NULL)
@@ -2720,7 +2738,8 @@ close_pb(PG_FUNCTION_ARGS)
 
 	result = (Point *) palloc(sizeof(Point));
 
-	box_closept_point(result, box, pt);
+	if (isnan(box_closept_point(result, box, pt)))
+		PG_RETURN_NULL();
 
 	PG_RETURN_POINT_P(result);
 }
@@ -2752,7 +2771,7 @@ close_sl(PG_FUNCTION_ARGS)
 
 	d1 = line_closept_point(NULL, line, &lseg->p[0]);
 	d2 = line_closept_point(NULL, line, &lseg->p[1]);
-	if (d1 < d2)
+	if (float8_lt(d1, d2))
 		*result = lseg->p[0];
 	else
 		*result = lseg->p[1];
@@ -2816,7 +2835,8 @@ close_ls(PG_FUNCTION_ARGS)
 
 	result = (Point *) palloc(sizeof(Point));
 
-	lseg_closept_line(result, lseg, line);
+	if (isnan(lseg_closept_line(result, lseg, line)))
+		PG_RETURN_NULL();
 
 	PG_RETURN_POINT_P(result);
 }
@@ -2831,11 +2851,11 @@ close_ls(PG_FUNCTION_ARGS)
 static float8
 box_closept_lseg(Point *result, BOX *box, LSEG *lseg)
 {
+	float8		dist,
+				d;
 	Point		point,
 				closept;
 	LSEG		bseg;
-	double		dist,
-				d;
 
 	if (box_interpt_lseg(result, box, lseg))
 		return 0.0;
@@ -2848,7 +2868,7 @@ box_closept_lseg(Point *result, BOX *box, LSEG *lseg)
 
 	statlseg_construct(&bseg, &box->high, &point);
 	d = lseg_closept_lseg(&closept, &bseg, lseg);
-	if (d < dist)
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		if (result != NULL)
@@ -2859,7 +2879,7 @@ box_closept_lseg(Point *result, BOX *box, LSEG *lseg)
 	point.y = box->low.y;
 	statlseg_construct(&bseg, &box->low, &point);
 	d = lseg_closept_lseg(&closept, &bseg, lseg);
-	if (d < dist)
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		if (result != NULL)
@@ -2868,7 +2888,7 @@ box_closept_lseg(Point *result, BOX *box, LSEG *lseg)
 
 	statlseg_construct(&bseg, &box->high, &point);
 	d = lseg_closept_lseg(&closept, &bseg, lseg);
-	if (d < dist)
+	if (float8_lt(d, dist))
 	{
 		dist = d;
 		if (result != NULL)
@@ -2887,7 +2907,8 @@ close_sb(PG_FUNCTION_ARGS)
 
 	result = (Point *) palloc(sizeof(Point));
 
-	box_closept_lseg(result, box, lseg);
+	if (isnan(box_closept_lseg(result, box, lseg)))
+		PG_RETURN_NULL();
 
 	PG_RETURN_POINT_P(result);
 }
@@ -2920,7 +2941,9 @@ close_lb(PG_FUNCTION_ARGS)
 static bool
 line_contain_point(LINE *line, Point *point)
 {
-	return FPzero(line->A * point->x + line->B * point->y + line->C);
+	return FPzero(float8_pl(float8_pl(float8_mul(line->A, point->x),
+									  float8_mul(line->B, point->y)),
+							line->C));
 }
 
 Datum
@@ -3001,7 +3024,7 @@ on_ppath(PG_FUNCTION_ARGS)
 	PATH	   *path = PG_GETARG_PATH_P(1);
 	int			i,
 				n;
-	double		a,
+	float8		a,
 				b;
 
 	/*-- OPEN --*/
@@ -3012,8 +3035,7 @@ on_ppath(PG_FUNCTION_ARGS)
 		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;
 		}
@@ -3099,10 +3121,10 @@ box_interpt_lseg(Point *result, BOX *box, LSEG *lseg)
 	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))
@@ -3209,7 +3231,7 @@ static void
 make_bound_box(POLYGON *poly)
 {
 	int			i;
-	double		x1,
+	float8		x1,
 				y1,
 				x2,
 				y2;
@@ -3220,13 +3242,13 @@ make_bound_box(POLYGON *poly)
 	y2 = y1 = poly->p[0].y;
 	for (i = 1; i < poly->npts; i++)
 	{
-		if (poly->p[i].x < x1)
+		if (float8_lt(poly->p[i].x, x1))
 			x1 = poly->p[i].x;
-		if (poly->p[i].x > x2)
+		if (float8_gt(poly->p[i].x, x2))
 			x2 = poly->p[i].x;
-		if (poly->p[i].y < y1)
+		if (float8_lt(poly->p[i].y, y1))
 			y1 = poly->p[i].y;
-		if (poly->p[i].y > y2)
+		if (float8_gt(poly->p[i].y, y2))
 			y2 = poly->p[i].y;
 	}
 
@@ -3739,8 +3761,8 @@ lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start)
 		 * 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);
 	}
@@ -3880,8 +3902,8 @@ static inline void
 point_add_point(Point *result, Point *pt1, Point *pt2)
 {
 	point_construct(result,
-					pt1->x + pt2->x,
-					pt1->y + pt2->y);
+					float8_pl(pt1->x, pt2->x),
+					float8_pl(pt1->y, pt2->y));
 }
 
 Datum
@@ -3903,8 +3925,8 @@ static inline void
 point_sub_point(Point *result, Point *pt1, Point *pt2)
 {
 	point_construct(result,
-					pt1->x - pt2->x,
-					pt1->y - pt2->y);
+					float8_mi(pt1->x, pt2->x),
+					float8_mi(pt1->y, pt2->y));
 }
 
 Datum
@@ -3926,8 +3948,10 @@ static inline void
 point_mul_point(Point *result, Point *pt1, Point *pt2)
 {
 	point_construct(result,
-					(pt1->x * pt2->x) - (pt1->y * pt2->y),
-					(pt1->x * pt2->y) + (pt1->y * pt2->x));
+					float8_mi(float8_mul(pt1->x, pt2->x),
+							  float8_mul(pt1->y, pt2->y)),
+					float8_pl(float8_mul(pt1->x, pt2->y),
+							  float8_mul(pt1->y, pt2->x)));
 }
 
 Datum
@@ -3948,18 +3972,15 @@ point_mul(PG_FUNCTION_ARGS)
 static inline void
 point_div_point(Point *result, Point *pt1, Point *pt2)
 {
-	double		div;
-
-	div = (pt2->x * pt2->x) + (pt2->y * pt2->y);
+	float8		div;
 
-	if (div == 0.0)
-		ereport(ERROR,
-				(errcode(ERRCODE_DIVISION_BY_ZERO),
-				 errmsg("division by zero")));
+	div = float8_pl(float8_mul(pt2->x, pt2->x), float8_mul(pt2->y, pt2->y));
 
 	point_construct(result,
-					((pt1->x * pt2->x) + (pt1->y * pt2->y)) / div,
-					((pt2->x * pt1->y) - (pt2->y * pt1->x)) / div);
+					float8_div(float8_pl(float8_mul(pt1->x, pt2->x),
+										 float8_mul(pt1->y, pt2->y)), div),
+					float8_div(float8_mi(float8_mul(pt1->y, pt2->x),
+										 float8_mul(pt1->x, pt2->y)), div));
 }
 
 Datum
@@ -4096,10 +4117,10 @@ boxes_bound_box(PG_FUNCTION_ARGS)
 
 	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);
 }
@@ -4419,7 +4440,8 @@ circle_in(PG_FUNCTION_ARGS)
 		s++;
 
 	circle->radius = single_decode(s, &s, "circle", str);
-	if (circle->radius < 0)
+	/* We have to accept NaN. */
+	if (circle->radius < 0.0)
 		ereport(ERROR,
 				(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
 				 errmsg("invalid input syntax for type %s: \"%s\"",
@@ -4486,7 +4508,8 @@ circle_recv(PG_FUNCTION_ARGS)
 	circle->center.y = pq_getmsgfloat8(buf);
 	circle->radius = pq_getmsgfloat8(buf);
 
-	if (circle->radius < 0)
+	/* We have to accept NaN. */
+	if (circle->radius < 0.0)
 		ereport(ERROR,
 				(errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
 				 errmsg("invalid radius in external \"circle\" value")));
@@ -4517,6 +4540,9 @@ circle_send(PG_FUNCTION_ARGS)
  *---------------------------------------------------------*/
 
 /*		circles identical?
+ *
+ * We consider NaNs values to be equal to each other to let those circles
+ * to be found.
  */
 Datum
 circle_same(PG_FUNCTION_ARGS)
@@ -4524,7 +4550,8 @@ circle_same(PG_FUNCTION_ARGS)
 	CIRCLE	   *circle1 = PG_GETARG_CIRCLE_P(0);
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 
-	PG_RETURN_BOOL(FPeq(circle1->radius, circle2->radius) &&
+	PG_RETURN_BOOL(((isnan(circle1->radius) && isnan(circle1->radius)) ||
+					FPeq(circle1->radius, circle2->radius)) &&
 				   point_eq_point(&circle1->center, &circle2->center));
 }
 
@@ -4537,7 +4564,7 @@ circle_overlap(PG_FUNCTION_ARGS)
 	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
@@ -4549,8 +4576,8 @@ 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?
@@ -4561,8 +4588,8 @@ 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?
@@ -4573,8 +4600,8 @@ 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
@@ -4586,8 +4613,8 @@ 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?
@@ -4599,7 +4626,7 @@ circle_contained(PG_FUNCTION_ARGS)
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 
 	PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
-						circle2->radius - circle1->radius));
+						float8_mi(circle2->radius, circle1->radius)));
 }
 
 /*		circle_contain	-		does circle1 contain circle2?
@@ -4611,7 +4638,7 @@ circle_contain(PG_FUNCTION_ARGS)
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 
 	PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
-						circle1->radius - circle2->radius));
+						float8_mi(circle1->radius, circle2->radius)));
 }
 
 
@@ -4623,8 +4650,8 @@ 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?
@@ -4635,8 +4662,8 @@ 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
@@ -4648,8 +4675,8 @@ 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
@@ -4661,8 +4688,8 @@ 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)));
 }
 
 
@@ -4775,7 +4802,7 @@ circle_mul_pt(PG_FUNCTION_ARGS)
 	result = (CIRCLE *) palloc(sizeof(CIRCLE));
 
 	point_mul_point(&result->center, &circle->center, point);
-	result->radius = circle->radius * HYPOT(point->x, point->y);
+	result->radius = float8_mul(circle->radius, HYPOT(point->x, point->y));
 
 	PG_RETURN_CIRCLE_P(result);
 }
@@ -4790,7 +4817,7 @@ circle_div_pt(PG_FUNCTION_ARGS)
 	result = (CIRCLE *) palloc(sizeof(CIRCLE));
 
 	point_div_point(&result->center, &circle->center, point);
-	result->radius = circle->radius / HYPOT(point->x, point->y);
+	result->radius = float8_div(circle->radius, HYPOT(point->x, point->y));
 
 	PG_RETURN_CIRCLE_P(result);
 }
@@ -4814,7 +4841,7 @@ 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));
 }
 
 
@@ -4839,10 +4866,11 @@ circle_distance(PG_FUNCTION_ARGS)
 	CIRCLE	   *circle2 = PG_GETARG_CIRCLE_P(1);
 	float8		result;
 
-	result = point_dt(&circle1->center, &circle2->center) -
-		(circle1->radius + circle2->radius);
-	if (result < 0)
-		result = 0;
+	result = float8_mi(point_dt(&circle1->center, &circle2->center),
+					   float8_pl(circle1->radius, circle2->radius));
+	if (result < 0.0)
+		result = 0.0;
+
 	PG_RETURN_FLOAT8(result);
 }
 
@@ -4852,7 +4880,7 @@ 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);
@@ -4864,7 +4892,7 @@ 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);
@@ -4881,9 +4909,11 @@ dist_pc(PG_FUNCTION_ARGS)
 	CIRCLE	   *circle = PG_GETARG_CIRCLE_P(1);
 	float8		result;
 
-	result = point_dt(point, &circle->center) - circle->radius;
-	if (result < 0)
-		result = 0;
+	result = float8_mi(point_dt(point, &circle->center),
+					   circle->radius);
+	if (result < 0.0)
+		result = 0.0;
+
 	PG_RETURN_FLOAT8(result);
 }
 
@@ -4897,9 +4927,10 @@ dist_cpoint(PG_FUNCTION_ARGS)
 	Point	   *point = PG_GETARG_POINT_P(1);
 	float8		result;
 
-	result = point_dt(point, &circle->center) - circle->radius;
-	if (result < 0)
-		result = 0;
+	result = float8_mi(point_dt(point, &circle->center), circle->radius);
+	if (result < 0.0)
+		result = 0.0;
+
 	PG_RETURN_FLOAT8(result);
 }
 
@@ -4921,10 +4952,10 @@ circle_center(PG_FUNCTION_ARGS)
 
 /*		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);
 }
 
 
@@ -4953,16 +4984,16 @@ 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);
 }
@@ -4978,8 +5009,8 @@ box_circle(PG_FUNCTION_ARGS)
 
 	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);
 
@@ -4996,8 +5027,8 @@ circle_poly(PG_FUNCTION_ARGS)
 	int			base_size,
 				size;
 	int			i;
-	double		angle;
-	double		anglestep;
+	float8		angle;
+	float8		anglestep;
 
 	if (FPzero(circle->radius))
 		ereport(ERROR,
@@ -5022,13 +5053,16 @@ circle_poly(PG_FUNCTION_ARGS)
 	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);
@@ -5057,12 +5091,13 @@ poly_to_circle(CIRCLE *result, POLYGON *poly)
 
 	for (i = 0; i < poly->npts; i++)
 		point_add_point(&result->center, &result->center, &poly->p[i]);
-	result->center.x /= poly->npts;
-	result->center.y /= poly->npts;
+	result->center.x = float8_div(result->center.x, poly->npts);
+	result->center.y = float8_div(result->center.y, poly->npts);
 
 	for (i = 0; i < poly->npts; i++)
-		result->radius += point_dt(&poly->p[i], &result->center);
-	result->radius /= poly->npts;
+		result->radius = float8_pl(result->radius,
+								   point_dt(&poly->p[i], &result->center));
+	result->radius = float8_div(result->radius, poly->npts);
 }
 
 Datum
@@ -5101,12 +5136,12 @@ poly_circle(PG_FUNCTION_ARGS)
 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;
@@ -5114,8 +5149,8 @@ point_inside(Point *p, int npts, Point *plist)
 	Assert(npts > 0);
 
 	/* compute first polygon point relative to single point */
-	x0 = plist[0].x - p->x;
-	y0 = plist[0].y - p->y;
+	x0 = float8_mi(plist[0].x, p->x);
+	y0 = float8_mi(plist[0].y, p->y);
 
 	prev_x = x0;
 	prev_y = y0;
@@ -5123,8 +5158,8 @@ point_inside(Point *p, int npts, Point *plist)
 	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)
@@ -5156,9 +5191,9 @@ point_inside(Point *p, int npts, Point *plist)
  */
 
 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))
@@ -5169,42 +5204,47 @@ lseg_crossing(double x, double y, double prev_x, double prev_y)
 		{						/* 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;
 		}
 	}
 }
@@ -5288,10 +5328,11 @@ plist_same(int npts, Point *p1, Point *p2)
  * case of hypot(inf,nan) results in INF, and not NAN.
  *-----------------------------------------------------------------------
  */
-double
-pg_hypot(double x, double y)
+float8
+pg_hypot(float8 x, float8 y)
 {
-	double		yx;
+	float8		yx,
+				result;
 
 	/* Handle INF and NaN properly */
 	if (isinf(x) || isinf(y))
@@ -5307,7 +5348,7 @@ pg_hypot(double x, double y)
 	/* Swap x and y if needed to make x the larger one */
 	if (x < y)
 	{
-		double		temp = x;
+		float8		temp = x;
 
 		x = y;
 		y = temp;
@@ -5323,5 +5364,9 @@ pg_hypot(double x, double y)
 
 	/* Determine the hypotenuse */
 	yx = y / x;
-	return x * sqrt(1.0 + (yx * yx));
+	result = x * sqrt(1.0 + (yx * yx));
+
+	check_float8_val(result, false, false);
+
+	return result;
 }
diff --git a/src/backend/utils/adt/geo_spgist.c b/src/backend/utils/adt/geo_spgist.c
index fea36f361a..4aff973ef3 100644
--- a/src/backend/utils/adt/geo_spgist.c
+++ b/src/backend/utils/adt/geo_spgist.c
@@ -84,14 +84,14 @@
  * 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;
@@ -100,8 +100,8 @@ compareDoubles(const void *a, const void *b)
 
 typedef struct
 {
-	double		low;
-	double		high;
+	float8		low;
+	float8		high;
 } Range;
 
 typedef struct
@@ -175,7 +175,7 @@ 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;
@@ -418,10 +418,10 @@ spg_box_quad_picksplit(PG_FUNCTION_ARGS)
 	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++)
@@ -434,10 +434,10 @@ spg_box_quad_picksplit(PG_FUNCTION_ARGS)
 		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;
 
diff --git a/src/include/utils/geo_decls.h b/src/include/utils/geo_decls.h
index 0e066894cd..9f8505804e 100644
--- a/src/include/utils/geo_decls.h
+++ b/src/include/utils/geo_decls.h
@@ -8,9 +8,6 @@
  *
  * 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
@@ -21,14 +18,14 @@
 #ifndef GEO_DECLS_H
 #define GEO_DECLS_H
 
-#include <math.h>
-
 #include "fmgr.h"
 
 /*--------------------------------------------------------------------
  * Useful floating point utilities and constants.
- *-------------------------------------------------------------------*/
-
+ *-------------------------------------------------------------------
+ *
+ * XXX: They are not NaN-aware.
+ */
 
 #define EPSILON					1.0E-06
 
@@ -57,7 +54,7 @@
  *-------------------------------------------------------------------*/
 typedef struct
 {
-	double		x,
+	float8		x,
 				y;
 } Point;
 
@@ -89,7 +86,7 @@ typedef struct
  *-------------------------------------------------------------------*/
 typedef struct
 {
-	double		A,
+	float8		A,
 				B,
 				C;
 } LINE;
@@ -124,7 +121,7 @@ typedef struct
 typedef struct
 {
 	Point		center;
-	double		radius;
+	float8		radius;
 } CIRCLE;
 
 /*
@@ -178,6 +175,6 @@ typedef struct
  * in geo_ops.c
  */
 
-extern double pg_hypot(double x, double y);
+extern float8 pg_hypot(float8 x, float8 y);
 
 #endif							/* GEO_DECLS_H */
-- 
2.13.6

