]> granicus.if.org Git - postgresql/blobdiff - src/backend/utils/adt/float.c
Update copyright for 2014
[postgresql] / src / backend / utils / adt / float.c
index 82c257c80157242fb27489db5805751ed3c8c007..e2a0812d194cf2ef3b2eac53f66da786500090c7 100644 (file)
  * float.c
  *       Functions for the built-in floating-point types.
  *
- * Portions Copyright (c) 1996-2003, PostgreSQL Global Development Group
+ * Portions Copyright (c) 1996-2014, PostgreSQL Global Development Group
  * Portions Copyright (c) 1994, Regents of the University of California
  *
  *
  * IDENTIFICATION
- *       $Header: /cvsroot/pgsql/src/backend/utils/adt/float.c,v 1.93 2003/08/04 02:40:04 momjian Exp $
+ *       src/backend/utils/adt/float.c
  *
  *-------------------------------------------------------------------------
  */
-/*----------
- * OLD COMMENTS
- *             Basic float4 ops:
- *              float4in, float4out, float4recv, float4send
- *              float4abs, float4um, float4up
- *             Basic float8 ops:
- *              float8in, float8out, float8recv, float8send
- *              float8abs, float8um, float8up
- *             Arithmetic operators:
- *              float4pl, float4mi, float4mul, float4div
- *              float8pl, float8mi, float8mul, float8div
- *             Comparison operators:
- *              float4eq, float4ne, float4lt, float4le, float4gt, float4ge, float4cmp
- *              float8eq, float8ne, float8lt, float8le, float8gt, float8ge, float8cmp
- *             Conversion routines:
- *              ftod, dtof, i4tod, dtoi4, i2tod, dtoi2, itof, ftoi, i2tof, ftoi2
- *
- *             Random float8 ops:
- *              dround, dtrunc, dsqrt, dcbrt, dpow, dexp, dlog1
- *             Arithmetic operators:
- *              float48pl, float48mi, float48mul, float48div
- *              float84pl, float84mi, float84mul, float84div
- *             Comparison operators:
- *              float48eq, float48ne, float48lt, float48le, float48gt, float48ge
- *              float84eq, float84ne, float84lt, float84le, float84gt, float84ge
- *
- *             (You can do the arithmetic and comparison stuff using conversion
- *              routines, but then you pay the overhead of invoking a separate
- *              conversion function...)
- *
- * XXX GLUESOME STUFF. FIX IT! -AY '94
- *
- *             Added some additional conversion routines and cleaned up
- *              a bit of the existing code. Need to change the error checking
- *              for calls to pow(), exp() since on some machines (my Linux box
- *              included) these routines do not set errno. - tgl 97/05/10
- *----------
- */
 #include "postgres.h"
 
 #include <ctype.h>
-#include <errno.h>
-#include <float.h>                             /* faked on sunos4 */
+#include <float.h>
 #include <math.h>
-
 #include <limits.h>
-/* for finite() on Solaris */
-#ifdef HAVE_IEEEFP_H
-#include <ieeefp.h>
-#endif
 
 #include "catalog/pg_type.h"
-#include "fmgr.h"
 #include "libpq/pqformat.h"
 #include "utils/array.h"
 #include "utils/builtins.h"
+#include "utils/sortsupport.h"
 
 
-#ifndef HAVE_CBRT
-static double cbrt(double x);
-#endif   /* HAVE_CBRT */
-
 #ifndef M_PI
 /* from my RH5.2 gcc math.h file - thomas 2000-04-03 */
 #define M_PI 3.14159265358979323846
 #endif
 
-#ifndef NAN
-#define NAN            (0.0/0.0)
-#endif
+/* 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};
 
-#ifndef SHRT_MAX
-#define SHRT_MAX 32767
-#endif
-#ifndef SHRT_MIN
-#define SHRT_MIN (-32768)
+#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
 
-/* ========== USER I/O ROUTINES ========== */
+/*
+ * 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)
 
 
-#define FLOAT4_MAX              FLT_MAX
-#define FLOAT4_MIN              FLT_MIN
-#define FLOAT8_MAX              DBL_MAX
-#define FLOAT8_MIN              DBL_MIN
+/* ========== USER I/O ROUTINES ========== */
 
 
 /* Configurable GUC parameter */
 int                    extra_float_digits = 0;         /* Added to DBL_DIG or FLT_DIG */
 
 
-static void CheckFloat4Val(double val);
-static void CheckFloat8Val(double val);
 static int     float4_cmp_internal(float4 a, float4 b);
 static int     float8_cmp_internal(float8 a, float8 b);
 
+#ifndef HAVE_CBRT
+/*
+ * Some machines (in particular, some versions of AIX) have an extern
+ * declaration for cbrt() in <math.h> but fail to provide the actual
+ * 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 */
+
 
 /*
- * check to see if a float4 val is outside of
- * the FLOAT4_MIN, FLOAT4_MAX bounds.
- *
- * raise an ereport warning if it is
-*/
-static void
-CheckFloat4Val(double val)
+ * 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
+}
+
+float
+get_float4_infinity(void)
 {
+#ifdef INFINITY
+       /* C99 standard way */
+       return (float) INFINITY;
+#else
+
        /*
-        * defining unsafe floats's will make float4 and float8 ops faster at
-        * the cost of safety, of course!
+        * 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.
         */
-#ifdef UNSAFE_FLOATS
-       return;
+       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
-       if (fabs(val) > FLOAT4_MAX)
-               ereport(ERROR,
-                               (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
-                                errmsg("float4 value out of range: overflow")));
-       if (val != 0.0 && fabs(val) < FLOAT4_MIN)
-               ereport(ERROR,
-                               (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
-                                errmsg("float4 value out of range: underflow")));
+       /* Assume we can get a NAN via zero divide */
+       return (double) (0.0 / 0.0);
+#endif
+}
 
-       return;
-#endif   /* UNSAFE_FLOATS */
+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
 }
 
+
 /*
- * check to see if a float8 val is outside of
- * the FLOAT8_MIN, FLOAT8_MAX bounds.
- *
- * raise an ereport error if it is
+ * 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.
  */
-static void
-CheckFloat8Val(double val)
+int
+is_infinite(double val)
 {
-       /*
-        * defining unsafe floats's will make float4 and float8 ops faster at
-        * the cost of safety, of course!
-        */
-#ifdef UNSAFE_FLOATS
-       return;
-#else
-       if (fabs(val) > FLOAT8_MAX)
-               ereport(ERROR,
-                               (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
-                                errmsg("float8 value out of range: overflow")));
-       if (val != 0.0 && fabs(val) < FLOAT8_MIN)
-               ereport(ERROR,
-                               (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
-                                errmsg("float8 value out of range: underflow")));
-#endif   /* UNSAFE_FLOATS */
+       int                     inf = isinf(val);
+
+       if (inf == 0)
+               return 0;
+       else if (val > 0)
+               return 1;
+       else
+               return -1;
 }
 
+
 /*
- *             float4in                - converts "num" to float
- *                                               restricted syntax:
- *                                               {<sp>} [+|-] {digit} [.{digit}] [<exp>]
- *                                               where <sp> is a space, digit is 0-9,
- *                                               <exp> is "e" or "E" followed by an integer.
+ *             float4in                - converts "num" to float4
  */
 Datum
 float4in(PG_FUNCTION_ARGS)
 {
        char       *num = PG_GETARG_CSTRING(0);
+       char       *orig_num;
        double          val;
        char       *endptr;
 
+       /*
+        * endptr points to the first character _after_ the sequence we recognized
+        * as a valid floating point number. orig_num points to the original input
+        * string.
+        */
+       orig_num = num;
+
+       /* skip leading whitespace */
+       while (*num != '\0' && isspace((unsigned char) *num))
+               num++;
+
+       /*
+        * Check for an empty-string input to begin with, to avoid the vagaries of
+        * strtod() on different platforms.
+        */
+       if (*num == '\0')
+               ereport(ERROR,
+                               (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
+                                errmsg("invalid input syntax for type real: \"%s\"",
+                                               orig_num)));
+
        errno = 0;
        val = strtod(num, &endptr);
-       if (*endptr != '\0')
+
+       /* did we not see anything that looks like a double? */
+       if (endptr == num || errno != 0)
        {
+               int                     save_errno = errno;
+
                /*
-                * XXX we should accept "Infinity" and "-Infinity" too, but what
-                * are the correct values to assign?  HUGE_VAL will provoke an
-                * error from CheckFloat4Val.
+                * C99 requires that strtod() accept NaN, [+-]Infinity, and [+-]Inf,
+                * but not all platforms support all of these (and some accept them
+                * but set ERANGE anyway...)  Therefore, we check for these inputs
+                * ourselves if strtod() fails.
+                *
+                * Note: C99 also requires hexadecimal input as well as some extended
+                * forms of NaN, but we consider these forms unportable and don't try
+                * to support them.  You can use 'em if your strtod() takes 'em.
                 */
-               if (strcasecmp(num, "NaN") == 0)
-                       val = NAN;
+               if (pg_strncasecmp(num, "NaN", 3) == 0)
+               {
+                       val = get_float4_nan();
+                       endptr = num + 3;
+               }
+               else if (pg_strncasecmp(num, "Infinity", 8) == 0)
+               {
+                       val = get_float4_infinity();
+                       endptr = num + 8;
+               }
+               else if (pg_strncasecmp(num, "+Infinity", 9) == 0)
+               {
+                       val = get_float4_infinity();
+                       endptr = num + 9;
+               }
+               else if (pg_strncasecmp(num, "-Infinity", 9) == 0)
+               {
+                       val = -get_float4_infinity();
+                       endptr = num + 9;
+               }
+               else if (pg_strncasecmp(num, "inf", 3) == 0)
+               {
+                       val = get_float4_infinity();
+                       endptr = num + 3;
+               }
+               else if (pg_strncasecmp(num, "+inf", 4) == 0)
+               {
+                       val = get_float4_infinity();
+                       endptr = num + 4;
+               }
+               else if (pg_strncasecmp(num, "-inf", 4) == 0)
+               {
+                       val = -get_float4_infinity();
+                       endptr = num + 4;
+               }
+               else if (save_errno == ERANGE)
+               {
+                       /*
+                        * Some platforms return ERANGE for denormalized numbers (those
+                        * that are not zero, but are too close to zero to have full
+                        * precision).  We'd prefer not to throw error for that, so try to
+                        * detect whether it's a "real" out-of-range condition by checking
+                        * to see if the result is zero or huge.
+                        */
+                       if (val == 0.0 || val >= HUGE_VAL || val <= -HUGE_VAL)
+                               ereport(ERROR,
+                                               (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+                                                errmsg("\"%s\" is out of range for type real",
+                                                               orig_num)));
+               }
                else
                        ereport(ERROR,
                                        (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
-                                        errmsg("invalid input syntax for float4: \"%s\"",
-                                                       num)));
+                                        errmsg("invalid input syntax for type real: \"%s\"",
+                                                       orig_num)));
        }
+#ifdef HAVE_BUGGY_SOLARIS_STRTOD
        else
        {
-               if (errno == ERANGE)
-                       ereport(ERROR,
-                                       (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
-                                        errmsg("\"%s\" is out of range for float4", num)));
+               /*
+                * Many versions of Solaris have a bug wherein strtod sets endptr to
+                * point one byte beyond the end of the string when given "inf" or
+                * "infinity".
+                */
+               if (endptr != num && endptr[-1] == '\0')
+                       endptr--;
        }
+#endif   /* HAVE_BUGGY_SOLARIS_STRTOD */
+
+       /* skip trailing whitespace */
+       while (*endptr != '\0' && isspace((unsigned char) *endptr))
+               endptr++;
+
+       /* if there is any junk left at the end of the string, bail out */
+       if (*endptr != '\0')
+               ereport(ERROR,
+                               (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
+                                errmsg("invalid input syntax for type real: \"%s\"",
+                                               orig_num)));
 
        /*
-        * if we get here, we have a legal double, still need to check to see
-        * if it's a legal float
+        * if we get here, we have a legal double, still need to check to see if
+        * it's a legal float4
         */
-       CheckFloat4Val(val);
+       CHECKFLOATVAL((float4) val, isinf(val), val == 0);
 
        PG_RETURN_FLOAT4((float4) val);
 }
@@ -226,22 +323,28 @@ float4out(PG_FUNCTION_ARGS)
 {
        float4          num = PG_GETARG_FLOAT4(0);
        char       *ascii = (char *) palloc(MAXFLOATWIDTH + 1);
-       int                     infflag;
-       int                     ndig;
 
        if (isnan(num))
                PG_RETURN_CSTRING(strcpy(ascii, "NaN"));
-       infflag = isinf(num);
-       if (infflag > 0)
-               PG_RETURN_CSTRING(strcpy(ascii, "Infinity"));
-       if (infflag < 0)
-               PG_RETURN_CSTRING(strcpy(ascii, "-Infinity"));
 
-       ndig = FLT_DIG + extra_float_digits;
-       if (ndig < 1)
-               ndig = 1;
-
-       sprintf(ascii, "%.*g", ndig, num);
+       switch (is_infinite(num))
+       {
+               case 1:
+                       strcpy(ascii, "Infinity");
+                       break;
+               case -1:
+                       strcpy(ascii, "-Infinity");
+                       break;
+               default:
+                       {
+                               int                     ndig = FLT_DIG + extra_float_digits;
+
+                               if (ndig < 1)
+                                       ndig = 1;
+
+                               snprintf(ascii, MAXFLOATWIDTH + 1, "%.*g", ndig, num);
+                       }
+       }
 
        PG_RETURN_CSTRING(ascii);
 }
@@ -273,43 +376,135 @@ float4send(PG_FUNCTION_ARGS)
 
 /*
  *             float8in                - converts "num" to float8
- *                                               restricted syntax:
- *                                               {<sp>} [+|-] {digit} [.{digit}] [<exp>]
- *                                               where <sp> is a space, digit is 0-9,
- *                                               <exp> is "e" or "E" followed by an integer.
  */
 Datum
 float8in(PG_FUNCTION_ARGS)
 {
        char       *num = PG_GETARG_CSTRING(0);
+       char       *orig_num;
        double          val;
        char       *endptr;
 
+       /*
+        * endptr points to the first character _after_ the sequence we recognized
+        * as a valid floating point number. orig_num points to the original input
+        * string.
+        */
+       orig_num = num;
+
+       /* skip leading whitespace */
+       while (*num != '\0' && isspace((unsigned char) *num))
+               num++;
+
+       /*
+        * Check for an empty-string input to begin with, to avoid the vagaries of
+        * strtod() on different platforms.
+        */
+       if (*num == '\0')
+               ereport(ERROR,
+                               (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
+                        errmsg("invalid input syntax for type double precision: \"%s\"",
+                                       orig_num)));
+
        errno = 0;
        val = strtod(num, &endptr);
-       if (*endptr != '\0')
+
+       /* did we not see anything that looks like a double? */
+       if (endptr == num || errno != 0)
        {
-               if (strcasecmp(num, "NaN") == 0)
-                       val = NAN;
-               else if (strcasecmp(num, "Infinity") == 0)
-                       val = HUGE_VAL;
-               else if (strcasecmp(num, "-Infinity") == 0)
-                       val = -HUGE_VAL;
+               int                     save_errno = errno;
+
+               /*
+                * C99 requires that strtod() accept NaN, [+-]Infinity, and [+-]Inf,
+                * but not all platforms support all of these (and some accept them
+                * but set ERANGE anyway...)  Therefore, we check for these inputs
+                * ourselves if strtod() fails.
+                *
+                * Note: C99 also requires hexadecimal input as well as some extended
+                * forms of NaN, but we consider these forms unportable and don't try
+                * to support them.  You can use 'em if your strtod() takes 'em.
+                */
+               if (pg_strncasecmp(num, "NaN", 3) == 0)
+               {
+                       val = get_float8_nan();
+                       endptr = num + 3;
+               }
+               else if (pg_strncasecmp(num, "Infinity", 8) == 0)
+               {
+                       val = get_float8_infinity();
+                       endptr = num + 8;
+               }
+               else if (pg_strncasecmp(num, "+Infinity", 9) == 0)
+               {
+                       val = get_float8_infinity();
+                       endptr = num + 9;
+               }
+               else if (pg_strncasecmp(num, "-Infinity", 9) == 0)
+               {
+                       val = -get_float8_infinity();
+                       endptr = num + 9;
+               }
+               else if (pg_strncasecmp(num, "inf", 3) == 0)
+               {
+                       val = get_float8_infinity();
+                       endptr = num + 3;
+               }
+               else if (pg_strncasecmp(num, "+inf", 4) == 0)
+               {
+                       val = get_float8_infinity();
+                       endptr = num + 4;
+               }
+               else if (pg_strncasecmp(num, "-inf", 4) == 0)
+               {
+                       val = -get_float8_infinity();
+                       endptr = num + 4;
+               }
+               else if (save_errno == ERANGE)
+               {
+                       /*
+                        * Some platforms return ERANGE for denormalized numbers (those
+                        * that are not zero, but are too close to zero to have full
+                        * precision).  We'd prefer not to throw error for that, so try to
+                        * detect whether it's a "real" out-of-range condition by checking
+                        * to see if the result is zero or huge.
+                        */
+                       if (val == 0.0 || val >= HUGE_VAL || val <= -HUGE_VAL)
+                               ereport(ERROR,
+                                               (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+                                  errmsg("\"%s\" is out of range for type double precision",
+                                                 orig_num)));
+               }
                else
                        ereport(ERROR,
                                        (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
-                                        errmsg("invalid input syntax for float8: \"%s\"",
-                                                       num)));
+                        errmsg("invalid input syntax for type double precision: \"%s\"",
+                                       orig_num)));
        }
+#ifdef HAVE_BUGGY_SOLARIS_STRTOD
        else
        {
-               if (errno == ERANGE)
-                       ereport(ERROR,
-                                       (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
-                                        errmsg("\"%s\" is out of range for float8", num)));
+               /*
+                * Many versions of Solaris have a bug wherein strtod sets endptr to
+                * point one byte beyond the end of the string when given "inf" or
+                * "infinity".
+                */
+               if (endptr != num && endptr[-1] == '\0')
+                       endptr--;
        }
+#endif   /* HAVE_BUGGY_SOLARIS_STRTOD */
+
+       /* skip trailing whitespace */
+       while (*endptr != '\0' && isspace((unsigned char) *endptr))
+               endptr++;
 
-       CheckFloat8Val(val);
+       /* if there is any junk left at the end of the string, bail out */
+       if (*endptr != '\0')
+               ereport(ERROR,
+                               (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
+                        errmsg("invalid input syntax for type double precision: \"%s\"",
+                                       orig_num)));
+
+       CHECKFLOATVAL(val, true, true);
 
        PG_RETURN_FLOAT8(val);
 }
@@ -323,22 +518,28 @@ float8out(PG_FUNCTION_ARGS)
 {
        float8          num = PG_GETARG_FLOAT8(0);
        char       *ascii = (char *) palloc(MAXDOUBLEWIDTH + 1);
-       int                     infflag;
-       int                     ndig;
 
        if (isnan(num))
                PG_RETURN_CSTRING(strcpy(ascii, "NaN"));
-       infflag = isinf(num);
-       if (infflag > 0)
-               PG_RETURN_CSTRING(strcpy(ascii, "Infinity"));
-       if (infflag < 0)
-               PG_RETURN_CSTRING(strcpy(ascii, "-Infinity"));
-
-       ndig = DBL_DIG + extra_float_digits;
-       if (ndig < 1)
-               ndig = 1;
 
-       sprintf(ascii, "%.*g", ndig, num);
+       switch (is_infinite(num))
+       {
+               case 1:
+                       strcpy(ascii, "Infinity");
+                       break;
+               case -1:
+                       strcpy(ascii, "-Infinity");
+                       break;
+               default:
+                       {
+                               int                     ndig = DBL_DIG + extra_float_digits;
+
+                               if (ndig < 1)
+                                       ndig = 1;
+
+                               snprintf(ascii, MAXDOUBLEWIDTH + 1, "%.*g", ndig, num);
+                       }
+       }
 
        PG_RETURN_CSTRING(ascii);
 }
@@ -396,8 +597,10 @@ Datum
 float4um(PG_FUNCTION_ARGS)
 {
        float4          arg1 = PG_GETARG_FLOAT4(0);
+       float4          result;
 
-       PG_RETURN_FLOAT4((float4) -arg1);
+       result = -arg1;
+       PG_RETURN_FLOAT4(result);
 }
 
 Datum
@@ -449,12 +652,8 @@ Datum
 float8abs(PG_FUNCTION_ARGS)
 {
        float8          arg1 = PG_GETARG_FLOAT8(0);
-       float8          result;
-
-       result = fabs(arg1);
 
-       CheckFloat8Val(result);
-       PG_RETURN_FLOAT8(result);
+       PG_RETURN_FLOAT8(fabs(arg1));
 }
 
 
@@ -467,9 +666,7 @@ float8um(PG_FUNCTION_ARGS)
        float8          arg1 = PG_GETARG_FLOAT8(0);
        float8          result;
 
-       result = ((arg1 != 0) ? -(arg1) : arg1);
-
-       CheckFloat8Val(result);
+       result = -arg1;
        PG_RETURN_FLOAT8(result);
 }
 
@@ -527,11 +724,19 @@ float4pl(PG_FUNCTION_ARGS)
 {
        float4          arg1 = PG_GETARG_FLOAT4(0);
        float4          arg2 = PG_GETARG_FLOAT4(1);
-       double          result;
+       float4          result;
 
        result = arg1 + arg2;
-       CheckFloat4Val(result);
-       PG_RETURN_FLOAT4((float4) result);
+
+       /*
+        * 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);
 }
 
 Datum
@@ -539,11 +744,11 @@ float4mi(PG_FUNCTION_ARGS)
 {
        float4          arg1 = PG_GETARG_FLOAT4(0);
        float4          arg2 = PG_GETARG_FLOAT4(1);
-       double          result;
+       float4          result;
 
        result = arg1 - arg2;
-       CheckFloat4Val(result);
-       PG_RETURN_FLOAT4((float4) result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
+       PG_RETURN_FLOAT4(result);
 }
 
 Datum
@@ -551,11 +756,12 @@ float4mul(PG_FUNCTION_ARGS)
 {
        float4          arg1 = PG_GETARG_FLOAT4(0);
        float4          arg2 = PG_GETARG_FLOAT4(1);
-       double          result;
+       float4          result;
 
        result = arg1 * arg2;
-       CheckFloat4Val(result);
-       PG_RETURN_FLOAT4((float4) result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2),
+                                 arg1 == 0 || arg2 == 0);
+       PG_RETURN_FLOAT4(result);
 }
 
 Datum
@@ -563,18 +769,17 @@ float4div(PG_FUNCTION_ARGS)
 {
        float4          arg1 = PG_GETARG_FLOAT4(0);
        float4          arg2 = PG_GETARG_FLOAT4(1);
-       double          result;
+       float4          result;
 
        if (arg2 == 0.0)
                ereport(ERROR,
                                (errcode(ERRCODE_DIVISION_BY_ZERO),
                                 errmsg("division by zero")));
 
-       /* Do division in float8, then check for overflow */
-       result = (float8) arg1 / (float8) arg2;
+       result = arg1 / arg2;
 
-       CheckFloat4Val(result);
-       PG_RETURN_FLOAT4((float4) result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
+       PG_RETURN_FLOAT4(result);
 }
 
 /*
@@ -592,7 +797,7 @@ float8pl(PG_FUNCTION_ARGS)
 
        result = arg1 + arg2;
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -605,7 +810,7 @@ float8mi(PG_FUNCTION_ARGS)
 
        result = arg1 - arg2;
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -618,7 +823,8 @@ float8mul(PG_FUNCTION_ARGS)
 
        result = arg1 * arg2;
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2),
+                                 arg1 == 0 || arg2 == 0);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -636,7 +842,7 @@ float8div(PG_FUNCTION_ARGS)
 
        result = arg1 / arg2;
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -654,9 +860,9 @@ static 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.
+        * 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))
        {
@@ -743,6 +949,24 @@ btfloat4cmp(PG_FUNCTION_ARGS)
        PG_RETURN_INT32(float4_cmp_internal(arg1, arg2));
 }
 
+static int
+btfloat4fastcmp(Datum x, Datum y, SortSupport ssup)
+{
+       float4          arg1 = DatumGetFloat4(x);
+       float4          arg2 = DatumGetFloat4(y);
+
+       return float4_cmp_internal(arg1, arg2);
+}
+
+Datum
+btfloat4sortsupport(PG_FUNCTION_ARGS)
+{
+       SortSupport ssup = (SortSupport) PG_GETARG_POINTER(0);
+
+       ssup->comparator = btfloat4fastcmp;
+       PG_RETURN_VOID();
+}
+
 /*
  *             float8{eq,ne,lt,le,gt,ge}               - float8/float8 comparison operations
  */
@@ -750,9 +974,9 @@ static 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.
+        * 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))
        {
@@ -839,6 +1063,44 @@ btfloat8cmp(PG_FUNCTION_ARGS)
        PG_RETURN_INT32(float8_cmp_internal(arg1, arg2));
 }
 
+static int
+btfloat8fastcmp(Datum x, Datum y, SortSupport ssup)
+{
+       float8          arg1 = DatumGetFloat8(x);
+       float8          arg2 = DatumGetFloat8(y);
+
+       return float8_cmp_internal(arg1, arg2);
+}
+
+Datum
+btfloat8sortsupport(PG_FUNCTION_ARGS)
+{
+       SortSupport ssup = (SortSupport) PG_GETARG_POINTER(0);
+
+       ssup->comparator = btfloat8fastcmp;
+       PG_RETURN_VOID();
+}
+
+Datum
+btfloat48cmp(PG_FUNCTION_ARGS)
+{
+       float4          arg1 = PG_GETARG_FLOAT4(0);
+       float8          arg2 = PG_GETARG_FLOAT8(1);
+
+       /* widen float4 to float8 and then compare */
+       PG_RETURN_INT32(float8_cmp_internal(arg1, arg2));
+}
+
+Datum
+btfloat84cmp(PG_FUNCTION_ARGS)
+{
+       float8          arg1 = PG_GETARG_FLOAT8(0);
+       float4          arg2 = PG_GETARG_FLOAT4(1);
+
+       /* widen float4 to float8 and then compare */
+       PG_RETURN_INT32(float8_cmp_internal(arg1, arg2));
+}
+
 
 /*
  *             ===================
@@ -866,7 +1128,7 @@ dtof(PG_FUNCTION_ARGS)
 {
        float8          num = PG_GETARG_FLOAT8(0);
 
-       CheckFloat4Val(num);
+       CHECKFLOATVAL((float4) num, isinf(num), num == 0);
 
        PG_RETURN_FLOAT4((float4) num);
 }
@@ -881,7 +1143,8 @@ dtoi4(PG_FUNCTION_ARGS)
        float8          num = PG_GETARG_FLOAT8(0);
        int32           result;
 
-       if ((num < INT_MIN) || (num > INT_MAX))
+       /* 'Inf' is handled by INT_MAX */
+       if (num < INT_MIN || num > INT_MAX || isnan(num))
                ereport(ERROR,
                                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                                 errmsg("integer out of range")));
@@ -898,15 +1161,13 @@ Datum
 dtoi2(PG_FUNCTION_ARGS)
 {
        float8          num = PG_GETARG_FLOAT8(0);
-       int16           result;
 
-       if ((num < SHRT_MIN) || (num > SHRT_MAX))
+       if (num < SHRT_MIN || num > SHRT_MAX || isnan(num))
                ereport(ERROR,
                                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
-                                errmsg("integer out of range")));
+                                errmsg("smallint out of range")));
 
-       result = (int16) rint(num);
-       PG_RETURN_INT16(result);
+       PG_RETURN_INT16((int16) rint(num));
 }
 
 
@@ -917,10 +1178,8 @@ Datum
 i4tod(PG_FUNCTION_ARGS)
 {
        int32           num = PG_GETARG_INT32(0);
-       float8          result;
 
-       result = num;
-       PG_RETURN_FLOAT8(result);
+       PG_RETURN_FLOAT8((float8) num);
 }
 
 
@@ -931,10 +1190,8 @@ Datum
 i2tod(PG_FUNCTION_ARGS)
 {
        int16           num = PG_GETARG_INT16(0);
-       float8          result;
 
-       result = num;
-       PG_RETURN_FLOAT8(result);
+       PG_RETURN_FLOAT8((float8) num);
 }
 
 
@@ -945,15 +1202,13 @@ Datum
 ftoi4(PG_FUNCTION_ARGS)
 {
        float4          num = PG_GETARG_FLOAT4(0);
-       int32           result;
 
-       if ((num < INT_MIN) || (num > INT_MAX))
+       if (num < INT_MIN || num > INT_MAX || isnan(num))
                ereport(ERROR,
                                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                                 errmsg("integer out of range")));
 
-       result = (int32) rint(num);
-       PG_RETURN_INT32(result);
+       PG_RETURN_INT32((int32) rint(num));
 }
 
 
@@ -964,29 +1219,25 @@ Datum
 ftoi2(PG_FUNCTION_ARGS)
 {
        float4          num = PG_GETARG_FLOAT4(0);
-       int16           result;
 
-       if ((num < SHRT_MIN) || (num > SHRT_MAX))
+       if (num < SHRT_MIN || num > SHRT_MAX || isnan(num))
                ereport(ERROR,
                                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
-                                errmsg("integer out of range")));
+                                errmsg("smallint out of range")));
 
-       result = (int16) rint(num);
-       PG_RETURN_INT16(result);
+       PG_RETURN_INT16((int16) rint(num));
 }
 
 
 /*
- *             i4tof                   - converts an int4 number to a float8 number
+ *             i4tof                   - converts an int4 number to a float4 number
  */
 Datum
 i4tof(PG_FUNCTION_ARGS)
 {
        int32           num = PG_GETARG_INT32(0);
-       float4          result;
 
-       result = num;
-       PG_RETURN_FLOAT4(result);
+       PG_RETURN_FLOAT4((float4) num);
 }
 
 
@@ -997,112 +1248,8 @@ Datum
 i2tof(PG_FUNCTION_ARGS)
 {
        int16           num = PG_GETARG_INT16(0);
-       float4          result;
-
-       result = num;
-       PG_RETURN_FLOAT4(result);
-}
-
-
-/*
- *             float8_text             - converts a float8 number to a text string
- */
-Datum
-float8_text(PG_FUNCTION_ARGS)
-{
-       float8          num = PG_GETARG_FLOAT8(0);
-       text       *result;
-       int                     len;
-       char       *str;
-
-       str = DatumGetCString(DirectFunctionCall1(float8out,
-                                                                                         Float8GetDatum(num)));
-
-       len = strlen(str) + VARHDRSZ;
-
-       result = (text *) palloc(len);
-
-       VARATT_SIZEP(result) = len;
-       memcpy(VARDATA(result), str, (len - VARHDRSZ));
-
-       pfree(str);
-
-       PG_RETURN_TEXT_P(result);
-}
 
-
-/*
- *             text_float8             - converts a text string to a float8 number
- */
-Datum
-text_float8(PG_FUNCTION_ARGS)
-{
-       text       *string = PG_GETARG_TEXT_P(0);
-       Datum           result;
-       int                     len;
-       char       *str;
-
-       len = (VARSIZE(string) - VARHDRSZ);
-       str = palloc(len + 1);
-       memcpy(str, VARDATA(string), len);
-       *(str + len) = '\0';
-
-       result = DirectFunctionCall1(float8in, CStringGetDatum(str));
-
-       pfree(str);
-
-       PG_RETURN_DATUM(result);
-}
-
-
-/*
- *             float4_text             - converts a float4 number to a text string
- */
-Datum
-float4_text(PG_FUNCTION_ARGS)
-{
-       float4          num = PG_GETARG_FLOAT4(0);
-       text       *result;
-       int                     len;
-       char       *str;
-
-       str = DatumGetCString(DirectFunctionCall1(float4out,
-                                                                                         Float4GetDatum(num)));
-
-       len = strlen(str) + VARHDRSZ;
-
-       result = (text *) palloc(len);
-
-       VARATT_SIZEP(result) = len;
-       memcpy(VARDATA(result), str, (len - VARHDRSZ));
-
-       pfree(str);
-
-       PG_RETURN_TEXT_P(result);
-}
-
-
-/*
- *             text_float4             - converts a text string to a float4 number
- */
-Datum
-text_float4(PG_FUNCTION_ARGS)
-{
-       text       *string = PG_GETARG_TEXT_P(0);
-       Datum           result;
-       int                     len;
-       char       *str;
-
-       len = (VARSIZE(string) - VARHDRSZ);
-       str = palloc(len + 1);
-       memcpy(str, VARDATA(string), len);
-       *(str + len) = '\0';
-
-       result = DirectFunctionCall1(float4in, CStringGetDatum(str));
-
-       pfree(str);
-
-       PG_RETURN_DATUM(result);
+       PG_RETURN_FLOAT4((float4) num);
 }
 
 
@@ -1119,11 +1266,8 @@ Datum
 dround(PG_FUNCTION_ARGS)
 {
        float8          arg1 = PG_GETARG_FLOAT8(0);
-       float8          result;
 
-       result = rint(arg1);
-
-       PG_RETURN_FLOAT8(result);
+       PG_RETURN_FLOAT8(rint(arg1));
 }
 
 /*
@@ -1204,12 +1348,12 @@ dsqrt(PG_FUNCTION_ARGS)
 
        if (arg1 < 0)
                ereport(ERROR,
-                               (errcode(ERRCODE_FLOATING_POINT_EXCEPTION),
+                               (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
                                 errmsg("cannot take square root of a negative number")));
 
        result = sqrt(arg1);
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1224,6 +1368,7 @@ dcbrt(PG_FUNCTION_ARGS)
        float8          result;
 
        result = cbrt(arg1);
+       CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1239,21 +1384,44 @@ dpow(PG_FUNCTION_ARGS)
        float8          result;
 
        /*
-        * We must check both for errno getting set and for a NaN result, in
-        * order to deal with the vagaries of different platforms...
+        * The SQL spec requires that we emit a particular SQLSTATE error code for
+        * certain error conditions.  Specifically, we don't return a
+        * divide-by-zero error code for 0 ^ -1.
+        */
+       if (arg1 == 0 && arg2 < 0)
+               ereport(ERROR,
+                               (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
+                                errmsg("zero raised to a negative power is undefined")));
+       if (arg1 < 0 && floor(arg2) != arg2)
+               ereport(ERROR,
+                               (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
+                                errmsg("a negative number raised to a non-integer power yields a complex result")));
+
+       /*
+        * pow() sets errno only on some platforms, depending on whether it
+        * follows _IEEE_, _POSIX_, _XOPEN_, or _SVID_, so we try to avoid using
+        * errno.  However, some platform/CPU combinations return errno == EDOM
+        * and result == Nan for negative arg1 and very large arg2 (they must be
+        * using something different from our floor() test to decide it's
+        * invalid).  Other platforms (HPPA) return errno == ERANGE and a large
+        * (HUGE_VAL) but finite result to signal overflow.
         */
        errno = 0;
        result = pow(arg1, arg2);
-       if (errno != 0
-#ifdef HAVE_FINITE
-               || !finite(result)
-#endif
-               )
-               ereport(ERROR,
-                               (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
-                                errmsg("result is out of range")));
+       if (errno == EDOM && isnan(result))
+       {
+               if ((fabs(arg1) > 1 && arg2 >= 0) || (fabs(arg1) < 1 && arg2 < 0))
+                       /* 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();
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1267,30 +1435,18 @@ dexp(PG_FUNCTION_ARGS)
        float8          arg1 = PG_GETARG_FLOAT8(0);
        float8          result;
 
-       /*
-        * We must check both for errno getting set and for a NaN result, in
-        * order to deal with the vagaries of different platforms. Also, a
-        * zero result implies unreported underflow.
-        */
        errno = 0;
        result = exp(arg1);
-       if (errno != 0 || result == 0.0
-#ifdef HAVE_FINITE
-               || !finite(result)
-#endif
-               )
-               ereport(ERROR,
-                               (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
-                                errmsg("result is out of range")));
+       if (errno == ERANGE && result != 0 && !isinf(result))
+               result = get_float8_infinity();
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1), false);
        PG_RETURN_FLOAT8(result);
 }
 
 
 /*
  *             dlog1                   - returns the natural logarithm of arg1
- *                                               ("dlog" is already a logging routine...)
  */
 Datum
 dlog1(PG_FUNCTION_ARGS)
@@ -1298,19 +1454,22 @@ dlog1(PG_FUNCTION_ARGS)
        float8          arg1 = PG_GETARG_FLOAT8(0);
        float8          result;
 
+       /*
+        * Emit particular SQLSTATE error codes for ln(). This is required by the
+        * SQL standard.
+        */
        if (arg1 == 0.0)
                ereport(ERROR,
-                               (errcode(ERRCODE_FLOATING_POINT_EXCEPTION),
-                                errmsg("cannot take log of zero")));
-
+                               (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
+                                errmsg("cannot take logarithm of zero")));
        if (arg1 < 0)
                ereport(ERROR,
-                               (errcode(ERRCODE_FLOATING_POINT_EXCEPTION),
-                                errmsg("cannot take log of a negative number")));
+                               (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
+                                errmsg("cannot take logarithm of a negative number")));
 
        result = log(arg1);
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1), arg1 == 1);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1324,19 +1483,23 @@ dlog10(PG_FUNCTION_ARGS)
        float8          arg1 = PG_GETARG_FLOAT8(0);
        float8          result;
 
+       /*
+        * Emit particular SQLSTATE error codes for log(). The SQL spec doesn't
+        * define log(), but it does define ln(), so it makes sense to emit the
+        * same error code for an analogous error condition.
+        */
        if (arg1 == 0.0)
                ereport(ERROR,
-                               (errcode(ERRCODE_FLOATING_POINT_EXCEPTION),
-                                errmsg("cannot take log of zero")));
-
+                               (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
+                                errmsg("cannot take logarithm of zero")));
        if (arg1 < 0)
                ereport(ERROR,
-                               (errcode(ERRCODE_FLOATING_POINT_EXCEPTION),
-                                errmsg("cannot take log of a negative number")));
+                               (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
+                                errmsg("cannot take logarithm of a negative number")));
 
        result = log10(arg1);
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1), arg1 == 1);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1350,18 +1513,18 @@ dacos(PG_FUNCTION_ARGS)
        float8          arg1 = PG_GETARG_FLOAT8(0);
        float8          result;
 
+       /*
+        * We use errno here because the trigonometric functions are cyclic and
+        * hard to check for underflow.
+        */
        errno = 0;
        result = acos(arg1);
-       if (errno != 0
-#ifdef HAVE_FINITE
-               || !finite(result)
-#endif
-               )
+       if (errno != 0)
                ereport(ERROR,
                                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                                 errmsg("input is out of range")));
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1), true);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1377,16 +1540,12 @@ dasin(PG_FUNCTION_ARGS)
 
        errno = 0;
        result = asin(arg1);
-       if (errno != 0
-#ifdef HAVE_FINITE
-               || !finite(result)
-#endif
-               )
+       if (errno != 0)
                ereport(ERROR,
                                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                                 errmsg("input is out of range")));
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1), true);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1402,16 +1561,12 @@ datan(PG_FUNCTION_ARGS)
 
        errno = 0;
        result = atan(arg1);
-       if (errno != 0
-#ifdef HAVE_FINITE
-               || !finite(result)
-#endif
-               )
+       if (errno != 0)
                ereport(ERROR,
                                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                                 errmsg("input is out of range")));
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1), true);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1428,16 +1583,12 @@ datan2(PG_FUNCTION_ARGS)
 
        errno = 0;
        result = atan2(arg1, arg2);
-       if (errno != 0
-#ifdef HAVE_FINITE
-               || !finite(result)
-#endif
-               )
+       if (errno != 0)
                ereport(ERROR,
                                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                                 errmsg("input is out of range")));
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1453,16 +1604,12 @@ dcos(PG_FUNCTION_ARGS)
 
        errno = 0;
        result = cos(arg1);
-       if (errno != 0
-#ifdef HAVE_FINITE
-               || !finite(result)
-#endif
-               )
+       if (errno != 0)
                ereport(ERROR,
                                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                                 errmsg("input is out of range")));
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1), true);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1478,17 +1625,13 @@ dcot(PG_FUNCTION_ARGS)
 
        errno = 0;
        result = tan(arg1);
-       if (errno != 0 || result == 0.0
-#ifdef HAVE_FINITE
-               || !finite(result)
-#endif
-               )
+       if (errno != 0)
                ereport(ERROR,
                                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                                 errmsg("input is out of range")));
 
        result = 1.0 / result;
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, true /* cotan(pi/2) == inf */ , true);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1504,16 +1647,12 @@ dsin(PG_FUNCTION_ARGS)
 
        errno = 0;
        result = sin(arg1);
-       if (errno != 0
-#ifdef HAVE_FINITE
-               || !finite(result)
-#endif
-               )
+       if (errno != 0)
                ereport(ERROR,
                                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                                 errmsg("input is out of range")));
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1), true);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1529,16 +1668,12 @@ dtan(PG_FUNCTION_ARGS)
 
        errno = 0;
        result = tan(arg1);
-       if (errno != 0
-#ifdef HAVE_FINITE
-               || !finite(result)
-#endif
-               )
+       if (errno != 0)
                ereport(ERROR,
                                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                                 errmsg("input is out of range")));
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, true /* tan(pi/2) == Inf */ , true);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1554,7 +1689,7 @@ degrees(PG_FUNCTION_ARGS)
 
        result = arg1 * (180.0 / M_PI);
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1580,7 +1715,7 @@ radians(PG_FUNCTION_ARGS)
 
        result = arg1 * (M_PI / 180.0);
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1593,8 +1728,8 @@ drandom(PG_FUNCTION_ARGS)
 {
        float8          result;
 
-       /* result 0.0-1.0 */
-       result = ((double) random()) / ((double) MAX_RANDOM_VALUE);
+       /* result [0.0 - 1.0) */
+       result = (double) random() / ((double) MAX_RANDOM_VALUE + 1);
 
        PG_RETURN_FLOAT8(result);
 }
@@ -1607,11 +1742,15 @@ Datum
 setseed(PG_FUNCTION_ARGS)
 {
        float8          seed = PG_GETARG_FLOAT8(0);
-       int                     iseed = (int) (seed * MAX_RANDOM_VALUE);
+       int                     iseed;
 
+       if (seed < -1 || seed > 1)
+               elog(ERROR, "setseed parameter %f out of range [-1,1]", seed);
+
+       iseed = (int) (seed * MAX_RANDOM_VALUE);
        srandom((unsigned int) iseed);
 
-       PG_RETURN_INT32(iseed);
+       PG_RETURN_VOID();
 }
 
 
@@ -1621,11 +1760,13 @@ setseed(PG_FUNCTION_ARGS)
  *             FLOAT AGGREGATE OPERATORS
  *             =========================
  *
- *             float8_accum    - accumulate for AVG(), STDDEV(), etc
- *             float4_accum    - same, but input data is float4
- *             float8_avg              - produce final result for float AVG()
- *             float8_variance - produce final result for float VARIANCE()
- *             float8_stddev   - produce final result for float STDDEV()
+ *             float8_accum            - accumulate for AVG(), variance aggregates, etc.
+ *             float4_accum            - same, but input data is float4
+ *             float8_avg                      - produce final result for float AVG()
+ *             float8_var_samp         - produce final result for float VAR_SAMP()
+ *             float8_var_pop          - produce final result for float VAR_POP()
+ *             float8_stddev_samp      - produce final result for float STDDEV_SAMP()
+ *             float8_stddev_pop       - produce final result for float STDDEV_POP()
  *
  * The transition datatype for all these aggregates is a 3-element array
  * of float8, holding the values N, sum(X), sum(X*X) in that order.
@@ -1637,17 +1778,18 @@ setseed(PG_FUNCTION_ARGS)
  */
 
 static float8 *
-check_float8_array(ArrayType *transarray, const char *caller)
+check_float8_array(ArrayType *transarray, const char *caller, int n)
 {
        /*
-        * We expect the input to be a 3-element float array; verify that. We
+        * We expect the input to be an N-element float array; verify that. We
         * don't need to use deconstruct_array() since the array data is just
-        * going to look like a C array of 3 float8 values.
+        * going to look like a C array of N float8 values.
         */
        if (ARR_NDIM(transarray) != 1 ||
-               ARR_DIMS(transarray)[0] != 3 ||
+               ARR_DIMS(transarray)[0] != n ||
+               ARR_HASNULL(transarray) ||
                ARR_ELEMTYPE(transarray) != FLOAT8OID)
-               elog(ERROR, "%s: expected 3-element float8 array", caller);
+               elog(ERROR, "%s: expected %d-element float8 array", caller, n);
        return (float8 *) ARR_DATA_PTR(transarray);
 }
 
@@ -1660,63 +1802,99 @@ float8_accum(PG_FUNCTION_ARGS)
        float8          N,
                                sumX,
                                sumX2;
-       Datum           transdatums[3];
-       ArrayType  *result;
 
-       transvalues = check_float8_array(transarray, "float8_accum");
+       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);
 
-       transdatums[0] = Float8GetDatumFast(N);
-       transdatums[1] = Float8GetDatumFast(sumX);
-       transdatums[2] = Float8GetDatumFast(sumX2);
+       /*
+        * 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;
 
-       result = construct_array(transdatums, 3,
-                                                        FLOAT8OID,
-                                                sizeof(float8), false /* float8 byval */ , 'd');
+               PG_RETURN_ARRAYTYPE_P(transarray);
+       }
+       else
+       {
+               Datum           transdatums[3];
+               ArrayType  *result;
+
+               transdatums[0] = Float8GetDatumFast(N);
+               transdatums[1] = Float8GetDatumFast(sumX);
+               transdatums[2] = Float8GetDatumFast(sumX2);
 
-       PG_RETURN_ARRAYTYPE_P(result);
+               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);
-       float4          newval4 = PG_GETARG_FLOAT4(1);
+
+       /* do computations as float8 */
+       float8          newval = PG_GETARG_FLOAT4(1);
        float8     *transvalues;
        float8          N,
                                sumX,
-                               sumX2,
-                               newval;
-       Datum           transdatums[3];
-       ArrayType  *result;
+                               sumX2;
 
-       transvalues = check_float8_array(transarray, "float4_accum");
+       transvalues = check_float8_array(transarray, "float4_accum", 3);
        N = transvalues[0];
        sumX = transvalues[1];
        sumX2 = transvalues[2];
 
-       /* Do arithmetic in float8 for best accuracy */
-       newval = newval4;
-
        N += 1.0;
        sumX += newval;
+       CHECKFLOATVAL(sumX, isinf(transvalues[1]) || isinf(newval), true);
        sumX2 += newval * newval;
+       CHECKFLOATVAL(sumX2, isinf(transvalues[2]) || isinf(newval), true);
 
-       transdatums[0] = Float8GetDatumFast(N);
-       transdatums[1] = Float8GetDatumFast(sumX);
-       transdatums[2] = Float8GetDatumFast(sumX2);
+       /*
+        * 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;
+
+               PG_RETURN_ARRAYTYPE_P(transarray);
+       }
+       else
+       {
+               Datum           transdatums[3];
+               ArrayType  *result;
+
+               transdatums[0] = Float8GetDatumFast(N);
+               transdatums[1] = Float8GetDatumFast(sumX);
+               transdatums[2] = Float8GetDatumFast(sumX2);
 
-       result = construct_array(transdatums, 3,
-                                                        FLOAT8OID,
-                                                sizeof(float8), false /* float8 byval */ , 'd');
+               result = construct_array(transdatums, 3,
+                                                                FLOAT8OID,
+                                                                sizeof(float8), FLOAT8PASSBYVAL, 'd');
 
-       PG_RETURN_ARRAYTYPE_P(result);
+               PG_RETURN_ARRAYTYPE_P(result);
+       }
 }
 
 Datum
@@ -1727,12 +1905,12 @@ float8_avg(PG_FUNCTION_ARGS)
        float8          N,
                                sumX;
 
-       transvalues = check_float8_array(transarray, "float8_avg");
+       transvalues = check_float8_array(transarray, "float8_avg", 3);
        N = transvalues[0];
        sumX = transvalues[1];
        /* ignore sumX2 */
 
-       /* SQL92 defines AVG of no values to be NULL */
+       /* SQL defines AVG of no values to be NULL */
        if (N == 0.0)
                PG_RETURN_NULL();
 
@@ -1740,7 +1918,36 @@ float8_avg(PG_FUNCTION_ARGS)
 }
 
 Datum
-float8_variance(PG_FUNCTION_ARGS)
+float8_var_pop(PG_FUNCTION_ARGS)
+{
+       ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
+       float8     *transvalues;
+       float8          N,
+                               sumX,
+                               sumX2,
+                               numerator;
+
+       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);
+
+       /* 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)
 {
        ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
        float8     *transvalues;
@@ -1749,7 +1956,7 @@ float8_variance(PG_FUNCTION_ARGS)
                                sumX2,
                                numerator;
 
-       transvalues = check_float8_array(transarray, "float8_variance");
+       transvalues = check_float8_array(transarray, "float8_var_samp", 3);
        N = transvalues[0];
        sumX = transvalues[1];
        sumX2 = transvalues[2];
@@ -1759,6 +1966,7 @@ float8_variance(PG_FUNCTION_ARGS)
                PG_RETURN_NULL();
 
        numerator = N * sumX2 - sumX * sumX;
+       CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true);
 
        /* Watch out for roundoff error producing a negative numerator */
        if (numerator <= 0.0)
@@ -1768,7 +1976,7 @@ float8_variance(PG_FUNCTION_ARGS)
 }
 
 Datum
-float8_stddev(PG_FUNCTION_ARGS)
+float8_stddev_pop(PG_FUNCTION_ARGS)
 {
        ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
        float8     *transvalues;
@@ -1777,7 +1985,36 @@ float8_stddev(PG_FUNCTION_ARGS)
                                sumX2,
                                numerator;
 
-       transvalues = check_float8_array(transarray, "float8_stddev");
+       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);
+
+       /* 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)
+{
+       ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
+       float8     *transvalues;
+       float8          N,
+                               sumX,
+                               sumX2,
+                               numerator;
+
+       transvalues = check_float8_array(transarray, "float8_stddev_samp", 3);
        N = transvalues[0];
        sumX = transvalues[1];
        sumX2 = transvalues[2];
@@ -1787,6 +2024,7 @@ float8_stddev(PG_FUNCTION_ARGS)
                PG_RETURN_NULL();
 
        numerator = N * sumX2 - sumX * sumX;
+       CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true);
 
        /* Watch out for roundoff error producing a negative numerator */
        if (numerator <= 0.0)
@@ -1795,6 +2033,429 @@ float8_stddev(PG_FUNCTION_ARGS)
        PG_RETURN_FLOAT8(sqrt(numerator / (N * (N - 1.0))));
 }
 
+/*
+ *             =========================
+ *             SQL2003 BINARY AGGREGATES
+ *             =========================
+ *
+ * The transition datatype for all these aggregates is a 6-element array of
+ * float8, holding the values N, sum(X), sum(X*X), sum(Y), sum(Y*Y), sum(X*Y)
+ * in that order.  Note that Y is the first argument to the aggregates!
+ *
+ * It might seem attractive to optimize this by having multiple accumulator
+ * functions that only calculate the sums actually needed.     But on most
+ * modern machines, a couple of extra floating-point multiplies will be
+ * insignificant compared to the other per-tuple overhead, so I've chosen
+ * to minimize code space instead.
+ */
+
+Datum
+float8_regr_accum(PG_FUNCTION_ARGS)
+{
+       ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
+       float8          newvalY = PG_GETARG_FLOAT8(1);
+       float8          newvalX = PG_GETARG_FLOAT8(2);
+       float8     *transvalues;
+       float8          N,
+                               sumX,
+                               sumX2,
+                               sumY,
+                               sumY2,
+                               sumXY;
+
+       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);
+       sumX2 += newvalX * newvalX;
+       CHECKFLOATVAL(sumX2, isinf(transvalues[2]) || isinf(newvalX), true);
+       sumY += newvalY;
+       CHECKFLOATVAL(sumY, isinf(transvalues[3]) || isinf(newvalY), true);
+       sumY2 += newvalY * newvalY;
+       CHECKFLOATVAL(sumY2, isinf(transvalues[4]) || isinf(newvalY), true);
+       sumXY += newvalX * newvalY;
+       CHECKFLOATVAL(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;
+               transvalues[2] = sumX2;
+               transvalues[3] = sumY;
+               transvalues[4] = sumY2;
+               transvalues[5] = sumXY;
+
+               PG_RETURN_ARRAYTYPE_P(transarray);
+       }
+       else
+       {
+               Datum           transdatums[6];
+               ArrayType  *result;
+
+               transdatums[0] = Float8GetDatumFast(N);
+               transdatums[1] = Float8GetDatumFast(sumX);
+               transdatums[2] = Float8GetDatumFast(sumX2);
+               transdatums[3] = Float8GetDatumFast(sumY);
+               transdatums[4] = Float8GetDatumFast(sumY2);
+               transdatums[5] = Float8GetDatumFast(sumXY);
+
+               result = construct_array(transdatums, 6,
+                                                                FLOAT8OID,
+                                                                sizeof(float8), FLOAT8PASSBYVAL, 'd');
+
+               PG_RETURN_ARRAYTYPE_P(result);
+       }
+}
+
+Datum
+float8_regr_sxx(PG_FUNCTION_ARGS)
+{
+       ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
+       float8     *transvalues;
+       float8          N,
+                               sumX,
+                               sumX2,
+                               numerator;
+
+       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);
+
+       /* 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)
+{
+       ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
+       float8     *transvalues;
+       float8          N,
+                               sumY,
+                               sumY2,
+                               numerator;
+
+       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);
+
+       /* 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)
+{
+       ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
+       float8     *transvalues;
+       float8          N,
+                               sumX,
+                               sumY,
+                               sumXY,
+                               numerator;
+
+       transvalues = check_float8_array(transarray, "float8_regr_sxy", 6);
+       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);
+
+       /* 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);
+       float8     *transvalues;
+       float8          N,
+                               sumX;
+
+       transvalues = check_float8_array(transarray, "float8_regr_avgx", 6);
+       N = transvalues[0];
+       sumX = transvalues[1];
+
+       /* if N is 0 we should return NULL */
+       if (N < 1.0)
+               PG_RETURN_NULL();
+
+       PG_RETURN_FLOAT8(sumX / N);
+}
+
+Datum
+float8_regr_avgy(PG_FUNCTION_ARGS)
+{
+       ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
+       float8     *transvalues;
+       float8          N,
+                               sumY;
+
+       transvalues = check_float8_array(transarray, "float8_regr_avgy", 6);
+       N = transvalues[0];
+       sumY = transvalues[3];
+
+       /* if N is 0 we should return NULL */
+       if (N < 1.0)
+               PG_RETURN_NULL();
+
+       PG_RETURN_FLOAT8(sumY / N);
+}
+
+Datum
+float8_covar_pop(PG_FUNCTION_ARGS)
+{
+       ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
+       float8     *transvalues;
+       float8          N,
+                               sumX,
+                               sumY,
+                               sumXY,
+                               numerator;
+
+       transvalues = check_float8_array(transarray, "float8_covar_pop", 6);
+       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);
+
+       PG_RETURN_FLOAT8(numerator / (N * N));
+}
+
+Datum
+float8_covar_samp(PG_FUNCTION_ARGS)
+{
+       ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
+       float8     *transvalues;
+       float8          N,
+                               sumX,
+                               sumY,
+                               sumXY,
+                               numerator;
+
+       transvalues = check_float8_array(transarray, "float8_covar_samp", 6);
+       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);
+
+       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,
+                               sumX,
+                               sumX2,
+                               sumY,
+                               sumY2,
+                               sumXY,
+                               numeratorX,
+                               numeratorY,
+                               numeratorXY;
+
+       transvalues = check_float8_array(transarray, "float8_corr", 6);
+       N = transvalues[0];
+       sumX = transvalues[1];
+       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);
+       numeratorY = N * sumY2 - sumY * sumY;
+       CHECKFLOATVAL(numeratorY, isinf(sumY2) || isinf(sumY), true);
+       numeratorXY = N * sumXY - sumX * sumY;
+       CHECKFLOATVAL(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);
+       float8     *transvalues;
+       float8          N,
+                               sumX,
+                               sumX2,
+                               sumY,
+                               sumY2,
+                               sumXY,
+                               numeratorX,
+                               numeratorY,
+                               numeratorXY;
+
+       transvalues = check_float8_array(transarray, "float8_regr_r2", 6);
+       N = transvalues[0];
+       sumX = transvalues[1];
+       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);
+       numeratorY = N * sumY2 - sumY * sumY;
+       CHECKFLOATVAL(numeratorY, isinf(sumY2) || isinf(sumY), true);
+       numeratorXY = N * sumXY - sumX * sumY;
+       CHECKFLOATVAL(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));
+}
+
+Datum
+float8_regr_slope(PG_FUNCTION_ARGS)
+{
+       ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
+       float8     *transvalues;
+       float8          N,
+                               sumX,
+                               sumX2,
+                               sumY,
+                               sumXY,
+                               numeratorX,
+                               numeratorXY;
+
+       transvalues = check_float8_array(transarray, "float8_regr_slope", 6);
+       N = transvalues[0];
+       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);
+       numeratorXY = N * sumXY - sumX * sumY;
+       CHECKFLOATVAL(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);
+       float8     *transvalues;
+       float8          N,
+                               sumX,
+                               sumX2,
+                               sumY,
+                               sumXY,
+                               numeratorX,
+                               numeratorXXY;
+
+       transvalues = check_float8_array(transarray, "float8_regr_intercept", 6);
+       N = transvalues[0];
+       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);
+       numeratorXXY = sumY * sumX2 - sumX * sumXY;
+       CHECKFLOATVAL(numeratorXXY, isinf(sumY) || isinf(sumX2) ||
+                                 isinf(sumX) || isinf(sumXY), true);
+       if (numeratorX <= 0)
+               PG_RETURN_NULL();
+
+       PG_RETURN_FLOAT8(numeratorXXY / numeratorX);
+}
+
 
 /*
  *             ====================================
@@ -1816,7 +2477,7 @@ float48pl(PG_FUNCTION_ARGS)
        float8          result;
 
        result = arg1 + arg2;
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1828,7 +2489,7 @@ float48mi(PG_FUNCTION_ARGS)
        float8          result;
 
        result = arg1 - arg2;
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1840,7 +2501,8 @@ float48mul(PG_FUNCTION_ARGS)
        float8          result;
 
        result = arg1 * arg2;
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2),
+                                 arg1 == 0 || arg2 == 0);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1857,7 +2519,7 @@ float48div(PG_FUNCTION_ARGS)
                                 errmsg("division by zero")));
 
        result = arg1 / arg2;
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1876,7 +2538,7 @@ float84pl(PG_FUNCTION_ARGS)
 
        result = arg1 + arg2;
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1889,7 +2551,7 @@ float84mi(PG_FUNCTION_ARGS)
 
        result = arg1 - arg2;
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1902,7 +2564,8 @@ float84mul(PG_FUNCTION_ARGS)
 
        result = arg1 * arg2;
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2),
+                                 arg1 == 0 || arg2 == 0);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -1920,7 +2583,7 @@ float84div(PG_FUNCTION_ARGS)
 
        result = arg1 / arg2;
 
-       CheckFloat8Val(result);
+       CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
        PG_RETURN_FLOAT8(result);
 }
 
@@ -2044,14 +2707,107 @@ float84ge(PG_FUNCTION_ARGS)
        PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) >= 0);
 }
 
+/*
+ * 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
+ * with the specified characteristics. An operand smaller than the
+ * lower bound is assigned to bucket 0. An operand greater than the
+ * upper bound is assigned to an additional bucket (with number
+ * count+1). We don't allow "NaN" for any of the float8 inputs, and we
+ * don't allow either of the histogram bounds to be +/- infinity.
+ */
+Datum
+width_bucket_float8(PG_FUNCTION_ARGS)
+{
+       float8          operand = PG_GETARG_FLOAT8(0);
+       float8          bound1 = PG_GETARG_FLOAT8(1);
+       float8          bound2 = PG_GETARG_FLOAT8(2);
+       int32           count = PG_GETARG_INT32(3);
+       int32           result;
+
+       if (count <= 0.0)
+               ereport(ERROR,
+                               (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
+                                errmsg("count must be greater than zero")));
+
+       if (isnan(operand) || isnan(bound1) || isnan(bound2))
+               ereport(ERROR,
+                               (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
+                        errmsg("operand, lower bound, and upper bound cannot be NaN")));
+
+       /* Note that we allow "operand" to be infinite */
+       if (isinf(bound1) || isinf(bound2))
+               ereport(ERROR,
+                               (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
+                                errmsg("lower and upper bounds must be finite")));
+
+       if (bound1 < bound2)
+       {
+               if (operand < bound1)
+                       result = 0;
+               else if (operand >= bound2)
+               {
+                       result = count + 1;
+                       /* check for overflow */
+                       if (result < count)
+                               ereport(ERROR,
+                                               (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+                                                errmsg("integer out of range")));
+               }
+               else
+                       result = ((float8) count * (operand - bound1) / (bound2 - bound1)) + 1;
+       }
+       else if (bound1 > bound2)
+       {
+               if (operand > bound1)
+                       result = 0;
+               else if (operand <= bound2)
+               {
+                       result = count + 1;
+                       /* check for overflow */
+                       if (result < count)
+                               ereport(ERROR,
+                                               (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
+                                                errmsg("integer out of range")));
+               }
+               else
+                       result = ((float8) count * (bound1 - operand) / (bound1 - bound2)) + 1;
+       }
+       else
+       {
+               ereport(ERROR,
+                               (errcode(ERRCODE_INVALID_ARGUMENT_FOR_WIDTH_BUCKET_FUNCTION),
+                                errmsg("lower bound cannot equal upper bound")));
+               result = 0;                             /* keep the compiler quiet */
+       }
+
+       PG_RETURN_INT32(result);
+}
+
 /* ========== PRIVATE ROUTINES ========== */
 
 #ifndef HAVE_CBRT
+
 static double
 cbrt(double x)
 {
        int                     isneg = (x < 0.0);
-       double          tmpres = pow(fabs(x), (double) 1.0 / (double) 3.0);
+       double          absx = fabs(x);
+       double          tmpres = pow(absx, (double) 1.0 / (double) 3.0);
+
+       /*
+        * The result is somewhat inaccurate --- not really pow()'s fault, as the
+        * exponent it's handed contains roundoff error.  We can improve the
+        * accuracy by doing one iteration of Newton's formula.  Beware of zero
+        * input however.
+        */
+       if (tmpres > 0.0)
+               tmpres -= (tmpres - absx / (tmpres * tmpres)) / (double) 3.0;
 
        return isneg ? -tmpres : tmpres;
 }