From 6e1458bba06d5382ab14fdf49b3201f766a62a00 Mon Sep 17 00:00:00 2001
From: Emre Hasegeli <emre@hasegeli.com>
Date: Sat, 28 May 2016 18:16:05 +0200
Subject: [PATCH 2/6] float-header-v11

Provide header file for built-in float datatypes

Even though, some datatypes under adt/ have separate header files,
most of the simple ones do not.  Their public functions were on
the builtins.h.  We would need to make more functions of floats public
to let the geometric types built on top of them.  This is a good
opportunity to make a separate header file for floats.

1acf7572554515b99ef6e783750aaea8777524ec made _cmp functions public
to solve NaN problem locally for GiST indexes.  This patch reworks it
in favour of a more extensive API.  Kevin Grittner suggested to design
the API using inline functions.  They are easier to use compared
to macros, and avoid double-evaluation hazards.
---
 contrib/btree_gin/btree_gin.c                 |   1 +
 contrib/btree_gist/btree_ts.c                 |   1 +
 contrib/cube/cube.c                           |   2 +-
 contrib/cube/cubeparse.y                      |   2 +-
 contrib/postgres_fdw/postgres_fdw.c           |   1 +
 src/backend/access/gist/gistget.c             |   2 +-
 src/backend/access/gist/gistproc.c            |  55 +--
 src/backend/access/gist/gistutil.c            |   2 +-
 src/backend/utils/adt/float.c                 | 589 ++++++--------------------
 src/backend/utils/adt/formatting.c            |   1 +
 src/backend/utils/adt/geo_ops.c               |   6 +-
 src/backend/utils/adt/geo_spgist.c            |   2 +-
 src/backend/utils/adt/numeric.c               |   1 +
 src/backend/utils/adt/rangetypes_gist.c       |   3 +-
 src/backend/utils/adt/rangetypes_selfuncs.c   |   3 +-
 src/backend/utils/adt/rangetypes_typanalyze.c |   3 +-
 src/backend/utils/adt/timestamp.c             |   1 +
 src/backend/utils/misc/guc.c                  |   1 +
 src/include/utils/builtins.h                  |  14 -
 src/include/utils/float.h                     | 376 ++++++++++++++++
 src/include/utils/geo_decls.h                 |   1 +
 21 files changed, 554 insertions(+), 513 deletions(-)
 create mode 100644 src/include/utils/float.h

diff --git a/contrib/btree_gin/btree_gin.c b/contrib/btree_gin/btree_gin.c
index 2473f79ca1..301caefd47 100644
--- a/contrib/btree_gin/btree_gin.c
+++ b/contrib/btree_gin/btree_gin.c
@@ -3,20 +3,21 @@
  */
 #include "postgres.h"
 
 #include <limits.h>
 
 #include "access/stratnum.h"
 #include "utils/builtins.h"
 #include "utils/bytea.h"
 #include "utils/cash.h"
 #include "utils/date.h"
+#include "utils/float.h"
 #include "utils/inet.h"
 #include "utils/numeric.h"
 #include "utils/timestamp.h"
 #include "utils/varbit.h"
 
 PG_MODULE_MAGIC;
 
 typedef struct QueryInfo
 {
 	StrategyNumber strategy;
diff --git a/contrib/btree_gist/btree_ts.c b/contrib/btree_gist/btree_ts.c
index 18740cad38..49d1849d88 100644
--- a/contrib/btree_gist/btree_ts.c
+++ b/contrib/btree_gist/btree_ts.c
@@ -2,20 +2,21 @@
  * contrib/btree_gist/btree_ts.c
  */
 #include "postgres.h"
 
 #include <limits.h>
 
 #include "btree_gist.h"
 #include "btree_utils_num.h"
 #include "utils/builtins.h"
 #include "utils/datetime.h"
+#include "utils/float.h"
 
 typedef struct
 {
 	Timestamp	lower;
 	Timestamp	upper;
 } tsKEY;
 
 /*
 ** timestamp ops
 */
diff --git a/contrib/cube/cube.c b/contrib/cube/cube.c
index 3e3d83323c..ed41110a17 100644
--- a/contrib/cube/cube.c
+++ b/contrib/cube/cube.c
@@ -7,21 +7,21 @@
 ******************************************************************************/
 
 #include "postgres.h"
 
 #include <float.h>
 #include <math.h>
 
 #include "access/gist.h"
 #include "access/stratnum.h"
 #include "utils/array.h"
-#include "utils/builtins.h"
+#include "utils/float.h"
 
 #include "cubedata.h"
 
 PG_MODULE_MAGIC;
 
 /*
  * Taken from the intarray contrib header
  */
 #define ARRPTR(x)  ( (double *) ARR_DATA_PTR(x) )
 #define ARRNELEMS(x)  ArrayGetNItems( ARR_NDIM(x), ARR_DIMS(x))
diff --git a/contrib/cube/cubeparse.y b/contrib/cube/cubeparse.y
index 1b65fa967c..deb2efdc0d 100644
--- a/contrib/cube/cubeparse.y
+++ b/contrib/cube/cubeparse.y
@@ -1,20 +1,20 @@
 %{
 /* contrib/cube/cubeparse.y */
 
 /* NdBox = [(lowerleft),(upperright)] */
 /* [(xLL(1)...xLL(N)),(xUR(1)...xUR(n))] */
 
 #include "postgres.h"
 
 #include "cubedata.h"
-#include "utils/builtins.h"
+#include "utils/float.h"
 
 /* All grammar constructs return strings */
 #define YYSTYPE char *
 
 /*
  * Bison doesn't allocate anything that needs to live across parser calls,
  * so we can easily have it use palloc instead of malloc.  This prevents
  * memory leaks if we error out during parsing.  Note this only works with
  * bison >= 2.0.  However, in bison 1.875 the default is to use alloca()
  * if possible, so there's not really much problem anyhow, at least if
diff --git a/contrib/postgres_fdw/postgres_fdw.c b/contrib/postgres_fdw/postgres_fdw.c
index 7992ba5852..80043badc3 100644
--- a/contrib/postgres_fdw/postgres_fdw.c
+++ b/contrib/postgres_fdw/postgres_fdw.c
@@ -28,20 +28,21 @@
 #include "optimizer/cost.h"
 #include "optimizer/clauses.h"
 #include "optimizer/pathnode.h"
 #include "optimizer/paths.h"
 #include "optimizer/planmain.h"
 #include "optimizer/restrictinfo.h"
 #include "optimizer/var.h"
 #include "optimizer/tlist.h"
 #include "parser/parsetree.h"
 #include "utils/builtins.h"
+#include "utils/float.h"
 #include "utils/guc.h"
 #include "utils/lsyscache.h"
 #include "utils/memutils.h"
 #include "utils/rel.h"
 #include "utils/sampling.h"
 #include "utils/selfuncs.h"
 
 PG_MODULE_MAGIC;
 
 /* Default CPU cost to start up a foreign query. */
diff --git a/src/backend/access/gist/gistget.c b/src/backend/access/gist/gistget.c
index ca21cf7047..de0bcabb55 100644
--- a/src/backend/access/gist/gistget.c
+++ b/src/backend/access/gist/gistget.c
@@ -13,21 +13,21 @@
  *-------------------------------------------------------------------------
  */
 #include "postgres.h"
 
 #include "access/gist_private.h"
 #include "access/relscan.h"
 #include "catalog/pg_type.h"
 #include "miscadmin.h"
 #include "pgstat.h"
 #include "lib/pairingheap.h"
-#include "utils/builtins.h"
+#include "utils/float.h"
 #include "utils/memutils.h"
 #include "utils/rel.h"
 
 /*
  * gistkillitems() -- set LP_DEAD state for items an indexscan caller has
  * told us were killed.
  *
  * We re-read page here, so it's important to check page LSN. If the page
  * has been modified since the last read (as determined by LSN), we cannot
  * flag any entries because it is possible that the old entry was vacuumed
diff --git a/src/backend/access/gist/gistproc.c b/src/backend/access/gist/gistproc.c
index 97e6dc9910..4bbfdadf9f 100644
--- a/src/backend/access/gist/gistproc.c
+++ b/src/backend/access/gist/gistproc.c
@@ -27,62 +27,53 @@
 
 
 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
 
-/* Convenience macros for NaN-aware comparisons */
-#define FLOAT8_EQ(a,b)	(float8_cmp_internal(a, b) == 0)
-#define FLOAT8_LT(a,b)	(float8_cmp_internal(a, b) < 0)
-#define FLOAT8_LE(a,b)	(float8_cmp_internal(a, b) <= 0)
-#define FLOAT8_GT(a,b)	(float8_cmp_internal(a, b) > 0)
-#define FLOAT8_GE(a,b)	(float8_cmp_internal(a, b) >= 0)
-#define FLOAT8_MAX(a,b)  (FLOAT8_GT(a, b) ? (a) : (b))
-#define FLOAT8_MIN(a,b)  (FLOAT8_LT(a, b) ? (a) : (b))
-
 
 /**************************************************
  * Box ops
  **************************************************/
 
 /*
  * Calculates union of two boxes, a and b. The result is stored in *n.
  */
 static void
 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);
+	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
 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))
+	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);
@@ -137,27 +128,27 @@ gist_box_consistent(PG_FUNCTION_ARGS)
 												 query,
 												 strategy));
 }
 
 /*
  * Increase BOX b to include addon.
  */
 static void
 adjustBox(BOX *b, const BOX *addon)
 {
-	if (FLOAT8_LT(b->high.x, addon->high.x))
+	if (float8_lt(b->high.x, addon->high.x))
 		b->high.x = addon->high.x;
-	if (FLOAT8_GT(b->low.x, addon->low.x))
+	if (float8_gt(b->low.x, addon->low.x))
 		b->low.x = addon->low.x;
-	if (FLOAT8_LT(b->high.y, addon->high.y))
+	if (float8_lt(b->high.y, addon->high.y))
 		b->high.y = addon->high.y;
-	if (FLOAT8_GT(b->low.y, addon->low.y))
+	if (float8_gt(b->low.y, addon->low.y))
 		b->low.y = addon->low.y;
 }
 
 /*
  * The GiST Union method for boxes
  *
  * returns the minimal bounding box that encloses all the entries in entryvec
  */
 Datum
 gist_box_union(PG_FUNCTION_ARGS)
@@ -609,36 +600,36 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
 		i1 = 0;
 		i2 = 0;
 		rightLower = intervalsLower[i1].lower;
 		leftUpper = intervalsUpper[i2].lower;
 		while (true)
 		{
 			/*
 			 * Find next lower bound of right group.
 			 */
 			while (i1 < nentries &&
-				   FLOAT8_EQ(rightLower, intervalsLower[i1].lower))
+				   float8_eq(rightLower, intervalsLower[i1].lower))
 			{
-				if (FLOAT8_LT(leftUpper, intervalsLower[i1].upper))
+				if (float8_lt(leftUpper, intervalsLower[i1].upper))
 					leftUpper = intervalsLower[i1].upper;
 				i1++;
 			}
 			if (i1 >= nentries)
 				break;
 			rightLower = intervalsLower[i1].lower;
 
 			/*
 			 * Find count of intervals which anyway should be placed to the
 			 * left group.
 			 */
 			while (i2 < nentries &&
-				   FLOAT8_LE(intervalsUpper[i2].upper, leftUpper))
+				   float8_le(intervalsUpper[i2].upper, leftUpper))
 				i2++;
 
 			/*
 			 * Consider found split.
 			 */
 			g_box_consider_split(&context, dim, rightLower, i1, leftUpper, i2);
 		}
 
 		/*
 		 * Iterate over upper bound of left group finding greatest possible
@@ -646,35 +637,35 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
 		 */
 		i1 = nentries - 1;
 		i2 = nentries - 1;
 		rightLower = intervalsLower[i1].upper;
 		leftUpper = intervalsUpper[i2].upper;
 		while (true)
 		{
 			/*
 			 * Find next upper bound of left group.
 			 */
-			while (i2 >= 0 && FLOAT8_EQ(leftUpper, intervalsUpper[i2].upper))
+			while (i2 >= 0 && float8_eq(leftUpper, intervalsUpper[i2].upper))
 			{
-				if (FLOAT8_GT(rightLower, intervalsUpper[i2].lower))
+				if (float8_gt(rightLower, intervalsUpper[i2].lower))
 					rightLower = intervalsUpper[i2].lower;
 				i2--;
 			}
 			if (i2 < 0)
 				break;
 			leftUpper = intervalsUpper[i2].upper;
 
 			/*
 			 * Find count of intervals which anyway should be placed to the
 			 * right group.
 			 */
-			while (i1 >= 0 && FLOAT8_GE(intervalsLower[i1].lower, rightLower))
+			while (i1 >= 0 && float8_ge(intervalsLower[i1].lower, rightLower))
 				i1--;
 
 			/*
 			 * Consider found split.
 			 */
 			g_box_consider_split(&context, dim,
 								 rightLower, i1 + 1, leftUpper, i2 + 1);
 		}
 	}
 
@@ -748,42 +739,42 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
 		{
 			lower = box->low.x;
 			upper = box->high.x;
 		}
 		else
 		{
 			lower = box->low.y;
 			upper = box->high.y;
 		}
 
-		if (FLOAT8_LE(upper, context.leftUpper))
+		if (float8_le(upper, context.leftUpper))
 		{
 			/* Fits to the left group */
-			if (FLOAT8_GE(lower, context.rightLower))
+			if (float8_ge(lower, context.rightLower))
 			{
 				/* Fits also to the right group, so "common entry" */
 				commonEntries[commonEntriesCount++].index = i;
 			}
 			else
 			{
 				/* Doesn't fit to the right group, so join to the left group */
 				PLACE_LEFT(box, i);
 			}
 		}
 		else
 		{
 			/*
 			 * Each entry should fit on either left or right group. Since this
 			 * entry didn't fit on the left group, it better fit in the right
 			 * group.
 			 */
-			Assert(FLOAT8_GE(lower, context.rightLower));
+			Assert(float8_ge(lower, context.rightLower));
 
 			/* Doesn't fit to the left group, so join to the right group */
 			PLACE_RIGHT(box, i);
 		}
 	}
 
 	/*
 	 * Distribute "common entries", if any.
 	 */
 	if (commonEntriesCount > 0)
@@ -853,24 +844,24 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
  * equivalent to box_same().
  */
 Datum
 gist_box_same(PG_FUNCTION_ARGS)
 {
 	BOX		   *b1 = PG_GETARG_BOX_P(0);
 	BOX		   *b2 = PG_GETARG_BOX_P(1);
 	bool	   *result = (bool *) PG_GETARG_POINTER(2);
 
 	if (b1 && b2)
-		*result = (FLOAT8_EQ(b1->low.x, b2->low.x) &&
-				   FLOAT8_EQ(b1->low.y, b2->low.y) &&
-				   FLOAT8_EQ(b1->high.x, b2->high.x) &&
-				   FLOAT8_EQ(b1->high.y, b2->high.y));
+		*result = (float8_eq(b1->low.x, b2->low.x) &&
+				   float8_eq(b1->low.y, b2->low.y) &&
+				   float8_eq(b1->high.x, b2->high.x) &&
+				   float8_eq(b1->high.y, b2->high.y));
 	else
 		*result = (b1 == NULL && b2 == NULL);
 	PG_RETURN_POINTER(result);
 }
 
 /*
  * Leaf-level consistency for boxes: just apply the query operator
  */
 static bool
 gist_box_leaf_consistent(BOX *key, BOX *query, StrategyNumber strategy)
diff --git a/src/backend/access/gist/gistutil.c b/src/backend/access/gist/gistutil.c
index 55cccd247a..7e43b8fc87 100644
--- a/src/backend/access/gist/gistutil.c
+++ b/src/backend/access/gist/gistutil.c
@@ -15,21 +15,21 @@
 
 #include <float.h>
 #include <math.h>
 
 #include "access/gist_private.h"
 #include "access/htup_details.h"
 #include "access/reloptions.h"
 #include "catalog/pg_opclass.h"
 #include "storage/indexfsm.h"
 #include "storage/lmgr.h"
-#include "utils/builtins.h"
+#include "utils/float.h"
 #include "utils/syscache.h"
 
 
 /*
  * Write itup vector to page, has no control of free space.
  */
 void
 gistfillbuffer(Page page, IndexTuple *itup, int len, OffsetNumber off)
 {
 	OffsetNumber l = InvalidOffsetNumber;
diff --git a/src/backend/utils/adt/float.c b/src/backend/utils/adt/float.c
index bc6a3e09b5..e4d097c6a0 100644
--- a/src/backend/utils/adt/float.c
+++ b/src/backend/utils/adt/float.c
@@ -16,62 +16,29 @@
 
 #include <ctype.h>
 #include <float.h>
 #include <math.h>
 #include <limits.h>
 
 #include "catalog/pg_type.h"
 #include "common/int.h"
 #include "libpq/pqformat.h"
 #include "utils/array.h"
-#include "utils/builtins.h"
+#include "utils/float.h"
+#include "utils/fmgrprotos.h"
 #include "utils/sortsupport.h"
 
 
-#ifndef M_PI
-/* from my RH5.2 gcc math.h file - thomas 2000-04-03 */
-#define M_PI 3.14159265358979323846
-#endif
-
-/* Radians per degree, a.k.a. PI / 180 */
-#define RADIANS_PER_DEGREE 0.0174532925199432957692
-
-/* Visual C++ etc lacks NAN, and won't accept 0.0/0.0.  NAN definition from
- * http://msdn.microsoft.com/library/default.asp?url=/library/en-us/vclang/html/vclrfNotNumberNANItems.asp
- */
-#if defined(WIN32) && !defined(NAN)
-static const uint32 nan[2] = {0xffffffff, 0x7fffffff};
-
-#define NAN (*(const double *) nan)
-#endif
-
 /* not sure what the following should be, but better to make it over-sufficient */
 #define MAXFLOATWIDTH	64
 #define MAXDOUBLEWIDTH	128
 
-/*
- * check to see if a float4/8 val has underflowed or overflowed
- */
-#define CHECKFLOATVAL(val, inf_is_valid, zero_is_valid)			\
-do {															\
-	if (isinf(val) && !(inf_is_valid))							\
-		ereport(ERROR,											\
-				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),	\
-		  errmsg("value out of range: overflow")));				\
-																\
-	if ((val) == 0.0 && !(zero_is_valid))						\
-		ereport(ERROR,											\
-				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),	\
-		 errmsg("value out of range: underflow")));				\
-} while(0)
-
-
 /* Configurable GUC parameter */
 int			extra_float_digits = 0; /* Added to DBL_DIG or FLT_DIG */
 
 /* Cached constants for degree-based trig functions */
 static bool degree_consts_set = false;
 static float8 sin_30 = 0;
 static float8 one_minus_cos_60 = 0;
 static float8 asin_0_5 = 0;
 static float8 acos_0_5 = 0;
 static float8 atan_1_0 = 0;
@@ -102,100 +69,20 @@ static void init_degree_constants(void);
  * function, which causes configure to not set HAVE_CBRT.  Furthermore,
  * their compilers spit up at the mismatch between extern declaration
  * and static definition.  We work around that here by the expedient
  * of a #define to make the actual name of the static function different.
  */
 #define cbrt my_cbrt
 static double cbrt(double x);
 #endif							/* HAVE_CBRT */
 
 
-/*
- * Routines to provide reasonably platform-independent handling of
- * infinity and NaN.  We assume that isinf() and isnan() are available
- * and work per spec.  (On some platforms, we have to supply our own;
- * see src/port.)  However, generating an Infinity or NaN in the first
- * place is less well standardized; pre-C99 systems tend not to have C99's
- * INFINITY and NAN macros.  We centralize our workarounds for this here.
- */
-
-double
-get_float8_infinity(void)
-{
-#ifdef INFINITY
-	/* C99 standard way */
-	return (double) INFINITY;
-#else
-
-	/*
-	 * On some platforms, HUGE_VAL is an infinity, elsewhere it's just the
-	 * largest normal double.  We assume forcing an overflow will get us a
-	 * true infinity.
-	 */
-	return (double) (HUGE_VAL * HUGE_VAL);
-#endif
-}
-
-/*
-* The funny placements of the two #pragmas is necessary because of a
-* long lived bug in the Microsoft compilers.
-* See http://support.microsoft.com/kb/120968/en-us for details
-*/
-#if (_MSC_VER >= 1800)
-#pragma warning(disable:4756)
-#endif
-float
-get_float4_infinity(void)
-{
-#ifdef INFINITY
-	/* C99 standard way */
-	return (float) INFINITY;
-#else
-#if (_MSC_VER >= 1800)
-#pragma warning(default:4756)
-#endif
-
-	/*
-	 * On some platforms, HUGE_VAL is an infinity, elsewhere it's just the
-	 * largest normal double.  We assume forcing an overflow will get us a
-	 * true infinity.
-	 */
-	return (float) (HUGE_VAL * HUGE_VAL);
-#endif
-}
-
-double
-get_float8_nan(void)
-{
-	/* (double) NAN doesn't work on some NetBSD/MIPS releases */
-#if defined(NAN) && !(defined(__NetBSD__) && defined(__mips__))
-	/* C99 standard way */
-	return (double) NAN;
-#else
-	/* Assume we can get a NAN via zero divide */
-	return (double) (0.0 / 0.0);
-#endif
-}
-
-float
-get_float4_nan(void)
-{
-#ifdef NAN
-	/* C99 standard way */
-	return (float) NAN;
-#else
-	/* Assume we can get a NAN via zero divide */
-	return (float) (0.0 / 0.0);
-#endif
-}
-
-
 /*
  * Returns -1 if 'val' represents negative infinity, 1 if 'val'
  * represents (positive) infinity, and 0 otherwise. On some platforms,
  * this is equivalent to the isinf() macro, but not everywhere: C99
  * does not specify that isinf() needs to distinguish between positive
  * and negative infinity.
  */
 int
 is_infinite(double val)
 {
@@ -340,21 +227,21 @@ float4in(PG_FUNCTION_ARGS)
 	if (*endptr != '\0')
 		ereport(ERROR,
 				(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
 				 errmsg("invalid input syntax for type %s: \"%s\"",
 						"real", orig_num)));
 
 	/*
 	 * if we get here, we have a legal double, still need to check to see if
 	 * it's a legal float4
 	 */
-	CHECKFLOATVAL((float4) val, isinf(val), val == 0);
+	check_float4_val((float4) val, isinf(val), val == 0);
 
 	PG_RETURN_FLOAT4((float4) val);
 }
 
 /*
  *		float4out		- converts a float4 number to a string
  *						  using a standard output format
  */
 Datum
 float4out(PG_FUNCTION_ARGS)
@@ -690,35 +577,35 @@ float4up(PG_FUNCTION_ARGS)
 	PG_RETURN_FLOAT4(arg);
 }
 
 Datum
 float4larger(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
 	float4		result;
 
-	if (float4_cmp_internal(arg1, arg2) > 0)
+	if (float4_gt(arg1, arg2))
 		result = arg1;
 	else
 		result = arg2;
 	PG_RETURN_FLOAT4(result);
 }
 
 Datum
 float4smaller(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
 	float4		result;
 
-	if (float4_cmp_internal(arg1, arg2) < 0)
+	if (float4_lt(arg1, arg2))
 		result = arg1;
 	else
 		result = arg2;
 	PG_RETURN_FLOAT4(result);
 }
 
 /*
  *		======================
  *		FLOAT8 BASE OPERATIONS
  *		======================
@@ -757,35 +644,35 @@ float8up(PG_FUNCTION_ARGS)
 	PG_RETURN_FLOAT8(arg);
 }
 
 Datum
 float8larger(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
 	float8		result;
 
-	if (float8_cmp_internal(arg1, arg2) > 0)
+	if (float8_gt(arg1, arg2))
 		result = arg1;
 	else
 		result = arg2;
 	PG_RETURN_FLOAT8(result);
 }
 
 Datum
 float8smaller(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
 	float8		result;
 
-	if (float8_cmp_internal(arg1, arg2) < 0)
+	if (float8_lt(arg1, arg2))
 		result = arg1;
 	else
 		result = arg2;
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		====================
  *		ARITHMETIC OPERATORS
@@ -796,234 +683,165 @@ float8smaller(PG_FUNCTION_ARGS)
  *		float4pl		- returns arg1 + arg2
  *		float4mi		- returns arg1 - arg2
  *		float4mul		- returns arg1 * arg2
  *		float4div		- returns arg1 / arg2
  */
 Datum
 float4pl(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
-	float4		result;
 
-	result = arg1 + arg2;
-
-	/*
-	 * There isn't any way to check for underflow of addition/subtraction
-	 * because numbers near the underflow value have already been rounded to
-	 * the point where we can't detect that the two values were originally
-	 * different, e.g. on x86, '1e-45'::float4 == '2e-45'::float4 ==
-	 * 1.4013e-45.
-	 */
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
-	PG_RETURN_FLOAT4(result);
+	PG_RETURN_FLOAT4(float4_pl(arg1, arg2));
 }
 
 Datum
 float4mi(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
-	float4		result;
 
-	result = arg1 - arg2;
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
-	PG_RETURN_FLOAT4(result);
+	PG_RETURN_FLOAT4(float4_mi(arg1, arg2));
 }
 
 Datum
 float4mul(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
-	float4		result;
 
-	result = arg1 * arg2;
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2),
-				  arg1 == 0 || arg2 == 0);
-	PG_RETURN_FLOAT4(result);
+	PG_RETURN_FLOAT4(float4_mul(arg1, arg2));
 }
 
 Datum
 float4div(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
-	float4		result;
 
-	if (arg2 == 0.0)
-		ereport(ERROR,
-				(errcode(ERRCODE_DIVISION_BY_ZERO),
-				 errmsg("division by zero")));
-
-	result = arg1 / arg2;
-
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
-	PG_RETURN_FLOAT4(result);
+	PG_RETURN_FLOAT4(float4_div(arg1, arg2));
 }
 
 /*
  *		float8pl		- returns arg1 + arg2
  *		float8mi		- returns arg1 - arg2
  *		float8mul		- returns arg1 * arg2
  *		float8div		- returns arg1 / arg2
  */
 Datum
 float8pl(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
-	float8		result;
 
-	result = arg1 + arg2;
-
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
-	PG_RETURN_FLOAT8(result);
+	PG_RETURN_FLOAT8(float8_pl(arg1, arg2));
 }
 
 Datum
 float8mi(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
-	float8		result;
 
-	result = arg1 - arg2;
-
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
-	PG_RETURN_FLOAT8(result);
+	PG_RETURN_FLOAT8(float8_mi(arg1, arg2));
 }
 
 Datum
 float8mul(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
-	float8		result;
 
-	result = arg1 * arg2;
-
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2),
-				  arg1 == 0 || arg2 == 0);
-	PG_RETURN_FLOAT8(result);
+	PG_RETURN_FLOAT8(float8_mul(arg1, arg2));
 }
 
 Datum
 float8div(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
-	float8		result;
 
-	if (arg2 == 0.0)
-		ereport(ERROR,
-				(errcode(ERRCODE_DIVISION_BY_ZERO),
-				 errmsg("division by zero")));
-
-	result = arg1 / arg2;
-
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
-	PG_RETURN_FLOAT8(result);
+	PG_RETURN_FLOAT8(float8_div(arg1, arg2));
 }
 
 
 /*
  *		====================
  *		COMPARISON OPERATORS
  *		====================
  */
 
 /*
  *		float4{eq,ne,lt,le,gt,ge}		- float4/float4 comparison operations
  */
 int
 float4_cmp_internal(float4 a, float4 b)
 {
-	/*
-	 * We consider all NANs to be equal and larger than any non-NAN. This is
-	 * somewhat arbitrary; the important thing is to have a consistent sort
-	 * order.
-	 */
-	if (isnan(a))
-	{
-		if (isnan(b))
-			return 0;			/* NAN = NAN */
-		else
-			return 1;			/* NAN > non-NAN */
-	}
-	else if (isnan(b))
-	{
-		return -1;				/* non-NAN < NAN */
-	}
-	else
-	{
-		if (a > b)
-			return 1;
-		else if (a < b)
-			return -1;
-		else
-			return 0;
-	}
+	if (float4_gt(a, b))
+		return 1;
+	if (float4_lt(a, b))
+		return -1;
+	return 0;
 }
 
 Datum
 float4eq(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
 
-	PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) == 0);
+	PG_RETURN_BOOL(float4_eq(arg1, arg2));
 }
 
 Datum
 float4ne(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
 
-	PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) != 0);
+	PG_RETURN_BOOL(float4_ne(arg1, arg2));
 }
 
 Datum
 float4lt(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
 
-	PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) < 0);
+	PG_RETURN_BOOL(float4_lt(arg1, arg2));
 }
 
 Datum
 float4le(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
 
-	PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) <= 0);
+	PG_RETURN_BOOL(float4_le(arg1, arg2));
 }
 
 Datum
 float4gt(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
 
-	PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) > 0);
+	PG_RETURN_BOOL(float4_gt(arg1, arg2));
 }
 
 Datum
 float4ge(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
 
-	PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) >= 0);
+	PG_RETURN_BOOL(float4_ge(arg1, arg2));
 }
 
 Datum
 btfloat4cmp(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
 
 	PG_RETURN_INT32(float4_cmp_internal(arg1, arg2));
 }
@@ -1045,99 +863,79 @@ btfloat4sortsupport(PG_FUNCTION_ARGS)
 	ssup->comparator = btfloat4fastcmp;
 	PG_RETURN_VOID();
 }
 
 /*
  *		float8{eq,ne,lt,le,gt,ge}		- float8/float8 comparison operations
  */
 int
 float8_cmp_internal(float8 a, float8 b)
 {
-	/*
-	 * We consider all NANs to be equal and larger than any non-NAN. This is
-	 * somewhat arbitrary; the important thing is to have a consistent sort
-	 * order.
-	 */
-	if (isnan(a))
-	{
-		if (isnan(b))
-			return 0;			/* NAN = NAN */
-		else
-			return 1;			/* NAN > non-NAN */
-	}
-	else if (isnan(b))
-	{
-		return -1;				/* non-NAN < NAN */
-	}
-	else
-	{
-		if (a > b)
-			return 1;
-		else if (a < b)
-			return -1;
-		else
-			return 0;
-	}
+	if (float8_gt(a, b))
+		return 1;
+	if (float8_lt(a, b))
+		return -1;
+	return 0;
 }
 
 Datum
 float8eq(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) == 0);
+	PG_RETURN_BOOL(float8_eq(arg1, arg2));
 }
 
 Datum
 float8ne(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) != 0);
+	PG_RETURN_BOOL(float8_ne(arg1, arg2));
 }
 
 Datum
 float8lt(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) < 0);
+	PG_RETURN_BOOL(float8_lt(arg1, arg2));
 }
 
 Datum
 float8le(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) <= 0);
+	PG_RETURN_BOOL(float8_le(arg1, arg2));
 }
 
 Datum
 float8gt(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) > 0);
+	PG_RETURN_BOOL(float8_gt(arg1, arg2));
 }
 
 Datum
 float8ge(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) >= 0);
+	PG_RETURN_BOOL(float8_ge(arg1, arg2));
 }
 
 Datum
 btfloat8cmp(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
 
 	PG_RETURN_INT32(float8_cmp_internal(arg1, arg2));
 }
@@ -1200,21 +998,21 @@ ftod(PG_FUNCTION_ARGS)
 
 
 /*
  *		dtof			- converts a float8 number to a float4 number
  */
 Datum
 dtof(PG_FUNCTION_ARGS)
 {
 	float8		num = PG_GETARG_FLOAT8(0);
 
-	CHECKFLOATVAL((float4) num, isinf(num), num == 0);
+	check_float4_val((float4) num, isinf(num), num == 0);
 
 	PG_RETURN_FLOAT4((float4) num);
 }
 
 
 /*
  *		dtoi4			- converts a float8 number to an int4 number
  */
 Datum
 dtoi4(PG_FUNCTION_ARGS)
@@ -1425,36 +1223,36 @@ dsqrt(PG_FUNCTION_ARGS)
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float8		result;
 
 	if (arg1 < 0)
 		ereport(ERROR,
 				(errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
 				 errmsg("cannot take square root of a negative number")));
 
 	result = sqrt(arg1);
 
-	CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
+	check_float8_val(result, isinf(arg1), arg1 == 0);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		dcbrt			- returns cube root of arg1
  */
 Datum
 dcbrt(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float8		result;
 
 	result = cbrt(arg1);
-	CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
+	check_float8_val(result, isinf(arg1), arg1 == 0);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		dpow			- returns pow(arg1,arg2)
  */
 Datum
 dpow(PG_FUNCTION_ARGS)
 {
@@ -1493,40 +1291,40 @@ dpow(PG_FUNCTION_ARGS)
 			/* The sign of Inf is not significant in this case. */
 			result = get_float8_infinity();
 		else if (fabs(arg1) != 1)
 			result = 0;
 		else
 			result = 1;
 	}
 	else if (errno == ERANGE && result != 0 && !isinf(result))
 		result = get_float8_infinity();
 
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
+	check_float8_val(result, isinf(arg1) || isinf(arg2), arg1 == 0);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		dexp			- returns the exponential function of arg1
  */
 Datum
 dexp(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float8		result;
 
 	errno = 0;
 	result = exp(arg1);
 	if (errno == ERANGE && result != 0 && !isinf(result))
 		result = get_float8_infinity();
 
-	CHECKFLOATVAL(result, isinf(arg1), false);
+	check_float8_val(result, isinf(arg1), false);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		dlog1			- returns the natural logarithm of arg1
  */
 Datum
 dlog1(PG_FUNCTION_ARGS)
 {
@@ -1541,21 +1339,21 @@ dlog1(PG_FUNCTION_ARGS)
 		ereport(ERROR,
 				(errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
 				 errmsg("cannot take logarithm of zero")));
 	if (arg1 < 0)
 		ereport(ERROR,
 				(errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
 				 errmsg("cannot take logarithm of a negative number")));
 
 	result = log(arg1);
 
-	CHECKFLOATVAL(result, isinf(arg1), arg1 == 1);
+	check_float8_val(result, isinf(arg1), arg1 == 1);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		dlog10			- returns the base 10 logarithm of arg1
  */
 Datum
 dlog10(PG_FUNCTION_ARGS)
 {
@@ -1571,21 +1369,21 @@ dlog10(PG_FUNCTION_ARGS)
 		ereport(ERROR,
 				(errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
 				 errmsg("cannot take logarithm of zero")));
 	if (arg1 < 0)
 		ereport(ERROR,
 				(errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
 				 errmsg("cannot take logarithm of a negative number")));
 
 	result = log10(arg1);
 
-	CHECKFLOATVAL(result, isinf(arg1), arg1 == 1);
+	check_float8_val(result, isinf(arg1), arg1 == 1);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		dacos			- returns the arccos of arg1 (radians)
  */
 Datum
 dacos(PG_FUNCTION_ARGS)
 {
@@ -1601,21 +1399,21 @@ dacos(PG_FUNCTION_ARGS)
 	 * range [-1, 1] to values in the range [0, Pi], so we should reject any
 	 * inputs outside that range and the result will always be finite.
 	 */
 	if (arg1 < -1.0 || arg1 > 1.0)
 		ereport(ERROR,
 				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
 				 errmsg("input is out of range")));
 
 	result = acos(arg1);
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		dasin			- returns the arcsin of arg1 (radians)
  */
 Datum
 dasin(PG_FUNCTION_ARGS)
 {
@@ -1631,21 +1429,21 @@ dasin(PG_FUNCTION_ARGS)
 	 * range [-1, 1] to values in the range [-Pi/2, Pi/2], so we should reject
 	 * any inputs outside that range and the result will always be finite.
 	 */
 	if (arg1 < -1.0 || arg1 > 1.0)
 		ereport(ERROR,
 				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
 				 errmsg("input is out of range")));
 
 	result = asin(arg1);
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		datan			- returns the arctan of arg1 (radians)
  */
 Datum
 datan(PG_FUNCTION_ARGS)
 {
@@ -1656,21 +1454,21 @@ datan(PG_FUNCTION_ARGS)
 	if (isnan(arg1))
 		PG_RETURN_FLOAT8(get_float8_nan());
 
 	/*
 	 * The principal branch of the inverse tangent function maps all inputs to
 	 * values in the range [-Pi/2, Pi/2], so the result should always be
 	 * finite, even if the input is infinite.
 	 */
 	result = atan(arg1);
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		atan2			- returns the arctan of arg1/arg2 (radians)
  */
 Datum
 datan2(PG_FUNCTION_ARGS)
 {
@@ -1681,21 +1479,21 @@ datan2(PG_FUNCTION_ARGS)
 	/* Per the POSIX spec, return NaN if either input is NaN */
 	if (isnan(arg1) || isnan(arg2))
 		PG_RETURN_FLOAT8(get_float8_nan());
 
 	/*
 	 * atan2 maps all inputs to values in the range [-Pi, Pi], so the result
 	 * should always be finite, even if the inputs are infinite.
 	 */
 	result = atan2(arg1, arg2);
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		dcos			- returns the cosine of arg1 (radians)
  */
 Datum
 dcos(PG_FUNCTION_ARGS)
 {
@@ -1721,21 +1519,21 @@ dcos(PG_FUNCTION_ARGS)
 	 * platform reports via errno, so also explicitly test for infinite
 	 * inputs.
 	 */
 	errno = 0;
 	result = cos(arg1);
 	if (errno != 0 || isinf(arg1))
 		ereport(ERROR,
 				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
 				 errmsg("input is out of range")));
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		dcot			- returns the cotangent of arg1 (radians)
  */
 Datum
 dcot(PG_FUNCTION_ARGS)
 {
@@ -1748,21 +1546,21 @@ dcot(PG_FUNCTION_ARGS)
 
 	/* Be sure to throw an error if the input is infinite --- see dcos() */
 	errno = 0;
 	result = tan(arg1);
 	if (errno != 0 || isinf(arg1))
 		ereport(ERROR,
 				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
 				 errmsg("input is out of range")));
 
 	result = 1.0 / result;
-	CHECKFLOATVAL(result, true /* cot(0) == Inf */ , true);
+	check_float8_val(result, true /* cot(0) == Inf */ , true);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		dsin			- returns the sine of arg1 (radians)
  */
 Datum
 dsin(PG_FUNCTION_ARGS)
 {
@@ -1774,21 +1572,21 @@ dsin(PG_FUNCTION_ARGS)
 		PG_RETURN_FLOAT8(get_float8_nan());
 
 	/* Be sure to throw an error if the input is infinite --- see dcos() */
 	errno = 0;
 	result = sin(arg1);
 	if (errno != 0 || isinf(arg1))
 		ereport(ERROR,
 				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
 				 errmsg("input is out of range")));
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		dtan			- returns the tangent of arg1 (radians)
  */
 Datum
 dtan(PG_FUNCTION_ARGS)
 {
@@ -1800,21 +1598,21 @@ dtan(PG_FUNCTION_ARGS)
 		PG_RETURN_FLOAT8(get_float8_nan());
 
 	/* Be sure to throw an error if the input is infinite --- see dcos() */
 	errno = 0;
 	result = tan(arg1);
 	if (errno != 0 || isinf(arg1))
 		ereport(ERROR,
 				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
 				 errmsg("input is out of range")));
 
-	CHECKFLOATVAL(result, true /* tan(pi/2) == Inf */ , true);
+	check_float8_val(result, true /* tan(pi/2) == Inf */ , true);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /* ========== DEGREE-BASED TRIGONOMETRIC FUNCTIONS ========== */
 
 
 /*
  * Initialize the cached constants declared at the head of this file
  * (sin_30 etc).  The fact that we need those at all, let alone need this
@@ -1952,21 +1750,21 @@ dacosd(PG_FUNCTION_ARGS)
 	if (arg1 < -1.0 || arg1 > 1.0)
 		ereport(ERROR,
 				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
 				 errmsg("input is out of range")));
 
 	if (arg1 >= 0.0)
 		result = acosd_q1(arg1);
 	else
 		result = 90.0 + asind_q1(-arg1);
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		dasind			- returns the arcsin of arg1 (degrees)
  */
 Datum
 dasind(PG_FUNCTION_ARGS)
 {
@@ -1987,21 +1785,21 @@ dasind(PG_FUNCTION_ARGS)
 	if (arg1 < -1.0 || arg1 > 1.0)
 		ereport(ERROR,
 				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
 				 errmsg("input is out of range")));
 
 	if (arg1 >= 0.0)
 		result = asind_q1(arg1);
 	else
 		result = -asind_q1(-arg1);
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		datand			- returns the arctan of arg1 (degrees)
  */
 Datum
 datand(PG_FUNCTION_ARGS)
 {
@@ -2017,21 +1815,21 @@ datand(PG_FUNCTION_ARGS)
 
 	/*
 	 * The principal branch of the inverse tangent function maps all inputs to
 	 * values in the range [-90, 90], so the result should always be finite,
 	 * even if the input is infinite.  Additionally, we take care to ensure
 	 * than when arg1 is 1, the result is exactly 45.
 	 */
 	atan_arg1 = atan(arg1);
 	result = (atan_arg1 / atan_1_0) * 45.0;
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		atan2d			- returns the arctan of arg1/arg2 (degrees)
  */
 Datum
 datan2d(PG_FUNCTION_ARGS)
 {
@@ -2051,21 +1849,21 @@ datan2d(PG_FUNCTION_ARGS)
 	 * result should always be finite, even if the inputs are infinite.
 	 *
 	 * Note: this coding assumes that atan(1.0) is a suitable scaling constant
 	 * to get an exact result from atan2().  This might well fail on us at
 	 * some point, requiring us to decide exactly what inputs we think we're
 	 * going to guarantee an exact result for.
 	 */
 	atan2_arg1_arg2 = atan2(arg1, arg2);
 	result = (atan2_arg1_arg2 / atan_1_0) * 45.0;
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		sind_0_to_30	- returns the sine of an angle that lies between 0 and
  *						  30 degrees.  This will return exactly 0 when x is 0,
  *						  and exactly 0.5 when x is 30 degrees.
  */
 static double
@@ -2172,21 +1970,21 @@ dcosd(PG_FUNCTION_ARGS)
 
 	if (arg1 > 90.0)
 	{
 		/* cosd(180-x) = -cosd(x) */
 		arg1 = 180.0 - arg1;
 		sign = -sign;
 	}
 
 	result = sign * cosd_q1(arg1);
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		dcotd			- returns the cotangent of arg1 (degrees)
  */
 Datum
 dcotd(PG_FUNCTION_ARGS)
 {
@@ -2237,21 +2035,21 @@ dcotd(PG_FUNCTION_ARGS)
 	result = sign * (cot_arg1 / cot_45);
 
 	/*
 	 * On some machines we get cotd(270) = minus zero, but this isn't always
 	 * true.  For portability, and because the user constituency for this
 	 * function probably doesn't want minus zero, force it to plain zero.
 	 */
 	if (result == 0.0)
 		result = 0.0;
 
-	CHECKFLOATVAL(result, true /* cotd(0) == Inf */ , true);
+	check_float8_val(result, true /* cotd(0) == Inf */ , true);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		dsind			- returns the sine of arg1 (degrees)
  */
 Datum
 dsind(PG_FUNCTION_ARGS)
 {
@@ -2291,21 +2089,21 @@ dsind(PG_FUNCTION_ARGS)
 	}
 
 	if (arg1 > 90.0)
 	{
 		/* sind(180-x) = sind(x) */
 		arg1 = 180.0 - arg1;
 	}
 
 	result = sign * sind_q1(arg1);
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		dtand			- returns the tangent of arg1 (degrees)
  */
 Datum
 dtand(PG_FUNCTION_ARGS)
 {
@@ -2356,64 +2154,56 @@ dtand(PG_FUNCTION_ARGS)
 	result = sign * (tan_arg1 / tan_45);
 
 	/*
 	 * On some machines we get tand(180) = minus zero, but this isn't always
 	 * true.  For portability, and because the user constituency for this
 	 * function probably doesn't want minus zero, force it to plain zero.
 	 */
 	if (result == 0.0)
 		result = 0.0;
 
-	CHECKFLOATVAL(result, true /* tand(90) == Inf */ , true);
+	check_float8_val(result, true /* tand(90) == Inf */ , true);
 	PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *		degrees		- returns degrees converted from radians
  */
 Datum
 degrees(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
-	float8		result;
 
-	result = arg1 / RADIANS_PER_DEGREE;
-
-	CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
-	PG_RETURN_FLOAT8(result);
+	PG_RETURN_FLOAT8(float8_div(arg1, RADIANS_PER_DEGREE));
 }
 
 
 /*
  *		dpi				- returns the constant PI
  */
 Datum
 dpi(PG_FUNCTION_ARGS)
 {
 	PG_RETURN_FLOAT8(M_PI);
 }
 
 
 /*
  *		radians		- returns radians converted from degrees
  */
 Datum
 radians(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
-	float8		result;
 
-	result = arg1 * RADIANS_PER_DEGREE;
-
-	CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
-	PG_RETURN_FLOAT8(result);
+	PG_RETURN_FLOAT8(float8_mul(arg1, RADIANS_PER_DEGREE));
 }
 
 
 /*
  *		drandom		- returns a random number
  */
 Datum
 drandom(PG_FUNCTION_ARGS)
 {
 	float8		result;
@@ -2491,144 +2281,105 @@ check_float8_array(ArrayType *transarray, const char *caller, int n)
  * This function is used only in two stage aggregation and
  * shouldn't be called outside aggregate context.
  */
 Datum
 float8_combine(PG_FUNCTION_ARGS)
 {
 	ArrayType  *transarray1 = PG_GETARG_ARRAYTYPE_P(0);
 	ArrayType  *transarray2 = PG_GETARG_ARRAYTYPE_P(1);
 	float8	   *transvalues1;
 	float8	   *transvalues2;
-	float8		N,
-				sumX,
-				sumX2;
 
 	if (!AggCheckCallContext(fcinfo, NULL))
 		elog(ERROR, "aggregate function called in non-aggregate context");
 
 	transvalues1 = check_float8_array(transarray1, "float8_combine", 3);
-	N = transvalues1[0];
-	sumX = transvalues1[1];
-	sumX2 = transvalues1[2];
-
 	transvalues2 = check_float8_array(transarray2, "float8_combine", 3);
 
-	N += transvalues2[0];
-	sumX += transvalues2[1];
-	CHECKFLOATVAL(sumX, isinf(transvalues1[1]) || isinf(transvalues2[1]),
-				  true);
-	sumX2 += transvalues2[2];
-	CHECKFLOATVAL(sumX2, isinf(transvalues1[2]) || isinf(transvalues2[2]),
-				  true);
-
-	transvalues1[0] = N;
-	transvalues1[1] = sumX;
-	transvalues1[2] = sumX2;
+	transvalues1[0] = transvalues1[0] + transvalues2[0];
+	transvalues1[1] = float8_pl(transvalues1[1], transvalues2[1]);
+	transvalues1[2] = float8_pl(transvalues1[2], transvalues2[2]);
 
 	PG_RETURN_ARRAYTYPE_P(transarray1);
 }
 
 Datum
 float8_accum(PG_FUNCTION_ARGS)
 {
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8		newval = PG_GETARG_FLOAT8(1);
 	float8	   *transvalues;
-	float8		N,
-				sumX,
-				sumX2;
 
 	transvalues = check_float8_array(transarray, "float8_accum", 3);
-	N = transvalues[0];
-	sumX = transvalues[1];
-	sumX2 = transvalues[2];
-
-	N += 1.0;
-	sumX += newval;
-	CHECKFLOATVAL(sumX, isinf(transvalues[1]) || isinf(newval), true);
-	sumX2 += newval * newval;
-	CHECKFLOATVAL(sumX2, isinf(transvalues[2]) || isinf(newval), true);
 
 	/*
 	 * If we're invoked as an aggregate, we can cheat and modify our first
 	 * parameter in-place to reduce palloc overhead. Otherwise we construct a
 	 * new array with the updated transition data and return it.
 	 */
 	if (AggCheckCallContext(fcinfo, NULL))
 	{
-		transvalues[0] = N;
-		transvalues[1] = sumX;
-		transvalues[2] = sumX2;
+		transvalues[0] = transvalues[0] + 1.0;
+		transvalues[1] = float8_pl(transvalues[1], newval);
+		transvalues[2] = float8_pl(transvalues[2], float8_mul(newval, newval));
 
 		PG_RETURN_ARRAYTYPE_P(transarray);
 	}
 	else
 	{
 		Datum		transdatums[3];
 		ArrayType  *result;
 
-		transdatums[0] = Float8GetDatumFast(N);
-		transdatums[1] = Float8GetDatumFast(sumX);
-		transdatums[2] = Float8GetDatumFast(sumX2);
+		transvalues[0] = transvalues[0] + 1.0;
+		transvalues[1] = float8_pl(transvalues[1], newval);
+		transvalues[2] = float8_pl(transvalues[2], float8_mul(newval, newval));
 
 		result = construct_array(transdatums, 3,
 								 FLOAT8OID,
 								 sizeof(float8), FLOAT8PASSBYVAL, 'd');
 
 		PG_RETURN_ARRAYTYPE_P(result);
 	}
 }
 
 Datum
 float4_accum(PG_FUNCTION_ARGS)
 {
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 
 	/* do computations as float8 */
 	float8		newval = PG_GETARG_FLOAT4(1);
 	float8	   *transvalues;
-	float8		N,
-				sumX,
-				sumX2;
 
 	transvalues = check_float8_array(transarray, "float4_accum", 3);
-	N = transvalues[0];
-	sumX = transvalues[1];
-	sumX2 = transvalues[2];
-
-	N += 1.0;
-	sumX += newval;
-	CHECKFLOATVAL(sumX, isinf(transvalues[1]) || isinf(newval), true);
-	sumX2 += newval * newval;
-	CHECKFLOATVAL(sumX2, isinf(transvalues[2]) || isinf(newval), true);
 
 	/*
 	 * If we're invoked as an aggregate, we can cheat and modify our first
 	 * parameter in-place to reduce palloc overhead. Otherwise we construct a
 	 * new array with the updated transition data and return it.
 	 */
 	if (AggCheckCallContext(fcinfo, NULL))
 	{
-		transvalues[0] = N;
-		transvalues[1] = sumX;
-		transvalues[2] = sumX2;
+		transvalues[0] = transvalues[0] + 1.0;
+		transvalues[1] = float8_pl(transvalues[1], newval);
+		transvalues[2] = float8_pl(transvalues[2], float8_mul(newval, newval));
 
 		PG_RETURN_ARRAYTYPE_P(transarray);
 	}
 	else
 	{
 		Datum		transdatums[3];
 		ArrayType  *result;
 
-		transdatums[0] = Float8GetDatumFast(N);
-		transdatums[1] = Float8GetDatumFast(sumX);
-		transdatums[2] = Float8GetDatumFast(sumX2);
+		transvalues[0] = transvalues[0] + 1.0;
+		transvalues[1] = float8_pl(transvalues[1], newval);
+		transvalues[2] = float8_pl(transvalues[2], float8_mul(newval, newval));
 
 		result = construct_array(transdatums, 3,
 								 FLOAT8OID,
 								 sizeof(float8), FLOAT8PASSBYVAL, 'd');
 
 		PG_RETURN_ARRAYTYPE_P(result);
 	}
 }
 
 Datum
@@ -2664,21 +2415,21 @@ float8_var_pop(PG_FUNCTION_ARGS)
 	transvalues = check_float8_array(transarray, "float8_var_pop", 3);
 	N = transvalues[0];
 	sumX = transvalues[1];
 	sumX2 = transvalues[2];
 
 	/* Population variance is undefined when N is 0, so return NULL */
 	if (N == 0.0)
 		PG_RETURN_NULL();
 
 	numerator = N * sumX2 - sumX * sumX;
-	CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true);
+	check_float8_val(numerator, isinf(sumX2) || isinf(sumX), true);
 
 	/* Watch out for roundoff error producing a negative numerator */
 	if (numerator <= 0.0)
 		PG_RETURN_FLOAT8(0.0);
 
 	PG_RETURN_FLOAT8(numerator / (N * N));
 }
 
 Datum
 float8_var_samp(PG_FUNCTION_ARGS)
@@ -2693,21 +2444,21 @@ float8_var_samp(PG_FUNCTION_ARGS)
 	transvalues = check_float8_array(transarray, "float8_var_samp", 3);
 	N = transvalues[0];
 	sumX = transvalues[1];
 	sumX2 = transvalues[2];
 
 	/* Sample variance is undefined when N is 0 or 1, so return NULL */
 	if (N <= 1.0)
 		PG_RETURN_NULL();
 
 	numerator = N * sumX2 - sumX * sumX;
-	CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true);
+	check_float8_val(numerator, isinf(sumX2) || isinf(sumX), true);
 
 	/* Watch out for roundoff error producing a negative numerator */
 	if (numerator <= 0.0)
 		PG_RETURN_FLOAT8(0.0);
 
 	PG_RETURN_FLOAT8(numerator / (N * (N - 1.0)));
 }
 
 Datum
 float8_stddev_pop(PG_FUNCTION_ARGS)
@@ -2722,21 +2473,21 @@ float8_stddev_pop(PG_FUNCTION_ARGS)
 	transvalues = check_float8_array(transarray, "float8_stddev_pop", 3);
 	N = transvalues[0];
 	sumX = transvalues[1];
 	sumX2 = transvalues[2];
 
 	/* Population stddev is undefined when N is 0, so return NULL */
 	if (N == 0.0)
 		PG_RETURN_NULL();
 
 	numerator = N * sumX2 - sumX * sumX;
-	CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true);
+	check_float8_val(numerator, isinf(sumX2) || isinf(sumX), true);
 
 	/* Watch out for roundoff error producing a negative numerator */
 	if (numerator <= 0.0)
 		PG_RETURN_FLOAT8(0.0);
 
 	PG_RETURN_FLOAT8(sqrt(numerator / (N * N)));
 }
 
 Datum
 float8_stddev_samp(PG_FUNCTION_ARGS)
@@ -2751,21 +2502,21 @@ float8_stddev_samp(PG_FUNCTION_ARGS)
 	transvalues = check_float8_array(transarray, "float8_stddev_samp", 3);
 	N = transvalues[0];
 	sumX = transvalues[1];
 	sumX2 = transvalues[2];
 
 	/* Sample stddev is undefined when N is 0 or 1, so return NULL */
 	if (N <= 1.0)
 		PG_RETURN_NULL();
 
 	numerator = N * sumX2 - sumX * sumX;
-	CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true);
+	check_float8_val(numerator, isinf(sumX2) || isinf(sumX), true);
 
 	/* Watch out for roundoff error producing a negative numerator */
 	if (numerator <= 0.0)
 		PG_RETURN_FLOAT8(0.0);
 
 	PG_RETURN_FLOAT8(sqrt(numerator / (N * (N - 1.0))));
 }
 
 /*
  *		=========================
@@ -2800,30 +2551,30 @@ float8_regr_accum(PG_FUNCTION_ARGS)
 	transvalues = check_float8_array(transarray, "float8_regr_accum", 6);
 	N = transvalues[0];
 	sumX = transvalues[1];
 	sumX2 = transvalues[2];
 	sumY = transvalues[3];
 	sumY2 = transvalues[4];
 	sumXY = transvalues[5];
 
 	N += 1.0;
 	sumX += newvalX;
-	CHECKFLOATVAL(sumX, isinf(transvalues[1]) || isinf(newvalX), true);
+	check_float8_val(sumX, isinf(transvalues[1]) || isinf(newvalX), true);
 	sumX2 += newvalX * newvalX;
-	CHECKFLOATVAL(sumX2, isinf(transvalues[2]) || isinf(newvalX), true);
+	check_float8_val(sumX2, isinf(transvalues[2]) || isinf(newvalX), true);
 	sumY += newvalY;
-	CHECKFLOATVAL(sumY, isinf(transvalues[3]) || isinf(newvalY), true);
+	check_float8_val(sumY, isinf(transvalues[3]) || isinf(newvalY), true);
 	sumY2 += newvalY * newvalY;
-	CHECKFLOATVAL(sumY2, isinf(transvalues[4]) || isinf(newvalY), true);
+	check_float8_val(sumY2, isinf(transvalues[4]) || isinf(newvalY), true);
 	sumXY += newvalX * newvalY;
-	CHECKFLOATVAL(sumXY, isinf(transvalues[5]) || isinf(newvalX) ||
-				  isinf(newvalY), true);
+	check_float8_val(sumXY, isinf(transvalues[5]) || isinf(newvalX) ||
+					 isinf(newvalY), true);
 
 	/*
 	 * If we're invoked as an aggregate, we can cheat and modify our first
 	 * parameter in-place to reduce palloc overhead. Otherwise we construct a
 	 * new array with the updated transition data and return it.
 	 */
 	if (AggCheckCallContext(fcinfo, NULL))
 	{
 		transvalues[0] = N;
 		transvalues[1] = sumX;
@@ -2862,63 +2613,33 @@ float8_regr_accum(PG_FUNCTION_ARGS)
  * This function is used only in two stage aggregation and
  * shouldn't be called outside aggregate context.
  */
 Datum
 float8_regr_combine(PG_FUNCTION_ARGS)
 {
 	ArrayType  *transarray1 = PG_GETARG_ARRAYTYPE_P(0);
 	ArrayType  *transarray2 = PG_GETARG_ARRAYTYPE_P(1);
 	float8	   *transvalues1;
 	float8	   *transvalues2;
-	float8		N,
-				sumX,
-				sumX2,
-				sumY,
-				sumY2,
-				sumXY;
 
 	if (!AggCheckCallContext(fcinfo, NULL))
 		elog(ERROR, "aggregate function called in non-aggregate context");
 
 	transvalues1 = check_float8_array(transarray1, "float8_regr_combine", 6);
-	N = transvalues1[0];
-	sumX = transvalues1[1];
-	sumX2 = transvalues1[2];
-	sumY = transvalues1[3];
-	sumY2 = transvalues1[4];
-	sumXY = transvalues1[5];
-
 	transvalues2 = check_float8_array(transarray2, "float8_regr_combine", 6);
 
-	N += transvalues2[0];
-	sumX += transvalues2[1];
-	CHECKFLOATVAL(sumX, isinf(transvalues1[1]) || isinf(transvalues2[1]),
-				  true);
-	sumX2 += transvalues2[2];
-	CHECKFLOATVAL(sumX2, isinf(transvalues1[2]) || isinf(transvalues2[2]),
-				  true);
-	sumY += transvalues2[3];
-	CHECKFLOATVAL(sumY, isinf(transvalues1[3]) || isinf(transvalues2[3]),
-				  true);
-	sumY2 += transvalues2[4];
-	CHECKFLOATVAL(sumY2, isinf(transvalues1[4]) || isinf(transvalues2[4]),
-				  true);
-	sumXY += transvalues2[5];
-	CHECKFLOATVAL(sumXY, isinf(transvalues1[5]) || isinf(transvalues2[5]),
-				  true);
-
-	transvalues1[0] = N;
-	transvalues1[1] = sumX;
-	transvalues1[2] = sumX2;
-	transvalues1[3] = sumY;
-	transvalues1[4] = sumY2;
-	transvalues1[5] = sumXY;
+	transvalues1[0] = transvalues1[0] + transvalues2[0];
+	transvalues1[1] = float8_pl(transvalues1[1], transvalues2[1]);
+	transvalues1[2] = float8_pl(transvalues1[2], transvalues2[2]);
+	transvalues1[3] = float8_pl(transvalues1[3], transvalues2[3]);
+	transvalues1[4] = float8_pl(transvalues1[4], transvalues2[4]);
+	transvalues1[5] = float8_pl(transvalues1[5], transvalues2[5]);
 
 	PG_RETURN_ARRAYTYPE_P(transarray1);
 }
 
 
 Datum
 float8_regr_sxx(PG_FUNCTION_ARGS)
 {
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
@@ -2930,21 +2651,21 @@ float8_regr_sxx(PG_FUNCTION_ARGS)
 	transvalues = check_float8_array(transarray, "float8_regr_sxx", 6);
 	N = transvalues[0];
 	sumX = transvalues[1];
 	sumX2 = transvalues[2];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
 	numerator = N * sumX2 - sumX * sumX;
-	CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true);
+	check_float8_val(numerator, isinf(sumX2) || isinf(sumX), true);
 
 	/* Watch out for roundoff error producing a negative numerator */
 	if (numerator <= 0.0)
 		PG_RETURN_FLOAT8(0.0);
 
 	PG_RETURN_FLOAT8(numerator / N);
 }
 
 Datum
 float8_regr_syy(PG_FUNCTION_ARGS)
@@ -2959,21 +2680,21 @@ float8_regr_syy(PG_FUNCTION_ARGS)
 	transvalues = check_float8_array(transarray, "float8_regr_syy", 6);
 	N = transvalues[0];
 	sumY = transvalues[3];
 	sumY2 = transvalues[4];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
 	numerator = N * sumY2 - sumY * sumY;
-	CHECKFLOATVAL(numerator, isinf(sumY2) || isinf(sumY), true);
+	check_float8_val(numerator, isinf(sumY2) || isinf(sumY), true);
 
 	/* Watch out for roundoff error producing a negative numerator */
 	if (numerator <= 0.0)
 		PG_RETURN_FLOAT8(0.0);
 
 	PG_RETURN_FLOAT8(numerator / N);
 }
 
 Datum
 float8_regr_sxy(PG_FUNCTION_ARGS)
@@ -2990,22 +2711,22 @@ float8_regr_sxy(PG_FUNCTION_ARGS)
 	N = transvalues[0];
 	sumX = transvalues[1];
 	sumY = transvalues[3];
 	sumXY = transvalues[5];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
 	numerator = N * sumXY - sumX * sumY;
-	CHECKFLOATVAL(numerator, isinf(sumXY) || isinf(sumX) ||
-				  isinf(sumY), true);
+	check_float8_val(numerator, isinf(sumXY) || isinf(sumX) ||
+					 isinf(sumY), true);
 
 	/* A negative result is valid here */
 
 	PG_RETURN_FLOAT8(numerator / N);
 }
 
 Datum
 float8_regr_avgx(PG_FUNCTION_ARGS)
 {
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
@@ -3058,22 +2779,22 @@ float8_covar_pop(PG_FUNCTION_ARGS)
 	N = transvalues[0];
 	sumX = transvalues[1];
 	sumY = transvalues[3];
 	sumXY = transvalues[5];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
 	numerator = N * sumXY - sumX * sumY;
-	CHECKFLOATVAL(numerator, isinf(sumXY) || isinf(sumX) ||
-				  isinf(sumY), true);
+	check_float8_val(numerator, isinf(sumXY) || isinf(sumX) ||
+					 isinf(sumY), true);
 
 	PG_RETURN_FLOAT8(numerator / (N * N));
 }
 
 Datum
 float8_covar_samp(PG_FUNCTION_ARGS)
 {
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
@@ -3086,22 +2807,22 @@ float8_covar_samp(PG_FUNCTION_ARGS)
 	N = transvalues[0];
 	sumX = transvalues[1];
 	sumY = transvalues[3];
 	sumXY = transvalues[5];
 
 	/* if N is <= 1 we should return NULL */
 	if (N < 2.0)
 		PG_RETURN_NULL();
 
 	numerator = N * sumXY - sumX * sumY;
-	CHECKFLOATVAL(numerator, isinf(sumXY) || isinf(sumX) ||
-				  isinf(sumY), true);
+	check_float8_val(numerator, isinf(sumXY) || isinf(sumX) ||
+					 isinf(sumY), true);
 
 	PG_RETURN_FLOAT8(numerator / (N * (N - 1.0)));
 }
 
 Datum
 float8_corr(PG_FUNCTION_ARGS)
 {
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
 	float8	   *transvalues;
 	float8		N,
@@ -3120,26 +2841,26 @@ float8_corr(PG_FUNCTION_ARGS)
 	sumX2 = transvalues[2];
 	sumY = transvalues[3];
 	sumY2 = transvalues[4];
 	sumXY = transvalues[5];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
 	numeratorX = N * sumX2 - sumX * sumX;
-	CHECKFLOATVAL(numeratorX, isinf(sumX2) || isinf(sumX), true);
+	check_float8_val(numeratorX, isinf(sumX2) || isinf(sumX), true);
 	numeratorY = N * sumY2 - sumY * sumY;
-	CHECKFLOATVAL(numeratorY, isinf(sumY2) || isinf(sumY), true);
+	check_float8_val(numeratorY, isinf(sumY2) || isinf(sumY), true);
 	numeratorXY = N * sumXY - sumX * sumY;
-	CHECKFLOATVAL(numeratorXY, isinf(sumXY) || isinf(sumX) ||
-				  isinf(sumY), true);
+	check_float8_val(numeratorXY, isinf(sumXY) || isinf(sumX) ||
+					 isinf(sumY), true);
 	if (numeratorX <= 0 || numeratorY <= 0)
 		PG_RETURN_NULL();
 
 	PG_RETURN_FLOAT8(numeratorXY / sqrt(numeratorX * numeratorY));
 }
 
 Datum
 float8_regr_r2(PG_FUNCTION_ARGS)
 {
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
@@ -3160,26 +2881,26 @@ float8_regr_r2(PG_FUNCTION_ARGS)
 	sumX2 = transvalues[2];
 	sumY = transvalues[3];
 	sumY2 = transvalues[4];
 	sumXY = transvalues[5];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
 	numeratorX = N * sumX2 - sumX * sumX;
-	CHECKFLOATVAL(numeratorX, isinf(sumX2) || isinf(sumX), true);
+	check_float8_val(numeratorX, isinf(sumX2) || isinf(sumX), true);
 	numeratorY = N * sumY2 - sumY * sumY;
-	CHECKFLOATVAL(numeratorY, isinf(sumY2) || isinf(sumY), true);
+	check_float8_val(numeratorY, isinf(sumY2) || isinf(sumY), true);
 	numeratorXY = N * sumXY - sumX * sumY;
-	CHECKFLOATVAL(numeratorXY, isinf(sumXY) || isinf(sumX) ||
-				  isinf(sumY), true);
+	check_float8_val(numeratorXY, isinf(sumXY) || isinf(sumX) ||
+					 isinf(sumY), true);
 	if (numeratorX <= 0)
 		PG_RETURN_NULL();
 	/* per spec, horizontal line produces 1.0 */
 	if (numeratorY <= 0)
 		PG_RETURN_FLOAT8(1.0);
 
 	PG_RETURN_FLOAT8((numeratorXY * numeratorXY) /
 					 (numeratorX * numeratorY));
 }
 
@@ -3201,24 +2922,24 @@ float8_regr_slope(PG_FUNCTION_ARGS)
 	sumX = transvalues[1];
 	sumX2 = transvalues[2];
 	sumY = transvalues[3];
 	sumXY = transvalues[5];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
 	numeratorX = N * sumX2 - sumX * sumX;
-	CHECKFLOATVAL(numeratorX, isinf(sumX2) || isinf(sumX), true);
+	check_float8_val(numeratorX, isinf(sumX2) || isinf(sumX), true);
 	numeratorXY = N * sumXY - sumX * sumY;
-	CHECKFLOATVAL(numeratorXY, isinf(sumXY) || isinf(sumX) ||
-				  isinf(sumY), true);
+	check_float8_val(numeratorXY, isinf(sumXY) || isinf(sumX) ||
+					 isinf(sumY), true);
 	if (numeratorX <= 0)
 		PG_RETURN_NULL();
 
 	PG_RETURN_FLOAT8(numeratorXY / numeratorX);
 }
 
 Datum
 float8_regr_intercept(PG_FUNCTION_ARGS)
 {
 	ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
@@ -3236,24 +2957,24 @@ float8_regr_intercept(PG_FUNCTION_ARGS)
 	sumX = transvalues[1];
 	sumX2 = transvalues[2];
 	sumY = transvalues[3];
 	sumXY = transvalues[5];
 
 	/* if N is 0 we should return NULL */
 	if (N < 1.0)
 		PG_RETURN_NULL();
 
 	numeratorX = N * sumX2 - sumX * sumX;
-	CHECKFLOATVAL(numeratorX, isinf(sumX2) || isinf(sumX), true);
+	check_float8_val(numeratorX, isinf(sumX2) || isinf(sumX), true);
 	numeratorXXY = sumY * sumX2 - sumX * sumXY;
-	CHECKFLOATVAL(numeratorXXY, isinf(sumY) || isinf(sumX2) ||
-				  isinf(sumX) || isinf(sumXY), true);
+	check_float8_val(numeratorXXY, isinf(sumY) || isinf(sumX2) ||
+					 isinf(sumX) || isinf(sumXY), true);
 	if (numeratorX <= 0)
 		PG_RETURN_NULL();
 
 	PG_RETURN_FLOAT8(numeratorXXY / numeratorX);
 }
 
 
 /*
  *		====================================
  *		MIXED-PRECISION ARITHMETIC OPERATORS
@@ -3264,251 +2985,211 @@ float8_regr_intercept(PG_FUNCTION_ARGS)
  *		float48pl		- returns arg1 + arg2
  *		float48mi		- returns arg1 - arg2
  *		float48mul		- returns arg1 * arg2
  *		float48div		- returns arg1 / arg2
  */
 Datum
 float48pl(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
-	float8		result;
 
-	result = arg1 + arg2;
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
-	PG_RETURN_FLOAT8(result);
+	PG_RETURN_FLOAT8(float8_pl((float8) arg1, arg2));
 }
 
 Datum
 float48mi(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
-	float8		result;
 
-	result = arg1 - arg2;
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
-	PG_RETURN_FLOAT8(result);
+	PG_RETURN_FLOAT8(float8_mi((float8) arg1, arg2));
 }
 
 Datum
 float48mul(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
-	float8		result;
 
-	result = arg1 * arg2;
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2),
-				  arg1 == 0 || arg2 == 0);
-	PG_RETURN_FLOAT8(result);
+	PG_RETURN_FLOAT8(float8_mul((float8) arg1, arg2));
 }
 
 Datum
 float48div(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
-	float8		result;
 
-	if (arg2 == 0.0)
-		ereport(ERROR,
-				(errcode(ERRCODE_DIVISION_BY_ZERO),
-				 errmsg("division by zero")));
-
-	result = arg1 / arg2;
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
-	PG_RETURN_FLOAT8(result);
+	PG_RETURN_FLOAT8(float8_div((float8) arg1, arg2));
 }
 
 /*
  *		float84pl		- returns arg1 + arg2
  *		float84mi		- returns arg1 - arg2
  *		float84mul		- returns arg1 * arg2
  *		float84div		- returns arg1 / arg2
  */
 Datum
 float84pl(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
-	float8		result;
 
-	result = arg1 + arg2;
-
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
-	PG_RETURN_FLOAT8(result);
+	PG_RETURN_FLOAT8(float8_pl(arg1, (float8) arg2));
 }
 
 Datum
 float84mi(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
-	float8		result;
 
-	result = arg1 - arg2;
-
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
-	PG_RETURN_FLOAT8(result);
+	PG_RETURN_FLOAT8(float8_mi(arg1, (float8) arg2));
 }
 
 Datum
 float84mul(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
-	float8		result;
 
-	result = arg1 * arg2;
-
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2),
-				  arg1 == 0 || arg2 == 0);
-	PG_RETURN_FLOAT8(result);
+	PG_RETURN_FLOAT8(float8_mul(arg1, (float8) arg2));
 }
 
 Datum
 float84div(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
-	float8		result;
 
-	if (arg2 == 0.0)
-		ereport(ERROR,
-				(errcode(ERRCODE_DIVISION_BY_ZERO),
-				 errmsg("division by zero")));
-
-	result = arg1 / arg2;
-
-	CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
-	PG_RETURN_FLOAT8(result);
+	PG_RETURN_FLOAT8(float8_div(arg1, (float8) arg2));
 }
 
 /*
  *		====================
  *		COMPARISON OPERATORS
  *		====================
  */
 
 /*
  *		float48{eq,ne,lt,le,gt,ge}		- float4/float8 comparison operations
  */
 Datum
 float48eq(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) == 0);
+	PG_RETURN_BOOL(float8_eq((float8) arg1, arg2));
 }
 
 Datum
 float48ne(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) != 0);
+	PG_RETURN_BOOL(float8_ne((float8) arg1, arg2));
 }
 
 Datum
 float48lt(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) < 0);
+	PG_RETURN_BOOL(float8_lt((float8) arg1, arg2));
 }
 
 Datum
 float48le(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) <= 0);
+	PG_RETURN_BOOL(float8_le((float8) arg1, arg2));
 }
 
 Datum
 float48gt(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) > 0);
+	PG_RETURN_BOOL(float8_gt((float8) arg1, arg2));
 }
 
 Datum
 float48ge(PG_FUNCTION_ARGS)
 {
 	float4		arg1 = PG_GETARG_FLOAT4(0);
 	float8		arg2 = PG_GETARG_FLOAT8(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) >= 0);
+	PG_RETURN_BOOL(float8_ge((float8) arg1, arg2));
 }
 
 /*
  *		float84{eq,ne,lt,le,gt,ge}		- float8/float4 comparison operations
  */
 Datum
 float84eq(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) == 0);
+	PG_RETURN_BOOL(float8_eq(arg1, (float8) arg2));
 }
 
 Datum
 float84ne(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) != 0);
+	PG_RETURN_BOOL(float8_ne(arg1, (float8) arg2));
 }
 
 Datum
 float84lt(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) < 0);
+	PG_RETURN_BOOL(float8_lt(arg1, (float8) arg2));
 }
 
 Datum
 float84le(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) <= 0);
+	PG_RETURN_BOOL(float8_le(arg1, (float8) arg2));
 }
 
 Datum
 float84gt(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) > 0);
+	PG_RETURN_BOOL(float8_gt(arg1, (float8) arg2));
 }
 
 Datum
 float84ge(PG_FUNCTION_ARGS)
 {
 	float8		arg1 = PG_GETARG_FLOAT8(0);
 	float4		arg2 = PG_GETARG_FLOAT4(1);
 
-	PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) >= 0);
+	PG_RETURN_BOOL(float8_ge(arg1, (float8) arg2));
 }
 
 /*
  * Implements the float8 version of the width_bucket() function
  * defined by SQL2003. See also width_bucket_numeric().
  *
  * 'bound1' and 'bound2' are the lower and upper bounds of the
  * histogram's range, respectively. 'count' is the number of buckets
  * in the histogram. width_bucket() returns an integer indicating the
  * bucket number that 'operand' belongs to in an equiwidth histogram
diff --git a/src/backend/utils/adt/formatting.c b/src/backend/utils/adt/formatting.c
index 0e30810ae4..593da9f481 100644
--- a/src/backend/utils/adt/formatting.c
+++ b/src/backend/utils/adt/formatting.c
@@ -84,20 +84,21 @@
 
 #ifdef USE_ICU
 #include <unicode/ustring.h>
 #endif
 
 #include "catalog/pg_collation.h"
 #include "mb/pg_wchar.h"
 #include "utils/builtins.h"
 #include "utils/date.h"
 #include "utils/datetime.h"
+#include "utils/float.h"
 #include "utils/formatting.h"
 #include "utils/int8.h"
 #include "utils/numeric.h"
 #include "utils/pg_locale.h"
 
 /* ----------
  * Routines type
  * ----------
  */
 #define DCH_TYPE		1		/* DATE-TIME version	*/
diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index fc5b68d757..16bda3087c 100644
--- a/src/backend/utils/adt/geo_ops.c
+++ b/src/backend/utils/adt/geo_ops.c
@@ -14,27 +14,23 @@
  */
 #include "postgres.h"
 
 #include <math.h>
 #include <limits.h>
 #include <float.h>
 #include <ctype.h>
 
 #include "libpq/pqformat.h"
 #include "miscadmin.h"
-#include "utils/builtins.h"
+#include "utils/fmgrprotos.h"
 #include "utils/geo_decls.h"
 
-#ifndef M_PI
-#define M_PI 3.14159265358979323846
-#endif
-
 
 /*
  * Internal routines
  */
 
 enum path_delim
 {
 	PATH_NONE, PATH_OPEN, PATH_CLOSED
 };
 
diff --git a/src/backend/utils/adt/geo_spgist.c b/src/backend/utils/adt/geo_spgist.c
index 0fe40499e0..096e59f817 100644
--- a/src/backend/utils/adt/geo_spgist.c
+++ b/src/backend/utils/adt/geo_spgist.c
@@ -69,21 +69,21 @@
  *			src/backend/utils/adt/geo_spgist.c
  *
  *-------------------------------------------------------------------------
  */
 
 #include "postgres.h"
 
 #include "access/spgist.h"
 #include "access/stratnum.h"
 #include "catalog/pg_type.h"
-#include "utils/builtins.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
  * of the index, not the correctness.
  */
 static int
diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
index a1792f0b01..5fbfdaa9de 100644
--- a/src/backend/utils/adt/numeric.c
+++ b/src/backend/utils/adt/numeric.c
@@ -29,20 +29,21 @@
 #include "access/hash.h"
 #include "catalog/pg_type.h"
 #include "common/int.h"
 #include "funcapi.h"
 #include "lib/hyperloglog.h"
 #include "libpq/pqformat.h"
 #include "miscadmin.h"
 #include "nodes/nodeFuncs.h"
 #include "utils/array.h"
 #include "utils/builtins.h"
+#include "utils/float.h"
 #include "utils/guc.h"
 #include "utils/int8.h"
 #include "utils/numeric.h"
 #include "utils/sortsupport.h"
 
 /* ----------
  * Uncomment the following to enable compilation of dump_numeric()
  * and dump_var() and to get a dump of any result produced by make_result().
  * ----------
 #define NUMERIC_DEBUG
diff --git a/src/backend/utils/adt/rangetypes_gist.c b/src/backend/utils/adt/rangetypes_gist.c
index 2fe61043d9..450479e31c 100644
--- a/src/backend/utils/adt/rangetypes_gist.c
+++ b/src/backend/utils/adt/rangetypes_gist.c
@@ -9,21 +9,22 @@
  *
  * IDENTIFICATION
  *	  src/backend/utils/adt/rangetypes_gist.c
  *
  *-------------------------------------------------------------------------
  */
 #include "postgres.h"
 
 #include "access/gist.h"
 #include "access/stratnum.h"
-#include "utils/builtins.h"
+#include "utils/float.h"
+#include "utils/fmgrprotos.h"
 #include "utils/datum.h"
 #include "utils/rangetypes.h"
 
 
 /*
  * Range class properties used to segregate different classes of ranges in
  * GiST.  Each unique combination of properties is a class.  CLS_EMPTY cannot
  * be combined with anything else.
  */
 #define CLS_NORMAL			0	/* Ordinary finite range (no bits set) */
diff --git a/src/backend/utils/adt/rangetypes_selfuncs.c b/src/backend/utils/adt/rangetypes_selfuncs.c
index 59761ecf34..f9a5117492 100644
--- a/src/backend/utils/adt/rangetypes_selfuncs.c
+++ b/src/backend/utils/adt/rangetypes_selfuncs.c
@@ -14,21 +14,22 @@
  *	  src/backend/utils/adt/rangetypes_selfuncs.c
  *
  *-------------------------------------------------------------------------
  */
 #include "postgres.h"
 
 #include "access/htup_details.h"
 #include "catalog/pg_operator.h"
 #include "catalog/pg_statistic.h"
 #include "catalog/pg_type.h"
-#include "utils/builtins.h"
+#include "utils/float.h"
+#include "utils/fmgrprotos.h"
 #include "utils/lsyscache.h"
 #include "utils/rangetypes.h"
 #include "utils/selfuncs.h"
 #include "utils/typcache.h"
 
 static double calc_rangesel(TypeCacheEntry *typcache, VariableStatData *vardata,
 			  RangeType *constval, Oid operator);
 static double default_range_selectivity(Oid operator);
 static double calc_hist_selectivity(TypeCacheEntry *typcache,
 					  VariableStatData *vardata, RangeType *constval,
diff --git a/src/backend/utils/adt/rangetypes_typanalyze.c b/src/backend/utils/adt/rangetypes_typanalyze.c
index f2c11f54bc..9c50e4c1be 100644
--- a/src/backend/utils/adt/rangetypes_typanalyze.c
+++ b/src/backend/utils/adt/rangetypes_typanalyze.c
@@ -19,21 +19,22 @@
  *
  * IDENTIFICATION
  *	  src/backend/utils/adt/rangetypes_typanalyze.c
  *
  *-------------------------------------------------------------------------
  */
 #include "postgres.h"
 
 #include "catalog/pg_operator.h"
 #include "commands/vacuum.h"
-#include "utils/builtins.h"
+#include "utils/float.h"
+#include "utils/fmgrprotos.h"
 #include "utils/lsyscache.h"
 #include "utils/rangetypes.h"
 
 static int	float8_qsort_cmp(const void *a1, const void *a2);
 static int	range_bound_qsort_cmp(const void *a1, const void *a2, void *arg);
 static void compute_range_stats(VacAttrStats *stats,
 					AnalyzeAttrFetchFunc fetchfunc, int samplerows, double totalrows);
 
 /*
  * range_typanalyze -- typanalyze function for range columns
diff --git a/src/backend/utils/adt/timestamp.c b/src/backend/utils/adt/timestamp.c
index e6a1eed191..f8ffa8aa4f 100644
--- a/src/backend/utils/adt/timestamp.c
+++ b/src/backend/utils/adt/timestamp.c
@@ -27,20 +27,21 @@
 #include "common/int128.h"
 #include "funcapi.h"
 #include "libpq/pqformat.h"
 #include "miscadmin.h"
 #include "nodes/makefuncs.h"
 #include "nodes/nodeFuncs.h"
 #include "parser/scansup.h"
 #include "utils/array.h"
 #include "utils/builtins.h"
 #include "utils/datetime.h"
+#include "utils/float.h"
 
 /*
  * gcc's -ffast-math switch breaks routines that expect exact results from
  * expressions like timeval / SECS_PER_HOUR, where timeval is double.
  */
 #ifdef __FAST_MATH__
 #error -ffast-math is known to break this code
 #endif
 
 #define SAMESIGN(a,b)	(((a) < 0) == ((b) < 0))
diff --git a/src/backend/utils/misc/guc.c b/src/backend/utils/misc/guc.c
index 72f6be329e..9530a53ecb 100644
--- a/src/backend/utils/misc/guc.c
+++ b/src/backend/utils/misc/guc.c
@@ -72,20 +72,21 @@
 #include "storage/fd.h"
 #include "storage/large_object.h"
 #include "storage/pg_shmem.h"
 #include "storage/proc.h"
 #include "storage/predicate.h"
 #include "tcop/tcopprot.h"
 #include "tsearch/ts_cache.h"
 #include "utils/builtins.h"
 #include "utils/bytea.h"
 #include "utils/guc_tables.h"
+#include "utils/float.h"
 #include "utils/memutils.h"
 #include "utils/pg_locale.h"
 #include "utils/plancache.h"
 #include "utils/portal.h"
 #include "utils/ps_status.h"
 #include "utils/rls.h"
 #include "utils/snapmgr.h"
 #include "utils/tzparser.h"
 #include "utils/varlena.h"
 #include "utils/xml.h"
diff --git a/src/include/utils/builtins.h b/src/include/utils/builtins.h
index 8bb57c5829..6432222a76 100644
--- a/src/include/utils/builtins.h
+++ b/src/include/utils/builtins.h
@@ -43,34 +43,20 @@ extern int	namestrcmp(Name name, const char *str);
 
 /* numutils.c */
 extern int32 pg_atoi(const char *s, int size, int c);
 extern void pg_itoa(int16 i, char *a);
 extern void pg_ltoa(int32 l, char *a);
 extern void pg_lltoa(int64 ll, char *a);
 extern char *pg_ltostr_zeropad(char *str, int32 value, int32 minwidth);
 extern char *pg_ltostr(char *str, int32 value);
 extern uint64 pg_strtouint64(const char *str, char **endptr, int base);
 
-/* float.c */
-extern PGDLLIMPORT int extra_float_digits;
-
-extern double get_float8_infinity(void);
-extern float get_float4_infinity(void);
-extern double get_float8_nan(void);
-extern float get_float4_nan(void);
-extern int	is_infinite(double val);
-extern double float8in_internal(char *num, char **endptr_p,
-				  const char *type_name, const char *orig_string);
-extern char *float8out_internal(double num);
-extern int	float4_cmp_internal(float4 a, float4 b);
-extern int	float8_cmp_internal(float8 a, float8 b);
-
 /* oid.c */
 extern oidvector *buildoidvector(const Oid *oids, int n);
 extern Oid	oidparse(Node *node);
 extern int	oid_cmp(const void *p1, const void *p2);
 
 /* regexp.c */
 extern char *regexp_fixed_prefix(text *text_re, bool case_insensitive,
 					Oid collation, bool *exact);
 
 /* ruleutils.c */
diff --git a/src/include/utils/float.h b/src/include/utils/float.h
new file mode 100644
index 0000000000..c9f5fd30ca
--- /dev/null
+++ b/src/include/utils/float.h
@@ -0,0 +1,376 @@
+/*-------------------------------------------------------------------------
+ *
+ * float.h
+ *	  Definitions for the built-in floating-point types
+ *
+ * Portions Copyright (c) 1996-2016, PostgreSQL Global Development Group
+ * Portions Copyright (c) 1994, Regents of the University of California
+ *
+ *
+ * IDENTIFICATION
+ *	  src/include/utils/float.c
+ *
+ *-------------------------------------------------------------------------
+ */
+#ifndef FLOAT_H
+#define FLOAT_H
+
+#include <math.h>
+
+#ifndef M_PI
+/* From my RH5.2 gcc math.h file - thomas 2000-04-03 */
+#define M_PI 3.14159265358979323846
+#endif
+
+/* Radians per degree, a.k.a. PI / 180 */
+#define RADIANS_PER_DEGREE 0.0174532925199432957692
+
+/* Visual C++ etc lacks NAN, and won't accept 0.0/0.0.  NAN definition from
+ * http://msdn.microsoft.com/library/default.asp?url=/library/en-us/vclang/html/vclrfNotNumberNANItems.asp
+ */
+#if defined(WIN32) && !defined(NAN)
+static const uint32 nan[2] = {0xffffffff, 0x7fffffff};
+
+#define NAN (*(const float8 *) nan)
+#endif
+
+extern PGDLLIMPORT int extra_float_digits;
+
+/*
+ * Utility functions in float.c
+ */
+extern int	is_infinite(float8 val);
+extern float8 float8in_internal(char *num, char **endptr_p,
+				  const char *type_name, const char *orig_string);
+extern char *float8out_internal(float8 num);
+extern int	float4_cmp_internal(float4 a, float4 b);
+extern int	float8_cmp_internal(float8 a, float8 b);
+
+/*
+ * Routines to provide reasonably platform-independent handling of
+ * infinity and NaN
+ *
+ * We assume that isinf() and isnan() are available and work per spec.
+ * (On some platforms, we have to supply our own; see src/port.)  However,
+ * generating an Infinity or NaN in the first place is less well standardized;
+ * pre-C99 systems tend not to have C99's INFINITY and NAN macros.  We
+ * centralize our workarounds for this here.
+ */
+
+/*
+ * The funny placements of the two #pragmas is necessary because of a
+ * long lived bug in the Microsoft compilers.
+ * See http://support.microsoft.com/kb/120968/en-us for details
+ */
+#if (_MSC_VER >= 1800)
+#pragma warning(disable:4756)
+#endif
+static inline float
+get_float4_infinity(void)
+{
+#ifdef INFINITY
+	/* C99 standard way */
+	return (float) INFINITY;
+#else
+#if (_MSC_VER >= 1800)
+#pragma warning(default:4756)
+#endif
+
+	/*
+	 * On some platforms, HUGE_VAL is an infinity, elsewhere it's just the
+	 * largest normal float8.  We assume forcing an overflow will get us a
+	 * true infinity.
+	 */
+	return (float) (HUGE_VAL * HUGE_VAL);
+#endif
+}
+
+static inline float8
+get_float8_infinity(void)
+{
+#ifdef INFINITY
+	/* C99 standard way */
+	return (float8) INFINITY;
+#else
+
+	/*
+	 * On some platforms, HUGE_VAL is an infinity, elsewhere it's just the
+	 * largest normal float8.  We assume forcing an overflow will get us a
+	 * true infinity.
+	 */
+	return (float8) (HUGE_VAL * HUGE_VAL);
+#endif
+}
+
+static inline float4
+get_float4_nan(void)
+{
+#ifdef NAN
+	/* C99 standard way */
+	return (float) NAN;
+#else
+	/* Assume we can get a NAN via zero divide */
+	return (float) (0.0 / 0.0);
+#endif
+}
+
+static inline float8
+get_float8_nan(void)
+{
+	/* (float8) NAN doesn't work on some NetBSD/MIPS releases */
+#if defined(NAN) && !(defined(__NetBSD__) && defined(__mips__))
+	/* C99 standard way */
+	return (float8) NAN;
+#else
+	/* Assume we can get a NAN via zero divide */
+	return (float8) (0.0 / 0.0);
+#endif
+}
+
+/*
+ * Checks to see if a float4/8 val has underflowed or overflowed
+ */
+
+static inline void
+check_float4_val(float4 val, bool inf_is_valid, bool zero_is_valid)
+{
+	if (isinf(val) && !inf_is_valid)
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("value out of range: overflow")));
+
+	if (val == 0.0 && !zero_is_valid)
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("value out of range: underflow")));
+}
+
+static inline void
+check_float8_val(float8 val, bool inf_is_valid, bool zero_is_valid)
+{
+	if (isinf(val) && !inf_is_valid)
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("value out of range: overflow")));
+
+	if (val == 0.0 && !zero_is_valid)
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("value out of range: underflow")));
+}
+
+/*
+ * Routines for operations with the checks above
+ *
+ * There isn't any way to check for underflow of addition/subtraction
+ * because numbers near the underflow value have already been rounded to
+ * the point where we can't detect that the two values were originally
+ * different, e.g. on x86, '1e-45'::float4 == '2e-45'::float4 ==
+ * 1.4013e-45.
+ */
+
+static inline float4
+float4_pl(float4 val1, float4 val2)
+{
+	float4		result;
+
+	result = val1 + val2;
+	check_float4_val(result, isinf(val1) || isinf(val2), true);
+
+	return result;
+}
+
+static inline float8
+float8_pl(float8 val1, float8 val2)
+{
+	float8		result;
+
+	result = val1 + val2;
+	check_float8_val(result, isinf(val1) || isinf(val2), true);
+
+	return result;
+}
+
+static inline float4
+float4_mi(float4 val1, float4 val2)
+{
+	float4		result;
+
+	result = val1 - val2;
+	check_float4_val(result, isinf(val1) || isinf(val2), true);
+
+	return result;
+}
+
+static inline float8
+float8_mi(float8 val1, float8 val2)
+{
+	float8		result;
+
+	result = val1 - val2;
+	check_float8_val(result, isinf(val1) || isinf(val2), true);
+
+	return result;
+}
+
+static inline float4
+float4_mul(float4 val1, float4 val2)
+{
+	float4		result;
+
+	result = val1 * val2;
+	check_float4_val(result, isinf(val1) || isinf(val2),
+					 val1 == 0.0f || val2 == 0.0f);
+
+	return result;
+}
+
+static inline float8
+float8_mul(float8 val1, float8 val2)
+{
+	float8		result;
+
+	result = val1 * val2;
+	check_float8_val(result, isinf(val1) || isinf(val2),
+					 val1 == 0.0 || val2 == 0.0);
+
+	return result;
+}
+
+static inline float4
+float4_div(float4 val1, float4 val2)
+{
+	float4		result;
+
+	if (val2 == 0.0f)
+		ereport(ERROR,
+				(errcode(ERRCODE_DIVISION_BY_ZERO),
+				 errmsg("division by zero")));
+
+	result = val1 / val2;
+	check_float4_val(result, isinf(val1) || isinf(val2), val1 == 0.0f);
+
+	return result;
+}
+
+static inline float8
+float8_div(float8 val1, float8 val2)
+{
+	float8		result;
+
+	if (val2 == 0.0)
+		ereport(ERROR,
+				(errcode(ERRCODE_DIVISION_BY_ZERO),
+				 errmsg("division by zero")));
+
+	result = val1 / val2;
+	check_float8_val(result, isinf(val1) || isinf(val2), val1 == 0.0);
+
+	return result;
+}
+
+/*
+ * Routines for NaN-aware comparisons
+ *
+ * We consider all NANs to be equal and larger than any non-NAN. This is
+ * somewhat arbitrary; the important thing is to have a consistent sort
+ * order.
+ */
+
+static inline bool
+float4_eq(float4 val1, float4 val2)
+{
+	return isnan(val1) ? isnan(val2) : !isnan(val2) && val1 == val2;
+}
+
+static inline bool
+float8_eq(float8 val1, float8 val2)
+{
+	return isnan(val1) ? isnan(val2) : !isnan(val2) && val1 == val2;
+}
+
+static inline bool
+float4_ne(float4 val1, float4 val2)
+{
+	return isnan(val1) ? !isnan(val2) : isnan(val2) || val1 != val2;
+}
+
+static inline bool
+float8_ne(float8 val1, float8 val2)
+{
+	return isnan(val1) ? !isnan(val2) : isnan(val2) || val1 != val2;
+}
+
+static inline bool
+float4_lt(float4 val1, float4 val2)
+{
+	return !isnan(val1) && (isnan(val2) || val1 < val2);
+}
+
+static inline bool
+float8_lt(float8 val1, float8 val2)
+{
+	return !isnan(val1) && (isnan(val2) || val1 < val2);
+}
+
+static inline bool
+float4_le(float4 val1, float4 val2)
+{
+	return isnan(val2) || (!isnan(val1) && val1 <= val2);
+}
+
+static inline bool
+float8_le(float8 val1, float8 val2)
+{
+	return isnan(val2) || (!isnan(val1) && val1 <= val2);
+}
+
+static inline bool
+float4_gt(float4 val1, float4 val2)
+{
+	return !isnan(val2) && (isnan(val1) || val1 > val2);
+}
+
+static inline bool
+float8_gt(float8 val1, float8 val2)
+{
+	return !isnan(val2) && (isnan(val1) || val1 > val2);
+}
+
+static inline bool
+float4_ge(float4 val1, float4 val2)
+{
+	return isnan(val1) || (!isnan(val2) && val1 >= val2);
+}
+
+static inline bool
+float8_ge(float8 val1, float8 val2)
+{
+	return isnan(val1) || (!isnan(val2) && val1 >= val2);
+}
+
+static inline float4
+float4_min(float4 val1, float4 val2)
+{
+	return float4_lt(val1, val2) ? val1 : val2;
+}
+
+static inline float8
+float8_min(float8 val1, float8 val2)
+{
+	return float8_lt(val1, val2) ? val1 : val2;
+}
+
+static inline float4
+float4_max(float4 val1, float4 val2)
+{
+	return float4_gt(val1, val2) ? val1 : val2;
+}
+
+static inline float8
+float8_max(float8 val1, float8 val2)
+{
+	return float8_gt(val1, val2) ? val1 : val2;
+}
+
+#endif							/* FLOAT_H */
diff --git a/src/include/utils/geo_decls.h b/src/include/utils/geo_decls.h
index 0e066894cd..aeb31515e4 100644
--- a/src/include/utils/geo_decls.h
+++ b/src/include/utils/geo_decls.h
@@ -17,20 +17,21 @@
  *		and data types. There are still some more to do. - tgl 97/04/19
  *
  *-------------------------------------------------------------------------
  */
 #ifndef GEO_DECLS_H
 #define GEO_DECLS_H
 
 #include <math.h>
 
 #include "fmgr.h"
+#include "utils/float.h"
 
 /*--------------------------------------------------------------------
  * Useful floating point utilities and constants.
  *-------------------------------------------------------------------*/
 
 
 #define EPSILON					1.0E-06
 
 #ifdef EPSILON
 #define FPzero(A)				(fabs(A) <= EPSILON)
-- 
2.14.3 (Apple Git-98)

