float.c

Go to the documentation of this file.

/*-------------------------------------------------------------------------
 *
 * float.c
 *    Functions 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/backend/utils/adt/float.c
 *
 *-------------------------------------------------------------------------
 */
#include "postgres.h"

#include <ctype.h>
#include <float.h>
#include <math.h>
#include <limits.h>

#include "catalog/pg_type.h"
#include "libpq/pqformat.h"
#include "utils/array.h"
#include "utils/builtins.h"
#include "utils/sortsupport.h"


#ifndef M_PI
/* from my RH5.2 gcc math.h file - thomas 2000-04-03 */
#define M_PI 3.14159265358979323846
#endif

/* Radians per degree, a.k.a. PI / 180 */
#define RADIANS_PER_DEGREE 0.0174532925199432957692

/* Visual C++ etc lacks NAN, and won't accept 0.0/0.0.  NAN definition from
 * http://msdn.microsoft.com/library/default.asp?url=/library/en-us/vclang/html/vclrfNotNumberNANItems.asp
 */
#if defined(WIN32) && !defined(NAN)
static const uint32 nan[2] = {0xffffffff, 0x7fffffff};

#define NAN (*(const double *) nan)
#endif

/* not sure what the following should be, but better to make it over-sufficient */
#define MAXFLOATWIDTH   64
#define MAXDOUBLEWIDTH  128

/*
 * check to see if a float4/8 val has underflowed or overflowed
 */
#define CHECKFLOATVAL(val, inf_is_valid, zero_is_valid)         \
do {                                                            \
    if (isinf(val) && !(inf_is_valid))                          \
        ereport(ERROR,                                          \
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),   \
          errmsg("value out of range: overflow")));             \
                                                                \
    if ((val) == 0.0 && !(zero_is_valid))                       \
        ereport(ERROR,                                          \
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),   \
         errmsg("value out of range: underflow")));             \
} while(0)


/* Configurable GUC parameter */
int         extra_float_digits = 0;     /* Added to DBL_DIG or FLT_DIG */

/* Cached constants for degree-based trig functions */
static bool degree_consts_set = false;
static float8 sin_30 = 0;
static float8 one_minus_cos_60 = 0;
static float8 asin_0_5 = 0;
static float8 acos_0_5 = 0;
static float8 atan_1_0 = 0;
static float8 tan_45 = 0;
static float8 cot_45 = 0;

/*
 * These are intentionally not static; don't "fix" them.  They will never
 * be referenced by other files, much less changed; but we don't want the
 * compiler to know that, else it might try to precompute expressions
 * involving them.  See comments for init_degree_constants().
 */
float8      degree_c_thirty = 30.0;
float8      degree_c_forty_five = 45.0;
float8      degree_c_sixty = 60.0;
float8      degree_c_one_half = 0.5;
float8      degree_c_one = 1.0;

/* Local function prototypes */
static double sind_q1(double x);
static double cosd_q1(double x);
static void init_degree_constants(void);

#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 */


/*
 * Routines to provide reasonably platform-independent handling of
 * infinity and NaN.  We assume that isinf() and isnan() are available
 * and work per spec.  (On some platforms, we have to supply our own;
 * see src/port.)  However, generating an Infinity or NaN in the first
 * place is less well standardized; pre-C99 systems tend not to have C99's
 * INFINITY and NAN macros.  We centralize our workarounds for this here.
 */

double
get_float8_infinity(void)
{
#ifdef INFINITY
    /* C99 standard way */
    return (double) INFINITY;
#else

    /*
     * On some platforms, HUGE_VAL is an infinity, elsewhere it's just the
     * largest normal double.  We assume forcing an overflow will get us a
     * true infinity.
     */
    return (double) (HUGE_VAL * HUGE_VAL);
#endif
}

/*
* The funny placements of the two #pragmas is necessary because of a
* long lived bug in the Microsoft compilers.
* See http://support.microsoft.com/kb/120968/en-us for details
*/
#if (_MSC_VER >= 1800)
#pragma warning(disable:4756)
#endif
float
get_float4_infinity(void)
{
#ifdef INFINITY
    /* C99 standard way */
    return (float) INFINITY;
#else
#if (_MSC_VER >= 1800)
#pragma warning(default:4756)
#endif

    /*
     * On some platforms, HUGE_VAL is an infinity, elsewhere it's just the
     * largest normal double.  We assume forcing an overflow will get us a
     * true infinity.
     */
    return (float) (HUGE_VAL * HUGE_VAL);
#endif
}

double
get_float8_nan(void)
{
    /* (double) NAN doesn't work on some NetBSD/MIPS releases */
#if defined(NAN) && !(defined(__NetBSD__) && defined(__mips__))
    /* C99 standard way */
    return (double) NAN;
#else
    /* Assume we can get a NAN via zero divide */
    return (double) (0.0 / 0.0);
#endif
}

float
get_float4_nan(void)
{
#ifdef NAN
    /* C99 standard way */
    return (float) NAN;
#else
    /* Assume we can get a NAN via zero divide */
    return (float) (0.0 / 0.0);
#endif
}


/*
 * Returns -1 if 'val' represents negative infinity, 1 if 'val'
 * represents (positive) infinity, and 0 otherwise. On some platforms,
 * this is equivalent to the isinf() macro, but not everywhere: C99
 * does not specify that isinf() needs to distinguish between positive
 * and negative infinity.
 */
int
is_infinite(double val)
{
    int         inf = isinf(val);

    if (inf == 0)
        return 0;
    else if (val > 0)
        return 1;
    else
        return -1;
}


/* ========== USER I/O ROUTINES ========== */


/*
 *      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);

    /* did we not see anything that looks like a double? */
    if (endptr == num || errno != 0)
    {
        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_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 type real: \"%s\"",
                            orig_num)));
    }
#ifdef HAVE_BUGGY_SOLARIS_STRTOD
    else
    {
        /*
         * 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 float4
     */
    CHECKFLOATVAL((float4) val, isinf(val), val == 0);

    PG_RETURN_FLOAT4((float4) val);
}

/*
 *      float4out       - converts a float4 number to a string
 *                        using a standard output format
 */
Datum
float4out(PG_FUNCTION_ARGS)
{
    float4      num = PG_GETARG_FLOAT4(0);
    char       *ascii = (char *) palloc(MAXFLOATWIDTH + 1);

    if (isnan(num))
        PG_RETURN_CSTRING(strcpy(ascii, "NaN"));

    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);
}

/*
 *      float4recv          - converts external binary format to float4
 */
Datum
float4recv(PG_FUNCTION_ARGS)
{
    StringInfo  buf = (StringInfo) PG_GETARG_POINTER(0);

    PG_RETURN_FLOAT4(pq_getmsgfloat4(buf));
}

/*
 *      float4send          - converts float4 to binary format
 */
Datum
float4send(PG_FUNCTION_ARGS)
{
    float4      num = PG_GETARG_FLOAT4(0);
    StringInfoData buf;

    pq_begintypsend(&buf);
    pq_sendfloat4(&buf, num);
    PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
}

/*
 *      float8in        - converts "num" to float8
 */
Datum
float8in(PG_FUNCTION_ARGS)
{
    char       *num = PG_GETARG_CSTRING(0);

    PG_RETURN_FLOAT8(float8in_internal(num, NULL, "double precision", num));
}

/*
 * float8in_internal - guts of float8in()
 *
 * This is exposed for use by functions that want a reasonably
 * platform-independent way of inputting doubles.  The behavior is
 * essentially like strtod + ereport on error, but note the following
 * differences:
 * 1. Both leading and trailing whitespace are skipped.
 * 2. If endptr_p is NULL, we throw error if there's trailing junk.
 * Otherwise, it's up to the caller to complain about trailing junk.
 * 3. In event of a syntax error, the report mentions the given type_name
 * and prints orig_string as the input; this is meant to support use of
 * this function with types such as "box" and "point", where what we are
 * parsing here is just a substring of orig_string.
 *
 * "num" could validly be declared "const char *", but that results in an
 * unreasonable amount of extra casting both here and in callers, so we don't.
 */
double
float8in_internal(char *num, char **endptr_p,
                  const char *type_name, const char *orig_string)
{
    double      val;
    char       *endptr;

    /* 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 %s: \"%s\"",
                        type_name, orig_string)));

    errno = 0;
    val = strtod(num, &endptr);

    /* did we not see anything that looks like a double? */
    if (endptr == num || errno != 0)
    {
        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.
             *
             * On error, we intentionally complain about double precision not
             * the given type name, and we print only the part of the string
             * that is the current number.
             */
            if (val == 0.0 || val >= HUGE_VAL || val <= -HUGE_VAL)
            {
                char       *errnumber = pstrdup(num);

                errnumber[endptr - num] = '\0';
                ereport(ERROR,
                        (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                   errmsg("\"%s\" is out of range for type double precision",
                          errnumber)));
            }
        }
        else
            ereport(ERROR,
                    (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
                     errmsg("invalid input syntax for type %s: \"%s\"",
                            type_name, orig_string)));
    }
#ifdef HAVE_BUGGY_SOLARIS_STRTOD
    else
    {
        /*
         * 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++;

    /* report stopping point if wanted, else complain if not end of string */
    if (endptr_p)
        *endptr_p = endptr;
    else if (*endptr != '\0')
        ereport(ERROR,
                (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
                 errmsg("invalid input syntax for type %s: \"%s\"",
                        type_name, orig_string)));

    return val;
}

/*
 *      float8out       - converts float8 number to a string
 *                        using a standard output format
 */
Datum
float8out(PG_FUNCTION_ARGS)
{
    float8      num = PG_GETARG_FLOAT8(0);

    PG_RETURN_CSTRING(float8out_internal(num));
}

/*
 * float8out_internal - guts of float8out()
 *
 * This is exposed for use by functions that want a reasonably
 * platform-independent way of outputting doubles.
 * The result is always palloc'd.
 */
char *
float8out_internal(double num)
{
    char       *ascii = (char *) palloc(MAXDOUBLEWIDTH + 1);

    if (isnan(num))
        return strcpy(ascii, "NaN");

    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);
            }
    }

    return ascii;
}

/*
 *      float8recv          - converts external binary format to float8
 */
Datum
float8recv(PG_FUNCTION_ARGS)
{
    StringInfo  buf = (StringInfo) PG_GETARG_POINTER(0);

    PG_RETURN_FLOAT8(pq_getmsgfloat8(buf));
}

/*
 *      float8send          - converts float8 to binary format
 */
Datum
float8send(PG_FUNCTION_ARGS)
{
    float8      num = PG_GETARG_FLOAT8(0);
    StringInfoData buf;

    pq_begintypsend(&buf);
    pq_sendfloat8(&buf, num);
    PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
}


/* ========== PUBLIC ROUTINES ========== */


/*
 *      ======================
 *      FLOAT4 BASE OPERATIONS
 *      ======================
 */

/*
 *      float4abs       - returns |arg1| (absolute value)
 */
Datum
float4abs(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);

    PG_RETURN_FLOAT4((float4) fabs(arg1));
}

/*
 *      float4um        - returns -arg1 (unary minus)
 */
Datum
float4um(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float4      result;

    result = -arg1;
    PG_RETURN_FLOAT4(result);
}

Datum
float4up(PG_FUNCTION_ARGS)
{
    float4      arg = PG_GETARG_FLOAT4(0);

    PG_RETURN_FLOAT4(arg);
}

Datum
float4larger(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);
    float4      result;

    if (float4_cmp_internal(arg1, arg2) > 0)
        result = arg1;
    else
        result = arg2;
    PG_RETURN_FLOAT4(result);
}

Datum
float4smaller(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);
    float4      result;

    if (float4_cmp_internal(arg1, arg2) < 0)
        result = arg1;
    else
        result = arg2;
    PG_RETURN_FLOAT4(result);
}

/*
 *      ======================
 *      FLOAT8 BASE OPERATIONS
 *      ======================
 */

/*
 *      float8abs       - returns |arg1| (absolute value)
 */
Datum
float8abs(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);

    PG_RETURN_FLOAT8(fabs(arg1));
}


/*
 *      float8um        - returns -arg1 (unary minus)
 */
Datum
float8um(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    result = -arg1;
    PG_RETURN_FLOAT8(result);
}

Datum
float8up(PG_FUNCTION_ARGS)
{
    float8      arg = PG_GETARG_FLOAT8(0);

    PG_RETURN_FLOAT8(arg);
}

Datum
float8larger(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);
    float8      result;

    if (float8_cmp_internal(arg1, arg2) > 0)
        result = arg1;
    else
        result = arg2;
    PG_RETURN_FLOAT8(result);
}

Datum
float8smaller(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);
    float8      result;

    if (float8_cmp_internal(arg1, arg2) < 0)
        result = arg1;
    else
        result = arg2;
    PG_RETURN_FLOAT8(result);
}


/*
 *      ====================
 *      ARITHMETIC OPERATORS
 *      ====================
 */

/*
 *      float4pl        - returns arg1 + arg2
 *      float4mi        - returns arg1 - arg2
 *      float4mul       - returns arg1 * arg2
 *      float4div       - returns arg1 / arg2
 */
Datum
float4pl(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);
    float4      result;

    result = arg1 + arg2;

    /*
     * There isn't any way to check for underflow of addition/subtraction
     * because numbers near the underflow value have already been rounded to
     * the point where we can't detect that the two values were originally
     * different, e.g. on x86, '1e-45'::float4 == '2e-45'::float4 ==
     * 1.4013e-45.
     */
    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
    PG_RETURN_FLOAT4(result);
}

Datum
float4mi(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);
    float4      result;

    result = arg1 - arg2;
    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
    PG_RETURN_FLOAT4(result);
}

Datum
float4mul(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);
    float4      result;

    result = arg1 * arg2;
    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2),
                  arg1 == 0 || arg2 == 0);
    PG_RETURN_FLOAT4(result);
}

Datum
float4div(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);
    float4      result;

    if (arg2 == 0.0)
        ereport(ERROR,
                (errcode(ERRCODE_DIVISION_BY_ZERO),
                 errmsg("division by zero")));

    result = arg1 / arg2;

    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
    PG_RETURN_FLOAT4(result);
}

/*
 *      float8pl        - returns arg1 + arg2
 *      float8mi        - returns arg1 - arg2
 *      float8mul       - returns arg1 * arg2
 *      float8div       - returns arg1 / arg2
 */
Datum
float8pl(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);
    float8      result;

    result = arg1 + arg2;

    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
    PG_RETURN_FLOAT8(result);
}

Datum
float8mi(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);
    float8      result;

    result = arg1 - arg2;

    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
    PG_RETURN_FLOAT8(result);
}

Datum
float8mul(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);
    float8      result;

    result = arg1 * arg2;

    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2),
                  arg1 == 0 || arg2 == 0);
    PG_RETURN_FLOAT8(result);
}

Datum
float8div(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);
    float8      result;

    if (arg2 == 0.0)
        ereport(ERROR,
                (errcode(ERRCODE_DIVISION_BY_ZERO),
                 errmsg("division by zero")));

    result = arg1 / arg2;

    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
    PG_RETURN_FLOAT8(result);
}


/*
 *      ====================
 *      COMPARISON OPERATORS
 *      ====================
 */

/*
 *      float4{eq,ne,lt,le,gt,ge}       - float4/float4 comparison operations
 */
int
float4_cmp_internal(float4 a, float4 b)
{
    /*
     * We consider all NANs to be equal and larger than any non-NAN. This is
     * somewhat arbitrary; the important thing is to have a consistent sort
     * order.
     */
    if (isnan(a))
    {
        if (isnan(b))
            return 0;           /* NAN = NAN */
        else
            return 1;           /* NAN > non-NAN */
    }
    else if (isnan(b))
    {
        return -1;              /* non-NAN < NAN */
    }
    else
    {
        if (a > b)
            return 1;
        else if (a < b)
            return -1;
        else
            return 0;
    }
}

Datum
float4eq(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);

    PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) == 0);
}

Datum
float4ne(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);

    PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) != 0);
}

Datum
float4lt(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);

    PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) < 0);
}

Datum
float4le(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);

    PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) <= 0);
}

Datum
float4gt(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);

    PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) > 0);
}

Datum
float4ge(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);

    PG_RETURN_BOOL(float4_cmp_internal(arg1, arg2) >= 0);
}

Datum
btfloat4cmp(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);

    PG_RETURN_INT32(float4_cmp_internal(arg1, arg2));
}

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
 */
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;
    }
}

Datum
float8eq(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) == 0);
}

Datum
float8ne(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) != 0);
}

Datum
float8lt(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) < 0);
}

Datum
float8le(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) <= 0);
}

Datum
float8gt(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) > 0);
}

Datum
float8ge(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) >= 0);
}

Datum
btfloat8cmp(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);

    PG_RETURN_INT32(float8_cmp_internal(arg1, arg2));
}

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));
}


/*
 *      ===================
 *      CONVERSION ROUTINES
 *      ===================
 */

/*
 *      ftod            - converts a float4 number to a float8 number
 */
Datum
ftod(PG_FUNCTION_ARGS)
{
    float4      num = PG_GETARG_FLOAT4(0);

    PG_RETURN_FLOAT8((float8) num);
}


/*
 *      dtof            - converts a float8 number to a float4 number
 */
Datum
dtof(PG_FUNCTION_ARGS)
{
    float8      num = PG_GETARG_FLOAT8(0);

    CHECKFLOATVAL((float4) num, isinf(num), num == 0);

    PG_RETURN_FLOAT4((float4) num);
}


/*
 *      dtoi4           - converts a float8 number to an int4 number
 */
Datum
dtoi4(PG_FUNCTION_ARGS)
{
    float8      num = PG_GETARG_FLOAT8(0);
    int32       result;

    /* '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")));

    result = (int32) rint(num);
    PG_RETURN_INT32(result);
}


/*
 *      dtoi2           - converts a float8 number to an int2 number
 */
Datum
dtoi2(PG_FUNCTION_ARGS)
{
    float8      num = PG_GETARG_FLOAT8(0);

    if (num < SHRT_MIN || num > SHRT_MAX || isnan(num))
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("smallint out of range")));

    PG_RETURN_INT16((int16) rint(num));
}


/*
 *      i4tod           - converts an int4 number to a float8 number
 */
Datum
i4tod(PG_FUNCTION_ARGS)
{
    int32       num = PG_GETARG_INT32(0);

    PG_RETURN_FLOAT8((float8) num);
}


/*
 *      i2tod           - converts an int2 number to a float8 number
 */
Datum
i2tod(PG_FUNCTION_ARGS)
{
    int16       num = PG_GETARG_INT16(0);

    PG_RETURN_FLOAT8((float8) num);
}


/*
 *      ftoi4           - converts a float4 number to an int4 number
 */
Datum
ftoi4(PG_FUNCTION_ARGS)
{
    float4      num = PG_GETARG_FLOAT4(0);

    if (num < INT_MIN || num > INT_MAX || isnan(num))
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("integer out of range")));

    PG_RETURN_INT32((int32) rint(num));
}


/*
 *      ftoi2           - converts a float4 number to an int2 number
 */
Datum
ftoi2(PG_FUNCTION_ARGS)
{
    float4      num = PG_GETARG_FLOAT4(0);

    if (num < SHRT_MIN || num > SHRT_MAX || isnan(num))
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("smallint out of range")));

    PG_RETURN_INT16((int16) rint(num));
}


/*
 *      i4tof           - converts an int4 number to a float4 number
 */
Datum
i4tof(PG_FUNCTION_ARGS)
{
    int32       num = PG_GETARG_INT32(0);

    PG_RETURN_FLOAT4((float4) num);
}


/*
 *      i2tof           - converts an int2 number to a float4 number
 */
Datum
i2tof(PG_FUNCTION_ARGS)
{
    int16       num = PG_GETARG_INT16(0);

    PG_RETURN_FLOAT4((float4) num);
}


/*
 *      =======================
 *      RANDOM FLOAT8 OPERATORS
 *      =======================
 */

/*
 *      dround          - returns   ROUND(arg1)
 */
Datum
dround(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);

    PG_RETURN_FLOAT8(rint(arg1));
}

/*
 *      dceil           - returns the smallest integer greater than or
 *                        equal to the specified float
 */
Datum
dceil(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);

    PG_RETURN_FLOAT8(ceil(arg1));
}

/*
 *      dfloor          - returns the largest integer lesser than or
 *                        equal to the specified float
 */
Datum
dfloor(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);

    PG_RETURN_FLOAT8(floor(arg1));
}

/*
 *      dsign           - returns -1 if the argument is less than 0, 0
 *                        if the argument is equal to 0, and 1 if the
 *                        argument is greater than zero.
 */
Datum
dsign(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    if (arg1 > 0)
        result = 1.0;
    else if (arg1 < 0)
        result = -1.0;
    else
        result = 0.0;

    PG_RETURN_FLOAT8(result);
}

/*
 *      dtrunc          - returns truncation-towards-zero of arg1,
 *                        arg1 >= 0 ... the greatest integer less
 *                                      than or equal to arg1
 *                        arg1 < 0  ... the least integer greater
 *                                      than or equal to arg1
 */
Datum
dtrunc(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    if (arg1 >= 0)
        result = floor(arg1);
    else
        result = -floor(-arg1);

    PG_RETURN_FLOAT8(result);
}


/*
 *      dsqrt           - returns square root of arg1
 */
Datum
dsqrt(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    if (arg1 < 0)
        ereport(ERROR,
                (errcode(ERRCODE_INVALID_ARGUMENT_FOR_POWER_FUNCTION),
                 errmsg("cannot take square root of a negative number")));

    result = sqrt(arg1);

    CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
    PG_RETURN_FLOAT8(result);
}


/*
 *      dcbrt           - returns cube root of arg1
 */
Datum
dcbrt(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    result = cbrt(arg1);
    CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
    PG_RETURN_FLOAT8(result);
}


/*
 *      dpow            - returns pow(arg1,arg2)
 */
Datum
dpow(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);
    float8      result;

    /*
     * 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 == 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();

    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
    PG_RETURN_FLOAT8(result);
}


/*
 *      dexp            - returns the exponential function of arg1
 */
Datum
dexp(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    errno = 0;
    result = exp(arg1);
    if (errno == ERANGE && result != 0 && !isinf(result))
        result = get_float8_infinity();

    CHECKFLOATVAL(result, isinf(arg1), false);
    PG_RETURN_FLOAT8(result);
}


/*
 *      dlog1           - returns the natural logarithm of arg1
 */
Datum
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_INVALID_ARGUMENT_FOR_LOG),
                 errmsg("cannot take logarithm of zero")));
    if (arg1 < 0)
        ereport(ERROR,
                (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
                 errmsg("cannot take logarithm of a negative number")));

    result = log(arg1);

    CHECKFLOATVAL(result, isinf(arg1), arg1 == 1);
    PG_RETURN_FLOAT8(result);
}


/*
 *      dlog10          - returns the base 10 logarithm of arg1
 */
Datum
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_INVALID_ARGUMENT_FOR_LOG),
                 errmsg("cannot take logarithm of zero")));
    if (arg1 < 0)
        ereport(ERROR,
                (errcode(ERRCODE_INVALID_ARGUMENT_FOR_LOG),
                 errmsg("cannot take logarithm of a negative number")));

    result = log10(arg1);

    CHECKFLOATVAL(result, isinf(arg1), arg1 == 1);
    PG_RETURN_FLOAT8(result);
}


/*
 *      dacos           - returns the arccos of arg1 (radians)
 */
Datum
dacos(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    /* Per the POSIX spec, return NaN if the input is NaN */
    if (isnan(arg1))
        PG_RETURN_FLOAT8(get_float8_nan());

    /*
     * The principal branch of the inverse cosine function maps values in the
     * range [-1, 1] to values in the range [0, Pi], so we should reject any
     * inputs outside that range and the result will always be finite.
     */
    if (arg1 < -1.0 || arg1 > 1.0)
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("input is out of range")));

    result = acos(arg1);

    CHECKFLOATVAL(result, false, true);
    PG_RETURN_FLOAT8(result);
}


/*
 *      dasin           - returns the arcsin of arg1 (radians)
 */
Datum
dasin(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    /* Per the POSIX spec, return NaN if the input is NaN */
    if (isnan(arg1))
        PG_RETURN_FLOAT8(get_float8_nan());

    /*
     * The principal branch of the inverse sine function maps values in the
     * range [-1, 1] to values in the range [-Pi/2, Pi/2], so we should reject
     * any inputs outside that range and the result will always be finite.
     */
    if (arg1 < -1.0 || arg1 > 1.0)
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("input is out of range")));

    result = asin(arg1);

    CHECKFLOATVAL(result, false, true);
    PG_RETURN_FLOAT8(result);
}


/*
 *      datan           - returns the arctan of arg1 (radians)
 */
Datum
datan(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    /* Per the POSIX spec, return NaN if the input is NaN */
    if (isnan(arg1))
        PG_RETURN_FLOAT8(get_float8_nan());

    /*
     * The principal branch of the inverse tangent function maps all inputs to
     * values in the range [-Pi/2, Pi/2], so the result should always be
     * finite, even if the input is infinite.
     */
    result = atan(arg1);

    CHECKFLOATVAL(result, false, true);
    PG_RETURN_FLOAT8(result);
}


/*
 *      atan2           - returns the arctan of arg1/arg2 (radians)
 */
Datum
datan2(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);
    float8      result;

    /* Per the POSIX spec, return NaN if either input is NaN */
    if (isnan(arg1) || isnan(arg2))
        PG_RETURN_FLOAT8(get_float8_nan());

    /*
     * atan2 maps all inputs to values in the range [-Pi, Pi], so the result
     * should always be finite, even if the inputs are infinite.
     */
    result = atan2(arg1, arg2);

    CHECKFLOATVAL(result, false, true);
    PG_RETURN_FLOAT8(result);
}


/*
 *      dcos            - returns the cosine of arg1 (radians)
 */
Datum
dcos(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    /* Per the POSIX spec, return NaN if the input is NaN */
    if (isnan(arg1))
        PG_RETURN_FLOAT8(get_float8_nan());

    /*
     * cos() is periodic and so theoretically can work for all finite inputs,
     * but some implementations may choose to throw error if the input is so
     * large that there are no significant digits in the result.  So we should
     * check for errors.  POSIX allows an error to be reported either via
     * errno or via fetestexcept(), but currently we only support checking
     * errno.  (fetestexcept() is rumored to report underflow unreasonably
     * early on some platforms, so it's not clear that believing it would be a
     * net improvement anyway.)
     *
     * For infinite inputs, POSIX specifies that the trigonometric functions
     * should return a domain error; but we won't notice that unless the
     * platform reports via errno, so also explicitly test for infinite
     * inputs.
     */
    errno = 0;
    result = cos(arg1);
    if (errno != 0 || isinf(arg1))
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("input is out of range")));

    CHECKFLOATVAL(result, false, true);
    PG_RETURN_FLOAT8(result);
}


/*
 *      dcot            - returns the cotangent of arg1 (radians)
 */
Datum
dcot(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    /* Per the POSIX spec, return NaN if the input is NaN */
    if (isnan(arg1))
        PG_RETURN_FLOAT8(get_float8_nan());

    /* Be sure to throw an error if the input is infinite --- see dcos() */
    errno = 0;
    result = tan(arg1);
    if (errno != 0 || isinf(arg1))
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("input is out of range")));

    result = 1.0 / result;
    CHECKFLOATVAL(result, true /* cot(0) == Inf */ , true);
    PG_RETURN_FLOAT8(result);
}


/*
 *      dsin            - returns the sine of arg1 (radians)
 */
Datum
dsin(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    /* Per the POSIX spec, return NaN if the input is NaN */
    if (isnan(arg1))
        PG_RETURN_FLOAT8(get_float8_nan());

    /* Be sure to throw an error if the input is infinite --- see dcos() */
    errno = 0;
    result = sin(arg1);
    if (errno != 0 || isinf(arg1))
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("input is out of range")));

    CHECKFLOATVAL(result, false, true);
    PG_RETURN_FLOAT8(result);
}


/*
 *      dtan            - returns the tangent of arg1 (radians)
 */
Datum
dtan(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    /* Per the POSIX spec, return NaN if the input is NaN */
    if (isnan(arg1))
        PG_RETURN_FLOAT8(get_float8_nan());

    /* Be sure to throw an error if the input is infinite --- see dcos() */
    errno = 0;
    result = tan(arg1);
    if (errno != 0 || isinf(arg1))
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("input is out of range")));

    CHECKFLOATVAL(result, true /* tan(pi/2) == Inf */ , true);
    PG_RETURN_FLOAT8(result);
}


/* ========== DEGREE-BASED TRIGONOMETRIC FUNCTIONS ========== */


/*
 * Initialize the cached constants declared at the head of this file
 * (sin_30 etc).  The fact that we need those at all, let alone need this
 * Rube-Goldberg-worthy method of initializing them, is because there are
 * compilers out there that will precompute expressions such as sin(constant)
 * using a sin() function different from what will be used at runtime.  If we
 * want exact results, we must ensure that none of the scaling constants used
 * in the degree-based trig functions are computed that way.  To do so, we
 * compute them from the variables degree_c_thirty etc, which are also really
 * constants, but the compiler cannot assume that.
 *
 * Other hazards we are trying to forestall with this kluge include the
 * possibility that compilers will rearrange the expressions, or compute
 * some intermediate results in registers wider than a standard double.
 *
 * In the places where we use these constants, the typical pattern is like
 *      volatile float8 sin_x = sin(x * RADIANS_PER_DEGREE);
 *      return (sin_x / sin_30);
 * where we hope to get a value of exactly 1.0 from the division when x = 30.
 * The volatile temporary variable is needed on machines with wide float
 * registers, to ensure that the result of sin(x) is rounded to double width
 * the same as the value of sin_30 has been.  Experimentation with gcc shows
 * that marking the temp variable volatile is necessary to make the store and
 * reload actually happen; hopefully the same trick works for other compilers.
 * (gcc's documentation suggests using the -ffloat-store compiler switch to
 * ensure this, but that is compiler-specific and it also pessimizes code in
 * many places where we don't care about this.)
 */
static void
init_degree_constants(void)
{
    sin_30 = sin(degree_c_thirty * RADIANS_PER_DEGREE);
    one_minus_cos_60 = 1.0 - cos(degree_c_sixty * RADIANS_PER_DEGREE);
    asin_0_5 = asin(degree_c_one_half);
    acos_0_5 = acos(degree_c_one_half);
    atan_1_0 = atan(degree_c_one);
    tan_45 = sind_q1(degree_c_forty_five) / cosd_q1(degree_c_forty_five);
    cot_45 = cosd_q1(degree_c_forty_five) / sind_q1(degree_c_forty_five);
    degree_consts_set = true;
}

#define INIT_DEGREE_CONSTANTS() \
do { \
    if (!degree_consts_set) \
        init_degree_constants(); \
} while(0)


/*
 *      asind_q1        - returns the inverse sine of x in degrees, for x in
 *                        the range [0, 1].  The result is an angle in the
 *                        first quadrant --- [0, 90] degrees.
 *
 *                        For the 3 special case inputs (0, 0.5 and 1), this
 *                        function will return exact values (0, 30 and 90
 *                        degrees respectively).
 */
static double
asind_q1(double x)
{
    /*
     * Stitch together inverse sine and cosine functions for the ranges [0,
     * 0.5] and (0.5, 1].  Each expression below is guaranteed to return
     * exactly 30 for x=0.5, so the result is a continuous monotonic function
     * over the full range.
     */
    if (x <= 0.5)
    {
        volatile float8 asin_x = asin(x);

        return (asin_x / asin_0_5) * 30.0;
    }
    else
    {
        volatile float8 acos_x = acos(x);

        return 90.0 - (acos_x / acos_0_5) * 60.0;
    }
}


/*
 *      acosd_q1        - returns the inverse cosine of x in degrees, for x in
 *                        the range [0, 1].  The result is an angle in the
 *                        first quadrant --- [0, 90] degrees.
 *
 *                        For the 3 special case inputs (0, 0.5 and 1), this
 *                        function will return exact values (0, 60 and 90
 *                        degrees respectively).
 */
static double
acosd_q1(double x)
{
    /*
     * Stitch together inverse sine and cosine functions for the ranges [0,
     * 0.5] and (0.5, 1].  Each expression below is guaranteed to return
     * exactly 60 for x=0.5, so the result is a continuous monotonic function
     * over the full range.
     */
    if (x <= 0.5)
    {
        volatile float8 asin_x = asin(x);

        return 90.0 - (asin_x / asin_0_5) * 30.0;
    }
    else
    {
        volatile float8 acos_x = acos(x);

        return (acos_x / acos_0_5) * 60.0;
    }
}


/*
 *      dacosd          - returns the arccos of arg1 (degrees)
 */
Datum
dacosd(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    /* Per the POSIX spec, return NaN if the input is NaN */
    if (isnan(arg1))
        PG_RETURN_FLOAT8(get_float8_nan());

    INIT_DEGREE_CONSTANTS();

    /*
     * The principal branch of the inverse cosine function maps values in the
     * range [-1, 1] to values in the range [0, 180], so we should reject any
     * inputs outside that range and the result will always be finite.
     */
    if (arg1 < -1.0 || arg1 > 1.0)
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("input is out of range")));

    if (arg1 >= 0.0)
        result = acosd_q1(arg1);
    else
        result = 90.0 + asind_q1(-arg1);

    CHECKFLOATVAL(result, false, true);
    PG_RETURN_FLOAT8(result);
}


/*
 *      dasind          - returns the arcsin of arg1 (degrees)
 */
Datum
dasind(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    /* Per the POSIX spec, return NaN if the input is NaN */
    if (isnan(arg1))
        PG_RETURN_FLOAT8(get_float8_nan());

    INIT_DEGREE_CONSTANTS();

    /*
     * The principal branch of the inverse sine function maps values in the
     * range [-1, 1] to values in the range [-90, 90], so we should reject any
     * inputs outside that range and the result will always be finite.
     */
    if (arg1 < -1.0 || arg1 > 1.0)
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("input is out of range")));

    if (arg1 >= 0.0)
        result = asind_q1(arg1);
    else
        result = -asind_q1(-arg1);

    CHECKFLOATVAL(result, false, true);
    PG_RETURN_FLOAT8(result);
}


/*
 *      datand          - returns the arctan of arg1 (degrees)
 */
Datum
datand(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;
    volatile float8 atan_arg1;

    /* Per the POSIX spec, return NaN if the input is NaN */
    if (isnan(arg1))
        PG_RETURN_FLOAT8(get_float8_nan());

    INIT_DEGREE_CONSTANTS();

    /*
     * The principal branch of the inverse tangent function maps all inputs to
     * values in the range [-90, 90], so the result should always be finite,
     * even if the input is infinite.  Additionally, we take care to ensure
     * than when arg1 is 1, the result is exactly 45.
     */
    atan_arg1 = atan(arg1);
    result = (atan_arg1 / atan_1_0) * 45.0;

    CHECKFLOATVAL(result, false, true);
    PG_RETURN_FLOAT8(result);
}


/*
 *      atan2d          - returns the arctan of arg1/arg2 (degrees)
 */
Datum
datan2d(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);
    float8      result;
    volatile float8 atan2_arg1_arg2;

    /* Per the POSIX spec, return NaN if either input is NaN */
    if (isnan(arg1) || isnan(arg2))
        PG_RETURN_FLOAT8(get_float8_nan());

    INIT_DEGREE_CONSTANTS();

    /*
     * atan2d maps all inputs to values in the range [-180, 180], so the
     * result should always be finite, even if the inputs are infinite.
     *
     * Note: this coding assumes that atan(1.0) is a suitable scaling constant
     * to get an exact result from atan2().  This might well fail on us at
     * some point, requiring us to decide exactly what inputs we think we're
     * going to guarantee an exact result for.
     */
    atan2_arg1_arg2 = atan2(arg1, arg2);
    result = (atan2_arg1_arg2 / atan_1_0) * 45.0;

    CHECKFLOATVAL(result, false, true);
    PG_RETURN_FLOAT8(result);
}


/*
 *      sind_0_to_30    - returns the sine of an angle that lies between 0 and
 *                        30 degrees.  This will return exactly 0 when x is 0,
 *                        and exactly 0.5 when x is 30 degrees.
 */
static double
sind_0_to_30(double x)
{
    volatile float8 sin_x = sin(x * RADIANS_PER_DEGREE);

    return (sin_x / sin_30) / 2.0;
}


/*
 *      cosd_0_to_60    - returns the cosine of an angle that lies between 0
 *                        and 60 degrees.  This will return exactly 1 when x
 *                        is 0, and exactly 0.5 when x is 60 degrees.
 */
static double
cosd_0_to_60(double x)
{
    volatile float8 one_minus_cos_x = 1.0 - cos(x * RADIANS_PER_DEGREE);

    return 1.0 - (one_minus_cos_x / one_minus_cos_60) / 2.0;
}


/*
 *      sind_q1         - returns the sine of an angle in the first quadrant
 *                        (0 to 90 degrees).
 */
static double
sind_q1(double x)
{
    /*
     * Stitch together the sine and cosine functions for the ranges [0, 30]
     * and (30, 90].  These guarantee to return exact answers at their
     * endpoints, so the overall result is a continuous monotonic function
     * that gives exact results when x = 0, 30 and 90 degrees.
     */
    if (x <= 30.0)
        return sind_0_to_30(x);
    else
        return cosd_0_to_60(90.0 - x);
}


/*
 *      cosd_q1         - returns the cosine of an angle in the first quadrant
 *                        (0 to 90 degrees).
 */
static double
cosd_q1(double x)
{
    /*
     * Stitch together the sine and cosine functions for the ranges [0, 60]
     * and (60, 90].  These guarantee to return exact answers at their
     * endpoints, so the overall result is a continuous monotonic function
     * that gives exact results when x = 0, 60 and 90 degrees.
     */
    if (x <= 60.0)
        return cosd_0_to_60(x);
    else
        return sind_0_to_30(90.0 - x);
}


/*
 *      dcosd           - returns the cosine of arg1 (degrees)
 */
Datum
dcosd(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;
    int         sign = 1;

    /*
     * Per the POSIX spec, return NaN if the input is NaN and throw an error
     * if the input is infinite.
     */
    if (isnan(arg1))
        PG_RETURN_FLOAT8(get_float8_nan());

    if (isinf(arg1))
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("input is out of range")));

    INIT_DEGREE_CONSTANTS();

    /* Reduce the range of the input to [0,90] degrees */
    arg1 = fmod(arg1, 360.0);

    if (arg1 < 0.0)
    {
        /* cosd(-x) = cosd(x) */
        arg1 = -arg1;
    }

    if (arg1 > 180.0)
    {
        /* cosd(360-x) = cosd(x) */
        arg1 = 360.0 - arg1;
    }

    if (arg1 > 90.0)
    {
        /* cosd(180-x) = -cosd(x) */
        arg1 = 180.0 - arg1;
        sign = -sign;
    }

    result = sign * cosd_q1(arg1);

    CHECKFLOATVAL(result, false, true);
    PG_RETURN_FLOAT8(result);
}


/*
 *      dcotd           - returns the cotangent of arg1 (degrees)
 */
Datum
dcotd(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;
    volatile float8 cot_arg1;
    int         sign = 1;

    /*
     * Per the POSIX spec, return NaN if the input is NaN and throw an error
     * if the input is infinite.
     */
    if (isnan(arg1))
        PG_RETURN_FLOAT8(get_float8_nan());

    if (isinf(arg1))
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("input is out of range")));

    INIT_DEGREE_CONSTANTS();

    /* Reduce the range of the input to [0,90] degrees */
    arg1 = fmod(arg1, 360.0);

    if (arg1 < 0.0)
    {
        /* cotd(-x) = -cotd(x) */
        arg1 = -arg1;
        sign = -sign;
    }

    if (arg1 > 180.0)
    {
        /* cotd(360-x) = -cotd(x) */
        arg1 = 360.0 - arg1;
        sign = -sign;
    }

    if (arg1 > 90.0)
    {
        /* cotd(180-x) = -cotd(x) */
        arg1 = 180.0 - arg1;
        sign = -sign;
    }

    cot_arg1 = cosd_q1(arg1) / sind_q1(arg1);
    result = sign * (cot_arg1 / cot_45);

    /*
     * On some machines we get cotd(270) = minus zero, but this isn't always
     * true.  For portability, and because the user constituency for this
     * function probably doesn't want minus zero, force it to plain zero.
     */
    if (result == 0.0)
        result = 0.0;

    CHECKFLOATVAL(result, true /* cotd(0) == Inf */ , true);
    PG_RETURN_FLOAT8(result);
}


/*
 *      dsind           - returns the sine of arg1 (degrees)
 */
Datum
dsind(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;
    int         sign = 1;

    /*
     * Per the POSIX spec, return NaN if the input is NaN and throw an error
     * if the input is infinite.
     */
    if (isnan(arg1))
        PG_RETURN_FLOAT8(get_float8_nan());

    if (isinf(arg1))
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("input is out of range")));

    INIT_DEGREE_CONSTANTS();

    /* Reduce the range of the input to [0,90] degrees */
    arg1 = fmod(arg1, 360.0);

    if (arg1 < 0.0)
    {
        /* sind(-x) = -sind(x) */
        arg1 = -arg1;
        sign = -sign;
    }

    if (arg1 > 180.0)
    {
        /* sind(360-x) = -sind(x) */
        arg1 = 360.0 - arg1;
        sign = -sign;
    }

    if (arg1 > 90.0)
    {
        /* sind(180-x) = sind(x) */
        arg1 = 180.0 - arg1;
    }

    result = sign * sind_q1(arg1);

    CHECKFLOATVAL(result, false, true);
    PG_RETURN_FLOAT8(result);
}


/*
 *      dtand           - returns the tangent of arg1 (degrees)
 */
Datum
dtand(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;
    volatile float8 tan_arg1;
    int         sign = 1;

    /*
     * Per the POSIX spec, return NaN if the input is NaN and throw an error
     * if the input is infinite.
     */
    if (isnan(arg1))
        PG_RETURN_FLOAT8(get_float8_nan());

    if (isinf(arg1))
        ereport(ERROR,
                (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
                 errmsg("input is out of range")));

    INIT_DEGREE_CONSTANTS();

    /* Reduce the range of the input to [0,90] degrees */
    arg1 = fmod(arg1, 360.0);

    if (arg1 < 0.0)
    {
        /* tand(-x) = -tand(x) */
        arg1 = -arg1;
        sign = -sign;
    }

    if (arg1 > 180.0)
    {
        /* tand(360-x) = -tand(x) */
        arg1 = 360.0 - arg1;
        sign = -sign;
    }

    if (arg1 > 90.0)
    {
        /* tand(180-x) = -tand(x) */
        arg1 = 180.0 - arg1;
        sign = -sign;
    }

    tan_arg1 = sind_q1(arg1) / cosd_q1(arg1);
    result = sign * (tan_arg1 / tan_45);

    /*
     * On some machines we get tand(180) = minus zero, but this isn't always
     * true.  For portability, and because the user constituency for this
     * function probably doesn't want minus zero, force it to plain zero.
     */
    if (result == 0.0)
        result = 0.0;

    CHECKFLOATVAL(result, true /* tand(90) == Inf */ , true);
    PG_RETURN_FLOAT8(result);
}


/*
 *      degrees     - returns degrees converted from radians
 */
Datum
degrees(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    result = arg1 / RADIANS_PER_DEGREE;

    CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
    PG_RETURN_FLOAT8(result);
}


/*
 *      dpi             - returns the constant PI
 */
Datum
dpi(PG_FUNCTION_ARGS)
{
    PG_RETURN_FLOAT8(M_PI);
}


/*
 *      radians     - returns radians converted from degrees
 */
Datum
radians(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float8      result;

    result = arg1 * RADIANS_PER_DEGREE;

    CHECKFLOATVAL(result, isinf(arg1), arg1 == 0);
    PG_RETURN_FLOAT8(result);
}


/*
 *      drandom     - returns a random number
 */
Datum
drandom(PG_FUNCTION_ARGS)
{
    float8      result;

    /* result [0.0 - 1.0) */
    result = (double) random() / ((double) MAX_RANDOM_VALUE + 1);

    PG_RETURN_FLOAT8(result);
}


/*
 *      setseed     - set seed for the random number generator
 */
Datum
setseed(PG_FUNCTION_ARGS)
{
    float8      seed = PG_GETARG_FLOAT8(0);
    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_VOID();
}



/*
 *      =========================
 *      FLOAT AGGREGATE OPERATORS
 *      =========================
 *
 *      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.
 *
 * Note that we represent N as a float to avoid having to build a special
 * datatype.  Given a reasonable floating-point implementation, there should
 * be no accuracy loss unless N exceeds 2 ^ 52 or so (by which time the
 * user will have doubtless lost interest anyway...)
 */

static float8 *
check_float8_array(ArrayType *transarray, const char *caller, int n)
{
    /*
     * 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 N float8 values.
     */
    if (ARR_NDIM(transarray) != 1 ||
        ARR_DIMS(transarray)[0] != n ||
        ARR_HASNULL(transarray) ||
        ARR_ELEMTYPE(transarray) != FLOAT8OID)
        elog(ERROR, "%s: expected %d-element float8 array", caller, n);
    return (float8 *) ARR_DATA_PTR(transarray);
}

/*
 * float8_combine
 *
 * An aggregate combine function used to combine two 3 fields
 * aggregate transition data into a single transition data.
 * This function is used only in two stage aggregation and
 * shouldn't be called outside aggregate context.
 */
Datum
float8_combine(PG_FUNCTION_ARGS)
{
    ArrayType  *transarray1 = PG_GETARG_ARRAYTYPE_P(0);
    ArrayType  *transarray2 = PG_GETARG_ARRAYTYPE_P(1);
    float8     *transvalues1;
    float8     *transvalues2;
    float8      N,
                sumX,
                sumX2;

    if (!AggCheckCallContext(fcinfo, NULL))
        elog(ERROR, "aggregate function called in non-aggregate context");

    transvalues1 = check_float8_array(transarray1, "float8_combine", 3);
    N = transvalues1[0];
    sumX = transvalues1[1];
    sumX2 = transvalues1[2];

    transvalues2 = check_float8_array(transarray2, "float8_combine", 3);

    N += transvalues2[0];
    sumX += transvalues2[1];
    CHECKFLOATVAL(sumX, isinf(transvalues1[1]) || isinf(transvalues2[1]),
                  true);
    sumX2 += transvalues2[2];
    CHECKFLOATVAL(sumX2, isinf(transvalues1[2]) || isinf(transvalues2[2]),
                  true);

    transvalues1[0] = N;
    transvalues1[1] = sumX;
    transvalues1[2] = sumX2;

    PG_RETURN_ARRAYTYPE_P(transarray1);
}

Datum
float8_accum(PG_FUNCTION_ARGS)
{
    ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
    float8      newval = PG_GETARG_FLOAT8(1);
    float8     *transvalues;
    float8      N,
                sumX,
                sumX2;

    transvalues = check_float8_array(transarray, "float8_accum", 3);
    N = transvalues[0];
    sumX = transvalues[1];
    sumX2 = transvalues[2];

    N += 1.0;
    sumX += newval;
    CHECKFLOATVAL(sumX, isinf(transvalues[1]) || isinf(newval), true);
    sumX2 += newval * newval;
    CHECKFLOATVAL(sumX2, isinf(transvalues[2]) || isinf(newval), true);

    /*
     * If we're invoked as an aggregate, we can cheat and modify our first
     * parameter in-place to reduce palloc overhead. Otherwise we construct a
     * new array with the updated transition data and return it.
     */
    if (AggCheckCallContext(fcinfo, NULL))
    {
        transvalues[0] = N;
        transvalues[1] = sumX;
        transvalues[2] = sumX2;

        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), FLOAT8PASSBYVAL, 'd');

        PG_RETURN_ARRAYTYPE_P(result);
    }
}

Datum
float4_accum(PG_FUNCTION_ARGS)
{
    ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);

    /* do computations as float8 */
    float8      newval = PG_GETARG_FLOAT4(1);
    float8     *transvalues;
    float8      N,
                sumX,
                sumX2;

    transvalues = check_float8_array(transarray, "float4_accum", 3);
    N = transvalues[0];
    sumX = transvalues[1];
    sumX2 = transvalues[2];

    N += 1.0;
    sumX += newval;
    CHECKFLOATVAL(sumX, isinf(transvalues[1]) || isinf(newval), true);
    sumX2 += newval * newval;
    CHECKFLOATVAL(sumX2, isinf(transvalues[2]) || isinf(newval), true);

    /*
     * If we're invoked as an aggregate, we can cheat and modify our first
     * parameter in-place to reduce palloc overhead. Otherwise we construct a
     * new array with the updated transition data and return it.
     */
    if (AggCheckCallContext(fcinfo, NULL))
    {
        transvalues[0] = N;
        transvalues[1] = sumX;
        transvalues[2] = sumX2;

        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), FLOAT8PASSBYVAL, 'd');

        PG_RETURN_ARRAYTYPE_P(result);
    }
}

Datum
float8_avg(PG_FUNCTION_ARGS)
{
    ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
    float8     *transvalues;
    float8      N,
                sumX;

    transvalues = check_float8_array(transarray, "float8_avg", 3);
    N = transvalues[0];
    sumX = transvalues[1];
    /* ignore sumX2 */

    /* SQL defines AVG of no values to be NULL */
    if (N == 0.0)
        PG_RETURN_NULL();

    PG_RETURN_FLOAT8(sumX / N);
}

Datum
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;
    float8      N,
                sumX,
                sumX2,
                numerator;

    transvalues = check_float8_array(transarray, "float8_var_samp", 3);
    N = transvalues[0];
    sumX = transvalues[1];
    sumX2 = transvalues[2];

    /* Sample variance is undefined when N is 0 or 1, so return NULL */
    if (N <= 1.0)
        PG_RETURN_NULL();

    numerator = N * sumX2 - sumX * sumX;
    CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true);

    /* Watch out for roundoff error producing a negative numerator */
    if (numerator <= 0.0)
        PG_RETURN_FLOAT8(0.0);

    PG_RETURN_FLOAT8(numerator / (N * (N - 1.0)));
}

Datum
float8_stddev_pop(PG_FUNCTION_ARGS)
{
    ArrayType  *transarray = PG_GETARG_ARRAYTYPE_P(0);
    float8     *transvalues;
    float8      N,
                sumX,
                sumX2,
                numerator;

    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];

    /* Sample stddev is undefined when N is 0 or 1, so return NULL */
    if (N <= 1.0)
        PG_RETURN_NULL();

    numerator = N * sumX2 - sumX * sumX;
    CHECKFLOATVAL(numerator, isinf(sumX2) || isinf(sumX), true);

    /* Watch out for roundoff error producing a negative numerator */
    if (numerator <= 0.0)
        PG_RETURN_FLOAT8(0.0);

    PG_RETURN_FLOAT8(sqrt(numerator / (N * (N - 1.0))));
}

/*
 *      =========================
 *      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);
    }
}

/*
 * float8_regr_combine
 *
 * An aggregate combine function used to combine two 6 fields
 * aggregate transition data into a single transition data.
 * This function is used only in two stage aggregation and
 * shouldn't be called outside aggregate context.
 */
Datum
float8_regr_combine(PG_FUNCTION_ARGS)
{
    ArrayType  *transarray1 = PG_GETARG_ARRAYTYPE_P(0);
    ArrayType  *transarray2 = PG_GETARG_ARRAYTYPE_P(1);
    float8     *transvalues1;
    float8     *transvalues2;
    float8      N,
                sumX,
                sumX2,
                sumY,
                sumY2,
                sumXY;

    if (!AggCheckCallContext(fcinfo, NULL))
        elog(ERROR, "aggregate function called in non-aggregate context");

    transvalues1 = check_float8_array(transarray1, "float8_regr_combine", 6);
    N = transvalues1[0];
    sumX = transvalues1[1];
    sumX2 = transvalues1[2];
    sumY = transvalues1[3];
    sumY2 = transvalues1[4];
    sumXY = transvalues1[5];

    transvalues2 = check_float8_array(transarray2, "float8_regr_combine", 6);

    N += transvalues2[0];
    sumX += transvalues2[1];
    CHECKFLOATVAL(sumX, isinf(transvalues1[1]) || isinf(transvalues2[1]),
                  true);
    sumX2 += transvalues2[2];
    CHECKFLOATVAL(sumX2, isinf(transvalues1[2]) || isinf(transvalues2[2]),
                  true);
    sumY += transvalues2[3];
    CHECKFLOATVAL(sumY, isinf(transvalues1[3]) || isinf(transvalues2[3]),
                  true);
    sumY2 += transvalues2[4];
    CHECKFLOATVAL(sumY2, isinf(transvalues1[4]) || isinf(transvalues2[4]),
                  true);
    sumXY += transvalues2[5];
    CHECKFLOATVAL(sumXY, isinf(transvalues1[5]) || isinf(transvalues2[5]),
                  true);

    transvalues1[0] = N;
    transvalues1[1] = sumX;
    transvalues1[2] = sumX2;
    transvalues1[3] = sumY;
    transvalues1[4] = sumY2;
    transvalues1[5] = sumXY;

    PG_RETURN_ARRAYTYPE_P(transarray1);
}


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);
}


/*
 *      ====================================
 *      MIXED-PRECISION ARITHMETIC OPERATORS
 *      ====================================
 */

/*
 *      float48pl       - returns arg1 + arg2
 *      float48mi       - returns arg1 - arg2
 *      float48mul      - returns arg1 * arg2
 *      float48div      - returns arg1 / arg2
 */
Datum
float48pl(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);
    float8      result;

    result = arg1 + arg2;
    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
    PG_RETURN_FLOAT8(result);
}

Datum
float48mi(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);
    float8      result;

    result = arg1 - arg2;
    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
    PG_RETURN_FLOAT8(result);
}

Datum
float48mul(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);
    float8      result;

    result = arg1 * arg2;
    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2),
                  arg1 == 0 || arg2 == 0);
    PG_RETURN_FLOAT8(result);
}

Datum
float48div(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);
    float8      result;

    if (arg2 == 0.0)
        ereport(ERROR,
                (errcode(ERRCODE_DIVISION_BY_ZERO),
                 errmsg("division by zero")));

    result = arg1 / arg2;
    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
    PG_RETURN_FLOAT8(result);
}

/*
 *      float84pl       - returns arg1 + arg2
 *      float84mi       - returns arg1 - arg2
 *      float84mul      - returns arg1 * arg2
 *      float84div      - returns arg1 / arg2
 */
Datum
float84pl(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);
    float8      result;

    result = arg1 + arg2;

    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
    PG_RETURN_FLOAT8(result);
}

Datum
float84mi(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);
    float8      result;

    result = arg1 - arg2;

    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), true);
    PG_RETURN_FLOAT8(result);
}

Datum
float84mul(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);
    float8      result;

    result = arg1 * arg2;

    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2),
                  arg1 == 0 || arg2 == 0);
    PG_RETURN_FLOAT8(result);
}

Datum
float84div(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);
    float8      result;

    if (arg2 == 0.0)
        ereport(ERROR,
                (errcode(ERRCODE_DIVISION_BY_ZERO),
                 errmsg("division by zero")));

    result = arg1 / arg2;

    CHECKFLOATVAL(result, isinf(arg1) || isinf(arg2), arg1 == 0);
    PG_RETURN_FLOAT8(result);
}

/*
 *      ====================
 *      COMPARISON OPERATORS
 *      ====================
 */

/*
 *      float48{eq,ne,lt,le,gt,ge}      - float4/float8 comparison operations
 */
Datum
float48eq(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) == 0);
}

Datum
float48ne(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) != 0);
}

Datum
float48lt(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) < 0);
}

Datum
float48le(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) <= 0);
}

Datum
float48gt(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) > 0);
}

Datum
float48ge(PG_FUNCTION_ARGS)
{
    float4      arg1 = PG_GETARG_FLOAT4(0);
    float8      arg2 = PG_GETARG_FLOAT8(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) >= 0);
}

/*
 *      float84{eq,ne,lt,le,gt,ge}      - float8/float4 comparison operations
 */
Datum
float84eq(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) == 0);
}

Datum
float84ne(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) != 0);
}

Datum
float84lt(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) < 0);
}

Datum
float84le(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) <= 0);
}

Datum
float84gt(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) > 0);
}

Datum
float84ge(PG_FUNCTION_ARGS)
{
    float8      arg1 = PG_GETARG_FLOAT8(0);
    float4      arg2 = PG_GETARG_FLOAT4(1);

    PG_RETURN_BOOL(float8_cmp_internal(arg1, arg2) >= 0);
}

/*
 * 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      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;
}

#endif   /* !HAVE_CBRT */