/*------------------------------------------------------------------------- * * float.c * Functions for the built-in floating-point types. * * Portions Copyright (c) 1996-2025, 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 #include #include #include #include "catalog/pg_type.h" #include "common/int.h" #include "common/shortest_dec.h" #include "libpq/pqformat.h" #include "utils/array.h" #include "utils/float.h" #include "utils/fmgrprotos.h" #include "utils/sortsupport.h" /* * Configurable GUC parameter * * If >0, use shortest-decimal format for output; this is both the default and * allows for compatibility with clients that explicitly set a value here to * get round-trip-accurate results. If 0 or less, then use the old, slow, * decimal rounding method. */ int extra_float_digits = 1; /* 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(). * * The additional extern declarations are to silence * -Wmissing-variable-declarations. */ extern float8 degree_c_thirty; extern float8 degree_c_forty_five; extern float8 degree_c_sixty; extern float8 degree_c_one_half; extern float8 degree_c_one; 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); /* * We use these out-of-line ereport() calls to report float overflow, * underflow, and zero-divide, because following our usual practice of * repeating them at each call site would lead to a lot of code bloat. * * This does mean that you don't get a useful error location indicator. */ pg_noinline void float_overflow_error(void) { ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("value out of range: overflow"))); } pg_noinline void float_underflow_error(void) { ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("value out of range: underflow"))); } pg_noinline void float_zero_divide_error(void) { ereport(ERROR, (errcode(ERRCODE_DIVISION_BY_ZERO), errmsg("division by zero"))); } /* * 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 * * Note that this code now uses strtof(), where it used to use strtod(). * * The motivation for using strtof() is to avoid a double-rounding problem: * for certain decimal inputs, if you round the input correctly to a double, * and then round the double to a float, the result is incorrect in that it * does not match the result of rounding the decimal value to float directly. * * One of the best examples is 7.038531e-26: * * 0xAE43FDp-107 = 7.03853069185120912085...e-26 * midpoint 7.03853100000000022281...e-26 * 0xAE43FEp-107 = 7.03853130814879132477...e-26 * * making 0xAE43FDp-107 the correct float result, but if you do the conversion * via a double, you get * * 0xAE43FD.7FFFFFF8p-107 = 7.03853099999999907487...e-26 * midpoint 7.03853099999999964884...e-26 * 0xAE43FD.80000000p-107 = 7.03853100000000022281...e-26 * 0xAE43FD.80000008p-107 = 7.03853100000000137076...e-26 * * so the value rounds to the double exactly on the midpoint between the two * nearest floats, and then rounding again to a float gives the incorrect * result of 0xAE43FEp-107. * */ Datum float4in(PG_FUNCTION_ARGS) { char *num = PG_GETARG_CSTRING(0); PG_RETURN_FLOAT4(float4in_internal(num, NULL, "real", num, fcinfo->context)); } /* * float4in_internal - guts of float4in() * * This is exposed for use by functions that want a reasonably * platform-independent way of inputting floats. The behavior is * essentially like strtof + ereturn on error. * * Uses the same API as float8in_internal below, so most of its * comments also apply here, except regarding use in geometric types. */ float4 float4in_internal(char *num, char **endptr_p, const char *type_name, const char *orig_string, struct Node *escontext) { float val; char *endptr; /* * endptr points to the first character _after_ the sequence we recognized * as a valid floating point number. orig_string points to the original * input string. */ /* 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') ereturn(escontext, 0, (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), errmsg("invalid input syntax for type %s: \"%s\"", type_name, orig_string))); errno = 0; val = strtof(num, &endptr); /* did we not see anything that looks like a double? */ if (endptr == num || errno != 0) { int save_errno = errno; /* * C99 requires that strtof() 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 strtof() 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 strtof() 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 || #if !defined(HUGE_VALF) isinf(val) #else (val >= HUGE_VALF || val <= -HUGE_VALF) #endif ) { /* see comments in float8in_internal for rationale */ char *errnumber = pstrdup(num); errnumber[endptr - num] = '\0'; ereturn(escontext, 0, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("\"%s\" is out of range for type real", errnumber))); } } else ereturn(escontext, 0, (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), errmsg("invalid input syntax for type %s: \"%s\"", type_name, orig_string))); } /* 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') ereturn(escontext, 0, (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), errmsg("invalid input syntax for type %s: \"%s\"", type_name, orig_string))); return 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(32); int ndig = FLT_DIG + extra_float_digits; if (extra_float_digits > 0) { float_to_shortest_decimal_buf(num, ascii); PG_RETURN_CSTRING(ascii); } (void) pg_strfromd(ascii, 32, 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, fcinfo->context)); } /* * 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 + ereturn on error, but note the following * differences: * 1. Both leading and trailing whitespace are skipped. * 2. If endptr_p is NULL, we report 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. * * If escontext points to an ErrorSaveContext node, that is filled instead * of throwing an error; the caller must check SOFT_ERROR_OCCURRED() * to detect errors. * * "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. */ float8 float8in_internal(char *num, char **endptr_p, const char *type_name, const char *orig_string, struct Node *escontext) { 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') ereturn(escontext, 0, (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'; ereturn(escontext, 0, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("\"%s\" is out of range for type double precision", errnumber))); } } else ereturn(escontext, 0, (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), errmsg("invalid input syntax for type %s: \"%s\"", type_name, orig_string))); } /* 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') ereturn(escontext, 0, (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(32); int ndig = DBL_DIG + extra_float_digits; if (extra_float_digits > 0) { double_to_shortest_decimal_buf(num, ascii); return ascii; } (void) pg_strfromd(ascii, 32, 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(fabsf(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_gt(arg1, arg2)) result = arg1; else result = arg2; PG_RETURN_FLOAT4(result); } Datum float4smaller(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float4 arg2 = PG_GETARG_FLOAT4(1); float4 result; if (float4_lt(arg1, arg2)) 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_gt(arg1, arg2)) result = arg1; else result = arg2; PG_RETURN_FLOAT8(result); } Datum float8smaller(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 arg2 = PG_GETARG_FLOAT8(1); float8 result; if (float8_lt(arg1, arg2)) 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); PG_RETURN_FLOAT4(float4_pl(arg1, arg2)); } Datum float4mi(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_FLOAT4(float4_mi(arg1, arg2)); } Datum float4mul(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_FLOAT4(float4_mul(arg1, arg2)); } Datum float4div(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_FLOAT4(float4_div(arg1, arg2)); } /* * float8pl - returns arg1 + arg2 * float8mi - returns arg1 - arg2 * float8mul - returns arg1 * arg2 * float8div - returns arg1 / arg2 */ Datum float8pl(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_FLOAT8(float8_pl(arg1, arg2)); } Datum float8mi(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_FLOAT8(float8_mi(arg1, arg2)); } Datum float8mul(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_FLOAT8(float8_mul(arg1, arg2)); } Datum float8div(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_FLOAT8(float8_div(arg1, arg2)); } /* * ==================== * COMPARISON OPERATORS * ==================== */ /* * float4{eq,ne,lt,le,gt,ge} - float4/float4 comparison operations */ int float4_cmp_internal(float4 a, float4 b) { if (float4_gt(a, b)) return 1; if (float4_lt(a, b)) return -1; return 0; } Datum float4eq(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_BOOL(float4_eq(arg1, arg2)); } Datum float4ne(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_BOOL(float4_ne(arg1, arg2)); } Datum float4lt(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_BOOL(float4_lt(arg1, arg2)); } Datum float4le(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_BOOL(float4_le(arg1, arg2)); } Datum float4gt(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_BOOL(float4_gt(arg1, arg2)); } Datum float4ge(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_BOOL(float4_ge(arg1, arg2)); } Datum btfloat4cmp(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_INT32(float4_cmp_internal(arg1, arg2)); } 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) { if (float8_gt(a, b)) return 1; if (float8_lt(a, b)) return -1; return 0; } Datum float8eq(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_BOOL(float8_eq(arg1, arg2)); } Datum float8ne(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_BOOL(float8_ne(arg1, arg2)); } Datum float8lt(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_BOOL(float8_lt(arg1, arg2)); } Datum float8le(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_BOOL(float8_le(arg1, arg2)); } Datum float8gt(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_BOOL(float8_gt(arg1, arg2)); } Datum float8ge(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_BOOL(float8_ge(arg1, arg2)); } Datum btfloat8cmp(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_INT32(float8_cmp_internal(arg1, arg2)); } 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)); } /* * in_range support function for float8. * * Note: we needn't supply a float8_float4 variant, as implicit coercion * of the offset value takes care of that scenario just as well. */ Datum in_range_float8_float8(PG_FUNCTION_ARGS) { float8 val = PG_GETARG_FLOAT8(0); float8 base = PG_GETARG_FLOAT8(1); float8 offset = PG_GETARG_FLOAT8(2); bool sub = PG_GETARG_BOOL(3); bool less = PG_GETARG_BOOL(4); float8 sum; /* * Reject negative or NaN offset. Negative is per spec, and NaN is * because appropriate semantics for that seem non-obvious. */ if (isnan(offset) || offset < 0) ereport(ERROR, (errcode(ERRCODE_INVALID_PRECEDING_OR_FOLLOWING_SIZE), errmsg("invalid preceding or following size in window function"))); /* * Deal with cases where val and/or base is NaN, following the rule that * NaN sorts after non-NaN (cf float8_cmp_internal). The offset cannot * affect the conclusion. */ if (isnan(val)) { if (isnan(base)) PG_RETURN_BOOL(true); /* NAN = NAN */ else PG_RETURN_BOOL(!less); /* NAN > non-NAN */ } else if (isnan(base)) { PG_RETURN_BOOL(less); /* non-NAN < NAN */ } /* * Deal with cases where both base and offset are infinite, and computing * base +/- offset would produce NaN. This corresponds to a window frame * whose boundary infinitely precedes +inf or infinitely follows -inf, * which is not well-defined. For consistency with other cases involving * infinities, such as the fact that +inf infinitely follows +inf, we * choose to assume that +inf infinitely precedes +inf and -inf infinitely * follows -inf, and therefore that all finite and infinite values are in * such a window frame. * * offset is known positive, so we need only check the sign of base in * this test. */ if (isinf(offset) && isinf(base) && (sub ? base > 0 : base < 0)) PG_RETURN_BOOL(true); /* * Otherwise it should be safe to compute base +/- offset. We trust the * FPU to cope if an input is +/-inf or the true sum would overflow, and * produce a suitably signed infinity, which will compare properly against * val whether or not that's infinity. */ if (sub) sum = base - offset; else sum = base + offset; if (less) PG_RETURN_BOOL(val <= sum); else PG_RETURN_BOOL(val >= sum); } /* * in_range support function for float4. * * We would need a float4_float8 variant in any case, so we supply that and * let implicit coercion take care of the float4_float4 case. */ Datum in_range_float4_float8(PG_FUNCTION_ARGS) { float4 val = PG_GETARG_FLOAT4(0); float4 base = PG_GETARG_FLOAT4(1); float8 offset = PG_GETARG_FLOAT8(2); bool sub = PG_GETARG_BOOL(3); bool less = PG_GETARG_BOOL(4); float8 sum; /* * Reject negative or NaN offset. Negative is per spec, and NaN is * because appropriate semantics for that seem non-obvious. */ if (isnan(offset) || offset < 0) ereport(ERROR, (errcode(ERRCODE_INVALID_PRECEDING_OR_FOLLOWING_SIZE), errmsg("invalid preceding or following size in window function"))); /* * Deal with cases where val and/or base is NaN, following the rule that * NaN sorts after non-NaN (cf float8_cmp_internal). The offset cannot * affect the conclusion. */ if (isnan(val)) { if (isnan(base)) PG_RETURN_BOOL(true); /* NAN = NAN */ else PG_RETURN_BOOL(!less); /* NAN > non-NAN */ } else if (isnan(base)) { PG_RETURN_BOOL(less); /* non-NAN < NAN */ } /* * Deal with cases where both base and offset are infinite, and computing * base +/- offset would produce NaN. This corresponds to a window frame * whose boundary infinitely precedes +inf or infinitely follows -inf, * which is not well-defined. For consistency with other cases involving * infinities, such as the fact that +inf infinitely follows +inf, we * choose to assume that +inf infinitely precedes +inf and -inf infinitely * follows -inf, and therefore that all finite and infinite values are in * such a window frame. * * offset is known positive, so we need only check the sign of base in * this test. */ if (isinf(offset) && isinf(base) && (sub ? base > 0 : base < 0)) PG_RETURN_BOOL(true); /* * Otherwise it should be safe to compute base +/- offset. We trust the * FPU to cope if an input is +/-inf or the true sum would overflow, and * produce a suitably signed infinity, which will compare properly against * val whether or not that's infinity. */ if (sub) sum = base - offset; else sum = base + offset; if (less) PG_RETURN_BOOL(val <= sum); else PG_RETURN_BOOL(val >= sum); } /* * =================== * 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); float4 result; result = (float4) num; if (unlikely(isinf(result)) && !isinf(num)) float_overflow_error(); if (unlikely(result == 0.0f) && num != 0.0) float_underflow_error(); PG_RETURN_FLOAT4(result); } /* * dtoi4 - converts a float8 number to an int4 number */ Datum dtoi4(PG_FUNCTION_ARGS) { float8 num = PG_GETARG_FLOAT8(0); /* * Get rid of any fractional part in the input. This is so we don't fail * on just-out-of-range values that would round into range. Note * assumption that rint() will pass through a NaN or Inf unchanged. */ num = rint(num); /* Range check */ if (unlikely(isnan(num) || !FLOAT8_FITS_IN_INT32(num))) ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("integer out of range"))); PG_RETURN_INT32((int32) num); } /* * dtoi2 - converts a float8 number to an int2 number */ Datum dtoi2(PG_FUNCTION_ARGS) { float8 num = PG_GETARG_FLOAT8(0); /* * Get rid of any fractional part in the input. This is so we don't fail * on just-out-of-range values that would round into range. Note * assumption that rint() will pass through a NaN or Inf unchanged. */ num = rint(num); /* Range check */ if (unlikely(isnan(num) || !FLOAT8_FITS_IN_INT16(num))) ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("smallint out of range"))); PG_RETURN_INT16((int16) 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); /* * Get rid of any fractional part in the input. This is so we don't fail * on just-out-of-range values that would round into range. Note * assumption that rint() will pass through a NaN or Inf unchanged. */ num = rint(num); /* Range check */ if (unlikely(isnan(num) || !FLOAT4_FITS_IN_INT32(num))) ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("integer out of range"))); PG_RETURN_INT32((int32) num); } /* * ftoi2 - converts a float4 number to an int2 number */ Datum ftoi2(PG_FUNCTION_ARGS) { float4 num = PG_GETARG_FLOAT4(0); /* * Get rid of any fractional part in the input. This is so we don't fail * on just-out-of-range values that would round into range. Note * assumption that rint() will pass through a NaN or Inf unchanged. */ num = rint(num); /* Range check */ if (unlikely(isnan(num) || !FLOAT4_FITS_IN_INT16(num))) ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("smallint out of range"))); PG_RETURN_INT16((int16) 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); if (unlikely(isinf(result)) && !isinf(arg1)) float_overflow_error(); if (unlikely(result == 0.0) && arg1 != 0.0) float_underflow_error(); 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); if (unlikely(isinf(result)) && !isinf(arg1)) float_overflow_error(); if (unlikely(result == 0.0) && arg1 != 0.0) float_underflow_error(); 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 POSIX spec says that NaN ^ 0 = 1, and 1 ^ NaN = 1, while all other * cases with NaN inputs yield NaN (with no error). Many older platforms * get one or more of these cases wrong, so deal with them via explicit * logic rather than trusting pow(3). */ if (isnan(arg1)) { if (isnan(arg2) || arg2 != 0.0) PG_RETURN_FLOAT8(get_float8_nan()); PG_RETURN_FLOAT8(1.0); } if (isnan(arg2)) { if (arg1 != 1.0) PG_RETURN_FLOAT8(get_float8_nan()); PG_RETURN_FLOAT8(1.0); } /* * 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"))); /* * We don't trust the platform's pow() to handle infinity cases per POSIX * spec either, so deal with those explicitly too. It's easier to handle * infinite y first, so that it doesn't matter if x is also infinite. */ if (isinf(arg2)) { float8 absx = fabs(arg1); if (absx == 1.0) result = 1.0; else if (arg2 > 0.0) /* y = +Inf */ { if (absx > 1.0) result = arg2; else result = 0.0; } else /* y = -Inf */ { if (absx > 1.0) result = 0.0; else result = -arg2; } } else if (isinf(arg1)) { if (arg2 == 0.0) result = 1.0; else if (arg1 > 0.0) /* x = +Inf */ { if (arg2 > 0.0) result = arg1; else result = 0.0; } else /* x = -Inf */ { /* * Per POSIX, the sign of the result depends on whether y is an * odd integer. Since x < 0, we already know from the previous * domain check that y is an integer. It is odd if y/2 is not * also an integer. */ float8 halfy = arg2 / 2; /* should be computed exactly */ bool yisoddinteger = (floor(halfy) != halfy); if (arg2 > 0.0) result = yisoddinteger ? arg1 : -arg1; else result = yisoddinteger ? -0.0 : 0.0; } } else { /* * pow() sets errno on only some platforms, depending on whether it * follows _IEEE_, _POSIX_, _XOPEN_, or _SVID_, so we must check both * errno and invalid output values. (We can't rely on just the * latter, either; some old platforms return a large-but-finite * HUGE_VAL when reporting overflow.) */ errno = 0; result = pow(arg1, arg2); if (errno == EDOM || isnan(result)) { /* * We handled all possible domain errors above, so this should be * impossible. However, old glibc versions on x86 have a bug that * causes them to fail this way for abs(y) greater than 2^63: * * https://sourceware.org/bugzilla/show_bug.cgi?id=3866 * * Hence, if we get here, assume y is finite but large (large * enough to be certainly even). The result should be 0 if x == 0, * 1.0 if abs(x) == 1.0, otherwise an overflow or underflow error. */ if (arg1 == 0.0) result = 0.0; /* we already verified y is positive */ else { float8 absx = fabs(arg1); if (absx == 1.0) result = 1.0; else if (arg2 >= 0.0 ? (absx > 1.0) : (absx < 1.0)) float_overflow_error(); else float_underflow_error(); } } else if (errno == ERANGE) { if (result != 0.0) float_overflow_error(); else float_underflow_error(); } else { if (unlikely(isinf(result))) float_overflow_error(); if (unlikely(result == 0.0) && arg1 != 0.0) float_underflow_error(); } } PG_RETURN_FLOAT8(result); } /* * dexp - returns the exponential function of arg1 */ Datum dexp(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 result; /* * Handle NaN and Inf cases explicitly. This avoids needing to assume * that the platform's exp() conforms to POSIX for these cases, and it * removes some edge cases for the overflow checks below. */ if (isnan(arg1)) result = arg1; else if (isinf(arg1)) { /* Per POSIX, exp(-Inf) is 0 */ result = (arg1 > 0.0) ? arg1 : 0; } else { /* * On some platforms, exp() will not set errno but just return Inf or * zero to report overflow/underflow; therefore, test both cases. */ errno = 0; result = exp(arg1); if (unlikely(errno == ERANGE)) { if (result != 0.0) float_overflow_error(); else float_underflow_error(); } else if (unlikely(isinf(result))) float_overflow_error(); else if (unlikely(result == 0.0)) float_underflow_error(); } 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); if (unlikely(isinf(result)) && !isinf(arg1)) float_overflow_error(); if (unlikely(result == 0.0) && arg1 != 1.0) float_underflow_error(); 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); if (unlikely(isinf(result)) && !isinf(arg1)) float_overflow_error(); if (unlikely(result == 0.0) && arg1 != 1.0) float_underflow_error(); 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); if (unlikely(isinf(result))) float_overflow_error(); 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); if (unlikely(isinf(result))) float_overflow_error(); 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); if (unlikely(isinf(result))) float_overflow_error(); 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); if (unlikely(isinf(result))) float_overflow_error(); 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"))); if (unlikely(isinf(result))) float_overflow_error(); 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; /* Not checking for overflow because cot(0) == Inf */ 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"))); if (unlikely(isinf(result))) float_overflow_error(); 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"))); /* Not checking for overflow because tan(pi/2) == Inf */ 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); if (unlikely(isinf(result))) float_overflow_error(); 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); if (unlikely(isinf(result))) float_overflow_error(); 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; if (unlikely(isinf(result))) float_overflow_error(); 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; if (unlikely(isinf(result))) float_overflow_error(); 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); if (unlikely(isinf(result))) float_overflow_error(); 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; /* Not checking for overflow because cotd(0) == Inf */ 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); if (unlikely(isinf(result))) float_overflow_error(); 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; /* Not checking for overflow because tand(90) == Inf */ PG_RETURN_FLOAT8(result); } /* * degrees - returns degrees converted from radians */ Datum degrees(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); PG_RETURN_FLOAT8(float8_div(arg1, RADIANS_PER_DEGREE)); } /* * dpi - returns the constant PI */ Datum dpi(PG_FUNCTION_ARGS) { PG_RETURN_FLOAT8(M_PI); } /* * radians - returns radians converted from degrees */ Datum radians(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); PG_RETURN_FLOAT8(float8_mul(arg1, RADIANS_PER_DEGREE)); } /* ========== HYPERBOLIC FUNCTIONS ========== */ /* * dsinh - returns the hyperbolic sine of arg1 */ Datum dsinh(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 result; errno = 0; result = sinh(arg1); /* * if an ERANGE error occurs, it means there is an overflow. For sinh, * the result should be either -infinity or infinity, depending on the * sign of arg1. */ if (errno == ERANGE) { if (arg1 < 0) result = -get_float8_infinity(); else result = get_float8_infinity(); } PG_RETURN_FLOAT8(result); } /* * dcosh - returns the hyperbolic cosine of arg1 */ Datum dcosh(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 result; errno = 0; result = cosh(arg1); /* * if an ERANGE error occurs, it means there is an overflow. As cosh is * always positive, it always means the result is positive infinity. */ if (errno == ERANGE) result = get_float8_infinity(); if (unlikely(result == 0.0)) float_underflow_error(); PG_RETURN_FLOAT8(result); } /* * dtanh - returns the hyperbolic tangent of arg1 */ Datum dtanh(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 result; /* * For tanh, we don't need an errno check because it never overflows. */ result = tanh(arg1); if (unlikely(isinf(result))) float_overflow_error(); PG_RETURN_FLOAT8(result); } /* * dasinh - returns the inverse hyperbolic sine of arg1 */ Datum dasinh(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 result; /* * For asinh, we don't need an errno check because it never overflows. */ result = asinh(arg1); PG_RETURN_FLOAT8(result); } /* * dacosh - returns the inverse hyperbolic cosine of arg1 */ Datum dacosh(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 result; /* * acosh is only defined for inputs >= 1.0. By checking this ourselves, * we need not worry about checking for an EDOM error, which is a good * thing because some implementations will report that for NaN. Otherwise, * no error is possible. */ if (arg1 < 1.0) ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("input is out of range"))); result = acosh(arg1); PG_RETURN_FLOAT8(result); } /* * datanh - returns the inverse hyperbolic tangent of arg1 */ Datum datanh(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 result; /* * atanh is only defined for inputs between -1 and 1. By checking this * ourselves, we need not worry about checking for an EDOM error, which is * a good thing because some implementations will report that for NaN. */ if (arg1 < -1.0 || arg1 > 1.0) ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("input is out of range"))); /* * Also handle the infinity cases ourselves; this is helpful because old * glibc versions may produce the wrong errno for this. All other inputs * cannot produce an error. */ if (arg1 == -1.0) result = -get_float8_infinity(); else if (arg1 == 1.0) result = get_float8_infinity(); else result = atanh(arg1); PG_RETURN_FLOAT8(result); } /* ========== ERROR FUNCTIONS ========== */ /* * derf - returns the error function: erf(arg1) */ Datum derf(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 result; /* * For erf, we don't need an errno check because it never overflows. */ result = erf(arg1); if (unlikely(isinf(result))) float_overflow_error(); PG_RETURN_FLOAT8(result); } /* * derfc - returns the complementary error function: 1 - erf(arg1) */ Datum derfc(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 result; /* * For erfc, we don't need an errno check because it never overflows. */ result = erfc(arg1); if (unlikely(isinf(result))) float_overflow_error(); PG_RETURN_FLOAT8(result); } /* ========== GAMMA FUNCTIONS ========== */ /* * dgamma - returns the gamma function of arg1 */ Datum dgamma(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 result; /* * Handle NaN and Inf cases explicitly. This simplifies the overflow * checks on platforms that do not set errno. */ if (isnan(arg1)) result = arg1; else if (isinf(arg1)) { /* Per POSIX, an input of -Inf causes a domain error */ if (arg1 < 0) { float_overflow_error(); result = get_float8_nan(); /* keep compiler quiet */ } else result = arg1; } else { /* * Note: the POSIX/C99 gamma function is called "tgamma", not "gamma". * * On some platforms, tgamma() will not set errno but just return Inf, * NaN, or zero to report overflow/underflow; therefore, test those * cases explicitly (note that, like the exponential function, the * gamma function has no zeros). */ errno = 0; result = tgamma(arg1); if (errno != 0 || isinf(result) || isnan(result)) { if (result != 0.0) float_overflow_error(); else float_underflow_error(); } else if (result == 0.0) float_underflow_error(); } PG_RETURN_FLOAT8(result); } /* * dlgamma - natural logarithm of absolute value of gamma of arg1 */ Datum dlgamma(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float8 result; /* * Note: lgamma may not be thread-safe because it may write to a global * variable signgam, which may not be thread-local. However, this doesn't * matter to us, since we don't use signgam. */ errno = 0; result = lgamma(arg1); /* * If an ERANGE error occurs, it means there was an overflow or a pole * error (which happens for zero and negative integer inputs). * * On some platforms, lgamma() will not set errno but just return infinity * to report overflow, but it should never underflow. */ if (errno == ERANGE || (isinf(result) && !isinf(arg1))) float_overflow_error(); PG_RETURN_FLOAT8(result); } /* * ========================= * 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 naive schoolbook implementation of these aggregates works by * accumulating sum(X) and sum(X^2). However, this approach suffers from * large rounding errors in the final computation of quantities like the * population variance (N*sum(X^2) - sum(X)^2) / N^2, since each of the * intermediate terms is potentially very large, while the difference is often * quite small. * * Instead we use the Youngs-Cramer algorithm [1] which works by accumulating * Sx=sum(X) and Sxx=sum((X-Sx/N)^2), using a numerically stable algorithm to * incrementally update those quantities. The final computations of each of * the aggregate values is then trivial and gives more accurate results (for * example, the population variance is just Sxx/N). This algorithm is also * fairly easy to generalize to allow parallel execution without loss of * precision (see, for example, [2]). For more details, and a comparison of * this with other algorithms, see [3]. * * The transition datatype for all these aggregates is a 3-element array * of float8, holding the values N, Sx, Sxx 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...) * * [1] Some Results Relevant to Choice of Sum and Sum-of-Product Algorithms, * E. A. Youngs and E. M. Cramer, Technometrics Vol 13, No 3, August 1971. * * [2] Updating Formulae and a Pairwise Algorithm for Computing Sample * Variances, T. F. Chan, G. H. Golub & R. J. LeVeque, COMPSTAT 1982. * * [3] Numerically Stable Parallel Computation of (Co-)Variance, Erich * Schubert and Michael Gertz, Proceedings of the 30th International * Conference on Scientific and Statistical Database Management, 2018. */ 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 N1, Sx1, Sxx1, N2, Sx2, Sxx2, tmp, N, Sx, Sxx; transvalues1 = check_float8_array(transarray1, "float8_combine", 3); transvalues2 = check_float8_array(transarray2, "float8_combine", 3); N1 = transvalues1[0]; Sx1 = transvalues1[1]; Sxx1 = transvalues1[2]; N2 = transvalues2[0]; Sx2 = transvalues2[1]; Sxx2 = transvalues2[2]; /*-------------------- * The transition values combine using a generalization of the * Youngs-Cramer algorithm as follows: * * N = N1 + N2 * Sx = Sx1 + Sx2 * Sxx = Sxx1 + Sxx2 + N1 * N2 * (Sx1/N1 - Sx2/N2)^2 / N; * * It's worth handling the special cases N1 = 0 and N2 = 0 separately * since those cases are trivial, and we then don't need to worry about * division-by-zero errors in the general case. *-------------------- */ if (N1 == 0.0) { N = N2; Sx = Sx2; Sxx = Sxx2; } else if (N2 == 0.0) { N = N1; Sx = Sx1; Sxx = Sxx1; } else { N = N1 + N2; Sx = float8_pl(Sx1, Sx2); tmp = Sx1 / N1 - Sx2 / N2; Sxx = Sxx1 + Sxx2 + N1 * N2 * tmp * tmp / N; if (unlikely(isinf(Sxx)) && !isinf(Sxx1) && !isinf(Sxx2)) float_overflow_error(); } /* * 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)) { transvalues1[0] = N; transvalues1[1] = Sx; transvalues1[2] = Sxx; PG_RETURN_ARRAYTYPE_P(transarray1); } else { Datum transdatums[3]; ArrayType *result; transdatums[0] = Float8GetDatumFast(N); transdatums[1] = Float8GetDatumFast(Sx); transdatums[2] = Float8GetDatumFast(Sxx); result = construct_array_builtin(transdatums, 3, FLOAT8OID); PG_RETURN_ARRAYTYPE_P(result); } } Datum float8_accum(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 newval = PG_GETARG_FLOAT8(1); float8 *transvalues; float8 N, Sx, Sxx, tmp; transvalues = check_float8_array(transarray, "float8_accum", 3); N = transvalues[0]; Sx = transvalues[1]; Sxx = transvalues[2]; /* * Use the Youngs-Cramer algorithm to incorporate the new value into the * transition values. */ N += 1.0; Sx += newval; if (transvalues[0] > 0.0) { tmp = newval * N - Sx; Sxx += tmp * tmp / (N * transvalues[0]); /* * Overflow check. We only report an overflow error when finite * inputs lead to infinite results. Note also that Sxx should be NaN * if any of the inputs are infinite, so we intentionally prevent Sxx * from becoming infinite. */ if (isinf(Sx) || isinf(Sxx)) { if (!isinf(transvalues[1]) && !isinf(newval)) float_overflow_error(); Sxx = get_float8_nan(); } } else { /* * At the first input, we normally can leave Sxx as 0. However, if * the first input is Inf or NaN, we'd better force Sxx to NaN; * otherwise we will falsely report variance zero when there are no * more inputs. */ if (isnan(newval) || isinf(newval)) Sxx = get_float8_nan(); } /* * 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] = Sx; transvalues[2] = Sxx; PG_RETURN_ARRAYTYPE_P(transarray); } else { Datum transdatums[3]; ArrayType *result; transdatums[0] = Float8GetDatumFast(N); transdatums[1] = Float8GetDatumFast(Sx); transdatums[2] = Float8GetDatumFast(Sxx); result = construct_array_builtin(transdatums, 3, FLOAT8OID); 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, Sx, Sxx, tmp; transvalues = check_float8_array(transarray, "float4_accum", 3); N = transvalues[0]; Sx = transvalues[1]; Sxx = transvalues[2]; /* * Use the Youngs-Cramer algorithm to incorporate the new value into the * transition values. */ N += 1.0; Sx += newval; if (transvalues[0] > 0.0) { tmp = newval * N - Sx; Sxx += tmp * tmp / (N * transvalues[0]); /* * Overflow check. We only report an overflow error when finite * inputs lead to infinite results. Note also that Sxx should be NaN * if any of the inputs are infinite, so we intentionally prevent Sxx * from becoming infinite. */ if (isinf(Sx) || isinf(Sxx)) { if (!isinf(transvalues[1]) && !isinf(newval)) float_overflow_error(); Sxx = get_float8_nan(); } } else { /* * At the first input, we normally can leave Sxx as 0. However, if * the first input is Inf or NaN, we'd better force Sxx to NaN; * otherwise we will falsely report variance zero when there are no * more inputs. */ if (isnan(newval) || isinf(newval)) Sxx = get_float8_nan(); } /* * 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] = Sx; transvalues[2] = Sxx; PG_RETURN_ARRAYTYPE_P(transarray); } else { Datum transdatums[3]; ArrayType *result; transdatums[0] = Float8GetDatumFast(N); transdatums[1] = Float8GetDatumFast(Sx); transdatums[2] = Float8GetDatumFast(Sxx); result = construct_array_builtin(transdatums, 3, FLOAT8OID); PG_RETURN_ARRAYTYPE_P(result); } } Datum float8_avg(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 *transvalues; float8 N, Sx; transvalues = check_float8_array(transarray, "float8_avg", 3); N = transvalues[0]; Sx = transvalues[1]; /* ignore Sxx */ /* SQL defines AVG of no values to be NULL */ if (N == 0.0) PG_RETURN_NULL(); PG_RETURN_FLOAT8(Sx / N); } Datum float8_var_pop(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 *transvalues; float8 N, Sxx; transvalues = check_float8_array(transarray, "float8_var_pop", 3); N = transvalues[0]; /* ignore Sx */ Sxx = transvalues[2]; /* Population variance is undefined when N is 0, so return NULL */ if (N == 0.0) PG_RETURN_NULL(); /* Note that Sxx is guaranteed to be non-negative */ PG_RETURN_FLOAT8(Sxx / N); } Datum float8_var_samp(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 *transvalues; float8 N, Sxx; transvalues = check_float8_array(transarray, "float8_var_samp", 3); N = transvalues[0]; /* ignore Sx */ Sxx = transvalues[2]; /* Sample variance is undefined when N is 0 or 1, so return NULL */ if (N <= 1.0) PG_RETURN_NULL(); /* Note that Sxx is guaranteed to be non-negative */ PG_RETURN_FLOAT8(Sxx / (N - 1.0)); } Datum float8_stddev_pop(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 *transvalues; float8 N, Sxx; transvalues = check_float8_array(transarray, "float8_stddev_pop", 3); N = transvalues[0]; /* ignore Sx */ Sxx = transvalues[2]; /* Population stddev is undefined when N is 0, so return NULL */ if (N == 0.0) PG_RETURN_NULL(); /* Note that Sxx is guaranteed to be non-negative */ PG_RETURN_FLOAT8(sqrt(Sxx / N)); } Datum float8_stddev_samp(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 *transvalues; float8 N, Sxx; transvalues = check_float8_array(transarray, "float8_stddev_samp", 3); N = transvalues[0]; /* ignore Sx */ Sxx = transvalues[2]; /* Sample stddev is undefined when N is 0 or 1, so return NULL */ if (N <= 1.0) PG_RETURN_NULL(); /* Note that Sxx is guaranteed to be non-negative */ PG_RETURN_FLOAT8(sqrt(Sxx / (N - 1.0))); } /* * ========================= * SQL2003 BINARY AGGREGATES * ========================= * * As with the preceding aggregates, we use the Youngs-Cramer algorithm to * reduce rounding errors in the aggregate final functions. * * The transition datatype for all these aggregates is a 6-element array of * float8, holding the values N, Sx=sum(X), Sxx=sum((X-Sx/N)^2), Sy=sum(Y), * Syy=sum((Y-Sy/N)^2), Sxy=sum((X-Sx/N)*(Y-Sy/N)) in that order. * * Note that Y is the first argument to all these 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, Sx, Sxx, Sy, Syy, Sxy, tmpX, tmpY, scale; transvalues = check_float8_array(transarray, "float8_regr_accum", 6); N = transvalues[0]; Sx = transvalues[1]; Sxx = transvalues[2]; Sy = transvalues[3]; Syy = transvalues[4]; Sxy = transvalues[5]; /* * Use the Youngs-Cramer algorithm to incorporate the new values into the * transition values. */ N += 1.0; Sx += newvalX; Sy += newvalY; if (transvalues[0] > 0.0) { tmpX = newvalX * N - Sx; tmpY = newvalY * N - Sy; scale = 1.0 / (N * transvalues[0]); Sxx += tmpX * tmpX * scale; Syy += tmpY * tmpY * scale; Sxy += tmpX * tmpY * scale; /* * Overflow check. We only report an overflow error when finite * inputs lead to infinite results. Note also that Sxx, Syy and Sxy * should be NaN if any of the relevant inputs are infinite, so we * intentionally prevent them from becoming infinite. */ if (isinf(Sx) || isinf(Sxx) || isinf(Sy) || isinf(Syy) || isinf(Sxy)) { if (((isinf(Sx) || isinf(Sxx)) && !isinf(transvalues[1]) && !isinf(newvalX)) || ((isinf(Sy) || isinf(Syy)) && !isinf(transvalues[3]) && !isinf(newvalY)) || (isinf(Sxy) && !isinf(transvalues[1]) && !isinf(newvalX) && !isinf(transvalues[3]) && !isinf(newvalY))) float_overflow_error(); if (isinf(Sxx)) Sxx = get_float8_nan(); if (isinf(Syy)) Syy = get_float8_nan(); if (isinf(Sxy)) Sxy = get_float8_nan(); } } else { /* * At the first input, we normally can leave Sxx et al as 0. However, * if the first input is Inf or NaN, we'd better force the dependent * sums to NaN; otherwise we will falsely report variance zero when * there are no more inputs. */ if (isnan(newvalX) || isinf(newvalX)) Sxx = Sxy = get_float8_nan(); if (isnan(newvalY) || isinf(newvalY)) Syy = Sxy = get_float8_nan(); } /* * 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] = Sx; transvalues[2] = Sxx; transvalues[3] = Sy; transvalues[4] = Syy; transvalues[5] = Sxy; PG_RETURN_ARRAYTYPE_P(transarray); } else { Datum transdatums[6]; ArrayType *result; transdatums[0] = Float8GetDatumFast(N); transdatums[1] = Float8GetDatumFast(Sx); transdatums[2] = Float8GetDatumFast(Sxx); transdatums[3] = Float8GetDatumFast(Sy); transdatums[4] = Float8GetDatumFast(Syy); transdatums[5] = Float8GetDatumFast(Sxy); result = construct_array_builtin(transdatums, 6, FLOAT8OID); 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 N1, Sx1, Sxx1, Sy1, Syy1, Sxy1, N2, Sx2, Sxx2, Sy2, Syy2, Sxy2, tmp1, tmp2, N, Sx, Sxx, Sy, Syy, Sxy; transvalues1 = check_float8_array(transarray1, "float8_regr_combine", 6); transvalues2 = check_float8_array(transarray2, "float8_regr_combine", 6); N1 = transvalues1[0]; Sx1 = transvalues1[1]; Sxx1 = transvalues1[2]; Sy1 = transvalues1[3]; Syy1 = transvalues1[4]; Sxy1 = transvalues1[5]; N2 = transvalues2[0]; Sx2 = transvalues2[1]; Sxx2 = transvalues2[2]; Sy2 = transvalues2[3]; Syy2 = transvalues2[4]; Sxy2 = transvalues2[5]; /*-------------------- * The transition values combine using a generalization of the * Youngs-Cramer algorithm as follows: * * N = N1 + N2 * Sx = Sx1 + Sx2 * Sxx = Sxx1 + Sxx2 + N1 * N2 * (Sx1/N1 - Sx2/N2)^2 / N * Sy = Sy1 + Sy2 * Syy = Syy1 + Syy2 + N1 * N2 * (Sy1/N1 - Sy2/N2)^2 / N * Sxy = Sxy1 + Sxy2 + N1 * N2 * (Sx1/N1 - Sx2/N2) * (Sy1/N1 - Sy2/N2) / N * * It's worth handling the special cases N1 = 0 and N2 = 0 separately * since those cases are trivial, and we then don't need to worry about * division-by-zero errors in the general case. *-------------------- */ if (N1 == 0.0) { N = N2; Sx = Sx2; Sxx = Sxx2; Sy = Sy2; Syy = Syy2; Sxy = Sxy2; } else if (N2 == 0.0) { N = N1; Sx = Sx1; Sxx = Sxx1; Sy = Sy1; Syy = Syy1; Sxy = Sxy1; } else { N = N1 + N2; Sx = float8_pl(Sx1, Sx2); tmp1 = Sx1 / N1 - Sx2 / N2; Sxx = Sxx1 + Sxx2 + N1 * N2 * tmp1 * tmp1 / N; if (unlikely(isinf(Sxx)) && !isinf(Sxx1) && !isinf(Sxx2)) float_overflow_error(); Sy = float8_pl(Sy1, Sy2); tmp2 = Sy1 / N1 - Sy2 / N2; Syy = Syy1 + Syy2 + N1 * N2 * tmp2 * tmp2 / N; if (unlikely(isinf(Syy)) && !isinf(Syy1) && !isinf(Syy2)) float_overflow_error(); Sxy = Sxy1 + Sxy2 + N1 * N2 * tmp1 * tmp2 / N; if (unlikely(isinf(Sxy)) && !isinf(Sxy1) && !isinf(Sxy2)) float_overflow_error(); } /* * 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)) { transvalues1[0] = N; transvalues1[1] = Sx; transvalues1[2] = Sxx; transvalues1[3] = Sy; transvalues1[4] = Syy; transvalues1[5] = Sxy; PG_RETURN_ARRAYTYPE_P(transarray1); } else { Datum transdatums[6]; ArrayType *result; transdatums[0] = Float8GetDatumFast(N); transdatums[1] = Float8GetDatumFast(Sx); transdatums[2] = Float8GetDatumFast(Sxx); transdatums[3] = Float8GetDatumFast(Sy); transdatums[4] = Float8GetDatumFast(Syy); transdatums[5] = Float8GetDatumFast(Sxy); result = construct_array_builtin(transdatums, 6, FLOAT8OID); PG_RETURN_ARRAYTYPE_P(result); } } Datum float8_regr_sxx(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 *transvalues; float8 N, Sxx; transvalues = check_float8_array(transarray, "float8_regr_sxx", 6); N = transvalues[0]; Sxx = transvalues[2]; /* if N is 0 we should return NULL */ if (N < 1.0) PG_RETURN_NULL(); /* Note that Sxx is guaranteed to be non-negative */ PG_RETURN_FLOAT8(Sxx); } Datum float8_regr_syy(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 *transvalues; float8 N, Syy; transvalues = check_float8_array(transarray, "float8_regr_syy", 6); N = transvalues[0]; Syy = transvalues[4]; /* if N is 0 we should return NULL */ if (N < 1.0) PG_RETURN_NULL(); /* Note that Syy is guaranteed to be non-negative */ PG_RETURN_FLOAT8(Syy); } Datum float8_regr_sxy(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 *transvalues; float8 N, Sxy; transvalues = check_float8_array(transarray, "float8_regr_sxy", 6); N = transvalues[0]; Sxy = transvalues[5]; /* if N is 0 we should return NULL */ if (N < 1.0) PG_RETURN_NULL(); /* A negative result is valid here */ PG_RETURN_FLOAT8(Sxy); } Datum float8_regr_avgx(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 *transvalues; float8 N, Sx; transvalues = check_float8_array(transarray, "float8_regr_avgx", 6); N = transvalues[0]; Sx = transvalues[1]; /* if N is 0 we should return NULL */ if (N < 1.0) PG_RETURN_NULL(); PG_RETURN_FLOAT8(Sx / N); } Datum float8_regr_avgy(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 *transvalues; float8 N, Sy; transvalues = check_float8_array(transarray, "float8_regr_avgy", 6); N = transvalues[0]; Sy = transvalues[3]; /* if N is 0 we should return NULL */ if (N < 1.0) PG_RETURN_NULL(); PG_RETURN_FLOAT8(Sy / N); } Datum float8_covar_pop(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 *transvalues; float8 N, Sxy; transvalues = check_float8_array(transarray, "float8_covar_pop", 6); N = transvalues[0]; Sxy = transvalues[5]; /* if N is 0 we should return NULL */ if (N < 1.0) PG_RETURN_NULL(); PG_RETURN_FLOAT8(Sxy / N); } Datum float8_covar_samp(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 *transvalues; float8 N, Sxy; transvalues = check_float8_array(transarray, "float8_covar_samp", 6); N = transvalues[0]; Sxy = transvalues[5]; /* if N is <= 1 we should return NULL */ if (N < 2.0) PG_RETURN_NULL(); PG_RETURN_FLOAT8(Sxy / (N - 1.0)); } Datum float8_corr(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 *transvalues; float8 N, Sxx, Syy, Sxy; transvalues = check_float8_array(transarray, "float8_corr", 6); N = transvalues[0]; Sxx = transvalues[2]; Syy = transvalues[4]; Sxy = transvalues[5]; /* if N is 0 we should return NULL */ if (N < 1.0) PG_RETURN_NULL(); /* Note that Sxx and Syy are guaranteed to be non-negative */ /* per spec, return NULL for horizontal and vertical lines */ if (Sxx == 0 || Syy == 0) PG_RETURN_NULL(); PG_RETURN_FLOAT8(Sxy / sqrt(Sxx * Syy)); } Datum float8_regr_r2(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 *transvalues; float8 N, Sxx, Syy, Sxy; transvalues = check_float8_array(transarray, "float8_regr_r2", 6); N = transvalues[0]; Sxx = transvalues[2]; Syy = transvalues[4]; Sxy = transvalues[5]; /* if N is 0 we should return NULL */ if (N < 1.0) PG_RETURN_NULL(); /* Note that Sxx and Syy are guaranteed to be non-negative */ /* per spec, return NULL for a vertical line */ if (Sxx == 0) PG_RETURN_NULL(); /* per spec, return 1.0 for a horizontal line */ if (Syy == 0) PG_RETURN_FLOAT8(1.0); PG_RETURN_FLOAT8((Sxy * Sxy) / (Sxx * Syy)); } Datum float8_regr_slope(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 *transvalues; float8 N, Sxx, Sxy; transvalues = check_float8_array(transarray, "float8_regr_slope", 6); N = transvalues[0]; Sxx = transvalues[2]; Sxy = transvalues[5]; /* if N is 0 we should return NULL */ if (N < 1.0) PG_RETURN_NULL(); /* Note that Sxx is guaranteed to be non-negative */ /* per spec, return NULL for a vertical line */ if (Sxx == 0) PG_RETURN_NULL(); PG_RETURN_FLOAT8(Sxy / Sxx); } Datum float8_regr_intercept(PG_FUNCTION_ARGS) { ArrayType *transarray = PG_GETARG_ARRAYTYPE_P(0); float8 *transvalues; float8 N, Sx, Sxx, Sy, Sxy; transvalues = check_float8_array(transarray, "float8_regr_intercept", 6); N = transvalues[0]; Sx = transvalues[1]; Sxx = transvalues[2]; Sy = transvalues[3]; Sxy = transvalues[5]; /* if N is 0 we should return NULL */ if (N < 1.0) PG_RETURN_NULL(); /* Note that Sxx is guaranteed to be non-negative */ /* per spec, return NULL for a vertical line */ if (Sxx == 0) PG_RETURN_NULL(); PG_RETURN_FLOAT8((Sy - Sx * Sxy / Sxx) / N); } /* * ==================================== * 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); PG_RETURN_FLOAT8(float8_pl((float8) arg1, arg2)); } Datum float48mi(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_FLOAT8(float8_mi((float8) arg1, arg2)); } Datum float48mul(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_FLOAT8(float8_mul((float8) arg1, arg2)); } Datum float48div(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_FLOAT8(float8_div((float8) arg1, arg2)); } /* * float84pl - returns arg1 + arg2 * float84mi - returns arg1 - arg2 * float84mul - returns arg1 * arg2 * float84div - returns arg1 / arg2 */ Datum float84pl(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_FLOAT8(float8_pl(arg1, (float8) arg2)); } Datum float84mi(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_FLOAT8(float8_mi(arg1, (float8) arg2)); } Datum float84mul(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_FLOAT8(float8_mul(arg1, (float8) arg2)); } Datum float84div(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_FLOAT8(float8_div(arg1, (float8) arg2)); } /* * ==================== * COMPARISON OPERATORS * ==================== */ /* * float48{eq,ne,lt,le,gt,ge} - float4/float8 comparison operations */ Datum float48eq(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_BOOL(float8_eq((float8) arg1, arg2)); } Datum float48ne(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_BOOL(float8_ne((float8) arg1, arg2)); } Datum float48lt(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_BOOL(float8_lt((float8) arg1, arg2)); } Datum float48le(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_BOOL(float8_le((float8) arg1, arg2)); } Datum float48gt(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_BOOL(float8_gt((float8) arg1, arg2)); } Datum float48ge(PG_FUNCTION_ARGS) { float4 arg1 = PG_GETARG_FLOAT4(0); float8 arg2 = PG_GETARG_FLOAT8(1); PG_RETURN_BOOL(float8_ge((float8) arg1, arg2)); } /* * float84{eq,ne,lt,le,gt,ge} - float8/float4 comparison operations */ Datum float84eq(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_BOOL(float8_eq(arg1, (float8) arg2)); } Datum float84ne(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_BOOL(float8_ne(arg1, (float8) arg2)); } Datum float84lt(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_BOOL(float8_lt(arg1, (float8) arg2)); } Datum float84le(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_BOOL(float8_le(arg1, (float8) arg2)); } Datum float84gt(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_BOOL(float8_gt(arg1, (float8) arg2)); } Datum float84ge(PG_FUNCTION_ARGS) { float8 arg1 = PG_GETARG_FLOAT8(0); float4 arg2 = PG_GETARG_FLOAT4(1); PG_RETURN_BOOL(float8_ge(arg1, (float8) arg2)); } /* * Implements the float8 version of the width_bucket() function * defined by SQL2003. See also width_bucket_numeric(). * * 'bound1' and 'bound2' are the lower and upper bounds of the * histogram's range, respectively. 'count' is the number of buckets * in the histogram. width_bucket() returns an integer indicating the * bucket number that 'operand' belongs to in an equiwidth histogram * with the specified characteristics. An operand smaller than the * lower bound is assigned to bucket 0. An operand greater than or equal * to 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) 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) { if (pg_add_s32_overflow(count, 1, &result)) ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("integer out of range"))); } else { if (!isinf(bound2 - bound1)) { /* The quotient is surely in [0,1], so this can't overflow */ result = count * ((operand - bound1) / (bound2 - bound1)); } else { /* * We get here if bound2 - bound1 overflows DBL_MAX. Since * both bounds are finite, their difference can't exceed twice * DBL_MAX; so we can perform the computation without overflow * by dividing all the inputs by 2. That should be exact too, * except in the case where a very small operand underflows to * zero, which would have negligible impact on the result * given such large bounds. */ result = count * ((operand / 2 - bound1 / 2) / (bound2 / 2 - bound1 / 2)); } /* The quotient could round to 1.0, which would be a lie */ if (result >= count) result = count - 1; /* Having done that, we can add 1 without fear of overflow */ result++; } } else if (bound1 > bound2) { if (operand > bound1) result = 0; else if (operand <= bound2) { if (pg_add_s32_overflow(count, 1, &result)) ereport(ERROR, (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), errmsg("integer out of range"))); } else { if (!isinf(bound1 - bound2)) result = count * ((bound1 - operand) / (bound1 - bound2)); else result = count * ((bound1 / 2 - operand / 2) / (bound1 / 2 - bound2 / 2)); if (result >= count) result = count - 1; result++; } } 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); }