From cd5678df60bd2a290264a2c5e5b3f7f4201d00f0 Mon Sep 17 00:00:00 2001
From: Tomas Vondra <tomas@2ndquadrant.com>
Date: Thu, 26 Jul 2018 14:39:29 +0200
Subject: [PATCH 2/6] 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.  The API uses inline functions
instead of macros, as they are easier to use compared to macros,
and avoid double-evaluation hazards.

Discussion: https://www.postgresql.org/message-id/CAE2gYzxF7-5djV6-cEvqQu-fNsnt%3DEqbOURx7ZDg%2BVv6ZMTWbg%40mail.gmail.com
---
 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            |  56 ++-
 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               |   7 +-
 src/backend/utils/adt/geo_spgist.c            |   3 +-
 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 ++++++++++++++++
 20 files changed, 556 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 b6d22d2b00..2ecf7a2d87 100644
--- a/contrib/btree_gin/btree_gin.c
+++ b/contrib/btree_gin/btree_gin.c
@@ -10,6 +10,7 @@
 #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"
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
@@ -9,6 +9,7 @@
 #include "btree_utils_num.h"
 #include "utils/builtins.h"
 #include "utils/datetime.h"
+#include "utils/float.h"
 
 typedef struct
 {
diff --git a/contrib/cube/cube.c b/contrib/cube/cube.c
index f02ac24ea1..dfa8465d74 100644
--- a/contrib/cube/cube.c
+++ b/contrib/cube/cube.c
@@ -13,7 +13,7 @@
 #include "access/gist.h"
 #include "access/stratnum.h"
 #include "utils/array.h"
-#include "utils/builtins.h"
+#include "utils/float.h"
 
 #include "cubedata.h"
 
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
@@ -7,7 +7,7 @@
 #include "postgres.h"
 
 #include "cubedata.h"
-#include "utils/builtins.h"
+#include "utils/float.h"
 
 /* All grammar constructs return strings */
 #define YYSTYPE char *
diff --git a/contrib/postgres_fdw/postgres_fdw.c b/contrib/postgres_fdw/postgres_fdw.c
index 5699252091..0803c30a48 100644
--- a/contrib/postgres_fdw/postgres_fdw.c
+++ b/contrib/postgres_fdw/postgres_fdw.c
@@ -35,6 +35,7 @@
 #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"
diff --git a/src/backend/access/gist/gistget.c b/src/backend/access/gist/gistget.c
index c4e8a3b913..ad07b9e63c 100644
--- a/src/backend/access/gist/gistget.c
+++ b/src/backend/access/gist/gistget.c
@@ -22,7 +22,7 @@
 #include "storage/predicate.h"
 #include "pgstat.h"
 #include "lib/pairingheap.h"
-#include "utils/builtins.h"
+#include "utils/float.h"
 #include "utils/memutils.h"
 #include "utils/rel.h"
 
diff --git a/src/backend/access/gist/gistproc.c b/src/backend/access/gist/gistproc.c
index 0536b318cc..d6ce5ccf6c 100644
--- a/src/backend/access/gist/gistproc.c
+++ b/src/backend/access/gist/gistproc.c
@@ -22,6 +22,7 @@
 #include "access/gist.h"
 #include "access/stratnum.h"
 #include "utils/builtins.h"
+#include "utils/float.h"
 #include "utils/geo_decls.h"
 
 
@@ -33,15 +34,6 @@ static bool rtree_internal_consistent(BOX *key, BOX *query,
 /* 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
@@ -53,10 +45,10 @@ static bool rtree_internal_consistent(BOX *key, BOX *query,
 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);
 }
 
 /*
@@ -73,8 +65,8 @@ size_box(const BOX *box)
 	 *
 	 * 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;
 
 	/*
@@ -143,13 +135,13 @@ gist_box_consistent(PG_FUNCTION_ARGS)
 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;
 }
 
@@ -615,9 +607,9 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
 			 * 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++;
 			}
@@ -630,7 +622,7 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
 			 * left group.
 			 */
 			while (i2 < nentries &&
-				   FLOAT8_LE(intervalsUpper[i2].upper, leftUpper))
+				   float8_le(intervalsUpper[i2].upper, leftUpper))
 				i2++;
 
 			/*
@@ -652,9 +644,9 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
 			/*
 			 * 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--;
 			}
@@ -666,7 +658,7 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
 			 * 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--;
 
 			/*
@@ -754,10 +746,10 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
 			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;
@@ -775,7 +767,7 @@ gist_box_picksplit(PG_FUNCTION_ARGS)
 			 * 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);
@@ -859,10 +851,10 @@ gist_box_same(PG_FUNCTION_ARGS)
 	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);
diff --git a/src/backend/access/gist/gistutil.c b/src/backend/access/gist/gistutil.c
index 12804c321c..dddfe0ae2c 100644
--- a/src/backend/access/gist/gistutil.c
+++ b/src/backend/access/gist/gistutil.c
@@ -21,7 +21,7 @@
 #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"
 
 
diff --git a/src/backend/utils/adt/float.c b/src/backend/utils/adt/float.c
index b86205b098..df35557b73 100644
--- a/src/backend/utils/adt/float.c
+++ b/src/backend/utils/adt/float.c
@@ -23,44 +23,11 @@
 #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
-
-/*
- * 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 */
 
@@ -106,86 +73,6 @@ static double cbrt(double x);
 
 
 /*
- * 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
@@ -343,7 +230,7 @@ float4in(PG_FUNCTION_ARGS)
 	 * 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);
 }
@@ -693,7 +580,7 @@ float4larger(PG_FUNCTION_ARGS)
 	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;
@@ -707,7 +594,7 @@ float4smaller(PG_FUNCTION_ARGS)
 	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;
@@ -760,7 +647,7 @@ float8larger(PG_FUNCTION_ARGS)
 	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;
@@ -774,7 +661,7 @@ float8smaller(PG_FUNCTION_ARGS)
 	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;
@@ -799,19 +686,8 @@ 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
@@ -819,11 +695,8 @@ 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
@@ -831,12 +704,8 @@ 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
@@ -844,17 +713,8 @@ 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));
 }
 
 /*
@@ -868,12 +728,8 @@ 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
@@ -881,12 +737,8 @@ 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
@@ -894,13 +746,8 @@ 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
@@ -908,17 +755,8 @@ 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));
 }
 
 
@@ -934,31 +772,11 @@ float8div(PG_FUNCTION_ARGS)
 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
@@ -967,7 +785,7 @@ 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
@@ -976,7 +794,7 @@ 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
@@ -985,7 +803,7 @@ 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
@@ -994,7 +812,7 @@ 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
@@ -1003,7 +821,7 @@ 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
@@ -1012,7 +830,7 @@ 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
@@ -1048,31 +866,11 @@ btfloat4sortsupport(PG_FUNCTION_ARGS)
 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
@@ -1081,7 +879,7 @@ 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
@@ -1090,7 +888,7 @@ 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
@@ -1099,7 +897,7 @@ 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
@@ -1108,7 +906,7 @@ 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
@@ -1117,7 +915,7 @@ 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
@@ -1126,7 +924,7 @@ 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
@@ -1341,7 +1139,7 @@ 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);
 }
@@ -1566,7 +1364,7 @@ dsqrt(PG_FUNCTION_ARGS)
 
 	result = sqrt(arg1);
 
-	CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
+	check_float8_val(result, isinf(arg1), arg1 == 0);
 	PG_RETURN_FLOAT8(result);
 }
 
@@ -1581,7 +1379,7 @@ dcbrt(PG_FUNCTION_ARGS)
 	float8		result;
 
 	result = cbrt(arg1);
-	CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
+	check_float8_val(result, isinf(arg1), arg1 == 0);
 	PG_RETURN_FLOAT8(result);
 }
 
@@ -1653,7 +1451,7 @@ dpow(PG_FUNCTION_ARGS)
 	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);
 }
 
@@ -1672,7 +1470,7 @@ dexp(PG_FUNCTION_ARGS)
 	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);
 }
 
@@ -1701,7 +1499,7 @@ dlog1(PG_FUNCTION_ARGS)
 
 	result = log(arg1);
 
-	CHECKFLOATVAL(result, isinf(arg1), arg1 == 1);
+	check_float8_val(result, isinf(arg1), arg1 == 1);
 	PG_RETURN_FLOAT8(result);
 }
 
@@ -1731,7 +1529,7 @@ dlog10(PG_FUNCTION_ARGS)
 
 	result = log10(arg1);
 
-	CHECKFLOATVAL(result, isinf(arg1), arg1 == 1);
+	check_float8_val(result, isinf(arg1), arg1 == 1);
 	PG_RETURN_FLOAT8(result);
 }
 
@@ -1761,7 +1559,7 @@ dacos(PG_FUNCTION_ARGS)
 
 	result = acos(arg1);
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
@@ -1791,7 +1589,7 @@ dasin(PG_FUNCTION_ARGS)
 
 	result = asin(arg1);
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
@@ -1816,7 +1614,7 @@ datan(PG_FUNCTION_ARGS)
 	 */
 	result = atan(arg1);
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
@@ -1841,7 +1639,7 @@ datan2(PG_FUNCTION_ARGS)
 	 */
 	result = atan2(arg1, arg2);
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
@@ -1881,7 +1679,7 @@ dcos(PG_FUNCTION_ARGS)
 				(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);
 }
 
@@ -1908,7 +1706,7 @@ dcot(PG_FUNCTION_ARGS)
 				 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);
 }
 
@@ -1934,7 +1732,7 @@ dsin(PG_FUNCTION_ARGS)
 				(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);
 }
 
@@ -1960,7 +1758,7 @@ dtan(PG_FUNCTION_ARGS)
 				(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);
 }
 
@@ -2112,7 +1910,7 @@ dacosd(PG_FUNCTION_ARGS)
 	else
 		result = 90.0 + asind_q1(-arg1);
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
@@ -2147,7 +1945,7 @@ dasind(PG_FUNCTION_ARGS)
 	else
 		result = -asind_q1(-arg1);
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
@@ -2177,7 +1975,7 @@ datand(PG_FUNCTION_ARGS)
 	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);
 }
 
@@ -2211,7 +2009,7 @@ datan2d(PG_FUNCTION_ARGS)
 	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);
 }
 
@@ -2332,7 +2130,7 @@ dcosd(PG_FUNCTION_ARGS)
 
 	result = sign * cosd_q1(arg1);
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
@@ -2397,7 +2195,7 @@ dcotd(PG_FUNCTION_ARGS)
 	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);
 }
 
@@ -2451,7 +2249,7 @@ dsind(PG_FUNCTION_ARGS)
 
 	result = sign * sind_q1(arg1);
 
-	CHECKFLOATVAL(result, false, true);
+	check_float8_val(result, false, true);
 	PG_RETURN_FLOAT8(result);
 }
 
@@ -2516,7 +2314,7 @@ dtand(PG_FUNCTION_ARGS)
 	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);
 }
 
@@ -2528,12 +2326,8 @@ 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));
 }
 
 
@@ -2554,12 +2348,8 @@ 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));
 }
 
 
@@ -2651,31 +2441,16 @@ float8_combine(PG_FUNCTION_ARGS)
 	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);
 }
@@ -2686,20 +2461,8 @@ 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
@@ -2708,9 +2471,9 @@ float8_accum(PG_FUNCTION_ARGS)
 	 */
 	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);
 	}
@@ -2719,9 +2482,9 @@ float8_accum(PG_FUNCTION_ARGS)
 		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,
@@ -2739,20 +2502,8 @@ float4_accum(PG_FUNCTION_ARGS)
 	/* 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
@@ -2761,9 +2512,9 @@ float4_accum(PG_FUNCTION_ARGS)
 	 */
 	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);
 	}
@@ -2772,9 +2523,9 @@ float4_accum(PG_FUNCTION_ARGS)
 		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,
@@ -2824,7 +2575,7 @@ float8_var_pop(PG_FUNCTION_ARGS)
 		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)
@@ -2853,7 +2604,7 @@ float8_var_samp(PG_FUNCTION_ARGS)
 		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)
@@ -2882,7 +2633,7 @@ float8_stddev_pop(PG_FUNCTION_ARGS)
 		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)
@@ -2911,7 +2662,7 @@ float8_stddev_samp(PG_FUNCTION_ARGS)
 		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)
@@ -2960,16 +2711,16 @@ float8_regr_accum(PG_FUNCTION_ARGS)
 
 	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
@@ -3022,49 +2773,19 @@ float8_regr_combine(PG_FUNCTION_ARGS)
 	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);
 }
@@ -3090,7 +2811,7 @@ float8_regr_sxx(PG_FUNCTION_ARGS)
 		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)
@@ -3119,7 +2840,7 @@ float8_regr_syy(PG_FUNCTION_ARGS)
 		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)
@@ -3150,8 +2871,8 @@ float8_regr_sxy(PG_FUNCTION_ARGS)
 		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 */
 
@@ -3218,8 +2939,8 @@ float8_covar_pop(PG_FUNCTION_ARGS)
 		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));
 }
@@ -3246,8 +2967,8 @@ float8_covar_samp(PG_FUNCTION_ARGS)
 		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)));
 }
@@ -3280,12 +3001,12 @@ float8_corr(PG_FUNCTION_ARGS)
 		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();
 
@@ -3320,12 +3041,12 @@ float8_regr_r2(PG_FUNCTION_ARGS)
 		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 */
@@ -3361,10 +3082,10 @@ float8_regr_slope(PG_FUNCTION_ARGS)
 		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();
 
@@ -3396,10 +3117,10 @@ float8_regr_intercept(PG_FUNCTION_ARGS)
 		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();
 
@@ -3424,11 +3145,8 @@ 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
@@ -3436,11 +3154,8 @@ 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
@@ -3448,12 +3163,8 @@ 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
@@ -3461,16 +3172,8 @@ 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));
 }
 
 /*
@@ -3484,12 +3187,8 @@ 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
@@ -3497,12 +3196,8 @@ 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
@@ -3510,13 +3205,8 @@ 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
@@ -3524,17 +3214,8 @@ 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));
 }
 
 /*
@@ -3552,7 +3233,7 @@ 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
@@ -3561,7 +3242,7 @@ 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
@@ -3570,7 +3251,7 @@ 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
@@ -3579,7 +3260,7 @@ 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
@@ -3588,7 +3269,7 @@ 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
@@ -3597,7 +3278,7 @@ 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));
 }
 
 /*
@@ -3609,7 +3290,7 @@ 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
@@ -3618,7 +3299,7 @@ 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
@@ -3627,7 +3308,7 @@ 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
@@ -3636,7 +3317,7 @@ 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
@@ -3645,7 +3326,7 @@ 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
@@ -3654,7 +3335,7 @@ 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));
 }
 
 /*
diff --git a/src/backend/utils/adt/formatting.c b/src/backend/utils/adt/formatting.c
index a345c65605..30696e3575 100644
--- a/src/backend/utils/adt/formatting.c
+++ b/src/backend/utils/adt/formatting.c
@@ -91,6 +91,7 @@
 #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"
diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index 3e7519015d..1fb2ff2603 100644
--- a/src/backend/utils/adt/geo_ops.c
+++ b/src/backend/utils/adt/geo_ops.c
@@ -21,13 +21,10 @@
 
 #include "libpq/pqformat.h"
 #include "miscadmin.h"
-#include "utils/builtins.h"
+#include "utils/float.h"
+#include "utils/fmgrprotos.h"
 #include "utils/geo_decls.h"
 
-#ifndef M_PI
-#define M_PI 3.14159265358979323846
-#endif
-
 
 /*
  * Internal routines
diff --git a/src/backend/utils/adt/geo_spgist.c b/src/backend/utils/adt/geo_spgist.c
index d7c807f8d9..fea36f361a 100644
--- a/src/backend/utils/adt/geo_spgist.c
+++ b/src/backend/utils/adt/geo_spgist.c
@@ -76,7 +76,8 @@
 #include "access/spgist.h"
 #include "access/stratnum.h"
 #include "catalog/pg_type.h"
-#include "utils/builtins.h"
+#include "utils/float.h"
+#include "utils/fmgrprotos.h"
 #include "utils/geo_decls.h"
 
 /*
diff --git a/src/backend/utils/adt/numeric.c b/src/backend/utils/adt/numeric.c
index 4d5042b278..444e575e1d 100644
--- a/src/backend/utils/adt/numeric.c
+++ b/src/backend/utils/adt/numeric.c
@@ -36,6 +36,7 @@
 #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"
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
@@ -16,7 +16,8 @@
 
 #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"
 
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
@@ -21,7 +21,8 @@
 #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"
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
@@ -26,7 +26,8 @@
 
 #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"
 
diff --git a/src/backend/utils/adt/timestamp.c b/src/backend/utils/adt/timestamp.c
index 9a481f6eb9..284e14dd73 100644
--- a/src/backend/utils/adt/timestamp.c
+++ b/src/backend/utils/adt/timestamp.c
@@ -33,6 +33,7 @@
 #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
diff --git a/src/backend/utils/misc/guc.c b/src/backend/utils/misc/guc.c
index a88ea6cfc9..c123de1a59 100644
--- a/src/backend/utils/misc/guc.c
+++ b/src/backend/utils/misc/guc.c
@@ -80,6 +80,7 @@
 #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"
diff --git a/src/include/utils/builtins.h b/src/include/utils/builtins.h
index 88a42b345c..6f55699912 100644
--- a/src/include/utils/builtins.h
+++ b/src/include/utils/builtins.h
@@ -52,20 +52,6 @@ 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);
diff --git a/src/include/utils/float.h b/src/include/utils/float.h
new file mode 100644
index 0000000000..0e8483c930
--- /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. */
+#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 float4
+get_float4_infinity(void)
+{
+#ifdef INFINITY
+	/* C99 standard way */
+	return (float4) 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 (float4) (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 (float4) NAN;
+#else
+	/* Assume we can get a NAN via zero divide */
+	return (float4) (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(const float4 val, const bool inf_is_valid,
+				 const bool zero_is_valid)
+{
+	if (!inf_is_valid && unlikely(isinf(val)))
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("value out of range: overflow")));
+
+	if (!zero_is_valid && unlikely(val == 0.0))
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("value out of range: underflow")));
+}
+
+static inline void
+check_float8_val(const float8 val, const bool inf_is_valid,
+				 const bool zero_is_valid)
+{
+	if (!inf_is_valid && unlikely(isinf(val)))
+		ereport(ERROR,
+				(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+				 errmsg("value out of range: overflow")));
+
+	if (!zero_is_valid && unlikely(val == 0.0))
+		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(const float4 val1, const float4 val2)
+{
+	float4		result;
+
+	result = val1 + val2;
+	check_float4_val(result, isinf(val1) || isinf(val2), true);
+
+	return result;
+}
+
+static inline float8
+float8_pl(const float8 val1, const float8 val2)
+{
+	float8		result;
+
+	result = val1 + val2;
+	check_float8_val(result, isinf(val1) || isinf(val2), true);
+
+	return result;
+}
+
+static inline float4
+float4_mi(const float4 val1, const float4 val2)
+{
+	float4		result;
+
+	result = val1 - val2;
+	check_float4_val(result, isinf(val1) || isinf(val2), true);
+
+	return result;
+}
+
+static inline float8
+float8_mi(const float8 val1, const float8 val2)
+{
+	float8		result;
+
+	result = val1 - val2;
+	check_float8_val(result, isinf(val1) || isinf(val2), true);
+
+	return result;
+}
+
+static inline float4
+float4_mul(const float4 val1, const 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(const float8 val1, const 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(const float4 val1, const 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(const float8 val1, const 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(const float4 val1, const float4 val2)
+{
+	return isnan(val1) ? isnan(val2) : !isnan(val2) && val1 == val2;
+}
+
+static inline bool
+float8_eq(const float8 val1, const float8 val2)
+{
+	return isnan(val1) ? isnan(val2) : !isnan(val2) && val1 == val2;
+}
+
+static inline bool
+float4_ne(const float4 val1, const float4 val2)
+{
+	return isnan(val1) ? !isnan(val2) : isnan(val2) || val1 != val2;
+}
+
+static inline bool
+float8_ne(const float8 val1, const float8 val2)
+{
+	return isnan(val1) ? !isnan(val2) : isnan(val2) || val1 != val2;
+}
+
+static inline bool
+float4_lt(const float4 val1, const float4 val2)
+{
+	return !isnan(val1) && (isnan(val2) || val1 < val2);
+}
+
+static inline bool
+float8_lt(const float8 val1, const float8 val2)
+{
+	return !isnan(val1) && (isnan(val2) || val1 < val2);
+}
+
+static inline bool
+float4_le(const float4 val1, const float4 val2)
+{
+	return isnan(val2) || (!isnan(val1) && val1 <= val2);
+}
+
+static inline bool
+float8_le(const float8 val1, const float8 val2)
+{
+	return isnan(val2) || (!isnan(val1) && val1 <= val2);
+}
+
+static inline bool
+float4_gt(const float4 val1, const float4 val2)
+{
+	return !isnan(val2) && (isnan(val1) || val1 > val2);
+}
+
+static inline bool
+float8_gt(const float8 val1, const float8 val2)
+{
+	return !isnan(val2) && (isnan(val1) || val1 > val2);
+}
+
+static inline bool
+float4_ge(const float4 val1, const float4 val2)
+{
+	return isnan(val1) || (!isnan(val2) && val1 >= val2);
+}
+
+static inline bool
+float8_ge(const float8 val1, const float8 val2)
+{
+	return isnan(val1) || (!isnan(val2) && val1 >= val2);
+}
+
+static inline float4
+float4_min(const float4 val1, const float4 val2)
+{
+	return float4_lt(val1, val2) ? val1 : val2;
+}
+
+static inline float8
+float8_min(const float8 val1, const float8 val2)
+{
+	return float8_lt(val1, val2) ? val1 : val2;
+}
+
+static inline float4
+float4_max(const float4 val1, const float4 val2)
+{
+	return float4_gt(val1, val2) ? val1 : val2;
+}
+
+static inline float8
+float8_max(const float8 val1, const float8 val2)
+{
+	return float8_gt(val1, val2) ? val1 : val2;
+}
+
+#endif							/* FLOAT_H */
-- 
2.13.6

