diff options
Diffstat (limited to 'src/backend/utils/adt/float.c')
| -rw-r--r-- | src/backend/utils/adt/float.c | 4080 |
1 files changed, 4080 insertions, 0 deletions
diff --git a/src/backend/utils/adt/float.c b/src/backend/utils/adt/float.c new file mode 100644 index 0000000..098bbb3 --- /dev/null +++ b/src/backend/utils/adt/float.c @@ -0,0 +1,4080 @@ +/*------------------------------------------------------------------------- + * + * float.c + * Functions for the built-in floating-point types. + * + * Portions Copyright (c) 1996-2021, 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 "common/int.h" +#include "common/shortest_dec.h" +#include "libpq/pqformat.h" +#include "miscadmin.h" +#include "utils/array.h" +#include "utils/float.h" +#include "utils/fmgrprotos.h" +#include "utils/sortsupport.h" +#include "utils/timestamp.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(). + */ +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; + +/* State for drandom() and setseed() */ +static bool drandom_seed_set = false; +static unsigned short drandom_seed[3] = {0, 0, 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); + char *orig_num; + float 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 %s: \"%s\"", + "real", orig_num))); + + 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. + * + * Use isinf() rather than HUGE_VALF on VS2013 because it + * generates a spurious overflow warning for -HUGE_VALF. Also use + * isinf() if HUGE_VALF is missing. + */ + if (val == 0.0 || +#if !defined(HUGE_VALF) || (defined(_MSC_VER) && (_MSC_VER < 1900)) + isinf(val) +#else + (val >= HUGE_VALF || val <= -HUGE_VALF) +#endif + ) + 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 %s: \"%s\"", + "real", orig_num))); + } + + /* 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 %s: \"%s\"", + "real", orig_num))); + + PG_RETURN_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(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)); +} + +/* Convenience macro: set *have_error flag (if provided) or throw error */ +#define RETURN_ERROR(throw_error, have_error) \ +do { \ + if (have_error) { \ + *have_error = true; \ + return 0.0; \ + } else { \ + throw_error; \ + } \ +} while (0) + +/* + * float8in_internal_opt_error - 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. + * + * When "*have_error" flag is provided, it's set instead of throwing an + * error. This is helpful when caller need to handle errors by itself. + */ +double +float8in_internal_opt_error(char *num, char **endptr_p, + const char *type_name, const char *orig_string, + bool *have_error) +{ + double val; + char *endptr; + + if (have_error) + *have_error = false; + + /* 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') + RETURN_ERROR(ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type %s: \"%s\"", + type_name, orig_string))), + have_error); + + 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'; + RETURN_ERROR(ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("\"%s\" is out of range for type double precision", + errnumber))), + have_error); + } + } + else + RETURN_ERROR(ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type " + "%s: \"%s\"", + type_name, orig_string))), + have_error); + } + + /* 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') + RETURN_ERROR(ereport(ERROR, + (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), + errmsg("invalid input syntax for type " + "%s: \"%s\"", + type_name, orig_string))), + have_error); + + return val; +} + +/* + * Interface to float8in_internal_opt_error() without "have_error" argument. + */ +double +float8in_internal(char *num, char **endptr_p, + const char *type_name, const char *orig_string) +{ + return float8in_internal_opt_error(num, endptr_p, type_name, + orig_string, NULL); +} + + +/* + * 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((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_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); +} + + +/* + * drandom - returns a random number + */ +Datum +drandom(PG_FUNCTION_ARGS) +{ + float8 result; + + /* Initialize random seed, if not done yet in this process */ + if (unlikely(!drandom_seed_set)) + { + /* + * If possible, initialize the seed using high-quality random bits. + * Should that fail for some reason, we fall back on a lower-quality + * seed based on current time and PID. + */ + if (!pg_strong_random(drandom_seed, sizeof(drandom_seed))) + { + TimestampTz now = GetCurrentTimestamp(); + uint64 iseed; + + /* Mix the PID with the most predictable bits of the timestamp */ + iseed = (uint64) now ^ ((uint64) MyProcPid << 32); + drandom_seed[0] = (unsigned short) iseed; + drandom_seed[1] = (unsigned short) (iseed >> 16); + drandom_seed[2] = (unsigned short) (iseed >> 32); + } + drandom_seed_set = true; + } + + /* pg_erand48 produces desired result range [0.0 - 1.0) */ + result = pg_erand48(drandom_seed); + + PG_RETURN_FLOAT8(result); +} + + +/* + * setseed - set seed for the random number generator + */ +Datum +setseed(PG_FUNCTION_ARGS) +{ + float8 seed = PG_GETARG_FLOAT8(0); + uint64 iseed; + + if (seed < -1 || seed > 1 || isnan(seed)) + ereport(ERROR, + (errcode(ERRCODE_INVALID_PARAMETER_VALUE), + errmsg("setseed parameter %g is out of allowed range [-1,1]", + seed))); + + /* Use sign bit + 47 fractional bits to fill drandom_seed[] */ + iseed = (int64) (seed * (float8) UINT64CONST(0x7FFFFFFFFFFF)); + drandom_seed[0] = (unsigned short) iseed; + drandom_seed[1] = (unsigned short) (iseed >> 16); + drandom_seed[2] = (unsigned short) (iseed >> 32); + drandom_seed_set = true; + + 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 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(transdatums, 3, + FLOAT8OID, + sizeof(float8), FLOAT8PASSBYVAL, TYPALIGN_DOUBLE); + + 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(transdatums, 3, + FLOAT8OID, + sizeof(float8), FLOAT8PASSBYVAL, TYPALIGN_DOUBLE); + + 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(transdatums, 3, + FLOAT8OID, + sizeof(float8), FLOAT8PASSBYVAL, TYPALIGN_DOUBLE); + + 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(transdatums, 6, + FLOAT8OID, + sizeof(float8), FLOAT8PASSBYVAL, TYPALIGN_DOUBLE); + + 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(transdatums, 6, + FLOAT8OID, + sizeof(float8), FLOAT8PASSBYVAL, TYPALIGN_DOUBLE); + + 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 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) + { + if (pg_add_s32_overflow(count, 1, &result)) + 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) + { + if (pg_add_s32_overflow(count, 1, &result)) + 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); +} |